Nonequilibrium axial charge production in expanding glasma flux tubes

Axial charge production at the early stage of heavy-ion collisions is investigated within the framework of real-time lattice simulations at leading order in QCD coupling. Starting from color glass condensate initial conditions, the time evolution of quantum quark fields under classical color gauge fields is computed on a lattice in longitudinally expanding geometry. We consider simple color charge distributions in Lorentz contracted nuclei that realize flux tube-like configurations of color fields carrying nonzero topological charge after a collision. By employing the Wilson fermion extended to the longitudinally expanding geometry, we demonstrate the realization of the axial anomaly on the real-time lattice.


I. INTRODUCTION
In relativistic heavy-ion collisions, CP-violating configurations of color gauge fields can be generated locally either by gauge field dynamics at the instant of a collision or sphaleron transitions at later times [1][2][3][4][5]. Quarks interacting with such gauge fields induce the imbalance of axial charge due to the quantum phenomenon of axial anomaly. In presence of a strong U(1) magnetic field, which may be generated in off-central collisions, the axial charge asymmetry can be converted to a flow of electric current along the magnetic field [6]. This phenomenon is called chiral magnetic effect (CME) [7][8][9]. Experimental searches for this novel phenomenon have been carried out at RHIC and the LHC [10][11][12], where a charge dependence of azimuthal correlations was measured [13]. However, the observation of the CME in heavy-ion collisions still remains inconclusive due to large backgrounds [14].
On the theory side, there have been numerous developments in the description of the transport phenomena associated with the CME based on the chiral kinetic theories [15][16][17][18][19][20] and the anomalous hydrodynamics [21][22][23][24][25][26]. To make predictions of observable consequences of the CME, some of these frameworks need the information of the axial charge distribution as an initial condition as well as * ntanji@ectstar.eu arXiv:1805.00775v2 [hep-ph] 19 Jul 2018 of the weak-coupling and strong-field approximation, the Yang-Mills equations couple to current induced by the quarks representing the effect of the backreaction [52]. In this study, we consider the weak-coupling and strong-field limit and thus neglect the backreaction.
The aim of this paper is to present formulation and numerical results for the axial charge production in the glasma gauge fields taking the expanding geometry specific to the early stage of heavy-ion collisions into account. To manifest the axial anomaly on a lattice, one has to take care of the fermion doubling problem [53]. In the context of the real-time lattice simulations, the Wilson fermion method has been successfully applied to the description of the axial anomaly in nonexpanding systems [47][48][49][54][55][56][57]. We will employ the Wilson fermion method that is extended to the expanding geometry, which was first introduced in Ref. [50]. Since the glasma gauge fields are produced as a consequence of the interactions between colliding two sheets of CGC, it is important for consistency to take the interactions of the quark fields with the CGC fields into account, i.e. to solve the Dirac equation under the CGC gauge fields. The Dirac equation is analytically solvable until the time right after a collision [43], and the solution that explicitly manifests the boost invariance of the system has been derived in Ref. [44]. We will employ this solution as an initial condition for the time evolution after a collision.
In the framework of the CGC, classical gauge fields are emitted from color charges that represent hard degrees of freedom in a nucleus, and the distributions of the color charges are treated as random variables [30][31][32]. When two nuclei collide, the longitudinal color fields are generated depending on the color charges of each nucleus. Since the two nuclei are causally separated before the collision and their color distributions are random, also the topological charge density F F has random nature: it fluctuates event by event, and in each event it has a random distribution in the transverse plane. As a first step to elucidate the axial charge production in the early stage of heavy-ion collisions, instead of the random color distributions, we consider fixed configurations of the color charges that realize simple flux tube-like configurations of color fields that has nonzero F F . By this setup, we aim at simulating the axial charge production in a domain where F F happens to be nonzero in a single collision event.
The paper is organized as follows. In Sec. II, the formulation of the Dirac field and the axial anomaly in the boost-invariantly expanding system is explained. In Sec. III, we first review the CGC initial conditions for the gauge fields and the quarks fields, and then we explicitly construct the color charge distribution that realizes the flux-tube structure of the glasma color field. After we discuss the formulation of the problem on the real-time lattice in Sec. IV, we present our numerical results in Sec. V. First, we consider uniform glasma fields by taking the limit of large flux-tube width and verify that the axial anomaly is correctly realized on the real-time lattice in the expanding geometry. Then we show results for the glasma flux tube configuration. Section VI is devoted to concluding remarks.

II. AXIAL ANOMALY IN THE BJORKEN FRAME
In the high energy limit of a heavy-ion collision, the system right after the collision shows boostinvariant expansion to the longitudinal direction, which can be conveniently described in terms of proper time τ and space-time rapidity η defined by as well as transverse coordinates x ⊥ = (x, y).
In the original rest frame, the vacuum expectation of the axial current density is given by where Ψ(x) denotes the quark field operator. The axial current obeys the Adler-Bell-Jackiw [58,59] anomaly equation where m denotes quark mass and the summation over the color indices a = 1, · · · , N 2 c − 1 is implied. 2 To respect the boost invariance, we compute all expectation values in the Bjorken frame that moves to the longitudinal direction with the local velocity of v z = z/t = tanh η. Moving to the Bjorken frame, the axial current is transformed as where 2 In this paper, we consider only one quark flavor. Since we neglect the backreaction, different flavors contribute to the axial anomaly just additively.
is the boost operator to the Bjorken frame for four-vectors. Here and in the following, quantities in the Bjorken frame are denoted with a hat . To solve the Dirac equation under boost-invariant background fields, it is convenient to treat the quark field operator boosted to the Bjorken frame, where e − η 2 γ 0 γ 3 is the boost operator to the Bjorken frame for spinors. The factor √ τ is just a convention to make the following equation simpler. The Dirac equation for the boosted field is with D µ = ∂ µ +igA µ being the covariant derivative [44]. Here and in the following, repeated indices i imply the summation over i = 1, 2. In terms of the boosted field operator, the axial current in the Bjorken frame (4) is simply rewritten as Such expectations of fermion operators can be expressed by fermion mode functions [54]. The mode functions are introduced by the mode expansion of the field operator where a p ⊥ ,ν,s,c and b p ⊥ ,ν,s,c are annihilation operators of a quark and an antiquark, respectively, having momentum (p ⊥ , ν), which is conjugate to (x ⊥ , η), spin s and color c. The superscripts + and − in the mode functions distinguish the positive and the negative energy solutions. By substituting Eq. (9) into (8), we find the expression In terms of the quantities in the Bjorken frame, the anomaly equation (3) is rewritten as We note that the two terms in the right hand side are Lorentz scalars. In boost-invariant background fields we consider in this study, the η-derivative term drops. The axial charge density per unit transverse area and unit space-time rapidity is related with j 0 5 as where we have introduced a shorthand notation for the pseudoscalar condensate, In deriving Eq. (13), we have assumed that the axial charge density is vanishing at τ = 0, which is the case for the CGC initial condition discussed in the next section.

III. CGC INITIAL CONDITIONS
In the CGC effective theory, hard degrees of freedom in a high energy nucleus are treated as classical sources of radiation, while soft degrees of freedom are described as classical gauge fields that couple to the hard sources via the Yang-Mills equations The classical sources of two colliding nuclei running with the speed of light are represented by a current where ρ (n) (x ⊥ ) (n = 1, 2) denote the color charge densities of the two nuclei in the transverse plane, and x ± are light-cone coordinates defined by x ± = (t ± z)/ √ 2. With the initial condition A µ = 0 at t → −∞, the Yang-Mills equations can be solved analytically up to the τ = 0 + surface 3 [60].
The solution at τ = 0 + in the Fock-Schwinger gauge A τ = 0 is where α i (n) are transverse pure gauges associated with gauge factors tudinal ones, This color field configuration in general carries nonzero topological charge density E a · B a , and hence can generate axial charge through quark production [2]. We emphasize that nonzero E a ·B a exists only after a collision. The CGC color fields localized on the light cones, x ± = 0, are only transverse ones and electric and magnetic fields are orthogonal to each other, E ⊥ B [2]. Therefore, the axial charge density is vanishing at the instant of a collision, τ = 0.
In the McLerran-Venugopalan (MV) model [61], the color charges ρ (n) (x ⊥ ) are assumed to be distributed randomly in the transverse plane according to a Gaussian probability distribution. In the present study, we consider a fixed configuration of ρ (n) (x ⊥ ) that corresponds to a flux tubelike configuration of the color electromagnetic fields in order to elucidate the nonequilibrium axial charge production in a simpler situation.
The leading order dynamics of fermions under strong gauge fields can be described by the Dirac equation for the fermion mode functions [51,54]. Under the CGC gauge fields, the Dirac equation can be solved analytically up to the τ = 0 + surface [43,44]. The mode solution at τ = 0 + with the initial condition of the negative-energy free spinor at t → −∞ is where M p = m 2 + p 2 ⊥ is the transverse mass, v s (p) is the negative-energy free spinor, and χ a (a = 1, · · · , N c ) are unit vectors in the color space [44]. 4 The Fourier transform of the gauge factors V n (p ⊥ ) are defined as 4 This expression is slightly changed from that given in Ref. [44]; the sign of the transverse momentum index is flipped, p ⊥ ↔ −p ⊥ , the integration variables are shifted as q ⊥ → q ⊥ + p ⊥ , and the overall normalizations differ by the factor √ 4π.
For the spinor satisfying v † s (p)v s (p ) = 2 p 2 + m 2 δ ss , the mode functions are normalized as A. Glasma flux tube We will construct the gauge factors V (n) such that the initial longitudinal electric and magnetic fields have localized x ⊥ dependences. In the following, we consider the color SU(2) theory for simplicity.
We suppose that each of the color sources ρ (n) (x ⊥ ) has only one color component, and write with real functions Θ n (x ⊥ ) and the Pauli matrices σ i . For these V (n) , the transverse gauge fields Since the color orientations of α (1) and α (2) are different, the electric and magnetic fields right after the collision given by (21) and (22) can be nonzero for these configurations. In a realistic situation of a heavy-ion collision, the color charges ρ (n) (x ⊥ ) are distributed randomly in the transverse plane. In that case, the color orientations of the pure gauges α (n) at a certain x ⊥ should be given randomly, and the value of E a ·B a right after the collision fluctuates in the transverse plane. In Eqs. (26), we have chosen one of specific color configurations that realize nonzero E a · B a right after the collision. We note that the color configurations of each nucleus, Eqs. (26), should be understood as ones given in a common gauge that is globally fixed. Until the instant of collision, the color orientations of each nucleus do not have any physical meaning because the two nuclei are causally separated and one can apply independent gauge transformations to them. However, it is not the case anymore after the collision. Since the color fields after the collision are generated by the interaction between the two nuclei, they should be computed in a common gauge.
We further assume that the functions Θ n (x ⊥ ) depend on x and y only through the combination of ξ n = x cos θ n + y sin θ n with real parameters θ n . Then, α i (n) can be written as where we have introduced For these α i (n) , the initial electric and magnetic fields read The topological charge density E a ·B a is nonzero when θ 1 − θ 2 = πn/2 (n: integers).
To gain a flux tube-like structure with a Gaussian profile for the electric and magnetic fields, we assume where Q n are parameters that have the mass dimension one, ∆ characterizes the width of a flux tube, and Erf(x) is the error function. This leads Then, the x ⊥ dependence of the electric and magnetic fields turns out to be a distorted Gaussian 5 , B. Uniform glasma By taking the limit of an infinitely wide flux tube, ∆ → ∞, we obtain a uniform color field configuration that carry nonzero topological charge. In this limit, the functions Θ n (ξ n ) become 5 Although the electric and magnetic fields are localized in the transverse plane (unless cos(θ1 − θ2) = ±1), the color charges ρ (n) that generate these fields have infinitely elongated structures to the directions along (x, y) = (sin θn, − cos θn): The electric and magnetic fields are induced in the overlapped region of the color charges ρ (1) and ρ (2) . Note that these color charge distributions are globally color neutral, d 2 x ⊥ ρ (n) (x ⊥ ) = 0. This condition is necessary for the inverse Laplacian in Eq. (20) being well-defined. and the electric and magnetic fields become uniform, By substituting V (n) = exp iQ n ξ n σn 2 into Eq. (23), we find the quark mode function at τ = 0 + for this background, where we have introduced two-dimensional vectors q n = ( 1 2 Q n cos θ n , 1 2 Q n sin θ n ). By using this expression, one can directly confirm that the axial charge density j 0 5 is vanishing at τ = 0 + . Nevertheless this expression is indicative of axial charge imbalance that is induced right after τ = 0 + , as the γ + and γ − projections are asymmetric when q 1 = q 2 .

IV. LATTICE FORMULATION
Since it is difficult to analytically solve the Dirac equation and the classical Yang-Mills equation for τ > 0 with the CGC initial conditions, we resort to numerical computations on the real-time lattice. The lattice discretization method we employ is the same as that used in Ref. [50]. We will review it in this section to make the paper self-contained, and also explain issues specific to the present study. The space coordinates (x ⊥ , η) are discretized into N ⊥ × N ⊥ × N η grids with spacings (a ⊥ , a ⊥ , a η ). The transverse and the longitudinal system size are L ⊥ = N ⊥ a ⊥ and L η = N η a η , respectively. The periodic boundary condition is imposed on all the fields.
On the spatial lattice, the gauge fields are represented by link variables U i , U η and electric fields E i , E η , where i = 1, 2 denotes the transverse directions. The link variables are related with the original gauge fields as The physical electric fields in the Bjorken frame are related with the lattice electric fields as The lattice version of the classical Yang-Mills equations in the expanding geometry are and where the subscript 'traceless' means The plaquettes variables U µ,ν (x) and U µ,−ν (x) are defined by and withμ representing the unit displacement in the µ direction on the lattice. To give a definition of magnetic fields, we further introduce other kinds of plaquettes, where T a is the generators of SU(N c ).

B. Quark sector
The Dirac equation for the mode functions on the lattice has almost the same form as that for the field operator in the continuum (7), The differences are the form of the covariant derivative and the addition of the Wilson term. To improve the convergence to the continuum limit, we employ the O(a 3 )-improved lattice derivatives [48,49]. The covariant derivative is given by with coefficients c 1 = 4/3 and c 2 = −1/6. As the spatial Wilson term extended to the expanding geometry, we employ where r ⊥ and r η are real parameters, and T is a quantity that has the dimension of time. In Ref. [50], T = τ 0 (initial time) was employed. In the present study, we use T = τ , which is essential to compute the early-time behavior of the axial charge production.
On the lattice with the periodic boundary condition, the plane wave factor is replaced as where n x,y,η are integers for the space coordinates; (x, y, η) = (a ⊥ n x , a ⊥ n y , a η n η ), while integers k x,y,ν specify the momentum modes. By the latter integers, fermion momenta are discretized as and The dispersion of these lattice momenta is illustrated in Fig. 1. The regions for integers k satisfying |k| > k max ≈ 0.286N correspond to fermion doubler modes. In the Wilson fermion method, the doublers are decoupled from the physical modes being made heavy by the Wilson term.
For the initial conditions of the quark mode functions, we replace (23) by where momenta q ⊥ = (q x , q y ) are associated with integers (j x , j y ). The transverse masses M p and M p+q in the above express contain the contribution from the Wilson term, where the Wilson mass m W depends on k x,y,ν and τ as The Fourier transform of the gauge factors is discretized as The most costly part in our numerical computation is solving the Dirac equation for the fermion mode functions. For general background gauge fields, the numerical cost to solve Eq. (52) is proportional to the square of the lattice size, (N 2 ⊥ N η ) 2 , since we have to solve the equations for all the modes and each mode function has dependence on the space coordinates. In the longitudinally expanding system, the longitudinal lattice size N η especially needs to be large to resolve the longitudinal momentum scales that rapidly vary as 1/τ , and the numerical cost becomes unacceptably expensive. One possible way to reduce the numerical cost is the use of the stochastic method for fermions [51,62]. However, this method is not suitable to the computations of local quantities which are not averaged over space, especially quantities related with axial anomaly, because of large statistical errors. In this study, therefore, we stick to the direct method of solving the Dirac equation for the mode functions, which does not involve statistical errors. Fortunately, the cost of the mode function method can be reduced by the factor of N η in boost-invariant backgrounds because the η-dependence of the mode functions is known to be e iνη .
To further reduce the numerical cost, we introduce cutoffs in the momentum space. As illustrated in Fig. 1, the lattice fermion modes contain unphysical doubler modes. A naive way to regulate the doubler modes in the mode function method is just to cut off these modes from the computation. This approach has been successfully employed in Ref. [45] and also in Ref. [50] being combined with the stochastic fermion method. However, this approach is not applicable to the computation of the axial anomaly because it amounts to introducing a cutoff for canonical momentum and thus breaks the gauge invariance [47]. Therefore, in the present study we employ the Wilson fermion, which amounts to introducing a cutoff for kinetic momentum. Once the doubler modes are suppressed by the Wilson term in a gauge-invariant way, we can introduce momentum cutoffs without affecting the axial anomaly. As depicted in Fig. 1, we put a cutoff k Λ for the momentum integers between k max and N/2, and excluded the modes for |k| > k Λ from the computation. We have explicitly confirmed for the uniform glasma configuration that this cutoff does not alter the results for the axial charge production as long as k Λ is not too close to k max . We typically choose a value of k Λ so that about 20% of modes is excluded in each dimension. By this, the total numerical cost becomes about half (0.8 3 ≈ 0.5).

C. Axial anomaly on the lattice
On the lattice, the axial anomaly is a nontrivial issue. If one uses a naively discretized fermion action, degenerated doubler modes appear as shown in Fig. 1 and the axial anomaly is not realized due to the cancellation among the doublers [53]. For the axial anomaly, one needs to eliminate the doublers to spoil this cancellation. As already discussed in the previous subsection, we employ the Wilson fermion method for this purpose. The Wilson term (54) introduced in the Dirac equation (52) makes the doublers as heavy as the lattice ultraviolet (UV) cutoff scale and decouples them from the dynamics. In the following, we explain how the axial anomaly is realized by the Wilson fermion.
Since we treat the time τ as a continuum variable, the definition of the time component of the axial current is the same as that in the continuum, The spatial components must be modified on the lattice as [49] j i This expression is valid also for i = 3 though the third component does not appear in the following equations due to the boost invariance. For this definition of the axial current, the anomaly equation with ∇ i denoting the backward difference The pseudoscalar condensateη(x) is represented by the mode functions as From the Dirac equation (52) that includes the Wilson term, one can derive the relation where w(x) stands for the expectation of a fermion operator involving iγ 5 and W , Comparing this equation with Eq. (64), we notice that the axial anomaly is realized by the Wilson fermion if This relation has been proven to hold in the continuum limit in the context of the Euclidean lattice gauge theory [63,64]. In the context of the real-time lattice computations, it has been numerically confirmed for nonexpanding systems [47][48][49]57].

D. Boundary condition
As already noted, we impose the periodic boundary condition (b.c.) on the spatial lattice. For the flux tube configuration introduced in Sec. III A, we need a special care for the periodicity in the transverse directions. In the gauge sector, the transverse link variables U i and the longitudinal electric field E η must satisfy the periodic b.c. at the initial time. Other components trivially satisfy the b.c. as they are vanishing then. Since the initial longitudinal electric field has a localized Gaussian profile, the field is vanishing at the boundaries and satisfies the periodic b.c. if the flux tube is located sufficiently away from the boundaries. In contrast, the initial transverse gauge fields have an elongated structure to the directions along (x, y) = (sin θ 1 , − cos θ 1 ) and (sin θ 2 , − cos θ 2 ) for the functions Θ n (x ⊥ ) given by (32). On a square transverse lattice and for θ 1 , θ 2 = πn/2 (n : integers), it is impossible that these gauge fields satisfy the periodic b.c. as long as there is only one flux tube. Therefore, we put two flux tubes in the transverse plane. In the following, we fix the angle parameters to be θ 1 = 0 and θ 2 = π/4, with which the initial electric and magnetic fields have the same field strength. For these angles, gauge configuration with two flux tubes located at (x 1 , y 1 ) = (L ⊥ /4, L ⊥ /4) and (x 2 , y 2 ) = (3L ⊥ /4, 3L ⊥ /4) satisfies the periodic b.c. as long as ∆ L ⊥ . This configuration is realized by the functions Also the quark mode functions (58) must satisfy the periodic b.c. If the factors V n (x ⊥ ) = exp [iΘ n (x ⊥ )σ n /2] are periodic, this requirement is fulfilled. For the angles θ 1 = 0 and θ 2 = π/4, the differences of Θ n (x, y) at the boundaries are for ∆ L ⊥ . Therefore, in order that the factors are periodic, the flux is quantized as √ πQ 1 ∆ = 2πn 1 , 1 2 √ πQ 2 ∆ = 2πn 2 , (n 1 , n 2 : integers). (77)

V. NUMERICAL RESULTS
In this section, we present numerical results of solving the Yang-Mills equations (41)(42)(43)(44) and the Dirac equation (52) on the lattice for the axial charge production in the longitudinally expanding geometry.

A. Uniform glasma
As a simple test for the real-time lattice computations of the axial charge production in the expanding geometry, we first consider the uniform glasma configuration introduced in Sec. III B.
For background fields which are uniform not only in the η-direction but also in the transverse directions, the numerical cost to solve the Dirac equation is significantly reduced since the spacedependence of the mode functions is completely known.
We denote the typical energy scale of the glasma by Q, and initialize the gauge fields by setting Q 1 = Q 2 = 2 1/4 Q, θ 1 = 0, and θ 2 = π/4. This choice of the parameters results in the initial gauge fields gA 1 (τ = 0) = 2 1/4 Q σ 1 2 and Since all these fields are uniform, they trivially satisfy the periodic boundary condition.
As the initial condition for the quark mode functions, we employ the expression (39) after replacing the plane wave factor and the momenta by corresponding lattice expressions. By the replacement (55), also the quark initial condition (39) satisfies the periodic boundary condition.
In actual numerical computations, we cannot take the initial time τ 0 to be exactly zero. Instead we take a small value of τ 0 as Qτ 0 = 10 −3 . We have confirmed that varying it between 5 · 10 −4 and In this and the following figures, all quantities are shown in dimensionless unit scaled by appropriate powers of Q. Furthermore, the factor g 2 is multiplied to the field strength to make it order one. 6 The result shown in Fig. 2 looks similar to that with the MV model initial condition first shown in Ref. [2]; Initially only the longitudinal components are nonzero. As the longitudinal components decrease in time, the transverse components are induced, and eventually all the components decay in time. We point out that the decay of the fields is in fact nontrivial for the uniform system. For example, if the initial field has only longitudinal electric component (which can be realized by e.g. θ 1 − θ 2 = 0), the field strength stays constant even in the longitudinally expanding system. 7 This is because the nonlinear terms in the Yang-Mills equations do not play any role for that field configuration. Therefore, the decay of the uniform glasma seen in Fig. 2 is caused by the interplay between the system expansion and the nonlinear interaction of the color fields.
As shown in Fig. 3, E a ·B a exhibits damped oscillation in time. Due to the oscillation, E a ·B a changes its sign leading to nonmonotonic behavior of the axial charge density as we will discuss below. Qτ In a uniform system, the transverse divergence term of the axial current disappears, and the anomaly relation leads Each term in this equation is plotted as a function of time in Fig. 4 for quark mass m/Q = 0.01. For such light quark mass, the pseudoscalar condensate term 2m τ 0 dτ η is negligible, and hence the remaining two terms must agree for the realization of the axial anomaly on the lattice. Indeed, the axial charge density and the time-integration of the topological charge density show an agreement, though we see some deviations especially at later time. This result demonstrates that the axial anomaly can be described by the Wilson fermion even in the longitudinally expanding geometry.
In Figs. 5 and 6, numerical results for the axial charge density computed with different lattice parameters are shown. Both for the longitudinal (Fig. 5) and the transverse (Fig. 6) lattice parameters, the results are nearly insensitive to the changes of either the UV cutoff scale 1/a and the infrared scale 1/L. So far, we have shown numerical results for light quark mass m/Q = 0.01, in which case the pseudoscalar condensate term is negligible. We now present in Fig. 7 the dependence of the axial charge density on the quark mass. In the computations with the quark masses m/Q = 0.3 and 0.5, we have used the technique of the Wilson parameter averaging [49]. For lighter masses, m/Q = 0.01 and 0.1, the curves are nearly overlapped indicating these quarks can be regarded as almost massless. In contrast, the results for heavier quarks, m/Q = 0.3 and 0.5, show significant deviations from those for the light quarks. This is because the pseudoscalar condensate term is comparable to other terms in the anomaly relation (82) for these masses. To illustrate it, we plot τ 0 τ E a · B a dτ indicates the realization of the axial anomaly. For this quark mass, the pseudoscalar condensate term is indeed as large as other terms especially at later times. However, at very early times, Qτ < ∼ 0.5, the rise of the pseudoscalar condensate term is slower than the other terms. This is the reason why the axial charge densities show little dependence on the quark masses at the early times in Fig. 7. At the later times, both of the axial charge density and the pseudoscalar condensate term show oscillation.
Interestingly, their oscillation phases are different, and therefore the pseudoscalar condensate does not always diminish the axial charge density in this oscillating background field.
Before closing this subsection, we show a rough estimate of the axial charge density in physical Qτ these values should not be taken too seriously because the color field configuration considered in the calculation is not so realistic; it is uniform in the transverse plane and the color directions of the fields are chosen such that E a ·B a is maximum for fixed energy density. As we see in the next subsection, the amount of axial charge density is reduced for inhomogeneous configurations because the spatial divergence term of the axial current takes some fraction in the anomaly relation. Also considering the color SU(3) theory instead of SU(2) employed in study would modify the results quantitatively.

B. Glasma flux tubes
In the previous subsection, we have confirmed that the axial anomaly on the real-time lattice in the expanding geometry can be described by using the Wilson fermion for the uniform glasma Next, we present numerical results of solving the Dirac equation under the glasma flux tube configuration for the axial charge production. First, let us look at the anomaly relation averaged over space, where the pseudoscalar condensate term is dropped out since it is negligible for m/Q = 0.01, and the transverse divergence term of the axial current is absent as it disappears by the space averaging. This relation is examined in Fig. 12, where the two terms in both sides of the equation are plotted separately as a function of time. Although there are noticeable deviations, the overall behaviors of the two curves roughly agree, demonstrating the realization of the axial anomaly in this transversally inhomogeneous system. Compared to the corresponding result in the uniform system shown in Fig. 4, what is remarkable is the nearly monotonic increase at later times Qτ > ∼ 3.
In the uniform glasma, E a · B a continues to oscillate even at the later times, and hence also its time integral does. In the flux tube configuration, E a ·B a decays faster, and this decaying behavior rather helps nonzero axial charge density remain at the later times. For example, if E a ·B a decays as 1/τ , the time-integral τ 0 τ E a ·B a dτ increases linearly in time. A similar observation, axial charges persist to be present after coherent gauge fields die out due to an instability, has been is trivially satisfied as we solve the lattice Dirac equation (52). Therefore, the agreement between the terms τ 0 τ E a ·B a dτ and τ 0 τ w dτ is the condition for the realization of the anomaly relation (13). Both at Qτ = 0.5 and 1, the overall behavior of the curves for these two terms roughly agree although at the most 20-40% of deviations are present. We infer that these local deviations are due to insufficient resolutions at UV scales. In an inhomogeneous system, resolutions at small scales in the coordinates space, or equivalently UV scales in the momentum space, is more important than in a uniform system. Since the Wilson term changes the UV sector of the theory, space-dependent quantities may be more affected by it. In our system, the situation is further complicated by the fact that the longitudinal momentum scales vary rapidly in time as 1/τ . We expect that the use of finer lattices would improve the accuracy. However, it is extraordinary challenging due to expensive numerical cost for the mode function method. Improvement of the computational method is desirable. We leave these issues for future investigations.
Having compromised with this accuracy of the anomaly relation in the present study, let us compare the axial charge density dN 5 d 2 x ⊥ dη and the the spatial divergence term of the axial current Qx τ 0 τ ∂ i j i 5 dτ . The latter term characterizes the outflow of axial charge. At the very early time, Qτ = 0.5, the outflow term τ 0 τ ∂ i j i 5 dτ is much smaller than other terms. In this case, one can compute the axial charge density directly from E a ·B a without solving the Dirac equation similarly to the uniform system. Already at Qτ = 1, however, the outflow term becomes comparable to other terms signaling the propagation of the axial charge in the transverse plane. When the outflow term is not negligible, one cannot predict the amount of axial charge anymore directly from the anomaly relation. It is necessary to compute the full quantum dynamics of the quark fields by solving the Dirac equation as we have done in this study.

VI. CONCLUSIONS
In this paper we have investigated the axial charge production in the early stage of heavy-ion collisions by using the real-time lattice simulation method for classical gauge fields and quantum quark fields. To consistently include the effects of the colliding nuclei on the evolution of the quark fields, the solution of the Dirac equation under the CGC gauge fields has been used for the initial condition of the numerical computation for the time evolution after the collision.
First, we considered the uniform configuration for the glasma gauge fields, and demonstrated that the Wilson fermion method generalized to the expanding geometry can correctly describe the axial anomaly on the lattice. In case of the uniform configuration, the color electromagnetic fields continue to coherently oscillate. Consequently, the axial charge density per unit rapidity exhibits oscillating behavior. Next, we computed the evolution of the glasma flux tubes. Due to the longitudinal expansion and the dynamics in the transverse plane, the color fields show decoherence at later times leading to the decay of the topological charge. The oscillating behavior of the axial charge density seen at earlier times is terminated by the decoherence. We find that the decay of the fields rather helps nonzero axial charge persists to be present at later times.
The present work may provide an important basis for future investigations, which include firstprinciples-based simulations of the CME in real-time along the lines of [48,49]. It is also important to consider more realistic configurations of the gauge fields in heavy-ion collisions. In that case, the net axial charge would be vanishing after space or event averaging, and one has to investigate fluctuations of axial charge and their possible connections to observables.