fermion.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 #ifndef O2SCL_FERMION_H
24 #define O2SCL_FERMION_H
25 
26 /** \file fermion.h
27  \brief File defining \ref o2scl::fermion
28 */
29 #include <string>
30 #include <iostream>
31 #include <fstream>
32 #include <cmath>
33 
34 // For gsl_sf_fermi_dirac_int()
35 #include <gsl/gsl_specfunc.h>
36 
37 #include <o2scl/constants.h>
38 #include <o2scl/funct.h>
39 #include <o2scl/root.h>
40 #include <o2scl/root_cern.h>
41 #include <o2scl/part.h>
42 
43 #ifndef DOXYGEN_NO_O2NS
44 namespace o2scl {
45 #endif
46 
47  /** \brief Fermion class
48 
49  This class adds two member data variables, \ref kf and \ref
50  del, for the Fermi momentum and the gap, respectively.
51  */
52  class fermion : public part {
53 
54  public:
55 
56  /// Fermi momentum
57  double kf;
58  /// Gap
59  double del;
60 
61  /// Create a fermion with mass \c mass and degeneracy \c dof.
62  fermion(double mass=0, double dof=0);
63 
64  virtual ~fermion() {
65  }
66 
67  /// Return string denoting type ("fermion")
68  virtual const char *type() { return "fermion"; }
69 
70  /// Copy constructor
71  fermion(const fermion &f) {
72  g=f.g;
73  m=f.m;
74  ms=f.ms;
75  n=f.n;
76  ed=f.ed;
77  pr=f.pr;
78  mu=f.mu;
79  en=f.en;
80  nu=f.nu;
83  kf=f.kf;
84  del=f.del;
85  }
86 
87  /// Copy construction with operator=()
88  fermion &operator=(const fermion &f) {
89  if (this!=&f) {
90  g=f.g;
91  m=f.m;
92  ms=f.ms;
93  n=f.n;
94  ed=f.ed;
95  pr=f.pr;
96  mu=f.mu;
97  en=f.en;
98  nu=f.nu;
101  kf=f.kf;
102  del=f.del;
103  }
104  return *this;
105  }
106 
107  };
108 
109  /** \brief Fermion properties at zero temperature
110 
111  This is a base class for the computation of fermionic statistics
112  at zero temperature. The more general case of finite temperature
113  is taken care of by \ref fermion_eval_thermo objects. The
114  primary functions are calc_mu_zerot() and calc_density_zerot()
115  which compute all the thermodynamic quantities as a function of
116  the chemical potential, or the density, respectively.
117 
118  \future Use hypot() and other more accurate functions for the
119  analytic expressions for the zero temperature integrals. [Progress
120  has been made, but there are probably other functions which may
121  break down for small but finite masses and temperatures]
122  */
124 
125  public:
126 
127  fermion_zerot() {
128  };
129 
130  virtual ~fermion_zerot() {
131  }
132 
133  /// \name Zero-temperature fermions
134  //@{
135  /** \brief Calculate the Fermi momentum from the density
136 
137  Uses the relation \f$ k_F = ( 6 \pi^2 n /g )^{1/3} \f$
138  */
139  void kf_from_density(fermion &f);
140 
141  /** \brief Energy density at T=0 from \ref fermion::kf and \ref part::ms
142 
143  Calculates the integral
144  \f[
145  \varepsilon = \frac{g}{2 \pi^2} \int_0^{k_F} k^2
146  \sqrt{k^2+m^{* 2}} d k
147  \f]
148  */
149  void energy_density_zerot(fermion &f);
150 
151  /** \brief Pressure at T=0 from \ref fermion::kf and \ref part::ms
152 
153  Calculates the integral
154  \f[
155  P=\frac{g}{6 \pi^2} \int_0^{k_F} \frac{k^4}{\sqrt{k^2+m^{* 2}}} d k
156  \f]
157  */
158  void pressure_zerot(fermion &f);
159 
160  /** \brief Zero temperature fermions from \ref part::mu or
161  \ref part::nu and \ref part::ms
162  */
163  virtual void calc_mu_zerot(fermion &f);
164 
165  /** \brief Zero temperature fermions from \ref part::n and \ref part::ms
166  */
167  virtual void calc_density_zerot(fermion &f);
168  //@}
169 
170  };
171 
172  /** \brief Fermion with finite-temperature thermodynamics
173  [abstract base]
174 
175  This is an abstract base for the computation of
176  finite-temperature fermionic statistics. Different children
177  (e.g. \ref fermion_eff and \ref fermion_rel) use different
178  techniques to computing the momentum integrations.
179 
180  Because massless fermions at finite temperature are much
181  simpler, there are separate member functions included in this
182  class to handle them. The functions massless_calc_density() and
183  massless_calc_mu() compute the thermodynamics of massless
184  fermions at finite temperature given the density or the chemical
185  potentials. The functions massless_pair_density() and
186  massless_pair_mu() perform the same task, but automatically
187  include antiparticles.
188 
189  The function massless_calc_density() uses a \ref root object to
190  solve for the chemical potential as a function of the density.
191  The default is an object of type root_cern. The function
192  massless_pair_density() does not need to use the \ref root
193  object because of the simplification afforded by the inclusion
194  of antiparticles.
195 
196  \future Create a Chebyshev approximation for inverting the
197  the Fermi functions for massless_calc_density() functions?
198  */
200 
201  public:
202 
204 
205  virtual ~fermion_eval_thermo() {
206  }
207 
208  /** \brief Non-degenerate expansion for fermions
209 
210  Attempts to evaulate thermodynamics of a non-degenerate
211  fermion. If the result is accurate to within the requested
212  precision, this function returns <tt>true</tt>, and otherwise
213  this function returns <tt>false</tt> and the values in stored
214  in the <tt>pr</tt>, <tt>n</tt>, <tt>en</tt>, and <tt>ed</tt>
215  field are meaningless.
216 
217  If \f$ \mu \f$ is negative and sufficiently far from zero,
218  then the thermodynamic quantities are smaller than the smallest
219  representable double-precision number. In this case,
220  this function will return <tt>true</tt> and report all
221  quantities as zero.
222 
223  Defining \f$ \psi \equiv (\mu-m)/T \f$, \f$ t \equiv T/m \f$,
224  and \f$ d \equiv g~m^4/(2 \pi^2) \f$ the pressure
225  in the non-degenerate limit (\f$ \psi \rightarrow - \infty \f$)
226  is (\ref Johns96)
227  \f[
228  P = d \sum_{n=1}^{\infty} P_n
229  \f]
230  where
231  \f[
232  P_n \equiv \left(-1\right)^{n+1} \left(\frac{t^2}{n^2}\right)
233  e^{n \left(\psi+1/t\right)} K_2 \left( \frac{n}{t} \right)
234  \f]
235  The density is then
236  \f[
237  n = d \sum_{n=1}^{\infty} \frac{n P_n}{T}
238  \f]
239  and the entropy density is
240  \f[
241  s = \frac{d}{m} \sum_{n=1}^{\infty} \left\{ \frac{2 P_n}{t}
242  -\frac{n P_n}{t^2}+
243  \frac{\left(-1\right)^{n+1}}{2 n}
244  e^{n \left(\psi+1/t\right)} \left[ K_1 \left( \frac{n}{t}
245  \right)+K_3 \left( \frac{n}{t} \right) \right]
246  \right\}
247  \f]
248 
249  This function is accurate over a wide range of conditions
250  when \f$ \psi < -4 \f$.
251 
252  The ratio of the nth term to the first term in the pressure
253  series is
254  \f[
255  R_n \equiv \frac{P_{n}}{P_{1}} = \frac{(-1)^{n+1}
256  e^{(n-1)(\psi+1/t)} K_2(n/t) }{n^2 K_2(1/t)}
257  \f]
258  This function currently uses 20 terms in the series and
259  immediately returns <tt>false</tt> if \f$ |R_{20}| \f$
260  is greater than <tt>prec</tt>
261 
262  In the nondegenerate and nonrelativistic (\f$ t \rightarrow 0
263  \f$) limit, the argument to the Bessel functions and the
264  exponential becomes too large. In this case, it's better to
265  use the expansions, e.g. for \f$ x \equiv n/t \rightarrow
266  \infty \f$,
267  \f[
268  \sqrt{\frac{2 x}{\pi}} e^{x} K_2(x) \approx
269  1 + \frac{3}{8 x} - \frac{15}{128 x^2} + ...
270  \f]
271  The current code currently goes up to \f$ x^{-12} \f$ in the
272  expansion, which is enough for the default precision of \f$
273  10^{-18} \f$ since \f$ (20/700)^{12} \sim 10^{-19} \f$.
274  */
275  virtual bool calc_mu_ndeg(fermion &f, double temper,
276  double prec=1.0e-18, bool inc_antip=false);
277 
278  /** \brief Degenerate expansion for fermions
279 
280  Attempts to evaulate thermodynamics of a degenerate fermion.
281  If the result is accurate to within the requested precision,
282  this function returns <tt>true</tt>, and otherwise this
283  function returns <tt>false</tt> and the values in stored in
284  the <tt>pr</tt>, <tt>n</tt>, <tt>en</tt>, and <tt>ed</tt>
285  field are meaningless.
286 
287  The pressure, density, and energy density, should be accurate
288  to the requested precision, but the first term in the series
289  expansion for the entropy is zero, so the entropy is one order
290  lower in accuracy.
291 
292  \future Make a function like this for dndm, dsdT, etc.
293  for fermion_deriv .
294  */
295  virtual bool calc_mu_deg(fermion &f, double temper,
296  double prec=1.0e-18);
297 
298  /** \brief Calculate properties as function of chemical potential
299  */
300  virtual void calc_mu(fermion &f, double temper)=0;
301 
302  /** \brief Calculate properties as function of density
303 
304  \note This function returns an integer value, in contrast to
305  \ref calc_mu(), because of the potential for non-convergence.
306  */
307  virtual int calc_density(fermion &f, double temper)=0;
308 
309  /** \brief Calculate properties with antiparticles as function of
310  chemical potential
311  */
312  virtual void pair_mu(fermion &f, double temper)=0;
313 
314  /** \brief Calculate properties with antiparticles as function of
315  density
316 
317  \note This function returns an integer value, in contrast to
318  \ref pair_mu(), because of the potential for non-convergence.
319  */
320  virtual int pair_density(fermion &f, double temper)=0;
321 
322  /// \name Massless fermions
323  //@{
324  /// Finite temperature massless fermions
325  virtual void massless_calc_mu(fermion &f, double temper);
326 
327  /// Finite temperature massless fermions
328  virtual void massless_calc_density(fermion &f, double temper);
329 
330  /** \brief Finite temperature massless fermions and antifermions
331  */
332  virtual void massless_pair_mu(fermion &f, double temper);
333 
334  /** \brief Finite temperature massless fermions and antifermions
335 
336  In the cases \f$ n^3 \gg T \f$ and \f$ T \gg n^3 \f$ ,
337  expansions are used instead of the exact formulas to avoid
338  loss of precision.
339 
340  In particular, using the parameter
341  \f[
342  \alpha = \frac{g^2 \pi^2 T^6}{243 n^2}
343  \f]
344  and defining the expression
345  \f[
346  \mathrm{cbt} = \alpha^{-1/6} \left( -1 + \sqrt{1+\alpha}\right)^{1/3}
347  \f]
348  we can write the chemical potential as
349  \f[
350  \mu = \frac{\pi T}{\sqrt{3}} \left(\frac{1}{\mathrm{cbt}} -
351  \mathrm{cbt} \right)
352  \f]
353 
354  These expressions, however, do not work well when \f$ \alpha
355  \f$ is very large or very small, so series expansions are
356  used whenever \f$ \alpha > 10^{4} \f$ or
357  \f$ \alpha < 3 \times 10^{-4} \f$. For small \f$ \alpha \f$,
358  \f[
359  \left(\frac{1}{\mathrm{cbt}} -
360  \mathrm{cbt} \right) \approx
361  \frac{2^{1/3}}{\alpha^{1/6}} -
362  \frac{\alpha^{1/6}}{2^{1/3}} +
363  \frac{\alpha^{5/6}}{6{\cdot}2^{2/3}} +
364  \frac{\alpha^{7/6}}{12{\cdot}2^{1/3}} -
365  \frac{\alpha^{11/6}}{18{\cdot}2^{2/3}} -
366  \frac{5 \alpha^{13/6}}{144{\cdot}2^{1/3}} +
367  \frac{77 \alpha^{17/6}}{2592{\cdot}2^{2/3}}
368  \f]
369  and for large \f$ \alpha \f$,
370  \f[
371  \left(\frac{1}{\mathrm{cbt}} -
372  \mathrm{cbt} \right) \approx
373  \frac{2}{3} \sqrt{\frac{1}{\alpha}} -
374  \frac{8}{81} \left(\frac{1}{\alpha}\right)^{3/2} +
375  \frac{32}{729} \left(\frac{1}{\alpha}\right)^{5/2}
376  \f]
377 
378  This approach works to within about 1 \part in \f$ 10^{12} \f$,
379  and is tested in <tt>fermion_ts.cpp</tt>.
380 
381  \future This could be improved by including more terms
382  in the expansions.
383  */
384  virtual void massless_pair_density(fermion &f, double temper);
385  //@}
386 
387  /** \brief Set the solver for use in massless_calc_density() */
389  massless_root=&rp;
390  return;
391  }
392 
393  /** \brief The default solver for massless_calc_density()
394 
395  We default to a solver of type root_cern here since we
396  don't have a bracket or a derivative.
397  */
399 
400  /// Return string denoting type ("fermion_eval_thermo")
401  virtual const char *type() { return "fermion_eval_thermo"; }
402 
403  /** \brief Test the thermodynamics of calc_density() and
404  calc_mu()
405 
406  This compares the approximation to the exact results over a
407  grid with \f$ T = \left\{10^{-2},1,10^{2}\right\} \f$, \f$
408  \log_{10} (m/T) = \left\{ -3,-2,-1,0,1,2,3\right\} \f$, and
409  \f$ \log_{10} \psi = \left\{ -3,-2,-1,0,1\right\} \f$, where
410  \f$ \psi \equiv \left(\mu-m\right)/T \f$ using
411  calc_density() and calc_mu(), with both inc_rest_mass
412  taking both values <tt>true</tt> and <tt>false</tt>.
413 
414  The <tt>verbose</tt> parameter controls the amount of
415  output, and \c fname is the filename for the file
416  <tt>fermion_cal.o2</tt>.
417 
418  \future Also calibrate massless fermions?
419  \future Convert into separate class?
420  */
421  virtual double calibrate(fermion &f, int verbose=0, std::string fname="");
422 
423 #ifndef DOXYGEN_NO_O2NS
424 
425  protected:
426 
427  /// A pointer to the solver for massless fermions
429 
430  /// Solve for the chemical potential for massless fermions
431  double massless_solve_fun(double x, fermion &f, double temper);
432 
433 #endif
434 
435  };
436 
437 #ifndef DOXYGEN_NO_O2NS
438 }
439 #endif
440 
441 #endif
fermion & operator=(const fermion &f)
Copy construction with operator=()
Definition: fermion.h:88
double nu
Effective chemical potential.
Definition: part.h:113
Fermion properties at zero temperature.
Definition: fermion.h:123
double n
Number density.
Definition: part.h:101
root_cern def_massless_root
The default solver for massless_calc_density()
Definition: fermion.h:398
Fermion class.
Definition: fermion.h:52
fermion(double mass=0, double dof=0)
Create a fermion with mass mass and degeneracy dof.
double en
Entropy density.
Definition: part.h:109
Fermion with finite-temperature thermodynamics [abstract base].
Definition: fermion.h:199
fermion(const fermion &f)
Copy constructor.
Definition: fermion.h:71
double m
Mass.
Definition: part.h:99
virtual const char * type()
Return string denoting type ("fermion_eval_thermo")
Definition: fermion.h:401
root * massless_root
A pointer to the solver for massless fermions.
Definition: fermion.h:428
double del
Gap.
Definition: fermion.h:59
Particle base class.
Definition: part.h:92
double mu
Chemical potential.
Definition: part.h:107
double pr
Pressure.
Definition: part.h:105
double ms
Effective mass (Dirac unless otherwise specified)
Definition: part.h:111
double kf
Fermi momentum.
Definition: fermion.h:57
double ed
Energy density.
Definition: part.h:103
bool non_interacting
True if the particle is non-interacting (default true)
Definition: part.h:119
bool inc_rest_mass
If true, include the mass in the energy density and chemical potential (default true) ...
Definition: part.h:117
double g
Degeneracy (e.g. spin and color if applicable)
Definition: part.h:97
void set_massless_root(root<> &rp)
Set the solver for use in massless_calc_density()
Definition: fermion.h:388
virtual const char * type()
Return string denoting type ("fermion")
Definition: fermion.h:68

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).