Virialized equation of state for warm and dense stellar plasmas in proto-neutron stars and Supernova matter

We present microscopic Molecular Dynamics simulations including the efficient Ewald sum procedure to study warm and dense stellar plasmas consisting of finite-size ion charges immerse in a relativistic neutralizing electron gas. For densities typical of Supernova matter and crust in a proto-neutron star, we select a representative single ion composition and obtain the virialized equation of state (vEoS). We scrutinize the finite-size and screening corrections to the Coulomb potential appearing in the virial coefficients 𝐵 2 , 𝐵 3 and 𝐵 4 as a function of temperature. In addition, we study the thermal heat capacity at constant volume, 𝐶 𝑉 , and the generalized Mayer’s relation i.e. the difference 𝐶 𝑃 − 𝐶 𝑉 with 𝐶 𝑃 being the heat capacity at constant pressure, obtaining clear features signaling the onset of the liquid-gas phase transition. Our findings show that microscopic simulations reproduce the discontinuity in 𝐶 𝑉 , whose value lies between that of idealized gas and crystallized configurations. We study the pressure isotherms marking the boundary of the metastable region before the gaseous transition takes place. The resulting vEoS displays a behaviour where effective virial coefficients include extra density dependence showing a generalized density-temperature form. As an application we parametrize pressure as a function of density and temperature under the form of an artificial neural network showing the potential of machine learning for future regression analysis in more refined multicomponent approaches. This is of interest to size the importance of these corrections in the liquid-gas phase transition in warm and dense plasma phases contributing to the cooling behaviour of early Supernova phases and proto-neutron stars.


INTRODUCTION
In nature there are several astrophysical scenarios where the relevant dynamics of matter can be described in terms of screened ion systems.Examples of them are low density warm Supernova matter or the outer layers of proto Neutron Stars (NSs).Ions are immerse in a neutralizing degenerate electron Fermi sea formed at densities beyond 10 6 g/cm 3 , as a result of electrons being stripped off atoms and screening the otherwise nude Coulomb interaction.While at sufficiently low temperatures, below the melting  <  melt , matter arranges in a bcc, fcc or hcp crystallized configuration, for densities below the neutron drip (Kozhberov & Potekhin 2021), at higher temperatures and/or densities matter adopts a fluid state undergoing for this a phase transition, see Lomba et al. (1994).Such plasmas display, naturally, a complex composition (Caplan et al. 2022), however typical approaches to their modelization have considered pure systems, i.e. a one component plasma (OCP) with a single ionic species for the sake of simplicity.However, in warm and dense stellar plasmas like the ones cited, a multicomponent (MCP) scenario is the most likely, see Fig. 1 in Yudin et al. (2018).The MCP system arises due to chemical and electrostatic equilibrium at given thermo-★ E-mail: david.barbag@usal.es† E-mail: albertus@usal.es‡ E-mail: mperezga@usal.esdynamical conditions from an existing previous composition late in the Supernova aftermath.Thus in the more realistic picture, the composition displays some variability with the ion population showing some spread around a most frequent ion.
Previous works at very low densities have provided the finitetemperature EoS in a model independent way under the form of the virial EoS.In Horowitz & Schwenk (2006) second virial coefficients were obtained for a system with neutrons, protons, and -particles along with experimental information on binding energies and phase shifts.O'Connor et al. (2007) and Mallik et al. (2008) added heavier species, obtaining as a result that they contribute substantially to the virial series for pressure.For light nuclei the quantum statistical model of Röpke et al. (2013) extended the quasiparticle virial expansion to include arbitrary point-like clusters.An application of this to H plasmas has been recently developed (Röpke 2023).
The link between the macroscopic expression given by the virialized equation of state (vEoS) and the microscopic world was made in the 1930s by Uhlenbeck and Beth when they showed that the virial coefficients could be given in terms of integrals involving the interaction potential energy (Uhlenbeck & Beth 1936).Besides, in realistic cooling MCP scenarios in the partially accreted inner crust at    ≲ 1 MeV ( B is the Boltzmann constant) heat capacity at constant volume is a function of density, mostly detemined by ions as shown in (Potekhin & Chabrier 2021).Thus understanding how the crust forms from an originally warm dense plasma in a proto NS is interesting by itself.Note the densities involved are those below the neutron drip,    ∼ 10 11.6 g/cm 3 , being this value somewhat model dependent.The crust formation in a newly born NS proceeds from a    ≳ 20 MeV initial state in the aftermath of gravitational collapse to cool down to    ≲ 0.1 MeV in days to months.
Here we focus on the properties of a plasma describing the warm and dense Supernova matter in the regime of nuclear statistical equilibrium (NSE) where photo-disintegration and radiative captures shape the nuclear population or in early proto NS outer layers.Previous studies (Ishizuka et al. 2003) show the critical temperature may be very high   ∼ 14 MeV for lepton fraction   = (0.3 − 0.4), which is the ratio expected in actual supernova explosions, see Oertel et al. (2017).These critical temperatures are much higher than those in neutrino-less Supernova matter and comparable to those in symmetric nuclear matter.The boiling points in those cases remain at   ∼ 1 MeV even at very low densities for all lepton fractions, see Dinh Thi, H. et al. (2023) for a recent calculation.From this the plasma composition is determined by Saha equations or extended NSE models involving nuclei, neutrons, protons, and  particles and when T cools down below some value, the composition frozens and is essentially unchanged, see a discussion in Raduta & Gulminelli (2019).Generally speaking, to model the cooling behaviour in this complex system one should rely on screened ionic degrees of freedom and solve (ideally in radial coordinate, ) for local (non-redshifted) temperature  (), in a locally flat space-time.Assuming no external mass transfer nor radiation the heat equation can be written as where ,   and  0 are the mass density, heat capacity per unit volume at constant pressure and thermal conductivity, respectively.  is the volumetric emissivity accounting for heat transfer.
In this work we perform Molecular Dynamics (MD) simulations for a low density ion system whose composition is compatible with those obtained in the literature for NS crusts, see Pearson et al. (2018) and Murarka et al. (2022), although we note there is some spread regarding the election of the single-ion species in the OCP description for same density range.In particular, we use the innermost ion of table 4 in Pearson et al. (2018) and focus on ion densities in the interval   ∈ [2, 7] × 10 −6 fm −3 within the OCP approximation with the ion species (, ) = (38, 128) having its charge spread under a finite-width gaussian distribution with zero spin.Accordingly, the temperature range that we will explore in our MD simulations is    ∈ [2, 11] MeV, corresponding to a range of Γ  ∈ [3.84, 32.06], so that crystallization does not happen and liquid and gas phases appear.The screening is set at  = 0.622.We will consider, relying on previous works Ishizuka et al. (2003), that dissolution of the nuclear clusterized entities in our neutral neutron rich plasma   ∼ 0.3 does not take place either and, in this spirit, we extend our analysis to the relatively higher T range for the sake of completion.In that work the authors discussed that for systems with   ≲ 0.35 boiling temperatures   ≲ 10 MeV.In the same line as shown in Dinh Thi, H. et al. (2023) clusterized degrees of freedom are still valid in the warm phases within the crust density range in the proto-NS, with compositions compatible to those in the inner crust of a catalyzed star.
The simulations we present were performed using our original code USALMD Gauss−Ewald using a multicore computing setup with efficient parallelization in Fortran+OpenMP, already presented in Barba-González et al. (2022).In our runs all positions and velocities of ions are followed for at least 10 5 fm/c to assure thermalization before extracting the thermodynamical properties of the system.The NVT ensemble is used, whereas we fix the volume by using periodic boundary conditions in an Ewald summation (Ewald 1921) paradigm calculation, and temperature is fixed by periodically rescaling the kinetic energy although its convergence to thermostats has been described in the past (Pérez-García 2006).
The paper is organized in the following manner.In Section (2) we present the framework to study the screened realistic ion system at finite temperature.We obtain the pressure  from the diagonal components of the general expression for the stress tensor arising from a Yukawa-like screened electrostatic interaction of finite-size ions with the Ewald construction.More in particular we obtain the pressure at finite temperature describing in detail the virial coefficients for the reference OCP and ion species (, ) we use.From our simulated cases and to illustrate the potential of machine learning methods we perform a regression analysis and present some useful expressions using artificial neural networks (ANN) in the appendix B for the density-temperature phase space we sample.In Section (3), we further characterize the system describing the liquid-gas phase transition it undergoes by using the Ehrenfest criterion based on derivatives of the thermodynamical potential.We focus on constant pressure and volume heat capacities   ,   along with the generalized Mayer's relationship.We discuss how the non-ideal screened interaction ion system along with finite size corrections distort the ideal gas result, that is safely recovered in the low density, high temperature regime.We discuss our findings comparing to some existing warm EoS in different frameworks and conclude in Section (4).

REALISTIC VIRIAL EQUATION OF STATE FOR THE WARM AND DENSE SCREENED PLASMA
Finite temperature ion plasmas at a particular number density   and temperature  have been characterized by the dilute Coulomb fluid theory based on the parameter Γ  =  2  involving point-like particles.Here  is the charge of the ions and  = (3/4  ) 1 3 is the average interparticle distance between them.However, at large electron mass densities   > 10 6 g/cm 3 typical in the NS outer crust it is expected that the ion plasma behaviour is better described by the Yukawa theory due to the electron screening of Coulomb interaction, as their degeneracy parameter i.e.Fermi energy to thermal energy ratio,  =  , /   ≫ 1, signals their degenerate contribution.The finite electron number density   =   yields the ion system electrically charge neutral, also meaning that there is an additional parameter  =   /, with   the Thomas-Fermi screening length, that shifts the crystallization point to lower temperatures in the Yukawa treatment Γ() > Γ  , see (Vaulina et al. 2002;Khrapak & Thomas 2015).As it follows from the previous definition of  and for the OCP approximation we use, its density dependence drops off thus remaining fixed in the density range under study.On a more general MCP system, however, this should not be the case.For  = 0 the Coulomb case is recovered.However, note that a more refined treatment using Jancovici screening instead of the usual Yukawa screening may distort this behaviour, for a discussion see Potekhin & Chabrier (2021).Importantly, this description assumes ions as point-like particles with no spatial spread of the charge.When more realistic distributions are considered, an additional parameter  = 1/ √ , where  is a measure of the particle finite-width charge spread, it results a less energetically bound system, which can even be structurally affected shifting the crystallization point, see Barba-González et al. (2022).
As we will discuss, the weaker interaction in Yukawa systems, and finite size effects to a less extent, impact the EoS (, ).For a OCP description based on the most frequent ion at given density in the cold NS outer crust see Fantina et al. (2020) or Dinh Thi, H. et al. (2023) for proto-neutron stars.Besides in the non-rigid crust approximation the low T behaviour of ions can impact the shear and bulk viscosities and mechanical properties with important consequences on thermal properties, bursts (Yakovlev et al. 2020) and the ability of the star to sustain oscillation modes (Fantina et al. 2018;Suleiman et al. 2022;Gusakov & Chugunov 2021;Samuelsson & Andersson 2007) just to cite some previous works.
As already presented in the context of realistic crust calculations in Barba-González et al. (2022) the Ewald sum procedure allows the very efficient calculation of high order correlations in systems with two-body interactions characterized by a non-zero spatial scale .Examples of this have been already discussed in the literature, see Salin & Caillol (2003) or Watanabe et al. (2003).
In our calculation we consider the electron-screened Coulomb potential created at distance |ì  − ì   | by the ith ion at position ì   i.e. the Debye-Hückel potential.In the static and point-like ion approximations the potential displays a Yukawa-like form where   ≡  is the Thomas-Fermi screening length and  , and   are the electron Fermi momentum and mass, respectively. is the fine structure constant.Electron number density is written as for a set of ions  = 1, ...  , being in general a multicomponent system.
In our scenario thermal quantum effects are constrained by the fact that they will not disturb the individual treatment of ions and collective of electrons as long as the thermal de Broglie wavelengths of the ions and electrons behave accordingly.Usually they are defined as , respectively.The quantum effects on ion motion are important either at   I ≳  I , where   ≡  or at  ≪  p , where  p ≡ ℏ p /  ≲ 10 8 K is the plasma temperature determined by the ion plasma frequency . Note that a fermion/boson ion idealized gas EoS will not be affected by such quantum corrections under this treatment, so that our MD treatment remains fully consistent for the  −  ranges under study.
Following Barba-González et al. (2022) we will be considering a gaussian shape charge distribution for ions in the form with   a characteristic width parameter depending on the ion (, ).
For details we refer to Appendix A. For the non point-like charges the ionic potential is obtained by solving the Poisson equation and is given by Salin & Caillol (2003), and used by Barba-González et al. (2022) in the following form with erfc the complementary error function.
From the previous we can use the efficient Ewald procedure and separate the potential energy expression  into short, long and selfinteraction contributions.We detail the technical steps in Appendix A. Thus one can write  =  short−range +  long−range −  self and similarly, the force over a particle  be decomposed as a result of twobody contributions ì These lengthy expressions include not only the real and screening corrections but also the finite-size ion gaussian charges within the Ewald summation frame.It is important to carefully implement the efficient sum, as the screening charge associated inverse square width  Ewald also controls how the summation is distributed between the short-range and long-range parts.It has been fixed ensuring that the minimum image convention is sufficient in real space (Holden et al. 2013).
The energy density  is obtained in the system as the contribution from the kinetic and potential parts where   is the Kronecker delta and ,  = 1, 2, 3 are tensorial indexes for the XYZ spatial coordinates.  =    3 is the number of ions considered in our system and  3 =  is the volume of the cubic simulation box.
In order to obtain the pressure  in the system from the Ewald construction we first obtain the general stress tensor Π tot  tailored specifically to use the Ewald procedure under the screened Coulomb potential created by the gaussian ion sources.Let us note that in the fluid phase the off-diagonal components are zero as there is no shear.
From the stress tensor diagonal components we obtain the corresponding pressure  = 1 3 Tr Π tot  .This second-rank symmetric tensor is obtained from contributions of the real and Fourier reciprocal space as follows Subsequently, it can be put into the form (see Appendix A)

𝛼𝛽
. As the interactions in real space may be written in a pairwise fashion, Π real  can be evaluated directly from the real forces by employing the virial expression The reciprocal part is obtained including the widths   and  Ewald associated to the gaussian charge spread and Ewald construction respectively.For point-like approaches see (Aguado & Madden 2003;Salin & Caillol 2003) coincident with the stress tensor calculations in Nosé & Klein (1983 This expression is obtained by taking the derivatives of the total energy with respect to generic deformations in the simulation box.The volume part of the tensor comes from the last term in  short−range (Eq.(A4) in the Appendix).It does not depend on the particle positions but it depends on volume, so that it adds to the stress tensor in the same way the reciprocal part in Eq. ( 10) does, as the derivative of volume with respect to deformations is finite.Then where  short−range is given by Eq. (A1).We now proceed to express thermodynamic pressure in a suitable way in the screened finite-size ion system.In the point-like approach the virial expansion series adopts the generic form (Kamerlingh Onnes 1902) with  1 ≡ 1 and   with  = 2, 3, 4....,  being the T-dependent virial coefficients up to kth order.In particular, the  = 2 term accounts for the contribution of two-body interactions and has the form while effective three-body interactions are approximately described using the coefficient and a similar expression for  4 ().Eqs. ( 13), ( 14) are expressed in terms of Mayer cluster integrals (Mayer & Montroll 1941) involving − 1 so that the larger |  | is, the more non-ideal the gas behaves.  with  ≥ 2 give the corrections to the ideal behaviour.Having obtained the ∼ O ( 3  ) EoS it is possible to obtain the virial coefficient  4 , and so on.So, from the  2 , it is possible to compute approximately the full correction due to interactions in the procedure called Kirkwood superposition approximation.In our work we will restrict to fourth order vEoS.Additional refinements are indeed possible, so that coefficients in the expansion yield information beyond a selected order in many-body correlations in the system, see Oertel et al. (2017) for a discussion.
In some works where interactions are relatively weak, a comparison is performed with respect to the ideal gas case internal energies and pressure,   ,   , see Shahzad & He (2012).Therefore real plasma thermodynamic properties are expressed in terms of reduced values i.e. /  and /  , so that typically a departure from unity shows the degree of changes with respect to ideal fluid.
The virial EoS as laid out in Eq. ( 12) represents the limit in which the degrees of freedom of the system are point-like particles.In our previous work (Barba-González et al. 2022) we showed the effect of introducing a gaussian ion charge spread in the crystal lattice formation.It is also expected that the EoS, and in turn the virial coefficients are tuned by this effect.The finite size of the particles is usually introduced by means of an excluded volume term ( −   ) in the Van der Waals approach.
In the low density and non-relativistic ion dynamics involved in this scenario, we seek to model the dependence on the finite width of the (, ) nucleus with charge spread  in terms of the adimensional parameter  = 1/ √ .Note  = 0 for point-like particles as the inverse width  → ∞.Each ion is modeled as a finite-size gaussian charge density distribution as in Eq. ( 5) where specifically , and following previous work (Barba-González et al.
In order to illustrate the systems under study in Fig.
( The plotting software we use is Paraview (Ahrens et al. 2005).We can see how in the more sparse system the weaker interactions will provide a closer behaviour to ideal, however it is expected that high order correlations due to density-dependent corrections will not be negligible.The depicted two crust ion densities in the OCP approximation used differ roughly an order of magnitude.There is a qualitative, visual, difference in the density profiles: whereas in the smaller density there are bigger voids and high density regions and in denser sample matter is more compactified in smaller clusters.It is important to remember that the MD simulation degrees of freedom are ions, and eventually as it cools down it approaches the crystallized phase, fulfilled at higher densities/lower temperatures than the ones shown.
In Fig.
(2) we show the EoS for both point-like and realistic ions from our MD simulations, as described in Section (1), i.e. reduced ionic pressure versus density for a range of temperatures and densities corresponding to Γ  ∈ [3.84, 32.06].The screening is set at  = 0.622.It is seen that the pressure tends to the ideal gas value when the ionic density is reduced, as expected from Eq. ( 12), for each set temperature.The smaller the temperature the lower the pressure goes, even as far as becoming negative for    ≲ 4 MeV.This behaviour has already been reported for less refined Yukawa OCP not including Ewald summation and gaussian particles, see, for example, Khrapak et al. (2015) and Shahzad & He (2012).
From these curves in Fig.
(2) we fit the virial expansion to our OCP data with ion (Z,A) = (38, 128) using least-squares and obtain the virial coefficients from either Eq. ( 12) or Eq. ( 15) for pointlike or gaussian ion runs.All the fits converge and are adequate for the representation of the numerical outputs within the gas regime, when the liquid-gas transition has already ocurred, in the temperature range    ∈ [7, 11] MeV (Γ  ∈ [3.84, 9.16]).The screening is set at  = 0.622.We show in Figs. ( 3), ( 4), ( 5) the virial coefficients obtained by this procedure, as a function of temperature.In the figure, coefficients   have been conveniently normalized by multiplying them by 1  −1 where  is a typical volume for our simulation box.It is seen that, as expected for the screened Thomas-Fermi interaction, two-body forces are dominant and so from Eqs. ( 13), ( 14), | 2 | > | 3 | > | 4 | for the full range of temperatures.This is also consistent with coefficients of increasing order  decreasing accordingly to add up into the convergent series in Eq. ( 15).In physical terms, this decline in the coefficients can be better understood when looking at the density profiles of the resultant fluid system with fixed (Z, A) evolved by our MD code, see Fig.

(1).
The effect of the ionic spread and density correlations is shown in Figs.
(3), ( 4) and ( 5) for  2 ,  3 ,  4 , respectively.We label the virial coefficients  ′  () for the gaussian case and   () for the point-like one (same expression as in Eq.( 15) when  = 0).They are perfectly compatible within statistical uncertainty for most of the temperature range considered.Together with them, we consider an effective virial coefficient that we define as  ,eff = (1 + )   ′  and plot it in shaded blue color band for all densities considered.The figures clearly show that density correlations along with finite ion size effects are included in the  parameter and play a key role so that the virial expansion coefficients get clearly shifted as interactions are smeared out when ion charge is spatially distributed.This fact, combined with the aforementioned compatibility of the coefficients, means that the complete information about the finite ionic spread is encoded in the factors (1 + )  that arise in each term of the virial expansion.This, in turn, supports the choice of parametrization in Eq. ( 15).Thus we not only show the existence of a measurable effect of the more realistic ions, but also are able to compute these effects in terms of the OCP virial expansion excluded volume factors encoded in  that can distort point-like values ∼ 25% to 100%.We emphasize that a fully consistent MCP system will yield a further ,  dependence confirming a robust effect already presented in Barba-González et al. (2022) along the whole range of densities of interest.Even a taylored A-dependent ion OCP system arising from energy density functionals will reveal interesting features regarding the combined effect of heavy species and composition (Barba-González et al. 2024).This is typically overlooked in current works, mostly restricting to nuclei with  ≲ 8.In line with this, even the presence of exotic species like the tetraneutron have been shown to largely distort the neutron rich ion fractions, and thus the EoS, see Pais et al. ( 2023).

LIQUID-GAS PHASE TRANSITION AND HEAT CAPACITIES IN THE VIRIALIZED EOS
The low density screened ion system we simulate in the range of temperatures and densities    ∈ [2, 11] MeV,   ∈ [2, 7] × 10 −6 fm −3 can be characterized by the isotherm pressure curves.This is equivalent to a range of temperatures and densities corresponding to [3.84, 32.06].The screening is set at  = 0.622.In Fig. ( 6) we plot the pressure as a function of simulation volumes using realistic gaussian charge distribution (solid) and point-like (dashed) ions.
Our set of MD simulations is used to perform an analysis of the characteristic features of the liquid-gas transition.Due to the nonlinear aspects present we choose the powerful techniques of Artificial Neural Networks (ANN), see appendix B. In brief, they consist in a set of layers i.e. input layer, hidden layers and output layer.The hidden layer can be more than one in number, thus referring to deep ANN with each layer consisting of  number of neurons.Each layer will be having an activation function associated with each of the neurons.The activation function is the function that is responsible for introducing non-linearity in the relationship, typically a sigmoid ().In our case, the output layer must contain a linear activation function as the output (pressure) has in principle no bounds.Each layer can also have regularizers associated with it.Regularizers are responsible for preventing overfitting.ANN consist of, first, a forward propagation phase i.e. is the process of multiplying weights with each data entry and adding them to the bias.Later, a backward propagation phase updates the weights in the model.Backward propagation requires an optimization function and a loss function, to be selected by the actual implementation of the ANN.
In appendix B we provide an optimized parametrization with two perceptron layers, each with two neurons capable of reproducing the non-linearity of the realistic warm ion gas virialized EoS with improved Ewald sums.The explicit form reads  = 0.1407320023 ×  2 ( ñI , T) + 10.61859989 ×10 −6 MeV fm 3 , (16) where  2 ( ñI , T) is the output of the second perceptron layer depending on reduced input data ñI , T. This parametrization produces a relative error of 1.91 % when predicting pressures within the sample provided in this work.Note that in this example illustrating the OCP results we have kept fixed the ion (, ) but without loss of generality this information can be included when constructing and training the ANN.This constitutes a first step and proof of concept of the powerful techniques of machine learning in this context.A more refined version for the full crust will be provided in Barba-González et al. (2024).
In Fig.
(2) the difference between the ideal case (dotted lines) and the more realistic case was already depicted under the form of reduced pressure values.From Fig. ( 6), we can see the pressure for the spread out charges is always bigger, as the system is less bound.Note this is consistent with our previous results (Barba-González et al. 2022) regarding the crystal lattices for realistic ions in OCP or MCP mixtures.Equally consistent is the fact that for bigger volumes (lower densities) the finite size effect has a reduced impact, justifying the  parametrization used for the virial expansion in Eq. ( 12) developed above.
We observe there is a clear sign change in the derivatives     that will characterize the thermodynamics as the system expands or contracts, that is a projection of the vEoS over the  −  plane.Note that in our setting the Gibbs energy is only dependent on ,  as the number of particles remains fixed for each simulation (no species conversion) and no external work is exerted on the box walls.In the pressure isotherms depicted we can clearly differentiate three regions, according to a critical temperature   , usually coined as subcritical ( <   ), critical ( ∼   ) and supercritical ( >   ).Regarding the OCP the Maxwell construction determines the coexistence region of liquid and gas, and the critical isotherm with  =   that in our setting is located   ≃ 6.5 MeV.Although in our OCP we do not consider variability in the ion species we expect some dependence on the charge spread (ion composition).
The thermal properties of the NVT system undergoing a phase transition are better understood when looking at the heat capacities, i.e. the thermodynamical potential response to energy variations with temperature.Heat capacities are defined, generically, as with   and   the heat capacities at constant pressure and volume, respectively.
In our screened system the thermal behaviour is determined by the impact of the electron component.This is in-built in the screened potential caused by effective forces driving the inter-particle interaction.Generically the two components that largely dominate the specific heats in the crust behave in a differentiated way.On the one hand the specific heat capacities   ,   for the ultra-relativistic electrons   ∼ 0 forming a highly degenerate electron Fermi gas will contribute to the heat capacity per unit volume (Potekhin & Chabrier 2021) at  ≪  F, with  F, =  , /  as It is usually assumed that in the outer crust the main contribution is given by that of the ions since electrons are highly degenerate and  , ≈  ,  .In an ideal classical crystal, the ion heat capacity is  I = 3 B while for ideal gas ∼ 3/2  .Quantum effects strongly reduce, however, this value for temperatures much lower than that  [5.82, 21.37]) and   = 1.1 × 10 −4 (orange line, Γ  ∈ [14.60, 53.53])for a common screening parameter  = 0.622.For the lower density case there is a discontinuity signaling the onset of the liquid-gas phase transition while this happens at higher  for the   = 1.1 × 10 −4 case, thus not shown.Realistic ion sizes (solid lines) shift   per ion as compared to point-like (dashed lines) ions.We include associated error bands corresponding to our simulations and fitting procedure, see text for details.associated to ion plasma frequency that we will not consider in our scenario.
As seen from the  isotherms shown in Fig. ( 6) the sign change in the derivative     marks the onset of a 1st-order transition.Namely for our ion species the critical isotherm lies approximately   ≃ 6.5 MeV, and has an associated flat profile of pressure as a function of V, i.e.     = 0. We also note that below and above the critical isotherm the ion system has a discontinuity   =    as depicted in Fig. ( 7) for an illustrative selected ion density   = 7 × 10 −6 fm −3 (blue lines).We have corroborated this by fitting for the slope of the () data in the two different temperature ranges, before and after the expected first order phase transition.We include associated error bands corresponding to the fitting procedure from our statistically distributed simulation energy values.For the sake of comparison we also plot the case with   = 1.1 × 10 −4 fm −3 (orange line) as it does not undergo any phase transition at this temperature range.We note that in both cases the value is largely increased ∼ 40% with respect to the ideal gas value ∼ 3/2  in a non-magnetized system.Realistic ion sizes (solid lines) shift   per ion as compared to point-like (dashed lines) ions less than ∼ 5%.
In addition to the ion heat capacity at constant volume we have also studied the ion heat capacity at constante pressure.Note that since we simulate our ion system under the NVT ensemble (where pressure is not longer constant) we calculate this quantity indirectly as deduced from the generalized Mayer relation   −   , a well known relation between the isochoric and isobaric heat capacities.The explicit form is In the above equation    is the inverse of the thermal expansion coefficient, and    is 1 times the inverse of the isothermal compressibility.The system that we are simulating, the cooling crust of a proto-NS or warm Supernova matter, is a non-ideal gaseous-liquid system, which at very high temperature is expected to reproduce the ideal gas Mayer's relationship.This is illustrated by substituting the virial expansion Eq. ( 15) into Eq.( 19).In this way for the realistic vEoS with gaussian ion corrections we obtain where we define  ′ 1 = 1.Note that for a point-like approximation in a non-interacting system we recover the usual ideal Mayer relation From the − behaviour in Fig. ( 6) of the critical isotherm limiting with associated vanishing of the derivative     = 0, together with the fact that the derivative     is finite, means that the difference between   and   from Eq. ( 19) diverges.The conclusion that it is   which does diverge can be arrived to by means of realizing that   vs  curves are smooth for any of the densities in the studied range, thus   always remains a finite value.This further confirms the phase transition as first-order.
In Fig. ( 8) we show the Mayer's relation for two densities   = 3 × 10 −6 fm −3 (blue line) and densities   = 7 × 10 −6 fm −3 (orange line) for finite size ion species (, ) = (38, 128).We can clearly see that divergence signaling the liquid-gas phase transition appearing at two different temperatures (larger as density grows) corresponding to Γ  ∼ 10 for both densities.The metastable region with phase coexistence is depicted to the left of each asymptote (dashed vertical line in the shaded regions).The negative values   −   must be understood as the presence of a coexisting liquid-gas phase.Note that in the finite temperature ion gas phase   / ∼ 8 being a factor nearly 4 larger than its isochoric counterpart   ∼ 4  .This fact may further impact the cooling behaviour of species in gas phases in the warm matter later forming the outer NS crust as derived from the heat equation Eq. ( 1).
In order to compare our findings in the OCP setting used in this work to other different warm EoS in the literature incorporating clusters, we have selected a canonical example of the Relativistic Mean Field (RMF) model calculation i.e.TM1e (Shen et al. 2020) including clusters up to  = 8.It is important to emphasize that in our scenario we deal with OCP systems, so that pressure and energy density are generated by a single species while in the TM1e this is done by a few of them arising from thermodynamic conditions previously settled.Typically, the heavy species formed near the liquid-gas boundary for Supernova matter, as presented in Ishizuka et al. (2003), are usually not considered in this scenario as the lepton rich matter (and later on neutrino-less matter) cools down.Our scenario agrees with that presented by Ishikuza et al as our transition arises at about ∼ 1.4 boil their boiling temperature at the same densities but computed in a very different framework.In these RMF models, due to the light nuclei involved in the calculations, see also Pais et al. (2023) in the same fashion, we expect that pressure obtained will be systematically higher.In the work of Murarka et al. (2022) using finite ion degrees of freedom in the OCP scenario, the NS outer crust EoS at  = 0 was constructed based on the predictions of a Machine Learning algorithm from AME2016 nuclear dataset immerse in a relativistic electron gas and lattice dynamics.As a NS cools as a function of temperature for densities   = 3 × 10 −6 (blue line) with Γ  ∈ [5.37, 24.17] and   = 7 × 10 −6 (orange line) corresponding to Γ  ∈ [7.12, 32.06] for finite size ion species ( , ) = (38, 128) corresponding to  = 0.622 The metastable region with liquid-gas coexistence is depicted to the left of each asymptote (dashed vertical line in the shaded regions).
down, a possible MCP distribution becomes more peaked so that at cristallization one can consider that the composition is frozen and unchanged to the  = 0 ground state.In this spirit we have used the composition for the NS outer crust computing warm EoS values in the density-temperature range we study.For the sake of discussion, we considered illustrative conditions near the liquid-gas phases at  = 10 MeV and neutron rich clusterized matter at   ∼ 0.3 were attained.For an ion density   = 2 × 10 −6 (  = 2.56 × 10 −4 fm −3 ) which gives a Coulomb parameter Γ  = 4.22 our OCP simulation yields  0,OCP = 1.0672×10 −5 MeV fm −3 while  0,OCP = 0.240310 MeVfm −3 .In the RMF value pressure is a factor ∼ 300 larger i.e.  Shen ∼ 300 0,OCP with a similar value of  Shen ∼  0,OCP .Instead for the cold OCP approach with degenerate electrons deduced by Murarka et al. (2022)  Murarka ∼ 50 0,OCP and similar value for  so that  Murarka ∼  0,OCP .In this last case the ion predicted is (, ) = (36, 116) at the fixed density value.Both warm RMF and cold NS outer crust EoS yield systematically larger pressure values when compared to our simulated samples.Thus we conclude the finite size ion degrees of freedom can not be overlooked in the SN matter or proto NS involving heavy ions as the full picture can not be understood without them.The determination of the liquid-gas transition is key to the conditions of heavy cluster formation, and temperatures up to a dozen MeV will play a role, thus having an impact on the outer layers of cooling proto NS.

CONCLUSIONS
We have performed microscopic Molecular Dynamics simulations with periodic boundary conditions of warm and neutral plasmas with spin-zero ions interacting by electron screened Coulomb potentials at finite temperature.We have used efficient Ewald sums for finite size gaussian ion charge distributions in the system.We have focused on the OCP description with single ion species characteristic of a given density in the crust.These conditions are of interest to describe low density Supernova matter and early phases of proto neutron stars.Although more refined treatments involving MCP systems are indeed possible we study and characterize this type of matter and the liquid-gas phase transition it undergoes and find distinctive features with our improved treatment.We compute numerical expressions for the energy density and stress tensor for screened gaussian ion charges including the additional screening  Ewald introduced by the Ewald summation method for the first time.By using this method low density thermodynamical quantities are accurately described to longer distances.In addition our system also displays a typical Thomas-Fermi screening length from the relativistic degenerate electron plasma that we effectively take at  = 0.
We dynamically follow the positions, velocities and accelerations from the prescribed set of conditions in a NVT ensemble that we thermalize before results are obtained in a temperature range where crystallization does not happen.For the fluid phases we describe a set of pressure isotherm curves   () as a function of volume where metastable coexisting (liquid-gas) and stable (gas) phases appear as T grows.We are able to determine differences among the pointlike and finite size ion treatments for the selected species so that when parametrizing the gas EoS clear differences appear in the virial coefficients up to 4th order, the meaningful range where we restrict our calculations.For the systems we explored these corrections to   can be as large as ∼ 25 − 100%.In particular we extract density and ion size dependence (by means of an extra parameter ) in terms of a 1 +  power expansion, so that at low densities we consider the excluded volume effects.In the usual treatments this parameter  = 0 as ions are considered point-like.In our OCP calculation we also find dependence of the ion mass although we can not further size its dependence leaving this for a future contribution.When comparing our findings the warm OCP yields systematically smaller pressures than RMF or clusterized matter using the same composition predicted for  = 0 NS outer crust.Although a more refined MCP system will improve this picture we expect this reduction is robust and due to finite-size ions.
As an illustrative example of usability we provide an ANN giving a rather efficient expression of the virial EoS in our OCP system (,   ) leaving for a future work a more detailed MCP calculation.As we can capture the dynamical transition liquid-gas in the MD treatment we are able to determine the energy and pressure derivatives involved in the heat capacity at constant volume and constant pressure by generalizing the Mayer relation for the Yukawa potential system with finite size particles for the first time.We find the ion heat capacity   is no longer T dependent with its value being ∼ 40% larger than the usual non-magnetized ideal gas value ∼ 3/2  .We can clearly see that the discontinous jump in the   signals the liquid-gas transition critical temperature in our simulated system.This confirms the first-order phase transition.From this we are able to determine   using the generalized Mayer relation accessible from data in our simulation.We find both ion heat capacities are different in the gaseous phase so that cooling behaviour of systems kept under constant V or P would be different.We foresee that with a MCP or OCP parametrization of crust densities Supernova matter or proto-NS early cooling phases could be affected under the usual paradigm of equivalent idealized ion   ,   values.This calculation constitutes an example of the complementary study of microscopic models based on energy functionals where particle dynamics can provide the effective description of screened ions with improved force and energy summed techniques.Early phases of Supernova matter where a distribution of nuclei with some heavy species and low mass nuclei are present can benefit from our virialized EoS at finite temperature.Preliminary comparisons with RMF and energy functional minimization from AME2016 datasets show a distintive feature, with those EoS pressure values over the ones estimated in the OCP scenario even is the liquid-gas boundary yields similar boiling temperatures.Composition is thus playing an important role and heavy species with  > 50 must be present in the consistent picture when including the proto NS.Note also that under dissolution or charge-changing reactions this description must incorporate different species as well.Even if OCP is most surely not realized in the polluted scenarios in the aftermath of a Supernova explossion we foresee our results remain valid for each species component with finite size ion charge (, ) at given density during its associated lifetime.

Figure 2 .
Figure 2. Reduced pressure, /   ≡ /  , as obtained from our MD simulations with screened Coulomb, Ewald sums and finite ion width, as a function of ionic density for different temperatures.We use OCP with (Z,A) = (38,128) for a range of temperatures and densities corresponding to Γ  ∈ [3.84, 32.06].The screening is set at  = 0.622.At low densities and high T ideal (expanding) behavior is recovered approaching unity, while the pressure reverses sign and becomes attractive for lower T. Solid color bands indicate gaussian ions, whereas dashed ones depict point-like particles.Band width represent standard deviation uncertainty from our MD simulations.

Figure 6 .
Figure 6.Pressure as a function of volume for OCP samples at different temperatures.We use ions ( , ) = (38, 128).Data are the same as in Fig. 2 with Γ  ∈ [3.84, 32.06] and  = 0.622.Ranges reflect uncertainty intervals corresponding to one standard deviation on the pressure.Solid and dashed ranges represent gaussian and point-like ions, while dotted lines are the ideal gas values  ∝ 1  .One can see a sign change in   as a function of temperature, going from positive to negative derivative as  grows.

Figure 7 .
Figure 7. Specific heat per ion at constant volume,    , for screened interactions using Ewald sum at densities   = 7 × 10 −6 (blue lines, Γ  ∈[5.82, 21.37]) and   = 1.1 × 10 −4 (orange line, Γ  ∈ [14.60, 53.53])for a common screening parameter  = 0.622.For the lower density case there is a discontinuity signaling the onset of the liquid-gas phase transition while this happens at higher  for the   = 1.1 × 10 −4 case, thus not shown.Realistic ion sizes (solid lines) shift   per ion as compared to point-like (dashed lines) ions.We include associated error bands corresponding to our simulations and fitting procedure, see text for details.