Holographic bottomonium formation in a cooling strong-interaction medium at finite baryon density

The shrinking of the bottomonium spectral function towards narrow quasi-particle states in a cooling strong-interaction medium at finite baryon density is followed within a holographic bottom-up model. The 5-dimensional Einstein-dilaton-Maxwell background is adjusted to lattice-QCD results of sound velocity and susceptibilities. The zero-temperature bottomonium spectral function is adjusted to experimental $\Upsilon$ ground-state mass and first radial excitations. At baryo-chemical potential $\mu_B = 0$, these two pillars let emerge the narrow quasi-particle state of the $\Upsilon$ ground state at a temperature of about 150 MeV. Excited states are consecutively formed at lower temperatures by about 10 (20) MeV for the $2S$ ($3S$) vector states. The baryon density, i.e. $\mu_B>0$, pulls that formation pattern to lower temperatures. At $\mu_B =$ 200 MeV, we find a shift by about 15 MeV.


I. INTRODUCTION
The observation of sequential bottomonium suppression [1][2][3][4][5] in relativistic heavy-ion collisions at LHC has sparked a series of dedicated investigations, e.g. [6][7][8][9][10][11][12][13][14]. Such heavyquark flavor degrees of freedom receive currently some interest as valuable probes of hot and dense strong-interaction matter produced in heavy-ion collisions at LHC energies. The information encoded, e.g. in heavy quarkonia (QQ = cc or bb) observables, supplements penetrating electromagnetic probes and hard (jet) probes and the rich flow observables, thus complementing each other in characterizing the dynamics of quarks and gluons up to the final hadronic states (cf. contributions in [15] for the state of the art). Heavy quarks emerge essentially in early, hard processes, that is, they witness the course of a heavy-ion collision -either as individual entities or subjects of dissociating and regenerating bound states.
Accordingly, the heavy-quark physics addresses such issues as charm (c,c) and bottom (b, b) dynamics related to transport coefficients [14,[16][17][18][19][20][21] in the rapidly evolving and highly anisotropic ambient quark-gluon medium [22,23] as well asQQ states as open quantum systems [24][25][26][27]. The wealth of experimental data from LHC, and also from RHIC, enables a tremendous refinement of our understanding of heavy-quark dynamics. For a recent survey on the quarkonium physics we refer the interested reader to [28].
The yields of various hadron species, light nuclei and anti-nuclei emerging from heavyion collisions at LHC energies are well described by the thermo-statistical hadronization model [29,30] over an interval of nine orders of magnitude. The final hadrons and nuclear clusters are determined by two parameters: the freeze-out temperature T f o ≈ 155 MeV and a freeze-out volume depending on the system size or centrality of the collision. Due to the near-perfect matter-antimatter symmetry at top LHC energies the baryo-chemical potential µ B is exceedingly small, µ B /T f o 1. While the authors of [31] see a delicate interplay of elastic and inelastic hadron reactions as governing principle of the hadro-chemical freezeout, it is argued in [30] that the freeze-out of color-neutral objects happens just in the demarcation region of hadron matter to quark-gluon plasma, i.e. confined vs. deconfined strong-interaction matter. In fact, lattice QCD results report a pseudo-critical temperature of T pc = (156 ± 1.5) MeV [32] and (158.0 ± 0.6) MeV [33] -values agreeing with the disappearance of the chiral condensates and the maximum of some susceptibilities. The key is the adjustment of physical quark masses and the use of 2+1 flavors [34,35], in short QCD 2+1 (phys). Details of the coincidence of deconfinement and chiral symmetry restoration are matter of debate [36]. Reference [37] advocates flavor-dependent freeze-out temperatures.
Note that at T pc no phase transition happens, rather the thermodynamics is characterized by a cross-over accompanied by a pronounced nearby minimum of the sound velocity. This situation continues to non-zero baryon density as long as the baryo-chemical potential µ B is small, µ B /T pc 1.
Among the tools for describing hadrons as composite strong-interaction systems is holography. Anchored in the famous AdS/CFT correspondence, holographic bottom-up approaches have facilitated a successful description of mass spectra, coupling strengths/decay constants etc. of various hadron species. While the direct link to QCD by a holographic QCD-dual or rigorous top-down formulations are still missing, one has to restrict the accessible observables to explore certain frameworks and scenarios. We consider here a framework which merges for the first time (i) QCD 2+1 (phys) thermodynamics described by a dynamical holographic gravity-dilaton-Maxwell background and (ii) holographic probe quarkonia.
We envisage a scenario which embodies QCD thermodynamics of QCD 2+1 (phys) and the emergence of hadron states at T c at the same time. One motivation of our work is the exploration of a holographic model which is in agreement with the above hadron phenomenology in heavy-ion collisions at LHC energies. Early holographic studies [38][39][40] to hadrons at finite temperatures faced the problem of meson melting at temperatures significantly below the deconfinement temperature T pc . Several proposals have been made [41][42][43] to find rescue avenues which accommodate hadrons at and below T pc . Otherwise, a series of holographic models of hadron melting without reference to realistic QCD thermodynamics, e.g. [44][45][46][47][48][49][50][51][52] -mostly with emphasis on quarkonium melting -, finds quarkonia states well above, at and below T pc in agreement with lattice QCD results [53][54][55][56]. It is therefore tempting to account for the proper QCD-related background.
In the temperature region T = O(T pc ), the impact of charm and bottom degrees of freedom on the quark-gluon-hadron thermodynamics is minor [57]. Thus, we consider quarkonia, in particular bottomonium, as test particles. We follow [58][59][60][61] and model the holographic background by a gravity-dilaton set-up supplemented by a Maxwell field [62,63], i.e. without adding further fundamental degrees of freedom to the dilaton. That is, the dilaton potential and its coupling to the Maxwell field are adjusted to QCD 2+1 (phys) lattice data. Our emphasis is here on the formation of bottomonium in a cooling strong-interaction environment.
Thereby, the bottomonium properties are described by a spectral function. The primary aim of the present paper is to study the impact of a finite baryon density of the stronginteraction medium, thus complementing [64,65]. Finite baryon effects become relevant at smaller beam energies, e.g. at RHIC, and are systematically accessible in the beam energy scans [66][67][68]. We restrict ourselves to equilibrium and leave non-equilibrium effects, e.g. [69,70], for future work.
Such µ B > 0 effects on holographic bottomonium spectroscopy have been considered, e.g.

II. BOTTOM-UP MODEL FOR QUARKONIA
In thermal equilibrium, the admixture of equilibrated heavy quarks in strong-interaction matter at T < 250 MeV is small [57,77]. Rather, initial hard parton interactions (essentially gluon fusion) create heavy quarks in heavy-ion collisions. Thus, heavy quark pairs serve as test particles and need not to be back-reacted. In particular, quarkonia constituents are decoupled from the ambient quark content, with the exception of the gluon component.
In a model with minimalistic field content one would prefer to keep the effective gravity-dilaton background (extended by the Maxwell field B for mimicking µ B > 0) to catch QCD thermodynamics and attribute to the test particles solely one vector field A. A U (1) gauge field A(z) in the bulk is supposed to be the dual of the vector meson current operatorQγ µ Q at the boundary. The string-frame action is accordingly where F A stands for the Abelian field-strength tensor of A and k V = Nc 24π 2 with number N c = 3 of colors.
In contrast to common previous practice, the background quantities g 5 (metric determinant) as well as φ (dilaton field) and B (Maxwell field) are universal for any test particle, therefore, G m encodes solely the essential properties of the respective test particle. We attribute the quarkonia masses to the considered test particle. Rather than including the heavy-quark masses explicitly, we encode them in the following manner in G m . From the the ansatz A µ = µ ϕ(z) exp{ip ν x ν } with µ, ν = 0, · · · , 3, which uniformly separates the z dependence of the gauge field by the bulk-to-boundary propagator ϕ for all components of A, and the constant polarization vector µ and gauges A z = 0 and ∂ µ A µ = 0, the equation of motion follows from the action (1) as which is cast in the form of a one-dimensional Schrödinger equation with the tortoise coor- by the transformation ψ(ξ) = ϕ(z(ξ)) exp{ 1 2 ξ 0 dz S(ξ)} and p µ p µ → m 2 n . One has to employ z(ξ) from solving ∂ ξ = (1/f )∂ z . The Schrödinger-equivalent potential in (3) is as a function of ξ(z) with A prime means the derivative w.r.t. the bulk coordinate z.
At T = 0, we have f = 1 and ξ = z, and m n in Eq. (3) is the quarkonium mass spectrum to be used as input. Therefore, the Schrödinger-equivalent potential U (z) must be chosen in such a manner to deliver the wanted values of m n . With given U (z), the Ricatti equation (4) must be solved for S, which in turn determines the heavy-quark mass-specific function G m (φ) via Eq. (5). This G m (φ) is assumed as independent of temperature and baryo-chemical potential, i.e. is ready for direct use at T > 0 and µ B > 0 as well.
In the described chain of operations for getting G m , the zero-temperature background quantities A(z) and φ(z) are needed. They are determined by the temperature independent dilation potential V (φ), which is adjusted to lattice-QCD thermodynamics data, briefly recalled in the next section.

BOTTOM-UP MODEL
We follow here closely the Einstein-dilaton-Maxwell (EdM) model of [75], see also [73,74,76]. The EdM action reads where R is the Einstein-Hilbert part, with warp function A and blackening function f , already used in Section II.
We relegate the field equations following from the action (6) in the coordinates (7) to Appendix A, but mention here the employed parameterizations and refer the interested reader to [75] for listings of the parameters a 1,··· ,10 , φ x , c 0,··· ,5 etc.
Despite the direct relation to an observabe, the location of min{v 2 s } is not so precisely constrained by lattice-QCD data as that of the maximum of chiral susceptibility which determines quite accurately T pc [32,33]. In so far, the curves T min{v 2 s } (µ B ) and T pc (µ B ) need not to coincide. The EdM model with these input data is then ready to transport the thermodynamic information from µ B = 0 to µ B > 0, thus uncovering the T -µ B plane. This is very much the spirit of the quasi-particle model [78,79], where a flow equation facilitates such a transport.

IV. SPECTRAL FUNCTIONS
The equation of motion (2) of ϕ can also be employed to compute quarkonia spectral functions, cf. [38, 47-49, 80, 81]. For ω 2 = p µ p µ > 0 fixed, the asymptotic boundary behavior facilitates two linearly independent solutions by considering the leading-order terms on both sides of the interval [0, z H ]. (i) For z → 0 + , one has the general solution ϕ(z → 0) → A(ω)ϕ 1 + B(ω)ϕ 2 , due to the AdS asymptotic at the boundary, with two ω-dependent complex constants A and B, and horizon, z → z − H , the asymptotic behavior of solutions of (2) is steered by the poles of 1/f and 1/f 2 . The two linearly independent solutions are ϕ where ϕ ± represent out-going and in-falling solutions, respectively. The general near-horizon solution is given by ϕ(z → z H ) → C(ω)ϕ ( z) + +D(ω)ϕ − (z), again with complex constants C and D which depend on ω. The side conditions for the bulk-to-boundary propagator are The corresponding retarded Green function G R of the dual current operatorQγ µ Q, defined within the framework of the holographic dictionary via a generating functional by with A 0 µ ≡ µ exp{ip ν x ν } for µ ∈ {1, 2, 3} [81]. The quantity S V, on-shell m denotes here the action (1) with the solution ϕ from (2). Finally, the spectral function ρ follows from It has the dimension of energy squared, suggesting to use L 2 ρ or ρ/ω 2 as convenient representations.

V. NUMERICAL RESULTS
The spectral function ρ(ω, T, µ B ) is accessible by numerical means by the following chain of operations: (i) solving the equations of motion (A1 -A4) following from the action (6) with boundary conditions (A5 -A10) for the background encoded in with the prescribed V (φ 0 ) from Eq. Our results are exhibited in Fig. 1 for the Υ meson. We emphasize that neither an explicit quark-mass dependence enters our approach (instead, quark masses are implicitly accounted via U 0 (z) for entering G m ) nor a confinement criterion (instead, narrow spectral functions as quasi-particle states are considered as confined J P C = 1 −−b b states). In so far, the emergence of such narrow quasi-particle ground states at T = O(T pc ) is astonishing.
The higher the excitation, the later the excited-state formation happens when considering Focusing on the crucial temperature region near T pc or T v 2 s , one observes how rapidly the ground state evolves towards a sharp quasi-particle within this narrow interval of T An adaptiveω grid with minimum spacing of 10 −8 L is employed. Note that feeding is not included.
emphasize that, at µ B = 0, the first and (weakly) the second excitations are identifiable as peak structures, in contrast to [52], where these excitations appear as molten, while the ground state persists up to high temperatures since it is kept by a narrow deep well potential.
To highlight the µ B dependence, we exhibit in Fig. 3  to follow the evolution of the spectral function. Figure 3 provides some guidance for that.
The relation of the spectral function to the resulting µ + µ − spectrum from Υ → µ + µ − may be elaborated as in previous studies, e.g. by superimposing the thermal yield (which needs a model of the space-time evolution of the fireball) and the post-freeze-out contribution (which is directly related to the Υ(nS) yields and feedings) and the various background sources. This is beyond the scope of our paper. Nevertheless, the emerging picture of our model (see increasing with (or independent of) radial quantum number n, in contrast to experimental data. This imperfection seems to continue to T > 0: the strengths of (i.e. yields from) excitations become too large. Therefore, we do not introduce a continuum or background subtraction, as discussed in [52], and leave further refinements (e.g. the options offered in Appendix B 3 w.r.t. decay constants) and feeding corrections to follow-up work.

VI. CONCLUSION AND SUMMARY
Using a bottom-up holographic model with minimalistic field content, we investigate the impact of the finite baryo-chemical potential µ B on bottomonium formation at temperatures in the order of the hadron-chemical freeze-out in relativistic heavy-ion collisions. The model has two pillars: (i) the vector meson part which employs the bottomonium masses of ground state and radial excitations as input to adjust a suitable Schrödinger-equivalent potential, and (ii) the Einstein-dilaton-Maxwell background which is adjusted independently to lattice-QCD thermodynamics (sound velocity and light-quark susceptibilities). The field content is as follows: (i) a bottomonium-specific function G m (φ), which encodes implicitly the b quark masses via the Schrödinger-equivalent potential U 0 (z) and is essential for the bulk-to-boundary propagator ϕ(z), and (ii) the gravity-dilaton-Maxwell part, determined dynamically by the dilaton potential V (φ) and the dilaton-Maxwell coupling F(φ). Since there is neither a confinement criterion nor a chiral condensate as order parameter in such an approach, we consider the shrinking of the spectral function ρ (determined by ϕ) to a narrow quasi-particle state as bottomonium formation in a cooling strong-interaction medium. Despite a simple two-parameter Schrödinger-equivalent potential U 0 , we find the bottomonium ground-state formation at about a temperature of 150 MeV at µ B = 0. Increasing µ B drops the formation temperature. Excited states are consecutively formed at lower temperatures, or the spectral strengths are not yet concentrated completely at the T = 0 quasi-particle energy at a given temperature. This fits well in the experimental observation that the Υ(2S) and Υ(3S) states are hardly identifiable in the di-lepton spectra in heavy-ion collisions ate LHC, while the ground state is clearly visible [1]. In contrast, Υ(1S, 2S, 3S) are clearly seen in proton-proton collisions at the same beam beam energies per nucleon.
Our approach assumes rapid thermalization and equilibration, since the cooling of the medium is handled as sequence of equilibrium states. Off-equilibrium phenomena up to dynamical freeze-out need to be considered in refined investigations. Highly desirable would be a closer contact to string theory to overcome the deployed phenomenological parameterizations steering our two pillars, background thermodynamics and vacuum mass spectrum.
Appendix A: Using the bulk coordinate z in the EdM model The equations of motion for the dilaton φ(z), the Maxwell field Φ(z) with the coupling function F(φ), following from the action (6), read in the coordinates (7) with warp factor exp{A(z)} and blackening function f (z) The leading-order initial conditions are (i) near boundary, i.e. z → 0 + , References [73][74][75][76] use essentially the coordinates originally employed in [62,63]  of U 0 (z). Here, we suppose that Eq. (B1) delivers a set of discrete eigenvaluesÊ n ≡ L 2 m 2 n , n = 0, 1, 2, · · · by the requirement of square-integrable solutions y n .
We emphasize again hat U 0 (z) is an independent input in our approach which determines  2) and (10) the spectral function. In so far, the choice of U 0 deserves some special attention.

Approximately uncovering the Υ mass spectrum
The famous soft-wall (SW) model [83] employs U SW 0 = U U V 0 +U IR 0 = 3 4 z −2 +(a/L) 2 (z/L) 2 with the leading-order asymptotic parts It has one free parameter, a, and, in general, cannot accommodate independently ground state mass and level spacing at the same time. Nevertheless, it delivers via L 2 m 2 n = 4a(n+1), n = 0, 1, 2 · · · , the Regge type mass spectrum -in [83] termed "linear confinement". Despite the imperfection, it has been used in [47,48] for an investigation of the termal behavior of the , step (i) would fix a, and b is obtained in step (ii). While such a two-step fine-tuning procedure looks promising, it could be hampered by a problem which we faced, e.g. in [64,65]: unfavorable parameterizations of U 0 (z) can lead to too low formation temperatures, such that at T pc quasi-particles are not yet formed, in contrast to the common understanding of hadron formation in relativistic heavy-ion collisions discussed in the introduction. The origin of the affair can be qualitatively explained within the transparent model U SW,b 0 . At finite temperatures, U 0 (z) → U T (z, z H (T )). Since To accommodate the ground state in such a potential, the IR turning point (t.p., in the spirit of WKB) z IR t.p. ≈ 2 √ a + b/a must obey z IR t.p. < z H ≈ 1/(4πT ). In other words, to allow for an "unmolested" state at given temperature T , the parameter a must be sufficiently large to get small z IR t.p. . This is the reasoning of considering the quark-mass effect encoded in m 0 as primary quantity and a less strict parameter adjustment for the level spacing, as deployed in the above parameter setting. (In stark contrast, [84] puts emphasis on the correct decay constants and is less restrictive to the bottomonium mass spectrum with the advantage of rather persistent states up to high temperatures T > T pc .) with constant parameters b U V,IR ,Ũ 0 , λ and scale setting parameter L. The options λ = 1 (dip related to discontinuity at x 0 ), λ → 1 (mimicking a Dirac delta dip at x 0 whenŨ 0 ∝ 1/(λ − 1)) and λ > 1 (box-like dip within x 0 · · · λx 0 ) w.r.t. fine-tuning of the mass spectrum are discussed in the Supplemental Material. We finish this essay by the expectation that the tendency of the µ B dependence of spectral functions is not obstructed by the details of approximately or accurately adjusting parameters of U 0 to the Υ mass spectrum.
a particular point in parameter space, we are interested in the systematic enabled by the ansatz (B4).
The solutions of the respective Schrödinger equations are as follows: where J 1 and Y 1 stand for the Bessel functions, and D ν are parabolic cylinder functions; x ≡ z/L from here on. 2 Given the asymptotic behavior (i) lim x→0 Y 1 (x) → −2/(π x) + · · · , the square-integrability requirement in x ∈ (0, x 0 ) facilitates c U V 2 = 0, and (ii) since with (i) W = √ 1 + 4α 2 for a = 0 and α 2 = 3/4, yielding (I), as well as (ii) α 2 = 0, yielding (III) (see [86] for the definition and properties of the confluent hypergeometric (Kummer) function of second kind U (a, b, x), also known as Tricomi function, and the associated Laguerre polynomial L ρ γ (x)). Square integrability for x ∈ (0, ∞) unravels the energy eigenvaluesÊ = a(4n+2 α 2 + 1 4 +2) which become these of the soft-wall model for α 2 = 3/4, i.e.Ê = 4a(n + 1). Adding a Dirac delta at x δ with strength −∆, the solutions (*) for x < x δ and x > x δ must be matched at x δ with the above condition for the discontinuity of y : The continuity condition of y at x δ determines one of the four integration constants of the general solution (*) with its two branches and the discontinuity condition of y fixes another. The remaining two constants are determined by the claim of square integrability and the related appropriate boundary conditions and the normalization of the solution in case of normalizability for all energy eigenvaluesÊ. In contrast to highly symmetric square-well/harmonic oscillator + Dirac delta models of [87,88], the ground state energy cannot be dialed independently of the excited states and their (non-uniform) level spacing. 3 Since for ν ≡Ê −a 2a = n, n = 0, 1, 2 · · · , the parabolic cylinder functions are related to Hermite polynomials H n via D n (x) = 2 −n/2 exp{−x 2 /4}H n (x/ √ 2), and Eq. (III) becomes is for an altitude-limited r.h.s. hard wall yieldingÊ ≈ (x /x 0 ) 2 with x being the zeroes of J 1 , J 1 (x ) = 0, = 1, 2, · · · , see dashed curves. In both cases, the respective energiesÊ n, must be smaller than the maxima of the altitude-limited hard walls:Ê n < L 2 U U V 0 (x − 0 ) or E < L 2 U IR 0 (x + 0 ).
Despite the noticeable change of the energy eigenvaluesÊ n with changing parameter x 0 (see solid and dotted curves in the left panel), the level spacing is less influenced, see solid and dotted curves in the middle panel. The level spacing, given be differences of energy eigenvaluesÊ n −Ê 0 = L 2 (m 2 n − m 2 0 ) (≈ 4a(n + 1) for x 0 < 1, see middle panel), can be related to the Υ(1S, 2S, 3S) mass-squared differences m 2 n − m 2 0 , 10.815 GeV 2 (n = 1) and 17.623 GeV 2 (n = 2), pointing either to a ≈ 0.1 (n = 1) or a ≈ 0.08 (n = 2), since L −1 = 5.148 GeV. Clearly, the choice λ = 1 only approximately accommodates the proper level spacing, as evidenced by two slightly different values of a.
In fact, the ratioÊ 2 −Ê 0 E 1 −Ê 0 (which is independent of L 2 , a and a shift by b) quantifies better the level spacing. It can be directly related to the experimental ratio The case λ → 1 means replacing theŨ 0 dip in Eq. (B4) for x ∈ [x 0 , λx 0 ] by a Dirac delta, −∆δ(x − x 0 ). As pointed out above and in footnote 2, the solutions y(Eq. (I)) and y(Eq. (III)) and their derivatives must be matched at x 0 to yield for discrete energy eigenvaluesÊ(x 0 , a, ∆) to be numerated by n, see dotted curves in Fig. 5.