fermion_deriv_rel.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_DERIV_REL_H
24 #define O2SCL_FERMION_DERIV_REL_H
25 
26 /** \file fermion_deriv_rel.h
27  \brief File defining \ref o2scl::fermion_deriv_rel_tl
28 */
29 
30 #include <string>
31 #include <iostream>
32 #include <fstream>
33 #include <cmath>
34 #include <o2scl/constants.h>
35 #include <o2scl/root_cern.h>
36 #include <o2scl/inte.h>
37 #include <o2scl/inte_qag_gsl.h>
38 #include <o2scl/inte_qagiu_gsl.h>
39 
40 #include <o2scl/part_deriv.h>
41 #include <o2scl/fermion_rel.h>
42 
43 #ifndef DOXYGEN_NO_O2NS
44 namespace o2scl {
45 #endif
46 
47  /** \brief Equation of state for a relativistic fermion
48 
49  \note This class only has preliminary support for
50  inc_rest_mass=true (more testing should be done, particularly
51  for the "pair" functions)
52 
53  This implements an equation of state for a relativistic fermion
54  using direct integration. After subtracting the rest mass from
55  the chemical potentials, the distribution function is
56  \f[
57  \left\{1+\exp[(\sqrt{k^2+m^{* 2}}-m-\nu)/T]\right\}^{-1}
58  \f]
59  where \f$ k \f$ is the momentum, \f$ \nu \f$ is the effective
60  chemical potential, \f$ m \f$ is the rest mass, and \f$ m^{*}
61  \f$ is the effective mass. For later use, we define \f$ E^{*} =
62  \sqrt{k^2 + m^{*2}} \f$ . The degeneracy parameter is
63  \f[
64  \psi=(\nu+(m-m^{*}))/T
65  \f]
66  For \f$ \psi \f$ greater than \ref deg_limit (degenerate
67  regime), a finite interval integrator is used and for \f$ \psi
68  \f$ less than \ref deg_limit (non-degenerate regime), an
69  integrator over the interval from \f$ [0,\infty) \f$ is
70  used. The upper limit on the degenerate integration is given by
71  the value of the momentum \f$ k \f$ which is the solution of
72  \f[
73  (\sqrt{k^2+m^{*,2}}-m-\nu)/T=\mathrm{f{l}imit}
74  \f]
75  which is
76  \f[
77  \sqrt{(m+{\cal L})^2-m^{*2}}
78  \f]
79  where \f$ {\cal L}\equiv\mathrm{f{l}imit}\times T+\nu \f$ .
80 
81  For the entropy integration, we set the lower limit
82  to
83  \f[
84  2 \sqrt{\nu^2+2 \nu m} - \mathrm{upper~limit}
85  \f]
86  since the only contribution to the entropy is at the Fermi surface.
87  \comment
88  I'm not sure, but I think this is an expression determined
89  from a small T taylor expansion of the argument of the
90  exponential.
91  \endcomment
92 
93  In the non-degenerate regime, we make the substitution \f$ u=k/T
94  \f$ to help ensure that the variable of integration scales
95  properly.
96 
97  Uncertainties are given in \ref unc.
98 
99  \b Evaluation \b of \b the \b derivatives
100 
101  The relevant
102  derivatives of the distribution function are
103  \f[
104  \frac{\partial f}{\partial T}=
105  f(1-f)\frac{E^{*}-m-\nu}{T^2}
106  \f]
107  \f[
108  \frac{\partial f}{\partial \nu}=
109  f(1-f)\frac{1}{T}
110  \f]
111  \f[
112  \frac{\partial f}{\partial k}=
113  -f(1-f)\frac{k}{E^{*} T}
114  \f]
115  \f[
116  \frac{\partial f}{\partial m^{*}}=
117  -f(1-f)\frac{m^{*}}{E^{*} T}
118  \f]
119 
120  We also need the derivative of the entropy integrand w.r.t. the
121  distribution function, which is
122  \f[
123  {\cal S}\equiv f \ln f +(1-f) \ln (1-f) \qquad
124  \frac{\partial {\cal S}}{\partial f} = \ln
125  \left(\frac{f}{1-f}\right) =
126  \left(\frac{\nu-E^{*}+m}{T}\right)
127  \f]
128  where the entropy density is
129  \f[
130  s = - \frac{g}{2 \pi^2} \int_0^{\infty} {\cal S} k^2 d k
131  \f]
132 
133  The derivatives can be integrated directly (\ref method = \ref
134  direct) or they may be converted to integrals over the
135  distribution function through an integration by parts (\ref
136  method = \ref by_parts)
137  \f[
138  \int_a^b f(k) \frac{d g(k)}{dk} dk = \left.f(k) g(k)\right|_{k=a}^{k=b}
139  - \int_a^b g(k) \frac{d f(k)}{dk} dk
140  \f]
141  using the distribution function for \f$ f(k) \f$ and 0 and
142  \f$ \infty \f$ as the limits, we have
143  \f[
144  \frac{g}{2 \pi^2} \int_0^{\infty} \frac{d g(k)}{dk} f dk =
145  \frac{g}{2 \pi^2} \int_0^{\infty} g(k) f (1-f) \frac{k}{E^{*} T} dk
146  \f]
147  as long as \f$ g(k) \f$ vanishes at \f$ k=0 \f$ .
148  Rewriting,
149  \f[
150  \frac{g}{2 \pi^2} \int_0^{\infty} h(k) f (1-f) dk =
151  \frac{g}{2 \pi^2} \int_0^{\infty} f \frac{T}{k}
152  \left[ h^{\prime} E^{*}-\frac{h E^{*}}{k}+\frac{h k}{E^{*}} \right] dk
153  \f]
154  as long as \f$ h(k)/k \f$ vanishes at \f$ k=0 \f$ .
155 
156  \b Explicit \b forms
157 
158  1) The derivative of the density wrt the chemical potential
159  \f[
160  \left(\frac{d n}{d \mu}\right)_T =
161  \frac{g}{2 \pi^2} \int_0^{\infty} \frac{k^2}{T} f (1-f) dk
162  \f]
163  Using \f$ h(k)=k^2/T \f$ we get
164  \f[
165  \left(\frac{d n}{d \mu}\right)_T =
166  \frac{g}{2 \pi^2} \int_0^{\infty}
167  \left(\frac{k^2+E^{*2}}{E^{*}}\right) f dk
168  \f]
169 
170  2) The derivative of the density wrt the temperature
171  \f[
172  \left(\frac{d n}{d T}\right)_{\mu} =
173  \frac{g}{2 \pi^2} \int_0^{\infty} \frac{k^2(E^{*}-m-\nu)}{T^2}
174  f (1-f) dk
175  \f]
176  Using \f$ h(k)=k^2(E^{*}-\nu)/T^2 \f$ we get
177  \f[
178  \left(\frac{d n}{d T}\right)_{\mu} =
179  \frac{g}{2 \pi^2} \int_0^{\infty} \frac{f}{T}
180  \left[2 k^2+E^{*2}-E^{*}\left(\nu+m\right)-
181  k^2 \left(\frac{\nu+m}{E^{*}}\right)\right] dk
182  \f]
183 
184  3) The derivative of the entropy wrt the chemical potential
185  \f[
186  \left(\frac{d s}{d \mu}\right)_T =
187  \frac{g}{2 \pi^2} \int_0^{\infty} k^2 f (1-f)
188  \frac{(E^{*}-m-\nu)}{T^2} dk
189  \f]
190  This verifies the Maxwell relation
191  \f[
192  \left(\frac{d s}{d \mu}\right)_T =
193  \left(\frac{d n}{d T}\right)_{\mu}
194  \f]
195 
196  4) The derivative of the entropy wrt the temperature
197  \f[
198  \left(\frac{d s}{d T}\right)_{\mu} =
199  \frac{g}{2 \pi^2} \int_0^{\infty} k^2 f (1-f)
200  \frac{(E^{*}-m-\nu)^2}{T^3} dk
201  \f]
202  Using \f$ h(k)=k^2 (E^{*}-\nu)^2/T^3 \f$
203  \f[
204  \left(\frac{d s}{d T}\right)_{\mu} =
205  \frac{g}{2 \pi^2} \int_0^{\infty} \frac{f(E^{*}-m-\nu)}{E^{*}T^2}
206  \left[E^{* 3}+3 E^{*} k^2- (E^{* 2}+k^2)(\nu+m)\right] d k
207  \f]
208 
209  5) The derivative of the density wrt the effective mass
210  \f[
211  \left(\frac{d n}{d m^{*}}\right)_{T,\mu} =
212  -\frac{g}{2 \pi^2} \int_0^{\infty}
213  \frac{k^2 m^{*}}{E^{*} T} f (1-f) dk
214  \f]
215  Using \f$ h(k)=-(k^2 m^{*})/(E^{*} T) \f$ we get
216  \f[
217  \left(\frac{d n}{d m^{*}}\right)_{T,\mu} =
218  -\frac{g}{2 \pi^2} \int_0^{\infty}
219  m^{*} f dk
220  \f]
221  \comment
222  This derivative may be written in terms of the
223  others
224  \f[
225  \left(\frac{d n}{d m^{*}}\right)_{T,\mu} = \frac{3 n}{m^{*}}
226  - \frac{T}{m^{*}}\left[ \left(\frac{d n}{d T}\right)_{\mu}
227  +\frac{\mu}{T} \left(\frac{d n}{d \mu}\right)_{T}
228  \right] - \left(\frac{d n}{d \mu}\right)_{T}
229  \f]
230  \endcomment
231 
232  \note The dsdT integration may fail if the system is
233  very degenerate. When method is byparts, the integral involves a
234  large cancellation between the regions from \f$ k \in (0,
235  \mathrm{ulimit/2}) \f$ and \f$ k \in (\mathrm{ulimit/2},
236  \mathrm{ulimit}) \f$. Switching to method=direct and setting the
237  lower limit to \f$ \mathrm{llimit} \f$, may help, but recent
238  testing on this gave negative values for dsdT. For very
239  degenerate systems, an expansion may be better than trying
240  to perform the integration. The value of the integrand
241  at k=0 also looks like it might be causing difficulties.
242 
243  \future The option err_nonconv=false is not really implemented
244  yet.
245 
246  \future The \ref pair_density() function is a bit slow because
247  it computes the non-derivative thermodynamic quantities
248  twice, and this could be improved.
249  */
250  template<class fp_t=double>
252 
253  public:
254 
255  /// Create a fermion with mass \c m and degeneracy \c g
257 
258  deg_limit=2.0;
259  upper_limit_fac=20.0;
260 
262  nit=&def_nit;
263  dit=&def_dit;
264 
267 
268  exp_limit=200.0;
269 
270  err_nonconv=true;
271 
272  last_method=0;
273  }
274 
275  virtual ~fermion_deriv_rel_tl() {
276  }
277 
278  /** \brief Limit of arguments of exponentials for Fermi functions
279  (default 200.0)
280  */
281  fp_t exp_limit;
282 
283  /** \brief The critical degeneracy at which to switch integration
284  techniques (default 2.0)
285  */
286  fp_t deg_limit;
287 
288  /** \brief The limit for the Fermi functions (default 20.0)
289 
290  fermion_deriv_rel will ignore corrections smaller than about
291  \f$ \exp(-\mathrm{f{l}imit}) \f$ .
292  */
294 
295  /// Storage for the most recently calculated uncertainties
297 
298  /// Object for computing non-derivative quantities
300 
301  /** \name Method of computing derivatives
302  */
303  //@{
304  /// Method (default is \ref automatic)
305  int method;
306  /// Automatically choose method
307  static const int automatic=0;
308  /// In the form containing \f$ f(1-f) \f$ .
309  static const int direct=1;
310  /// Integrate by parts
311  static const int by_parts=2;
312  //@}
313 
314  /** \brief An integer indicating the last numerical method used
315 
316  The function \ref calc_mu() sets this integer to a two-digit
317  number. It is equal to 10 times the value reported by \ref
318  o2scl::fermion_rel::calc_mu_tlate() plus a value from the list
319  below corresponding to the method used for the derivatives
320  - 1: nondegenerate expansion
321  - 2: degenerate expansion
322  - 3: nondegenerate integrand, using \ref by_parts for
323  \ref method
324  - 4: nondegenerate integrand, using user-specified value
325  for \ref method
326  - 5: degenerate integrand, using \ref direct
327  - 6: degenerate integrand, using \ref by_parts
328  - 7: degenerate integrand, using user-specified
329  value for \ref method
330 
331  The function \ref nu_from_n() sets this value equal to
332  100 times the value reported by
333  \ref o2scl::fermion_rel_tl::nu_from_n_tlate() .
334 
335  The function \ref calc_density() sets this value equal to the
336  value from \ref o2scl::fermion_deriv_rel_tl::nu_from_n() plus the
337  value from \ref o2scl::fermion_deriv_rel_tl::calc_mu() .
338 
339  */
341 
342  /** \brief If true, call the error handler when convergence
343  fails (default true)
344  */
346 
347  /** \brief Calculate properties as function of chemical potential
348  */
349  virtual int calc_mu(fermion_deriv &f, fp_t temper) {
350 
351  fr.calc_mu_tlate<fermion_deriv>(f,temper);
353 
354  int iret;
355 
356  if (temper<=0.0) {
357  O2SCL_ERR("T=0 not implemented in fermion_deriv_rel().",exc_eunimpl);
358  }
359 
360  if (f.non_interacting==true) { f.nu=f.mu; f.ms=f.m; }
361 
362  fp_t prefac=f.g/2.0/this->pi2;
363 
364  // Compute the degeneracy parameter
365 
366  bool deg=false;
367  fp_t psi;
368  if (f.inc_rest_mass) psi=(f.nu-f.ms)/temper;
369  else psi=(f.nu+f.m-f.ms)/temper;
370  if (psi>deg_limit) deg=true;
371 
372  // Try the non-degenerate expansion if psi is small enough
373  if (psi<-4.0) {
374  bool acc=this->calc_mu_ndeg(f,temper,1.0e-14);
375  if (acc) {
376  unc.n=f.n*1.0e-14;
377  unc.ed=f.ed*1.0e-14;
378  unc.pr=f.pr*1.0e-14;
379  unc.en=f.en*1.0e-14;
380  unc.dndT=f.dndT*1.0e-14;
381  unc.dsdT=f.dsdT*1.0e-14;
382  unc.dndmu=f.dndmu*1.0e-14;
383  last_method+=1;
384  return 0;
385  }
386  }
387 
388  // Try the degenerate expansion if psi is large enough
389  if (psi>20.0) {
390  bool acc=this->calc_mu_deg(f,temper,1.0e-14);
391  if (acc) {
392  unc.n=f.n*1.0e-14;
393  unc.ed=f.ed*1.0e-14;
394  unc.pr=f.pr*1.0e-14;
395  unc.en=f.en*1.0e-14;
396  unc.dndT=f.dndT*1.0e-14;
397  unc.dsdT=f.dsdT*1.0e-14;
398  unc.dndmu=f.dndmu*1.0e-14;
399  last_method+=2;
400  return 0;
401  }
402  }
403 
404  if (deg==false) {
405 
406  // Set integration method
407  if (method==automatic) {
409  last_method+=3;
410  } else {
412  last_method+=4;
413  }
414 
415  // The non-degenerate case
416 
417  funct density_T_fun_f=
418  std::bind(std::mem_fn<fp_t(fp_t,fermion_deriv &,fp_t)>
420  this,std::placeholders::_1,std::ref(f),temper);
421  iret=nit->integ_err(density_T_fun_f,0.0,0.0,f.dndT,unc.dndT);
422  if (iret!=0) {
423  O2SCL_ERR2("dndT integration (ndeg) failed in ",
424  "fermion_deriv_rel::calc_mu().",
425  exc_efailed);
426  }
427  f.dndT*=prefac;
428  unc.dndT*=prefac;
429 
430  funct density_mu_fun_f=
431  std::bind(std::mem_fn<fp_t(fp_t,fermion_deriv &,fp_t)>
433  this,std::placeholders::_1,std::ref(f),temper);
434  iret=nit->integ_err(density_mu_fun_f,0.0,0.0,f.dndmu,unc.dndmu);
435  if (iret!=0) {
436  O2SCL_ERR2("dndmu integration (ndeg) failed in ",
437  "fermion_deriv_rel::calc_mu().",
438  exc_efailed);
439  }
440  f.dndmu*=prefac;
441  unc.dndmu*=prefac;
442 
443  funct entropy_T_fun_f=
444  std::bind(std::mem_fn<fp_t(fp_t,fermion_deriv &,fp_t)>
446  this,std::placeholders::_1,std::ref(f),temper);
447  iret=nit->integ_err(entropy_T_fun_f,0.0,0.0,f.dsdT,unc.dsdT);
448  if (iret!=0) {
449  O2SCL_ERR2("dsdT integration (ndeg) failed in ",
450  "fermion_deriv_rel_tl<fp_t>::calc_mu().",exc_efailed);
451  }
452  f.dsdT*=prefac;
453  unc.dsdT*=prefac;
454 
455  } else {
456 
457  // Compute the upper limit for degenerate integrals
458 
459  fp_t arg;
460  if (f.inc_rest_mass) {
461  arg=pow(upper_limit_fac*temper+f.nu,2.0)-f.ms*f.ms;
462  } else {
463  arg=pow(upper_limit_fac*temper+f.nu+f.m,2.0)-f.ms*f.ms;
464  }
465  // Set to zero to avoid uninit'ed var. warnings
466  fp_t ul=0.0;
467  if (arg>0.0) {
468  ul=sqrt(arg);
469  } else {
470  O2SCL_ERR2("Zero density in degenerate limit in fermion_deriv_rel::",
471  "calc_mu(). Variable deg_limit set improperly?",
472  exc_efailed);
473  }
474 
475  // Compute the lower limit for the entropy and derivative
476  // integrations
477 
478  fp_t ll;
479  if (f.inc_rest_mass) {
480  arg=pow(-upper_limit_fac*temper+f.nu,2.0)-f.ms*f.ms;
481  if (arg>0.0 && (f.ms-f.nu)/temper<-upper_limit_fac) {
482  ll=sqrt(arg);
483  } else {
484  ll=-1.0;
485  }
486  } else {
487  arg=pow(-upper_limit_fac*temper+f.nu+f.m,2.0)-f.ms*f.ms;
488  if (arg>0.0 && (f.ms-f.nu-f.m)/temper<-upper_limit_fac) {
489  ll=sqrt(arg);
490  } else {
491  ll=-1.0;
492  }
493  }
494 
495  // Set integration method
496  if (method==automatic) {
497  if ((!f.inc_rest_mass && (f.nu+f.m-f.ms)/temper>1.0e3) ||
498  (f.inc_rest_mass && (f.nu-f.ms)/temper>1.0e3)) {
500  last_method+=5;
501  } else {
503  last_method+=6;
504  }
505  } else {
507  last_method+=7;
508  }
509 
510  funct deg_density_mu_fun_f=
511  std::bind(std::mem_fn<fp_t(fp_t,fermion_deriv &,fp_t)>
513  this,std::placeholders::_1,std::ref(f),temper);
514  if (intl_method==direct && ll>0.0) {
515  iret=dit->integ_err(deg_density_mu_fun_f,ll,ul,
516  f.dndmu,unc.dndmu);
517  } else {
518  iret=dit->integ_err(deg_density_mu_fun_f,0.0,ul,
519  f.dndmu,unc.dndmu);
520  }
521  if (iret!=0) {
522  O2SCL_ERR2("dndmu integration (deg) failed in ",
523  "fermion_deriv_rel_tl<fp_t>::calc_mu().",exc_efailed);
524  }
525  f.dndmu*=prefac;
526  unc.dndmu*=prefac;
527 
528  funct deg_density_T_fun_f=std::bind
529  (std::mem_fn<fp_t(fp_t,fermion_deriv &,fp_t)>
531  this,std::placeholders::_1,std::ref(f),temper);
532  if (intl_method==direct && ll>0.0) {
533  iret=dit->integ_err(deg_density_T_fun_f,ll,ul,f.dndT,unc.dndT);
534  } else {
535  iret=dit->integ_err(deg_density_T_fun_f,0.0,ul,f.dndT,unc.dndT);
536  }
537  if (iret!=0) {
538  O2SCL_ERR2("dndT integration (deg) failed in ",
539  "fermion_deriv_rel_tl<fp_t>::calc_mu().",exc_efailed);
540  }
541  f.dndT*=prefac;
542  unc.dndT*=prefac;
543 
544  funct deg_entropy_T_fun_f=std::bind
545  (std::mem_fn<fp_t(fp_t,fermion_deriv &,fp_t)>
547  this,std::placeholders::_1,std::ref(f),temper);
548  if (intl_method==direct && ll>0.0) {
549  iret=dit->integ_err(deg_entropy_T_fun_f,ll,ul,f.dsdT,unc.dsdT);
550  } else {
551  iret=dit->integ_err(deg_entropy_T_fun_f,0.0,ul,f.dsdT,unc.dsdT);
552  }
553  if (iret!=0) {
554  O2SCL_ERR2("dsdT integration (deg) failed in ",
555  "fermion_deriv_rel_tl<fp_t>::calc_mu().",exc_efailed);
556  }
557  f.dsdT*=prefac;
558  unc.dsdT*=prefac;
559 
560  }
561 
562  if (!std::isfinite(f.en)) {
563  O2SCL_ERR2("Entropy not finite in ",
564  "fermion_deriv_rel_tl<fp_t>::calc_mu().",exc_efailed);
565  }
566  f.pr=-f.ed+temper*f.en+f.nu*f.n;
567  // Pressure uncertainties are not computed
568  unc.pr=0.0;
569 
570  return 0;
571  }
572 
573  /** \brief Calculate properties as function of density
574  */
575  virtual int calc_density(fermion_deriv &f, fp_t temper) {
576 
577  if (f.non_interacting==true) { f.ms=f.m; f.nu=f.mu; }
578 
579  nu_from_n(f,temper);
580  int lm=last_method;
581 
582  if (f.non_interacting) { f.mu=f.nu; }
583 
584  calc_mu(f,temper);
585  last_method+=lm;
586 
587  return 0;
588  }
589 
590  /** \brief Calculate properties with antiparticles as function of
591  chemical potential
592  */
593  virtual int pair_mu(fermion_deriv &f, fp_t temper) {
594  if (f.non_interacting) { f.nu=f.mu; f.ms=f.m; }
595 
596  fermion_deriv antip(f.ms,f.g);
597  f.anti(antip);
598 
599  calc_mu(f,temper);
600  int lm=last_method*100;
601  calc_mu(antip,temper);
602  last_method+=lm;
603 
604  f.n-=antip.n;
605  f.pr+=antip.pr;
606  if (f.inc_rest_mass) {
607  f.ed+=antip.ed;
608  } else {
609  f.ed=f.ed+antip.ed+2.0*antip.n*f.m;
610  }
611  f.en+=antip.en;
612  f.dsdT+=antip.dsdT;
613  f.dndT-=antip.dndT;
614  f.dndmu+=antip.dndmu;
615 
616  return 0;
617  }
618 
619  /** \brief Calculate properties with antiparticles as function of
620  density
621  */
622  virtual int pair_density(fermion_deriv &f, fp_t temper) {
623  int ret=fr.pair_density_tlate<fermion_deriv>(f,temper);
624  pair_mu(f,temper);
625  return 0;
626  }
627 
628  /// Calculate effective chemical potential from density
629  virtual int nu_from_n(fermion_deriv &f, fp_t temper) {
630  int ret=fr.nu_from_n_tlate<fermion_deriv>(f,temper);
632  return ret;
633  }
634 
635  /** \brief Set inte objects
636 
637  The first integrator is used for non-degenerate integration
638  and should integrate from 0 to \f$ \infty \f$ (like \ref
639  o2scl::inte_qagiu_gsl). The second integrator is for the
640  degenerate case, and should integrate between two finite
641  values.
642  */
643  void set_inte(inte<> &unit, inte<> &udit) {
644  nit=&unit;
645  dit=&udit;
646  return;
647  }
648 
649  /** \brief Set the solver for use in calculating the chemical
650  potential from the density */
652  density_root=&rp;
653  return;
654  }
655 
656  /// The default integrator for the non-degenerate regime
658 
659  /// The default integrator for the degenerate regime
661 
662  /// The default solver for npen_density() and pair_density()
664 
665  /// Return string denoting type ("fermion_deriv_rel")
666  virtual const char *type() { return "fermion_deriv_rel"; };
667 
668  protected:
669 
670 #ifndef DOXYGEN_NO_O2NS
671 
672  /// The internal integration method
673  int intl_method;
674 
675  /// The integrator for non-degenerate fermions
677 
678  /// The integrator for degenerate fermions
680 
681  /// The solver for calc_density() and pair_density()
683 
684  /** \name The integrands, as a function of \f$ u=k/T \f$, for
685  non-degenerate integrals
686  */
687  //@{
688  fp_t density_T_fun(fp_t u, fermion_deriv &f, fp_t T) {
689  fp_t k=u*(T), E, ret;
690  if (f.inc_rest_mass) {
691  E=gsl_hypot(k,f.ms);
692  if (intl_method==direct) {
693  fp_t ff=fermi_function(E,f.nu,T,exp_limit);
694  ret=k*k*(E-f.nu)/T*ff*(1.0-ff);
695  } else {
696  ret=(2.0*k*k/T+E*E/T-E*(f.nu)/T-k*k*(f.nu)/T/E)*
697  T*fermi_function(E,f.nu,T,exp_limit);
698  }
699  } else {
700  E=gsl_hypot(k,f.ms);
701  if (intl_method==direct) {
702  E-=f.m;
703  fp_t ff=fermi_function(E,f.nu,T,exp_limit);
704  ret=k*k*(E-f.nu)/T*ff*(1.0-ff);
705  } else {
706  ret=(2.0*k*k/T+E*E/T-E*(f.nu+f.m)/T-k*k*(f.nu+f.m)/T/E)*
707  T*fermi_function(E-f.m,f.nu,T,exp_limit);
708  }
709  }
710  return ret;
711  }
712 
713  fp_t density_mu_fun(fp_t u, fermion_deriv &f, fp_t T) {
714  fp_t k=u*(T), E, ret;
715  if (f.inc_rest_mass) {
716  E=gsl_hypot(k,f.ms);
717  if (intl_method==direct) {
718  fp_t ff=fermi_function(E,f.nu,T,exp_limit);
719  ret=k*k*ff*(1.0-ff);
720  } else {
721  ret=T*(E*E+k*k)/E*fermi_function(E,f.nu,T,exp_limit);
722  }
723  } else {
724  E=gsl_hypot(k,f.ms);
725  if (intl_method==direct) {
726  E-=f.m;
727  fp_t ff=fermi_function(E,f.nu,T,exp_limit);
728  ret=k*k*ff*(1.0-ff);
729  } else {
730  ret=T*(E*E+k*k)/E*fermi_function(E-f.m,f.nu,T,exp_limit);
731  }
732  }
733  return ret;
734  }
735 
736  fp_t entropy_T_fun(fp_t u, fermion_deriv &f, fp_t T) {
737  fp_t k=u*T, E, ret;
738  if (f.inc_rest_mass) {
739  E=gsl_hypot(k,f.ms);
740  if (intl_method==direct) {
741  fp_t ff=fermi_function(E,f.nu,T,exp_limit);
742  ret=T*k*k*ff*(1.0-ff)*pow(E-f.nu,2.0)/pow(T,3.0);
743  } else {
744  ret=(E-f.nu)/E/T*
745  (pow(E,3.0)+3.0*E*k*k-(E*E+k*k)*(f.nu))*
746  fermi_function(E,f.nu,T,exp_limit);
747  }
748  } else {
749  E=gsl_hypot(k,f.ms);
750  if (intl_method==direct) {
751  E-=f.m;
752  fp_t ff=fermi_function(E,f.nu,T,exp_limit);
753  ret=T*k*k*ff*(1.0-ff)*pow(E-f.nu,2.0)/pow(T,3.0);
754  } else {
755  ret=(E-f.m-f.nu)/E/T*
756  (pow(E,3.0)+3.0*E*k*k-(E*E+k*k)*(f.nu+f.m))*
757  fermi_function(E-f.m,f.nu,T,exp_limit);
758  }
759  }
760  return ret;
761  }
762 
763  fp_t density_ms_fun(fp_t u, fermion_deriv &f, fp_t T) {
764  fp_t k=u*T, E, ret;
765  if (f.inc_rest_mass) {
766  E=gsl_hypot(k,f.ms);
767  if (intl_method==direct) {
768  fp_t ff=fermi_function(E,f.nu,T,exp_limit);
769  ret=-k*k*f.ms/(E)/T*ff*(1.0-ff);
770  } else {
771  ret=-f.ms*fermi_function(E,f.nu,T,exp_limit);
772  }
773  } else {
774  E=gsl_hypot(k,f.ms);
775  if (intl_method==direct) {
776  E-=f.m;
777  fp_t ff=fermi_function(E,f.nu,T,exp_limit);
778  ret=-k*k*f.ms/(E+f.m)/T*ff*(1.0-ff);
779  } else {
780  ret=-f.ms*fermi_function(E-f.m,f.nu,T,exp_limit);
781  }
782  }
783  ret*=T;
784  return ret;
785  }
786 
787  //@}
788 
789  /** \name The integrands, as a function of momentum, for the
790  degenerate integrals
791  */
792  //@{
793  fp_t deg_density_T_fun(fp_t k, fermion_deriv &f, fp_t T) {
794  fp_t E, ret;
795  if (f.inc_rest_mass) {
796  E=gsl_hypot(k,f.ms);
797  if (intl_method==direct) {
798  fp_t ff=fermi_function(E,f.nu,T,exp_limit);
799  ret=k*k*(E-f.nu)/T/T*ff*(1.0-ff);
800  } else {
801  ret=(2.0*k*k/T+E*E/T-E*(f.nu)/T-k*k*(f.nu)/T/E)*
802  fermi_function(E,f.nu,T,exp_limit);
803  }
804  } else {
805  E=gsl_hypot(k,f.ms);
806  if (intl_method==direct) {
807  E-=f.m;
808  fp_t ff=fermi_function(E,f.nu,T,exp_limit);
809  ret=k*k*(E-f.nu)/T/T*ff*(1.0-ff);
810  } else {
811  ret=(2.0*k*k/T+E*E/T-E*(f.nu+f.m)/T-k*k*(f.nu+f.m)/T/E)*
812  fermi_function(E-f.m,f.nu,T,exp_limit);
813  }
814  }
815  return ret;
816  }
817 
818  fp_t deg_density_mu_fun(fp_t k, fermion_deriv &f, fp_t T) {
819  fp_t E, ret;
820  if (f.inc_rest_mass) {
821  E=gsl_hypot(k,f.ms);
822  if (intl_method==direct) {
823  fp_t ff=fermi_function(E,f.nu,T,exp_limit);
824  ret=k*k/T*ff*(1.0-ff);
825  } else {
826  ret=(E*E+k*k)/E*fermi_function(E,f.nu,T,exp_limit);
827  }
828  } else {
829  E=gsl_hypot(k,f.ms);
830  if (intl_method==direct) {
831  E-=f.m;
832  fp_t ff=fermi_function(E,f.nu,T,exp_limit);
833  ret=k*k/T*ff*(1.0-ff);
834  } else {
835  ret=(E*E+k*k)/E*fermi_function(E-f.m,f.nu,T,exp_limit);
836  }
837  }
838  return ret;
839  }
840 
841  fp_t deg_entropy_T_fun(fp_t k, fermion_deriv &f, fp_t T) {
842  fp_t E, ret;
843  E=gsl_hypot(k,f.ms);
844  if (f.inc_rest_mass) {
845  fp_t ff=fermi_function(E,f.nu,T,exp_limit);
846  if (intl_method==direct) {
847  ret=k*k*ff*(1.0-ff)*pow(E-f.nu,2.0)/pow(T,3.0);
848  } else {
849  ret=(E-f.nu)/E/T/T*
850  (pow(E,3.0)+3.0*E*k*k-(E*E+k*k)*f.nu)*ff;
851  }
852  } else {
853  fp_t ff=fermi_function(E-f.m,f.nu,T,exp_limit);
854  if (intl_method==direct) {
855  ret=k*k*ff*(1.0-ff)*pow(E-f.nu-f.m,2.0)/pow(T,3.0);
856  } else {
857  ret=(E-f.m-f.nu)/E/T/T*
858  (pow(E,3.0)+3.0*E*k*k-(E*E+k*k)*(f.nu+f.m))*ff;
859  }
860  }
861  return ret;
862  }
863 
864  fp_t deg_density_ms_fun(fp_t k, fermion_deriv &f, fp_t T) {
865  fp_t E, ret;
866  if (f.inc_rest_mass) {
867  E=gsl_hypot(k,f.ms);
868  if (intl_method==direct) {
869  fp_t ff=fermi_function(E,f.nu,T,exp_limit);
870  ret=-k*k*f.ms/E/T*ff*(1.0-ff);
871  } else {
872  ret=-f.ms*fermi_function(E,f.nu,T,exp_limit);
873  }
874  } else {
875  E=gsl_hypot(k,f.ms);
876  if (intl_method==direct) {
877  E-=f.m;
878  fp_t ff=fermi_function(E,f.nu,T,exp_limit);
879  ret=-k*k*f.ms/(E+f.m)/T*ff*(1.0-ff);
880  } else {
881  ret=-f.ms*fermi_function(E-f.m,f.nu,T,exp_limit);
882  }
883  }
884  return ret;
885  }
886 
887  //@}
888 
889 #endif
890 
891  };
892 
893  /** \brief Double-precision version of
894  \ref o2scl::fermion_deriv_rel_tl
895  */
897 
898 #ifndef DOXYGEN_NO_O2NS
899 }
900 #endif
901 
902 #endif
o2scl::part_deriv_press_tl::dsdT
fp_t dsdT
Derivative of entropy density with respect to temperature.
Definition: part_deriv.h:135
o2scl::fermion_deriv_rel_tl::def_dit
inte_qag_gsl def_dit
The default integrator for the degenerate regime.
Definition: fermion_deriv_rel.h:660
o2scl::fermion_deriv_rel_tl::upper_limit_fac
fp_t upper_limit_fac
The limit for the Fermi functions (default 20.0)
Definition: fermion_deriv_rel.h:293
o2scl::fermion_deriv_thermo_tl
Compute properties of a fermion including derivatives [abstract base].
Definition: part_deriv.h:712
o2scl::fermion_rel_tl::calc_mu_tlate
void calc_mu_tlate(fermion_t &f, fp_t temper)
Calculate properties as function of chemical potential (template version)
Definition: fermion_rel.h:541
o2scl::inte_qag_gsl
o2scl::part_deriv_press_tl::dndT
fp_t dndT
Derivative of number density with respect to temperature.
Definition: part_deriv.h:132
o2scl::part_tl::mu
fp_t mu
Chemical potential.
Definition: part.h:118
o2scl::part_tl::ms
fp_t ms
Effective mass (Dirac unless otherwise specified)
Definition: part.h:122
o2scl::fermion_deriv_rel_tl::type
virtual const char * type()
Return string denoting type ("fermion_deriv_rel")
Definition: fermion_deriv_rel.h:666
exc_efailed
exc_efailed
O2SCL_ERR2
#define O2SCL_ERR2(d, d2, n)
o2scl::fermion_deriv_rel_tl::method
int method
Method (default is automatic)
Definition: fermion_deriv_rel.h:305
o2scl::part_tl::n
fp_t n
Number density.
Definition: part.h:112
o2scl::fermion_deriv_thermo_tl< double >::calc_mu_ndeg
virtual bool calc_mu_ndeg(fermion_deriv &f, double temper, double prec, bool inc_antip=false)
Calculate properties as a function of chemical potential using a nondegenerate expansion.
Definition: part_deriv.h:837
o2scl::fermion_deriv_rel_tl::unc
fermion_deriv unc
Storage for the most recently calculated uncertainties.
Definition: fermion_deriv_rel.h:296
o2scl::fermion_deriv_rel_tl::fermion_deriv_rel_tl
fermion_deriv_rel_tl()
Create a fermion with mass m and degeneracy g.
Definition: fermion_deriv_rel.h:256
o2scl::part_tl::non_interacting
bool non_interacting
True if the particle is non-interacting (default true)
Definition: part.h:130
o2scl::fermion_deriv_rel_tl::pair_density
virtual int pair_density(fermion_deriv &f, fp_t temper)
Calculate properties with antiparticles as function of density.
Definition: fermion_deriv_rel.h:622
o2scl::fermion_deriv_rel_tl::intl_method
int intl_method
The internal integration method.
Definition: fermion_deriv_rel.h:666
o2scl::fermion_deriv_rel_tl::density_root
root * density_root
The solver for calc_density() and pair_density()
Definition: fermion_deriv_rel.h:682
o2scl::fermion_rel_tl::pair_density_tlate
int pair_density_tlate(fermion_t &f, fp_t temper)
Calculate thermodynamic properties with antiparticles from the density (template version)
Definition: fermion_rel.h:981
o2scl::fermion_deriv_rel_tl::fr
fermion_rel fr
Object for computing non-derivative quantities.
Definition: fermion_deriv_rel.h:299
exc_eunimpl
exc_eunimpl
o2scl::part_tl::en
fp_t en
Entropy density.
Definition: part.h:120
o2scl::fermion_deriv_rel_tl::nit
inte * nit
The integrator for non-degenerate fermions.
Definition: fermion_deriv_rel.h:676
o2scl::fermion_deriv_rel_tl::nu_from_n
virtual int nu_from_n(fermion_deriv &f, fp_t temper)
Calculate effective chemical potential from density.
Definition: fermion_deriv_rel.h:629
o2scl::fermion_deriv_rel_tl::set_density_root
void set_density_root(root<> &rp)
Set the solver for use in calculating the chemical potential from the density.
Definition: fermion_deriv_rel.h:651
o2scl::root
o2scl::inte_qagiu_gsl
o2scl::fermion_deriv_rel_tl::automatic
static const int automatic
Automatically choose method.
Definition: fermion_deriv_rel.h:307
o2scl::part_tl::g
fp_t g
Degeneracy (e.g. spin and color if applicable)
Definition: part.h:108
o2scl::part_tl::ed
fp_t ed
Energy density.
Definition: part.h:114
o2scl::fermion_rel_tl< double >
o2scl::fermion_deriv_rel_tl::set_inte
void set_inte(inte<> &unit, inte<> &udit)
Set inte objects.
Definition: fermion_deriv_rel.h:643
fermi_function
double fermi_function(double E, double mu, double T, double limit=40.0)
o2scl::fermion_deriv_rel_tl::deg_limit
fp_t deg_limit
The critical degeneracy at which to switch integration techniques (default 2.0)
Definition: fermion_deriv_rel.h:286
o2scl::fermion_deriv_rel_tl::direct
static const int direct
In the form containing .
Definition: fermion_deriv_rel.h:309
o2scl::fermion_deriv_tl< double >
o2scl::part_tl::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::part_deriv_press_tl::dndmu
fp_t dndmu
Derivative of number density with respect to chemical potential.
Definition: part_deriv.h:129
o2scl::fermion_deriv_thermo_tl< double >::pi2
double pi2
Desc.
Definition: part_deriv.h:726
o2scl::fermion_deriv_rel_tl::by_parts
static const int by_parts
Integrate by parts.
Definition: fermion_deriv_rel.h:311
o2scl::part_tl::m
fp_t m
Mass.
Definition: part.h:110
o2scl::fermion_deriv_rel_tl::pair_mu
virtual int pair_mu(fermion_deriv &f, fp_t temper)
Calculate properties with antiparticles as function of chemical potential.
Definition: fermion_deriv_rel.h:593
o2scl::fermion_deriv_rel_tl::calc_mu
virtual int calc_mu(fermion_deriv &f, fp_t temper)
Calculate properties as function of chemical potential.
Definition: fermion_deriv_rel.h:349
o2scl::part_tl::pr
fp_t pr
Pressure.
Definition: part.h:116
o2scl::fermion_rel_tl::last_method
int last_method
An integer indicating the last numerical method used.
Definition: fermion_rel.h:383
o2scl::fermion_rel_tl::nu_from_n_tlate
int nu_from_n_tlate(fermion_t &f, fp_t temper)
Calculate the chemical potential from the density (template version)
Definition: fermion_rel.h:391
o2scl::fermion_deriv_thermo_tl< double >::calc_mu_deg
virtual bool calc_mu_deg(fermion_deriv &f, double temper, double prec)
Calculate properties as a function of chemical potential using a degenerate expansion.
Definition: part_deriv.h:767
o2scl::part_tl::anti
virtual void anti(part_tl &ax)
Make ap an anti-particle with the same mass and degeneracy.
Definition: part.h:205
o2scl::root_cern
o2scl::fermion_deriv_rel_tl::def_density_root
root_cern def_density_root
The default solver for npen_density() and pair_density()
Definition: fermion_deriv_rel.h:663
o2scl::fermion_deriv_rel_tl::def_nit
inte_qagiu_gsl def_nit
The default integrator for the non-degenerate regime.
Definition: fermion_deriv_rel.h:657
o2scl::inte
o2scl::fermion_deriv_rel_tl::err_nonconv
bool err_nonconv
If true, call the error handler when convergence fails (default true)
Definition: fermion_deriv_rel.h:345
o2scl::fermion_deriv_rel_tl
Equation of state for a relativistic fermion.
Definition: fermion_deriv_rel.h:251
O2SCL_ERR
#define O2SCL_ERR(d, n)
o2scl::funct
std::function< double(double)> funct
o2scl::fermion_deriv_rel_tl::calc_density
virtual int calc_density(fermion_deriv &f, fp_t temper)
Calculate properties as function of density.
Definition: fermion_deriv_rel.h:575
o2scl::fermion_deriv_rel_tl::last_method
int last_method
An integer indicating the last numerical method used.
Definition: fermion_deriv_rel.h:340
o2scl::part_tl::nu
fp_t nu
Effective chemical potential.
Definition: part.h:124
o2scl::fermion_deriv_rel_tl::exp_limit
fp_t exp_limit
Limit of arguments of exponentials for Fermi functions (default 200.0)
Definition: fermion_deriv_rel.h:281
psi
const double psi
o2scl::inte::integ_err
virtual int integ_err(func_t &func, fp_t a, fp_t b, fp_t &res, fp_t &err)=0
o2scl::fermion_deriv_rel_tl::dit
inte * dit
The integrator for degenerate fermions.
Definition: fermion_deriv_rel.h:679

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