Holographic Charged Fluid with Chiral Electric Separation Effect

Hydrodynamics with both vector and axial currents is under study within a holographic model, consisting of canonical $U(1)_V\times U(1)_A$ gauge fields in an asymptotically AdS$_5$ black brane. When gravitational back-reaction is taken into account, the chiral electric separation effect (CESE), namely the generation of an axial current as the response to an external electric field, is realized naturally. Via fluid/gravity correspondence, all the first order transport coefficients in the hydrodynamic constitutive relations are evaluated analytically: they are functions of vector chemical potential $\mu$, axial chemical potential $\mu_5$ and the fluid's temperature $T$. Apart from the proportionality factor $\mu\mu_5$, the CESE conductivity is found to be dependent on the dimensionless quantities $\mu/T$ and $\mu_5/T$ nontrivially. As a complementary study, frequency-dependent transport phenomena are revealed through linear response analysis, demonstrating perfect agreement with the results obtained from fluid/gravity correspondence.


Introduction
Fluid dynamics is an effective low energy description of most interacting systems at finite temperature. Within such a hydrodynamic approximation, the entire dynamics of a microscopic theory is reduced to that of macroscopically conserved currents, such as the stress-energy tensor and charge current operators computed in a locally near equilibrium thermal state. An essential ingredient of any fluid dynamics is the constitutive relations which relate the macroscopically conserved currents to the hydrodynamic variables, such as fluid velocity and charge densities, and to external forces like external electromagnetic fields. Derivative expansion in the hydrodynamic variables and external forces accounts for deviations from thermal equilibrium. At each order, the derivative expansion is fixed by thermodynamics and symmetries, up to a finite number of transport coefficients such as viscosity and diffusion coefficients. These transport coefficients are not calculable from hydrodynamics itself but have to be deduced experimentally or computed from the underlying microscopic theory.
For relativistic fluid dynamics, the stress-energy tensor is conveniently parameterized as where u µ , E, P are the velocity, energy density and pressure of the fluid, and P µν = u µ u ν + η µν is the projection tensor. Up to the first order in the derivative expansion, the viscous component π µν takes the form, where η, ζ are the shear viscosity and bulk viscosity, respectively. Throughout this work we will take the Landau-Lifshitz frame so that u µ π µν = 0. For the fluid with conserved charges, one also needs to specify constitutive relations for the associated currents. Indeed, the charge transport properties are found to be useful in probing the structure and dynamics of matter. A well-known example is the Ohm's law J i ≡ ψ γ i ψ = σE i , which states the generation of an electric current in response to an external electric field for a normal conducting media. In recent years, exploration of other possible electric current generation, particularly for a system with charged chiral fermions, has attracted much interest. It turns out that the celebrated microscopic chiral anomaly induces fascinating anomalous transport phenomena that break the space parity symmetry. One such example is the chiral magnetic effect (CME) [1][2][3][4]: generation of an electric current directed along an externally applied magnetic field, J = ξ B B. The existence of CME relies on chirality imbalance between left-and right-handed chiral fermions, usually parameterized by an axial chemical potential µ 5 . For a rotating hydrodynamic flow, chiral anomaly induces a vector current along the fluid vorticity, J = 1 2 ξ ∇ × u, which is called the chiral vortical effect (CVE) [5][6][7].
On the other hand, an axial current J i 5 ≡ ψ γ i γ 5 ψ also exists for a system with charged chiral fermions. In fact, via the chiral anomaly effect, an axial current is generated along an external magnetic field, J 5 = ξ 5B B, which is referred to as chiral separation effect (CSE) [8,9]. Interestingly, the interplay of CME and CSE predicts a gapless wave mode propagating along the magnetic field, called chiral magnetic wave (CMW) [10]. In heavy-ion collisions, the CMW induces an electric quadrupole moment of the quark-gluon plasma [11], leaving experimentally observable effects [11][12][13][14]. It is important to stress that chiral anomaly-induced transport phenomena summarized above are non-dissipative [7], that is they do not contribute to entropy production. While observable signatures predicted by the CME have not been conclusively detected in heavy-ion collisions [15][16][17][18], the CME may explain a large negative magnetoresistance observed in Dirac and Weyl semi-metals [19][20][21]. We recommend recent reviews [22][23][24][25][26][27] and references therein on the subject of anomalous transport phenomena.
The separation of chiral charge could also be induced by an external electric field, J 5 = σ 5e E, which is called chiral electric separation effect (CESE) [28,29]. However, it is important to emphasize that the CESE does not originate from chiral anomaly but is simply the result of conduction of a chiral many-body environment [28]. On very general grounds, the CESE exists only when both vector and axial chemical potentials are nonzero. Based on Kubo formula, the CESE conductivity has been computed for thermal quantum electrodynamics (QED) [28] and quantum chromodynamics (QCD) plasmas [29] at high-temperature regime up to leading-log order in gauge couplings. In the strongly-coupled regime, the CESE was studied in [30,31] by using the holographic QCD model of [32,33]. Recently, the CESE conductivity is also computed within the framework of kinetic theory [34][35][36][37] under the relaxation time approximation (RTA). If the CESE conductivity is parameterized as σ 5e = µµ 5χA , then all these studies indicate thatχ A depends on µ, µ 5 weakly. In analogy with CMW, the CESE, combined with the Ohm's law, results in a new gapless wave mode propagating along the electric field, which is called density wave [28] or chiral electric wave (CEW) [31].
The transport phenomena reviewed above could be summarized into hydrodynamic constitutive relations for vector and axial currents. In the Landau-Lifshitz frame where u µ J µ = −ρ and u µ J µ 5 = −ρ 5 , we are ended up with [38][39][40], (1.4) where ρ, ρ 5 are vector and axial charge densities. The external electromagnetic fields E µ , B µ and the fluid's vorticity ω µ are (1. 5) In (1.3)(1.4), the Wiedemann-Franz law has been used to relate thermal conductivity to electrical conductivity. Theσ 5e -andσ e -terms are relevant to chiral charge diffusions, while ξ 5 -term is an axial analogue of CVE. Time/space evolution of the system is determined by solving the hydrodynamic equations where the axial current is not conserved due to chiral anomaly effect and C denotes the anomaly coefficient. ξ, ξ 5 , ξ B , ξ 5B -terms are chiral anomaly-induced and correspond to nondissipative transport phenomena, particularly making zero entropy production. In contrast, σ e , σ 5e ,σ e ,σ 5e are dissipative transport coefficients and have to be determined by the underlying microscopic theory and will be the focus of present work. We would like to point out that studies in references [28][29][30][31] were carried out under the decoupling limit (or probe limit), where the vector and axial currents were taken as decoupled from the stress-energy tensor of the fluid. While the decoupling limit allows deriving a simple Kubo formula for CESE conductivity, it does result in loss of some important physical contents, such as a pole in DC electrical conductivity due to spatial translation invariance [41][42][43]. In the probe limit, the authors of [30,31] found a nonzero CESE conductivity for the holographic QCD model, which consists of probe D8/D8-branes in the background geometry of N c stacks of D4-brane. A crucial point of [30,31] in realizing CESE is about the explicit interaction between the vector and axial gauge fields in the bulk, which originates from the nonlinear Dirac-Born-Infeld (DBI) action on the world-volumes of D8/D8. Indeed, it was found that the CESE conductivity does vanish in a simple holographic setup with canonical U (1) V × U (1) A Maxwell fields probing the fixed Schwarzschild-AdS 5 background [44][45][46].
In this work, we reconsider the holographic U (1) V × U (1) A model and take into account the gravitational back-reaction effect on the gauge sectors of the bulk theory. While the anomaly coefficient C is fixed by the microscopic theory, we will take it as a free parameter and turn it off for the present study. Equivalently, the bulk Chern-Simons actions will be neglected in the holographic setup of [44][45][46]. Note that there is no explicit interaction between U (1) V and U (1) A gauge fields in the bulk action. However, as will be clear later, the gravitational back-reaction will result in a coupling between perturbations of U (1) V and U (1) A bulk fields. This is the mechanism of generating CESE in such a simple holographic model. In the next subsection, we will summarize the results and make the comparison with previous works. The remaining sections supplement all calculation details.

Summary of the results
Except for those chiral anomaly-induced terms, we will re-derive the first order constitutive relations (1.1),(1.3),(1.4) via fluid/gravity correspondence [47]. For a specific holographic model to be introduced in Section 2, we analytically compute all dissipative transport coefficients. The calculation details will be presented in Section 3. The shear and bulk viscosities in (1.2) take universal values of holographic conformal fluids (holographically described by two-derivative Einstein gravity), where s is the entropy density of the dual fluid. The dissipative transport coefficients in (1.3) and (1.4) could be presented in several different forms. In the first form, they are where σ Q and σ 0 are given by It is important to stress that σ e , σ 5e andσ e ,σ 5e are intrinsic conductivities, which are usually referred to as "quantum critical" [48] or "incoherent" [49,50] conductivities in holographic framework. We have reserved the notation σ Q to the intrinsic electrical conductivity for the single charge case, since σ e = σ Q once ρ 5 = 0 or µ 5 = 0. In the probe limit where µ = µ 5 = 0, we have σ e =σ e = σ 0 = πT , but σ 5e ,σ 5e vanish.
In (1.8), the relation σ 5e =σ 5e originates from Onsager reciprocal relation for a timereversal symmetric system. Additionally, there are mirror symmetric relations σ e ↔σ e and σ 5e ↔σ 5e under the exchange µ ↔ µ 5 , which are specific to our holographic model 1 . Thanks to these symmetric relations, in the discussions below we will focus on σ e , σ 5e . As will be clear in Section 4, these symmetric relations for DC transport coefficients are extendable to associated alternating current (AC) conductivities in (4.13). While the second law of thermodynamics, i.e., non-negativeness of divergence for the entropy current, could not fix the values of dissipative transport coefficients, it does set constraints for them [29,30], which are obviously satisfied by our results (1.8)(1.9). From the dual fluid's point of view, it is more natural to re-express σ e ,σ e , σ 5e ,σ 5e in terms of chemical potentials and temperature. As a result, (1.8) are turned into (1.12) In (1.11) we defined χ A to measure deviation of the CESE conductivity σ 5e from the conjectured universal factor Tμμ 5 [28,29]. We also defined χ V to measure deviation of the Ohmic conductivity σ e from its probe limit. In the high temperature (small chemical potential) limit or low temperature (large chemical potential) limit, the conductivities (1.11) behave as, σ e πT, σ 5e −πTμμ 5 ,μ,μ 5 1, σ e πTμ 2 5 6(μ 2 +μ 2 5 ) , σ 5e − πTμμ 5 6(μ 2 +μ 2 5 ) ,μ,μ 5 1.
For illustration, in Figure 1 we plot σ e /(πT ) and −σ 5e /(πT ) as functions of dimensionless chemical potentialsμ andμ 5 . As seen from Figure 1, the axial chemical potentialμ 5 has the effect of enhancing σ e while the vector chemical potentialμ diminishes it. On the other hand, the CESE conductivity σ 5e depends on bothμ andμ 5 non-trivially. Particularly, the deviation factor χ A is a monotonic decaying function ofμ andμ 5 . In the high-temperature regime whereμ,μ 5 1, the factor χ A could be thought of as unity, which is obviously violated in the lower temperature limit. In Figure 2 we show density plots for the deviation factors χ V and χ A .
Below we compare our results with those obtained within various other models [28][29][30][31], see Table 1. The first two calculations in Table 1 were based on perturbative thermal QED and   QCD, respectively. To leading-log order in gauge couplings, the conductivities were computed using the Kubo formula in [28,29]. The rest calculations in Table 1 are based on specific holographic models: Sakai-Sugimoto (S-S) model in [30] versus holographic U (1) V × U (1) A model of present work. As shown in Table 1, the results obtained from S-S model show the behavior of the pre-factors in σ e , σ 5e ∝ N 2 c g 2 YM [30], which is quite different from those in weakly coupled regime [28,29]. While in our case, the conductivities depend on the bulk gauge couplings q V , q A . If we employ the top-down setups in [52,53], they could be related to the parameters of boundary theory by 1/q 2 , where N c and N f are the numbers of the colors and flavors of the dual theory.
Although these dissipative transport coefficients are model dependent, in the high temperature regime, they do share some universal features, such as the linear dependence in T (for σ e ) and Tμμ 5 (for σ 5e ). In addition, the Ohmic conductivity σ e has to be positive [28], but the sign of CESE conductivity σ 5e could not be fixed by the second law of thermodynamics [28]. Particularly, σ 5e was found to be positive for the models in [28][29][30], but it is shown to be negative from present study (1.11). As shown in Table 1, from weak to strong coupling, the conductivities σ e , σ 5e show different dependence on the gauge coupling of boundary gauge theory in various models.  13.0 20.5 e 3 ln(1/e) 14.5 Note: In this table, e, g c , g YM are the gauge couplings of QED, QCD, the dual SU (N c ) gauge theory, respectively. For the calculations in QGP with two light quarks (u, d), Q e = Diag(2/3, −1/3), Q V and Q A are the vector and axial charge matrices in flavor space. For the S-S model, the results above are read off from relevant numerical plots around T 0.2GeV [30]. Actually, in the high-temperature regime where the chemical potentials are suppressed, σ e ∝ T 2 /m KK and m KK is the Kaluza-Klein mass in the S-S model [31,51]. For the results of present study, we have restored the bulk gauge couplings q 2 V and q 2 A , which are related to parameters of boundary theory by 1/q 2 [52,53].
From Table 1, in the high temperature regime that T µ, µ 5 , the conductivities in holographic U (1) A × U (1) V model share the same temperature scalings as those in the QED plasma or QGP. The Ohmic conductivity depends linearly on the temperature, σ e ∝ T , while the CESE conductivity depends inversely on the temperature, σ 5e ∝ µµ 5 T . It will become more evident from dimensional analysis as implemented in [29] within relativistic kinetic theory. In the high temperature regime, the conductivities can be expanded in terms of the ratios T , where f (T ) is a function of the temperature and d 02 , d 20 , d 11 are constants [29]. We know that the conductivities have the same dimension as the temperature or chemical potentials, So, if there are no extra dimensional physical quantities in the model, it is valid to conclude that f (T ) T , which also supports the general scaling behaviors as summarized in Table 1.
On the universal scaling behaviors in the high temperature regime, for a physical explanation, it is nature to assume that the interaction between the charges in the fluid will become weaker in the higher temperature. The Ohmic conductivity σ e describes the response of the charged vector current to the external electrical field, σ e J/ E. When the temperature is increased, σ e will be enhanced because the moving of the charges in the fluid will become easier. However, the CESE conductivity σ 5e measures the response of the axial current to the external electrical field, σ 5e J 5 / E. When temperature is increased, the interaction between the axial current and electrical field will become weaker and σ 5e will decrease due to the additional interaction factor μ T μ 5 T . While the CESE is a non-anomalous transport phenomenon, it may induce phenomenological consequences in heavy-ion collisions, namely the net charge distribution and correlation patterns in Cu+Au collisions as discussed in [28]. Admittedly, this should be a mixture due to CESE, CME, and CSE. However, CESE, along with other robust anomalous transport phenomena, is masked by various backgrounds in heavy-ion collisions, making it very difficult to pin down, not even to explore its properties. On the other hand, the exotic topological states of metal, such as Dirac and Weyl semi-metals, provide an experimental playground to study potential observable effects of CESE and other anomalous transports in a controllable way. See [27,[35][36][37] for recent progress on this topic, as well as the relevant investigations from holographic models [54][55][56][57][58][59][60][61][62].
Below we would like to rewrite the currents (1.3) and (1.4) in a linear response form, in which the electric field E i and thermal gradient ∇ i T are taken as external sources. In other words, the fluid velocity u i will be eliminated using the current conservation law in (1.6). Here, the chemical potentials µ, µ 5 will be taken as constant. Consequently, where ∂ t → −iω is used. Then, the currents (1.3) and (1.4) turn into whereσ e is presented in (1.8). Physically, σ (ω) is the low-frequency limit of the Ohmic electrical conductivity, and σ 5(ω) measures the chiral electrical separation effect. Obviously, in the probe limit, σ (ω) and σ 5(ω) are dominated by their intrinsic parts σ e and σ 5e . α (ω) and α 5(ω) are the thermoelectric conductivities for vector and axial currents. The heat current is where the transport coefficients arē κ (ω) is the low-frequency limit of thermal conductivity, and it is fully determined by the intrinsic conductivities σ e , σ 5e in (1.8). Here we would like to emphasize that σ (ω) , σ 5(ω) , σ (ω) , α (ω) , α 5(ω) ,ᾱ (ω) ,κ (ω) are different from the intrinsic conductivities σ e , σ 5e ,σ e ,σ 5e : while the former could be directly read off from associated Kubo formulas, the latter are useful in parameterizing the hydrodynamic constitutive relations (1.3)(1.4). In the probe limit, the differences between them are absent.
Note the appearance of i/ω pieces in the physical conductivities σ (ω) and σ 5(ω) . Via Kramers-Kronig relations [63,64], this i/ω immediately implies the existence of a δ(ω) in real parts of σ (ω) and σ 5(ω) , as is required by translation invariance. However, in holographic models with spatial translational symmetry, it is hard to reveal the delta peak in DC conductivity analytically. Once translational symmetry is broken, the DC conductivities become finite with the delta peak removed. Within the relaxation time approximation, the momentum dissipation effect would result in the replacement where τ corresponds to momentum relaxation time. Now all physical conductivities become finite and, in particular, they are split into two parts: the "coherent" pieces (related to the momentum dissipation) and the "inherent" ones (the universal pieces). Admittedly, a more rigorous treatment of momentum dissipation along the line [41-43, 65, 66] would be useful in clarifying physical meanings of these transport coefficients, and we will address this elsewhere.
The remaining sections are structured as follows. In Section 2, we present the holographic model. In Section 3, with anomalous terms neglected, we re-derive the first-order constitutive relations (1.1) (1.3) (1.4) using the fluid/gravity correspondence, and analytically compute all dissipative transport coefficients. In Section 4, we obtain AC conductivities through linear response analysis. Section 5 contains the conclusion and discussions. Two appendices A and B provide further calculation details.

Holographic Model for Fluid with
We consider (4 + 1)-dimensional Einstein gravity with a negative cosmological constant Λ = −6/L 2 in the bulk, which is coupled to U (1) V × U (1) A gauge fields (see, e.g., [44,45]). The total bulk action is The field strengths of the bulk U (1) gauge fields are defined as In our notations, A M andÃ M denote the vector and axial bulk gauge fields, which are dual to vector and axial currents J µ , J µ 5 of the boundary conformal field theory (CFT), respectively. The Gibbons-Hawking-York term S K is where K = γ µν K µν . γ µν is the induced metric on the boundary hypersurface Σ defined by the equation r = r c , and K µν is the extrinsic curvature tensor on Σ, where L n is the Lie derivative along the unit normal vector n M of the hypersurface Σ: The counter-term action S ct is [53, 67-71] where R is the Ricci scalar of the induced metric γ µν . Note minimal subtraction scheme has been utilized in writing down the counter-terms for bulk Maxwell fields. F µν andF µν in (2.6) are the projections of bulk field strengths F M N andF M N onto the hypersurface Σ. We would like to stress once again that the possible Chern-Simons termsÃ ∧ F ∧ F and A ∧F ∧F have been ignored in (2.1), which amounts to switching off the chiral anomaly effect in the dual boundary theory [44,45]. Indeed, in order for the chiral anomaly to take effect, the dual plasma should either be exposed to a magnetic environment or rotate. Thus, with an electric field as the only external source, all the anomalous transports in (1.3)(1.4) vanish accidentally.
Variation of total bulk action S M in (2.1) with respect to bulk metric g M N gives rise to the Einstein equation, The bulk stress-energy tensor has been denoted as T bulk M N andT bulk M N , which should not be confused with that of the dual boundary theory. q andq measure the strength of backreaction of bulk gauge fields on the bulk geometry. The Maxwell equations for A M andÃ M are, According to Anti-de Sitter/Conformal Field Theory (AdS/CFT) correspondence [72][73][74], the expectation values of stress-energy tensor and currents on the boundary theory are defined as In terms of bulk fields, (2.10) turns into where the counter-term which will affect the study of Section 4 beyond first-order transport coefficients. In (2.11)-(2.13), G µν is the Einstein tensor associated with the induced metric γ µν , and D is the covariant derivative operator compatible with γ µν . The bulk equations (2.7)-(2.9) can be classified into dynamical components and constraint ones. Moreover, the constraint components correspond to conservation laws for boundary stress-energy tensor and currents: να are external vector and axial electromagnetic field strengths, respectively. In what follows, we will work under the convention 16πG 5 = 1 and L = 1. The bulk gauge couplings q V , q A will be absorbed into redefinitions of bulk gauge fields A M → q V A M andÃ M → q AÃ M . The boundary CFT in thermal equilibrium corresponds to a homogeneous solution of the bulk theory (2.1). We assume the presence of finite vector and axial charge densities. Consequently, the homogeneous solution of the bulk theory is the AdS 5 black brane with two charges, where M, Q,Q are constant parameters of the bulk theory. In (2.16), r h is the largest root for f (r) = 0, defining the location of event horizon of the two-charges AdS 5 black brane. The Hawking temperature, identified as the temperature of dual boundary field theory, is which should be non-negative, setting constraints on the combination (Q 2 +Q 2 ). At the horizon with a constant time section, the line element degenerates into ds 2 | Horizon = r 2 h δ ij dx i dx j , with i, j = 1, 2, 3. Thus, the entropy density of the black brane (2.16), which will be identified as that of the dual CFT, turns out to be where in the second equality we made use of the normalization convention 16πG 5 = 1. The subscript " 0 " in (2.17) (2.18) is to emphasize that they are constant, as compared to temperature field T (x α ) in Section 3. Finally, via (2.11)-(2.13), the dual stress-energy tensor and currents of the boundary theory are where the subscript " (0) " is to mark that these quantities are associated to a state in thermal equilibrium. If we make the following identifications,

First order hydrodynamics from fluid/gravity correspondence
In this section, we construct the first-order hydrodynamics dual to the bulk theory (2.1) via the fluid/gravity correspondence [47,75].

Set Up the Fluid/Gravity Calculations
In this subsection, we set up the stage for performing fluid/gravity calculations for the bulk theory of (2.1). Following the standard procedure of fluid/gravity correspondence [47,75], we make a Lorenz boost transformation for the static solution (2.16) along the boundary coordinates L µ ν is the Lorentz boost matrix, which has been parameterized via a four-velocity u µ . Note the four-velocity u µ is a time-like unit vector obeying η µν u µ u ν = −1, which leads to u 0 = − √ 1 + u 2 . After the transformation (3.1), the homogenous solution (2.16) turns into where a constant vector field A One key procedure of fluid/gravity correspondence is to promote the constant parameters in (3.2) to arbitrary functions of boundary coordinates [7,47,75,76], where we have chosen the frame that u i = 0 and A ( ) µ = 0 at the origin. Moreover, a formal parameter ε is introduced to count the number of derivatives in the expansion, and eventually will be set to unity. All calculations in this section will be accurate up to the first-order in the derivative expansion. Consequently, up to order O(ε 1 ) the promoted metric and gauge fields become where As explained below (3.3), in order to satisfy the equations of motion (2.7)-(2.9), suitable corrections must be added on top of (3.5). Up to the first-order in the derivative expansion, where α ij (r) is a traceless symmetric tensor of rank two under SO(3) rotational symmetry between the spatial directions x i . In the parameterizing of the corrections (3.7), we choose the following gauge convention Plugging the total bulk metric and gauge fields into (2.7)-(2.9) results in a system of ordinary differential equations for those corrections in (3.7). We need to specify suitable boundary conditions in order to fully determine the corrections. The first type of boundary condition is the requirement of asymptotic AdS, which fixes the large r behavior for the corrections, The second type of boundary condition is the regularity requirement for all the corrections in (3.7), h(r), k(r), j i (r), α ij (r), a t (r), a i (r),ã t (r),ã i (r) are regular over r ∈ [r h , +∞), (3.11) which turns out to be effective at the event horizon r = r h . The remaining ambiguity of determining the corrections of (3.7) will be fixed by the frame convention. We will work in Landau-Lifshitz frame so that where E, ρ, ρ 5 are the energy density, vector and axial charge densities of the fluid, respectively, Note the identification made in (3.13) is the promotion from their in-equilibrium counterparts (2.20). Up to the first-order in the derivative expansion, the Landau-Lifshitz frame conditions (3.12) turn into constraints on the form of boundary stress-energy tensor and currents, 14) For the sake of later presentation, we rewrite the expressions (2.11)-(2.13) in terms of those corrections in (3.7), as well as (3.20) The limit of r → ∞ is assumed implicitly in expressions above (3.16)-(3.20), and we have also ignored the terms that will be explicitly vanishing as r → ∞.
While the bulk corrections will be constructed around x µ = 0, they do contain enough information to write down the total bulk metric and gauge fields about any point, valid up to the first-order in the derivative expansion. Instead of following this approach of [47], we will compute the boundary stress-energy tensor and currents using thus-constructed solutions, via the formulas (3.16)- (3.20). Eventually, we will lift up thus-obtained constitutive relations into a covariant form.

First-Order Charged Fluid: CESE and Other Conductivities
Following Section 3.1, it is straightforward to solve the bulk equations (2.7)-(2.9) and obtain the corrections in (3.7). In what follows, we summarize the final results and leave all the calculation details in Appendix A.1. We present by grouping them into different sectors under SO(3) symmetry of the boundary spatial directions.
In the scalar sector, In the vector sector, thanks to going beyond the probe limit, dynamical equations for j i , a i ,ã i are coupled, see (A.13)-(A.15). It is exactly this coupling that makes CESE nonvanish, as opposed to the probe limit [45,46]. Since the final solutions in the vector sector are very lengthy, we record their near boundary expansions only, The tensor sector is the most simple one Now it is direct to compute the boundary stress-energy tensor and currents by substituting where the projection tensor P µν and shear tensor σ µν are defined as E is the energy density and P is the pressure of the fluid, which satisfy E = 3P and The shear viscosity η and bulk viscosity ζ are In what follows we outline the strategy of implementing this transformation but defer technical details to Appendix A.2.
The vector and axial charge densities are promoted versions of (2.20): In the fluid/gravity correspondence, chemical potentials are defined as Note independent transport coefficients are σ e and σ 5e , which are "quantum critical" or "incoherent" conductivities of hydrodynamics as in [48][49][50]. This is partly due to the timereversal symmetry of our holographic model. Physically, bulk quantities in (3.34)(3.35) should be eliminated in favor of thermodynamic variables of the boundary fluid. There are several ways of presenting the results. When discussing single charge limit or probe limit, we find it more transparent to split the conductivities (3.34)(3.35) into two parts, represented by σ Q and σ 0 , as flashed in (1.8)(1.9). On the other hand, in comparison with relevant works [28][29][30][31], we express σ e and σ 5e as functions of chemical potentials µ, µ 5 and fluid temperature T , as summarized in (1.11). We also recast the hydrodynamic constitutive relations (1.3)(1.4) into linear response form, see (1.15) (1.16). Then, we give a clarification for the differences between physical conductivities directly read off from Kubo formulas and intrinsic ones parameterizing hydrodynamic constitutive relations.

Holographic AC Conductivities from Linear Response Analysis
As a complementary study of Section 3, in this section, we reveal some transport phenomena of the holographic model (2.1) through linear response analysis. We focus on the vector and axial currents generated by the external vector and axial electric fields and thermal gradient. We assume these external sources are weak in amplitudes and oscillate in time only.

Black Brane Fluctuations and Conductivity Matrix
To study linear response transports within the holographic framework, we follow standard procedure and perturb the homogeneous black brane (2. 16) g M N = g Black brane fluctuations (4.1) could be classified into decoupled sectors [80] according to their transformation properties under the remaining symmetry group SO(3). For the purpose of computing electrical and thermal conductivities, we consider the helicity one sector only where is a formal parameter marking the linearization. The calculations below will be accurate up to O( 1 ), as required for linear response study.
Fluctuations (4.3) satisfy a system of partial differential equations, whose derivation is presented in Appendix B, see (B.8)-(B.12). While there is no explicit interaction term between A andÃ in the bulk action (2.1), fluctuations δA i and δÃ i do interact via δg ti . As in Section 3, this is exactly our point that the CESE can be realized by going beyond probe limit in a simple holographic model.
We proceed by deriving the compact forms of stress-energy tensor and currents of the boundary theory. To this end, we solve (B.8)-(B.11) near conformal boundary r = ∞, where the non-normalizable modes δh The rest normalizable modes δa (2) i , δã ti , (4.10) where the constraint relation (4.7) has been used to simplify T ti . The heat current Q i is defined as where T ti could be computed from T ti by a weakly curved boundary metric η µν + δh ti . Thus, the currents (4.10)-(4.12), as the linear response to external sources E i , E 5 i , ∇ i T , can be summarized in a compact matrix form where (B.4)-(B.6) have been used. The electrical conductivities are . (4.14) The thermoelectric conductivities are Finally, the thermal conductivity is As seen from (4.14)(4.15)(4.16), only σ and σ 5 are independent:σ could be extracted from σ via the exchange µ ↔ µ 5 ; all remaining conductivities are determined by their combinations. σ andσ are the Ohmic electrical conductivities for vector current and axial current, and σ 5 corresponds to the chiral electric separation effect whileσ 5 is its vector analogue. The CESE conductivity σ 5 is invariant under the exchange of µ and µ 5 . In (4.15), α and α 5 are the thermoelectric conductivities of generating vector and axial currents, respectively. Thanks to time reversal symmetry, Onsager reciprocal relationsᾱ = α andᾱ 5 = α 5 do hold, and the conductivity matrix in (4.13) is symmetric. In (4.16),κ is the heat conductivity. In the numerator ofκ, we have made the replacement E → E + P = 4M as done in [81][82][83]. This added P is actually a "contact term" [84] due to translation invariance. Then, it is consistent with (1.21) obtained by fluid/gravity calculations.
Note the appearance of i/ω in the conductivities (4.14)(4.15)(4.16). From the Kramers-Kronig relation [63,64], this means there must be a delta function δ(ω) in real parts of all the conductivities. However, in holographic models with spatial translational invariance, it is not easy to track the delta-peak directly. For consistency, we just added this δ(ω) as underlined terms above [63,64].

AC Conductivities: Numerical Plots
In this subsection, we present numerical results for the AC conductivities σ and σ 5 while depositing more technical details in Appendix B. Our results for frequency dependence of Ohmic conductivity σ are plotted in Figures 3, 4, 5.   Figure 3 is about the plot of σ/(πT ) as a function of dimensionless frequencyω when eitherμ = 0 orμ 5 = 0. As seen from left panels of Figure 3, whenμ = 0, i.e., neglecting the back-reaction effect from the vector field in the bulk, the low-frequency limit of σ is always finite, which is in consistent with the analytical result (4.17). For representative values of µ 5 chosen by us, increasing the axial chemical potentialμ 5 results in reasonably profound modification for the infrared behaviors of σ. From the right-bottom panel of Figure 3, it is clear that, onceμ = 0 (i.e. beyond probe limit), there appears divergent behavior for imaginary part of σ: Im[σ] ∼ω −1 , which is also in agreement with analytical result in (4.17).  Via Kramers-Kronig relation, this 1/ω-behavior in Im[σ] means that Re[σ] ∼ δ(ω). Moreover, the strength of back-reaction due to vector field in the bulk is also reflected by different curves in the right-bottom panel in Figure 3. Figure 4 is to further explore the effects of chemical potentials on σ. In the infrared regime of ω, Im[σ] is dominated by the divergent behavior ∼ 1/ω, which implies a delta-peak δ(ω) for Re [σ]. While we could not display this δ(ω) in numerical plots, we find the finite piece in Re[σ] and particularly that lim ω→0 Re[σ] = σ e . Let us focus on the infrared regime (roughly with 0 <ω 2) of Re[σ]. Our observation is that increasingμ will diminish Re[σ] while increasingμ 5 will result in enhancement of Re[σ]. This is exactly consistent with the DC limit σ e in (1.11), which has been plotted in the left panel of Figure 1. The behavior for largerω is supposed to be controlled by UV conformal symmetry. To confirm this observation, in Figure 5 we plot Ohmic conductivity σ, as a function ofω, for larger chemical potentials µ andμ 5 . Roughly, whenω 6, Re[σ] becomes insensitive to change of chemical potentials. Figure 6 is to show the frequency-dependence of the CESE conductivity σ 5 and the generalized deviation factor χ A ≡ σ 5 /(μμ 5 πT ) extending (1.11) to AC case. Given that σ 5 is symmetric under the exchange of µ and µ 5 , it is legitimate to constrain to either µ ≥ µ 5 or µ ≤ µ 5 without loss of generality. Just like σ displayed in Figures 3, 4, 5, Im[σ 5 ] shows diverging behavior (i.e.∼ 1/ω) near ω = 0, which, via Kramers-Kronig relation, indicates Re[σ] ∼ δ(ω). Once away from ω = 0, the imaginary parts approach zero soon, while the real parts evolve in a more profound fashion as ω is increased. Particularly, we observe a damped oscillating behavior for Re[σ 5 ], where the asymptotic regime is achieved aroundω = 5.5 for µ = 0. This is roughly in agreement with the numerical results of [31], although we have the different mechanism of generating CESE. Moreover, from Figure 6 we observe that increasing µ orμ 5 would delay the achievement of the asymptotic regime.  Figure 6. Left: CESE conductivity σ 5 /(πT ) (Real and Imaginary parts) as a function of dimensionless frequencyω = ω/(πT ). Right: The deviation factor σ 5 /(μμ 5 πT ) (Real and Imaginary parts) as a function ofω. As in Figure 3, the divergent behavior (∼ 1/ω) in Im[σ 5 ] implies the presence of a delta-peak atω = 0 in Re[σ 5 ], which we could not display numerically.

Conclusion and Discussions
In this work, we explored transport properties of strongly coupled matter, which is holographically described by (4 + 1)-dimensional Einstein gravity, coupled to U (1) V × U (1) A gauge fields, in the asymptotic AdS 5 black brane. Our main finding is the nonzero CESE conductivity when the gravitational back-reaction effect is taken into account, see (1.11) for the hydrodynamic limit and (4.14) for its extension to an AC conductivity. We confirmed our results with two complementary studies -fluid/gravity calculations versus linear response analysis. Within the former framework, we constructed the first-order constitutive relations for stress-energy tensor, vector and axial currents for the holographic matter in the long wavelength and low frequency limit. Following the linear response approach, we revealed the frequency-dependence of Ohmic, CESE and thermoelectric conductivities.
As the second task, we clarified the relations between the dissipative transport coefficients in the hydrodynamic constitutive relations (1.3)(1.4) and those appearing in the conductivity matrix (4.13). While the "intrinsic" conductivities σ e , σ 5e etc. are widely used in the framework of fluid dynamics, the physical observable are indeed those appearing in the conductivity matrix (4.13). Indeed, when the hydrodynamic description is reformulated into the linear response form, we find perfect agreement between these two different approaches, see (1.17)(1.18) (1.20) and (4.17).
Since we have turned off the possible Chern-Simons terms in our holographic model (2.1), it would be interesting to check if the CESE conductivity will be corrected by the chiral anomaly. From recent works [45,46,86], the anomalous corrections to normal transport coefficients start from the second order in the derivative expansion. A study along the line of [45,46,86] will be helpful in clarifying this issue.

A Technical Details in the Fluid/Gravity Calculations
In this appendix, we collect some computational details and useful relations in the fluid/gravity calculations of Section 3.

A.1 Solving the Bulk Equations
Following the standard procedure of fluid/gravity correspondence, we first solve constraint equations to derive relations among fluid-dynamical variables. Practically, we find it more convenient to consider certain combinations of constraint and dynamical equations. Below is the listing of solutions to constraint equations, The asymptotic condition (3.10) requires C h 1 = 0. By Landau-Lifshitz frame condition (3.14), the integration constant C h 2 will also be fixed to zero. Therefore, h(r) = 0. From the time components of Maxwell equations W t = 0 andW t = 0, where C a 1 , C a 2 , Cã 1 , Cã 2 are integration constants to be determined by boundary conditions (3.10)-(3.12). First, nonzero C a 1 , Cã 1 correspond to non-normalizable modes and would violate the asymptotic requirement (3.10). So, C a 1 = Cã 1 = 0. The Landau-Lifshitz frame conditions (3.15) require C a 2 , Cã 2 to vanish, i.e. C a 2 = Cã 2 = 0. The remaining equation of the scalar sector is g rr E rr + g rt E tr = 0, which is solved by A nonzero C k 1 would cause the T tt computed from (3.16) to be in contradiction with the Landau-Lifshitz frame convention (3.14). So, C k 1 = 0. Therefore, all the integration constants in the scalar sector have to be set to zero. The solutions in the scalar sector are summarized as below k(r) = 2 3 II. Vector Sector. -Now we consider the helicity one sector, which consists of a i (r),ã i (r), j i (r) and turns out to be more involved. First consider the Maxwell equations W i = 0 andW i = 0 which are dynamical equations for a i (r) andã i (r), but coupled to j i (r). The Einstein equation W ri = 0 corresponds to dynamical equation for j i (r), which couples to a i (r) andã i (r). Our strategy of solving the coupled differential equations (A.13)-(A.15) is to get rid of j i (r) and derive decoupled differential equations for suitably combined variables from a i (r),ã i (r). To this end, we first integrate over r once in the equation (A.15). As a result, where the integration constant is fixed by the Landau-Lifshitz frame convention, i.e., in the near boundary expansion for j i (r) the constant should be zero in order to be consistent with (3.14). The combinations Q 0 × (A.13) +Q 0 × (A.14), andQ 0 × (A.13) − Q 0 × (A.14) give rise to Then, substituting (A.16) into (A.17) yields The equation (A.18) can be solved via direct integration over r, where the lower limit of the inner integral is fixed by regularity (3.11) at the unperturbed horizon r = r h and the upper limit of the outer integral is fixed by asymptotic requirement (3.10). However, the equation (A. 19) is more complicated and cannot be solved by directly integrating over r. Indeed, the homogeneous version of (A.19) has two linearly independent solutions given by H 1 (r)=r 5 f 0 (r) and H 2 (r)=H 1 (r) ∞ r dr r 3 f 0 (r)H 1 (r) 2 −1 , where the second one H 2 (r) is obtained by the Liouville formula. Then, one could proceed to solve (A.19) by using the method of variation of parameters. In practical calculations, we make a coordinate transformation by u = r h /r and perform the solving of (A. 19) in Mathematica. The regularity at the unperturbed horizon and asymptotic requirement completely fix both integration constants. Since the final solution looks quite complicated, we only record the large r behavior for Q 0 a i (r) +Q 0ãi (r), Therefore, the near boundary behaviors for a i (r) andã i (r) are With large r behaviors of a i (r),ã i (r) at hand, the equation (A.16) could be solved near the boundary r = ∞, yielding III. Tensor Sector. -Finally, the tensor equation W ij − 1 3 δ ij W kk = 0 gives the dynamical equation for α ij (r): which can be solved by direct integration over r. The final solution for α ij (r) is

A.2 Useful Relations in Deriving Dual Currents
In this appendix, we summarize some formulas that are quite lengthy but useful towards deriving non-anomalous parts of the currents in (1. where ∂ t u i will be replaced via the constraint relation (A.4) Meanwhile, the derivative ∂ i M would be replaced by which could be inverted to express ∂ i r h , ∂ i Q and ∂ iQ in terms of ∂ i T , ∂ i (µ/T ) and ∂ i (µ 5 /T ).

(B.32)
We numerically solve (B.32) within the pseudospectral collocation method. The boundary conditions for S + , S − could be straightforwardly obtained from those for S 1 , S 2 . For the sake of performing numerical calculations within the spectral collocation method, we summarize the boundary conditions as equalities