Inverse Anisotropic Catalysis in Holographic QCD

We investigate the effects of anisotropy on the chiral condensate in a holographic model of QCD with a fully backreacted quark sector at vanishing chemical potential. The high temperature deconfined phase is therefore a neutral and anisotropic plasma showing different pressure gradients along different spatial directions, similar to the state produced in noncentral heavy-ion collisions. We find that the chiral transition occurs at a lower temperature in the presence of anisotropy. Equivalently, we find that anisotropy acts destructively on the chiral condensate near the transition temperature. These are precisely the same footprints as the"inverse magnetic catalysis"i.e. the destruction of the condensate with increasing magnetic field observed earlier on the lattice, in effective field theory models and in holography. Based on our findings we suggest, in accordance with the conjecture of [1], that the cause for the inverse magnetic catalysis may be the anisotropy caused by the presence of the magnetic field instead of the charge dynamics created by it. We conclude that the weakening of the chiral condensate due to anisotropy is more general than that due to a magnetic field and we coin the former"inverse anisotropic catalysis". Finally, we observe that any amount of anisotropy changes the IR physics substantially: the geometry is $\text{AdS}_4 \times \mathbb{R}$ up to small corrections, confinement is present only up to a certain scale, and the particles acquire finite widths.


Introduction
Understanding all corners of the phase diagram of quantum chromodynamics (QCD) is a major focus of current research. Besides theoretical curiosity, studying QCD matter in extreme conditions is crucial in many physical situations ranging from the ultra-relativistic heavy-ion collision experiments at RHIC and LHC, to the core of neutron stars and magnetars, and to early cosmology [2][3][4][5][6][7][8][9].
According to our current understanding, colliding heavy ions create a strongly-coupled deconfined plasma state known as the quark-gluon plasma (QGP) that behaves almost as a perfect fluid [10][11][12][13][14][15][16]. In the event of off-central, i.e. with nonvanishing impact parameter, collisions the plasma is highly anisotropic 1 and is typically created in the presence of a strong magnetic field, which can reach up to eB{m 2 π " 5´10 [20][21][22][23][24][25][26]. The interplay between strong magnetic fields, strong interactions and finite temperature has been studied extensively in the literature, and is known to lead to rich phenomenology, see [27][28][29] for reviews. One of the most surprising effects in this context is the phenomenon known as Inverse Magnetic Catalysis (IMC). IMC refers to the observation that, at temperatures of the order of " 150 MeV, the presence of a strong external magnetic field has a destructive effect on the chiral condensate xqqy [30][31][32][33]. This phenomenon was first observed on the lattice, and cannot be explained by standard perturbative calculations. In fact, perturbative QCD predicts the exact opposite effect, dubbed as Magnetic Catalysis (MC) [34][35][36]. The intuition behind the Magnetic Catalysis is that, in the presence of a strong magnetic field, charged particles freeze in their lowest Landau level, effectively reducing the dimensionality to p1`1q. Since the IR dynamics in gauge theories in lower dimensions is much stronger in comparison to their higher dimensional counterparts, this leads to a strengthening of the condensate and catalysis of chiral symmetry breaking [29]. The lattice results of [30][31][32][33] on the other hand indicate that the inverse effect, i.e. weakening of the condensate arises again from the strongly coupled dynamics around the deconfinement temperature at stronger magnetic fields.
The exact mechanism that leads to IMC remains elusive today. One compelling idea [37,38] based on lattice calculations, is that IMC arises due to competition between the "valence" and the "sea" quarks in the quark propagator 2 : where ZpBq is the path integral without the propagator, the trace and the determinant are taken over the spin and the momentum space and { DpBq includes coupling of fermions both to the external magnetic field and to the gluons A a µ . For a magnetic field in the x 3 -direction The "valence" contribution arises from the quark operators inside the path integral (1.1) i.e. from the trace. The effect of B through this contribution always tend to catalyze the condensate simply because B increases the spectral density of the zero energy modes of the Dirac operator. The "sea" contribution, on the other hand, comes from the determinant that describes fluctuations around the gluon path integral. The B dependence of this contribution 1 Anisotropy is present even in the central collisions, thanks to fluctuations in the initial shape of the participant nuclei. This can be inferred from the fact that the elliptic flow parameter v2 is typically nonvanishing even at vanishing centrality, see e.g. [17][18][19]. 2 For simplicity, we consider only one fermion flavor with mass m, but the idea applies more generally.
suppresses the condensate around the deconfinement temperature [37,38]. Another idea is based on a competition between the total magnetic field dependence of the quark propagator and of the QCD coupling constant when RG scale is taken at B [29].
In this paper, we put forward and test an alternative idea: inverse magnetic catalysis results from the anisotropy in space-time caused by the presence of the external magnetic field. Imagine, as a result of some mechanism, the SOp3q rotational symmetry in space is broken down to SOp2q. This may result from a constant external magnetic field as in the lattice QCD studies above, from unequal values of spin-spin interaction constants in different directions as in certain spin models, or from the deformed geometry one puts the system on as in the anisotropic QGP state produced in the off-central collisions. In these cases one can describe anisotropy in the effective action by introducing a flat but anisotropic background metric g µν " η ab e a µ e b ν , e a µ " δ a µ for µ " 0, 1, 2, where we have chosen the anisotropy in the z direction and denoted it by W 0 . One may try to set W 0 to zero by a rescaling x 3 Ñ e´W 0 x 3 but this is not an RG invariant operation hence invalid in the full quantum theory 3 . Let us for simplicity assume that there is no external magnetic field, hence the only source of anisotropy is W 0 . The covariant derivative for the fermions now becomes 4 instead of (1.2). Note that we can separate the W 0 dependent and independent parts, just as we did with the B dependent and independent parts in (1.2). Now one can separate the contribution of W 0 in the condensate expectation value into valence vs sea quarks simply by replacing { DpBq by { DpW 0 q in (1.1): It is then tempting to ask whether one obtains a dependence of the condensate on W 0 similar to the dependence on B above. In particular, it is tempting to ask whether one observes IMC solely due to W 0 . It is not easy to answer this question directly in field theory as IMC is supposed to arise from strong coupling dynamics. One complication that is immediate to see arises from renormalization. In the renormalized theory the bare coupling W 0 will be replaced by an RG-scale dependent anisotropy parameter W pµ RG q. In this paper we will instead answer to the aforementioned question affirmatively using the techniques of the gauge-gravity duality that are suitable for studying the effects of anisotropy in the full non-perturbative system. The gauge/gravity duality is established as a powerful theoretical tool to study especially the qualitative aspects of a large class of strongly-coupled gauge theories in a completely nonperturbative manner. Several works have already approached the problem of the dependence of the condensate on the magnetic field in the holographic context, including [39][40][41][42][43][44][45][46][47][48][49][50][51]. In [45], in particular, the authors gave a heuristic explanation of the IMC inspired by the aforementioned competition between the valence v.s the sea quarks but translated in the gravity language. Just as in (1.1), there are two contributions that can be separately recognized [45] in the gravitational description as well. The first one, "valence" comes from explicit dependence of the open string tachyon equation of motion on B, that is the bulk field dual to the condensate, while the second, "sea" refers to an indirect effect coming from the backreaction of B on the geometry. The authors of [45] pointed out that is natural to identify the former explicit dependence with the valence, and the latter, implicit dependence with the sea contributions, respectively. If true, then, it would imply that the backreaction contribution is responsible for the IMC.
In this paper we consider a holographic QCD theory with no external magnetic field but with anisotropy. One way to introduce anisotropy to the system is to turn on a relevant (or marginal) operator that (i) depends explicitly on one of the spatial directions and (ii) couples only to the color degrees of freedom. Indeed, this kind of deformation has been previously considered in the context of holography, e.g. in [1,52,53] for massless quark flavors. In these papers the authors considered a θ-parameter (which sources the pseudo-scalar operator Tr F^F ) that depends linearly in one of the spatial directions, θpxq " a x 3 , as a way to introduce anisotropy into the system. Now let us see that the field theory with this spatially dependent θ term and massless quarks can also be put in a form similar to (1.2), hence the expectation value of the quark condensate can again be split into the valence and the sea parts as above. Consider the generating function of QCD both with a nontrivial θ term and an external axial gauge field A 5 : Dq DA a e´ş LrA a ,qs`A 5¨J 5`θ Tr‹F^F (1 .6) where LrA a , qs is the Lagrangian for the massless QCD and J 5 is the anomalous chiral current. We do not turn on an external electric gauge field for simplicity. Let us call the anomaly coefficient c a i.e. we have the non-conservation equation This generating function enjoys invariance under the generalized chiral transformation 5 Therefore a nontrivial space dependent θ term, θ " ax 3 can be set to zero by turning on an external axial gauge field A 5,µ " a{c a δ 3 µ . This means that the action expectation value of the quark condensate in the theory with θ " ax 3 and A 5 " 0, which we considered in the previous paragraph, can be written as This is again in the form 1.2 where the anisotropy enters the inverse propagator linearly and the contribution of a to the quark condensate can be divided into the valence and the sea parts as above.
Anisotropic confining gauge theories were revisited in the holographic approach recently in [1]. One of the important lessons of this paper was that, in the color sector, the anisotropic deformation reduces the confinement-deconfinement phase transition temperature. Since the effects of B on confinement are qualitatively the same, this result reinforced the intuition of [45], and led to the conjecture that anisotropy by itself could explain the phenomenon of IMC. 6 In order to study the behavior of the chiral condensate xqqy in the presence of anisotropy we have to consider an extra flavor sector on top of the model in [1]. Alternatively, we can introduce the same anisotropic deformation in the models originally considered in [45], at zero magnetic field. The difference between these two approaches boils down to the choice of potentials for the dilaton field. We choose to do the latter, because the choice of potentials is better motivated than in the former models. In this case, the color sector of the theory is taken to be Improved Holographic QCD (IHQCD) [54,55]. This is a bottom-up Einstein-Dilaton theory with a specific potential for the dilaton, which mimics many of the phenomenological signatures of QCD. On top of this theory, we also consider a flavor sector based on a pair of space filling D4´D4 branes [56,57]. However, since flavor physics is suppressed in the large N c limit, one must consider an appropriate limit in order to properly take into account the backreaction of flavors. Specifically, one must take both N c Ñ 8 and N f Ñ 8, while keeping their ratio x " N f {N c fixed. This is known as the Veneziano limit, and defines the V-QCD model [58] which is the model we use as the holographic dual of QCD in this paper.
The paper is organized as follows. In section 2 we start by giving a brief overview of the model, discussing in detail the color and flavor sectors mentioned above, as well as presenting the relevant equations of motion and constraints. In section 3 we discuss the IR asymptotics in detail and show, in particular, the drastic effects induced by the anisotropic deformation. In section 4 we solve numerically the equations of motion and find the relevant anisotropic black brane solutions. We also study the thermodynamics of the models by working out the free energy in the canonical ensemble and discussing in detail the role of the anisotropic deformation. In section 5 we compute various observables of physical interest. First, we devote our attention to study the chiral condensate, which was the original motivation of the paper. In addition, we study the meson and glueball spectra, quark-antiquark potential and entanglement entropy, to further characterize the behavior of the new IR fixed points. We close in 6 with a discussion of our results and some outlook.

Holographic setup
The holographic model we will consider has two parts, the gluon sector and the flavor sector. The gluon sector is based on the so-called improved holographic QCD model (IHQCD) [54,55], and the flavor sector is defined in terms of a generalized tachyon Dirac-Born-Infeld action arising from a pair of space filling D4´D4 branes [56,57]. The two actions fully backreact in the Veneziano limit, which defines the V-QCD model [58]: The gluon sector contains a finite set of bulk fields dual to relevant or marginal operators that dominate the dynamics in the IR. Among these we have the stress-energy tensor T µν , which is dual to the metric g µν , the glueball operator Tr F 2 dual to the dilaton λ and a pseudo-scalar operator Tr F^F dual to the axion χ. The latter operator is introduced in order to break isotropy as in [1,52,53]. Notice that proper implementation of the QCD axial anomaly would require coupling of the axion to the flavor sector which is of the leading order in the Veneziano limit [59,60]. As we will only use the axion to break the isotropy, considering such couplings is not necessary and we omit them for simplicity. Finally, the flavor sector includes an additional field, the tachyon τ , which is dual to the quark bilinear operatorqq.
Constraints to the potential functions and couplings in the action from various sources have been discussed in detail in earlier literature [54,55,[58][59][60][61][62][63][64][65]. In the current study, the coupling Zpλq between the dilaton and the axion is taken from [62,66] while the other potentials are taken from [63,67]. Explicitly, the potentials are given by

6)
Zpλq " 1`λ 4 10 , Notice that we have set the overall constant in Zpλq (Z 0 in the notation of [62]) to unity, since it can be reabsorbed in the normalization of χ.
Our Ansatz for the metric and other bulk fields is the following: This Ansatz automatically satisfy the equations of motion for the axion χ, while introducing anisotropy in the x 3 direction. Moreover the dependence of metric field W on the holographic coordinate r precisely corresponds to the dependence of the renormalized anisotropy parameter on the RG scale discussed in the introduction. The Einstein equations that follow from the action are: There is also a first order constraint, which is given by (2.12) The equation of motion for the dilaton is Finally, the equation of motion for the tachyon is (2.14) We will solve this set of equations numerically. To do this, we make use of a scale symmetry present in the equations of motion: This corresponds to the scale symmetry of the action on the field theory side. The background solution has a nontrivial dependence on the bulk coordinate r, which breaks this symmetry and introduces an energy scale Λ, analogous to Λ QCD . The precise definition of Λ will be given in section 4. Before discussing the numerical solutions, we will study in detail the IR structure of these equations and obtain analytic solutions in this regime.

IR behavior
Let us then discuss the asymptotic geometry and RG flow at zero temperature (i.e., for f " 1) in the IR. As it turns out, turning on any nonzero anisotropic parameter a changes the IR structure drastically. Motivated by the fact that in the isotropic case the quarks are either asymptotically decoupled in the IR in the chirally broken phase or affect the gluon dynamics only trivially in the symmetric phase [58,59], we start by considering the system without quarks, i.e., taking the limit x Ñ 0 above. It is useful to write the equations of motion in another form. We define λ " e φ , dr dA e A " q "´e p , and Ă W " W`A. In terms of them, the three independent equations of motion are where dots denote the derivatives with respect to A.

AdS 4 IR fixed point in the chirally symmetric phase
We start by discussing the exact fixed point solutions which are realized in the chirally symmetric phase in the zero temperature limit. First we notice that the above equations of motion seem to admit an exact fixed point solution, determined by the equations e 2p˚V g pλ˚q " 9 , 3a 2 Zpλ˚q " 2V g pλ˚q (3.4) and with Ă W a constant which we set to zero (as it can be absorbed in a). This fixed point (without additional requirements) is however not realized as an endpoint of any holographic RG flow. The reason can be seen as follows. The second order dilaton EoM may be written as The term in the latter square brackets vanishes at the fixed point (because also 9 Ă W " 0) but the term in the first square brackets does not. After dividing by 9 φ, the first two terms already imply that : φ is finite. So even if the EoM is satisfied when 9 φ vanishes exactly, any small perturbation will lead to a sizable : φ and therefore fast deviation from the fixed point. In particular, the fixed point cannot be reached asymptotically for A Ñ˘8.
If in addition to the conditions (3.4) we impose (again taking Ă W " 0) the issue with the flow of the dilaton is gone and the fixed point is stable and physical. Notice that this condition can also be written as so that the flow equations (3.13)-(3.17) also have trivial solutions. The additional condition cannot however be satisfied for a generic a, but only for some specific value which we denote by a˚. As one can check, for the potentials (2.4)-(2.6) there is no such solution (excluding the runaway solution at φ˚" 8).
The situation is different if we consider the backreaction of the flavors in the chirally symmetric phase, τ " 0. In this case the EoMs for the glue are obtained by replacing V g by the effective potential V eff " V g´x V f 0 [58]. As it turns out, after the replacement and for the potentials specified above, a nontrivial fixed point solution pφ˚, p˚, a˚q does exist for x À 1. As we show below, this fixed point is indeed realized the IR limit of the T " 0 RG flows in the symmetric phase for x " 1{3.
As p˚is fixed, the resulting geometry is AdS 4ˆR : e Aprq {A 1 prq "´e p˚i s solved by e A " e p˚{ r, and the warp factor of dx 2 3 is e Ă W " e A`W " const. as we pointed out above.

Rolling IR fixed point in the tachyonic phase
The low temperature geometries in the chirally broken phase (and also for the pure glue case x " 0) have an interesting structure which is drastically different from that of the isotropic solutions. In order to analyze them, we start by studying variations around the exact fixed point discussed above, which leads to a "slow roll" behavior as we will now demonstrate. Without loss of generality, we may take Ă W pA " 0q " 0 and discuss the evolution of the system near A " 0. The "slow roll" will be driven by Eq. (3.2). We substitute the following Ansatz in the equations of motion: where all coefficients are taken to be small and the fixed point values may depend on them, p˚"p˚pC i q andφ˚"φ˚pC i q. As it turns out, the EoMs are satisfied to first order in the coefficients C i if the proportionality Zpφq 9 V g pφq 3 holds at least for the exponential terms in φ, in consistency with (3.7). The slowly rolling functions ppAq and φpAq satisfy the equations up to corrections OpC 2 i q. This implies in particular that C p "´C W {2 and that C φ " C W V g pφ˚q{V 1 g pφ˚q. The conditions forp˚andφ˚are obtained by setting A " 0 and they are corrected by OpC W q terms with respect to (3.4). Thereforep˚andφ˚approach the fixed point values p˚and φ˚defined by (3.4) as C W Ñ 0. The value of C W can be related to the deviation from the law Zpφq 9 V g pφq 3 but this requires considering higher order corrections in C W which we will not do here. Instead, we will present a more precise way of discussing the IR flow as follows.
The IR flow will actually determined by an fixed point which is slowly moving due to the flow of, say Ă W . In order to guarantee that the flow stays at the fixed point, it is sufficient, to a very good precision, to simply set all second derivatives with respect to A to zero. Notice that (3.3) then can be solved for Ă W 1 , and the system can be rearranged to read For these to lead to a consistent flow (rather than being satisfied only at a single point), we require that their derivatives are also satisfied, imposing the approximation that second derivative with respect to A are set to zero. This leads to two new equations which are equivalent to the above equations at certain values of φ and q. This then determines the flow. Since the potentials are possibly complicated functions of φ, it is convenient to solve a (or actually ae´Ă W ) and treat φ as a parameter. That is, the above three equations and their derivatives lead to five independent equations which we solve for φ 1 , p 1 , Ă W 1 , p, and ae´Ă W . There is a trivial solution without flow given by Eqs. (3.4). The nontrivial solution is given by Remarks: • The equations (3.13)-(3.15) are consistent with the relations of the coefficients C i obtained above in Sec. 3.1 when Zpφq is roughly 9 V g pφq 3 , i.e., d dφ log Zpφq´3 d dφ log V g pφq is small. Now the value of C W is fixed and can be read from (3.14).
• For the flow to present an acceptable IR asymptotics we need such potentials that φ 1 ă 0 in (3.13).
• Interestingly, inserting to the flow equations the exactly exponential potentials V " 8{3γφ reproduces the exact scaling solution of [1] even when deviation from Zpφq 9 V g pφq 3 is sizable. This is because corrections to the flow involve higher order logarithmic derivatives which vanish for exactly exponential potentials.

Behavior at asymptotically large λ
The potentials in Eqs. (2.4)-(2.7) are such that the IR solutions are very precisely described by the above flow equations because the proportionality Zpφq 9 V g pφq 3 is only violated by logarithmic corrections.  We first discuss the asymptotic IR behavior of the geometry which can be solved analytically by using the flow equations. This has been worked out for more generic potentials and higher order corrections in [68]. In order to find the asymptotics, we substitute the potentials in the flow equations, i.e., we take V g 9 V IR e 4φ{3 ? φ and Z 9 Z IR e 4φ . We obtain up to corrections suppressed by 1{φ. Integrating these equations, we obtain the asymptotics for the metric factors and φ: (3.20) These formulas describe an approximate AdS 4ˆR geometry with multiplicative corrections of the form e # ? log r . It is instructive to write down the string frame metric: Notice the cancellation of the square roots in the warp factor, after which the first term is the AdS 4 metric with multiplicative logarithmic corrections. equations for the pure glue case (x " 0) and various values of the asymmetry parameter a.

Recall that the IR (UV) is at
The exact numerical solutions of the full equations of motion are given by the solid curves whereas the flows from Eq. (3.13) are given as the dashed curves. We took the initial conditions for the latter flows from the exact solutions at A "´10. When a{Λ ! 1, the structure of the solution is as follows.
• For A " 0, the geometry has the same UV asymptotics as in the absence of asymmetry [54,55], i.e., the RG flow is given by φ "´log A.
• For A ! 0 and when φ ! φ˚there is an intermediate regime where the background follows the "standard" IR asymptotics φ " 3 2 A of the a " 0 case [54,55]. Here the fixed point value φ˚is given by Eqs. (3.4): for small a we have that φ˚"´3plog aq{4.
• For A ! 0 and φ Á φ˚the flow is the "slow roll" described by Eq. (3.13). As A decreases, after a short transition regime, the exact numerical solution is practically indistinguishable from that given by Eq. (3.13).
When a " Op1q the intermediate region is absent. As a Ñ 0, the region with the intermediate behavior grows. The "slow roll" region is pushed deeper and deeper in the IR and disappears in the limit. Therefore the geometry approaches the a " 0 thermal gas geometry, but the convergence is not uniform. In Fig. 1 (right) we compare the RG flow of the coupling given by Eq. (3.13) and to the asymptotic formula in Eq. (3.20) for various "boundary conditions" near A " 0, as determined by the value of a. The color coding for the values of a is the same as in the left hand plot. We notice that the analytic asymptotic formula works reliably only at extremely high values of´A.
Next, we make the following consistency checks on the discussion above. In Appendix A we verify that for the potentials used in this article the tachyon indeed diverges fast enough in the IR for it to decouple so that the above analysis is consistent even for the tachyonic solutions at x ą 0.
We also plot the string frame Ricci scalar which governs the higher order stringy corrections in Fig. 2. We see that |R s | remains small in the IR in our region of interest so that the corrections are small and our approach is consistent. Notice however that at extremely deep in the IR the Ricci scalar grows as shown on the right hand plot and eventually blows up. This behavior is due to the logarithmic term in the AdS part of (3.21). Therefore the zero temperature geometries corresponding to the flow would receive large stringy corrections. Notice also that the limit of a Ñ 0 appears nonuniform as one can see from the left hand plot in Fig. 2: as a decreases the amplitude of R s in the IR grows but the "slow roll" regime, where R s is sizable, is pushed deeper in the IR.
Finally, it is important to understand what the divergence of the Ricci scalar means for the finite temperature solutions which we consider below. Since the IR geometry is close to AdS 4 , small black holes will have temperatures T " 1{r h " e A h where the subscript "h" refers to the value at the horizon. Therefore the finite temperature geometries are consistent, that is free of stringy corrections, down to exponentially small temperatures: for example the Ricci scalar reaches the value´2 around A "´4000 which gives the temperature T {Λ " 10´1 700 . We conclude that we can trust the finite temperature solutions down to an extremely small temperature scale, which for any practical application can be taken to be zero.

Thermodynamics
To investigate thermodynamics of the system we consider the black brane solutions with a nontrivial blackening factor f " f prq in (2.9) which vanishes at some horizon r h , f pr h q " 0. r h then parametrizes the different black-brane solutions with different temperatures T " f 1 pr h q{4π. We solve the equations of motion (2.10), (2.13) and (2.14) numerically by shooting from the horizon towards the boundary. More specifically, we set the data of the nonnormalizable modes according to We shoot from the horizon until the desired value of m q is reached -the sources for the other fields can then be set to be those of (4.1) by using symmetries of the solution [45,58]. In particular, the symmetry Eq. (2.15) can be used to set 7 the value of the source Λ of the dilaton field, which characterizes the energy scale of the solution (in analogy to Λ QCD on the field theory side). From each solution we read off the value of the normalizable modes and extract from them the expectation values of the operators dual to the various bulk fields. For simplicity we will focus only on solutions with m q " 0. Other values of m q are straightforward to obtain but the extraction of the quark condensate is numerically more demanding in these cases.
In order to study the thermodynamics of our solutions we also need to compute the free energy F in the canonical ensemble. In the holographic description, this amounts to evaluating the on-shell Euclidean action (2.1), appended by the standard Gibbons-Hawking and counterterms [55,64]. In practice, however, we compute the background subtracted free energy directly by integrating the first law of thermodynamics while keeping all sources fixed. Before proceeding further it is worth noting that more generally, given the bulk fields of our gravitational action, the Ward identity for the boundary stress-energy tensor reads: where O α is the operator dual to bulk field α. From the form of our Ansatz (2.9) we know that xO χ y " 0 for our solutions, so the third term in (4.3) is generally absent.
In Figure 3 we plot the free energy F {Λ 4 as a function of T {Λ for various values of a{Λ in the case where x " 0 so the flavor physics is suppressed. Other values of x behave qualitatively similarly 8 so we will omit their discussion for simplicity. We observe the following general behavior: i) for small T {Λ there are always two solutions. First, there is the ground state heated up to temperature T , which is referred to as the thermal gas solution. The gravitational background dual to this state is obtained from the black brane solution (2.9) by sending f prq Ñ 1 and compactifying the Euclidean time direction with period 1{T . We take this solution as our reference background for the free energy computation so it corresponds to the horizontal axis F " 0 in the figure. 9 Second, we also observe a black hole solution which, in fact, dominates the ensemble at small T {Λ (for non-zero a). The fact that there is a black hole branch that dominates at small T {Λ is in stark contrast with the standard, isotropic models of holographic QCD, for which the thermal gas solution is the dominant one at low temperatures. 10 Indeed, this behavior at low T {Λ is crucially related to the fact that, for any non-zero a, the IR structure of the theory is drastically modified by the anisotropic deformation, approaching a nearly AdS 4 geometry in the deep IR (with logarithmic corrections). At very low temperatures, T {Λ ! 1, the free energy scales as F 9´T 3 while the entropy density scales as s 9 T 2 . ii) The free energy exhibits a swallow tail behavior at intermediate values of T {Λ. In this range of temperatures, besides the thermal gas solution, there are 8 As we shall show below, for x ą 0 there are some additional features which do not affect the main points discussed here. 9 States dual to black hole solutions have free energies of order OpN 2 c q, while the thermal gas solution has a free energy of order OpN 0 c q. 10 This is in analogy to the results for the canonical ensemble of charged black holes, where the anisotropy parameter a plays the role of the charge [69][70][71]. The swallow tail behavior in these cases is associated to a Van der Waals liquid-gas phase transition. up to three black hole solutions with free energies of order OpN 2 c q. Two of these solutions are sub-dominant in the ensemble, and the third one dominates in this regime. iii) At large enough temperatures there is only one black hole solution -in addition to the thermal gas solution-which dominates the ensemble. Since the geometry is asymptotically AdS 5 the free energy and entropy density at very high temperatures, T {Λ " 1, scale as F 9´T 4 and s 9 T 3 , respectively. Altogether, if one follows the dependence of the dominant phase as a function of T , one finds a single first order phase transition connecting two different branches of black hole solutions, I and II as shown in Figure 3, that dominate at small and large temperatures, respectively.
In the limit a Ñ 0 the free-energy of the black hole solutions I in figure 3 apparently vanishes which suggests that they approach the thermal gas solution. This is indeed the case: the horizon is pushed deeper and deeper in the IR as a decreases, and disappears from the limiting solution. Therefore the thermodynamics approaches smoothly that of the standard symmetric (a " 0) IHQCD [72,73] and the black hole phase I is replaced by the thermal gas solution. This is also consistent with the analysis of the zero temperature IR geometry in the same limit in section 3.2.2.
In Figure 4 we plot phase diagrams for some representative values of x " N f {N c . Specifically, the values that we consider are the following: x " 0, x " 1{3, x " 2{3 and x " 1, respectively. In these diagrams, the black solid lines correspond to a first order transition between two black hole solutions. Such a transition is present for the x " 0 case (discussed in the previous paragraph) regardless the value of the anisotropic parameter a{Λ. We notice that for x " 1{3 and x " 2{3 this transition eventually becomes second order at large enough anisotropies, while for x " 1, one can distinguish two transitions, one of first order (which disappears at large anisotropies) and another one of second order. The second order transitions are shown as black dashed lines. Curiously, the second order line eventually becomes first order again for large values of a only in the case x " 1. If one were to plot the analog of figure 3 for x ą 0, the swallowtail structure would disappear at the value of a for which the first order transition disappears. Notice that the results at a " 0 are consistent with the earlier findings [45,48,63,67,74,75] in the model with the same choices for the potentials. We also distinguish between chirally symmetric and chirally broken phases in figure 4, depending on the value of the quark condensate xqqy (whether it is zero or non-zero, respectively). For the case x " 0 the chiral symmetry is determined through the solutions to the DBI action in the probe limit, x Ñ 0`. See section 5.1 for a more detailed discussion on the quark condensate and the role of the anisotropy. The IR behavior at finite x is generally dominated by the slow rolling tachyonic phase discussed in section 3.2. The only exception is the x " 1{3 case, for an intermediate regime of a{Λ, roughly between a{Λ " 1 and 3. In this range of anisotropies, the low temperature regime realizes the chirally symmetric IR fixed point discussed in section 3.1. Finally, following the criterion discussed in [55,76], we indicate whether or not the quark-antiquark potential (for quarks separated in the longitudinal or transverse directions) has a linear behavior, indicating confinement in both directions, anisotropic confinement or deconfinement. These different phases are separated by black dotted lines, indicating a crossover. We observe that as x is increased, the value of a up to which we have confinement-according to this criterion-decreases. We discuss the explicit calculation of the quark-antiquark potential in section 5.3.
An interesting comparison can now be made between figure 4 and the phase diagrams of [45]. In the latter case, the chiral transition temperature first decreases as a function of the magnetic field B, and then at very large values of B, it starts increasing again. In figure 4, this behavior is qualitatively similar, with the anisotropy parameter taking the role of the magnetic field. Note indeed that we see the chiral transition temperature first decrease as a function of a, then increase beyond some large enough a. This draws an interesting parallel between the effect of a magnetic field on the chiral transition and that of an anisotropy. In the following subsection we make this analogy more concrete by studying the analog of the magnetic susceptibility. Further in section 5.1 we demonstrate that the chiral condensate also behaves in the same way in the presence of an anisotropy as it does under the influence of magnetic field. In particular, we find an effect analogous to the 'inverse magnetic catalysis' [30,32].

Anisotropic susceptibility
In analogy to magnetization in a theory with nontrivial external magnetic field, we define the response to the anisotropy parameter a as M a "´B F Ba .
In the absence of a better name, we call this object the "anisotropization". Let us also define which, at a " 0, coincides with the usual definition of a susceptibility. 11 Holographically, this 'anisotropic susceptibility' can be computed as follows. Because the free energy is given by the on-shell action, it is enough to consider the derivative of the explicit a dependence in (2.1). Substituting χ " ax 3 we find that This integral is however UV divergent, and needs to be renormalized [78]. We discuss the renormalization procedure in Appendix B. The result obtained by evaluating the renormalized integral numerically for x " 0 can be seen in Fig. 5. Notice that χ a is not well-defined in the confined phase (I) when a " 0, as the integral (4.4) is IR divergent. We see that χ a decreases without limit as a Ñ 0 at small temperatures as demonstrated in the bottom plot. Notice that the susceptibility contains a scheme-dependent piece 9c 1`a 2 c 2 (see Appendix B), but the divergence is scheme-independent.
As was done in [45,47], one can then derive the following relation for the transition temperature of a first order transition as a function of a: where ∆χ a and ∆s are the differences between the two phases of the susceptibility and the entropy density, respectively. From this equation we can relate the sign of the slope to the jumps of the derivatives of χ a and s. As seen from Fig. 5, for x " 0 and 0 ď a{Λ ď 1 the jump ∆χ a is positive which agrees with T c decreasing with a{Λ in this region in Fig. 4. We have checked that at large values of a{Λ, where T c increases with a (see Fig. 4), ∆χ a also has the opposite sign.

Chiral condensate and inverse anisotropic catalysis
Once a numerical solution for the metric and other bulk fields is obtained, with boundary conditions as in (4.1), we can extract various observables of interest. In this section we will start by studying the chiral condensate xqqy, which can be computed from the normalizable mode of the tachyon. More specifically, for m q " 0, we obtain that the leading behavior of the tachyon near the boundary is given by: with b 0 " 1 3 p11´2xq and γ 0 " 3 2 . In Figure 6 we show the dependence of the chiral condensate xqqy{Λ 3 on T {Λ for various choices of a{Λ and x. The x " 0 results here are obtained by solving the tachyon equation in the probe limit.  Fig. 3). Bottom: χ a at zero temperature as a function of a{Λ.
One can see that for temperatures above the chiral transition the condensate vanishes, and below the transition it is nonzero. The order of the phase transitions can also be seen by looking at whether the condensate is continuous across the transition. For x " 0, we observe that the chiral condensate decreases with increasing a for all temperatures. For finite x, we observe two effects: The condensate first decreases with a, and then at larger a, it increases again. 12 In order to study this behavior in more detail, it is convenient to define ΣpT, aq " xqqypT, aq xqqyp0, 0q , ∆ΣpT, aq " ΣpT, aq´ΣpT, 0q.
This combination is analogous to the quantity which has been computed on the lattice at finite magnetic field (see for example [32]). In figure 7, ∆Σ is plotted as a function of a for different temperatures, and for x " 1. We see first a decrease of the condensate with a for small a followed by an increase at larger a. The decrease of the condensate is most pronounced for temperatures close to the chiral transition. This is in close analogy to the analysis of [45], which is the same model as in this work, setting a " 0, and including a magnetic field in the flavor action. There, it was found that for small temperatures, the condensate always increases with a magnetic field B, while for temperatures close to the chiral transition, there is first a decrease of the condensate with B, and then an increase. 13 The decrease of the condensate (inverse magnetic catalysis) with a magnetic field was discovered on the lattice [30,32]. The analogous behavior we observe here is therefore evidence for the claim made in [1], namely that a possible cause for the inverse magnetic catalysis is the anisotropy, which can be induced by the magnetic field as in the lattice studies or explicitly as we do in this paper. What we see is that just the presence of anisotropy has the same effect as the magnetic field. We thus call this behavior "inverse anisotropic catalysis".

Particle spectra
We observe from the phase diagram, figure 4, that anisotropy may destroy confinement even at relatively low values of a{Λ. In this section we analyze this phenomenon further by discussing how the meson and glueball states (at zero temperature) melt as a is increased. Dissociation of mesons due to anisotropy has also been studied in other holographic models [79,80]. We will discuss the particle spectrum in the helicity two glueball tower which is relatively easy to analyze. We have also computed that the spectral functions in other sectors, including all flavored states, and helicity one flavor singlet states, and find a similar behavior to that of the helicity two glueballs. The fluctuations in these additional sectors are presented in Appendix C. 13 In [45], there is also a direct coupling of the magnetic field to the condensate, in addition to the backreaction through the geometry. This direct effect was found to always increase the condensate, in agreement with lattice results [37]. An analogous direct coupling is absent in this work.
In order to analyze the helicity two glueballs we write an Ansatz for the perturbation of the metric as δg 12 " δg 21 " e 2Aprq e ikµx µ hprq where we the sum in the plane wave term goes over the time and space indices. By boosting in the x 1,2 directions (and assuming that |ω| ą |q K |) we can transform the momentum to the form q µ " pω, 0, 0, qq. We find the following equation for the fluctuations in this frame: The Schrodinger form is obtained by defining hprq " e´3 Aprq{2´W prq{2 ψprq: Inserting here the asymptotic IR solution (3.20) from above, i.e., AdS 4ˆR with logarithmic corrections, we find that V S prq " 2{r 2 as r Ñ 8, indicating that the spectrum is continuous. That is, the glueballs have finite widths even at zero temperature, signaling the possibility of the decay to the AdS 4 vacuum. Spectral density of the helicity two glueballs can be obtained from the correlator of the energy momentum tensor T 12 , which is determined (in momentum space) through the UV coefficient of the IR-regular solution to (5.3). The UV expansion of the properly normalized solution is given by where Gpω, qq is the correlator. For positive ω, definition of IR regularity is obtained via analytic continuation from the upper half of the complex ω-plane which picks up the solution with the IR asymptotics ψprq 9 e iωr (i.e. the sign in the exponent is plus) so that hprq 9 e´3 Aprq{2´W prq{2 e iωr . Notice moreover that the higher order corrections to the source term in (5.6) are real for real ω. This is the case because then the fluctuation equation has real coefficients, and the only source for complex behavior is the IR boundary condition which only affects the integration constants (i.e., the correlator) in the UV expansion. Therefore Im hprq " Im Gpω, qq r 4 for the solution with the normalization h Ñ 1 in the UV. The coefficient of the r 4 term here is the spectral density which can be therefore extracted unambiguously. We show the resulting spectral density, normalized by ω 4 , for x " 0 and for three choices for the values of a in Fig. 8. As expected, the glueballs have finite widths even though we are working at zero temperature. At a{Λ " 10´5 the first 3-4 states are clearly peaked. As a is increased, the widths of the glueballs grow, and for a{Λ " 1 the lowest glueball is already very  wide while the other peaks have melted away. Despite this the spectral density is still heavily suppressed in the regime corresponding to the mass gap at a " 0, i.e., for 0 ă ω{Λ À 4. We also note that the peak arising at ω " 0 is due to the IR behavior Im Gpω, 0q 9 ω 3 which reflects the IR geometry of the system being approximately AdS 4 .
For higher values of x the spectral density is qualitatively similar to x " 0. We show the dependence of the width (more precisely the full width at half maximum) of the lowest glueball mode on a for x " 0, 1{3, 2{3 and 1 in Fig. 9. The widths first increase when x grows from 0 to 1{3 but then decrease as x is further increased.

Quark-antiquark potential
In the absence of anisotropy in the low-temperature phase, our model is known to show linear confinement, i.e. the quark-antiquark potential V grows linearly with the separation L of the test quark and the antiquark for large enough L. A potential of such a shape prevents colored particles from escaping their bound states and becoming free. On the other hand, in section 5.2 the glueball states were shown to melt at zero temperature in the presence of anisotropy indicating that quarks can in fact escape their bound states in an anisotropic vacuum state. In this section, we investigate the mechanism behind this decay in some more detail by studying quark-antiquark potential at zero temperature. See also [81,82] for the analysis of the potentials in an anisotropic setup in slightly simpler models.
In holography, the quark-antiquark potential is a sum over saddle points, which has several terms contributing [83]. In principle, one should take all of these into account, but as a first approximation, we work in the α 1 Ñ 0 limit, where only the one with the smallest action contributes. Also, we choose to neglect the graviton exchange contribution that plays an important role by maintaining smoothness of the Polyakov loop two-point function in L [83], because this contribution does not have a qualitative effect on our results. The remaining terms can be computed by evaluating the on-shell Nambu-Goto action for a static string in the 5D spacetime described by the string frame metric, with the endpoints of the string attached to the AdS boundary [84,85]. Following [55,76], we find that at large L, the quark-antiquark potential parallel to the anisotropy V grows linearly with L if A S`W {2 has a minimum, with A S " A`2 3 log λ the string frame scale factor. Similarly, we find that the quark-antiquark potential perpendicular to the anisotropy V K grows linearly with L at large L if A S has a minimum. In principle, one could also look at the potential in arbitrary directions, but for lack of symmetry this is more difficult, and will be left for future work.
The presence of a minimum in A S resp. A S`W {2 described above is the criterion used to label a phase as confining in the phase diagram (Fig. 4). However, while this indeed gives some indication of where in the phase diagram there is confinement, it does not give the full picture. This can be seen in figure 8, where we see that the first glueball is already no longer absolutely stable in regions that satisfy the confinement criterion used in figure 4. In figure 8, we also see that the first glueball peak is still very narrow for anisotropies orders of magnitude larger than the anisotropy for which the minimum of A S p`W {2q disappears. To address this discrepancy, one has to compute not just the large L behavior of V , but the entire function.
Instead of computing V as a function of L, we can obtain both V and L as functions of  Figure 10. Quark-antiquark potential in the direction parallel to the anisotropy V , and in the direction perpendicular to the anisotropy V K , for x " 1. Note that each curve has multiple branches, and that every curve has a branch which is identically zero. Also, for a{Λ " 10´6, 10´5, the linear branch of V continues to infinity, while for V K this only happens for a{Λ " 10´6. Note that not just the branch with the smallest action is plotted. In the approximation we consider the quark-antiquark potential at separation L is given by the lowest branch.
the worldsheet turning point in the bulk r F in the following way [55,76]: with T f the string tension. The last term in (5.8) is the UV regulator which in our convention equals (twice) the action of a straight string hanging from the boundary down to the IR. Similarly, V K can be obtained using the same formulas with W set to 0. Results of these computations for x " 1 are shown in figure 10.
One immediate observation is that V has multiple branches. For each L, the smallest V corresponds to the globally stable branch. Note also that there always exists a zero branch, corresponding to two detached portions of string which have fallen into the deep IR [83]. The existence of this zero branch is due to the IR behavior which is, as we have seen, qualitatively different from the case without anisotropy. In the isotropic case the string falling into the deep IR has a diverging on-shell NG action, whereas in the anisotropic case the on-shell action is zero. 14 The existence of the zero branch is important, because it means that even though for small enough nonzero anisotropy A S has a minimum, indicating linear confinement according to [55], the branch of solutions corresponding to this case may be globally unstable, meaning that bound states can decay. We remark that our potentials at nonzero a are similar to those contructed in part by deep learning methods in [86], but the regime where the potential is linear is much shorter.
One can now also explain why the first glueball peak is still very narrow at anisotropies for which A S does not exhibit a minimum. The lowest lying bound state decaying to free particles corresponds to the ends of the string starting close together and ending up freely moving a large distance apart, so that the zero branch is the stable one. In almost 15 all cases we examined, the melting of the first glueball peak occurs when the nontrivial branches of V K disappear. Therefore we conjecture that the presence of nontrivial structure makes the worldsheet corresponding to glueball decay have a large action S decay associated to crossing this nontrivial structure. In other words, the nontrivial structure in V K prevents the decay of the first glueball. This would then contribute to a lifetime τ 9 exppS decay q, which, if S decay is large enough, will make the excitations very long-lived. To support this conjecture, we have to find the corresponding time-dependent worldsheet solutions. This is a hard problem and we leave it for future work.
Finally, we would like to comment that the computation described above assumes that the string lies in the pX r , X 3 q-plane when computing V , and in the pX r , X 1 q-plane when computing V K . In principle, while these solutions indeed solve the equations of motion, they could be unstable in the neglected directions. This would lead to yet more branches of solutions. Again, we leave exploration of these possibilities to future.

Entanglement entropy
Another interesting observable probing the IR structure of the zero temperature geometry is the entanglement entropy. The quark-antiquark potential is determined by a minimal length of a string stretching in the bulk, where the length is computed in the string frame, whereas the entanglement entropy arises from a similar minimization procedure for a higher dimensional surface in the Einstein frame [87]. We compute the entanglement entropy for two different regions of the boundary: • A region defined by 0 ă x 3 ă L. Note that the x 3 direction is parallel to the direction of anisotropy. We therefore denote the entanglement entropy of this region by S E, .
• A region defined by 0 ă x 1 ă L. Since x 1 is perpendicular to the direction of anisotropy, we denote the entanglement entropy of this region by S E,K .
14 It is exactly zero because of the way the on-shell action is renormalized. This renormalization has no physical effect. 15 For x " 1{3 the nontrivial structure of V never goes away, but glueballs still melt. Note that this is also the value of x for which we have the different structure in the phase diagram.
The problem of determining the entanglement entropy reduces to finding the minimal area of a three dimensional spatial surface in the Einstein frame metric stretching between x 3 " 0 and x 3 " L for the first region, and between x 1 " 0 and x 1 " L for the second region. Subtracting the UV divergence in the same way as was done for the quark-antiquark potential, one finds both S E, and L as functions of the turning point in the bulk r F : 16 10) with A V an infinite factor that arises because the region is spatially infinite in two dimensions. The analogous formulas for S E,K read The result of this computation for x " 1 is shown in figure 11. Note that in the figure also the unstable branches are shown. The entanglement entropy corresponds to the lowest branch.
One can see that the results for values of a{Λ up to about 0.1 are similar to isotropic results (see, e.g., [88,89]). When the RT surface probes the deep IR, so near where in our normalization the entanglement entropy vanishes, the result depends heavily on the orientation of the entangling region even small values of a{Λ. In particular, for S E, , the curve of the connected surfaces always reaches pL " 0, S E, " 0q, where it connects to the branch of the disconnected solution which have S E, " 0 for all values of L in our normalization. In the case of S E,K , however, the point L " 0 " S E,K is not reached for any nonzero values of a{Λ. The swallowtail structure present at zero anisotropy quickly gets smaller, and eventually vanishes entirely as a{Λ grows.
When a{Λ is roughly Op1q the result starts being very different depending on the orientation of the region. For S E, , the entanglement entropy crosses zero for a smaller value of L, while S E,K has the opposite behavior. The former behavior is similar to what has been observed in the presence of a magnetic field in a simpler setup [90]. Note that the curves at large a{Λ look qualitatively similar to what happens with the quark-antiquark potential. The most striking difference between the entanglement entropy and the quark-antiquark potential is that the entanglement entropy is almost independent of a for small enough a, while the quark-antiquark potential is very sensitive to even small values of a. We can explain this difference with the observation that at small a, the combination 2A S`W appearing in (5.9) either has a minimum, or is close to having one (i.e. 2A S`W is almost flat). Because of this the result is very sensitive to the boundary conditions string and tiny details of the geometry, leading to the observed strong dependence on a{Λ even for a{Λ ! 1. Moreover, we note that the potentials have nontrivial structure up to relatively large value of L " 10{Λ for tiny values of a{Λ, demonstrating that the string indeed probes the deep IR geometry.
The combination 4A`W corresponding to 2A S`W for the entanglement entropies in (5.11) and (5.13) though, is much steeper. We notice that the nontrivial behavior for the entanglement entropy takes place at smaller values of LΛ as compared to the quark-antiquark potential. For the entanglement entropy the characteristic scale is L " 1{Λ, which is to be expected in the absence of nontrivial dynamics and other scales. Therefore the characteristic RT surface remains relatively close to the boundary, r À 1{Λ, and is insensitive to the modified IR geometry due to small amounts of anisotropy (see Sec. 3.2.2 for the discussion on the geometry).
At first glance, one would suspect that this computation suffers from the same potential caveat as the quark-antiquark potential, namely the possibility that the string is unstable towards twisting in the bulk. However, since in this case the 'string' is just a slice of the Ryu-Takayanagi surface, such an instability cannot occur, because it would always amount to an RT surface diffeomorphism.

Conclusions and outlook
In this article we carried out a detailed analysis of anisotropic effects in improved holographic QCD including a fully backreacted quark sector in the Veneziano limit (V-QCD). The anisotropy was sourced by an axion profile with a linear dependence on one of the spatial coordinates, which is equivalent to deforming the theory by a space-dependent θ term. In turn, the action and the solution still remain translationally invariant along the relevant direction. We found numerical solutions of the equations of motion representing a family of anisotropic black branes and a thermal gas solution that is obtained from the black branes by sending the horizon to the deep interior of the geometry.
The latter thermal gas solution corresponds to the vacuum (zero temperature and entropy), heated up to a temperature T , of the corresponding anisotropic gauge theory at strong coupling. First we studied the RG flow in this vacuum state and showed that the IR end point is a scale invariant (almost) fixed point with approximate conformal symmetry SOp2, 3qˆR.
In the bulk picture this corresponds to a deep interior geometry AdS 4ˆR up to logarithmic corrections. This behavior is very similar to the IR theory obtained by deforming N " 4 super Yang-Mills with a magnetic field at strong coupling [91,92] which results in AdS 3ˆR 2 in the bulk. In that case the understanding in the field theory was that the fermions could effectively move only along the B direction at very low energies, due to Landau quantization in the transverse directions, and this gives rise to a CFT 2 which explains AdS 3 in the bulk. In the present case, we may try a similar understanding that directly follows from the picture 17 of the partition function (1.9) and (1.10). Modulo coupling to the non-Abelian gauge fields, precisely this field theory is used to study the Weyl semimetals at weak coupling [93][94][95]. On the other hand, it is quite reasonable to expect that integrating out the gauge fields produces an effective mass term Mqq for the Weyl quarks, just like in the NJL models. For values of M ą a{c a the Weyl semimetals are in the insulator phase as momentum in the x 3 -direction becomes gapped [93][94][95]. This then would restrict the motion of the quarks to the transverse plane and in the far IR would give rise to an approximate CFT 3 which would explain the AdS 4 factor. This interesting IR fixed point determines many of the interesting aspects we observe in this paper. We would like to remark that this IR geometry is present even for an infinitesimally small a and is drastically different than the a " 0 isotropic case. This drastic change was found to be due to a subtle competition between the potentials V g (which determines the behavior of the pure glue geometry) and Z (which controls the effect of the axion), and leads to a number of geometrical imprints in the IR.
Next we studied the thermodynamics and the phase diagram at finite temperature and found a rich structure with competition between confined/deconfined as well as chirally broken/symmetric phases as a function of T and a. This is summarized in figure 4. In particular 17 It should be noted however that our holographic model does not precisely correspond to this field theory. This is clear from the fact that the well-defined external symmetry (1.8) should correspond to a bulk dual written in terms ofÃ5`dχ whereÃ5 is a bulk gauge field dual to the axial U(1). We indeed have room for this in the V-QCD theory but decided not to include the dynamics ofÃ5 for simplicity.
we found that anisotropy can both deconfine the theory and destroy the condensate depending on the choice of x " N f {N c . Finally, we calculated certain observables of interest, such as the chiral condensate, particle spectra, quark-antiquark potentials and entanglement entropy which help investigate the effect of anisotropy both in the vacuum and the plasma states.
Perhaps our most important result concerns the behavior of the chiral condensate as a function of anisotropy. In [45], it was conjectured that the backreaction of a magnetic field onto the geometry was responsible for inverse magnetic catalysis, while the direct coupling was responsible for magnetic catalysis. In the present work we consider an electric neutral plasma, yet observe a very similar pattern of (inverse) magnetic catalysis generated by a instead of B. In particular both in the behavior of the chiral transition temperature, as well as the chiral condensate, the similarity is striking. We called the destruction of the chiral condensate by anisotropy, that we observed at strong coupling, the "inverse anisotropic catalysis". These observations lead us to the conclusion that the usual inverse magnetic catalysis that is observed in the presence of a magnetic field may equally well be caused by the anisotropy brought in the system by it. This is remarkable because, if we extrapolate it to QCD, it would mean that inverse magnetic catalysis could be recast entirely in terms of "inverse anisotropic catalysis". It would be extremely interesting to directly check this proposal on the lattice; see more on this below.
Another interesting result in the same context was the identification of a universal order parameter for (inverse) anisotropic catalysis. In order to do so we defined the quantity "anisotropic susceptibility" in analogy to the magnetic susceptibility. This susceptibility was then identified as a natural order parameter, similar to what was done in [45,47] for the magnetic case. We then obtained a relation between the critical temperature of the first order transition and the anisotropic susceptibility in Eq. (4.5).
One (possibly alarming) property of the IR geometry was the divergence of the string frame Ricci scalar, as seen in Fig. 2. The divergence is, however, extremely weak: it is due to corrections " log log r in the asymptotic IR geometry, and therefore practically absent in any anisotropic finite temperature solution. At zero temperature for any positive a, though, the deep IR is expected to receive stringy corrections which are absent at vanishing a. This indicates that the limits a Ñ 0 and α 1 Ñ 0 do not commute. Another related effect is the IR divergence of the anisotropic susceptibility for the anisotropic ground state, discussed in Sec. 4.1.
We also identified several physical consequences of the novel geometry in the IR: first, is the fact that the thermal gas solutions, which represent the confined isotropic phase, is replaced by a branch of small black holes. This is similar to what happens generally in the canonical ensemble of charged black holes, where the anisotropy parameter a plays the role of the charge. This has been seen explicitly in the model used in the current paper [67], where the resulting IR geometry was found to be AdS 2ˆR 3 , as well as in other works [69][70][71]. As a result, we find that linear confinement persists but only up to a certain length scale. This is evident from the behavior of the quark-antiquark potential, which shows signs of instability at large enough separations. Related to this, we observe that mesons and glueballs indeed fail to be absolutely stable but instead develop narrow widths, indicating the possibility of decaying to the AdS 4 vacuum.
There are some open questions and extensions of our work that are worth exploring: • Interplay with magnetic field. A natural future extension of the present work would be to introduce a magnetic field on top of the anisotropic background and to study the interplay between magnetic and anisotropic effects. This extension has a rich parameter space as there is the possibility of introducing an angle between the magnetic field and the direction singled out by the anisotropy.
• Potentials and universality. Here we used the potentials (in particular V g and Z) motivated by earlier work, which produce physics qualitatively similar to QCD. There is however still some freedom in choosing these potentials even taking into account all known constraints. In particular, the function Z could be modified by logarithmic corrections in λ in the IR, and the effects due to this could be studied. It could also be interesting to make completely different choices for the potentials and study the universality of our results. For instance, a Z which grows slower in the IR (e.g. Z " const.): in this case the drastic change in the IR structure between the a " 0 and a ą 0 solutions would be absent. Studying this could therefore be interesting even though this choice seems less motivated by the comparison to QCD.
• String embeddings. In the computation of the quark-antiquark potential, there is the possibility that the solutions for the string worldsheet that were found are unstable towards twisting in the bulk. Future work could study this possible instability. A more difficult problem to tackle is the physical significance of the unstable branches of the quark-antiquark potential. To investigate this, one would have to study time-dependent solutions of the string worldsheet. This is an inherently difficult problem.
• Thermalization and isotropization. Another interesting future project would be to study the time-dependence of the model more in general, and to address questions such as thermalization and isotropization of the QGP in the presence of fully backreacted quarks and anisotropic deformations (e.g. a space-dependent θ term and/or a magnetic field). As a first step, one could compute the full quasinormal mode spectrum of the system and transport coefficients. Some steps in this direction were recently taken in [96].
• Field theory and lattice QCD. We initiated a study of the problem directly in field theory in the Introduction, see e.g. equations (1.9) and (1.10). It is tempting to carry out the calculation of the quark condensate perturbatively using this setting. However, in analogy with an external magnetic field [29] we do not expect the Inverse Anisotropic Catalysis phenomenon to be present at weak coupling. Nevertheless this calculation may give us hints towards a better understanding of the phenomenon. Another idea is to study the problem in an effective field theory such as the NJL, which in fact provided important insights for the IMC [40]. Finally, it would be very interesting to directly check our proposal on the lattice. Anisotropy on a lattice can be introduced by assigning a different number of lattice points in one spatial direction than the directions transverse to it. This setting will simply correspond to the discretization of the system in (1.4). This is of course not easily done as said, because introducing fermions on an anisotropic lattice turns out to be a challenging problem, see e.g. [97]. It may be easier to consider the picture given by (1.10) instead. Since this is essentially a system of Weyl semimetal with the Dirac cone separated in the left and the right parts in the momentum space by an axial gauge field, and since it is possible to study Weyl semimetals on the lattice, we expect this picture to be more suitable for the lattice studies. Yet, one has to check directly if this action is plagued by a fermion sign problem. Finally, we note that the effect of anisotropy on the critical temperature can already be checked in a pure glue setting as in [1], which, on the lattice, is computationally easier than a model which includes fermions.
We hope to come back to these points in the near future.
We further insert the rough approximation e´A » r{q where q is constant. We obtain The solution in the simplest approximation is therefore For the IR asymptotics of section 3.2, q 2 V g " const. Moreover we have chosen the potentials such that V g κ " const at large λ, we find that q 2 {pκpλqq " const and the tachyon therefore obeys a power law in the IR. This is enough for the tachyon to decouple and for the assumptions we made above to be valid. Notice however that it may be enough to modify the subleading logarithmic corrections to the IR asymptotics of, say, κ to change this conclusion.

B Holographic renormalization
In this appendix we write down the counterterms needed to regularize the free energy and the anisotropic susceptibility χ a . In the case of flat boundary metric, the required counterterms are [78] S ct "´M 3 N 2 where γ is the boundary metric and we also used the fact that the only field depending on the spatial coordinates is χ. The various functions are defined as follows. First, U pλq is the superpotential. Up to a choice of scheme, we can choose any solution of the superpotential equation (Eq. (1.6) in [54]). In this work however we only use explicit counterterms to renormalize χ a (while the free energy is obtained by integrating the first law of thermodynamics) and its renormalization is independent of U pλq. The function Θpλq can be found in general as an integral defined in terms of the superpotential and the function Zpλq [78]. The function cpλq cancels a remaining logarithmic divergence and only its expansion in the UV is needed. For our purposes it is sufficient to note that the following counterterms for the susceptibility are equivalent to the general prescription: where the integral arises from Θpλq,Ã andλ give the solution at a " 0 and T " 0, we introduced a UV cutoff , and For the numerical evaluation of χ a we chose r 0 such thatÃpr 0 q " 0 in units where " 1.

C Fluctuation equations
We give here the fluctuation equations for two additional modes apart from the helicity two glueballs discussed in the main text. We have checked numerically that the spectrum of these sectors behaves qualitatively similarly as the helicity two glueballs.

C.1 Flavor nonsinglet mesons
The flavor nonsinglet mesons (states in the adjoint representation of the unbroken SU pN f q V ) are decoupled from the glueballs. The spatial asymmetry is mediated to the meson sector (only) by the metric, and therefore the changes with respect to the isotropic case [59,98] are minor. In particular, as it turns out, the asymmetry does not lead to the mixing of any of the helicity zero states in our setup. For example, the fluctuation equation for the rho mesons is 1 V f pλ, τ q wpλ, τ q 2 e A˘W G B r`Vf pλ, τ q wpλ, τ q 2 e A˘W G´1 B r ψ V˘``ω 2´q2 e´2 W˘ψ V " 0 (C.1) with plus (minus) signs in the exponents for helicity one (zero) states. Here ψ V prq is the radial wave function for the fluctuations.

C.2 Helicity one glueballs
The helicity one states also turn out to be simple. For B " 0 there are no background gauge fields, and since the action is quadratic in gauge fields, gauge field fluctuations decouple from the metric. The metric fluctuations, in the gauge where δg µr " 0, are δg it and δg i3 with i " 1, 2. Constraints arising from the ir components of the Einstein equations eliminate two of these. Moreover the states with positive and negative helicities decouple. Therefore the remaining two physical fluctuations are decoupled and satisfy the same equations. Defining (considering the positive helicity for example) δg 13`i δg 23 " e 2Aprq e´i ωt`iqx 3 h 3 prq , δg 1t`i δg 2t " e 2Aprq e´i ωt`iqx 3 h t prq , (C.2) and choosing the invariant combination ζprq " ωh 3 prq`qh t prq, the fluctuation equation reads ζ 2 prq`˜3A 1 prq´ω 2`q2 e´2 W prq ω 2´q2 e´2 W prq W 1 prq¸ζ 1 prq`´ω 2´q2 e´2 W prq¯ζ prq " 0 .

(C.3)
At q " 0 we see that the only change with respect to the equation for the helicity two glueballs is the opposite sign in the term 9 W 1 prq.