fermion_deriv_rel.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2018, 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  \future This class may need more corrections to ensure
99  quantities like \f$ \sqrt{k^2+m^{*2}}-m \f$ are computed
100  accurately when \f$ m^{*}\approx m \ll k \f$ .
101 
102  \todo Call error handler if inc_rest_mass is true or update
103  to properly treat the case when inc_rest_mass is true.
104 
105  \b Evaluation \b of \b the \b derivatives
106 
107  The relevant
108  derivatives of the distribution function are
109  \f[
110  \frac{\partial f}{\partial T}=
111  f(1-f)\frac{E^{*}-m-\nu}{T^2}
112  \f]
113  \f[
114  \frac{\partial f}{\partial \nu}=
115  f(1-f)\frac{1}{T}
116  \f]
117  \f[
118  \frac{\partial f}{\partial k}=
119  -f(1-f)\frac{k}{E^{*} T}
120  \f]
121  \f[
122  \frac{\partial f}{\partial m^{*}}=
123  -f(1-f)\frac{m^{*}}{E^{*} T}
124  \f]
125 
126  We also need the derivative of the entropy integrand w.r.t. the
127  distribution function, which is
128  \f[
129  {\cal S}\equiv f \ln f +(1-f) \ln (1-f) \qquad
130  \frac{\partial {\cal S}}{\partial f} = \ln
131  \left(\frac{f}{1-f}\right) =
132  \left(\frac{\nu-E^{*}+m}{T}\right)
133  \f]
134  where the entropy density is
135  \f[
136  s = - \frac{g}{2 \pi^2} \int_0^{\infty} {\cal S} k^2 d k
137  \f]
138 
139  The derivatives can be integrated directly (\ref method = \ref
140  direct) or they may be converted to integrals over the
141  distribution function through an integration by parts (\ref
142  method = \ref by_parts)
143  \f[
144  \int_a^b f(k) \frac{d g(k)}{dk} dk = \left.f(k) g(k)\right|_{k=a}^{k=b}
145  - \int_a^b g(k) \frac{d f(k)}{dk} dk
146  \f]
147  using the distribution function for \f$ f(k) \f$ and 0 and
148  \f$ \infty \f$ as the limits, we have
149  \f[
150  \frac{g}{2 \pi^2} \int_0^{\infty} \frac{d g(k)}{dk} f dk =
151  \frac{g}{2 \pi^2} \int_0^{\infty} g(k) f (1-f) \frac{k}{E^{*} T} dk
152  \f]
153  as long as \f$ g(k) \f$ vanishes at \f$ k=0 \f$ .
154  Rewriting,
155  \f[
156  \frac{g}{2 \pi^2} \int_0^{\infty} h(k) f (1-f) dk =
157  \frac{g}{2 \pi^2} \int_0^{\infty} f \frac{T}{k}
158  \left[ h^{\prime} E^{*}-\frac{h E^{*}}{k}+\frac{h k}{E^{*}} \right] dk
159  \f]
160  as long as \f$ h(k)/k \f$ vanishes at \f$ k=0 \f$ .
161 
162  \b Explicit \b forms
163 
164  1) The derivative of the density wrt the chemical potential
165  \f[
166  \left(\frac{d n}{d \mu}\right)_T =
167  \frac{g}{2 \pi^2} \int_0^{\infty} \frac{k^2}{T} f (1-f) dk
168  \f]
169  Using \f$ h(k)=k^2/T \f$ we get
170  \f[
171  \left(\frac{d n}{d \mu}\right)_T =
172  \frac{g}{2 \pi^2} \int_0^{\infty}
173  \left(\frac{k^2+E^{*2}}{E^{*}}\right) f dk
174  \f]
175 
176  2) The derivative of the density wrt the temperature
177  \f[
178  \left(\frac{d n}{d T}\right)_{\mu} =
179  \frac{g}{2 \pi^2} \int_0^{\infty} \frac{k^2(E^{*}-m-\nu)}{T^2}
180  f (1-f) dk
181  \f]
182  Using \f$ h(k)=k^2(E^{*}-\nu)/T^2 \f$ we get
183  \f[
184  \left(\frac{d n}{d T}\right)_{\mu} =
185  \frac{g}{2 \pi^2} \int_0^{\infty} \frac{f}{T}
186  \left[2 k^2+E^{*2}-E^{*}\left(\nu+m\right)-
187  k^2 \left(\frac{\nu+m}{E^{*}}\right)\right] dk
188  \f]
189 
190  3) The derivative of the entropy wrt the chemical potential
191  \f[
192  \left(\frac{d s}{d \mu}\right)_T =
193  \frac{g}{2 \pi^2} \int_0^{\infty} k^2 f (1-f)
194  \frac{(E^{*}-m-\nu)}{T^2} dk
195  \f]
196  This verifies the Maxwell relation
197  \f[
198  \left(\frac{d s}{d \mu}\right)_T =
199  \left(\frac{d n}{d T}\right)_{\mu}
200  \f]
201 
202  4) The derivative of the entropy wrt the temperature
203  \f[
204  \left(\frac{d s}{d T}\right)_{\mu} =
205  \frac{g}{2 \pi^2} \int_0^{\infty} k^2 f (1-f)
206  \frac{(E^{*}-m-\nu)^2}{T^3} dk
207  \f]
208  Using \f$ h(k)=k^2 (E^{*}-\nu)^2/T^3 \f$
209  \f[
210  \left(\frac{d s}{d T}\right)_{\mu} =
211  \frac{g}{2 \pi^2} \int_0^{\infty} \frac{f(E^{*}-m-\nu)}{E^{*}T^2}
212  \left[E^{* 3}+3 E^{*} k^2- (E^{* 2}+k^2)(\nu+m)\right] d k
213  \f]
214 
215  5) The derivative of the density wrt the effective mass
216  \f[
217  \left(\frac{d n}{d m^{*}}\right)_{T,\mu} =
218  -\frac{g}{2 \pi^2} \int_0^{\infty}
219  \frac{k^2 m^{*}}{E^{*} T} f (1-f) dk
220  \f]
221  Using \f$ h(k)=-(k^2 m^{*})/(E^{*} T) \f$ we get
222  \f[
223  \left(\frac{d n}{d m^{*}}\right)_{T,\mu} =
224  -\frac{g}{2 \pi^2} \int_0^{\infty}
225  m^{*} f dk
226  \f]
227  \comment
228  This derivative may be written in terms of the
229  others
230  \f[
231  \left(\frac{d n}{d m^{*}}\right)_{T,\mu} = \frac{3 n}{m^{*}}
232  - \frac{T}{m^{*}}\left[ \left(\frac{d n}{d T}\right)_{\mu}
233  +\frac{\mu}{T} \left(\frac{d n}{d \mu}\right)_{T}
234  \right] - \left(\frac{d n}{d \mu}\right)_{T}
235  \f]
236  \endcomment
237 
238  \note The dsdT integration may fail if the system is
239  very degenerate. When method is byparts, the integral involves a
240  large cancellation between the regions from \f$ k \in (0,
241  \mathrm{ulimit/2}) \f$ and \f$ k \in (\mathrm{ulimit/2},
242  \mathrm{ulimit}) \f$. Switching to method=direct and setting the
243  lower limit to \f$ \mathrm{llimit} \f$, may help, but recent
244  testing on this gave negative values for dsdT. For very
245  degenerate systems, an expansion may be better than trying
246  to perform the integration. The value of the integrand
247  at k=0 also looks like it might be causing difficulties.
248 
249  \todo err_nonconv=false not implemented yet.
250 
251  \future It might be worth coding up direct differentiation, or
252  differentiating the eff results, as these may succeed more
253  generally.
254 
255  \future This class will have difficulty with extremely degenerate
256  or extremely non-degnerate systems. Fix this.
257 
258  \future Create a more intelligent method for dealing with bad
259  initial guesses for the chemical potential in calc_density().
260  */
262 
263  public:
264 
265  /// Create a fermion with mass \c m and degeneracy \c g
267 
268  virtual ~fermion_deriv_rel();
269 
270  /** \brief Limit of arguments of exponentials for Fermi functions
271  (default 200.0)
272  */
273  double exp_limit;
274 
275  /** \brief The critical degeneracy at which to switch integration
276  techniques (default 2.0)
277  */
278  double deg_limit;
279 
280  /** \brief The limit for the Fermi functions (default 20.0)
281 
282  fermion_deriv_rel will ignore corrections smaller than about
283  \f$ \exp(-\mathrm{f{l}imit}) \f$ .
284  */
286 
287  /// Storage for the most recently calculated uncertainties
289 
290  /** \name Method of computing derivatives
291  */
292  //@{
293  /// Method (default is \ref automatic)
294  int method;
295  /// Automatically choose method
296  static const int automatic=0;
297  /// In the form containing \f$ f(1-f) \f$ .
298  static const int direct=1;
299  /// Integrate by parts
300  static const int by_parts=2;
301  //@}
302 
303  /** \brief If true, call the error handler when convergence
304  fails (default true)
305  */
307 
308  /** \brief Calculate properties as function of chemical potential
309  */
310  virtual int calc_mu(fermion_deriv &f, double temper);
311 
312  /** \brief Calculate properties as function of density
313  */
314  virtual int calc_density(fermion_deriv &f, double temper);
315 
316  /** \brief Calculate properties with antiparticles as function of
317  chemical potential
318  */
319  virtual int pair_mu(fermion_deriv &f, double temper);
320 
321  /** \brief Calculate properties with antiparticles as function of
322  density
323  */
324  virtual int pair_density(fermion_deriv &f, double temper);
325 
326  /// Calculate effective chemical potential from density
327  virtual int nu_from_n(fermion_deriv &f, double temper);
328 
329  /** \brief Set inte objects
330 
331  The first integrator is used for non-degenerate integration
332  and should integrate from 0 to \f$ \infty \f$ (like \ref
333  o2scl::inte_qagiu_gsl). The second integrator is for the
334  degenerate case, and should integrate between two finite
335  values.
336  */
337  void set_inte(inte<> &unit, inte<> &udit);
338 
339  /** \brief Set the solver for use in calculating the chemical
340  potential from the density */
342  density_root=&rp;
343  return;
344  }
345 
346  /// The default integrator for the non-degenerate regime
348 
349  /// The default integrator for the degenerate regime
351 
352  /// The default solver for npen_density() and pair_density()
354 
355  /// Return string denoting type ("fermion_deriv_rel")
356  virtual const char *type() { return "fermion_deriv_rel"; };
357 
358  /** \brief Calibrate with more accurate tabulated results
359 
360  This compares the approximation to the exact results over a
361  grid with \f$ T = \left\{10^{-2},1,10^{2}\right\} \f$, \f$
362  \log_{10} (m/T) = \left\{ -3,-2,-1,0,1,2,3\right\} \f$, and
363  \f$ \log_{10} \psi = \left\{ -3,-2,-1,0,1\right\} \f$, where
364  \f$ \psi \equiv \left(\mu-m\right)/T \f$ using
365  calc_density() and calc_mu(), with both inc_rest_mass
366  taking both values <tt>true</tt> and <tt>false</tt>.
367 
368  The <tt>verbose</tt> parameter controls the amount of
369  output, and \c fname is the filename for the file
370  <tt>fermion_cal.o2</tt>.
371  */
372  double deriv_calibrate(fermion_deriv &f, int verbose,
373  std::string fname="");
374 
375  protected:
376 
377 #ifndef DOXYGEN_NO_O2NS
378 
379  /// The internal integration method
381 
382  /// The integrator for non-degenerate fermions
384 
385  /// The integrator for degenerate fermions
387 
388  /// The solver for calc_density() and pair_density()
390 
391  /** \name The integrands, as a function of \f$ u=k/T \f$, for
392  non-degenerate integrals
393  */
394  //@{
395  double density_fun(double u, fermion_deriv &f, double T);
396  double energy_fun(double u, fermion_deriv &f, double T);
397  double entropy_fun(double u, fermion_deriv &f, double T);
398  double density_T_fun(double k, fermion_deriv &f, double T);
399  double density_mu_fun(double k, fermion_deriv &f, double T);
400  double entropy_T_fun(double k, fermion_deriv &f, double T);
401  double density_ms_fun(double k, fermion_deriv &f, double T);
402  //@}
403 
404  /** \name The integrands, as a function of momentum, for the
405  degenerate integrals
406  */
407  //@{
408  double deg_density_fun(double u, fermion_deriv &f, double T);
409  double deg_energy_fun(double u, fermion_deriv &f, double T);
410  double deg_entropy_fun(double u, fermion_deriv &f, double T);
411  double deg_density_T_fun(double k, fermion_deriv &f, double T);
412  double deg_density_mu_fun(double k, fermion_deriv &f, double T);
413  double deg_entropy_T_fun(double k, fermion_deriv &f, double T);
414  double deg_density_ms_fun(double k, fermion_deriv &f, double T);
415  //@}
416 
417  /** \brief Solve for the chemical potential from the density
418  for calc_density()
419  */
420  double solve_fun(double x, fermion_deriv &f, double T);
421 
422  /** \brief Solve for the chemical potential from the density
423  for pair_density()
424  */
425  double pair_fun(double x, fermion_deriv &f, double T);
426 
427 #endif
428 
429  };
430 
431 #ifndef DOXYGEN_NO_O2NS
432 }
433 #endif
434 
435 #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:134
Compute properties of a fermion including derivatives [abstract base].
Definition: part_deriv.h:199
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).