Holographic quenches and anomalous transport

We study the response of the chiral magnetic effect due to continuous quenches induced by time dependent electric fields within holography. Concretely, we consider a holographic model with dual chiral anomaly and compute the electric current parallel to a constant, homogeneous magnetic field and a time dependent electric field in the probe approximation. We explicitly solve the PDEs by means of pseudospectral methods in spatial and time directions and study the transition to an universal"fast"quench response. Moreover, we compute the amplitudes, i.e.,~residues of the quasi normal modes, by solving the (ODE) Laplace transformed equations. We investigate the possibility of considering the asymptotic growth rate of the amplitudes as a well defined notion of initial time scale for linearized systems. Finally, we highlight the existence of Landau level resonances in the electrical conductivity parallel to a magnetic field at finite frequency and show explicitly that these only appear in presence of the anomaly. We show that the existence of these resonances induces, among others, a long-lived AC electric current once the electric field is switched off.


Introduction
In the last decade, there has been an increasing interest in macroscopic anomaly-driven effects. To date, several macroscopic phenomena have been shown to be related to the existence of quantum anomalies in the underlying microscopic theory. This includes, among others, the chiral magnetic (vortical) effect [1,2] which consists of the generation of an electromagnetic current in the direction of an external magnetic field (vortex), a massless propagating mode known as the chiral magnetic wave [3], and the enhancement of the electric conductivity in presence of parallel electromagnetic fields known as negative magnetoresistance [4].
From the phenomenological point of view, the need for massless fermions strongly restricts the variety of experiments where these effects might be verified. To date the two main test grounds for anomalous transport are the state of matter known as quark gluon plasma (QGP) and Weyl/Dirac semimetals. In heavy ion collisions (HIC) the QGP is generated for ∼ 10 −24 s. In HIC very strong and transient magnetic fields are generated. Due to the high temperatures reached in the QGP, it is reasonable to neglect the quark masses. Moreover, experimental evidence [5] shows that the QGP is, despite the high temperatures, strongly coupled.
The second kind of systems in consideration is of a very different nature. Weyl (Dirac) semimetals are solid state crystals. These present effective Weyl (Dirac) fermionic degrees of freedom in concrete points of the Brioullin zone known as Weyl nodes. The main difference between Dirac and Weyl semimetals is that in the former the nodes for opposite chiralities coincide in momentum space while in the latter these are separated. Remarkably, one of the aforementioned effects, the "negative magnetoresistance", has already been measured in several Weyl and Dirac semimetals [6][7][8][9]. Although it is usually said that Weyl semimetals are the 3D graphene, it is not yet clear whether the former are strongly coupled too. Therefore it is of special interest to study these systems from a strong coupling perspective. In this sense the gauge/gravity duality [10] (for textbooks see [11][12][13]) is a specially adequate theoretical tool to study these processes.
A complete study of anomalous transport effects in both HIC and condensed matter systems requires to address the question of how the system behaves out of equilibrium. Several studies in this direction include computations for fixed axial charge and dynamic magnetic fields at weak coupling to leading order in α s [14,15] and at strong coupling via holography [16,17]. The fate of the chiral magnetic effect (CME) in presence of timedependent external electric and magnetic fields has been studied at weak coupling for concrete configurations of the electromagnetic fields [18,19]. Moreover, the authors of [20] studied the evolution of the CME during thermalization in holography. It is desirable to have a complete description of anomalous transport for realistic systems.
Most of the past effort in this direction had its focus on the phenomenology of the QGP in HIC. However, the appearance of Weyl/Dirac semimetals has affected the order of relevance of the different features to be addressed. An example of this is the axial charge time-dependence. In the field of QGP and HIC it has been argued [20][21][22] that the estimated decay rate of axial charge is of the order of the lifetime of the plasma and therefore as a first approximation the charge can be taken to be constant. In "chiral" condensed matter system, however, the story is quite different.
Indeed, it is known that in these systems, the chiral magnetic current vanishes in equilibrium even at very large temperatures [23]. This is due to the fact that the Weyl nodes are filled up to the same energy (the Fermi energy) and that the Nielsen-Ninomiya theorem [24] applies. For materials with the Weyl nodes placed at the same energy, one can generate a net chiral magnetic current by inducing a finite axial chemical potential. This is achieved by the non-conservation of the axial charge density in the presence of parallel electric and magnetic fields. Therefore, the mechanism studied here can be used to generate the CME in Weyl semimetals. 1 Motivated by these ideas, we focus in this paper on the time evolution of the CME in presence of time dependent axial charge. In order to introduce axial charge in the system in a time dependent fashion, we consider the anomaly expression in presence of external electromagnetic fields and in absence of axial currentṡ (1.1) This equation tells us that the evolution of the axial charge is determined by the time dependence of E · B. Concretely, we will consider a static magnetic field and a dynamic electric field.
The possible strongly coupled nature of Dirac/Weyl semimetals motivates us to apply holography. This is not, however, the only reason to consider holography for this topic. As extensively discussed in the literature, holography is nicely suited to perform real time computations which, by means of the holographic dictionary, reduce to solving the bulk system of PDEs. As shown later, our case can be reduced to one hyperbolic PDE describing a quenched current operator in presence of a magnetic field. This connects the initial motivation of this work to the topic of (quantum) quenches in holography.
In recent years, there has been a lot of interest in studying (quantum) quenches due to new experimental results. From the theoretical point of view, we have to compute the (quantum) evolution of an isolated system in the presence of a time-dependent parameter in the Hamiltonian. As explained above, for strongly coupled quantum field theories, gauge/gravity dualities are an ideal playground to tackle such questions. One remarkable outcome of previous studies of quenches in holography is the appearance of universal scaling behaviours in the "early time" response [26] for fast quenches. We explore these ideas within our system and determine how the anomaly affects them.
In order to study all these features, it is in principle necessary to solve the initialboundary value problem i.e. solving PDEs. Indeed we explicitly solve the equations by means of pseudospectral methods. However, for linear systems there is a way out of this. Based on several works from the numerical GR community in flat space [27][28][29], we are able to numerically compute the residue of the QNMs of the system by just solving the ODEs that result from Laplace transforming the PDEs. In other words, one does not need to solve the explicit time dependent problem to obtain the residues. As shown in [29], generically, this information allows to determine the full response only from an initial time scale τ on. τ is obtained from the amplitudes and therefore depends not only on the system under consideration but on the initial and boundary conditions as well. When this applies, it means that the knowledge of all QNM and their residues is not enough to determine the response of the system at initial times. We explore the applicability of these ideas in AdS and give a physical interpretation to τ from the dual theory point of view.
The paper is structured as follows. In section 2, we present and discuss our holographic model and the ansatz that we use for the different fields. This section includes some discussion regarding the 1-point function renormalization, gauge fixing and the definition of the current. In section 3, we adapt the spectral decomposition analysis developed for asymptotically flat Schwarzschild in [29] to AdS-Schwarzschild. Concretely, we explicit compute the amplitudes of the QNMs in our system for different boundary data. In section 4, we focus on the initial time behaviour of the current. First we look at the dependence of τ on the different parameters of the system. After this, following [26], we look at the response of the current for fast quenches and its dependence on the anomaly coefficient and the magnetic field. Finally, in section 5 we study the late time response where we identify QNMs approaching the real axis with increasing magnetic field. This modes are responsible for a long lived oscillatory late time behaviour of the current and for a resonant behaviour under oscillatory sources. By computing the QNMs with backreaction, we argue that this is indeed a consequence of the axial anomaly in our class of holographic models. Moreover, we relate the QNMs to resonances due to Landau levels. We finish this work in section 6 with a discussion of the results and some hints on possible future directions.

The model
In this work we consider a U (1) × U (1) model, which has been previously applied to study anomalous transport in holography [30][31][32]. None of these studies, however, has considered the effect that dynamical generation of axial charge may have in the response. Our aim here is to compute time dependent setups for the CME in holography. Concretely we focus on the evolution of the effect in presence of electric field driven axial charge generation. This is in contrast to the usual "by hand" implementation of net charge via (axial) chemical potential [16,30].
The 5-dimensional matter Lagrangian consists of two photon fields (A µ , V µ ) coupled by a Chern-Simons (CS) term in the bulk withε ijkl being the epsilon symbol in the boundary theory and Latin indices running from 0 to 3. These Ward identities correspond to an abelian anomaly. The Bardeen counter term has been chosen such that only one of the currents is not conserved. It is therefore customary to identify J V and J A with the dual vector and axial currents, respectively. This allows us to meaningfully consider external "electromagnetic" fields and study the evolution of the system in such backgrounds. Moreover, the relative factor between both CS terms is chosen so that the relative factor in the r.h.s. of equation (2.2) coincides with that of the abelian part of the anomaly in QCD [33]. As a first step towards solving the full non-linear system we consider the probe limit of the system. This limits the validity of the results. We comment on this later, when the results are discussed. As background metric we consider Schwarzschild-AdS 5 . In infalling Eddington-Finkelstein coordinates the metric reads Note that the conformal boundary is located at ρ = 0 and without loss of generality we have rescaled the coordinates such that the AdS radius is L = 1 and the black-hole horizon is at ρ = 1. The equations of motion for the gauge fields read

Setup
As explained in the introduction, in order to control the dynamic generation of axial charge we need a time dependent E · B. We keep the magnetic field static and just set a time dependent electric field. Therefore, we consider the following ansatz with boundary conditionV 6) and the usual regularity conditions for A v (v, ρ) at the horizon. Here, the time derivate ∂ v is represented by the dot. This corresponds to a time independent magnetic field and time dependent electric field both homogeneous and pointing in the z-direction in the boundary theory. Due to the anomaly, a non-trivial time component of the axial gauge field is needed to get equations for V z and V y consistent with the choices in equation (2.5). The equations for this ansatz read with the notation for the radial derivative ∂ ρ . Integrating equation (2.9) in time one obtains Substituting this result back into equation (2.7), one finds C 1 (ρ) = ρC. We will later see in equation (2.19) that C is just a gauge shift for V z and will fix it to C = 0. We can then reduce the system of PDEs to a single hyperbolic PDE for the z−component of the vector field The asymptotic boundary expansions read

12)
In order to compute the dual correlators we need to specify boundary terms. We have with γ the induced metric on the boundary. The first term is an infinite counter-term fixed by the theory. In addition, for convenience in our numerical calculations, we fix the renormalization scheme by introducing the second (finite) gauge invariant counter term. We make this choice to avoid explicit contributions of the logarithmic coefficient in (2.12) to the one point function which in our scheme is given by where the subscript stands for "consistent" current. From previous studies of the CME in holography it is known that this version of the current can be problematic [16,30] due to the difficulty to distinguish between the chemical potential and the dual coupling associated to A v in certain gauge. It is now well understood [30] that imposing a vanishing A v at the horizon gives rise to a vanishing consistent current in the time independent limit. One way to circumvent these issues is to compute the covariant current instead In Eddington-Finkelstein coordinates, however, regularity does not imply a vanishing gauge field at the horizon. Therefore we are free to choose a gauge in which the field vanishes at the boundary A v (ρ = 0) = A 0 = 0. With this choice the covariant and consistent currents are indistinguishable From now on we stick to this choice and omit the "vector" superscript. Then, our problem reduces to computing the coefficient of the normalisable mode of the V z field. To this end, we must solve the dynamical equation (2.11) with the boundary condition (2.6) and additional initial data. In order to construct the initial data, we assume that our system is initially in equilibrium, i.e. described by a stationary configuration. Mathematically, this assumption implies that all time derivatives of the boundary data V 0 (v) must vanish at v = 0. In practice, this condition is hardly realised. Nonetheless, as we are going to mention afterwards, we choose specific quenches profiles for the function V 0 (v) in such a way that d j dv j V 0 (0) ∼ 0. The time-independent equation (2.11) can be rewritten in terms of a new variable u ≡ ρ 2 as an inhomogeneous Legendre equation with regular solutions Here, P l are the l-th Legendre functions of the first kind. The constant C is just a shift in the field. It becomes clear that A v in equation (2.10) is C independent. We will from now on gauge fix it to C = 0. The constant C 2 determines the initial value of the current and it can be related to the initial value of the axial chemical potential Since the adiabatic regime implies that the equilibrium formula is valid at any time, the adiabatic response is given by the time independent solution. The asymptotic expansion of the stationary solutions in equation (2.19) read with Γ the Euler gamma functions. Therefore, in the adiabatic limit Moreover, for κB 1 it asymptotes as It is possible to understand the above equation in terms of the effect of the Chiral Magnetic Wave. If a system displays a massless mode with a finite residue R 0 , the response to an externally applied electric field In our situation, R 0 = −12κB for large κB and thus at low frequencies the D.C. conductivity scales as In order to compute the current outside the adiabatic limit, we make use of the asymptotic expansion (2.12) and introduce an auxiliary field U (v, ρ) via Substituting equation (2.24) into the original equation (2.11) we obtain (2.25) with λ = 12 κB. By doing this, we automatically incorporate the information about the boundary condition (2.6) into the inhomogeneity S(v, ρ) given by

26) and
Therefore, equation (2.25) is to be solved as an initial value problem with The most natural way to solve (2.25) is to numerically integrate it in time. In this work, the time evolution is performed with the fully spectral code 2 introduced in [34] (the numerical details are described in appendix A). We use this highly accurate numerical method to obtain the full response of the current operator. In particular, it provides us with a reference solution, against which we can compare the results from the strategy based on the spectral decomposition of the solution U (v, ρ) in terms of the quasi-normal modes.

QNM amplitudes from Laplace analysis in AdS
As stated in the introduction, solving the PDE is not the only way to determine the response of the current. Given the spectral decomposition for the dual operator solving the problem reduces to obtaining the eigenfrequencies ω n and the amplitudes A n for given initial and boundary data. It is well known how to obtain the quasi-normal frequencies in holography [35]. However this is not the case for the amplitudes A n . Of course a possibility is to fit the data obtained from the PDE to the previous formula (see for example [36]). It is, nevertheless, possible to compute A n directly without need of explicitly solving the PDEs. Although this topic has a long history in the community of numerical general relativity [27,28] only recently [29] has it been shown how to compute these amplitudes in the case of asymptotically flat Schwarzschild black hole for rather generic initial data. More specifically, the authors formulate the wave-equation describing the propagation of fields on the Schwarzschild background in terms of a particular coordinate system, where the surfaces of constant time are spacelike hypersurfaces that penetrate the black-hole horizon and extend to future null infinity. By using the standard framework provided by the Laplace-transformation, they develop a semi-analytical algorithm to obtain "eigenvalues" and "eigenvectors" (related only to the wave-equation in question) and amplitudes (related to the particular initial data being used). Most importantly, they introduce a well defined time scale τ for which the spectral decomposition in the form (3.1) is valid, based on the growth rate of such amplitudes.
In this section, we provide the first steps towards the application of these techniques in asymptotically AdS and we focus on the role that the τ may have in the dual theory.

Laplace transformation
We consider the compact form of the dynamical equation (2.25) for the function U (v, ρ) with given initial data U in (ρ) specified in (2.28). Here, α(ρ) and β(ρ) are differential operators acting on the radial coordinate ρ, which can be easily read from (2.25). The former is second order, whereas the latter is first order. As mentioned, the source term S(v, ρ) contains information about the boundary data V 0 (v) and is generically written in the form Here, N S = 4 and the functions a i (ρ) are summarised in equations (2.27). Applying the Laplace transformation 3Ū and taking into account that for a generic function f (v, ρ) As observed in the previous section, ideally the boundary data is such that all the time derivatives vanish initially. A condition that is only approximatively realised in practice.

Quasi-normal amplitudes
Quasi-normal modes are the complex s n −values for which φ n (ρ) is a regular solution to the homogeneous equation This equation is solved in Mathematica within the framework of the "Generalised Eigenvalue problem" (see appendix A.1 for further details). As discussed in [29], the functionŪ (ρ; s) has poles on the quasi-normal modes. We introduce the decompositionŪ Besides, we identify the differential operator on the l.h.s of equation (3.6) as A(s) = α+sβ and rewrite it in the form Then, equation (3.6) in the limit s → s n gives as regularity condition At second order, we obtain equation (3.12) is to be solved simultaneously forW (ρ) and η n with the normalization conditionW (ρ 0 ) =W 0 (for further details see appendix A.2).
Then, the standard recipe to obtain the spectral representation of the solution U (v, ρ) is as follows: 1. one expresses the time-dependent solution U (v, ρ) as the inverse Laplace-transformation; 2. the integration path of the inverse Laplace-transformation in the complex s−plane is deformed in order to include the QNMs (see [29] for further details). By doing so, we obtain with the last term corresponding to the integral along the half-circle C as |s| → ∞; 3. for asymptotically flat spacetimes, it is argued in [29] that the validity of (3.13) depends on the behaviour of |η n φ n (ρ)| for large n. It is shown that generically |η n φ n (ρ)| ∼ e −τ Re(sn) . For v < τ , the functionÛ (ρ; s) diverges exponentially as |s| → ∞. On the other hand, for v > τ , the contribution from the integral C is zero and we obtain the desired spectral representation.

Numerical evidences for the recipe in AdS-Schwarzschild
Here, we present numerical evidence showing that the recipe discussed above may also work for asymptotically AdS spacetimes. Using (2.17) we identify 4 A n = −2 η n φ n (0). (3.14) First, for the quenches to be presented in (4.1), we calculate the quasi-normal amplitudes A n using equation (3.12) and we observe that they indeed behave as Here ω n are the Fourier frequencies. We expect 5 this behaviour to be valid for any quench which approaches a stationary configuration as v → ∞. In particular, we show in the left panel of figure 1 an example for the Gaussian quench with κB = 2 and Λ = 0.1, for which we read the time τ = 3.0.
Second, we check that the sum in (3.1) converges for v > τ , based on In figure 1 we also show the comparison of the spectral analysis to the numerical data obtained by solving the PDE for this specific quench. As mentioned before, here we obtain τ = 3.0 and the time evolution based on the spectral decomposition matches the numerical evolution accordingly.
In practice, since we can only take a finite number of quasi-normal modes into account, we can distinguish three regimes in the plot. The first is v < τ where the spectral decomposition diverges and does not resemble the current response at all. There is a second, intermediate, regime τ < v 4 where the fit is not completely accurate. In order to obtain a perfect fit in this region, we should also take into account higher QNMs. Finally for v 4 the spectral decomposition fits perfectly the numerical data.
The regime v < τ requires further investigations. It is still not clear to us whether the inclusion of all higher QNMs and possibly an analytic continuation can cure the deviation between the numerical time evolution and the spectral decomposition observed in figure 1. We leave further investigation in this topic for future work. To summarise, τ defines the time from which an explicit solving of the time evolution can be substituted by an analysis based on the spectral decomposition i.e. on the quasinormal modes of the system. It depends not only on the system under consideration but on the concrete initial and boundary conditions. Moreover, there is in principle no bound on the value of this quantity. This means that for specific systems or initial/boundary conditions it might be possible to get τ < 0.
In [37], it is argued that τ can be understood as the time it takes for (a compact support) initial data to propagate to a point whose lightcone contains all the initial data information. So far, we have not found a similar intuition in AdS. It is worth mentioning that [37] assumes that the initial data is a function of compact support and does not consider the possibility of a boundary with "ingoing" data.
In the following, we investigate the dependence of τ on the width of the quenches and on the effective anomaly-magnetic field parameter κB. With this we would like to explore the possibility of considering τ as a well defined notion of "initial time" for linear (or linearised) systems.

Initial time
In the previous section we reviewed how to compute the residues of the QNMs of the system and discussed the mutual growth rate τ . As emphasised before, this time scale is fixed not only by the system under consideration but by the initial and boundary data as well. In this section we study the dependence of τ on the different quenches and on κB and we explore the possibility of identifying it with the out of equilibrium or "initial response" time of the system.
Concretely we have computed τ as a function of the anomaly parameter κB for quenches of the form where v i fixes the "centre of the source", and of Λ, which is inverse to the width of the signal. This quantity should be though of as the abruptness of the quench, so that Λ 1 implies a very fast quench. Let us remark that the main qualitative difference between the two sources is that Gaussian sources introduce no final net axial charge in the system while hyperbolic tangents do. This can be easily understood by considering that the total axial charge introduced in the system goes like 6 As in previous studies of holographic quenches [26], we expect the quench scale Λ to play an important role in the response of the system. For non-compact support functions like the ones we are considering, it is necessary to set some reference time to compare τ with. The most natural quantity to consider in our case is the temporal distance to the centre of the source: τ − v i . In figures 2 and 3 we show the dependence of τ − v i on Λ for several values of κB for Gaussian and tanh quenches respectively. For Gaussian quenches (see figure 2) we find a linear dependence on Λ −1 . In particular, for κB = 0 we observe no dependence on this quench parameter. In the r.h.s. of figure 2 we show the dependence of the slope (τ − v i )/Λ −1 with κB. As shown in the figure the data fits to a line with slope −4κB. The fact that this slope is negative means that for the same quench, the initial time is smaller the higher κB. Hence, the system is closer to the adiabatic response the stronger the anomaly term. In other words, the relaxation time of the system is smaller for larger κB.
For tanh quenches (see figure 3) we find again a linear dependence on Λ −1 for small Λ (Λ 10). In this case, contrary to what we observed in the Gaussian quench, κB only shifts the lines and does not affect the slope. For larger Λ, the time scale τ − v i shows a non-linear dependence on Λ −1 . Although qualitatively different, the behaviour of τ for tanh quenches is still compatible with the intuition of a smaller relaxation times the higher κB.
Summarising, in all cases the data fits to a linear dependence on Λ −1 in the regime of small enough Λ. That τ − v i is generically negative and increases for faster quenches seems natural, since for very slow quenches one gets closer to the adiabatic regime and one expects to find τ − v i → −∞. By looking at a different quantity, obtained within the initial time v < τ , we will confirm this intuition in next section.

Universality in fast quenches
In the previous subsection we have seen how τ depends on Λ and κB. As already argued this provides us with some notion of "initial" response of the system. In this subsection we investigate the initial response regime defined as v < τ . Let us recall that for this part of the response one has to, by definition, explicitly integrate equation (2.25) in time.
Our basic motivation comes from the studies of quantum quenches in [26,38]. In these papers a universal response was found for certain kind of smooth, fast enough quenches in generic CFT's. Concretely the authors of [26] looked at a free scalar holographic model. Among others, it was shown that the time τ ex , at which the operator deviates 5% from the adiabatic response, follows a simple universal law for fast quenches. We perform an analogous analysis in our system.
Before we proceed to explain the details, some comments are in order. A first observation is that our system is linear, in contrast with that of [26]. As we will see this does not prevent the appearance of an universal behaviour for fast quenches. In addition to this, we emphasise that our main interest is to determine how the anomaly affects the universal behaviour. Therefore we focus on the dependence of the universality on the value of our anomaly/magnetic parameter κB. Moreover, we focus on Gaussian sources henceforth. To give a qualitative idea of the system response, we show the generic behaviour of the current in figure 4. In the l.h.s. we fix κB and plot the current for two qualitatively different values of the width Λ of Gaussian. The qualitative difference between a "fast" quench (blue) and a "slow" quench (orange) is apparent: in the slow case the current approaches the adiabatic behaviour, given by equation (2.22), depicted in black. The fast quench shows a more complicated structure and deviates significantly from its corresponding adiabatic behaviour. In the r.h.s. we fix the width and show the response for different values of the magnetic field. Again, one can observe a qualitative difference; the lower κB the bigger the deviation from the adiabatic response. This simple qualitative observations seem consistent with those made in the previous section. We would like now to see whether the system shows an universal behaviour for fast quenches. This will shed some light concerning the relaxation time. In order to characterise this we have found it useful to look at the following quantity with J z (v 1 ) and J z (v 2 ) the first minima/maxima of the current, respectively. We have computed δ as a function of Λ for several values of the κB. Our results are summarised in figure 5. The plot shows that fast quenches indeed show an universal, B independent, behaviour. Concretely we find and a fitting error for α of the order 10 −3 . As one can observe in figure 5, as we increase κB, we need faster quenches, i.e. higher Λ, to get to the universal regime. Despite of the smoothness of the transition, it is possible to define a critical Λ using an interpolating function for the (logarithmic) data. It is natural to identify Λ C with the absolute maximum of the second derivative of the interpolating function. We get same results using either splines or Hermite polynomials as interpolants and for several interpolation orders. In the left panel of figure 5 we have exaggerated the corresponding points. In the right panel we show the values of Λ C against κB. By fitting this data to a curve of the form a + bx c we find a = −0.08, b = 4.33 and c = 0.96 with fitting errors 0.19, 0.19, 0.07 respectively. We conclude that the transition to the "universal", fast quench, regime is delayed linearly with κB. This suggests that the effective relaxation time of the system is inversely proportional to κB. This is in agreement with the τ behaviour for tanh quenches displayed in figure 3 where we see that the nonlinear ("fast") behaviour appears for higher Λ the higher κB.

Landau resonances
Let us now look at the response of the current for late times. Concretely, we consider times bigger than the transition time τ as defined in equation (3.15). An interesting result in this regime is depicted in figure 6. Here we show the behaviour of the current for a fixed Gaussian source and several values of the magnetic field. The decay rate of the late time oscillations decreases as we increase κB. In particular, for κB 1 this relaxation time is very small compared with all other scales in the system. This signals the existence of a QNM approaching the real axis for increasing κB. We have computed the late time current for a variety of quenches (Gaussian, tanh, oscillatory) with analogous behaviour; the only qualitative difference we have found is the value around which the current oscillates and that is given by the total amount of axial charge induced in the system (see equation (4.2)). There is, in fact, an analytic argument to see that indeed a QNM must approach the real axis with increasing κB by writing the perturbation equation for the relevant field component as a Schrödinger-type equation (see appendix B for further discussion).
Next we want to check whether this behaviour is indeed an anomaly driven phenomenon, at least for a class of holographic models. The main obstacle we face is that in the probe limit the magnetic field always appears in combination with the Chern-Simons coupling κ in the equations. Therefore, in this limit it is not possible to disentangle anomaly from magnetic field effects. Moreover, the probe approximation only makes sense for small values of the fields and we are now interested in the fate of the QNMs as we increase B. This motivates us to compute the QNM spectrum of the system including backreaction. Details on the computation can be found in appendix C. Our main goal is to address the following questions: How do the lowest QNMs behave for increasing B? Is any mode crossing to the upper half plane Im[ω] > 0? What is the role of the axial anomaly in the dynamics of the QNMs?
We have looked at the dependence of the lowest QNMs on κ and B at zero momentum. The main message is summarized in the l.h.s. of figure 7. Here we show the behaviour of the three lowest QNMs as one increases the magnetic field for κ = 0 and κ = 1. The difference is apparent. For κ = 1 the modes approach the real axis and the real part increases with B. Conversely, if we set κ = 0 the modes approach the imaginary axis. There is an intermediate region, for approximately 0 < κ 1/2. To illustrate this we show the fate of the lowest QNM for κ ∈ {0.1, 0.2, 0.34, 0.5, 1} in the r.h.s. of figure 7. Here we see that in the low κ regime the imaginary part of the frequency begins decreasing again once the magnetic field is large enough. We found similar results for higher QNMs. This implies that for low values of κ no resonances are to be found, independently of the value of B. For κ 1/2 we find the mode approaching monotonically the real axis in the regime of B allowed by our numerics. For these values the QNMs approach the axis up to very small absolute values (Im[ω] ∼ −10 −6 for κ = 1/2) of the imaginary part and resonances are found. In all cases the approach is faster the higher κ. We didn't find any mode crossing to the upper half plane for any of the values considered.
In addition to this, it is clarifying to look at the behaviour of the real part of the frequency. By fitting it to a function of the form Let us make a remark. In principle, for large values of κ we cannot exclude the possibility of the modes going away from the real axis again for very large values of B. If this is the case, then the resonances would be restricted to a certain regime of the magnetic field.
Next, we inquire the dependence of b (see equation (5.1)) on κ. For the data used in figure 8 we find b ∼ κ 0.510 with a fitting error of 10 −3 for the exponent. Therefore, for high enough κ and large enough magnetic field we conjecture that Re[ω] ∼ √ κB.
Finally, the slow decay rate of the current is not the only effect that "almost normal" modes have in the current. One can consider the possibility of directly exciting these modes by switching on an oscillatory source with appropriate frequency. In figure 9 we show the behaviour of the current within the probe limit for a source of the form It corresponds to an oscillatory electric field with an amplitude envelope 7 chosen for numerical convenience. As expected, the current response is qualitatively similar to that of a damped driven oscillator. For frequencies close to the resonant frequency ω c (fixed by κB) the amplitude is dramatically enhanced and several characteristic frequencies appear. At the critical value of ω the amplitude is maximally amplified and it asymptotes a certain value. The higher κB the stronger the resonance.

Discussion
In this paper we have considered quenches in an holographic system with a U (1) × U (1) anomaly. Concretely, we looked at the electric current generated parallel to a magnetic field known as chiral magnetic current and its response to a parallel time-dependent external electric field. Our motivations were twofold. On the one hand, as argued in the introduction, it is necessary to address the question of how anomalous transport behaves out of equilibrium and how it reacts to the dynamical evolution of the axial charge. On the other hand, we were interested in how the axial anomaly affects previous results [26] regarding the universal response to fast quenches. As a first step, we have only considered the probe limit of the theory. The probe approximation limits the validity of the results to field configurations with a small stress-energy tensor compared to the temperature of the black hole. As a consequence, this approximation is only justified for sufficiently "small" magnetic fields and the sufficiently "slow" quenches. We have already made progress in the task of including backreaction into the system and will present our results in a follow up paper. However, we remark that one of our main results, the Landau resonances for large magnetic fields, were obtained with backreaction.
In the probe limit, our setup yields in a linear hyperbolic PDE, which has to be solved as a initial/boundary value problem. We explicitly evolve this equation in time with a fully spectral method [34]. Moreover, we also consider a different strategy to obtain the solution to the problem. Based on the discussion from [27][28][29] we read the so called quasinormal amplitudes out of the initial/boundary data and express the solution as a spectral decomposition based on the quasi-normal modes. As discussed in [29,37] for asymptotically flat spacetimes, this analysis is expected to be valid for a time scale v > τ , where τ is related to the growth rate of the amplitudes. We observed that the same analysis works in our asymptotically AdS system. Moreover, based on [39], we expect this approach to be valid in generic asymptotically AdS spacetimes.
We have investigated the behaviour of τ for Gaussian and tanh quenches with the idea of associating this quantity to a well defined notion of initial time. Our results are positive, although some unsolved questions remain. We find that the qualitative dependence of τ on the abruptness of the quench and κB matches the expectations: τ becomes smaller close to the adiabatic regime or for higher κB.
However, it is still to be investigated how τ depends on the precise quench protocol and whether τ can be used to identify a generic classification for the quenches. In particular, the Gaussian quench showed a rather unexpected behaviour at κB = 0, namely the time scale τ appeared to be independent of the width. This issues deserve deeper analysis and we leave it to future work.
In addition, we have studied the behaviour of the dual current operator for fast quenches (motivated by [26]) where a universal response regime was found in a holographic model with a dual scalar operator. Our results show that for fast quenches the quantity δ defined in equation (4.3) shows a universal behaviour: it is independent of the magnetic field and it scales with a fixed exponent with the width of the quench. Moreover, we have shown that this universal regime is suppressed the higher κB, i.e. one needs to perform faster quenches to reach the universal regime the higher κB.
Our results for both τ and δ indicate that the relaxation time of the system is smaller with increasing κB. This is coherent with the notion of κB being a coupling of the current to the electric field 8 .
The last main point of this work is the emphasis on the dynamics of the current for late times. By explicitly computing the fate of the current after quenches in the electric field, we have shown that oscillatory, long-lived currents are produced in presence of the anomaly and strong enough magnetic fields in our holographic model. By computing the QNMs of the system with backreaction we have been able to show that this is a consequence of the presence of Landau levels 9 . This has been done by looking at the dependence of the lowest QNMs on κ and B. We have observed that for non-zero κ the real part of the frequency goes as ∼ √ B for large B, pointing towards the appearance of Landau levels for κ > 0. As already discussed in the main text, the fact that the characteristic √ B behaviour is absent for κ = 0 indicates that no dual fermionic d.o.f. are charged under the symmetries in consideration when the Chern-Simons term is not present in the bulk, i.e. when there is no dual anomaly.
Whether resonances are to be found in the conductivity depends not only on the magnetic field, but on κ as well. As we saw, the imaginary part of the quasi normal frequencies asymptotes the real axis only if κ 1/2. It is tempting to speculate whether this is related to the existence of a quantum critical point for κ > 1/2 studied in [42][43][44].
In addition, we explored the idea of directly exciting the Landau resonances by means of an oscillatory external electric field finding a resonant pattern as expected 10 . Moreover it would be interesting to check whether long lived oscillating currents can be generated by means of an electric quench in Weyl semimetals. A more realistic treatment of this would require to introduce some ("intrinsic") mechanism for axial charge relaxation [31]. We propose to check this statements in the weak coupling regime. Regarding this resonances, it is worth mentioning that the system shows several similarities to holographic p-wave superconductors. Both present dissipationless transport driven by a broken symmetry and both cases [46] present resonances at finite frequencies in the anisotropic direction. It can be interesting to explore this analogy in more depth.
Let us finally mention that initially we had Weyl/Dirac semimetals in the back of our heads, which is described at strong coupling by [47]. However, we have considered a simpler model which we think captures the core properties we wanted to focus on. It would be interesting to make an analogous study with the holographic model [47]. development in the application of spectral methods along the time direction in dynamical scenarios [34].
Here, we discuss the main features of the method with emphasis on its application for eigenvalue problems and for systems of differential equations with extra parameters. Then, we review the algorithm for the fully spectral code introduced in [34].

A.1 The eigenvalue problem: quasi-normal modes
For convenience, let us reproduce here the general representation of the ordinary differential equation characterising the QNM as we introduced in (3.8) (for a fixed background) or (C.9) (for the backreacted system) This equation is to be solved in the domain ρ ∈ [0, 1]. Due to the choice of the ingoing Eddington-Finkelstein coordinates to the background metric, we obtain only a linear term in the complex parameter s n characterising the QNM 11 . Furthermore, thanks to this coordinate choice, the surfaces of constant time penetrate the black-hole horizon and therefore the ingoing boundary condition on the horizon is automatically realised here by the geometry, i.e., they correspond to the regularity condition of (3.8) at ρ = 1. In the same way, equation (3.8) provides us with the proper regularity conditions at ρ = 0, so no further information is needed for the solution of equation (3.8).
In order to discretise equation (3.8), we fix a numerical resolution N ρ and we introduce a Chebyshev-Lobatto grid in the ρ direction via Then, we recall that α(ρ) and β(ρ) are (second/first order, respectively) differential operators acting on the radial coordinate ρ. Hence, a discrete representationα andβ is obtained by substituting the differential operators ∂ ρ and ∂ 2 ρρ by the discrete Chebyshev-Lobatto spectral differentiation matricesD ρ andD 2 ρρ =D ρ ·D ρ , whose expression can be found in [49]. In its discrete form, equation (3.8) has the structure of a generalised eigenvalue problem and both eigenvalues s n and eigenvectors φ n are easily obtained with Mathematica.
First thing to observe is that we obtain a total of N ρ + 1 eigenvalues and eigenvectors. Yet, not all of them are trustful numerical solution. It is crucial that one performs a convergence test to determine which of the obtained solutions are stable and converge to fixed value as one increases the resolution N ρ . Our empirical observation shows that for a given N ρ , the first n QNM ∼ N ρ QNMs correspond to a reliable numerical solution. In order to keep high-accurate solutions, we fix N ρ = 300 and considered only the first n QNM = 12 quasi-normal modes in this work. Secondly, we would like to stress that we are also interested in the eigenfunctions φ n , since they are needed in the calculation of the QNM amplitudes η n . Each component of φ n corresponds to the value of the function φ(ρ) at the grid point, i.e., φ j n ≡ φ(ρ j ). Notice that the eigenfunctions are uniquely defined up to a normalisation constant. Thus, we work with the conveniently rescale quantity 12 φ norm n = φ n /φ Nρ n which give us φ norm n (0) = 1.

A.2 Equation with a free parameter: QNM amplitudes
We now proceed with the discussion of the solution of equation (3.12) (which we reproduce here once more for convenience) We remind that A(s n ) = α[φ n ] + s n β[φ n ]. Moreover, we have simplified the source term on the r.h.s. with the introduction ofQ = β[U in ] −S.
Our objective is to numerically solve equation (3.12) for both W (ρ) and η n . We also make use of spectral methods for this task so, in principle, the discretisation procedure in terms of the Chebyshev-Lobatto grid points (A.1) and the discrete Chebyshev-Lobatto spectral differentiation matrices is the same as described in A.1.
Note, however, that we have now a total of N ρ + 2 unknown variables, which we collect into a single vector (with W i ≡ W (ρ i )) The discrete version of equation (3.12) provides us only with N ρ + 1 equations, though. The extra condition needed to complete the system is obtained after observing that the solutionW (ρ) is not unique. AssumingW a (ρ) is solution to equation (3.12), then,W b (ρ) = W a (ρ) + Cφ n (ρ) also satisfies (3.12). Therefore, we must impose an extra normalisation condition to W (ρ), which here we fix as Taking into account the discretised spectral representation of equation (3.12) together with the condition (A.3), and introducing the notation β φ =β · φ, we end up with the linear systemM · X = q, whose components can be explicitly expressed as We observe that η n does not depend on the normalisation W (ρ 0 ) = W 0 . While the algorithm used to solve equation (3.12) is essentially the one described above, we must face a caveat introduced by the logarithmic terms presented in the source term Q(ρ) [see equations (2.27) and (3.7)], which leads to an algebraic convergence rate of the spectral scheme. In fact, for generically non-vanishing boundary data V 0 (v), the sourcē Q(ρ) is merely C 0 due to the term ∼ ρ log(ρ) and in practice, one faces difficulty in finding the amplitudes η n already for n 2.
To overcome this problem, we introduce a new coordinate [52] z ∈ [0, 1] via ρ = e 1−1/z , (A.5) which allows one to map the problematic C k−1 terms ∼ ρ k log(ρ) into the C ∞ expressions ∼ (1 − 1/z)e (k−k/z) . By rewriting 13 equation (3.12) in terms of z and discretising the system with Chebyshev-Lobatto grid points z i together with differentiation matricesD z , we obtain a spectral convergence rate 14 for the solution W(z) = W (ρ(z)) and therefore a much more stable scheme for finding the amplitudes η n . A second remark regarding the efficiency of the code is related to the presence of rather huge numbers (∼ 10 30 −10 50 ). Such values are a consequence of the Laplace-transformation of the boundary functionV 0 (s n ), which contributes significantly to the source functionQ(ρ) (see (3.7)). Therefore, it is also convenient to rescale the vectors in (A.4) by q =V 0 (s n ) p and X =V 0 (s n ) Y and solve the systemM · Y = p.
In order to obtain the desired high accuracy for the first n QNM = 12 quasi-normal modes amplitudes η n , we set the resolution to N z = 600. It is true that both resolutions N ρ = 300 (for the solution of the eigenvalue problem -sec. A.1) and N z = 600 are quite extreme for codes based on spectral methods. We mention however, that for a given physical parameter κB the generalised eigenvalue problem and the inversion of the matrixM must be performed only once. The results can be saved and applied afterwards several times for any r.h.s q. Still, a more sophisticated and efficient approach based on the work [29] together with a systematic study of the QNM amplitudes in a more generic context is planned to be presented in a forthcoming article.

A.3 Hyperbolic equation: fully spectral code
We end this section with a brief discussion of the fully spectral code used to solve the equation (2.25) in terms of the coordinates {v, ρ}. A detailed description of the method can be found in [34]. Let us start by assuming we are looking for a solution in a generic time interval v ∈ [v a , v b ]. Given the initial data U a (ρ) at v = v a , we introduce the auxiliary For prescribed numerical resolution N ρ and N v , we work with the Chebyshev-Lobatto grid in the ρ-direction given in equation (A.1) and with the Chebyshev-Radau collocation points in the v-direction The function values V ki = V are stored in a vector X T = (V ki ) k=0...Nv, i=0...Nρ , from which Chebyshev coefficients c kj of the field P (v, ρ) are computed by inverting the equations 13 The point z = 0 must be treated with care. In this limit, one obtains W,z(0) = 0 and this property should be explicitly implemented when constructing the matrixM . 14 Here, spectral convergence means that the convergence rate is faster than algebraic. However, the rate does not decay exponentially since W(z) is not analytic at z = 0.
with T j (µ) = cos[j arccos(µ)] the Chebyshev polynomials of the first kind. After calculating spectral approximations of the fields' derivatives, equation (2.25) yields a linear algebraic system, which is solved with the iterative BiCGStab method. Furthermore, we also provide the BiCGstab method with a pre-conditioner based on a Singly Implicitly Diagonally Runge-Kutta (SDIRK) method [34].
After obtaining the solution U (v, ρ) for v ∈ [v a , v b ], the values U (v b , ρ) at the upper time boundary v b serve as initial data for a subsequent time domain v ∈ [v b , v c ]. This procedure allows to divide the whole time interval v ∈ [0, v final ] into smaller sub-intervals, with the very first one being [0, v a ]. The size of each time interval ∆v can be chosen according to the quench profile. If ∆ is a characteristic time length for a given quench, then we fix ∆v = ∆/4. Furthermore, we set the resolution in the time direction to N v ∼ 25. The radial direction requires a higher resolution (N ρ ∼ 100) due to the presence of logarithmic terms in the source function S(v, ρ). ) for κB = 0 against r. The curve grows monotonically with increasing r. Right: Effective potential for κB = 0.1, 1, 2 (redblue) against r. The κB term becomes relevant close to the horizon and forms a higher barrier for increasing κB.
Imposing this expression to match with (B.1) fixes S and H up to an unimportant constant and leads to The explicit expression effective potential is V eff (κ 2 B 2 ; r) = r 4 − 1 576 (κB) 2 + 3r 4 + 5 4r 6 . (B.7) Note that lim r→∞ V eff = ∞, whereas V eff (r = 1) = 0. We plot this function in Figure 10. As explained in [53] the κB = 0 case is compatible with a conductivity displaying a continuous spectrum due to the fact that the potential is unbounded. However, as κB increases the effective potential develops a barrier close to the horizon, creating more bound states as κB grows bigger. For κB → ∞ the barrier is infinite and we have a binding potential: there is only a discrete set of allowed frequencies and the QNMs do not display an imaginary part because they cannot reach the black hole due to the infinite barrier.

C Quasi Normal Modes
As explained in the main text, in order to be able to vary κ and B independently we have to take backreaction into account (see e.g. [54]). Since the magnetic field breaks rotational invariance we choose the following ansatz for the metric ds 2 = 1 ρ 2 −U (ρ) dv 2 − 2 dv dρ + W (ρ) 2 (dx 2 + dy 2 ) + H(ρ) 2 dz 2 . (C.1) The equations of motion for the background in the so called trace reduced form read R µν = −4 g µν + τ 2 − 1 6 g µν F αβ F αβ + H αβ H αβ + g αβ (F µα F νβ + H µα H νβ ) . (C.2) Due to diffeomorphism invariance we can set the black hole horizon to ρ = 1. The metric function U (ρ) is the blackening factor and hence has to vanish at the horizon U (1) = 0. Before solving the equations of motion we first consider the asymptotic expansions. Imposing asymptotically AdS where we set the term linear in ρ in equation (C.4) (and therefore in equations (C.5) and (C.6)) to zero in order to fix the remaining diffeomorphisms. Furthermore we can fix 2τ 2 = 1 since it appears in the equations always as a product with the magnetic field B.
Using the ansatz V m = V y (x) = Bx and A m = 0 we can solve the equations of motion by a spectral method. To do so we first rescale our functions properly using the expansions equation (C.4)-(C.6) and solve for given κ and B equation (C.2) for the rescaled functions.
With the background solutions at hand we can look at the fluctuations (c.f. [55]). Since the factor proportional to the Chern-Simons coupling is metric independent the metric fluctuations will decouple from the fluctuations of the gauge fields. Furthermore the (1) and (2)  We solve this equation by means of pseudospectral methods imposing the non-normalizable mode to zero. The QNMs we are interested in converge slowly for high values of the magnetic field. To ensure the convergence we improve the spectral solution by a coordinate mapping of the radial variable ρ →ρ 2 which moves the gridpoints in the direction of the boundary localised atρ = 0. Furthermore we use a grid with N ρ = 200.