Estimate for the bulk viscosity of strongly coupled quark matter using perturbative QCD and holography

Modern hydrodynamic simulations of core-collapse supernovae and neutron-star mergers require knowledge not only of the equilibrium properties of strongly interacting matter, but also of the system’s response to perturbations, encoded in various transport coefficients. Using perturbative and holographic tools, we derive here an improved weak-coupling and a new strong-coupling result for the most important transport coefficient of unpaired quark matter, its bulk viscosity. These results are combined in a simple analytic pocket formula for the quantity that is rooted in perturbative Quantum Chromodynamics at high densities but takes into account nonperturbative holographic input at neutron-star densities, where the system is strongly coupled. This expression can be used in the modeling of unpaired quark matter at astrophysically relevant temperatures and densities.


INTRODUCTION
During the last ten years, neutron stars (NSs) and their binary mergers -observable through both electromagnetic and gravitational waves (GW) [1, 2] -have established themselves as the leading laboratory for dense Quantum Chromodynamics (QCD) matter.While the observable properties of single quiescent NSs and even the inspiral parts of NS mergers are mostly determined by the equation of state (EoS) of the constituent matter, the ringdown phase of a NS merger constitutes a considerably more complicated out-of-equilibrium system.In preparation for the eventual observation of a ringdown GW signal, extensive hydrodynamic simulations of NS mergers are currently being carried out, with one crucial challenge being to correctly account for energy dissipation and transport in NS matter [3].
Among the different transport coefficients, the bulk viscosity ζ, which quantifies energy dissipation during a rapid compression or expansion of matter, stands out as particularly important [4][5][6][7][8][9][10][11][12].For isolated NSs, it affects the emission of continuous GWs [13], expected to be detectable in next-generation GW observatories such as the Einstein Telescope [14] and Cosmic Explorer [15], and determines the maximal rotation frequencies of pulsars in a temperature-dependent fashion, giving rise to the so-called r-mode stability window in the 1-100 keV range [16][17][18] (for a review of NS oscillatory modes, see [19]).In NS mergers, the bulk viscosity on the other hand provides damping for density oscillations, affecting both the inspiral [20] and post-merger dynamics, of which the latter involves temperatures up to tens of MeVs.The bulk viscosity may indeed leave a detectable imprint on the post-merger GW waveform [21][22][23][24][25][26], the magnitude of which is however still under discussion [27].
The dominant contribution to the bulk viscosity comes about when weak interactions cannot keep pace with the compression rate, leading to deviations from beta equilibrium and a non-equilibrium contribution to the pressure, against which work can be done.This effect peaks when the timescales of macroscopic oscillations and microscopic flavor-changing rates match.In the nuclear matter phase, the value of ζ depends on multiple factors, such as whether direct Urca processes are allowed or if hyperons or Cooper pairing between nucleons are present, each affecting in particular the temperature scale where ζ reaches its maximal value (see ref. [28] for a review).
The first milliseconds of a binary NS merger are known to involve baryon densities up to several nuclear saturation densities n sat ≈ 0.16/fm 3 as well as temperatures up to several tens of MeV's (see, e.g., [31]).Such conditions may also lead to the creation of deconfined QM [32][33][34][35][36], F r e e q u a r k s , 5 n s a t FIG. 1.The bulk viscosity ζ of NS matter, evaluated at rotation frequency ω = 2π × 1 kHz and given as a function of T for a baryon density nB ≈ 5nsat.The uncertainty bands of the holographic results are assessed via their matching to QCD: The D3-D7 result is matched to pQCD quark densities including their uncertainty bands, while the uncertainty of the V-QCD result is estimated by varying the parameters of the model within limits set by lattice-QCD results.Finally, nuclear and hyperonic matter results (labeled "Nucl."and "Hyperons") from Refs.[24,29] are shown for comparison.Note that for technical reasons, the V-QCD result is shown for 7nsat and the hyperonic one for 4.5nsat.We observe that our QM results always peak within the r-mode stability window 1-100 keV, but are strongly suppressed at the O(10 MeV) temperatures involved in NS mergers.This may, however, be related to the absence of quark pairing in our setup (see [30] for a counterexample in the Color-Flavor-Locked case).
the transport properties of which differ significantly from those of nuclear matter.While the value of the QM bulk viscosity is expected to strongly depend on the presence and details of quark pairing, differences between various partially paired configurations are expected to be smaller than between quark and nuclear matter [28].This makes the bulk viscosity an interesting quantity for tracking the possible creation of QM during mergers.Despite the phenomenological importance of the bulk viscosity, our ability to predict its behavior remains limited owing to the unavailability of controlled firstprinciples quantum-field-theory methods at NS densities.The leading first-principles tools include perturbative QCD (pQCD), available only at very high densities (see, e.g., [37][38][39]), and holography, which describes the strong-coupling limit of a class of QCD-like theories [40][41][42][43][44].For QM, leading-order perturbative results for several transport coefficients were derived some thirty years ago [45,46] and improved to next-to-leading order (NLO) later [47,48], whereas at strong coupling, the shear viscosity and the electrical and thermal conductivities were first evaluated only recently in two holographic models [49,50].For the bulk viscosity, only the minuscule purely QCD contribution has been considered in recent literature [49,51], but for the dominant contribution stemming from an interplay between the electroweak and strong sectors, no strong-coupling prediction is currently available at all.
In this work, we derive state-of-the-art results for the thermodynamic response of QM to a change in its flavor content, thus providing novel predictions for the bulk viscosity.We do so using both perturbative and holographic methods, and in particular derive the first strongcoupling predictions for the quantity.Our results are applicable for unpaired QM and serve as a starting point for any partially unpaired phase [28].
The main result of our work is shown in fig. 1, where we display the bulk viscosity of NS matter as a function of temperature for a baryon density of roughly 5n sat .For QM, we include results corresponding to the freetheory limit, evaluated at a fixed strange quark mass (m s = 93.4MeV), as well as our two holographic models, D3-D7 and V-QCD, but not the pQCD result, which is not under quantitative control at intermediate densities.
For the confined phase, we display results corresponding to both nuclear [24] and hyperonic [29] matter.As we discuss in detail in the remaining sections of this letter, our results paint a consistent picture of the behavior of the QM bulk viscosity that displays a stark qualitative difference to that witnessed in the confined phases of QCD.Furthermore, we observe that for astrophysically relevant densities and temperatures, nearly all temperature dependence in the QM result originates from the flavor-changing interactions.For our D3-D7 computation, this leads to a simple analytic result for ζ, given in eq. ( 5) below, that we suggest for use as an approximation for the bulk viscosity of unpaired QM in future phenomenological applications.

SETUP
For unpaired three-flavor QM in the neutrinotransparent regime, the leading contribution to the bulk viscosity arises from W -boson exchange in the process u + d ←→ u + s.Outside beta equilibrium, i.e., when the d and s quark chemical potentials differ µ d ̸ = µ s , the quark densities n d and n s change with rates proportional to an electroweak rate λ 1 [52][53][54], so that Neglecting quark masses, the leading low-T contribution to the rate becomes [54] 1 where G F is the Fermi constant and θ c the Cabibbo angle.The quartic prefactor on the right-hand side represents the only known O(α s ) correction to the rate, which is moreover logarithmically enhanced at low temperatures as it originates from a so-called non-Fermi-liquid (nFL) contribution to the specific heat of QM [56] (see also [57,58]).As discussed in detail around fig. 7 of the Supplemental Material, this correction allows us to gauge the importance of the (partially unknown) O(α s ) corrections to the rate: for σ = 0, the result reduces to the leading-order rate, while for σ ≡ 4α s /(9π) and Λ ≈ 0.158 s one recovers the result derived in [56].
While the unknown QCD corrections to the rate may be sizable, we note that the qualitative behavior of the rate likely remains the same at strong coupling: In holography, the QCD contribution to the rate, replacing the leading-order multiplicative factor µ 5 d T2 above, is available from the convolution of two flavor-current correlators.For these correlators, calculations at non-zero quark densities in the D3-D7 model show a linear dependence on the temperature at low frequencies [59][60][61], consistent with the formula we use.Furthermore, the normalization of the correlators depends on the number of colors and flavors but not on the 't Hooft coupling, thus keeping the rate constant in the strong-coupling limit.
A study of energy dissipation during a compressiondecompression cycle near beta equilibrium connects ζ to various susceptibilities χ ij ≡ ∂ 2 p/∂µ i ∂µ j and reaction rates (see Supplemental Material sec.A). 2 If we only take into account the u + d ←→ u + s process [48], this leads to where the coefficients A 1 and C 1 , determined by various susceptibilities and quark densities, are found in eqs.(38) and (39) of the Supplemental Material and ω denotes the angular frequency of density oscillations (see [19] for discussion).3 The combination of susceptibilities appearing in eq. ( 38) vanishes if the d and s quarks are degenerate in mass -a fact most easily verified if (38) is given in terms of the inverse susceptibility matrix (see Supplemental Material for details).This implies that a nonzero strange quark mass must be implemented in both the weak-and strong-coupling setups, which we briefly introduce below.

METHODS
In this section, we review our perturbative and holographic determinations of the susceptibilities that enter eq.(3).In both calculations, we treat electrons as non-interacting and (numerically) solve the corresponding chemical potential µ e from the charge neutrality condition 2n u /3 − n d /3 − n s /3 = n e = T 2 µ e /3 − µ 3 e /(3π 2 ).Together with the beta-equilibrium conditions µ s = µ d , µ u = µ d − µ e , this allows us to obtain ζ in terms of µ d , T , ω.Finally, our results will depend on the parameter X ≡ Λ/(2µ d ) which parametrizes our results' dependence on the unphysical renormalization scale Λ in the MS scheme.It appears directly in our pQCD results and indirectly in the D3-D7 ones, where it enters through the high-density matching of the model to pQCD.

Perturbative QCD
For vanishing quark masses, the perturbative pressure of deconfined unpaired QCD matter is known up to order α 5/2 s at nonzero temperatures and densities [62,63] and up to partial O(α 3 s ) in the T = 0 limit [39,64,65].Up to the highest fully known order O(α 5/2 s ), the result can be split into two distinct terms corresponding to contributions from the hard and soft momentum scales, which for µ ≫ T are of order µ and α 1/2 s µ, respectively.We treat the additional mass-dependent contribution to the pressure p m within the mass-expansion scheme of [66], where m s is formally treated as a quantity of O(α 1/2 s µ) and the light quark masses are neglected.This mass expansion is performed to O(m4 s ) and up to a combined O(α 5/2 s ) (for the full mass-dependence at T = 0, see [37]).For the value of the s quark mass, we use the physical MS renormalized value m s ≈ 93.4 MeV [67].We have confirmed that additionally including nonzero m u and m d terms would lead to a vanishingly small effect, while the chemical potentials realized in NSs are not large enough to allow for heavier quarks.For the soft contribution, evaluated in the massless limit, we furthermore use an analytic small-T /µ expansion derived in [63] that is valid for T ≲ 100 MeV.Mass corrections to this result start at O(α 3 s ) and can therefore be neglected.The perturbative pressure described above can be readily differentiated to obtain predictions for the coefficients A 1 , C 1 and eventually for ζ as functions of the three quark chemical potentials and the renormalization scale parameter X.The results constitute lengthy closedform expressions in terms of standard special functions and their derivatives, allowing for inexpensive evaluation of the necessary quantities.

Holography
The D3-D7 model [68] is the holographic dual of N = 4 SU(N c ) super Yang-Mills theory with N f copies of N = 2 hypermultiplets in the quenched approximation N f /N c ≪ 1.It consists of N f probe D7-branes embedded in the AdS 5 × S 5 spacetime, while baryon charge is introduced by turning on an electric field on the D7-branes [69,70] and temperature by modifying the geometry to that of a black brane.Following [71], we extrapolate the model to the physically relevant N c = N f = 3 and fix α s ≈ 0.285 so that the pressure matches the Stefan-Boltzmann value at high density, extending the model's validity towards higher densities.Although the field content of the model differs from that of QCD, we note that the thermodynamic coefficients A 1 and C 1 , obtained through chemicalpotential derivatives of the pressure, are highly insensitive the additional fields in the D3-D7 model.
At vanishing temperature, the pressure of the D3-D7 model takes the simple form [71,72] where M i are the constituent quark masses that we fix by equating quark densities with pQCD at µ d = 1 GeV and varying X ∈ [1/2, 2].Doing so, we obtain M u ∈ (522.5, 434.6)MeV, M d ∈ (526.4,435.9)MeV, and M s ∈ (541.8,450.1)MeV, within this interval in X.In what follows, in addition to estimating uncertainties by matching to pQCD at different values of X, we also vary this matching density within we finally compute the pressure numerically, following methods introduced in [69,70].
The other holographic model we use is V-QCD [73], which is a bottom-up model tuned to reproduce QCD physics as closely as possible (see, e.g., the reviews [42,43,74]).It combines the improved holographic QCD model for pure Yang-Mills theory [75,76] to a description of flavors introduced via tachyonic brane actions [77][78][79], featuring, e.g., a running α s as reviewed in the Supplemental Material.Given that quarks are treated as unquenched (N f /N c ∼ 1) in V-QCD, the model should capture their physics more realistically than the D3-D7 model.Indeed, V-QCD by construction agrees with various qualitative properties of QCD (such as confinement and asymptotic freedom), and its parameters are fitted to data, including lattice results for the pressure [80,81] and baryon number susceptibilities [81] at µ = 0.The model is consistent with all known astrophysical observations in the NS-matter regime [82,83], but eventually becomes inconsistent with pQCD at high densities [84].
In this paper, we otherwise follow the treatment of the above V-QCD papers but relax the assumption of exact chiral symmetry in the QM phase by turning on a nonzero strange quark mass, thus extending the prescription of [85].The corresponding mass parameter of the model is fixed by demanding that the masses of kaons and η mesons are well reproduced in the vacuum (see Supplemental Material and refs.[86][87][88][89] for details).We find that this procedure under-predicts the dependencies of quark number susceptibilities on the strange quark mass at zero µ and high T , where the results can be benchmarked against lattice data [90].This leads us to expect that this model similarly under-predicts the effects of the strange quark mass in physical quantities at high densities.
Finally, we quantify the underlying uncertainty of our results by allowing the V-QCD parameters vary within limits set by the lattice QCD fit in the chirally symmetric phase [81,91], but otherwise follow the computational strategy of [92] in determining the quantities appearing in eq.(3).In both holographic setups, the variation procedure we perform thus corresponds to the uncertainties associated with the respective matching procedures.

RESULTS
Our main result for the bulk viscosity of unpaired QM is displayed in fig. 1.It highlights a qualitative contrast between the behavior of ζ in the confined and deconfined phases of QCD, with the more suppressed QM results peaking at lower temperatures, and in addition demonstrates the important effect of interaction corrections in the latter case.Consistently with our expectations for quantities that vanish in the degenerate-mass limit, V-QCD appears to predict somewhat lower values for ζ than our other methods, but nevertheless retains the same qualitative features.
A closer inspection of our results reveals a number of interesting further findings.Explicit calculations show that in all three approaches, the bulk viscosity is insensitive to the T -dependence originating from the coefficients A 1 and C 1 of eq.(3).As demonstrated in fig.6 of Supplemental Material, to a good accuracy we can indeed set T = 0 in these functions and only keep the Tdependence of the electroweak rate λ 1 in eq.(2).Another universal characteristic that all our results exhibit is an approximate quartic dependence on the strange quark mass which has been noted before in [55].
While the full ζ depends on the rate λ 1 , we may construct physical features of the bulk viscosity that are  38) and ( 39) and are independent of the oscillation frequency ω and the electroweak rate λ1, as indicated by the expressions on the right vertical axes.The error bars in these panels capture the variation of model parameters in the different models as described in the main text.
sensitive only to QCD input.For example, the peak value of the viscosity, ζ peak ≡ ζ(T peak ), and its rescaled zero-frequency limit λ 1 ζ(ω = 0) that corresponds to the DC bulk viscosity entering the Israel-Stewart theory [12,93,94] are completely insensitive to the electroweak rate and can be fully extracted from the coefficients A 1 and C 1 in eqs.( 38)- (39).These two quantities are shown in fig.2, where we observe a good agreement between our pQCD and D3-D7 results for densities where both predictions are available, while V-QCD again appears to underestimate the quantities (see discussion in Supplemental Material).Setting T = 0 in A 1 and C 1 , we find that the D3-D7 calculation leads to a remarkably simple analytic formula as a function of where we have defined We stress that for the M i in this formula, one should use the constituent quark-mass ranges listed below eq.( 4), leading to the uncertainty ranges visible in fig. 1.To express this as function of n B for the small temperatures of relevance to BNS mergers, one can further use the T = 0 pressure in eq. ( 4) to numerically relate n B to µ d in beta equilibrium.
Returning finally to the bulk viscosity itself, we note that it is straightforward to compare our NNLO pQCD results to lower perturbative orders, as shown in fig. 4 of the Supplemental Material.We find that the difference between the NLO and NNLO results is non-negligible even at 40n sat , and that the results diverge rapidly at lower densities, making extrapolation to the NS realm impossible.While the naive free quark expression can, in principle, be extrapolated to low densities, it completely fails to take into account the effects of interactions, which become increasingly important much before the hadronic phase is eventually reached.For phenomenological purposes, the compact D3-D7 bulk viscosity of eq. ( 5) is on the other hand appealing as it is rooted in pQCD but takes into account the strongly coupled nature of the theory at low densities.To this end, despite its limitations discussed above, we recommend the use of this result in the modeling of dense unpaired QM at astrophysically relevant densities and temperatures, and similarly expect the the V-QCD result to provide a reasonable lower bound for the bulk viscosity.
An important limitation of our present approach is finally related to the fact that the pairing channel and the magnitude of the superconducting gap in low-and moderate-density QM remains unknown (though see [95] for a recent model-independent study bounding the gap at high densities).To obtain estimates for the bulk viscosity in various pairing channels, corrections to both the electroweak rate in eq.(2) and to the thermodynamic functions entering through eqs.( 38)-( 39) should be separately considered.While the latter are expected to be subleading, the former may be substantial given that the contribution of gapped quark modes to the reaction rate is exponentially suppressed.While the detailed evaluation of these corrections is left for future work, we note that the electroweak rate receives O(α s ) QCD corrections even in the unpaired phase, some of which are presently known [56].Their effect is studied in fig.7 of the Supplemental Material, where we observe that, in agreement with the λ 1 -independence of ζ peak , they primarily simply shift the peak of the viscosity to lower temperatures.dation of Korea (NRF) funded by the Korean government (MSIT) (grant number 2021R1A2C1010834).T.G. has been supported in part by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) project-ID 279384907-SFB 1245, by the State of Hesse within the Research Cluster ELEMENTS (pro-jectID 500/10.006),and by the ERC Advanced Grant "JETSET: Launching, propagation and emission of relativistic jets from binary mergers and across mass scales" (Grant No. 884631).C.H. is partially supported by the AEI and the MCIU through the Spanish grant PID2021-123021NB-I00 and by FICYT through the Asturian grant SV-PA-21-AYUD/2021/52177.N.J., R.P., and A.V have been supported in part by the Research Council of Finland grants no.322507, 345070, 347499, 353772, and 354533 as well as by the ERC Consolidator grant no.725369.S.S. acknowledges support of the DFG cluster of excellence ORIGINS funded by the DFG under Germany's Excellence Strategy -EXC-2094-390783311. * jesus.cruz@correo.nucleares.unam.mx† gorda@itp.uni-frankfurt.de‡ hoyoscarlos@uniovi.es§ niko.jokela@helsinki.fi¶ matti.jarvinen@apctp.org∥ aleksi.kurkela@uis.no* * risto.paatelainen@helsinki.fi† † saga.saeppi@tum.de‡ ‡ aleksi.vuorinen@helsinki.fi

Supplemental Material
In the three sections of the Supplemental Material, we go through several details of our calculations and results that provide additional context to the main text.These include a compact derivation of the main formulas for the bulk viscosity, eqs.(3)- (39), that we include here for completeness; details of how we have supplemented the V-QCD setup with a nonzero s quark mass parameter; and a quantitative analysis and comparison of our three independent results for ζ.

A: In-equilibrium processes and bulk viscosity
We review first the processes taking place in near-equilibrium unpaired quark matter and provide a relatively detailed derivation of the bulk viscosity formulas used in the main text.It follows the arguments introduced in [4, 6] and reviewed in [28], but presents more details and extends their analysis to include non-diagonal flavor susceptibilities.
A.1: Beta equilibrium, charge neutrality, and production rates The relevant flavor-changing electroweak processes produced by the exchange of a W -boson are which lead to the chemical equilibration conditions (we assume no neutrino trapping, µ νe = 0) In addition, local charge neutrality implies the following relation for the densities When the system is taken out of equilibrium, the densities of various particle species change according to [55] δn u = δn e , δn d + δn s = −δn e (9) so that the baryon number and electric charge densities do not change, For small deviations from beta equilibrium, the rates of change of the densities through the electroweak processes are proportional to where λ i stand for the electroweak reaction rates for processes in eq. ( 6) enumerated left-to-right.Since the baryon number and the electric charge are not affected by the above processes, we may use the baryon or quark density and the following partial fractions as thermodynamic variables: In this case, the reaction rates to consider are A.2: Formula for the bulk viscosity Changes in volume due to a radial pulsation of the system also change the baryon density.Considering the simple case of a small homogeneous oscillation of period τ = 2π/ω, we can parameterize the deviation of the quark density n q from its equilibrium value n 0 q as ∆n q , fulfilling n q ≈ n 0 q + ∆n q sin(ωt) , which can alternatively be seen as a change in the specific volume V q = 1/n q , such that dV q dt ≈ −ω ∆n q (n 0 q ) 2 cos(ωt) . ( The oscillatory radial pulsation leads to the dissipation of energy Ė due to a nonzero bulk viscosity of the medium, ζ.Averaging over an oscillation period, we obtain where v is the velocity of the fluid due to the pulsation.Using the continuity equation of baryon charge, we can further write the divergence of the velocity in the form so that the dissipated energy becomes At the same time, the dissipated energy can be seen to result from the irreversibility of the compression-decompression cycle, for which the work performed in an infinitesimal change of the specific volume equals dW q = p dV q with p the pressure.Averaging over a period, we again obtain where a comparison with eq. ( 15) shows that nonzero contributions come from terms in p(t) proportional to cos(ωt).Let us next inspect the pressure, which is a function of chemical potentials µ a in the grand canonical ensemble.A shift of these variables translates into where the partial derivatives are taken while leaving other thermodynamic variables fixed, namely the temperature and other chemical potentials.Denoting the matrix of susceptibilities as it is then straightforward to derive the following relation between the rates from (12): Using ( 12), the pressure becomes in turn Next, we use known relations between the partial fractions of ( 12) to convert eq. ( 13) into a set of differential equations for the chemical potentials.Defining the chemical potentials describing the deviation from the beta equilibrium we obtain (α, β = 1, 2; a = u, d, s, e) where the coefficients read M a s,e = (χ −1 ) a s,e , and In terms of µ 1 and µ 2 , the change in the pressure (23) then becomes At this point, we can reduce the problem to an algebraic system by expanding the chemical potentials in sine and cosine µ e = µ 0 e + ω∆n q [c e cos(ωt) + s e sin(ωt)] , which lead to the set of equations (Γ = 1, 2, s, e and β = 1, 2) Introducing the pressure ( 28) in (19) with the solutions we have written, we obtain now and comparing with ( 16), we find a simple expression for the bulk viscosity, Defining finally a matrix with components the solution for the cosine coefficients in (30) becomes A.3: Non-leptonic limit To conclude the derivation of the bulk viscosity formulas, let us finally briefly study how the result simplifies is the leptonic processes are neglected by setting λ 2 = λ 3 = 0 and only the non-leptonic processes u + d ←→ u + s retained.Using the symmetry of the susceptibility matrix, we now define the coefficients whereby the matrix M in (33) simplifies to Solving for the coefficients in (34) and introducing the result in the bulk viscosity formula (32), we obtain i.e., eq. ( 3) from the main text.Further inverting the susceptibility matrix, we see that the coefficients A 1 and C 1 acquire the forms where H p = det χ ij is the Hessian determinant of the pressure p.These results include the effects of off-diagonal susceptibilities, which to the best of our knowledge have not been included in the determination of ζ for unpaired quark matter before.If we finally assume that the susceptibilities are diagonal, χ b a = χ aa δ b a , the coefficients take the simple form in which case re-inverting the susceptibility matrix becomes trivial.

B: Details on the V-QCD calculation
Given the sensitivity of the bulk viscosity on the quark masses, we need to generalize the V-QCD model in a way that includes quark flavors with non-degenerate masses.While the flavor dependence is implicitly included in the model as defined in earlier literature [73,92], the implications of flavor dependence have not been properly studied before.Here, we briefly review the computation of the QM equation of state and Nambu-Goldstone boson masses in the flavored case.We stress, however, that the flavor-dependent V-QCD extends the flavor-independent case in a rather simple fashion: to the most part, one just needs to sum over the contributions from each flavor.
The fields of the holographic dictionary relevant for the present setup are the following: • The dilaton ϕ, which is a scalar field dual to F 2 , where F is the field strength operator of the gluon field.The source corresponding to this operator is the gauge coupling of QCD.
• The tachyon T ij , where i, j = 1 . . .N f are flavor indices, which is dual to the operator ψi ψ j with ψ i the i-quark field in QCD.The source corresponding to this operator is the quark mass matrix.
• The (vectorial) bulk gauge fields A ij µ , which are dual to the currents ψi γ µ ψ j .The sources for the temporal components include the chemical potentials of different quark flavors.
The relevant terms in the model action are on the other hand given by where S g is the action of the improved holographic QCD (IHQCD) model [75,76].It describes the gluon sector of the theory via five-dimensional Einstein-dilaton gravity where the Planck mass M p and the potential V g need to be determined by comparing to QCD data [80].The flavor action on the other hand takes the form of a generalized Dirac-Born-Infeld action [77,78] where we assume that all flavor-dependent fields are diagonal: Here, the capital Latin indices run over all five dimensions whereas the Greek indices denote the Lorentz indices in four dimensions.
The potentials V g , V f , κ, and w in the above actions can be determined by comparing to QCD data [80,81,86,87,91].In this article, we use the potential sets 5b, 7a, and 8b fitted to lattice thermodynamics in [81,91].We have also checked that the potentials constructed in [87], which additionally use data for meson masses and decay constants, produce similar results.
We use the metric ansatz with the AdS 5 boundary located at r = 0.As we are interested in the equilibrium thermodynamics of homogeneous phases, the metric as well as the background values for the gauge fields and scalars are assumed to be functions of r only.We only consider black-hole backgrounds, which have a horizon at some r = r h where the blackening factor vanishes, f (r h ) = 0.The temperature and the entropy density are then given by the surface gravity and area of the horizon: We use the standard radial gauge A (i) r = 0.In order to study the model at nonzero quark chemical potentials, we turn on temporal components Φ i ≡ A (i) t for the gauge fields.The variables Φ i are then cyclic: the action only depends on them via the derivatives Φ ′ i ≡ ∂ r Φ i .Consequently, the gauge fields can be integrated out leading to [92] ni where n i is the number of quarks with flavor i.The total quark number and baryon number become then respectively.The gauge fields must vanish at the horizon, Φ i (r h ) = 0, while the quark chemical potentials are given in terms of their boundary values, Finally, we also need to control the flavor-dependent quark masses.They are given as the sources of the tachyon fields, i.e., the coefficients m i in the boundary expansion [73] where the AdS radius ℓ and the energy scale Λ are determined by the boundary behavior of the metric and the dilaton In these formulas γ, ℓ, and v 1 are found explicitly in terms of the near-boundary expansions of the potentials V g , V f , and κ, but Λ needs to be determined numerically [73] and is typically around Λ QCD .Note that m i do not correspond to the physical quark masses but differ from those defined, e.g., in chiral perturbation theory by multiplicative constants.
After introducing the numerical methods we employ in the V-QCD calculation below, we will briefly discuss the determination of the m s parameter by fitting the masses of pseudoscalar mesons.

B.1: Numerical method
In order to solve the backgrounds numerically, we first fix N c = N f = 3 and m u = m d = 0 and choose as the potentials the sets 5b, 7a, and 8b from [81,91].As we are studying the deconfined phase without spontaneous chiral symmetry breaking, the fields τ u and τ d vanish identically.We then proceed to solving the Einstein equations and the equation for τ s numerically by shooting from the horizon towards the boundary.Since the entropy in ( 46) is directly given in term of horizon quantities, it turns out that the natural input charges are the ratios ñi = 4πn i s = ni e 3A(r h ) ; (52) see [92] for details.The input values to the computation are then the horizon values of the scalars ϕ h = ϕ(r h ), τ sh = τ s (r h ), and the three charges ñu , ñd , and ñs .After the solution has been found, we can compute m s , s, T , µ u , µ d , and µ s as functions of the input parameters using the formulas given above.The horizon value of e ϕ can be identified as the 't Hooft coupling in this setup [75], i.e., α s = e ϕ h /(12π) for N c = 3.The value of τ sh is determined for each solution by fitting the strange quark mass to the pseudoscalar meson masses as we discuss below.As above, the relations between quark chemical potentials are solved by requiring charge neutrality and β-equilibrium after adding a free electron gas on top of the strongly coupled model.The susceptibilities χ ij are then found as numerical derivatives of the various quark number densities with respect to the chemical potentials.Here, we use a brute force method where we vary all input parameters except for τ sh on a small four-dimensional grid around a chosen point, and obtain the partial derivatives by fitting the results.

B.2: Determining the strange quark mass parameter
Finally, let us briefly discuss how the strange quark mass parameter m s is determined by fitting the masses of the pseudo-Goldstone bosons in QCD.For this we need to consider the ansatz where the phases π a are small fluctuations and t a are the generators of SU(3).The fluctuation analysis is then a generalization of the flavor-independent analysis performed in [88,89] (To be precise, one also needs to fluctuate the gauge fields, the longitudinal components of which mix with the pseudoscalars).This generalization is largely straightforward, but there is a small complication: when the masses of different quark flavors are unequal, the background and the fluctuations involve matrices in flavor space that do not necessarily commute.The general prescription of the tachyonic Dirac-Born-Infeld actions for such cases is not known, but for the purposes of this article it is enough to adopt a simple prescription: we arrange the background matrices (functions of the background tachyon field) in the fluctuation action to the left and the fluctuation wave functions (combinations of π a t a ) to the right in the single trace appearing in the action for the fluctuations.We focus on the case of equal u and d quark masses, in which case τ u = τ d ≡ τ ud , and the standard pion, kaon, and η states in π a t a form a diagonal basis for the fluctuations.In this case, the ambiguity in defining the fluctuations will only affect the kaons, as pions and the η commute with the background.
To proceed, we write the usual plane wave Ansatz for the fluctuations, arriving at the following fluctuation equations for the wave functions where there is no sum over a, and The index i takes two values, ud and s, and the coefficients c ia are determined such that Tr diag(τ ud , τ ud , τ s )t 2 a = c uda τ ud + c sa τ s .
If we adopt a normalization Tr [t a t b ] = δ ab /2, the coefficients are given by In order to fit the meson masses we construct the vacuum, i.e., the zero-temperature solution of V-QCD (see [73]) at finite quark masses and search for normalizable solutions to (55) with lowest mass −q 2 in each sector.We vary the quark masses determined through (50) and search for the best fit to the experimental values of the pion, kaon, and η masses (cf.Table I) finding agreement with the experimental values within a few percent.In the analysis of the bulk viscosity, we set m ud to zero for simplicity, and use the values of m s from this table.
As we are working in the deconfined phase at finite temperature and density, it would be better to fit the strange quark mass to some observables in this phase rather than to the masses of pseudoscalar mesons.Indeed, it has been found that simultaneously fitting the model using data from both the confined phase (such as meson masses) and the deconfined phase (such as lattice results for the EoS at zero density and high temperature) is challenging and easily leads to tensions in the fit parameters [87].In the deconfined phase, there is lattice data for quark number susceptibilities at high temperatures and small densities, showing strong dependence on the strange quark mass [90].However, fitting the strange-quark-mass dependence using this data in combination with the flavored version of the V-QCD model described above does not lead to a good fit.This is not surprising: The quark mass dependence of the flavored V-QCD model has not been compared in detail to such data in the literature so far, and there may be additional parameters that one must include in the holographic setup.
Nevertheless, we are able to make a general observation on the nature of the quark-mass corrections.The magnitude of the strange-quark-mass dependence of the V-QCD quark number susceptibilities significantly underestimates the effect seen in lattice data when using the value of the strange quark mass that is fitted to the kaon and η masses in Table I.While the lattice data describes QCD only at low densities, this observation strongly suggests that the quarkmass dependence remains underestimated at high densities.This is in agreement with our results for the quantities that depend on the parameter A 1 -vanishing at zero quark mass -in fig. 2 of the main text: While the density dependence of the V-QCD result appears to be very similar to that of the pQCD result, the overall normalization of the former is clearly lower.Furthermore, the difference is so large that no smooth interpolation between the results is possible.We are planning to carry out a detailed analysis of the fit to the quark number susceptibilities in future work.

C: Analyzing the behavior of the results
To keep the figures in the main text as readable as possible, we have relegated a finer analysis of our results to this section.Below, we in turn investigate the limiting behavior of the bulk viscosity at high densities, compare our NNLO pQCD result with a lower-order ones, and study the validity of various approximations concerning, e.g., the inclusion of only diagonal susceptibilities or only a part of the full temperature dependence of the result.At the end of the section, we also study the effect of the nFL correction to the electroweak rate, visible as the prefactor in eq. ( 2), on the bulk viscosity.This study gauges the stability of our results with respect to a future addition of currently unknown O(α s ) corrections to the electroweak rate.
C.1: High-density limit and comparison with lower-order results Starting from the high-density limit, fig. 3 shows our D3-D7 and pQCD results as functions of temperature for a range of baryon densities above those displayed in fig. 1.We observe a consistent pattern, where the peak value of the viscosity decreases and moves to smaller temperatures with increasing density, although in the D3-D7 case the peak value appears to ultimately saturate.Interestingly, below the peak temperature the value of the bulk viscosity increases as a function of density in both cases, but this behavior reverses after the peak.The density dependence is also seen to become stronger at higher temperatures, with the viscosity bands being quite tightly bundled together at lower values of T especially in the pQCD result.
For pQCD, we are also able to compare our results with lower-order ones.This is done in fig.4, which shows our NNLO result together with a calculation where one only uses terms up to O(α s ) in the weak-coupling expansion of the thermodynamic functions, including the mass correction of O(m 2 s ) ("NLO"), as well as the leading-order (LO) result for non-interacting quarks.The NLO result is very similar to an existing NLO computation [48], with small discrepancies arising due to the mass expansion scheme we use as well as the lack of a running coupling in [48].One distinctive feature of these results is a clear increase in the renormalization-scale-related uncertainty at lower temperatures, which can be traced back to the logarithmically enhanced nFL corrections to the electroweak rate, visible in eq.(2).

C.2: The accuracy of various approximations
In the main text, the bulk viscosity is evaluated using the full susceptibility matrix, which is diagonal in the D3-D7 model due to the quarks being treated in the quenched approximation but contains off-diagonal elements in the pQCD and V-QCD calculations.As we can see from fig. 5, which displays the latter two results in the diagonal approximation, the difference is, however, modest.For the pQCD calculation in particular, using the simplified eq. ( 40) leads to a result that is numerically remarkably close to the full one.Nevertheless, we use the full expressions elsewhere in this letter, as the computational gain from using the diagonal approximation is quite small.In V-QCD, the difference of the two results is more significant, with the off-diagonal terms increasing the bulk viscosity below the peak and decreasing it at higher temperatures.
Next, we briefly inspect the validity of our earlier claim that the temperature dependence of the bulk viscosity originates almost entirely from the electroweak rate in eq.(2), while ζ is nearly independent of the T -dependence of the quark densities and susceptibilities.This comparison is performed in fig.6, from which we indeed see that the full and approximate results are accurate up to O(10%) corrections both in the D3-D7 and pQCD calculations for all T ≤ 100 MeV.As argued before, this implies that we may evaluate the pressure and its derivatives at a fixed (small) temperature T , varying only the chemical potentials.
Finally, let us briefly look into the stability of our results with respect to potentially sizable QCD corrections to the electroweak rate λ 1 .This is most straightforwardly achieved by comparing bulk-viscosity results evaluated with and without the nFL correction to eq. (2).Precisely this is done in fig.7, where we observe, in accordance with the discussion of main text around fig. 2, that the peak value and shape of the results stay intact, but the value of the peak temperature shifts due to the change in λ 1 .While the qualitative contrast between our QM results and the confined-phase curves in fig. 1 is left unchanged, this clearly highlights the importance of determining the remaining O(α s ) corrections to electroweak rate.
e r o n s , 4 .5 n s a t FIG. 2. A comparison of the values of two quantities charac-terizing the bulk viscosity: its zero-frequency limit and peak value, ζ(ω = 0) and ζ peak .These quantities are multiplied by different factors so that they depend only on A1 and C1 in eqs.(38) and(39) and are independent of the oscillation frequency ω and the electroweak rate λ1, as indicated by the expressions on the right vertical axes.The error bars in these panels capture the variation of model parameters in the different models as described in the main text.

FIG. 3 .FIG. 4 .
FIG.3.The bulk viscosity computed in the D3-D7 model and in pQCD as a function of T over a large range of densities.

FIG. 5 .
FIG.5.The bulk viscosity computed in pQCD (left) and in V-QCD (right), showing this time results obtained with the complete susceptibility matrix and its diagonal approximation.In the left panel, the inset shows a zoomed-in bulk viscosity near the peak, emphasizing the small value of the difference.

TABLE I .
Fit results for the quark and meson masses in the flavor-dependent holographic model.