Class fermion_rel2_tl (o2scl)

O2scl : Class List

template<class fermion_t = fermion_tl<double>, class fd_inte_t = class o2scl::fermi_dirac_integ_gsl, class be_inte_t = o2scl::bessel_K_exp_integ_gsl, class root_t = root_cern<>, class func_t = funct, class fp_t = double>
class fermion_rel2_tl : public o2scl::fermion_thermo_tl<fermion_tl<double>, class o2scl::fermi_dirac_integ_gsl, o2scl::bessel_K_exp_integ_gsl, root_cern<>, funct, double>, private o2scl::fermion_rel2_integ<double>

Equation of state for a relativistic fermion.

This class computes the thermodynamics of a relativistic fermion either as a function of the density or the chemical potential. It employs direct integration, using two different integrators for the degenerate and non-degenerate regimes. The default integrators are inte_qag_gsl (for degenerate fermions) and inte_qagiu_gsl (for non-degenerate fermions). For the functions calc_mu() and calc_density(), if the temperature argument is less than or equal to zero, the functions fermion_zerot::calc_mu_zerot() and fermion_zerot::calc_density_zerot() will be used to compute the result.

Template types:

  • fermion_t: the type of particle object which will be passed to the methods defined by this class

  • fd_inte_t: the integration object for massless fermions

  • be_inte_t: the integration object for the nondegenerate limit

  • inte_t: the generic integration object

  • density_root_t: the solver for fixed densities

  • root_t: the solver for massless fermions

  • func_t: the function object type

  • fp_t: the floating point type

Note that density_root_t and root_t are almost always the same type.

Degeneracy parameter:

Define the degeneracy parameter

\[ \psi=(\nu-m^{*})/T \]
where \( \nu \) is the effective chemical potential (including the rest mass) and \( m^{*} \) is the effective mass. For \( \psi \) smaller than min_psi, the non-degenerate expansion in fermion_thermo::calc_mu_ndeg() is attempted first. If that fails, then integration is used. For \( \psi \) greater than deg_limit (degenerate regime), a finite interval integrator is used and for \( \psi \) less than deg_limit (non-degenerate regime), an integrator over the interval from \( [0,\infty) \) is used. In the case where part::inc_rest_mass is false, the degeneracy parameter is
\[ \psi=(\nu+m-m^{*})/T \]

Integration limits:

The upper limit on the degenerate integration is given by

\[ \mathrm{upper~limit} = \sqrt{{\cal L}^2-m^{*,2}} \]
where \( {\cal L}\equiv u T+\nu \) and \( u \) is fermion_rel2::upper_limit_fac . In the case where part::inc_rest_mass is false, the result is
\[ \mathrm{upper~limit} = \sqrt{(m+{\cal L})^2-m^{*2}} \]

The entropy is only significant at the Fermi surface, thus in the degenerate case, the lower limit of the entropy integral can be given be determined by the value of \( k \) which solves

\[ - u = \frac{\sqrt{k^2+m^{* 2}}-\nu}{T} \]
The solution is
\[ \mathrm{lower~limit} = \sqrt{(-u T+{\nu})^2-m^{*,2}} \]
but this solution is only valid if \( (m^{*}-\nu)/T < -u \). In the case where part::inc_rest_mass is false, the result is
\[ \mathrm{lower~limit} = \sqrt{(-u T + m +\nu)^2-m^{*,2}} \]
which is valid if \( (m^{*}-\nu - m)/T < -u \).

Entropy integrand:

In the degenerate regime, the entropy integrand

\[ - k^2 \left[ f \log f + \left(1-f\right) \log \left(1-f \right) \right] \]
where \( f \) is the fermionic distribution function can lose precision when \( (E^{*} - \nu)/T \) is negative and sufficiently large in absolute magnitude. Thus when \( (E^{*} - \nu)/T < S \) where \( S \) is stored in deg_entropy_fac (default is -30), the integrand is written as
\[ -k^2 \left( E/T-\nu/T \right) e^{E/T-\nu/T} \, . \]
If \( (E - \nu)/T < S \) is less than -1 times exp_limit (e.g. less than -200), then the entropy integrand is assumed to be zero.

Non-degenerate integrands:

The integrands in the non-degenerate regime are written in a dimensionless form, by defining \( u \) with the relation \( p = \sqrt{\left(T u + m^{*}\right)^2-m^{* 2}} \), \( y \equiv \nu/ T \), and \( \eta \equiv m^{*}/T \). Then, \( p/T=\sqrt{u^2+2 \eta u} \) \( E/T = \mathrm{mx+u} \) and \( p/T^2 dp = 2(\eta+u) du \) The density integrand is

\[ \left(\eta+u\right) \sqrt{u^2+2 (\eta) u} \left(\frac{e^{y}}{e^{\eta+u}+e^{y}}\right) \, , \]
the energy integrand is
\[ \left(\eta+u\right)^2 \sqrt{u^2+2 (\eta) u} \left(\frac{e^{y}}{e^{\eta+u}+e^{y}}\right) \, , \]
and the entropy integrand is
\[ \left(\eta+u\right) \sqrt{u^2+2 (\eta) u} \left(t_1+t_2\right) \, , \]
where
\[\begin{split}\begin{eqnarray*} t_1 &=& \log \left(1+e^{y-\eta-u}\right)/ \left(1+e^{y-\eta-u}\right) \nonumber \\ t_2 &=& \log \left(1+e^{\eta+u-y}\right)/ \left(1+e^{\eta+u-y}\right) \, . \end{eqnarray*}\end{split}\]

Accuracy:

The default settings for for this class give an accuracy of at least 1 part in \( 10^6 \) (and frequently better than this).

When the integrators provide numerical uncertainties, these uncertainties are stored in unc. In the case of calc_density() and pair_density(), the uncertainty from the numerical accuracy of the solver is not included. There is also a relatively small inaccuracy due to the mathematical evaluation of the integrands which is not included in unc.)

One can improve the accuracy to within 1 part in \( 10^{10} \) using

 fermion_rel2 rf(1.0,2.0); rf.upper_limit_fac=40.0;
rf.dit->tol_abs=1.0e-13; rf.dit->tol_rel=1.0e-13;
rf.nit->tol_abs=1.0e-13; rf.nit->tol_rel=1.0e-13;
rf.density_root->tol_rel=1.0e-10; 
which decreases the both the relative and absolute tolerances for both the degenerate and non-degenerate integrators and improves the accuracy of the solver which determines the chemical potential from the density. Of course if these tolerances are too small, the calculation may fail.

Todo

In class fermion_rel_tl:

  • Future: I had to remove the shared_ptr stuff because the default algorithm types don’t support multiprecision, but it might be nice to restore the shared_ptr mechanism somehow.

  • Future: The expressions which appear in in the integrand functions density_fun(), etc. could likely be improved, especially in the case where ref o2scl::part::inc_rest_mass is <tt>false</tt>. There should not be a need to check if <tt>ret</tt> is finite.

  • Future: It appears this class doesn’t compute the uncertainty in the chemical potential or density with calc_density(). This could be fixed.

  • Future: I’d like to change the lower limit on the entropy integration, but the value in the code at the moment (stored in <tt>ll</tt>) makes bm_part2.cpp worse.

  • Future: The function pair_mu() should set the antiparticle integrators as done in fermion_deriv_rel.

Numerical parameters

bool err_nonconv

If true, call the error handler when convergence fails (default true)

fp_t min_psi

The smallest value of \( (\mu-m)/T \) for which integration is used (default -4.0)

fp_t deg_limit

The critical degeneracy at which to switch integration techniques (default 2.0)

fp_t upper_limit_fac

The factor for the degenerate upper limits (default 20.0)

fp_t deg_entropy_fac

A factor for the degenerate entropy integration (default 30.0)

int verbose

Verbosity parameter (default 0)

bool use_expansions

If true, use expansions for extreme conditions (default true)

fp_t tol_expan

Tolerance for expansions (default \( 10^{-14} \))

bool verify_ti

If true, verify the thermodynamic identity (default false)

fp_t therm_ident

Value for verifying the thermodynamic identity.

o2scl::root_brent_gsl<func_t, fp_t> alt_solver

Alternate solver.

fermion_t unc

Storage for the uncertainty.

root<func_t, func_t, fp_t> *density_root

The solver for calc_density()

density_root_t def_density_root

The default solver for the chemical potential given the density.

int last_method

An integer indicating the last numerical method used.

In all functions

  • 0: no previous calculation or last calculation failed

In nu_from_n():

  • 1: default solver

  • 2: default solver with smaller tolerances

  • 3: bracketing solver

In calc_mu():

  • 4: non-degenerate expansion

  • 5: degenerate expansion

  • 6: exact integration, non-degenerate integrands

  • 7: exact integration, degenerate integrands, lower limit on entropy integration

  • 8: exact integration, degenerate integrands, full entropy integration

  • 9: T=0 result

In calc_density(), the integer is a two-digit number. The first digit (1 to 3) is the method used by nu_from_n() and the second digit is one of

  • 1: nondegenerate expansion

  • 2: degenerate expansion

  • 3: exact integration, non-degenerate integrands

  • 4: exact integration, degenerate integrands, lower limit on entropy integration

  • 5: exact integration, degenerate integrands, full entropy integration If calc_density() uses the T=0 code, then last_method is 40.

In pair_mu(), the integer is a three-digit number. The third digit is always 0 (to ensure a value of last_method which is unique from the other values reported from other functions as described above). The first digit is the method used for particles from calc_mu() above and the second digit is the method used for antiparticles.

In pair_density(), the integer is a four-digit number. The first digit is from the list below and the remaining three digits, if nonzero, are from pair_mu().

  • 1: T=0 result

  • 2: default solver

  • 3: bracketing solver

  • 4: default solver with smaller tolerances

  • 5: default solver with smaller tolerances in log units

  • 6: bracketing in log units

std::string last_method_s

String detailing last method used.

inline fermion_rel_tl()

Create a fermion with mass m and degeneracy g.

inline virtual ~fermion_rel_tl()
inline virtual const char *type()

Return string denoting type (“fermion_rel”)

Template versions of base functions

inline int nu_from_n(fermion_t &f, fp_t temper)

Calculate the chemical potential from the density (template version)

inline virtual int calc_mu(fermion_t &f, fp_t temper)

Calculate properties as function of chemical potential.

inline virtual int calc_density(fermion_t &f, fp_t temper)

Calculate properties as function of density.

This function uses the current value of nu (or mu if the particle is non interacting) for an initial guess to solve for the chemical potential. If this guess is too small, then this function may fail.

\verbatim embed:rst

.. todo::

In function fermion_rel_tl::calc_density()

- Future: There is still quite a bit of code duplication
between this function and \ref calc_mu() .

\endverbatim

inline virtual void pair_mu(fermion_t &f, fp_t temper)

Calculate properties with antiparticles as function of chemical potential.

inline virtual int pair_density(fermion_t &f, fp_t temper)

Calculate thermodynamic properties with antiparticles from the density.

\verbatim embed:rst

.. todo::

In function pair_density():

- This actually works for negative densities some of the
time, but the solver probably doesn't work as well there and
we need to document the density expectations for this 
function.

\endverbatim
inline fp_t solve_fun(fp_t x, fermion_t &f, fp_t T)

Solve for the chemical potential given the density.

inline fp_t pair_fun(fp_t x, fp_t density_match, fermion_t &f, fp_t T, bool log_mode)

Solve for the chemical potential given the density with antiparticles.

\verbatim embed:rst

.. todo::

In function fermion_rel_tl::calc_density()

- Future: Particles and antiparticles have different
degeneracy factors, so we separately use the expansions one
at a time. It is probably better to separately generate a
new expansion function which automatically handles the sum
of particles and antiparticles.

\endverbatim