eos_had_skyrme.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2020, 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_skyrme.h
24  \brief File defining \ref o2scl::eos_had_skyrme
25 */
26 #ifndef O2SCL_SKYRME_EOS_H
27 #define O2SCL_SKYRME_EOS_H
28 
29 #include <iostream>
30 #include <string>
31 #include <cmath>
32 
33 #include <o2scl/constants.h>
34 #include <o2scl/mroot.h>
35 #include <o2scl/eos_had_base.h>
36 #include <o2scl/part.h>
37 #include <o2scl/fermion_nonrel.h>
38 #include <o2scl/fermion_deriv_nr.h>
39 
40 #ifndef DOXYGEN_NO_O2NS
41 namespace o2scl {
42 #endif
43 
44  /** \brief Skyrme hadronic equation of state
45 
46  Equation of state of nucleonic matter based on
47  the Skryme interaction from \ref Skyrme59 .
48 
49  \b Hamiltonian
50 
51  The Hamiltonian is defined (using the notation of
52  \ref Steiner05b ) as
53  \f[
54  {\cal H} =
55  {\cal H}_{k1} +
56  {\cal H}_{k2} +
57  {\cal H}_{k3} +
58  {\cal H}_{p1} +
59  {\cal H}_{p2} +
60  {\cal H}_{p3} +
61  {\cal H}_{g1} +
62  {\cal H}_{g2} \, .
63  \f]
64 
65  The kinetic terms are:
66  \f[
67  {\cal H}_{k1} = \frac{\tau_n}{2 m_n} +
68  \frac{\tau_p}{2 m_p}
69  \f]
70  \f[
71  {\cal H}_{k2} =
72  n \left(\tau_n + \tau_p \right) \left[ \frac{t_1}{4}
73  \left( 1 + \frac{x_1}{2} \right)
74  + \frac{t_2}{4} \left( 1 + \frac{x_2}{2} \right) \right]
75  \f]
76  \f[
77  {\cal H}_{k3} =
78  \left( \tau_n n_n + \tau_p n_p \right) \left[ \frac{t_2}{4}
79  \left( \frac{1}{2} + x_2 \right)
80  - \frac{t_1}{4} \left( \frac{1}{2} + x_1 \right) \right]
81  \f]
82  where \f$ \tau_i \f$ is the Fermi gas energy density
83  of particle \f$ i \f$ .
84 
85  The potential terms are:
86  \f[
87  {\cal H}_{p1} =
88  \frac{t_0}{2}
89  \left[ \left( 1 + \frac{x_0}{2} \right) n^2 -
90  \left( {\textstyle \frac{1}{2}} + x_0 \right)
91  \left( n_n^2 + n_p^2 \right) \right]
92  \f]
93  \f[
94  {\cal H}_{p2} =
95  \frac{a t_3}{6} \left[ \left( 1 + \frac{x_3}{2} \right) n^{\alpha}
96  n_n n_p + 2^{\alpha-2} \left(1 - x_3\right)
97  \left(n_n^{\alpha+2} + n_p^{\alpha+2}\right) \right]
98  \f]
99  \f[
100  {\cal H}_{p3} =
101  \frac{b t_3}{12} \left[ \left(1 + \frac{x_3}{2} \right) n^{\alpha+2} -
102  \left(\frac{1}{2} + x_3 \right) n^{\alpha}
103  \left( n_n^2+n_p^2 \right) \right]
104  \f]
105 
106  \b Gradient \b and \b Spin-Orbit \b Terms
107 
108  The gradient terms are displayed here for completeness even though
109  they are not computed in the code:
110  \f[
111  {\cal H}_{g1} =
112  \frac{3}{32} \left[ t_1 \left(1 - x_1 \right) -
113  t_2 \left(1 + x_2 \right) \right] \left[ \left( \nabla n_n\right)^2 +
114  \left( \nabla n_p \right)^2 \right]
115  \f]
116  \f[
117  {\cal H}_{g2} =
118  \frac{1}{8} \left[ 3 t_1 \left( 1 +
119  \frac{x_1}{2} \right) - t_2 \left(1 + \frac{x_2}{2} \right) \right]
120  \nabla n_n \nabla n_p
121  \f]
122 
123  The values \f$ a=0, b=1 \f$ give the standard definition of the
124  Skyrme Hamiltonian \ref Skyrme59, while \f$a=1, b=0\f$ contains
125  the modifications suggested by \ref Onsi94.
126 
127  The spin-orbit term is (following \ref Steiner05)
128  \f[
129  {\cal H}_{J} = -\frac{W_0}{2} \left( n_n \vec{\nabla} \cdot
130  \vec{J}_n + n_p \vec{\nabla} \cdot \vec{J}_p + n \vec{\nabla}
131  \cdot \vec{J} \right) + \frac{t_1}{16} \left(\vec{J}_n^2 +
132  \vec{J}_p^2 - x_1 \vec{J}^2\right) - \frac{t_2}{16}
133  \left(\vec{J}_n^2 + \vec{J}_p^2 + x_2 \vec{J}^2\right)
134  \f]
135  where sometimes the \f$ \vec{J}^2 \f$ terms are not included.
136  Alternatively, one can separate the isoscalar and isovector
137  parts in the first term
138  \f[
139  {\cal H}_{J} = - b_4 n \vec{\nabla} \cdot \vec{J} -
140  b_4^{\prime} n_n \vec{\nabla} \cdot \vec{J}_n -
141  b_4^{\prime} n_p \vec{\nabla} \cdot \vec{J}_p
142  \f]
143  then the earlier Skyrme interactions have \f$ b_4 =
144  b_4^{\prime} = W_0/2 \f$. For example, for SLy4,
145  \f$ b_4 = b_4^{\prime} = W_0/2 = 61.5~\mathrm{MeV} \f$.
146 
147  Three quantities are defined in \ref Steiner05b for
148  use in computing the properties of matter at saturation
149  \f[
150  t_3^{\prime} = \left(a + b\right) t_3 \, ,
151  \f]
152  \f[
153  C = \frac{3 }{10 m} \left( \frac{3 \pi^2 }
154  {2} \right)^{2/3} \, ,
155  \f]
156  and
157  \f[
158  \beta = \frac{M}{2} \left[ \frac{1}{4}
159  \left( 3 t_1 + 5 t_2 \right)
160  + t_2 x_2 \right] \, . \\
161  \f]
162 
163  \b Units
164 
165  Quantities which have units containing powers of energy are
166  divided by \f$\hbar c\f$ to ensure all quantities are in units
167  of \f$ \mathrm{fm} \f$. The \f$x_i\f$ and \f$\alpha\f$ are
168  unitless, while the original units of the \f$t_i\f$ are:
169  - \f$t_0\f$ - \f$\mathrm{MeV}\f$ \f$\mathrm{fm}^3\f$
170  - \f$t_1\f$ - \f$\mathrm{MeV}\f$ \f$\mathrm{fm}^5\f$
171  - \f$t_2\f$ - \f$\mathrm{MeV}\f$ \f$\mathrm{fm}^5\f$
172  - \f$t_3\f$ - \f$\mathrm{MeV}\f$ \f$\mathrm{fm}^{3(1+\alpha)}\f$
173 
174  These are stored internally with units of:
175  - \f$t_0\f$ - \f$\mathrm{fm}^2\f$
176  - \f$t_1\f$ - \f$\mathrm{fm}^4\f$
177  - \f$t_2\f$ - \f$\mathrm{fm}^4\f$
178  - \f$t_3\f$ - \f$\mathrm{fm}^{2+3 \alpha}\f$
179 
180  \b Other \b Notes
181 
182  The functions for the usual saturation properties are based
183  partly on \ref Brack85.
184 
185  Models are taken from the references: \ref Bartel79, \ref
186  Beiner75, \ref Chabanat95, \ref Chabanat97, \ref Danielewicz09,
187  \ref Dobaczewski94, \ref Dutta86, \ref Friedrich86, \ref Onsi94,
188  \ref Reinhard95, and \ref Tondeur84, and \ref VanGiai81 .
189 
190  The variables \f$ \nu_n \f$ and \f$ \nu_p \f$ contain the
191  expressions \f$ (-\mu_n+V_n)/T \f$ and \f$ (-\mu_p+V_p)/T \f$
192  respectively, where \f$ V \f$ is the potential part of the
193  single particle energy for particle i (i.e. the derivative of
194  the Hamiltonian w.r.t. density while energy density held
195  constant). Equivalently, \f$ \nu_n\f$ is just \f$ -k_{F_n}^2/ 2
196  m_n^{*} \f$.
197 
198  \note The finite temperature code does not attempt to include
199  antiparticles and uses \ref
200  o2scl::fermion_nonrel_tl::calc_density(). At finite temperature,
201  pure neutron matter implies a zero proton number density which
202  means the proton chemical potential is \f$ - \infty \f$ and thus
203  set to the additive inverse of the return value of the function
204  <tt>numeric_limits<double>::infinity()</tt> The case of pure
205  proton matter is handled similarly. Negative densities result in
206  calling the error handler.
207 
208  Skyrme models are loaded using \ref o2scl_hdf::skyrme_load() .
209  The full list is given in the \o2 repository in
210  <tt>o2scl/data/o2scl/skdata/model_list</tt>.
211 
212  \b Todos
213 
214  \todo
215  - Convert W0 to b4 and b4p everywhere
216  - Remove use of mnuc in calparfun()?
217  - Update reference list.
218 
219  \future
220  - This EOS typically converges very well. One exception seems
221  to be using <tt>calc_temp_p()</tt> at very low densities. I have
222  had problems, for example, with <tt>mun=5.0, mup=6.5</tt>
223  at <tt>T=1.0/197.33</tt>.
224 
225  */
227 
228  protected:
229 
230  /// \name EOS helper functions
231  //@{
232  /** \brief Compute the coefficients of the potential energy which
233  have unique dependence on the densities
234  */
235  void hamiltonian_coeffs(double &ham1, double &ham2,
236  double &ham3, double &ham4,
237  double &ham5, double &ham6);
238 
239  /** \brief Compute the base thermodynamic quantities
240 
241  This function computes the energy density, pressure,
242  entropy, and chemical potentials.
243  */
244  template<class fermion_t>
245  void base_thermo
246  (fermion_t &ne, fermion_t &pr, double ltemper, thermo &locth,
247  double term, double term2, double ham1, double ham2,
248  double ham3, double ham4, double ham5, double ham6) {
249 
250  double nb=ne.n+pr.n;
251  double na=pow(nb,alpha);
252  double npa=pow(pr.n,alpha);
253  double nna=pow(ne.n,alpha);
254 
255  double ham=ne.ed+pr.ed+ham1*nb*nb+ham2*(ne.n*ne.n+pr.n*pr.n)+
256  ham3*na*ne.n*pr.n+ham4*(nna*ne.n*ne.n+npa*pr.n*pr.n)+
257  ham5*nb*nb*na+ham6*(ne.n*ne.n+pr.n*pr.n)*na;
258 
259  double gn, gp;
260  if (ne.inc_rest_mass) {
261  gn=2.0*ne.ms*(ne.ed-ne.n*ne.m);
262  } else {
263  gn=2.0*ne.ms*ne.ed;
264  }
265  if (pr.inc_rest_mass) {
266  gp=2.0*pr.ms*(pr.ed-pr.n*pr.m);
267  } else {
268  gp=2.0*pr.ms*pr.ed;
269  }
270 
271  // Variables dhdn{n,p} are the partial derivatives of the
272  // Hamiltonian wrt the neutron and proton densities
273  double common=2.0*ham1*nb+ham5*(2.0+alpha)*nb*na;
274  double dhdnn=common+2.0*ham2*ne.n+ham3*na*pr.n*(alpha*ne.n/nb+1.0)+
275  ham4*(nna*ne.n*(2.0+alpha))+
276  ham6*(2.0*ne.n*na+(ne.n*ne.n+pr.n*pr.n)*alpha*na/nb);
277  double dhdnp=common+2.0*ham2*pr.n+ham3*na*ne.n*(alpha*pr.n/nb+1.0)+
278  ham4*(npa*pr.n*(2.0+alpha))+
279  ham6*(2.0*pr.n*na+(ne.n*ne.n+pr.n*pr.n)*alpha*na/nb);
280 
281  // Compute the chemical potentials
282  ne.mu=ne.nu+dhdnn+(gn+gp)*term+gn*term2;
283  pr.mu=pr.nu+dhdnp+(gn+gp)*term+gp*term2;
284 
285  // Thermodynamics
286  locth.ed=ham;
287  locth.en=ne.en+pr.en;
288  locth.pr=ltemper*locth.en+ne.mu*ne.n+pr.mu*pr.n-locth.ed;
289 
290  return;
291  }
292 
293  /** \brief Compute second derivatives of the free energy
294  */
295  template<class fermion_t>
296  void second_deriv
297  (fermion_t &ne, fermion_t &pr, double ltemper, thermo &locth,
298  thermo_np_deriv_helm &locthd, double term, double term2,
299  double ham1, double ham2, double ham3,
300  double ham4, double ham5, double ham6) {
301 
302  double nb=ne.n+pr.n;
303  double na=pow(nb,alpha);
304  double npa=pow(pr.n,alpha);
305  double nna=pow(ne.n,alpha);
306 
307  // Second derivatives of the potential energy
308 
309  double opatpa=(1.0+alpha)*(2.0+alpha);
310  double common2=2.0*ham1+2.0*ham2;
311  double dhdnn2=common2+ham4*nna*opatpa+ham5*na*opatpa+
312  ham3*(2.0*pr.n/nb*na*alpha+
313  ne.n*pr.n*na/nb/nb*(-1.0+alpha)*alpha)+
314  ham6*(2.0*na+4.0*ne.n/nb*na*alpha+na/nb/nb*
315  (ne.n*ne.n+pr.n*pr.n)*(-1.0+alpha)*alpha);
316  double dhdnp2=common2+ham4*npa*opatpa+ham5*na*opatpa+
317  ham3*(2.0*ne.n/nb*na*alpha+
318  ne.n*pr.n*na/nb/nb*(-1.0+alpha)*alpha)+
319  ham6*(2.0*na+4.0*pr.n/nb*na*alpha+na/nb/nb*
320  (ne.n*ne.n+pr.n*pr.n)*(-1.0+alpha)*alpha);
321  double dhdnndnp=2.0*ham1+na/nb/nb*
322  (ham5*nb*nb*opatpa+ham6*alpha*
323  (4.0*ne.n*pr.n+ne.n*ne.n*(1.0+alpha)+pr.n*pr.n*(1.0+alpha))+
324  ham3*(ne.n*ne.n*(1.0+alpha)+pr.n*pr.n*(1.0+alpha)+
325  ne.n*pr.n*(2.0+alpha+alpha*alpha)));
326 
327  // For the kinetic part, convert from the (mu,T)
328  // to (n,T) representation
329 
330  double n_dsdT_f=0.0, p_dsdT_f=0.0;
331  double n_dmudT_f=0.0, p_dmudT_f=0.0;
332  double n_dmudn_f=0.0, p_dmudn_f=0.0;
333  ne.deriv_f(n_dmudn_f,n_dmudT_f,n_dsdT_f);
334  pr.deriv_f(p_dmudn_f,p_dmudT_f,p_dsdT_f);
335 
336  // The product 10 mstar^2 epsilon
337 
338  double hn, hp;
339  if (ne.inc_rest_mass) {
340  hn=10.0*ne.ms*(ne.ed-ne.n*ne.m)*ne.ms;
341  } else {
342  hn=10.0*ne.ms*ne.ed*ne.ms;
343  }
344  if (pr.inc_rest_mass) {
345  hp=10.0*pr.ms*(pr.ed-pr.n*pr.m)*pr.ms;
346  } else {
347  hp=10.0*pr.ms*pr.ed*pr.ms;
348  }
349 
350  // Now combine to compute the six derivatives
351 
352  locthd.dsdT=n_dsdT_f+p_dsdT_f;
353 
354  locthd.dmundT=2.0*ltemper*ne.ms*(term+term2)*n_dsdT_f+
355  2.0*ltemper*pr.ms*term*p_dsdT_f+n_dmudT_f;
356 
357  locthd.dmupdT=2.0*ltemper*pr.ms*(term+term2)*p_dsdT_f+
358  2.0*ltemper*ne.ms*term*n_dsdT_f+p_dmudT_f;
359 
360  locthd.dmundnn=-(term+term2)*(term+term2)*hn-term*term*hp+
361  pow(1.0+3.0*(term+term2)*ne.n*ne.ms,2.0)*n_dmudn_f+
362  9.0*term*term*pr.n*pr.ms*pr.n*pr.ms*p_dmudn_f+dhdnn2;
363 
364  locthd.dmupdnp=-(term+term2)*(term+term2)*hp-term*term*hn+
365  pow(1.0+3.0*(term+term2)*pr.n*pr.ms,2.0)*p_dmudn_f+
366  9.0*term*term*ne.n*ne.ms*ne.n*ne.ms*n_dmudn_f+dhdnp2;
367 
368  locthd.dmudn_mixed=-term*(term+term2)*(hn+hp)+
369  3.0*term*ne.ms*ne.n*(1.0+3.0*(term+term2)*ne.n*ne.ms)*n_dmudn_f+
370  3.0*term*pr.ms*pr.n*(1.0+3.0*(term+term2)*pr.n*pr.ms)*p_dmudn_f+
371  dhdnndnp;
372 
373  return;
374  }
375 
376  /** Check the EOS input
377 
378  This function checks that
379  - the densities and temperature are finite and positive
380  - the spin denegeracies are correct
381  - the masses are sensible
382  - the values of 'non_interacting' are false
383  - the alpha parameter is positive
384  - the temperature is not negative
385  */
386  template<class fermion_t>
387  void check_input(fermion_t &ne, fermion_t &pr, double T) {
388 
389  if (!std::isfinite(ne.n) || !std::isfinite(pr.n) ||
390  !std::isfinite(T)) {
391  O2SCL_ERR2("Nucleon densities or temperature not finite in ",
392  "eos_had_skyrme::check_input().",exc_einval);
393  }
394  if (ne.n<0.0 || pr.n<0.0) {
395  std::string str=((std::string)"Nucleon densities negative, n_n=")+
396  std::to_string(ne.n)+", n_p="+std::to_string(pr.n)+", in "+
397  "eos_had_skyrme::check_input().";
398  O2SCL_ERR(str.c_str(),exc_einval);
399  }
400  if (fabs(ne.g-2.0)>1.0e-10 || fabs(pr.g-2.0)>1.0e-10) {
401  O2SCL_ERR((((std::string)"Neutron (")+std::to_string(ne.g)+
402  ") or proton ("+std::to_string(pr.g)+") spin deg"+
403  "eneracies wrong in "+
404  "eos_had_skyrme::check_input().").c_str(),
405  exc_einval);
406  }
407  if (fabs(ne.m-4.5)>1.0 || fabs(pr.m-4.5)>1.0) {
408  O2SCL_ERR((((std::string)"Neutron (")+std::to_string(ne.m)+
409  ") or proton ("+std::to_string(pr.m)+") masses wrong "+
410  "in eos_had_skyrme::check_input().").c_str(),
411  exc_einval);
412  }
413  if (ne.non_interacting==true || pr.non_interacting==true) {
414  O2SCL_ERR2("Neutron or protons non-interacting in ",
415  "eos_had_skyrme::check_input().",exc_einval);
416  }
417  if (alpha<=0.0) {
418  O2SCL_ERR2("Parameter alpha negative in ",
419  "eos_had_skyrme::calc_e().",exc_einval);
420  }
421  if (T<0.0) {
422  O2SCL_ERR2("Temperature negative in ",
423  "eos_had_skyrme::calc_e().",exc_einval);
424  }
425 
426  return;
427  }
428  //@}
429 
430  public:
431 
432  /// \name Basic usage
433  //@{
434  /// Create a blank Skyrme EOS
435  eos_had_skyrme();
436 
437  /// Destructor
438  virtual ~eos_had_skyrme() {};
439 
440  /** \brief Equation of state as a function of densities
441  at finite temperature
442  */
443  virtual int calc_temp_e(fermion &ne, fermion &pr, double temper,
444  thermo &th);
445 
446  /** \brief Equation of state as a function of the densities at
447  finite temperature and including second derivatives
448  */
449  virtual int calc_deriv_temp_e(fermion_deriv &ne, fermion_deriv &pr,
450  double temper, thermo &th,
451  thermo_np_deriv_helm &thd);
452 
453  /** \brief Equation of state as a function of densities at
454  zero temperature
455  */
456  virtual int calc_e(fermion &ne, fermion &pr, thermo &lt);
457 
458  /// Return string denoting type ("eos_had_skyrme")
459  virtual const char *type() { return "eos_had_skyrme"; }
460  //@}
461 
462  /// \name Basic Skyrme model parameters
463  //@{
464  double t0, t1, t2, t3, x0, x1, x2, x3, alpha, a, b;
465  //@}
466 
467  /// \name Other parameters
468  //@{
469  /** \brief Spin-orbit splitting (in \f$ \mathrm{fm}^{-1} \f$)
470 
471  This is unused, but included for possible future use and
472  present in the internally stored models.
473  */
474  double W0;
475 
476  /// Isoscalar spin-orbit term (in \f$ \mathrm{fm}^{-1} \f$)
477  double b4;
478 
479  /// Isovector spin-orbit term (in \f$ \mathrm{fm}^{-1} \f$)
480  double b4p;
481 
482  /// Bibliographic reference
483  std::string reference;
484 
485  /** \brief Use eos_had_base methods for saturation properties
486 
487  This can be set to true to check the difference between
488  the exact expressions and the numerical values from
489  class eos_had_base.
490  */
492  //@}
493 
494  /** \name Saturation properties
495 
496  These calculate the various saturation properties exactly from
497  the parameters at any density. These routines often assume that
498  the neutron and proton masses are equal.
499  */
500  //@{
501  /** \brief Calculate binding energy
502 
503  \f[
504  \frac{E}{A} = C n_B^{2/3} \left( 1 + \beta n_B \right) +
505  \frac{3 t_0}{8} n_B + \frac{t_3^{\prime}}{16} n_B^{\alpha+1}
506  \f]
507  */
508  virtual double feoa(double nb);
509 
510  /** \brief Calculate effective mass
511 
512  \f[
513  M^{*}/M = \left(1+ \beta n_B \right)^{-1} \\
514  \f]
515  */
516  virtual double fmsom(double nb);
517 
518  /** \brief Calculate compressibility in nuclear (isospin-symmetric
519  matter)
520 
521  \f[
522  K = 10 C n_B^{2/3} + \frac{27}{4} t_0 n_B + 40 C \beta n_B^{5/3} +
523  \frac{9 t_3^{\prime}}{16}
524  \alpha \left( \alpha+1 \right) n_B^{1 + \alpha} +
525  \frac{9 t_3^{\prime}}{8} \left( \alpha+1 \right) n_B^{1 + \alpha}
526  \f]
527  */
528  virtual double fcomp_nuc(double nb);
529 
530  /** \brief Calculate symmetry energy
531 
532  If pf=0.5, then the exact expression below is used.
533  Otherwise, the method from class eos_had_base is used.
534 
535  \f[
536  E_{sym} = \frac{5}{9} C n^{2/3} + \frac{10 C m}{3}
537  \left[ \frac{t_2}{6} \left(1 + \frac{5}{4} x_2 \right) -
538  \frac{1}{8} t_1 x_1 \right] n^{5/3}
539  - \frac{t_3^{\prime}}{24}
540  \left({\textstyle \frac{1}{2}} + x_3 \right) n^{1+\alpha} -
541  \frac{t_0}{4} \left( {\textstyle \frac{1}{2}} + x_0 \right) n
542  \f]
543  */
544  virtual double fesym(double nb, double alpha=0.0);
545 
546  /** \brief Skewness in nuclear (isospin-symetric) matter
547 
548  \f[
549  2 C n_B^{2/3} \left(9-5/M^{*}/M\right)+
550  \frac{27 t_3^{\prime}}{16} n^{1+\alpha} \alpha
551  \left(\alpha^2-1\right)
552  \f]
553  */
554  virtual double fkprime_nuc(double nb);
555  //@}
556 
557  /// \name Compute and test Landau parameters
558  //@{
559  /** \brief Check the Landau parameters for instabilities
560 
561  This returns zero if there are no instabilities.
562  */
563  int check_landau(double nb, double m);
564 
565  /** \brief Calculate the Landau parameters for nuclear matter
566 
567  Given \c n0 and \c m, this calculates the Landau parameters in
568  nuclear matter as given in \ref Margueron02
569 
570  \todo This needs to be checked.
571 
572  (Checked once on 11/05/03)
573  */
574  void landau_nuclear(double n0, double m,
575  double &f0, double &g0, double &f0p,
576  double &g0p, double &f1, double &g1,
577  double &f1p, double &g1p);
578 
579  /** \brief Calculate the Landau parameters for neutron matter
580 
581  Given 'n0' and 'm', this calculates the Landau parameters in
582  neutron matter as given in \ref Margueron02
583 
584  \todo This needs to be checked
585 
586  (Checked once on 11/05/03)
587  */
588  void landau_neutron(double n0, double m, double &f0, double &g0,
589  double &f1, double &g1);
590 
591  //@}
592 
593  /// \name Other functions
594  //@{
595  /** \brief Evaluate the effective masses for neutrons and
596  protons
597  */
598  template<class fermion_t>
599  void eff_mass(fermion_t &ne, fermion_t &pr, double &term,
600  double &term2) {
601 
602  // Landau effective masses
603  double nb=ne.n+pr.n;
604  term=0.25*(t1*(1.0+x1/2.0)+t2*(1.0+x2/2.0));
605  term2=0.25*(t2*(0.5+x2)-t1*(0.5+x1));
606  ne.ms=ne.m/(1.0+2.0*(nb*term+ne.n*term2)*ne.m);
607  pr.ms=pr.m/(1.0+2.0*(nb*term+pr.n*term2)*pr.m);
608  return;
609  }
610 
611  /** \brief Calculate \f$ t_0,t_1,t_2,t_3 \f$ and \f$ \alpha \f$ from
612  the saturation properties.
613 
614  In nuclear matter:
615 
616  \f$ E_b=E_b(n_0,M^{*},t_0,t_3,\alpha) \f$ \n
617  \f$ P=P(n_0,M^{*},t_0,t_3,\alpha) \f$ \n
618  \f$ K=K(n_0,M^{*},t_3,\alpha) \f$
619  (the \f$ t_0 \f$ dependence vanishes) \n
620  \f$ M^{*}=M^{*}(n_0,t_1,t_2,x_2) \f$
621  (the \f$ x_1 \f$ dependence cancels), \n
622  \f$ E_{sym}=E_{sym}(x_0,x_1,x_2,x_3,t_0,t_1,t_2,t_3,\alpha) \f$
623 
624  To fix the couplings from the saturation properties, we take
625  \f$ n_0, M^{*}, E_b, K \f$ as inputs, and we can fix \f$
626  t_0,t_3,\alpha \f$ from the first three relations, then use
627  \f$ M^{*}, E_b \f$ to fix \f$ t_2 \f$ and \f$ t_1 \f$. The
628  separation into two solution steps should make for better
629  convergence. All of the x's are free parameters and should be
630  set before the function call.
631 
632  The arguments \c gt0, \c gt3, \c galpha, \c gt1, and \c gt2
633  are used as initial guesses for skyme_eos::t0, eos_had_skyrme::t3,
634  eos_had_skyrme::alpha, eos_had_skyrme::t1, and eos_had_skyrme::t2
635  respectively.
636 
637  \todo Does this work for both 'a' and 'b' non-zero?
638 
639  \todo Compare to similar formulas in \ref Margueron02
640  */
641  int calpar(double gt0=-10.0, double gt3=70.0, double galpha=0.2,
642  double gt1=2.0, double gt2=-1.0);
643 
644  // Unfinished.
645  /* \brief
646  From \ref Margueron02
647  */
648  // int calpar_new(double m);
649 
650  /** \brief Set using alternate parameterization
651 
652  From \ref Bender03 as in, e.g. \ref Kortelainen14
653  \f{eqnarray*}
654  C^{\rho \rho}_{00} &=& 3 t_0/8 \nonumber \\
655  C^{\rho \rho}_{10} &=& -t_0/4 \left(\frac{1}{2}+x_0 \right)
656  \nonumber \\
657  C^{\rho \rho}_{0D} &=& t_3/16 \nonumber \\
658  C^{\rho \rho}_{1D} &=& -t_3/24 \left(\frac{1}{2}+x_3\right)
659  \nonumber \\
660  C^{\rho \tau}_{0} &=& 3 t_1/16+t_2/4\left(\frac{5}{4}+x_2\right)
661  \nonumber \\
662  C^{\rho \tau}_{1} &=& -t_1/8 \left(\frac{1}{2}+x_1\right) +
663  t_2/8 \left(\frac{1}{2}+x_2\right) \nonumber \\
664  C^{\rho \Delta \rho}_{0} &=& -9/64 t_1+t_2/16
665  \left(\frac{5}{4}+x_2\right) \nonumber \\
666  C^{\rho \Delta \rho}_{1} &=& 3/32 t_1 \left(\frac{1}{2}+x_1\right) +
667  t_2/32 \left(\frac{1}{2}+x_2\right) \nonumber \\
668  C^{\rho \nabla J}_{0} &=& -b_4 -b_4^{\prime}/2 \nonumber \\
669  C^{\rho \nabla J}_{1} &=& -b_4^{\prime}/2
670  \f}
671 
672  The parameters should have the following units
673  - <tt>Crr00</tt>: \f$ \mathrm{fm}^2 \f$
674  - <tt>Crr10</tt>: \f$ \mathrm{fm}^2 \f$
675  - <tt>Crr0D</tt>: \f$ \mathrm{fm}^{3 \alpha+2} \f$
676  - <tt>Crr1D</tt>: \f$ \mathrm{fm}^{3 \alpha+2} \f$
677  - <tt>Crt0</tt>: \f$ \mathrm{fm}^4 \f$
678  - <tt>Crt1</tt>: \f$ \mathrm{fm}^4 \f$
679  - <tt>CrDr0</tt>: \f$ \mathrm{fm}^4 \f$
680  - <tt>CrDr1</tt>: \f$ \mathrm{fm}^4 \f$
681  - <tt>CrnJ0</tt>: \f$ \mathrm{fm}^{-1} \f$
682  - <tt>CrnJ1</tt>: \f$ \mathrm{fm}^{-1} \f$
683  - <tt>alpha2</tt>: unitless
684 
685  \todo These expressions are not exactly the same
686  as those in \ref Bender03, so I need to find out why
687  and make this more clear.
688  */
689  void alt_params_set
690  (double Crr00, double Crr10, double Crr0D, double Crr1D,
691  double Crt0, double Crt1, double CrDr0, double CrDr1,
692  double CrnJ0, double CrnJ1, double alpha2);
693 
694  /** \brief Get alternate parameterization
695 
696  The parameters will have the following units
697  - <tt>Crr00</tt>: \f$ \mathrm{fm}^2 \f$
698  - <tt>Crr10</tt>: \f$ \mathrm{fm}^2 \f$
699  - <tt>Crr0D</tt>: \f$ \mathrm{fm}^{3 \alpha+2} \f$
700  - <tt>Crr1D</tt>: \f$ \mathrm{fm}^{3 \alpha+2} \f$
701  - <tt>Crt0</tt>: \f$ \mathrm{fm}^4 \f$
702  - <tt>Crt1</tt>: \f$ \mathrm{fm}^4 \f$
703  - <tt>CrDr0</tt>: \f$ \mathrm{fm}^4 \f$
704  - <tt>CrDr1</tt>: \f$ \mathrm{fm}^4 \f$
705  - <tt>CrnJ0</tt>: \f$ \mathrm{fm}^{-1} \f$
706  - <tt>CrnJ1</tt>: \f$ \mathrm{fm}^{-1} \f$
707  - <tt>alpha2</tt>: unitless
708 
709  See \ref eos_had_skyrme::alt_params_set().
710  */
711  void alt_params_get
712  (double &Crr00, double &Crr10, double &Crr0D, double &Crr1D,
713  double &Crt0, double &Crt1, double &CrDr0, double &CrDr1,
714  double &CrnJ0, double &CrnJ1, double &alpha2);
715 
716  /** \brief Use the specified saturation properties and couplings
717  and the function \ref alt_params_set() to set the
718  Skyrme coefficients
719 
720  This function uses the relations in \ref Kortelainen10 .
721  The parameters should have the following units
722  - <tt>n0</tt>: \f$ \mathrm{fm}^{-3} \f$
723  - <tt>EoA</tt>: \f$ \mathrm{fm}^{-1} \f$
724  - <tt>K</tt>: \f$ \mathrm{fm}^{-1} \f$
725  - <tt>Ms_star</tt>: unitless
726  - <tt>a</tt>: \f$ \mathrm{fm}^{-1} \f$
727  - <tt>L</tt>: \f$ \mathrm{fm}^{-1} \f$
728  - <tt>Mv_star</tt>: unitless
729  - <tt>CrDr0</tt>: \f$ \mathrm{fm}^{-3} \f$
730  - <tt>CrDr1</tt>: \f$ \mathrm{fm}^{-3} \f$
731  - <tt>CrnJ0</tt>: \f$ \mathrm{fm}^{-3} \f$
732  - <tt>CrnJ1</tt>: \f$ \mathrm{fm}^{-3} \f$
733 
734  \ref Kortelainen10 assumed equal neutron and proton masses, so
735  this function uses \f$ \hbar^2/(2m) = \hbar^2/(m_n+m_p) \f$
736  and the neutron and proton masses in \ref
737  eos_had_base::def_neutron and \ref eos_had_base::def_proton,
738  respectively. To obtain the results in the original paper, set
739  neutron and proton masses to ensure that \f$ \hbar^2/(2m) =
740  20.73553~\mathrm{MeV}~\mathrm{fm}^2 \f$ .
741  */
743  (double n0, double EoA, double K, double Ms_star, double a, double L,
744  double Mv_star, double CrDr0, double CrDr1, double CrnJ0, double CrnJ1);
745  //@}
746 
747  /// \name Particle classes
748  //@{
749  /// Thermodynamics of non-relativistic fermions
750  fermion_nonrel nrf;
751 
752  /// Thermodynamics of non-relativistic fermions with derivatives
753  fermion_deriv_nr nrfd;
754  //@}
755 
756 #ifndef DOXYGEN_NO_O2NS
757 
758  protected:
759 
760  /// \name Functions and parameters for calpar()
761  //@{
762  int calparfun(size_t nv, const ubvector &x, ubvector &y);
763  int calparfun2(size_t nv, const ubvector &x, ubvector &y);
764  double fixn0, fixeoa, fixesym, fixcomp, fixmsom;
765  //@}
766 
767 #endif
768 
769  };
770 
771 #ifndef DOXYGEN_NO_O2NS
772 }
773 #endif
774 
775 #endif
o2scl::eos_had_skyrme::alt_params_set
void alt_params_set(double Crr00, double Crr10, double Crr0D, double Crr1D, double Crt0, double Crt1, double CrDr0, double CrDr1, double CrnJ0, double CrnJ1, double alpha2)
Set using alternate parameterization.
thermo_tl< double >::ed
double ed
o2scl::eos_had_skyrme::parent_method
bool parent_method
Use eos_had_base methods for saturation properties.
Definition: eos_had_skyrme.h:491
boost::numeric::ublas::vector< double >
o2scl::eos_had_skyrme::type
virtual const char * type()
Return string denoting type ("eos_had_skyrme")
Definition: eos_had_skyrme.h:459
O2SCL_ERR2
#define O2SCL_ERR2(d, d2, n)
o2scl::eos_had_temp_eden_base
A hadronic EOS at finite temperature based on a function of the densities [abstract base].
Definition: eos_had_base.h:1186
o2scl::eos_had_skyrme::b4p
double b4p
Isovector spin-orbit term (in )
Definition: eos_had_skyrme.h:480
o2scl::eos_had_skyrme::fmsom
virtual double fmsom(double nb)
Calculate effective mass.
thermo_tl< double >::en
double en
thermo_tl< double >
o2scl::eos_had_skyrme::~eos_had_skyrme
virtual ~eos_had_skyrme()
Destructor.
Definition: eos_had_skyrme.h:438
o2scl::eos_had_base::n0
double n0
Saturation density in .
Definition: eos_had_base.h:346
o2scl::eos_had_skyrme::check_input
void check_input(fermion_t &ne, fermion_t &pr, double T)
Definition: eos_had_skyrme.h:387
o2scl::eos_had_skyrme::calc_deriv_temp_e
virtual int calc_deriv_temp_e(fermion_deriv &ne, fermion_deriv &pr, double temper, thermo &th, thermo_np_deriv_helm &thd)
Equation of state as a function of the densities at finite temperature and including second derivativ...
o2scl::eos_had_skyrme::eff_mass
void eff_mass(fermion_t &ne, fermion_t &pr, double &term, double &term2)
Evaluate the effective masses for neutrons and protons.
Definition: eos_had_skyrme.h:599
o2scl::eos_had_skyrme::alt_params_get
void alt_params_get(double &Crr00, double &Crr10, double &Crr0D, double &Crr1D, double &Crt0, double &Crt1, double &CrDr0, double &CrDr1, double &CrnJ0, double &CrnJ1, double &alpha2)
Get alternate parameterization.
o2scl::eos_had_skyrme::hamiltonian_coeffs
void hamiltonian_coeffs(double &ham1, double &ham2, double &ham3, double &ham4, double &ham5, double &ham6)
Compute the coefficients of the potential energy which have unique dependence on the densities.
o2scl::eos_had_skyrme::calpar
int calpar(double gt0=-10.0, double gt3=70.0, double galpha=0.2, double gt1=2.0, double gt2=-1.0)
Calculate and from the saturation properties.
o2scl::eos_had_skyrme::b4
double b4
Isoscalar spin-orbit term (in )
Definition: eos_had_skyrme.h:477
o2scl::eos_had_skyrme::calc_temp_e
virtual int calc_temp_e(fermion &ne, fermion &pr, double temper, thermo &th)
Equation of state as a function of densities at finite temperature.
exc_einval
exc_einval
o2scl::eos_had_skyrme::nrf
fermion_nonrel nrf
Thermodynamics of non-relativistic fermions.
Definition: eos_had_skyrme.h:750
thermo_tl< double >::pr
double pr
o2scl::eos_had_skyrme::landau_neutron
void landau_neutron(double n0, double m, double &f0, double &g0, double &f1, double &g1)
Calculate the Landau parameters for neutron matter.
o2scl::eos_had_skyrme::nrfd
fermion_deriv_nr nrfd
Thermodynamics of non-relativistic fermions with derivatives.
Definition: eos_had_skyrme.h:753
o2scl::eos_had_skyrme
Skyrme hadronic equation of state.
Definition: eos_had_skyrme.h:226
o2scl::eos_had_skyrme::landau_nuclear
void landau_nuclear(double n0, double m, double &f0, double &g0, double &f0p, double &g0p, double &f1, double &g1, double &f1p, double &g1p)
Calculate the Landau parameters for nuclear matter.
o2scl::eos_had_skyrme::fkprime_nuc
virtual double fkprime_nuc(double nb)
Skewness in nuclear (isospin-symetric) matter.
o2scl::eos_had_skyrme::eos_had_skyrme
eos_had_skyrme()
Create a blank Skyrme EOS.
O2SCL_ERR
#define O2SCL_ERR(d, n)
o2scl::eos_had_skyrme::second_deriv
void second_deriv(fermion_t &ne, fermion_t &pr, double ltemper, thermo &locth, thermo_np_deriv_helm &locthd, double term, double term2, double ham1, double ham2, double ham3, double ham4, double ham5, double ham6)
Compute second derivatives of the free energy.
Definition: eos_had_skyrme.h:297
o2scl::eos_had_skyrme::feoa
virtual double feoa(double nb)
Calculate binding energy.
o2scl::eos_had_skyrme::reference
std::string reference
Bibliographic reference.
Definition: eos_had_skyrme.h:483
o2scl::eos_had_skyrme::base_thermo
void base_thermo(fermion_t &ne, fermion_t &pr, double ltemper, thermo &locth, double term, double term2, double ham1, double ham2, double ham3, double ham4, double ham5, double ham6)
Compute the base thermodynamic quantities.
Definition: eos_had_skyrme.h:246
o2scl::eos_had_skyrme::fcomp_nuc
virtual double fcomp_nuc(double nb)
Calculate compressibility in nuclear (isospin-symmetric matter)
o2scl::eos_had_skyrme::W0
double W0
Spin-orbit splitting (in )
Definition: eos_had_skyrme.h:474
o2scl::eos_had_skyrme::calc_e
virtual int calc_e(fermion &ne, fermion &pr, thermo &lt)
Equation of state as a function of densities at zero temperature.
o2scl::eos_had_skyrme::alt_params_saturation
void alt_params_saturation(double n0, double EoA, double K, double Ms_star, double a, double L, double Mv_star, double CrDr0, double CrDr1, double CrnJ0, double CrnJ1)
Use the specified saturation properties and couplings and the function alt_params_set() to set the Sk...
o2scl::eos_had_skyrme::check_landau
int check_landau(double nb, double m)
Check the Landau parameters for instabilities.
o2scl::eos_had_skyrme::fesym
virtual double fesym(double nb, double alpha=0.0)
Calculate symmetry energy.

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