fermion.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2020, 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_tl
28 */
29 #include <string>
30 #include <iostream>
31 #include <fstream>
32 #include <cmath>
33 
34 #include <boost/math/constants/constants.hpp>
35 #include <boost/math/special_functions/bessel.hpp>
36 
37 // For gsl_sf_fermi_dirac_int()
38 #include <gsl/gsl_specfunc.h>
39 
40 #include <o2scl/constants.h>
41 #include <o2scl/funct.h>
42 #include <o2scl/root.h>
43 #include <o2scl/root_cern.h>
44 #include <o2scl/part.h>
45 
46 #ifndef DOXYGEN_NO_O2NS
47 namespace o2scl {
48 #endif
49 
50  /** \brief Fermion class
51 
52  This class adds two member data variables, \ref kf and \ref
53  del, for the Fermi momentum and the gap, respectively.
54  */
55  template<class fp_t=double> class fermion_tl : public part_tl<fp_t> {
56 
57  public:
58 
59  /// Fermi momentum
60  fp_t kf;
61  /// Gap
62  fp_t del;
63 
64  /// Create a fermion with mass \c mass and degeneracy \c dof.
65  fermion_tl(fp_t mass=0, fp_t dof=0) : part_tl<fp_t>(mass,dof) {
66  kf=0.0;
67  del=0.0;
68  }
69 
70  virtual ~fermion_tl() {
71  }
72 
73  /// Return string denoting type ("fermion_tl")
74  virtual const char *type() { return "fermion_tl"; }
75 
76  /// Copy constructor
77  fermion_tl(const fermion_tl &f) {
78  this->g=f.g;
79  this->m=f.m;
80  this->ms=f.ms;
81  this->n=f.n;
82  this->ed=f.ed;
83  this->pr=f.pr;
84  this->mu=f.mu;
85  this->en=f.en;
86  this->nu=f.nu;
89  kf=f.kf;
90  del=f.del;
91  }
92 
93  /// Copy construction with operator=()
95  if (this!=&f) {
96  this->g=f.g;
97  this->m=f.m;
98  this->ms=f.ms;
99  this->n=f.n;
100  this->ed=f.ed;
101  this->pr=f.pr;
102  this->mu=f.mu;
103  this->en=f.en;
104  this->nu=f.nu;
107  kf=f.kf;
108  del=f.del;
109  }
110  return *this;
111  }
112 
113  };
114 
115  typedef fermion_tl<double> fermion;
116 
117  /** \brief Fermion properties at zero temperature
118 
119  This is a base class for the computation of fermionic statistics
120  at zero temperature. The more general case of finite temperature
121  is taken care of by \ref fermion_thermo_tl objects. The
122  primary functions are calc_mu_zerot() and calc_density_zerot()
123  which compute all the thermodynamic quantities as a function of
124  the chemical potential, or the density, respectively.
125 
126  \future Use hypot() and other more accurate functions for the
127  analytic expressions for the zero temperature integrals. [Progress
128  has been made, but there are probably other functions which may
129  break down for small but finite masses and temperatures]
130  */
131  template<class fp_t=double> class fermion_zerot_tl {
132 
133  protected:
134 
135  /// Desc
136  fp_t pi;
137 
138  /// Desc
139  fp_t pi2;
140 
141  public:
142 
143  fermion_zerot_tl() {
144  pi=boost::math::constants::pi<fp_t>();
145  pi2=boost::math::constants::pi_sqr<fp_t>();
146  };
147 
148  virtual ~fermion_zerot_tl() {
149  }
150 
151  /// \name Zero-temperature fermions
152  //@{
153  /** \brief Calculate the Fermi momentum from the density
154 
155  Uses the relation \f$ k_F = ( 6 \pi^2 n /g )^{1/3} \f$
156  */
158  f.kf=cbrt(6.0*pi2/f.g*f.n);
159  return;
160  }
161 
162  /** \brief Energy density at T=0 from \ref o2scl::fermion_tl::kf and
163  \ref o2scl::part_tl::ms
164 
165  Calculates the integral
166  \f[
167  \varepsilon = \frac{g}{2 \pi^2} \int_0^{k_F} k^2
168  \sqrt{k^2+m^{* 2}} d k
169  \f]
170  */
172  fp_t r,efs;
173  if (f.kf>0.0) {
174  if (f.ms<=0.0) {
175  f.ed=f.g*(pow(f.kf,4.0)/8.0/pi2);
176  } else {
177  efs=gsl_hypot(f.kf,f.ms);
178  r=(f.kf+efs)/f.ms;
179  f.ed=f.g/16.0/pi2*(2.0*f.kf*pow(efs,3.0)-f.kf*efs*f.ms*f.ms
180  -pow(f.ms,4.0)*log(r));
181  }
182  } else {
183  f.ed=0.0;
184  }
185  return;
186  }
187 
188  /** \brief Pressure at T=0 from \ref o2scl::fermion_tl::kf and
189  \ref o2scl::part_tl::ms
190 
191  Calculates the integral
192  \f[
193  P=\frac{g}{6 \pi^2} \int_0^{k_F} \frac{k^4}{\sqrt{k^2+m^{* 2}}} d k
194  \f]
195  */
197  fp_t r,efs;
198  if (f.kf>0.0) {
199  if (f.ms<=0.0) {
200  f.pr=f.g*(pow(f.kf,4.0)/24.0/pi2);
201  } else {
202  efs=gsl_hypot(f.kf,f.ms);
203  r=(f.kf+efs)/f.ms;
204  f.pr=f.g/48.0/pi2*(2.0*efs*pow(f.kf,3.0)-
205  3.0*f.kf*efs*f.ms*f.ms
206  +3.0*pow(f.ms,4.0)*log(r));
207  }
208  } else {
209  f.pr=0.0;
210  }
211  return;
212  }
213 
214  /** \brief Zero temperature fermions from \ref o2scl::part_tl::mu or
215  \ref o2scl::part_tl::nu and \ref o2scl::part_tl::ms
216  */
217  virtual void calc_mu_zerot(fermion_tl<fp_t> &f) {
218  bool nulessthan0;
219  if (f.non_interacting) { f.nu=f.mu; f.ms=f.m; }
220 
221  if (f.inc_rest_mass) {
222  if (f.nu>f.ms) {
223  nulessthan0=false;
224  f.kf=sqrt(f.nu*f.nu-f.ms*f.ms);
225  } else {
226  nulessthan0=false;
227  f.kf=0.0;
228  }
229  } else {
230  fp_t nupm=f.nu+f.m;
231  if ((nupm)>f.ms) {
232  nulessthan0=false;
233  f.kf=sqrt(nupm*nupm-f.ms*f.ms);
234  } else {
235  nulessthan0=false;
236  f.kf=0.0;
237  }
238  }
239  f.n=f.g/6.0/pi2*pow(f.kf,3.0);
241  pressure_zerot(f);
242  f.en=0.0;
243  if (!f.inc_rest_mass) f.ed-=f.n*f.m;
244  if (nulessthan0==true) {
245  f.n*=-1.0;
246  f.kf*=-1.0;
247  }
248 
249  return;
250  }
251 
252  /** \brief Zero temperature fermions from \ref o2scl::part_tl::n
253  and \ref o2scl::part_tl::ms
254  */
256  if (f.non_interacting) { f.ms=f.m; }
257 
258  f.kf=cbrt(6.0*pi2/f.g*f.n);
259  f.nu=gsl_hypot(f.kf,f.ms);
261  pressure_zerot(f);
262  f.en=0.0;
263 
264  if (!f.inc_rest_mass) {
265  f.nu-=f.m;
266  f.ed-=f.n*f.m;
267  }
268 
269  if (f.non_interacting) { f.mu=f.nu; }
270 
271  return;
272  }
273  //@}
274 
275  };
276 
277  /** \brief Double-precision version of \ref o2scl::fermion_zerot_tl
278  */
280 
281  /** \brief Fermion with finite-temperature thermodynamics
282  [abstract base]
283 
284  This is an abstract base for the computation of
285  finite-temperature fermionic statistics. Different children
286  (e.g. \ref fermion_eff and \ref fermion_rel_tl) use different
287  techniques to computing the momentum integrations.
288 
289  Because massless fermions at finite temperature are much
290  simpler, there are separate member functions included in this
291  class to handle them. The functions massless_calc_density() and
292  massless_calc_mu() compute the thermodynamics of massless
293  fermions at finite temperature given the density or the chemical
294  potentials. The functions massless_pair_density() and
295  massless_pair_mu() perform the same task, but automatically
296  include antiparticles.
297 
298  The function massless_calc_density() uses a \ref root object to
299  solve for the chemical potential as a function of the density.
300  The default is an object of type root_cern. The function
301  massless_pair_density() does not need to use the \ref root
302  object because of the simplification afforded by the inclusion
303  of antiparticles.
304 
305  \future Create a Chebyshev approximation for inverting the
306  the Fermi functions for massless_calc_density() functions?
307  */
308  template<class fp_t=double>
309  class fermion_thermo_tl : public fermion_zerot_tl<fp_t> {
310 
311  public:
312 
315  }
316 
317  virtual ~fermion_thermo_tl() {
318  }
319 
320  /** \brief Calculate thermodynamic properties from the chemical
321  potential using a nondegenerate expansion
322  */
323  template<class fermion_t>
324  bool calc_mu_ndeg_tlate(fermion_t &f, fp_t temper,
325  fp_t prec, bool inc_antip) {
326 
327  if (f.non_interacting==true) { f.nu=f.mu; f.ms=f.m; }
328 
329  // Compute psi and tt
330  fp_t psi, psi_num;
331  if (f.inc_rest_mass) {
332  psi_num=f.nu-f.ms;
333  } else {
334  psi_num=f.nu+f.m-f.ms;
335  }
336  psi=psi_num/temper;
337  fp_t tt=temper/f.ms;
338 
339  // Return false immediately if we're degenerate
340  if (inc_antip==false && psi>-1.0) return false;
341 
342  // Prefactor 'd' in Johns96
343  fp_t prefac=f.g/2.0/this->pi2*pow(f.ms,4.0);
344 
345  // One term is always used, so only values of max_term greater than
346  // 0 are useful.
347  static const size_t max_term=200;
348 
349  // Maximum argument for exponential
350  // fp_t log_dbl_max=709.78;
351 
352  // Return zero if psi+1/t is too small
353  if (psi+1.0/tt<-700.0) {
354  f.n=0.0;
355  f.ed=0.0;
356  f.pr=0.0;
357  f.en=0.0;
358  return true;
359  }
360 
361  // -----------------------------------------------------
362  // Return early if the last term is going to be too large.
363 
364  // Ratio of last term to first term in the pressure expansion
365  fp_t rat;
366  fp_t dj1=((fp_t)max_term), jot1=max_term/tt;
367  fp_t dj2=1.0, jot2=1.0/tt;
368 
369  // fp_t Kn21=boost::math::cyl_bessel_k(2,jot1)*exp(jot1);
370  // fp_t Kn22=boost::math::cyl_bessel_k(2,jot2)*exp(jot2);
371  //std::cout << boost::math::cyl_bessel_k(2,jot1)*exp(jot1) << std::endl;
372  //std::cout << gsl_sf_bessel_Kn_scaled(2.0,jot1) << std::endl;
373  //exit(-1);
374 
375  if (inc_antip==false) {
376  rat=exp(dj1*psi)/jot1/jot1*gsl_sf_bessel_Kn_scaled(2.0,jot1);
377  rat/=exp(dj2*psi)/jot2/jot2*gsl_sf_bessel_Kn_scaled(2.0,jot2);
378  } else {
379  if (f.inc_rest_mass) {
380  rat=exp(-jot1)*2.0*cosh(dj1*f.nu/temper)/jot1/jot1*
381  gsl_sf_bessel_Kn_scaled(2.0,jot1);
382  rat/=exp(-jot2)*2.0*cosh(dj2*f.nu/temper)/jot2/jot2*
383  gsl_sf_bessel_Kn_scaled(2.0,jot2);
384  } else {
385  rat=exp(-jot1)*2.0*cosh(dj1*(f.nu+f.m)/temper)/jot1/jot1*
386  gsl_sf_bessel_Kn_scaled(2.0,jot1);
387  rat/=exp(-jot2)*2.0*cosh(dj2*(f.nu+f.m)/temper)/jot2/jot2*
388  gsl_sf_bessel_Kn_scaled(2.0,jot2);
389  }
390  }
391 
392  // If the ratio between the last term and the first term is
393  // not small enough, return false
394  if (std::isfinite(rat) && rat>prec) {
395  return false;
396  }
397 
398  fp_t first_term=0.0;
399  f.pr=0.0;
400  f.n=0.0;
401  f.en=0.0;
402 
403  for(size_t j=1;j<=max_term;j++) {
404 
405  fp_t pterm, nterm, enterm;
406 
407  ndeg_terms(j,tt,psi*tt,f.ms,f.inc_rest_mass,inc_antip,
408  pterm,nterm,enterm);
409 
410  if (j==1) first_term=pterm;
411  f.pr+=pterm;
412  f.n+=nterm;
413  f.en+=enterm;
414 
415  // If the first term is zero, then the rest of the terms
416  // will be zero so just return early
417  if (first_term==0.0) {
418  f.pr=0.0;
419  f.n=0.0;
420  f.ed=0.0;
421  f.en=0.0;
422  return true;
423  }
424 
425  // Stop if the last term is sufficiently small compared to
426  // the first term
427  if (j>1 && fabs(pterm)<prec*fabs(first_term)) {
428  f.pr*=prefac;
429  f.n*=prefac;
430  f.en*=prefac;
431  f.ed=-f.pr+f.nu*f.n+temper*f.en;
432  return true;
433  }
434 
435  // End of 'for(size_t j=1;j<=max_term;j++)'
436  }
437 
438  // We failed to add enough terms, so return false
439  return false;
440  }
441 
442  /** \brief Calculate thermodynamic properties from the chemical
443  potential using a degenerate expansion
444  */
445  template<class fermion_t>
446  bool calc_mu_deg_tlate(fermion_t &f, fp_t temper,
447  fp_t prec) {
448 
449  // Handle the zero-temperature limit
450  if (temper==0.0) {
451  this->calc_mu_zerot(f);
452  return true;
453  }
454 
455  // Double check to ensure T and mass are positive
456  if (temper<0.0 || f.ms<0.0) {
457  O2SCL_ERR2("Temperature or mass negative in fermion_thermo",
458  "::calc_mu_deg_tlate().",exc_einval);
459  }
460 
461  if (f.non_interacting==true) { f.nu=f.mu; f.ms=f.m; }
462 
463  // Compute psi and tt
464  fp_t psi;
465  if (f.inc_rest_mass) psi=(f.nu-f.ms)/temper;
466  else psi=(f.nu+f.m-f.ms)/temper;
467  fp_t tt=temper/f.ms;
468 
469  // Return false immediately psi<0 where the expressions below
470  // don't work because of the square roots
471  if (psi<0.0) return false;
472 
473  // Prefactor 'd' in Johns96
474  fp_t prefac=f.g/2.0/this->pi2*pow(f.ms,4.0);
475 
476  // Define x = psi * t = (mu/m - 1) and related values
477  fp_t x=psi*tt;
478  fp_t sx=sqrt(x);
479  fp_t s2x=sqrt(2.0+x);
480  fp_t x2=x*x;
481  fp_t x3=x2*x;
482  fp_t x4=x2*x2;
483 
484  // Evaluate the first and last term for the pressure
485  fp_t pterm1;
486  if (x>1.0e-5) {
487  pterm1=(x*(1.0+x)*(2.0+x)*(-3.0+2.0*x*(2.0+x))+6.0*sx*s2x*
488  log((sx+s2x)/sqrt(2.0)))/24.0/sx/s2x;
489  } else {
490  pterm1=x2*sx*(29568.0+15840.0*x+1540.0*x2-105.0*x3)/55440.0/sqrt(2.0);
491  }
492  fp_t pterm4=-31.0*pow(this->pi*tt,6.0)/1008.0*(1.0+x)*
493  sx*s2x/pow(x*(2.0+x),4.0);
494 
495  // Check if we're going to succeed
496  if (fabs(pterm4)/fabs(pterm1)>prec) {
497  return false;
498  }
499 
500  // First order density term (first order entropy term is zero)
501  fp_t nterm1=sx*s2x*x*(2.0+x)/3.0/f.ms;
502 
503  // Second order terms
504  fp_t pterm2=tt*tt*this->pi2/6.0*(1.0+x)*sx*s2x;
505  fp_t nterm2=tt*tt*this->pi2/6.0*(1.0+4.0*x+2.0*x2)/
506  f.ms/sx/s2x;
507  fp_t enterm2=tt*this->pi2/3.0*(1.0+x)*sx*s2x/f.ms;
508 
509  // Third order terms
510  fp_t pterm3=7.0*pow(this->pi*tt,4.0)/360.0*(1.0+x)*
511  (-1.0+4.0*x+2.0*x2)/pow(x*(2.0+x),1.5);
512  fp_t nterm3=7.0*pow(this->pi*tt,4.0)/120.0/sx/s2x/
513  x2/(x+2.0)/(x+2.0)/f.ms;
514  fp_t enterm3=7.0*pow(this->pi*tt,4.0)/tt/90.0*(1.0+x)*
515  (-1.0+4.0*x+2.0*x2)/f.ms/sx/s2x/x/(x+2.0);
516 
517  // Fourth order terms for density and entropy
518  fp_t nterm4=31.0*pow(this->pi*tt,6.0)/1008.0*sx*s2x*
519  (7.0+12.0*x+6.0*x2)/f.ms/pow(x*(2.0+x),5.0);
520  fp_t enterm4=-31.0*pow(this->pi*tt,6.0)/tt/168.0*sx*s2x*
521  (1.0+x)/pow(x*(2.0+x),4.0);
522 
523  // Add up all the terms
524  f.pr=prefac*(pterm1+pterm2+pterm3+pterm4);
525  f.n=prefac*(nterm1+nterm2+nterm3+nterm4);
526  f.en=prefac*(enterm2+enterm3+enterm4);
527  f.ed=-f.pr+f.nu*f.n+temper*f.en;
528 
529  return true;
530  }
531 
532  /** \brief Non-degenerate expansion for fermions
533 
534  Attempts to evaluate thermodynamics of a non-degenerate
535  fermion. If the result is accurate to within the requested
536  precision, this function returns <tt>true</tt>, and otherwise
537  this function returns <tt>false</tt> and the values in stored
538  in the <tt>pr</tt>, <tt>n</tt>, <tt>en</tt>, and <tt>ed</tt>
539  field are meaningless.
540 
541  If \f$ \mu \f$ is negative and sufficiently far from zero,
542  then the thermodynamic quantities are smaller than the smallest
543  representable double-precision number. In this case,
544  this function will return <tt>true</tt> and report all
545  quantities as zero.
546 
547  Defining \f$ \psi \equiv (\mu-m)/T \f$, \f$ t \equiv T/m \f$,
548  and \f$ d \equiv g~m^4/(2 \pi^2) \f$ the pressure
549  in the non-degenerate limit (\f$ \psi \rightarrow - \infty \f$)
550  is (\ref Johns96)
551  \f[
552  P = d \sum_{n=1}^{\infty} P_n
553  \f]
554  where
555  \f[
556  P_n \equiv \left(-1\right)^{n+1} \left(\frac{t^2}{n^2}\right)
557  e^{n \left(\psi+1/t\right)} K_2 \left( \frac{n}{t} \right)
558  \f]
559  The density is then
560  \f[
561  n = d \sum_{n=1}^{\infty} \frac{n P_n}{T}
562  \f]
563  and the entropy density is
564  \f[
565  s = \frac{d}{m} \sum_{n=1}^{\infty} \left\{ \frac{2 P_n}{t}
566  -\frac{n P_n}{t^2}+
567  \frac{\left(-1\right)^{n+1}}{2 n}
568  e^{n \left(\psi+1/t\right)} \left[ K_1 \left( \frac{n}{t}
569  \right)+K_3 \left( \frac{n}{t} \right) \right]
570  \right\}
571  \f]
572 
573  This function is accurate over a wide range of conditions
574  when \f$ \psi < -4 \f$.
575 
576  The ratio of the nth term to the first term in the pressure
577  series is
578  \f[
579  R_n \equiv \frac{P_{n}}{P_{1}} = \frac{(-1)^{n+1}
580  e^{(n-1)(\psi+1/t)} K_2(n/t) }{n^2 K_2(1/t)}
581  \f]
582  This function currently uses 20 terms in the series and
583  immediately returns <tt>false</tt> if \f$ |R_{20}| \f$
584  is greater than <tt>prec</tt>
585 
586  In the nondegenerate and nonrelativistic (\f$ t \rightarrow 0
587  \f$) limit, the argument to the Bessel functions and the
588  exponential becomes too large. In this case, it's better to
589  use the expansions, e.g. for \f$ x \equiv n/t \rightarrow
590  \infty \f$,
591  \f[
592  \sqrt{\frac{2 x}{\pi}} e^{x} K_2(x) \approx
593  1 + \frac{3}{8 x} - \frac{15}{128 x^2} + ...
594  \f]
595  The current code currently goes up to \f$ x^{-12} \f$ in the
596  expansion, which is enough for the default precision of \f$
597  10^{-18} \f$ since \f$ (20/700)^{12} \sim 10^{-19} \f$.
598  */
599  virtual bool calc_mu_ndeg(fermion &f, fp_t temper,
600  fp_t prec=1.0e-18, bool inc_antip=false) {
601  return calc_mu_ndeg_tlate<fermion>(f,temper,prec,inc_antip);
602  }
603 
604  /** \brief Degenerate expansion for fermions
605 
606  Attempts to evaulate thermodynamics of a degenerate fermion.
607  If the result is accurate to within the requested precision,
608  this function returns <tt>true</tt>, and otherwise this
609  function returns <tt>false</tt> and the values in stored in
610  the <tt>pr</tt>, <tt>n</tt>, <tt>en</tt>, and <tt>ed</tt>
611  field are meaningless.
612 
613  The pressure, density, and energy density, should be accurate
614  to the requested precision, but the first term in the series
615  expansion for the entropy is zero, so the entropy is one order
616  lower in accuracy.
617 
618  \future Make a function like this for dndm, dsdT, etc.
619  for fermion_deriv .
620  */
621  virtual bool calc_mu_deg(fermion &f, fp_t temper,
622  fp_t prec=1.0e-18) {
623  return calc_mu_deg_tlate<fermion>(f,temper,prec);
624  }
625 
626  /** \brief Calculate properties as function of chemical potential
627  */
628  virtual void calc_mu(fermion &f, fp_t temper)=0;
629 
630  /** \brief Calculate properties as function of density
631 
632  \note This function returns an integer value, in contrast to
633  \ref calc_mu(), because of the potential for non-convergence.
634  */
635  virtual int calc_density(fermion &f, fp_t temper)=0;
636 
637  /** \brief Calculate properties with antiparticles as function of
638  chemical potential
639  */
640  virtual void pair_mu(fermion &f, fp_t temper)=0;
641 
642  /** \brief Calculate properties with antiparticles as function of
643  density
644 
645  \note This function returns an integer value, in contrast to
646  \ref pair_mu(), because of the potential for non-convergence.
647  */
648  virtual int pair_density(fermion &f, fp_t temper)=0;
649 
650  /// \name Massless fermions
651  //@{
652  /// Finite temperature massless fermions
653  virtual void massless_calc_mu(fermion &f, fp_t temper) {
654 
655  fp_t fm2, fm3;
656 
657  if (f.non_interacting) { f.nu=f.mu; }
658 
659  fm2=gsl_sf_fermi_dirac_int(2,f.nu/temper);
660  fm3=gsl_sf_fermi_dirac_int(3,f.nu/temper);
661 
662  f.n=f.g/this->pi2*pow(temper,3.0)*fm2;
663  f.ed=f.g*3.0/this->pi2*pow(temper,4.0)*fm3;
664  f.pr=f.ed/3.0;
665  f.en=(f.ed+f.pr-f.n*f.nu)/temper;
666 
667  return;
668  }
669 
670  /// Finite temperature massless fermions
671  virtual void massless_calc_density(fermion &f, fp_t temper) {
672  fp_t x, T=temper;
673 
674  x=f.ms+temper;
675  funct mf2=std::bind(std::mem_fn<fp_t(fp_t,fermion &,fp_t)>
677  this,std::placeholders::_1,std::ref(f),temper);
678  massless_root->solve(x,mf2);
679  f.nu=x;
680 
681  massless_calc_mu(f,temper);
682 
683  // If the particle is non-interacting, then need to set
684  // mu=nu to get the entropy right
685  if (f.non_interacting) { f.mu=f.nu; }
686 
687  return;
688  }
689 
690  /** \brief Finite temperature massless fermions and antifermions
691  */
692  virtual void massless_pair_mu(fermion &f, fp_t temper) {
693  fp_t pitmu, pitmu2, nu2;
694 
695  if (f.non_interacting) { f.nu=f.mu; f.ms=f.m; }
696  if (f.nu==0.0) {
697  f.n=0.0;
698  f.ed=f.g/8.0/this->pi2*7.0/15.0*
699  pow(this->pi*temper,4.0);
700  f.pr=f.ed/3.0;
701  f.en=(f.ed+f.pr-f.n*f.mu)/temper;
702  } else {
703  nu2=f.nu*f.nu;
704  pitmu=this->pi*temper/f.nu;
705  pitmu2=pitmu*pitmu;
706  f.n=f.g*f.nu*nu2/6.0/this->pi2*(1.0+pitmu2);
707  f.ed=f.g*nu2*nu2/8.0/this->pi2*(1.0+2.0*pitmu2+
708  7.0/15.0*pitmu2*pitmu2);
709  f.pr=f.ed/3.0;
710  f.en=(f.ed+f.pr-f.n*f.mu)/temper;
711 
712  // Might the following work better for the energy density?
713  // pit=pi*temper;
714  // pit2=pit*pit;
715  // ed=g/8.0/pi2*(nu2*nu2+2.0*pit2*nu2+7.0/15.0*pit2*pit2);
716 
717  }
718 
719  return;
720  }
721 
722  /** \brief Finite temperature massless fermions and antifermions
723 
724  In the cases \f$ n^3 \gg T \f$ and \f$ T \gg n^3 \f$ ,
725  expansions are used instead of the exact formulas to avoid
726  loss of precision.
727 
728  In particular, using the parameter
729  \f[
730  \alpha = \frac{g^2 \pi^2 T^6}{243 n^2}
731  \f]
732  and defining the expression
733  \f[
734  \mathrm{cbt} = \alpha^{-1/6} \left( -1 + \sqrt{1+\alpha}\right)^{1/3}
735  \f]
736  we can write the chemical potential as
737  \f[
738  \mu = \frac{\pi T}{\sqrt{3}} \left(\frac{1}{\mathrm{cbt}} -
739  \mathrm{cbt} \right)
740  \f]
741 
742  These expressions, however, do not work well when \f$ \alpha
743  \f$ is very large or very small, so series expansions are
744  used whenever \f$ \alpha > 10^{4} \f$ or
745  \f$ \alpha < 3 \times 10^{-4} \f$. For small \f$ \alpha \f$,
746  \f[
747  \left(\frac{1}{\mathrm{cbt}} -
748  \mathrm{cbt} \right) \approx
749  \frac{2^{1/3}}{\alpha^{1/6}} -
750  \frac{\alpha^{1/6}}{2^{1/3}} +
751  \frac{\alpha^{5/6}}{6{\cdot}2^{2/3}} +
752  \frac{\alpha^{7/6}}{12{\cdot}2^{1/3}} -
753  \frac{\alpha^{11/6}}{18{\cdot}2^{2/3}} -
754  \frac{5 \alpha^{13/6}}{144{\cdot}2^{1/3}} +
755  \frac{77 \alpha^{17/6}}{2592{\cdot}2^{2/3}}
756  \f]
757  and for large \f$ \alpha \f$,
758  \f[
759  \left(\frac{1}{\mathrm{cbt}} -
760  \mathrm{cbt} \right) \approx
761  \frac{2}{3} \sqrt{\frac{1}{\alpha}} -
762  \frac{8}{81} \left(\frac{1}{\alpha}\right)^{3/2} +
763  \frac{32}{729} \left(\frac{1}{\alpha}\right)^{5/2}
764  \f]
765 
766  This approach works to within about 1 \part in \f$ 10^{12} \f$,
767  and is tested in <tt>fermion_ts.cpp</tt>.
768 
769  \future This could be improved by including more terms
770  in the expansions.
771  */
772  virtual void massless_pair_density(fermion &f, fp_t temper) {
773 
774  fp_t t2=temper*temper,pitmu,pitmu2,nu2;
775  fp_t cbt, alpha, two13, alpha16;
776 
777  if (f.non_interacting) { f.ms=f.m; }
778  if (f.n<=0.0) {
779  f.nu=0.0;
780  f.ed=f.g/8.0/this->pi2*7.0/15.0*pow(this->pi*temper,4.0);
781  f.pr=f.ed/3.0;
782  } else {
783  alpha=f.g*f.g*this->pi2*t2*t2*t2/243.0/f.n/f.n;
784  if (alpha>1.0e4) {
785  f.nu=(2.0/3.0/sqrt(alpha)-8.0/81.0/pow(alpha,1.5)+
786  32.0/729.0/pow(alpha,2.5))*this->pi*temper/sqrt(3.0);
787  } else if (alpha<3.0e-4) {
788  two13=cbrt(2.0);
789  alpha16=pow(alpha,1.0/6.0);
790  f.nu=(two13/alpha16-alpha16/two13+alpha/alpha16/6.0/two13/two13
791  +alpha*alpha16/12.0/two13-alpha*alpha/alpha16/18.0/two13/two13-
792  5.0*alpha*alpha*alpha16/144.0/two13+
793  77.0/2592.0*alpha*alpha*alpha/alpha16/two13/two13)*
794  this->pi*temper/sqrt(3.0);
795  } else {
796  cbt=pow(-1.0+sqrt(1.0+alpha),1.0/3.0)/pow(alpha,1.0/6.0);
797  f.nu=this->pi*temper/sqrt(3.0)*(1.0/cbt-cbt);
798  }
799  pitmu=this->pi*temper/f.nu;
800  pitmu2=pitmu*pitmu;
801  nu2=f.nu*f.nu;
802  f.ed=f.g*nu2*nu2/8.0/this->pi2*
803  (1.0+2.0*pitmu2+7.0/15.0*pitmu2*pitmu2);
804  f.pr=f.ed/3.0;
805 
806  if (!std::isfinite(f.nu)) {
807  std::string str="Chemical potential not finite ("+dtos(f.nu)+
808  ") in fermion::massless_pair_density().";
809  O2SCL_ERR(str.c_str(),exc_efailed);
810  }
811  }
812 
813  if (f.non_interacting) { f.mu=f.nu; }
814  f.en=(f.ed+f.pr-f.n*f.nu)/temper;
815 
816  return;
817  }
818  //@}
819 
820  /** \brief Set the solver for use in massless_calc_density() */
822  massless_root=&rp;
823  return;
824  }
825 
826  /** \brief The default solver for massless_calc_density()
827 
828  We default to a solver of type root_cern here since we
829  don't have a bracket or a derivative.
830  */
832 
833  /// Return string denoting type ("fermion_thermo")
834  virtual const char *type() { return "fermion_thermo"; }
835 
836  /** \brief Compute a term in the nondegenerate expansion
837  */
838  void ndeg_terms(size_t j, fp_t tt,
839  fp_t xx, fp_t m, bool inc_rest_mass,
840  bool inc_antip, fp_t &pterm, fp_t &nterm,
841  fp_t &enterm) {
842 
843  fp_t dj=((fp_t)j);
844  fp_t jot=dj/tt;
845 
846  // fp_t Kn1=boost::math::cyl_bessel_k(1,jot1)*exp(jot);
847  // fp_t Kn2=boost::math::cyl_bessel_k(2,jot1)*exp(jot);
848  // fp_t Kn3=boost::math::cyl_bessel_k(2,jot1)*exp(jot);
849 
850  if (inc_antip==false) {
851  pterm=exp(jot*xx)/jot/jot*gsl_sf_bessel_Kn_scaled(2.0,jot);
852  if (j%2==0) pterm*=-1.0;
853  nterm=pterm*jot/m;
854  fp_t enterm1=(4.0*tt-dj*xx-dj)/dj/tt*nterm;
855  fp_t enterm2=exp(jot*xx)/dj*gsl_sf_bessel_Kn_scaled(1.0,jot)/m;
856  if (j%2==0) {
857  enterm=enterm1-enterm2;
858  } else {
859  enterm=enterm1+enterm2;
860  }
861  } else {
862  pterm=exp(-jot)*2.0*cosh(jot*(xx+1.0)/tt)/jot/jot*
863  gsl_sf_bessel_Kn_scaled(2.0,jot);
864  if (j%2==0) pterm*=-1.0;
865  nterm=pterm*tanh(jot*(xx+1.0))*jot;
866  fp_t enterm1=-(1.0+xx)/tt*nterm/m;
867  fp_t enterm2=2.0*exp(-jot*xx)/dj*cosh(jot*(xx+1.0))*
868  gsl_sf_bessel_Kn_scaled(3.0,jot)/m;
869  if (j%2==0) {
870  enterm=enterm1-enterm2;
871  } else {
872  enterm=enterm1+enterm2;
873  }
874  }
875 
876  return;
877  }
878 
879 #ifndef DOXYGEN_NO_O2NS
880 
881  protected:
882 
883  /// A pointer to the solver for massless fermions
885 
886  /// Solve for the chemical potential for massless fermions
887  fp_t massless_solve_fun(fp_t x, fermion &f, fp_t temper) {
888  fp_t fm2=gsl_sf_fermi_dirac_int(2,x/(temper));
889  return f.g*pow(temper,3.0)*fm2/this->pi2/f.n-1.0;
890  }
891 
892 #endif
893 
894  };
895 
896  /** \brief Double-precision version of \ref o2scl::fermion_thermo_tl
897  */
899 
900 #ifndef DOXYGEN_NO_O2NS
901 }
902 #endif
903 
904 #endif
o2scl::part_tl
Particle base class.
Definition: part.h:103
o2scl::fermion_thermo_tl::calc_mu_deg_tlate
bool calc_mu_deg_tlate(fermion_t &f, fp_t temper, fp_t prec)
Calculate thermodynamic properties from the chemical potential using a degenerate expansion.
Definition: fermion.h:446
o2scl::fermion_tl::operator=
fermion_tl & operator=(const fermion_tl &f)
Copy construction with operator=()
Definition: fermion.h:94
o2scl::fermion_zerot_tl::energy_density_zerot
void energy_density_zerot(fermion_tl< fp_t > &f)
Energy density at T=0 from o2scl::fermion_tl::kf and o2scl::part_tl::ms.
Definition: fermion.h:171
o2scl::fermion_thermo_tl::massless_calc_density
virtual void massless_calc_density(fermion &f, fp_t temper)
Finite temperature massless fermions.
Definition: fermion.h:671
o2scl::fermion_thermo_tl::calc_mu
virtual void calc_mu(fermion &f, fp_t temper)=0
Calculate properties as function of chemical potential.
o2scl::root::solve
virtual int solve(fp_t &x, func_t &func)=0
dtos
std::string dtos(const fp_t &x, int prec=6, bool auto_prec=false)
o2scl::fermion_zerot_tl::pi
fp_t pi
Desc.
Definition: fermion.h:136
o2scl::part_tl< double >::mu
double mu
Chemical potential.
Definition: part.h:118
o2scl::part_tl< double >::ms
double ms
Effective mass (Dirac unless otherwise specified)
Definition: part.h:122
exc_efailed
exc_efailed
O2SCL_ERR2
#define O2SCL_ERR2(d, d2, n)
o2scl::part_tl< double >::n
double n
Number density.
Definition: part.h:112
o2scl::fermion_zerot_tl::calc_density_zerot
virtual void calc_density_zerot(fermion_tl< fp_t > &f)
Zero temperature fermions from o2scl::part_tl::n and o2scl::part_tl::ms.
Definition: fermion.h:255
o2scl::fermion_thermo_tl::pair_mu
virtual void pair_mu(fermion &f, fp_t temper)=0
Calculate properties with antiparticles as function of chemical potential.
o2scl::part_tl< double >::non_interacting
bool non_interacting
True if the particle is non-interacting (default true)
Definition: part.h:130
o2scl::fermion_thermo_tl::calc_mu_ndeg_tlate
bool calc_mu_ndeg_tlate(fermion_t &f, fp_t temper, fp_t prec, bool inc_antip)
Calculate thermodynamic properties from the chemical potential using a nondegenerate expansion.
Definition: fermion.h:324
o2scl::fermion_thermo_tl::type
virtual const char * type()
Return string denoting type ("fermion_thermo")
Definition: fermion.h:834
o2scl::fermion_thermo_tl::pair_density
virtual int pair_density(fermion &f, fp_t temper)=0
Calculate properties with antiparticles as function of density.
o2scl::part_tl< double >::en
double en
Entropy density.
Definition: part.h:120
o2scl::fermion_thermo_tl::massless_root
root * massless_root
A pointer to the solver for massless fermions.
Definition: fermion.h:884
o2scl::fermion_thermo_tl::massless_pair_density
virtual void massless_pair_density(fermion &f, fp_t temper)
Finite temperature massless fermions and antifermions.
Definition: fermion.h:772
o2scl::root
o2scl::fermion_thermo_tl::set_massless_root
void set_massless_root(root<> &rp)
Set the solver for use in massless_calc_density()
Definition: fermion.h:821
o2scl::part_tl< double >::g
double g
Degeneracy (e.g. spin and color if applicable)
Definition: part.h:108
o2scl::fermion_thermo_tl::calc_mu_ndeg
virtual bool calc_mu_ndeg(fermion &f, fp_t temper, fp_t prec=1.0e-18, bool inc_antip=false)
Non-degenerate expansion for fermions.
Definition: fermion.h:599
o2scl::part_tl< double >::ed
double ed
Energy density.
Definition: part.h:114
o2scl::fermion_zerot_tl
Fermion properties at zero temperature.
Definition: fermion.h:131
o2scl::fermion_tl
Fermion class.
Definition: fermion.h:55
exc_einval
exc_einval
o2scl::fermion_zerot_tl::kf_from_density
void kf_from_density(fermion_tl< fp_t > &f)
Calculate the Fermi momentum from the density.
Definition: fermion.h:157
o2scl::part_tl< double >::inc_rest_mass
bool inc_rest_mass
If true, include the mass in the energy density and chemical potential (default true)
Definition: part.h:128
o2scl::fermion_thermo_tl::calc_density
virtual int calc_density(fermion &f, fp_t temper)=0
Calculate properties as function of density.
o2scl::part_tl< double >::m
double m
Mass.
Definition: part.h:110
o2scl::fermion_zerot_tl::pressure_zerot
void pressure_zerot(fermion_tl< fp_t > &f)
Pressure at T=0 from o2scl::fermion_tl::kf and o2scl::part_tl::ms.
Definition: fermion.h:196
o2scl::fermion_thermo_tl::massless_calc_mu
virtual void massless_calc_mu(fermion &f, fp_t temper)
Finite temperature massless fermions.
Definition: fermion.h:653
o2scl::part_tl< double >::pr
double pr
Pressure.
Definition: part.h:116
o2scl::fermion_tl::kf
fp_t kf
Fermi momentum.
Definition: fermion.h:60
o2scl::fermion_tl::type
virtual const char * type()
Return string denoting type ("fermion_tl")
Definition: fermion.h:74
o2scl::root_cern
o2scl::fermion_thermo_tl::calc_mu_deg
virtual bool calc_mu_deg(fermion &f, fp_t temper, fp_t prec=1.0e-18)
Degenerate expansion for fermions.
Definition: fermion.h:621
o2scl::fermion_tl::del
fp_t del
Gap.
Definition: fermion.h:62
O2SCL_ERR
#define O2SCL_ERR(d, n)
o2scl::funct
std::function< double(double)> funct
o2scl::fermion_tl::fermion_tl
fermion_tl(fp_t mass=0, fp_t dof=0)
Create a fermion with mass mass and degeneracy dof.
Definition: fermion.h:65
o2scl::fermion_thermo_tl::ndeg_terms
void ndeg_terms(size_t j, fp_t tt, fp_t xx, fp_t m, bool inc_rest_mass, bool inc_antip, fp_t &pterm, fp_t &nterm, fp_t &enterm)
Compute a term in the nondegenerate expansion.
Definition: fermion.h:838
o2scl::fermion_thermo_tl
Fermion with finite-temperature thermodynamics [abstract base].
Definition: fermion.h:309
o2scl::fermion_zerot_tl::calc_mu_zerot
virtual void calc_mu_zerot(fermion_tl< fp_t > &f)
Zero temperature fermions from o2scl::part_tl::mu or o2scl::part_tl::nu and o2scl::part_tl::ms.
Definition: fermion.h:217
o2scl::fermion_thermo_tl::massless_solve_fun
fp_t massless_solve_fun(fp_t x, fermion &f, fp_t temper)
Solve for the chemical potential for massless fermions.
Definition: fermion.h:887
o2scl::fermion_thermo_tl::massless_pair_mu
virtual void massless_pair_mu(fermion &f, fp_t temper)
Finite temperature massless fermions and antifermions.
Definition: fermion.h:692
o2scl::fermion_thermo_tl::def_massless_root
root_cern def_massless_root
The default solver for massless_calc_density()
Definition: fermion.h:831
o2scl::part_tl< double >::nu
double nu
Effective chemical potential.
Definition: part.h:124
psi
const double psi
o2scl::fermion_zerot_tl::pi2
fp_t pi2
Desc.
Definition: fermion.h:139
o2scl::fermion_tl::fermion_tl
fermion_tl(const fermion_tl &f)
Copy constructor.
Definition: fermion.h:77

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