nucmass_ldrop.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 nucmass_ldrop.h
24  \brief File defining \ref o2scl::nucmass_ldrop
25 */
26 #ifndef LDROP_MASS_H
27 #define LDROP_MASS_H
28 
29 #include <cmath>
30 #include <string>
31 #include <map>
32 #include <o2scl/nucleus.h>
33 #include <o2scl/nucmass.h>
34 #include <o2scl/constants.h>
35 #include <o2scl/eos_had_base.h>
36 #include <o2scl/eos_had_rmf.h>
37 #include <o2scl/fermion_eff.h>
38 #include <o2scl/inte_qagiu_gsl.h>
39 
40 #ifndef DOXYGEN_NO_O2NS
41 namespace o2scl {
42 #endif
43 
44  /** \brief Simple liquid drop mass formula
45 
46  Includes bulk \part plus surface and Coulomb (no pairing)
47  without neutron skin and without any isospin contribution
48  to the surface energy.
49 
50  The NL4 EOS is loaded by default.
51 
52  \warning This class sets part::inc_rest_mass to true
53  for the particle objects specified in set_n_and_p().
54 
55  \hline
56 
57  <b>Central densities</b>
58 
59  Given a saturation density, \f$ n_0 \f$ and a transition
60  density, \f$ n_t \f$, we set \f$ I = 1 - 2 Z/A \f$, and then
61  assume \f$ \delta = I \f$. We further assume that the
62  isospin-asymmetric saturation density \f$ n_L \f$ is
63  \f[
64  n_L = n_0 + n_1 \delta^2
65  \f]
66  and then we can compute \f$ n_{p} = (1 - \delta)/2 n_L \f$ and
67  \f$ n_{n} = (1 + \delta)/2 n_L \f$ .
68 
69  Note that \f$ \delta = I \f$ implies no neutron skin. A neutron
70  skin occurs when \f$ \delta < I \f$, and \f$ \delta = 0 \f$
71  implies a "maximum skin size" which is occurs when no extra
72  neutrons are in center and all extra neutrons are located in the
73  skin, i.e. \f$ N_{\mathrm{skin}} = N-Z \f$.
74 
75  <b>Nuclear radii</b>
76 
77  The neutron and proton radii are determined from the
78  central densities with
79  \f{eqnarray*}
80  R_n &=& \left( \frac{3 N}{4 \pi n_n} \right)^{1/3} \nonumber \\
81  R_n &=& \left( \frac{3 Z}{4 \pi n_p} \right)^{1/3}
82  \f}
83 
84  <b>Bulk energy contribution</b>
85 
86  The bulk binding energy contribution ( \f$ \sim -16 \f$
87  MeV per nucleon) and the symmetry energy are computing using the
88  hadronic EOS (either \ref def_had_eos or the EOS specified in
89  the most recent call to set_eos_had_temp_base() ). The bulk
90  energy per baryon is
91  \f[
92  E_{\mathrm{bulk}}/A = \frac{\hbar c}{n_{L} }
93  \left[\varepsilon(n_n,n_p) - n_n m_n - n_p m_p \right]
94  \f]
95 
96  <b>Surface energy contribution</b>
97 
98  The surface energy density is (\ref Ravenhall83)
99  \f[
100  \varepsilon = \frac{\chi d \sigma}{R}
101  \f]
102  where \f$ \sigma \f$ is the surface tension. The factor \f$ \chi
103  \f$ is typically taken care of by the caller, so we ignore it
104  for now. To compute the surface energy per baryon, we divide by
105  the baryon density, \f$ n_n + n_p \f$. We can rewrite this
106  \f[
107  E_{\mathrm{surf}} = \frac{3 \sigma}{n_n + n_p}
108  \left[ \frac{3 A}{ 4 (n_n+n_p) \pi}
109  \right]^{-1/3}
110  \f]
111  or
112  \f[
113  E_{\mathrm{surf}} = \frac{\sigma}{n_L}
114  \left(\frac{36 \pi n_L}{A} \right)^{1/3}
115  \f]
116  where the surface tension \f$ \sigma \f$ (in MeV) is given in
117  \ref surften.
118 
119  Taking a typical value, \f$ \sigma =1.1~\mathrm{MeV} \f$ and
120  \f$ n_L = 0.16~\mathrm{fm}^{-3} \f$, gives the standard result,
121  \f$ E_{\mathrm{surf}}/A = 18~\mathrm{MeV}~A^{-1/3} \f$.
122 
123  <b>Coulomb energy contribution</b>
124 
125  The Coulomb energy density (\ref Ravenhall83) is
126  \f[
127  \varepsilon_{\mathrm{Coul}} = \frac{4 \pi}{5} n_p^2 e^2 R_p^2
128  \f]
129  The energy per baryon is
130  \f[
131  E_{\mathrm{Coul}}/A = \frac{4 \pi}{5 n_L} n_p^2 e^2 R_p^2
132  \f]
133  This is the expression used in the code, except for a prefactor
134  \ref coul_coeff which is a fit parameter and should be close to
135  unity.
136 
137  Assuming \f$ R_p = R \f$
138  and \f$ Z = \frac{4 \pi}{3} R^3 n_p \f$
139  and \f$ R = \left[ 3 A / (4 \pi n_L) \right]^{1/3} \f$
140  gives
141  \f[
142  E_{\mathrm{Coul}}/A = \frac{6^{2/3}}{5}
143  \pi^{1/3} e^2 n_L^{1/3} \frac{Z^2}{A^{4/3}}
144  \f]
145  and taking \f$ n_L = 0.16~\mathrm{fm}^{-3} \f$ and
146  \f$ e^2 = \hbar c/137 \f$ gives the standard result
147  \f[
148  E_{\mathrm{Coul}}/A = 0.76~\mathrm{MeV}~Z^2 A^{-4/3}
149  \f]
150 
151  \todo 12/4/14: This doesn't gracefully handle negative values of
152  n0 and n1 as then the neutron and proton densities become
153  negative. This needs to be addressed. For now, there is a
154  fix at line 246 in nucmass_ldrop.cpp .
155 
156  \hline
157 
158  <b>References</b>
159 
160  Designed for \ref Steiner08 based on \ref Lattimer85 and
161  \ref Lattimer91 .
162 
163  \hline
164  */
166 
167  public:
168 
169  nucmass_ldrop();
170 
171  /// \name Input parameters
172  //@{
173  /// Density asymmetry (default 0)
174  double n1;
175 
176  /** \brief Saturation density ( The default is \f$ 0.16
177  \mathrm{fm}^{-3} \f$)
178  */
179  double n0;
180 
181  /// Surface tension in MeV (default 1.1 MeV)
182  double surften;
183 
184  /// Coulomb coefficient (default 1.0)
185  double coul_coeff;
186  //@}
187 
188  /// \name Output quantities
189  //@{
190  /// Internal average neutron density
191  double nn;
192 
193  /// Internal average proton density
194  double np;
195 
196  /// Neutron radius
197  double Rn;
198 
199  /// Proton radius
200  double Rp;
201 
202  /// Bulk \part of energy
203  double bulk;
204 
205  /// Surface \part of energy
206  double surf;
207 
208  /// Coulomb \part of energy
209  double coul;
210  //@}
211 
212  /** \brief Given \c Z and \c N, return the mass excess in MeV
213 
214  \comment
215  We don't need to implement mass_excess() for integers because
216  that's done in the parent nucmass_cont.
217  \endcomment
218  */
219  virtual double mass_excess_d(double Z, double N);
220 
221  /// Given \c Z and \c N, return the mass excess in MeV
222  virtual double mass_excess(int Z, int N) {
223  return mass_excess_d(Z,N);
224  }
225 
226  /** \brief Given \c Z and \c N, return the binding energy in MeV
227 
228  This function is currently independent of \c npout, \c nnout,
229  and \c chi.
230  */
231  virtual double drip_binding_energy_d(double Z, double N,
232  double npout, double nnout,
233  double chi, double T);
234 
235  /// \name EOS and particle parameters
236  //@{
237  /// Change the base hadronic EOS
239  heos=&uhe;
240  return 0;
241  }
242 
243  /// The default hadronic EOS
245 
246  /// Change neutron and proton objects
247  void set_n_and_p(fermion &un, fermion &up) {
248  n=&un;
249  p=&up;
250  return;
251  }
252 
253  /// Default neutron
255 
256  /// Default proton
258  //@}
259 
260  /// \name Fitting functions
261  //@{
262  /// Fix parameters from an array for fitting
263  virtual int fit_fun(size_t nv, const ubvector &x);
264 
265  /// Fill array with guess from present values for fitting
266  virtual int guess_fun(size_t nv, ubvector &x);
267  //@}
268 
269  /// Return the type, \c "nucmass_ldrop".
270  virtual const char *type() { return "nucmass_ldrop"; }
271 
272  /// Energy and pressure
274 
275 #ifndef DOXYGEN_INTERNAL
276 
277  protected:
278 
279  /// Pointer to neutron
281  /// Pointer to proton
283  /// The base EOS for bulk matter
285 
286 #endif
287 
288  };
289 
290  /** \brief More advanced liquid drop model
291 
292  In addition to the physics in \ref nucmass_ldrop, this includes
293  corrections for
294  - finite temperature
295  - neutron skin
296  - an isospin-dependent surface energy
297  - decrease in the Coulomb energy from external protons
298 
299  \note The input parameter T should be given in units of inverse
300  Fermis -- this is a bit unusual since the binding energy is
301  returned in MeV, but we keep it for now.
302 
303  <b>Bulk energy</b>
304 
305  The central densities and radii, \f$ n_n, n_p, R_n, R_p \f$
306  are all determined in the same way as \ref nucmass_ldrop,
307  except that now \f$ \delta \equiv I \zeta \f$, where
308  \f$ \zeta \f$ is stored in \ref doi . Note that this
309  means \f$ N > Z~\mathrm{iff}~R_n>R_p \f$.
310 
311  If \ref new_skin_mode is false, then the bulk energy is
312  also computed as in \ref nucmass_ldrop. Otherwise, the
313  number of nucleons in the core is computed with
314  \f{eqnarray*}
315  A_{\mathrm{core}} = Z (n_n+n_p)/n_p~\mathrm{for}~N\geq Z \\
316  A_{\mathrm{core}} = N (n_n+n_p)/n_p~\mathrm{for}~Z>N \\
317  \f}
318  and \f$ A_{\mathrm{skin}} = A - A_{\mathrm{core}} \f$.
319  The core contribution to the bulk energy is
320  \f[
321  E_{\mathrm{core}}/A = \left(\frac{A_{\mathrm{core}}}{A}\right)
322  \frac{\hbar c}{n_{L} }
323  \left[\varepsilon(n_n,n_p) - n_n m_n - n_p m_p \right]
324  \f]
325  then the skin contribution is
326  \f[
327  E_{\mathrm{skin}}/A = \left(\frac{A_{\mathrm{skin}}}{A}\right)
328  \frac{\hbar c}{n_{L} }
329  \left[\varepsilon(n_n,0) - n_n m_n \right]~\mathrm{for}~N>Z
330  \f]
331  and
332  \f[
333  E_{\mathrm{skin}}/A = \left(\frac{A_{\mathrm{skin}}}{A}\right)
334  \frac{\hbar c}{n_{L} }
335  \left[\varepsilon(0,n_p) - n_p m_p \right]~\mathrm{for}~Z>N
336  \f]
337 
338  <b>Surface energy</b>
339 
340  If \ref full_surface is false, then the surface energy is
341  just that from \ref nucmass_ldrop , with an extra factor
342  for the surface symmetry energy
343  \f[
344  E_{\mathrm{surf}} = \frac{\sigma}{n_L}
345  \left(\frac{36 \pi n_L}{A} \right)^{1/3}
346  \left( 1- \sigma_{\delta} \delta^2 \right)
347  \f]
348  where \f$ \sigma_{\delta} \f$ is unitless and stored in \ref ss.
349 
350  If \ref full_surface is true, then the surface energy is modified
351  by a cubic dependence for the medium and contains finite temperature
352  corrections.
353 
354  <b>Coulomb energy</b>
355 
356  The Coulomb energy density (\ref Ravenhall83) is
357  \f[
358  \varepsilon = 2 \pi e^2 R_p^2 n_p^2 f_d(\chi_p)
359  \f]
360  where the function \f$ f_d(\chi) \f$ is
361  \f[
362  f_d(\chi_p) = \frac{1}{(d+2)}
363  \left[ \frac{2}{(d-2)} \left( 1 - \frac{d}{2}
364  \chi_p^{(1-2/d)} \right) + \chi_p \right]
365  \f]
366 
367  This class takes \f$ d=3 \f$ .
368 
369  <b>Todos and Future</b>
370 
371  \todo This is based on LPRL, but it's a little different in
372  Lattimer and Swesty. I should document what the difference is.
373 
374  \todo The testing could be updated.
375 
376  \future Add translational energy?
377 
378  \future Remove excluded volume correction and compute nuclear
379  mass relative to the gas rather than relative to the vacuum.
380 
381  \future In principle, Tc should be self-consistently determined
382  from the EOS.
383 
384  \future Does this work if the nucleus is "inside-out"?
385 
386  \comment
387  \todo The choice of nn and np from n0 and n1 is very closely
388  related to FRDM (\ref Moller95). We should document this here.
389 
390  I've taken this out, because it appears to me that Moller '95
391  actually set this term (n1 = -3 L/K) to zero.
392  \endcomment
393 
394  \comment
395 
396  \hline
397 
398  Excluded volume and \ref rel_vacuum:
399 
400  Typically in a single-nucleus EOS with a neutron drip
401  (ignoring translational degrees of freedom for the nucleus)
402  \f[
403  f = n_N m_N + (1-\chi_n) f_{n,\mathrm{drip}}
404  \f]
405  where
406  \f[
407  m_N = \frac{A}{n_n+n_p}(f - n_n m_n - n_p m_p)
408  \f]
409  Since \f$ n_N = 3/(4 \pi R_{\mathrm{ws}}^3) \f$, and
410  \f$ \chi_n = (R_n/R_{\mathrm{ws}})^3 \f$, this is
411  \f[
412  f = \frac{3}{4 \pi R_{\mathrm{ws}}^3}
413  \left[ m_N - f_{n,\mathrm{drip}} \frac{4 \pi}{3} R_n^3 \right]
414  + f_{n,\mathrm{drip}}
415  \f]
416 
417  \endcomment
418 
419  \hline
420 
421  <b>References</b>
422 
423  Designed in \ref Steiner08 and \ref Souza09 based in part
424  on \ref Lattimer85 and \ref Lattimer91 .
425 
426  \hline
427  */
429 
430  public:
431 
433 
434  /// Return the type, \c "nucmass_ldrop_skin".
435  virtual const char *type() { return "nucmass_ldrop_skin"; }
436 
437  /// Fix parameters from an array for fitting
438  virtual int fit_fun(size_t nv, const ubvector &x);
439 
440  /// Fill array with guess from present values for fitting
441  virtual int guess_fun(size_t nv, ubvector &x);
442 
443  /** \brief If true, properly fix the surface for the pure neutron
444  matter limit (default true)
445  */
447 
448  /** \brief If true, separately compute the skin for the bulk energy
449  (default false)
450  */
452 
453  /// Ratio of \f$ \delta/I \f$ (default 0.8).
454  double doi;
455 
456  /// Surface symmetry energy (default 0.5)
457  double ss;
458 
459  /// \name Input parameters for temperature dependence
460  //@{
461  /// Exponent (default 1.25)
462  double pp;
463 
464  /// Coefficient (default 0.935)
465  double a0;
466 
467  /// Coefficient (default -5.1)
468  double a2;
469 
470  /// Coefficient (default -1.1)
471  double a4;
472  //@}
473 
474  /** \brief If true, define the nuclear mass relative to the vacuum
475  (default true)
476  */
478 
479  /** \brief The critical temperature of isospin-symmetric matter in
480  \f$ fm^{-1} \f$ (default \f$ 20.085/(\hbar c)\f$.)
481  */
482  double Tchalf;
483 
484  /** \brief Return the free binding energy of a \nucleus in a many-body
485  environment
486  */
487  virtual double drip_binding_energy_d(double Z, double N,
488  double npout, double nnout,
489  double chi, double T);
490  };
491 
492  /** \brief Liquid drop model with pairing
493 
494  This class adds a pairing correction
495  \f[
496  E_{\mathrm{pair}}/A = - \frac{\zeta}{2 A^{3/2}}
497  \left[
498  \cos(Z \pi) + \cos(N \pi)
499  \right]
500  \f]
501  where \f$ \zeta \f$ is stored in \ref Epair. The trigonometric
502  functions give the expected result for integer values of
503  N and Z.
504  */
506 
507  public:
508 
509  /// Return the type, \c "nucmass_ldrop_pair".
510  virtual const char *type() { return "nucmass_ldrop_pair"; }
511 
513  nfit=7;
514  Epair=13.0;
515  }
516 
517  /// Fix parameters from an array for fitting
518  virtual int fit_fun(size_t nv, const ubvector &x);
519 
520  /// Fill array with guess from present values for fitting
521  virtual int guess_fun(size_t nv, ubvector &x);
522 
523  /// Pairing energy coefficient (default 13 MeV)
524  double Epair;
525 
526  /// Most recently computed pairing energy per baryon
527  double pair;
528 
529  /** \brief Return the free binding energy of a \nucleus in a many-body
530  environment
531  */
532  virtual double drip_binding_energy_d
533  (double Z, double N, double npout, double nnout,
534  double chi, double T);
535 
536  };
537 
538 #ifndef DOXYGEN_NO_O2NS
539 }
540 #endif
541 
542 #endif
eos_had_rmf def_had_eos
The default hadronic EOS.
double surf
Surface part part of energy.
double pp
Exponent (default 1.25)
fermion def_neutron
Default neutron.
More advanced liquid drop model.
bool rel_vacuum
If true, define the nuclear mass relative to the vacuum (default true)
virtual const char * type()
Return the type, "nucmass_ldrop_skin".
double doi
Ratio of (default 0.8).
bool full_surface
If true, properly fix the surface for the pure neutron matter limit (default true) ...
double pair
Most recently computed pairing energy per baryon.
double n0
Saturation density ( The default is )
void set_n_and_p(fermion &un, fermion &up)
Change neutron and proton objects.
Liquid drop model with pairing.
virtual const char * type()
Return the type, "nucmass_ldrop_pair".
double Epair
Pairing energy coefficient (default 13 MeV)
fermion * n
Pointer to neutron.
A finite temperature hadronic EOS [abstract base].
Definition: eos_had_base.h:962
double coul
Coulomb part part of energy.
double ss
Surface symmetry energy (default 0.5)
bool new_skin_mode
If true, separately compute the skin for the bulk energy (default false)
Relativistic mean field theory EOS.
Definition: eos_had_rmf.h:295
double n1
Density asymmetry (default 0)
virtual double mass_excess_d(double Z, double N)
Given Z and N, return the mass excess in MeV.
double surften
Surface tension in MeV (default 1.1 MeV)
Simple liquid drop mass formula.
eos_had_temp_base * heos
The base EOS for bulk matter.
double bulk
Bulk part part of energy.
fermion * p
Pointer to proton.
double coul_coeff
Coulomb coefficient (default 1.0)
virtual int fit_fun(size_t nv, const ubvector &x)
Fix parameters from an array for fitting.
fermion def_proton
Default proton.
double a2
Coefficient (default -5.1)
double Rp
Proton radius.
double a0
Coefficient (default 0.935)
double Rn
Neutron radius.
thermo th
Energy and pressure.
virtual double mass_excess(int Z, int N)
Given Z and N, return the mass excess in MeV.
virtual int guess_fun(size_t nv, ubvector &x)
Fill array with guess from present values for fitting.
double np
Internal average proton density.
int set_eos_had_temp_base(eos_had_temp_base &uhe)
Change the base hadronic EOS.
double a4
Coefficient (default -1.1)
virtual double drip_binding_energy_d(double Z, double N, double npout, double nnout, double chi, double T)
Given Z and N, return the binding energy in MeV.
double nn
Internal average neutron density.
double Tchalf
The critical temperature of isospin-symmetric matter in (default .)
virtual const char * type()
Return the type, "nucmass_ldrop".

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