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

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