Holographic quark matter with colour superconductivity and a stiff equation of state for compact stars

We present a holographic model of QCD with a first order chiral restoration phase transition with chemical potential, mu. The first order behaviour follows from allowing a discontinuity in the dual description as the quarks are integrated out below their constituent mass. The model predicts a deconfined yet massive quark phase at intermediate densities (350 MeV<mu<500 MeV), above the nuclear density phase, which has a very stiff equation of state and a speed of sound close to one. We also include a holographic description of a colour superconducting condensate in the chirally restored vacuum and study the resulting equation of state. They provides a well behaved first order transition from the deconfined massive quark phase at very high density (mu>500 MeV). We solve the Tolman-Oppenheimer-Volkoff equations with the resulting equations of state and find stable hybrid stars with quark cores. We compute the tidal deformability for these hybrid stars and show they are consistent with LIGO/Virgo data on a neutron star collision. Our holographic model shows that quark matter could be present at the core of such compact stars.


I INTRODUCTION
There is a growing literature [1][2][3][4][5][6][7][8] attempting to use holography [9] to describe the equation of state (EoS) of deconfined quark matter to determine whether it can play a role in neutron star cores. In the first such paper [1] the exact results at finite density [10,11] for the D3/probe D7 system [12], which is dual to quark multiplets in a supersymmetric theory, were applied in this context. This model described a possible deconfined but massive dense quark phase. The resulting EoS was not stiff enough to support quark cores in compact stars such as neutron stars or hybrid stars.
Recently we adapted this system to a model that included a running anomalous dimension, γ, for the quark condensate through an effective dilaton profile [4]. That running allowed a description of the chiral restoration transition away from the deconfined massive phase. The transition occurs dependent upon whether the Brietenlohner Freedman (BF) bound [13] is violated for a scalar in the bulk dual to the chiral condensate. We parameterized the running of γ so we could control the derivative of the running at the BF bound violation point. For all choices of the derivative the model displayed a second order transition from the chirally broken phase to the chirally symmetric phase as µ increased. Generically the EoS stiffened relative to the pure D3/D7 case with the speed of sound squared in the material peaking, dependent on the chosen derivative, at 0.6c 2 for chemical potentials of order the chiral restoration point. Even this was not sufficient to convincingly support large mass neutron stars.
In this paper we present a related model in which the chirally broken phase resists transition to the chirally restored phase leading to a first order chiral transition. Prior to the transition the EoS is even stiffer than our previous example with a speed of sound rising close to the speed of light. The key extra ingredient relative to our previous work is that we have allowed a discontinuity in the holographic description at the scale of the IR constituent quark mass, where one might expect the quarks to be integrated out of the running dynamics of the gauge fields. This seems rather natural and it is interesting that our first sensible attempt has led to a very different transition and a much stiffer material.
Our model remains based on the Dirac Born Infeld (DBI) action of a D7 brane in AdS 5 with a scalar field describing the chiral condensate and a U(1) gauge field for the chemical potential. In the spirit of the model in [14] (motivated by [15,16]) we then include a mass term for the scalar by hand that allows us to include the running γ. We input this form from the perturbative QCD running result for γ (allowing it to naively extend to the non-perturbative regime). When the scalar mass passes through the BF bound a chiral condensate is induced. This philosophy has been used before in [14] to successfully describe T = µ = 0 QCD and strongly coupled theories beyond the Standard Model. Previously such models have been able to dodge the question of the deep IR physics below where the quarks become on mass shell as a result of the formation of their constituent mass. Here at finite µ an explicit description is needed for the IR regime -we propose a simple completion where the anomalous dimension is switched off in this regime and sensible physics results.
We note that this modelling does not explicitly include confinement -the philosophy is that confinement is a property of the pure Yang Mills theory at scales below the IR constituent quark mass and whilst implied in the model is not directly included. We assume that as soon as quark density switches on confinement is lost in the plasma. We find a phase with quark density and chiral symmetry breaking -we refer to this phase as a deconfined massive quark phase. It is this proposed phase that has a stiff EoS and could play a role in hybrid star cores.
We do not attempt to describe the nuclear physics phase of QCD here. Holography lies close to large N c where baryons are infinitely massive so it is perhaps not a good starting point (there are good attempts to describe the nuclear phase holographically -for example [2,[17][18][19] -and one might seek to include those descriptions in the future). Instead we simply take results from the nuclear physics/neutron star literature [20] that provide three possible varying stiffness EoS for the nuclear phase. Above 308MeV this phase takes over from the vacuum of the holographic theory. In our model though we find that at yet higher density the deconfined massive quark phase becomes the true vacuum. Then, the first order transition to the chirally symmetric phase occurs and the stiffness of the EoS falls off sharply.
The EoS can be inserted into the Tolman-Oppenheimer-Volkoff (TOV) equations to seek pressure versus radius relations inside compact stars (see for example [21,22]). Conditions for stability are discussed in [23,24] which we will review below. Here we show, by solving the TOV equations, that we do find stable stars with the deconfined massive quark phase in the core. This is our first interesting result.
The second task we undertake is to include a holographic description of a possible colour superconducting state [25] above the first order chiral restoration transition. Traditionally the holography community has declared describing colour superconductivity as very hard because the naive coloured qq order parameter is not gauge invariant and suppressed at large N c [26]. One would need to describe the breaking of the colour group and this remains a tricky issue (see [27] for recent work in this direction). However, in [28] we proposed to be more cavalier at a phenomenological level and simply allow the inclusion of gauge non-invariant operators and neglect their colour symmetry breaking effects in the dynamics. This was motivated by the idea that the coloured density of quarks and monopoles (associated with confinement) are already likely to have given Debye masses [29] to the gluons before the Cooper pairs form. In this spirit we include a new scalar field dual to the Cooper Pair in analogy to the scalar describing the chiral condensate. Since the Cooper pair carries net baryon charge it couples directly to the U(1) gauge field and the chemical potential itself generates a BF bound violation that can trigger a superconducting phase [32,33].
We construct a bottom up model of the superconducting condensation in the chiral restored phase -the transition between these two phases is second order at a finite µ (the value is dependent on the precise coupling used but typical before the true first order transition from the deconfined massive quark phase). We can describe both a two flavour condensate that breaks the colour symmetry to SU (2) or a three flavour model with a colour flavour locked condensate (although we stress we neglect the impact on the glue sector). The presence of the condensate(s) increases the pressure of the chirally restored phase and pushes the speed of sound squared up. However, for sensible choices of parameters that give a condensate at a scale in the 10s to 100s MeV, c 2 s ≤ 0.5 and this phase is not stiff enough to serve as the core of a stable hybrid star. Nevertheless the increase in pressure makes the first order transition from the deconfined massive quark phase to the now superconducting chirally restored phase occur at a lower µ. This serves to complete our model since the speed of sound rises above one in the deconfined massive phase if allowed to persist to too high µ.
Thus we present a holographic model that describes a chirally broken vacuum at µ = 0. We allow a first order transition to a nuclear phase from 308 MeV, then the holographic model predicts a first order transition to a stiff deconfined massive quark phase above about 350 MeV before a final first order transition to a superconducting chirally restored phase at around 500 MeV. The EoS still supports the hybrid stars we have previously discussed. This at least encourages experimental studies of neutron stars to seek such exotic hybrid stars.
In our final section we also seek to challenge our models of hybrid stars with the recent data from LIGO/Virgo on a neutron star collision [35] which has been used to provide constraints on the tidal deformability parameter,λ (tid) , of neutron stars near 1.4 M . We briefly review how to computeλ (tid) from the TOV equations [36][37][38] and then compute for example equations of state where we have predicted hybrid stars in the appropriate mass range. We find the predictions lie within the allowed region but they could be tested as future data accumulates. This paper is organised as follows: In section II we review our models of the QCD phases and their EoS -here we review the base D3/D7 model, the nuclear EoS we use, the holographic model of the deconfined massive phase and the chirally restored vacuum, plus finally we add colour superconducting condensates to the chirally restored phase. In section III we then solve the TOV equations to find M versus R relations for hybrid stars. In Section IV we compute the tidal deformability parameter and compare ot the LIGO/Virgo data. In Section V we summarize and conclude.

II DESCRIPTIONS OF QCD PHASES WITH µ
In this section we will work through the descriptions we use of the µ = 0 chirally broken vacuum, the nuclear physics phase, a deconfined massive quark phase and a high temperature chirally restored phase with colour superconductivity. All of these descriptions are holographic except for the nuclear phase. Since our holographic models are inspired by the D3/probe D7 model we review that briefly first. Note that we do not include temperature in any of these discussions since neutron star cores are likely cool. In principle though one could straightforwardly include temperature holographically by allowing a black hole spacetime [9].
IIA Review of the base D3/D7 probe model At strong 't Hooft coupling λ and in the large N c limit, the dual description of the N = 4 SU (N c ) SY M theory is given by a classical type IIB SUGRA in an AdS 5 × S 5 space time [9]. The flavor sector can be introduced as N f D7 branes extended along the AdS 5 geometry and warping an S 3 sphere [12]. At zero temperature the metric background is where (t, x) are the gauge theory coordinates, the ρ and Ω 3 are on the D7 brane world volume and two transverse directions to the D7 brane are χ and Ω 1 . The energy scale of the boundary theory corresponds to the radial direction r 2 = ρ 2 + χ 2 . The AdS radius is denoted by R. To study quarks, consider a D7 probe brane in the background geometry in a quenched approximation when N f << N c . There is also a U (1) gauge field A a where a = 0, 1, .., 7 runs over the world volume of the D7 brane. The DBI action for the probe D7 branes is here T D7 is the D7 brane tension, the world volume coordinates are ξ a , the induced metric is denoted as g ab and F ab is the world volume U (1) gauge field. The two fields we will concentrate on are the gauge field A t (ρ), with the field strength F ρt = A t (ρ), and the embedding scalar field χ(ρ) corresponding to the transverse direction to the D7 brane. They satisfy the following action at zero temperature where N 7 = N f T D7 V 3 and V 3 = 2π 2 is the volume of the unit S 3 on the D7 brane. In the AdS/CFT dictionary α 2 = 1 λ where λ is the 't Hooft coupling, T D7 = (2π) −7 α −4 and N 7 = N f Ncλ 16π 4 . We divide both sides of the D7 brane action in (3) by the volume of boundary space time R 3,1 and henceforth work with action densities.
The holographic interpretation is that the two constants of integration for χ are the quark mass m and the quark condensate c. For A t we have the chemical potential µ and the density d.
It is helpful to rescale the factor of 2πα into A t . The constants of integration in the solutions µ and d are then on a footing with the constants in χ (m and c) since to move from distances in AdS to energy scales in the field theory one multiplies by 1/2πα . Formally one needs to set a scale in the theory by picking for example the IR quark mass' value -after expressing all physical quantities as ratios of this setting scale the 2πα factors then cancel in the ratios.
Thus it is useful to work with the action (we will reinstate the overall factor of N 7 shortly) There are two constants of the motion d, c because only derivatives of A t (ρ) and χ(ρ) appear in the Lagrangian. Thus and the equation following from ∂L ∂A t One can solve for A t (ρ) and χ (ρ) in terms of the constants c and d which can be integrated analytically. The solution is [10] with B an incomplete Beta function. For c = d = 0 one finds constant solutions for A t (ρ) and χ(ρ). To obtain the physical solutions for the system with density one needs d 2 − c 2 positive. Note these solutions have χ = A t = 0 at ρ = 0 so "spike" out of the origin. The density of quarks are D3/D7 strings that pull the D7 to the origin.
The action density evaluated on the solutions is (reinstating N 7 ) One needs a regulator S 0 to obtain a finite value for the density action S. Numerically, we consider a cutoff at Λ, a factor of 20 above the IR quark mass, and subtract S 0 = − N7Λ 4 4 . Then one can define the thermodynamic density free energy by the renormalized action as F = − (S − S 0 ) and as a result an analytic form for the density free energy [10] To match the asymptotic UV form known from QCD one can pick λ = 3π/η 3 so that: We will use N f = N c = 3.
At even infintessimally small temperature, this theory is deconfined. The phase therefore describes a vacuum with a density of quarks of mass m. This mass, identified with the constituent quark mass, must be put in by hand and there is no chiral symmetry breaking mechanism. The EoS, which relates the pressure P to the energy density E is found from Here the pressure is too small and the quark interiors of stars can not support neutron stars -see [1].
IIB The µ = 0 chiral symmetry breaking phase Our base motivation here and in our previous paper [4] is to include the QCD running of the gauge coupling and chiral symmetry breaking into the D3/D7 system to see if the EoS stiffens. In [4] we added the running as a ρ dependent dilaton prefactor to the action (3). As described in [16], the crucial role the dilaton plays is to provide a running anomalous dimension for the quark bilinear operator which displays as a ρ dependent mass for the field χ (after expanding the dilaton). In the previous paper we used a set of functions that ran from γ = 0 in the UV through the critical γ = 1 in the IR (with varying derivative at this point) which indeed triggered chiral symmetry breaking. These models all showed a second order transition from the chirally broken to the chirally symmetric phase. We found the EoS stiffened around the transition so that the speed of sound became as large as 0.6c yet this was still not stiff enough to support hybrid stars.
Here we will take what appears only a slightly different approach, which is to not introduce running through a dilaton factor but directly through a ρ dependent mass term for χ. This approach has been taken previously in [14]. Our original motivation for this was that we wanted to add a field for a colour superconducting order parameter in sympathy with χ but we didn't want that field to experience the same running as χ which an overall dilaton factor would introduce. We will see that this ansatz can lead to a yet stiffer EoS. Thus we take the Lagrangian at µ = 0 If we write the additional mass term ∆m 2 purely as a function of ρ then were this term to lead to a violation of the BF bound in some range of small ρ then the instability would exist however large χ were to grow. Therefore we identify the RG scale in this term with ρ 2 + χ 2 (this naturally happens in the D7 probe action where this quantity is the radial distance in the background space).
When ∆m 2 = 0, near the boundary which corresponds to the UV, the solution is given asymptotically by χ(ρ) = m + c/ρ 2 , with c = qq of dimension three and m, the mass, of dimension one (note χ and ρ have dimension one). For non-zero ∆m 2 , the solution takes the form Here γ is precisely the anomalous dimension of the quark mass. The BF bound below which an instability occurs is given by ∆m 2 = −1 when γ = 1. To directly control the running of the dimension (which is our goal) it is best to allow ∆m 2 to have ρ dependence at the level of the equation of motion. This effectively neglects a term in the equation of motion ρχ 2 ∂ ∂L ∆m 2 -this term would in any case only be large when ∆m 2 is varying fast at the BF bound violation point.
We will fix the functional form of ∆m 2 using the one loop running of the gauge coupling in QCD with N f = 3 flavours transforming in the fundamental representation. This is found by solving: with Q the renormalization group scale.
The one loop result for the anomalous dimension of the quark mass is We stress that using the perturbative result outside the perturbative regime is a sensible but non-rigorous, phenomenological parametrization of the running.
We will identify the RG scale, Q with the AdS radial parameter ρ 2 + χ 2 in our model. Working perturbatively from the AdS result m 2 = ∆(∆ − 4) we have To find numerical solutions for χ's vacuum configuration, we need an IR boundary condition. In top down models χ (0) = 0 is the condition for a regular solution [12]. In previous papers [14] using this model this condition has been replaced by the very similar on mass shell condition Here χ 0 is the IR value of the quark mass where the on-shell condition ρ = χ is realized. Thus one shoots out form the 45 • line in the χ − ρ plane to find the value of χ 0 that gives the desired quark mass at some UV value -here we will require that the mass vanishes in the UV. The resulting solution for χ is shown in red in Figure 1. Note here and henceforth we will use χ 0 , the IR quark mass, to set the scales in the theory rather than Λ QCD in the one loop running but there remains just one scale introduced via Λ QCD .
In previous studies the solution above the on-mass shell point has been sufficient -bound states masses can be determined by looking at fluctuations of this solution. Now though we wish to compute the action of this configuration. There are two complications. Firstly since we imposed ∆m 2 at the level of the equations of motion we have neglected one term dependent on the derivative of ∆m 2 in the equation of motion so it is inconsistent to then use ∆m 2 directly in the action. Here though this is a small error since ∆m 2 only has a large derivative in a very small region of ρ and we will neglect this error.
Secondly we have no solution below ρ = χ 0 yet the chirally symmetric solution χ = 0, which we will want to compare the action of our solution to, extends all the way to ρ = 0. This problem will become worse below when we allow solutions with density where one expects the solution for χ to "spike" to the origin of the χ − ρ plane -with the current boundary conditions we will lose all of this part of the solution. Our resolution of this issue here is pragmatic, based on simply obtaining sensible looking solutions in the region interior to the 45 • line. We will set ∆M 2 = 0 in the region χ 2 + ρ 2 < 2χ 2 0 . The solutions of the equations of motion are then just those of the base D3/D7 probe system in (6). Thus for example at d = 0 they are the solutions χ = m. We will require the solution to match (χ and χ ) to our exterior solution on the χ 2 + ρ 2 = 2χ 2 0 circle. Thus we extend the solution in Figure 1 into the IR with the blue solution shown. These solutions are now a sensible approximation to the forms found for χ in complete D3/D7 models with chiral symmetry breaking.
Even now there remains an ambiguity as to the constant prefactor between our UV and IR action pieces. We will keep this ambiguity as a multiplier k IR on the IR action. In the philosophy of this modelling we assume that chiral symmetry breaking occurs at a higher RG scale than confinement. Below χ 0 the quarks should integrate out from the dynamics leaving the pure glue theory to provide confinement. Since we don't include this dynamics our IR action is likely out by a constant factor. The k IR choice is one way to include this factor in the dynamics.
Thus the solution for χ in Figure 1 is our description of the µ = 0 vacuum of QCD. We will use the action of this configuration(for a given choice of k IR ) to set our zero of potential energy. In fact this state will persist until quark density switches on at µ = χ 0 (a scale that is naturally of order 330 MeV in QCD -one third of the proton mass). However, before that point we must allow for a density of nucleons to set in.

IIC Nuclear phase
At small chemical potentials the nuclear transition in QCD is well understood: the confined, chirally broken vacuum is empty until a chemical potential of µ = 308.55 MeV when there is a first order phase transition to nuclear matter. This transition is well studied and the nuclear matter EoS has been explored in [20]. There the authors combined observations of a 1.97 solar mass neutron star with effective field theory (EFT) to construct the EoS, extrapolating with a constrained piecewise polytropic form. Here holography is probably least able to help -given its origin at infinite N c , baryons are naturally very heavy and far from the QCD limit. Thus, following several other authors [1][2][3], we will simply use the results of [20] to model the nuclear phase. Note there have been attempts to study the QCD nuclear phase  holographically, for example in [2,[17][18][19], but this will not be our focus in this paper.
Three ansatz for the EoS (soft, medium and stiff) are presented in Table 5 of [20] -they give the energy density and pressure for different densities. We have encoded their data as a Mathematica fitting polynomial for the analysis below and we plot these in Figure 2.

IID The dense quark phases
We next consider the (separate) transitions associated with the onset of quark density and to a chirally symmetric quark phase in our holographic model. We allow for a quark density by including a U(1) gauge field in addition to the action (11) Here A t has UV asymptotic solution µ + d/ρ 2 where d is the density. We apply a Legendre transformation to obtain the action in terms of the density d Then the equations of motion are: Note in the first equation we have again suppressed the term ρχ 2 ∂ ∂L ∆m 2 .
We solve the equations of motion in two steps for each d.
We divide the space using the IR value of the quark The blue curve is the solution below the constituent mass scale and the red is that above.
mass at d = 0, χ 0 . We obtain solutions for χ(ρ) and A t (ρ) in two intervals: first for 0 ≤ ρ 2 + χ 2 ≤ 2χ 2 0 , which we call the region below the constituent mass and the other 2χ 2 0 ≤ ρ 2 + χ 2 ≤ Λ 2 U V , with Λ U V a large UV cut off, which we call the region above the constituent mass.
Below the constituent mass we fixed ∆m 2 to zero. For a given d we shoot from the origin of the χ − ρ plane with different gradients for χ. We solve until we reach the surface χ 2 + ρ 2 = 2χ 2 0 when we read off χ and χ plus A t . These numerical solutions can also be checked against the analytic form in (6). External to the circle we use the running ∆m 2 from the QCD perturbative running and match the initial conditions provided from the interior on the circle. We then seek amongst those solutions the one that shoots to a zero UV quark mass. Then in the UV we can read off the value of the chemical potential from the A t solution. We repeat this for each value of d.
The results are shown in Figure 3. The chirally broken phase exhibits a second order transition where density switches on. This behaviour is controlled by the low ρ phase with ∆m 2 = 0 -it is just the transition of the N = 2 model where a spike grows from the origin of the ρ − χ plane connecting to the flat embedding. The exterior region (in red) plays no role initially. As d increases the model resists returning towards the χ = 0 chirally symmetric phase with the maximum value of χ even increasing. This is the phase we call the deconfined massive quark phase.
After d = 0.554χ 3 0 there are no non-trivial solutions that have a zero UV mass so by this value of d a transformation to the chiral restored phase must have occurred (this puts some constraints on the parameter k IR as we will see).
We compute the free energy by obtaining the on-shell action for each value of d. The integration follows the same separation into the regions above and below the constituent mass.We weight the IR piece's action by the parameter k IR . One must be careful when splitting into sub-regions that any counter term is the same for each computation. We normalize so that the vacuum energy is zero for the d = 0 embedding as previously discussed.
We show some example plots of the pressure (minus the free energy) against µ for various k IR in Figure 4. For k IR = 1 and 2 the system does not make sense. The chirally broken state ceases to exist before it stops being the true vacuum. On the other hand for k IR = 0.575 the system is more sensible -the chirally restored vacuum becomes preferred and the chirally broken state becomes metastable before it ceases to exist. This provides a sensible description of a first order chiral restoration transition. We also show the case k IR = 0.1 where the transition occurs at lower µ.  and k IR = 2 (orange). The lower dark blue line corresponds to the chirally restored phase (L = 0 phase) which asymptotes to 1/3.
Our expectation is that in the sensible systems the chirally broken phase is rather stiff. It is resisting the transition to the chirally restored phase. A good test of this is to determine the speed of sound squared, c 2 s (which is simply ∂P/∂E)-see Figure 5. We show the results for the four values of k IR in Figure 4 and also for the chirally symmetric phase χ = 0. We plot for values of µ above the transition where density switches on. The c 2 s in the chirally symmetric phase is 1/3. The speed of sound in the chirally broken phase though rises much higher and even passes through the speed of light c = 1. Note that for the cases of k IR = 1, 2 the speed of sound has rather strange behaviour including a turning pointthis suggests again that these choices of k IR do not make physical sense. The cases we will continue with between k IR = 0.1 − 0.575 have monotonic rising behaviour. For the moment we will allow the speed of sound to lie greater than one but will return to address this issue when we add a colour superconducting condensate to the chirally restored phase.
We can next set χ 0 = 330 MeV and compare the free energy of these phases to the nuclear phase's free energy. We do this in Figure 6  In the second plot in Figure 6 we show the variation of c 2 s with µ for the case of the soft nuclear EoS (the EoS is piecewise constructed so there are discontinuities -we just quote these from [20]) and the k IR = 0.575 case for the chirally broken and chirally symmetric vacua. The vertical dotted lines show where the phase transitions between phases occur. For the moment c 2 s rises above one before the final transition to the chirally restored phase. In the next section we will show that, by modifying the chirally restored phase by including a colour superconducting condensate, the transition away from the deconfined massive quark phase can occur earlier removing the region with c 2 s > 1.

IIE Colour superconducting phases
There has been considerable speculation in recent years that there may be a colour superconducting phase of high density QCD [25]. In the presence of a Fermi surface and any attractive interaction the formation of a di-quark condensate is expected [30,31]. In the two flavour theory the spinless condensate is in the fundamental representation of colour SU(3) and a single coloured qq bilinear condenses. In the three flavour theory a colour flavour locking (CFL) state is expected to form with three qq bilinears non-zero.
Holographically it has been shown that the presence of a chemical potential, through a dual gauge field A t , causes a charged scalar's mass to be driven through the BF bound and cause condensation [32,33]. Thus baryon number charged operators such as the qq bilinears would be expected to condense. The holographic dual is formally a description of gauge invariant operators and so it has proven hard to describe superconducting operators which are colour charged and should break the gauge group. However, in [28] we proposed that phenomenologically one can be more relaxed about this constraint. In a quark gluon plasma near a confining region of the phase diagram one expects a plasma of quarks but also potentially a plasma of colour magnetic monopoles that will play a part in the confinement mechanism. If these are present then the electric and magnetic gluon fields will all already have a Debye mass [29] and the gauged nature of colour will be blurred. We proposed to simply neglect the back reaction of the coloured condensates on the gluons but use holography to describe the condensation mechanism and to compute the vacuum energy. In this spirit we will include the superconducting phase into our holographic model for neutron stars.
We will describe each condensing qq operator by a scalar field ψ i that we introduce into the holographic model in analogy to the chiral condensate field χ (both are dimension 3 scalars). In addition though because the qq operator carries baryon number it will couple directly to the baryon number U(1) gauge field in the bulk. Thus we propose the action We have trialled actions where ψ i enter the square root term but have not been able to make them give sensible profiles for ψ i particularly because in the deep IR the square root approaches zero. Here the action for ψ i is the kinetic term emerging from the DBI action in the expansion where all fields are small and with the derivative promoted to a covariant derivative. This is intended in the same spirit as ∆m 2 is added, being the leading term for χ when aspects of the metric or dilaton contribute to its running γ.
Q is the quark number charge of the qq bilinear (which we set to 2). The final issue though is that we must match the running coupling strength, G[ρ], of the U(1) gauge field. In principle one should match this as the coupling runs from the perturbative regime but it may not be appropriate to just use the one loop running for α(ρ). In addition to that running one also expects this coupling to run logarithmically as one approaches the Fermi surface -see [31] for example. In addition the strength of the attraction in different channels depends on group theory factors which could suppress the Cooper pair condensation coupling by as a much as an order of magnitude [28]. Together these effects could, at larger µ change the coupling further. We will therefore take G to have the form where κ is a free parameter we will vary. Note allowing this choice enables us to find a wider set of solutions than were found in [34].
We do not expect both χ and ψ to condense together. For example, Lagrangian terms we could include such as |φ| 2 |ψ| 2 would tend to fight against any BF bound violation for one field if the other field condenses. We will therefore just concentrate on Cooper pair formation in the chirally symmetric χ = 0 phase to see how its presence effects that phase.
The equations of motion are We solve the equations between a large UV cut off where ψ i ∼ J i + O i /ρ 2 with J i a source and O i the Cooper pair vev, and the IR scale √ 2χ 0 where the running has become strong enough to cause χ condensation at µ = 0. Now we need suitable IR boundary conditions. It is not clear what to pick although any none extreme choice give similar behaviour -we pick ψ i = −ψ/χ 0 which has the same proportionality as the usual holographic superconducting case where the embedding ends on a black hole. Thus we can now set ψ( √ 2χ 0 ) to find solutions that asymptote to J = 0.
For A t we use the N = 2 theory A t at χ = 0 for various d and use the values of the solutions at ρ = √ 2χ 0 to set boundary conditions for A t externally.
Finally we note that the number of ψ i fields is easily dealt with. If there are N such degenerate fields then in the lower equation of (22) there is simply a factor of N -it can be absorbed into the normalization of ψ i . Since the top equation in (22) is linear in ψ i this rescaling does not change the solution. Similarly at the level of the action (20) rescaling ψ 2 i by 1/N whilst summing over N copies leaves the action invariant. Thus the difference between the theory of the two flavour condensate (where there is one ψ i ) and the colour flavour locked phase (with three) is just a rescaling of the condensate by a factor of √ 3. We will therefore restrict to one ψ i for the analysis to come but the free energy/pressure analysis is the same for the colour flavour locked phase.
The proof of all this construction is whether we obtain sensible phenomenology for the Cooper pair formation.
In Figure 7 we plot some example embeddings for κ = 1 and varying µ. We indeed find profiles that asymptote to J = 0 at each µ and where the gap size grows with µ. We plot O against µ for varying κ in Figure 8. Now we can see that there is a second order transition (from the chirally symmetric vacuum) to the colour superconducting phase with µ. Note that this will not be a physical transition because at lower scales the deconfined massive phase is preferred to the χ = 0 state -the true transition will be a first order transition from the deconfined massive phase to the superconducting phase. For κ = 10,  at µ 1.5χ 0 (approximately 500 MeV) where this first order transition to this phase will occur, the condensate's scale is of order χ 3 0 -(330 MeV) 3 -which is possibly large relative to the supposed gap scale although it serves as a sensible upper possible case. For κ = 1 the condensate only switches on close to µ 1.5χ 0 and can be an order of magnitude smaller which is again a sensible lower estimate for the condensate's value. As µ increases in all cases the condensate grows in rough proportion to µ.
We plot the pressure (minus the free energy) of the solutions in Figure 9. We see that the pressure of the phase is raised depending on the size of κ. We will adjust κ to move the transition to the deconfined phase (now with colour superconductivity) of Figure 6 to lower µ so that the speed of sound in the deconfined massive quark phase never rises above one. We can make the transition occur just before the speed of sound passes through 1 with κ = 0.85 for the case of k IR = 0.575; κ = 0.89 for k IR = 0.35 and finally κ = 0.94 for the case of k IR = 0.1 To display this graphically we again set χ 0 = 330MeV and compare the free energy of these phases to the nuclear phases' free energy and to the deconfined massive quark phase of the previous sections. We do this in Figure 10 (top left) showing the cases k IR = 0.1, 0.35 and 0.575. As seen before the transition to the nuclear phase occurs at 308 MeV. At 330 MeV the deconfined massive quark phase's pressure begins to rise. The k IR = 0.575 curve rapidly becomes the true vacuum relative to even the least stiff nuclear phase. The case k IR = 0.1 only becomes the true vacuum relative to the stiffest nuclear equation of state. For intermediate k IR one can achieve curves between these limits -for example k IR = 0.35 grows to dominate the medium and stiffest nuclear curves but does not replace the soft nuclear curve. Now though we also include the chirally restored vacuum with colour superconductivity curves for κ = 0.94, 0.89 and 0.85. They rise sharply in pressure and become the true vacuum in the range µ =450-500 MeV.
In the remaining plots in Figure 10 we show again the variation of c 2 s with µ in a number of these scenarios.
For example, in the upper right we show c 2 s for the soft nuclear EoS, for the k IR = 0.575 case for the chirally broken, and for the κ = 0.85 case of the colour superconducting vacuum. The vertical dotted lines show where the phase transitions between phases occur. We have tuned κ so that the transition to the colour superconducting phase occurs just before the speed of sound passes through 1 (in order to have the stiffest EoS we can). The inclusion of colour superconductivity in the chirally restored phase does raise c 2 s but only a little to around 0.4c 2 . This will not be sufficient to support neutron stars if this material forms the core as we will see in the next section. The crucial role colour superconductivity is playing here is to reduce the critical µ for the transition from the deconfined massive quark phase to ensure c 2 s doesn't rise above 1.
The lower two figures in Figure 10 show a variety of scenarios for the medium and stiff nuclear EoS. In all cases there are, with increasing µ, the phases: chiral broken; nuclear; massive deconfined quark; chirally restored with colour superconductivity. In each case c 2 s rises close to 0.7-0.8 in the nuclear phase then to close to 1 in the deconfined massive quark phase.  Figure 10: Transitions from the nuclear phase to the deconfined phase and then to superconducting phase. The different colors represent the cases of massive chirally broken phase with k IR = 0.575 (blue), k IR = 0.35 (magenta) and with k IR = 0.1 (purple) as in Figure 6, whereas the superconducting cases are κ = 0.85 (brown), κ = 0.89 (gray) and κ = 0.94 (black). The dashed vertical lines represent the transition from nuclear to chirally broken quark matter and from chirally broken quark matter to the superconducting phase.

III NEUTRON STAR MASS-RADIUS RELATIONS
The mass-radius relation for neutron stars is determined by the EoS of the neutron/quark matter. One solves the Tolman-Oppenheimer-Volkov (TOV) equations (see for example [21,22]) which are the relativistic equations that model hydrostatic equilibrium inside the stars. G N is Newton's constant. Here M and P are the mass and pressure in the star as a function of radius r. To integrate the equations we need to input the EoS E(P ), as well as the central pressure P c = P (r = 0) as initial condition, and the output are the mass M (r) and pressure P (r) of the corresponding star. The radius R of the star will be the value of r at which the pressure vanishes. Then varying the initial condition P c as a parameter we can construct a curve for the mass of the star M = M (r = R) against R.
It is useful to place the TOV equations in their dimensionless form: Where r = r 0 ξ, M = m 0 y(ξ), P = p 0 p(ξ), E = 0 e(ξ), A = as is sensible in the context of the nuclear equation of state discussed above; this choice then fixes the rest of our scale parameters.
One can make a radial perturbation about a solution. In terms of the mass vs radius curve one increases the value of the central density E c whilst keeping the same mass. If ∂M (Ec) ∂Ec > 0 then the corresponding equilibrium solution for this new configuration has a higher mass and therefore there is a deficit of mass. The gravitational force then needs to be balanced by increasing the central pressure. The forces acting on the matter in the star will therefore act to return the new configuration toward its original unperturbed state. However for the case in which ∂M (Ec) ∂Ec ≤ 0, if the star is perturbed, the forces acting on the perturbed star will act to drive it further from its original point in the mass vs radius curve. Therefore the condition for stability is given by As mentioned in [24] we can also determine the stability of a star from the mass vs radius curve using the Bardeen, Thorne and Meltzer (BTM) criteria [23] which established a simple formulation to know if all its radial modes are stable: 1. At each extremum where the M (R) curve rotates counter-clockwise with increasing central pressure, one radial stable mode becomes unstable.
2. At each extremum where the M (R) curve rotates clockwise with increasing central pressure, one unstable radial stable mode becomes stable.
We now perform these calculations for the EoS we obtained in Figure 10. To summarize our model has three parameters: χ 0 which is the IR quark mass that we have set to 330MeV; k IR which must lie below 0.575 to ensure there is a sensible transition from the deconfined massive quark phase to the chirally restored phase; and κ which determines the strength of the superconducting interaction which we have used to move the chiral restoration transition so that c 2 s is never greater than one.
We present the mass radius relations for neutron stars that we obtain in Figure 11. At low mass the stars are entirely neutron matter and the M −R plot is determined by the nuclear EoS -these are the green, orange and red lines depending on the choice of nuclear EoS. In each case though we now propose a transition to hybrid stars with deconfined massive quarks in the core. These are the blue, magneta or purple lines departing from the nuclear curves -the stable parts of these curves are marked by dashed lines. Finally there is a transition in the centre of the star to chirally symmetric quark matter leading to unstable stars. These branches angle off the deconfined massive curves down to the left -the higher example of this branch in each case is the chirally restored phase without superconductivity whilst the lower example has a superconducting condensate present.  Figure 11. Mass vs radius curves for the case of k IR = 0.575 (blue), k IR = 0.35 (magenta) and k IR = 0.1 (purple). The curves leaving the green/red/orange nuclear EoS prediction is the transition to a quark phase from Figure 10. The case with k IR = 0.35 only has a transition from medium (orange) and stiff nuclear matter (red) and the case with k IR = 0.1 only has a transition from stiff nuclear matter. The stable branch where c 2 s ≤ 1 is indicated in dashed cyan/pink . The transition to a superconducting state for κ = 0.85(brown), κ = 0.89(gray) and κ = 0.94(black) just before the speed of sound goes beyond 1 is also shown. Let's look at the top plot in detail as an example with the softest nuclear EoS. Here we only had the case k IR = 0.575 where the deconfined massive quark phase became the vacuum. To invoke a transition to the superconducting phase when the c 2 s has just risen to one we set κ = 0.85 -see the top two plots of Figure 10. The top plot of Figure 11 shows the resulting stars. Where the curve is green the nuclear phase only plays a role -the star is neutrons to the core. The blue line marks where the star has begun to have a deconfined massive quark phase in its core. If we did not include the superconducting phase but instead as in Figure 6 transitioned to the chirally restored χ = 0 phase this branch extends to the highest point. After the core of the star experiences the transition to the chirally restored phase the stars become unstable -this is the sharp transition to the blue dashed line that angles down to the left. The region of these stars which satisfy the stability criteria above and have c 2 s < 1 are marked by the dashed section of line. Finally if we allow a transition to the superconducting phase rather than the χ = ψ = 0 phase then we obtain the brown line -these stars with superconducting cores are again unstable but now the transition to them occurs at c 2 s = 1 in the deconfined massive quark phase, leaving a fully sensible picture of the dynamics at all µ. This EoS does not support neutron stars as high in mass as the observational ∼2 solar mass limit so is presumably not a good description of QCD.
In the central figure of Figure 11 we show example cases using the medium stiffness EoS for the nuclear phase.
Here there are deconfined massive quark phases for lower k IR and we show the cases of 0.575 and 0.35. The plots show the same structure and elements as for the top plot as we have described. Again there are stars with quark cores but still reaching two solar masses is a struggle.
Finally in the bottom picture we show three cases for the stiffest nuclear EoS -k IR = 0.575, 0.35 and 0.1. Between the last two of these values we find solutions with deconfined massive quark cores and a upper most mass for stable stars between 2 and 2.5 solar masses. This is a considerable success. We have taken sensible phenomenological holographic models of the QCD EoS and shown that such stars can exist within sensible choices of parameters. This leads credence to the idea that quark cores can exist in neutron stars and hopefully encourages study for signals of such cores in gravitational wave signals from neutron star collisions. It is interesting to look at the structure of stars in this range. In Figure 12 we plot the pressure vs radius profiles of the stars for the case of the stiffest nuclear equation of state and k IR = 0.35, colouring the radial regions that are nuclear matter, deconfined massive quarks and the colour superconducting phase (stars with superconducting cores are unstable). We see that the quark core of some stars can be quite substantial.

IV LIGO CONSTRAINTS FOR TIDAL DEFORMABILITES
It is expected that in a colliding binary system of two neutron stars, the tidal forces between the two objects would have a measurable effect in the gravitational wave signal that could be observed using gravitational wave detectors. In [35] the measurement of this effect was reported as a limit given for the tidal deformabilities of the two stars involved in the merger. To calculate the tidal deformability for specific solutions of the TOV equations that represent a neutron star, we follow references [36][37][38].
The tidal deformation between neutron stars in a binary system connects the EoS, that describe the matter inside neutron stars, to the gravitational wave emission during the inspiral. It has been shown that a small tidal signature arises in the inspiral below 400 Hz [39]. This signature amounts to a phase correction which can be described in terms of a single EoS dependent tidal deformability parameterλ (tid) , which is the ratio of each star's induced quadrupole moment to the tidal field of its companion in the binary system. The parameterλ (tid) depends on the EoS via both the neutron star radius R and mass M , and a dimensionless quantity k (tid) 2 called the Love number: A quick summary of this formalism is: one considers a static, spherically symmetric star of mass M placed in a time-independent external quadrupolar tidal field E ij .
In response, the star will develop a quadrupole moment Q ij . In the star's local rest frame, for large values of the radial coordinate r, the metric coefficient g tt is given by [40]: where n i = x i /r. This expansion defines the traceless tensors E ij and Q ij . To linear order, the induced quadrupole will be of the form Thus Q ij and E ij are defined as the coefficients in an asymptotic expansion of the metric at large distances from the star.
The perturbation to the metric can be expanded in spherical harmonics. If one allows just the l = 2 harmonic as a perturbation with a fixed spin axis then E and Q can be written in terms of Y 20 . In the metric Y 20 is then multiplied by a function of r, H(r). Now one solves the Einstein equations. A first degree differential equation is obtained for the radial function of the spherical harmonics H(r) and, when solved together with the TOV equations, can give the value of the tidal deformabilitȳ λ (tid) from the following expression for the Love number κ  Figure 13. The dimensionless tidal deformability as a function of the mass (in units of solar masses) for a holographic quark, and nuclear equations of state. The soft nuclear phase (green) has a transition to a chirally broken quark matter phase with k IR = 0.575 (blue). The medium nuclear phase (orange) has a transition to a chirally broken quark matter phase with k IR = 0.35 (magenta). The LIGO/Virgo upper bound ofλ (tid) = 800 at 1.4M is indicated by the horizontal line.  Figure 14. The tidal deformabilitiesλ i obtained for two stars with masses corresponding to those involved in the binary Neutron Star merger observed by LIGO and Virgo [35], corresponding to masses m 1 ∈ [1.36, 1.60]M and m 2 ∈ [1.17, 1.36]M (low-spin prior). The curves stand for the corresponding quark matter phases displayed in Fig. 13; the chirally broken quark matter phase has k IR = 0.575 (blue) and the chirally broken quark matter phase has k IR = 0.35 (magenta). The black curve is a sketch of the 90% experimental bound contour given in figure 5 of [35].
correspond to the two stars in the collision. To compare our results to these values, we show in Figure 14 how our example EoSs relate to these contours. The curves are generated by independently determining the tidal deformabilities for each of the stars involved in the merger.
To describe the binary system we take a chirp mass given by equation (33) where m 1 and m 2 are the masses of the components of the two body system. We can solve for m 1 in terms of m 2 for a fixed value of M. Then we can vary m 2 as a parameter and use the relation we have for the tidal deformability as a function of the mass of the star to obtain a relation between the tidal deformability for one of the two stars of the binary system with respect to the second one. This means that the the two stars involved in the binary neutron star merger correspond to masses m 1 ∈ [1.36, 1.60]M and m 2 ∈ [1.17, 1.36]M . We observe that the two cases considered in Figure 15 fit inside the 90% probability contour given in [35]. It is reasonable to hope that as more events are recorded that the bounds will begin to probe our models.

V DISCUSSION
In this paper we have presented a holographic model of quark dynamics that provides a fully self consistent equation of state for quark matter and that allows for the existence of hybrid stars. The model has a sequence of transitions with µ: the chirally broken vacuum first order transitions to a nuclear density phase above 308 MeV (we have not described this transition holographically); there is then a first order transition to a phase with a density of chirally broken deconfined massive quarks; finally there is a first order transition to a chiral restored but colour superconducting phase. The speed of sound in each of the nuclear and deconfined massive quark phases grows with µ until the transition to the next phase (see Figure  10 for examples). Thus phases that resist second order transitions generate stiff matter.
The holographic model of the quark dynamics has two regimes above and below the IR constituent quark mass.
Above that scale the model is an AdS scalar dual to the chiral condensate with a radially dependent mass set by the running of the anomalous dimension γ in the gauge theory. When this running violates the BF bound in AdS a chiral condensate forms. At scales below the IR constituent quark mass we have needed to introduce a distinct description in the regime where the quarks should be integrated out -we have simply turned off the running mass in this very low energy regime which seems to provide a sensible IR completion of the model. The, very natural, discontinuity is though what drives the deconfinement transition to be first order (rather than the second order transitions we saw in smoother descriptions in [4]).
In the chirally restored phase we have also introduced an AdS scalar to describe a colour superconducting Cooper pair condensate. Solutions with this present somewhat raise the pressure of the chirally restored vacuum.
Our model has three parameters: χ 0 the IR quark mass that we have set to 330 MeV and sets the scale of the theory (this scale is formally introduced through the running coupling); k IR a parameter that weights the relative contributions to the action of the two regimes above and below the constituent quark mass; and κ that controls the size of the interaction that triggers the superconducting condensate. The description of the nuclear phase taken from [20] also contains a range of EoS adding essentially an extra parameter through that choice. We have used the freedom of these parameters to construct a consistent set of phase transitions. The resulting equations of state include stiff nuclear and then deconfined massive quark phases.
We have solved the TOV equations using these equations of state to compute the structure of neutron/hybrid stars. We find hybrid stars with deconfined massive quark cores for stars in the 1-2.5 solar mass range depending on the precise choice of parameters. The superconducting chirally restored phase is not stiff enough to support stars for sensible values of the condensate.
We also compared our stable solutions with observations made by the LIGO and Virgo collaboration of the tidal deformabilities obtained from the detection of gravitational waves from a binary neutron star inspiral. We observed an agreement with the data reported by LIGO and Virgo as the two cases we considered fit inside the 90% probability contour.
We conclude that sensible holographic models can provide support to the idea that quark matter can be present at the cores of the most massive neutron stars observed. Our introduction of a deconfined massive quark phase requires a separation between chiral symmetry breaking and confinement at high density but this is far from impossible. Overall then this is an exciting conclusion that hopefully motivates further astrophysical and gravitational wave analysis of neutron stars and their collisions.