eos_had_rmf.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 /** \file eos_had_rmf.h
24  \brief File defining \ref o2scl::eos_had_rmf
25 */
26 #ifndef O2SCL_RMF_EOS_H
27 #define O2SCL_RMF_EOS_H
28 
29 #include <string>
30 #include <cmath>
31 #include <o2scl/lib_settings.h>
32 #include <o2scl/constants.h>
33 #include <o2scl/mm_funct.h>
34 
35 #include <o2scl/part.h>
36 #include <o2scl/eos_had_base.h>
37 #include <o2scl/fermion.h>
38 
39 #ifndef DOXYGEN_NO_O2NS
40 namespace o2scl {
41 #endif
42 
43  /** \brief Relativistic mean field theory EOS
44 
45  This class computes the properties of nucleonic matter using a
46  mean-field approximation to a field-theoretical model.
47 
48  Before sending neutrons and protons to these member functions,
49  the masses should be set to and the degeneracy factor should be
50  set to 2. Some models which can be loaded using
51  <tt>o2scl_hdf::rmf_load()</tt> expect that the neutron and
52  proton masses are set to the value stored in \ref mnuc.
53 
54  \note Since this EOS uses the effective masses and chemical
55  potentials in the \ref o2scl::part class, the values of
56  <tt>o2scl::part::non_interacting</tt> for neutrons and protons
57  are set to false in many of the functions.
58 
59  \note Matter at two different densities can have the same
60  chemical potentials, so the behavior of the function \ref
61  o2scl::eos_had_rmf::calc_temp_p() is ambiguous. This arises
62  because the field equations have more than one solution for a
63  specified chemical potential. Internally, \ref
64  o2scl::eos_had_rmf::calc_temp_p() either uses the initial guess
65  specified by a call to \ref o2scl::eos_had_rmf::set_fields(), or
66  uses hard-coded initial guess values typical for saturation
67  densities. In order to ensure that the user gets the desired
68  solution to the field equations, it may be necessary to specify
69  a sufficiently accurate initial guess. There is no ambiguity in
70  the behavior of \ref o2scl::eos_had_rmf::calc_eq_temp_p(),
71  however.
72 
73  \note This class can fail to solve the meson field equations or
74  fail to solve for the nucleon densities. By default the error
75  handler is called when this happens. If \ref
76  eos_had_base::err_nonconv is false, then functions which don't
77  converge (which also return <tt>int</tt>) will return a non-zero
78  value. Note that the solvers (in \ref def_sat_mroot and \ref
79  o2scl::eos_had_base::def_mroot) also has its own data member
80  indicating how to handle nonconvergence \ref
81  o2scl::mroot::err_nonconv which is separate.
82 
83  \comment
84  AWS, 11/17/13: It is not clear that this is entirely necessary
85  as almost all the CONV_ERR calls in eos_had_rmf.cpp are due to calls
86  to solvers. It could be that then err_nonconv can be removed and
87  all the eos_had_rmf functions just always directly return any
88  nonzero values they get from solvers. One nice thing about the
89  explicit CONV_ERR calls in eos_had_rmf.cpp is that it makes the code
90  easier to read. In any case err_nonconv should probably be
91  pushed up to eos_had_base.
92  AWS: 6/14/17: Moved err_nonconv up to eos_had_base now because
93  I need it in eos_had_skyrme and it's needed in eos_had_eden_base
94  anyway, so it might as well be in the parent.
95  \endcomment
96 
97  \hline
98  \b Background
99 
100  The full Lagragian can be written as a sum of several terms
101  \f[
102  {\cal L} = {\cal L}_{\mathrm{Dirac}} + {\cal L}_{\sigma} +
103  {\cal L}_{\omega} + {\cal L}_{\rho} + {\cal L}_{\mathrm{int}} \, .
104  \f]
105 
106  The part for the nucleon fields is
107  \f[
108  {\cal L}_{\mathrm{Dirac}} =
109  \bar{\Psi}_i \left[ i {{\partial}\!\!\!{\slash}} -
110  g_{\omega} {{\omega}\!\!\!{\slash}} - \frac{g_{\rho}}{2}
111  {{\vec{\rho}}\!\!\!{\slash}}
112  \vec{\tau} - M_i + g_{\sigma} \sigma -
113  e q_i A\!\!\!{\slash} \right] \Psi_i
114  \f]
115  where \f$ \Psi \f$ is the nucleon field and \f$ \sigma, \omega
116  \f$ and \f$ \rho \f$ are the meson fields. The meson masses
117  are \f$ m_{\sigma}, m_{\omega} \f$ and \f$ m_{\rho}
118  \f$ and meson-nucleon
119  couplings are \f$ g_{\sigma}, g_{\omega} \f$ and \f$ g_{\rho}
120  \f$ . The couplings \c cs, \c cw, and \c cr are related to \f$
121  g_{\sigma}, g_{\omega} \f$ and \f$ g_{\rho} \f$ by
122  \f[
123  c_{\sigma} = g_{\sigma}/m_{\sigma} \quad
124  c_{\omega} = g_{\omega}/m_{\omega} \quad \mathrm{and} \quad
125  c_{\rho} = g_{\rho}/m_{\rho}
126  \f]
127  The nucleon masses are in \f$ M_i \f$ and stored in
128  <tt>part::m</tt> and \f$ q_i \f$ just represents the charge (1
129  for protons and 0 for neutrons). The Coulomb field, \f$ A_{\mu}
130  \f$, is ignored in this class, but used in \ref
131  o2scl::nucleus_rmf.
132 
133  The part for the \f$ \sigma \f$ field is
134  \f[
135  {\cal L}_{\sigma} =
136  {\textstyle \frac{1}{2}} \left( \partial_{\mu} \sigma \right)^2
137  - {\textstyle \frac{1}{2}} m^2_{\sigma} \sigma^2
138  - \frac{b M}{3} \left( g_{\sigma} \sigma\right)^3
139  - \frac{c}{4} \left( g_{\sigma} \sigma\right)^4 \, .
140  \f]
141  where \f$ m_{\sigma} \f$ is the meson mass,
142  \f$ b \f$ and \f$ c \f$ are unitless couplings and
143  \f$ M \f$ is a dimensionful scale, ususally taken to be
144  939 MeV (which need not be equal to \f$ M_i \f$ above).
145  The coefficients \f$ b \f$ and \f$ c \f$ are related to the somewhat
146  standard \f$ \kappa \f$ and \f$ \lambda \f$ by:
147  \f[
148  \kappa=2 M b \quad \lambda=6 c;
149  \f]
150 
151  The part for the \f$ \omega \f$ field is
152  \f[
153  {\cal L}_{\omega} =
154  - {\textstyle \frac{1}{4}} f_{\mu \nu} f^{\mu \nu}
155  + {\textstyle \frac{1}{2}} m^2_{\omega}\omega^{\mu}\omega_{\mu}
156  + \frac{\zeta}{24} g_{\omega}^4 \left(\omega^\mu \omega_\mu\right)^2
157  \f]
158  where \f$ m_{\omega} \f$ is the meson mass.
159 
160  The part for the \f$ \rho \f$ field is
161  \f[
162  {\cal L}_{\rho} =
163  - {\textstyle \frac{1}{4}} \vec{B}_{\mu \nu} \cdot \vec{B}^{\mu \nu}
164  + {\textstyle \frac{1}{2}} m^2_{\rho} \vec{\rho}^{~\mu} \cdot
165  \vec{\rho}_{~\mu}
166  + \frac{\xi}{24} g_{\rho}^4 \left(\vec{\rho}^{~\mu}\right) \cdot
167  \vec{\rho}_{~\mu}
168  \f]
169 
170  Finally, additional meson interactions are
171  \f[
172  {\cal L}_{\mathrm{int}} =
173  g_{\rho}^2 f (\sigma, \omega) \vec{\rho}^{~\mu} \cdot
174  \vec{\rho}_{~\mu} \nonumber \\
175  \f]
176  The function \f$ f \f$ is the coefficient of \f$ g_r^2 \rho^2 \f$
177  \f$ f(\sigma,\omega) = b_1 \omega^2 + b_2 \omega^4 + b_3 \omega^6 +
178  a_1 \sigma + a_2 \sigma^2 + a_3 \sigma^3 + a_4 \sigma^4 +
179  a_5 \sigma^5 + a_6 \sigma^6 \f$
180  where the notation from \ref Horowitz01 is:
181  \f$ f(\sigma,\omega) = \lambda_4 g_s^2 \sigma^2 +
182  \lambda_v g_w^2 \omega^2 \f$
183  This implies \f$ b_1=\lambda_v g_w^2 \f$ and
184  \f$ a_2=\lambda_4 g_s^2 \f$
185 
186  The couplings, \c cs, \c cw, and \c cr all have units of \f$
187  \mathrm{fm} \f$, and the couplings \c b, \c c, \c zeta and \c xi are
188  unitless. The additional couplings from \ref Steiner05b, \f$ a_i
189  \f$ have units of \f$ \mathrm{fm}^{(i-2)} \f$ and the couplings
190  \f$ b_j \f$ have units of \f$ \mathrm{fm}^{(2j-2)} \f$ .
191 
192  When the variable \ref zm_mode is true, the effective mass is
193  fixed using the approach of \ref Zimanyi90 .
194 
195  The expressions for the energy densities are often simplified in
196  the literature using the field equations. These expressions are
197  not used in this code since they are only applicable in infinite
198  matter where the field equations hold, and are not suitable for
199  use in applications (such as to finite nuclei in \ref
200  o2scl::nucleus_rmf) where the spatial derivatives of the fields
201  are non-zero. Notice that in the proper expressions for the
202  energy density the similarity between terms in the pressure up
203  to a sign. This procedure allows one to verify the thermodynamic
204  identity even if the field equations are not solved and allows
205  the user to add gradient terms to the energy density and
206  pressure.
207 
208  See also \ref Muller96 .
209 
210  \hline
211  \b Field \b equations
212 
213  The field equations are:
214  \f[
215  0 = m_{\sigma}^2 \sigma - g_{\sigma} \left( n_{s n} + n_{s p} \right)
216  + b M g_{\sigma}^3 \sigma^2 + c g_{\sigma}^4 \sigma^3 -
217  g_{\rho}^2 \rho^2 \frac{\partial f}{\partial \sigma}
218  \f]
219  \f[
220  0 = m_{\omega}^2 \omega - g_{\omega} \left(n_n+n_p\right)
221  + \frac{\zeta}{6} g_{\omega}^4 \omega^3 + g_{\rho}^2 \rho^2
222  \frac{\partial f}{\partial \omega}
223  \f]
224  \f[
225  0 = m_{\rho}^2 \rho + \frac{1}{2} g_{\rho} \left(n_n-n_p\right)
226  + 2 g_{\rho}^2 \rho f + \frac{\xi}{6} g_{\rho}^4 \rho^3
227  \f]
228 
229  \hline
230  \b Saturation \b properties
231 
232  Defining
233  \f[
234  U(\sigma)=\frac{1}{2} m_\sigma^2\sigma^2+\frac{b M}{3}(g_\sigma\sigma)^3
235  +\frac{c}{4}(g_\sigma\sigma)^4\;,
236  \f]
237  the binding energy per particle in symmetric matter at equilibrium
238  is given by
239  \f[
240  \frac{E}{A} = \frac{1}{n_0} \left[U(\sigma_0)+
241  \frac{1}{2} m_\omega\omega_0^2+
242  \frac{\zeta}{8}(g_\omega\omega_0)^4+\frac{2}{\pi^2}
243  \int\limits_0^{k_F} dk k^2\sqrt{k^2+M^{*2}} \right]
244  \f]
245  where the Dirac
246  effective mass is \f$ M^{*}_i = M_i - g_{\sigma}\sigma_0 \f$ .
247  The compressibility is given by
248  \f[
249  K=9\frac{g_\omega^2}{m_\omega^2}n_0+3\frac{k_F^2}{E_F^*}
250  -9n_0\frac{M^{*2}}{E_F^{*2}}\left[\left(\frac{1}{g_\sigma^2}
251  \frac{\partial^2}{\partial\sigma_0^2}+\frac{3}{g_\sigma M^*}
252  \frac{\partial}{\partial\sigma_0}\right)
253  U(\sigma_0)-3\frac{n_0}{E_F^*}\right]^{-1}\;.
254  \f]
255  The symmetry energy of bulk matter is given by
256  \f[
257  E_{sym} = \frac{k_F^2}{6 E_F^{*}} + \frac{ n }
258  {8 \left(g_{\rho}^2/m_{\rho}^2 + 2 f (\sigma_0, \omega_0)
259  \right)} \, .
260  \f]
261 
262  In the above equations, the subscipt \f$ 0 \f$ denotes the mean
263  field values of \f$ \sigma \f$ and \f$ \omega \f$ . For the case
264  \f$ f=0 \f$ , the symmetry energy varies linearly with the density at
265  large densities. The function \f$ f \f$ permits variations in the
266  density dependence of the symmetry energy above nuclear matter
267  density.
268 
269  \hline
270 
271  \todo
272  - The functions fcomp_fields(), fkprime_fields(), and fesym_fields()
273  are not quite correct if the neutron and proton masses are different.
274  For this reason, they are currently unused by saturation().
275  - The fix_saturation() and calc_cr() functions use mnuc, and should
276  be modified to allow different neutron and proton masses.
277  - Check the formulas in the "Background" section
278  - Make sure that this class properly handles particles for which
279  inc_rest_mass is true/false
280  - The calc_e() function fails to converge at lower densities.
281  See the testing code which has trouble with NL3 and RAPR
282 
283  \future
284  - Finish putting the err_nonconv system into calc_p(),
285  calc_temp_e() and fix_saturation(), etc.
286  - It might be nice to remove explicit reference to the meson
287  masses in functions which only compute nuclear matter since they
288  are unnecessary. This might, however, demand redefining some of
289  the couplings.
290  - Fix calc_p() to be better at guessing
291  - The number of couplings is getting large, maybe new
292  organization is required.
293  - Overload eos_had_base::fcomp() with an exact version
294  - It would be nice to analytically compute the Jacobian
295  of the field equations for the solver
296 
297  */
299 
300  public:
301 
302  /// \name Other data members
303  //@{
304  /** \brief The number of separate calls to the solver
305  that the <tt>calc_e</tt> functions take (default 20)
306 
307  Values larger than about \f$ 10^4 \f$ are probably
308  not useful.
309  */
310  size_t calc_e_steps;
311 
312  /** \brief If true, solve for relative densities rather than
313  absolute densities (default false)
314 
315  Setting this to true makes \ref calc_temp_e() and \ref
316  calc_e() more accurate at low densities.
317  */
319 
320  /// Modifies method of calculating effective masses (default false)
321  bool zm_mode;
322 
323  /** \brief Verbosity parameter
324 
325  If this is greater than zero, then some functions report
326  on their progress.
327  - The function \ref saturation() reports progress towards
328  computing the properties of nuclear matter near saturation.
329  - The functions \ref calc_e() and \ref calc_temp_e() report
330  progress on solving for matter at a fixed density.
331  */
332  int verbose;
333 
334  /** \brief If true, throw exceptions when the function calc_e()
335  does not converge (default true)
336  */
338  //@}
339 
340  /// \name Masses
341  //@{
342  /** \brief The scale \f$ M \f$
343 
344  This need not be exactly equal to the neutron or proton mass,
345  but provides the scale for the coupling \c b.
346  */
347  double mnuc;
348 
349  /// \f$ \sigma \f$ mass (in \f$ \mathrm{fm}^{-1} \f$ )
350  double ms;
351 
352  /// \f$ \omega \f$ mass (in \f$ \mathrm{fm}^{-1} \f$ )
353  double mw;
354 
355  /// \f$ \rho \f$ mass (in \f$ \mathrm{fm}^{-1} \f$ )
356  double mr;
357 
358  //@}
359 
360  /// \name Standard couplings (including nonlinear sigma terms)
361  //@{
362  double cs, cw, cr, b, c;
363  //@}
364 
365  /// \name Quartic terms for omega and rho.
366  //@{
367  double zeta, xi;
368  //@}
369 
370  /// \name Additional isovector couplings
371  //@{
372  double a1, a2, a3, a4, a5, a6, b1, b2, b3;
373  //@}
374 
375  eos_had_rmf();
376 
377  /* \brief Load parameters for model named 'model'
378 
379  Presently accepted values from file rmfdata/model_list:
380  \include rmfdata/model_list
381 
382  In these files, the nucleon and meson masses are by default
383  specified in MeV, and cs, cw, and cr are given in fm. The
384  parameters b and c are both unitless. If the bool 'oakstyle' is
385  true, then load() assumes that gs, gw, and gr have been given
386  where gs and gw are as usual, but gr is a factor of two smaller
387  than usual, and g2 and g3 have been given where g2 = -b M gs^3
388  and g3 = c gs^4. If tokistyle is true, then it is additionally
389  assumed that c3 is given where c3=zeta/6*gw^4.
390 
391  If \c external is true, then model is the filename (relative
392  to the current directory) of the file containing the model
393  parameters. Otherwise, the model is assumed to be present in
394  the \o2 library data directory.
395  */
396  //int load(std::string model, bool external=false);
397 
398  /// \name Compute EOS
399  //@{
400  /** \brief Equation of state as a function of density
401 
402  Initial guesses for the chemical potentials are taken
403  from the user-given values. Initial guesses for the fields
404  can be set by set_fields(), or default values will be used.
405  After the call to calc_e(), the final values of the fields
406  can be accessed through get_fields().
407 
408  This is a little more robust than the standard version
409  in the parent \ref eos_had_base.
410 
411  \future Improve the operation of this function when the
412  proton density is zero.
413  */
414  virtual int calc_e(fermion &ne, fermion &pr, thermo &lth);
415 
416  /** \brief Equation of state as a function of chemical potential
417 
418  Solves for the field equations automatically.
419 
420  \future It may be possible to make the solver for the
421  field equations more robust
422  */
423  virtual int calc_p(fermion &ne, fermion &pr, thermo &lth);
424 
425  /** \brief Equation of state and meson field equations
426  as a function of chemical potentials
427 
428  This calculates the pressure and energy density as a function of
429  \f$ \mu_n,\mu_p,\sigma,\omega,rho \f$ . When the field equations
430  have been solved, \c f1, \c f2, and \c f3 are all zero.
431 
432  The thermodynamic identity is satisfied even when the field
433  equations are not solved.
434 
435  \future Probably best to have f1, f2, and f3 scaled
436  in some sensible way, i.e. scaled to the fields?
437  */
438  virtual int calc_eq_p(fermion &neu, fermion &p, double sig,
439  double ome, double rho, double &f1,
440  double &f2, double &f3, thermo &th);
441 
442  /** \brief Equation of state and meson field equations as a
443  function of chemical potentials at finite temperature
444 
445  Analogous to \ref calc_eq_p() except at finite temperature.
446  */
447  virtual int calc_eq_temp_p(fermion &ne, fermion &pr, double temper,
448  double sig, double ome, double rho, double &f1,
449  double &f2, double &f3, thermo &th);
450 
451  /** \brief Equation of state as a function of chemical potential
452 
453  Solves for the field equations automatically.
454  */
455  virtual int calc_temp_p(fermion &ne, fermion &pr, double T,
456  thermo &lth);
457 
458  /** \brief Equation of state as a function of densities at
459  finite temperature
460  */
461  int calc_temp_e(fermion &ne, fermion &pr, double T,
462  thermo &lth);
463  //@}
464 
465  /// \name Saturation properties
466  //@{
467  /** \brief Calculate cs, cw, cr, b, and c from the saturation
468  properties
469 
470  Note that the meson masses and \ref mnuc must be specified
471  before calling this function.
472 
473  This function does not give correct results when bool zm_mode
474  is true.
475 
476  \c guess_cs, \c guess_cw, \c guess_b, and \c guess_c are
477  initial guesses for \c cs, \c cw, \c b, and \c c respectively.
478 
479  \todo
480  - Fix this for zm_mode=true
481  - Ensure solver is more robust
482 
483  */
484  int fix_saturation(double guess_cs=4.0, double guess_cw=3.0,
485  double guess_b=0.001, double guess_c=-0.001);
486 
487  /** \brief Calculate properties of nuclear matter at the
488  saturation density
489 
490  This function first constructs an initial guess, increasing
491  the chemical potentials if required to ensure the neutron and
492  proton densities are finite, and then uses \ref
493  eos_had_rmf::sat_mroot to solve the field equations and ensure
494  that the neutron and proton densities are equal and the
495  pressure is zero. The quantities \ref eos_had_base::n0, \ref
496  eos_had_base::eoa, and \ref eos_had_base::msom can be computed
497  directly, and the compressibility, the skewness, and the
498  symmetry energy are computed using the functions
499  fkprime_fields() and fesym_fields(). This function overrides
500  the generic version in \ref eos_had_base.
501 
502  If \ref verbose is greater than zero, then then this function
503  reports details on the initial iterations to get the initial
504  guess for the solver.
505  */
506  virtual int saturation();
507 
508  /** \brief Calculate symmetry energy assuming the field
509  equations have already been solved
510 
511  This may only work at saturation density and may assume
512  equal neutron and proton masses.
513  */
514  double fesym_fields(double sig, double ome, double nb);
515 
516  /** \brief Calculate the compressibility assuming the field
517  equations have already been solved
518 
519  This may only work at saturation density and may assume
520  equal neutron and proton masses.
521  */
522  double fcomp_fields(double sig, double ome, double nb);
523 
524  /** \brief Calculate compressibilty and \c kprime assuming the field
525  equations have already been solved
526 
527  This may only work at saturation density and may assume
528  equal neutron and proton masses.
529 
530  \todo This function, \ref o2scl::eos_had_rmf::fkprime_fields() is
531  currently untested.
532  */
533  void fkprime_fields(double sig, double ome, double nb,
534  double &k, double &kprime);
535  //@}
536 
537  /// \name Fields and field equations
538  //@{
539  /** \brief A function for solving the field equations
540 
541  The values <tt>x[0], x[1]</tt>, and <tt>x[2]</tt> should be
542  set to \f$ \sigma, \omega \f$ , and \f$ \rho \f$ on input (in
543  \f$ \mathrm{fm}^{-1} \f$ ) and on exit, <tt>y[0], y[1]</tt>
544  and <tt>y[2]</tt> contain the field equations and are zero
545  when the field equations have been solved.
546  */
547  int field_eqs(size_t nv, const ubvector &x, ubvector &y);
548 
549  /** \brief A function for solving the field equations at finite
550  temperature
551 
552  The values <tt>x[0], x[1]</tt>, and <tt>x[2]</tt> should be
553  set to \f$ \sigma, \omega \f$ , and \f$ \rho \f$ on input (in
554  \f$ \mathrm{fm}^{-1} \f$ ) and on exit, <tt>y[0], y[1]</tt>
555  and <tt>y[2]</tt> contain the field equations and are zero
556  when the field equations have been solved.
557  */
558  int field_eqsT(size_t nv, const ubvector &x, ubvector &y);
559 
560  /** \brief Set a guess for the fields for the next call to calc_e(),
561  calc_p(), or saturation()
562  */
563  virtual int set_fields(double sig, double ome, double lrho) {
564  sigma=sig;
565  omega=ome;
566  rho=lrho;
567  guess_set=true;
568  return 0;
569  }
570 
571  /** \brief Return the most recent values of the meson fields
572 
573  This returns the most recent values of the meson fields set by
574  a call to \ref saturation(), \ref calc_e(), or
575  \ref calc_p(fermion &, fermion &, thermo &).
576  */
577  int get_fields(double &sig, double &ome, double &lrho) {
578  sig=sigma;
579  ome=omega;
580  lrho=rho;
581  return 0;
582  }
583  //@}
584 
585  /// Return string denoting type ("eos_had_rmf")
586  virtual const char *type() { return "eos_had_rmf"; }
587 
588  /// \name Solver
589  //@{
590  /** \brief Set class mroot object for use calculating saturation density
591  */
593  sat_mroot=&mrx;
594  return 0;
595  }
596 
597  /** \brief The default solver for calculating the saturation
598  density
599 
600  Used by fn0() (which is called by saturation()) to solve
601  saturation_matter_e() (1 variable).
602  */
604  //@}
605 
606  /// \name Functions dealing with naturalness
607  //@{
608  /** \brief Set the coefficients of a eos_had_rmf object to their
609  limits from naturalness
610 
611  As given in \ref Muller96 .
612 
613  The definition of the vector-isovector field and coupling
614  matches what is done here. Compare the Lagrangian above
615  with Eq. 10 from the reference.
616 
617  The following couplings should all be of the same
618  size:
619  \f[
620  \frac{1}{2 c_s^2 M^2}, \frac{1}{2 c_v^2 M^2}
621  \frac{1}{8 c_{\rho}^2 M^2},~\mathrm{and}~\frac{
622  \bar{a}_{ijk} M^{i+2 j+2 k-4}}{2^{2 k}}
623  \f]
624  which are equivalent to
625  \f[
626  \frac{m_s^2}{2 g_s^2 M^2}, \frac{m_v^2}{2 g_v^2 M^2}
627  \frac{m_{\rho}^2}{8 g_{\rho}^2 M^2},~\mathrm{and}~\frac{
628  a_{ijk} M^{i+2 j+2 k-4}}{g_s^i g_v^{2 j}
629  g_{\rho}^{2 k} 2^{2 k}}
630  \f]
631 
632  The connection the \f$ a_{ijk} \f$ 's and the coefficients
633  that are used here is
634  \f{eqnarray*}
635  \frac{b M}{3} g_{\sigma}^3 \sigma^3 &=& a_{300}~\sigma^3
636  \nonumber \\
637  \frac{c}{4} g_{\sigma}^4 \sigma^4 &=& a_{400}~\sigma^4
638  \nonumber \\
639  \frac{\zeta}{24} g_{\omega}^4 \omega^4 &=& a_{020}~\omega^4
640  \nonumber \\
641  \frac{\xi}{24} g_{\rho}^4 \rho^4 &=& a_{002}~\rho^4
642  \nonumber \\
643  b_1 g_{\rho}^2 \omega^2 \rho^2 &=& a_{011}~\omega^2 \rho^2
644  \nonumber \\
645  b_2 g_{\rho}^2 \omega^4 \rho^2 &=& a_{021}~\omega^4 \rho^2
646  \nonumber \\
647  b_3 g_{\rho}^2 \omega^6 \rho^2 &=& a_{031}~\omega^6 \rho^2
648  \nonumber \\
649  a_1 g_{\rho}^2 \sigma^1 \rho^2 &=& a_{101}~\sigma^1 \rho^2
650  \nonumber \\
651  a_2 g_{\rho}^2 \sigma^2 \rho^2 &=& a_{201}~\sigma^2 \rho^2
652  \nonumber \\
653  a_3 g_{\rho}^2 \sigma^3 \rho^2 &=& a_{301}~\sigma^3 \rho^2
654  \nonumber \\
655  a_4 g_{\rho}^2 \sigma^4 \rho^2 &=& a_{401}~\sigma^4 \rho^2
656  \nonumber \\
657  a_5 g_{\rho}^2 \sigma^5 \rho^2 &=& a_{501}~\sigma^5 \rho^2
658  \nonumber \\
659  a_6 g_{\rho}^2 \sigma^6 \rho^2 &=& a_{601}~\sigma^6 \rho^2
660  \nonumber
661  \f}
662 
663  Note that Muller and Serot use the notation
664  \f[
665  \frac{\bar{\kappa} g_s^3 }{2} = \frac{\kappa}{2} = b M
666  g_s^3 \qquad \mathrm{and} \qquad
667  \frac{\bar{\lambda} g_s^4}{6} = \frac{\lambda}{6}
668  = c g_s^4
669  \f]
670  which differs slightly from the "standard" notation above.
671 
672  We need to compare the values of
673  \f{eqnarray*}
674  &\frac{m_s^2}{2 g_s^2 M^2}, \frac{m_v^2}{2 g_v^2 M^2}
675  \frac{m_{\rho}^2}{8 g_{\rho}^2 M^2},\frac{b}{3},
676  \frac{c}{4}
677  &
678  \nonumber \\
679  &\frac{\zeta}{24}, \frac{\xi}{384},
680  \frac{b_1}{4 g_{\omega}^2},
681  \frac{b_2 M^2}{4 g_{\omega}^4},
682  \frac{b_3 M^4}{4 g_{\omega}^6},
683  \frac{a_1}{4 g_{\sigma} M},&
684  \nonumber \\
685  &\frac{a_2}{4 g_{\sigma}^2},
686  \frac{a_3 M}{4 g_{\sigma}^3},
687  \frac{a_4 M^2}{4 g_{\sigma}^4},
688  \frac{a_5 M^3}{4 g_{\sigma}^5},~\mathrm{and}~\frac{a_6 M^4}
689  {4 g_{\sigma}^6}\, .&
690  \f}
691 
692  These values are stored in the variables cs, cw, cr, b, c,
693  zeta, xi, b1, etc. in the specified \ref eos_had_rmf object. All
694  of the numbers should be around 0.001 or 0.002.
695 
696  For the scale \f$ M \f$, \ref mnuc is used.
697 
698  \todo I may have ignored some signs in the above, which are
699  unimportant for this application, but it would be good to fix
700  them for posterity.
701 
702  */
704 
705  double gs=cs*ms;
706  double gw=cw*mw;
707  double gr=cr*mr;
708 
709  re.cs=0.5/cs/cs/mnuc/mnuc;
710  re.cw=0.5/cw/cw/mnuc/mnuc;
711  re.cr=0.125/cr/cr/mnuc/mnuc;
712  re.b=b/3.0;
713  re.c=c/4.0;
714 
715  re.zeta=zeta/24.0;
716  re.xi=xi/384.0;
717 
718  re.b1=b1/gw/gw/4.0;
719  re.b2=b2/pow(gw,4.0)/4.0*mnuc*mnuc;
720  re.b3=b3/pow(gw,6.0)/4.0*pow(mnuc,4.0);
721 
722  re.a1=a1/gs/4.0/mnuc;
723  re.a2=a2/pow(gs,2.0)/4.0;
724  re.a3=a3/pow(gs,3.0)/4.0*mnuc;
725  re.a4=a4/pow(gs,4.0)/4.0*mnuc*mnuc;
726  re.a5=a5/pow(gs,5.0)/4.0*pow(mnuc,3.0);
727  re.a6=a6/pow(gs,6.0)/4.0*pow(mnuc,4.0);
728 
729  return;
730  }
731 
732  /** \brief Provide the maximum values of the couplings assuming
733  a limit on naturalness
734 
735  The limits for the couplings are function of the nucleon and
736  meson masses, except for the limits on \c b, \c c, \c zeta,
737  and \c xi which are independent of the masses because of the
738  way that these four couplings are defined.
739  */
740  void naturalness_limits(double value, eos_had_rmf &re) {
741 
742  double gs=cs*ms;
743  double gw=cw*mw;
744  double gr=cr*mr;
745 
746  re.cs=value*2.0*mnuc*mnuc;
747  re.cw=value*2.0*mnuc*mnuc;
748  re.cr=value*8.0*mnuc*mnuc;
749  re.b=value*3.0;
750  re.c=value*4.0;
751 
752  re.zeta=value*24.0;
753  re.xi=value*384.0;
754 
755  re.b1=value*gw*gw*4.0;
756  re.b2=value*pow(gw,4.0)*4.0/mnuc/mnuc;
757  re.b3=value*pow(gw,6.0)*4.0/pow(mnuc,4.0);
758 
759  re.a1=value*gs*4.0*mnuc;
760  re.a2=value*pow(gs,2.0)*4.0;
761  re.a3=value*pow(gs,3.0)*4.0/mnuc;
762  re.a4=value*pow(gs,4.0)*4.0/mnuc/mnuc;
763  re.a5=value*pow(gs,5.0)*4.0/pow(mnuc,3.0);
764  re.a6=value*pow(gs,6.0)*4.0/pow(mnuc,4.0);
765 
766  return;
767  }
768  //@}
769 
770 #ifndef DOXYGEN_INTERNAL
771 
772  protected:
773 
774  /** \brief Temporary baryon density
775  */
776  double n_baryon;
777 
778  /** \brief Temporary charge density
779 
780  \future Should use eos_had_base::proton_frac instead?
781  */
782  double n_charge;
783 
784  /// \name The meson fields
785  //@{
786  double sigma, omega, rho;
787  //@}
788 
789  /// Temperature for solving field equations at finite temperature
790  double fe_temp;
791 
792  /// For calc_e(), if true, then solve for neutron matter
794 
795  /// For calc_e(), if true, then solve for proton matter
797 
798  /// True if a guess for the fields has been given
799  bool guess_set;
800 
801  /// The solver to compute saturation properties
803 
804  /// The function for fix_saturation()
805  int fix_saturation_fun(size_t nv, const ubvector &x, ubvector &y);
806 
807  /// Compute matter at zero pressure (for saturation())
808  virtual int zero_pressure(size_t nv, const ubvector &ex,
809  ubvector &ey);
810 
811  /// The function for calc_e()
812  virtual int calc_e_solve_fun(size_t nv, const ubvector &ex,
813  ubvector &ey);
814 
815  /// The function for calc_temp_e()
816  virtual int calc_temp_e_solve_fun(size_t nv, const ubvector &ex,
817  ubvector &ey);
818 
819  /** \brief Calculate the \c cr coupling given \c sig and \c ome
820  at the density 'nb'.
821 
822  Used by fix_saturation().
823  */
824  int calc_cr(double sig, double ome, double nb);
825 
826  /// Temperature storage for calc_temp_e()
827  double ce_temp;
828 
829 #endif
830 
831  };
832 
833 #ifndef DOXYGEN_NO_O2NS
834 }
835 #endif
836 
837 #endif
mroot< mm_funct, ubvector, jac_funct > * sat_mroot
The solver to compute saturation properties.
Definition: eos_had_rmf.h:802
virtual int calc_eq_temp_p(fermion &ne, fermion &pr, double temper, double sig, double ome, double rho, double &f1, double &f2, double &f3, thermo &th)
Equation of state and meson field equations as a function of chemical potentials at finite temperatur...
virtual int calc_p(fermion &ne, fermion &pr, thermo &lth)
Equation of state as a function of chemical potential.
double ms
mass (in )
Definition: eos_had_rmf.h:350
void check_naturalness(eos_had_rmf &re)
Set the coefficients of a eos_had_rmf object to their limits from naturalness.
Definition: eos_had_rmf.h:703
double mnuc
The scale .
Definition: eos_had_rmf.h:347
bool guess_set
True if a guess for the fields has been given.
Definition: eos_had_rmf.h:799
mroot_hybrids def_sat_mroot
The default solver for calculating the saturation density.
Definition: eos_had_rmf.h:603
virtual int set_fields(double sig, double ome, double lrho)
Set a guess for the fields for the next call to calc_e(), calc_p(), or saturation() ...
Definition: eos_had_rmf.h:563
bool ce_prot_matter
For calc_e(), if true, then solve for proton matter.
Definition: eos_had_rmf.h:796
double fesym_fields(double sig, double ome, double nb)
Calculate symmetry energy assuming the field equations have already been solved.
virtual int calc_eq_p(fermion &neu, fermion &p, double sig, double ome, double rho, double &f1, double &f2, double &f3, thermo &th)
Equation of state and meson field equations as a function of chemical potentials. ...
double mr
mass (in )
Definition: eos_had_rmf.h:356
double fe_temp
Temperature for solving field equations at finite temperature.
Definition: eos_had_rmf.h:790
bool zm_mode
Modifies method of calculating effective masses (default false)
Definition: eos_had_rmf.h:321
double n_baryon
Temporary baryon density.
Definition: eos_had_rmf.h:776
Relativistic mean field theory EOS.
Definition: eos_had_rmf.h:298
A hadronic EOS at finite temperature based on a function of the chemical potentials [abstract base]...
size_t calc_e_steps
The number of separate calls to the solver that the calc_e functions take (default 20) ...
Definition: eos_had_rmf.h:310
virtual int saturation()
Calculate properties of nuclear matter at the saturation density.
bool ce_neut_matter
For calc_e(), if true, then solve for neutron matter.
Definition: eos_had_rmf.h:793
bool calc_e_relative
If true, solve for relative densities rather than absolute densities (default false) ...
Definition: eos_had_rmf.h:318
double mw
mass (in )
Definition: eos_had_rmf.h:353
virtual int calc_e(fermion &ne, fermion &pr, thermo &lth)
Equation of state as a function of density.
int fix_saturation(double guess_cs=4.0, double guess_cw=3.0, double guess_b=0.001, double guess_c=-0.001)
Calculate cs, cw, cr, b, and c from the saturation properties.
void fkprime_fields(double sig, double ome, double nb, double &k, double &kprime)
Calculate compressibilty and kprime assuming the field equations have already been solved...
virtual int zero_pressure(size_t nv, const ubvector &ex, ubvector &ey)
Compute matter at zero pressure (for saturation())
int calc_temp_e(fermion &ne, fermion &pr, double T, thermo &lth)
Equation of state as a function of densities at finite temperature.
virtual int calc_temp_e_solve_fun(size_t nv, const ubvector &ex, ubvector &ey)
The function for calc_temp_e()
double ce_temp
Temperature storage for calc_temp_e()
Definition: eos_had_rmf.h:827
virtual int calc_temp_p(fermion &ne, fermion &pr, double T, thermo &lth)
Equation of state as a function of chemical potential.
virtual int calc_e_solve_fun(size_t nv, const ubvector &ex, ubvector &ey)
The function for calc_e()
double kprime
Skewness in .
Definition: eos_had_base.h:352
bool err_nonconv
If true, throw exceptions when the function calc_e() does not converge (default true) ...
Definition: eos_had_rmf.h:337
int calc_cr(double sig, double ome, double nb)
Calculate the cr coupling given sig and ome at the density &#39;nb&#39;.
double fcomp_fields(double sig, double ome, double nb)
Calculate the compressibility assuming the field equations have already been solved.
int field_eqsT(size_t nv, const ubvector &x, ubvector &y)
A function for solving the field equations at finite temperature.
int field_eqs(size_t nv, const ubvector &x, ubvector &y)
A function for solving the field equations.
int fix_saturation_fun(size_t nv, const ubvector &x, ubvector &y)
The function for fix_saturation()
virtual int set_sat_mroot(mroot< mm_funct, ubvector, jac_funct > &mrx)
Set class mroot object for use calculating saturation density.
Definition: eos_had_rmf.h:592
virtual const char * type()
Return string denoting type ("eos_had_rmf")
Definition: eos_had_rmf.h:586
double n_charge
Temporary charge density.
Definition: eos_had_rmf.h:782
void naturalness_limits(double value, eos_had_rmf &re)
Provide the maximum values of the couplings assuming a limit on naturalness.
Definition: eos_had_rmf.h:740
int get_fields(double &sig, double &ome, double &lrho)
Return the most recent values of the meson fields.
Definition: eos_had_rmf.h:577
int verbose
Verbosity parameter.
Definition: eos_had_rmf.h:332

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