Anomalous quartic term in the expansion of the symmetry energy

TL;DR

This paper investigates the quartic term in the symmetry energy within relativistic mean field theory, revealing it can become very large and questioning the validity of the standard even power expansion due to potential non-analytic contributions.

Contribution

It introduces the possibility of significant non-analytic contributions to the symmetry energy, challenging the conventional polynomial expansion approach.

Findings

Quartic term can reach very large values in the model.

Standard even power expansion may be inadequate.

Non-analytic contributions might be important.

Abstract

The quartic term in the framework of relativistic mean field theory with inclusion of scalar meson interactions is investigated. It is shown that the quartic term in the asymmetric expansion of nuclear matter energy may reach very large values. This makes the even power expansion of asymmetry questionable and suggests possible non-analytic contributions to the energy of matter.

Tables

Table 1. Table 1: Model parameters for isoscalar and isovector sectors.

isoscalar sector	$C_{\sigma }^2=11\left[{\mathrm{fm}}^2\right]$
	$C_{\omega }^2=6.48\left[{\mathrm{fm}}^2\right]$
	$b=0.054$
	$c=-0.0057$
isovector sector	standard RMF theory
	$C_{\delta }^2\left[{\mathrm{fm}}^2\right]$	$C_{\rho }^2\left[{\mathrm{fm}}^2\right]$	$L\left[\mathrm{MeV}\right]$
	1.0	7.1	91.8
	2.5	12.3	103.7
	3.25	14.8	110.1
	3.5	15.7	112.3
	linear $\sigma$-$\delta$ meson interaction ($\alpha =1$)
	$g_{\alpha =1}=-0.009{\mathrm{fm}}^{-1}$
	$C_{\delta }^2\left[{\mathrm{fm}}^2\right]$	$C_{\rho }^2\left[{\mathrm{fm}}^2\right]$	$L\left[\mathrm{MeV}\right]$
	1.0	7.4	88.8
	2.5	14.8	77.6
	3.25	19.6	57.2
	3.5	21.3	46.8
	quadratic $\sigma$-$\delta$ meson interaction ($\alpha =2$)
	$g_{\alpha =2}=-0.004$
	$C_{\delta }^2\left[{\mathrm{fm}}^2\right]$	$C_{\rho }^2\left[{\mathrm{fm}}^2\right]$	$L\left[\mathrm{MeV}\right]$
	1.0	7.3	88.3
	2.5	13.6	77.6
	3.25	17.1	62.3
	3.5	18.4	55.2

Equations

$$\begin{array}{rl}\epsilon (n,x)=& \epsilon (n,1/2)+S_2(n)(1-2x{)}^2+\\ & S_4(n)(1-2x{)}^4+\mathcal{O}((1-2x{)}^6).\end{array}$$

$$S_{2i}(n)={\frac{1}{2^{2i}\cdot (2i)!}\frac{{\partial }^{2i}\left(\epsilon \left(n,x\right)/n\right)}{{\partial x}^{2i}}|}_{x=\frac{1}{2}}.$$

$${\epsilon }_{kin}=({\int }_0^{k_p}+{\int }_0^{k_n})\sqrt{k^2+m^2}\frac{k^{2}dk}{{\pi }^2},$$

$$S_2^{kin}=\frac{k^2}{6\sqrt{k^2+m^2}}$$

$$S_4^{kin}=\frac{k^2\left(10k^4+11k^2{m}^2+4m^4\right)}{648{\left(k^2+m^2\right)}^{5/2}}$$

$$S_2^{kin}=\frac{1}{3}{E}_0^{(nr)}$$

$$S_4^{kin}=\frac{1}{81}{E}_0^{(nr)}$$

$$\mathcal{L}={\mathcal{L}}_0+{\mathcal{L}}_{int}+{\mathcal{L}}_{\sigma \delta }.$$

$${\mathcal{L}}_{int}=g_{\sigma }\sigma \bar{\psi }\psi -g_{\omega }{\omega }_{\mu }\bar{\psi }{\gamma }^{\mu }\psi -\frac{1}{2}{g}_{\rho }{\vec{\rho }}_{\mu }\bar{\psi }{\gamma }^{\mu }\vec{\tau }\psi +g_{\delta }\vec{\delta }\bar{\psi }\vec{\tau }\psi -U(\sigma ),$$

$${\mathcal{L}}_{\sigma \delta }={\tilde{g}}_{\alpha }{\sigma }^{\alpha }{\vec{\delta }}^2,\alpha =1,2.$$

$$m_{\sigma }^2\sigma =g_{\sigma }\left(n_p^s+n_n^s\right)-\frac{\partial U}{\partial \sigma }-{\tilde{g}}_{\alpha }\alpha {\sigma }^{\alpha -1}{\delta }^2,$$

$$m_{\delta }^2{\delta }_3=g_{\delta }\left(n_p^s-n_n^s\right)-2{\tilde{g}}_{\alpha }{\sigma }^{\alpha }\delta .$$

$$n_i=\frac{2}{{(2\pi )}^3}{\int }_0^{k_i}{d}^{3}k,i=p,n$$

$$n_i^s=\frac{2}{{(2\pi )}^3}{\int }_0^{k_i}\frac{m_i^{\ast }}{\sqrt{k^2+m_i^{\ast 2}}}{d}^{3}k,i=p,n$$

$$\left\{\begin{array}[]{ll}m_{p}^{\ast}=m-g_{\sigma}{\sigma}-g_{\delta}{{\delta}},\\ m_{n}^{\ast}=m-g_{\sigma}{\sigma}+g_{\delta}{{\delta}},\end{array}\right.$$

$$\begin{array}{rl}\epsilon =& \frac{1}{2}{m}_{\sigma }^2{\sigma }^2+\frac{1}{2}{m}_{\omega }^2{\omega }^2+\frac{1}{2}{m}_{\rho }^2{\rho }^2+\frac{1}{2}{m}_{\delta }^2{\delta }^2+U(\sigma )\\ & +g_{\alpha }{\sigma }^{\alpha }{\delta }^2+\sum _{i=n,p}{\int }_0^{k_i}{k}^2\sqrt{k^2+m_i^{\ast 2}}\frac{dk}{{\pi }^2}.\end{array}$$

$$S_4(n)={\frac{1}{384}\frac{{\partial }^4\left(\epsilon \left(n,x\right)/n\right)}{{\partial x}^4}|}_{x=\frac{1}{2}}.$$

$$S_2(n)=\frac{k^2}{6E}+\frac{1}{8}n{C}_{\rho }^2-C_{\delta }^2\frac{k^3{m}^{\ast 2}}{3E^2{\pi }^2{f}_{\delta }},$$

$$f_{\delta }=1+C_{\delta }^2\mathcal{A}+8C_{\delta }^2{\sigma }^{\alpha }{g}_{\alpha },$$

$$\mathcal{A}(k,m^{\ast })=\frac{4}{(2\pi {)}^3}{\int }_0^k\frac{k^2{d}^{3}k}{(k^2+m^{\ast 2}{)}^{3/2}}=\frac{3n_s}{m^{\ast }}-\frac{3n}{E}.$$

$$\begin{array}{rl}{S}_4& (n)=\frac{k^2\left(10k^4+11k^2{m}^{\ast 2}+4m^{\ast 4}\right)}{648E^5}-g_{\alpha }^2\frac{8{\alpha }^2{m}^{\ast 4}{n}^3{\sigma }^{2(\alpha -1)}{C}_{\delta }^8}{E^4{f}_{\delta }^4{f}_{\sigma }}-\\ & g_{\alpha }\frac{2\alpha m^{\ast 2}{n}^2{\sigma }^{\alpha -1}{C}_{\delta }^4\left(6E{k}^2{m}^{\ast }n{C}_{\delta }^2{f}_{\delta }+9n{C}_{\delta }^4\left(3k^2{m}^{\ast }n+2m^{\ast 3}n-2E^3{n}_s\right)+E^2{k}^2{m}^{\ast }{f}_{\delta }^2\right)}{3E^7{f}_{\delta }^4{f}_{\sigma }}-\\ & \frac{n}{216E^{11}{k}^4{f}_{\delta }^4{f}_{\sigma }}(-3E^5{k}^8{m}^{\ast 2}{f}_{\delta }^4-4E^4{k}^6{m}^{\ast 2}{C}_{\delta }^2{f}_{\delta }^3\left(E{f}_{\sigma }\left(4k^2+m^{\ast 2}\right)+9k^{2}n\right)-\\ & 18E^2{k}^4{m}^{\ast }n{C}_{\delta }^4{f}_{\delta }^2\left(m^{\ast }\left(k^2{f}_{\sigma }\left(k^2-2m^{\ast 2}\right)\left(k^2+m^{\ast 2}\right)+3En\left(3E^4+2k^4-4k^2{m}^{\ast 2}-3m^{\ast 4}\right)\right)-6E^4{k}^2{n}_s\right)+\\ & 108E^2{k}^4{m}^{\ast }{n}^2{C}_{\delta }^6{f}_{\delta }\left(k^2{m}^{\ast 3}\left(E{f}_{\sigma }+12n\right)+6E^3{k}^2{n}_s-9E^4{m}^{\ast }n+9m^{\ast 5}n\right)+27n^2{C}_{\delta }^8(-36E^7{k}^2{n}_s\left(k^2{n}_s-3E{m}^{\ast }n\right)+\\ & m^{\ast }\left(f_{\sigma }\left(7k^8{m}^{\ast 3}n+2E^7{k}^4{n}_s-3E^8{k}^2{m}^{\ast }n\right)-9E^4{m}^{\ast }n\left(4k^2{m}^{\ast }\left(\left(4k^2+3m^{\ast 2}\right)n_s-6E{m}^{\ast }n\right)+9E^{5}n\right)\right)\\ & {\left(m^{\ast }\right)}^{6}n\left(f_{\sigma }\left(14k^6+10k^4{m}^{\ast 2}+3k^2{m}^{\ast 4}\right)-9En\left(-18E^4+16k^4+24k^2{m}^{\ast 2}+9m^{\ast 4}\right)\right)))\end{array}$$

$$f_{\sigma }=\mathcal{A}+\frac{1}{C_{\sigma }^2}+U^{\prime \prime }(\sigma ),$$

Full text

I Introduction

The fundamental quantity for describing neutron star matter is the energy density $\varepsilon(n,x)$ which is a function of the baryon number $n$ and proton fraction $x=\frac{n_{p}}{n}$. The proton fraction is the number of protons per baryon. The dependence of the energy of nuclear matter on the isospin asymmetry, represented by $x$, can be described by the quantity called the symmetry energy. It plays an important role in many aspects of nuclear physics and astrophysics, such as the mass-radius relation and the composition of the crust. There are many scientific papers based on both theoretical calculations and experimental measurements, which focus on the energy symmetry behavior and its dependence on the density. Although it is relatively well determined at the nuclear saturation density $n_{0}=0.16\rm~{}fm^{-3}$, little is known about its behavior at densities much higher than $n_{0}$. Moreover, the asymmetry dependence of $\varepsilon$ is also poorly known. It can be quantitatively measured by expansion in powers of $1-2x$. The isospin symmetry ensures that only even powers are present:

[TABLE]

As per Taylor’s theorem, the coefficients of Eq. (1) for order $2i$ take the following forms:

[TABLE]

The first term in Eq. (1) describes symmetric nuclear matter, which can be determined relatively well around $n_{0}$ using the compressibility coefficient. The rest of this equation deals with asymmetric nuclear matter (ANM). Although details of ANM are still far from certain, great progress has occurred in recent years, particularly in the determination of the quadratic term $S_{2}$. The behavior of $S_{2}$, as a function of the density, is better understood, especially at low energy densities, due to the determination of the symmetry energy slope at $n_{0}$ (for a review see Horowitz:2014bja ). The higher orders of this equation are still poorly known. There is widespread conviction that the next term in Eq. (1), called the quartic term is much smaller than the quadratic one and hence is usually neglected. Ignoring all terms of higher than the second order is called the parabolic approximation and is in common use. Interestingly, the quartic term may have some influence on neutron star structure and behavior, such as the density at which the core-crust transition appears or the threshold density for the URCA cooling process Steiner:2006bx ; Xu:2009vi ; Xu:2008vz ; Seif:2013tja .

The omission of the quartic term is usually justified by the comparison of quadratic and quartic terms coming from the kinetic energy contribution in the absence of any interactions

[TABLE]

then

[TABLE]

where $m$ is the nucleon mass and $k$ is the Fermi momentum for symmetric matter $k=(\frac{3\pi^{2}}{2}n)^{1/3}$. Their nonrelativistic versions take a simpler form,

[TABLE]

and here $E_{0}^{(nr)}$ is the nonrelativistic Fermi energy $E_{0}^{(nr)}=\frac{k^{2}}{2m}$. At the saturation point, the values of $S_{2}^{kin}(n_{0})$ and $S_{4}^{kin}(n_{0})$ are $12\rm~{}MeV$ and $0.5\rm~{}MeV$ respectively. However, recently presented considerations based on a high momentum tail in the distribution function Cai:2015mac indicate that the kinetic energy contribution could be much higher. From Eqs. (6) and (7) it can be seen that they scale as $n^{2/3}$, and at very high densities, typical for a neutron star center, $S_{2}^{kin}$ and $S_{4}^{kin}$ do not exceed $50\rm~{}MeV$ and $3\rm~{}MeV$, respectively. The inclusion of the effective nucleon mass $m^{*}$ (which decreases with the density, and so increases these kinetic energy terms slightly), does not significantly change the situation. Comparison of the kinetic contribution to $S_{2}$ with its experimental value at saturation point $S_{2}(n_{0})=31.2~{}\rm~{}MeV$ indicates a large role for interactions. Moreover, the interaction contribution usually has a stronger density dependence and becomes more significant as the density increases.

In the relativistic mean field (RMF) framework, the quartic term has been analyzed in Lee:1998zzd and Cai:2011zn and appears to be small, typically, one order of magnitude smaller than $S_{2}$. Other theoretical predictions for $S_{4}$ based on the Hartree-Fock theory Chen:2009wv ; Xu:2010kf ; Pu:2017kjx ; Liu:2018far lead to similar conclusions.

Some attempts were made to estimate the fourth order contribution to the symmetry energy at saturation density by fitting to the nuclear binding energy Jiang:2014jya ; Nandi:2016xnh ; Wang:2017itv or double beta decay Wan:2018lvn . Results for the quartic term were unexpectedly high, as large as 20 MeV. Such uncertainty concerning the value of the quartic term gained additional complexity from recent calculations based on the chiral effective field theory. In Ref. Kaiser:2015qia , it was shown that the chiral expansion introduces a logarithmic asymmetry dependence $(1-2x)^{4}\ln|1-2x|$ which makes the expansion given by Eq. (1) questionable.

In our work we have focused on the quartic term determination and its density dependence within the framework of the RMF model, enriched through scalar meson interactions, e.g., by inclusion of the $\sigma\!-\!\delta$ crossing term Zabari:2019ukk . We show that the expansion given by Eq. (1) must be treated with caution, because the quartic term appears to be unexpectedly large.

II model

The Lagrangian $\mathcal{L}$ in the RMF framework includes nucleons and the four meson fields $\sigma$, $\omega$, $\rho$, and $\delta$. It is given by

[TABLE]

Additionally, we have taken into consideration the interaction between the scalar mesons $\sigma$ and $\delta$ denoted in Eq.(8) as $L_{\sigma\delta}$. The first component $\mathcal{L}_{0}$ includes the standard components for non-interacting nucleons described by a bispinor field ${\psi_{p},\psi_{n}}$ along with mesons $\sigma$, $\omega$, $\rho$, and $\delta$. Expressions for those can be found in Zabari:2019ukk . The second part of Eq. (8), where the coupling constants appear, is given by

[TABLE]

where $U\left(\sigma\right)=\frac{1}{3}bm{(g_{\sigma}\sigma)}^{3}+\frac{1}{4}c{(g_{\sigma}\sigma)}^{4}$ is the self-interaction of the $\sigma$ meson. Finally, the third component, which represents a new interaction between two scalar mesons, has the following form

[TABLE]

We have considered two cases of this interaction: linear $\alpha=1$ and quadratic $\alpha=2$. Furthermore, this interaction has an influence on the field equations for the mean meson fields: $\sigma$ and $\delta_{3}$ (in the following we use the notation $\delta_{3}\equiv\delta$), coming from the minimum energy condition

[TABLE]

[TABLE]

The field equations for $\omega$ and $\rho$ remain unaffected, following the standard RMF model. Here, the following notation has been used:

[TABLE]

[TABLE]

[TABLE]

where $i=p,n$ stands for proton and neutron, $n=n_{p}+n_{n}$ is the baryon density, $n_{i}^{s}$ is the scalar density, $k_{i}=\left(3\ \pi^{2}n_{i}\right)^{\frac{1}{3}}$ denotes the Fermi momenta, and $m_{i}^{*}$ are the effective masses. Finally, the energy density has the form:

[TABLE]

All relevant parameters of the model are listed in Table 1. They are determined by the properties of matter at the saturation point. The constants in the isoscalar sector, $C_{\sigma},C_{\omega},b$, and $c$, are chosen to obtain the binding energy and compressibility of matter at $n_{0}$, $(\varepsilon/n-m)|_{n_{0}}=-16$ MeV, and $K=9n_{0}^{2}\left.\frac{\partial(\varepsilon/n)}{\partial n^{2}}\right|_{n_{0}}=230$ MeV, respectively.

In the isovector sector we fit the coupling constants $C_{\rho},C_{\delta},g_{\alpha}$ to two quantities given by the experimental data: the symmetry energy $S_{2}(n_{0})=30$ MeV and its slope $L=3n_{0}\left.\frac{\partial S_{2}}{\partial n}\right|_{n_{0}}$. For vanishing $g_{\alpha}$ the constants $C_{\rho},C_{\delta}$ are correlated by the value of $S_{2}(n_{0})=30$ MeV, as shown in Kubis:1997ew , and the slope $L$ gets a value that is too large. Recent measurements suggest that the slope value is 60 MeV with an uncertainty $\pm 30$ MeV Oertel:2016bki . Such low values of $L$ can be obtained in our model only by selecting a negative cross-term coupling $g_{\alpha}$. For linear and quadratic coupling we used the values $g_{1}=-0.009~{}\rm fm^{-1}$ and $g_{2}=-0.004$. The remaining two couplings $C_{\rho},C_{\delta}$ are uniquely determined by $S_{2}(n_{0})$ and $L$. The relations between measured quantities and coupling constants are listed in Table 1. A more detailed discussion of the model parameters is given in another work Zabari:2019ukk .

III Quartic term

The quartic term $S_{4}(n)$, which represents the fourth order of the Taylor expansion of Eq. (1), is defined as the fourth derivative of the energy density $\varepsilon(n,x)$ with respect to the proton fraction $x$:

[TABLE]

The computation can be made using the energy density $\varepsilon$ specified in Eq. (16) which is characterized by meson fields $\sigma$, $\omega$, $\rho$, and $\delta$ or, alternatively, it can be expressed by the effective masses of protons and neutrons given by Eq. (15). The energy density is minimized with respect to the meson fields, which corresponds to vanishing partial derivatives $\frac{\partial\varepsilon}{\partial\sigma}=0$ and $\frac{\partial\varepsilon}{\partial\delta}=0$ when $n$ and $x$ are kept constant, which leads to the equations of motion, (11) and (12). The minimum-energy condition means that the meson fields become functions of the densities $\sigma(n,x),\delta(n,x)$. The second derivative $\frac{\partial^{2}\varepsilon(n,x)}{\partial x^{2}}$, which corresponds to $S_{2}(n)$, holds the first derivative of the $\delta$ and $\sigma$ meson fields: $\frac{\partial\delta}{\partial x},\frac{\partial\sigma}{\partial x}$, while the fourth derivative, which corresponds to $S_{4}(n)$, includes field derivatives up to third order. In general, the $n$-th derivative of $\varepsilon$ requires the $(n-1)$-th derivative of the meson fields. For the isoscalar field $\sigma$, only even order derivatives remain at $x=1/2$ whereas for the isovector field $\delta$, only odd order derivatives remain. Summing up, for $S_{2}$ one needs to know $\frac{\partial\delta}{\partial x}$ and $S_{4}$ requires $\frac{\partial\delta}{\partial x},\frac{\partial^{3}\delta}{\partial x^{3}},\frac{\partial^{2}\sigma}{\partial x^{2}}$, which can be obtained by the differentiation of the equations of motion, Eqs. (11) and (12). In the following, we assume $k_{p}=k_{n}=k$ and $m_{p}^{\ast}=m_{n}^{\ast}=m^{\ast}$ for symmetric nuclear matter. Here, $k(n)={(\frac{3}{2}\pi^{2}n)}^{1/3}$. The symmetry energy has the form

[TABLE]

where $E=E(k,m^{*})=\sqrt{k^{2}+m^{\ast 2}}$ and $f_{\delta}$ is a function depending on the nucleon-$\delta$ coupling constant and the $\sigma$-$\delta$ meson interaction coupling constant. Its form is given as

[TABLE]

where the $\sigma\!-\!\delta$ coupling $\tilde{g}_{\alpha}$, appearing in Eqs. (11) and (12), is replaced here by a combination of coupling constants $g_{\alpha}=\displaystyle\frac{\tilde{g}_{\alpha}}{4g_{\sigma}^{\alpha}g_{\delta}^{2}}$, which allows one to simplify the following formulas. The quantity $\cal A$ is a function of scalar density $n_{s}=n_{p}^{s}+n_{n}^{s}$ and vector density $n=n_{p}+n_{n}$ at $x=1/2$

[TABLE]

The derivation of the quartic term is more laborious but attainable:

[TABLE]

Its complexity is associated with the higher derivatives of the meson fields $\delta$ and $\sigma$. One may note the appearance of $f_{\sigma}$

[TABLE]

which includes couplings from the isoscalar sector due to the appearance of the isoscalar field derivative as explained previously.

Here we show the results for the quartic term density dependence for different $C_{\delta}$ coupling. At the saturation point $S_{4}(n_{0})$, the values remain between 0.44 and 0.65 MeV for the linear model and between 0.49 and 0.61 MeV for the quadratic model. $C_{\delta}$ ranges within 0 and 3.5 $\rm fm^{2}$. These are typical values which one may encounter in RMF models. However, as the density increases the behavior of $S_{4}$ changes dramatically. As one can observe in Fig. 1, at densities of a few times $n_{0}$, for the quadratic model, $S_{4}$ is anomalously large, reaching over a thousand MeV. In the linear model the values of the quartic term at these densities oscillate between -100 and 100 MeV. In order to check whether such large values are the result of erroneous derivations, we numerically calculated the fourth derivative of the energy expressed by the polynomial-interpolated function up to sixth order. We obtained complete agreement of the numerical derivative with the analytical one given by Eq.(21). One must remember that those unexpectedly large values of the quartic term appear for the most appropriate values of $C_{\delta}^{2}=3.5~{}\rm fm^{2}$ which corresponds to the most likely symmetry energy slope $L$ around 50 MeV. When $S_{4}$ is much greater than the energy difference between symmetric and neutron matter, one may suspect that the expansion given by Eq. (1) is no longer valid. Indeed, in Fig. 2 the second and fourth order expansions of the energy for pure neutron matter are compared with the exact result. For the quadratic model especially, the subsequent correction in the energy expansion is greater than the difference between symmetric and pure neutron matter. It can be seen already that the second order expansion (with parabolic approximation) differs greatly from the exact result and adding the fourth order term does not improve it.

For further evidence we have analyzed the dependence of the energy on the proton fraction with a fixed density, shown in Fig. 3. At low densities the second and fourth order expansions work very well, they both overlap with the exact solution, independently of the coupling value $C_{\delta}$. However, as shown in Fig. 1, at higher densities, expansion (1) diverges significantly with increasing asymmetry.

IV Discussion

In this work the expansion of energy in powers of isospin asymmetry up to fourth order was studied for nuclear models with a $\sigma\mathchar 45\relax\delta$ interaction term, Eq.(10). Two types of couplings were considered: linear, $\alpha=1$; and quadratic, $\alpha=2$. We focused on the density dependence of the quartic term, which appeared to be unexpectedly large at densities a bit higher than $n_{0}$ for both kinds of interactions. In the case of quadratic coupling the quartic term is so large that it makes the power series completely useless. When subsequent terms in the power expansion at a given point of a function are rapidly increasing it could be a signal that the function is not analytical at this point. Such a function can have regular derivatives only up to some order. For example, for $f(\beta)=\beta^{9/2}+(1-\beta^{2})^{8}$ we get the expansion around 0, $f(\beta)\approx 1-8\beta^{2}+28\beta^{4}$, where the subsequent derivatives rapidly increase. Up to fourth order the function seems to be analytical, however the fifth derivative of that function does not exist. It is difficult to say whether the case is the same in our model. In principle, it would be necessary to know all higher derivatives. Even the sixth derivative is highly complex. Nonanalyticity of the expansion of the asymmetry was recently signaled in Kaiser:2015qia , where it was obtained from chiral perturbation theory with the inclusion of an iterated $1\pi$ exchange carried out by the contact interaction. It would be very interesting if, similarly, such nonanalyticity in the energy could be obtained in the RMF only after the inclusion of a specific scalar meson interaction. We are aware that our final conclusion here requires further investigation.