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

TL;DR

This paper develops a gauge-invariant lattice formulation for Abelian gauge theories coupled with a shift-symmetric axion-like field, enabling accurate real-time simulations of topological effects relevant in particle physics and cosmology.

Contribution

It introduces a lattice implementation that preserves gauge invariance and shift symmetry, with a novel topological density definition and an iterative scheme for inhomogeneous fields.

Findings

Lattice formulation reproduces continuum limit to order (dx_^2)

Defines a topological density with a lattice total derivative representation

Enables simulation of real-time dynamics of topological charge and axion interactions

Abstract

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.

Equations

$$Q={\partial }_{\mu }{K}^{\mu },$$

$$Q={\Delta }_{\mu }^{+}{K}^{\mu },$$

$$S$$

$$D_{\mu }{D}^{\mu }\phi$$

$${\partial }_{\nu }{F}^{\mu \nu }-\frac{e^2}{4{\pi }^2}\frac{a}{M}{\partial }_{\nu }{\tilde{F}}^{\mu \nu }$$

$${\partial }_0{\partial }_{0}a-c_s^2{\partial }_i{\partial }_{i}a$$

$$j_{\mu }=(\rho ,\vec{J})\equiv (2\mathrm{Im}\{{\varphi }^{\ast }\dot{\varphi }\},2\mathrm{Im}\{{\varphi }^{\ast }\vec{D}\varphi \})$$

$$D_o{D}_o\phi -\vec{D}\vec{D}\phi$$

$$\dot{\vec{E}}+\vec{\nabla }\times \vec{B}+\frac{e^2}{4{\pi }^2}\frac{a}{M}(\dot{\vec{B}}-\vec{\nabla }\times \vec{E})$$

$$\vec{\nabla }\vec{E}+\frac{e^2}{4{\pi }^2}\frac{a}{M}\vec{\nabla }\vec{B}$$

$$\ddot{a}-c_s^2{\vec{\nabla }}^{2}a$$

$$\vec{\nabla }a=0.$$

$$a\equiv \alpha M,\dot{a}\equiv \mu M,$$

$$S$$

$$D_o{D}_o\phi -D_j{D}_j\phi$$

$$\dot{\vec{E}}+\vec{\nabla }\times \vec{B}$$

$$\vec{\nabla }\vec{E}$$

$$\dot{\mu }$$

$$\dot{\mu }=\frac{3}{{\pi }^2}\frac{1}{T^2}\frac{1}{V}{\int }_{{}_{V\to \infty }}{d}^{3}x\vec{E}\cdot \vec{B}.$$

$$M=T(\mathrm{temperature}),$$

$$c_s^2=12.$$

$$\displaystyle\begin{array}[]{rclr}\Delta_{\mu}^{\pm}\phi\equiv{\pm 1\over dx}(\phi_{\pm\mu}-\phi)&\rightarrow&(\partial_{\mu}\varphi)(x)+O(dx^{2})&\,,~{}x\equiv(n\pm{1\over 2}\hat{\mu})dx\,,\vspace*{0.2cm}\\ (D_{\mu}^{\pm}\varphi)\equiv{\pm 1\over dx}(U_{\pm\mu}\varphi_{\pm\mu}-\varphi)&\rightarrow&(D_{\mu}\varphi)(x)\left(1-{1\over 2}idxA_{\mu}(x)\right)&,~{}x\equiv(n\pm{1\over 2}\hat{\mu})dx\,,\end{array}$$

$$\phi (n)⟶e^{+i\beta (n)}\phi (n),A_{\mu }(n+\frac{1}{2}\hat{\mu })⟶A_{\mu }(n+\frac{1}{2}\hat{\mu })+{\Delta }_{\mu }^{+}\beta (n+\frac{1}{2}\hat{\mu }),$$

$$U_{\pm \mu ,n}⟶e^{i\beta }{U}_{\pm \mu ,n}{e}^{-i{\beta }_{\pm \mu }},D_{\mu }^{\pm }\phi ⟶e^{i\beta }{D}_{\mu }^{\pm }\phi .$$

$$S_{\mathrm{AH}}^L$$

$${\Delta }_o^{-}{\Delta }_o^{+}\phi -\sum _i{D}_i^{-}{D}_i^{+}\phi +V_{,{\phi }^{\ast }}=0$$

$${\Delta }_o^{-}{\Delta }_o^{+}(A_i)-\sum _j\left({\Delta }_j^{-}{\Delta }_j^{+}(A_i)-{\Delta }_i^{+}{\Delta }_j^{-}(A_j)\right)=2e^2\mathrm{Im}\{{\phi }^{\ast }{D}_i\phi \}$$

$$\sum _i{\Delta }_i^{-}{\Delta }_o^{+}(A_i)=2e^2\mathrm{Im}\{{\phi }^{\ast }{\Delta }_{+o}\phi \}(\mathrm{Gauss}\mathrm{Constraint})$$

$$\displaystyle|D_{\mu}^{+}\varphi|^{2}(l)\big{|}_{{l\equiv n+{\hat{\mu}\over 2}}}$$

$$\displaystyle\left(\Delta_{\mu}^{+}A_{\nu}-\Delta_{\nu}^{+}A_{\mu}\right)^{2}(l)\big{|}_{{l\equiv n+{\hat{\mu}\over 2}+{\hat{\nu}\over 2}}}$$

$$\displaystyle(D_{\mu}^{-}D_{\mu}^{+}\varphi)(l)\big{|}_{l\equiv n}$$

$$\displaystyle\Delta_{o}^{-}\Delta_{o}^{+}A_{i}(l)\big{|}_{l\equiv n+{\hat{i}\over 2}}$$

$$\displaystyle\Delta_{j}^{-}\Delta_{j}^{+}A_{i}(l)\big{|}_{l\equiv n+{\hat{i}\over 2}}$$

$$\displaystyle\Delta_{i}^{+}\Delta_{j}^{-}A_{j}(l)\big{|}_{l\equiv n+{\hat{i}\over 2}}$$

$$\displaystyle{\rm Im}\{\varphi^{*}D_{i}^{+}\varphi\}(l)\big{|}_{l\equiv n+{\hat{i}\over 2}}$$

$$S_{ac}\equiv \frac{1}{4{\pi }^2}\int d^{4}x\alpha \vec{E}\vec{B}.$$

$$S_{ac}^{L(1)}\propto \sum _{\vec{n},n_o}\alpha \sum _i{E}_i{B}_i=\sum _{\vec{n},n_o}\alpha \sum _i({\Delta }_o^{+}{A}_i-{\Delta }_i^{+}{A}_o){ϵ}_{ijk}{\Delta }_j^{+}{A}_k.$$

Full text

1 Introduction

Real-time evolution of classical fields has many applications in different areas of high energy physics and cosmology. These include the creation and evolution of topological defects in the early Universe Grigoriev:1988bd ; Vincent:1997cx ; Rajantie:1999mp ; Davis:2000kv ; Hindmarsh:2001vp ; Hindmarsh:2000kd ; Bevis:2006mj ; Figueroa:2012kw ; Daverio:2015nva ; Hindmarsh:2017qff , non-perturbative investigations of hot sphaleron transitions related to fermion number non-conservation in the electroweak theory Grigoriev:1989je ; Ambjorn:1990pu ; Ambjorn:1990wn ; Moore:1998swa ; Moore:1999fs ; GarciaBellido:1999sv ; Moore:2000mx ; GarciaBellido:2003wd ; Tranberg:2003gi ; Tranberg:2006dg ; DOnofrio:2014rug , the analysis of inflationary preheating Khlebnikov:1996mc ; Prokopec:1996rr ; Felder:2000hj ; Rajantie:2000nj ; GarciaBellido:2002aj ; DiazGil:2007dy ; DiazGil:2008tf ; Figueroa:2015rqa ; Enqvist:2015sua ; Figueroa:2016wxr , the generation of cosmological perturbations Bassett:1998wg ; Bassett:1999mt ; Bassett:1999ta ; Finelli:2000ya ; Chambers:2007se ; Bond:2009xx and gravitational waves Khlebnikov:1997di ; Easther:2006gt ; GarciaBellido:2007dg ; GarciaBellido:2007af ; Dufaux:2007pt ; Dufaux:2008dn ; Figueroa:2011ye ; Bethke:2013aba ; Bethke:2013vca ; Figueroa:2016ojl during preheating, and different aspects of quark-gluon plasma, see e.g. Attems:2016ged and references therein. The classical approximation to quantum dynamics should work when the relevant distance scales of the problem exceed considerably the typical quantum distances, at finite temperatures $T$ given by $1/T$ (distance between particles) and $1/(gT)$ (Debye screening length, where $g$ is the gauge coupling).

The numerical procedures for Abelian or non-Abelian scalar-gauge theories, such as the Standard Model is well developed, see e.g. Ambjorn:1990pu 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 Bodeker:1999gx ). 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 $Pontryagin$ density, $Q=F_{\mu\nu}\tilde{F}^{\mu\nu}$, where the dual of the field strength is defined as usual by $\tilde{F}^{\mu\nu}\equiv{1\over 2}\epsilon^{{\mu\nu}{\alpha\beta}}F_{{\alpha\beta}}$, with $\epsilon^{{\mu\nu}{\alpha\beta}}$ being the completely anti-symmetric tensor in four dimensions, with $\epsilon^{0123}\equiv 1$. The most interesting examples contain an axion field coupled linearly to $Q$ as $a(x)F_{\mu\nu}\tilde{F}^{\mu\nu}$, and the theories with non-zero chemical potential $\mu$ for chiral fermions, leading to an effective bosonic Hamiltonian containing the Chern-Simons term for the gauge fields $\mu n_{\rm cs}$, with $n_{\rm cs}\propto\int 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 space-time can be expressed via an integral over the boundary. To put it in other words, in continuum one can write

[TABLE]

where $K^{\mu}$ is the Chern-Simons current. The lattice analogue of this relation is

[TABLE]

where $\Delta_{\mu}^{+}$ is the lattice difference in positive direction of $\mu$ 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). For example, it is Eq. (1) which makes the axion mass $m_{a}$ to be zero in all orders of perturbation theory ($m_{a}\neq 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) is still approximate, but as precise as possible Moore:1999fs .

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. (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 (2) on the lattice Luscher:1981zq ; Woit:1983tq . The applications of our formulation can include, for example, the study of the late stage of inflation and preheating in axion-inflation models Freese:1990rb ; Sorbo:2011rz ; Cook:2011hg ; Barnaby:2011qe ; Linde:2012bt ; Pajer:2013fsa ; Adshead:2015pva ; Adshead:2016iae , the clarification of the role of Standard Model hypercharge group U(1) in baryogenesis and in magnetic field generation Joyce:1997uy ; Giovannini:1997eg ; Long:2013tha ; Long:2016uez , the modeling the different aspects of chiral magnetic effects Vilenkin:1980fu ; Joyce:1997uy ; Alekseev:1998ds ; Fukushima:2008xe , or the study of the problem of chiral fermionic charge evolution in high temperature electrodynamics Boyarsky:2011uy ; Boyarsky:2012ex . An accompanying paper Figueroa:2017tbf 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)\tilde{F}_{\mu\nu}F^{\mu\nu}$, 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,{\bf x})$ in Sect. 4.2. In Sect. 5 we finally discuss the lattice formulation of the Chern-Simons number $n_{\rm cs}\propto\int d^{4}x\,Q$ based on the lattice version(s) of $Q=\tilde{F}_{\mu\nu}F^{\mu\nu}$ 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. (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-time111We choose a metric signature $\eta_{\mu\nu}=\eta^{\mu\nu}=(-,+,+,+)$., in the presence of an axion-type field $a(x)$ linearly coupled to the $Pontryagin$ density $F_{\mu\nu}\tilde{F}^{\mu\nu}$ of a U(1) gauge field,

[TABLE]

We consider a ’Higgs’ sector as $\mathcal{L}_{\varphi}=(D_{0}\varphi)^{*}(D_{0}\varphi)-(\vec{D}\varphi)^{*}(\vec{D}\varphi)+V(\varphi^{*}\varphi)$, with $\varphi=(\varphi_{1}+i\varphi_{2})/\sqrt{2}$ a $U(1)$ charged field [with $\varphi_{i}\in\Re$], $A_{\mu}=(\phi,\vec{A})$ the gauge field, $e^{2}$ the gauge coupling strength, and $D_{\mu}\equiv\partial_{\mu}-iA_{\mu}$ the covariant derivative. The field strength and its dual are defined as usual by $F_{\mu\nu}\equiv\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ and $\tilde{F}^{\mu\nu}\equiv{1\over 2}\epsilon^{{\mu\nu}{\alpha\beta}}F_{{\alpha\beta}}$. Given our choice of $A_{\mu}\equiv(\phi,\vec{A})$, we define the electric and magnetic fields as $E^{i}=E_{i}=\dot{A}_{i}-\partial_{i}\phi$ and $B^{i}=B_{i}=\epsilon_{ijk}\partial_{j}A_{k}$. This leads to the relations $E^{i}\equiv F_{0i}=-F^{0i}$, $B^{i}\equiv{1\over 2}\epsilon^{ijk}F_{jk}$, so that $F_{{\mu\nu}}=(\delta_{\mu 0}\delta_{\nu i}-\delta_{\mu i}\delta_{\nu 0})E^{i}+(\delta_{\mu i}\delta_{\nu j}-\delta_{\mu j}\delta_{\nu i})\epsilon_{ijk}B^{k}$. Given our metric signature, we obtain $F_{{\mu\nu}}\tilde{F}^{{\mu\nu}}=+4\vec{E}\vec{B}$, and arrive at the final vectorial expressions of Eq. (3).

Lagrangian Eq. (3) 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_{s}^{2}$ as a free parameter, as this will be convenient for us later on. Action Eq. (3) is invariant under the transformations $\varphi(x)\rightarrow e^{+i\beta(x)}\varphi(x)$, $A_{\mu}(x)\rightarrow A_{\mu}(x)+\partial_{\mu}\beta(x)$, with $\beta(x)\in\Re$ and $e^{i\beta(x)}\in U(1)$. Varying the action, one obtains the equations of motion (EOM)

[TABLE]

where the (unit-charge) current is defined as $j^{\mu}=2{\rm Im}\{\varphi^{*}D^{\mu}\varphi\}$, so that

[TABLE]

Eqs. (4)-(6) can be rewritten in a vectorial form as

[TABLE]

with Eq. (10) representing the Gauss constraint in the presence of an axion.

As mentioned in Section 1, it is well known that the Pontryagin density represents a topological term, as it can be written as a total derivative $F_{\mu\nu}\tilde{F}^{\mu\nu}=\partial_{\mu}K^{\mu}$. This is reflected by the Bianchi identities, i.e. the term $\partial_{\nu}\tilde{F}^{{\mu\nu}}=0$ in Eq. (5), or equivalently its vectorial counterparts $(\dot{\vec{B}}-\vec{\nabla}\times E)=0$ in Eq. (9), and $\vec{\nabla}\vec{B}=0$ in Eq. (10). Those terms simply represent vanishing contributions in the EOM, so it is customary to remove them. It is nonetheless convenient for us to keep such terms in the EOM, despite their null contribution. The reason for this will become clear, however, only in Sec. 4, after we introduce the lattice discretization scheme(s) for the action Eq. (3).

For the time being, let us simply note now that due to the topological nature of the $F_{\mu\nu}\tilde{F}^{\mu\nu}$ operator, action Eq. (3) is also (’topologically’) invariant under $a(x)\rightarrow a(x)+C$, with $C$ an arbitrary constant. This is reflected in the fact that the linear coupling of $a(x)$ to the Pontryagin density $\propto\int d^{4}x\,a(x)\,F_{\mu\nu}\tilde{F}^{\mu\nu}$, represents in reality a derivative coupling: thanks to the total derivative nature of $F_{\mu\nu}\tilde{F}^{\mu\nu}=\partial_{\mu}K^{\mu}$, after integration by parts, we obtain $\propto\int d^{4}x\,K^{\mu}\partial_{\mu}a(x)$. Once the Bianchi identities are considered, the terms proportional to $a(x)$ in Eq. (5) [equivalently in Eqs. (9), (10)] disappear, as it must, for a derivative coupling. As we will see later in Sec. 4, the lattice equivalent of the Bianchi identities, and hence the topological nature of the lattice equivalent of $F_{\mu\nu}\tilde{F}^{\mu\nu}$, depend crucially on the lattice discretization scheme. The lattice equivalent terms to $\propto\vec{\nabla}\vec{B}$, $\propto(\dot{\vec{B}}-\vec{\nabla}\times E)$ in the discrete EOM, are actually not granted to vanish by default. Achieving such a goal will represent, in fact, a guiding principle towards the construction of a correct lattice formulation of $F_{\mu\nu}\tilde{F}^{\mu\nu}$, admitting a total derivative representation on the lattice.

Let us remark that the $shift$ symmetry enjoyed by $a(x)$, together with the linear nature of its coupling to the Pontryagin density, is at the heart of the original introduction of the axion as a solution for the strong CP problem222Of course the QCD axion is rather coupled to Tr $\tilde{G}_{\mu\nu}G^{\mu\nu}$, where $G_{\mu\nu}$ is the gluon field strength.. Once $a(x)$ is considered as a dynamical field, the interaction $a(x)F_{\mu\nu}\tilde{F}^{\mu\nu}$ in Eq. (3) leads naturally to extra contributions in the EOM. Assuming a fully space-time dependent axion field, naturally leads to a contribution in Eq. (5) as $\propto\tilde{F}^{{\mu\nu}}\partial_{\nu}a$, or equivalently by the vectorial counterparts $\propto(-{\dot{a}}\vec{B}+{\vec{\nabla}a}\times\vec{E})$ and $\propto{\vec{\nabla}a}\cdot\vec{B}$ in Eq. (9) and Eq. (10), respectively. Besides, the Pontryagin density $F_{\mu\nu}\tilde{F}^{\mu\nu}$ acts as a source for $a(x)$ in the $rhs$ of Eq. (11). 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. (3) with $a(x)=0$.

2.1 Mimicking a chemical potential

Interestingly, through an adequate interpretation of field variables and parameters, Eq. (3) can be mapped into the description of a gauge theory with a chemical potential $\mu$ for chiral fermions (for more details see Figueroa:2017tbf ). Starting from Eq. (3) 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

[TABLE]

We then introduce the following convenient ’dimensionally reduced’ variables

[TABLE]

so that Eq. (3) is reduced to

[TABLE]

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_{s}^{2}$, see Eq. (21). We will describe first, however, the new dynamical equations that follow from minimizing the new re-written action.

Varying Eq. (14) we obtain the equations of motion of the system, which we write directly in a vectorial form as

[TABLE]

Once again, let us note that the terms $\propto\alpha(\dot{\vec{B}}-\vec{\nabla}\times E)$ in Eq. (9) and $\propto\alpha\vec{\nabla}\vec{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_{s}^{2}$ to appropriate values, so that the set of Eqs. (15)-(18) properly describe an Abelian gauge theory with the chemical potential $\mu$ for chiral fermionic charge. In particular, the EOM of a chemical potential follows from anomaly equation Adler:1969gk ; Bell:1969ts and in our case has the form Boyarsky:2012ex

[TABLE]

In light of Eq. (18), we must identify the mass scale with the temperature of the system,

[TABLE]

and fix the dimensionless parameter to

[TABLE]

Let us note that in the original action Eq. (3), $c_{s}^{2}$ 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_{s}^{2}$ represents simply a parameter in the theory. One should not conclude therefore that the value given by Eq. (21) represents a super-luminal axion speed of propagation, as in the chemical potential context where such value is determined, $c_{s}^{2}$ does not represent, in first place, a propagation speed. Eqs. (15)-(18) together with Eqs. (20)-(21), describe appropriately an Abelian gauge with chemical potential $\mu$.

3 Part I. Lattice formulation of Abelian gauge theories to order $\mathcal{O}(dx_{\mu}^{2})$, Scalar-Electrodynamics

In this section we briefly summarize the basics of the lattice formulation of a gauge theory. We focus, for convenience, in the non-compact formulation of an Abelian gauge theory (for related discussion see, e.g. Kajantie:1998rz ). 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.

3.1 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},\vec{n})=(n_{o},n_{1},n_{2},n_{3})$ displaced in the $\mu-$direction by one unit lattice spacing, $n+\hat{\mu}$, will be often referred simply as $n+\mu$ or by $+\mu$. For example, $\varphi_{+\mu}\equiv\varphi(n+\hat{\mu})$. Components of gauge fields live in between lattice sites in the direction of the component, so $A_{\mu}\equiv A_{\mu}(n+{1\over 2}\hat{\mu})$, $A_{\mu,+\nu}\equiv A_{\mu}(n+{1\over 2}\hat{\mu}+\hat{\nu})$, etc. For simplicity of the notation, we will refer to both the lattice spacing $\Delta x$ and the time step $\Delta t$, simply as $dx$, so if we write e.g. $+\hat{\mu}dx$, this should be interpreted as a time step advancement $+\Delta t$ if $\mu=0$, or as a unit displacement $+\Delta x$ in a given spatial direction $\mu=1,2$ or $3$. We will loosely speak of the lattice spacing order $\mathcal{O}(dx)$, independently or whether we are referring to $\mathcal{O}(\Delta t)$ or $\mathcal{O}(\Delta x)$.

We define a lattice link, as usual, like $U_{\mu}\equiv U_{\mu}(n+{1\over 2}\hat{\mu})\equiv e^{-i\int_{x(n)}^{x(n+\hat{\mu})}A_{\mu}(x^{\prime})dx^{\prime\mu}}$ $\simeq e^{-idx^{\mu}A_{\mu}(n+{1\over 2}\hat{\mu})}$. We also define $U_{-\mu}\equiv U_{\mu,-\mu}^{*}\equiv U_{\mu}^{*}(n-{1\over 2}\hat{\mu})\simeq e^{+idx^{\mu}A_{\mu}(n-{1\over 2}\hat{\mu})}$. Forward (+) and backward (-), ordinary and covariant derivatives, are defined in the lattice by

[TABLE]

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

[TABLE]

with $\beta$ an arbitrary function, so that the links and covariant derivatives transform as

[TABLE]

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

[TABLE]

from where the lattice gauge invariance (based on the transformations defined above) is rather explicit. We will refer to Eq. (27) as the $Abelian-Higgs$ (AH) lattice action. As we will show in Section 3.2, this action reproduces to order $\mathcal{O}(dx^{2})$ the continuum action Eq. (3) in the absence of an axion field ($a(x)=0$). Varying Eq. (27) with respect to the different fields, we obtain the lattice equivalent of the dynamical equations, which read (taking the Coulomb Gauge $A_{o}=0\leftrightarrow U_{o}=1$)

[TABLE]

3.2 Recovering the continuum limit to order $\mathcal{O}(dx^{2})$

Let us note that, given a continuum action $S=\int dt\,d^{3}x\,\mathcal{L}_{C}$, with lagrangian given by the sum of various operators $\mathcal{L}_{C}=\sum_{p}\mathcal{O}_{C}^{(p)}$, we try to emulate the same physical system by defining a lattice action $S=\Delta t\Delta x^{3}\sum_{\vec{n},n_{o}}\mathcal{L}_{L}$, with the lattice lagrangian given by the sum of lattice operators $\mathcal{L}_{L}=\sum_{p}\mathcal{O}_{L}^{(p)}$. Each operator $\mathcal{O}_{L}^{(p)}$, when expanded around the lattice site where it is defined, reproduces the continuum operator $\mathcal{O}_{C}^{(p)}$ to some $n$-th order in the lattice spacing, $\mathcal{O}_{L}^{(p)}\rightarrow\mathcal{O}_{C}^{(p)}+\mathcal{O}(dx^{n})$. For consistency we require that each and every operator $\mathcal{O}_{L}^{(p)}$ reproduces the corresponding continuum term $\mathcal{O}_{L}^{(p)}$ to the same $n$-th order.

Let us inspect in detail action Eq. (27). Each term in the action reproduces the equivalent continuum terms to second order in $dx$, i.e. $\mathcal{O}_{\rm L}^{(p)}\rightarrow\mathcal{O}_{\rm C}^{(p)}+\mathcal{O}(dx^{2}),~{}\forall p$. In order to verify this, it is crucial to make a suitable interpretation of the lattice sites where each operator naturally ’lives’. For instance, let us consider a lattice operator mimicking the kinetic term of a scalar field, $\mathcal{O}_{L}^{(k)}={1\over 2}(\Delta_{\mu}^{+}\phi)^{2}$. If we were to expand $\Delta_{\mu}^{+}\varphi$ around the position $x=n\cdot dx$, we would obtain $(\Delta_{\mu}^{+}\phi)(n)\rightarrow(\partial_{\mu}\phi)({\vec{x}})+\mathcal{O}(dx)$. Consequently $\mathcal{O}_{L}^{(k)}\rightarrow\mathcal{O}_{C}^{(k)}+\mathcal{O}(dx)$, so the continuum limit is only reproduced to linear order. The natural site where $\mathcal{O}_{L}^{(k)}$ really lives is however, not $x=ndx$, but rather in $x_{\mu/2}\equiv(n+{1\over 2}\hat{\mu})dx$, as it involves a finite difference of the field evaluated (with equal weight) at both $n$ and $n+\hat{\mu}$ lattice sites. If we expand $\mathcal{O}_{L}^{(k)}$ around ${x}_{\mu/2}$, then we obtain $(\Delta_{\mu}\phi)(n+{1\over 2}\hat{\mu})\rightarrow(\partial_{\mu}\phi)({x}_{\mu/2})+\mathcal{O}(dx^{2})$. 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 $\mathcal{O}_{L}^{(p)}$ in the lattice lagrangian $\mathcal{L}_{L}=\sum_{p}\mathcal{O}_{L}^{(p)}$, may (possibly) live at different lattice sites. There is however no contradiction in this result, for as long as each operator $\mathcal{O}_{L}^{(p)}$ 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. (27), 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. (27) reproduces the continuum to second order in $dx$, i.e. $\mathcal{O}_{\rm L}^{(p)}\rightarrow\mathcal{O}_{\rm C}^{(p)}+\mathcal{O}(dx^{2}),~{}\forall 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 $\vec{x}\equiv ldx$ the physical coordinate where the continuum limit is reproduced. A simple expansion leads to

[TABLE]

Let us turn our attention now into the lattice EOM Eqs. (28)-(30). Let us begin by the EOM of the charged scalar field, which involve the terms $V_{,\varphi^{*}}$ and $D_{\mu}^{-}D_{\mu}^{+}\varphi$. We can expand each term around $l=n$, as the operator $D_{\mu}^{+}\varphi$ lives naturally at $n+{1\over 2}\hat{\mu}$, but $D_{\mu}^{-}(D_{\mu}^{+}\varphi)$ makes it live back to $l=(n+{1\over 2}\hat{\mu})-{1\over 2}\hat{\mu}=n$,

[TABLE]

Consistently, the term $V_{,\varphi^{*}}$ in Eq. (30) [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. (30) reproduces correctly the corresponding continuum Eq. (15) up to $\mathcal{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+{\hat{j}\over 2}$ for the dynamical equation Eq. (29). This can be easily shown by expanding each term of this equation

[TABLE]

A similar analysis can be done expanding the terms of Eq. (28) [Gauss constraint] around their natural lattice site $l=n+{1\over 2}\hat{0}$, leading to the same conclusion: when interpreting that the lattice operators involved in the discrete equations live at $l=n+{1\over 2}\hat{0}$ [Gauss constraint Eq. (28)] or $l=n+{1\over 2}\hat{i}$ [Dynamical Eq. (29)], the discrete EOM of the gauge fields reproduce correctly the continuum equations Eqs. (16), (17) [in the absence of an axion] up to $\mathcal{O}(dx^{2})$ corrections.

4 Part II. Lattice formulation of Abelian gauge theories to order $\mathcal{O}(dx_{\mu}^{2})$, 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, $S_{ac}\equiv{1\over(4\pi)^{2}}\int d^{4}x{a\over M}F_{{\mu\nu}}\tilde{F}^{{\mu\nu}}$. 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 $\alpha(x)=\alpha(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 $\alpha(x)=\alpha(t,{\bf 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 $\alpha\equiv{a\over M}$, $\mu\equiv\dot{\alpha}$ defined in Eq. (13), 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_{\mu}\equiv(\phi,\vec{A})$, we find $F_{{\mu\nu}}\tilde{F}^{{\mu\nu}}=+4\vec{E}\vec{B}$, so that we can write the above continuum action as

[TABLE]

We will refer to the interaction described by Eq. (38) as an $axionic-coupling$. Our main aim now is to formulate a latttice version of the continuum action Eq. (38), from which to derive a discrete version of the EOM in the continuum Eqs. (8)-(11) [or Eqs. (15)-(18) in the case of a homogeneous axion mimicking a chemical potential]. In light of the EOM in the continuum, we can foresee three possible problems arising when formulating a lattice version of Eq. (38):

1. i)

The terms $\alpha(-\dot{\vec{B}}+\vec{\nabla}\times\vec{E})$ and $\alpha\vec{\nabla}{\vec{B}}$ in Eqs. (16), (17), vanish in the continuum thanks to the Bianchi identities $\partial_{\nu}\tilde{F}^{\mu\nu}=0$, which are equivalent to $\dot{\vec{B}}=\vec{\nabla}\times\vec{E}$ and $\vec{\nabla}{\vec{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=\tilde{F}_{\mu\nu}\tilde{F}^{\mu\nu}$ is topological admitting a total derivative representation $Q=\Delta_{\mu}^{+}K^{\mu}$. 2. ii)

Assuming that a correct version of the Bianchi identities follows naturally from a given topological lattice formulation of ${F}_{\mu\nu}\tilde{F}^{\mu\nu}$, another problem may arise. The terms $\propto\mu\vec{B}$ and $\propto\vec{E}\vec{B}$ in Eqs. (9), (11), indicate possible obstructions to achieving an explicit scheme to solve iteratively the set of lattice coupled equations reproducing the continuum Eqs. (8)-(11). 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. 3. iii)

When considering a fully inhomogeneous axion-like field $a(x)$, terms proportional to spatial variations $\propto\nabla a(x)$ appear in the EOM. In particular, the term $\nabla a\times\vec{E}$ in the $rhs$ of electric field evolution equation Eq. (9) $\dot{E}_{i}=[...]+{e^{2}\over 4\pi^{2}}(\nabla a\times\vec{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 $i$th 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 $\nabla a\times\vec{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 $\alpha\,\vec{E}\vec{B}$, determining their ’appropriateness’ based on the ability of each lattice formulation to address the previous criteria $i)-iii)$.

4.1 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 $\alpha(x)=\alpha(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 $\mu=\dot{\alpha}$. 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,{\bf x})$, which we will address in Sect. 4.2, building up from our previous findings on the homogeneous case.

4.1.1 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, $E_{i}\equiv(\Delta_{o}^{+}A_{i}-\Delta_{i}^{+}A_{o})$ and $B_{i}\equiv\epsilon_{ijk}\Delta_{j}^{+}A_{k}$, like

[TABLE]

A priori, this looks like benign lattice operator, since it describes the continuum action up to order $\mathcal{O}(dx^{2})$, as we showed already in Sect. 3: each term of $F_{{\mu\nu}}\equiv(\Delta_{\mu}^{+}A_{\nu}-\Delta_{\nu}^{+}A_{\nu})$ individually reproduces the continuum expression up to to order $\mathcal{O}(dx^{2})$, when interpreting that $F_{{\mu\nu}}$ lives at $n+{1\over 2}\hat{\mu}+{1\over 2}\hat{\nu}$. When varying this lattice action with respect $A_{i}$, we obtain (recall that we are assuming now $\alpha(x)=\alpha(t)$ homogeneous) a term $\propto\alpha[\Delta_{o}^{-}B_{i}-(\nabla^{-}\times\vec{E})_{i}]$. When varying with respect to $A_{o}$, we obtain a term $\propto\alpha\sum\Delta_{i}^{-}B_{i}$. As expected, these terms resemble the continuum analogues $\propto\alpha(\partial_{o}{\vec{B}}-\vec{\nabla}\times\vec{E})$ and $\propto\alpha\vec{\nabla}{\vec{B}}$ in Eqs. (16), (17). However, contrary to the continuum analogues, they do not vanish. The reason is simple, the correct discrete version of $\vec{\nabla}{\vec{B}}$ is rather $\sum_{i}\Delta_{i}^{+}B_{i}=0$, simply because the lattice magnetic field was defined in terms of forward derivatives, $B_{i}\equiv\epsilon_{ijk}\Delta_{j}^{+}A_{k}$. To obtain the desired result, one needs to take the divergence over $B_{i}$ with a forward derivative: as $\Delta_{i}^{+}$ and $\Delta_{j}^{+}$ commute, $\Delta_{i}^{+}\Delta_{j}^{+}A_{k}$ is symmetric in $i\leftrightarrow j$, and hence its contraction with $\epsilon_{ijk}$ vanishes, $\sum_{i}\Delta_{i}^{+}B_{i}=\epsilon_{ijk}\Delta_{i}^{+}\Delta_{j}^{+}A_{k}=0$. However, the variation of Eq. (39) naturally led to $\sum\Delta_{i}^{-}B_{i}=\epsilon_{ijk}\Delta_{i}^{-}\Delta_{j}^{+}A_{k}\neq 0$, which does not vanish because $\Delta_{i}^{-}\Delta_{j}^{+}A_{k}$ is not symmetric in $i\leftrightarrow j$, given that $\Delta_{i}^{-}$ and $\Delta_{j}^{+}$ do not commute. Similarly, the term $[\Delta_{o}^{-}B_{i}-(\vec{\nabla}^{-}\times\vec{E})_{i}\neq 0$] does not vanish. The appropriate version in the lattice of the Bianchi identity should rather be built as $[(\vec{\Delta}^{+}\times\vec{E})_{i}-\Delta^{+}_{o}B_{i}]=0$, as one can easily check that $\Delta_{o}^{+}B_{i}\equiv\epsilon_{ijk}\Delta_{j}^{+}\Delta_{o}^{+}A_{k}$ = $\epsilon_{ijk}\Delta_{j}^{+}(\Delta_{o}^{+}A_{k}-\Delta_{k}^{+}A_{o})$ $\equiv(\vec{\Delta}^{+}\times\vec{E})_{i}$. Variation of the lattice action Eq. (39) produced instead the expression $[(\vec{\nabla}^{-}\times\vec{E})_{i}-\Delta_{o}^{-}B_{i}]$ = $\epsilon_{ijk}(\Delta_{j}^{-}\Delta_{o}^{+}-\Delta_{o}^{-}\Delta_{j}^{+})A_{k}\neq 0$, which does not vanish, simply because $\Delta_{j}^{-}\Delta_{o}^{+}$ and $\Delta_{o}^{-}\Delta_{j}^{+}$ do not commute.

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. (39) as a lattice operator is actually not consistent. Note that even though it reproduces the continuum result to order $\mathcal{O}(dx^{2})$, it consists however in the product of three fields, $\alpha,E_{i}$ and $B_{i}$ that live, not only in different lattice sites, also at different time steps. Whereas $\alpha$ lives at $(n_{o},\vec{n})$, $E_{i}$ lives at $(n_{o}+{1\over 2},\vec{n}+{1\over 2}\hat{i})$ and $B_{i}$ lives at $(n_{o},\vec{n}+{1\over 2}\hat{j}+{1\over 2}\hat{k})$. 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. (39) 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

[TABLE]

and note that each of these fields can be expressed as

[TABLE]

For convenience we also define

[TABLE]

which reproduce the continuum as

[TABLE]

A ’symmetrized’ operator that reproduces the continuum expression of $S_{ac}$ at $l=(\vec{n},n_{o}+{1\over 2})$ to order $\mathcal{O}(dx^{2})$, can be easily proposed based on the above expressions,

[TABLE]

Varying Eq. (4.1.1) with respect $A_{i}$, produces a term in the discrete EOM of the gauge field as $\alpha\left(\sum_{j,k}\epsilon_{ijk}(\Delta^{+}_{j}+\Delta^{-}_{j})E^{(8)}_{k}\right.$ $-\left.(\Delta_{o}^{+}+\Delta_{o}^{-})B_{i}^{(8)}\right)$, which resembles the continuum term $\alpha(\epsilon_{ijk}\partial_{j}E_{k}-\dot{B}_{i})$ in Eq. (16). Whereas the latter vanishes thanks to the Bianchi identity in the continuum $\vec{\nabla}\times\vec{E}=\dot{\vec{B}}$, it can be shown, with a bit of algebra, that a lattice version of this identity holds as $\sum_{j,k}\epsilon_{ijk}(\Delta^{+}_{j}+\Delta^{-}_{j})E^{(8)}_{k}$ = $(\Delta_{o}^{+}+\Delta_{o}^{-})B_{i}^{(8)}$. This implies that the new term encountered in the gauge field discrete EOM just vanishes. Similarly, when varying Eq. (4.1.1) with respect $A_{o}$, we produce a new term $\alpha\sum_{i}(2+dt\Delta^{+}_{t})\Delta_{i}^{-}B_{i}^{(8)}$, which again resembles the term in the continuum $\alpha\partial_{i}B_{i}$ in Eq. (17), which vanishes due to the Bianchi identity $\vec{\nabla}\vec{B}=0$. With a bit of algebra, it can be shown that the analogous lattice identity reads $\Delta_{i}^{-}B_{i}^{(8)}=(\Delta_{i}^{+}+\Delta_{i}^{-})B_{i}^{(4)}=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.1.1) represents a good lattice candidate from which to derive [together with action Eq. (27)] 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.

4.1.2 Explicit scheme for real time evolution

The operator Eq. (4.1.1) proposed to represent an axionic coupling, exhibits various features: $i)$ it reproduces correctly the continuum term to order $\mathcal{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. (15)-(18) in the continuum. In fact, varying Eq. (4.1.1) with respect to the gauge fields, generates a term $\propto{1\over 2}\left(\mu_{-0}B_{i}^{(8)}+\mu B_{i,+0}^{(8)}\right)$, which reduces correctly to the continuum term $\mu\vec{B}$ in Eq. (16), to order $\mathcal{O}(dx^{2})$. At the same time, varying Eq. (4.1.1) with respect to $\alpha(t)$, generates a term $\propto{1\over 2}(E_{i}^{(2)}+E_{i,-0}^{(2)})B_{i}^{(4)}$ sourcing the chemical potential, which again reduces correctly to the continuum source term $\propto\vec{E}\vec{B}$ in Eq. (18), to order $\mathcal{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 $\mu$ and $A_{i,+0}$ (to obtain $B_{i,+0}^{(8)}$), and at the same time to find $A_{i,+0}$ and $\mu$ we need $E_{i}$. One would need to express the term $\mu B_{i,+0}^{(8)}$ in the gauge field EOM in terms of $E_{i}$, and then solve for $E_{i}$, but this would complicate the equations unnecessary333Actually, it is not even feasible, in principle, to solve this system, as e.g. $B_{1,+0}^{(8)}$ depends $A_{2,+\hat{1}+\hat{0}}$ and $A_{3,+\hat{1}+\hat{0}}$, which depend respectively on $E_{2,+1}$ and $E_{3,+1}$. So at the end, the equation for the $\hat{1}$-component of the electric field depends on $E_{1}$ but also on $E_{2,+1}$ and $E_{3,+1}$, and analogously for the other electric field component equations. Therefore, one cannot just simply solve for the $E_{i}$ 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.. A simpler solution consists in choosing another lattice operator that upon variation over the gauge fields, prevents the duplication of the $\sim\mu B_{i}$ terms at the two different times. The duplication of these terms originated from Eq. (4.1.1) due to the fact that the electric field $E_{i}^{(4)}$ multiplying $\alpha B_{i}^{(4)}$ in the operator, is equivalent to the sum of electric fields at two time steps, $E_{i}^{(4)}\equiv{1\over 2}(E_{i}^{(2)}+E_{i,-0}^{(2)})$. Hence, when varying with respect to $A_{i}$, two terms of the form $\sim\mu B_{i}$ are generated, $\left(\mu_{-0}B_{i}^{(8)}+\mu B_{i,+0}^{(8)}\right)$, 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

[TABLE]

This operator fails however in ’symmetrizing’ the expression for $\vec{E}\cdot\vec{B}$ around a common time: as $B_{i}^{(4)}$ lives at integer times $n_{o}$ whereas $E_{i}^{(2)}$ lives at time semi-integer timEs $(n_{o}+{1\over 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_{ac}^{L(1)}$ Eq. (39).

A better solution is found if the lattice operator $\sim\alpha(E\cdot B)$ is built such that each element $\alpha$, $\vec{E}$ and $\vec{B}$, live separately in a given common space-time site. In other words, contrary to the operator $S_{ac}^{L(3)}$ Eq. (4.1.2), where the electric and magnetic fields lived at different time steps separated away by half step ${dt\over 2}$, now both $\vec{E}$ and $\vec{B}$ must live at the same time step. As the natural time where electric fields live is $(n_{o}+{1\over 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,

[TABLE]

where $\equiv$ 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 $\vec{\nabla}\vec{B}$ and $(\dot{\vec{B}}-\vec{\nabla}\times\vec{E})$ in the discrete equations of motion due to the lattice Bianchi identities, and generate the following terms [recall that here we still consider $\alpha(x)=\alpha(t)$],

[TABLE]

From here we see that $S_{ac}^{L(6)}$ Eq. (4.1.2) does not allow for an explicit scheme of iteration444In this occasion this occurs because even though there is a single electric field $E_{i}^{(2)}$ [hence defined at its natural time $(n_{o}+1/2)dt$], the latter is multiplied by $(\alpha+\alpha_{+0})$ with the axion-like field living at two different times., similarly as what it happened with $S_{ac}^{L(2)}$ Eq. (4.1.1). Only actions $S_{ac}^{L(4)}$ Eq. (52) and $S_{ac}^{L(5)}$ Eq. (53) allow allow for an explicit scheme of iteration, and in fact produce essentially an equivalent set of coupled discrete EOM.

In conclusion, although all operators $S_{ac}^{L(1)}$-$S_{ac}^{L(6)}$ reproduce correctly the continuum gauge-axionic interaction to order $\mathcal{O}(dx^{2})$, $S_{ac}^{L(1)}$ Eq. (39) and $S_{ac}^{L(3)}$ Eq. (4.1.2) fail to generate vanishing terms in the discrete EOM from the lattice Bianchi identities, whereas $S_{ac}^{L(2)}$ Eq. (4.1.1) and $S_{ac}^{L(6)}$ Eq. (4.1.2) produce correctly null terms due to the lattice Bianchi identities, but fail to generate explicit iterative schemes for solving the set of coupled finite difference EOM. Only $S_{ac}^{L(4)}$ Eq. (52) and $S_{ac}^{L(5)}$ Eq. (53), produce correctly vanishing terms due to the lattice Bianchi identities, while maintaining an explicit iterative scheme for solving the set of lattice EOM. Actually, the functional form of $S_{ac}^{L(4)}$ Eq. (52) [equivalently $S_{ac}^{L(5)}$ Eq. (53)] indicates us something important: the natural time steps where an shift-symmetric field lives are also semi-integer times, and not integer steps as for ordinary scalar fields. Pseudo-scalar fields (independently of whether they are homogeneous as we have just required so far) must live in semi-integer times in the lattice. In light of this, we can interpret $S_{ac}^{L(4)}$ Eq. (52) and $S_{ac}^{L(5)}$ Eq. (53) 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

[TABLE]

where it is important to note that both $\alpha$ and $A_{o}$ live at space-time sites $(n_{o}+{1\over 2},\vec{n})$, whereas $A_{i}$ live at $(n_{o},\vec{n}+{\hat{i}\over 2})$. Consequently, $\mu\equiv\Delta_{o}^{-}\alpha$ lives naturally at $(n_{o},n)$, $E_{i}$ at $(n_{o}+{1\over 2},\vec{n}+{1\over 2}\hat{i})$, $E_{i}^{(2)}$ at $(n_{o}+{1\over 2},\vec{n})$, whereas $B_{i}$ lives at $(n_{o},n+{1\over 2}\hat{j}+{1\over 2}\hat{k})$, and $B_{i}^{(4)}$ at $(n_{o},\vec{n})$.

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},\vec{l})$, is a crucial aspect for determining the right functional form the lattice operator. As we saw, however, this is not enough, as $S_{ac}^{L(2)}$ Eq. (4.1.1) verifies that, with $(l_{o},\vec{l})=(n_{o},\vec{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_{i}^{(2)}$, which lives in integer lattice sites, but semi-integer times, i.e. $(l_{o},\vec{l})=(n_{o}+{1\over 2},\vec{n})$. This implies, correspondingly, that the appropriate magnetic field representation must be ${1\over 2}(B_{i}^{(4)}+B_{i,+0}^{(4)})$, so that it also lives in $(l_{o},\vec{l})=(n_{o}+{1\over 2},\vec{n})$. The linear part in the (pseudo-)scalar field $\alpha$, can then be made meaningfully defined only if $\alpha$ also lives at $(l_{o},\vec{l})=(n_{o}+{1\over 2},\vec{n})$.

The final set of equations of motion obtained from varying $S_{\rm AH}^{L}+S_{{\alpha(t)}}^{L}$ [Eq. (27) + Eq. (60)], mimicking a system with chemical potential $\mu\equiv\Delta_{o}^{-}\alpha$ at finite temperature $T$, where recall that $M=T$ and $c_{s}^{2}=12$, are (we use the Coulomb gauge $A_{o}=0$, so that $U_{o}=1$)

[TABLE]

Let us emphasize that, only thanks to the fact that $\alpha$ lives at $(l_{o},l)=(n_{o}+{1\over 2},n)$, so that $\mu$ lives at $(l_{o},\vec{l})=(n_{o},\vec{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. (15)-(18), up to order $\mathcal{O}(dx^{2})$.

4.2 Abelian gauge theory with an inhomogeneous axion

Let us just recall action Eq. (3) for a fully inhomegeneous axion-like field $a(x)$,

[TABLE]

with $\mathcal{L}_{\varphi}\equiv(D_{0}\varphi)^{*}(D_{0}\varphi)-(\vec{D}\varphi)^{*}(\vec{D}\varphi)+V(\varphi^{*}\varphi)$ characterizing the Higgs sector. Variation of the action produces the EOM (imposing already the Bianchi identities)

[TABLE]

Various differences arise with respect the set of Eqs. (15)-(18) describing an homogeneous field $a(x)=a(t)$. First, the gauge field EOM Eq. (71) includes now a term $\propto{\vec{\nabla}a}\times\vec{E}$. Secondly, the Gauss law Eq. (72) includes a term $\propto{\vec{\nabla}a}\cdot\vec{B}$. Third, the axion dynamics follows a (sourced) wave equation, where $c_{s}^{2}$ 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. (69) and the corresponding set of Eqs. (70)-(73), 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.

4.2.1 Implicit scheme for real time evolution

In order to find a lattice formulation for the interaction $a(x)\tilde{F}_{\mu\nu}F^{\mu\nu}\propto a\vec{E}\cdot\vec{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_{ac}^{L(1)}$-$S_{ac}^{L(6)}$ proposed in Sec. 4.1, as we know they all reproduce correctly the continuum interaction $a(x)\tilde{F}_{\mu\nu}F^{\mu\nu}$ to order $\mathcal{O}(dx^{2})$. The spatial dependence of $a(x)$ does not change the fact that $S_{ac}^{L(1)}$ Eq. (39) and $S_{ac}^{L(3)}$ Eq. (4.1.2), fail again to generate vanishing terms in the discrete EOM due to the lattice Bianchi identities. We also find that $S_{ac}^{L(2)}$ Eq. (4.1.1) and $S_{ac}^{L(6)}$ Eq. (4.1.2) produce correctly null terms due to the lattice Bianchi identities, but fail again to generate explicit iterative schemes for solving the set of coupled finite different EOM. The expressions $S_{ac}^{L(4)}$ Eq. (52) and $S_{ac}^{L(5)}$ Eq. (53) produce correctly vanishing terms in the gauge field discrete EOM due to the lattice Bianchi identities. However, contrary to the homogeous $a(x)=a(t)$ case, $S_{ac}^{L(4)}$-$S_{ac}^{L(5)}$ do not lead now to an explicit iterative scheme for solving the set of lattice EOM when $a(x)=a(t,{\bf x})$ is inhomogeous.

Let us see this more in detail, considering $S_{ac}^{L(4)}$ Eq. (52) as the lattice representation of the axion-gauge field interaction [$S_{ac}^{L(5)}$ Eq. (53) is essentially equivalent]. We write the final lattice action mimicking to order $O(dx^{2})$ the continuum action Eq. (69) as

[TABLE]

Varying this action, we obtain a set of finite difference coupled equations,

[TABLE]

which reproduce to order $\mathcal{O}(dx^{2})$ the set of continuum Eqs. (70)-(73), when expanding each equation around its natural lattice site. Note that the natural sites ascribed to each discrete EOM in Eqs. (83) coincide with those listed in the $rhs$ of Eqs. (68), so we do not repeat them here. A simple inspection of the $lhs$ of each equation suffices anyways to determine these sites, e.g. knowing that the lattice representation of the electric field $E_{i}=\Delta_{o}^{+}A_{i}$ lives at $l=(n_{o}+{1\over 2},\vec{n}+{1\over 2}\hat{i})$, then equation $E_{i}\equiv[...]$ must be expanded around $x=((n_{o}+{1\over 2})dt,(\vec{n}+{1\over 2}\hat{i})dx)$, $\Delta_{o}^{-}E_{i}=[...]$ around $x=(n_{o}dt,(\vec{n}+{1\over 2}\hat{i})dx)$, etc.

Unfortunately, the set of finite difference coupled Eqs. (83) cannot be solved with an explicit scheme. In fact, these equations can only be formally solved by a non-local solution, even though we started from a local Lagrangian Eq. (4.2.1). The reason for this is the following. Say we consider the $\hat{1}$-component of the curl product involving the electric field in the $rhs$ of the gauge field EOM $\Delta_{o}^{-}E_{1}=[...]$, i.e. $\sim[(\Delta a)\times\vec{E}^{(2)}]_{1}$. This term depends on $E_{2,+1}$ and $E_{3,+1}$, so the equation to update the $\hat{1}$-component of the electric field depends on $E_{2,+1}$ and $E_{3,+1}$, and analogously for the other electric field component equations. In other words, there is a first-neighbour coupling of the electric field components. This simply makes impossible to solve explicitly for the $E_{i}$ components at a given lattice site (in a given time step), as we would need to know the electric field components at all other lattice sites (at the new time step).

The origin of this problem lies of course in the form of the interaction $\sim a\vec{E}\vec{B}$. In the moment $a(x)=a(\vec{x},t)$ is inhomogeneous, this brings up the term $\vec{\nabla}a\times\vec{E}$ in the $rhs$ of the gauge field continuum EOM Eq. (71), $\dot{\vec{E}}=[...]$. As in the lattice $E_{i}$ lives at $\vec{n}+{1\over 2}\hat{i}$, the equivalent lattice expression representing the continuum term $\sim(\vec{\nabla}a\times\vec{E})$ must be evaluated at both $\vec{n}$ and $\vec{n}+\hat{i}$. This explains in fact the $(2+dx\Delta^{+}_{i})$ operation in $rhs$ of the equivalent discrete EOM $\Delta_{o}^{-}E_{i}=[...]$ within Eqs. (83). This turns non-local the set of coupled equations in finite differences Eqs. (83), hence making it unfeasible to solve them iteratively by an explicit scheme at each lattice site.

In conclusion, it does not seem possible to find a set of discrete equations reproducing to order $\mathcal{O}(dx^{2})$ the dynamics of an Abelian gauge theory in the presence of a general axion-like field $a(x)=a(t,{\bf 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

[TABLE]

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.2.1) as

[TABLE]

so that using this approximation, we obtain an approximate solution to Eq. (4.2.1) as

[TABLE]

The solution for the updated electric field found this way makes Eq. (85) to reproduce the continuum Eq. (71) to order $\mathcal{O}(dt)$. We can however build now an iterative solution as

[TABLE]

so that by successive iterations we approach closer an closer to the correct solution to Eq. (4.2.1)

[TABLE]

Of course, it is enough, in principle, to iterate just two times, so that $E_{i,+{\hat{0}\over 2}}\simeq\left.E_{i,+{\hat{0}\over 2}}\right|_{2}$ solves Eq. (4.2.1) reproducing the continuum Eq. (71) to order $\mathcal{O}(dt^{2})$. On the other hand, the Gauss law within the set of Eqs. (83), should be exact (up to computer machine precision) as long as $E_{i}$ is the exact solution to Eq. (4.2.1). However, as we are now approximating the electric field at each time step as $E_{i}\simeq\left.E_{i}\right|_{n}$, it is not clear a priori the accuracy attained in the (now approximated) Gauss law

[TABLE]

particularly after only $n=2$ iterations. Even though successive solutions $\left.E_{i}\right|_{n}$ with increasingly larger $n$, can never be better than order $\mathcal{O}(dt^{2})$ with respect $\left.E_{i}\right|_{2}$, it might very well be the case that, in order to fulfill the Gauss law with sufficient precision, $\left.E_{i}\right|_{n}$ is required to an order $n\gg 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. axion-inflation Freese:1990rb ; Sorbo:2011rz ; Cook:2011hg ; Barnaby:2011qe ; Linde:2012bt ; Pajer:2013fsa ; Adshead:2015pva ; Adshead:2016iae , gauge fields are largely excited due to their axionic-coupling to a shift-symmetric 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. (83) simply using Eq. (4.2.1) to solve for the electric fields, i.e. with $E_{i}\simeq\left.E_{i}\right|_{1}$, and yet maintain a good accuracy close to $\mathcal{O}(dt^{2})$. The reason for this is that even though the approximated terms in the $rhs$ of Eq. (4.2.1) have a reduced accuracy of $\mathcal{O}(dt)$ instead of $\mathcal{O}(dt^{2})$, they are also suppressed by $\vec{\nabla}\alpha$. Thus, these terms may be negligible in the dynamics, allowing to solve iteratively the set of Eqs. (83) together with Eq. (4.2.1) to advance the electric fields, yet with $\mathcal{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.

5 Lattice Chern-Simons number(s)

Let us now consider the definition of the Chern-Simons number in the continuum theory

[TABLE]

In the lattice, we can define an equivalent quantity describing this topological number, by considering some lattice representation of $\vec{E}\cdot\vec{B}$, and substituting the space-time integral by finite sums, $\int d^{4}x\rightarrow\Delta t\Delta x^{3}\sum_{n_{o},\vec{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 axionic-coupling. For example, we want that the lattice representations of $\vec{E}$ and $\vec{B}$ live at the same space-time site, so that we expand the lattice version of $\vec{E}\cdot\vec{B}$ around a common site, in order to reproduce the continuum limit to a given order $\mathcal{O}(dx^{n})$. Obviously, this prevent us from simply using our ordinary representation of the lattice electric and magnetic fields, $E_{i}\equiv(\Delta_{o}^{+}A_{i}-\Delta_{i}^{+}A_{o})$ and $B_{i}\equiv(\Delta_{j}^{+}A_{k}-\Delta_{k}^{+}A_{j})$, as these expressions reproduce their continuous analogues to order $\mathcal{O}(dx^{2})$ only when interpreting that $E_{i}$ and $B_{i}$ live at different space-time sites, $l=(n_{o}+{1\over 2},\vec{n})$ and $l=(n_{o},\vec{n}+{1\over 2}\hat{j}+{1\over 2}\hat{k})$, respectively. We could consider to expand these lattice fields around the common space-time lattice site $(n_{o},\vec{n})$, but then the topological density built like that, $\sum_{i}E_{i}B_{i}=\sum_{i,j,k}\epsilon_{ijk}(\Delta_{o}^{+}A_{i}-\Delta_{i}^{+}A_{o})\Delta_{j}^{+}A_{k}$, would only reproduce the continuum limit to linear order $\mathcal{O}(dx)$. As we have already derived the lattice EOM reproducing the system continuum dynamics up to order $\mathcal{O}(dx^{2})$, we should clearly aim now for a description of the Chern-Simons density to (at least) order $\mathcal{O}(dx^{2})$.

Given our experience in Sect. 4 in building up an axionic-type coupling to gauge fields up to order $\mathcal{O}(dx^{2})$, a natural candidate to describe the Chern-Simons density at a common space-time site, is naturally given by $\sum_{i}E_{i}^{(4)}B_{i}^{(4)}$, which reproduces the continuum density $\vec{E}\cdot\vec{B}$ around $l=(n_{o},n)$ to order $\mathcal{O}(dx^{2})$. Therefore, we propose as a lattice candidate to describe the Chern-Simons number in the lattice, the following expression

[TABLE]

A well known identity in the continuum is

[TABLE]

where the additive constant is simply given by the initial value

[TABLE]

We can easily demonstrate this identity using integration by parts and the fact that fields vanish at infinity: $\int d^{4}x~{}\vec{E}\cdot\vec{B}$ = $\int d^{4}x\,(\dot{A}_{i}-\partial_{i}\phi)\epsilon_{ijk}\partial_{j}A_{k}$ = $-\int d^{4}x\,A_{k}\epsilon_{ijk}\partial_{j}(\dot{A}_{i}-\partial_{i}\phi)$ = $\int d^{4}x\,A_{k}\dot{B}_{k}$ = \int d^{3}x\,A_{i}B_{i}\big{|}_{0}^{t}-\int d^{4}x~{}B_{k}\dot{A}_{k} = \int d^{3}x\,A_{i}B_{i}\big{|}_{0}^{t} $-\int d^{4}x\,E_{i}B_{i}$ + $\int d^{4}x\,\epsilon_{ijk}\partial_{i}\phi\partial_{j}A_{k}$, where the last term is equal to $-\int d^{4}xA_{k}\epsilon_{ijk}\partial_{i}\partial_{j}\phi=0$, and hence 2\int_{0}^{t}\,dt\int d^{3}x\,\vec{E}\cdot\vec{B}=\int d^{3}x\,\vec{A}\cdot\vec{B}\big{|}_{0}^{t}. Therefore, a good starting point to check whether Eq. (93) describes well a Chern-Simons number, would be to see if it verifies the analogous property to Eq. (94) 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. (68).

Let us define

[TABLE]

which naturally lives at $l=(n_{o},\vec{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

[TABLE]

Essentially, inside a lattice sum $\sum_{\vec{n}}\sum_{i}$, one can always exchange $A_{i,+p\hat{0}}^{(2)}B_{i,+q\hat{0}}^{(4)}$ by $A_{i,+q\hat{0}}^{(2)}B_{i,+p\hat{0}}^{(4)}$, precisely thanks to the periodicity of the boundary conditions. Thanks to this property, it is easy to check that

[TABLE]

where the sum over time steps has been explicitly divided in $p$ steps, and the additive constant is given by the fields’ initial value $\mathcal{D}_{o}\equiv-{1\over 2\Delta t}\sum_{n}\sum_{i}A_{i,-\hat{0}}^{(2)}B_{i}^{(4)}$. We have just demonstrated the following lattice identity

[TABLE]

with

[TABLE]

Our lattice expression for the Chern-Simons number Eq. (93) successfully passes the first check, as it verifies the identity Eq. (99), which clearly represents the analogue to the continuum identity Eq. (94).

5.1 Chern-Simons number and chemical potential

Let us consider now an Abelian gauge theory with chemical potential $\mu$. As a consequence of integrating in time the EOM for the chemical potential Eq. (18), and using the Chern-Simons number definition Eq. (92), the following relation follows in the continuum

[TABLE]

Therefore, a relevant test we can impose over our lattice implementation of the Chern-Simons number Eq. (99), is to see whether it verifies some lattice version of the relation Eq. (101). From the discrete evolution equation for the chemical potential, see Eqs. (68), we see that the chemical potential after $p$ time steps, is given by

[TABLE]

where we denote the initial chemical potential value as $\mu_{o}\equiv\mu(n_{o}=0)$, and the lattice volume as $V_{L}\equiv(Ndx)^{3}$, with $N$ is the number of lattice sites per dimension. Eq. (102) has lead us in fact to define a new lattice Chern-Simons number as

[TABLE]

with initial constant

[TABLE]

Notice that in order to arrive at the final expression of Eq. (103) of the form $\sim\sum_{\vec{n}}AB$ from the initial expression of the form $\sim\sum_{n_{o},\vec{n}}EB$, we have applied again the trick expressed by Eq. (97). We conclude that in order to correctly mimic the continuum relation Eq. (101), we are forced to define a second Chern-Simons number by Eq. (103), which is similar to that of Eq. (99), but still different. As a matter of fact, we have already shown in Eq. (103) that this new Chern-Simons verifies a lattice analogue to the continuum relation Eq. (94) $\int dx^{4}~{}\vec{E}\cdot\vec{B}={1\over 2}\int dx^{3}\vec{A}\cdot\vec{B}~{}+~{}\mathcal{C}_{o}$. A simple check also shows that $n_{\rm cs}^{L(2)}$ is built from a Chern-Simons lattice density which reproduces exactly the continuum expression $\vec{E}\cdot\vec{B}$ up to order $\mathcal{O}(dx^{2})$.

A Chern-Simons number given by Eq. (103) is therefore, not only equally valid as Eq. (99), but also preferable. Contrary to Eq. (99), it is exactly linearly proportional to the chemical potential at arbitrary times, see Eq. (102), as expected in order to mimic the continuum relation Eq. (101). Besides, as our final choice for the axion-gauge interaction is based on the lagrangian density $E_{i}^{(2)}(B_{i}^{(4)}+B_{i,+0}^{(4)})$, see Eq. (4.2.1), there is no surprise whatsoever that such expression [and not $E_{i}^{(4)}B_{i}^{(4)}$ on which Eq. (99) is based on] defines the natural Chern-Simons number of the system: it is such density the one naturally sourcing the discrete EOM of the chemical potential, see Eqs. (68).

Does this mean that our expression $n_{\rm cs}^{L(1)}$ Eq. (93) is actually incorrect ? Actually not really, as the different between $n_{\rm cs}^{L(2)}$ Eq. (103) and $n_{\rm cs}^{L(1)}$ Eq. (93) can be exactly determined. With a bit of algebra, we find

[TABLE]

Note that we have used the fact that the $\left(A_{i,+p\hat{0}}^{(2)}-A_{i,+(p-1)\hat{0}}^{(2)}\right)$ is simply555By construction $\left(A_{i,+p\hat{0}}^{(2)}-A_{i,+(p-1)\hat{0}}^{(2)}\right)\equiv dtE_{i}((p-{1\over 2})\hat{0},\vec{n}+{1\over 2}\hat{i})$, so even though it changes in time, this change is not cumulative, and represents always a small perturbation. of order $\mathcal{O}(dt)$, whereas the initial constant $-dtE_{i}^{(2)}B_{i}^{(4)}$ is, of course, of order $\mathcal{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 $\mathcal{O}(dt)$ perturbation. We have in fact numerically checked the relation Eq. (5.1) 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_{\rm cs}^{L(2)}$ Eq. (103) sources exactly the chemical potential [whereas $n_{\rm cs}^{L(1)}$ Eq. (99) does not], we propose to work only with $n_{\rm cs}^{L(2)}$.

Let us also note that using the trick expressed by Eq. (97), 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. (60), like

[TABLE]

where we have discarded additive constants in the last expression of Eq. (106), 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. (106), we also obtain identical EOM as in Eq. (68), as it should be. Writing $S_{\alpha(t)}^{L}$ as in the last expression in Eq. (106), is perhaps the most natural thing to do in the presence of a chemical potential (homogeneous axion), as $\mu\equiv\Delta_{o}^{-}\alpha$ is simply given by the time derivative of an auxiliar variable $\alpha$, but the latter plays no dynamical role. The last expression in Eq. (106) eliminates precisely the presence of $\alpha$ in the action, and maintains only terms involving explicitly $\mu=\Delta_{o}^{-}\alpha$.

5.2 Derivative representation of Pontryagin density $\tilde{F}_{\mu\nu}F^{\mu\nu}$

Our former demonstration(s) that the Chern-Simons number(s) can be written in either form $\sim\sum_{n_{o},\vec{n}}\vec{E}\vec{B}$ or as $\sim\sum_{\vec{n}}\vec{A}\vec{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\equiv\tilde{F}_{\mu\nu}F^{\mu\nu}=\partial_{\mu}K^{\mu}$, with $K^{\mu}$ the Chern-Simons current, we expect that in the discrete some analogous relation may exist, so that $Q=\Delta_{\mu}^{+}K^{\mu}$. In order to see this, let us focus on the Chern-Simons number $n_{\rm cs}^{L(2)}$ (analogous derivations can be applied to $n_{\rm cs}^{L(1)}$). If we redo the steps in Eq. (103), we immediately realize that the lattice Chern-Simons number $n_{\rm cs}^{L(2)}$ can be re-written like

[TABLE]

where

[TABLE]

and

[TABLE]

Note that this does not mean that we can locally substitute the expression for the lattice Pontryagin density $Q_{L}$ by $\Delta_{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.

[TABLE]

If we also sum over the time history of the system, we then arrive immediately at

[TABLE]

Using the definition of $K^{o}$ in Eq. (109), we see that Eq. (5.2) coincides exactly, as it should, with the final expression of Eq. (103).

Let us find the spatial $K^{i}$ components of the lattice representation of the Chern-Simons current. Starting from the definition of $Q_{L}$ in Eq. (108), we can write

[TABLE]

and immediately recognize the first term in the $rhs$ of Eq. (112) as equal to $\Delta_{0}^{+}K^{0}$. On the other hand, the second term in the $rhs$ of Eq. (112), can be re-written like

[TABLE]

where we have identified

[TABLE]

We have demonstrated therefore that we can write

[TABLE]

as a local identity at every lattice site, with $K^{\mu}$ defined by components, $K^{0}$ given by Eq. (109) and $K^{i}$ given by Eq. (114). Given that we consider periodic boundary conditions, let us note that it must be true that $\sum_{\vec{n}}\sum_{i}\Delta_{i}^{+}K^{i}=0$, no matter the expression for $K^{i}$. Therefore, the identity in Eq. (110) will not change in any case, if we just substitute $\Delta_{o}^{+}K^{o}$ by $\sum_{\mu}\Delta_{\mu}^{+}K^{\mu}$ in the $rhs$ of such equation,

[TABLE]

This implies that the expression for the Chern-Simons number Eq. (5.2), given only in terms of $K^{0}$ (and not $K^{i}$’s), is of course unchanged,

[TABLE]

as it should, independently of the expression Eq. (114) we found for $K^{i}$.

We have found, as promised, a lattice expression for the Pontryagin density, Eq. (115), that admits a total derivative representation, with $K^{0}$ and $K^{i}$ given by Eqs. (109), (114), respectively.

5.3 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 Kajantie:1998rz , 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 $\hat{z}$ direction, $\vec{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.

[TABLE]

This condition creates a flux of magnitude $\Phi_{\rm flux}=\equiv 2\pi n_{\rm mag}$ orthogonal to the $xy$-plane of the lattice, with area $A\equiv(Ndx)^{2}$,

[TABLE]

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 Kajantie:1998rz for details. We need therefore to take $n_{\rm mag}$ as a positive integer, with $n_{\rm mag}=0$ simply representing the absence of magnetic field.

The ’twisted’ boundary condition for $A_{1}$ Eq. (118) implies that whenever we want to calculate a magnetic field $B_{3}(n_{1},n_{2},n_{3})$ at the location $(n_{1},n_{2},n_{3})=(1,N-1,n_{3})$, instead of simply computing $(\Delta_{1}^{+}A_{2}-\Delta_{2}^{+}A_{1})$, we should really make the calculation

[TABLE]

as this is equivalent to assume that the first component of the gauge field makes a ’discrete jump’ $A_{1}(1,0,n_{3})=A_{1}(1,N,n_{3})+{2\pi n_{\rm mag}\over dx}$, just as required by Eq. (118).

Let us now see how the Chern-Simons number Eq. (103) is modified in the presence of a background magnetic field, introduced in the lattice through the condition Eq. (118). Let us begin by noticing that the relation Eq. (97), used previously to relate the $\sim EB$ and $\sim AB$ expressions of the CS number [see e.g. Eq. (99) or Eq. (103)], does not hold anymore, as it relies on the (now lost) periodic boundary conditions of the gauge fields. When the gauge field follows instead the twisted boundary condition Eq. (118), Eq. (97) needs to be modified. Given some $p,q$ integer numbers, we obtain now

[TABLE]

where $\mathcal{M}_{p}$ represents a magnetic correction given by

[TABLE]

Essentially, inside a lattice sum $\sum_{\vec{n}}\sum_{i}$, one can still substitute $A_{i,+p\hat{0}}^{(2)}B_{i,+q\hat{0}}^{(4)}$ by $A_{i,+q\hat{0}}^{(2)}B_{i,+p\hat{0}}^{(4)}$, as long as outside the sum one compensates by adding the magnetic correction $\mathcal{M}_{p}$. This correction is precisely a manifestation of the loss of periodic boundary conditions Eq. (118) for the gauge field component $A_{1}$.

In light of Eq. (121), the expression of the lattice Chern-Simons number Eq. (103), changes in the presence of an external magnetic field to

[TABLE]

with $\mathcal{D}_{o}^{(2)}$ the same constant as in Eq. (104), and $\mathcal{M}_{0},\mathcal{M}_{p}$ given by Eq. (123), representing respectively initial and final (after $p$ time steps) magnetic field corrections. Equivalently,

[TABLE]

6 Summary and Discussion

In this paper we have derived a lattice representation of an axionic-interaction $a(x)\tilde{F}_{\mu\nu}F^{\mu\nu}$, presenting step by step the necessary ingredients to achieve a formulation that $i)$ reproduces the continuum limit to order $\mathcal{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. (68). We generalized our results afterwards, in Sect. 4.2, to the case of a fully inhomogeneous axion $a(x)=a(t,{\bf x})$. We showed that the set of discrete lattice Eqs. (83) do not admit a simple local explicit solution (while preserving the $\mathcal{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_{\rm cs}\propto\int d^{4}x\,Q$ (Sect. 5) based on the lattice version(s) of $Q=\tilde{F}_{\mu\nu}F^{\mu\nu}$ developed in Sect. 4. We put special care in the need to achieve a lattice formulation that admits a total derivative representation for $Q=\Delta_{\mu}^{+}K^{\mu}$. We showed explicitly that such total derivative representation exists, and provided the expression for the $K^{\mu}$ components, see Eqs. (109), (114). 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 Figueroa:2017tbf 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 Figueroa:2017tbf 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\neq 0$ is similar to that in non-Abelian gauge theories. The diffusion rate is related to the rate of chiral charge non-conservation, 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 $\mu$ 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 non-linear dynamics in axion-inflation models Freese:1990rb ; Sorbo:2011rz ; Cook:2011hg ; Barnaby:2011qe ; Linde:2012bt ; Pajer:2013fsa . 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 Linde:2012bt . 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 $\mathcal{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 Adshead:2015pva ; Adshead:2016iae . This represents a complicate non-linear evolution stage after inflation, for which our lattice formalism can be particularly suitable.

Acknowledgements

We thank Dmitry Levkov and Alberto Ramos for clarifications about iterative schemes of evolution to solve implicit chain of discrete equations of motion. This work was supported by the ERC-AdG-2015 grant 694896. The work of M.S. was supported partially by the Swiss National Science Foundation.