.. _fermi-dirac-notes: ======================= Fermi-Dirac Integrals ======================= The :mod:`mmf.math.special.fermi` module provides efficient ways of calculating the following generalized Fermi-Dirac integrals: .. math:: \begin{aligned} \eta &= \beta(\mu-m), & \theta &= \frac{1}{\beta m}, \end{aligned}\\ \begin{aligned} F_{k}(\eta, \theta) &= \int_{0}^{\infty} \d{t}\; \frac{t^{k}\sqrt{1+ \theta t/2}} {1+e^{t-\eta}},\\ F^{(a)}_{k}(\beta\mu, \beta m) &= F_{k}(\eta, \theta) - F_{k}(\eta - \beta\mu, \theta),\\ &= \int_{0}^{\infty} \d{t}\; \frac{t^{k}\sqrt{1+ \theta t/2}\sinh(\eta + \theta^{-1})} {\cosh(t + \theta^{-1}) + \cosh(\eta + \theta^{-1})} \end{aligned} Following [Natarajan:1993rz]_, we note that the key to evaluating these is to use the analytic properties. The integrands are well behaved except near the poles at $t = \beta(\pm \mu - m) \pm i\pi$. First we note the asymptotic behaviour of the integrand: .. math:: \frac{t^{k}\sqrt{1+ \theta t/2}\sinh(\eta + \theta^{-1})} {\cosh(t + \theta^{-1}) + \cosh(\eta + \theta^{-1})} \rightarrow \sqrt{2\theta}\sinh(\eta + \theta^{-1})e^{-\theta^{-1}}t^{k+1/2}e^{-t}. For $\tau \gg 1$: .. math:: \int_{\tau}^{\infty} t^{a}e^{-t} \d{t} \sim \tau^{a} e^{-\tau} Thus, if we want the integral accurate to an absolute accuracy of $\varepsilon$, then we can neglect the tail for .. math:: t^{k+1/2}e^{-t} < \frac{\varepsilon e^{\theta^{-1}}} {\sqrt{2\theta}\sinh(\eta + \theta^{-1})} \approx \sqrt{\frac{2}{\theta}}\epsilon e^{-\eta},\\ t > -\left(k+\tfrac{1}{2}\right) \operatorname{W}_{-1}\frac{-c}{k + \tfrac{1}{2}} The integrands have the form .. math:: ct^{k}\sqrt{1 + \tfrac{\theta}{2}t} \left( \frac{\cdots} {\left[\cosh(t + \beta m) + \cosh(\eta + \beta m)\right]^b} \right) The integrand is peaked about $t_0$: .. math:: t_0 \sim k + \tfrac{1}{2} + {\rm LambertW}\left( (k + \tfrac{1}{2})e^{\eta - k - \tfrac{1}{2}} \right) \approx k + \tfrac{1}{2} + \eta Here is a plot of various integrands .. plot:: :include-source: import mmf.math def f(t, k, eta, theta, b, c): beta_m = 1.0/theta beta_mu = eta + beta_m if 1 == b: if 0 == c: num = np.sinh(beta_mu) else: num = np.cosh(beta_mu) + np.exp(-t-beta_m) elif 2 == b: if 0 == c: num = np.sinh(t + beta_m)*np.sinh(beta_mu) else: num = 1 + np.cosh(beta_mu)*np.cosh(t+beta_m) res = (t**k*np.sqrt(1 + t/beta_m/2)* num/(np.cosh(t+beta_m) + np.cosh(beta_mu))**b) return res t_0 = np.linspace(0,5,100) for eta in [0.01,0.5,1.0,2.0,5,100]: for theta in [0.01, 1.0, 100.0]: for b in [1,2]: for c in [0,1]: for k in [0.5,1.5,2.5]: a = k + 0.5 t0 = a + mmf.math.LambertW(a*np.exp(eta - a)) t1 = a + eta f0 = a*t0**(k-1)*np.sqrt(1 + t0*theta/2)*np.exp(eta - t0) f0 = f(t0, k, eta, theta, b, c) plt.plot(t_0, f(t_0*t0, k, eta, theta, b, c)/f0, '-b') plt.plot(t_0, f(t_0*t1, k, eta, theta, b, c)/f0, ':y') plt.xlabel("t/t_0") plt.ylabel("integrand/integrand(t_0)") As shown in [Schwartz:1969if]_, the truncation error of the Trapezoidal method is of order $\exp(-2\pi w/h)$ where $w$ is the distance to the nearest pole and $h$ is the step-size. The poles in the integrand are at must have a step-size of $h \lesssim -\ln(\epsilon)/2\pi d$ As suggested in [Natarajan:1993rz]_, given that the pole is located off the real axis about $\eta$, we break the interval into two regions $[0,\eta]$ and $[\eta, 60\eta]$ and use the IMT method which approximates the integral .. math:: \int_{0}^{1}\d{x}\; f(x) \approx \sum_{k=1}^{N-1} f(x_k)x'_{k} where .. math:: x_k = \int_{0}^{t_{k}} x'(t)\d{t},\\ t_k = \frac{k}{N},\\ x'(t) = \frac{1}{Q}\exp\left(-\frac{1}{t} - \frac{1}{1-t}\right),\\ Q = \int_{0}^{1}\exp\left(-\frac{1}{t} - \frac{1}{1-t}\right). This is known to improve the analytic properties at the endpoints. We using successively more points until the desired accuracy is achieved. Here are some examples of the transformed functions as seen by the trapezoidal method. The yellow curves are for the interval $[0, k+ \tfrac{1}{2} + \eta]$ while the red curves are for the region $[k+ \tfrac{1}{2} + \eta, 60\eta]$. .. plot:: :include-source: import mmf.math.integrate.integrate_1d.imt imt = mmf.math.integrate.integrate_1d.imt.IMT() def f(t, k, eta, theta, b, c): beta_m = 1.0/theta beta_mu = eta + beta_m if 1 == b: if 0 == c: num = np.sinh(beta_mu) else: num = np.cosh(beta_mu) + np.exp(-t-beta_m) elif 2 == b: if 0 == c: num = np.sinh(t + beta_m)*np.sinh(beta_mu) else: num = 1 + np.cosh(beta_mu)*np.cosh(t+beta_m) res = (t**k*np.sqrt(1 + t/beta_m/2)* num/(np.cosh(t+beta_m) + np.cosh(beta_mu))**b) return res t = np.linspace(0,1,50) for b in [1, 2]: for c in [0,1]: plt.subplot(2,2,b+2*c) if b == 1 and c == 0: plt.title("n") if b == 1 and c == 1: plt.title("n_s") if b == 2 and c == 0: plt.title("dn_dm") if b == 2 and c == 1: plt.title("dn_dmu") for eta in [0.01,1.0,100]: for theta in [0.01, 1.0, 100.0]: for k in [0.5,1.5,2.5]: a = k + 0.5 #t0 = a + mmf.math.LambertW(a*np.exp(eta - a)) t0 = a + eta f_ = lambda x: f(x, k, eta, theta, b, c) phi0 = imt.phi(f_, 0, t0)(t) phi1 = imt.phi(f_, t0, 60*eta)(t) print t0/eta phi_max = max(phi0.max(), phi1.max()) plt.plot(t, phi0/phi_max, '--y') plt.plot(t, phi1/phi_max, '--r') trapezoidal rule on the transformed .. [Natarajan:1993rz] A. Natarajan and N. Mohan Kumar, `"On the numerical evaluation of the generalised Fermi-Dirac integrals" `_, Computer Physics Communications 76(1), 48-50 (1993) .. [Schwartz:1969if] Charles Schwartz, `"Numerical integration of analytic functions" `_, Journal of Computational Physics 4(1), 19--29 (1969) Physical Quantities =================== Here we present the relevant physical quantities for a single species of free fermion with mass $m$ and chemical potential $\mu$, and dispersion $E_k = \sqrt{k^2 + m^2}$. Everything is derived from the pressure: .. math:: \begin{aligned} P &= \frac{1}{\beta}\int\dbar^{3}{\vect{k}}\;\left[ \ln\left(1+e^{-\beta(E_k-\mu)}\right) + \ln\left(1+e^{-\beta(E_k+\mu)}\right) \right],\\ &= \frac{\sqrt{2}m^4}{3\pi^2}\theta^{5/2}\left[ F_{3/2} + \frac{\theta}{2}F_{5/2}\right] \end{aligned} The second term is the antiparticle contribution which may or may not be included. The standard Fermi-Dirac integrals do not include it, but we will since it is almost always more physical and inexpensive to include. The relevant physical quantities follow from differentiating: .. math:: \begin{aligned} n &= \pdiff{P}{\mu} \geq 0,\\ n^{(s)} &= -\pdiff{P}{m} \geq 0,\\ n_{,m} &= \pdiff{n}{m} = -n^{(s)}_{,\mu} = -\pdiff{n^{(s)}}{\mu} = \frac{\partial^{2}P}{\partial\mu \partial m},\\ n_{,\mu} &= \pdiff{n}{\mu} = \pdiff[2]{P}{\mu} \geq 0,\\ n^{(s)}_{,m} &= \pdiff{n^{(s)}}{m} = -\pdiff[2]{P}{m}. \end{aligned} Here are the expressions: .. math:: \begin{aligned} P &= \frac{(2\beta m)^{3/2}}{6\pi^2 \beta^4} \int_{0}^{\infty}\d{t}\;\sqrt{1 + \tfrac{\theta}{2}t} \left(t^{3/2} + \tfrac{\theta}{2} t^{5/2}\right) \left(\frac{\cosh(\beta\mu) + e^{-(t+\beta m)}} {\cosh(t + \beta m) + \cosh(\eta + \beta m)} \right),\\ n &= \frac{(\beta m)^{3/2}}{\sqrt{2}\pi^2 \beta^3} \int_{0}^{\infty}\d{t}\;\sqrt{1 + \tfrac{\theta}{2}t} \left(t^{1/2} + \theta t^{3/2}\right) \left(\frac{\sinh(\beta\mu)} {\cosh(t + \beta m) + \cosh(\eta + \beta m)} \right),\\ n^{(s)} &= \frac{(\beta m)^{3/2}}{\sqrt{2}\pi^2 \beta^3} \int_{0}^{\infty}\d{t}\;\sqrt{1 + \tfrac{\theta}{2}t} t^{1/2} \left(\frac{\cosh(\beta\mu) + e^{-(t+\beta m)}} {\cosh(t + \beta m) + \cosh(\eta + \beta m)} \right),\\ n_{,\mu} &= \frac{(\beta m)^{3/2}}{\sqrt{2}\pi^2 \beta^2} \int_{0}^{\infty}\d{t}\;\sqrt{1 + \tfrac{\theta}{2}t} \left(t^{1/2} + \theta t^{3/2}\right) \left(\frac{1 + \cosh(t + \beta m)\cosh(\beta\mu)} {\left[\cosh(t + \beta m) + \cosh(\eta + \beta m)\right]^2} \right),\\ n^{(s)}_{,\mu} = - n_{,m} &= \frac{(\beta m)^{3/2}}{\sqrt{2}\pi^2 \beta^2} \int_{0}^{\infty}\d{t}\;\sqrt{1 + \tfrac{\theta}{2}t} t^{1/2} \left(\frac{\sinh(t + \beta m)\sinh(\beta\mu)} {\left[\cosh(t + \beta m) + \cosh(\eta + \beta m)\right]^2} \right),\\ \end{aligned} In each case, the behaviour of the integrand is dominated by the denominator which is minimized at $t=\eta$. Here are some missing intermediate steps: .. math:: \begin{aligned} n_{\beta}(\mu, m) &= \int \dbar^{3}{k}\; \frac{1}{1+e^{\beta(\sqrt{k^2 + m^2} - \mu)}} = \frac{1}{2\pi^2}\int_{0}^{\infty}\d{k}\; \frac{1}{1+e^{\beta(\sqrt{k^2 + m^2} - \mu)}} = \left(\frac{m}{2\beta}\right)^{3/2}\frac{2}{\pi^2} \int_{0}^{\infty}\d{t}\; \frac{E}{m}\frac{\sqrt{t}\sqrt{1+\theta t/2}} {1+e^{t-\eta}}. \end{aligned}