Dynamical Phase Transition in Neutron Stars

One of the most exciting topics in astrophysics is the study of compact objects (white dwarfs, neutron stars [NSs], and black holes). NSs particularly have proven to be an excellent astrophysical laboratory to study the properties of matter under extreme conditions of high densities and low temperatures (see, e.g., Weber 1999; Glendenning 2000). It is a complementary approach to the study of high-temperature relativistic heavy-ion collisions, which explores the high-temperature low-density regime.

Baade & Zwicky (1934) first gave the hypothesis that compact stars are made of neutrons. The verification of this conjecture was only established after the discovery of pulsar (Hewish et al. 1968) and connecting them with rotating NS (Gold 1968). At the theoretical end, Tolman, Oppenheimer, and Volkoff were the first to solve the hydrostatic equilibrium configuration of NSs within the framework of general relativity, also known as TOV equation. To close the problem and to calculate the mass and radius of NS, one needs an equation of state (EoS) that describes the physical properties of matter at those densities. The search for a proper EoS has motivated many researchers to study purely nuclear NS matter, commonly known as nuclear matter (NM).

It is known that the supernova matter is composed of relativistic neutrons, protons, electrons, and degenerate neutrinos. It is also expected that the fundamental component of NS would have these relativistic particles, along with some fraction of muons in them. They can even have massive hyperons and some condensate in them. Since the discovery of the first pulsar, there have been numerous studies and constant improvement while describing the matter that makes up NSs. The studies include a large number of many-body calculations where the nucleons interact via scalar, vector, and pseudo-vector mesons. The maximum mass of such NSs lies in the range between 1.8 and 2.5 M_⊙ (Fuchs 2008; Li & Schulze 2008). These calculations are compatible with the measured value 1.97 ± 0.04 M_⊙ for the mass of the pulsar PSR J16142230 (Demorest et al. 2010; Antonidis et al. 2013). The presence of heavy hyperons tends to make the EoS softer, thereby lowering the maximum mass, which is not compatible with new mass measurement. However, some recent relativistic mean-field calculations suggest that NS can support 2 M_⊙ stars, even after including hyperons.

The study of the evolution of NSs has revealed that NSs are nothing but the remnant stars after supernova explosion, the catastrophic gravitational collapse of massive stars (mass M_⊙ > 8 M_⊙) (Haensel et al. 2007), at the end of their evolution. They have a radius of the order of 10 km and a mass of around 2 M_⊙ (obtained by solving the TOV equation). The central core of a NS has a density of about few times 10¹⁵ gm/cc, and the surface magnetic field ranges from 10⁸ to 10¹⁵ G. At such extreme densities in the central core, the nuclear matter is no longer the stable ground state. It is prone to convert to three-flavor quark matter (QM; consisting of up, down, and strange quarks), which is the ground state at such densities. The strange quark matter (SQM) conjecture (Itoh 1970; Bodmer 1971; Witten 1984) was supported by model calculations (Alcock et al. 1986). The most straightforward model that describes the properties of QM at such high densities is the MIT bag model (Chodos et al. 1974). New refined models based on results from recent experiments in laboratories have also been proposed (Nambu & Jona-Lasinio 1961; Ghosh et al. 2006, 2008). However, the quark sector is still not well understood, as the nature of strong interaction at extreme condition remains a challenge.

Thus, NM at the core of NSs (high density) is likely to be unstable against QM and will eventually decay. The problem of astrophysical phase transition (PT) is a primary focus of this article. The PT from ordinary nuclear matter to QM can take place at the cores of NSs where density is expected to have a value of few times nuclear density. Compact stars, therefore, can be made of either nuclear matter or QM. Stars that contain only nuclear matter are called NS, and stars having some quantity of deconfined QM (which can be either two-flavor QM or stable SQM) in them are called quark stars (QSs). The size of the core depends on the critical density for the quark-hadron PT and the EoS describing the matter phases.

It is most likely that the PT happens due to some sudden density fluctuation at the center of the star, which induces a sufficient density and pressure discontinuity. As the discontinuity (usually a shock) propagates outward, a combustion process starts, resulting in the PT of NM to QM. Recent studies have explored different PT happening in hot proto-NSs (Bombaci et al. 2004; Drago et al. 2004; Mintz et al. 2010; Gulminelli et al. 2013) and cold NS with hyperons (Iida & Sato 1998). Major rearrangement in the star interior is expected due to the PT, which would subsequently release a large amount of energy, yielding some detectable consequences like gamma-ray bursts (Bombaci & Datta 2000; Berezhiani et al. 2003; Drago et al. 2004; Mallick & Sahu 2014), changes in the cooling rate (Sedrakian 2013), and the gravitational wave (GW) emission (Lin et al. 2006; Abdikamalov et al. 2009). If detected, all these signals could give valuable insight to constrain the properties of matter in NSs.

There has been a lot of studies of PT in the literature; however, the exact nature remains uncertain. In the literature one can find two very different scenarios: (i) the PT is a slow deflagration process and never a detonation (Drago et al. 2007; Mallick & Singh 2017), and (ii) the PT from confined to deconfined matter is a fast detonation-like process, which lasts about 1 ms (Bhattacharyya et al. 2006, 2007; Mishustin et al. 2015). If the process is quick burning and very violent (detonation), there could be compelling GW signals coming from them, which could be detected at least in the second or third generation of VIRGO and LIGO GW detectors (Lin et al. 2006; Abdikamalov et al. 2009). The earliest calculation (Olinto 1987) assumed the conversion to proceed via a slow combustion, where the conversion process depends strongly on the temperature of the star. Later, Horvath & Benvenuto (2013) studied the stability of the conversion process and concluded that under the influence of gravity, the conversion process could become unstable, leading to fast detonation. A relativistic calculation was done (Cho et al. 1996) to determine the nature of the conversion process, employing the energy–momentum conservation and baryon number conservation (also known as the Rankine–Hugoniot condition). In most of this calculation, the dynamical evolution of the shock wave was not considered, and the combustion process was categorized by employing the energy–momentum and baryon number conservation across the shock discontinuity.

Very simplistic studies of dynamics of the PT were carried out first by Lin et al. (2006). They assumed that the PT occurs instantaneously at the star core that gets converted. The outer matter remains in the nuclear phase; a mixed phase region separates the inner and outer surface. The QM core has a radius of about 5 km. The typical timescale for the conversion to happen is of the order of 0.5 ms assuming that the sound speed in the star quark core is about 0.3–0.5c. They also carried out a simplistic calculation of the GW emission, due to such a PT induced collapse. Employing such a model in Newtonian hydrodynamics, they predicted that the gravitational strain h is of the order of (3–15) × 10⁻²³ for a source at a distance of 10 Mpc. The energy carried by such GW is in the range (0.3–2.8) × 10⁵¹ erg. Abdikamalov et al. (2009) improved similar calculation for an axis-symmetric star, taking GR into account. They also introduced a finite timescale for the PT instead of instantaneous PT. The GW strain h for their calculation lies in the range (0.1–2.4) × 10⁻²³ for different models considered. The energy output of the GW during the first 50 ms of the evolution is in the range between 10⁴⁸ and 10⁵⁰ erg. They inferred that the detection possibility of such GW signals is quite low for first-generation or, in that sense, second-generation detectors, but entirely possible for third-generation detectors, especially for the GW signals coming from Virgo cluster.

Herzog & Ropke (2011) did another calculation of 3D hydrodynamic simulation of the combustion of NS to QS. In their work, they used more realistic EoS both for the quark and hadronic matter, instead of the polytropic EoS used by previous authors. They assume that PT proceeds via slow deflagration process from the center to the surface, converting NM to QM. However, the combustion velocity or the mean propagation velocity of the flame is described by the turbulent burning speed, which results from instabilities of the flaming front. Their hydrodynamic front stops somewhere inside the star, and the outer layers are still in the hadronic phase.

In this paper, we address the dynamical PT scenario, where a shock propagation induces the PT. The PT from NS to QS is presumably the first order PT. The phase transformation is usually assumed to begin near the center of a star, where the density is at least few times the critical density (Abdikamalov et al. 2009). Several processes can trigger the PT: spin-down of star, accretion of matter on the stellar surface (Alcock et al. 1986), or just cooling. Such PT could happen in newly born stars or could occur even in old NS that is accreting matter from its companion. In our study we assume a sudden density rise at the center of the star. This can happen when there is sudden spin-up (Chubarian et al. 2000) or spin-down (Glendenning et al. 1997; Spyrou & Stergioulas 2002) of the star. Due to sudden change in their spin and to conserve angular momentum, the density at the center of a star changes abruptly, thereby creating a density discontinuity near the star core. The formation of millisecond pulsars is thought to be a process where an old NS gains angular momentum by accretion from a binary companion, and the amount of spin-up is related directly to the amount of mass accreted (Burderi et al. 1999). If the process is true, then we would have a considerable population of rapidly rotating NSs that are quite massive (and hence have high central densities). Recent observation also seems to confirm this assumption, as there are quite a significant number of millisecond pulsars, which indeed are massive (Freire et al. 2008) and thus can be the candidates of those stars that are undergoing (or have undergone) a PT of the type under discussion.

We study the spread of the shock from the center to the surface, and the point of the contact discontinuity that propagates with the shock is the point where the NM deconfines to QM. The shock propagation is solved using the hydrodynamic Euler equations. The EoS of the NM and QM is obtained by fitting piecewise polytrope to the actual EoS. The density fluctuation at the center starts the whole process, and there are pressure and density discontinuity at the core of the star. Depending on the density of the NS (precisely that of the nuclear matter) at any particular radius, the combusted QM density at that radius is obtained by employing the energy–momentum and baryon number conservation (see Mallick & Singh 2017 for a detailed discussion). The calculated QM density serves as the initial condition for the hydrodynamic differential equations. With time, the shock wave propagates (the discontinuity) outward to the periphery of the star, converting NM (ahead of the shock) to QM (behind the shock). The hydrodynamic equations govern the evolution of the shock wave along the star. We ensure that the point of contact discontinuity is the point where the actual PT from NM to QM is taking place. As the shock propagates outward to the low-density region, the discontinuity reduces. Also, if there is a formation of the mixed phase, it will dissolve this sharp discontinuity. The smoothing of the discontinuity will depend on how much time is needed for the mixed phase region to grow. It is quite possible that the shock can travel a considerable distance before it smooths out. Therefore, after propagating some distance in the star, the PT will stop. Usually this happens somewhere inside the star, depending on the EoS we choose. We thereby have a star that has a quark core and a hadronic outer layer. Such stars are known as hybrid stars (HSs). During the conversion process, the stars suffer significant structural transformation and contracts considerably, which gives rise to GWs.

The paper is arranged as follows. In Section 2 we will be discussing our numerical techniques and the codes used to generate our results. In Section 3 we present the EoS that were used to model the compact stars, while in Section 4 we show our results for the PT of NS to QS. Finally, in Section 5 we summarize our findings and conclude them.

In our present work, we have used the open source GR1D code (O'Connor & Ott 2009), which is a spherically symmetric general-relativistic hydrodynamics code. In this section we summarize the formulation of curvature and hydrodynamics equations on which the GR1D simulation is based. After presenting the essential equations, we also describe the numerical methods that have been used in GR1D code.

We have made changes to the GR1D code to suit our needs, such as to include the profiles obtained by solving the TOV equation as initial configuration, to locate the position of shock during the time evolution, and to permit a change in the EoS in the region surpassed by shock to mimic PT.

2.1. Curvature Equations

GR1D uses spherically symmetric metric in RGPS coordinates (Gourgoulhon 1991; Romero et al. 1996). The metric is given by

$$\begin{eqnarray}&&{{ds}}^{2}=-\alpha {(r,t)}^{2}{{dt}}^{2}+X{(r,t)}^{2}{{dr}}^{2}+{r}^{2}d{\theta }^{2}+{r}^{2}{\sin }^{2}\theta d{\phi }^{2},\end{eqnarray} \tag{ 1 }$$

where $\alpha (r,t)=\exp [{\rm{\Phi }}(r,t)]$, with Φ(r, t) being the metric potential and $X(r,t)={\left(1-\tfrac{2m(r,t)}{r}\right)}^{-1/2}$ with m(r, t) being the gravitational mass at the radial distance r.

We assume the matter to be a perfect fluid, described by a mass current density J^μ = ρu^μ and stress energy tensor ${T}^{\mu \nu }=\rho {{hu}}^{\mu }{u}^{\nu }+{g}^{\mu \nu }P$. In the equations, ρ is the rest mass density, P is the fluid pressure, and h = 1 + e + P/ρ is the specific enthalpy, with e being specific internal energy. ${u}^{\mu }=(W/\alpha ,{Wv}/X,0,0)$ is the fluid 4-velocity. $W=\sqrt{\tfrac{1}{1-{v}^{2}}}$ denotes the Lorentz factor and v is the physical radial velocity.

The expression for gravitational mass can be derived from Hamiltonian constraint equation, and the expression for metric potential can be derived from the momentum constraints. It comes out to be

$$\begin{eqnarray}&&m(r,t)=4\pi {\int }_{0}^{r}(\rho {{hW}}^{2}-P)r{{\prime} }^{2}{dr}^{\prime} \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}{\rm{\Phi }}(r,t) & = & {\int }_{0}^{r}{X}^{2}\left[\displaystyle \frac{m(r^{\prime} ,t)}{r{{\prime} }^{2}}+4\pi r^{\prime} (\rho {{hW}}^{2}{v}^{2}+P)\right]\\ & & \times {dr}^{\prime} +{{\rm{\Phi }}}_{0}.\end{eqnarray} \tag{ 3 }$$

Φ₀ is obtained by matching the Φ(r, t) at the star surface (r = R_*) to the metric potential of the Schwarzschild metric given by

$$\begin{eqnarray}&&{\rm{\Phi }}({R}_{* },t)=\mathrm{ln}[\alpha ({R}_{* },t)]=\displaystyle \frac{1}{2}\mathrm{ln}\left[1-\displaystyle \frac{2m({R}_{* },t)}{{R}_{* }}\right].\end{eqnarray} \tag{ 4 }$$

2.2. Evolution Equations

The GR hydrodynamics equation used in GR1D is based on flux conservative Valencia formulation (Banyuls et al. 1997; Font et al. 2000; Font 2008) with modifications for spherically symmetric flow (Romero et al. 1996). The evolution equations for matter current density are obtained from the continuity equation

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }{J}^{\mu }=0\end{eqnarray} \tag{ 5 }$$

and for the matter fields from the local conservation laws for stress energy tensor

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{\mu }{T}^{\mu \nu }=0.\end{eqnarray} \tag{ 6 }$$

In the coordinate frame of GR1D where u^μ = (W/α, Wv/X, 0, 0), the evolution of rest mass density is given by

$$\begin{eqnarray}&&{\partial }_{t}(D)+\displaystyle \frac{1}{{r}^{2}}{\partial }_{r}\left(\displaystyle \frac{\alpha {r}^{2}}{X}{Dv}\right)=0,\end{eqnarray} \tag{ 7 }$$

where D is the conserved variable given by D = XρW.

The momentum evolution is given by

$$\begin{eqnarray}&&{\partial }_{t}({S}^{r})+\displaystyle \frac{1}{{r}^{2}}{\partial }_{r}\left[\displaystyle \frac{\alpha {r}^{2}}{X}({S}^{r}v+P)\right]=\alpha X[({S}^{r}v-\tau -D)\\ &&\,\left.\times \left(8\pi {rP}+\displaystyle \frac{m}{{r}^{2}}\right)+\displaystyle \frac{{Pm}}{{r}^{2}}+\displaystyle \frac{2P}{{X}^{2}r}\right]+\alpha W({{vQ}}_{E}^{0}+{Q}_{M}^{0}),\end{eqnarray} \tag{ 8 }$$

where S^r is the conserved variable given by S^r = ρhW²v and conserved variable τ = ρhW² − P − D.

The energy evolution equation is given by

$$\begin{eqnarray}&&{\partial }_{t}(\tau )+\displaystyle \frac{1}{{r}^{2}}{\partial }_{r}\left[\displaystyle \frac{\alpha {r}^{2}}{X}({S}^{r}-{vD})\right]=\alpha W({Q}_{E}^{0}+{{vQ}}_{M}^{0}),\end{eqnarray} \tag{ 9 }$$

where ${Q}_{E}^{0}$ and ${Q}_{M}^{0}$ denotes the energy and momentum source terms. The conserved variables are function of primitive variables ρ, e, v, and P.

This set of evolution equations can be written as

$$\begin{eqnarray}&&{\partial }_{t}{\boldsymbol{U}}+\displaystyle \frac{1}{{r}^{2}}{\partial }_{r}\left[\displaystyle \frac{\alpha {r}^{2}}{X}{\boldsymbol{F}}\right]={\boldsymbol{S}},\end{eqnarray} \tag{ 10 }$$

where ${\boldsymbol{U}}$ is the set of conserved variables, ${\boldsymbol{F}}$ is their flux vector, and the ${\boldsymbol{S}}$ vector contains gravitational and geometric sources. They can expressed as

$$\begin{eqnarray}&&{\boldsymbol{U}}=[D,{S}^{r},\tau ]\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{\boldsymbol{F}}=[{Dv},{S}^{r}v+P,{S}^{r}-{Dv}]\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}{\boldsymbol{S}} & = & \phantom{\rule{}{3.25ex}}[0,\alpha X({S}^{r}v-\tau -D)\left(8\pi {rP}+\displaystyle \frac{m}{{r}^{2}}\right)+\alpha {PX}\displaystyle \frac{m}{{r}^{2}}\\ & & +\displaystyle \frac{2\alpha P}{{Xr}}+\alpha W({{vQ}}_{E}^{0}+{Q}_{M}^{0}),\alpha W({Q}_{E}^{0}+{{vQ}}_{M}^{0})\phantom{\rule{}{3.25ex}}].\end{eqnarray} \tag{ 13 }$$

2.3. Numerical Method

To develop the numerical code, we first discretized the space. Then, applying the method of lines (MoL, Hyman 1976), and by using the standard second-order or third-order Runge–Kutta, time integration of conserved variables is carried out. The spatial discretization is done by finite volume approach (Romero et al. 1996; Font 2008). The variables are defined at cell centers i and are interpolated at cell interfaces. At the cell interfaces, the inter-cell fluxes are evaluated. To interpolate, the third-order piecewise-parabolic method (PPM; Colella & Woodward 1984) is used to smoothen parts of flow. The primitive variables are interpolated, and then corresponding conserved variables are found at the cell interfaces. Also, piecewise-constant reconstruction is used, and to avoid oscillations near origin in the innermost three to five zones specifically, the piecewise linear rebuilding (van Leer 1977) is used.

After the conserved variables are evaluated at the cell interfaces, physical interface fluxes ${{\boldsymbol{F}}}_{i+1/2}$ are found with the HLLE Reimann solver. The right-hand side flux update term for ${{\boldsymbol{U}}}_{i}$ comes out to be

$$\begin{eqnarray}{\mathrm{RHS}}_{i} & = & -\displaystyle \frac{1}{{r}_{i}^{2}{\rm{\Delta }}{r}_{i}}\left[\displaystyle \frac{{\alpha }_{i+1/2}{r}_{i+1/2}^{2}}{{X}_{i+1/2}}{{\boldsymbol{F}}}_{i+1/2}\right.\\ & & \left.-\displaystyle \frac{{\alpha }_{i-1/2}{r}_{i-1/2}^{2}}{{X}_{i-1/2}}{{\boldsymbol{F}}}_{i-1/2}\right].\end{eqnarray} \tag{ 14 }$$

Gravitational and geometrical sources are not taken into consideration in flux evaluation but are coupled with the MoL integration. Once the conserved variables are updated, the primitive variables needed for the next time step need to be found. Since the primitive variables cannot be written algebraically in terms of conserved variables (Font et al. 2000), the primitive variables are found iteratively with an initial guess using P_old (from previous time step) with help of the expressions

$$\begin{eqnarray}&&v=\displaystyle \frac{{S}^{r}}{\tau +D+{P}_{\mathrm{old}}}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\rho =\displaystyle \frac{D}{{XW}}\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&e=\displaystyle \frac{\tau +D+{P}_{\mathrm{old}}(1-{W}^{2})}{\rho {W}^{2}}-1,\end{eqnarray} \tag{ 17 }$$

where the X is calculated from the conserved variables and W is calculated using the value of v.

The EoS is used to find the pressure (with the density obtained from above expression). This process is iterated using the Newton–Raphson method till the pressure value between the successive iteration step differs by 10⁻¹⁰. Thus, in each time step, we find the v, ρ, e at various cells (space discretized blocks) and finally the P using the iteration process.

2.4. Modifications in GR1D

NSs is assumed to contain mainly neutrons, protons, electrons with other baryons, and leptons in a small fraction. The sigma, omega, and rho mesons also form the constituent of NSs. In literature, different authors have considered different EoS, depending upon the constituent particles they would like to discuss.

3.1. PLZ EOS

The PLZ parameter setting is chosen for describing the properties of matter in the nucleonic form. The corresponding Lagrangian is given in the following form (Boguta & Bodmer 1977; Serot & Walecka 1986; Glendenning & Moszkowski 1991; ${\hslash }=c=1$):

$$\begin{eqnarray}{{ \mathcal L }}_{H} & = & \displaystyle \sum _{n}{\bar{\psi }}_{n}\left[{\gamma }_{\mu }(i{\partial }^{\mu }-{g}_{\omega n}{\omega }^{\mu }-\displaystyle \frac{1}{2}{g}_{\rho n}{\boldsymbol{\tau }}.{{\boldsymbol{\rho }}}^{\mu })\right.\\ & & -({m}_{n}-{g}_{\sigma n}\sigma )\phantom{\rule{}{3.25ex}}]{\psi }_{n}\\ & & +\displaystyle \frac{1}{2}({\partial }_{\mu }\sigma {\partial }^{\mu }\sigma -{m}_{\sigma }^{2}{\sigma }^{2})-\displaystyle \frac{1}{3}b{\sigma }^{3}-\displaystyle \frac{1}{4}c{\sigma }^{4}\\ & & -\displaystyle \frac{1}{4}{\omega }_{\mu \nu }{\omega }^{\mu \nu }\\ & & +\displaystyle \frac{1}{2}{m}_{\omega }^{2}{\omega }_{\mu }{\omega }^{\mu }-\displaystyle \frac{1}{4}{{\boldsymbol{\rho }}}_{\mu \nu }.{{\boldsymbol{\rho }}}^{\mu \nu }+\displaystyle \frac{1}{2}{m}_{\rho }^{2}{{\boldsymbol{\rho }}}_{\mu }.{{\boldsymbol{\rho }}}^{\mu }\\ & & +\displaystyle \sum _{l}{\bar{\psi }}_{l}[i{\gamma }_{\mu }{\partial }^{\mu }-{m}_{l}]{\psi }_{l}.\end{eqnarray} \tag{ 18 }$$

The baryonic matter is assumed to contain nucleons (n) and leptons (l = e^±, μ^±). The leptons are assumed to be non-interacting, whereas the baryons interact with the scalar σ, the isoscalar-vector ω_μ, and the isovector–vector ρ_μ mesons. The fundamental properties of NM and those of finite nuclei are used to fit the adjustable parameter of the model. We can use different fitting parameters to generate different EoS. In the present problem, we have used a PLZ Reinhard (1988) parameter set to construct our star.

3.2. Bag Model EOS

To describe the QM, we use a simplistic MIT bag model (Chodos et al. 1974). Such a simple model is enough for our calculation, because we are not concerned about the microscopic nature of the matter but instead the macroscopic properties—namely the mass, radius, pressure, and so on. Recent calculation has shown that inclusion of the quark interaction in this model makes it possible to satisfy the present heavy pulsar mass bound. The potential for this model can be described as

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{Q}=\displaystyle \sum _{i}{{\rm{\Omega }}}_{i}+\displaystyle \frac{{\mu }^{4}}{108{\pi }^{2}}(1-{a}_{4})+B,\end{eqnarray} \tag{ 19 }$$

where Ω_i is the potential for species i and B is the bag constant. The model has only quarks and leptons. μ signifies the baryon chemical potential, and a₄ is the interaction parameter among the quarks, varied between 1 (no interaction) and 0 (full interaction).

Due to the first order PT happening at the core of the star, nucleons are deconfined to mostly up (u) and down (d) quarks. The density at the core of stars is such that it contains mostly u and d quarks. The masses of the u and d quarks are 5 and 10 MeV, respectively. We choose the values of B^1/4 = 140 MeV and a₄ = 0.5. Later, the excess of d quarks may undergo weak decay to form s quarks; however, with our choice of parameter setting, this two-flavor QM is also quite stable with respect to the nuclear matter.

3.3. Polytropic EoS

Employing such complex EoS in the hydrodynamic calculation is a herculean task. Usually the hydrodynamic calculation is performed using simple ideal gas EoS. However, they are not suitable for modeling NSs. For some astrophysics calculation, polytropic EoS are often used. In our estimation, we will use such polytropes. Dealing with the ideal Fermi gas model containing only electrons or neutrons, it comes out that the degeneracy pressure P follows a power-law dependence on mass density ρ, given by P = Kρ^γ. This relation between pressure and density is called a polytropic EoS. Here K and γ are constants. In the literature, polytropic EoS has been used to construct the profiles of a white dwarf or NS with proper choice of K and γ (Shapiro & Teukolsky 1983). Though polytropic EoS is very simplistic in some situations, it can give a good description of the star or parts of the star. To replicate a real EoS profile (obtained by using relativistic mean-field calculation), it is better suited to map the EoS with more than one polytrope, with different γ. Finding the right ones and constructing a proper, thermodynamic consistent, EoS with them is difficult, and also the result will be not unique. However, another way of constructing polytropic EoS was given by Bonazzola et al. (1993). The polytrope can be described as

$$\begin{eqnarray}\epsilon (n) & = & {m}_{B}n+\displaystyle \frac{k{\epsilon }_{0}}{\gamma -1}{\left(\displaystyle \frac{n}{{n}_{0}}\right)}^{\gamma }\\ p(n) & = & k{\epsilon }_{0}{\left(\displaystyle \frac{n}{{n}_{0}}\right)}^{\gamma },\end{eqnarray} \tag{ 20 }$$

where n is the baryon number density, is the energy density, and p is the pressure. k and γ are dimensionless parameters, m_B is the mean baryon mass m_B = 931.2 MeV, and n₀ and ₀ are arbitrary number and energy densities. For NS matter, values around nuclear density are appropriate, so, for example, n₀ = 0.14–1.7 fm⁻³ and ${\epsilon }_{0}={m}_{B}{n}_{0}\mathrm{MeV}\,{\mathrm{fm}}^{3}$. With those values, suitable k are around 0.05 and γ 2.45.

The EoS profiles are plotted in Figure 1(a). We see that PLZ is quite stiffer than the quark EoS. All the EoS can very well be described with more than one polytrope. Fitting them with a single polytrope is a bit difficult for the whole density range, but this can be done for a specific density range for which we are interested. The problem with fitting the nuclear EoS with more than one polytrope is that during the actual hydrodynamic simulation, the points where we match the two polytropes give rise to substantial numerical noise. It can interfere with the real shock wave, and possibly cloud some vital physics. Therefore, in our work, we have used a single polytrope to describe nuclear and QM. In Figure 1 we have also shown the fitted polytropes, and we find that they match quite well with the actual nuclear matter EoS for a broad density range. The polytropic EoS can generate 2.32 solar mass NS (similar to that of PLZ EoS), with γ = 2452, n₀ = 0.148, and k = 0.05.

Zoom In Zoom Out Reset image size

In Figure 1(b) we have plotted an M-R curve for the above discussed EoS. The PLZ EoS generates stars with a maximum mass of around 2.32 M_⊙ with a radius of 11.5 km. The fitted polytrope maximum mass and radius are reasonably close to that. The M-R sequence with QM EoS generates stars with a maximum mass close to 2 M_⊙, however with a much smaller radius of 10.5 km. Such stars are composed entirely of QM and are not altogether stable against NM. Therefore, a pure QS may not be entirely stable, and its outer region may convert into NM. Such conditions can give rise to stars that have a quark core surrounded by hadronic outer layers known as HS. In our calculation, we will use the quark (two-flavor) and PLZ EoS for the actual PT scenario.

In solving the hydrodynamic equation, we have only used the polytropic ones. Solving the hydrodynamic equations with nuclear and even simple quark EoS is quite a difficult task, where we may have to change the whole formulation. However, the ultimate physics does not change when we are studying the overall PT proceeding as a shock front and looking at it macroscopically. We are mimicking the general nature of pressure and energy as a polytrope, and macroscopically it remains the same. The microscopic nature may also play an important role; however, for the moment we are more interested in the macroscopic behavior.

The PT is brought about by a sudden density fluctuation at the center. Somewhere near the core (for our case we are taking it to be at about 1 km), the density fluctuation creates a shock like discontinuity. On the right side of the collapse, there is unburnt NM, whereas on the left side there is burnt QM. The shock then propagates outward to the surface with time. We study the shock propagation by solving the hydrodynamic equations, and find the spatial and temporal evolution of the collapse. As the shock propagates outward, it converts unburnt NM burn to QM, thereby enabling a PT. The initial conditions of the shock are fixed from the jump conditions (energy–momentum and baryon number conservation) that were previously discussed by us in Mallick and Singh (2017) and Mishustin et al. (2015).

We have modified a well-established hydrodynamic GR1D code for our calculation of astrophysical PT from NS to QS. The stepping stone of this problem is the old Sod shock tube problem. Therefore, as an initial check, we first show the results of the Sod shock tube problem. The problem is solved for an ideal gas EoS given by

$$\begin{eqnarray}&&p=(\gamma -1)\rho e,\end{eqnarray} \tag{ 21 }$$

where p is the pressure, γ is the adiabatic index, and e is the internal energy density. The problem is solved in the usual one dimension, with a discontinuity separating the two sides. The discontinuity is kept at x = 0.5, and the γ is taken to be 1.6. On the left side of the discontinuity, the fluid is at higher pressure and density, whereas on the right side the fluid is at much lower pressure and density. With time, the discontinuity propagates to the right, giving rise to shock and rarefaction waves simultaneously. This Sod shock tube problem is one of the oldest test problems for the Riemann problem where we solve the Euler's equation. The left and right side of quantities are

$$\begin{eqnarray}&&{\rho }_{l}=10\hspace{0.2cm}{p}_{l}=13.2\hspace{0.2cm}{v}_{l}=0,\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{\rho }_{r}=1\hspace{0.2cm}{p}_{r}=0\hspace{0.2cm}{v}_{r}=0.\end{eqnarray} \tag{ 23 }$$

In Figure 2, we show the time evolution of the density and pressure. The evolution of the shock from the initial configuration is shown for two non-zero time instances, t = 0.3 μ and t = 0.8 μ. As expected, we find that the shock discontinuity has propagated in the right, whereas the rarefaction wave has propagated to the left. Such behavior is found for both the density and pressure curves. The velocity of the shocked matter is shown in Figure 3, where we find that only the shocked matter (matter affected by the shock) and rarefaction wave have attained some non-zero velocity. The unshocked matter velocity remains zero. All the curves match with the usual shock tube problem.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the propagation of a pure shock wave in an NS, we first solve the TOV equation for a given EoS with a particular value of central density. The solution of the TOV provides the pressure and density as a function of the radius of the star. Using this as an input, the hydrodynamic equations are solved with the given pressure and density function. A relativistic code is needed to solve the problem because the EoS of the NSs are usually relativistic, and it also avoids the superluminal velocity problem. The problem is solved for a spherically symmetric star, and the discontinuity is kept near the center of the star, where the density is maximum. In our problem we have assumed a reflecting boundary condition at both ends of the star. If the initial discontinuity is kept very near the center, the rarefaction wave going backward will get reflected and interfere with the forward propagating shock wave and may hide some important physics that we want to study. However, we will show that keeping it at a reasonable distance does not change the qualitative nature of our results. We have used a polytrope that mimics the PLZ EoS model and have constructed a star with a central density of 2.4 times nuclear density, which yields a star of mass 1.79 M_⊙ and radius of 13.2 km. In Figure 4 we plot the density (Figure 4(a)) and pressure (Figure 4(b)) as a function of radial distance. The density and pressure fall smoothly from the center to the surface for an unshocked star. We have kept the discontinuity at a distance of 1.05 km from the center of the star. We have shown five time slice (at four different times, including initial configuration) plots of the density and pressure evolution as a function of the radial distance. The initial discontinuity is given for the density, and the pressure discontinuity is obtained from the polytrope. The initial matter velocities on either side of the shock are kept zero. We see that with time the discontinuity proceeds toward the periphery of the star from the high-density to the low-density region and its strength gets reduced. Pressure plot also suggests similar behavior.

As the shock wave travels outward, it carries a lot of material with it, which thereby reduces the central density of the star. Such a feature is clear from the figure. If the shock can propagate through the whole star, and at the surface if it still has considerable strength, then it would expel some of the matter from the star interior, similar to that of supernovae explosion. However, gravity should play a huge role toward the restructuring of the whole star as the shock propagates. The effect of gravity has not been taken into account in our calculation, as it would be a very challenging task. The role of gravity and shock initial strength would play a decisive role toward restructuring of the star while the shock wave propagates.

In Figure 5(a) we plot the velocity of the shocked and unshocked matter. At an initial time before the shock propagates, all the matter velocities are taken to be zero. As the collapse spreads, the speed of both shocked and unshocked matter gets changed. With time, the velocity of the matter grows as the shock propagates outward. The velocity profile also takes the shape of discontinuity, and the maximum speed is obtained at the shock boundary. The velocity of the unshocked matter also gets changed as the density profile of the star changes due to the propagation of the discontinuity (because we are dealing with a very dense closed system). The finite velocity at the shock surface indicates that the shock would not stop at the star boundary and is likely to be expelled out of the star. Gravity plays a decisive role in determining how much matter, if any, should be ejected from the star. The finite velocity of the shock at the boundary indicates that the density and pressure would also rise to a limited value at the surface. However, it is not huge. In our problem, we, thus have fixed the pressure to be zero at the surface.

Zoom In Zoom Out Reset image size

The location of the shock discontinuity at each instance is known from the density and pressure plots. Differentiating the shock location numerically with time, we calculate the velocity of the shock. In Figure 5(b) we plot the shock velocity with time. Initially, the shock velocity starts with 0.68 times the light speed (c) and rises with time. The final shock velocity is about 0.94 c. Such high shock propagation indicated that the shock could propagate the entire star in about 50 microseconds (μs) .

We have constrained that the pressure should be zero at the surface. However, such a strict boundary condition will not affect the shock propagation much, as has been shown in Figure 6. In Figure 6(a) we plot the pressure evolution as the shock propagates outward and we see that the pressure at the boundary is not zero. This is because of the boundary condition and also of the fact that we are dealing with the matter at very high density. However, when we plot the shock location (r_s) as a function of time, we find that shock location does not change much. This also means that the shock velocity will also not differ much.

Zoom In Zoom Out Reset image size

The two-dimensional plots of the density evolution are shown in Figure 7. These 2D plots are obtained by rotating the 1D plots by 360° about the radial plane. The color coordinates show the intensity or the value of the density. The red color is for maximum density, and the purple for minimum density. As the shock propagates, the discontinuity of the density follows it.

Till now we have calculated the shock propagation in a star without involving PT. We currently study a scenario where a shock propagation is such that it induces a PT. This can be achieved numerically if we ensure that the shocked material has the properties of QM and the unshocked material has nuclear properties, and they obey the energy–momentum and baryon number conservation at the shock boundary. This can be done once we locate the position of the shock discontinuity, which is done by our code. The code finds the shock boundary separating the two phases and then calculates the matter properties on either side of the shock, following the conservation condition. As the shock wave propagates, it converts NM to QM; however, after some distance as the density decreases, the shock wave starts to lose its strength. Also, if there is a formation of the mixed phase, it will dissolve this sharp discontinuity. The smoothing of the discontinuity will depend on how much time is needed for the mixed phase region to grow, whereas our propagation of the shock is very fast. It is quite possible that it can travel a considerable distance before it smooths out. Therefore, after propagating some distance in the star, the PT will stop and eventually the shock wave may die out. We should mention that we have not done the exact calculation on how long it would take for the sharp discontinuity to dissolve, due to the formation of the mixed phase. Such a calculation may be important but involves much complexity. For a simplicity, we do not consider such a case in our study. The lowest density of the QM EoS that we have used is about 1.8 times the nuclear density, which corresponds to a distance of about 7 km in the star. Therefore, after a distance of 7 km, the shock evolves without bringing about a PT. Therefore we have an HS, a star with outer NM, and a quark core. We have kept the initial shock discontinuity at about 1.05 km from the center of the star. As soon as the shock starts to propagate, the post-shock quantities are replaced by the quark EoS, whereas the pre-shock quantiles retain their values. With time, the shock propagates outward, keeping such criteria intact. Therefore, we obtain a PT from nuclear matter to QM.

In Figure 8 we have shown the evolution of the density and pressure with time (displayed for five time instances). The density and pressure of this plot differ from the previous growth in the fact that there is a generation of sharp peak discontinuity at the shock boundary. The peak discontinuity is present as the shock front propagates outward. The PT from NM to QM continues till 7 km from the center of the star. To travel 7 km, the shock takes about 25.8 μs. However, beyond 7 km the shock wave only propagates through the star without converting NM to QM. From Figure 8 it is also clear that beyond 7 km, the peak suddenly disappears and only the shock discontinuity propagates.

Figure 9(a) shows the evolution of matter velocities with distance and time. As done in our earlier study, we set the matter velocity initially as zero. At a later time, the matter velocities peak up both in the shocked and unshocked matter. The peaking of speed is a manifestation of the fact that we are dealing with a dense closed system and also that the PT is a very fast process. The matter velocities for the shock-induced PT are faster than ordinary shocks. In Figure 9(b) we plot the speed of the shock with time. The initial shock velocity is about $0.7c$, and with time as the collapse propagates outward, the shock velocity rises and becomes 0.8 c at 7 km. As the PT stops in the star and only normal shock propagates, the velocity falls off suddenly (the dip around 25 μs). The velocity then grows very slowly and again reaches 0.94 c at the star boundary. Due to a sudden decrease in the shock velocity around 7 km, the time taken for the shock to reach the surface is slightly larger.

Zoom In Zoom Out Reset image size

In Figure 10 we plot the 2D plot of the density evolution for a shock-induced PT. The peak-like feature is also visible in this scenario. The lowering of the density at the center, as the shock propagates out, is also observed. As the front propagates outward, the sharp discontinuity disappears.

In our calculation, we have assumed the initial matter velocity of the shock to be zero. However, it may happen that the initial speed is not zero but a finite number. However, at the core of the star, the initial velocity can never be very enormous, and it would not change our results much.

In Figure 11 we show the evolution of ρ as a function of radial distance when the density discontinuity is kept at a distance of 400 m. From the figure, we can see that just after the formation of the discontinuity, the shock wave travels forward and the rarefaction wave travels back. As we have a reflecting condition at the boundary, the rarefaction wave gets reflected and starts to interfere with the forward progressing shock wave (the green dashed curve). The shock wave thereby starts losing its strength as its progress forward and after some time becomes very thin. This may not be suitable for our study. Therefore, we have kept our shock at such a distance that the rarefaction wave cannot interfere with it and the shock progresses smoothly.

Zoom In Zoom Out Reset image size

In Figure 12(a) we show the evolution of ρ as a function of radial distance when the discontinuity is a bit smaller. There is not an appreciable difference in the evolution of the front; however, with less discontinuity the front velocity is smaller, as expected (as shown in Figure 12(b)). The velocity of the shock decreases with a decrease in the initial discontinuity, and it takes a bit longer to travel through the whole star. For much lower discontinuity, it becomes difficult in our code to locate the point of PT and therefore the evolution of the PT. Hence we have done our studies with amplitudes of discontinuity that are 1.33 and 1.25 times the central density of the star.

Zoom In Zoom Out Reset image size

The combustion process of burning NM to QM is an exothermic process. Therefore, although the NM may have zero temperature, the burned QM may have a finite temperature. However, by such combustion, the star is not expected to heat up much, and the maximum temperature rise may be of the order of 10–20 MeV (Mishustin et al. 2015). In our calculation, we have calculated the QM EoS at a temperature of 15 MeV. As our QM EoS is quite stiff, such change in temperature does not change the QM EoS considerably. As can be seen from Figure 13, there is not much change in the overall density and/or pressure evolution of the shock, and the quantitative change in the PT timescale is also quite small. (It now takes 51.07 s to travel across the whole star). This is because the shock velocity for the finite temperature QM EoS is smaller than the zero temperature. This is expected because some amount of energy is going into heat energy, which manifests itself in the form of a finite temperature of the QM. As the available energy for the shock propagation is now less, the shock evolution is a bit slower.

Zoom In Zoom Out Reset image size

In this work, we have studied the dynamical evolution of the shock front as it travels from the core to the surface of a star. The shock front generates at the center of a star due to sudden density and pressure fluctuation. The density fluctuation can be initiated by a number of reasons, ranging from simple cooling to sudden braking of the star. Once the shock front generates at such high density, it propagates outward, converting NM to QM. We have studied this PT from NM to QM with our hydrodynamic code. We assumed that the shock generates near the center of the star at a distance around 1 km. In our study, we have employed PLZ EoS to model the nuclear matter and simple MIT bag model with strong coupling to describe the QM. We have assumed that the PT or the burning of matter happens at the shock boundary. In this work, we have studied the PT using the relativistic hydrodynamic equations and have not taken gravity into account.

We have modeled the PT of NS by modifying a Sod-type problem. The problem has been modified to incorporate the star profile, such as the density and pressure variation, as a function of radius (instead of "x"). Next, we have located the position of the shock discontinuity, and by our assumption, this is the place where PT happens. Once we have located the shock discontinuity, we have ensured that the post-shock burnt matter should have a different EoS than the pre-shock unburnt matter obeying the conservation conditions. As time evolves and as the shock propagates outward, this is repeated for every small time step, thereby imitating the PT scenario.

It is difficult to employ real nuclear and QM EoS in the hydrodynamic equations to describe matter properties. One of the ways out is using piecewise polytropic EoS. However, the higher the number of piecewise polytropes, the higher the chance of numerical fluctuations. To avoid such difficulty, we have employed a single polytrope to describe matter properties. The single polytrope cannot precisely replicate all stars that are described by using nuclear or quark EoS. However, we made sure that the stars which would be studied for PT scenario are closely represented by the polytropic stars that we employ.

In the literature, we have seen that the shock properties are affected when the matter velocities on either side of the front get very large (close to that of c). However, for practical purposes at such high density, it is difficult for the matter to obtain such high velocity, and therefore we have kept the initial speed to be zero at either side of the front. We are assuming a spherically symmetric non-rotating star. Studying the shock propagation, we find that the matter velocity takes non-zero values at places where the discontinuity has still not arrived. Such non-zero speed is because of the boundary condition, and also because we are dealing with dense confined matter. However, this does not change the pressure and density profile much, and so we have kept the surface pressure as zero at all times.

In the actual scenario, the real pressure may not be zero as the shock reaches the surface, and thereby it may expel some matter into space. However, such a situation should be complemented with gravity, which must be taken into account when we study the PT. Gravity is the galvanizing force that would play a huge role while restructuring of the star takes place after the propagation of the shock front. In our simple picture, we have not considered gravity, and we wish to address such detailed context soon.

Studying the shock propagation, we find that the shock strength decreases as the shock propagates to the boundary. However, the velocity of the shock increases as it travels outward to the low-density region. The shock takes about 50 μs to travel through the entire star. However, when we are studying the PT brought about by a shock, we see a very narrow peak-like structure as the shock is generated. The sharp discontinuity smooths out with time (and space). The PT in NS takes about 50 μs. This result is quite remarkable and different from any other calculation done earlier. In most of the previous works, the time of PT is of the order of milliseconds. Such significant change in the shock propagation time can have considerable significance in determining the observational outcome of NSs, like GW and gamma-ray bursts. Such small PT time would generate a stronger GW, however, lasting for about 50 microseconds. This is different from any other GW signals coming from other processes like black hole or NS mergers, where the GW last for longer times. This signal would be accompanied by other gamma-ray and neutrino signals, which would come from the PT of NM to QM. The timing of the gamma-ray and neutrino signal would also differ from any other signal coming from another process, due to the time difference. If such signals could be detected, then we can conclude that the PT in astrophysical scenarios is a real event and there can be QSs along with NSs.

We are in the process of refining our calculation in various aspects, such as employing real nuclear matter EoS in our problem and also employing some piecewise polytropes and trying to reduce numerical fluctuation at the fitting polytrope boundary. We are also trying to incorporate gravity in our problem, which would galvanize the matter. By incorporating gravity, we can even have situation where as the shock propagates, the radius of the star changes, as there is PT from NM to QM. The emergence of Gibbs phase and the smoothing of the discontinuity with time is another aspect we wish to study. All such calculation are in our immediate agenda.

R.M. would like to thank SERB, the Government of India, for monetary support in the form of Ramanujan Fellowship and Early Career Research Award. R.M. and R.P. would like to thank IISER Bhopal for providing all the research and infrastructure facilities.