23 #ifndef O2SCL_FERMION_H
24 #define O2SCL_FERMION_H
34 #include <boost/math/constants/constants.hpp>
35 #include <boost/math/special_functions/bessel.hpp>
38 #include <gsl/gsl_specfunc.h>
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>
46 #ifndef DOXYGEN_NO_O2NS
74 virtual const char *
type() {
return "fermion_tl"; }
115 typedef fermion_tl<double> fermion;
144 pi=boost::math::constants::pi<fp_t>();
145 pi2=boost::math::constants::pi_sqr<fp_t>();
177 efs=gsl_hypot(f.
kf,f.
ms);
180 -pow(f.
ms,4.0)*log(r));
202 efs=gsl_hypot(f.
kf,f.
ms);
204 f.
pr=f.
g/48.0/
pi2*(2.0*efs*pow(f.
kf,3.0)-
206 +3.0*pow(f.
ms,4.0)*log(r));
233 f.
kf=sqrt(nupm*nupm-f.
ms*f.
ms);
244 if (nulessthan0==
true) {
259 f.
nu=gsl_hypot(f.
kf,f.
ms);
308 template<
class fp_t=
double>
323 template<
class fermion_t>
325 fp_t prec,
bool inc_antip) {
327 if (f.non_interacting==
true) { f.nu=f.mu; f.ms=f.m; }
331 if (f.inc_rest_mass) {
334 psi_num=f.nu+f.m-f.ms;
340 if (inc_antip==
false &&
psi>-1.0)
return false;
343 fp_t prefac=f.g/2.0/this->
pi2*pow(f.ms,4.0);
347 static const size_t max_term=200;
353 if (
psi+1.0/tt<-700.0) {
366 fp_t dj1=((fp_t)max_term), jot1=max_term/tt;
367 fp_t dj2=1.0, jot2=1.0/tt;
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);
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);
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);
394 if (std::isfinite(rat) && rat>prec) {
403 for(
size_t j=1;j<=max_term;j++) {
405 fp_t pterm, nterm, enterm;
410 if (j==1) first_term=pterm;
417 if (first_term==0.0) {
427 if (j>1 && fabs(pterm)<prec*fabs(first_term)) {
431 f.ed=-f.pr+f.nu*f.n+temper*f.en;
445 template<
class fermion_t>
456 if (temper<0.0 || f.ms<0.0) {
457 O2SCL_ERR2(
"Temperature or mass negative in fermion_thermo",
461 if (f.non_interacting==
true) { f.nu=f.mu; f.ms=f.m; }
465 if (f.inc_rest_mass)
psi=(f.nu-f.ms)/temper;
466 else psi=(f.nu+f.m-f.ms)/temper;
471 if (
psi<0.0)
return false;
474 fp_t prefac=f.g/2.0/this->
pi2*pow(f.ms,4.0);
479 fp_t s2x=sqrt(2.0+x);
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;
490 pterm1=x2*sx*(29568.0+15840.0*x+1540.0*x2-105.0*x3)/55440.0/sqrt(2.0);
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);
496 if (fabs(pterm4)/fabs(pterm1)>prec) {
501 fp_t nterm1=sx*s2x*x*(2.0+x)/3.0/f.ms;
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)/
507 fp_t enterm2=tt*this->pi2/3.0*(1.0+x)*sx*s2x/f.ms;
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);
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);
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;
600 fp_t prec=1.0e-18,
bool inc_antip=
false) {
601 return calc_mu_ndeg_tlate<fermion>(f,temper,prec,inc_antip);
623 return calc_mu_deg_tlate<fermion>(f,temper,prec);
659 fm2=gsl_sf_fermi_dirac_int(2,f.
nu/temper);
660 fm3=gsl_sf_fermi_dirac_int(3,f.
nu/temper);
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;
675 funct mf2=std::bind(std::mem_fn<fp_t(fp_t,
fermion &,fp_t)>
677 this,std::placeholders::_1,std::ref(f),temper);
693 fp_t pitmu, pitmu2, nu2;
698 f.
ed=f.
g/8.0/this->
pi2*7.0/15.0*
699 pow(this->
pi*temper,4.0);
704 pitmu=this->
pi*temper/f.
nu;
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);
774 fp_t t2=temper*temper,pitmu,pitmu2,nu2;
775 fp_t cbt, alpha, two13, alpha16;
780 f.
ed=f.
g/8.0/this->
pi2*7.0/15.0*pow(this->
pi*temper,4.0);
783 alpha=f.
g*f.
g*this->
pi2*t2*t2*t2/243.0/f.
n/f.
n;
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) {
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);
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);
799 pitmu=this->
pi*temper/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);
806 if (!std::isfinite(f.
nu)) {
807 std::string str=
"Chemical potential not finite ("+
dtos(f.
nu)+
808 ") in fermion::massless_pair_density().";
834 virtual const char *
type() {
return "fermion_thermo"; }
839 fp_t xx, fp_t m,
bool inc_rest_mass,
840 bool inc_antip, fp_t &pterm, fp_t &nterm,
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;
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;
857 enterm=enterm1-enterm2;
859 enterm=enterm1+enterm2;
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;
870 enterm=enterm1-enterm2;
872 enterm=enterm1+enterm2;
879 #ifndef DOXYGEN_NO_O2NS
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;
900 #ifndef DOXYGEN_NO_O2NS