fermion_rel.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_REL_FERMION_H
24 #define O2SCL_REL_FERMION_H
25 
26 /** \file fermion_rel.h
27  \brief File defining \ref o2scl::fermion_rel
28 */
29 
30 #include <string>
31 #include <iostream>
32 #include <fstream>
33 #include <cmath>
34 #include <o2scl/constants.h>
35 #include <o2scl/mroot.h>
36 #include <o2scl/inte.h>
37 #include <o2scl/fermion.h>
38 
39 #ifndef DOXYGEN_NO_O2NS
40 namespace o2scl {
41 #endif
42 
43  /** \brief Equation of state for a relativistic fermion
44 
45  This class computes the thermodynamics of a relativistic fermion
46  either as a function of the density or the chemical potential.
47  It employs direct integration, using two different integrators
48  for the degenerate and non-degenerate regimes. The default
49  integrators are inte_qag_gsl (for degenerate fermions) and
50  inte_qagiu_gsl (for non-degenerate fermions). For the functions
51  calc_mu() and calc_density(), if the temperature argument is
52  less than or equal to zero, the functions \ref
53  fermion_zerot::calc_mu_zerot() and \ref
54  fermion_zerot::calc_density_zerot() will be used to compute the
55  result.
56 
57  \hline
58  <b>Degeneracy parameter:</b>
59 
60  Define the degeneracy parameter
61  \f[
62  \psi=(\nu-m^{*})/T
63  \f]
64  where \f$ \nu \f$ is the effective chemical potential (including
65  the rest mass) and \f$
66  m^{*} \f$ is the effective mass. For \f$ \psi \f$ smaller than
67  \ref min_psi, the non-degenerate expansion in \ref
68  fermion_eval_thermo::calc_mu_ndeg() is attempted first. If that
69  fails, then integration is used. For \f$ \psi \f$ greater than
70  \ref deg_limit (degenerate regime), a finite interval integrator
71  is used and for \f$ \psi \f$ less than \ref deg_limit
72  (non-degenerate regime), an integrator over the interval from
73  \f$ [0,\infty) \f$ is used. In the case where \ref
74  part::inc_rest_mass is false, the degeneracy parameter is
75  \f[
76  \psi=(\nu+m-m^{*})/T
77  \f]
78 
79  <b>Integration limits:</b>
80 
81  The upper limit on the degenerate integration is given by
82  \f[
83  \mathrm{upper~limit} = \sqrt{{\cal L}^2-m^{*,2}}
84  \f]
85  where \f$ {\cal L}\equiv u T+\nu \f$ and \f$ u \f$ is \ref
86  fermion_rel::upper_limit_fac . In the case where \ref
87  part::inc_rest_mass is false, the result is
88  \f[
89  \mathrm{upper~limit} = \sqrt{(m+{\cal L})^2-m^{*2}}
90  \f]
91 
92  The entropy is only significant at the Fermi surface, thus
93  in the degenerate case, the lower limit of the entropy
94  integral can be given be determined by the value of \f$ k \f$
95  which solves
96  \f[
97  - u = \frac{\sqrt{k^2+m^{* 2}}-\nu}{T}
98  \f]
99  The solution is
100  \f[
101  \mathrm{lower~limit} = \sqrt{(-u T+{\nu})^2-m^{*,2}}
102  \f]
103  but this solution is only valid if \f$ (m^{*}-\nu)/T < -u \f$.
104  In the case where part::inc_rest_mass is false, the result is
105  \f[
106  \mathrm{lower~limit} = \sqrt{(-u T + m +\nu)^2-m^{*,2}}
107  \f]
108  which is valid if \f$ (m^{*}-\nu - m)/T < -u \f$.
109 
110  <b>Entropy integrand:</b>
111 
112  In the degenerate regime, the entropy integrand
113  \f[
114  - k^2 \left[ f \log f + \left(1-f\right) \log
115  \left(1-f \right) \right]
116  \f]
117  where \f$ f \f$ is the fermionic distribution function can lose
118  precision when \f$ (E^{*} - \nu)/T \f$ is negative and
119  sufficiently large in absolute magnitude. Thus when \f$ (E^{*} -
120  \nu)/T < S \f$ where \f$ S \f$ is stored in \ref deg_entropy_fac
121  (default is -30), the integrand is written as
122  \f[
123  -k^2 \left( E/T-\nu/T \right) e^{E/T-\nu/T} \, .
124  \f]
125  If \f$ (E - \nu)/T < S \f$ is less than -1 times \ref exp_limit
126  (e.g. less than -200), then the entropy integrand is assumed
127  to be zero.
128 
129  <b>Non-degenerate integrands:</b>
130 
131  \comment
132  It's not at all clear that this dimensionless form is more
133  accurate than other potential alternatives. On the other hand,
134  it seems that the uncertainties in the integrations are larger
135  than the errors made by the integrand at present.
136  \endcomment
137  The integrands in the non-degenerate regime are written
138  in a dimensionless form, by defining \f$ u \f$ with
139  the relation
140  \f$ p = \sqrt{\left(T u + m^{*}\right)^2-m^{* 2}} \f$,
141  \f$ y \equiv \nu/ T \f$, and
142  \f$ \mathrm{mx} \equiv m^{*}/T \f$.
143  The density integrand is
144  \f[
145  \left(\mathrm{mx}+u\right) \sqrt{u^2+2 (\mathrm{mx}) u}
146  \left(\frac{e^{y}}{e^{\mathrm{mx}+u}+e^{y}}\right) \, ,
147  \f]
148  the energy integrand is
149  \f[
150  \left(\mathrm{mx}+u\right)^2 \sqrt{u^2+2 (\mathrm{mx}) u}
151  \left(\frac{e^{y}}{e^{\mathrm{mx}+u}+e^{y}}\right) \, ,
152  \f]
153  and the entropy integrand is
154  \f[
155  \left(\mathrm{mx}+u\right) \sqrt{u^2+2 (\mathrm{mx}) u}
156  \left(t_1+t_2\right) \, ,
157  \f]
158  where
159  \f{eqnarray*}
160  t_1 &=& \log \left(1+e^{y-\mathrm{mx}-u}\right)/
161  \left(1+e^{y-\mathrm{mx}-u}\right) \nonumber \\
162  t_2 &=& \log \left(1+e^{\mathrm{mx}+u-y}\right)/
163  \left(1+e^{\mathrm{mx}+u-y}\right) \, .
164  \f}
165 
166  \hline
167  <b>Accuracy:</b>
168 
169  The default settings for for this class give an accuracy of at
170  least 1 part in \f$ 10^6 \f$ (and frequently better than this).
171 
172  When the integrators provide numerical uncertainties, these
173  uncertainties are stored in \ref unc. In the case of
174  calc_density() and pair_density(), the uncertainty from the
175  numerical accuracy of the solver is not included. (There is also
176  a relatively small inaccuracy due to the mathematical evaluation
177  of the integrands which is not included in \ref unc.)
178 
179  One can improve the accuracy to within 1 part in \f$ 10^{10} \f$
180  using
181  \code
182  fermion_rel rf(1.0,2.0);
183  rf.upper_limit_fac=40.0;
184  rf.dit->tol_abs=1.0e-13;
185  rf.dit->tol_rel=1.0e-13;
186  rf.nit->tol_abs=1.0e-13;
187  rf.nit->tol_rel=1.0e-13;
188  rf.density_root->tol_rel=1.0e-10;
189  \endcode
190  which decreases the both the relative and absolute tolerances
191  for both the degenerate and non-degenerate integrators and
192  improves the accuracy of the solver which determines the
193  chemical potential from the density. Of course if these
194  tolerances are too small, the calculation may fail.
195 
196  \hline
197  <b>Todos:</b>
198 
199  \future The expressions which appear in in the integrand
200  functions density_fun(), etc. could likely be improved,
201  especially in the case where \ref o2scl::part::inc_rest_mass is
202  <tt>false</tt>. There should not be a need to check if
203  <tt>ret</tt> is finite.
204 
205  \future It appears this class doesn't compute the uncertainty in
206  the chemical potential or density with calc_density(). This
207  could be fixed.
208 
209  \future I'd like to change the lower limit on the entropy
210  integration, but the value in the code at the moment (stored
211  in <tt>ll</tt>) makes bm_part2.cpp worse.
212 
213  \future The function pair_mu() should set the antiparticle
214  integrators as done in fermion_deriv_rel.
215  */
217 
218  public:
219 
220  /// \name Numerical parameters
221  //@{
222  /** \brief If true, call the error handler when convergence
223  fails (default true)
224  */
226 
227  /** \brief The smallest value of \f$ (\mu-m)/T \f$ for which
228  integration is used
229  */
230  double min_psi;
231 
232  /** \brief The critical degeneracy at which to switch integration
233  techniques (default 2)
234  */
235  double deg_limit;
236 
237  /** \brief The limit for exponentials to ensure integrals are finite
238  (default 200)
239  */
240  double exp_limit;
241 
242  /// The factor for the degenerate upper limits (default 20)
244 
245  /// A factor for the degenerate entropy integration (default 30)
247  //@}
248 
249  /// Storage for the uncertainty
251 
252  /// If true, use expansions for extreme conditions (default true)
254 
255  /// Create a fermion with mass \c m and degeneracy \c g
256  fermion_rel();
257 
258  virtual ~fermion_rel();
259 
260  /** \brief Calculate properties as function of chemical potential
261  */
262  virtual void calc_mu(fermion &f, double temper);
263 
264  /** \brief Calculate properties as function of density
265 
266  This function uses the current value of \c nu (or \c mu if the
267  particle is non interacting) for an initial guess to solve for
268  the chemical potential. If this guess is too small, then this
269  function may fail.
270  */
271  virtual int calc_density(fermion &f, double temper);
272 
273  /** \brief Calculate properties with antiparticles as function of
274  chemical potential
275  */
276  virtual void pair_mu(fermion &f, double temper);
277 
278  /** \brief Calculate properties with antiparticles as function of
279  density
280  */
281  virtual int pair_density(fermion &f, double temper);
282 
283  /** \brief Calculate effective chemical potential from density
284 
285  \future This function might be improved by generating a
286  bracket for a bracketing solver, rather than \ref
287  o2scl::root_cern which is the default for \ref
288  o2scl::fermion_rel::density_root.
289  */
290  virtual int nu_from_n(fermion &f, double temper);
291 
292  /// The non-degenerate integrator
293  std::shared_ptr<inte<> > nit;
294 
295  /// The degenerate integrator
296  std::shared_ptr<inte<> > dit;
297 
298  /// The solver for calc_density()
299  std::shared_ptr<root<> > density_root;
300 
301  /// Return string denoting type ("fermion_rel")
302  virtual const char *type() { return "fermion_rel"; }
303 
304  protected:
305 
306 #ifndef DOXYGEN_INTERNAL
307 
308  /// The integrand for the density for non-degenerate fermions
309  double density_fun(double u, fermion &f, double T);
310 
311  /// The integrand for the energy density for non-degenerate fermions
312  double energy_fun(double u, fermion &f, double T);
313 
314  /// The integrand for the entropy density for non-degenerate fermions
315  double entropy_fun(double u, fermion &f, double T);
316 
317  /// The integrand for the density for degenerate fermions
318  double deg_density_fun(double u, fermion &f, double T);
319 
320  /// The integrand for the energy density for degenerate fermions
321  double deg_energy_fun(double u, fermion &f, double T);
322 
323  /// The integrand for the entropy density for degenerate fermions
324  double deg_entropy_fun(double u, fermion &f, double T);
325 
326  /// Solve for the chemical potential given the density
327  double solve_fun(double x, fermion &f, double T);
328 
329  /** \brief Solve for the chemical potential given the density
330  with antiparticles
331 
332  \future Particles and antiparticles have different degeneracy
333  factors, so we separately use the expansions one at a time. It
334  is probably better to separately generate a new expansion
335  function which automatically handles the sum of particles and
336  antiparticles.
337  */
338  double pair_fun(double x, fermion &f, double T, bool log_mode);
339 
340 #endif
341 
342  };
343 
344 #ifndef DOXYGEN_NO_O2NS
345 }
346 #endif
347 
348 #endif
double upper_limit_fac
The factor for the degenerate upper limits (default 20)
Definition: fermion_rel.h:243
bool err_nonconv
If true, call the error handler when convergence fails (default true)
Definition: fermion_rel.h:225
double deg_entropy_fac
A factor for the degenerate entropy integration (default 30)
Definition: fermion_rel.h:246
double exp_limit
The limit for exponentials to ensure integrals are finite (default 200)
Definition: fermion_rel.h:240
bool use_expansions
If true, use expansions for extreme conditions (default true)
Definition: fermion_rel.h:253
double min_psi
The smallest value of for which integration is used.
Definition: fermion_rel.h:230
Fermion class.
Definition: fermion.h:52
std::shared_ptr< inte<> > nit
The non-degenerate integrator.
Definition: fermion_rel.h:293
virtual void pair_mu(fermion &f, double temper)
Calculate properties with antiparticles as function of chemical potential.
Fermion with finite-temperature thermodynamics [abstract base].
Definition: fermion.h:199
double pair_fun(double x, fermion &f, double T, bool log_mode)
Solve for the chemical potential given the density with antiparticles.
virtual const char * type()
Return string denoting type ("fermion_rel")
Definition: fermion_rel.h:302
std::shared_ptr< inte<> > dit
The degenerate integrator.
Definition: fermion_rel.h:296
double entropy_fun(double u, fermion &f, double T)
The integrand for the entropy density for non-degenerate fermions.
std::shared_ptr< root<> > density_root
The solver for calc_density()
Definition: fermion_rel.h:299
fermion unc
Storage for the uncertainty.
Definition: fermion_rel.h:250
double deg_limit
The critical degeneracy at which to switch integration techniques (default 2)
Definition: fermion_rel.h:235
Equation of state for a relativistic fermion.
Definition: fermion_rel.h:216
virtual int nu_from_n(fermion &f, double temper)
Calculate effective chemical potential from density.
double deg_density_fun(double u, fermion &f, double T)
The integrand for the density for degenerate fermions.
double solve_fun(double x, fermion &f, double T)
Solve for the chemical potential given the density.
double deg_entropy_fun(double u, fermion &f, double T)
The integrand for the entropy density for degenerate fermions.
double energy_fun(double u, fermion &f, double T)
The integrand for the energy density for non-degenerate fermions.
virtual int calc_density(fermion &f, double temper)
Calculate properties as function of density.
fermion_rel()
Create a fermion with mass m and degeneracy g.
virtual void calc_mu(fermion &f, double temper)
Calculate properties as function of chemical potential.
double deg_energy_fun(double u, fermion &f, double T)
The integrand for the energy density for degenerate fermions.
virtual int pair_density(fermion &f, double temper)
Calculate properties with antiparticles as function of density.
double density_fun(double u, fermion &f, double T)
The integrand for the density for non-degenerate fermions.

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