Neutron stars and phase diagram in a hard-wall AdS/QCD model

We study the phase diagram of (large-$N_c$) QCD using a simplistic holographic hard-wall model with a dynamical scalar field and a homogeneous Ansatz representing a smeared instanton/baryon density. The resulting phase diagram is qualitatively consistent with expectations, including a mesonic, baryonic, quarkyonic, and quark-gluon plasma phase. As in other holographic models, we also find a baryonic popcorn transition, which appears at large chemical potential as a crossover. We then evaluate the nuclear matter equation of state, which turns out to be rather stiff with a large peaked sound velocity above the conformal limit, construct corresponding neutron stars using the TOV equations, and finally use a full numerical gravity/hydrodynamics computation to extract the gravitational wave signal of neutron star mergers.


Introduction
Determining the phase diagram of QCD, the theory governing the constituents of nuclear matter at the fundamental level, is one of the most important issues in modern theoretical physics, with applications to the cosmology of the early universe, neutron stars, hadronic physics, cosmic rays, heavy-ion collision experiments and more. The confinement of quarks and gluons into hadrons and in particular baryons is an essential feature of nature, which allows the low-energy physics to be described in terms of atoms, molecules and electrons with suitable theories at those scales. This fact is strongly tied with asymptotic freedom of QCD, which permits to describe the high-energy processes with satisfactory precision using perturbation theory. Unfortunately, ordinary matter is governed by the low-energy limit of QCD, where the theory runs into strong coupling and is tremendously hard to study theoretically. Brute-force numerical computations of QCD in the version where spacetime is discretized on a lattice is so far the most reliable source of our knowledge about the phase diagram of QCD. Unfortunately, lattice QCD works well only for small or vanishing chemical potentials, due to a technical problem known as the sign problem, which is a severe obstacle to the convergence of the numerical computations [1][2][3].
In 1997 an alternative toolbox became available to theoretical physicists by the discovery of the AdS/CFT correspondence by Maldacena [4]. In the original and most precise version, it relates a particular conformal field theory (CFT), namely 4-dimensional N = 4 super-Yang-Mills theory in the limit of large color number N c , to supergravity in 5-dimensional anti-de Sitter (AdS) space. The exciting property about the correspondence or duality is that the CFT at strong coupling is mapped to gravity at weak coupling, which means that perturbation theory in the bulk in a theory with one extra dimension can be used to studying the field theory living on the boundary of AdS at strong coupling. Conformal field theories do not possess particles as we know them, but describe operators, their fall-offs, and their anomalous dimensions. A large volume of work in the past 25 years, however, has shown that various deformations of the correspondence are possible, which has led to enormous progress in the understanding of strongly coupled gauge theories, see e.g. [5,6] for reviews.
There are basically two ways to use the holographic correspondence to studying QCD, which are top-down or bottom-up, both of which are often referred to by now as holographic QCD. The top-down category of holographic models are based on constructions in string theory that give rise to QCD-like theories in the low-energy limit. The most successful of those is arguably the Witten-Sakai-Sugimoto model [7,8], where flavor and color branes geometrically describe the theory of strong interactions with chiral quarks resulting from strings with one end point on a flavor brane and one end point on the color branes, whereas gluons have both end points on the color branes and are thus adjoint fields. In the confined phase, the flavor branes and anti-branes corresponding to left and right-handed quarks merge in the cigar-shaped bulk geometry and thus geometrically explain nonabelian chiral symmetry breaking.
A simpler holographic setup is inspired by the same type of geometric construction and a particularly simple model considered originally by Polchinski and Strassler where the AdS space is cut off at a finite value of the radial coordinate [9][10][11][12]. Meson and baryons were then subsequently added to the model, providing the basis of perhaps the simplest possible phenomenologically viable holographic description of QCD and hadron physics [13][14][15][16][17], which we shall take as basis for our study. In contrast to the Witten-Sakai-Sugimoto model, where deconfinement and chiral symmetry restoration can be set up to occur at different temperatures [18,19], the simplest version of the hard-wall model has them locked together. In order to enrich the phase structure of the theory also in the hard-wall models, we implement a so-called double hard-wall condition, where the wall for the gluons is kept at the hard-wall, but allowing for the flavor gauge and scalar fields to end on a second hard-wall, thus giving the possibility of separately adjusting the deconfinement transition.
As is well-known, the low-energy limit of the strong interactions is sufficiently well described by only two flavors of quarks, namely the up and the down quarks -enough for describing proton, neutrons and pions. The flavor symmetry of this low-energy sector is thus SU (2), which in 4-dimensional spacetime possesses a topological degree, which counts the number of instantons. In the holographic context, the instantons in the bulk are identified as the baryons of the theory.
In this paper, we are interested in the part of the QCD phase diagram where the chemical potential is finite to large. For this reason, we employ a homogeneous Ansatz, describing the instantons in the approximation suitable for large densities or equivalently finite/large chemical potential. By large densities, we mean larger than the saturation density of nuclear matter but below the asymptotic, perturbative regime. The phase diagram we find is rather rich, phenomenologically, and is consistent with common lore for the different phases and their approximate placement in the diagram [20,21]. Indeed we find mesonic, baryonic and quarkyonic (baryonic popcorn [22]) phases as we increase the chemical potential for small but finite temperatures. For large temperatures, the model loses the hard-wall as it moves behind the black hole horizon in the bulk geometry and the theory is thus in the chiral symmetry restoration phase. We nevertheless find a new speculative phase transition at large temperatures for relatively large chemical potential, where a kind of quarkyonic phase appears in our model -in the sense that deconfined quarks and baryons coexist -reminiscent of the quarkyonic phase found in the Witten-Sakai-Sugimoto model in [23]. We use the model in the baryonic phase to calculate the equation of state, which for a range of densities of order of a few times the saturation density for nuclear matter is a difficult region for other types of models to make predictions for. Indeed, traditional nuclear models are trustable only at lower densities and perturbative QCD is trustable only at much larger densities than can be realized in neutron stars. Since already simple hard-wall holographic QCD models have proven surprisingly successful in describing lowenergy hadron physics [6,24], it is interesting to explore what they are predicting for the equation of state in exactly such range of densities. The equation of state we find in our simple model is rather stiff compared to other holographic models in the literature 1 [27][28][29]; at any rate, the behavior of the speed of sound having a large peak before asymptotically tending towards the conformal value c 2 CFT = 1/3 is a rather welcome feature. Whereas perturbative quark matter that arises at sufficiently high baryonic densities has a speed of sound below the conformal value, a stiffer equation of state is required for the density ranges of neutron stars' cores to reach the observed masses of at least 2M , so that a peak structure at densities higher than saturation is what is expected [30].
While our treatment does not take care of nuclear matter at lower energies and would need substantial refinements for including the outer regions of a neutron star, we also use the obtained equation of state to simulate a neutron star merger of two neutron stars of 1.4 solar masses each at a separation distance of 45 kilometers and calculate the gravitational wave spectrum using a full-fledged numerical gravity and hydrodynamics code.
The paper is organized as follows. In sec. 2, we set up the action of the model and introduce our notation and conventions, whereas in sec. 3, we adapt the model to the case of finite temperature, introduce a generalization of the double hard-wall as well as a dynamical scalar. In sec. 4, we discuss baryonic matter in the confined and deconfined phases and finally, in sec. 5, we calculate the entire phase diagram of the model. In sec. 6, we consider neutron stars and their gravitational waves. Finally, we conclude the paper with a discussion and an outlook in sec. 7. We have delegated details of the speed of sound at large densities, the single instanton used to calibrate the model, boundary conditions and the stress-energy tensor to the appendices A, B, C and D, respectively.

Action and notation
Largely following the setup and notation of [31,32], the background geometry is taken to be that of a slice of AdS 5 ending at bulk coordinate z 0 with curvature scale indicated by L. We will work in the mostly minus convention (when we do not Wick rotate) so that the metric assumes the form: The flavor field content is given by the presence of two U (2) gauge vectors, L M (x, z), R M (x, z), dual to the left and right-handed quark currents and a bi-fundamental complex scalar Φ(x, z) dual to the order parameter of chiral symmetry breaking. The action can be divided into three main contributions: a gauge part containing Yang-Mills-like terms, a Chern-Simons term to account for flavor anomalies, and a piece containing kinetic and 1 See also [25,26] for recent reviews on holographic modeling of neutron stars.
interaction terms for the scalar. The minimal action reads [31,32]:  Requiring that the vector-vector correlator has the same asymptotic short-distance behavior as in QCD fixes M 5 L = N c /(12π 2 ) (in [14] this parameter is denoted by 1/g 2 5 ). To break chiral symmetry, the following boundary conditions are imposed in [31,32] on the IR brane: As we will see, the homogeneous Ansatz we wish to use cannot lead to a finite baryon density if both the above boundary conditions are employed. Changing boundary conditions can in principle introduce terms via the Chern-Simons action that are not present in 4D QCD: we discuss this aspect in appendix C, while for now we just keep Eq. (2.9), and will use another IR boundary condition on the L i , R i fields to fix baryon number density.
The scalar field Φ, upon imposition of the boundary conditions acquires a nontrivial vacuum as given in [32]: Here it was chosen to have the boundary condition as which is appropriate for the mesonic sector with both M q , ξ > 0 but, as shown in [31], there also exists a sector with the minus sign: it is in this sector that the single baryon configuration exists. However, the homogeneous nuclear matter Ansatz will not need this sign change: we will also discuss this topic in Appendix C.
As can be seen, in the usual scenario in which L 2 M 2 Φ = −3 (that is, for α = 1), the vacuum of the scalar reduces to the configuration where we see the explicit dependence on the infrared boundary condition of the parameter multiplying the cubic term (holographically dual to the chiral condensate).

Infrared boundary condition from a boundary action
In the previous sections we treated ξ as a free parameter that is determined just by fitting some phenomenological data. As noted, in the chiral limit M q = 0 the parameter ξ plays the role of the chiral condensate σ. If the quark mass is nonvanishing, then it can simply be regarded as the IR boundary value of the scalar field, dual to a combination of quark mass and chiral condensate.
However, following [33] the parameter ξ can be thought of as originating from a variational principle, by minimizing the energy with the addition of an IR localized boundary term to allow for a nonvanishing ξ: Let us now analyze what happens in the simpler scenario in which M 2 Φ L 2 = −3, M q = 0, α = 1; with these parameters the vacuum of the scalar field is given by so that in the vacuum, ξ is stabilized to due to the IR potential. Also the value ξ = 0 makes the energy stationary, but it is a local maximum.
The parameters λ, m b determine this value, but usually only ξ 0 is relevant, so it is simply fitted, and the former parameters are traded away. However, we note that in nuclear matter the situation can and will change, as the coupling between the scalar and the gauge fields can introduce corrections to the stable value of ξ, potentially restoring chiral symmetry for certain combinations of the other parameters (chemical potential µ and temperature T ): we will discuss this possibility in subsequent sections, after introducing suitable Ansätze to account for the presence of baryonic matter in the bulk.
3 Hard-wall model at finite temperature 3

.1 Hawking-Page transition
As shown in [34], the deconfinement transition happens via a Hawking-Page transition from the cutoff thermal AdS geometry to the AdS black hole geometry described by the metric (not yet continued to Euclidean signature) The temperature of the dual theory is determined by the periodicity of the euclidean time coordinate 0 ≤ τ < 1/T , which is unconstrained in the thermal AdS case, but fixed in the black hole case by the regularity of the near-horizon solution For a critical value z c of the horizon position, and thus for a certain temperature T d , the black hole phase becomes energetically favored, hence the deconfinement transition. This transition depends only on the cutoff scale z 0 : Here z 0 is the cutoff for the geometry: we labeled it z IR in the previous section, here we introduce this new notation in order to not make confusion in the next section, so for the moment we have two names for two quantities, z IR (the cutoff for the flavor fields) and z 0 (the cutoff for the geometry), that are usually taken to coincide, see Fig. 1.
This evaluation also neglects the role of the vacuum configuration of the scalar Φ: its vacuum energy changes the critical horizon coordinate to favor the confined phase and introduces a dependence on z IR in the expression for z c : however the trace of the scalar is proportional to N f so it would be of the same order as the backreaction of the fields on the geometry, a correction we are neglecting in the present work.

Double hard-wall model
One issue with the usual hard wall approach is that, as can be seen by the previous section on deconfinement, the critical value z c of the black hole horizon at which the phase transition happens is always lower than the cutoff z 0 = z IR . The horizon of the black hole, as soon as it appears, hides the IR brane on which the spontaneous breaking of chiral symmetry is realized by the boundary conditions of the gauge and scalar fields: this implies that as soon as the theory undergoes a deconfining transition, it also unavoidably restores chiral symmetry. However, in general these phase transitions can be widely separated [35].
To have a richer phase diagram with separate confined, deconfined, and chiral symmetry restored phases, we postulate that the geometrical cutoff and the gauge cutoff do not coincide, so we keep the notation z IR for the cutoff of the flavor gauge fields propagation which we permit to be smaller than z 0 where the geometry ends.
This kind of construction looks artificial, but an analogous configuration is actually realized in the top-down holographic QCD model of Witten-Sakai-Sugimoto when the flavor D-branes are not chosen to be antipodal: in this configuration, the flavor gauge fields that live on the branes also do not exist for all values of the radial coordinate. Since the Hawking-Page transition is sensitive only to z 0 , to have T d = T c we need to have Figure 1: On the left: a sketch of our extended hard-wall model, with two separate cutoffs, z IR for the flavor gauge fields and z 0 for the bulk geometry formed by confined gluons. On the right: the model after the Hawking-Page transition. The value of z h determines the temperature, deconfinement happens when z h < z c , while chiral symmetry is restored when z h < z IR .
Since our critical temperature corresponds to a critical z c = 2 −1/4 z 0 , having a higher value of z 0 than in [34] (since their z 0 corresponds now to our z IR ) means ending up with a lower temperature for the Hawking-Page transition departing even more from a phenomenologically accurate value, if z IR is fixed by the mass of the rho meson to be as is usually done. An alternative way to look at this choice without invoking a "double wall" structure, is simply to take the deconfinement temperature to be a free parameter, disregarding the dynamics of the deconfining Hawking-Page transition and instead imposing the change in geometry by hand: we will take this approach in drawing the phase diagram, plotting temperatures as ratios of the critical temperature T c = (πz IR ) −1 .
This choice will not affect the main results of this paper, namely the equation of state of cold nuclear matter and the properties of neutron stars, whose typical temperatures are far below the scale of T c and T d , for which we shall consider the model at T = 0. The benefit of having a setup where the deconfinement temperature T d is different from the critical temperature T c rather lies in reproducing the analogous freedom that is found in the Witten-Sakai-Sugimoto model with non-antipodal D-branes, thereby permitting to see whether any of the features of the latter can be reproduced in our much simpler model.

Scalar thermal vacuum
We have already shown in Eq. (2.13) the vacuum configuration of the scalar for nonvanishing M q and ξ in the truncated AdS 5 geometry. Since we want to study the theory also in the deconfined phase, we need to find the vacuum configuration also for the thermal black hole geometry, for the case in which the black hole horizon lies beyond the IR-brane.
Parametrizing again Φ = v(z)1, the equation of motion for the scalar alone is modified by the presence of blackening factors f (z) as keeping the same boundary conditions provided in [32]. The presence of the blackening factors modifies the configuration (2.13) by the introduction of hypergeometric functions, so that a general solution reads with C 1 , C 2 integration constants (but potentially dependent on T through z h ) to be fixed by the boundary conditions. Here we provide the full configuration for our boundary conditions: we employ the following shorthand notation Then the vacuum of the scalar in the deconfined background is given by: with the constant C being

Homogeneous baryonic matter
From the holographic point of view baryons are solitonic configurations in the flavor gauge fields: to build a many-instantons configuration accounting for interactions and minimization of the energy with respect to the moduli of the solitons is extremely challenging, so usually an approximation scheme is necessary to extract some qualitative result from the model.
Another possibility is to take the fields to only depend on z to begin with, instead of starting with individual baryons and then performing some averaging. A configuration like this 2 would, of course, be far from reality when densities are low, but as the density increases the system tends to look homogeneous up to very short distances in R 3 .
The benefit of using this approximation is threefold: first, it reduces the number of variables in our equations of motion, secondly it decouples the nonabelian part of the scalar field Φ from the baryonic matter, hence allowing us to only turn on its U (1) term, and lastly it allows us to easily solve the equations of motion numerically without fixing the configuration of certain fields, so that the resulting field configuration will be a true classical solution of the theory.
The gauge fields and the scalar are taken to depend only on z, so as to reproduce a distribution of a large number of baryons smeared in flat three-dimensional space. This requirement alone, together with the condition (2.9) for the gauge fields, and with M q = m1 for the scalar, enforces the following form [32]: The imposition of R i = −L i is also necessary in order to obtain a nonvanishing baryon number. The field strengths are then given by: From now on we shall simplify our equations by setting (z IR = 1 corresponds to measuring all dimensionful physical quantities in some units essentially fixed by the rho mass, to be fitted later, whereas L = z IR = 1 can be achieved by appropriate rescalings.) Given the Ansatz in the previous paragraph, we can extract a one-dimensional effective Lagrangian density from the model. It is given by From this Lagrangian density we can obtain the grand-canonical potential via the holographic correspondence as Equations of motion The set of equations of motion for the z−dependent fields of the homogeneous Ansatz are easily obtained as: The boundary conditions in the UV descend from and are thus translated into Dirichlet conditions for the fields of the homogeneous Ansatz: We also need to specify conditions for the IR boundary: for the fields Φ(z) and a 0 (z) we still make use of the original boundary conditions, that is As mentioned earlier, the only IR condition we have imposed for the field H(z) is automatically satisfied by construction. We still need one more boundary condition for the field H(z) to be able to integrate the differential equation: this freedom reflects the freedom that we have to choose an arbitrary baryon number for the constructed Ansatz. The baryon number is obtained as the topological charge Since H(0) = 0, the value of B can be fixed by H(1): since we are working with an infinite volume it is more meaningful to trade the global topological charge B with its density d, obtaining the following boundary conditions for H(1) at a given value of d, and for the other fields ; a 0 (z IR = 1) = 0 ; ω 0 = ξ. Note that the field a 0 appears only through its derivative: we will thus be able to determine its shape up to an additive constant. This is not unexpected, and is in agreement with the holographic interpretation of such constant, that is, the chemical potential µ B : it can easily be understood by inspection of the Chern-Simons action, which reduces to a coupling between µ B and the baryon number (density). We are thus free to choose this value, as we are free to choose the chemical potential we are working at. Furthermore, since Eq. (4.14) contains the baryon number density (which reduces to a boundary term upon integration), it can be immediately integrated and solved for a 0 (z): . (4.24) The integration constant C can be fixed using the Neumann boundary condition at z = z IR = 1 for a 0 (z), and then the Dirichlet condition for H(z): So after all we can solve the two equations (4.13), (4.15) (using (4.24)) to find H(z), ω(z) and then plugging H(z) into Eq. (4.24) together with the boundary condition a 0 (0) = µ to find a 0 (z). Note however that with this notation, the integration constant µ doesn't give the physical chemical potential for baryons µ B , but is proportional to it via which will be important to take into account when fitting the free parameters to obtain phenomenologically meaningful quantities.

Confined phase
To move to the finite temperature theory we perform a Wick rotation as t = −iτ : because of this, the time components of vector fields will be defined as L 0 = iL τ . Keeping the same homogeneous Ansatz as that of Eq.
Consistency of the Chern-Simons term with the other ones requires the field a τ to be imaginary, so we rewrite it as a τ = −i a 0 with a 0 (z) a real field. We then define the euclidean Lagrangian L E as L E = −iL. With these prescriptions we end up with the very same form of the Lagrangian as in the zero temperature scenario, and thus with the same equations of motion: no temperature induced corrections to the baryonic matter appear in this particular phase. By substituting the numerical solution of the equations of motion for this phase into the action, we can compute the grand canonical potential.

Deconfined phase
When the geometry undergoes a transition to the black hole background, the effective onedimensional Lagrangian (and correspondingly the equations of motion) change, including additional blackening factors in some terms: leading to the following set of equations: Boundary conditions can be imposed as in the T = 0 scenario. This time z h (hence temperature) enters directly into the equations of motion, so the corresponding nuclear matter solution will have temperature induced corrections: this translates in the phase diagram T vs µ into the phase transition line between mesonic phase and nuclear matter phase no longer being a straight vertical line (a feature already seen for pointlike instantons). Numerical evaluation confirms this and also establishes the transition line remaining of first order. To show how this happens, let us start from a simplified picture, in which the value of ξ does not affect the shape of H(z): if this was the case, then the only relevant additional term in the stabilization of ξ would be given by the one that explicitly contains ω 0 (z) Upon integration over z and exploiting the assumption that ∂ ξ H = 0, we find the following stable value for ξ that minimizes the energy: where λ is the self coupling of ξ, see eq. (2.17). The integral J depends on the density and behaves roughly as d 2/3 , but in general it is sufficient to note that it is a monotonic increasing function of the baryon density. We note that in nuclear matter we then have indeed a dynamical chiral condensate. More interestingly, the additional term can in principle allow for a restoration of chiral symmetry: if the right-hand side of Eq. (5.2) becomes negative, this signals that only the solution ξ = 0 survives, and it becomes a local minimum. This would establish a critical density at which chiral symmetry is restored even at zero temperature.
In reality, the function H(z) also depends weakly on ξ and so does J , so the behavior of the energy as a function of d, ξ has to be evaluated numerically.
To perform the numerical analysis, we need to find a scheme to fit the free parameters of the model. At this stage, the free parameters are L, m b , λ. The parameter m b can be traded for the value of ξ since it is obtained by minimizing the energy. The most common fit for L has it chosen to be L −1 = 320MeV, to reproduce the ρ−meson mass, while the parameter ξ in the vacuum at zero temperature should be set as 3 ξ = 4 √ 2: however in holographic models the fits of the mesonic sector are often at odds with observables in the baryonic sector. Since in this work we are mainly concerned with the physics of baryons, we choose a fit that reproduces baryonic observables relevant to the construction of the phase diagram and neutron stars, keeping only λ as a free parameter. We choose our parameters to reproduce the mass of the baryon and the onset of the nuclear matter phase, independently of the chosen λ (neglecting binding energy): These two values are in principle not independent as they must coincide, however since we employ the homogeneous Ansatz, the system loses track of the mass of the single baryon. In order to determine the mass of a single baryon, we construct a single instanton holographically in the theory, see App. B.
If we plot the free energy (see Fig. 2 for its complete behavior for a particular choice of µ) as a function of d, ξ, we immediately note that the only two local minima lie on the two axes d = 0, ξ = 0, so that the two phases of the theory are one in which there are no baryons and the chiral symmetry is spontaneously broken (d = 0 and finite ξ), and one in which there are baryons and the chiral symmetry is restored (finite d and ξ = 0). Which one of the two local minima is the global one depends on the chemical potential µ.
To study the phases of this system we numerically find those local minima (checking that this structure of local minima holds true for each µ), compare them, and determine which one is the global minimum of the function at each given value of µ (note that the value of the minimum on the d = 0 axis is independent on µ, but the one on the ξ = 0 axis is not). For the purpose of this work we have chosen a nominal value of λ = 0.002, a value close to the small λ limit which gives the best results for the mass-radius relation of neutron stars among this class of calibrations.

Chiral symmetry restoration at finite temperature
We introduce temperature dependence by employing the deconfining geometry of the AdS black hole. Naively chiral symmetry is automatically restored as soon as the horizon of the black hole reaches z h = z IR = 1 and hides the boundary conditions for the scalar field Φ(z IR ) = ξ1. However, before this happens, the true minimum of the energy can in principle be realized for ξ = 0 even at lower temperatures, depending on the µ and T dependence of the energy density. Because of the black hole metric, the IR boundary term of the action is modified as We find that this is indeed the case, and that the chiral restoration transition once again coincides with a baryonic onset, but this time the critical chemical potential µ on (T ) is a strictly decreasing function of the temperature. Interestingly, for T = T c , the curve ends at a finite value of µ on (T c ): this way our phase diagram shows a triple point as expected for large-N c QCD.
However, depending on the choice of the parameters, the equilibrium value of ξ reaches zero for temperatures lower than T c in a second order phase transition: with our choices the critical temperatures for this dynamical chiral restoration is T = 0.939T c , for the parameter set that fits the baryonic onset and baryon mass, and T = 0.9983T c (barely visible in Fig. 6) for the one that gives more realistic neutron stars. This is clearly an artifact of the model, but may be interpreted as a trend towards a weakening of the first-order phase transition around and below the value of µ at the triple point, which in real QCD should actually be a critical endpoint. Colors from blue to red indicate an increase in baryon density, and the continuous development of a peak in the distribution away from the infrared wall can be clearly seen, signaling a crossover to a quarkyonic phase.

Quarkyonic phase
We now turn our attention to the high µ region: the equilibrium density d is a monotonically increasing function of µ, so this means also analyzing the solutions of the equations of motion at high density. As we increase the density, we observe a continuous deformation of the baryon number density distribution, as it changes from being peaked on the infrared brane, where its boundary condition encodes the value of d, to a configuration where it develops a second peak at a distance from the hard-wall which is monotonically increasing with d (see Fig. 3 for a set of density profiles where the peak development is manifest). This appearance of a peak in the distribution at a finite distance from the hard wall (in general, from the position towards which the baryons are pulled by the gravitational background) is reminiscent of the "popcorn transition" in the non-homogeneous scenario. It was argued that the distribution in the holographic direction should be related to a spectrum of energies for the condensed baryons [37], so it is dual to a Fermi sphere. However the baryon at this level is still a classical object, neither fermion nor boson, so the fermionic nature has to come from quarks: this transition then indicates the onset of a quarkyonic phase for cold and dense nuclear matter [22]. The nuclear homogeneous matter in our model exhibits the same feature, performing a (continuous) "popcorn transition", where the distribution of baryon charge density is gradually pushed away from the infrared brane. This continuous distribution is in contrast to the case of individual layers of baryons considered in [22] where the position modulus in the holographic direction is determined dynamically.
Because this is not a transition that we determine by comparing order parameters, it should be regarded as a crossover 4 , with the exact value of µ at which it is happening being set by our choice of the identification of when we consider this "popcorn transition" to have happened. Note that this would not be the case if we were working within the single-instanton approximation, since in that case we would be able to determine whether a single layer or a double layer of instantons is favored at a given density. The choice we employed in this work to draw a crossover line in the phase diagram is to declare the transition to happen when the baryon density distribution is no longer a monotonic function of z, that is, when it develops a local maximum far from the infrared wall. However, the presence of this phase has no impact on the physics of neutron stars in our model, which we study below, since the high densities at which this crossover happens are never reached in their cores.
In [22] the baryonic popcorn transition was proposed as a true phase transition from nuclear matter to quarkyonic matter, whereas in our case these phases are strictly speaking one. In fact, a continuity between the two has been suggested by the lattice study of [39]. A characteristic property of quarkyonic matter has been argued in [30] to be a steep rise of the speed of sound after baryon on-set to values of order one, above the conformal value c 2 s = 1/3. As we shall show next, this indeed happens in our model quickly after baryon on-set, already before our analog of the popcorn transition, with implications for neutron stars like in [30]. The dashed sections represent unstable branches, never realized since at those chemical potentials the mesonic phase is energetically favored. The colors correspond (as will throughout all this article) to the two fit choices we employed, that we call fit A (blue) and fit B (red) in the next section. Where the curves coincide for the two fits, we draw a single purple line. The saturation density is approximately at d A 0.2 for the fit A, and at d B 0.39 for the fit B. The behavior is not monotonic and it lies above the "sound barrier" value of 1/3, reaching a maximum close to double this value, enabling the possibility of highly stiff equations of state.

Speed of sound
In early studies of holographic models, it was conjectured that the conformal value of the speed of sound presents an upper bound [40] for such models, but this was disproved later by holographic models at finite density [27,29,[41][42][43][44][45]. 5 In fact, the comparatively simple hard-wall AdS/QCD model that we are studying here also provides such a counterexample, which makes it possible to satisfy the empirical evidence from neutron stars that call for a speed of sound above the conformal value [47][48][49].
In the T = 0 limit, the squared speed of sound can be obtained from the simplified formula: In Fig. (4) we plot the squared speed of sound in the baryonic phase with the parameter choices described in the previous section. As can be seen the speed of sound rapidly increases above the sound barrier of c 2 s = 1/3, to reach almost double this value, then slowly decreases. This fits well the expectations of [30] for a quarkyonic nature of nuclear matter, which in our model is in fact continuously connected to the baryonic popcorn phase. If the latter is interpreted as a quarkyonic phase [22], nuclear matter in our model may be viewed as always having a quarkyonic nature. However, also in the Witten-Sakai-Sugimoto model, where the quarkyonic phase obtained in [23] is separated from the nuclear matter phase, a speed of sound peaking to values above the conformal one has been found in the ordinary nuclear matter phase [43] (albeit not as high as in our model).
Our numerical analysis indicates that for large µ, the squared speed of sound decreases continuously without falling below 1/3. In Appendix A we argue that asymptotically the conformal value is reached, in contrast to the Witten-Sakai-Sugimoto model, where the asymptotic value of the squared speed of sound is given by 2/5 [46]. At any rate, eventually perturbative QCD should take over, where the conformal value is reached from below [50].

Above the critical temperature
When z h reaches z IR , all boundary conditions that encode spontaneous chiral symmetry breaking remain suddenly hidden beyond a black hole horizon, effectively keeping chiral symmetry restored in the theory for any value of µ. Note that an explicit breaking (absent in this work) would be introduced via a quark mass which appears in the UV boundary condition for the scalar, so the explicit breaking cannot be undone (as expected). From these considerations follows that we can always assume the following boundary conditions for the scalar: and since its equation of motion is homogeneous, we can always trivially solve it for The presence of the black hole horizon forces us to choose Dirichlet boundary conditions for the fields L 0 , R 0 in the infrared: the usual recipe is to set A 0 (z IR ) = 0; however, we immediately see that this choice is not unique, since the value of the gauge field is not gauge invariant. We could set this value to any constant, modifying the interpretation of the value A 0 (z U V ) as the chemical potential. To avoid such redefinition, we will choose: Note that now the value of µ will determine the value of the derivative a 0 (z), so it will enter the dynamics of the fields.
Starting from the easier scenario of d = 0 (expected to hold for low chemical potential) we can see how a quark density arises: the equation of motion for a 0 in integrated form reads −2M 5 a(z) a 0 (z) = k, (5.9) with k being an integration constant. If we were to impose the usual boundary condition a 0 (z IR ) = 0, then we would conclude that k has to vanish for every temperature. However, the IR boundary is now hidden behind the horizon, where we have to impose a new boundary condition. If we take as the boundary condition the vanishing of the field on the horizon (while retaining the UV boundary condition fixing µ), then we end up with We can see from the holographic dictionary that k is then associated to the baryon density, a role played before by the baryonic matter density from H(z IR ): in this phase we have no baryon matter explicitly built into the fields (in fact, this density does not contribute to the topological charge), so this contribution can be interpreted as coming from free quarks -hence we identify the matter in this phase with a quark-gluon plasma. This construction opens up a possibility for the presence of a quarkyonic phase along the lines of [23], where in the top-down model of Witten-Sakai-Sugimoto, baryons are included forming a symmetrical layer far from the branes' joining point. For us, this amounts to relaxing the assumption that solitons are located on the IR brane (that is inaccessible in the chirally restored phase) and include both their contribution and that of free quarks in the equation of motion for the field a 0 .
To do so, we turn on the SU (2) fields with boundary conditions (i.e. adding "hair" to the AdS black hole) and we solve again numerically the coupled system of H, a 0 , choosing d that minimizes the grand canonical potential. The numerical evaluation leads to the conclusion that a phase with both baryons and free quarks arises at high chemical potential for every temperature above T c , with the transition curve being almost a straight line in the µ-T plane (see Figs. 5 and 6). The transition is of first order, and the baryonic matter appears always with its distribution peaked at a finite distance from the horizon, signaling a quarkyonic nature of this phase, which is marked "quarkyonic * " in Figs. 5 and 6. 6 While Fig. 5 shows the qualitative structure in the single-wall case (z 0 = z IR ), where there is no temperature dependence for T < T d , Fig. 6 shows the numerical results obtained for one arbitrary choice of T d < T c and the resulting temperature dependence of the otherwise straight lines separating mesonic and baryonic phases and the baryonic popcorn crossover. The two sets of choices for the free parameters we employed here (and that we It should be noted that the phase diagram that we have obtained here does not include backreactions of the flavor physics on the geometry, which is always AdS. These backreactions will certainly become increasingly important as the baryon density is increased such that the phase transition at T /T c = 1 acquires a nontrivial dependence on µ. The expectation is that T c (µ) decreases with increasing µ. If so, the phase diagram that we The green line represents a(n arbitrary in our setup) deconfining transition, below which every transition line is perfectly vertical (there is no temperature dependence). The blue and red colors represent respectively two choices of parameter sets: the blue one has parameters fitted to reproduce the correct chemical potential at the baryonic onset and the correct baryon mass from the hard-wall soliton (see App. B), while the red one has parameters chosen to have the equation of state lie predominantly in a phenomenologically relevant region, so as to obtain more realistic neutron stars. The solid curved lines for T /T c < 1 represent the corresponding first order transitions to nuclear matter and chiral restoration, while the dotted lines mark crossovers to a regime similar to the baryonic popcorn phase of [22] that is identified with a quarkyonic phase. The horizontal dashed red and blue lines represent a second-order transition to chiral restoration within the mesonic phase, with the colored regions identifying phases with broken chiral symmetry. Above the critical temperature, a phase of quark-gluon plasma exists, while for large enough chemical potentials, and as a function of temperature, baryons also appear in this phase, realizing a quarkyonic phase of a fundamentally different nature than the one below T c with coexistence of deconfined quarks and baryons.
have obtained in Fig. 6 in terms of T /T c vs µ could be mapped to one in terms of T vs µ where the horizontal lines bend down at large µ, and if its structure remains intact it could be interpreted as predicting a phase transition inside the baryonic phase, where the latter does not need to extend to high absolute temperatures but could in fact be restricted to the low-temperature domain.

Neutron stars
In the following we shall take the above results for the low-temperature equation of state at face value and consider the consequences for a toy neutron star formed by the corresponding baryonic matter. Unlike the holographic models reviewed in [25,26] we refrain from matching to realistic equation of states for conventional nuclear matter at moderate densities, thereby completely neglecting any effects from the crust of a neutron star. The modeling of a phenomenologically viable neutron star crust cannot be done purely within the homogeneous Ansatz and will require additional input, such as the physics of individual baryons and the presence of leptons, or simply a hybrid matching with other phenomenological equations of state. Also, the nontrivial phase diagram that we have obtained above for temperatures above the deconfinement temperature will play no role here.

Tolman-Oppenheimer-Volkov equations
The equations are solved by using as a boundary condition P (r = 0) = P 0 for a range of values of P 0 , and the radius of the neutron star obtained, is defined as the value of R for which P (R) = 0. Different values of P 0 will produce results of R, m that define a curve in the m-R plot to be compared with observational data. The TOV equation is dependent on the equation of state of the baryonic matter in the model, and since our quantities at equilibrium and T = 0 are solely dependent on µ, we effectively solve for µ by minimizing the grand potential to obtain d(µ), and feeding it to P, ε so as to obtain P (µ, d(µ)), ε(µ, d(µ)).
Since we work in the grand canonical ensemble, we use the holographic dictionary as before to identify the grand potential with (minus) the on-shell action. Then we can use the relations: In agreement with the identification of the grand potential as Ω = E −μN : note that we have an additional factor of 1/2 that consistently stems from the definition of µ we employed. With all these identifications we can plot the P -ε curve that defines our equa- For the blue lines (fit A), we fitted the baryonic mass and the chemical potential at baryon onset to the phenomenological values, leaving us with one free parameter which is then chosen to reproduce the best possible mass-radius relation, which in this case is quite far from what is expected, leading to quite heavy stars. For the red plots (fit B) all parameters are instead chosen to obtain a good compromise between highest mass, radius for stars of 1.4M (whose most recent measurement by NICER is represented by the gray interval), and tidal deformability (see Fig. 8). We see that in both cases the measured radius is not completely compatible with our predictions: however, the effect of a crust is expected to increase the radius of neutron stars, refining the precision of the "phenomenological" set of parameters.
tion of state, and the m-R curve that arises from the solutions to the TOV equations. The identification of P and ε with the dual holographic quantities can also be done by constructing the stress-energy tensor T i j and reading off T 0 0 and T i i from there. In general, the hydrodynamic pressure found this way is different from the thermodynamical one, but they coincide at equilibrium, as we have verified numerically (see Appendix D for this derivation of P, ε).
In Fig. 7 we plot the equation of state and the mass-radius relations for both the fit of the parameters according to (5.13), and an alternate scenario (5.14) where we do not fit any observable quantity, but we choose parameters such that the equation of state falls as much as possible within the window allowed by observational data, and such as to obtain more realistic mass-radius relations. With fit A (5.13), we find that the speed of sound inside stable neutron stars never reaches the maximum allowed by the equation of state, peaking instead at a value of c 2 s = 0.645 at the center of the most massive stable star, while the polytropic index γ = d ln P/d ln ε reaches a minimum value of γ = 1.5; with fit B (5.14) instead the allowed speed of sound values extend beyond the peak in the curve c 2 s (µ), so that the maximum value reached corresponds to the global maximum c 2 s = 0.659, while the polytropic index inside stable neutron stars reaches a lower minimum value of γ = 1.31. Similarly to the results obtained in Ref. [28] with the Witten-Sakai-Sugimoto model, we thus find that at the center of neutron stars γ drops below the value 1.75 used in Ref. [52] as criterion for the presence of quark matter, while within the holographic model the nuclear matter phase persists.  Figure 8: Tidal deformabilities as a function of star mass in the two parameter choices considered. Color coding of the plots is as in Fig. 7. We see that both fits lead to acceptable results, with the "phenomenological" one being extremely close to the lower end of the allowed region. Effects of a crust are expected to increase tidal deformabilities, pushing this last fit to more realistic predictions.
Another important property of neutron stars is the tidal deformability: in Fig. 8 we plot the tidal deformability Λ as a function of the star's mass: we see that by choosing the more phenomenologically accurate equation of state, we obtain the least likely tidal deformability for a star with mass of 1.4M , at the lower boundary of the allowed values. The equation of state obtained by fitting the baryon mass and onset instead produces a tidal deformability well within the experimental bounds.

Neutron star mergers and gravitational waves
We will now use the equation of state (EoS) found in the previous sections and tested by the TOV equations in a full numerical gravity simulation, giving rise to a gravitational wave signal, in principle detectable at the LIGO and VIRGO experiments.
The setup that we will consider is a binary of two neutron stars, each of mass 1.4 solar masses (M ) with an initial separation of 45 km. The initial configuration is made with LORENE [53] using first init bin and then coal with suitable initial enthalpies (large enough for creating an initial configuration of the desired mass. We perform the numerical gravity simulations using WhiskyTHC [54][55][56] using the method of lines with a third order Runge-Kutta (RK3) evolution, a hybrid EoS mode in the "thorn" EOS Thermal with the holographic EoS input as the cold EoS and the temperature dependence given by a simple Gamma law with coefficient Γ th = 1.75 [57]: this accounts for shock-heating effects during the merger [58] by adding a thermal component to the equation of state with P th = Γ th ρ(ε total − ε c ) = Γ th ρε th where ρ = m b d is the baryon mass density, m b is the baryon rest mass and d is the number density. The total pressure is then P from the cold equation of state plus the thermal component P th , i.e. P total = P + P th . The result for the choice of parameters fitted to the baryon mass and baryon onset is shown in Fig. 9, whereas the result for the "phenomenological" parameter choice is shown in Fig. 10, both for 1.4M + 1.4M neutron stars as initial condition.

Conclusions and outlook
We studied the phase diagram of a holographic "hard-wall" model of QCD, identifying baryonic and quarkyonic phase transitions for arbitrary temperatures, with a nontrivial temperature dependence arising from a generalization to separate infrared walls for gluonic and quark degrees of freedom. To do so, we employed the approximation of homogeneous baryonic matter, which holds true for sufficiently high densities, while we expect corrections to this picture at low densities to arise by treating baryons as individual topological solitons. At high densities we found an analogue of the baryonic popcorn transition observed before in [22] in the Witten-Sakai-Sugimoto model, which we interpret as a continuous crossover to a quarkyonic phase. Another quarkyonic phase emerged at high densities above the critical temperature for chiral symmetry restoration, which could be taken as an indication that in a more refined treatment with backreacted geometry and µ-dependent critical temperature there is a larger quarkyonic regime with embedded phase transition.
In the low-temperature phases, we have evaluated the resulting equation of state with a view towards modeling cold nuclear matter at high densities. Independently of the fit chosen for the free parameters in our model, the resulting speed of sound at zero temperature in this baryonic matter turned out to be rather high, quickly rising above the conformal value of c 2 s = 1/3 and reaching a maximum of approximately c 2 s = 0.66 before monotonically decreasing with further increases of the chemical potential. A steep rise of the speed of sound to values of order one has been argued in [30] to be associated with a quarkyonic nature of baryonic matter, which is consistent with the continuous connection between the baryonic and quarkyonic (baryonic popcorn) phases in our model.  neutron stars: to do so, we employed two different sets of parameter choices, one that reproduces the correct baryon mass and the correct baryon onset, while the other gives a good compromise for the observable properties of neutron stars. The resulting neutron stars turn out to be quite compact, with maximum masses higher than what is expected, with the radius for stars of mass 1.4M being slightly outside (smaller than) the allowed interval, and with tidal deformability for the same star exactly on the lower bound of the allowed interval. However, we have to bear in mind an important limitation in our setup: the homogeneous Ansatz is expected to only describe accurately the core of neutron stars, while we do not have a description for a phase with more diluted solitons, which would provide the stars with a softer crust. The effects of the presence of a crust are expected to be an increase in the radius of stars (at fixed mass), and an increase in the tidal deformability. To put this intuition to the test, we extended our holographic equation of state with the tabulated "SLy4" one [59] and built again neutron stars with the phenomenological set of parameters. As a result, the radius of the 1.4M star increased to over 12.3 km, while its tidal deformability increased over the allowed interval, to a value of about 775. This result opens up the possibility to obtain a better phenomenological fit by slightly changing parameters in order to obtain a lower maximum mass, which would in turn also lower the tidal deformability and the radius. We leave the modeling of a suitable crust phase for neutron stars within the same holographic model, and thus a more refined computation of neutron stars' properties, for a future work.

Appendices
A The speed of sound in the limit of large densities In this appendix we will consider taking the limit of large chemical potential, which corresponds to large densities. We will do the considerations for T = 0 and in the baryonic phase, so ξ = 0. The equation of motion (4.15) for ω 0 is trivial for ξ = 0 and hence ω 0 = 0 everywhere. H has to obey the boundary conditions (4.18) and (4.23), which for large d becomes very energetically expensive. We conjecture, based on numerical experience, that the limit of the function H(z) tends to a step function for µ → ∞: and for large µ, we have that d ∝ µ α with α ≥ 2. Integrating the Chern-Simons term by parts, we have the contributing terms in the Lagrangian: Integrating over z we obtain where we have defined the regularized integral where (d) parametrizes our ignorance about how the true solution approaches the step function in the limit of d → ∞. Minimizing the energy with respect to d yields where in the second line, we have assumed a power-law behavior of (d) ∝ d −γ . The solution to the above equation in terms of µ 3 is where I is logarithmically dependent on d. In the limit d → ∞ the speed of sound becomes: and the conformal value of the speed of sound squared is thus obtained.

B The baryon as an instanton
In this appendix, we will calculate the mass of the baryon using the Ansatz for a single instanton, following [31,32,60], giving rise to the reduced energy functionals where we have defined the complex fields the covariant derivative Dμ = ∂μ − iAμ, the abelian field strength Aμν = ∂μAν − ∂νAμ and finally, the barred indices,μ,ν run over r, z with Euclidean metric.

The equations of motion relevant for the instanton are
Dμ(a(z)Dμφ) + a(z) the Lorentz gauge condition is the boundary condition in the infrared (z = z IR ) are and in the ultraviolet (z = z U V ) they are whereas at spatial infinity, the boundary condition for the single instanton are and at the origin we have By using suitable initial conditions and a numerical scheme sometimes called arrested Newton flow, we solve the PDEs numerically and integrate the energy functionals (B.6) and (B.7) to find the energy of the instanton, which is identified with the mass of the baryon (nucleon). In Fig. 11 we show an example of the instanton solution in the bulk for ξ = 4, M 5 = Nc 12π 2 , N c = 3 and M 2 Φ = −3. For numerical convenience, we have defined a new scalar field which is defined as being relative to the vacuum solution (up to a factor a range of ξs is shown in Fig. 12. We use the data of Fig. 12 to calibrate the model by fitting the energy of the instanton to the mass of the proton/neutron (we work in the limit of unbroken isospin): where L is the curvature scale in the metric as given in Eq. (2.1).

C Boundary conditions and baryon number
Here we discuss our choice of boundary conditions for the homogeneous Ansatz. As noted in [32], a gauge transformation g L,R = exp(iα L,R ) that does not reduce to unity at the UV brane correctly reproduces the QCD global anomaly, but also carries an infrared term which is not present in QCD and which should be canceled for the holographic model to be consistently describing QCD.
The action is gauge invariant except for the Chern-Simons term, whose variation amounts to: We want to cancel the infrared term, which can be done by the following choice of boundary conditions: The second one of these conditions is automatically satisfied by the homogeneous Ansatz, while the first is in general not satisfied for µ = i: If we want to describe baryonic matter within this Ansatz, we have to give up on this choice of boundary condition, as it is incompatible with H 3 (z IR ) = 4π 2 d. However, this does not mean that the infrared variation of the Chern-Simons term cannot vanish. Let us take a look atω 1 4 ;ω The conditions (C.3) enforce the vanishing of the variation by cancellation between the L and R terms: another possibility is that they both vanish independently which is precisely what happens for the homogeneous Ansatz. The first term in Eq. (C.4) amounts to Since the form lives on a four-dimensional surface orthogonal to z, then actually α i = µ i . But from the homogeneous Ansatz we know that the only dependence of the fields is on z, and only the 0-th component of the abelian field is turned on, so that: hence the first term vanishes. Similarly, the second term is proportional to: and because of anti-symmetry of the wedge product it will select the only term with all indices different, which then vanishes because of the presence of F 0k . In the end, eachω 1 4 (α, A) independently vanishes and the unwanted anomaly is still canceled with our choice of boundary conditions in the homogeneous Ansatz.

D Pressure and energy density from stress-energy tensor
Here we want to compute the pressure and the energy density of the baryonic phase, starting from the stress-energy tensor, verifying that the definitions match the relations with the grand-canonical potential. We start by recalling the definition of the stress-energy tensor T αβ of matter fields in general relativity: where L matter indicates the Lagrangian density of matter fields, not including the √ −g factor (as it is part of the integration measure). In our model, there are two kinds of action terms for matter fields: the ones from L g , L Φ which depend on the metric, and the Chern-Simons one from L CS , which is topological and hence does not contribute to T αβ by definition (it can and will still contribute via the equations of motion). In the baryonic phase, and in the chiral limit we are considering, the scalar field vanishes because of the minimum energy boundary condition ξ = 0, hence L Φ also vanishes, so for us holds the following identification: where now the indices are to be lowered with the full curved metric g αβ , as opposed to with η αβ as was the case in Eq. (2.3). To perform the computation it is useful to use the following form of the stress-energy tensor: αβ + T We want to compute the energy density E and the hydrodynamic pressure P of the model: Let us begin by computing the energy density, starting from the first term of Eq. (D.3), remembering that with the homogeneous Ansatz, only the field strengths L ij , L iz and L 0z are turned on: The integral term can be immediately seen to give rise to −L CS after making use of the equation of motion for the field a 0 . The boundary term vanishes in the infrared, but despite the warp factor approaching zero on the ultraviolet boundary, there also a 0 vanishes. We can then make use of the ultraviolet expansion of a 0 to compute the contribution of this term: The second term in the definition of the stress-energy tensor is proportional to L matter , and since this differs from L g only by a warp factor, it is easy to check that once we raise one index and include the factor of √ −g we obtain the simple contribution: where L on−shell is the total on-shell Lagrangian integrated over z.
We now turn to compute the hydrodynamic pressure: again we divide the computation into two parts, this time starting with the trivial one T p(2) j = −g pi g ij L matter = −δ p j L matter , (D.11) which upon taking the trace and integrating becomes: The first term is instead more complicated, again requiring to rely on integration by parts and the equations of motion: differentiating with respect to g ij removes the abelian fields from the computation (with the homogeneous Ansatz only L 0z is turned on), so we find 2 ∂L matter ∂g ij = −2M 5 Tr (L ik L jk ) g kk + Tr (L iz L jz ) g zz + {L ↔ R} = −4a −2 (z)M 5 Tr H 4 4 ika jkb τ a τ b η kk + Tr H 2 4 τ i τ j η zz = −a −2 (z)M 5 2H 4 ika jka η kk + 2H 2 δ ij η zz . (D.13) As before, we now raise one index, take the trace: and then integrate, using integration by parts: As a first step, we recognize in the integral term the equation of motion for H(z), so we plug it into our expression and then we integrate again by parts, to remove the derivative of the abelian field and make the presence of the Chern-Simons Lagrangian manifest: We are left with the boundary term to evaluate: it is straightforward that the two terms evaluated on the UV boundary individually vanish, so that this is truly a purely infrared term. When evaluated on a solution of the equations of motion that is also in thermodynamic equilibrium, this boundary term is found to vanish (an evaluation we did numerically). Putting all the results together, we can then write the following identifications, where with L on−shell we mean the Lagrangian density integrated along z (that is, the Lagrangian density from the point of view of the boundary theory): The holographic duality also identifies the grand canonical potential Ω with the on-shell Lagrangian, so that the usual thermodynamic relation for homogeneous matter P V = −Ω (D. 18) correctly emerges. The relation between energy density and pressure also is correctly reproduced: Note that to find all these thermodynamical results, we had to plug a thermodynamically stable configuration in the boundary term that enters the definition of P : this is however not unexpected, since thermodynamical and hydrodynamical pressure are in principle not the same quantity, but they are required to coincide at equilibrium. The stress-energy tensor computes the hydrodynamical pressure, hence the need to make use of the condition of thermodynamical equilibrium to connect with relations from thermodynamics.