fermion_deriv_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_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
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 
42 #ifndef DOXYGEN_NO_O2NS
43 namespace o2scl {
44 #endif
45 
46  /** \brief Equation of state for a relativistic fermion
47 
48  \note This class only has preliminary support for
49  inc_rest_mass=true (more testing should be done, particularly
50  for the "pair" functions)
51 
52  This implements an equation of state for a relativistic fermion
53  using direct integration. After subtracting the rest mass from
54  the chemical potentials, the distribution function is
55  \f[
56  \left\{1+\exp[(\sqrt{k^2+m^{* 2}}-m-\nu)/T]\right\}^{-1}
57  \f]
58  where \f$ k \f$ is the momentum, \f$ \nu \f$ is the effective
59  chemical potential, \f$ m \f$ is the rest mass, and \f$ m^{*}
60  \f$ is the effective mass. For later use, we define \f$ E^{*} =
61  \sqrt{k^2 + m^{*2}} \f$ . The degeneracy parameter is
62  \f[
63  \psi=(\nu+(m-m^{*}))/T
64  \f]
65  For \f$ \psi \f$ greater than \ref deg_limit (degenerate
66  regime), a finite interval integrator is used and for \f$ \psi
67  \f$ less than \ref deg_limit (non-degenerate regime), an
68  integrator over the interval from \f$ [0,\infty) \f$ is
69  used. The upper limit on the degenerate integration is given by
70  the value of the momentum \f$ k \f$ which is the solution of
71  \f[
72  (\sqrt{k^2+m^{*,2}}-m-\nu)/T=\mathrm{f{l}imit}
73  \f]
74  which is
75  \f[
76  \sqrt{(m+{\cal L})^2-m^{*2}}
77  \f]
78  where \f$ {\cal L}\equiv\mathrm{f{l}imit}\times T+\nu \f$ .
79 
80  For the entropy integration, we set the lower limit
81  to
82  \f[
83  2 \sqrt{\nu^2+2 \nu m} - \mathrm{upper~limit}
84  \f]
85  since the only contribution to the entropy is at the Fermi surface.
86  \comment
87  I'm not sure, but I think this is an expression determined
88  from a small T taylor expansion of the argument of the
89  exponential.
90  \endcomment
91 
92  In the non-degenerate regime, we make the substitution \f$ u=k/T
93  \f$ to help ensure that the variable of integration scales
94  properly.
95 
96  Uncertainties are given in \ref unc.
97 
98  \todo This needs to be corrected to calculate \f$ \sqrt{k^2+m^{*
99  2}}-m \f$ gracefully when \f$ m^{*}\approx m \ll k \f$ .
100  \todo Call error handler if inc_rest_mass is true or update
101  to properly treat the case when inc_rest_mass is true.
102 
103  \b Evaluation \b of \b the \b derivatives
104 
105  The relevant
106  derivatives of the distribution function are
107  \f[
108  \frac{\partial f}{\partial T}=
109  f(1-f)\frac{E^{*}-m-\nu}{T^2}
110  \f]
111  \f[
112  \frac{\partial f}{\partial \nu}=
113  f(1-f)\frac{1}{T}
114  \f]
115  \f[
116  \frac{\partial f}{\partial k}=
117  -f(1-f)\frac{k}{E^{*} T}
118  \f]
119  \f[
120  \frac{\partial f}{\partial m^{*}}=
121  -f(1-f)\frac{m^{*}}{E^{*} T}
122  \f]
123 
124  We also need the derivative of the entropy integrand w.r.t. the
125  distribution function, which is
126  \f[
127  {\cal S}\equiv f \ln f +(1-f) \ln (1-f) \qquad
128  \frac{\partial {\cal S}}{\partial f} = \ln
129  \left(\frac{f}{1-f}\right) =
130  \left(\frac{\nu-E^{*}+m}{T}\right)
131  \f]
132  where the entropy density is
133  \f[
134  s = - \frac{g}{2 \pi^2} \int_0^{\infty} {\cal S} k^2 d k
135  \f]
136 
137  The derivatives can be integrated directly (\ref method = \ref
138  direct) or they may be converted to integrals over the
139  distribution function through an integration by parts (\ref
140  method = \ref by_parts)
141  \f[
142  \int_a^b f(k) \frac{d g(k)}{dk} dk = \left.f(k) g(k)\right|_{k=a}^{k=b}
143  - \int_a^b g(k) \frac{d f(k)}{dk} dk
144  \f]
145  using the distribution function for \f$ f(k) \f$ and 0 and
146  \f$ \infty \f$ as the limits, we have
147  \f[
148  \frac{g}{2 \pi^2} \int_0^{\infty} \frac{d g(k)}{dk} f dk =
149  \frac{g}{2 \pi^2} \int_0^{\infty} g(k) f (1-f) \frac{k}{E^{*} T} dk
150  \f]
151  as long as \f$ g(k) \f$ vanishes at \f$ k=0 \f$ .
152  Rewriting,
153  \f[
154  \frac{g}{2 \pi^2} \int_0^{\infty} h(k) f (1-f) dk =
155  \frac{g}{2 \pi^2} \int_0^{\infty} f \frac{T}{k}
156  \left[ h^{\prime} E^{*}-\frac{h E^{*}}{k}+\frac{h k}{E^{*}} \right] dk
157  \f]
158  as long as \f$ h(k)/k \f$ vanishes at \f$ k=0 \f$ .
159 
160  \b Explicit \b forms
161 
162  1) The derivative of the density wrt the chemical potential
163  \f[
164  \left(\frac{d n}{d \mu}\right)_T =
165  \frac{g}{2 \pi^2} \int_0^{\infty} \frac{k^2}{T} f (1-f) dk
166  \f]
167  Using \f$ h(k)=k^2/T \f$ we get
168  \f[
169  \left(\frac{d n}{d \mu}\right)_T =
170  \frac{g}{2 \pi^2} \int_0^{\infty}
171  \left(\frac{k^2+E^{*2}}{E^{*}}\right) f dk
172  \f]
173 
174  2) The derivative of the density wrt the temperature
175  \f[
176  \left(\frac{d n}{d T}\right)_{\mu} =
177  \frac{g}{2 \pi^2} \int_0^{\infty} \frac{k^2(E^{*}-m-\nu)}{T^2}
178  f (1-f) dk
179  \f]
180  Using \f$ h(k)=k^2(E^{*}-\nu)/T^2 \f$ we get
181  \f[
182  \left(\frac{d n}{d T}\right)_{\mu} =
183  \frac{g}{2 \pi^2} \int_0^{\infty} \frac{f}{T}
184  \left[2 k^2+E^{*2}-E^{*}\left(\nu+m\right)-
185  k^2 \left(\frac{\nu+m}{E^{*}}\right)\right] dk
186  \f]
187 
188  3) The derivative of the entropy wrt the chemical potential
189  \f[
190  \left(\frac{d s}{d \mu}\right)_T =
191  \frac{g}{2 \pi^2} \int_0^{\infty} k^2 f (1-f)
192  \frac{(E^{*}-m-\nu)}{T^2} dk
193  \f]
194  This verifies the Maxwell relation
195  \f[
196  \left(\frac{d s}{d \mu}\right)_T =
197  \left(\frac{d n}{d T}\right)_{\mu}
198  \f]
199 
200  4) The derivative of the entropy wrt the temperature
201  \f[
202  \left(\frac{d s}{d T}\right)_{\mu} =
203  \frac{g}{2 \pi^2} \int_0^{\infty} k^2 f (1-f)
204  \frac{(E^{*}-m-\nu)^2}{T^3} dk
205  \f]
206  Using \f$ h(k)=k^2 (E^{*}-\nu)^2/T^3 \f$
207  \f[
208  \left(\frac{d s}{d T}\right)_{\mu} =
209  \frac{g}{2 \pi^2} \int_0^{\infty} \frac{f(E^{*}-m-\nu)}{E^{*}T^2}
210  \left[E^{* 3}+3 E^{*} k^2- (E^{* 2}+k^2)(\nu+m)\right] d k
211  \f]
212 
213  5) The derivative of the density wrt the effective mass
214  \f[
215  \left(\frac{d n}{d m^{*}}\right)_{T,\mu} =
216  -\frac{g}{2 \pi^2} \int_0^{\infty}
217  \frac{k^2 m^{*}}{E^{*} T} f (1-f) dk
218  \f]
219  Using \f$ h(k)=-(k^2 m^{*})/(E^{*} T) \f$ we get
220  \f[
221  \left(\frac{d n}{d m^{*}}\right)_{T,\mu} =
222  -\frac{g}{2 \pi^2} \int_0^{\infty}
223  m^{*} f dk
224  \f]
225  \comment
226  This derivative may be written in terms of the
227  others
228  \f[
229  \left(\frac{d n}{d m^{*}}\right)_{T,\mu} = \frac{3 n}{m^{*}}
230  - \frac{T}{m^{*}}\left[ \left(\frac{d n}{d T}\right)_{\mu}
231  +\frac{\mu}{T} \left(\frac{d n}{d \mu}\right)_{T}
232  \right] - \left(\frac{d n}{d \mu}\right)_{T}
233  \f]
234  \endcomment
235 
236  \note The dsdT integration may fail if the system is
237  very degenerate. When method is byparts, the integral involves a
238  large cancellation between the regions from \f$ k \in (0,
239  \mathrm{ulimit/2}) \f$ and \f$ k \in (\mathrm{ulimit/2},
240  \mathrm{ulimit}) \f$. Switching to method=direct and setting the
241  lower limit to \f$ \mathrm{llimit} \f$, may help, but recent
242  testing on this gave negative values for dsdT. For very
243  degenerate systems, an expansion may be better than trying
244  to perform the integration. The value of the integrand
245  at k=0 also looks like it might be causing difficulties.
246 
247  \todo err_nonconv=false not implemented yet.
248 
249  \future It might be worth coding up direct differentiation, or
250  differentiating the eff results, as these may succeed more
251  generally.
252 
253  \future This class will have difficulty with extremely degenerate
254  or extremely non-degnerate systems. Fix this.
255 
256  \future Create a more intelligent method for dealing with bad
257  initial guesses for the chemical potential in calc_density().
258  */
260 
261  public:
262 
263  /// Create a fermion with mass \c m and degeneracy \c g
265 
266  virtual ~fermion_deriv_rel();
267 
268  /** \brief Limit of arguments of exponentials for Fermi functions
269  (default 200.0)
270  */
271  double exp_limit;
272 
273  /** \brief The critical degeneracy at which to switch integration
274  techniques (default 2.0)
275  */
276  double deg_limit;
277 
278  /** \brief The limit for the Fermi functions (default 20.0)
279 
280  fermion_deriv_rel will ignore corrections smaller than about
281  \f$ \exp(-\mathrm{f{l}imit}) \f$ .
282  */
284 
285  /// Storage for the most recently calculated uncertainties
287 
288  /** \name Method of computing derivatives
289  */
290  //@{
291  /// Method (default is \ref automatic)
292  int method;
293  /// Automatically choose method
294  static const int automatic=0;
295  /// In the form containing \f$ f(1-f) \f$ .
296  static const int direct=1;
297  /// Integrate by parts
298  static const int by_parts=2;
299  //@}
300 
301  /** \brief If true, call the error handler when convergence
302  fails (default true)
303  */
305 
306  /** \brief Calculate properties as function of chemical potential
307  */
308  virtual int calc_mu(fermion_deriv &f, double temper);
309 
310  /** \brief Calculate properties as function of density
311  */
312  virtual int calc_density(fermion_deriv &f, double temper);
313 
314  /** \brief Calculate properties with antiparticles as function of
315  chemical potential
316  */
317  virtual int pair_mu(fermion_deriv &f, double temper);
318 
319  /** \brief Calculate properties with antiparticles as function of
320  density
321  */
322  virtual int pair_density(fermion_deriv &f, double temper);
323 
324  /// Calculate effective chemical potential from density
325  virtual int nu_from_n(fermion_deriv &f, double temper);
326 
327  /** \brief Set inte objects
328 
329  The first integrator is used for non-degenerate integration
330  and should integrate from 0 to \f$ \infty \f$ (like \ref
331  o2scl::inte_qagiu_gsl). The second integrator is for the
332  degenerate case, and should integrate between two finite
333  values.
334  */
335  void set_inte(inte<> &unit, inte<> &udit);
336 
337  /** \brief Set the solver for use in calculating the chemical
338  potential from the density */
340  density_root=&rp;
341  return;
342  }
343 
344  /// The default integrator for the non-degenerate regime
346 
347  /// The default integrator for the degenerate regime
349 
350  /// The default solver for npen_density() and pair_density()
352 
353  /// Return string denoting type ("fermion_deriv_rel")
354  virtual const char *type() { return "fermion_deriv_rel"; };
355 
356  /// Calibrate with more accurate tabulated results
357  double deriv_calibrate(fermion_deriv &f, int verbose,
358  std::string fname="");
359 
360  protected:
361 
362 #ifndef DOXYGEN_NO_O2NS
363 
364  /// The internal integration method
366 
367  /// The integrator for non-degenerate fermions
369 
370  /// The integrator for degenerate fermions
372 
373  /// The solver for calc_density() and pair_density()
375 
376  /** \name The integrands, as a function of \f$ u=k/T \f$, for
377  non-degenerate integrals
378  */
379  //@{
380  double density_fun(double u, fermion_deriv &f, double T);
381  double energy_fun(double u, fermion_deriv &f, double T);
382  double entropy_fun(double u, fermion_deriv &f, double T);
383  double density_T_fun(double k, fermion_deriv &f, double T);
384  double density_mu_fun(double k, fermion_deriv &f, double T);
385  double entropy_T_fun(double k, fermion_deriv &f, double T);
386  double density_ms_fun(double k, fermion_deriv &f, double T);
387  //@}
388 
389  /** \name The integrands, as a function of momentum, for the
390  degenerate integrals
391  */
392  //@{
393  double deg_density_fun(double u, fermion_deriv &f, double T);
394  double deg_energy_fun(double u, fermion_deriv &f, double T);
395  double deg_entropy_fun(double u, fermion_deriv &f, double T);
396  double deg_density_T_fun(double k, fermion_deriv &f, double T);
397  double deg_density_mu_fun(double k, fermion_deriv &f, double T);
398  double deg_entropy_T_fun(double k, fermion_deriv &f, double T);
399  double deg_density_ms_fun(double k, fermion_deriv &f, double T);
400  //@}
401 
402  /** \brief Solve for the chemical potential from the density
403  for calc_density()
404  */
405  double solve_fun(double x, fermion_deriv &f, double T);
406 
407  /** \brief Solve for the chemical potential from the density
408  for pair_density()
409  */
410  double pair_fun(double x, fermion_deriv &f, double T);
411 
412 #endif
413 
414  };
415 
416 #ifndef DOXYGEN_NO_O2NS
417 }
418 #endif
419 
420 #endif
inte * nit
The integrator for non-degenerate fermions.
void set_inte(inte<> &unit, inte<> &udit)
Set inte objects.
fermion_deriv unc
Storage for the most recently calculated uncertainties.
root * density_root
The solver for calc_density() and pair_density()
inte * dit
The integrator for degenerate fermions.
bool err_nonconv
If true, call the error handler when convergence fails (default true)
static const int by_parts
Integrate by parts.
double upper_limit_fac
The limit for the Fermi functions (default 20.0)
virtual int pair_mu(fermion_deriv &f, double temper)
Calculate properties with antiparticles as function of chemical potential.
A fermion with derivative information.
Definition: part_deriv.h:345
Compute properties of a fermion including derivatives [abstract base].
Definition: part_deriv.h:410
root_cern def_density_root
The default solver for npen_density() and pair_density()
inte_qag_gsl def_dit
The default integrator for the degenerate regime.
int intl_method
The internal integration method.
inte_qagiu_gsl def_nit
The default integrator for the non-degenerate regime.
fermion_deriv_rel()
Create a fermion with mass m and degeneracy g.
double deg_limit
The critical degeneracy at which to switch integration techniques (default 2.0)
static const int automatic
Automatically choose method.
virtual int nu_from_n(fermion_deriv &f, double temper)
Calculate effective chemical potential from density.
static const int direct
In the form containing .
virtual int calc_mu(fermion_deriv &f, double temper)
Calculate properties as function of chemical potential.
int method
Method (default is automatic)
double solve_fun(double x, fermion_deriv &f, double T)
Solve for the chemical potential from the density for calc_density()
void set_density_root(root<> &rp)
Set the solver for use in calculating the chemical potential from the density.
virtual const char * type()
Return string denoting type ("fermion_deriv_rel")
double exp_limit
Limit of arguments of exponentials for Fermi functions (default 200.0)
Equation of state for a relativistic fermion.
virtual int calc_density(fermion_deriv &f, double temper)
Calculate properties as function of density.
virtual int pair_density(fermion_deriv &f, double temper)
Calculate properties with antiparticles as function of density.
double deriv_calibrate(fermion_deriv &f, int verbose, std::string fname="")
Calibrate with more accurate tabulated results.
double pair_fun(double x, fermion_deriv &f, double T)
Solve for the chemical potential from the density for pair_density()

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