Lattice implementation of Abelian gauge theories with Chern-Simons number and an axion field

Real time evolution of classical gauge fields is relevant for a number of applications in particle physics and cosmology, ranging from the early Universe to dynamics of quark-gluon plasma. We present a lattice formulation of the interaction between a $shift$-symmetric field and some $U(1)$ gauge sector, $a(x)\tilde{F}_{\mu\nu}F^{\mu\nu}$, reproducing the continuum limit to order $\mathcal{O}(dx_\mu^2)$ and obeying the following properties: (i) the system is gauge invariant and (ii) shift symmetry is exact on the lattice. For this end we construct a definition of the {\it topological number density} $Q = \tilde{F}_{\mu\nu}F^{\mu\nu}$ that admits a lattice total derivative representation $Q = \Delta_\mu^+ K^\mu$, reproducing to order $\mathcal{O}(dx_\mu^2)$ the continuum expression $Q = \partial_\mu K^\mu \propto \vec E \cdot \vec B$. If we consider a homogeneous field $a(x) = a(t)$, the system can be mapped into an Abelian gauge theory with Hamiltonian containing a Chern-Simons term for the gauge fields. This allow us to study in an accompanying paper the real time dynamics of fermion number non-conservation (or chirality breaking) in Abelian gauge theories at finite temperature. When $a(x) = a(\vec x,t)$ is inhomogeneous, the set of lattice equations of motion do not admit however a simple explicit local solution (while preserving an $\mathcal{O}(dx_\mu^2)$ accuracy). We discuss an iterative scheme allowing to overcome this difficulty.

The numerical procedures for Abelian or non-Abelian scalar-gauge theories, such as the Standard Model is well developed, see e.g. [12] for first studies in 3+1 dimensions. In the simplest realisation it consists in the following steps (for inclusion in numerical simulations of hard thermal loops and Langevin-like dynamics see [49]). One uses the standard lattice formulation of the gauge theory action with exact lattice gauge invariance and Minkowski signature of the metric. The variation of the action with respect to the gauge fields (living at the links of the lattice) and scalar fields (living at the lattice sites) produce a system of equations in which the dynamical variables at time slice t n are expressed via those at two preceding times t n−1 and t n−2 . So, giving an initial condition at t 0 and t 1 (the Cauchy problem) allows to follow the evolution of the system and address all sorts of questions one is interested in. The initial conditions are chosen depending on the physical system under consideration: for sphaleron transitions they are taken from an equilibrium ensemble at some temperature, for preheating from knowing the spectrum of initial fluctuations of the fields after inflation, etc.
It is very important that the procedure discussed above is gauge invariant. The variation of the action with respect to the zero component of the gauge field gives the lattice Gauss constraint, which is exactly conserved during the (lattice) time evolution. The gauge invariance insures the stability of numerics and keeps an (approximate) energy conservation; features that are necessary for a continuum interpretations of the results. It is an empirical fact that the formulations that are not gauge invariant in discrete space-time (and only invariant when the lattice spacing and time step go to zero) are plagued with different non-physical numerical instabilities.
For a number of applications the simplest gauge-scalar actions should be extended by an addition of the pieces that contain the so called topological term or P ontryagin density, Q = F µνF µν , where the dual of the field strength is defined as usual byF µν ≡ 1 2 µναβ F αβ , with µναβ being the completely anti-symmetric tensor in four dimensions, with 0123 ≡ 1. The most interesting examples contain an axion field coupled linearly to Q as a(x)F µνF µν , and the theories with non-zero chemical potential µ for chiral fermions, leading to an effective bosonic Hamiltonian containing the Chern-Simons term for the gauge fields µn cs , with n cs ∝ d 4 x Q.
To investigate the time evolution of these systems one is faced with the following problem. The realisation of Q on the lattice should be done in such a way that the continuum topological properties of Q hold: the integral of Q over the volume of spacetime can be expressed via an integral over the boundary. To put it in other words, in continuum one can write where K µ is the Chern-Simons current. The lattice analogue of this relation is where ∆ + µ is the lattice difference in positive direction of µ axis (more details are provided in Sections 3,4 and 5). If the lattice (gauge-invariant) definition of Q does not satisfy this property, we will get unwanted lattice artifacts and the continuum extrapolation would be difficult. There are several interesting quantities which are very sensitive to the topological property (1.1). For example, it is Eq. (1.1) which makes the axion mass m a to be zero in all orders of perturbation theory (m a = 0 being a non-perturbative phenomenon). The extraction of the non-Abelian sphaleron rate in the symmetric phase of the electroweak theory from diffusion of Chern-Simons number requires a careful construction of the lattice version of Q in which the property (1.1) is still approximate, but as precise as possible [15].
The aim of the present paper is to set up a lattice gauge invariant formulation for real time simulations of Abelian gauge theories, respecting the topological property Eq. (1.2) of Q exactly. As for non-Abelian theories, the mission seems to be impossible due to well known difficulties of defining Q obeying the property (1.2) on the lattice [50,51]. The applications of our formulation can include, for example, the study of the late stage of inflation and preheating in axion-inflation models [52][53][54][55][56][57][58][59], the clarification of the role of Standard Model hypercharge group U(1) in baryogenesis and in magnetic field generation [60][61][62][63], the modeling the different aspects of chiral magnetic effects [60,[64][65][66], or the study of the problem of chiral fermionic charge evolution in high temperature electrodynamics [67,68]. An accompanying paper [69] is devoted to the last subject.
This paper is organised as follows. In Sect. 2 we review the classical equations of motion in the continuum of an Abelian gauge theory with an axion field. We also discuss how to map that system into an Abelian gauge theory with chemical potential. In Sect. 3 we review the essence of the non-compact lattice formulation of an Abelian gauge theory, so that we set notation and conventions for the following sections. In Sect. 4 we build a lattice implementation of an axionic-interaction a(x)F µν F µν , deriving step by step the necessary ingredients to achieve a formulation consistent with the (lattice version) of the Bianchi identities, and solvable by an iterative scheme of evolution. We first consider in Sect. 4.1 the case of a homogeneous axion a(x) = a(t), and later generalize to a fully inhomogeneous axion a(x) = a(t, x) in Sect. 4.2. In Sect. 5 we finally discuss the lattice formulation of the Chern-Simons number n cs ∝ d 4 x Q based on the lattice version(s) of Q =F µν F µν developed in Sect. 4. We put special care in the need to achieve a lattice formulation that admits a total derivative representation for Q as in Eq. (1.2). In Sect. 6 we summarize our results and discuss some of the potential applications in particle physics and cosmology.
2 Abelian gauge theory with an axion. Theory in the continuum Let us begin by considering the action of an Abelian gauge theory in flat space-time 1 , in the presence of an axion-type field a(x) linearly coupled to the P ontryagin density F µνF µν of a U(1) gauge field, We consider a 'Higgs' sector as the gauge field, e 2 the gauge coupling strength, and D µ ≡ ∂ µ − iA µ the covariant derivative. The field strength and its dual are defined as usual by F µν ≡ ∂ µ A ν − ∂ ν A µ andF µν ≡ 1 2 µναβ F αβ . Given our choice of A µ ≡ (φ, A), we define the electric and magnetic fields as Given our metric signature, we obtain F µνF µν = +4 E B, and arrive at the final vectorial expressions of Eq. (2.1).
Lagrangian Eq. (2.1) describes a system of scalar electro-dynamics in the presence of an axion-like field a(x), with M some mass scale undetermined at this point. Note that we maintain explicitly the speed of propagation of the axion c 2 s as a free parameter, as this will be convenient for us later on. Action Eq.
Varying the action, one obtains the equations of motion (EOM) where the (unit-charge) current is defined as j µ = 2Im{ϕ * D µ ϕ}, so that Eqs. (2.2)-(2.4) can be rewritten in a vectorial form as . The promotion of a(x) into a dynamical field can therefore affect notably the dynamics of the system with respect to standard scalar electrodynamics described by Eq. (2.1) with a(x) = 0.

Mimicking a chemical potential
Interestingly, through an adequate interpretation of field variables and parameters, Eq. (2.1) can be mapped into the description of a gauge theory with a chemical potential µ for chiral fermions (for more details see [69]). Starting from Eq. (2.1) will allow us, through a formal trick, to bring up a Lagrangian formulation into this problem. In order to see this, let us begin by demanding that a(x) = a(t) is a spatially homogeneous field, so that We then introduce the following convenient 'dimensionally reduced' variables a ≡ αM ,ȧ ≡ µM , (2.11) so that Eq. (2.1) is reduced to As we will see next, the requisite to describe a gauge theory at high temperatures and in the presence of a chemical potential, will fix the mass scale M and parameter c 2 s , see Eq. (2.19). We will describe first, however, the new dynamical equations that follow from minimizing the new re-written action.
Varying Eq. (2.12) we obtain the equations of motion of the system, which we write directly in a vectorial form as Once again, let us note that the terms ∝ α(˙ B − ∇ × E) in Eq. (2.7) and ∝ α ∇ B, which vanish in the continuum, are only maintained in the above equations for later convenience when discretising the system in Sect. 4. We can now fix the mass scale M and the parameter c 2 s to appropriate values, so that the set of Eqs. (2.13)-(2.16) properly describe an Abelian gauge theory with the chemical potential µ for chiral fermionic charge. In particular, the EOM of a chemical potential follows from anomaly equation [70,71] and in our case has the form [68]  Let us note that in the original action Eq. (2.1), c 2 s represents the speed of propagation of the axion field. However, as we considering the axion now as a homogeneous field a(x) = a(t), c 2 s represents simply a parameter in the theory. One should not conclude therefore that the value given by Eq. related discussion see, e.g. [72]). We just intent to set notation and basic concepts, which we will be used later on when discussing the lattice formulation of an Abelian gauge theory with an axion. A reader already familiar with lattice gauge invariant techniques can skip this section and jump directly into Sect. 4.

Non-compact formulation of Scalar-Electrodynamics
Let us, first of all, set some notation. We will not consider summation over repeated indices, as this can lead to confusion. A lattice point n = (n o , n) = (n o , n 1 , n 2 , n 3 ) displaced in the µ−direction by one unit lattice spacing, n +μ, will be often referred simply as n + µ or by +µ. For example, ϕ +µ ≡ ϕ(n +μ). Components of gauge fields live in between lattice sites in the direction of the component, so A µ ≡ A µ (n + 1 2μ ), A µ,+ν ≡ A µ (n + 1 2μ +ν), etc. For simplicity of the notation, we will refer to both the lattice spacing ∆x and the time step ∆t, simply as dx, so if we write e.g. +μdx, this should be interpreted as a time step advancement +∆t if µ = 0, or as a unit displacement +∆x in a given spatial direction µ = 1, 2 or 3. We will loosely speak of the lattice spacing order O(dx), independently or whether we are referring to O(∆t) or O(∆x).
We define a lattice link, as usual, like ) . Forward (+) and backward (-), ordinary and covariant derivatives, are defined in the lattice by where we have indicated the order in the lattice spacing to which one recovers the continuum limit, as well as the natural space-time location in the continuum the derivatives live. A lattice gauge transformation under U (1) corresponds to ϕ(n) −→ e +iβ(n) ϕ(n) , with β an arbitrary function, so that the links and covariant derivatives transform as We will use these transformation rules to build a gauge invariant lattice action in the following. In this section we ignore the axion field, so for the time being we just consider a lattice action for scalar-electrodynamics only. Using a non-compact formulation, we can write from where the lattice gauge invariance (based on the transformations defined above) is rather explicit. We will refer to Eq. (3.4) as the Abelian − Higgs (AH) lattice action. As we will show in Section 3.2, this action reproduces to order O(dx 2 ) the continuum action Eq. (2.1) in the absence of an axion field (a(x) = 0). Varying Eq. (3.4) with respect to the different fields, we obtain the lattice equivalent of the dynamical equations, which read (taking the Coulomb Gauge 3.2 Recovering the continuum limit to order O(dx 2 ) Let us note that, given a continuum action S = dt d 3 x L C , with lagrangian given by the sum of various operators L C = p O (p) C , we try to emulate the same physical system by defining a lattice action S = ∆t∆x 3 n,no L L , with the lattice lagrangian given by the sum of lattice operators L L = p O If we were to expand ∆ + µ ϕ around the position x = n · dx, we would obtain C +O(dx), so the continuum limit is only reproduced to linear order. The natural site where O (k) L really lives is however, not x = ndx, but rather in x µ/2 ≡ (n + 1 2μ )dx, as it involves a finite difference of the field evaluated (with equal weight) at both n and n +μ lattice sites. If we expand O (k) . Therefore, in order to analyze the continuum limit of a lattice operator, one must first recognize which is the natural lattice site where it lives, and only then expand around this site. This will lead to the fact that each operator O L is interpreted as living in its natural lattice site, and reproducing the continuum limit to a common n-th order. Once this is imposed, the EOM obtained from varying a lattice action satisfying these two requisites, are guaranteed to reproduce the continuum EOM to the same order. Besides, each lattice EOM will live naturally in a well defined lattice site common to all the terms involved in the discrete equation, which determined the site around which to compute the continuum limit. Let us note, however, that the lattice site where different equations live, do not need to be the same.
Let us check all this with the terms of the action Eq. (3.4), as this will serve as a good training for our task in Sect. 4, building an appropriate lattice operator for an axionic interaction. As mentioned before, each term in action Eq. (3.4) reproduces the continuum to second order in dx, i.e. O (p) , ∀p, for as long as each term is expanded around its natural lattice site. Let us check this term by term (recall there is no implicit sum over repeated indexes), denoting by l the natural lattice site of each operator, and by x ≡ ldx the physical coordinate where the continuum limit is reproduced. A simple expansion leads to Let us turn our attention now into the lattice EOM Eqs. (3.5)-(3.7). Let us begin by the EOM of the charged scalar field, which involve the terms V ,ϕ * and D − µ D + µ ϕ. We can expand each term around l = n, as the operator D + µ ϕ lives naturally at n + 1 2μ , but D − µ (D + µ ϕ) makes it live back to l = (n + 1 2μ ) − 1 2μ = n, Consistently, the term V ,ϕ * in Eq. (3.7) [which reproduces the continuum to any order in dx], lives naturally at the same lattice site l = n. Therefore, when interpreting that the lattice operators involved in the discrete EOM of the charged scalar live at l = n, Eq. (3.7) reproduces correctly the corresponding continuum Eq. (2.13) up to O(dx 2 ) corrections. The EOM of the gauge fields contain more terms, but they all live consistently at the same lattice site, e.g. at l = n +ĵ 2 for the dynamical equation Eq. (3.6). This can be easily shown by expanding each term of this equation A similar analysis can be done expanding the terms of Eq.

Axionic-coupling
Let us now turn our discussion into the formulation of a proper lattice equivalent for the continuum interaction between the gauge fields and an axion, In this section we aim to a general formulation of an Abelian gauge theory with an axion field. For simplicity we will first start dealing with the case of a homogeneous axion α(x) = α(t) in Sect. 4.1. Our findings in Sect. 4.1 will be actually applicable as well to the case of a fully inhomogeneous axion α(x) = α(t, x), but as the latter introduces further complications, we postpone the discussion about a fully spatially-dependent axion for Sect. 4.2.
We will introduce the dimensionally reduced variables α ≡ a M , µ ≡α defined in Eq. (2.11), so we prevent this way having to drag the scale M along our derivations. Given our choice of metric signature (−, +, +, +) and gauge field representation A µ ≡ (φ, A), we find F µνF µν = +4 E B, so that we can write the above continuum action as We will refer to the interaction described by Eq. (4.1) as an axionic − coupling. Our main aim now is to formulate a latttice version of the continuum action Eq. i) The terms α(−˙ B + ∇ × E) and α ∇ B in Eqs. (2.14), (2.15), vanish in the continuum thanks to the Bianchi identities ∂ νF µν = 0, which are equivalent to˙ B = ∇ × E and ∇ B = 0. The equivalent terms in the discrete EOM are however not granted to vanish, as this depends on the lattice representation of the electric and magnetic fields, and on the choice of lattice derivatives. It is therefore crucial that we find a lattice representation of S ac so that the equivalent discrete terms in the lattice EOM vanish identically (or at least to the same order in the lattice spacing to which the discrete EOM reproduce the continuum). In other words we seek a lattice formulation of S ac so that the lattice expression of Q =F µνF µν is topological admitting a total derivative ii) Assuming that a correct version of the Bianchi identities follows naturally from a given topological lattice formulation of F µνF µν , another problem may arise. The terms ∝ µ B and ∝ E B in Eqs. (2.7), (2.9), indicate possible obstructions to achieving an explicit scheme to solve iteratively the set of lattice coupled equations reproducing the continuum Eqs. (2.6)-(2.9). Even though an implicit scheme for finite difference coupled equations can be solved by non-linear numerical methods (applied at every lattice site), this makes the continuum limit less transparent, and results typically in a computationally more expensive procedure (if not unfeasible). Therefore, achieving a simple explicit scheme for solving iteratively the set of lattice coupled equations that will mimic the continuum EOM, will be a strong requisite, unless we prove that such a scheme cannot be developed.
iii) When considering a fully inhomogeneous axion-like field a(x), terms proportional to spatial variations ∝ ∇a(x) appear in the EOM. In particular, the term ∇a × E in the rhs of electric field evolution equation Eq. (2.7)Ė i = [...] + e 2 4π 2 (∇a × E) i , introduces a 'mixing' of the E i component, naturally living in the i-th direction, with the components E j , E k , naturally living along the transverse directions to the ith axis. As electric field components live naturally in between lattice sites (at the links), this will imply a mixture of orthogonal links. Some spatial averaging over neighboring position along the i-th axis will be needed, to force the term ∇a × E (in the rhs of the equation) to live at the same location where E i lives (in the lhs of the equation). This will create a non-local interaction, possibly preventing the development of an explicit iterative scheme to solve the resulting set of finite difference coupled equations.
In the following we will investigate various lattice versions of the interaction α E B, determining their 'appropriateness' based on the ability of each lattice formulation to address the previous criteria i) − iii).

Abelian gauge theory with a homogeneous axion
We will start considering in this section the simplest case of a gauge theory with a homogeneous axion α(x) = α(t). As discussed in Sect. 2, this problem can be mapped into the description of a gauge theory in the presence of a chemical potential µ =α. Both in Sect. 4.1.1 and 4.1.2 we will stick to a(x) = a(t), simply to make more transparent the discussion about the importance of achieving a good lattice representation of the Bianchi identities, as well of an explicit iterative scheme to solve the set of coupled lattice EOM. The conclusions that will be reached in Sects. 4.1.1, 4.1.2 will be equally applicable to the case of a fully inhomogeneous axion a(x) = a(t, x), which we will address in Sect. 4.2, building up from our previous findings on the homogeneous case.

Lattice formulation of the Bianchi identities
Let us first of all take the simplest possible approach, and attempt to describe S ac using the lattice definition of electric and magnetic fields introduced in Section 3, A priori, this looks like benign lattice operator, since it describes the continuum action up to order O(dx 2 ), as we showed already in Sect. 3: each term of individually reproduces the continuum expression up to to order O(dx 2 ), when interpreting that F µν lives at n + 1 2μ + 1 2ν . When varying this lattice action with respect A i , we obtain (recall that we are assuming now α( However, contrary to the continuum analogues, they do not vanish. The reason is simple, the correct discrete version of ∇ B is rather i ∆ + i B i = 0, simply because the lattice magnetic field was defined in terms of forward derivatives, To obtain the desired result, one needs to take the divergence over B i with a forward derivative: as ∆ + i and ∆ + j commute, does not vanish. The appropriate version in the lattice of the Bianchi identity should rather be built as As anticipated, generating the appropriated vanishing terms in the discrete EOM (due to the lattice version of the Bianchi identities), is not automatically granted. The problem arises because our choice of Eq. (4.2) as a lattice operator is actually not consistent. Note that even though it reproduces the continuum result to order O(dx 2 ), it consists however in the product of three fields, α, E i and B i that live, not only in different lattice sites, also at different time steps. Whereas α lives at (n o , n), E i lives at (n o + 1 2 , n + 1 2î ) and B i lives at (n o , n + 1 2ĵ + 1 2k ). To make consistent an action formed by the sum of several lattice operators, let us recall the rule we already discussed in Sec. 3, that all operators must reproduce the continuum limit to the same order, when expended around their natural site. It is therefore implicit in that statement, that each operator must have a well defined natural site where to live. This is precisely the reason why the previous operator Eq. (4.2) is inconsistent, as there is no natural site ascribed to it. The solution passes trough "symmetrizing" the factors in the operator, so that the factors built up from different fields, live nonetheless at the same lattice site. Let us define and note that each of these fields can be expressed as For convenience we also define which reproduce the continuum as A 'symmetrized' operator that reproduces the continuum expression of S ac at l = ( n, n o + 1 2 ) to order O(dx 2 ), can be easily proposed based on the above expressions, Varying Eq. (4.13) with respect A i , produces a term in the discrete EOM of the gauge field as α j,k ijk ( , which resembles the continuum term α( ijk ∂ j E k −Ḃ i ) in Eq. (2.14). Whereas the latter vanishes thanks to the Bianchi identity in the continuum ∇ × E =˙ B, it can be shown, with a bit of algebra, that a lattice version of this identity holds as j,k ijk ( i . This implies that the new term encountered in the gauge field discrete EOM just vanishes. Similarly, when varying Eq. (4.13) with respect A o , we produce a new term i , which again resembles the term in the continuum α∂ i B i in Eq. (2.15), which vanishes due to the Bianchi identity ∇ B = 0. With a bit of algebra, it can be shown that the analogous lattice identity reads ∆ − i B = 0, so that the new term encountered in the analogous discrete Gauss law, just vanishes. This completes the proof that the new operator Eq. (4.13) represents a good lattice candidate from which to derive [together with action Eq. (3.4)] a set of coupled finite different equations reproducing correctly the functional form of continuum EOM. As we anticipated, there is however another problem yet to be circumvented, related to the solubility of a set of coupled finite different equations.

Explicit scheme for real time evolution
The operator Eq. (4.13) proposed to represent an axionic coupling, exhibits various features: i) it reproduces correctly the continuum term to order O(dx 2 ), and ii) it reproduces correctly a lattice version of the Bianchi identities, so that the discrete EOM reproduce correctly the functional form of the dynamical Eqs. (2.13)-(2.16) in the continuum. In fact, varying Eq. (4.13) with respect to the gauge fields, generates a term ∝ 1 2 µ −0 B i + µB (8) i,+0 , which reduces correctly to the continuum term µ B in Eq. (2.14), to order O(dx 2 ). At the same time, varying Eq. (4.13) with respect to α(t), generates a term ∝ 1 2 (E i + E i,−0 )B (4) i sourcing the chemical potential, which again reduces correctly to the continuum source term ∝ E B in Eq. (2.16), to order O(dx 2 ). The set of coupled discrete equations one obtains, cannot be put however in an explicit iterative scheme, because in order to find E i we need µ and A i,+0 (to obtain B (8) i,+0 ), and at the same time to find A i,+0 and µ we need E i . One would need to express the term µB (8) i,+0 in the gauge field EOM in terms of E i , and then solve for E i , but this would complicate the equations unnecessary 3 . A simpler solution consists in choosing another lattice operator that upon variation over the gauge fields, prevents the duplication of the ∼ µB i terms at the two different times. The duplication of these terms originated from Eq. (4.13) due to the fact that the electric field E i,−0 ). Hence, when varying with respect 3 Actually, it is not even feasible, in principle, to solve this system, as e.g. B 1,+0 depends A 2,+1+0 and A 3,+1+0 , which depend respectively on E2,+1 and E3,+1. So at the end, the equation for the1-component of the electric field depends on E1 but also on E2,+1 and E3,+1, and analogously for the other electric field component equations. Therefore, one cannot just simply solve for the Ei components at a given time step and lattice site, as we would need to know the electric field at all other lattice sites. We will actually also encounter this problem in Sec. 4.2 when solving the system for an inhomogenous axion-like field, so we postpone any further discussion on this till then.
to A i , two terms of the form ∼ µB i are generated, µ −0 B (8) i + µB (8) i,+0 , each evaluated at a different time.
It is easy to build some lattice operator that verifies the last requisite, i.e. not involving the sum of two gauge field conjugate momenta (i.e. electric fields) at different times, like for example This operator fails however in 'symmetrizing' the expression for E · B around a common time: as B i lives at integer times n o whereas E (2) i lives at time semi-integer timEs (n o + 1 2 ), the lattice equivalent to the terms that should be vanishing due to the discretized Bianchi identities, will fail to vanish, as it happened already with the operator S L(1) ac

Eq. (4.2).
A better solution is found if the lattice operator ∼ α(E · B) is built such that each element α, E and B, live separately in a given common space-time site. In other words, contrary to the operator S

L(3) ac
Eq. (4.14), where the electric and magnetic fields lived at different time steps separated away by half step dt 2 , now both E and B must live at the same time step. As the natural time where electric fields live is (n o + 1 2 )dt, and given that we do not want to sum over electric fields at different times, the natural option will be to make the magnetic field to live in the same semi-integer time. There are 3 options for this, where ≡ should be understood as an equivalence modulo an additive constant shift. Varying each of these action terms with respect to the gauge fields, we obtain correct vanishing versions of ∇ B and (˙ B− ∇× E) in the discrete equations of motion due to the lattice Bianchi identities, and generate the following terms [recall that here we still consider α(x) = α(t)], i,−0 (B (4) i + B Eq. (4.16) as equivalent descriptions of a final suitable action. Let us notice as well that, as a consequence of this, the kinetic term of the axion should be also defined differently compared to ordinary scalar fields, since we want to obtain a finite difference evaluated at integer times.
Putting all together, the discretized action from where the dynamical effects of the presence of a homogeneous axion to be derived reads where it is important to note that both α and A o live at space-time sites (n o + 1 2 , n), whereas A i live at (n o , n +î 2 ). Consequently, µ ≡ ∆ − o α lives naturally at (n o , n), E i at (n o + 1 2 , n + 1 2î ), E (2) i at (n o + 1 2 , n), whereas B i lives at (n o , n + 1 2ĵ + 1 2k ), and B (4) i at (n o , n). 4 In this occasion this occurs because even though there is a single electric field E (2) i [hence defined at its natural time (no + 1/2)dt], the latter is multiplied by (α + α+0) with the axion-like field living at two different times.
It is perhaps relevant to stress that, at the end, the need for all the factors multiplying within a given operators living at the same space-time site (l o , l), is a crucial aspect for determining the right functional form the lattice operator. As we saw, however, this is not enough, as S

L(2) ac
Eq. (4.13) verifies that, with (l o , l) = (n o , n), but still has problems. One must first recognize the natural space-time site where the operator lives. As we do not want to have the sum of electric fields at different times in the operator, the appropriate representation for the electric field is E (2) i , which lives in integer lattice sites, but semiinteger times, i.e. (l o , l) = (n o + 1 2 , n). This implies, correspondingly, that the appropriate magnetic field representation must be 1 2 (B Site Let us emphasize that, only thanks to the fact that α lives at (l o , l) = (n o + 1 2 , n), so that µ lives at (l o , l) = (n o , n), we can make real sense of these discrete equations. Only thanks to this interpretation, the terms within each of the above equations live at a given common natural space-time site (specifically indicated above in the rhs of each equation), around which we can expand each equation and reproduce the continuum analogue Eqs. (2.13)-(2.16), up to order O(dx 2 ).

Abelian gauge theory with an inhomogeneous axion
Let us just recall action Eq. (2.1) for a fully inhomegeneous axion-like field a(x), characterizing the Higgs sector. Variation of the action produces the EOM (imposing already the Bianchi identities) Various differences arise with respect the set of Eqs. (2.13)-(2.16) describing an homogeneous field a(x) = a(t). First, the gauge field EOM Eq. (4.23) includes now a term ∝ ∇a × E. Secondly, the Gauss law Eq. (4.24) includes a term ∝ ∇a · B. Third, the axion dynamics follows a (sourced) wave equation, where c 2 s now plays the role of a real propagation speed. In other words, there are new terms affecting the dynamics (except in the higgs sector) due to the spatial dependence of a(x). Action Eq. (4.21) and the corresponding set of Eqs. (4.22)-(4.25), cannot be used anymore to describe the problem of an Abelian gauge theory with a chemical potential. They rather describe an Abelian gauge theory in the presence to a fully inhomogeneous axion field. We shall understand that from now on we are dealing with such scenario.

Implicit scheme for real time evolution
In order to find a lattice formulation for the interaction a(x)F µν F µν ∝ a E · B where a(x) is now an inhomogeneous field, we can proceed in the same manner as in Sect. 4.1, when a(x) = a(t) was simply considered a homogeneous field. The same considerations apply now in order to find an appropriate lattice representation of the gauge-axion interaction. We can thus survey the same lattice implementations S Eq. (4.14), fail again to generate vanishing terms in the discrete EOM due to the lattice Bianchi identities. We also find that S Varying this action, we obtain a set of finite difference coupled equations,  In conclusion, it does not seem possible to find a set of discrete equations reproducing to order O(dx 2 ) the dynamics of an Abelian gauge theory in the presence of a general axion-like field a(x) = a(t, x), and at the same time being solvable by an explicit scheme. An approximate way around this difficulty can however be obtained by an implicit scheme as follows. Let us write the discrete equation evolving the electric fields, but writing down only the terms that prevented us from achieving an explicit scheme where [...] represents all the terms in the rhs of the discrete equation which involve only the amplitude (or spatial gradients) of the gauge field A i , but not the electric field. Now let us suppose that we approximate the electric field term in the rhs of Eq. (4.28) as so that using this approximation, we obtain an approximate solution to Eq. (4.28) as The solution for the updated electric field found this way makes Eq. (4.29) to reproduce the continuum Eq. (4.23) to order O(dt). We can however build now an iterative solution as (4.31) so that by successive iterations we approach closer an closer to the correct solution to Eq. (4.28) Of course, it is enough, in principle, to iterate just two times, so that E i,+0 2 E i,+0 2 2 solves Eq. (4.28) reproducing the continuum Eq. (4.23) to order O(dt 2 ). On the other hand, the Gauss law within the set of Eqs. (4.27), should be exact (up to computer machine precision) as long as E i is the exact solution to Eq. (4.28). However, as we are now approximating the electric field at each time step as E i E i | n , it is not clear a priori the accuracy attained in the (now approximated) Gauss law particularly after only n = 2 iterations. Even though successive solutions E i | n with increasingly larger n, can never be better than order O(dt 2 ) with respect E i | 2 , it might very well be the case that, in order to fulfill the Gauss law with sufficient precision, E i | n is required to an order n 2. We have not investigated explicitly this aspect in simulations, as this will depend most likely on the specific scenario under study. We leave therefore this check as a future task to be considered when applying our formalism into specific scenarios where a time-dependent and fully-inhomogeneous axion may play a central role.
As a last comment, let us note that in relevant cosmological scenarios like e.g. axioninflation [52][53][54][55][56][57][58][59], gauge fields are largely excited due to their axionic-coupling to a shiftsymmetric field a(x) which plays the role of the inflaton, and hence is (mostly) homogeneous. Only when the gauge fields are largely excited towards the end of inflation or during preheating, will they back react into the axion field, breaking its (classical) homogeneity. For most of the dynamics the axion-inflaton field remains therefore almost homogeneous. It is therefore conceivable that one may solve the system of Eqs. (4.27) simply using Eq. (4.30) to solve for the electric fields, i.e. with E i E i | 1 , and yet maintain a good accuracy close to O(dt 2 ). The reason for this is that even though the approximated terms in the rhs of Eq. (4.30) have a reduced accuracy of O(dt) instead of O(dt 2 ), they are also suppressed by ∇α. Thus, these terms may be negligible in the dynamics, allowing to solve iteratively the set of Eqs. (4.27) together with Eq. (4.30) to advance the electric fields, yet with O(dt 2 ) accuracy. As only dedicated simulations can resolve this issue, we leave as future work the test of this circumstance within these scenarios.

Lattice Chern-Simons number(s)
Let us now consider the definition of the Chern-Simons number in the continuum theory (5.1) In the lattice, we can define an equivalent quantity describing this topological number, by considering some lattice representation of E · B, and substituting the space-time integral by finite sums, d 4 x → ∆t∆x 3 no, n . In order to do this, we just need to follow a similar logic as in Sect. 4, when we disccused the different lattice representations of an axioniccoupling. For example, we want that the lattice representations of E and B live at the same space-time site, so that we expand the lattice version of E · B around a common site, in order to reproduce the continuum limit to a given order O(dx n ). Obviously, this prevent us from simply using our ordinary representation of the lattice electric and magnetic fields, , as these expressions reproduce their continuous analogues to order O(dx 2 ) only when interpreting that E i and B i live at different space-time sites, l = (n o + 1 2 , n) and l = (n o , n + 1 2ĵ + 1 2k ), respectively. We could consider to expand these lattice fields around the common space-time lattice site (n o , n), but then the topological density built like that, would only reproduce the continuum limit to linear order O(dx). As we have already derived the lattice EOM reproducing the system continuum dynamics up to order O(dx 2 ), we should clearly aim now for a description of the Chern-Simons density to (at least) order O(dx 2 ).
Given our experience in Sect. 4 in building up an axionic-type coupling to gauge fields up to order O(dx 2 ), a natural candidate to describe the Chern-Simons density at a common space-time site, is naturally given by i E (4) i B (4) i , which reproduces the continuum density E · B around l = (n o , n) to order O(dx 2 ). Therefore, we propose as a lattice candidate to describe the Chern-Simons number in the lattice, the following expression A well known identity in the continuum is where the additive constant is simply given by the initial value We can easily demonstrate this identity using integration by parts and the fact that fields vanish at infinity: Therefore, a good starting point to check whether Eq. (5.2) describes well a Chern-Simons number, would be to see if it verifies the analogous property to Eq. (5.3) in the lattice. Of course, in the lattice we cannot represent an infinite volume, but rather we take periodic boundary conditions to mimic this. Hence, our demonstration should rely on the use of periodic boundary conditions, which in any case we have implicitly assumed in previous derivations like e.g. the discrete EOM Eqs. (4.20).
Let us define A which naturally lives at l = (n o , n). Let p, q to be non-negative integer numbers. We can now observe that due to the periodic boundary conditions, the following property holds Essentially, inside a lattice sum n i , one can always exchange A (2) i,+p0 B (4) i,+q0 by A (2) i,+q0 B (4) i,+p0 , precisely thanks to the periodicity of the boundary conditions. Thanks to this property, it is easy to check that where the sum over time steps has been explicitly divided in p steps, and the additive constant is given by the fields' initial value i . We have just demonstrated the following lattice identity Our lattice expression for the Chern-Simons number Eq. (5.2) successfully passes the first check, as it verifies the identity Eq. (5.8), which clearly represents the analogue to the continuum identity Eq. (5.3).

Chern-Simons number and chemical potential
Let us consider now an Abelian gauge theory with chemical potential µ. As a consequence of integrating in time the EOM for the chemical potential Eq. (2.16), and using the Chern-Simons number definition Eq. (5.1), the following relation follows in the continuum Therefore, a relevant test we can impose over our lattice implementation of the Chern-Simons number Eq. (5.8), is to see whether it verifies some lattice version of the relation Eq. (5.10). From the discrete evolution equation for the chemical potential, see Eqs. (4.20), we see that the chemical potential after p time steps, is given by cs , (5.11) where we denote the initial chemical potential value as µ o ≡ µ(n o = 0), and the lattice volume as V L ≡ (N dx) 3 , with N is the number of lattice sites per dimension. Eq. (5.11) has lead us in fact to define a new lattice Chern-Simons number as i .  Eq. (5.2) can be exactly determined. With a bit of algebra, we find Note that we have used the fact that the A i,+p0 − A (2) i,+(p−1)0 is simply 5 of order O(dt), 5 By construction A (2) i,+p0 − A (2) i,+(p−1)0 ≡ dtEi((p − 1 2 )0, n + 1 2î ), so even though it changes in time, this change is not cumulative, and represents always a small perturbation.
whereas the initial constant −dtE (2) i B (4) i is, of course, of order O(dt) by construction. Therefore, we see that both lattice Chern-Simon numbers 'track' each other in time, and their relative (density) difference is always a small O(dt) perturbation. We have in fact numerically checked the relation Eq. (5.14) in lattice simulations, and found that it is verified exactly (to machine precision). In conclusion, even though both lattice Chern-Simons numbers represent a valid description, given that n cs .
Let us also note that using the trick expressed by Eq. (5.6), we can also write the lattice axionic-interaction term in the case of an Abelian gauge theory with chemical potential, in an alternative way to Eq. (4.19), like where we have discarded additive constants in the last expression of Eq. (5.15), as they do no contribute to the EOM when varying the action. One can check, of course, that from variations of the last expression in Eq. (5.15), we also obtain identical EOM as in Eq. (4.20), as it should be. Writing S L α(t) as in the last expression in Eq. (5.15), is perhaps the most natural thing to do in the presence of a chemical potential (homogeneous axion), as µ ≡ ∆ − o α is simply given by the time derivative of an auxiliar variable α, but the latter plays no dynamical role. The last expression in Eq. (5.15) eliminates precisely the presence of α in the action, and maintains only terms involving explicitly µ = ∆ − o α.

Derivative representation of Pontryagin densityF µν F µν
Our former demonstration(s) that the Chern-Simons number(s) can be written in either form ∼ no, n E B or as ∼ n A B, indicates that we can find a derivative representation of our lattice expressions of the Pontragyn density. Whereas in the continuum we can write the identity Q ≡F µν F µν = ∂ µ K µ , with K µ the Chern-Simons current, we expect that in the discrete some analogous relation may exist, so that Q = ∆ + µ K µ . In order to see this, let us focus on the Chern-Simons number n L(2) cs (analogous derivations can be applied to n L(1) cs ). If we redo the steps in Eq. (5.12), we immediately realize that the lattice Chern-Simons number n L(2) cs can be re-written like and Note that this does not mean that we can locally substitute the expression for the lattice Pontryagin density Q L by ∆ + o K o , like if it was an identity at every lattice site. However, whenever summing over the lattice volume, we can do such a replacement locally inside the argument of the lattice sum, i.e.
If we also sum over the time history of the system, we then arrive immediately at , (5.21) and immediately recognize the first term in the rhs of Eq. (5.21) as equal to ∆ + 0 K 0 . On the other hand, the second term in the rhs of Eq. (5.21), can be re-written like where we have identified We have demonstrated therefore that we can write as a local identity at every lattice site, with K µ defined by components, K 0 given by Eq. (5.18) and K i given by Eq. (5.23). Given that we consider periodic boundary conditions, let us note that it must be true that n i ∆ + i K i = 0, no matter the expression for K i . Therefore, the identity in Eq. (5.19) will not change in any case, if we just substitute This implies that the expression for the Chern-Simons number Eq. (5.20), given only in terms of K 0 (and not K i 's), is of course unchanged, as it should, independently of the expression Eq. (5.23) we found for K i . We have found, as promised, a lattice expression for the Pontryagin density, Eq. (5.24), that admits a total derivative representation, with K 0 and K i given by Eqs. (5.18), (5.23), respectively.

Chern-Simons number in the presence of a magnetic field
Let us now turn our attention into the case where a background magnetic field is present in the system. Following [72], we can introduce a magnetic flux in the lattice by demanding that the boundary conditions of the gauge fields A i are not periodic. Without loss of generality, we can consider an external magnetic field in theẑ direction, B = (0, 0, B). Such magnetic field can be introduced in the lattice, by demanding that only the component A 1 (n 1 , n 2 , n 3 ) of the gauge field is 'aperiodic' at a given x-site, say n 1 = 1, at the boundary of the lattice y-axis, independently of its location within the z-axis, i.e. dx (A j (n 1 , 0, n 3 ) − A j (n 1 , N, n z )) = 2πn mag δ 1j δ 1n 1 . (5.27) This condition creates a flux of magnitude Φ flux =≡ 2πn mag orthogonal to the xy-plane of the lattice, with area A ≡ (N dx) 2 , In principle, any flux can be generated. However, quantization of this flux is required for maintaining a periodic action in the non-compact formulation of a gauge theory, see [72] for details. We need therefore to take n mag as a positive integer, with n mag = 0 simply representing the absence of magnetic field.

Summary and Discussion
In this paper we have derived a lattice representation of an axionic-interaction a(x)F µν F µν , presenting step by step the necessary ingredients to achieve a formulation that i) reproduces the continuum limit to order O(dx 2 ), ii) it is consistent with the (lattice version of the) Bianchi identities, and iii) it is solvable by an iterative scheme of evolution. We first considered in Sect. 4.1 the case of a homogeneous axion a(x) = a(t), deriving a lattice representation of the axion-gauge interaction that leads to a set of discrete couple equations solvable by an explicit local iterative scheme, see Eqs. (4.20). We generalized our results afterwards, in Sect. 4.2, to the case of a fully inhomogeneous axion a(x) = a(t, x). We showed that the set of discrete lattice Eqs. (4.27) do not admit a simple local explicit solution (while preserving the O(dx 2 ) difference with respect to the continuum). We have proposed an implicit scheme to overcome this difficulty. We have also introduced consistent lattice formulation(s) of the Chern-Simons number n cs ∝ d 4 x Q (Sect. 5) based on the lattice version(s) of Q =F µν F µν developed in Sect. 4. We put special care in the need to achieve a lattice formulation that admits a total derivative representation for Q = ∆ + µ K µ . We showed explicitly that such total derivative representation exists, and provided the expression for the K µ components, see Eqs. (5.18), (5.23). Finally, we derived the analogous lattice expressions for the Chern-Simons number in the presence of an external magnetic field.
A number of potential applications of our lattice formalism has been already mentioned in the Introduction. In particular, in an accompanying paper [69] we study a number of questions one can address in high temperature electrodynamics with non-zero background magnetic field and fermion chemical potential. This theory was formulated on the lattice in Section 4.1. We investigate in [69] the random walk of the topological charge and show that it has a diffusive behavior in the presence of magnetic field B. This indicates that the mechanism for fermionic number non-conservation for B = 0 is similar to that in non-Abelian gauge theories. The diffusion rate is related to the rate of chiral charge nonconservation, and we present new results concerning the value of this rate. Our formulation allows us to study the dynamics of instabilities in the presence of non-zero µ and elucidate the role of thermal fluctuations of the gauge fields. We find several interesting behaviors that we did not expect a priori.
In addition, let us note that our formalism can be useful as well for the study of the nonlinear dynamics in axion-inflation models [52][53][54][55][56][57]. In these scenarios, gauge fields coupled to a pseudo-scalar inflaton, are excited to high occupation states. Towards the last stages of inflation the system becomes non-linear: the gauge fields are so largely excited, that they significantly back-react into the inflaton dynamics, affecting also inflationary expansion rate. The large amplification of the gauge fields leads to a very efficient generation of gravitational waves and scalar density perturbations, both with non-Gaussian statistics. If the amplitude of the scalar perturbations sourced by the gauge fields is too large, primordial black holes (PBHs) in excess to the current bounds may appear [56]. The details of both the gravitational waves and scalar perturbations (possibly leading to PBHs), depend very sensitively on the late non-linear stages of inflation, where analytical techniques can only provide an order of magnitude estimation of the dynamics. Our formalism, however, is suitable for solving the complicate non-linear dynamics numerically on a lattice, deviating from the continuum (classical) theory only to order O(dx 2 ). Besides, in axion-inflation scenarios, preheating is driven by the so called tachyonic resonance of the gauge fields, which occurs precisely due to the axionic-coupling with the inflaton [58,59]. This represents a complicate non-linear evolution stage after inflation, for which our lattice formalism can be particularly suitable.