Phase diagram of 4D field theories with chiral anomaly from holography

Within gauge/gravity duality, we study the class of four dimensional CFTs with chiral anomaly described by Einstein-Maxwell-Chern-Simons theory in five dimensions. In particular we determine the phase diagram at finite temperature, chemical potential and magnetic field. At high temperatures the solution is given by an electrically and magnetically charged AdS Reissner-Nordstroem black brane. For sufficiently large Chern-Simons coupling and at sufficiently low temperatures and small magnetic fields, we find a new phase with helical order, breaking translational invariance spontaneously. For the Chern-Simons couplings studied, the phase transition is second order with mean field exponents. Since the entropy density vanishes in the limit of zero temperature we are confident that this is the true ground state which is the holographic version of a chiral magnetic spiral.


Introduction
One of the amazing developments emerging from the research in string theory is the idea of a gauge/gravity duality [1]. Remarkably, the duality relates the strongly coupled regime of gauge theories to the weakly coupled regime of the dual string theory or (super-)gravity. Consequently, it has become a powerful tool to study strongly interacting systems by using a conjectured dual weakly coupled gravitational theory. At present, holographic descriptions of non-perturbative phenomena include, among other applications, condensed matter physics, high energy physics and quark-gluon plasma. 1 Many of these real world systems of interest involve a finite chemical potential and strongly-coupled degrees of freedom. However, only a few reliable methods exist to compute physical observables for these systems, with real-time physics being particularly difficult to study.

JHEP03(2016)164
Using the framework of gauge/gravity duality it is possible to study strongly coupled conformal field theories at finite temperature and finite charge density. For example, the simplest bottom-up holographic model of strongly-interacting matter is Einstein-Maxwell theory with negative cosmological constant. The dual field theory is some CFT with a global U(1) symmetry. The asymptotically AdS Reissner-Nordstroem black brane solution describes thermal equilibrium states with a finite U(1) charge density.
If the field theory spacetime is even, and the global current is anomalous, the dual gravitational theory also contains a Chern-Simons term for the U(1) gauge field. In a series of papers [9][10][11][12][13][14], D'Hoker and Kraus constructed the electrically and magnetically charged black brane solutions in five dimensional Einstein-Maxwell-Chern-Simons theory which are dual to strongly coupled four dimensional conformal field theories with chiral anomaly at finite temperature, chemical potential and magnetic field. The phase diagram of these field theories exhibits interesting features, such as a quantum critical point.
The instability of the asymptotically Reissner-Nordstroem black brane solution has attracted much attention due to its relevance to quantum phase transitions in the dual strongly interacting quantum field theory at finite density. For example, in the presence of (charged) scalar fields or non-abelian vector fields new phases were studied which are reminiscent of s-wave and p-wave superfluids [15][16][17][18]. Moreover, for large enough chiral anomaly, and for low temperatures new spatially modulated phases were found [19][20][21] at zero magnetic field.
In this paper, we consider the class of strongly coupled four dimensional CFTs with chiral anomaly whose gravitational dual description is given in terms of Einstein-Maxwell-Chern-Simons theory. 2 We study thermal equilibrium states for finite temperature, chemical potential and magnetic field and determine the phase diagram. Depending on the coefficient of the chiral anomaly we find a spatially modulated phase 3 extending the results of [19] to non-zero magnetic fields. In particular, the quantum critical point [11,12] is hidden within this new phase.
The new spatially modulated phase discussed in this paper may be viewed as a holographic version of a chiral spiral [35], 4 although technically speaking we only have one anomalous current in contrast 5 to QCD. Moreover, the interplay between the quantum critical point and spatially modulated phases is also observed in certain meta-magnetic materials in condensed matter physics. In particular these materials may have a quantum critical point due to the meta-magnetic phase transition which is hidden behind a nematic phase (e.g. see [38]).
The remainder of the paper is organised as follows. In section 2 we summarize the holographic setup used here. First, we present our coordinate ansatz in the gravitational theory which exhibits Bianchi VII 0 symmetry implying that the corresponding equations 2 It would be interesting to add scalar fields along the lines of [22] and investigate the interplay between the quantum critical point as well as spatially modulated and s-wave superfluid phases. 3 Spatially modulated phases in the presence of magnetic fields were also discussed in [23][24][25][26][27][28][29][30][31][32][33][34]. 4 For other work on holographic chiral spirals see [36,37]. 5 Hence, if we compare to QCD, the anomalous current may be identified with the axial current. Moreover, µ should be viewed as an axial chemical potential, and B as an axial magnetic field. of motion are ordinary differential equations. Then, we briefly state the asymptotic expansions close to the horizon and conformal boundary and discuss how to extract the thermodynamic observables of the dual CFT from the gravitational theory.
In section 3 we numerically construct asymptotically AdS 5 black brane solutions with non-trivial electric charge density and magnetic field breaking translational invariance spontaneously. In particular, we determine the phase diagram at finite temperature, chemical potential and magnetic field. Moreover, we characterise the new phase by identifying the order parameters and critical exponents close the phase transition. Finally we explicitly show that the entropy density vanishes in the limit of zero temperature.
In section 4 we summarise the results of the paper and conclude with a few remarks on possible future directions. More details concerning the equations of motion, thermodynamics, special cases and numerics are given in the appendices to the paper. Note that in appendix D some techniques are presented to improve the numerical accuracy, specifically at low temperatures, which may be relevant also for other holographic setups.

Holographic setup
To describe a strongly coupled four-dimensional field theory with chiral anomaly within the framework of gauge/gravity duality, we consider a simple toy model on the gravity side. Let us collect the minimal features of this toy model. First, it has to contain a metric (with components g mn ), and a U(1) gauge field A = A m dx m . Second, the spacetime should be asymptotically AdS 5 , i.e. using coordinates (x m ) = (x µ , z), where x µ may be identified with the field theory coordinates 6 and z is the radial coordinate of AdS 5 . The line-element ds 2 reads ds 2 = L 2 z 2 dz 2 + η µν dx µ dx ν , (2.1) in the limit z → 0. Here, L is the radius of AdS 5 . The metric g mn is dual to the energymomentum tensor T cf t µν of the corresponding four-dimensional CFT, while the gauge field A = A m dx m is dual to the current J cf t µ on the field theory side. Due to the chiral anomaly of the CFT the current J cf t µ is not conserved, i.e.
where µνρσ denotes the totally antisymmetric tensor of four dimensional Minkowski spacetime which we normalise such that 0123 = 1. Moreover,F is the field strength tensor of an external fieldÃ which may be viewed as a source term for the current J cf t . The chiral anomaly with its coefficient γ is dual to a Chern-Simons term on the gravity side, which is another vital ingredient of the gravitational toy model. To summarize, the action of the gravitational toy model (with the features highlighted above) reads 3) 6 Moreover, x 0 = t and x i are the spatial coordinates of the field theory.

JHEP03(2016)164
where g = det g mn . Moreover, F = dA is the field strength tensor of the U(1) gauge field 7 A and 2κ 2 ≡ 16πG 5 with the five-dimensional gravitational constant G 5 . For γ = 2/ √ 3, the action (2.3) coincides with the bosonic part of minimal gauged supergravity in five dimensions and hence is a consistent truncation of the most general class of type IIB supergravity in ten dimensions or supergravity in eleven dimensions which are dual to N = 1 superconformal field theories, see e.g. [39][40][41][42]. However, in this paper γ is treated as a free parameter and we study the phase diagram as a function of γ. This action has to be supplemented by boundary terms [43][44][45] of the form Here, h µν is the metric induced by g mn on the conformal boundary of AdS 5 . The extrinsic curvature K mn is given by where ∇ is the covariant derivative and n m are the components of the outward pointing normal vector of the boundary ∂M. Moreover, K is the trace of the extrinsic curvature with respect to the metric at the boundary. The first term in (2.4) is the standard Gibbons-Hawking term which is needed for the well-posedness of the variational principle. The other terms remove divergencies and hence are required for the proper renormalisation of various physical quantities [43][44][45]. In our case, the conformal boundary of AdS 5 is flat Minkowski space and hence the Ricci scalar associated with the boundary metric, R(h), vanishes. Note that the last term in (2.4) is not invariant under diffeomorphisms. As we will see explicitly, this term in the boundary action is needed to remove the divergence associated with the trace anomaly on the field theory side. To keep notations compact we set 2κ 2 = 1 as well as L = 1 from now on. The equations of motion associated with the action (2.3) read for the metric, as well as for the gauge field. Equivalently, we can rewrite (2.8) as Note that the two gauge fields A andÃ are closely related but are not identical. In particular, the gauge field A lives in the five-dimensional curved spacetime, while the external fieldÃ is dual to the current J cf t and hence is defined on the four-dimensional Minkowski spacetime of the dual field theory.

JHEP03(2016)164
where˜ mnopq is the totally antisymmetric Levi-Civita symbol in five spacetime dimensions with˜ t123z = 1. Moreover, the field strength tensor has to satisfy the Bianchi identitiy dF = 0. For example, a solution to the equations of motion (2.7) and (2.8) is given by the electrically charged, asymptotically AdS 5 Reissner-Nordstroem (RN) black brane with metric and gauge field given by Herefrom we get the field strength tensor The function u(z) and E(z) are given by leading to the electric field e(z) = ρ z. (2.13) The parameter ρ = 2 µ is related to the density of the dual field theory, i.e. J cf t t = −ρ. In a series of papers [9][10][11][12][13] (see [14] for a review) the electrically charged RN black brane in asymptotically AdS 5 spacetime was generalized to allow for a constant non-vanishing magnetic field B. 8 Without loss of generality, we can assume that the constant magnetic field B is aligned in x 3 direction. As reviewed in the introduction, and as explicitly reproduced in appendix C.2, the solution exhibits interesting features, such as a quantum critical point.
Here, we study particular instabilities against spatial modulation for the electrically and magnetically charged Reissner-Nordstroem black brane. As shown in [20,21] for γ > γ c ≈ 1.158, the electrically charged AdS-RN black brane is unstable against spatial modulation below a critical temperature, suggesting that the system is in a spatially modulated phase in which the current acquires a helical order. The corresponding backreacted solution for zero magnetic field was presented in [19]. In this paper we find numerical evidence that the helical structure at zero magnetic field persists at finite magnetic field, at least in some part of the phase diagram. In particular we construct electrically and magnetically charged black branes with (reduced) Bianchi VII 0 symmetry, which give rise to the helical order in the currents and energy-momentum tensor.
The Bianchi VII 0 symmetry is manifest using the one-forms ω i defined by 14) 8 We sometimes refer to this solution as the charged magnetic brane solution. Note that ω 1 ∧ ω 2 = dx 1 ∧ dx 2 as well as dω 1 = k ω 2 ∧ ω 3 , dω 2 = −k ω 1 ∧ ω 3 and dω 3 = 0. The meaning of the differential forms ω i , or to be precise their dual tangent vectors, is apparent from figure 1. Using the differential forms (2.14) we assume 9 that the helical structure is parallel to the magnetic field B, i.e. ω 3 is aligned along the magnetic field. ω 1 and ω 2 span the plane of the spatial directions perpendicular to the magnetic fields, i.e. the (x 1 ,x 2 )-plane. In this paper, we determine the phase diagram at finite magnetic field, chemical potential and temperature. In particular, we use the following ansatz for the metric and for the field strength tensor F = dA. Note that B has to be independent of z in order to satisfy the Bianchi identitiy dF = 0. The ansatz specified by (2.15) and (2.16) respects the Bianchi VII 0 symmetry mentioned above. The field strength tensor may be obtained from a gauge field A of the form where A t (z) is related to the electric field e(z) by E (z) = e(z). Moreover, P (z) = p(z), where denotes the derivative with respect to z. Note that the ansatz (2.15), and (2.16) (or (2.17) repsectively) generalises [11] and [19]: for α(z) = 1, b(z) = Q(z) = 0 as well as g(z) = 0 and k = 0 we obtain the original D'Hoker

JHEP03(2016)164
and Kraus background [11], while in the limiting case c(z) = 0 as well as B = 0 and g(z) = p(z) = 0 we obtain a setup equivalent to Donos and Gauntlett [19] as discussed in appendix C.1.
Inserting the ansatz (2.15) and (2.16) into the Einstein equations (2.7), we obtain nine differential equations, i.e. seven second order differential equations for the metric functions u(z), v(z), w(z), q(z), c(z), α(z) and g(z) as well as two constraints which we denote by CON 1 and CON 2 and contain only first derivatives of the metric fields. Moreover, there are three independent equations of motion (2.8) for the gauge fields. While b(z) has to satisfy a second order differential equation, the fields p(z) and e(z) satisfy first order equations of motion which can be recasted in the form The full form of the remaining equations of motion is not very enlightening. For these reasons we do not display them here (see appendix A for more details). Moreover, we explicitly checked that the two constraints are consistent. To be precise, using the equations of motion of the metric and gauge fields we showed that the constraints satisfy the following differential equations ∂ z (CON 1 (z)) + f (z) CON 1 (z) = 0 ∂ z (CON 2 (z)) +f (z) CON 2 (z) = 0 (2.20) for some functions f (z) andf (z). These differential equations for CON i (z) can be formally solved in terms of exponential functions. Hence, solving the constraints CON i (z 0 ) = 0 for some z = z 0 guarantees that the constraints are satisfied for all z.

Asymptotic expansions
In order to solve the equations of motion, we first consider the asymptotic expansion of the metric and gauge fields close to the horizon and the conformal boundary of the spacetime. For instance, imposing asymptotically AdS (see discussion in appendix A) we obtain for the metric functions while the gauge field functions read Using diffeomorphisms we can shift the horizon to z = 1. The event horizon condition imposes that u(1) = 0. Then, it follows from the regularity conditions that c(1) = q(1) = 0. Therefore, the expansion around z = 1 assumes the following structure for the metric functions Furthermore, we also impose regulartiy for A t at the horizon, i.e. E(1) = 0. In turn, we obtain for the gauge field functions The boldface letters in (2.21)-(2.24) denote quantities that are not determined by the expansion, i.e., their values are obtained only after a global solution is found. In the ansatz for the gauge field (2.17) the functions E(z) and P (z) appear. By integrating (2.18) we can determine P (z) and E(z). For γB = 0, the function E(z) is given by

JHEP03(2016)164
Since we identify E(0) with the chemical potential µ, i.e. A t (0) = −E(0) = µ, we can use eqs. (2.19) and (2.25) as well as the asymptotic expansion (2.22) to obtain Similarly, we can solve the equation (2.18) to determine P (z). Note that in this case, we only have to demand regularity at the horizon and hence we cannot fix P (1). However, we only want to study systems where we do not allow for a source term of the operator dual to P (z) and hence we have to demand P (0) = 0.

Thermodynamics
Next we describe how to extract thermodynamic information from our solutions which describe thermal equilibrium states in the dual CFT. We will work in the grand canonical ensemble in which the chemical potential µ is fixed. In order to analyse the thermodynamical properties of the black brane solution we have to analytically continue to a Euclidean time τ = t (E) by a Wick rotation of the form τ = i t.
Since the metric and the vector field should be real in Euclidean signature, we also have to introduce q (E) = −i q, c (E) = −i c as well as e (E) = −i e. The leading order coefficients for q (E) , c (E) and e (E) at the horizon, see (2.23) and (2.24), are denoted byq 1(E) ,c 1(E) and e 0(E) respectively. Hence the Euclidean metric and the field strength tensor of the gauge field near the horizon z = 1 are given by as well as In the last lines of (2.27) and of (2.28) we kept only the leading terms in the near horizon limit. The temperature of the black brane solution can be deduced by demanding regularity of the Euclidean metric (2.27). In particular, we find that the temperature is given by The entropy S is given by the Bekenstein-Hawking entropy of the black brane. Due to the infinite area of the event horizon it is more convenient to work with the entropy density s.
Since we work in units where 2κ 2 = 16πG 5 = 1, the entropy density is given by

JHEP03(2016)164
This entropy density can be also deduced from the grand canonical potential Ω. In AdS/CFT, the grand canonical potential is identified with T times the on-shell bulk action in Euclidean signature. We thus analytically continue to Euclidean signature and compactify the Euclidean time direction τ with period 1/T . In order to determine the total Euclidean action S (E)tot , we first perform the Wick-rotation on the action S as given by (2.3) (including its boundary terms (2.4)), and denote the result byS grav andS bdy respectively. Then the corresponding Euclidean actions S (E)grav and S (E)bdy are given by S (E)grav = −iS grav and S (E)bdy = −iS bdy . Note that S (E)tot = −S tot . We will use the action S tot in Minkowski signature from now on. Hence the grand-canonical potential is given by where o.s. indicates that we have to evaluate the total action on-shell. The total Euclidean on-shell action is displayed in appendix B. Using the boundary and horizon expansions (2.21)-(2.24), we obtain for the density of the grand canonical potential (which we also denote by Ω to keep the notation simple) Due to the standard recipes of AdS/CFT, in particular the formula 10 (see [44]) we can extract the energy-momentum tensor of the dual conformal field theory. The nonvanishing components of the energy-momentum tensor are given by 10 From now on, we drop the superscript cft of the energy momentum tensor T cf t µν and of the current J cf t µ to keep the notation simple. Moreover, recall that we set 2κ 2 ≡ 1 from the beginning.

JHEP03(2016)164
In particular, the trace of the energy momentum tensor is given by Similarly, we can also read off the expectation value of the current in the dual field theory using the relation [10] For our ansatz, the non-zero components of the current are given by Hence we can write Ω in the form where we have assumed that u (0) < 0 and henceū 1 > 0 which is the case for our numerical results. Moreover, the charge density reads J t = ρ, while the energy density is given by Note that for γ = 0, the grand canonical potential (2.39) reduces to its standard form Ω = U −sT −µ J t . If both γ and B are non-vanishing we obtain an additional contribution to the grand canonical potential due to the chiral anomaly. Moreover, combining (2.38) and (2.26), we obtain a relation between J x 3 and µ of the form This is precisely the chiral magnetic effect. Due to the chiral anomaly we also expect to find [46]

The magnetic helical black brane
We present now the results describing our magnetic helical black brane solution. Appendix D provides more details on the numerical techniques employed here. Due to the invariance of the metric (2.15) under scale transformationx m = λ x m , we express the results in terms of dimensionless quantities normalised by µ. Since under the scale transformation we haveμ = λ µ, the relevant physical observables arē As described in [19], forB = 0 one expects to construct the spatially modulated black brane solutions provided the Chern-Simons coupling be γ >  we focus ourselves on the results with γ = 1.5. Besides, as a generalisation of the particular results from [19], we also comment on some specific features of the case γ = 1.7.
We first address the question in which region of the parameter space {k,B} we expect new solutions. Then we construct these solutions and single out the thermodynamically preferred ones. Following this we discuss thermodynamic properties of these solutions, with particular emphasis on the behaviour of near the critical temperature and in the low temperature limit.

The phase boundary
The magnetic helical black brane solution is described by the existence of a function b(z) = 0. As described in appendix C.2, in the limiting case b(z) = 0, one obtains the electrically charged Reissner-Nordstroem black brane withB = 0 or the electrically and magnetically charged brane forB = 0. The boundary between the two regimes is therefore naturally defined as the region for which b(z) ≈ 0. 11 In the left panel of figure 2, the highest temperatureT for which the magnetic helical black brane exists is plotted as a function ofk andB. In other words, above the surface only the RN black brane exists while below the surface both the RN black brane as well as the magnetic helical black brane solution coexist. Projecting the surface onto theB = 0 plane we reproduce 12 the expected profile showed in [19] and the contour plot in the right panel of the same figure highlights the isothermal curves on the (k,B)-plane.
To appreciate this property, in figure 3 we restrict ourselves to the (T ,k)-plane with B = constant. The left panel corresponds to the same results as in the previous figure, i.e. with γ = 1.5. It becomes evident that the critical temperatureT It is worth mentioning the equivalent results for γ = 1.7, depicted in the right panel of figure 3. For valuesB 0.274, we observe that the magnetic helical phase lies entirely 11 Numerically, the boundary is characterised by b2 ≈ 10 −9 as introduced in (2.21), see discussion in appendix D. 12 The results shown in figure 2 are for γ = 1.5 while the results explicitly shown in [19] are for γ = 1.7. within a closed curve. In the particular example withB = 0.275 displayed here, phase transitions occur at bothT C ≈ 0.03909 andT C ≈ 0.022061.
After identifying the critical temperature in the (T ,k)-plane, we study the dependence ofT C on the magnetic fieldB and present the results in figure 4. Here again, it is evident that there exists a valueB 0 asT C → 0, which limits the region where the magnetic helical solution is expected to be found. For the particular examples treated here, these values areB 0 ≈ 0.279 (γ = 1.5) andB 0 ≈ 0.274 (γ = 1.7). We also identify in the same figure the quantum critical pointB C as found in [11] (see also the discussion in appendix C.2). In particular, the critical values areB C ≈ 0.185 andB C ≈ 0.220 for γ = 1.7 and γ = 1.5, respectively. It is interesting to notice thatB C lies within the new phase region, meaning that the phase transition should occur before the system reaches the quantum critical point.
The important question that arises now is what happens to the system as we lower the temperature and move inside the new phase along curves of constantB. Of particular interests is the regionB <B C and the behaviour of the entropys in the low temperature regime. We address this issue and discuss further details about the thermodynamics in the next section.

Thermodynamic results
For fixedB and for fixed temperatureT we construct the solutions for different values ofk. The solution corresponding to the physical state minimizes the grand canonical potential Ω(k). The corresponding value fork minimising the grand canonical potential is denoted byk * . For fixedB we repeat this procedure for smaller temperaturesT and hence obtain a trajectoryk * (T ) in the (T ,k)-plane of thermodynamically preferred solutions. This trajectory is shown in figure 5 for the valuesB = 0.200 <B C andB = 0.250 >B C . In both cases, note that when lowering the temperatureT , the wave-numberk * (T ) decreases and hence the pitchp * = 2π/k * increases. Along such trajectories of thermodynamically preferred solutions, we evaluate the observables derived in section 2.2 and compare them to the corresponding values from the charged magnetic solution. In all the following figures, a continuous line represents a result within the new magnetic helical phase, whereas the dashed lines depict the results of [11] (see appendix C.2 for more details). First, in figure 6 we present the entropy densitys as a function ofT . Let us first concentrate on the dahsed lines corresponding to the charged magnetic black brane. ForB <B C the entropy goes to a non-vanishing constant asT → 0 in agreement with the results of [11]. However, for the new helical magnetic black brane construct in this paper, we observe thats → 0 asT → 0 regardless of the value ofB. Due to the vanishing entropy density we are confident that the magnetic helical black brane is dual to the true ground state of the CFT. Moreover, for fixedB the entropy is continuous close to the phase transition, i.e. forT T C (B). Hence the phase transition is second order.
Next, we turn to the non-vanishing components of the energy-momentum tensor T µν and the current J µ of the dual field theory. Figure 7 depicts the components T tt , In all cases there are expected small deviation between the helical magnetic black brane and the charge magnetic black brane. Moreover, from eq. (2.40) and the normalisation (3.1) it is clear that J x 3 = −γB is a constant. Furthermore, we also confirm that (2.41) holds numerically, i.e. in terms of dimensionless quantities T tx 3 = 1 2 γB. Note that the relation (2.41) is also satisfied for the charged magnetic black brane [11] which we explicitly demonstrate in appendix C.2. Finally, we display the results for the components T tω 1 , T ω 1 ω 1 − T ω 2 ω 2 , T ω 1 x 3 and J ω 1 . As one can see in the left panels of figure 8, these components vanish asT →T C .

Summary and outlook
We studied strongly coupled four dimensional CFTs with chiral anomaly whose gravitational description is given in terms of Einstein-Maxwell-Chern-Simons theory in asymptotically AdS spacetime. In particular, we investigated the phase diagram at finite temperature, chemical potential and magnetic field and found a new spatially modulated phase for low temperatures and small magnetic fields. This is dual to asymptotically AdS 5 black brane solutions with non-trivial electric charge density and magnetic field spontaneously acquiring a helical current which we construct numerically. Note that this helical phase is supposed to exist only for large enough coefficient γ of the chiral anomaly. In this paper we presented results mainly for γ = 1.5.
The new spatially modulated phase has interesting features. First, the quantum critical point at B = B C is hidden in this new phase, at least for the values of γ studied here. Second, the phase transition is second order with mean field exponents. Third, our numerical results indicate that the entropy density vanishes in the limit of zero temperature supporting our speculations that this is the true ground state of the system. Fourth, we have extracted how the wave-number of the helical structure changes with the parameters of the phase diagram. Finally, during the course of this work, as a side product we developed new numerical techniques which may be also useful for other holographic systems.

JHEP03(2016)164
It will be very interesting to further analyse the system by addressing questions such as: does the phase diagram change for other values of the chiral anomaly coefficient? Are the states with helical structure aligned to the magnetic field really thermodynamically favoured? Answering this question will require to solve partial differential equations on the gravity side.
Moreover, it will be worthwile to explore if there exists a simple relation between the location of the quantum critical point, given by B = B C , and the location of the phase boundary B = B 0 at zero temperature. Note that both, the new phase and the quantum critical point, are controlled by the chiral anomaly coefficient and hence such a relation may exist although it is not obvious in terms of the dual field theory.
Another future direction is to study transport coefficients and quasi-normal modes within the new phase extending the results of [47]. Finally, the results presented here can be generalised to models with an anomaly structure closer to the one of QCD and Weyl semimetals. In particular, the interplay between a magnetic field and chemical potentials for the vector/axial charge may be interesting.

Acknowledgments
We are very grateful to Marcus Ansorg, Amadeo Jimenez Alba and Sebastian Möckel for valuable discussions and for insightful comments on the draft. JL and RM acknowledge financial support by Deutsche Forschungsgemeinschaft (DFG) GRK 1523/2. RM was supported by CNPq under the programme "Ciência sem Fronteiras". MA would like to thank KITPC and the organizers of the workshop "Holographic duality for condensed matter systems", as well as University of Vienna and the organizers of the workshop "Vienna Central European Seminar" for their hospitality.

A Equations of motion
In this section we give some details on the equations of motions and, in particular, how we treat them as a boundary value problem. For convenience, let us reproduce here eqs. (2.7) and (2.9) and define them as Furthermore, we complete the one-forms (2.14) with ω 0 = dt and ω 4 = dz and express the equations of motion in terms of the tetrad basis ω a = ω a m dx m , i.e., we look specifically at 13 13 In the tetrad basis, the equations do not present any trigonometric term related to cos(k x3) or sin(k x3).

JHEP03(2016)164
The non-zeros components {E 00 , E 01 , E 03 , E 11 , E 13 , E 22 , E 24 , E 33 , E 44 } and {M 1234 , M 1245 , M 2345 } form a system of 12 ordinary differential equations (ODE) for our 10 field variables: the components of the metric (7 functions) and gauge field (3 functions). In spite of being overdetermined, this system of ODE is consistent as already shown in the main text. Therefore, we must solve 10 out of 12 equations and ensure that the remaining 2 are satisfied for at least one value of z. Yet, we must assure that the chosen equations are independent of each other. The first point to notice is that the second derivatives appearing in each of the Maxwell-Chern-Simons equations involve only one of each gauge field function. In other words, the equations M 1234 , M 1245 and M 2345 can be straightforward regarded individually as equations for E(z), P (z) and b(z), respectively. Next, we observe that E 24 contains only first order derivatives. Finally, one can work the second derivates out of the remaining E ab and obtain equations for each one of the metric fields u(z), c(z), v(z), w(z), α(z), g(z) or q(z). This procedure leaves us with 7 second order ODEs and one additional first order ODE (apart from E 24 ).
By sorting out the second derivatives, we can associate for each one of the fields its respective ODE. In this way, we need boundary values at both z = 0 and z = 1. The only exception is the function g(z), for which we work with the first order ODE E 24 and therefore we are only allowed to fix the value at one of the surfaces. To exemplify the structure of the system of equations, let us collect the field variables and the equations of motion into the vector notation The boundary values are a mixture of regularity conditions imposed by the equations of motion and physical assumptions insuring the surfaces z = 0 and z = 1 to represent the AdS boundary and the event horizon, respectively. For example, at z = 1 the horizon condition tells us that u(1) = 0. Moreover due to regularity, we have to impose E(1) = 0. By imposing such conditions on the remaining equations, we are left with regularity conditions involving the value of fields and their first derivatives at z = 1 (Robin boundary conditions). In some specific cases, the conditions are rather simple and reduce to q(1) = 0, c(1) = 0.
The next step is to study the asymptotic expansion around the AdS boundary z = 0. In its most generic form, the expansions read u(z) = 1 + u 2 z 2 + u 4 z 4 + O(z 6 ) + z 4 ln(z) û 4 + O(z 2 ) as well as The quantities in boldface are free parameters, which can not be determined by the series expansion. All the other terms are fixed by them if one considers all equations of motion (including here the two first order differential equations). For example, we find that Asymptotically AdS solutions require v 0 = 1, w 0 = 1, α 0 = 1, c 0 = 0, g 0 = 0, q 0 = 0 , and thus b 0 = 0 from (A.6). In our numerics we demand u 1 = 0 in order to fix all remaining diffeomorphisms. For the gauge field functions, we fix the chemical potential 14 E 0 = −µ. Moreover, we do not allow for source term for operators dual to P (z). Hence we impose P 0 = 0. With this conditions, the expansions around z = 0 assume the much simpler form given by (2.21). Once the solution is available, the thermodynamic observables (see appendix B) require the knowledge of some coefficients related to higher derivatives, such as u 4 , w 4 , a 4 , c 4 and 14 From the ODE point of view, we could also prescribe E2 = ρ/2 instead of E0.
-20 -JHEP03(2016)164 q 4 . Not only do we lose accuracy by calculating them numerically, but there are also some cases in which the derivative might not even exist due to the presence of terms z 4 ln(z). In order to get access to all the needed coefficients with a reliable high accuracy, we incorporate the boundary conditions into our variables and introduce auxiliary fields via

B More on thermodynamics
Calculating thermodynamic properties requires the evaluation of the on-shell action, which in principle is an integral over the numerically determined solutions. It is more enlightening to have an analytic expression in which only boundary values of the solution have to be inserted. In this appendix we show how to rewrite (part of) the on-shell Lagrangian as a total derivative. This is usually also used to get Smarr Type formulas, for example see [48,49]. We present here systematically how to do this step by step. First by taking the trace of the Einstein equations, R mn = −4g mn + T mn , R = −20 + T , we can reformulate the Einstein-Maxwell part of the Lagrangian as The equations of motion (B.1) may be written in the following form whereT n m = F mo F no /2. Hence we obtain for L EM In order to rewrite √ −g L EM as a total derivative we have to massage R t t andT t t .

JHEP03(2016)164
Let us start with R t t by using the identity R m n ξ n = ∇ m ∇ n ξ n for an arbitrary Killing vector ξ (with components ξ m ). If the metric depends only on the radial coordinate z the identity can be expressed as a total derivative In order to rewrite √ −gT t t as a total derivative we analyse the Maxwell equations in the following systematic way: the F ∧ F term has a rather simple structure while the term d * F is a sum of total derivatives of z and x 3 The relevant parts of * F are those, without either dz or dx 3 . This gives three independent equations of motion. By comparing to the remainderT t t given by we see that the Maxwell equation with legs dx 1 ∧ dx 2 ∧ dx 3 ∧ dz corresponds to equation relating p(z) and P * (z) in (2.18). For convenience, we re-write it here in the form Now, in order to rewrite −2 √ −gT t t as a total derivative, we multiply M 1d (z) by E(z) and substract again the additional term of ∂ z (E(z) M 1d (z)). 15 P * (z) was first introduced in eq. (2.19).

JHEP03(2016)164
Finally we also have to consider the Chern-Simons term plus the remaining term E(z) R 1 (z), 14) The absent terms in (B.14) are proportional to x 2 cos(kx 3 ) and vanish when we integrate over any symmetric interval with respect to x 2 . The expression (B.14) can be reformulated in terms of the following two total derivatives ∂ ∂z plus the remaining term 1 3 B γ E(z) p(z). Collecting everything we end up with the action density s, defined by S = vol(R 3,1 ) s, Finally, we have to insert the boundary and horizon expansions (2.21)-(2.24) The divergent parts are canceled by appropriate counterterms given by (2.4). Hence the final result for the on-shell action density reads Note that the final result still contains an integral and hence the Lagrangian does not seem to reduce to a total derivative. However, the integral expression is actually a boundary term in x 2 direction. We checked this explicitly by computing the Noether charges along the lines of [50,51].

C Special cases
In this section, we look at particular limits of our system of equations in order recover the results already available in the literature. We begin with the discussion of the special case B = 0 corresponding to the helical black brane solutions studied by Donos and Gauntlett in [19]. Afterwards, we comment on the case B = 0 and its quantum critical point observed in the studies of the charged magnetic branes solutions by D'Hoker and Kraus in [11]. In our coordinates, the helical black brane solution [19] corresponds to the choice of In this case, the line element (2.15), the field strength tensor (2.16) and the gauge field (2.17) read Note that our ansatz differs slightly from the one presented in [19]: first, we compactify r by r = 1/z with z ∈ [0, 1], second we re-label the one-forms: ω 1 here = ω 2 DG , ω 2 here = ω 3 DG and ω 3 here = ω 1 DG and third we use a slightly different metric ansatz which is closer to the one used by D'Hoker and Kraus (see appendix C.2). In particular, while in [19] the metric components in terms of z satisfy z 4 g tt g zz = −f (z) 2 , we choose an ansatz such that z 4 g tt g zz = −1 and z 4 g ω 2 ω 2 g ω 1 ω 1 = v(z) 4 . Finally, the authors of [19] use scaling freedom of the coordinates to set µ = 1, while the coordinate location of the event horizon is not known a priori. In contrast we fix the horizon to be located at z = 1, and hence we are not allowed to set µ to one. Our numerics pass an important check: we can reproduce all their results down to temperatures of orderT ≈ 10 −5 . The authors of [19] reported some difficulties in studying the behaviour of the solutions in the regime of very low temperatures,T → 0. In this regime, some functions develop strong gradients around the horizon. An example is depicted in figure 9 for the results shown in section 3.2 (γ = 1.5). As we drop the temperature, the functionq(z) becomes steeper around z = 1. Numerically, the solution must be obtained either by a massive increase in the number of grid points or by the development of specific techniques adapted to this drawback. In this paper, we use the so called analytical mesh-refinements [52,53], described in appendix D to circumvent these problems.

C.2 Quantum critical point
We now turn our attention to the charged magnetic solution constructed by D'Hoker and Kraus [11]. In our conventions, it corresponds to the choice of q(z) = g(z) = b(z) = 0 and α(z) = 1, which leads to 16 16 One must be careful with the definition of F and A within the action. Our Sgrav, given by eq. (2.3), coincides with Donos and Gauntlett [19], which in turn differs from the one used by D'Hoker and Kraus (DK) [11]. It is crucial to keep in mind that F here = 2FDK and A here = 2ADK. The relevant quantities µ, ρ and B pick up this factor of 2 accordingly.  Figure 9. Functionq(z) normalised by its value on z = 1. In the regimeT → 0, the function develops very strong gradients around the horizon. The numerical solution requires higher resolution and, eventually, specific techniques must be employed (such as the analytical mesh-refinement -see appendix D).
Moreover, note that D'Hoker and Kraus normalise all their dimensionful quantities by the charge density ρ = e (0). For example, the dimensionless magnetic fieldB is given bŷ , whereas we consider the dimensionless magnetic fieldB given byB = B µ 2 . As usual, the charge density ρ and the chemical potential µ are thermodynamically conjugate. In other words, one may consider the chemical potential to be a function of the charge density or vice-versa. The choice of perspective, i.e. canonical ensemble versus grand canonical ensemble, does not affect the dynamics of the field theory.
Due to different normalisations used we have to be careful when using results from [12]. For example, we first have to translate the locationB C of the quantum critical point quoted in [12]. The corresponding value ofB C in our notation is related toB C byB C = γ 2B3 C . For instance, for γ = 1.7 the quantum critical point is located atB C ≈ 0.400 which corresponds toB C ≈ 0.185 in our notation. For γ = 1.5 one obtainsB C ≈ 0.461 or B C ≈ 0.220, respectively. The location of the quantum critical point is shown in the phase diagram figure 4. To check our numerics and the correct normalisation we determine the behaviour of the entropy densitys as a function ofT close to the quantum critical point B ≈B C . We display in figure 10 the entropys as a function ofT for γ = 1.5. Needless to say that we reproduce the results from [11].
A final comment concerning the behaviour of T tx 3 given the different normalisations.
As numerically checked alongB = constant, is also a constant.
However, alongB, the left panel of figure 11 shows that neither T tx 3 nor T tx 3 = T  Figure 10. Entropys as a function of the temperatureT for the charged magnetic brane. For γ = 1.5 the quantum critical point is located atB C ≈ 0.220 and the entropy decays ass ∼T 1/3 . ForB = 0.215 <B C the entropy goes to a constant in the low temperature regime, whereas for B = 0.225 >B C one hass ∼T . For high temperatures we find the expected behaviours ∼T 3 . All results are in agreement with [11]. should still hold. We explicitly checked the validity of this expression in the right panel of figure 11 for γ = 1.5 along the critical valueB = 0.461. In particular, the inset displays the difference 1 − 8 c 4 γ B µ 2 , which is limited only by the numerical round-off error.

D Numerical details
As mentioned before, the system of equations (A.3) expressed in terms of the auxiliary variables (A.7) is solved numerically by means of a spectral method. In this section, we elaborate further on this topic and we give more details on the numerics involved.

JHEP03(2016)164
Spectral methods are best applied to differential equations whose solutions are known to be analytic. In such a case, the error coming from the numerics decays exponentially as one increases the grid resolution. On the other hand, the presence of logarithmic terms spoils this properties, rendering a merely algebraic convergence rate. Yet, the introduction of the auxiliary variables (A.7) removed the leading z 4 ln(z) terms and we verify that our solutions show typically a rather efficient convergence rate, even for large values ofB.
In addition to the high accuracy, spectral methods are also flexible enough to deal with other unknown parameter apart from the field functions. In order to fix these parameters, one needs to specify extra conditions together with the equations of motion. As a global scheme, the method makes no distinction between the unknown functions and parameters and solves the system of all variables at once. In the context of gauge/gravity dualities, it is possible to address the low temperature behaviour within spectral methods. In addition, spectral methods also allow to include singular points and hence the equations of motion can be solved in the whole domain. More details on spectral methods can be found in [54,55] as well as in the specialised reviews [56,57].

D.1 Spectral methods
Let n fields be the total number of unknown functions X I (with I = 1, . . . , n fields ) defined on the domain z ∈ [0, 1]. Let us also assume the existence of n par unknown real parameters χ A (with A = 1, . . . , n par ). Moreover, we consider the following system of equations F I (X I , X I , X I , z; χ A ) = 0 for 0 < z < 1 , F I 1 (X I , X I ; χ A ) = 0 for z = 1 , Φ A (X I , X I , X I ; χ A ) = 0 .

(D.1)
Here, X I and X I are, respectively, the first and second derivative of the functions X I (z). F I 0 and F I 1 are the boundary conditions, whereas F I represent the equations of motions. Finally, Φ A stand for the extra conditions that fix the unknown parameters.
In order to solve numerically the system of equations (D.1), we first provide a numerical resolution N and expand the functions X I (z) as In the expansion above, the basis functions are the Chebyshev polynomials of first kind T k (ξ) = cos[k arccos(ξ)], ξ ∈ [−1, 1], while R I (z) are the residual functions. To determine the Chebyshev coefficients c I k , we specify a set of grid points z i (with i = 0, . . . , N ) and impose that the residual function vanishes exactly at the grid points, i.e., R(z i ) = 0. In other words, at the grid points, the unknown functions are given exactly by the spectral representation The Chebyshev coefficients are then obtained after the inversion of (D.3). From the c I k , we can obtain the coefficients c I k and c I k describing the spectral representation of the derivatives X I (z) and X I (z) as in eq. (D.2). In this paper, we work with the Chebyshev-Lobatto grid points Now we can combine the function values X I i := X I (z i ) and the parameters χ A into the single vector X of length n total = (N + 1)n Fields + n par from which we can also form the vectors X and X representing the discrete spectral derivatives with respect to z. These vectors are finally used to evaluate the system of equations (D.1) at the grid points (D.4), giving us a non-linear algebraic system to be solved for the components of the vector X. The solution of the algebraic system (D.6) is obtained by a Newton-Raphson method, i.e. given an initial guess X 0 , the solution is iteratively approximated by The inversion of the Jacobian matrixĴ( X) = ∂ F /∂ X is performed with a LU decomposition. One can show that the Newton-Rapshon scheme always converges, providing the initial guess X 0 is sufficiently close to a solution.

D.2 Numerical solution: the initial guess
In this work, we are looking for the numerical solution of the n Fields = 10 metric and gauge field functions. IfB = 0, the electrically charged Reissner-Nordstroem black brane (2.10) is always a solution, regardless ofk and γ. Similarly, forB = 0 one always obtains charged magnetic black brane as trivial solution. The interesting cases are those values of γ andk for which a non-trivial solution with b(z) = 0 and q(z) = 0 also exists. Unfortunately, our first experiences showed that, for given boundary valueẼ(0) = −µ, the Newton-Raphson method always converged to a trivial solution of the non-linear algebraic system (D.6), regardless of the initial guess X 0 . In order to obtain the new solution describing the condensed phase, an extremely careful fine-tuning to the initial guess X 0 seems to be needed.
To solve this issue, we consider the parameter µ as a further unknown variable (thus n par = 1) and impose the extra conditionb(0) = b 2 = ε. By fixing a value ε = 0, we enforce that the Newton-Raphson scheme will necessarily converge to the non-trivial solution.
For a givenk, γ, we start withB ∼ 0, and ε ∼ 0. Then, the new solution should be just a small perturbation of the AdS Reissner-Nordstroem spacetime. Therefore, our -28 -

JHEP03(2016)164
initial-guess constitute of eqs. (2.12)-(2.13), with the slight modificationb(z) = ε. Besides, we must also provide an initial-guess for the variable µ. In all our experiments, the value µ = 2 was sufficient for the convergence of the Newton-Raphson scheme. Once a solution is available, we can use it as initial-guess for a modified set parameters γ,B,k, ε .
From the numerical point of view, fixing γ,B,k, ε is an efficient method to find the non-trivial solution. However, from the physical perspective, a system with a constant temperatureT , i.e. specified by γ,B,k,T is what one really wants to describe. Note that, as an alternative to the extra conditionb(0) = b 2 = ε, we can indeed imposeT =constant. This corresponds to looking for the value of µ leading to a solution with a fixed valuẽ u(1). Unfortunately, this approach does not guarantee that the method will give us the non-trivial solution. Depending on how far the initial-guess is from the trivial solution, the Newton-Rapshon scheme might converge to the charged magnetic black brane solution (withb(z) = 0). Therefore, for a fixed {γ,B}, our algorithm is a combination of both possibilities and can be divided into three stages: I) Phase boundary: we set ε = 10 −9 and scan the values ofk ∈ [k 0 ,k 1 ] to get the phase boundary. 17 With the knowledge ofT (k) we find the point {k C ,T C } for which T C =T (k C ) is at a maximum.
II) Condensed phase: we then keepk =k C fixed and find new solutions inside the phase by slowly increasing the values of ε until a given 0 <T 0 <T C is achieved. Typically, we setT 0 = 0.95T C .
III) Constant temperature: with the solution inside the phase provided by step II as initial guess, we no longer need to keep ε fixed. We now solve at surfaces of constant T ∈ [T min ,T 0 ] ∪ [T 0 ,T C ] in a given intervalk ∈ [k 0 (T ),k 1 (T )] and find the physical statek * (T ) for which the grand canonical potentialΩ * =Ω(k * ) is at a minimum.
To illustrate the solutions, we show in figures 12 the results for the metric and gauge field functions with γ = 1.5 for the following two configurations: giving respectively µ = 2.172 and µ = 2.982 for the chemical potential.

D.3 Numerical error
Finally, we discuss the accuracy of our method. Given a a high resolution N max , we consider a reference solution X I (z; N max ), ρ(N max ) and define, for a lower resolution N < N max , the numerical error of the solution X I (z; N ), µ(N ) by    . ForB = 0 we obtain the expected exponential convergence rate. ForB = 0 one expects a merely algebraic decay due to the logarithmic terms. Its contribution, however, enters only within the machine precision. Figure 13 displays the error for the configurations mentioned in the previous section. The left panel shows the caseB = 0 and we see the typical exponential convergence provided by the spectral method. The right panel depicts the caseB = 0.275. Despite the presence of logarithmic terms, the convergence rate is very efficient and we do not observe a significant influence of an algebraic decay within the machine limits imposed by round off errors.
Even though the method provides a high accuracy solution for moderate values ofT , we note that the small temperature regime requires a massive increase in resolution. This feature becomes evident in figure 14, where we compare the convergence rate (for instance for µ) at different temperatures. As already mentioned in appendix C.1 and illustrated in figure 9, the main reason is the presence of strong gradients around the horizon z = 1. One technique to deal with such strong gradients is the so called analytical mesh-refinement, described in [52,53]. It consists of mapping the coordinate z ∈ [0, 1] into ζ ∈ [0, 1] via By choosing an adequate parameter λ, the mapping increases the number of grid points around z = 1 and smoothen out the solution. In our case we set λ = b 2 ln(T ) . In figure 15 one sees the significant improvement of the convergence rate, specially in the casē B = 0. ForB = 0, the method is still effective at low temperatures, but it also intensifies the algebraic decay rate introduced by logarithmic terms. Even though the analytical mesh-refinement provides the necessary tools to study the low temperature limit within the scope of this work, we note that there are still possibilities for enhancing the accuracy of theB = 0 case. 18 18 An option is to use a multi-domain code with another coordinate map to remove the logarithmic terms [58].   Figure 15. Numerical error for low temperatures with the analytical mesh-refinement. Left panel: forB = 0 the convergence rate is significantly improved and we achieve the machine precision with a moderate number of grid points. Right panel: forB = 0 we still obtain a better convergence rate, though the method worsens the algebraic decay due to the logarithmic terms.

JHEP03(2016)164
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.