Low Density Neutron Star Matter with Quantum Molecular Dynamics: The Role of Isovector Interactions

Abstract

The effect of isospin-dependent nuclear forces on the inner crust of neutron stars is modeled within the framework of Quantum Molecular Dynamics (QMD). To successfully control the density dependence of the symmetry energy of neutron-star matter below nuclear saturation density, a mixed vector-isovector potential is introduced. This approach is inspired by the baryon density and isospin density-dependent repulsive Skyrme force of asymmetric nuclear matter. In isospin-asymmetric nuclear matter, the system shows nucleation, as nucleons are arranged into shapes resembling nuclear pasta. The dependence of clusterization in the system on the isospin properties is also explored by calculating two-point correlation functions. We show that, as compared to previous results that did not involve such mixed interaction terms, the energy symmetry slope L is successfully controlled by varying the corresponding coupling strength. Nevertheless, the effect of changing the slope of the nuclear symmetry energy L on the crust-core transition density does not seem significant. To the knowledge of the authors, this is the first implementation of such a coupling in a QMD model for isospin asymmetric matter, which is relevant to the inner crust of neutron and proto-neutron stars.

1. Introduction

Matter in neutron stars presents the largest densities achieved in the Universe, making their equation of state (EOS) hard to determine. Seeking the EOS of neutron-star matter (NSM) is a flourishing field of interest due to the presence of neutron rich matter with magnetic fields that can be larger than ${10}^{12}$ G with the possibility of exotic particles, and a phase transition to deconfined quark matter. The crust of a neutron star contains nuclei embedded in a sea of electrons. As the density increases from the surface of the neutron star towards its core, these nuclei undergo a neutronization process, eventually reaching a state of high neutron to proton asymmetry, which is followed by a transition to uniform nuclear matter at the core. Since matter above nuclear saturation density is unattainable in terrestrial conditions (except in heavy ion collisions with larger temperatures), neutron stars are considered to hold the key to the mysteries of dense nuclear matter.

Several approaches have been employed to study the properties of nuclear matter in the context of neutron stars. One of the prominent methods is Quantum Molecular Dynamics (QMD), which allows for the incorporation of competing nuclear forces of attraction and repulsion in dynamical simulations. QMD as a framework for simulating heavy-ion collisions was proposed by J. Aichelin and H. Stöcker [1]. Until then, nuclear matter simulations were only possible microscopically through one-body models, such as the Vlaslov–Uehling–Uhlenbeck (VUU) theory, and macroscopically by fluid dynamical models [2,3]. QMD combines classical molecular dynamics with quantum corrections, the most important of which is the Pauli principle. Peilert et al. [4] used QMD for the first time to simulate clustering in nuclear matter at sub-saturation densities. They performed uniform nuclear matter simulations with nucleons, which were sampled only in momentum space, for the density range $0<\rho <2{\rho }_0$ (where ${\rho }_0$ is the nuclear saturation density). These were then compared with simulations where nucleons were free to move in position space, showing a decrease in binding energy per nucleon ($E/A$), for the latter case, of about 8 MeV towards a more bound system for sub-saturation densities at a near-zero temperature. In the same work, the authors also took snapshots of simulated nuclear matter for different mean densities below ${\rho }_0$, which was useful to visualize clustered matter at $\rho =0.1{\rho }_0$, but did not help deduce the properties of single clusters (unless a computationally expensive time average of many simulations could be done).

Later, results for sub-saturation density nuclear matter at zero temperature were published by Maruyama et al. [5], where the the number of nucleons was significantly expanded (by ≈4 times) in the simulated infinite nuclear matter system. In addition to partially observing transient shapes like holes, slabs, and cylinders in clustered nuclear matter, they also extended the calculations to asymmetric nuclear matter, and obtained similar clusterization effects. This is necessary to evaluate the properties of NSM, which is highly asymmetric at saturation and sub-saturation densities. Further improvements to NSM simulations were made by Watanabe et al. [6] by implementing larger relaxation time scales and analyses of spatial distribution of nucleons. In a similar analysis, utilizing the Indiana University Molecular Dynamics framework, Sagert et al. [7] have shown nuclear pasta through similar 3D Skyrme Hartree–Fock (SHF) simulations. Recently, Schramm and Nandi [8] studied the asymmetry dependence of the transition density from asymmetric to homogeneous nuclear matter in the inner crust using QMD.

In this article, the asymmetry dependence shown by R. Nandi and S. Schramm [8] is modified to have better control on the symmetry energy slope (L). The inspiration is taken from the coupling of omega ($\omega$) and rho ($\rho$) meson fields in the Relativistic Mean-Field (RMF) theory. The model is first applied to isospin chains of finite nuclei, and then to nuclear matter at ${\rho }_0$. Symmetry energy at saturation density is re-evaluated along with its slope L. The primary aim of this work is to successfully control the density dependence of symmetry energy, and of pure neutron matter, by calibrating the $\omega -\rho$ type coupling according to established constraints. The expected clustering of nuclear matter at densities ≈$0.1{\rho }_0$ is also addressed.

The structure of the article is as follows: the general formalism is outlined in Section 2. Then a study of parameters of different strengths of the $\omega -\rho$ coupling in elucidated in Section 3. The conclusions are presented in Section 4, along with an outlook for the model under study.

2. Formalism

2.1. The Canonical Formalism: Hamiltonian and Equations of Motion

Quantum Molecular Dynamics (QMD) is a model used to accomplish dynamical simulations of nuclear matter by incorporating correlation effects between the constituents of the simulated N-body system. Peilert et al. [4] studied non-uniformities that give rise to clustering in nuclear matter. A model based on QMD for heavy-ion collisions through an N-body approach was proposed as early as the late 1980s (Aichelin and Stöcker [1]). The reader can refer to Ref. [9] for a thorough review of the method and its theoretical background. A brief insight into the working of QMD and the relevance to this project is provided in this section based on a review by Maruyama et al. [10].

In a Classical Molecular Dynamics (CMD) simulation of nucleons, particles are simulated as solid elastic spheres, and their motion is governed by Newton’s equations of motion. Inter-particle potentials quantify the force experienced by a particle, given the positions of other particles. QMD introduces quantum behavior to the system of nucleons by including the following modifications:

(a)

In QMD for nuclear matter, a nucleon is represented by a fixed-width Gaussian wavepacket in the form of a single particle wave function

$$\psi \left({\mathbf{r}}_{\mathbf{i}}\right)=\frac{1}{{\left(2\pi C_W\right)}^{3/4}}exp\left(-\frac{{(\mathbf{r}-{\mathbf{R}}_{\mathbf{i}})}^2}{4C_W}+i\mathbf{r}·{\mathbf{P}}_{\mathbf{i}}\right),$$

(1)

with ${\mathbf{R}}_{\mathbf{i}}$ and ${\mathbf{P}}_{\mathbf{i}}$ as the centers of position and momentum of the wave packet, respectively. $C_W$ denotes the width of the wave packet. The motion of the wave packet or ‘nucleon’ is determined by forces derived from inter-particle potentials in the QMD Hamiltonian. The total wave function of the N-nucleon system is obtained through a direct product given by

$$\mathsf{\Psi }\left(\mathbf{r}\right)=\prod _i^N\psi \left({\mathbf{r}}_{\mathbf{i}}\right),$$

(2)

(b)

The nucleon wavefunctions are not anti-symmetrized to explicitly manifest fermionic characteristics. As a result, the energy states violate the Pauli principle, as they all attain minimum energy. This problem was addressed phenomenologically (see the review in Ref. [10] for further references) by mimicking the Pauli principle through a repulsive 2-body potential called the Pauli potential ($V_{Pauli}$). The potential effectively repels nucleons with the same spin and isospin from coming close in phase space, since it is a function of both distance in coordinate and momentum space. In the ground state, nucleons have non-zero momentum values and do not all exist in the lowest energy state.

The Hamiltonian of the nucleon-nucleon interaction is given by Ref. [11]

$$\mathcal{H}=K+V_{\mathit{Pauli}}+V_{\mathit{Skyrme}}+V_{\mathit{sym}}+V_{\mathit{MD}}+V_{\mathit{Coul}},$$

(3)

where K is the kinetic energy, $V_{Pauli}$ is the Pauli potential, $V_{Skyrme}$ is the potential similar to Skyrme like interactions, $V_{sym}$ is the isospin dependent potential, $V_{MD}$ is the momentum dependent potential, and $V_{Coul}$ is the Coulomb potential. The expressions for the potential and kinetic terms are

$$\begin{array}{cc}\hfill K& \hfill =\sum _i\frac{{\mathbf{P}}_{\mathbf{i}}^{\mathbf{2}}}{2m_i},\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill V_{\mathrm{Pauli}}& \hfill =\frac{C_{\mathrm{P}}}{2}{\left(\frac{1}{q_0{p}_0}\right)}^3\sum _{i,j(\ne i)}\exp \left[-\frac{{\left({\mathbf{R}}_{\mathbf{i}}-{\mathbf{R}}_{\mathbf{j}}\right)}^2}{2q_0^2}-\frac{{\left({\mathbf{P}}_{\mathbf{i}}-{\mathbf{P}}_{\mathbf{j}}\right)}^2}{2p_0^2}\right]{\delta }_{{\tau }_i{\tau }_j}{\delta }_{{\sigma }_i{\sigma }_j},\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill V_{\mathrm{Skyrme}}& \hfill =\frac{\alpha }{2{\rho }_0}\sum _{i,j(\ne i)}{\rho }_{ij}+\frac{\beta }{(1+\theta ){\rho }_0^{\theta }}\sum _i{\left[\sum _{j(\ne i)}{\tilde{\rho }}_{ij}\right]}^{\theta },\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill V_{\mathrm{Sym}}& \hfill =\frac{C_{\mathrm{s}}}{2{\rho }_0}\sum _{i,j(\ne i)}\left(1-2\left|{\tau }_i-{\tau }_j\right|\right){\rho }_{ij},\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill V_{\mathrm{MD}}& \hfill =\frac{C_{\mathrm{ex}}^{\left(1\right)}}{2{\rho }_0}\sum _{i,j(\ne i)}\frac{1}{1+{\left[\frac{{\mathbf{P}}_{\mathbf{i}}-{\mathbf{P}}_{\mathbf{j}}}{{\mu }_1}\right]}^2}{\rho }_{ij}+\frac{C_{\mathrm{ex}}^{\left(2\right)}}{2{\rho }_0}\sum _{i,j(\ne i)}\frac{1}{1+{\left[\frac{{\mathbf{P}}_{\mathbf{i}}-{\mathbf{P}}_{\mathbf{j}}}{{\mu }_2}\right]}^2}{\rho }_{ij},\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill V_{\mathrm{Coul}}& \hfill =C_{coul}\frac{e^2}{2}\sum _{i,j(\ne i)}\left({\tau }_i+\frac{1}{2}\right)\left({\tau }_j+\frac{1}{2}\right)\int \int d^3\mathbf{r}{d}^3{\mathbf{r}}^{\prime }\frac{1}{\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|}{\rho }_i\left(\mathbf{r}\right){\rho }_j\left({\mathbf{r}}^{\prime }\right),\hfill \end{array}$$

(9)

where the nucleon mass, spin, and isospin are represented by $m_i$, ${\sigma }_i$ and ${\tau }_i$, respectively. The values of the parameters are listed in

Table 1. Parameter set for nucleon-nucleon interaction (values from Ref. [5] parameterized to reproduce properties of the ground states of the finite nuclei and saturation properties of the nuclear matter). The parameters are optimized to give $E/A$ ≈ −16 MeV for symmetric nuclear matter at saturation ${\rho }_0$.

Parameter	Value
$C_P$ (MeV)	207
$p_0$ (MeV/c)	120
$q_0$ (MeV)	1.644
$\alpha$ (MeV)	−92.86
$\beta$ (MeV)	169.28
$\theta$	1.33333
$C_{ex}^{\left(1\right)}$ (MeV)	−258.54
$C_{ex}^{\left(2\right)}$ (MeV)	375.6
${\mu }_1$ (fm${}^{-1}$)	2.35
${\mu }_2$ (fm${}^{-1}$)	0.4
${\rho }_0$ (fm${}^{-3}$)	0.165
$C_S$ (MeV)	25
$C_W$ (fm${}^2$)	2.1
$C_{coul}$	0 or 1

.

The overlap between single nucleon densities ${\rho }_{ij}$ and ${\tilde{\rho }}_{ij}$, which depends on positions ${\mathbf{R}}_{\mathbf{i}}\mathrm{and}{\mathbf{R}}_{\mathbf{j}}$, is calculated as

$${\rho }_{ij}\equiv \int d^3\mathbf{r}{\rho }_i\left(\mathbf{r}\right){\rho }_j\left(\mathbf{r}\right),{\tilde{\rho }}_{ij}\equiv \int d^3\mathbf{r}{\tilde{\rho }}_i\left(\mathbf{r}\right){\tilde{\rho }}_j\left(\mathbf{r}\right),$$

(10)

where the single nucleon densities are given by

$$\begin{array}{ll}{\rho }_i\left(\mathbf{r}\right)& ={\left|{\psi }_i\left(\mathbf{r}\right)\right|}^2=\frac{1}{{\left(2\pi C_W\right)}^{3/2}}\exp \left[-\frac{{\left(\mathbf{r}-{\mathbf{R}}_i\right)}^2}{2C_W}\right],\\ {\tilde{\rho }}_i\left(\mathbf{r}\right)& =\frac{1}{{\left(2\pi {\tilde{C}}_W\right)}^{3/2}}\exp \left[-\frac{{\left(\mathbf{r}-{\mathbf{R}}_{\mathbf{i}}\right)}^2}{2{\tilde{C}}_W}\right],\end{array}$$

(11)

along with the modified width

$${\tilde{C}}_W=\frac{1}{2}{(1+\theta )}^{1/\theta }{C}_W,$$

(12)

which is calculated in this form to incorporate the effect of density-dependent term in Equation (6) (see Section II.B. of Ref. [5] for details).

2.2. Vector-Isovector Interaction Formalism

As the numerical model to simulate nuclear matter in conditions pertaining to neutron star crusts has been outlined above, we now move on to the introduction of a nucleon-nucleon interaction potential based on the RMF $\omega -\rho$ vector interaction.

C. J. Horowitz and J. Piekarewicz [12,13] added isoscalar-isovector coupling terms to the non-linear Lagrangian for nuclear matter, and achieved softening of symmetry energy to control the neutron skin thickness in ²08Pb. They introduced a RMF Lagrangian density, where the interaction part has the following terms:

$$\begin{array}{l}{\mathcal{L}}_{\mathrm{int}}=\stackrel{-}{\psi }\left[g_s\varphi -\left(g_{\mathrm{v}}{V}_{\mu }+\frac{g_{\rho }}{2}\tau ·{\mathbf{b}}_{\mu }+\frac{e}{2}\left(1+{\tau }_3\right)A_{\mu }\right){\gamma }^{\mu }\right]\psi \\ -\frac{\kappa }{3!}{\left(g_{\mathrm{s}}\varphi \right)}^3-\frac{\lambda }{4!}{\left(g_{\mathrm{s}}\varphi \right)}^4+\frac{\zeta }{4!}{g}_{\mathrm{v}}^4{\left(V_{\mu }{V}^{\mu }\right)}^2+\frac{\xi }{4!}{g}_{\rho }^4{\left({\mathbf{b}}_{\mu }·{\mathbf{b}}^{\mu }\right)}^2\\ +g_{\rho }^2{\mathbf{b}}_{\mu }·{\mathbf{b}}^{\mu }\left[{\mathsf{\Lambda }}_4{g}_{\mathrm{s}}^2{\varphi }^2+{\mathsf{\Lambda }}_{\mathrm{v}}{g}_{\mathrm{v}}^2{V}_{\mu }{V}^{\mu }\right],\end{array}$$

(13)

where $\psi$ and $\stackrel{-}{\psi }$ are the baryon and conjugate baryon fields, respectively. V represents the isoscalar $\omega$ meson field, $\varphi$ represents the isoscalar-scalar $\sigma$ meson field, isovector $\mathbf{b}$ is the $\rho$-meson field, and the photon is denoted by A. $g_v$, $g_s$, and $g_{\rho }$ are the respective coupling constants. A similar Lagrangian with a non-linear $\omega -\rho$ interaction term is employed by F. Grill, H. Pais et al. [14] to study the effect of the symmetry energy slope parameter, L, on the profile of the neutron star crust within a Thomas–Fermi formalism.

Note that a softening of the symmetry energy around saturation can also be achieved through the use of density dependent couplings (See Figure 4 and the right panel of Figure 2 of Ref. [15]).

According to the RMF framework, the equation of motion for the $\omega$-meson field takes the form,

$$m_{\omega }^2〈V_0〉-\sum _{B=n,p}{g}_v{\rho }_B+\frac{\zeta }{3!}{g}_{\mathrm{v}}^4{〈V_0〉}^3+2g_{\rho }^2{\mathsf{\Lambda }}_{\mathrm{v}}{g}_{\mathrm{v}}^2{〈b_0〉}^2〈V_0〉=0,$$

(14)

and similarly for the the $\rho$-meson field:

$$m_{\rho }^2〈b_0〉-\sum _{B=n,p}{g}_{\rho }({\rho }_p-{\rho }_n)+\frac{\xi }{3!}{g}_{\rho }^4{〈b_0〉}^3+2g_{\rho }^2{\mathsf{\Lambda }}_4{g}_s^2〈b_0〉{〈{\varphi }_0〉}^2+2g_{\rho }^2{\mathsf{\Lambda }}_{\mathrm{v}}{g}_{\mathrm{v}}^2〈b_0〉{〈V_0〉}^2=0.$$

(15)

From Equations (14) and (15), it is clear that the mean $\omega$ and $\rho$ meson fields depend on the baryon density ${\rho }_B$ and isospin density ${\rho }_I={\rho }_p-{\rho }_n$, respectively (linearly, if we ignore higher-order contributions). Equation (13) shows how the mixed coupling $\omega -\rho$ potential part of the Lagrangian density depends quadratically on the $\omega$ and $\rho$ meson fields, from which we can conclude its dependency to be $\sim {\rho }_B^2{\rho }_I^2$. Let us approximate it for our QMD model, motivated by the density-dependent repulsive Skyrme potential as in Equation (6) with a term quadratic in both the ${\rho }_B$, and in ${\rho }_I$ as

$$V_{\omega \rho }=\frac{C_{\omega \rho }}{5{\rho }_0^4}\sum _{i,k}<{\rho }_i{>}^2<\widetilde{{\rho }_k}{>}^2,$$

(16)

where $<{\rho }_i>$ and $<\widetilde{{\rho }_k}>$ are the averaged ${\rho }_B$ and ${\rho }_I$ respectively, with the following expressions:

$$\begin{array}{ccc}\hfill <{\rho }_i>& =& \hfill \sum _{j(\ne i)}{\rho }_{ij}=\sum _{j(\ne i)}\frac{e^{-{({\mathbf{R}}_{\mathbf{i}}-{\mathbf{R}}_{\mathbf{j}})}^2/4C_W}}{{\left(4\pi C_W\right)}^{3/2}}\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill <\widetilde{{\rho }_k}>& =& \hfill \sum _{l(\ne k)}{c}_{kl}{\rho }_{kl}=\sum _{l(\ne k)}(1-2|{\tau }_k-{\tau }_l\left|\right)\frac{e^{-{({\mathbf{R}}_{\mathbf{k}}-{\mathbf{R}}_{\mathbf{l}})}^2/4C_W}}{{\left(4\pi C_W\right)}^{3/2}}.\hfill \end{array}$$

(18)

The summation needs to be calculated before squaring in Equation (16). A similar calculation has already been made for the repulsive part of the Skyrme potential.

The components of force for the $\omega -\rho$ term can be derived from the potential

$$\begin{array}{lll}-f_m^x& =& \frac{\partial V_{\omega \rho }}{\partial X_m}\\ & =& \frac{2C_{\omega \rho }}{5{\rho }_0^4}\sum _{j,k}[\left(<{\rho }_m>+<{\rho }_j>\right)\frac{X_m-X_j}{2L}{\rho }_{mj}<\widetilde{{\rho }_k}{>}^2\\ & & +\left(<\widetilde{{\rho }_m}>+<\widetilde{{\rho }_j}>\right)<{\rho }_k{>}^2\frac{X_m-X_j}{2L}{c}_{mj}{\rho }_{mj}]\\ & =& \frac{2C_{\omega \rho }}{5{\rho }_0^4}[\sum _k<\widetilde{{\rho }_k}{>}^2\left\{\sum _j{\rho }_{mj}\frac{d{X}_{mj}}{2L}(<{\rho }_m>+<{\rho }_j>)\right\}\\ & & +\sum _k<{\rho }_k{>}^2\left\{\sum _j{C}_{mj}{\rho }_{mj}\frac{d{X}_{mj}}{2L}(<{\tilde{\rho }}_m>+<{\tilde{\rho }}_j>)\right\}]\end{array}$$

(19)

where $X_m$ and $X_j$ are the x-coordinates of the centers of the positions of $m-$th and $j-$the particles, respectively.

2.3. Modeling of Infinite Systems: Achieving the Ground State Configuration

Different methods can be employed to achieve the ground state configuration of nuclear matter for a given density or temperature. Peilert et al. [4] calculated $E/A$ values for finite nuclei, and subsequently studied infinite nuclear matter using a version of the QMD model. They found that nuclear matter simulated at temperatures near $T=0$ MeV showed clustering among nucleons at sub-saturation densities. Later, Maruyama et al. [5] employed QMD to study the dynamical evolution of nuclear matter into pasta phases. In this work, we follow the method employed by Maruyama et al., obtaining the energy-minimum configuration of nuclear matter by distributing nucleons randomly in phase space, and then cooling down the system to achieve the minimum energy state of the system. This allows for arbitrary nuclear shapes and incorporates thermal fluctuations, giving an insight into the formation process of such structures.

To achieve equilibrium in the nuclear matter system, we use the following equations of motion along with damping factors ${\xi }_R$ and ${\xi }_P$:

$$\begin{array}{cc}\hfill \stackrel{.}{{\mathbf{R}}_{\mathbf{i}}}& \hfill =\frac{\partial \mathcal{H}}{\partial {\mathbf{P}}_{\mathbf{i}}}-{\xi }_R\frac{\partial \mathcal{H}}{\partial {\mathbf{R}}_{\mathbf{i}}},\hfill \\ \hfill \stackrel{.}{{\mathbf{P}}_{\mathbf{i}}}& \hfill =-\frac{\partial \mathcal{H}}{\partial {\mathbf{R}}_{\mathbf{i}}}-{\xi }_P\frac{\partial \mathcal{H}}{\partial {\mathbf{P}}_{\mathbf{i}}},\hfill \end{array}$$

(20)

where $\mathcal{H}$ is given by Equation (3) and the factors ${\xi }_R$ and ${\xi }_P$ are adjusted according to the relaxation time scale, with a fixed value of either 0 or −0.1.

The system is cooled from an initial temperature maintained by the Nosé–Hoover thermostat. The thermostat introduces additional coordinates and velocities in the Hamiltonian of the system in order to mimic a thermal bath in contact with the system. The extended Hamiltonian ${\mathcal{H}}_{Nose}$ appears as

$$\begin{array}{cc}\hfill {\mathcal{H}}_{\mathit{Nose}}& \hfill =\sum _{i=1}^N\frac{{\mathbf{P}}_{\mathbf{i}}^{\mathbf{2}}}{2m_i}+\mathcal{U}\left(\left\{{\mathbf{R}}_{\mathbf{i}}\right\},\left\{{\mathbf{P}}_{\mathbf{i}\mathbf{j}}\right\}\right)+\frac{s{p}_s^2}{2Q}+g\frac{lns}{\beta }\hfill \\ & \hfill =\mathcal{H}+\frac{s{p}_s^2}{2Q}+g\frac{lns}{\beta },\hfill \end{array}$$

(21)

where s is the additional dynamical variable for time scaling, $p_s$ is the momentum conjugate to s, $\mathcal{U}\left(\left\{{\mathbf{R}}_i\right\},\left\{{\mathbf{P}}_i\right\}\right)=\mathcal{H}-K$ is the potential which depends on both positions and momenta, Q is the thermal inertial parameter corresponding to a coupling constant between the system and thermostat, g is a parameter to be determined as $3N$ by a condition for generating the canonical ensemble in the classical molecular dynamics simulations, and $\beta$ is defined as $\beta \equiv 1/\left(k_{\mathrm{B}}{T}_{\mathrm{set}}\right)$ [16,17]. The energy of the nuclear matter system is not conserved, but ${\mathcal{H}}_{Nose}$ is. The most important variables here are $\beta =1/k_B{T}_{set}$, where $T_{set}$ is the desired input temperature, and $Q\approx {10}^8$ MeV ${(\mathrm{fm}/\mathrm{c})}^2$. More details can be found in Refs. [17,18] and sources therein.

At sub-saturation densities, local minima may take place around the actual global energy-minimum value of the ground state that the damping coefficients lead to if not chosen carefully. The simulation results should be checked to avoid local energy minima by repeating the cooling procedure.

3. Results

3.1. Simulation Procedure

Using the theoretical framework established in the previous chapters, the QMD simulation of a system of neutrons and protons is carried out. The final temperature after cooling down from a finite temperature was set to 0, so as to imitate the conditions in a neutron star’s inner crust.

A cubic box confines the nucleons. The size of the box is determined by the number of nucleons N and average density ${\rho }_{av}$. Periodic boundary conditions are imposed and the motion of nucleons is imitated across 26 cells surrounding the central primitive cell. The value of N is set to 1024, such that for homogeneous symmetric nuclear matter the number for protons/neutrons with spin up and protons/neutrons with spin down is equal (proton fraction $Y_p$ = 0.5 with 512 particles each of protons and neutrons). Hence, there is no magnetic polarization. Electrons are treated as a uniform background gas that makes the system charge neutral.

The nucleons are initially distributed randomly in phase space. The system is brought to thermal equilibrium at T = 20 MeV for about 1000 fm/c. The system, initially kept at a constant temperature by the Nosé–Hoover thermostat, is slowly cooled down in accordance with the equations of motion (Equation (20)), until the temperature is 0. To attain the ground state configuration, the simulation requires about 1–2 days of computation time to reach ${10}^4$ fm/c when carried out on the Goethe-HLR CPU cluster at Goethe-University Frankfurt. The computer code for the simulations in this project was first used for QMD calculations in Ref. [8].

The set of values for the parameters used in interaction potentials constituting the Hamiltonian Equation (3) is given in Table 1. Additionally, the set of values for the coefficient $C_{\omega \rho }$ of $V_{\omega \rho }$ in Equation (16) are listed in

Table 2. Optimized values for coefficient of $V_{\omega \rho }$.

Set	${\mathit{C}}_{\mathit{\omega }\mathit{\rho }}$ (MeV)
I	0.02
II	0.01
III	0.005
IV	−0.01
V	−0.02

.

3.2. Finite Nuclei

We first calculate the binding energies of ground states of a number of finite nuclei and their isotopes. Five different values of the coefficient $C_{\omega \rho }$ are tested. All five reproduce the trend of binding energies per nucleon of various nuclear isotopes, as can be seen in Figure 1. Individual simulated energy values ($E_{calc}$) deviate from the experimental ($E_{exp}$) counterparts [19] by less than 10% in all cases. Considering a reasonable expectation of accuracy within the QMD model employed in this paper, there is a minor spread in the calculated values. It is clear that varying $C_{\omega \rho }$ does not have a significant impact on the binding energies per nucleon of finite nuclei, which can be explained by the non-dependence of symmetry energy in a finite nucleus to its slope L, and the fact that it rather depends on other parameters: the symmetry energy coefficient at saturation density, ratio of the surface symmetry coefficient to the volume symmetry coefficient, surface stiffness and obviously the mass number of the nucleus (see Refs. [20,21] and sources therein.) The nuclei chosen are heavy with Z larger than 40. Good results for lighter nuclei are not expected, based on the results in Figure 4 of Ref. [5]. For each isotope family, three nuclei are selected with $Y_P$ ranging from 0.3 to 0.5 to analyze the effect of isospin dependent interactions.

There is an anomaly in the form of binding energies per nucleon being about 0.65 MeV too deep compared with experimental values for all nuclei. Given the realistically achievable accuracy within a molecular dynamics approach, this deviation is acceptable. Nevertheless, the model reproduces the overall trends of the binding energies of various nuclei reasonably well, for all values of $C_{\omega \rho }$.

3.3. Pure Neutron Matter

An important final test of the model is the examination of the behavior of a pure neutron gas at nuclear and sub-nuclear densities. The energy per nucleon ${(E/N)}_n$ of pure neutron matter affects the densities at which NSM becomes uniform.

For this case, the same system is adapted to simulate nuclear matter with $Y_P=0.0$, i.e., 1024 neutrons in the primitive cell without protons. The results for pure neutron matter simulations for nuclear and sub-nuclear densities are shown in Figure 2. The density dependence of neutron matter (or the neutron matter EoS) is crucial, as $E/N$ is an input in the calculation of the symmetry energy.

In Figure 2, the neutron matter EoSs for different $C_{\omega \rho }$ from the QMD model can be compared with two other non-linear RMF models (IUFSU [23] and FSUgold [24]), which also include the $\omega -\rho$ coupling. The shaded area shows the results from Chiral EFT [22], providing robust theoretical constraints for neutron-matter equations of state. For $C_{\omega \rho }$ = 0.02, the EoS indicates a bit too much repulsion around the nuclear saturation density. For all values of $C_{\omega \rho }$ at low densities, binding is weaker than expected. In spite of these issues, all parameter sets with different strengths of the coefficient $C_{\omega \rho }$ appear to be in good qualitative agreement with the constraints.

The slope of the symmetry energy L and the pure neutron matter EoS are related, shown by Equation (19) in Ref. [11]:

$$L=3{\rho }_0\frac{\partial }{\partial {\rho }_n}{\left(\frac{{\epsilon }_n}{{\rho }_n}\right)}_{{\rho }_0},$$

(22)

where the energy density of pure neutron matter is given by ${\epsilon }_n$. The slope of the neutron-matter EoS decreases as $C_{\omega \rho }$ is lowered, which is consistent with the trend of the L values in

Table 3. Symmetry energies and corresponding slope values (parabolic approximation).

Set	${\mathit{C}}_{\mathit{\omega }\mathit{\rho }}$ (MeV)	$\mathit{S}(\mathit{\rho }$)	L
I	0.02	37.40	135.26
II	0.01	35.63	102.71
III	0.005	34.72	100.41
IV	−0.01	32.23	66.38
III	−0.02	30.52	48.32

. Therefore, varying the slope (and by extension the strength of the $\omega -\rho$ interaction) has a direct impact on the densities at which neutrons drip out of nuclei, and consequently on the nuclear pasta phases in NSM.

3.4. Determination of Symmetry Energy and Slope Parameter

In a free fermion gas of nucleons, the expression for energy per particle is

$$\frac{E}{A}\approx \frac{E}{A}({\beta }_{asy}=0)+E_{sym}{\beta }_{asy}^2+...,$$

(23)

where ${\beta }_{asy}$ defined as

$${\beta }_{asy}=\frac{n_n-n_p}{n_n+n_p}=\frac{N-Z}{A},$$

(24)

or in terms of proton and neutron densities ${\rho }_p$ and ${\rho }_n$,

$${\beta }_{asy}=\frac{{\rho }_n-{\rho }_p}{\rho }.$$

(25)

For an initial determination of $E_{sym}$ and L at saturation density, a parabolic approximation is applied, such that only the lowest-order non-vanishing term in ${\beta }_{asy}$ is retained. Rewriting the equation with the approximation gives

$$\frac{E}{A}=\frac{E}{A}({\beta }_{asy}=0)+S\left(\rho \right){\beta }_{asy}^2,$$

(26)

where $\frac{E}{A}({\beta }_{asy}=0)$ = ${(E/A)}_0$ is the energy per nucleon of symmetric matter, and $S\left(\rho \right)$ is the nuclear symmetry energy. Keeping the Coulomb interaction switched off, the simulation is run for many values of $Y_P$ at ${\rho }_0$ for all $C_{\omega \rho }$ in Table 2. The values for $E/A$ are fitted in Equation (26), and $S\left({\rho }_0\right)$ is obtained as a fit parameter from the plot of energies per nucleon shown in Figure 3.

The slope parameter L quantifies the density dependence of the symmetry energy, which can be used to practically calculate the possible L values as [25]

$$L=3{\rho }_0\frac{S\left(1.1{\rho }_0\right)-S\left(0.9{\rho }_0\right)}{1.1{\rho }_0-0.9{\rho }_0}.$$

(27)

Here, $S\left(1.1{\rho }_0\right)$ and $S\left(0.9{\rho }_0\right)$ are determined with the same procedure as for $S\left({\rho }_0\right)$ described above. The obtained values for $S\left(\rho \right)$ and L are listed in Table 3.

The calculations discussed above can be improved. Chen et al. [26] suggested that the description of the nuclear matter EoS can be made better by improving on the parabolic approximation. Through a systematic study of isospin dependence of saturation properties of asymmetric nuclear matter, it was concluded that the parabolic approximation produces good results for ${\beta }_{asy}^2\le$ 0.1, but for higher asymmetries the quartic term should also be included. In this work, where higher isospin asymmetries are simulated, the fit using the function in Equation (26), as can be seen in Figure 3, is not satisfactory. The slope values for ${\beta }_{asy}^2>$ 0.1 can therefore be modified by adding a quartic term to Equation (26), which now expands to

$$\frac{E}{A}=\frac{E}{A}({\beta }_{asy}=0)+S_{\left(2\right)}\left(\rho \right){\beta }_{asy}^2+S_{\left(4\right)}\left(\rho \right){\beta }_{asy}^4,$$

(28)

where $S_{\left(2\right)}\left(\rho \right)$=$S\left(\rho \right)$ and $S_{\left(4\right)}\left(\rho \right)$ is the fourth order term of nuclear symmetry energy. The binding energies for different ${\beta }_{asy}$ values are fitted to the Equation (28) and $S_{\left(2\right)}\left(\rho \right)$=$S\left(\rho \right)$ and $S_{\left(4\right)}\left(\rho \right)$ are obtained as fitting parameters. A better fit for energy per nucleon is achieved, as shown in Figure 4. The updated values for symmetry energy and slope are listed in

Table 4. Symmetry energies and corresponding slope values (quartic approximation).

Set	${\mathit{C}}_{\mathit{\omega }\mathit{\rho }}$ (MeV)	${\mathit{S}}_{\left(2\right)}\left(\mathit{\rho }\right)$	${\mathit{S}}_{\left(4\right)}\left(\mathit{\rho }\right)$	${\mathit{L}}_{\left(2\right)}$	${\mathit{L}}_{\left(4\right)}$
I	0.02	32.81	6.56	135.52	−0.38
II	0.01	30.93	6.72	71.88	44.11
III	0.005	30.18	6.59	99.52	1.30
IV	−0.01	27.56	6.68	61.32	7.24
V	−0.02	25.88	6.63	49.32	−1.43

. Figure 5 shows the density dependence of symmetry energy for 5 different parameter sets. The difference between results obtained using different parameter sets increases with density due to the quadratic dependence of $V_{\omega -\rho }$ on baryon and isospin densities, being very small for densities below $0.5{\rho }_0$.

Comparing Table 3 and Table 4, the effect of adding a fourth order term is a decrease in the symmetry energy in all cases. However, the decrease in L is not straightforward and is only seen for Sets II, III, and IV in the quartic case (as compared to the parabolically approximated case). It is clear in Table 4 that, with a decrease in $C_{\omega \rho }$, L can be lowered to optimal values for Sets II, III, VI, and V. However, the L for Set III was expected to be lower than that for Set II, in agreement with the trend of decreasing values going from Sets I to V. It appears that this value is indeed much closer to previous results in Ref. [8], where a symmetry energy of ≈29 MeV is associated with an L≈ 92 MeV. This can be interpreted in terms of the strength of the $\omega -\rho$ interaction energy being too low for Set III, which causes the prediction to agree closely with previous results that excluded it. Note that the values of symmetry energy and slope for set I are high when compared with experimental values [27,28,29].

3.5. Nucleon Distributions

Nuclear clustering cannot only occur in NSM, but also for more isospin symmetric matter as it undergoes the liquid-gas transition. Such matter can be studied, for example, in high energy nuclear collisions. In order to bridge the gap between such studies and NSM, the proper isospin dependence of the existence and occurrence of the liquid-gas phase separation needs to be understood. In the following we will show how our model can be used to study the occurrence of clustering of nuclear matter for nuclear matter with proton fractions between $0.3<Y_e<0.5.$

The nucleon distribution of nuclear matter at $T=0$ can be visualized in the simulation box. At every grid point, the density contribution of each nucleon is added to produce a density map. Since the Coulomb interaction is included, clusterization of nucleons in a lattice-like structure is clearly seen in the system.

In Figure 6, Figure 7 and Figure 8, clusterization is seen in the system at $T=0$. At every grid point in the simulation box, the density contribution of each nucleon gaussian wave packet is added to calculate the density map. In Figure 7, an increased $C_{\omega \rho }$ decreases the density of clusters, so they seem to break up into smaller clusters of lower densities. This implies that neutrons will drip out at lower densities for $C_{\omega \rho }$ = 0.02 than for 0.01. As the density is increased 3-folds (shown in Figure 8), the density map morphs into a more interesting structure.

3.6. Transition from Clustered to Uniform Nuclear Matter

Long-range correlations between nucleons can determine the density at which a liquid-gas phase transition occurs. To this end, a useful tool to analyse the spatial distribution of nucleons is the two-point density fluctuation correlation function ${\xi }_{NN}$ for nucleon density fluctuations defined as [8]:

$${\xi }_{NN}=〈{\mathsf{\Delta }}_N\left(\mathbf{x}\right){\mathsf{\Delta }}_N(\mathbf{x}+\mathbf{r})〉.$$

(29)

Here, the average denoted by $\langle ...\rangle$ is taken over the position $\mathrm{x}$ and in the direction of $\mathrm{r}$. The fluctuation ${\mathsf{\Delta }}_N\left(\mathrm{x}\right)$ of the nucleon density field ${\rho }_N\left(\mathrm{x}\right)$ is defined as

$${\mathsf{\Delta }}_N=\frac{{\rho }_N\left(\mathbf{x}\right)-{\rho }_{\mathrm{av}}}{{\rho }_{\mathrm{av}}},$$

(30)

where ${\rho }_{av}$ is the average density of the simulation box. Two-point correlation functions for $Y_P=0.3$ and $0.5$ (for $C_{\omega \rho }$ = 0.01, 0.005, and −0.01) are plotted in Figure 9. In all cases, an increase in density decreases the amplitude of ${\xi }_{NN}$, indicating a smoother nucleon density distribution. Correlations are highest near the origin as the nucleons have the strongest influence on their nearest neighbors. This also indicates clusterization at low densities. A negative value of ${\xi }_{NN}$ at a given r implies anti-clustering or regularity, which means the point at that r has a density lower than the average density of the simulation box.

All curves at densities higher than $0.8{\rho }_0$ are almost flat-lined at ${\xi }_{NN}=0$, indicating uniform matter above $0.8{\rho }_0$. Clear trends in the variation of cluster size and densities with $C_{\omega \rho }$ and L could not be deduced.

When nuclear matter is uniform, the two-point correlation vanishes. At a certain average density, the long-range correlations suddenly disappear (instead of gradually), indicating the density turning to uniform matter through a first-order phase transition, which corresponds to the liquid-gas transition. Similar conclusions were obtained in previous studies [6,8]. For all cases of $C_{\omega \rho }$ in Figure 9, more data is needed to find out the point of transition although the transition from asymmetric to uniform matter seems to occur between $\rho /{\rho }_0$=0.6 and 0.8.

Note that it has been shown using relativistic mean field models that effects on the slope of the symmetry energy induced by an additional $\omega -\rho$ interaction affect the crust-core transition: a smaller slope reproduces a larger onset [30,31,32]. Similar results were found within the Brueckner–Hartree–Fock approach [33] and in a detailed study involving different approaches [34].

4. Conclusions

The conditions in the inner crust of neutron stars have been simulated within a Quantum Molecular Dynamics (QMD) approach with periodic boundary conditions to imitate infinite uniform nuclear matter. The nucleon-nucleon interaction Hamiltonian for QMD was successively developed in earlier works by Aichelin and Stöcker [1], Peilert et al. [4], and Watanabe et al. [6] and consists of effective interaction potentials that take into account the Pauli principle, the Yukawa interaction, Coulomb interaction, and density dependent terms. In the current project, an isospin-dependent potential term to take into account the repulsion from interaction of omega and rho mesons has been implemented in the QMD Hamiltonian. The idea to include a mesonic self-interaction in nuclear matter calculations is not new. It was first introduced in an attempt to reduce the neutron-skin thickness of ²08Pb [12]. In a work proposing the IUFSU effective interaction [23], it was shown that increasing the $\omega -\rho$ coupling constant softens the EoS of nuclear matter at around saturation density and that the density dependence of symmetry energy is highly sensitive to it. This was done within a model based on a relativistic effective field theory. Later, it was shown to improve the radius and tidal deformability of neutron stars, leading RMF models to be in better agreement with observations [35].

The new $\omega -\rho$-inspired term in the QMD model Hamiltonian in this work is inspired from the density dependent repulsive Skyrme force and depends on the baryon density and isospin density of asymmetric nuclear matter. A few values for the coefficient of the $\omega -\rho$ potential were tested, which resulted in very different behavior of the symmetry energy and its slope L. First, the values and trends of binding energies per nucleon of ground states of several nuclear isotopes were reasonably reproduced compared with experimental values. Simulations for pure neutron matter resulted in a density dependent behavior that is largely similar for all $C_{\omega \rho }$ (coefficient of $\omega -\rho$ meson field interaction) below nuclear saturation density and is in good qualitative agreement within constraints from Chiral EFT. Around ${\rho }_0$, we can see a large divergence in the trend of $E/N$ and a decrease in maximum energies as $C_{\omega \rho }$ is lowered. The numerical data obtained by simulating asymmetric nuclear matter was fitted to the energy per nucleon expanded as Taylor series, keeping both the lowest and the second highest order term. The approximation in the second-order term, also called the parabolic approximation, gave a trend of symmetry energies that decreases as the coefficient of the $\omega -\rho$ potential also decreases. The corresponding slopes L exhibit a similar trend, although four of five tested parameter sets produced L values within established constraints [27,28,29]. The higher-order approximation, which is necessary to obtain a better fitting of the data to energy per nucleon, further reduced the symmetry energies for the same coefficient values whereas there are more variations in the corresponding L values amidst a general decrease. This behavior requires the inclusion of data for higher and in-between values of proton fractions to further improve fitting and obtain better symmetry energy and slope values.

The dependence of clusterization in the system, due to the nuclear liquid-gas transition, on the isospin properties was also explored by calculating two-point correlation functions. Although a detailed study of the structure of inhomogenous phases could not be accomplished due to inaccurate Coulomb energies, a visualization of the simulated system shows interesting pasta-like shapes. The transition from inhomogeneous to uniform matter is evaluated using a two-point density fluctuation correlation function and points to a first-order phase transition. Only a small change was observed in the effect of varying L on the transition density, which cannot be deemed as significant. The properties of the mixed phase with the newly integrated $\omega -\rho$-inspired interaction can be studied in a similar fashion to a work conducted earlier in Ref. [36] with the aim of giving a better range for the critical end-point of the liquid-gas phase transition in dense nuclear matter. The analysis of two-point density fluctuation correlations also reveals the size of clumps of nucleons in the system. If the evolution of clump sizes is tracked with respect to time, one can deduce where density fluctuations are amplified enough to have matter separate into domains of high and low densities forming a coexisting phase. The growth of instabilities or fluctuations point to a region of negative compressibility in the phase diagram of nuclear matter, where at a constant temperature an increase of density results in a decrease in pressure [37]. This region is called the spinodal region. Therefore, further study is needed to shed more light on the nuclear phase diagram. Steinheimer et al. [38,39] have conducted detailed analyses on experimental signals of the expected phase transition at large baryon densities and identifying spinodal clumping in high energy nuclear collisions. Studies of temperature, pressure, and time evolution of density fluctuations in nuclear matter are outside the scope of this project but is an interesting prospect for the future.

QMD has an advantage over other types of models in the possibility to track the trajectory of nucleons and study the non-averaged properties of clusters in nuclear matter at inner crust densities, unlike mean-field approaches. Implementing the $\omega -\rho$ interaction in a QMD model is an important step in efforts to constrain the density dependence of symmetry energy and at the same time observe effects of this interaction on the structure of nuclear matter within a dynamical framework. Pressure can be calculated using the simulated data to obtain the full equation of state, and subsequently a M-R curve for the model used in this work. We also aim to find a way to reconcile the model with the causality of sound speed, which is ensured in Relativistic Mean Field models but can be problematic in microscopic simulations. An exciting prospect is finite temperature calculations to check for phases of hot nuclear matter at sub-saturation densities, which is relevant for proto-neutron stars.