Hybrid Stars using the Quark-Meson Coupling and Proper Time NJL Models

Hybrid Stars using the Quark-Meson Coupling and Proper Time NJL Models

D. L. Whittenbury    H. H. Matevosyan    A. W. Thomas CSSM and ARC Centre of Excellence for Particle Physics at the Terascale,
School of Physical Sciences, University of Adelaide, Adelaide SA 5005, Australia
July 14, 2019
Abstract
Background

At high density deconfinement of hadronic matter may occur leading to quark matter. The immense densities reached in the inner core of massive neutron stars may be sufficient to facilitate the transition.

Purpose

To investigate a crossover transition between two phenomenological models which epitomise QCD in two different regimes, while incorporating the influence of quark degrees of freedom in both.

Method

We use the Hartree-Fock quark-meson coupling model and the proper time regularised three flavour NJL model to describe hadronic and quark matter, respectively. Hybrid equations of state are obtained by interpolating the energy density as a function of total baryonic density and calculating the pressure.

Results

Equations of state for hadronic, quark and hybrid matter and the resulting mass versus radius curves for hybrid stars are shown, as well as other relevant physical quantities such as species fractions and the speed of sound in matter.

Conclusions

The observations of massive neutron stars can certainly be explained within such a construction. However, the so-called thermodynamic correction arising from an interpolation method can have a considerable impact on the equation of state. The interpolation dependency of and physical meaning behind such corrections require further study.

I Introduction

It has long been thought that the densities reached inside the inner core of neutron stars may be sufficient to produce a phase transition from hadronic matter to deconfined quark matter forming hybrid stars. A transition to the chirally restored phase of QCD is also thought to occur at high density. It is unknown if these two transitions coincide and the form they may take. In understanding these transitions one would ideally like to use QCD directly, but this is currently too difficult. One typically resorts to using two phenomenological models which epitomise the key features of QCD in the two asymptotic regions of its phase diagram—one in the low density region modelling hadronic matter and the other modelling quark matter at intermediate to high density—and then construct a phase transition between the two. This means that the dissociation of hadrons does not occur naturally and is dependent on how we construct the transition.

Various models have been used to describe each phase and different constructions in characterising the process of deconfinement have been investigated, such as the Maxwell Rosenhauer and Staubo (1991); Klahn et al. (2007); Blaschke et al. (2008); Agrawal (2010); Blaschke et al. (2010); Bonanno and Sedrakian (2012); Lenzi and Lugones (2012); Logoteta et al. (2013); Plumari et al. (2013); Benić (2014); Klahn and Fischer (2015) and Gibbs Glendenning (1992); Schertler et al. (1999); Bentz et al. (2003); Lawley et al. (2006); Carroll et al. (2009); Grunfeld et al. (2008); Xu et al. (2010); Shao (2011); Logoteta et al. (2012a, b); Orsaria et al. (2013); Hell and Weise (2014); Miyatsu et al. (2015) constructions describing first order transitions, as well as interpolation/percolation constructions where the transition is taken to be of crossover type Asakawa and Hatsuda (1997); Masuda et al. (2013a, b); Hell and Weise (2014); Alvarez Castillo et al. (2014); Kojo et al. (2015). The calculated properties of hybrid stars are considerably influenced by the choice of models and the type of construction used to describe the transition. In this work we will investigate transitions of the second or crossover type. This kind of transition was, for example, recently investigated by Masuda et al Masuda et al. (2013a, b) using several hadronic models and the three momentum regularised NJL model. They concluded that massive hybrid stars () could be produced using a percolation picture provided the quark matter EoS was stiff enough and the transition occurred at moderately low density (), thus providing a possible reconciliation of exotic degrees of freedom with the observations of Demorest et al. Demorest et al. (2010) and Antoniadis et al. Antoniadis et al. (2013). We now examine this possibility using the Hartree–Fock QMC model Whittenbury et al. (2014) and both the proper–time and three momentum regularised versions of the NJL model.

In using these two models, quark degrees of freedom influence both regions, with the latter also exemplifying chiral symmetry breaking. With both models employing quark degrees of freedom it is hoped that they will be more reliable in the transition region. In this region, where hadrons and quarks are expected to coexist, it is likely that the quark substructure of hadrons would play an important role and their interaction with the external quarks strong. The QMC model has the advantage over models which employ point-like descriptions of hadrons by modelling the baryons as MIT bags. By incorporating the change of this internal structure in-medium it naturally incorporates many-body forces Guichon and Thomas (2004); Guichon et al. (2006), which are in fact crucial for nuclear saturation.

Within the QMC model the in-medium changes of the baryon masses are calculated through the bag equations and then parametrised as functions of the scalar field as

 M∗B=MB−wσBgσN¯σ+d2~wσB(gσN¯σ)2 , (1.1)

(where the weightings and simply allow the use of a unique coupling to nucleons ). Using this parametrisation and a corresponding density dependent coupling, we can solve for the equation of state in the same standard way as the Walecka model Serot and Walecka (1986), that is at the hadronic level. In this way the sub-structure of the baryons is entirely contained in the mass parametrisation. We use the parametrisation given in Ref. Guichon et al. (2008), which includes the effects of one gluon exchange.

The Nambu-Jona-Lasinio (NJL) model is an effective model of the strong interaction introduced in 1961 Nambu and Jona-Lasinio (1961a, b). Originally the model was formulated in terms of nucleons as a local effective interaction inspired by the BCS theory of superconductivity Bardeen et al. (1957a, b) (At the time of the model’s conception quarks were yet to be discovered.). Later it was redeveloped in terms of quarks to model low and intermediate energy QCD. A number of reviews are available on the NJL model and its applications  Vogl and Weise (1991); Klevansky (1992); Hatsuda and Kunihiro (1994); Ripka (1996, 1997); Buballa (2005). Its ability to model the breaking of chiral symmetry dynamically by spontaneously generating mass has made it very popular.

The NJL model assumes that at low energy scales the gluons acquire a large effective mass and can be integrated out to a good approximation, leaving a local (contact) four Fermi interaction between quarks. Upon integrating out gluons the local colour symmetry of QCD is reduced to a global symmetry and confinement is lost. This shortcoming of the NJL model will not be an issue when modelling quark matter at high density, where the quarks will be considered to be deconfined.

Moreover, the NJL model is not renormalizable. However, unrenormalizable theories are still useful and information can be obtained through the process of regularisation, whereby a cut-off is introduced setting the scale of the model. There are many regularisation schemes in common use. Here we will use Schwinger’s proper time regularisation (PTR) Klevansky (1992); Ripka (1997) and make comparisons to the simple, non-covariant three momentum regularisation (TMR).

In Ref. Whittenbury et al. (2014) we extended the QMC model by performing a Hartree–Fock calculation including the full vertex structure for the vector mesons. This extension only alters the exchange contribution, including not only the Dirac vector term, as was done in  Rikovska-Stone et al. (2007), but also the Pauli tensor term. These terms were already included within the QMC model by Krein et al. Krein et al. (1999) for symmetric nuclear matter and more recently by Miyatsu et al. Miyatsu et al. (2012). We generalised the work of Krein et al. Krein et al. (1999) by evaluating the full exchange terms for all octet baryons and adding them, as additional contributions, to the energy density. A consequence of this increased level of sophistication is that, if we insist on using the hyperon couplings predicted in the simple QMC model, with no coupling to the strange quarks, the hyperon is no longer bound.

The present line of research complements the recent work of Refs. Masuda et al. (2013a, b); Hell and Weise (2014); Alvarez Castillo et al. (2014); Kojo et al. (2015), which also considered deconfinement in neutron stars as a crossover rather than a bona-fide phase transition. The novelty of this work is that we incorporate the influence of quark degrees of freedom in the hadronic phase. Moreover, we use a covariantly regularised NJL model for modelling deconfined quark matter. This work also complements Ref. Miyatsu et al. (2015), which investigated hybrid star matter using a different variation of the Hartree-Fock QMC model and a bag model, where deconfinement was modelled as a first order transition using a Gibbs construction. The present version of the QMC model differs from Miyatsu et al. (2015) as we restrict the exchanged mesons to , , , and ; we use couplings as derived within the model and treat contact terms differently.

This paper is organised as follows. The Hartree-Fock QMC model is presented in Sec. II and then used in Sec. III to study hadronic matter in generalised beta-equilibrium. The proper time regularised NJL model is used to study three flavour quark matter in Sec. IV. In Sec. V, using these two models a hybrid equation of state is obtained assuming a faux crossover construction. In Sec. VI, we summarise the results obtained and draw our conclusions.

Ii Hadronic model with quark degrees of freedom: The QMC model

In this section we present the formalism used for hadronic matter in beta-equilibrium with leptons using the Hartree-Fock QMC model Whittenbury et al. (2014). We will briefly review the derivation of the equations given in Refs. Whittenbury et al. (2014); Whittenbury (2015) and the approximations used therein. All our parameters are fixed at saturation density in Symmetric (N=Z) Nuclear Matter. We then extrapolate the model to investigate high density matter in generalised beta-equilibrium (GBEM), which is relevant to neutron stars.

ii.1 The Lagrangian and Hamiltonian density

In our calculations we consider only the spin- octet baryons. These baryons interact via the exchange of mesons which couple directly to the quarks. The exchanged mesons included are the scalar-isoscalar (), vector-isoscalar (), vector-isovector (), and pseudo-vector-isovector () bosons. These mesons only couple with the light quarks by the phenomenological OZI rule. We include the full vertex structure for the vector mesons, that is, we include both the Dirac and Pauli terms.

The QMC Lagrangian density used in this work is given by a combination of baryon and meson components

 L=∑BLB+∑mLm , (2.1)

for the octet of baryons and selected mesons with the individual Lagrangian densities

 LB = ¯ΨB(iγμ∂μ−MB+gσB(σ)σ (2.2) −gωBγμωμ−fωB2MNσμν∂μων −gρBγμt⋅ρμ−fρB2MNσμνt⋅∂μρν −gA2fπχBBγμγ5τ⋅∂μπ)ΨB ,
 ∑mLm = 12(∂μσ∂μσ−m2σσ2)−14ΩμνΩμν+12m2ωωμωμ (2.3) −14Rμν⋅Rμν+12m2ρρμ⋅ρμ +12(∂μπ⋅∂μπ−m2ππ⋅π) ,

for which the vector meson field strength tensors are and . The , denote meson-baryon coupling constants for the mesons. The quantities in bold are vectors in isospin space with isospin matrices denoted by and isospin Pauli matrices by . For nucleons and cascade particles . The pion-baryon interaction used here is assumed to be described by an SU(3) invariant Lagrangian with the mixing parameter  Rikovska-Stone et al. (2007) from which the hyperon-pion coupling constants can be given in terms of the pion nucleon coupling, de Swart (1963); Rikovska-Stone et al. (2007); Massot and Chanfray (2008), where and are the axial vector coupling and the pion decay constant, respectively.

From the Lagrangian given in Eq. (2.1) we obtain through the Euler-Lagrange equations a system of coupled non-linear partial differential equations for the quantum fields. This is a difficult system of equations to solve and to make the problem tractable approximations are applied. Static, no sea and mean field approximations are typically used and are implemented here. However, when we consider the NJL model we include the Dirac sea of negative energy states.

From the Hamiltonian density the EoS of nuclear matter can be calculated. The Hamiltonian density in the static approximation is

 H=∫d3r⎧⎨⎩K +∑m∈{σ,ω,ρ,π}Hm⎫⎬⎭, (2.4)

where we have decomposed it into its baryon and meson components as

 K = ∑B¯ΨB[−i→γ⋅→∇+MB−gσB(σ)σ]ΨB, (2.5) Hσ = 12→∇σ⋅→∇σ+12m2σσ2, (2.6) Hω = ∑B¯ΨB[gωBγμωμ−fωB2MNσμi∂iωμ]ΨB (2.7) −12[→∇ωμ⋅→∇ωμ+(→∇⋅→ω)2+m2ωωμωμ], Hρ = ∑B¯ΨB[gρBγμt⋅ρ μ−fρB2MNσμit⋅∂iρ μ]ΨB (2.8) −12[→∇ρμ⋅→∇ρμ+(→∇⋅→ρ)2+m2ρρμ⋅ρμ], Hπ = −∑B¯ΨB[gA2fπχBBγ5τ⋅(→γ⋅→∇)π]ΨB (2.9) +12→∇π⋅→∇π+12m2ππ⋅π.

ii.2 Hartree-Fock approximation

To perform the Hartree-Fock approximation we follow Refs. Guichon et al. (2006); Rikovska-Stone et al. (2007); Massot and Chanfray (2008); Hu et al. (2010) by considering each meson field to be decomposed into two parts, a mean field part and a fluctuation part , such that and solve the equations of motion order by order. The fluctuation terms are to be considered small with respect to the mean field contribution, the exception to this being the and meson fluctuations. In this fashion, the meson equation of motion is decomposed according to

 (−→∇2 +m2σ)(¯σ+δσ) =−∂K∂σ=−∂K∂σ(¯σ)−δσ∂2K∂σ2(¯σ)−….

The terms on the r.h.s of Eq. (II.2) are expanded about their ground state expectation values (denoted by ). For the first term, we have

 ∂K∂¯σ≡∂K∂σ(¯σ) = ⟨∂K∂¯σ⟩+δ[∂K∂¯σ] (2.11) = (2.12)

and similarly for the second and so on. We now proceed to solve the meson equation of motion order by order. At the mean field or Hartree level we obtain

 (−→∇2+m2σ)¯σ=−⟨∂K∂¯σ⟩=∑B(−∂M∗B∂¯σ⟨¯ΨBΨB⟩) (2.13)

and at the Fock level

 (−→∇2+m2σ)δσ = −δσ[⟨∂2K∂¯σ2⟩+(∂2K∂¯σ2−⟨∂2K∂¯σ2⟩)] = −∂K∂¯σ+⟨∂K∂¯σ⟩−δσ∂2K∂¯σ2, (2.15)

where to this order ()

 ∂2K∂¯σ2−−−−→⟨∂2K∂¯σ2⟩. (2.16)

The fluctuation equation of motion becomes

 (−→∇2+m2σ)δσ=−∂K∂¯σ+⟨∂K∂¯σ⟩−δσ⟨∂2K∂¯σ2⟩. (2.17)

Eq. (2.17) can be re-expressed in terms of an in-medium -meson mass and the fluctuation of the scalar baryon current as

 (−→∇2+m∗ 2σ)δσ=∑B−∂M∗B∂¯σ(¯ΨBΨB−⟨¯ΨBΨB⟩), (2.18)

where

 m∗ 2σ=m2σ+⟨∂2K∂¯σ2⟩=m2σ+∑B∂2M∗B∂¯σ2⟨¯ΨBΨB⟩. (2.19)

This in-medium meson mass is only relevant to the fluctuating part and does not appear in the mean field portion of the meson’s equation of motion. This in-medium modification due to the baryons internal structure was included in Refs. Guichon et al. (2006); Rikovska-Stone et al. (2007); Massot and Chanfray (2008), but we will omit it here. We are neglecting this in-medium modification as we are approximating the Fock terms in the static approximation, omitting all other meson retardation effects and implementing a crude method of subtracting the contact terms that arise in the Fock terms. For these reasons it is reasonable to disregard it and use the free meson mass in the Fock term, thereby treating it in the same manner as the other mesons.

The expectation value of the field is given by

 ¯σ=−1m2σ⟨∂K∂¯σ⟩=−1m2σ∑B∂M∗B∂¯σ⟨¯ΨBΨB⟩, (2.20)

which is then determined numerically. Krein et al. Krein et al. (1999) also considered an additional correction involving the mean scalar field generated by the Fock terms. This can be done by considering the energy density as a functional and requiring it to be thermodynamically consistent, meaning that the total energy density, , is minimised with respect to . This amounts to Eq. (2.20) plus an additional term because of the dependence of the Fock contribution to the energy density on .

The fluctuation of the field can now be written in terms of the meson’s Green function as

 δσ(→r) = ∫d3r′ Δσ(→r−→r ′)(−∂K∂¯σ+⟨∂K∂¯σ⟩)(→r ′) (2.21) = ∑B∫d3r′d3q(2π)3 ei→q⋅(→r−→r ′)Δσ(→q) [−∂M∗B∂¯σ(¯ΨBΨB−⟨¯ΨBΨB⟩)(→r ′)].

By considering the meson fields decomposed into a mean field part and a fluctuating part it can be seen in Eq. (2.18) and (II.2), that this is related to a similar decomposition of the baryon currents. We introduce the following notation:

 ¯ΨB~ΓαBΨB = ⟨¯ΨB~ΓαBΨB⟩ = (2.24)

where denotes one of the meson-baryon interactions, appearing in the Lagrangian (Eq. (2.2)), written in terms of the appropriate Lorentz and isospin structures.

The analogous expressions for the remaining mesons follow in same manner. Here we summarise the equations of motion for all mesons decomposed into mean field and fluctuation components. For the meson mean fields we have

 ¯σ = −1m2σ∑B∂M∗B∂¯σ⟨¯ΨBΨB⟩ , (2.25a) ¯ω = 1m2ω∑BgωB⟨Ψ†BΨB⟩ , (2.25b) ¯ρ = 1m2ρ∑BgρB⟨Ψ†Bt3BΨB⟩ , (2.25c)

and . Each meson field fluctuation can be condensed to

 δϕ(→r)=∑B∫d3r′ Δϕ(→r−→r′)δ(¯ΨB~ΓϕBΨB)(→r′), (2.26)

where is the Yukawa propagator for the meson .

The decomposition of the meson fields also occurs in the Hamiltonian. For example, the baryon kinetic and the meson terms amount to

 ∫d3r [K+Hσ] (2.28) = ∫d3r[K(¯σ)+12δσ(∂K∂¯σ−⟨∂K∂¯σ⟩)+12m2σ¯σ2] = ∫d3r[∑B¯ΨB[−i→γ⋅→∇+MB−gσB(¯σ)¯σ]ΨB +12∑B∂M∗B∂¯σδσ(→r)δ(¯ΨBΨB)(→r)+12m2σ¯σ2].

The first term in the last line is the meson’s Fock term contribution to the energy. To proceed further, we must first explain how the in-medium Dirac equation for the baryons is solved in the Hartree-Fock approximation. This is presented in Sec. II.3. In the process of doing this, we Fourier transform to momentum space, where the energy density of nuclear matter is more easily evaluated.

ii.3 The in-medium Dirac equation

The in-medium Dirac equation, for a baryon in nuclear matter, can be written as

 (⧸p−Mi−Σi(p))ui(p,s)=0, (2.29)

where is the self-energy of the baryon. From parity conservation and translational, rotational and time reversal invariance, the self-energy can be decomposed into three scalar functions in the nuclear matter rest frame Serot and Walecka (1986),

 Σi(p)=Σsi(p)+γ0Σ0i(p)+→γ⋅→pΣvi(p). (2.30)

The functions , and are the scalar, temporal vector and spatial vector components of the self-energy, respectively.

If we introduce the following effective quantities

 M∗i(p) = Mi+Σsi(p), (2.31) E∗i(p) = Ei(p)+Σ0i(p)=√→p∗2+M∗2i, (2.32) →p∗ = →p(1+Σvi(p)), (2.33)

the Dirac equation can be written in a form which is formally equivalent to the Dirac equation in vacuum. Therefore, as in vacuum, the positive energy solution to the in-medium Dirac equation is

 (2.34)

where are Pauli spinors and we have used the normalisation convention for the spinor Serot and Walecka (1986).

From fully self-consistent calculations performed using QHD Serot and Walecka (1986); Bouyssy et al. (1987) and QMC Krein et al. (1999) models, it is known that the scalar and temporal vector self-energy components are approximately momentum independent and the spatial vector component is very small. Therefore, we proceed by carrying out the self-consistency approximately, as in Ref. Whittenbury (2015); Whittenbury et al. (2014), where we ignore these small contributions. The self-energy then has a form identical to the usual mean-field (Hartree) result and the small Fock corrections can be included by requiring thermodynamic consistency.

The no sea approximation is used, i.e., the negative energy states of the baryons are ignored. Therefore, the in-medium propagators for baryons propagating on-shell in the nuclear matter rest frame are given entirely by the Dirac portion of the baryon propagator Serot and Walecka (1986)

 Giαβ(p)=iπE∗i(p)(⧸p∗+M∗i)αβδ(p0−E(p))Θ(pF,i−|→p|), (2.35)

where is the Fermi momentum of baryon .

With the above approximations and definitions, the Fock contribution can be evaluated. For the meson Fock contribution to the energy density, we obtain

 ϵFσ = 12∑i(∂M∗i∂¯σ)2∫|→p|≤pF,id3p(2π)3∫|→p′|≤p′F,id3p′(2π)3 (2.36) Δσ(→q)Tr[(⧸p∗+M∗i)(⧸p′∗+M∗i)]4E∗i(p)E∗i(p′).

As can be seen in Eq. (2.36), there is an additional scalar dependence in this Fock term, which appears after explicit evaluation. A correction to the mean scalar field can easily be included numerically. This is a small contribution, as illustrated by the scenarios labelled “Fock ” in Ref. Whittenbury (2015); Whittenbury et al. (2014), and so we neglect it here. The Fock contributions of the other mesons are straightforwardly obtained in the same manner. They too have an additional dependence on the scalar mean field.

ii.4 Hartree-Fock equation of state

Here, to describe the baryon-meson interaction, we also introduce form factors, because of the extended nature of the baryons, by

 gαB−−−−→gαBFα(k2). (2.37)

The , , and form factors are all taken to have the dipole form with the same cut-off . We explored values of the cut-off mass in the range 0.9 – 1.3 GeV in Ref. Whittenbury et al. (2014). Clearly, these form factors are only of concern for the Fock terms as these allow for a finite momentum transfer, whereas Hartree contributions do not.

Within the QMC model, the hadronic energy density is the sum of the of the baryonic energy density in nuclear matter which is

 ϵB=2(2π)3∑B∫|→p|≤pFd3p √→p2+M∗2B, (2.38)

and the mesonic energy density . This can be divided into two parts, the Hartree and the Fock contribution. The total mesonic energy density is given by , where the Hartree and Fock components of the mesonic energy density are given respectively by

 (2.39)

where refers to the mean value of meson field and

 ϵF = (2.40) ∫|→p|≤pF∫|→p′|≤pF′d3pd3p′ ΞαBB′ ,

where . and , which arise from symmetry considerations, are given in Ref. Rikovska-Stone et al. (2007) and , is explained below. The additional contribution from the Fock terms to the scalar self-consistency significantly increases computation time and was shown to make only a small change in our results in Ref. Whittenbury et al. (2014). For this reason we neglect its correction to the mean field. The vector meson mean fields are given by

 ¯ω=∑igωim2ωρvi and ¯ρ=∑igρim2ρt3iρvi, (2.41)

where is the number density for baryon .

For , the integrand has the form

 ΞmBB′=12∑s,s′|¯uB′(p′,s′)ΓmBuB(p,s)|2Δm(→k) , (2.42)

where is the Yukawa propagator for meson with momentum and are the baryon spinors. In the above integrands we expand to isolate the momentum independent pieces and simply subtract contact terms as in Ref. Whittenbury et al. (2014). We emphasise here the importance of subtraction of the momentum independent piece, which when transformed to configuration space corresponds to a delta function.The removal of the contact terms is a common procedure because they represent very short range correlations between the baryons. To keep them in this model would not be consistent as it treats the baryons as clusters of quarks and not as point-like objects. Moreover, it is also required because of the suppression of the relative wave function at short distance originating from the repulsive hard core.

ii.5 QMC model parameters

The model dependence on all parameters was discussed in detail in Ref. Whittenbury et al. (2014). There are just the three main adjustable coupling constants, which control the coupling of the mesons to the two lightest quarks, , , and for ( for all mesons ). In addition, one has the meson masses, the value of the cut-off parameter appearing in the dipole form factors needed to evaluate the Fock terms and finally the bag radius of the free nucleon. The , , and couplings to the quarks are constrained to reproduce the standard empirical properties of symmetric (=) nuclear matter; the saturation density = 0.16 fm, the binding energy per nucleon at saturation of  MeV as well as the asymmetry energy coefficient  MeV Rikovska-Stone et al. (2007).

The , and meson masses are set to their experimental values. Whereas, the value of the mesons mass is taken to be  MeV Guichon et al. (1996, 2006); Whittenbury et al. (2014). The form factor cut-off mass, , controls the strength of the Fock terms. In Ref. Whittenbury et al. (2014) the preferred value was  0.9 GeV. Here we consider only two values, the preferred 0.9 GeV and 1.3 GeV. The latter was chosen as it produces an overly stiff EoS. All the other coupling constants in the expression for the total energy density are calculated within the QMC model or determined from symmetry considerations without further need for adjustable parameters. In particular, the tensor couplings are determined from experimental magnetic moments.

Iii Generalised beta-equilibrium matter and neutron stars

As the density of hadronic matter increases beyond saturation density nuclei dissolve to form an interacting system of nucleons and leptons. If this system survives longer than the time scale of weak interactions, 10 s, it is able to reach equilibrium with respect to beta decay and its inverse . When the total baryonic density reaches about 2 – 3 and because baryons obey the Pauli exclusion principle, it becomes energetically more favourable to create a slow and more massive hyperon, rather than another energetic nucleon. A generalised beta equilibrium develops with respect to all reactions involving either weak or strong interactions, that leads to the lowest energy state. Only two quantities are conserved in GBEM—the total charge (zero in stars) and total baryon number. Strangeness is conserved only on the time scale of strong interaction, 10 s, and lepton number is conserved only on the time-scale of tens of seconds, because of the diffusion of neutrinos out of the star Glendenning (1997).

We supplement the QMC model developed so far with non-interacting leptons and proceed to study matter in generalised beta-equilibrium (allowing for hyperons as well as nucleons). Only electrons and muons are considered, as tau leptons are too massive to be found in neutron stars. As we are considering old neutron stars, neutrinos are assumed to have radiated out of the star, so they can also be neglected. For the lepton masses we use their experimental values Beringer et al. (2012). The corresponding lepton energy density and number density are given by the usual formulas for a degenerate Fermi gas Glendenning (1997). The total energy density of the GBEM is then given by the sum of the hadron and lepton energy densities, . Similarly the total pressure is the sum of the hadron and lepton pressures and can be calculated as

 Ptotal=ρ2∂∂ρ(ϵtotalρ)=∑i,ℓμiρi−ϵtotal. (3.1)

To describe GBEM, we need to determine the lowest energy state under the two constraints of baryon number conservation and charge neutrality. For this we use the standard method of Lagrange multipliers. The equilibrium configuration of the system is then determined variationally.

We solve the following system of equations:

 0 = μi+Biλ+νQi , (3.2) 0 = μℓ−ν , (3.3) 0 = ∑iBiρvi−ρ , (3.4) 0 = ∑iBiρviQi+∑ℓρℓQℓ , (3.5)

where the baryon number, , is 0 or 1 and the charge number, , is 0 or 1. This system of equations is solved to obtain the number densities for each particle ( and ), , as well as the Lagrange multipliers (, ).

At Hartree–Fock level, the following formulas to numerically evaluate the chemical potentials, must be used to ensure we encapsulate the Fock contribution to the energy densities correctly

 μi=∂ϵtotal∂ρvi ,μℓ=∂ϵℓ∂ρℓ=√kℓ2F+m2ℓ. (3.6)

The EoS of GBEM is not valid in the outer regions (crust) of the star, where nuclei and nuclear processes become dominant. Following the customary procedure, we introduce a smooth transition between our EoS in GBEM and the standard low density EoS of Baym, Pethick and Sutherland (BPS) Baym et al. (1971) at low density. In order to calculate neutron star properties, such as the total gravitational mass, , within the stellar radius , we solve the TOV equations Oppenheimer and Volkoff (1939); Tolman (1939, 1934) for hydrostatic equilibrium of spherically symmetric (non-rotating) matter.

iii.1 Hadronic matter numerical results

Quark-meson coupling model numerical results were discussed in detail in Ref. Whittenbury et al. (2014). Here, we simply present the numerical results for two model variations,  GeV and the overly stiff  GeV. Figure 1 shows the species fraction for hadronic matter in generalised beta-equilibrium. The corresponding EoS are depicted in Fig. 2 .

Iv Quark model: Proper time regularised NJL model

In this section, we introduce the proper time regularised three flavour Nambu–Jona-Lasinio (NJL) model and use it to study quark matter.

We use an NJL Lagrangian which is inspired by one gluon exchange and proceed to calculate the corresponding effective potential in the mean field approximation. One gluon exchange can be approximated by constructing a current-current interaction using the conserved colour current . Performing Fierz transformations on this interaction Lagrangian in the exchanged -channel and adding it to the original, one obtains an interaction Lagrangian consisting of colour singlet and octet terms. As we are only interested in colour neutral matter, we ignore the colour octet terms. The Lagrangian density we investigate is then:

 LNJL = ¯ψ(i⧸∂−^m0)ψ (4.1) + GSN2F−1∑a=0[(¯ψλaψ)2+(¯ψiγ5λaψ)2] − GVN2F−1∑a=0[(¯ψγμλaψ)2+(¯ψγμγ5λaψ)2],

where diag. We used to denote the Gell-Mann matrices in colour space, whereas in Eq. (4.1) we use to label them in flavour space. Through the application of the Fierz transformations, one finds that the vector coupling is simply related to the scalar coupling, i.e., . However, in practice it can be constrained by some physical quantity such as a vector meson mass. We will treat the vector coupling as a free parameter, varying it from zero up to equality with the scalar coupling, in order to understand its effect on the quark equation of state.

The dynamically generated constituent quark masses in this NJL model in the MFA are then given by , where is the current quark mass of flavour . In vacuum, the explicit proper time regularised expression is

 Mi=mi+3GSMiπ2∫∞1Λ2UVdτ1τ2e−τM2i. (4.2)

At finite density Eq. (4.2) gains an additional correction.

iv.1 NJL model parameters

The NJL model has essentially five model parameters, the UV cut-off ; the scalar coupling ; the light and strange constituent quark masses or equivalently their current quark masses (there is a one-to-one relationship between constituent and current quark masses by Eq. (4.2)); and finally the vector coupling, , which is treated as a free parameter. As is common practice, we use pion phenomenology to fit the NJL model parameters. We constrain our model parameters in two ways. In parameter set PS1, we take as input the constituent quark masses  MeV and  MeV, the pion’s mass  MeV and its decay constant  MeV. By requiring these values for and , the UV cut-off, , is constrained to be  MeV. Then using the correspondence of the pole in the quark-anti-quark t-matrix in the pseudoscalar channel to the physical pion mass , along with the values and the scalar coupling is found to be  GeV. Finally, using , and the current quark masses can be calculated from Eq. (4.2)

In our first parameter set, PS1, the calculated current quark mass is  MeV, larger than the values typically used in the three momentum regularised versions of the model. As an additional test of sensitivity of the parameters to our fitting procedure we take instead the current quark masses as input,  MeV and  MeV, the pion’s mass  MeV and its decay constant  MeV. By following a similar procedure as above we determine the other parameters and calculate the constituent quark mass. This leads to a new and substantially different parameter set, PS2, with the constituent quark mass considerably lower. When fitting our model parameters, we are enforcing a scale in our model. With this in mind we should compare and choose the parameter set which is both consistent with hadron phenomenology (enforced through the above mentioned fitting procedures) and also favourable for modelling high density quark matter. We will compare the proper time regularised model with the three momentum regularised model Masuda et al. (2013b).

iv.2 At finite density

At finite density we have conservation of baryon number and associated chemical potentials. To handle this, an extra term is added to our NJL Lagrangian Eq. (4.1),

 LNJL→LNJL+¯ψ^μγ0ψ, (4.3)

where is the chemical potential matrix given by diag,,.

The inverse quark propagators in momentum space are now of the form

 S−1i(p)=(p0+~μi)γ0−→p⋅→γ−Mi (4.4)

for each flavour , where we have defined the reduced chemical potential

 ~μi=μi−4GVρvi≡μi−4GV⟨ψ†iψi⟩. (4.5)

Using standard methods Ripka (1997); Swanson (1992) the effective potential evaluated at finite density in the MFA is

 VNJLMF({Mi},{μi}) (4.6) = 2iNc∑i∈{u,d,s}∫d4k(2π)4 Log[k2−M2i+iϵk2−M2i0+iϵ] +∑i∈{u,d,s}(Mi−mi)28GS−∑i∈{u,d,s}(Mi0−mi)28GS − 2NC∑i∈{u,d,s}∫d3p(2π)3 Θ(~μi−Ep,i)(~μi−Ep,i) −∑i∈{u,d,s}(~μi−μi)28GV,

where is the vacuum value of the constituent quark mass and for flavour . In Eq. (4.6) we have subtracted a constant defining the effective potential, and hence the pressure (), to be zero in vacuum. With this definition of the effective potential and by the Gibbs-Duhem relation, the energy density also vanishes in vacuum. This is a common definition of the model as it is only defined up to a constant. However, some authors choose to exploit this degree of freedom by introducing a “Bag” constant, allowing the model to have non-zero values for the above thermodynamic variables in vacuum. The first two lines of Eq. (4.6) contain terms which are divergent and must be regularised. We choose to regularise using Schwinger’s covariant proper time method Klevansky (1992). The choice of the regularisation procedure is a defining decision of any NJL model.

After the analytic continuation to Euclidean space the stationary condition for the path integral translates to the condition that the mean field effective potential is determined at a global minimum and must therefore satisfy

 ∂VNJLMF∂Mi=0 and ∂2VNJLMF∂M2i≥0∀i∈{u,d,s}. (4.7)

The gap equation therefore acquires an additional contribution at finite density, which acts to restore chiral symmetry.

iv.3 Flavour independent vector interaction

We anticipate that the vector interaction is important (as is well known, see Ref. Klahn et al. (2007)) and that the strength and type of this interaction is crucial for a realistic description of quark matter. For this reason we introduce an alternative “simplified” vector interaction which is flavour independent, such that

 Lv=−gV(¯ψγμψ)2. (4.8)

This form of vector interaction has been used in many NJL studies. Particularly interesting are those that use it to produce high mass neutron stars if the coupling is large enough, see for example Refs. Masuda et al. (2013b); Hell and Weise (2014). With the vector interaction given by Eq. (4.8), rather than the flavour dependent interaction in Eq. (4.1), both the reduction in the chemical potentials and the contribution to the effective potential must be recalculated.

In symmetric two flavour quark matter ( and ) these two interactions are equivalent but differ otherwise. In asymmetric two flavour quark matter they will differ and there should be a substantial difference when strange quarks are present. In three flavour symmetric quark matter () the additional cross terms for the flavour independent interaction could give a substantial increase in pressure coming from the vector contribution. Of course, each of the quark chemical potentials will be reduced by the same amount, determined by the total quark density, as opposed to the flavour dependent interaction, where each quark’s chemical potential is only reduced by its own density. Consequently, the type of vector interaction could be important in the description of hybrid and quark stars, particularly when strange quarks are involved.

iv.4 Quarks in beta-equilibrium

Thermal equilibrium of quarks and leptons with respect to the weak and strong interactions, under the constraints of charge and baryon number conservation, is described by the following system of equations:

 23ρvu−13(ρvd+ρvs)−ρve−ρvμ = 0 (4.9) ρ−13(ρvu+ρvd+ρvs) = 0 (4.10) μd−μu−μe = 0 (4.11) μd−μs = 0 (4.12) μμ−μe = 0 (4.13)

The relations imposed in Eqs. (4.114.12) are between the thermodynamic chemical potentials and not the reduced chemical potentials. In terms of the individual quark densities (), the total quark density () and the total baryonic density () are defined as

 ρ≡ρtot3≡13∑i∈{u,d,s}ρvi. (4.14)

In the limit of zero vector coupling, the individual number density of each quark species is related to the respective chemical potential by

 ρvi=(piF)3π2=(−M2i+μ2i)3/2π2GV≠0−−−→(−M2i+~μ2i)3/2π2, (4.15)

where is the quark Fermi momentum. For non-zero vector coupling, the chemical potential in Eq. (4.15) is replaced with its reduced counterpart. The lepton chemical potentials are once again given by Eq. (3.6). This system of equations combined with the three gap equations, is then solved to determine the particle content and thermodynamic behaviour of three flavour quark matter in beta-equilibrium with leptons.

The pressure of quark matter is calculated from the thermodynamic relation

 P=−Vtotal=−VNJLMF({Mi},{μi})−Vl({μl}), (4.16)

where is the effective potential contribution of the non-interacting leptons. This gives the same pressure contribution as in Sec. II. The energy density is obtained from the following formula

 ϵtotal=Vtotal+∑i∈{u,d,s,e,μ}μiρvi, (4.17)

where in the second term is the un-reduced thermodynamic chemical potential.