Cool baryon and quark matter in holographic QCD

TL;DR

This paper develops a holographic QCD model that captures both baryonic and quark matter phases at low temperature and high density, revealing phase transitions and a stiff equation of state.

Contribution

It introduces a bottom-up holographic model incorporating baryons and quark matter with backreaction, providing new insights into dense QCD phases.

Findings

Identifies first-order phase transitions between vacuum, baryon, and quark matter.

Shows the baryonic phase has a stiff equation of state with high speed of sound.

Finds the second baryon modeling approach aligns well with phenomenological expectations.

Abstract

We establish a holographic bottom-up model which covers both the baryonic and quark matter phases in cold and dense QCD. This is obtained by including the baryons using simple approximation schemes in the V-QCD model, which also includes the backreaction of the quark matter to the dynamics of pure Yang-Mills. We examine two approaches for homogeneous baryon matter: baryons as a thin layer of noninteracting matter in the holographic bulk, and baryons with a homogeneous bulk gauge field. We find that the second approach exhibits phenomenologically reasonable features. At zero temperature, the vacuum, baryon, and quark matter phases are separated by strongly first order transitions as the chemical potential varies. The equation of state in the baryonic phase is found to be stiff, i.e., the speed of sound clearly exceeds the value $c_s^2=1/3$ of conformal plasmas at high baryon densities.

Equations

$$N_c\to \infty ,N_f\to \infty ,N_f/N_c\equiv x_f\mathrm{fixed},g^2{N}_c\mathrm{fixed}.$$

$$S_{\text{V-QCD}}=S_{\mathrm{glue}}+S_{\mathrm{DBI}}+S_{\mathrm{CS}},$$

$$S_{\mathrm{glue}}=M^3{N}_c^2\int d^{5}x\sqrt{-\det g}\left[R-\frac{4}{3}\frac{(\partial \lambda {)}^2}{{\lambda }^2}+V_g(\lambda )\right],$$

$$S_{\mathrm{DBI}}=-\frac{1}{2}{M}^3{N}_c\mathbb{T}r\int d^{5}x\left(V_f(\lambda ,T^{†}T)\sqrt{-\det {\mathbf{A}}^{(L)}}+V_f(\lambda ,T{T}^{†})\sqrt{-\det {\mathbf{A}}^{(R)}}\right),$$

$${\mathbf{A}}_{MN}^{(L)}$$

$${\mathbf{A}}_{MN}^{(R)}$$

$$D_{M}T={\partial }_{M}T+iT{A}_M^L-i{A}_M^{R}T.$$

$$F^{(L/R)}=d{A}_{L/R}-i{A}_{L/R}\wedge A_{L/R}.$$

$$A_L\to V_L{A}_L{V}_L^{†}-id{V}_L{V}_L^{†},A_R\to V_R{A}_R{V}_R^{†}-id{V}_R{V}_R^{†},$$

$$T\to V_{R}T{V}_L^{†},T^{†}\to V_L{T}^{†}{V}_R^{†},$$

$$T=\tau (r){\mathbb{I}}_{N_f},$$

$$V_f(\lambda ,T{T}^{†})=V_{f0}(\lambda )e^{-a{\tau }^2}.$$

$$d{s}^2=e^{2A(r)}(-f(r)d{t}^2+d{\mathbf{x}}^2+f(r{)}^{-1}d{r}^2).$$

$$S_{\mathrm{DBI}}^{(0)}=-M^3{N}_c{N}_f\int d^{5}x{V}_{f0}(\lambda )e^{-{\tau }^2}\sqrt{-\det g}\sqrt{1+e^{-2A}f\kappa (\lambda )({\tau }^{\prime }{)}^2-e^{-4A}w(\lambda {)}^2({\Phi }^{\prime }{)}^2}.$$

$$\int d^{4}x\mathbb{T}r\left(F_{rt}^{(L)}+F_{rt}^{(R)}\right)=0.$$

$${\mathbf{A}}_{MN}^{(L)}=g_{MN}$$

$$+\frac{\kappa (\lambda ){\tau }^2}{2}\left(A_M{A}_N+A_N{A}_M\right),$$

$${\tilde{g}}_{MN}\equiv g_{MN}+\kappa (\lambda ){\delta }_M^r{\delta }_N^r({\tau }^{\prime }{)}^2+w(\lambda )({\delta }_M^r{\delta }_N^t-{\delta }_M^t{\delta }_N^r){\Phi }^{\prime }$$

$${\left({\tilde{g}}^{-1}\right)}^{MP}{\mathbf{A}}_{PN}^{(L)}={\delta }_N^M+w(\lambda ){\left({\tilde{g}}^{-1}\right)}^{MP}{F}_{PN}^{(L)}+\frac{\kappa (\lambda ){\tau }^2}{2}{\left({\tilde{g}}^{-1}\right)}^{MP}\left(A_P{A}_N+A_N{A}_P\right).$$

$$\sqrt{-\det {\mathbf{A}}^{(L)}}\simeq \sqrt{-\det \tilde{g}}$$

$$\displaystyle+\frac{w(\lambda)^{2}}{8}\left(\left(\tilde{g}^{-1}\right)^{MN}F^{(L)}_{NM}\right)^{2}-\frac{w(\lambda)^{2}}{4}\left(\tilde{g}^{-1}\right)^{MN}F^{(L)}_{NP}\left(\tilde{g}^{-1}\right)^{PQ}F^{(L)}_{QM}\Bigg{]}$$

$${\left({\tilde{g}}^{-1}\right)}^{MN}{F}_{NM}^{(L)}=2{\Xi }^{-1}{e}^{-4A}w(\lambda ){\Phi }^{\prime }{F}_{rt}^{(L)},$$

$$\Xi =\frac{\det \tilde{g}}{\det g}=1+e^{-2A}f\kappa (\lambda )({\tau }^{\prime }{)}^2-e^{-4A}w(\lambda {)}^2({\Phi }^{\prime }{)}^2$$

$$({\tilde{g}}^{-1}{)}_s=e^{-2A}\mathrm{diag}\left(-f^{-1}{\Xi }^{-1}(1+e^{-2A}f\kappa (\lambda )({\tau }^{\prime }{)}^2),1,1,1,f{\Xi }^{-1}\right),$$

$$\displaystyle\sqrt{-\det{\bf A}^{(L)}}\simeq\sqrt{-\det\tilde{g}}\Bigg{[}1$$

$$\displaystyle-\frac{w(\lambda)^{2}}{4}\left(\tilde{g}^{-1}\right)_{s}^{MN}F^{(L)}_{NP}\left(\tilde{g}^{-1}\right)_{s}^{PQ}F^{(L)}_{QM}\Bigg{]}\,.$$

$$S_{\mathrm{DBI}}^{(1)}$$

$$\displaystyle\phantom{=}-\frac{w(\lambda)^{2}}{8}\left(\tilde{g}^{-1}\right)_{s}^{MN}\left(\tilde{g}^{-1}\right)_{s}^{PQ}\,\mathbb{T}r\left(F^{(L)}_{NP}F^{(L)}_{QM}+F^{(R)}_{NP}F^{(R)}_{QM}\right)\Bigg{]}\,,$$

$$S_{\mathrm{CS}}=\frac{i{N}_c}{4{\pi }^2}\int {\Omega }_5,$$

$$\begin{split}\Omega_{5}&=\frac{1}{6}\mathrm{Tr}\,e^{-\tau^{2}}\!\left\{-iA_{L}\wedge F^{(L)}\wedge F^{(L)}+\frac{1}{2}A_{L}\wedge A_{L}\wedge A_{L}\wedge F^{(L)}+\frac{i}{10}A_{L}\wedge A_{L}\wedge A_{L}\wedge A_{L}\wedge A_{L}\right.\\ &+iA_{R}\wedge F^{(R)}\wedge F^{(R)}-\frac{1}{2}A_{R}\wedge A_{R}\wedge A_{R}\wedge F^{(R)}-\frac{i}{10}A_{R}\wedge A_{R}\wedge A_{R}\wedge A_{R}\wedge A_{R}\\ &+\tau^{2}\Big{[}iA_{L}\wedge F^{(R)}\wedge F^{(R)}-iA_{R}\wedge F^{(L)}\wedge F^{(L)}+\frac{i}{2}(A_{L}\!-\!A_{R})\wedge(F^{(L)}\wedge F^{(R)}+F^{(R)}\wedge F^{(L)})\\ &+\frac{1}{2}A_{L}\wedge A_{L}\wedge A_{L}\wedge F^{(L)}-\frac{1}{2}A_{R}\wedge A_{R}\wedge A_{R}\wedge F^{(R)}+\frac{i}{10}A_{L}\wedge A_{L}\wedge A_{L}\wedge A_{L}\wedge A_{L}\\ &-\frac{i}{10}A_{R}\wedge A_{R}\wedge A_{R}\wedge A_{R}\wedge A_{R}\Big{]}\\ &+i\tau^{3}\,d\tau\wedge\Big{[}(A_{L}\wedge A_{R}-A_{R}\wedge A_{L})\wedge(F^{(L)}+F^{(R)})+iA_{L}\wedge A_{L}\wedge A_{L}\wedge A_{R}\\ &-\frac{i}{2}A_{L}\wedge A_{R}\wedge A_{L}\wedge A_{R}+iA_{L}\wedge A_{R}\wedge A_{R}\wedge A_{R}\Big{]}\\ &\left.+\frac{i}{20}\tau^{4}(A_{L}-A_{R})\wedge(A_{L}-A_{R})\wedge(A_{L}-A_{R})\wedge(A_{L}-A_{R})\wedge(A_{L}-A_{R})\right\}\end{split}$$

$$12{\tilde{\Omega }}_5=12{\Omega }_5$$

$$\displaystyle\ -iA_{R}\wedge A_{R}\wedge A_{R})\big{]}+i\,\mathbb{T}r\,d\big{[}\,e^{-\tau^{2}}\tau^{2}\Phi dt\wedge(-2A_{L}\wedge F^{(L)}-6A_{L}\wedge F^{(R)}$$

$$\displaystyle\ +iA_{L}\wedge A_{L}\wedge A_{L}+6A_{R}\wedge F^{(L)}+2A_{R}\wedge F^{(R)}-iA_{R}\wedge A_{R}\wedge A_{R})\big{]}\,.$$

$$\widetilde{\Omega}_{5}=\Omega_{5}\big{|}_{\Phi=0}+\frac{1}{6}\Phi dt\wedge H_{4}^{(\Phi)}\,,$$

Full text

††preprint: KUNS-2752

1 Introduction and summary

Studying dense matter in QCD has turned out to be a hard problem with many unresolved questions remaining Brambilla:2014jmp . The main theoretical tool, perturbation theory, applies only to asymptotically high densities where QCD becomes a free theory Kurkela:2009gj . Effective methods such as chiral perturbation theory are useful to describe nuclear matter at low densities Machleidt:2011zz ; Tews:2012fj . But the ranges where these methods can be trusted leave a wide gap at intermediate densities where reliable and accurate approaches are not available. Moreover, first principles lattice field theory methods only work at small values of the baryon chemical potential due to the famous sign problem, and cannot be used to study properties of cold QCD matter. Consequently, even basic observables such as the equation of state of cold matter have significant uncertainties after applying known theoretical and experimental constraints (see, e.g., Annala:2017llu ). The intermediate densities, where uncertainties are at their largest, are physically relevant: they contain the phase transition (or possibly a crossover, or several transitions) between the nuclear matter and quark matter phases, and the densities of neutron star cores are known to lie within this region.

In the absence of applicable first-principles methods, model computations can give useful information about the properties of strongly interacting QCD matter in the regime of the transition between the baryon and quark matter phases. In this article, we will study this regime by using gauge/gravity duality. One of the weaknesses in this approach is that no exact gravity dual for QCD is known, and typically the models available in the literature have similar features as QCD but fail to reproduce in detail the thermodynamics of QCD, for example. However, recently progress towards more realistic and reliable modeling of QCD has been made, which motivates us to apply these models to cold and dense QCD matter. We will use one of the most realistic holographic models available (V-QCD).

V-QCD is a class of holographic models for QCD, obtained through a fusion Jarvinen:2011qe of two frameworks: improved holographic QCD (IHQCD) Gursoy:2007cb ; Gursoy:2007er ; Gursoy:2008bu ; Gursoy:2008za ; Gursoy:2009jd for the gluon sector, and a setup based on Sen-like tachyonic Dirac-Born-Infeld (DBI) actions for the quark sector Bigazzi:2005md ; Casero:2007ae . The former framework is inspired by five dimensional noncritical string theory, and the latter is obtained by introducing a pair of space filling branes in this background. In the Veneziano limit where one takes both the number of colors and number of flavors to infinity keeping their ratio $x_{f}\equiv N_{f}/N_{c}$ fixed, the two sectors are fully backreacted as one expects for ordinary QCD (with three colors and 2-3 light flavors). The model is not derived from string theory strictly: in the end one switches to bottom-up approach because on the one hand the results do not match precisely with known QCD phenomenology, and on the other hand the stringy derivation cannot be made exact (in particular due to working in the Veneziano limit). Therefore one generalizes the action to contain certain potential functions, which are then chosen to agree with qualitative QCD features and/or fitted to lattice and experimental data, in a rough analogy to effective field theory.

The thermodynamics of V-QCD has been studied in earlier work Alho:2012mh ; Alho:2013hsa ; Alho:2015zua and shown to agree qualitatively with several known properties of QCD, such as the main features of the phase diagram as a function of temperature and chemical potential. After comparison with lattice data for the equation of state at small chemical potential, the model was shown to produce an equation of state for cold and dense quark matter which agrees with known experimental and theoretical constraints for QCD Jokela:2018ers . A remaining major task in order to establish a model including the basic features of cold QCD is the inclusion of baryon physics in V-QCD. In this article, we will take the first steps in this direction.

Baryons are introduced as solitonic configurations in holographic models for large-$N_{c}$ QCD. In top-down construction, a D-brane joining $N_{c}$ open strings gives a baryon vertex Witten:1998xy ; Gross:1998gk and provides a baryon number through Chern-Simons (CS) terms. Following this notion, baryons have been considered in effective holographic theories, especially in the Witten-Sakai-Sugimoto (WSS) model Witten:1998zw ; Sakai:2004cn ; Sakai:2005yt . First, an approximation by a small size instanton was used to introduce a baryon Hong:2007kx ; Hata:2007mb ; Hong:2007dq ; Hashimoto:2008zw ; Kim:2008pw , and then this has been generalized to include contributions beyond the instanton approximation Cherman:2009gb ; Cherman:2011ve ; Bolognesi:2013nja ; Rozali:2013fna . A baryon solution has been constructed also as in bottom-up AdS/QCD Pomarol:2008aa ; Panico:2008it . Multi-baryon solutions are also studied in the WSS model Kaplunovsky:2012gb ; deBoer:2012ij ; Kaplunovsky:2015zsa ; Preis:2016fsp ; BitaghsirFadafan:2018uzs . For dense baryonic matter in the QCD phase diagram, homogeneous approximations have been utilized in the WSS model Bergman:2007wp ; Rozali:2007rx ; Ghoroku:2012am ; Li:2015uea ; Elliot-Ripley:2016uwb as well as another approach based on probe branes Gwak:2012ht ; Evans:2012cx . With these approximations, baryonic matter phases can be realized in the low temperature high-density region in the phase diagram.

In this article, we carry out the first study of baryons in V-QCD. We restrict to approximations where the baryon configurations are homogeneous in spatial directions and which hence simplifies the analysis considerably. We work in an isospin symmetric setup and also neglect the effects due to light quark masses. We adopt two approaches:

• •

The first, given in Sec. 3, is to introduce a nondynamical thin layer of baryons localized in the holographic coordinate. This approach is essentially equivalent to treating the baryons as point-like sources, which is the picture arising in the WSS model at large coupling Bergman:2007wp ; Rozali:2007rx .

• •

The second, slightly more advanced method given in Sec. 4, employs a homogeneous ansatz with $SU(2)$ flavor symmetry for the spatial components of the non-Abelian flavor gauge field, sourcing baryon density through the CS coupling. The region where the solution is highly inhomogeneous in the bulk is modeled as a discontinuity of the homogeneous baryon field.

We summarize the main results from both these approaches in the following. The main message is that the first method works unsatisfactorily, but the second one provides reasonable results.

First, we consider the approach with a thin layer of baryons of Sec. 3. We also include the full backreaction of this layer to the five dimensional gravitational background. We make the following observations:

• •

Stabilizing the layer of baryon matter turns out to be hard as it has the tendency of decaying by falling in the IR (deep in the bulk) where the approximations made in this approach also break down. We can choose the potentials such that the layer stays near the boundary of the five dimensional space, but as it turns out, this leads to another problem: the obtained phase diagram (even in the absence of baryons) is at odds with QCD phenomenology. In particular, confinement can be obtained only at very small chemical potentials.

• •

If we choose the action so that the baryons are present, baryonic phases appear in the expected region of the phase diagram: at low temperatures and intermediate chemical potentials. Consequently, we obtain phases where the charge is sourced in part by the baryons and in part by quark matter. We, however, also obtain a chirally symmetric baryonic phase which looks exotic from a QCD intuition. The phase structure obtained in this approach is summarized in Fig. 1.

It is apparent that the problems observed within this approach are weaknesses of the approach rather than the V-QCD model. This motivates us to consider an improved method.

We discuss the setup with a homogeneous non-Abelian bulk field in Sec. 4, and demonstrate that this indeed improves the results of Sec. 3 in several ways. We consider this baryon ansatz in the probe limit, i.e., on top of a fixed gravitational background. The main results from this approach, and therefore the main results of this article, are the following:

• •

The approach is seen to capture the coupling of the baryons to another bulk field (the tachyon) which is dual to the $\bar{q}q$ operator and therefore controls chiral symmetry and its breaking. It is the coupling to the tachyon, which was missing in the simpler approach of Sec. 3, that prevents the baryons from falling in the IR.

• •

The phase diagram has the expected structure: baryons dominate at low temperatures and intermediate chemical potentials, between the confined vacuum phase (dominant at low chemical potentials) and deconfined quark gluon plasma (QGP) phase (dominant at high chemical potentials). All phase transitions between these phases are of first order. The phase diagram for this approach is shown in Fig. 2. Notice that the thermodynamics in the vacuum and in the baryon phase is independent of the temperature, but nontrivial temperature dependence is included in the QGP phase.

• •

As the density is increased, the equation of state in the baryonic phase becomes stiff, and the speed of sound rises well above the conformal value $c_{s}=1/\sqrt{3}$. This is interesting because with stiff equations of state it is easier to pass the constraints set by observations of masses and deformability of neutron stars. The basic picture is therefore the following: the nuclear (quark) matter has a stiff (soft) equation of state, and the latent heat at the baryon to QGP transition is sizable. This agrees with the earlier analysis in V-QCD Jokela:2018ers where polytropic interpolations were used to model the baryonic phase.

Readers interested in these main results can safely skip Sec. 3, as the discussion of Sec. 4 can be followed independently.

Another important result arises as a by-product of the problem found in Sec. 3. Namely, having the correct confinement properties of the phase diagram sets a previously unknown constraint to the V-QCD models. This constraint is actually completely independent of baryon physics. That is, the coupling of the gauge fields in the DBI action of the flavor sector (function $w(\lambda)$ defined below) is constrained to agree (up to small corrections) the IR behavior predicted by string theory, which complements similar results for the other coupling functions of the model found in the literature Gursoy:2007cb ; Gursoy:2007er ; Jarvinen:2011qe ; Arean:2012mq ; Arean:2013tja ; Arean:2016hcs . This is discussed in detail in Appendix A.

The paper is organized as follows. First, in Sec. 2, we introduce the V-QCD model. In particular, we work out the CS terms, which are essential for computing the baryon physics. In sections 3 and 4, we consider the thin layer and homogeneous non-Abelian approaches for the baryons, respectively. Discussion and outlook are given in Sec. 5. Additional details on the computations are given in the appendices.

2 V-QCD and setup for baryon physics

2.1 The holographic action

We start by reviewing the two basic building blocks of V-QCD. First, improved holographic QCD Gursoy:2007cb ; Gursoy:2007er ; Gursoy:2008bu ; Gursoy:2008za ; Gursoy:2009jd gives the description of the dynamics of gluons. It is a bottom-up model for pure Yang-Mills motivated by noncritical string theory. Second, the flavor sector is introduced through a tachyonic DBI action, inspired by a space filling $D4-\overline{D4}$ configuration Bigazzi:2005md ; Casero:2007ae . The two sectors are fully backreacted in the Veneziano limit:

[TABLE]

Such backreacted models (V-QCD) were constructed in Jarvinen:2011qe , and these are the models we discuss in this article. A similar setup was considered in the probe limit in Iatrakis:2010zf ; Iatrakis:2010jb .

The relevant part of the dictionary is the following:

• •

The dilaton field $\lambda=e^{\phi}$ is dual to the $\mathrm{Tr}\ F^{2}$ operator and therefore sources the ’t Hooft coupling in Yang-Mills theory. This is the only field from the IHQCD sector which we will consider in this article.

• •

The tachyon field $T^{ij}$ is dual111More precisely, the duality is defined through the boundary Lagrangian $\propto\bar{q}T(1+\gamma_{5})q/2+\bar{q}T^{\dagger}(1-\gamma_{5})q/2$ Arean:2016hcs . to the quark bilinear $\bar{q}^{i}q^{j}$ and sources the quark mass matrix. It arises from the strings stretching between the $D4$ and $\overline{D4}$ branes.

• •

The left and right handed gauge fields $\left(A_{L/R}^{\mu}\right)^{ij}$ living on the branes are dual to the left and right handed currents $\bar{q}^{i}\gamma^{\mu}(1\pm\gamma_{5})q^{j}/2$.

The action of the full model consist of several terms:222For finite temperature studies one must also include the appropriate Gibbons-Hawking term which will be specified below.

[TABLE]

where the first term is the action of IHQCD and the other two terms describe the dynamics of the flavor branes. We will first discuss the first two terms, and the relevant pieces of the CS action $S_{\mathrm{CS}}$ will be given in Sec. 2.4.

The action for the gluon dynamics is given by

[TABLE]

where $M$ is the five dimensional Planck mass. The dilaton potential will be chosen appropriately to mimic the physics of QCD.

The full flavored DBI action of the model reads

[TABLE]

with the radicands defined through

[TABLE]

and the covariant derivative given by

[TABLE]

Our convention for the field strengths is such that

[TABLE]

The fields $A_{L}$, $A_{R}$ and $T$ are $N_{f}\times N_{f}$ matrices in the flavor space, and $\mathbb{T}r$ denotes the trace over flavor indices. Notice that the full non-Abelian DBI action is not known, and typically a symmetrized trace prescription Tseytlin:1997csa is assumed. The first few corrections as a series in $F$ are know precisely Refolli:2001df ; Koerber:2001uu ; Grasso:2002wb ; Keurentjes:2004tu . In this article we will only consider non-Abelian configurations using the first nontrivial term in the expansion on top of an Abelian background, in which case ambiguities in the prescription are absent, and it is enough to use a standard trace in (4). Under the left and right $U(N_{f})$ gauge transformations the fields transform as

[TABLE]

with $V_{L}\,V_{L}^{\dagger}=\mathbb{I}_{N_{f}}=V_{R}\,V_{R}^{\dagger}$.

We will make the simplifying assumption that the couplings $w$ and $\kappa$ depend on $\lambda$ only. Moreover we consider backgrounds where the tachyon is flavor independent:

[TABLE]

and use a Sen-like tachyon potential

[TABLE]

where we take $a$ to be a constant. Its value can be absorbed into the normalization of the tachyon, so we will set $a=1$ from now on. As for the metric, our ansatz reads

[TABLE]

In order to determine the model completely, one needs to also specify the potential functions $V_{g}(\lambda)$, $V_{f0}(\lambda)$, $\kappa(\lambda)$, and $w(\lambda)$. The general idea is the following: the asymptotics of most of these functions both in the UV ($\lambda\to 0$) and in the IR ($\lambda\to\infty$) are tightly constrained by agreement with QCD. In more detail, constraints, e.g., from confinement Gursoy:2007cb ; Gursoy:2007er , consistency of the backgrounds, linear glueball/meson trajectories Gursoy:2007er ; Jarvinen:2011qe ; Arean:2013tja , and regularity of the model at finite $\theta$-angle Arean:2016hcs set the power and in some cases the subleading logarithmic term of the IR behavior for the functions (for example $V_{g}\sim\lambda^{4/3}\sqrt{\log\lambda}$ as $\lambda\to\infty$).

In the UV, i.e., in the weak coupling regime, holographic models are generally not reliable. However, to set the best boundary conditions for the IR physics, we choose the UV behavior of the functions to agree with QCD perturbation theory: as usual we require that the correct UV dimensions of the various operators are reproduced, but in addition we require agreement with asymptotic freedom Gursoy:2007cb ; Gursoy:2007er , with RG flow of the quark mass Jarvinen:2011qe , and behavior at large quark mass Jarvinen:2015ofa . Interestingly, this is obtained if all the functions go to constants in the UV with perturbative corrections in $\lambda$.

In the intermediate region, $\lambda=\mathcal{O}(1)$, the remaining degrees of freedom in the potentials need to be fitted to QCD data from experiments and from lattice computations. This has been considered for IHQCD in Gursoy:2009jd and started for full V-QCD in Jokela:2018ers .

For the baryon physics a particularly important choice is that of the function $w(\lambda)$. We will discuss the choice in more detail below. The explicit choices of potentials which we use in this article are given in Appendix B.

2.2 Thermodynamics in the absence of baryons

We first discuss the physics and the phase diagram in the absence of baryons. Then we only have a vectorial flavor singlet gauge field $A_{L}=A_{R}=\mathbb{I}_{N_{f}}\Phi(r)dt$ giving nonzero charge density and chemical potential. Inserting the expressions for the fields and the potentials the DBI action evaluates to

[TABLE]

The thermodynamics of IHQCD has been studied in Gursoy:2008bu ; Gursoy:2008za , and thermodynamics in the V-QCD setup has been discussed at zero $\mu$ in Alho:2012mh ; Alho:2015zua and at nonzero $\mu$ in Alho:2013hsa ; Jokela:2018ers at zero quark mass. Quark mass effects have been considered in Jarvinen:2015ofa . The V-QCD action has two classes of solutions which either have a horizon or not. The “thermal gas” solutions without a horizon extend from the UV boundary to the IR singularity which is of the “good” kind Gubser:2000nd . They have trivial thermodynamics: the pressure is zero, whereas the pressure is nontrivial and $\mathcal{O}(N_{c}^{2})$ for the black hole solutions with a horizon. The different scalings of the pressure with $N_{c}^{2}$ can be interpreted as an order parameter for confinement Alho:2012mh ; Alho:2015zua .

Both thermal gas (TG) and black hole (BH) solutions have two further variants which either have or do not have a scalar tachyon hair, i.e., nonzero bulk condensate of the field $\tau$. The (non)existence of tachyon hair determines through the dictionary whether chiral symmetry is broken or not. Therefore there are in total four phases. Studies with choices of the potentials $V_{g}$, $V_{f}$, $\kappa$, and $w$ that reproduce various features of QCD Alho:2012mh ; Alho:2013hsa ; Jokela:2018ers have shown that three of these solutions may be dominant for values of $x_{f}=\mathcal{O}\left(1\right)$ relevant for ordinary QCD:

• •

The tachyonic thermal gas solution, which is identified as the description of the confined and chirally broken phase in QCD. This phase appears at low values of the temperature and chemical potential, i.e., for $T\lesssim\Lambda$ and $\mu\lesssim\Lambda$ where $\Lambda$ is the characteristic energy scale of the model which we will define precisely below in Sec. 3.2.

• •

The tachyonless black hole solutions, which are identified as the chirally symmetric deconfined quark matter phase in QCD. This phase dominates for large values of temperature or the chemical potential.

• •

Tachyonic black hole solutions, which describe an intermediate deconfined but chirally broken phase. Depending on the precise choice of the potentials this phase may or may not be present in the phase diagram. As it turns out, fits to lattice data disfavor its existence Jokela:2018ers . For the potentials which will be used in this article, it is subdominant for all values of $T$ and $\mu$ and therefore does not appear in the phase diagram.

The confinement/deconfinement phase transition, which is realized as a Hawking-Page phase transition is always first order. It is possible that stringy loop corrections, which we shall not consider in this article, turn the first order transition into a higher order transition or a crossover at low values of $\mu$ Alho:2015zua . The chiral transition may also be of second order, if the intermediate phase exists so that it is separated from the confinement/deconfinement transition, but otherwise it is of first order. A rather similar phase structure has been found also in models based on a D3/D7 brane system, see Evans:2010iy ; Evans:2010hi .

As it turns out, requiring the phase diagram to have the desired phase structure leads to a new nontrivial constraint for the potentials. More precisely, requiring that the TG phase extends to nonzero $\mu$ constrains the IR asymptotics of the function $w(\lambda)$. We discuss this in detail in Appendix A. The result is that in the IR asymptotics, $w(\lambda)\sim\lambda^{-w_{p}}$ as $\lambda\to\infty$, we must have $w_{p}\geq 4/3$. Consideration of the meson spectra (in particular the splitting between vector and axial vector mesons) sets $w_{p}\leq 4/3$ Arean:2013tja , which pins down $w_{p}=4/3$ as the only remaining possibility. This choice was indeed used in Alho:2013hsa ; Jokela:2018ers . The result complements earlier findings for the leading IR asymptotics as follows: Results on confinement and on glueball trajectories fix $V_{g}\sim\lambda^{4/3}$ Gursoy:2007cb ; Gursoy:2007er . Meson trajectories and regularity of chirally broken solutions set $\kappa\sim\lambda^{-4/3}$ Arean:2013tja . Regularity of solutions with chiral symmetry and at finite $\theta$-angle, and agreement with QCD phase diagram as a function of $x_{f}=N_{f}/N_{c}$ require $4/3\leq v_{p}\leq 10/3$ in $V_{f0}\sim\lambda^{v_{p}}$. Moreover, complete regular solutions could be found numerically only for $v_{p}\lesssim 3$ Jarvinen:2011qe ; Arean:2016hcs . These results hold up to logarithmic terms in $\lambda$ (which are also fixed for $V_{g}$ and $\kappa$ Gursoy:2007cb ; Gursoy:2007er ; Arean:2013tja ). Interestingly, after including the result of Appendix A, the power laws for $V_{g}$, $\kappa$, and $w$, determined by matching with QCD agree with expectation from noncritical string theory in the Einstein frame (even though this was not required in the fit), see Arean:2013tja . In addition, the string theory result $v_{p}=7/3$ also lies in the acceptable range.

The detailed comparison of the full V-QCD model with lattice data was initiated in Jokela:2018ers . This was done by fitting the thermodynamics of the model in the deconfined, chirally symmetric phase to lattice data near $\mu=0$ in full QCD with 2+1 flavors. One of the main results after the fit was, that the intermediate phase with deconfinement and broken chiral symmetry, was absent (as we already mentioned above). Another important result was that the extrapolated equation of state (EoS) in the deconfined QGP phase from $\mu=0$ and finite $T$ to $T=0$ and finite $\mu$ was typically in agreement with theoretical and observational bounds. The backgrounds resulting from this fit will be the starting point for this article, on top of which the baryon dynamics will be added. Extrapolated EoSs for QCD on the $(\mu,T)$-plane have also been considered earlier using holography DeWolfe:2010he ; DeWolfe:2011ts ; Knaute:2017opk ; Critelli:2017oub and field theory Rebhan:2003wn ; Peshier:2005pp ; Ratti:2005jh .

2.3 DBI action for small non-Abelian gauge fields

Then we include the baryonic terms assuming that the amplitude of the soliton is so small that it can be treated as a small perturbation on top of the background. That is, with slight abuse of notation, we replace $A_{L/R}\to\mathbb{I}_{N_{f}}\Phi dt+A_{L/R}$ and treat $A_{L}$ and $A_{R}$ (but not $\Phi$) as small perturbations. We will specify below what exactly are the leading nontrivial terms in the expansion that we will consider.

This division of the gauge field into $\Phi$ and the flavor singlet part of $A_{L}+A_{R}$ is however not well-defined for generic baryon fields $A_{L/R}$. For our purposes it is enough to fix this by requiring that the soliton part satisfies

[TABLE]

We go on developing (4) as a series at small gauge fields. We note that

[TABLE]

where $A=A_{L}-A_{R}$. A similar identity holds for ${\bf A}^{(R)}$. The last two terms in (2.3) capture the contribution from the soliton and are treated as small perturbations. We define the effective metric as

[TABLE]

so that

[TABLE]

Taking the determinant and rearranging the terms, we find

[TABLE]

where we only included the terms up to quadratic order in $F^{(L/R)}$, corresponding to the expansion in $\alpha^{\prime}$, and the leading nontrivial additional term appearing due to the presence of the background tachyon field. We further notice that

[TABLE]

where

[TABLE]

since only the antisymmetric terms in $\tilde{g}^{-1}$ contribute. Let us denote by $(\tilde{g}^{-1})_{s}$ the remaining diagonal and symmetric terms:

[TABLE]

where the indexes are ordered as in the expressions for the metric above: $(t,x_{1},x_{2},x_{3},r)$. In the last two terms of (2.3) the contributions from the antisymmetric terms exactly cancel. Putting these observations together,

[TABLE]

We are now ready to write down the leading term of the DBI action in the flavored gauge fields:

[TABLE]

where we also included the terms arising from ${\bf A}^{(R)}$ and used (13). Notice that up to quadratic order in the gauge fields the DBI action is unambiguous: the result is independent of the order of the (non-Abelian) fields. For higher order terms a specific prescription (e.g., the symmetrized trace) would need to be chosen.

2.4 Chern-Simons terms

The CS terms determine how the solitons source baryonic charge. These terms depend on a CP-odd potential $V_{a}(\lambda,\tau)$ Arean:2016hcs which must satisfy certain requirements: the normalization in the UV ($\lambda=0=\tau$) must reproduce the correct axial anomaly and perturbative corrections in $\lambda$ must vanish due to the perturbative nonrenormalization of the anomaly. In principle we could work with a generic CP-odd potential $V_{a}(\lambda,\tau)$ but we choose the string motivated ansatz $V_{a}(\lambda,\tau)=e^{-b\tau^{2}}$. The inclusion of the constant $b$ reflects the findings of Arean:2016hcs : in order for the model to have regular IR solutions in the presence of a finite $\theta$-angle, the contributions from the CS terms had to vanish faster than those coming from the DBI. The easiest way to arrange this is to take $b>1$ in the CS action. In this section we will set $b=1$ for notational simplicity. It will be reintroduced later by rescaling $\tau$.

We compute here explicitly the coupling of $\Phi$ to the instanton density arising from these terms. The relevant CS term is given by Casero:2007ae

[TABLE]

where

[TABLE]

with the understanding that the contributions from the Abelian field $\Phi$ are included in the gauge fields here. The normalization of this term is consistent with the QCD flavor anomalies Casero:2007ae .

In order to extract the coupling between the solitonic components and $\Phi$ explicitly, we substituting $A_{L/R}\to\Phi dt+A_{L/R}$ in (24) and collect the coupling terms. Recall however that $\Omega_{5}$ is well defined only up to total derivatives. As it turns out, it is convenient to first modify the definition of $\Omega_{5}$ by adding the following total derivative terms

[TABLE]

Then we find that

[TABLE]

where

[TABLE]

and we have used the fact that $d\tau\wedge d\Phi=0$ since both fields are assumed to depend on $r$ only. We note that $H_{4}^{(\Phi)}$ is closed, $dH_{4}^{(\Phi)}=0$, and exact:

[TABLE]

The total charge density is defined as

[TABLE]

where we used the $\Phi$ equation of motion. Therefore the baryon charge is given by the coupling to $\Phi$ in the CS action:

[TABLE]

where $N_{b}$ is the total baryon number. We will compute this explicitly below within the approaches considered in this article.

3 Baryons as a thin layer of noninteracting bulk matter

The first approach we consider is to include baryons as a layer of noninteracting solitons. The layer is located in equilibrium at a finite nonzero value of the bulk and assumed to have a zero width in the holographic direction. This setup is similar to the approach in the WSS model where the baryons were treated as point-like sources in the limit of large coupling Bergman:2007wp ; Rozali:2007rx . When comparing to the WSS model it is useful to recall that the dynamics of chiral symmetry breaking can be discussed in terms of tachyon condensation as we did for V-QCD in Sec. 2 Bergman:2007pm ; Dhar:2007bz ; Dhar:2008um ; Jokela:2009tk . Notice however that in our model there will not be a limit (similar to the large coupling limit in the WSS model) in which the sizes of the solitons are suppressed. Since our approach in this section requires the extent of the baryons to be zero in the holographic direction, it should be considered as a rough approximation. Notice however that as we are neglecting the interactions between the solitons and our background solution is independent of the spacetime coordinates, the sizes of the solitons in spatial directions are irrelevant. The easiest approach is to consider the solitons to be of zero size in this direction also. We will consider another approach in Sec. 4 which will take the effects due to the finite size and interactions into account at least partially.

3.1 Setup

In order to establish the thermodynamics in the setup, we need to compute the mass of a single soliton (integral of the expanded DBI action). As discussed above, we will essentially treat the soliton as point-like. We first consider the simplest approach, where the tachyon field $\tau$ is completely ignored – this is a good approximation if the soliton is located very close to the boundary. We will argue how the tachyon effects can be taken into account later. In this approximation, we obtain

[TABLE]

Notice that the DBI action still involves the nontrivial effective metric $\left(\tilde{g}^{-1}\right)_{s}$. In order to simplify the analysis, we can rescale the coordinates and the gauge fields. Since the soliton is localized in $r$ and spatial coordinates we may rescale them but not the time. We define

[TABLE]

These rescalings were chosen such that the CS term remains invariant and the factors of $\left(\tilde{g}^{-1}\right)_{s}$ can be absorbed in the determinant of the rescaled metric:

[TABLE]

where the metric $\hat{g}$ is conformally flat,

[TABLE]

The result is similar in form to what has been found in probe brane models. Because the soliton is assumed to be localized in the $r$-direction, the result boils down to the Yang-Mills action in flat space where the solution (the BPST instanton) is known. The action may be evaluated as

[TABLE]

where we reinstated the unrescaled coordinate $r$. The location of the baryon $r_{b}$ will be determined by minimizing the action as we will show below. For a soliton corresponding to an antibaryon the sign of the CS term is opposite.

There is however no obvious reason (as we shall demonstrate below) why the soliton should stabilize very close to the UV boundary in our model. Therefore the tachyon dependence should not be discarded. We will discuss how they affect the computation starting with the CS term.

In the CS action the tachyon dependent terms are given in (2.4) and (2.4). For a soliton localized in the $r$-direction the CS term can be written as

[TABLE]

We expect that the integral here is quantized in units of $24\pi^{2}i$ also when the soliton is not close to the boundary, and the result therefore is the same as in (35). Since the integral couples simply to $\Phi(r_{b})$ this is consistent with the baryons carrying a fixed charge. Indeed as the form $H_{4}^{(\Phi)}$ is exact, the integral becomes a boundary term, which suggests that the quantization can be read off by inserting the asymptotic form of the soliton solution in this expression. However, in the absence of an explicit soliton solution which would take into account the coupling to the tachyon, we have not been able to prove this.

This result implies in particular that the integrals over the various terms in (2.4) on the soliton solution will need to grow large in order to compensate the factor $e^{-\tau^{2}}$ in the expression of this form if the soliton is located deep in the IR. That is, the amplitude and/or size of the soliton needs to grow large. Therefore, it is essential that for the approximations done in this section to work, the baryon is not located very deep in the IR.

Without better control of the soliton solution, it is hard to evaluate its contribution to the DBI action, i.e., the soliton mass, in the presence of the tachyon corrections. The main addition due to the tachyon is the factor $e^{-\tau^{2}}$ in the potential of the DBI term, see (2.3). The quantization argument of the CS term suggests that the contribution of the soliton grows if it is moved towards the IR such that it roughly cancels this term. Therefore our best guess for the effects of the tachyon is that they are absent at least when the soliton is not too deep in the IR, that is, we will also use the expressions in (35) in the presence of the coupling to the tachyon. We remind that we will consider a different approach in Sec. 4 which will capture the coupling to the tachyon.

In the WSS model the baryon action is obtained through a D4 action, with the brane wrapping the $S^{4}$ of the geometry, or equivalently by considering an expansion of the D8 actions at small gauge fields Sakai:2004cn ; Rozali:2007rx . Doing a simple minded mapping of this approach to our model, the solitonic solutions should correspond roughly to adding a D0 brane in the configuration. Indeed, noticing that $\sqrt{f}e^{A}=\sqrt{-g_{tt}}$, the first term in (35) takes the form of an action for a D0 brane sitting at the location of the baryon (with a certain $\lambda$ dependent potential).

The final action for a baryon gas with constant density is then obtained by a “convolution” which amounts to integrating the above actions $\int d^{3}x\,n_{b}$ to the actions in (35):

[TABLE]

3.2 Equations of motion and boundary conditions

The complete action of the model is given by the sum of the terms in (3), (12), and (37) in the current approach. We define the bulk charge density as

[TABLE]

The equation of motion for $\Phi$ implies

[TABLE]

The thermal gas solutions extend form $r=0$ (UV boundary) to $r=\infty$ (IR singularity)333Actually, for the choice of action for which we carry out numerical analysis below, the TG solutions turn out to be subdominant.. For these configurations all the charge originates from the baryons, and therefore $\rho(r)=0$ for $r>r_{b}$. Consequently, $\rho(r)=N_{c}n_{b}\equiv\rho_{b}$ for $0<r<r_{b}$. The black hole solutions extend from the boundary ($r=0$) to a horizon at some value $r=r_{h}$ of the bulk coordinate. For them, part of the charge may be hidden behind the horizon. Then $\rho$ is nonzero everywhere, and constant except for the discontinuity at the baryon location: $\rho(r)=\rho_{h}$ for $r_{b}<r<r_{h}$ and $\rho(r)=\rho_{h}+\rho_{b}\equiv\varrho$ for $0<r<r_{b}$.

The gauge field is obtained by inverting (38),

[TABLE]

where

[TABLE]

and integrating over $r$. Here we defined the normalized density as $\rho=M^{3}N_{c}N_{f}\hat{\rho}$. We will discuss below how the constant of integration is fixed.

When $r\neq r_{b}$, the other equations of motion are

[TABLE]

where we already inserted the solution for $\Phi^{\prime}$ from (40). For the thermal gas solutions, $f=\mathrm{const.}$ for $r>r_{b}$. At $r=r_{b}$ we will also need the junction conditions for the various functions which are derived in Appendix C.

The solutions satisfy the following boundary conditions in the UV:

[TABLE]

with $b_{0}=\frac{1}{24\pi^{2}}(11-2x_{f})$ and $\gamma_{0}=\frac{3}{16\pi^{2}}$. The boundary condition for $\lambda$ can here be taken as the definition of $\Lambda$. Dimensionful quantities are computed relative to the energy scale defined by $\Lambda$, and putting $\Lambda$ to a physically reasonable value, one can match with QCD. In all numerical examples considered in this article, we set the quark mass $m_{q}$ to zero.

In the IR, if a black hole is present, the solutions satisfy regularity conditions on the horizon. The thermal gas solution, which has no black hole, is obtained by taking a black hole solution and letting the horizon area approach zero. This is in accordance with the requirement that the IR singularity should be of the “good” type, i.e. it should be possible to clock the singularity with a horizon of infinitesimal area.

3.3 Location of the baryon

In order to determine the location of the baryon, because each baryon carries a fixed charge, we need to study the system in canonical ensemble Bergman:2007wp ; Rozali:2007rx . The Legendre transformed zeroth order action reads

[TABLE]

where the integrand has a discontinuity at $r=r_{b}$ and we also included the CS contribution. In order to determine $r_{b}$, we need to minimize the full action

[TABLE]

varying $r_{b}$. As we keep all sources fixed it is enough to evaluate the derivative with respect to the explicit dependence on $r_{b}$. After the Legendre transformation, this appears in the source term and through the discontinuity of $\hat{\rho}$. Taking the derivative gives the condition (see Appendix C)

[TABLE]

where $g(r_{b}^{\pm})\equiv\lim_{\epsilon\to 0+}\,g(r\pm\epsilon)$ and the averaged derivatives are defined by $\langle g^{\prime}\rangle\equiv(g^{\prime}(r^{+}_{b})+g^{\prime}(r^{-}_{b}))/2=\lim_{\epsilon\to 0+}(g^{\prime}(r_{b}+\epsilon)+g^{\prime}(r_{b}-\epsilon))/2$.

In the limit $\hat{\rho}\to 0$ the first term in (3.3) vanishes as $\propto\hat{\rho}^{2}$, $f$ tends to one, and the derivatives become continuous. Therefore, the condition becomes

[TABLE]

That is, the term defining the soliton mass in (35) should have a minimum in the bulk, and the minimum should be quite close to the boundary: otherwise the baryon will fall deep in the IR in the regime where the coupling of the soliton to the tachyon field becomes important. Our approximation, where the tachyon is essentially neglected will fail in this region.

In order to realize the minimum of the baryon mass at finite $r$ we need to choose the potential $w(\lambda)$ differently from earlier literature (see Arean:2013tja ; Arean:2016hcs ; Jokela:2018ers ). The easiest way444In principle one may find nontrivial solutions to (50) without this property, but numerical studies show that one needs to introduce potentials (e.g. $w(\lambda)$) with peculiar structures, which are likely to cause other issues. to guarantee a minimum is to require that the combination in the square brackets of (50) diverges in the IR (because it does diverge in the UV). If $V_{f0}\sim\lambda^{v_{p}}$ and $w\sim\lambda^{-w_{p}}$, this means that

[TABLE]

Because we need to choose $v_{p}\leq 10/3$ Arean:2013tja we have that $w_{p}\leq 4/3$. In practice, however, the bound is tighter: we have not been able to construct numerically regular backgrounds for potentials with $v_{p}\gtrsim 3$. Here we use $v_{p}=2$, so that $w_{p}\leq 2/3$. We will present some results from the numerical analysis of the model for such potentials with $w_{p}=2/3$ below.

However, as we discussed in Sec. 2.2 and in Appendix A, the phase diagram for this choice of $w_{p}$ has an undesired structure even in the absence of baryons: the thermal gas phase, which is identified as the confined chirally broken phase in QCD, becomes subdominant and is replaced by a phase with “tiny” black hole solutions. The requirement of the TG phase to be dominant leads to $w_{p}\geq 4/3$ which is in contradiction with (51) (inserting $v_{p}\lesssim 3$ which is in turn required for regular backgrounds).

This is, however, a problem with the approximations done in this section rather than a problem in the model: in Sec. 4 we will demonstrate that the coupling of the soliton to the tachyon (which is basically not included in the thin layer approximation scheme) will prevent the soliton from falling in the IR. Therefore, after the coupling to the tachyon has been added, it is actually possible to choose $w(\lambda)$ in the same way as in earlier literature and as required by the analysis of Appendix A. Then a physically reasonable phase diagram with both baryons and a TG phase can be obtained.

3.4 Thermodynamics

The grand potential is given by

[TABLE]

where the source term is

[TABLE]

the Gibbons-Hawking term $S_{\mathrm{GH}}$ is given in Appendix C, and the minus sign in (52) appears because we wrote our actions in Lorentzian signature.

In order to establish the thermodynamics of the system and check its consistency, we need to determine the integration constant in the definition of $\mu$. We reproduce here the basic arguments and delegate details to Appendix C.

As it turns out, if we choose a gauge where $\mu=\Phi(0)$, we need to require that the source term vanishes on-shell by setting

[TABLE]

This fixes the integration constant and ensures that the variation of the free energy follows the first law: if the source vanishes, its variation vanishes as well, and then the variation of the bulk term gives the first law by the standard calculation.

For black hole solutions, (54) can also be derived as the equilibrium condition for moving charge from the baryons to behind the horizon. In order to see this, consider the variation of the Legendre transformed action (3.3) with respect to the charge behind the horizon $\hat{\rho}_{h}$ such that the total charge $\hat{\varrho}=\hat{\rho}_{h}+\hat{\rho}_{b}$ is fixed. This gives

[TABLE]

where we used the EoM (38) on the second line. Therefore, requiring the variation to vanish results in (54).

Naively one might think that the discontinuities in the bulk profile induced by the point like source may also lead to nonzero terms at $r=r_{b}$ which could violate the first law. However, by replacing the delta distribution of the source by a smooth approximation, we see that no such terms can arise: the variation of the Lagrangian is still a total derivative and only boundary terms at $r=0$ are generated. In the point-like limit the junction conditions for the bulk fields at $r=r_{b}$, which are given above, guarantee that all contributions from $r=r_{b}$ in the variation of the grand potential cancel. We show this explicitly in Appendix C.

By using (38), for the thermal gas solution the chemical potential is then given by

[TABLE]

For black hole solutions, (54) sets another constraint which will fix the ratio of the charge of the baryons to the total charge. For the black holes we simply have the formula

[TABLE]

The free energy in each phase can then be obtained by integrating

[TABLE]

where $\varrho$ is the boundary value of the charge density, $\varrho=\rho(0)$, for the TG solutions the entropy vanishes, and for the BH solutions the entropy and the temperature are given by standard BH thermodynamics (see also Alho:2012mh ; Alho:2013hsa ).

3.5 Phase diagram

Using the results of the previous sections, we can compute the phase diagram. The non-baryonic solutions are constructed as in Alho:2012mh ; Alho:2015zua ; Alho:2013hsa . The baryonic solutions are constructed in the same way, i.e. by shooting from the horizon, with the difference that when it is encountered that equation (3.3) is satisfied, a discontinuity is inserted in accordance with the junction conditions derived in Appendix C.1. We choose $x_{f}=1$ and only consider solutions with vanishing quark mass. Both the non-baryonic and the baryonic solutions are then used to integrate the first law as described in Sec. 3.4.

This procedure allows one to compute the phase diagram given a set of potentials. The potentials which are used in the following are given in Appendix B. They have been obtained by fitting to lattice data in the vicinity of $\mu=0$ Jokela:2018ers . Note that in the thin layer approximation, we need to satisfy the bound (51) in order to have stable baryon matter. Consequently we choose a $w$-function with the asymptotics $w\sim\lambda^{-2/3}$ in the IR. The explicit choice is given in (101).555Note that in Appendix B, there are two sets of potentials where the most significant difference between the two choices for $w(\lambda)$ in (101) and in (105) is given by the IR asymptotics. We will use the second set later on. Also note that the $w$-potential was, by tuning the stabilization point, used to set the baryon mass (given by $-S_{\mathrm{DBI}}^{(1)}$ in (35)) to roughly the correct value.

As we pointed out in Sec. 2.2 and in Sec. 3.3, the choice for (101) also has an unintended consequence. Denoting $w\sim\lambda^{-w_{p}}$ for the power $w_{p}$ of the $w$-potential in the IR, one can show that the thermal gas phase is always subdominant at nonzero values of the chemical potential unless $w_{p}\geq 4/3$. This means that a limitation of this approximation is that it is incompatible with having a thermal gas phase at small but nonzero $T$ and $\mu$.666Note that exactly at $\mu=0$, these potentials still display a thermal gas phase, because as $\mu\rightarrow 0$, the size of the black hole shrinks to zero and the geometry approaches that of the thermal gas solution. The reason for this is discussed in more detail in Appendix A. We will use the $w$-potential given by (101) despite the issue with the TG phase, because for this approximation to yield a non-trivial result, we need to have stable baryons.

Moreover this choice of $w$-function destroys the Silver Blaze property of QCD (see Armoni:2014wla ): the pressure is no longer independent of the chemical potential at zero temperature and at nonzero chemical potentials up to a critical value. This happens because the thermal gas depends nontrivially on the chemical potential. It also means that, because the thermal gas phase is replaced by a small black hole, confinement properties are altered. In holography one can investigate confinement properties from the behavior of a string suspended from two points on the boundary. In particular, one can calculate the Wilson loop by computing the on-shell action of such a string Rey:1998ik ; Maldacena:1998im ; Kinar:1998vq ; Bak:2007fk . In this work, we define a phase as confining if the Wilson loop one computes in this way has a branch exhibiting an area law. The crucial difference between the confinement properties of a thermal gas and the confinement properties of the solutions one obtains in this section is that even if the Wilson loop has such a branch, these can correspond to string solutions which are unstable and subdominant to disconnected strings extending to the black hole from the boundary. Note that this instability is a separate effect from the “usual” breaking of the string in QCD due to a pair creation of light quarks. This pair creation effect is not included in the classical string computation we discuss here. Therefore the instability which we observe here is indeed undesirable for QCD. For a more detailed discussion of confinement in a similar geometry, we refer to Gursoy:2018ydr . Note that thermodynamic properties like the equation of state are much less affected by this change in geometry, as they do not probe the deep IR.

There are three order parameters by which we label the phases in the phase diagram:

• •

Chiral symmetry: If the chiral condensate $\langle\bar{q}q\rangle$ is zero, a phase displays chiral symmetry. If it is nonzero, chiral symmetry is broken.

• •

Appearance of baryons: If the $U(1)$ charge associated to quark number777Note that quark number is equal to baryon number divided by the number of colors. is all located behind a black hole horizon, that phase does not have baryons. If, however, a part of the baryon number charge is located at some point in the bulk, we say that a phase is baryonic. Note that even with some of the charge located outside the black hole, part of the charge always remains behind the black hole, making the baryonic phases a sort of mixture of baryons and quark matter.

• •

Confinement: A phase is called confining if the Wilson loop exhibits an area law. As was discussed before, this does not imply confinement in the usual sense. Instead, more weakly, it implies that the properties one usually associate with confinement, like a mass gap and the inability to pull apart two quarks, are only approximately satisfied.

By comparing these properties, there are in principle 8 phases we can have in the phase diagram. Of these, 5 are realized by the model. These are the 4 possibilities for deconfined phases, plus a confined, chirally broken nonbaryonic phase. Comparing free energies of these 5 phases, one obtains the phase diagram in Fig. 1.

Note that as one expects at $\mu=0$, there is a first order phase transition between a chirally symmetric, deconfined QGP phase, and a chirally broken confined phase. As $\mu$ is increased at fixed $T$, the chirally broken phase becomes deconfined as well. Since the geometry changes smoothly as confinement is lost, the transition is a crossover. Also note that the first order transition extending from $\mu=0$ ends in a critical point, turning into a second order phase transition just before the baryonic phases appear. Another critical point can be found between the chirally symmetric deconfined QGP phase and the chirally symmetric deconfined baryonic phase.

Now let us discuss the baryonic phases. The most important observation is their location in the phase diagram. In particular, note that the chirally broken baryonic phase appears at roughly $\mu=270\,\text{MeV}$. This number should be roughly equal to the baryon mass over the number of flavors, which for QCD is (up to the binding energy of nuclear matter) $m_{\text{proton}}/3\approx 313\,\text{MeV}$. This, while not quantitatively the same, is qualitatively in the right range. In particular, one can note that while the potentials were chosen to reproduce the correct baryon mass, the fact that the location of the transition indeed appears in the appropriate location is a non-trivial observation. Another thing to note is that the baryonic phases disappear above a finite value of $\mu$.888The first order transition between the chirally symmetric baryonic phase and the chirally symmetric nonbaryonic phase has a large numerical uncertainty. However, the chirally symmetric baryonic solutions stop existing altogether at a finite value of $\mu$, so even if we cannot pinpoint to great accuracy where the chirally symmetric baryonic phase becomes thermodynamically disfavored, we can say with certainty that this will happen at some finite value of $\mu$. A last observation is that the properties of the chirally symmetric baryonic phase seem somewhat contradictory. On the one hand, chiral symmetry is exact, meaning that there is no mechanism by which the quarks can gain mass. On the other hand, these massless quarks form bound states in the form of baryons. This phase could perhaps be studied in future work by studying the excitations of the theory in that region of the phase diagram.

It is clear that this approximation has its shortcomings, the most serious of which is that from QCD we expect the confined solutions to be described by a thermal gas, while in this approximation we obtain a phase which is not confining in the usual sense. In the next section, we take a different approach, in which these problems are not present.

4 Baryons from a homogeneous bulk gauge field

The approximation considered in the previous section obviously has some shortcomings, as we already pointed out. There is no reason to expect the baryonic soliton to be small in V-QCD, unlike in the WSS model where the size of the soliton $\sim 1/\sqrt{\lambda}$ goes to zero in the limit of strong coupling. Moreover, we basically neglected the tachyon, forcing us to choose a specific kind of potentials which keep the soliton close to the UV boundary where the tachyon is small. This choice of potentials was seen to cause problems with confinement properties. Furthermore, as it turns out, such potentials are slightly disfavored by the fit to lattice data Jokela:2018ers . When the baryon is no longer close to the boundary, the amplitude and/or size of the soliton must actually grow in order to account for the suppression due to the exponential factor $e^{-\tau^{2}}$ in order to keep the baryon charge fixed (see (36) and (2.4)).

If the soliton becomes sizable in the IR, configurations with a high density of solitons, dual to dense baryonic matter in QCD, may be described better by a homogeneous non-Abelian gauge field configuration than separate solitons. This is what we will attempt here. Interestingly, as it turns out, the approach will closely resemble the approximation with a thin noninteracting layer carried out above even though the starting point is completely different. A similar approach has been suggested in the context of the WSS model in Rozali:2007rx and further developed in Li:2015uea (see also Elliot-Ripley:2016uwb for a slightly different setup). We will treat the baryons as probe on top of the TG background in this section. This makes sense as the DBI action discussed above is known precisely only to leading nontrivial order in the non-Abelian field strengths of the solitons.

The basic idea of the setup is as follows. We consider a system with a high density of baryons (comparable to the saturation density on the QCD side) on top of the TG background and divide the space in to three regions in the $r$-direction:

Region close to the boundary, $r\ll r_{c}$, where $r_{c}$ is roughly the location of the soliton “centers”. At high density, the configuration in this region is assumed to be well approximated by a homogeneous baryon field. 2. 2.

The region in the middle, $r\sim r_{c}$. In this region, the configuration is highly inhomogeneous and nontrivial. 3. 3.

The region in the IR, $r\gg r_{c}$. In this region the baryon field is again taken to be homogeneous.

The idea is then that when the baryon density is high, the second, inhomogeneous region is not important for the main features of the phase diagram, and may be ignored. Its effect on the solutions is modeled through a discontinuity of the baryon field, as we shall discuss below. This is an uncontrolled approximation, but as we shall see, the results are encouraging.

4.1 Setup

We will only consider $SU(2)$ solitons in the thermal gas background here. We will use the first order series approximation to the DBI action rather than the full DBI (which is not known for non-Abelian fields). Our ansatz

[TABLE]

respects chiral symmetry and parity Pomarol:2007kr ; Pomarol:2008aa .

The ansatz (59) immediately leads to an issue which has also been observed in the WSS model Rozali:2007rx . Namely inserting it in the expression of the baryon charge gives

[TABLE]

where we reinstated the coefficient $b$. If $h(r)$ is a smooth function this evaluates to a boundary term. Both the UV and IR contributions however vanish: the diverging tachyon sets the action to zero in the IR, and the baryon field $h(r)$ vanishes in the UV due to boundary conditions. Therefore the baryon density is zero.

A nonzero baryon density can however arise from a discontinuity of the function $h(r)$. As discussed above, we use such an abrupt discontinuity of the function $h(r)$ to model the intermediate inhomogeneous regime, which then gives rise to a nonzero baryon density.

We then work out the action for the homogeneous ansatz. The DBI term in (2.3) simplifies to

[TABLE]

with $\Xi$ given in (19). Since the charge is only sourced by the baryons, which are treated as probes, we could also work with only the leading perturbation due to the charge in $\Xi$. We choose however to keep the nonlinear dependence on the charge here.

Importantly, the homogeneous ansatz satisfies the consistency condition (13). Putting this expression together with the DBI action in the absence of solitons and the CS term then gives the total action

[TABLE]

As we argued above, this action should be only trusted away from the vicinity of $r=r_{c}$ where $h$ is discontinuous. In particular we ignore the singular contributions which arise form the derivative of the discontinuities.999For example, for the CS term, this means effectively adding the term S_{\mathrm{Disc}}=\frac{2N_{c}}{\pi^{2}}\int d^{4}x\,\Phi\,e^{-b\,\tau^{2}}(1-2b\,\tau^{2})\ \mathrm{Disc}\,h^{3}\big{|}_{r=r_{c}} which cancels the singular contribution. In general, the prescription which we will use, amounts to interpreting the integrals as

[TABLE]

From the action (4.1), which depends on $\Phi^{\prime}$ through $\Xi$, we derive the charge density

[TABLE]

where we abbreviated

[TABLE]

The $\Phi$ equation of motion implies

[TABLE]

except for $r=r_{c}$ where $h$ is discontinuous. Our prescription stipulates that there are no $\delta$-function contributions at this point, so $\rho$ is continuous. We obtain that

[TABLE]

where $\varrho$ is the density at the boundary, i.e., the physical baryon density. The continuity of $\rho$ implies that it is given in terms of the discontinuity of $h$ as

[TABLE]

where the discontinuity is defined as $\mathrm{Disc}\,g(r)\equiv\lim_{\epsilon\to 0+}\left(g(r+\epsilon)-g(r-\epsilon)\right)$.

4.2 Location of solitons and consistency of thermodynamics

In this subsection we discuss the minimization of the action in particular to determine the location of the discontinuity $r_{c}$. Before going to the precise analysis, we point out how the main features arise from the action we wrote down above. First, as in Sec. 3 we need to work at fixed baryonic charge, rather than chemical potential. Therefore, $\varrho$ in (68) is kept fixed. This means that $\mathrm{Disc}\,h^{3}$ must diverge if $r_{c}$ is taken to deep in the IR (where the tachyon diverges) and also at the point where $2b\,\tau(r_{c})^{2}=1$. This necessarily means that the DBI action and consequently also the free energy diverge at these values of $r_{c}$. In particular, the coupling to the tachyon therefore prevents the baryon from falling in the IR.

Moreover, since (68) gives roughly $\varrho\sim h^{3}$, we expect that the free energy behaves as $F\sim h^{2}\sim\varrho^{2/3}$ at small $\varrho$ and as $F\sim h^{4}\sim\varrho^{4/3}$ at larger $\varrho$. Taking the derivative with respect to $\varrho$, we obtain that $\mu\sim\varrho^{-1/3}$ at small $\varrho$ and $\mu\sim\varrho^{1/3}$ at larger $\varrho$. This indicates that $\mu$ has a minimum and there is a first order phase transition between the empty TG phase and the baryonic phase. Moreover, in the stable phase (larger $\varrho$) the zero temperature speed of sound obeys $c_{s}^{2}=d\log\varrho/d\log\mu\approx 1/3$, i.e., it is roughly given by the conformal value. We will confirm numerically below that the transition is of first order and that the speed of sound is close to the conformal value.

We then move on to the precise analysis. In order to study the system at fixed charge, we perform a Legendre transformation

[TABLE]

Notice that (using our prescription at $r=r_{c}$)

[TABLE]

Therefore

[TABLE]

In order to write down the final expression for the Legendre transformed action we need to eliminate $\Phi^{\prime}$ from the Lagrangian density by inverting (64). When doing this one should recall that we are working in an expansion of the original DBI action at small amplitudes of $F^{(L/R)}$, i.e., at small $\alpha^{\prime}$. Therefore we will consistently ignore all such terms (typically higher powers of $h$) which correspond to higher order terms in the expansion of the DBI action. Doing this gives us the expression

[TABLE]

and allows us to cast the final action in a relatively simple form:

[TABLE]

where

[TABLE]

In order to find the value of $r_{c}$, we need to minimize $\widetilde{S}_{h}$ at fixed $\varrho$. In the numerical analysis carried out below, we will simply do this by numerically minimizing of the action (which is finite for our ansatz so that no regularization is needed) rather than by using the equilibrium conditions explicitly. It is however instructive to compute the conditions.

There are two contributions to the equilibrium condition for $r_{c}$. The first one arises from the variation of $h$ which is necessary as the relation (68) needs to remain satisfied. The second one arises from the discontinuity of the Lagrangian density.101010Recall that due to our prescription that delta functions at $r=r_{c}$ are neglected the two pieces of the Lagrangian for $r<r_{c}$ and for $r>r_{c}$ need to be treated independently. The variation of the action due to changing $r_{c}$ at the discontinuity needs to vanish. The condition can be written as

[TABLE]

evaluated on-shell at $r=r_{c}$, i.e., substituting the regular solution for $h$. Here $\widetilde{L}_{h}$ is the Lagrangian density of the Legendre transformed action $\widetilde{S}_{h}$. This condition is the analogue of the stability condition (3.3) in the thin layer approximation.

The remaining condition is the analogue of the condition (54) and therefore should fix the normalization of the Abelian gauge field $\Phi$. Similarly as in the case of (54), there are two ways to derive it. First, we may require that there are only boundary contributions to the variation of the action. Second, we may replace the thermal gas geometry by a tiny black hole and require that the baryon charges behind the horizon and at $r=r_{c}$ are at equilibrium.

We discuss first the variation of the action. We may keep $r_{c}$ fixed because assuming (75) its variation does not contribute. After this, the solution only depends on the parameter $\varrho$. It appears explicitly in the solution (67) and affects indirectly the profile of $h$. The explicit dependence on $\varrho$ evaluates to an integral of $\Phi^{\prime}$ after using the relation (64). The variation due to the change in $h$ also gives a localized term. We obtain

[TABLE]

The extra terms in this differential vanish if

[TABLE]

We then consider the second condition. Taking the amount of charge behind the horizon to be $\rho_{h}$, the solution (67) is modified to read

[TABLE]

with

[TABLE]

We then require the variation of the action with respect to $\rho_{h}$ near $\rho_{h}=0$ to vanish. The contribution due to the explicit dependence on $\rho_{h}$ again evaluates to an integral of $\Phi^{\prime}$. We obtain

[TABLE]

The condition is therefore very similar to (77), but the sign is opposite. We notice that it is the difference $\varrho-\rho_{h}$ which sources the discontinuity of $h$. In the limit of small density, $\rho$ can be neglected in the equation of motion of $h$ and then $h$ depends on the density only through the discontinuity. Therefore. in this limit, the variations of $\varrho$ and $\rho_{h}$ have exactly opposite effect on $h$. Consequently, the two conditions are equivalent.

At finite density, it appears that the conditions (77) and (80) differ. In the current setup we simply choose to satisfy the first condition, as we are anyhow neglecting backreaction of the charge and higher order corrections in $h$ which are important at larger charge densities. Therefore, within the range of consistency of our approach, the conditions are equivalent, and we will simply ignore the disagreement at high density. Alternatively, the system could be stabilized toward exchanging charge between the bulk and $r=\infty$ by adding some extra charge at $r=r_{c}$ (and naturally also by adding it at $r=\infty$).

4.3 Asymptotic behavior

We then discuss the asymptotics of the field $h$. It is straightforward to solve the asymptotics in the UV. After inserting the standard UV behavior in the action (73) and solving for $h$, we find the asymptotics typical for gauge fields: $h\simeq C_{1}+C_{2}r^{2}$. We require that non-Abelian sources vanish, and therefore $C_{1}=0$, but $C_{2}$ remains as a free parameter.

The IR behavior is more complicated. First, we notice that for $b>1$ in (67) and because the tachyon diverges faster than $r$ in the IR, the factor $\rho/(V_{\rho}we^{-2A})$ tends to zero in the IR unless $h$ diverges exponentially fast. We will assume that this is the case and analyze the resulting EoM, which is obtained after neglecting the terms involving $\rho^{2}$ in the action:

[TABLE]

First of all, $h=0$ is an exact solution to the full action. There is also a family of regular solutions

[TABLE]

which vanish in the IR. Here $C_{\tau}$ is a positive constant which can be expressed in terms of the IR expansion parameters of the potentials. One may check numerically that there are no other asymptotic solutions (for which tachyon would diverge so fast that the condition of $\rho/(V_{\rho}we^{-2A})$ vanishing in the IR would be violated).

The IR solutions in principle bring in one more parameter $h_{0}$ which also needs to be determined by minimizing the action. We have however checked that in the numerical analysis below the value of $h_{0}$ at the minimum is consistent with zero, i.e., in practice the solution in the IR region ($r>r_{c}$) simply vanishes.

4.4 Phase diagram and the equation of state

We work in the probe limit: we first construct the TG background solutions Jarvinen:2011qe in the absence of baryons, setting the quark mass to zero. The equations of motion and boundary conditions for the background fields $A$, $f$, $\lambda$ and $\tau$ are the same as those described in Sec. 3.2. We then insert the background in (73) and solve the resulting equations of motion for the baryon field. We explicitly integrate (73) to obtain the free energy in the canonical ensemble, where we subtract $\tilde{S}_{h}$ evaluated for $h=0$.111111This subtraction makes sure that the grand potential is exactly equal to the Legendre transform of the free energy, not just up to an additive constant. We then need to minimize this free energy over the free parameter that we left as boundary condition, namely $C_{2}$ as defined in the previous subsection. Subsequently, we do a Legendre transform to obtain the grand potential, which is used to compute the phase diagram, and which is equal to minus the pressure. Note that in all the subsequent results, we set $N_{c}=3$.

The precise choice of the various potential functions which we use in this section is given in Appendix B. Notice in particular that we choose the asymptotics of the $w$-function such that $w\sim\lambda^{-4/3}$ in the IR, with the explicit choice of the function given in (105). As we pointed out above in Sec. 3.3 and proved in Appendix A, we must choose the asymptotics $w\sim\lambda^{-w_{p}}$ with $w_{p}\geq 4/3$ in order to have a stable TG phase at small, nonzero $T$ and $\mu$, which is in turn necessary for the approach of this section to work. Moreover notice that in practice it is impossible to use the same sets of potentials in the two approximation schemes considered in this article: in the thin layer approximation of Sec. 3 we were forced to use $w\sim\lambda^{-w_{p}}$ with $w_{p}\leq 2/3$ in order to stabilize the baryon phases.

In the final results, we also compare the free energy of the TG and baryonic TG phases to the chirally symmetric, baryonless BH solutions at $x_{f}=1$, which are constructed as in Alho:2013hsa . The latter solutions therefore model the QGP, and the equation of state in this phase is that anchored to lattice data in Jokela:2018ers since the choice of potentials is the same. The comparison of the free energies between the baryonic and QGP phases however comes with uncertainties, which mostly arise from the normalization of the baryonic free energy. This is due to the approximations done in this section and the fact that the baryonic ansatz assumes $SU(2)$ flavor symmetry whereas the QGP equation of state was fitted to data with 2+1 dynamical quarks. In particular our result for the pressure in the baryonic phase may be somewhat low because of the simple approach.

The phase diagram one obtains is shown in Fig. 2. It features the following three phases, corresponding to the three kinds of numerical solutions we discussed above:

• •

A confining phase with broken chiral symmetry. In contrast to the confining phase obtained in the thin layer approximation, this phase is a thermal gas phase, which implies that, for instance, glueballs are absolutely stable, not just long-lived like they were in the thin layer approximation.

• •

A confining phase with broken chiral symmetry and condensed baryons, i.e., the new ingredient from the approach of this section.

• •

A deconfined QGP phase with chiral symmetry.

The pressure in the first phase in our approach is constant Alho:2013hsa , and it is independent of the temperature also in the baryonic phase. Therefore the baryon-vacuum transition line is exactly vertical. Temperature dependence in these phases would arise from stringy loop effects Alho:2015zua . In the QGP phase, nontrivial temperature effects are included, and the result for the pressure can be seen as an extrapolation from the fit to lattice data around $\mu=0$ as we have explained above Alho:2013hsa ; Jokela:2018ers .

The location of this transition depends on the parameter $b$ described above. In the results described in this section, $b=10$ was chosen to have the transition at approximately $\mu=313\,\text{MeV}\approx m_{\text{proton}}/3$ (where we ignored the small binding energy of nuclear matter). Another thing to note is that, unlike within a similar approximation in the WSS model121212A phase transition at a large value of the chemical potential is obtained in the WSS model in approaches using interacting solitons Li:2015uea ; BitaghsirFadafan:2018uzs . Li:2015uea , the baryonic phase does not survive to arbitrarily large values of the chemical potential. Lastly, because in the thermal gas phase temperature dependence is suppressed, the baryonic phase transition does not end at a critical point, but instead continues until the QGP phase becomes dominant.

In Fig. 3, the pressure at $T=0$ is plotted. The phase transitions are clearly visible. The latent heat at the vacuum to nuclear matter transition is $\Delta\epsilon\approx 51\ \mathrm{MeV}\,\mathrm{fm}^{-3}$, and at the nuclear matter to quark matter transition $\Delta\epsilon\approx 687\ \mathrm{MeV}\,\mathrm{fm}^{-3}$. In principle, this pressure could be used as an equation of state to solve the Tolman-Oppenheimer-Volkov equations, thereby obtaining a mass-radius relation for neutron stars. (See, for instance, Kim:2011da ; Ghoroku:2013gja for the application in the WSS model.) We leave this for future work.

The quark number density is shown in Fig. 4. From the fact that the number density jumps across the phase transitions, one can clearly see that the phase transitions are indeed first order. The baryon number density (which is obtained from the quark number density by dividing by $N_{c}=3$) at the end of the transition from the vacuum to nuclear matter is $n_{b}\approx 0.056\ \mathrm{fm}^{-3}$, i.e., $n_{b}\approx 0.35n_{s}$ where $n_{s}\approx 0.16\ \mathrm{fm}^{-3}$ is the nuclear saturation density. This indicates that the nuclear matter number densities and pressures as a function of the chemical potential lie below the band of nuclear matter EoS (obtained by extrapolating the low density equations of state to higher densities Hebeler:2013nza ; Annala:2017llu ). This may be at least in part due to the uncertainties in the normalization of the baryon action (73) and in the rough approximations done in this section. We comment more on this in Sec. 5. We also note that the baryon to quark matter transition is strongly first order: the density jumps at the transition roughly by a factor $4.5$.

Figure 5 shows the location of the discontinuity in the bulk. The soliton stays roughly in the same place for different values of the chemical potential. Also, in units of $\Lambda_{\text{QCD}}$, $r_{c}$ is roughly $\mathcal{O}(1)$, as one would expect.

The isothermal speed of sound is shown in Fig. 6. The vacuum (TG) phase does not have a speed of sound, as it has no pressure and no energy density. At the vacuum to baryon matter transition the speed of sound immediately jumps to a value, which is larger than expected for nuclear matter in this regime. Such a deviation is not surprising because the homogeneous approximation used here is expected to be reliable only at higher densities. The change of the speed of sound at the baryon to quark matter transition is relatively large, indicating a jump from a hard to soft equation of state at the strongly first order transition. Noteworthy is that the speed of sound is above the conformal value $c_{s}^{2}=1/3$ in two intervals of the chemical potentials. For large values of the chemical potential, it approaches the conformal value from above. Notice that it is likely that speeds of sound above the conformal value are necessary in order to pass the astrophysical constraints from observations of neutron stars and their mergers (see, e.g., Bedaque:2014sqa ; Tews:2018kmu ). Interestingly, the baryonic equation of state in the WSS model has been seen to also violate the conformal bound clearly BitaghsirFadafan:2018uzs . For other work towards realizing stiff phases in holographic models for dense QCD see Hoyos:2016cob ; Ecker:2017fyh .

Lastly, the adiabatic index, which is defined as $\Gamma=\mathrm{d}\log p/\mathrm{d}\log n$, is plotted in Fig. 7. Similarly as the speed of sound, this quantity measures the stiffness of the equation of state, and a piecewise constant $\Gamma$ is often used in polytropic realizations of the equation of state. Its order of magnitude in the baryonic phase is comparable with what is expected from nuclear matter models (see, for example, Douchin:2001sv ; Hebeler:2013nza ).

5 Discussion and Outlook

In this work, we explored how baryonic physics can be included in the V-QCD model. Baryons play a large role in the cold dense matter region of the QCD phase diagram. Neutron stars consist of matter located in this region of the phase diagram, which is one of the main motivations for this work. In holography, baryons are dual to solitons sitting in the bulk. Ideally, one would obtain these solitons explicitly. This is unfortunately challenging especially in bottom-up models with potentials that are phenomenologically matched to QCD. In this work, we tried to get around this problem by two different approximations: one in which the baryons are approximated as a thin layer of noninteracting solitons, and one in which they are described by a homogeneous non-Abelian gauge field configuration. These approximations are rough, and likely to miss some features of real baryon dynamics, but still we obtained several encouraging results.

Let us first discuss the thin layer approximation studied in Sec. 3. This approximation has baryonic matter appearing in roughly the expected region in the phase diagram in Fig. 1. Also, at large values of the chemical potential, the baryonic matter eventually makes way for a non-baryonic deconfined phase. However, this approximation also has some unrealistic aspects, such as the appearance of a chirally symmetric baryonic phase. The most serious of these issues is the stability of the location of the solitons. In particular, the solitons fall into the deep IR unless in the asymptotic IR behavior of $w\sim\lambda^{-w_{p}}$, we choose $w_{p}\leq 2/3$. Such a power law has other consequences though, the most important of which is that at finite chemical potential the thermal gas is replaced by a small black hole phase which approaches a thermal gas solution as $\mu\rightarrow 0$. This means that at finite chemical potential, the theory is only confining up to a certain distance scale, and that glueballs are no longer stable, but instead have a very long lifetime.

This issue is not present in the other approach we studied in Sec. 4, namely the homogeneous approximation. In this approach, the solitons naturally stabilize at a coordinate distance from the boundary of order $\mathcal{O}(1)$. In fact, as we discussed in Sec. 2.2, in order to satisfy phenomenological constraints the power law in the $w$ asymptotics has to be exactly equal to $-4/3$.

In the homogeneous approximation, the phase diagram (see Fig. 2) is qualitatively similar to what one would expect (within the limitations of the approach). For small temperatures and small chemical potentials, there is a thermal gas phase, which displays confinement, a linear glueball spectrum and broken chiral symmetry. At larger temperatures, there is a phase transition to a deconfined phase where chiral symmetry is restored. If instead of increasing the temperature one increases chemical potential, a phase of non-zero baryon density appears, which is still confining and has broken chiral symmetry. Increasing the chemical potential further, the baryonic phase gives way to a deconfined chirally symmetric plasma phase.

There is a hint that there is something going on in that region: at large chemical potential, as one takes $T\rightarrow 0$, the geometry becomes asymptotically $\text{AdS}_{2}$. This can be understood as a signal of a quantum critical regime Alho:2013hsa , which may be subject to instabilities. In this region the entropy is finite even at zero temperature. This could indicate that there is an operator missing in our bulk effective field theory, and the end point of the instability induced by this operator is another phase not included in the present analysis. It is tempting to interpret this as an instability towards an exotic, i.e., color superconducting phase in QCD. It would be interesting to understand if there are connections to the recent work on color superconductivity in other models Faedo:2017aoe ; BitaghsirFadafan:2018iqr ; Faedo:2018fjw .

In addition to the phase diagram, we computed several thermodynamical observables. The most striking aspect of these is that the speed of sound clearly exceeds the value of conformal theories. This result thus joins the list of examples where this happens Hoyos:2016cob ; Anabalon:2017eri ; Ecker:2017fyh . From the point of view of neutron star physics, this result is not unexpected though, since experimental evidence seems to indicate that speeds of sound above the conformal value are necessary Bedaque:2014sqa ; Tews:2018kmu .

For future work, it would be very interesting to apply our results to neutron star physics, following the ideas of Hoyos:2016zke ; Annala:2017tqz ; Jokela:2018ers . In particular, one could use the equation of state as input for the Tolman-Oppenheimer-Volkov equations to obtain a mass-radius relation for non-rotating neutron stars. Also, one could use the equation of state to simulate neutron star mergers. Interestingly, our equation of state includes a first order transition between the baryon and quark matter phases at low densities which could be easily reachable in merger events. Actually, the densities which we obtained in this model are even too low to be consistent with extrapolations of equations of state from nuclear matter and the bound of the maximal neutron star mass from Shapiro delay measurements. This value is affected, among other things, by the normalization of the baryon action (73). This normalization contains uncertainties due to the inhomogeneities neglected in Sec. 4, and because we restricted to the ansatz with $SU(2)$ flavor symmetry which is also a rough model in the Veneziano limit. Improving on these issues may increase the pressure in the baryonic phase pushing the transition from baryons to quarks to higher chemical potentials. Note that such a change is also likely to make the excess over the conformal value for the speed of sound more drastic, see Fig. 6. We could also simply treat the normalization of the action (3.3) as an additional fit parameter that could be determined e.g. by comparing to the nuclear saturation density as discussed in Sec. 4.4.

In light of neutron stars, it would also be interesting to expand the model to study the equation of state for neutron stars in more detail. Currently, the V-QCD model that we are considering has a flavor sector consisting of $N_{f}$ identical massless quarks, all with charge $+1$. In the interior of neutron stars, however, the number of neutrons is expected to be greater than the number of protons. In principle, in the V-QCD model one could also differentiate between different quark species, giving each a different mass and couplings to external gauge fields.131313Note that to incorporate the effect of multiple quark species onto the plasma phase, one would have to include the backreaction of each quark species. This would require that for each species $i$, there are $N_{f,i}$ identical copies. This would then allow to set the charges of the different species appropriately, as well as to introduce an isospin chemical potential, which could force an imbalance between the quark species to mimic that in a neutron star.

One could also, since neutron stars are known to have strong magnetic fields, study the effect of a magnetic field on the baryonic matter. Note that in Drwenski:2015sha ; Gursoy:2016ofp ; Gursoy:2017wzz , the effects of a magnetic field on the plasma phase was already studied. One could in principle then do a similar analysis as was done in this work to include the effect of the magnetic field on the baryonic matter too. This has been studied in the WSS model in Preis:2011sp . Another interesting extension would be to study the transport properties Gursoy:2009kk ; Iatrakis:2014txa of dense nuclear matter. Dissipation plays a role even in isolated neutron stars, and (at least) bulk viscosity is expected to be relevant in neutron star mergers Alford:2017rxf .

Another thing that would be interesting is to explore different choices for the CP-odd potential. In this work, we chose $V_{a}(\lambda,\tau)=\exp(-b\tau^{2})$ with $b=10$, but in principle one could even choose a different functional form entirely, possibly including dependence on $\lambda$. This will require establishing the CS terms of Sec. 2.4 for the more general choice of potential. It would be interesting to study the effects of different choices on the baryon physics. One could then try to match the potential to match nuclear matter models, or even the properties of neutron stars once more astronomical observations become available.

Another possible improvement to the approximations employed in this work would be to try to include the effect of backreaction of the homogeneous bulk gauge field $h$ onto the geometry. This is in principle possible, but it complicates the numerical analysis, therefore we leave this for future work. A similar possible improvement would be to include more terms in the DBI expansion using the results in the literature Refolli:2001df ; Koerber:2001uu ; Grasso:2002wb ; Keurentjes:2004tu . The full non-Abelian DBI action is not known, but one may use the action in (4) with the standard trace or by using the symmetrized trace prescription. Moreover, one could check whether the related approximation scheme suggested in Elliot-Ripley:2016uwb , which is essentially based on applying the homogeneous assumption directly to the field strengths $F$ rather than the gauge fields, leads to any relevant changes in the results.

A final thing that would be interesting, if it can be done, is to attempt to obtain the soliton solutions explicitly. This would certainly be worthwhile, for the following reason. As was mentioned in the beginning of Sec. 4, the homogeneous approximation implicitly assumes a high density of baryons. This means that the approximation is not well suited to the low densities one expects near the vacuum to nuclear matter phase transition. In this regime description in terms of weakly interacting solitons could be more appropriate. If soliton interactions can be taken into account, they can then describe both low densities and high densities equally well.

Acknowledgments

We thank for discussions F. Bigazzi, C. Ecker, U. Gürsoy, C. Hoyos, I. Iatrakis, N. Jokela, E. Kiritsis, A. Krikun, M. Lippert, F. Nitti, W. van der Schee, A. Schmitt, J. Sonnenschein, and J. Tarrio. We also thank C. Ecker, N. Jokela, E. Kiritsis, W. van der Schee, and A. Schmitt for comments on the draft. This work is partially supported by the Netherlands Organisation for Scientific Research (NWO) under the VIDI grant 680-47-518, and the Delta-Institute for Theoretical Physics ($\Delta$-ITP), both funded by the Dutch Ministry of Education, Culture and Science (OCW). The work of T. I. is supported by JSPS KAKENHI Grant Number 18H01214.

Appendix A Asymptotics of $w(\lambda)$ and the phase structure

As it turns out the IR asymptotics of the coupling function $w(\lambda)$ of the gauge field may affect drastically the phase diagram of the model, which leads to a constraint for the IR behavior of this function which is completely independent of the presence of baryons. The solutions to analyze in order to see this are the small temperature near-extremal tachyonic black holes. In particular, as it turns out, it is important to compute what is the chemical potential in such backgrounds in the limit of small black holes. In this limit the chemical potential may either tend to zero or to infinity. The former option means that tachyonic black holes exist at arbitrary small chemical potentials and temperatures, and they will also dominate over the thermal gas solutions (which do not have a black hole horizon). Therefore the phase diagram has an undesired structure where the thermal gas phase (identified as the chirally broken confined phase of QCD in Alho:2012mh ; Alho:2013hsa ) is subdominant for all positive $\mu$. The dominant black holes are very small (as we shall see below) and the backgrounds are close to the thermal gas, so that the change in thermodynamics with respect to the thermal gas solution is tiny. However, the dominance of the black holes mean also that the Silver Blaze property of QCD (see Armoni:2014wla ) is lost. As the result differs qualitatively from that obtained from QCD, the dominance of the black holes at small $T$ and $\mu$ is clearly disfavored.

We shall then analyze for which potentials this phenomenon occurs. Many of the details will be based on numerical analysis, and we will not present rigorous proofs. We assume that the IR behavior of the potentials is

[TABLE]

up to logarithmic corrections in $\lambda$. The parameters must satisfy the following constraints: $4/3\leq v_{p}\leq 10/3$ and $w_{p}\leq 4/3$ Jarvinen:2011qe ; Arean:2013tja . Further we will take the logarithmic corrections to $V_{g}$ and $\kappa$ to be those singled out in Gursoy:2007cb ; Gursoy:2007er ; Jarvinen:2011qe ; Arean:2013tja . After this, our potentials satisfy all qualitative IR constraints established in earlier work. The potentials used in the numerical analysis of this article (given in Appendix B) belong to this class of potentials.

Let us first recall what is the “standard” behavior of the background in the IR in the absence of charge and black hole horizon, i.e., for the thermal gas solutions, for the specified potentials. The IR behavior in the gluon sector is given by Gursoy:2007cb ; Gursoy:2007er

[TABLE]

For the tachyon we have

[TABLE]

where $C_{\tau}>1$ is a known constant.

As we shall show, the structure of the relevant small BH backgrounds has the following scaling regimes:

• •

For $r\Lambda\ll 1$, the background follows the “standard” UV behavior established in Gursoy:2007cb ; Gursoy:2007er ; Jarvinen:2011qe . The UV regime will be irrelevant for the thermodynamics.

• •

For $r_{*}\Lambda\gg r\Lambda\gg 1$, where $r_{*}$ will be specified below, the background follows the “standard” IR behavior.

• •

For $r_{h}\gg r\gg r_{*}$, the metric will follow the standard IR behavior but the tachyon will be frozen. Here $r_{h}$ is the value of $r$ at the horizon.

The UV regime will be irrelevant for our analysis, so we restrict to the IR regime $r/\Lambda\gg 1$. The backreaction of the flavor to the glue is determined by the effective potential Jarvinen:2011qe ; Alho:2013hsa :

[TABLE]

where

[TABLE]

For the relevant solutions the charge $\hat{\rho}/\Lambda^{3}$ will be tiny so that for $r\sim 1/\Lambda$ the charge is decoupled, i.e., $K\ll 1$. Therefore we will be considering the limit $\hat{\rho}\to 0$. As $r$ grows, the solution behaves as the thermal gas and follows the asymptotics given in (84) and in (85). In particular the growth of the tachyon decouples the flavor from the glue. Because of the exponential dependence of $K$ on the tachyon, its growth will result in $K$ being of the order of one at some location, which we mark by $r_{*}$. Since the tachyon dependence dominates in $K$, we have roughly $\hat{\rho}\sim e^{-\tau(r_{*})^{2}}$. Consequently,

[TABLE]

As $r$ grows further, $K$ will also grow. This is seen as follows. At large $K$ the tachyon EoM in (42) becomes

[TABLE]

so that the term in the square brackets is a constant. We note that $f\kappa(\lambda)/w(\lambda)$ is suppressed in the IR: this is clear for $w_{p}<4/3$, and must also hold for $w_{p}=4/3$ for the spectrum to agree with QCD Arean:2013tja . But $\tau^{\prime}/G$ is bounded in no matter how fast the tachyon grows. Therefore the constant must be zero, and the regular solution for the tachyon is simply the constant solution. Then using the asymptotics (84), $K$ grows towards the IR if $w_{p}>v_{p}-2$. This equation is satisfied for the potentials we used in this article, and we will also assume this in the analysis below.141414Our results will be valid also when this inequality is not satisfied, but the analysis will be more complicated.

When $K\gg 1$, the effective potential reads

[TABLE]

Since $w_{p}>v_{p}-2>-2/3$, the second term will grow towards the IR and eventually become comparable to $V_{g}$. This marks the regime where the regular solution must develop a horizon: we may verify numerically that solutions with horizon at even higher values of $r$ do not exist. The solution in the vicinity of the horizon is not analytically tractable, and consequently we cannot estimate the temperature of the solutions analytically, but we may check numerically that as $r_{h}$ approaches its highest value $T/\Lambda$ tends to zero. From the requirement of the vanishing of (90) we obtain that

[TABLE]

This shows that we have the hierarchy assumed above, $r_{h}/\Lambda\gg r_{*}/\Lambda\gg 1$, in the limit $\hat{\rho}\to 0$.

It remains to check the behavior of the chemical potential for these solutions. It is given by

[TABLE]

The integrand $-\Phi^{\prime}(r)$ increases fast with $r$ in the regime $r_{*}\Lambda\gg r\Lambda\gg 1$ due to the tachyon dependence. When $r_{h}\gg r\gg r_{*}$ we find that

[TABLE]

where we used the fact that for the potentials and tachyon asymptotics specified above, $G$ has mild dependence on $r$ which we ignore here. For $w_{p}<4/3$ we find that $-\Phi^{\prime}(r)$ decreases with $r$ in this region. Therefore the integral (92) is dominated at $r\sim r_{*}$, and we find

[TABLE]

up to subdominant multiplicative corrections. For $w_{p}<4/3$ (but not for $w_{p}=4/3$) and given (88), we therefore find that $\mu\to 0$ as $\hat{\rho}\to 0$. Combining this with the numerical observation that these black hole solutions may also have an arbitrarily small temperature, we conclude that tiny tachyonic and charged black holes exist at arbitrarily small temperatures and chemical potentials.

Since $r_{*}\to\infty$ as $\hat{\rho}\to 0$, the background approaches the thermal gas background pointwise in this limit. As the IR contributions to the grand potential are suppressed by the behavior of the metric and/or the tachyon, the value of grand potential for the black holes approaches that of the thermal gas. As $\hat{\rho}>0$ from the first law of thermodynamics it immediately follows that the black holes are dominant over the thermal gas solutions when $w_{p}<4/3$ (as was also seen numerically in Fig. 1).

Therefore we conclude that a reasonable phase diagram can only be obtained for $w_{p}=4/3$: whether $\mu$ tends to zero or infinity in the zero charge limit depends on subleading corrections to (94). Numerically we have verified in this article and in Alho:2013hsa that for $w_{p}=4/3$ there are indeed choices of potentials which give the desired phase structure (i.e., no tachyonic black holes at small $T$ and $\mu$).

Appendix B Choice of potentials in the V-QCD action

We used two sets of potentials in this article, one set for the thin layer approximation in Sec. 3, and another for the homogeneous baryons of Sec. 4. Notice that due to the constraints for the asymptotics of $w(\lambda)$ explained in Sec. 3.3, we cannot use the same set of potentials with both approaches. The sets of potentials were chosen to satisfy various asymptotic constraints (see, e.g., Arean:2013tja ; Arean:2016hcs ) and fitted to lattice data for QCD thermodynamics as explained in Jokela:2018ers . The latter set equals the set 7a in Jokela:2018ers up to small changes in the potential $\kappa(\lambda)$ which leave the thermodynamics of the deconfined phase unaffected. For the former set, where we chose the function $w$ to have the asymptotics $w\sim\lambda^{-2/3}$ as we discussed in the main text, the fit is somewhat worse than for the latter set. That is, the lattice data favors the asymptotics $w\sim\lambda^{-4/3}$.

Notice that while the number of free parameters is high, the fit to lattice data is “stiff”: after fixing the asymptotic behavior of the potentials, observables mostly depend very mildly on the details of the potentials in the regime $\lambda=\mathcal{O}(1)$. It is in this regime where most of the remaining freedom is, and it is controlled by the fit parameters discussed here.

For the potentials $V_{g}$, $V_{f0}$, and $\kappa$ we used the following ansatz in both cases:

[TABLE]

Here the UV parameters were given by

[TABLE]

where we set $x_{f}=N_{f}/N_{c}=1$, and chose $W_{0}=0$ for the first set (thin layer of noninteracting baryons) and $W_{0}=2.5$ for the second set (homogeneous field). The parameters for the glue potentials were chosen to be

[TABLE]

for both sets of potentials. They were obtained by fitting the lattice data Panero:2009tv for the thermodynamics of pure Yang-Mills theory Alho:2015zua ; Jokela:2018ers .

The remaining parameters (mostly governing the potentials in the IR) and the function $w(\lambda)$ are different for the two potential sets.

The first set, used in Sec. 3, has the following ansatz for $w(\lambda)$

[TABLE]

where the UV behavior was chosen to produce qualitatively correct thermodynamics at large $T$ and small $\mu/T$ Jokela:2018ers , and the IR asymptotics $w\sim\lambda^{-2/3}$ was chosen such that the baryon remains close to the UV boundary. The full set of UV and/or normalization parameters is

[TABLE]

Moreover, when $W_{0}=0$ the AdS radius $\ell=1$. The IR parameters are

[TABLE]

These parameters were chosen such that the EoS at $\mu=0$ and its first nonzero cumulant compare relatively well with the lattice data. As we remarked above, due to the requirement $w\sim\lambda^{-2/3}$, the fit is a bit worse than for the potentials considered in Jokela:2018ers (and therefore also worse than for the other set given below).

The other set, used in Sec. 4, is the “potentials 7a” constructed in Jokela:2018ers up to small modifications in the choice of $\kappa(\lambda)$ which do not affect the thermodynamics of the deconfined quark matter phase. We choose

[TABLE]

The UV parameters are given by

[TABLE]

with the AdS radius $\ell=\left(1-2.5/12\right)^{-1/2}$. The IR parameters are

[TABLE]

Appendix C Discontinuities and junction conditions for the thin layer approximation

In this Appendix we consider several technical details of the thin layer approximation of Sec. 3.

C.1 Junction conditions

As we are backreacting the baryons to the metric, we need the Israel junction conditions for the bulk fields at the location of the baryon. We will consider here the general case where part of the charge sits behind the horizon and there is additional charge due to the baryon. That is, $\hat{\rho}=\hat{\rho}_{h}$ for $r>r_{b}$ and $\hat{\rho}=\hat{\rho}_{h}+\hat{\rho}_{b}$ for $0<r<r_{b}$. Since the source is independent of the tachyon, its junction condition is obtained by requiring the continuity of the term in the square brackets in (42). Denoting $f(r_{b}^{\pm})=\lim_{r\to r_{b}\pm}f(r)$, we find

[TABLE]

where $K(r_{b}^{-})=(\hat{\rho}_{h}+\hat{\rho}_{b})^{2}K(r_{b}^{+})/\hat{\rho}_{h}^{2}$. The source term however does depend on the metric and the dilaton, so its contribution to the Einstein equations needs to be considered. We notice that the dependence on the metric is through the factor $\sqrt{-g_{tt}}=\sqrt{f}e^{A}$. Consequently, the constraint equation (43) is unchanged, but source terms are generated in the other equations:

[TABLE]

These equations imply the following discontinuities of the derivatives

[TABLE]

where

[TABLE]

The discontinuity of $\lambda$ can be read from its second order equation, which reads

[TABLE]

where

[TABLE]

This implies that

[TABLE]

C.2 Baryon location

Having solved the discontinuities we then derive the equilibrium condition (3.3). In order to compute the variation of the first term in (3.3), we first need to analyze the discontinuity of the tachyon more carefully. For clarity we replace the $\delta$ function by a continuous estimate $\delta_{\epsilon}(r-r_{b})$ and a continuous step function satisfying $\theta_{\epsilon}^{\prime}(r-r_{b})=\delta_{\epsilon}(r-r_{b})$. Then we can write that

[TABLE]

However, the only property of $K$ that we will need below is the fact that near the discontinuity

[TABLE]

since the terms where the derivative operates on the step function dominate.

The tachyon behavior near the discontinuity becomes

[TABLE]

The discontinuous behavior in the (derivative of the) first term of (3.3) is included in

[TABLE]

Notice that the discontinuity of the tachyon in $G(r)$ is a reaction to the explicit $r_{b}$ dependence of the charge density and the derivative should not act in this discontinuity. Integrating over $r$ yields

[TABLE]

The second term in (3.3) can be treated in a more straightforward way, noticing first that the continuous estimates for the discontinuities of $\lambda^{\prime}$, $A^{\prime}$, and $f^{\prime}$ are simply obtained by replacing $\delta\mapsto\delta_{\epsilon}$ (which is the case because the coefficients of the $\delta$-functions in (C.1) and (C.1) are continuous). After partial integration, the contribution from the second term is

[TABLE]

Taking into account the $\epsilon$ regularization, one immediately obtains the expression in (3.3) as the $r$-integral leads to the averages of the derivatives over the discontinuity.

As a consistency check, we show that the equilibrium condition follows from the consistency of the equations of motion. Namely, the constraint equation (43) can be written as

[TABLE]

where

[TABLE]

Moreover, using

[TABLE]

when evaluating the discontinuity of the left hand side in (124), we obtain

[TABLE]

which is the same condition as (3.3).

C.3 Variation of the on-shell action

Finally, we check explicitly that the contributions to the first law of thermodynamics from $r=r_{b}$ are absent. To this end we consider a generic variation of the on-shell action around a saddle-point configuration. As we already pointed out, the source term does not contribute to the variation thanks to the condition (54). As usual, the variation of the bulk term becomes a total derivative. In our case, the relevant variation terms read

[TABLE]

where we readily ignored the tachyon term which will not contribute at $r=r_{b}$. In order to make the variation of the gravitational action well-behaved we also need to add the Gibbons-Hawking term

[TABLE]

Adding the variation of this term, the result reads

[TABLE]

where the subscript “disc” refers to the fact that we are only keeping the terms which are potentially discontinuous at $r=r_{b}$.

It suffices to show that the expression in the square brackets of (131) is in fact continuous at $r=r_{b}$. To do this, we need to consider the variation of (54). This leads to two kind of contributions: one due to the variation of $r_{b}$ and the other due to the variation of the fields. The former was actually already computed above in Sec. C.2 using the Legendre transformed action. In order to see this explicitly, we again consider a continuous estimate of the delta function, so that the variation of the left hand side of (54) is interpreted as (omitting the trivial factor $\delta r_{b}$)

[TABLE]

Inserting here the solution of $\Phi^{\prime}$ from (38) and $G(r)$ from (119) we find that

[TABLE]

The variation of the right hand side of (54) is computed as in (C.2) so the variations of $r_{b}$ cancel after imposing (3.3). Therefore we are left with the variations of the fields in (54), which lead to

[TABLE]

We then find for the discontinuity in (131)

[TABLE]

Collecting the results, we have shown that the variation of the on-shell action only receives contributions at the boundary.

As a final remark, we notice that the Gibbons-Hawking term

[TABLE]

also contains a term localized at $r=r_{b}$ as $L_{\mathrm{GH}}$ is discontinuous at this point.