# Axion electrodynamics: Green’s functions, zero-point energy and optical
activity

Amedeo M. Favitta  Dipartimento di Fisica e Chimica—Emilio Segré, Università
degli Studi di Palermo, Via Archirafi 36, I-90123 Palermo, Italy  Iver H.
Brevik  Department of Energy and Process Engineering, Norwegian University of
Science and Technology, N-7491 Trondheim, Norway  Moshe M. Chaichian
Department of Physics, University of Helsinki and Helsinki Institute of
Physics, P.O. Box 64, FI-00014 University of Helsinki, Finland

(15/06/2023
To To appear in the Annals of Physics)

###### Abstract

Starting from the theory of Axion Electrodynamics, we work out the axionic
modifications to the electromagnetic Casimir energy using the Green’s function
method, both when the axion field is initially assumed purely time-dependent
and when the axion field configuration is a static domain wall, so purely
space-dependent. For the first case it means that the oscillating axion
background is taken to resemble an axion fluid at rest in a conventional
Casimir setup with two infinite parallel conducting plates, while in the
second case we evaluate the radiation pressure acting on an axion domain wall.
We extend previous theories in order to include finite temperatures. Various
applications are discussed. i) We review the theory of Axion Electrodynamics
and particularly the energy-momentum conservation in a linear dielectric and
magnetic material. We treat this last aspect by extending former results by
Brevik and Chaichian (2022) and Patkos (2022). ii) Adopting the model of the
oscillating axion background we discuss the axion-induced modifications to the
Casimir force between two parallel plates by using the Green’s function
method. iii) We calculate the radiation pressure acting on an axion domain
wall at finite temperature $T$. Our results for an oscillating axion field and
a domain wall are also useful for condensed matter physics, where some
topological materials, ”axionic topological insulators”, interact with the
electromagnetic field with a Chern-Simons interaction, like the one in Axion
Electrodynamics, and there are experimental systems analogous to time-
dependent axion fields and domain walls as the ones showed e.g. by Jiang, Q.
D., & Wilczek, F. (2019) and Fukushima et al. (2019). iv) We compare our
results, where we assume time-dependent or space-dependent axion
configurations, with the discussion of the optical activity of Axion
Electrodynamics by Sikivie (2021) and Carrol et al. (1990). We also make
comparisons with the properties of known materials, such as optically active,
chiral media and the Faraday effect.

##  I Introduction

Axions, the hypothetical particles suggested by Peccei and Quinn [1, 2, 3] in
order to solve the strong CP problem in QCD [4, 5, 6], have attracted a strong
theoretical interest because they are among some of the best known candidates
for dark matter, as first suggested in Refs.[7, 8, 9]. This topic is one of
the most important problems in cosmology, and is discussed lively in the
current literature [10, 2, 11, 12, 13, 14, 15, 16]. The axion is still a
hypothetical and undetected particle and there have been many efforts to
detect it. A most promising way to detect axions is based on its interaction
with the electromagnetic field. Axion Electrodynamics is a natural extension
of classical electrodynamics, where the axion field is coupled to the
electromagnetic field, and this interaction between the photonic field and the
axion field makes it possible to detect the axions experimentally [14, 17]. It
is worth noticing that the interaction between the electromagnetic field and
the axion field leads to modifications of the dispersion relation and the
zero-point energy of the electromagnetic field, similarly to chiral QED,
another interesting extension of quantum electrodynamics that is investigated
in Refs. [18, 19]. These last aspects are discussed in our work along with the
calculation of the axionic modified Green’s functions, where we find new
results, and compare them to the results of Ref. [20]. This work is also
related to earlier investigations in Refs. [11, 21, 22, 16, 23].

We present the basics of Axion Electrodynamics and introduce our notation in
Section II. Energy-momentum balance is considered in Section III. We develop
the basics of the Green’s function method in Axion Electrodynamics in Section
IV and use it in practice for time-dependent axion backgrounds, namely with
constant time derivative and oscillating behaviour in Section V. There we
calculate the Casimir force between two perfectly conducting parallel plates
in two time-dependent backgrounds. We treat space-dependent axion backgrounds,
namely dependent on just the space coordinate $z$, in Sections VI and VII. In
Section VI we develop the optical properties of a toy model of the axion
domain wall. The electromagnetic dispersion relations in such an axion
background are given in Subsection VI.2. General features of the Green’s
function method for domain walls are introduced in Section VII.

The Casimir force on a domain wall is obtained in Subsection VII.1. Subsection
VII.2 also deals with the calculation of the temperature-dependent
electromagnetic radiation pressure acting on axion domain walls and we exploit
our results for applications in Axion Cosmology and Condensed Matter Physics.
Section VIII compares our results with formerly known phenomena in Classical
Electrodynamics, such as the Faraday effect and optical activity in chiral
media, and former results on the optical activity in Axion Electrodynamics.
Our work on the Green’s functions for the two given types of axion backgrounds
of experimental interest are relevant for observing an axion field in the
current Universe. This is seen in the long wavelength approximation, i.e. we
assume the experimental set-up to have spatial dimension $L$ much smaller that
the de Broglie wavelength of the local axion field, so that we take
$a(t)=a_{0}\sin(\omega_{a}t)$ . If a further approximation of taking the time-
derivative of the axion field to be constant is just a pedagogical first
approximation for QCD axion, it is exact for the axion field in Weyl
semimetals. The axion field is mostly treated as a fixed background field,
since this is true at first order in the perturbation theory because there is
no back reaction of the electromagnetic ﬁeld onto axions in such an
approximation, partly analogous to what happen in the linear approximation for
weak gravitational perturbations in General Relativity [24]. We discuss this
point in Section III, where we also partly discuss the axion dynamics when
interacting with the electromagnetic field.

##  II Basic elements of Axion Electrodynamics

###  II.1 Lagrangian density of Axion Electrodynamics

In the following we adopt the metric convention $\eta_{\alpha\beta}={\rm
diag}(+1,-1,-1,-1)$ and the measure system of natural units system such that:
$\hbar=c=4\pi\epsilon_{0}=1$. The remaining unit is choosen to be energy,
measured in $\text{\,}\mathrm{e}\mathrm{V}$. We consider a pseudoscalar axion
field $a=a({\bf r},t)=a(x)$ present in the entire Universe, having a two-
photon interaction with the electromagnetic field. The total Lagrangian
density $\mathcal{L}$ describing the interaction between electromagnetic field
and axion field in the Minkowski space inside a linear dielectric and magnetic
material with dielectric permittivity $\varepsilon$ and magnetic permeability
$\mu$ is [16, 14, 11, 12, 15]:

$\mathcal{L}=-\frac{1}{4}F^{\alpha\beta}H_{\alpha\beta}+\mathcal{L}_{a}-J^{\mu}A_{\mu}+\frac{1}{4}g_{a\gamma\gamma}a(x)F_{\mu\nu}\tilde{F}^{\mu\nu},$
(1)

$\tilde{F}^{\alpha\beta}=\frac{1}{2}\epsilon^{\alpha\beta\gamma\delta}F_{\gamma\delta}$
is the dual Faraday tensor , $\epsilon^{\alpha\beta\gamma\delta}$ is the
totally antisymmetric symbol with $\epsilon^{0123}=1$, and $m_{a}$ is the mass
of the axion. The tensor $H$ is:

$H_{\alpha\beta}=\begin{pmatrix}0&-D_{x}&-D_{y}&-D_{z}\\\
D_{x}&0&H_{z}&-H_{y}\\\ D_{y}&-H_{z}&0&H_{x}\\\ D_{z}&H_{y}&-H_{x}&0\\\
\end{pmatrix},$ (2)

where $\vec{D}$ and $\vec{H}$ are the usual displacement electric and magnetic
fields that are related to electric $\vec{E}$ and magnetic $\vec{B}$ fields by
$\vec{D}=\varepsilon\vec{E}$ and $\vec{B}=\mu\vec{H}$ for a linear dielectric
and magnetic material. $J^{\mu}$ is an external classical electrical 4-current
that we present here for purpose of completeness. Finally,
$g_{a\gamma\gamma}=g_{\gamma}\frac{\alpha}{\pi}\frac{1}{f_{a}}$ where $\alpha$
is the usual fine structure constant and $f_{a}$ is the axion decay constant
whose value is only insufficiently known. Astrophysical and experimental
bounds indicate that $f_{a}\sim
10^{9}-10^{12}\,$\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$ [14].

In most cases in the following, we will use the non-dimensional quantities
$\Theta(x)=g_{a\gamma\gamma}a(x)$ and $\theta(x)=a(x)/f_{a}$. The axion mass
is related to the decay constant $f_{a}$ by [5]

$m_{a}\simeq\frac{\sqrt{m_{u}m_{d}}}{m_{u}+m_{d}}\frac{f_{\pi}m_{\pi}}{f_{a}}\simeq
6\times
10^{-6}\;$\text{\,}\mathrm{e}\mathrm{V}$\left(\frac{10^{12}\;$\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$}{f_{a}}\right),$
(3)

where $m_{\pi}\simeq$140\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$ is the pion
mass and $f_{\pi}\simeq$93\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$ is the
pion decay constant. $g_{\gamma}$ is a model-dependent constant of order unity
[14].

###  II.2 QCD Axion Lagrangian, domain walls and long wavelength approximation

$\mathcal{L}_{a}$ is the ”free Lagrangian” of the axion field, that in the
case of High Energy Physics keeps care of the effective interaction with
gluonic fields that gives rise to its mass. Here we neglect the interaction
with the gravitational field. For our purposes it is [12]:

$\mathcal{L}_{a}=\frac{1}{2}\partial_{\mu}a\,\partial^{\mu}a-V_{QCD}(a,T).$
(4)

The form of the potential $V_{QCD}(a,T)$ is a nontrivial point for finite
temperatures $T$ [12]. Using the dilute gas approximation one may adopt the
following form when $T=0$ and $T\gg\Lambda_{QCD}$:

$V_{QCD}(a)=f^{2}_{a}m^{2}_{a}(T)\left[1-\cos{\left(\frac{a}{f_{a}}\right)}\right].$
(5)

It is of interest for the present work that the field theory of an axion field
interacting with the gluonic fields by Eq. (5) has a $Z_{N_{DW}}$ symmetry
subgroup if $N_{DW}>1$, with associated $N_{DW}$ degenerate vacua, of the
Peccei-Quinn symmetry that leaves the potential $V_{QCD}(a)$ invariant [2,
11], namely:

$\theta(x)\rightarrow\theta(x)+2\pi k,\qquad k=0,1,...,N_{DW}-1.$ (6)

This discrete symmetry is spontaneously broken when the axion field acquires a
vacuum expectation value, which happens near to the temperature
$T\sim\Lambda_{QCD}$, leading to the formation of domain walls [12].
Generally, a domain wall is a two-dimensional topological defect separating
two domains of a field and are indeed associated to the spontaneous breaking
of a discrete symmetry [25, 11].

The $N_{DW}$ number is called ”domain wall number” in the literature. As
treated in several references, e.g. Refs. [11, 12, 21, 26], the creation of
domain walls is very sensitive for the evolution of the axion field during the
earliest stages of the Universe, overall for scenarios where $N_{DW}>1$ and
the spontaneous breaking of the Peccei-Quinn symmetry happens after inflation.
The scenario described by the formerly cited references seems to lead to an
early Universe dominated by axion domain walls, with an energy density much
bigger than the critical one in the current Universe.

In the following (namely Section (VII.2)) we consider the possibility whether
electromagnetic radiation could significantly contribute to the decay of
domain walls and we enforce what is claimed in Ref. [21], i.e. there cannot be
any significant contribution. In the following we describe a laboratory QCD
axion field by assuming a coherent time-oscillating axion background, namely
$\Theta(t)=\Theta_{0}sin(\omega_{a}t)$. The latter is a very good
approximation when the typical spatial dimension of the system $L$ is much
smaller than the de Broglie wavelength of the axion field (typically of the
order of $10^{-3}\;m_{a}^{-1}$ if we take $v\sim 10^{-3}$ and so
$\omega_{a}\sim m_{a}$ cf., Millar et al. [15]).

###  II.3 Axion-Modified Maxwell equations

It is straightforward [15] to see that the interaction term in Eq. (1) can be
written as an interaction Lagrangian of the form $-J^{\nu}A_{\nu}$ by taking

$J^{\nu}=g_{a\gamma\gamma}\tilde{F}^{\mu\nu}\partial_{\mu}a=\frac{1}{2}g_{a\gamma\gamma}\partial_{\mu}a\;\epsilon^{\mu\nu\rho\sigma}\partial_{\rho}A_{\sigma}.$
(7)

The 4-current generated by axions is thus $J_{a}=(\rho_{a},\vec{J}_{a})$,
where

$\begin{split}\rho_{a}&=g_{a\gamma\gamma}\vec{B}\cdot\nabla a,\\\
\vec{J}_{a}&=g_{a\gamma\gamma}\nabla
a\wedge\vec{E}-g_{a\gamma\gamma}\dot{a}\vec{B}.\end{split}$ (8)

The continuity equation $\dot{\rho}_{a}+\nabla\cdot\vec{J}_{a}=0$ expresses
the same as $\partial_{\mu}j_{a}^{\mu}=0$. Another straightforward way to see
this is by observing from the structure of Eq. (1) that it is possible to
redefine the tensor $H$ as:

$H^{a}_{\mu\nu}=H_{\mu\nu}-g_{a\gamma\gamma}a(x)\tilde{F}_{\mu\nu},$ (9)

in order to rewrite Eq. (1) with the interaction term
$\frac{1}{4}g_{a\gamma\gamma}a\tilde{F}_{\mu\nu}F^{\mu\nu}$ included inside
the ’free term’ $-\frac{1}{4}F^{\mu\nu}H_{\mu\nu}$. We then find that we can
use different constitutive relations for $\vec{D}$ and $\vec{H}$ in Axion
Electrodynamics:

$\displaystyle\vec{D}_{a}=\varepsilon\vec{E}+\vec{P}_{a},$ (10a)
$\displaystyle\vec{B}=\mu\vec{H}_{a}+\vec{M}_{a},$ (10b)

where we have defined the axion-induced polarization vector
$\vec{P}_{a}=-g_{a\gamma\gamma}a\vec{B}$ and magnetization
$\vec{M}_{a}=-\mu\,g_{a\gamma\gamma}a\vec{E}$, which can be interpreted as
polarization and magnetization vectors that are induced by the interaction
with the axion and lead to effective ”bound charge and current densities”
$\rho_{a}=-\nabla\cdot\vec{P}_{a}$ and
$\vec{J}_{a}=\nabla\wedge\vec{M}_{a}+\dot{P}_{a}$. This means that the
interaction between the axion and the electromagnetic field generates an
effective polarization and magnetization, analogously to gravitational
backgrounds.

The Euler-Lagrange equations associated with the Lagrangian density (1) are

$\displaystyle\ddot{a}-\nabla^{2}a+V^{\prime}(a)+g_{a\gamma\gamma}\vec{E}\cdot\vec{B}$
$\displaystyle=0,$ (11a)
$\displaystyle\nabla\cdot\vec{E}=\rho+g_{a\gamma\gamma}\vec{B}\cdot\nabla a,$
(11b)
$\displaystyle\nabla\wedge\vec{B}=\vec{J}+\dot{\vec{E}}-g_{a\gamma\gamma}\dot{a}\vec{B}+g_{a\gamma\gamma}\nabla
a\wedge\vec{E},$ (11c) $\displaystyle\nabla\cdot\vec{B}=0,$ (11d)
$\displaystyle\nabla\wedge\vec{E}=-\dot{\vec{B}}.$ (11e)

Equation (11a) is the Sine-Gordon equation for the axion field interacting
with the electromagnetic field, while equations (11b) and (11c) are the
modified Maxwell equations with sources and the equations (11d) and (11e) are
the ”constraint Maxwell equations” that are unchanged because they follow from
the Bianchi identities $\partial_{\mu}\tilde{F}^{\mu\nu}=0$.

The fields in (10a) satisfy the following forms of Maxwell equations:

$\displaystyle\nabla\cdot\vec{D}_{a}$ $\displaystyle=\rho,$ (12a)
$\displaystyle\nabla\wedge\vec{H}_{a}$
$\displaystyle=\vec{J}+\dot{\vec{D}}_{a},$ (12b)
$\displaystyle\nabla\cdot\vec{B}$ $\displaystyle=0,$ (12c)
$\displaystyle\nabla\wedge\vec{E}$ $\displaystyle=-\dot{\vec{B}}.$ (12d)

###  II.4 Axion Electrodynamics in Condensed Matter

As explained in Refs.[27, 28], the electromagnetic properties of a Weyl
semimetal can be described by a Chern-Simons theory with Lagrangian (1), where
the effective axion field is

$\Theta(x)=b_{0}t-\vec{b}\cdot\vec{r}=b_{\mu}x^{\mu}.$ (13)

The quantities $b_{\mu}$ are of relevant physical meaning for the material,
since they are related to the energy shift $b_{0}$ and the momentum shift
$\vec{b}$ of the specific Weyl point of the material and whose Hamiltonian is

$h_{W}(\vec{k})=b_{0}+v\vec{\sigma}\cdot(\vec{k}-\vec{b}).$ (14)

A Weyl semimetal constitutes an interesting phase of topological quantum
matter with fascinating physical properties such as protected surface states
and the electromagnetic response.

In the following, in order to consider an order of magnitude for $b_{\mu}$ in
the case of time-dependent axion field and space-dependent one, we consider
$b_{0}\sim 10^{-1}\,$\text{\,}\mathrm{e}\mathrm{V}$$ and $|\vec{b}|\sim
10\,$\text{\,}\mathrm{e}\mathrm{V}$$. In the first case, the time derivative
of the axion field is constant, so it is easier to deal with since e.g. the
Green’s function is dependent on $t-t^{\prime}$, and can be regarded as an
approximation of the second case when we consider times
$0<t\ll\omega_{a}^{-1}$. The second case is a coherently oscillating axion
field, which can appropriately approximate a current laboratory axion field
when the typical spatial dimension $L$ of the experimental system is much
smaller than the de Broglie wavelength of the axion field (see Refs. [15,
14]).

###  II.5 Electromagnetic 4-potential in Axion Electrodynamics

As in usual electrodynamics, the differential equations for the 4-potential
$A_{\mu}(x)=(\phi,\vec{A})$ depend on the choice of gauge. We can in general
write the electric and magnetic fields as

$\vec{E}=-\nabla\phi-\frac{\partial\vec{A}}{\partial
t}\qquad\vec{B}=\nabla\wedge\vec{A}.$ (15)

Consequently, we have in the vacuum:

$\displaystyle\nabla^{2}\Phi+\partial_{t}\left(\nabla\cdot\vec{A}\right)$
$\displaystyle=-\rho_{a},$ (16)
$\displaystyle(\partial_{t}^{2}-\nabla^{2})\vec{A}+\nabla\left(\nabla\cdot\vec{A}+\partial_{t}\Phi\right)$
$\displaystyle=\vec{J}_{a}.$ (17)

In the following we will use mainly the Coulomb gauge (radiation gauge if
$\rho_{a}=0$) $\nabla\cdot\vec{A}=0$ and the temporal gauge $A_{0}=\Phi=0.$ In
the case of temporal gauge we have:

$\Box\vec{A}+\nabla(\nabla\cdot\vec{A})=-\dot{\Theta}\;\nabla\wedge\vec{A}-\nabla\Theta\wedge\frac{\partial\vec{A}}{\partial
t}.$ (18)

In the case of the radiation gauge with a time-dependent axion field we get
the equations

$(\varepsilon\mu\partial_{t}^{2}-\nabla^{2})\vec{A}=-\mu\dot{\Theta}\;\nabla\wedge\vec{A}.$
(19)

##  III Energy-momentum balance with fixed and dynamical axion background

It is then easy to notice from the Lagrangian (1) that the total energy-
momentum balance reads (with $J^{\mu}=0$):

$\partial_{\mu}T^{\mu\nu}=0,$ (20)

where we have that $T^{\mu\nu}$ is simply the sum of the free Minkowski
stress-energy tensor
$T^{M\mu\nu}=-\frac{1}{4}\eta^{\mu\nu}F_{\alpha\beta}H^{\alpha\beta}+F^{\mu}_{\rho}H^{\nu\rho}$
and the free axion one
$T^{\mu\nu}_{a}=\partial^{\mu}a\,\partial^{\nu}a-\eta^{\mu\nu}(\frac{1}{2}\partial_{\rho}a\,\partial^{\rho}a-\frac{1}{2}m_{a}^{2}a^{2})$.
This simple result is a consequence of the fact that in this case the system
formed by the electromagnetic field and the axion field is treated as a closed
system (this is not rigorously exact since the potential $V(a)$ comes from the
effective interaction with the gluonic fields). However, if we take the two
fields alone such a relation as Eq.(20) is not valid since they are open
systems interacting between each other. If we consider the electromagnetic
field alone, one gets according to Ref. [16]:

$\partial_{\nu}T^{M\nu}_{\mu}=-f^{M}_{\mu},$ (21)

where $f^{M}_{\mu}=(f^{M}_{0},\vec{f}^{M})$ whose spatial components are the
components of Abraham’s force density
$\vec{f}^{A}=(\varepsilon\mu-1)\frac{\partial}{\partial
t}\Big{(}\vec{E}\wedge\vec{H}\Big{)}-(\vec{E}\cdot\vec{B})\nabla\Theta$ and
$f^{M}_{0}=-g_{a\gamma\gamma}\dot{a}\vec{E}\cdot\vec{B}$.

##  IV Green’s functions of Axion Electrodynamics

We first develop the Green’s function approach in Axion Electrodynamics, using
the temporal gauge. We start from Eq. (18) with $\varepsilon=\mu=1$ and define
the kernel $\vec{G}_{ij}(x,x^{\prime})$ such that:

$A_{i}(x)=\int d^{4}x^{\prime}\;G_{ij}(x,x^{\prime})\,J^{j}(x^{\prime}).$ (22)

Due to causality, the variable $t^{\prime}$ is only integrated over
$t^{\prime}\leq t$. From Eq. (18) we find

$\begin{split}\Box
G_{ij}(x,x^{\prime})+\partial_{i}(\partial_{k}G_{kj}(x,x^{\prime}))-\dot{\Theta}(x)\epsilon_{ikl}\partial_{k}G_{lj}(x,x^{\prime})-\epsilon_{ikl}\partial_{k}\Theta(x)\;\partial_{t}G_{lj}(x,x^{\prime})=\\\
=\delta_{ij}\delta^{(4)}(x^{\mu}-x^{\prime\mu}).\end{split}$ (23)

The $G$ kernel (22) is equal to $iD_{ij}(x^{\mu},y^{\mu})$ where $D$ is the
retarded Green’s function of the vector potential $\vec{A}$, defined as

$iD^{R}_{ij}(x^{\mu},y^{\mu})=\begin{cases}\langle{A_{i}(x)A_{j}(y)-A_{j}(x)A_{i}(y)}\rangle\;\;\;\textit{if}\;\;\;\;x^{0}-y^{0}>0,\\\
0\;\;\;\;\rm{otherwise},\end{cases}$ (24)

can be derived by the Schwinger-Dyson equations [6] for Axion Electrodynamics
which is indeed Eq. (23) with the temporal gauge. As shown in Refs. [29, 30]
and [31], the calculation of the retarded Green’s function (24) is useful in
order to get the significant two-point physical averages between electric and
magnetic fields. One can easily find that with this choice of gauge (valid
also in radiation gauge):

$\displaystyle\langle{E_{i}(x)E_{j}(x^{\prime})}\rangle$
$\displaystyle=\partial_{t}\partial_{t^{\prime}}\langle{A_{i}(x)A_{j}(x^{\prime})}\rangle,$
(25a) $\displaystyle\langle{B_{i}(x)B_{j}(x^{\prime})}\rangle$
$\displaystyle=\operatorname{curl}_{il}\operatorname{curl}^{\prime}_{jm}\langle{A_{l}(x)A_{m}(x^{\prime})}\rangle,$
(25b) $\displaystyle\langle{B_{i}(x)E_{j}(x^{\prime})}\rangle$
$\displaystyle=-\operatorname{curl}_{il}\partial_{t^{\prime}}\langle{A_{l}(x)A_{j}(x^{\prime})}\rangle,$
(25c)

and similarly for other two-point functions. We also use in the following
Fourier transforms (FT) of the two-point functions, e.g.

$\langle{A_{\alpha}(x)A_{\beta}(y)}\rangle_{\omega}=\int_{-\infty}^{+\infty}dt\,e^{i\omega
t}\langle{A_{\alpha}(t,\vec{x})A_{\beta}(0,\vec{y})}\rangle,$ (26)

and similarly for FTs, from which the relations (25) can be used to find the
corrispondent FTs for two-point functions involving electric and magnetic
fields.
For example, if the Green’s function is a function of the variable
$t-t^{\prime},x-x^{\prime}\,\,\text{and}\,\,y-y^{\prime}$ we can define the
reduced Green’s function $g_{ij}(z,z^{\prime},\omega,\tilde{k})$:

$G_{ij}(x,x^{\prime})=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\int_{-\infty}^{\infty}\frac{d^{2}\tilde{k}}{(2\pi)^{2}}e^{-i\omega(t-t^{\prime})}e^{i\tilde{k}\cdot(\vec{r}-\vec{r}^{\prime})}g_{ij}(z,z^{\prime},\omega,\tilde{k}),$
(27)

where we use the notation $\tilde{k}=(k_{x},k_{y})$. We also adopt the useful
quantity $\kappa^{2}=\omega^{2}-k_{x}^{2}-k_{y}^{2}.$ In the following we are
more focused on the Green’s function in the vacuum for calculating Casimir
forces for our systems of interest. From this the following stress tensor
component, important in a Casimir context, can be found straightforwardly
($\varepsilon=\mu=1$):

$\begin{split}\langle{T^{M}_{zz}}\rangle_{\omega,\tilde{k}}=\frac{1}{2i}\left[-\kappa^{2}g_{zz}+(\omega^{2}-k_{y}^{2})g_{xx}+(\omega^{2}-k_{x}^{2})g_{yy}+ik_{y}\left(\partial_{z}g_{yz}-\partial_{z}^{\prime}g_{zy}\right)\right.\\\
\left.+ik_{x}\left(\partial_{z}g_{xz}-\partial_{z}^{\prime}g_{zx}\right)+k_{x}k_{y}(g_{xy}+g_{yx})+\partial_{z}\partial_{z}^{\prime}(g_{xx}+g_{yy})\right].\end{split}$
(28)

For our purpose of evaluating zero-point energies we are also using the
retarded Green’s function with $\mu\neq 1$ in order to find the vacuum one by
using the fluctuation-dissipation theorem [32]:

$\langle{A_{\alpha}(x)A_{\beta}(y)}\rangle_{\omega}=\frac{i}{2}\left[D^{R}_{\alpha\beta}(\omega,\vec{x},\vec{y})-D^{R*}_{\beta\alpha}(\omega,\vec{y},\vec{x})\right]$
(29)

and then make the limit $\mu\rightarrow 1+i0$. Notice that in the case of
temporal gauge, the Green’s function is a $3\times 3$ tensor and not a
$4\times 4$ one. If the axion-generated charge density $\rho_{a}=0$, that is
the case when the axion field is only a function of time $t$, it is convenient
to adopt the radiation gauge and it is what we adopt in the following section.
We adopt the temporal gauge for space-dependent axion fields from Section
(VII) on.

##  V Purely time-dependent axion field

Here we observe that the only difference between the Green’s function in
temporal gauge and radiation gauge (with an axion field depending only on
time) relies on the resulting equation for Green’s function $G_{kj}(x,y)$ that
becomes from Eq. (19), if we assume a homogenous medium with electrical
permittivity $\epsilon=1$ and magnetic permeability $\mu$:

$\left[\delta_{ik}\Box(\mu)-\mu\dot{\Theta}(t)\epsilon_{lik}\nabla_{l}\right]G_{kj}(x,y)=\mu\delta_{ij}\delta^{(4)}(x^{\mu}-y^{\mu}),$
(30)

where the operator $\Box(\mu)=\mu\partial_{t}^{2}-\nabla^{2}$. Since the time
derivative of the axion field can be time-dependent, the solution of Eq.(30)
is generally dependent on $\vec{r}-\vec{r^{\prime}}$ and $t-t^{\prime}$, but
explicitly also on $t$. It is not trivial to find solutions with a generic
time dependence of the axion field. In the following we work out two specific
cases of interest we mentioned in Section II, i.e. $\Theta=\Theta_{0}t$ and
$\Theta=\Theta_{0}\sin(\omega_{a}t)$.

###  V.1 Constant axion time derivative $\dot{\Theta}$

If we consider $\dot{\Theta}$ as constant in time and uniform in space we can
take, similarly to [20], the FT of the Green’s function and get from Eq. (30)
the following equation:

$\left[\delta_{ik}(-\mu\omega^{2}+k^{2})+i\mu\dot{\Theta}\epsilon_{lik}k_{l}\right]\tilde{G}_{kj}(\vec{k},\omega)=\mu\delta_{ij}.$
(31)

The calculations in this case are of conceptual interest and can be applied in
a similar manner in a physically more interesting case, such as the
oscillating axion field. We will come back to this topic in the next section.

In order to get the Green’s function with vacuum surroundings we can put a
magnetic permeability $\mu=1+i\epsilon(\omega)$, where $\epsilon(\omega)$ is
an infinitesimal of the first order, and we evaluate the limit for
$\epsilon\rightarrow 0$ in the proper order, as analogously done in Ref. [32]
with an electrical permeability $\varepsilon(\omega)=1+i\epsilon(\omega)$ to
get the quantum mechanical result for the Green’s function. A magnetic
permeability of the form $\mu=\mu_{r}+i\mu^{\prime}$ can be found, e.g. for
real media where magnetic viscosity is present (see e.g. Ref. [33]):

$\mu=\Re{(\mu)}+i\epsilon\frac{\omega}{\omega_{0}},$ (32)

which is analogous to the case of electric permittivity, where we have the
imaginary part $\Im{(\mu)}>0$ when $\omega>0$. This also justifies the
presence of $\operatorname{sgn}(\omega)$.

From Eq. (31) we get, similarly to Ref. [20], the retarded Green’s function FT
as the inverse operator of the expression in square brackets on the left side
in the mentioned limit of
$\mu=1+i\epsilon\operatorname{sgn}(\omega)\rightarrow
1+i0\operatorname{sgn}(\omega)$:

$\tilde{G}_{jk}(\vec{k},\omega)=(\omega^{2}-|\vec{k}|^{2})\tilde{A}(\omega,\vec{k},\vec{\beta})\delta_{jk}+i\tilde{A}(\omega,\vec{k},\vec{\beta})\epsilon_{jkl}\beta_{l}+\frac{1}{\omega^{2}-|\vec{k}|^{2}+i0\operatorname{sgn}(\omega)}\tilde{A}(\omega,\vec{k},\vec{\beta})\beta_{j}\beta_{k},$
(33)

where $\vec{\beta}$ is a vector with components $\beta_{j}=\dot{\Theta}k_{j}$,
whose squared module is $\beta^{2}=\dot{\Theta}^{2}|\vec{k}|^{2}$ and

$\tilde{A}(\omega,\vec{k},\vec{\beta})=\frac{1}{(\omega^{2}-|\vec{k}|^{2})^{2}-\beta^{2}+i0\operatorname{sgn}(\omega)}.$
(34)

The ordinary electrodynamics limits can be obtained by setting
$\dot{\Theta}=0$. Those Fourier transforms we got can be helpful for two
applications we develop in the following, namely for evaluating the spectral
energy density $\rho_{em}(\omega)$ of the electromagnetic field at temperature
$T$, the Casimir force between two parallel conducting plates in such a fixed
axionic background and a discussion on dispersion relations.

####  V.1.1 Spectral energy density $\rho_{em}(\omega)$ of electromagnetic
field

In order to get the zero-point spectral energy density we calculate the FT
$\tilde{G}(\vec{r},\omega)$ of the Green’s function (33). When calculating
energy density we need the correlation functions in the limit
$\vec{r}\rightarrow 0$. The usual spectral energy density of the
electromagnetic field is

$\rho(\omega,\vec{r})d\omega=\frac{1}{2}\left[2\langle{\vec{E}^{2}(\vec{r})}\rangle_{\omega,T}+2\langle{\vec{B}^{2}(\vec{r})}\rangle_{\omega,T}\right]\frac{d\omega}{2\pi},$
(35)

where the factors 2 inside the brackets are inserted because we follow the
definition of spectral densities given in Ref. [32]: average spectral values
are defined as integrals in $\omega$ from $-\infty$ to $+\infty$. We also
exploit the relation:

$\frac{1}{x\pm i0}=\mathcal{P}{\frac{1}{x}}\mp i\pi\delta(x).$ (36)

We expand the spectral energy density in a perturbative series:

$\rho(\omega,\vec{r})=\rho^{(0)}(\omega,\vec{r})+\rho^{(1)}(\omega,\vec{r})+\rho^{(2)}(\omega,\vec{r})+...$
(37)

according to the perturbative factor, that we can take here to be
$\dot{\Theta}$. We can develop a perturbative expansion for the Green’s
function:

$\tilde{G}_{jk}(\vec{k},\omega)\sim\tilde{G}^{(0)}_{jk}(\vec{k},\omega)+\tilde{G}^{(1)}_{jk}(\vec{k},\omega)+\tilde{G}^{(2)}_{jk}(\vec{k},\omega),$
(38)

where

$\displaystyle\tilde{G}^{(0)}_{jk}(\vec{k},\omega)$
$\displaystyle=\frac{1}{\omega^{2}-k^{2}+i0\operatorname{sgn}{\omega}},$ (39a)
$\displaystyle\tilde{G}^{(1)}_{jk}(\vec{k},\omega)$
$\displaystyle=\frac{i}{(\omega^{2}-k^{2})^{2}+i0\operatorname{sgn}{\omega}}\epsilon_{jkl}\beta_{l}=\frac{i}{(\omega^{2}-k^{2})^{2}+i0\operatorname{sgn}{\omega}}\epsilon_{jkl}(\dot{\Theta}k_{l}),$
(39b) $\displaystyle\tilde{G}^{(2)}_{jk}(\vec{k},\omega)$
$\displaystyle=\frac{\beta^{2}}{(\omega^{2}-k^{2})^{3}+i0\operatorname{sgn}{\omega}}\left[\delta_{jk}+\frac{\epsilon_{jkl}\beta_{l}}{\omega^{2}-k^{2}}\right]+\frac{1}{(\omega^{2}-k^{2})^{3}+i0\operatorname{sgn}{\omega}}\dot{\Theta}^{2}k_{j}k_{k}.$
(39c)

As shown in Ref. [32] it is easy to obtain (by integrating over the domain of
$\vec{k}$) that the spectral density associated with the zero order term
$\rho^{(0)}(\omega)$ is

$\rho^{(0)}(\omega)d\omega=d\omega\;\omega^{2}\left[\frac{1}{2}\omega+\frac{\omega}{e^{\frac{\omega}{T}}-1}\right].$
(40)

This is the familiar result from QED: the zero-point energy of electromagnetic
field at a given temperature $T$ is the sum of the zero-point energy at the
zero Kelvin temperature and the black-body radiation energy at the given
temperature $T$. The first order term $\rho^{(1)}(\omega)$ is trivially zero,
while we have a contribution of order
$\mathcal{O}(g^{2}_{\alpha\gamma\gamma})$:

$\rho^{(2)}(\omega)\;d\omega=d\omega\;\omega^{2}\left[\frac{1}{2}\frac{\dot{\Theta}^{2}}{8\omega}+\frac{1}{e^{\frac{\omega}{T}}-1}\frac{\dot{\Theta}^{2}}{8\omega}\right].$
(41)

We notice in Eq. (41) the frequency dependence of the spectral energy density,
from which we can find the spectral emissivity using Planck’s law. The first
term on the right hand side is proportional to the frequency. As
$\dot{\Theta}$ can be regarded small for physical purposes, this frequency
dependence suggests that this can be physically interpreted as a blue noise.

The second term is similar to the first one as regards the frequency
dependence, but it has a Bose-Einstein temperature-dependent weight. For
$T\gg\omega$ it is term going as $\simeq T$, so it can be regarded as a flat
noise, and surely is bigger than the zero-temperature term. For $T\ll\omega$
the Bose-Einstein weight would be approximately $e^{-\frac{\omega}{T}}$, so it
gives an exponential noise.

####  V.1.2 Casimir force between parallel plates with time-increasing axion
field

Here we evaluate the axion modifications to the usual Casimir force between
two parallel perfectly conducting plates (see Refs. [30, 6] for a treatment
without an axion background). Our result is comparable to the one obtained by
Ref. [34], while improving their result.

We can evaluate the zero-point energy $u_{em}(L,\dot{\Theta})$ per unit
transverse area by evaluating the following expression:

$u_{em}(L,\dot{\Theta})=\frac{1}{2}\sum_{\pm}\sum_{n}\int\frac{d^{2}k}{(2\pi)^{2}}\sqrt{|\vec{k}|^{2}+\frac{n^{2}\pi^{2}}{L^{2}}\pm\dot{\Theta}\sqrt{|\vec{k}|^{2}+\frac{n^{2}\pi^{2}}{L^{2}}}}.$
(42)

This could be expected by a Casimir approach, namely by physical qualitative
means given by the dispersion relations . The expression (42) can be
demonstrated overall by a Green’s function approach. We assume the $z$-axis to
be the direction normal to the plates. We need to solve the Green’s function
equations with the proper boundary conditions. Since we assume to have two
perfect conducting plates we have boundary conditions for electric field and
magnetic field at $z=0,L$:

$B_{z}=0\qquad\vec{E}_{//}=0,$ (43)

where $\vec{E}_{//}$ is the parallel component of the electric field, from
which we get the correspondent ones for the Green’s function:

$\displaystyle g_{ij}(x-x^{\prime},y-y^{\prime},z,z^{\prime})|_{z=0,z=L}$
$\displaystyle=0\quad i\neq z,j\neq z,$ (44a)
$\displaystyle\partial_{z}g_{zz}(x-x^{\prime},y-y^{\prime},z,z^{\prime})|_{z=0,z=L}$
$\displaystyle=0.$ (44b)

These conditions, along with translational invariance along x and y
directions, lets us write the Green’s function in the following Fourier form
for $0<z<L$ and $0<z^{\prime}<L$:

$\displaystyle G_{ij}(x,x^{\prime})$
$\displaystyle=\int_{-\infty}^{+\infty}\frac{d\omega}{2\pi}\int_{-\infty}^{\infty}\frac{d^{2}\tilde{k}}{(2\pi)^{2}}e^{i\tilde{k}\cdot(\tilde{r}-\tilde{r}^{\prime})}\frac{2}{L}\sum_{n=0}^{\infty}\sin{\left(\frac{n\pi}{L}z\right)}\sin{\left(\frac{n\pi}{L}z^{\prime}\right)}\tilde{g}_{ij}(\tilde{k},n)\quad
i\neq z,j\neq z,$ (45a) $\displaystyle G_{zz}(x,x^{\prime})$
$\displaystyle=\int_{-\infty}^{+\infty}\frac{d\omega}{2\pi}\int_{-\infty}^{\infty}\frac{d^{2}\tilde{k}}{(2\pi)^{2}}e^{i\tilde{k}\cdot(\tilde{r}-\tilde{r}^{\prime})}\frac{2}{L}\sum_{n=0}^{\infty}\iota(n)\cos{\left(\frac{n\pi}{L}z\right)}\cos{\left(\frac{n\pi}{L}z^{\prime}\right)}\tilde{g}_{zz}(\tilde{k},n),$
(45b)

where $\tilde{k}=(k_{x},k_{y})$, $\tilde{r}=(x,y)$ and $\iota(n)=1$ for $n>0$,
while $\iota(n=0)=1/2$.
The reduced function $\tilde{G}(\tilde{k},n)$ consequently satisfies the
equations (31), with solution Eq. (33) ,where $k_{z}=\frac{n\pi}{L}$. We show
now two ways of demonstrating Eq. (42), then we calculate the Casimir force.

##### Energy density per transversal area from $\Im{(\mu)}\rightarrow 0$

Here we adopt the form (33) for the Green’s function. If we exploit the
relation (36) and use the Green’s function (33) we get

$\frac{1}{2}\langle{\vec{E}^{2}+\vec{B}^{2}}\rangle_{\omega,k_{x},k_{y},n}=\pi\frac{\omega}{2}\sum_{\pm}\left[\delta\left(\omega-\sqrt{|k|^{2}\pm\dot{\Theta}|k|}\right)-\delta\left(\omega+\sqrt{|k|^{2}\pm\dot{\Theta}|k|}\right)\right],$
(46)

where $|k|^{2}=k_{x}^{2}+k_{y}^{2}+\left(\frac{n\pi}{L}\right)^{2}$. We then
exploit the properties of the Dirac delta distribution by integrating in the
$\frac{\omega}{2\pi}$ domain from $-\infty$ to $+\infty$ and we then get

$\frac{1}{2}\langle{\vec{E}^{2}+\vec{B}^{2}}\rangle_{k_{x},k_{y},n}=\frac{1}{2}\sum_{\pm}\sqrt{k_{x}^{2}+k_{y}^{2}+\frac{n^{2}\pi^{2}}{L^{2}}\pm\dot{\Theta}\sqrt{k_{x}^{2}+k_{y}^{2}+\frac{n^{2}\pi^{2}}{L^{2}}}}.$
(47)

We also mention that $\langle{\vec{A}\cdot\vec{B}}\rangle_{k_{x},k_{y},n}$ and
$\langle{\vec{E}\cdot\vec{B}}\rangle_{k_{x},k_{y},n}$ are trivially equal to
zero. Consequently, we have
$\langle{T^{00}_{em}}\rangle_{k_{x},k_{y},n}=\frac{1}{2}\langle{\vec{E}^{2}+\vec{B}^{2}}\rangle_{k_{x},k_{y},n}$,
from which it is easy to get the energy per transverse area unit to be equal
to the expression (42).

##### Casimir force by the reduced Green’s function in the case $\mu=1$

In general, the Casimir force is the physical force between dielectric or
metallic surfaces (and in different configurations), caused by quantum
fluctuations of the electromagnetic field. From a general perspective, the
Casimir effect can be treated by different methods.

It can be derived from statistical equilibrium thermodynamics, as an
application of the fluctuation-dissipation theorem, which relates the two-
point functions of the electromagnetic field to the imaginary part of retarded
Green function, as discussed in Landau-Lifshitz book on statistical physics
[32] and mentioned before, as well as in several articles, for instance
Ref.[35]. Alternatively , the effect can be seen upon as a manifestation of
quantum field theory; cf. Milton’s book [30].

In both cases, the essential ingredient required in practical calculations is
the difference between the normal Maxwell stress components $T_{zz}$ on the
two sides of a boundary. We explicitly calculate the expression (28) with the
Green’s function solving the equation (31) in the form (45) when exactly
$\mu=1$. It is easy to get:

$\langle{T_{zz}}\rangle_{\omega,k_{x},k_{y}}|_{z=z^{\prime}=0,L}=-2i\sum_{n=1}^{+\infty}\frac{n^{2}\pi^{2}}{L^{3}}\frac{\kappa^{2}-\frac{n^{2}\pi^{2}}{L^{2}}}{(\kappa^{2}-\frac{n^{2}\pi^{2}}{L^{2}})^{2}-\dot{\Theta}^{2}(k_{x}^{2}+k_{y}^{2}+\frac{n^{2}\pi^{2}}{L^{2}})},$
(48)

We remind that $\kappa^{2}=\omega^{2}-k_{x}^{2}-k_{y}^{2}$. This series can be
treated by noticing that the $n$-term of the series can be written as:

$\frac{\frac{n^{2}\pi^{2}}{L^{2}}(\kappa^{2}-\frac{n^{2}\pi^{2}}{L^{2}})}{(\kappa^{2}-\frac{n^{2}\pi^{2}}{L^{2}})^{2}-\dot{\Theta}^{2}(k_{x}^{2}+k_{y}^{2}+\frac{n^{2}\pi^{2}}{L^{2}})}=\frac{1}{2}\sum_{\pm}\frac{\frac{n^{2}\pi^{2}}{L^{2}}}{\kappa^{2}-\frac{n^{2}\pi^{2}}{L^{2}}\pm\dot{\Theta}\sqrt{k_{x}^{2}+k_{y}^{2}+\frac{n^{2}\pi^{2}}{L^{2}}}},$
(49)

and one can perform a standard complex frequency rotation
$\omega\rightarrow\zeta=i\omega$ , similarly to what done in Ref. [30], and
integrate over $\zeta/{2\pi}$, obtaining:

$\langle{T_{zz}}\rangle_{k_{x},k_{y}}|_{z=z^{\prime}=0,L}=\sum_{\pm}\sum_{n=1}^{+\infty}\frac{n^{2}\pi^{2}}{L^{3}}\frac{1}{\sqrt{k_{x}^{2}+k_{y}^{2}+\frac{n^{2}\pi^{2}}{L^{2}}\pm\dot{\Theta}\sqrt{k_{x}^{2}+k_{y}^{2}+\frac{n^{2}\pi^{2}}{L^{2}}}}}.$
(50)

It is easy to see from Eqs. (48,(49) and (50) that, also up to contact terms,
we are dealing with divergent expressions. This is because the physical
quantity we observe (Casimir force) is the discontinuity of $T_{zz}$, not the
single values on the two sides of the interface(at $z=0^{+}$ as in this case).
We then need to consider $T_{zz}$ at $z=0^{-}$.

In order to do that, one needs to solve Eq. (30) with $z<0$ and $z^{\prime}<0$
(the case with $z>L$ gives the same results from symmetry) with the boundary
conditions (44) at z=0 and $g\sim e^{-ikz}$ for $z\rightarrow-\infty$. This
can be done in a similar fashion to what done in Stakgold’s book [36], by
noticing that, accordingly to the boundary condition at z=0, the sine (or
cosine) Fourier transform satisfies Eq.(31).

We can then solve the equations, as done in Ref. [36], and get the following
$T_{zz}$ expression:

$\langle{T_{zz}}\rangle_{\omega,k_{x},k_{y}}|_{z=z^{\prime}=0^{-}}=-2i\int_{0}^{+\infty}k_{z}^{2}\frac{\kappa^{2}-k_{z}^{2}}{(\kappa^{2}-k_{z}^{2})^{2}-\dot{\Theta}^{2}(k_{x}^{2}+k_{y}^{2}+k_{z}^{2})}dk_{z}.$
(51)

We manipulate this expression in a similar fashion to the one at
$0^{+}$,integrate over $\zeta$ and omit contact terms. We then get the Casimir
force as the discontinuity of $T_{zz}$:

$f(L,\dot{\Theta})=\int_{-\infty}^{\infty}\frac{dk_{x}}{2\pi}\int_{-\infty}^{\infty}\frac{dk_{y}}{2\pi}\frac{1}{L}\sum_{n=1}^{+\infty}\sum_{\pm}\left[\frac{\frac{n^{2}\pi^{2}}{L^{2}}\pm\dot{\Theta}\frac{n^{2}\pi^{2}}{L^{2}\sqrt{k_{x}^{2}+k_{y}^{2}+\frac{n^{2}\pi^{2}}{L^{2}}}}}{\sqrt{|\vec{k}|^{2}+n^{2}\pi^{2}/L^{2}\pm\dot{\Theta}\sqrt{|\vec{k}|^{2}+\frac{n^{2}\pi^{2}}{L^{2}}}}}\right].$
(52)

It is easy to notice that the following formal relation holds:

$f(L,\dot{\Theta})=-\frac{\partial u_{em}(L,\dot{\Theta})}{\partial L},$ (53)

confirming our results.

##### Calculation of the Casimir force

In order to evaluate the Casimir force deriving from the expression (42) it is
convenient to adopt the zeta function regularization method and extending it
to a dimension $d\neq 2$ (see e.g. Milton [30] and Brevik [37]). We employ the
Schwinger proper-time representation for the square root:

$u_{em}=\frac{1}{2}\sum_{\pm}\sum_{n}\int\frac{d^{d}k}{(2\pi)^{d}}\int_{0}^{+\infty}\frac{dt}{t}t^{-\frac{1}{2}}e^{-t\left(k^{2}+n^{2}\pi^{2}/L^{2}\pm\dot{\Theta}\sqrt{|\vec{k}|^{2}+\frac{n^{2}\pi^{2}}{L^{2}}}\right)}\frac{1}{\Gamma(-\frac{1}{2})},$
(54)

It is not easy to evaluate it exactly, so we will calculate it perturbatively
up to second order in $\dot{\Theta}$. It is possible to write the Taylor
series of the exponential in Eq. (54) up to the second order in $\dot{\Theta}$
and sum up in the two polarizations explicitly:

$u_{em}\simeq\sum_{n}\int\frac{d^{d}k}{(2\pi)^{d}}\int_{0}^{+\infty}\frac{dt}{t}t^{-\frac{1}{2}}e^{-t\left(k^{2}+n^{2}\pi^{2}/L^{2}\right)}\left[2+t^{2}\dot{\Theta}^{2}\left(k^{2}+\frac{n^{2}\pi^{2}}{L^{2}}\right)\right]\frac{1}{\Gamma(-\frac{1}{2})}.$
(55)

As it is usually done in this kind of calculations, we carry out the Gaussian
integration over k, use the Euler representation of Gamma function, carry out
the sum over $n$ via the definition of Riemann zeta function and we get:

$u_{em}\simeq u_{0}(L,d)+u_{ax}(L,d,\dot{\Theta}),$ (56)

where

$\displaystyle u_{0}(L,d)$
$\displaystyle=-\frac{1}{2\sqrt{\pi}}\frac{1}{(4\pi)^{d/2}}\left(\frac{\pi}{L}\right)^{d+1}\Gamma\left(-\frac{d+1}{2}\right)\zeta(-d-1),$
(57a) $\displaystyle u_{a}(L,d,\dot{\Theta})$
$\displaystyle=-\frac{\dot{\Theta}^{2}}{2\sqrt{\pi}}\frac{1}{(4\pi)^{d/2}}\left(\frac{\pi}{L}\right)^{d-1}\left[\Gamma\left(\frac{3-d}{2}\right)\zeta(-d-1)+\frac{d}{2}\,\Gamma\left(\frac{1-d}{2}\right)\zeta(-d+1)\right].$
(57b)

The first term $u_{0}$ is the familiar result of Quantum Electrodynamics,
while $u_{a}$ is the additional second-order axion contribution. For $d=2$
they are equal to

$\displaystyle u_{0}(L,d=2)=-\frac{\pi^{2}}{720}\frac{1}{L^{3}}\,\,\,,\qquad
u_{a}(L,\dot{\Theta},d=2)=-\dot{\Theta}^{2}\frac{7}{320}\frac{1}{L}\,\,\,,$
(58a)

from which:

$\displaystyle f_{0}(L)=-\frac{\pi^{2}}{240}\frac{1}{L^{4}}\,\,\,,\qquad
f_{a}(L,\dot{\Theta})=-\frac{7\dot{\Theta}^{2}}{320}\frac{1}{L^{2}}\,\,\,.$
(59a)

It is worth noticing that we get a second-order axion-induced modification of
the Casimir energy per transverse area going as $\sim\frac{1}{L}$, similarly
to what is obtained in Ref. [34], although our result differs by a
multiplicative factor of $0.63$ (our result is lower). There are some points
to notice about the possibilities of physical application of these results.

The approximation $\dot{\Theta}\sim$ constant for the fundamental physics
axion field could be applicable to approximate a coherently oscillating axion
field for very small times $t\ll\omega_{a}^{-1}$, meaning that the theory is a
reasonable approximation for frequencies $\omega\gg\omega_{a}$. It could then
be used when $\frac{1}{L}\gg\omega_{a}$. This means that the treatment of the
Casimir force between two conducting plates with such an axion background is
relevant for experimental purposes because the assumption for the perturbative
approach is $L\ll(\Theta_{0}\,\omega_{a})^{-1}$, which is a stricter
requirement than the long wavelength approximation, since for an axion being
the main component of Dark Matter we need $\Theta_{0}\sim 10^{-19}$ [37, 38]
at cosmological scales.

If we take then the case of High Energy Physics axion at cosmological scales,
$\dot{\Theta}\sim\Theta_{0}\omega_{a}$, and a comparison between the terms in
Eq. (58) tells us that $u_{a}$ is comparable with $u_{0}$ when
$L_{a}\sim(\Theta_{0}\,\omega_{a})^{-1}\sim
10^{14}\left(\frac{m_{a}}{10^{-2}$\text{\,}\mathrm{e}\mathrm{V}$}\right)$
meters, so $u_{a}$ would be very small for usual Casimir experimental devices,
if there are no significant local deviations in the axion field. However, this
axionic term could be significant in the case of topological insulators.
Namely, if we take $b_{0}\sim 10^{-1}\,$\text{\,}\mathrm{e}\mathrm{V}$$, this
corresponds to $\dot{\Theta}\sim 10^{-1}\,$\text{\,}\mathrm{e}\mathrm{V}$$, so
now $L_{a}\sim$2\text{\,}\mu\mathrm{m}$$.

####  V.1.3 Dispersion relations

Here we briefly discuss the dispersion relations for the case of a constant
$\dot{\Theta}$. From the Green’s function it is easy to find the 4 poles to
which they are associated with the dispersion relations:

$\omega_{\pm}=\sqrt{|\vec{k}|^{2}\pm\dot{\Theta}|\vec{k}|}\qquad\omega_{0}=|\vec{k}|.$
(60)

This last observation can be understood in the following way: one firstly
consider the Fourier equations of vector potential for a plane wave
propagating along the direction of the $z$-axis:

$\displaystyle\begin{split}(-\omega^{2}+k_{z}^{2})\mathcal{A}_{x}&=\dot{\Theta}k_{z}\mathcal{A}_{y},\\\
(-\omega^{2}+k_{z}^{2})\mathcal{A}_{y}&=-\dot{\Theta}k_{z}\mathcal{A}_{x},\\\
(-\omega^{2}+k_{z}^{2})\mathcal{A}_{z}&=0.\end{split}$ (61)

This assumption on the propagation direction does not lose generality, since
the treated system in the vacuum is isotropic and it is possible to find an
inertial reference frame, where the electromagnetic wave propagates along the
$z$-axis. It is then easy to demonstrate the dispersion relation
$\omega_{\pm}$ of (60) for real transverse photons. Indeed, if we define the
fields as $\mathcal{A}_{\pm}=\mathcal{A}_{x}\pm i\mathcal{A}_{y}$, one gets
the dispersion relations $\omega_{\pm}$ with $k_{x}=k_{y}=0$ . One then gets
the expression for $\omega_{\pm}$ in Eq. (60) by using the relativistic
invariance of the magnitude of the wave 4-vector. The physical meaning of this
is interesting: a left-circular polarized wave has a different frequency than
a right-circular one with the same $\vec{k}$. This is due to the optical
activity of the vacuum in Axion Electrodynamics as we shall discuss in Section
VIII.

We also notice that a solution for the equation of $\mathcal{A}_{z}$, the one
with the dispersion relation $\omega=|\vec{k}|$ and corresponding to
longitudinal photons, is simply zero, so that the dispersion relation
corresponds to virtual photons.

###  V.2 High frequency approximation for the Green’s function

Before working out the case of an oscillating axion field (by perturbative
methods), we can get some understanding of the behavior by considering the
High frequency approximation in that case. Start from Eq. (23), substitute the
form $\Theta(t)=\Theta_{0}\sin{(\omega_{a}t)}$, and calculate the FT:

$(-\omega^{2}+|\vec{k}|^{2})\tilde{G}_{ij}(\omega,\vec{k})-i\epsilon_{lik}k_{l}\frac{\Theta_{0}\omega_{a}}{2}\left[\tilde{G}_{ij}(\omega+\omega_{a},\vec{k})+\tilde{G}_{ij}(\omega-\omega_{a},\vec{k})\right]=\delta_{ij}.$
(62)

We have here used the modulation property of the Fourier Transform in time.
The High frequency approximation consists on assuming $\omega\gg\omega_{a}$.
In this way, up to first order in $\omega_{a}$, the Green’s function satisfies
the following equation,

$(-\omega^{2}+|\vec{k}|^{2})\tilde{G}_{ij}(\omega,\vec{k})-i\epsilon_{lik}k_{l}\Theta_{0}\omega_{a}\tilde{G}_{kj}(\omega,\vec{k})=\delta_{ij}.$
(63)

It has the same form as the preceding equation and has accordingly the same
solutions, with just the substitution
$\dot{\Theta}\rightarrow\Theta_{0}\omega_{a}$. The interpretation of Eq. (41)
terms are still valid in this case.

###  V.3 Case with $\Theta(t)=\Theta_{0}\,sin(\omega_{a}t)$: Spectral energy
density and production of real photons

Now we treat a case of physical interest as mentioned before. In order to
treat this case properly we need to work out Eq. (30) in more detail.
We assume that the wave vector is directed along a fixed direction, without
losing generality. Actually, the axion field is taken to be dependent on time
but spatially homogeneous, so it is isotropic. We make a FT in space and put
$\vec{k}=|\vec{k}|\hat{e}_{z}$,

$\begin{cases}(\partial_{t}^{2}+|\vec{k}|^{2})\tilde{G}_{xx}(\vec{k},t,t^{\prime})+i\dot{\Theta}k_{z}\tilde{G}_{zx}(\vec{k},t,t^{\prime})=\delta(t-t^{\prime}),\\\
(\partial_{t}^{2}+|\vec{k}|^{2})\tilde{G}_{yy}(\vec{k},t,t^{\prime})+i\dot{\Theta}k_{z}\tilde{G}_{xy}(\vec{k},t,t^{\prime})=\delta(t-t^{\prime}),\\\
(\partial_{t}^{2}+|\vec{k}|^{2})\tilde{G}_{zz}(\vec{k},t,t^{\prime})=\delta(t-t^{\prime}),\\\
(\partial_{t}^{2}+|\vec{k}|^{2})\tilde{G}_{xy}+i\dot{\Theta}k_{z}\tilde{G}_{yy}=0,\\\
(\partial_{t}^{2}+|\vec{k}|^{2})\tilde{G}_{yx}-i\dot{\Theta}k_{z}\tilde{G}_{xx}=0,\\\
(\partial_{t}^{2}+|\vec{k}|^{2})\tilde{G}_{zy}+i\dot{\Theta}k_{z}\tilde{G}_{xy}=0.\end{cases}$
(64)

We will solve this coupled differential equation system by a perturbative
approach. We expand the formalism up to second order in $\Theta_{0}$ in the
following way.
Consider $G_{xy}$ and $G_{yy}$, whose equations are coupled between themselves
only. To order zero we get

$\begin{split}G_{yy}(\vec{k},t,t^{\prime})&=H(t-t^{\prime})\;\times\frac{1}{|\vec{k}|}\,\sin\left[|\vec{k}|(t-t^{\prime})\right],\\\
G_{xy}(\vec{k},t,t^{\prime})&=0,\end{split}$ (65)

according to the initial conditions, where the function $H(t-t^{\prime})$ is
the usual Heaviside function. It is then easy to observe that the first-order
correction $G^{(1)}_{yy}$ is zero because it satisfies the same equation as
$G^{(0)}_{xy}$ with the same boundary conditions. This is not true for
$G^{(1)}_{xy}$ because $G^{(0)}_{yy}$ is not zero:

$(\partial_{t}^{2}+|\vec{k}|^{2})G^{(1)}_{xy}+i\dot{\Theta}k_{z}G^{(0)}_{yy}=0.$
(66)

The solution of this equation has the form of a convolution integral,

$\Theta_{0}G^{(1)}_{xy}=ik_{z}\int_{-\infty}^{\infty}d\omega\;G_{yy}(t,\omega)\mathcal{F}(\dot{\Theta}G^{(0)}_{yy})(\omega),$
(67)

where $\mathcal{F}(\dot{\Theta}G^{(0)}_{yy})(\omega)$ is the FT of
$\dot{\Theta}G^{(0)}_{yy}$ in $\omega$. It now becomes necessary to know the
specific time dependence of the axion field in order to evaluate the Fourier
Transform of the product of the time derivative of the axion field and the
Green’s function,

$\mathcal{F}(\dot{\Theta}G^{(0)}_{yy})(\omega)=\Theta_{0}\omega_{a}\int
dt^{\prime}e^{i\omega
t^{\prime}}\cos{(\omega_{a}t^{\prime})}G^{(0)}_{yy}(t^{\prime})=\frac{\Theta_{0}\omega_{a}}{2}\left[\frac{1}{(\omega-\omega_{a})^{2}-|\vec{k}|^{2}}+\frac{1}{(\omega+\omega_{a})^{2}-|\vec{k}|^{2}}\right].$
(68)

Analogously, we obtain:

$\tilde{G}^{(1)}_{xy}=i|\vec{k}|\frac{\Theta_{0}\omega_{a}}{2}\frac{1}{\omega^{2}-|\vec{k}|^{2}}\left[\frac{1}{(\omega-\omega_{a})^{2}-|\vec{k}|^{2}}+\frac{1}{(\omega+\omega_{a})^{2}-|\vec{k}|^{2}}\right],$
(69)
$\tilde{G}^{(2)}_{yy}=\tilde{G}^{(2)}_{xx}=-|\vec{k}|^{2}\frac{\Theta_{0}^{2}\omega_{a}^{2}}{4}\frac{1}{\omega^{2}-|\vec{k}|^{2}}\left[\frac{1}{(\omega-\omega_{a})^{2}-|\vec{k}|^{2}}+\frac{1}{(\omega+\omega_{a})^{2}-|\vec{k}|^{2}}\right]^{2}.$
(70)

This is a result compatible with what could have been expected from the HF
approximation discussed above. It is now understandable that according to the
results above, a useful ansatz for the expansion is

$\tilde{G}_{jk}(\vec{k},\omega)\sim\tilde{G}^{(0)}_{jk}(\vec{k},\omega)+\tilde{G}^{(1)}_{jk}(\vec{k},\omega)+\tilde{G}^{(2)}_{jk}(\vec{k},\omega),$
(71)

where

$\displaystyle\tilde{G}^{(0)}_{jk}(\vec{k},\omega)$
$\displaystyle=\frac{1}{(\omega^{2}-k^{2})},$ (72a)
$\displaystyle\tilde{G}^{(1)}_{jk}(\vec{k},\omega)$
$\displaystyle=\frac{i}{(\omega^{2}-k^{2})^{2}}\left[\frac{1}{(\omega-\omega_{a})^{2}-|\vec{k}|^{2}}+\frac{1}{(\omega+\omega_{a})^{2}-|\vec{k}|^{2}}\right]\epsilon_{jkl}\beta_{l},$
(72b) $\displaystyle\tilde{G}^{(2)}_{jk}(\vec{k},\omega)$
$\displaystyle=\frac{\beta^{2}}{(\omega^{2}-k^{2})}\left[\frac{1}{(\omega-\omega_{a})^{2}-|\vec{k}|^{2}}+\frac{1}{(\omega+\omega_{a})^{2}-|\vec{k}|^{2}}\right]\left[\delta_{jk}+\frac{\epsilon_{jkl}\beta_{l}}{\omega^{2}-k^{2}}\right]+$
(72c) $\displaystyle+$
$\displaystyle\frac{1}{(\omega^{2}-k^{2})}\left[\frac{1}{(\omega-\omega_{a})^{2}-|\vec{k}|^{2}}+\frac{1}{(\omega+\omega_{a})^{2}-|\vec{k}|^{2}}\right]^{2}\Theta_{0}^{2}\omega_{a}^{2}k_{j}k_{k}.$
(72d)

It is worth mentioning that the results for axion oscillating background are
compatible with models of axion echo, treated recently in Ref. [39]. The terms
needed in that context can be worked out as done before with a constant time
derivative.

Figure 1: Graph of the function f in Eq. (74) with $\Theta_{0}=1$. Notice that
for $\omega>\omega_{a}$ we have the expected behaviour from Eq. (41), whose
infrared cut-off is the peak at $\omega=\omega_{a}$, associated with the
limiting case of decay of one axion to only one detectable photon.
Furthermore, we have a peak at $\omega=\frac{\omega_{a}}{2}$, a result
compatible with the axion echo phenomenon [39], so the physical process behind
is the decay of an axion of energy $\omega_{a}$ to two photons with energy
$\frac{\omega_{a}}{2}$.

It is also possible in this case to evaluate the second-order contribution to
energy density:

$\rho^{(2)}(\omega)\;d\omega=d\omega\;\omega^{2}\left[\frac{1}{2}f(\omega,\Theta_{0},\omega_{a})+\frac{1}{e^{\omega/T}-1}f(\omega,\Theta_{0},\omega_{a})\right],$
(73)

where

$f(\omega,\Theta_{0},\omega_{a})=\frac{\Theta_{0}^{2}\omega_{a}^{2}\omega(8\omega^{4}+\omega_{a}^{4})}{2(\omega+\omega_{a})|\omega-\omega_{a}|(\omega_{a}^{2}-4\omega^{2})^{2}}.$
(74)

A numerical plot of the function $f$ is displayed in Figure 1: from it one
could notice how the validity of an expression of Casimir force between two
conducting plates analogous to the case of time increasing axion field is true
for $L^{-1}\gg\omega_{a}$, as mentioned before. For bigger L one appreciates
the deviations from that case and consider that the behaviour of the
modification to energy density going as $1/\omega$ needs to be corrected with
$1/|\omega-\omega_{a}|$ and has a lower cut-off at $\omega=\omega_{a}$.
Furthermore, one needs to take care of the additional contribution of ”axion
echo” [39] at $\omega=\frac{\omega_{a}}{2}$.

This last contribution has a straightforward physical interpretation :it is
associated to the production of virtual photons with frequency
$\omega=\frac{\omega_{a}}{2}$ from the decay of an axion, while the peak at
$\omega=\omega_{a}$ is associated to a decay of one axion to only one
detectable photon.

###  V.4 A possible way to boost axionic Casimir effect in a a current
Universe experimental set-up

In this subsection we will consider a dielectric system containing two
interfaces separating media of refractive indices $n_{1}$ and $n_{2}$ inside
an uniform magnetic field $B_{e}$ and we assume that the media are elastic so
that medium 2 can be ”turned back” and glued to the left side of medium 1.
This is a ring-formed system that we have been treated in Ref. [40], but we
extend the theory treated there.

Figure 2: Geometry and notation of the closed string of interest. Figure from
Ref. [40]

.

Figure 2 shows the configuration and here we treat the most fundamental
results from [40] in order to develop here our new results relative to axion-
induced radiation pressure.
$\sigma$ is the length coordinate along the string, such that the two
junctions are at $\sigma=0$ and $\sigma=L_{I}$. The total length of the string
is $L=L_{I}+L_{II}$ and the junctions $\sigma=0$ and $\sigma=L$ are
overlapping. The string is lying in the $xy$ plane, and a strong uniform
magnetic field ${\bf B}_{e}$ is applied in the $z$ direction. In Ref. [40] we
evaluated the stationary oscillations of the electromagnetic oscillations in
the string. If $E_{I}(\sigma,t)$ and $E_{II}(\sigma,t)$ are the electric
fields in the two regions, we have in complex representation

$\displaystyle E_{I}(\sigma,t)=\xi_{I}e^{in_{1}\omega\sigma-i\omega
t}+\eta_{I}e^{-in_{1}\omega\sigma-i\omega t},$ (75a) $\displaystyle
E_{II}(\sigma,t)=\xi_{II}e^{in_{2}\omega(\sigma-L_{I})-i\omega
t}+\eta_{II}e^{-in_{2}\omega(\sigma-L_{I})-i\omega t},$ (75b)

where $\xi_{I},\eta_{I},\xi_{II},\eta_{II}$ are constants. Analogously, using
the relationship $H=nE$, one gets the magnetic fields. One could treat such a
system by considering a generic $\omega$, however, as we discussed in Ref.
[40], the axion-generated oscillations would be suppressed by a factor
$\delta(\omega-\omega_{a})$, so the only relevant frequency is $\omega_{a}$
for our purposes. Then, we take $\omega=\omega_{a}$ in the following. By
imposing the boundary conditions at the junctions to be closed, we found it
for particular cases and here we deal with a particular case of interest,
namely when $n_{2}$ becomes large in comparison to $n_{1}$ so that the ratio
$x\rightarrow 0$ and we consider the case where the length $L_{I}\rightarrow
0$, corresponding to a kind of point defect sitting on an otherwise uniform
string. As $\delta_{1}\rightarrow 0$ in this case, we see that

$\xi_{I}=\frac{E_{0}}{2\varepsilon_{1}},\quad x\rightarrow 0,L_{I}\rightarrow
0,$ (76)

which is a real quantity. We also have for this last case:

$\displaystyle\eta_{I}=\xi_{I},$ (77) $\displaystyle\eta_{II}=\xi_{II}=0.$
(78)

We have then an electromagnetic energy density $u_{1def}$ inside the point
defect:

$u_{1def}=\frac{E_{0}^{2}}{4\epsilon_{1}}.$ (79)

It is then worth mentioning the similarities and differences with the system
of a single dielectric plate treated by Millar et al. [15]. Indeed our system
is like the one in [15], but when the two outward dielectric media are closed
to form a ring. That suggestes the possibility that the case of a closed ring
with $N$ multiple dielectric plates could be of interest, since we could
expect to have an enhancement of $\sim N^{2}$ for $u_{1def}$ as in the setup
of [15].

The treatment of a ring with multiple regions (e.g. multiple point defects) is
surely very complicated analitically and would need a numerical treatment.
This can be subject of future works, along with other possible geometries,
such a cylindrical one with multiple conducting or dielectric plates, as
mentioned in [40]. Here we just give an order of magnitude of the expected
radiation pressure $p$. If we consider the behaviour of $p$ going as $\sim
N^{2}$, we have:

$p\sim N^{2}\frac{E_{0}^{2}}{4\epsilon_{1}}.$ (80)

If we take from the system in Ref. [15] with the following orientative values
$B_{e}\sim$10\text{\,}\mathrm{T}$$, $\Theta_{0}\sim 10^{-21}$,
$\epsilon_{1}\sim 1$, $N\sim 80$, one gets:

$p\sim 3.2\times
10^{-13}\,\frac{1}{\epsilon_{1}}\,\left(\frac{\Theta_{0}}{10^{-21}}\right)^{2}\,\left(\frac{N}{80}\right)^{2}\left(\frac{B_{e}}{$10\text{\,}\mathrm{T}$}\right)^{2}\,$\text{\,}\mathrm{P}\mathrm{a}$\,.$
(81)

It is worth noticing that for the orientative values we get a pressure that is
just one order of magnitude lower than the lowest pressures gotten in
laboratory (see e.g. Refs. [41, 42]), but it would be possible to get or
overcome those values of pressure with a bigger number $N$.

##  VI Optical properties of a toy model and axion domain walls

The purpose of this subsection is to study the propagation of an
electromagnetic wave in the presence of a toy model for the axion domain wall.
The toy wall is taken to be a planar sheet in which the axion field increases
linearly with the longitudinal coordinate $z$ between $z=0$ and $z=L$ and it
is constant and uniform outside. We assume a static wall, so $a=a(z)$ only,
thus no time dependence. As mentioned above, we first consider a space-
dependent axion field with planar symmetry:

$\Theta(z)=\begin{cases}0\;\;\;\textit{if}\;\;\;z<0,\\\
\frac{\Theta_{0}}{L}z\;\;\;\textit{if}\;\;\;0<z<L,\\\
\Theta_{0}\;\;\;\textit{if}\;\;\;z>L.\end{cases}$ (82)

The model was treated earlier in Ref. [37]: we will here make some
corrections. It serves as a toy model for ’localized’ axion configurations
since space-dependent axion configurations have typical lenghts $L\sim
m_{a}^{-1}$.

Figure 3: Explicit graph of the static axion domain wall configuration
$\Theta(z)=g_{a\gamma\gamma}a(z)$ where $a(z)=\gamma_{0}(z)+\alpha_{0}(z)$ is
reported in solid line in blue, evaluated as the ones in Ref. [21].We reported
also the correspondent values for $\gamma_{0}(z)$ and $\alpha_{0}(z)$.

We will now investigate the axion modifications of the electromagnetic wave in
the system. We mention in passing that possible other approximation could be:

$\Theta(z)=\frac{\Theta_{0}}{2}\left[\tanh\left(\frac{z}{L}\right)-1\right].$
(83)

As shown in Ref. [21] the ’effective field’ in Axion-Modified Maxwell
equations for an actual cosmological static axion domain wall gets two
contributions: one $\alpha_{0}(z)$ from the phase of Peccei-Quinn field, i.e.
the axion field and $\gamma_{0}(z)$ from the neutral pion field, since the
pion has an interaction term to electromagnetic field analogous to Eq. (1).
The two contributions are related because of the equilibrium condition for a
static wall:

$\tan{\gamma_{0}}=-\xi\tan{\alpha_{0}},$ (84)

where $\xi=\frac{m_{u}-m_{d}}{m_{u}+m_{d}}$, where $m_{u}$ is the mass of the
up quark and $m_{d}$ is the mass of down quark. In our numerical evaluation
here we used the value $\xi=0.3$ as done by Ref. [21] and we evaluated
$\alpha_{0}(z)$ by solving numerically the equation:

$\frac{1}{m_{a}^{2}}\frac{d^{2}\alpha_{0}}{dz^{2}}=\frac{\sin{\alpha_{0}}\cos{\alpha_{0}}}{\sqrt{\cos^{2}{\alpha_{0}}+\xi^{2}\sin^{2}{\alpha_{0}}}},$
(85)

obtained by the same Ref. [21]. We plot one of the possible solutions of Eq.
(85) in Figure 3. It is worth noticing that it is a ’localized’ axion
configuration whose typical lenght is of the order of $m_{a}^{-1}$.

###  VI.1 Exact calculations for toy model

As above, we consider an axion configuration that is dependent on $z$ only. At
first, we assume $\beta(z)=\partial_{z}\Theta(z)$ to be an arbitrary function
of $z$. The incident electromagnetic wave is propagating along the $z$-axis
and is transversely polarized (the same configuration was assumed in Ref.
[21]). This configuration is useful because a domain wall is invariant under
Lorentz boosts parallel to the wall surface, so it is always possible to find
a reference frame where the wave is incident normally, as also noticed in the
same [21].
This simplifies the calculations to the first order in the 4-potential because
the effective axion charge density is zero. We then consider the following
equation:

$\Box\vec{A}=-\nabla\Theta\wedge\frac{\partial\vec{A}}{\partial t},$ (86)

and search for an exact solution. We calculate the FT of this equation in time
$t$ and in the $x$ and $y$ coordinates, and get the following expressions for
the FT components of the vector potential:

$\begin{split}(-\partial_{z}^{2}-\kappa^{2})\mathcal{A}_{x}&=i\omega\beta(z)\mathcal{A}_{y},\\\
(-\partial_{z}^{2}-\kappa^{2})\mathcal{A}_{y}&=-i\omega\beta(z)\mathcal{A}_{x},\\\
(-\partial_{z}^{2}-\kappa^{2})\mathcal{A}_{z}&=0.\end{split}$ (87)

When $k_{x}=k_{y}=0$, we have $\kappa=\omega$. Equations (87) imply that the
gauge fields defined as $\mathcal{A}_{\pm}=\mathcal{A}_{x}\pm
i\mathcal{A}_{y}$ satisfy:

$\begin{split}(-\partial_{z}^{2}-\kappa^{2})\mathcal{A}_{\pm}&=\pm\omega\beta(z)\mathcal{A}_{\pm}.\end{split}$
(88)

From now on we assume to deal with the toy model configuration (82), so
$\beta(z)=0$ outside the wall, while $\beta(z)=\frac{\Theta_{0}}{L}=\beta_{z}$
inside the wall. We notice that the axion configuration (82) has not well-
defined values of the z derivative at $z=0,L$ since the derivative has a
Heaviside behaviour in such points. However, this is not problematic for
treating Eq. (86) since it just tells us that the second z-derivative of
vector potential has a Heaviside behaviour, but vector potential and its first
z-derivative are continuous.

Indeed, in the following calculations we basically solve Eq. (86) for the
regions $z<0$, $z>L$ and $0<z<L$ (where $\beta(z)$ is well-defined) and then
impose continuity of vector potential and its first z-derivative. This last
procedure would be also good in the case one considers axion configurations
like (83) where there would not be Heaviside-like behaviours.

Moreover, we consider a plane wave propagating along z. Then, the equations
(88) can be interpreted as follows. For the field $\mathcal{A}_{-}$ we have a
usual harmonic oscillator with frequency
$\alpha_{-}=\sqrt{\omega^{2}+\omega\beta(z)}$ inside (and outside) the domain
wall, while for the field $\mathcal{A}_{+}$ we have a damped wave inside the
domain wall if $\omega>\beta_{z}$. The gauge fields $\mathcal{A}_{\pm}$
correspond to left/right circular polarizations, so the basic properties are
consequently different for a left circular polarized wave and a right circular
polarized one. We now evaluate the reflected component of the incident
orthogonal wave. Recall the geometry: the leftmost region is $z<0$, the slab
region is $0<z<L$, and the rightmost region is $z>L$. The continuity
conditions at the walls are that the $A_{\pm}$ fields as well as their $z$
derivatives are continuous. Thus in the leftmost region,

$\mathcal{A}_{\pm}(z)_{\text{left}}=A_{\pm}e^{-i\kappa z}+B_{\pm}e^{i\kappa
z},$ (89)

where the $B$ term is the reflected wave. The wave transmitted outside the
slab is

$\mathcal{A}_{\pm}(z)_{\text{right}}=C_{\pm}e^{-i\kappa z}.$ (90)

To simplify the notation, we write $\alpha$ instead of $\alpha_{\pm}$, and we
consider $\alpha$ to be the frequency inside the domain wall:

$A_{\pm}(z)_{\text{in}}=A^{\prime}_{\pm}e^{-i\alpha
z}+B^{\prime}_{\pm}e^{i\alpha z}.$ (91)

where $A^{\prime}_{\pm}$ is the wave transmitted inside the wall. In order to
evaluate the reflection coefficient (defined as
$R_{\pm}=\frac{B_{\pm}}{A_{\pm}})$ and the internal transmission coefficient
(defined as $T_{\pm}=\frac{A^{\prime}_{\pm}}{A_{\pm}}$), we need to impose the
continuity conditions on the surfaces $z=0$ and $z=L$:

$\begin{cases}A+B=A^{\prime}+B^{\prime}&\text{Continuity of
$\mathcal{A}_{\pm}$ at z=0},\\\ \omega A-\omega B=\alpha A^{\prime}-\alpha
B^{\prime}&\text{Continuity of first derivatives of $\mathcal{A}_{\pm}$ at
z=0},\\\ e^{-i\alpha L}A^{\prime}+e^{i\alpha L}B^{\prime}=Ce^{-i\omega
L}&\text{Continuity of $\mathcal{A}_{\pm}$ at z=L},\\\ -i\alpha e^{-i\alpha
L}A^{\prime}+i\alpha e^{i\alpha L}B^{\prime}=-i\omega Ce^{-i\omega
L}&\text{Continuity of first derivatives of $\mathcal{A}_{\pm}$ at
z=L}.\end{cases}$ (92)

From it we get

$\displaystyle T$ $\displaystyle=\frac{2e^{2i\alpha
L}\omega(\alpha+\omega)}{(-1+e^{2i\alpha L})\alpha^{2}+(-1+e^{2i\alpha
L})\omega^{2}+2(1+e^{2i\alpha L})\alpha\omega},$ (93a) $\displaystyle R$
$\displaystyle=-\frac{(-1+e^{2i\alpha
L})(\alpha^{2}-\omega^{2})}{(-1+e^{2i\alpha L})\alpha^{2}+(-1+e^{2i\alpha
L})\omega^{2}+2(1+e^{2i\alpha L})\alpha\omega}.$ (93b)

Some fundamental properties of the reflection coefficient are shown in Figs. 4
and 5. Modules and phases are given, where we have taken the unit frequency to
be $1/L$ and $\Theta_{0}=10^{-3}$. We notice that there is most reflection,
for both modes, when $\omega=0$ and $\omega=\frac{\Theta_{0}}{L}$ but there is
still reflection at higher frequencies at intervals of
$\omega=\frac{\Theta_{0}}{L}$.

Figure 4: Plot of reflection coefficient modulus for + (red) and -
(blue)polarizations. Figure 5: Plot of reflection coefficient phase for +
(orange) and - (blue) polarizations.

###  VI.2 Dispersion relation

In the previous section we obtained the following dispersion relation for the
toy model (82):

$\omega^{2}-k_{z}^{2}=\mp\omega\beta(z),$ (94)

implying

$\omega_{\pm}=\sqrt{\left(k_{z}^{2}+\frac{1}{4}\beta^{2}(z)\right)}\pm\frac{1}{2}\beta(z).$
(95)

This is the same dispersion relation as obtained in [37]. However, this is not
the general dispersion relation because, as discussed in the previous section,
it applies to a single linear polarization only in its reference frame where
the wave vector lies along z direction. Applying an inverse Lorentz boost (so
in x and y directions) we can generalize the relation to take the following
form, similar to the one obtained in Ref. [23]:

$\omega=\sqrt{k_{x}^{2}+k_{y}^{2}+\left(\sqrt{\left(k_{z}^{2}+\frac{1}{4}\beta^{2}_{z}\right)}\pm\frac{1}{2}\beta_{z}\right)^{2}}.$
(96)

##  VII Space-dependent axion field. Casimir aspects

We now consider the simplest case of a planar symmetry in which the axion
field is constant in time but depends on the longitudinal coordinate z only.
The gradient of the field is thus directed along the $z$-axis; we will denote
it $\nabla_{z}\Theta=\beta(z)$. We get then the following equations for the
reduced Green’s function in component form:

$\displaystyle-(\partial_{z}^{2}+\omega^{2}-k_{y}^{2})g_{xx}+ik_{x}(\partial_{z}g_{zx}+ik_{y}g_{yx})-i\omega\beta(z)g_{yx}=-\delta(z-z^{\prime}),$
(97a)
$\displaystyle-(\partial_{z}^{2}+\omega^{2}-k_{x}^{2})g_{yy}+ik_{y}(\partial_{z}g_{zy}+ik_{x}g_{xy})+i\omega\beta(z)g_{xy}=-\delta(z-z^{\prime}),$
(97b)
$\displaystyle-\kappa^{2}g_{zz}+\partial_{z}(ik_{x}g_{xz}+ik_{y}g_{yz})=-\delta(z-z^{\prime}),$
(97c)

and

$\displaystyle-(\partial_{z}^{2}+\omega^{2}-k_{x}^{2})g_{yx}+ik_{y}(ik_{x}g_{xx}+\partial_{z}g_{zx})+i\omega\beta(z)g_{xx}=0,$
(98a)
$\displaystyle-(\partial_{z}^{2}+\omega^{2}-k_{y}^{2})g_{xy}+ik_{x}(ik_{y}g_{yy}+\partial_{z}g_{zy})-i\omega\beta(z)g_{yy}=0$
(98b)
$\displaystyle-(\partial_{z}^{2}+\omega^{2}-k_{y}^{2})g_{xz}++ik_{x}(ik_{y}g_{yz}+\partial_{z}g_{zz})=0,$
(98c)
$\displaystyle-\kappa^{2}g_{zi}+\partial_{z}(ik_{y}g_{yi}+ik_{x}g_{xi})=0\text{
if $i\neq z$},$ (98d)

as well as other analogous expressions. We here develop briefly the case of
the domain wall.

###  VII.1 Casimir force on domain wall

We can manipulate the equations (97) and (98) and the other equations by some
preliminar manipulations. Combining (97 a),(98 a) and (98) we can get the two
following coupled equations for $g_{xx}$ and $g_{yy}$:

$\begin{split}\omega^{2}\left(\partial_{z}^{2}+\kappa^{2}+i\beta(z)\frac{k_{x}k_{y}}{\omega}\right)g_{xx}+i\omega\beta(z)(\omega^{2}-k_{x}^{2})g_{yx}&=(\omega^{2}-k_{x}^{2})\delta(z-z^{\prime}),\\\
\omega^{2}\left(\partial_{z}^{2}+\kappa^{2}-i\beta(z)\frac{k_{x}k_{y}}{\omega}\right)g_{yx}-i\omega\beta(z)(\omega^{2}-k_{y}^{2})g_{xx}&=-k_{x}k_{y}\delta(z-z^{\prime}).\end{split}$
(99)

It is convenient to define $g^{e}_{ij}=-\omega^{2}g_{ij}$ and we notice it is
also convenient to define the functions

$g_{\pm}(z,z^{\prime})=g^{e}_{xx}(z,z^{\prime})+\gamma_{\pm}(\tilde{k},\omega)g^{e}_{yx}(z,z^{\prime}),$
(100)

where

$\gamma_{\pm}(\tilde{k},\omega)=\frac{k_{x}k_{y}\pm
i\omega\kappa}{\omega^{2}-k_{y}^{2}}.$ (101)

We then get the following simpler equations:

$\begin{split}\left[\partial_{z}^{2}+\kappa^{2}\pm\beta(z)\kappa\right]g_{\pm}=\left[\frac{\omega^{2}\kappa^{2}}{\omega^{2}-k_{y}^{2}}\mp\frac{i\omega\kappa
k_{x}k_{y}}{\omega^{2}-k_{y}^{2}}\right]\delta(z-z^{\prime}).\end{split}$
(102)

From solving Eq. (102) it is then straightforward to find $g_{xx}$ and
$g_{yx}$ and from these two can find $g_{zx}$.
It is also easy to find an analogous111It is trivial that it is analogous
because the system is invariant under rotation around $z$-axis result with
$g_{yy}$, $g_{xy}$ and $g_{zy}$:

$\begin{split}\left[\partial_{z}^{2}+\kappa^{2}\pm\beta(z)\kappa\right]\tilde{g}_{\pm}=\left[\frac{\omega^{2}\kappa^{2}}{\omega^{2}-k_{x}^{2}}\mp\frac{i\omega\kappa
k_{x}k_{y}}{\omega^{2}-k_{x}^{2}}\right]\delta(z-z^{\prime}),\end{split}$
(103)

where

$\tilde{g}_{\pm}(z,z^{\prime})=g^{e}_{yy}(z,z^{\prime})+\tilde{\gamma}_{\pm}(\tilde{k},\omega)g^{e}_{xy}(z,z^{\prime}).$
(104)

We notice then that for expressions such as $\tilde{g}_{\pm}(z,z^{\prime})$
and $g_{\pm}(z,z^{\prime})$ we need from the continuity of Green’s function
and their equation to have for $z,z^{\prime}<0$:

$g_{\pm}(z,z^{\prime})\propto\frac{1}{2\kappa}(e^{-i\kappa|z-z^{\prime}|}+R_{\pm}(\kappa)e^{i\kappa(z+z^{\prime})}),$
(105)

where $R_{\pm}(\kappa)$ are exactly the reflection coefficients (93). Then, if
we calculate the zz component of Energy-momentum tensor by expression (28), we
can get after a lot of calculations the temperature-dependent Casimir force
per unit area:

$f(T,L)=-T\sum_{\pm}\sum_{m}^{+\infty\,^{\prime}}\int_{\zeta_{m}}^{+\infty}\kappa^{2}d\kappa\,\frac{|R_{\pm}(\kappa)|^{2}e^{-2L\sqrt{q^{2}\pm
q\beta_{z}}}}{1-|R_{\pm}(\kappa)|^{2}e^{-2L\sqrt{q^{2}\pm q\beta_{z}}}}.$
(106)

In the following Subsection, we discuss its high temperature limit ($L^{-1}\ll
T$), that can be obtained by considering the first term $m=0$ in the sum and
it is equal to:

$\begin{split}f(T,L)=2\int\frac{d^{3}k}{(2\pi)^{3}}\left[R^{2}_{+}(k_{z})+R^{2}_{-}(k_{z})\right]\frac{k^{2}_{z}}{\omega}\frac{1}{e^{\beta\omega}-1}H(k_{z})=\\\
=2\int\frac{d^{2}k}{(2\pi)^{3}}\int_{0}^{+\infty}dk_{z}\left[R^{2}_{+}(k_{z})+R^{2}_{-}(k_{z})\right]\frac{k^{2}_{z}}{\omega}\frac{1}{e^{\beta\omega}-1},\end{split}$
(107)

It can be demonstrated by calculating the last integral in the variables
$k_{x}$ and $k_{y}$. However, it is very hard to develop the calculations with
the exact explicit expressions (93) of $R^{2}_{+}(\kappa)$. In the following,
we approximate the reflection coefficients $R_{\pm}(\kappa)$ in a similar way
as in Ref. [21] and also based on our results for the toy model. For
simplicity we define the useful quantity $m_{L}=\frac{\Theta_{0}}{L}$.

The reflection coefficients can be approximated to be
$R^{2}_{+}(k_{z})\sim\zeta m_{L}\,\delta(k_{z})$ and
$R^{2}_{-}(k_{z})\sim\zeta\,m_{L}\,\left[\delta(k_{z}+m_{L})+\delta(k_{z}-m_{L})\right]$
where $\zeta$ is a numerical factor such that fits better the behaviour of the
reflection coefficients. This can be justified by noticing that, although the
exact behaviours of the square modula of $R_{\pm}(k_{z})$ are complicated, it
is clear that a crucial frequency for our system is the same $m_{L}$, and by
the graphs in Figures 4 and 5 and a qualitative description of them, they can
be approximated to be
$R^{2}_{+}(k_{z})\sim\frac{\zeta}{\pi}\Theta_{0}^{2}m_{L}^{2}\frac{\sin^{2}{(k_{z}/m_{L})}}{k^{2}_{z}}$
and
$R^{2}_{-}(k_{z})\sim\frac{\zeta}{\pi}\Theta_{0}^{2}\frac{\sin^{2}{(k_{z}/m_{L}-1)}}{(k_{z}/m_{L}-1)^{2}}$.
Then one exploits the distributional relation:

$\frac{\sin^{2}{(\epsilon
t)}}{t^{2}}\sim\pi\epsilon\delta(t)\qquad\text{for}\qquad\epsilon\rightarrow
0,$ (108)

and we can get the expressions if $\epsilon=(m_{L})^{-1}$.

The use of the Dirac delta limit is only useful for simplifying the
calculations, however, an exact calculation would correspond to taking care of
the secondary peaks of the graphs in Figure 4 along with the detailed
structure of them and also of the main peak (the last one is namely $k_{z}=0$
for $R_{-}(k_{z})$ and $k_{z}=m_{L}$ for $R_{+}(k_{z})$) just modify the
numerical factor $\zeta$. The Dirac delta approximation then takes care of the
main contribution coming from the main peak of Figure 4, approximating it as a
strict peak.

###  VII.2 Evaluation of the temperature-dependent electromagnetic pressure on
domain wall

Figure 6: Visual representation useful for the calculation of radiation
pressure acting on an axion domain wall

As we mentioned previously, the calculation of the Temperature dependent
Casimir pressure acting on an axion domain wall is not immediate. We can
anyway get a first physical idea by getting the expression (107), making a
calculation similar to the one made by Huang and Sikivie [21] and to the one
made in Ref. [22], that also highlights the physical meaning and let us extend
the result to non-static walls.
Our results can be first applied to QCD axion domain walls, but also to ALP
domain walls [22] and overall strong topological insulators [23, 27, 43].
Namely, the last cited application can be the same system treated in Ref.
[23], if we consider the toy model configuration for the axion field (82), but
without the presence of two conducting plates, so the comparison between our
results and the ones in the cited reference is straightforward.

The idea is basically to calculate kinetically the radiation pressure acting
on the toy model (82) at temperature $T$, whose results can be also indicative
for cosmological walls, as we discussed before. If we have a circularly
polarized electromagnetic wave incident on the axion domain wall this is
partly reflected as discussed in Section VI. We refer to Figure 5 for the
geometry of the system.

If the incident wave (that can be left-circularly or right-circularly
polarized) has momentum density that is equal to its energy density
$u_{em,\pm}(\omega)$ it is partly reflected back with momentum density that is
equal to $R^{2}_{\pm}(k_{z})\,u_{em,\pm}(\omega)$ and partly transmitted with
momentum density equal to $T^{2}_{\pm}(k_{z})\,u_{em,\pm}(\omega)$, depending
on its polarization.

We preliminarly observe that, as also visible from what gotten for the toy
model, the reflection and transmission coefficients can only depend on
$k_{z}$, since the motion parallel to the domain wall cannot affect the
dynamics as the domain wall is invariant for parallel boosts. It then
experiences a variation of momentum equal to $\delta t\,\delta
A\,\left[2R^{2}_{\pm}(k_{z})\cos{\theta}^{2}\right]u_{em,\pm}(\omega)$ where
$\delta A$ is the differential area and $\delta t$ is the differential time.
That expression is consequence of the relation
$R^{2}_{\pm}(k_{z})+T^{2}_{\pm}(k_{z})=1$, so
$1+R^{2}_{\pm}(k_{z})-T^{2}_{\pm}(k_{z})=2R^{2}_{\pm}(k_{z})$. If we divide by
$\delta A$ and $\delta t$ we get the pressure exerted from that mode with that
specific polarization.

Summing on the two polarization and all frequencies, we get the following
pressure $\mathcal{P}_{L}$ acting on the left of the static wall, similarly to
what is obtained in Ref. [22]:

$\begin{split}\mathcal{P}_{L}=\mathcal{P}=2\int\frac{d^{3}k}{(2\pi)^{3}}\left[R^{2}_{+}(k_{z})+R^{2}_{-}(k_{z})\right]\frac{k^{2}_{z}}{\omega}\frac{1}{e^{\beta\omega}-1}H(k_{z})=\\\
=2\int\frac{d^{2}k}{(2\pi)^{3}}\int_{0}^{+\infty}dk_{z}\left[R^{2}_{+}(k_{z})+R^{2}_{-}(k_{z})\right]\frac{k^{2}_{z}}{\omega}\frac{1}{e^{\beta\omega}-1},\end{split}$
(109)

which is equal to the pressure $\mathcal{P}_{R}$ acting on the right of the
wall because a static domain wall is left-right symmetric. $H$ is the
Heaviside function, taking care of just the contribution of photons coming
from left and here $\omega=\sqrt{k^{2}_{x}+k^{2}_{y}+k^{2}_{z}}$, so we
exploited the relation $k_{z}=\omega\cos{\theta}$.

The calculation of $\mathcal{P}$ is not trivial because it is dependent on the
details of the axion configuration, from which reflection coefficients depend.
For the case of the toy model we can evaluate it by using the obtained
explicit expressions. In this way, the integral (109) becomes explicitly as:

$\mathcal{P}_{L}=2m_{L}^{3}\int\frac{d^{2}k}{(2\pi)^{3}}\frac{1}{\sqrt{k^{2}+m_{L}^{2}}}\frac{1}{e^{\beta\sqrt{k^{2}+m_{L}^{2}}}-1},$
(110)

since one can use the parity property of the integrand under the change of
variable $k_{z}\rightarrow-k_{z}$ and the defining property of the Dirac delta
$\int_{-\infty}^{\infty}dxf(x)\,\delta(x-x_{0})=f(x_{0})$. This integral can
be evaluated by noticing that the integrand is dependent on the variable
$k=\sqrt{k_{x}^{2}+k_{y}^{2}}$, on the physical quantities $m_{L}$ and
$\beta$, so we can use polar coordinates and making the substitution
$k_{i}\rightarrow k_{i}/m_{L}$ to get an integral just dependent on $\beta
m_{L}$:

$\mathcal{P}_{L}=\frac{m_{L}^{4}}{\pi^{2}}\int_{0}^{+\infty}dk^{\prime}\,\frac{k^{\prime}}{\sqrt{k^{\prime
2}+1}}\frac{1}{e^{\beta m_{L}\sqrt{k^{\prime 2}+1}}-1}.$ (111)

This integral can be evaluated by noticing that a primitive function
(indefinite) integral $F(x)$ of the function

$f(x)=\frac{x}{\sqrt{x^{2}+1}}\frac{1}{e^{\alpha\sqrt{x^{2}+1}}-1}$ (112)

is trivially the function
$F(x)=\alpha^{-1}\ln{(e^{\alpha\sqrt{x^{2}+1}}-1)}-\sqrt{x^{2}+1}$.
Consequently, it can be easily found that

$\mathcal{P}=\frac{m_{L}^{4}}{\pi^{2}}\left[1-\frac{T}{m_{L}}\ln{(e^{\beta
m_{L}}-1)}\right].$ (113)

We get then expressions agreeing partly with those mentioned in Ref. [21]:

$\displaystyle\mathcal{P}=\frac{m_{L}^{3}}{\pi^{2}}\,\,Te^{-\frac{m_{L}}{T}}\,\,\,\,\,\,\,\,\,\,\text{for
$T\ll m_{L}$},$ (114a)
$\displaystyle\mathcal{P}=\frac{m_{L}^{3}}{\pi^{2}}\,\,T\,\ln{\left(\frac{T}{m_{L}}\right)}\,\,\,\,\,\,\,\,\,\text{for
$T\gg m_{L}$}.$ (114b)

The first limiting case (114a) can be found by adopting the asymptotic
expression $\ln{(x-1)}\sim\ln{x}-\frac{1}{x}$ for $x\rightarrow+\infty$, where
$x$ is in our case $e^{\beta m_{L}}$. The second case (114b) can be obtained
by using $\ln{(e^{x}-1)}=\ln{x}+\frac{x}{2}+\mathcal{O}(x^{4})$ for
$x\rightarrow 0$ and leaving the more relevant term for $T\gg m_{L}$, i.e. the
logarithm $\ln{x}$.
We notice how our result (114b) is different from that claimed in Ref. [21],
since here we get a behaviour going like $\sim m_{L}^{3}T\ln(T/m_{L})$, that
is proportional to $\Theta_{0}^{3}$, which is different from a behaviour
$\sim\Theta_{0}^{2}m_{a}^{2}T^{2}$. This is a relevant difference and it is a
consequence of the different dependence on temperature $T$ and the observation
of the relevance of the threeshold frequency $m_{L}$, instead of $L^{-1}$,
where the reflection coefficients are significant.

It is also of interest to treat the case of a domain wall moving at constant
velocity $v$, namely along the direction of the $z$-axis, relatively to the
reference frame where the electromagnetic radiation emission is an isotropic
blackbody one at temperature $T$. In such a case, expressions for
$\mathcal{P}_{L}$ and $\mathcal{P}_{R}$ are similar to the ones in the static
case $v=0$, but there are two main differences. The first one is that we do
not expect $\mathcal{P}_{L}$ and $\mathcal{P}_{R}$ to be equal since the left-
right symmetry is broken in such a case, while the second one relies on the
need of taking care of the Doppler effect. Consequently we get, similarly to
Ref. [22],

$\mathcal{P}_{L,R}=\frac{m_{L}^{4}}{\pi^{2}}\int_{0}^{+\infty}dk^{\prime}\,\frac{k^{\prime}}{\sqrt{k^{\prime
2}+1}}\frac{1}{e^{\gamma(v)\beta m_{L}(\sqrt{k^{\prime 2}+1}\pm v)}-1}.$ (115)

where the $+$ is valid for $\mathcal{P}_{L}$ while $-$ for $\mathcal{P}_{R}$.
This integral can be calculated by noticing that the integrand has a primitive
function (indefinite integral) that is analogous to the one of static case.
One then gets the exact solution:

$\mathcal{P}_{L,R}=\frac{1}{\pi^{2}}m_{L}^{4}e^{\mp\gamma(v)\beta
m_{L}v}\left[1-\frac{T}{m_{L}}\ln{(e^{\beta m_{L}(1\pm v)}-1)}\right].$ (116)

We highlight, as mentioned before, that the pressures in the limit $\beta
m_{a}\gg 1$ can be evaluated by just substituting the Bose-Einstein
distribution with the limiting Boltzmann factor $e^{-\gamma(v)\beta
m_{a}(\sqrt{k^{\prime 2}+1}\pm v)}$, which is indeed a good approximation for
$\beta m_{a}\gg 1$, since it is surely much bigger than 1. This substitution
is also true in the ultrarelativistic limit $v\simeq 1$ for the same reason.
We have for such limiting cases ($\beta m_{a}\gg 1$ and/or $v\simeq 1$):

$\mathcal{P}_{L,R}\propto m_{L}^{3}T\sqrt{1-v^{2}}e^{-\gamma(v)\beta
m_{L}(1\pm v)}$ (117)

The non-relativistic regime $v\ll 1$ can be well approximated by the solution
(113) with $v=0$.

####  VII.2.1 Application to Axion Cosmology

In the case of interest for Axion cosmology $T\gg m_{L}$, if we consider a
static wall, we adopt the expression (114b) for such a case. This result gives
us a pressure:

$P_{em}^{0}\sim
10^{-27}\left(\frac{\Theta_{0}}{10^{-3}}\right)^{3}\left(\frac{m_{a}}{10^{-2}\,$\text{\,}\mathrm{e}\mathrm{V}$}\right)^{3}\left(\frac{T}{$1\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$}\right)\left[1+\ln{\left(10^{-14}\frac{T}{m_{a}}\right)}\right]\,$\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$^{4}.$
(118)

This very low pressure can be compared for example to gravitational pressure.
Indeed, according to Ref. [21] in the $N_{DW}=1$ case, a closed axion domain
wall oscillates and radiates gravitational energy and, if one neglects other
contributions, the ’equation of motion’ is:

$\frac{d\sigma_{s}}{dt}=-P_{grav}=-G\sigma_{s}^{2},$ (119)

where $G$ is the gravitational constant and $\sigma_{s}$ is the surface energy
density of the domain wall. According to the values reported in Ref. [11] and
the dynamics described by Eq. (119), we expect a life-time $\tau_{grav}$:

$\tau_{grav}\sim(G\sigma_{s})^{-1}\sim 2\times
10^{6}\,$\text{\,}\mathrm{s}$\left(\frac{10^{9}\,$\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$}{f_{a}}\right),$
(120)

and a pressure $P_{grav}$:

$P_{grav}\sim G\sigma^{2}_{s}\sim
0.8375\left(\frac{f_{a}}{10^{9}\,$\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$}\right)^{2}$\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$^{4}.$
(121)

This last pressure is much bigger than $P_{em}^{0}$ for every reasonable axion
mass, and with a temperature $T\simeq\Lambda_{QCD}$ it is much worse for an
actual domain wall since, according to Refs.[11, 26], two axion walls repulse
each other with acceleration $a=2\pi G\sigma_{s}$, so they reach a speed $c$
into a time $\tau\sim(2\pi G\sigma_{s})^{-1}$. That means that, if we
qualitatively estimate the electromagnetic radiation by considering
approximately the velocity-dependent pressure (gotten with zero acceleration)
we get fainter pressures that approach to zero when $v\rightarrow 1$. We can
analogously get a partial life-time associated to the electromagnetic
radiation pressure (118) that is

$\tau_{a\gamma\gamma}=\frac{\sigma_{s}}{P^{0}_{em}}\sim 35.75\times
10^{24}\left(\frac{\Theta_{0}}{10^{-3}}\right)^{-3}\left(\frac{m_{a}}{10^{-2}\,$\text{\,}\mathrm{e}\mathrm{V}$}\right)^{-3}\left(\frac{T}{$1\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$}\right)^{-1}\left[1+\ln{\left(10^{-14}\frac{T}{m_{a}}\right)}\right]^{-1}\,$\text{\,}\mathrm{s}$.$
(122)

This value is useful for comparisons with what one gets for $N_{DW}>1$
scenarios. In fact, as mentioned before, the cosmological scenario is
problematic for $N_{DW}>1$ where there is the possibility of a wall-dominated
Universe with an energy density inconsistent with observations. As shown e.g.
by Ref. [11], the typical times for such disaster events would be the times
$\tau_{2}$, where domain wall form, and the age of the Universe $\tau_{DW}$
where domain walls would dominate when $N_{DW}>1$.

The comparison of the times $\tau_{DW}$ and $\tau_{2}$ in the former reference
with the very high decay life-time in Eq. (122), that is far bigger than the
current age of the Universe, suggests us that the electromagnetic radiation
pressure is irrelevant for the evolution of cosmological axion domain walls.

##  VIII Comparison with Faraday Effect and chiral media and about the optical
activity in Axion Electrodynamics

It is worth mentioning that the mode splittings are analogous to those
encountered in the Faraday effect, i.e., the polarization rotation that is
proportional to the longitudinal strong magnetic field. For treatises of this
effect, cf., for instance, [44, 45]. Moreover, there is also a strong
connection with the Casimir polarization rotation observed in a chiral medium
when there is a strong transverse magnetic field present. Cf. [46, 47]. This
effect has been treated for Axion Electrodynamics, however with fixed axion
field that holds constant space and time derivatives, in Refs. [14] and [48].
It is indeed easy to notice how our theory with fixed axion background is
equivalent to a modified Electrodynamics with an additional Chern-Simons term:

$\mathcal{L}_{a}=-\frac{1}{2}p_{\alpha}A_{\beta}\tilde{F}^{\alpha\beta},$
(123)

if $p_{\alpha}=g_{a\gamma\gamma}\partial_{\alpha}a(x)$, regardless if axion
derivatives are constant and uniform or not.
The theory of this modified electrodynamics when $p_{\alpha}$ are constants is
very well treated in Ref. [48]. Here we mention the similarities and the
differences between our results and the ones in [48, 14], in order to
highlight what we get when considering more physically interesting axion
configurations, such as a coherent oscillating axion field and an axion domain
wall. As found in Refs. [48, 14], an Axion Electrodynamics theory with fixed
axion field that holds constant and uniform derivatives is characterized by
having an optically active vacuum

$\frac{d\Phi}{dt}=\frac{1}{2}g_{a\gamma\gamma}\sqrt{\frac{\mu}{\varepsilon}}\left(\dot{a}+\frac{\omega}{k^{2}}\vec{k}\cdot\nabla
a\right).$ (124)

This comes from the dispersion relation:

$\omega_{\pm}=\frac{|\vec{k}|}{\sqrt{\varepsilon\mu}}\pm\frac{1}{2}g_{a\gamma\gamma}\sqrt{\frac{\mu}{\varepsilon}}\left(\dot{a}+\frac{\omega}{|\vec{k}|^{2}}\vec{k}\cdot\nabla
a\right)+\mathcal{O}(g_{a\gamma\gamma}^{2}),$ (125)

which is compatible with the dispersion relations obtained in Subsections
(V.1) and (VI.2). However, this is not enough to grasp all the properties of
real axion backgrounds, such as an oscillating one and an axion domain wall.
These two are characterized by a typical frequency (in the case
$a(t)=a_{0}\sin(\omega_{a}t)$ it is the axion frequency $\omega_{a}$, while
for the toy model of domain wall it is $m_{L}=\frac{\Theta_{0}}{L}$). The
result is valid in the limit of frequencies and wave numbers much bigger than
the typical frequency, because they correspond to modes with very low
wavelengths or periods, so they do not ’appreciate’ the space-time variations
of the derivatives of axion field. It is worth noticing that the deviations we
got for the two systems: a domain wall is characterized by having a more
significant reflection coefficient at $k_{z}=m_{L}=\frac{\Theta_{0}}{L}$,
which means a very steady variation of the polarization plane near to the
interface, if compared to the variation of the axion field along the light
trail (that is arbitrarily small for an arbitrary small neighborhood of the
interface). Similarly, it happens with $a(t)=a_{0}\sin(\omega_{a}t)$ that, in
the forementioned phenomenon of axion echo , we have production of a fainter
radiation when inputting an electromagnetic wave with frequency
$\omega=\frac{\omega_{a}}{2}$, that lets the polarization plane rotate with a
behaviour not described by the theory with constant axion derivatives.
Furthermore, when generating a strong magnetic field there is production of
faint photons with frequency $\omega=\omega_{a}$ inside an electromagnetic
cavity (see Refs. [14, 15, 40]), so there is a phenomenon more similar to the
Faraday effect. Those arguments show how the validity of ’achromaticity of the
optical activity of Axion Electrodynamics’ is only valid in the regime of
’high frequencies’.

##  IX Conclusions

In this article we have treated the theory of zero-point energy modifications
in Axion Electrodynamics with a focus on physical systems of interest for the
High Energy Physics axion and the effective axion field in topological
insulators. We reviewed the theory of Axion Electrodynamics and particularly
the energy-momentum conservation in a linear dielectric and magnetic material.
Adopting the model of the oscillating axion background, we discussed the
axion-induced modifications to the Casimir force between two parallel plates
by using the Green’s function method and suggested a way of enhancing axion-
induced Casimir force through a closed string of dielectric point defects. We
calculated the radiation pressure acting on a axion domain wall at temperature
$T$. Our results for an oscillating axion field and a domain wall are useful
for ”axionic topological insulators” and experimental systems analogous to the
ones showed e.g. by Jiang, Q. D., & Wilczek, F. (2019) [46] and Fukushima et
al. (2019) [23]. Finally, we compared our results, where we assume time-
dependent or space-dependent axion configurations, with the discussion of the
optical activity of Axion Electrodynamics by Sikivie (2021) [14] and Carrol et
al. (1990) [48] and showed the ’cromaticity’ of the vacuum in Axion
Electrodynamics.

Note: During the revision of our paper, a recent work [49] has appeared on
arXiv, where the Authors discuss also finite-temperature Casimir force in
Axion Electrodynamics for other cases of interest.

## Acknowledgments

We thank the Erasmus project, the Norwegian University of Science and
Technology, and the University of Palermo, for the support of the Erasmus
Traineeship for Amedeo Maria Favitta. Part of this work was performed in
connection with his Master Thesis. We are grateful to Yuri N. Obukhov, Markku
Oksanen and Roberto Passante for valuable informations and several
suggestions.

## References

  * Peccei and Quinn [1977] R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977).
  * Peccei [2008] R. D. Peccei, The Strong CP Problem and Axions, in  _Axions: Theory, Cosmology, and Experimental Searches_, edited by M. Kuster, G. Raffelt, and B. Beltrán (Springer Berlin Heidelberg, 2008).
  * Chadha-Day  _et al._ [2021] F. Chadha-Day, J. Ellis, and D. J. E. Marsh, Sci. Adv. 8 (2021).
  * Bigi  _et al._ [2021] I. I. Bigi, G. Ricciardi, and M. Pallavicini,  _New Era for CP Asymmetries: Axions and Rare Decays of Hadrons and Leptons_ (World Scientific, 2021).
  * Weinberg [1978] S. Weinberg, Phys. Rev. Lett. 40, 223 (1978).
  * Schwartz [2014] M. D. Schwartz,  _Quantum field theory and the standard model_ (Cambridge University Press, 2014).
  * Preskill  _et al._ [1983] J. Preskill, M. B. Wise, and F. Wilczek, Phys. Lett. B 120, 127 (1983).
  * Abbott and Sikivie [1983] L. F. Abbott and P. Sikivie, Phys. Lett. B 120, 133 (1983).
  * Dine and Fischler [1983] M. Dine and W. Fischler, Phys. Lett. B 120, 137 (1983).
  * Amruth  _et al._ [lack] A. Amruth, T. Broadhurst, J. Lim, M. Oguri, G. F. Smoot, J. M. Diego, E. Leung, R. Emami, J. Li, T. Chiueh, _et al._ , Nature Astronomy , 1 (2023).
  * Sikivie [2008] P. Sikivie, Axion cosmology, in  _Axions: Theory, Cosmology, and Experimental Searches_ , edited by M. Kuster, G. Raffelt, and B. Beltrán (Springer Berlin Heidelberg, Berlin, Heidelberg, 2008) pp. 19–50.
  * Hiramatsu  _et al._ [2012] T. Hiramatsu, M. Kawasaki, K. Saikawa, and T. Sekiguchi, Phys. Rev. D 85, 105020 (2012).
  * Kawasaki  _et al._ [2015] M. Kawasaki, K. Saikawa, and T. Sekiguchi, Phys. Rev. D 91, 065014 (2015).
  * Sikivie [2021] P. Sikivie, Rev. Mod. Phys. 93, 015004 (2021).
  * Millar  _et al._ [2017] A. J. Millar, G. G. Raffelt, J. Redondo, and F. D. Steffen, Journal of Cosmology and Astroparticle Physics 01, 061 (2017).
  * Brevik and Chaichian [2022a] I. H. Brevik and M. M. Chaichian, International Journal of Modern Physics A 37, 2250151 (2022a), https://doi.org/10.1142/S0217751X22501512 .
  * Braggio  _et al._ [2022] C. Braggio, R. Calabrese, G. Carugno, G. Fiscelli, M. Guarise, A. Khanbekyan, A. Noto, R. Passante, L. Rizzuto, G. Ruoso, and L. Tomassetti, Applied Sciences 12, 6492 (2022).
  * Canfora  _et al._ [2022] F. Canfora, D. Dudal, T. Oosthuyse, P. Pais, and L. Rosa, Journal of High Energy Physics 2022, 1 (2022).
  * Oosthuyse and Dudal [lack] T. Oosthuyse and D. Dudal, Interplay between chiral medium and perfect electromagnetic conductor plates: repulsive vs. attractive casimir force transitions (2023), arXiv:2301.12870 [hep-th] .
  * Adshead  _et al._ [2020] P. Adshead, P. Draper, and B. Lillard, Phys. Rev. D 102, 123011 (2020).
  * Huang and Sikivie [1985] M. C. Huang and P. Sikivie, Phys. Rev. D 32, 1560 (1985).
  * Blasi  _et al._ [lack] S. Blasi, A. Mariotti, A. Rase, A. Sevrin, and K. Turbang, Journal of Cosmology and Astroparticle Physics 2023 (04), 008.
  * Fukushima  _et al._ [lack] K. Fukushima, S. Imaki, and Z. Qiu, Phys. Rev. D 100, 045013 (2019).
  * Kiefer [2009] C. Kiefer, Quantum gravity (general) and applications, in  _Compendium of Quantum Physics_ , edited by D. Greenberger, K. Hentschel, and F. Weinert (Springer Berlin Heidelberg, Berlin, Heidelberg, 2009) pp. 565–572.
  * Vilenkin and Shellard [1994] A. Vilenkin and E. P. S. Shellard,  _Cosmic strings and other topological defects_ (Cambridge University Press, 1994).
  * Ipser and Sikivie [1984] J. Ipser and P. Sikivie, Phys. Rev. D 30, 712 (1984).
  * Nenno  _et al._ [2020] D. M. Nenno, C. A. Garcia, J. Gooth, C. Felser, and P. Narang, Nature Reviews Physics 2, 682 (2020).
  * Zyuzin and Burkov [2012] A. A. Zyuzin and A. A. Burkov, Phys. Rev. B 86, 115133 (2012).
  * Landau and Lifshitz [2013] L. D. Landau and E. Lifshitz,  _Electrodynamics of continuous media_ (Elsevier, 2013).
  * Milton [1999] K. A. Milton, The Casimir Effect: Physical Manifestations of Zero Point Energy (World Scientific, 1999).
  * Birkeland and Brevik [2007] O. J. Birkeland and I. Brevik, Physical Review E 76, 066605 (2007).
  * Landau  _et al._ [1987] L. Landau, E. Lifshitz, and L. Pitaevskij,  _Statistical physics, part 2: Theory of the condensed state_ (Pergamon Press, 1987).
  * Ku  _et al._ [1997] W. Ku, F. Ge, and J. Zhu, Journal of Physics D: Applied Physics 30, 1842 (1997).
  * Kharlanov and Zhukovsky [2010] O. G. Kharlanov and V. C. Zhukovsky, Phys. Rev. D 81, 025015 (2010).
  * Brevik  _et al._ [2022] I. Brevik, B. Shapiro, and M. G. Silveirinha, International Journal of Modern Physics A 37, 2241012 (2022), https://doi.org/10.1142/S0217751X22410123 .
  * Stakgold and Holst [lack] I. Stakgold and M. J. Holst,  _Green’s functions and boundary value problems_ (John Wiley & Sons, 2011).
  * Brevik [2021] I. Brevik, Universe 7, 133 (2021).
  * Brevik and Chaichian [2022b] I. H. Brevik and M. M. Chaichian, The European Physical Journal C 82, 202 (2022b).
  * Arza and Sikivie [2019] A. Arza and P. Sikivie, Phys. Rev. Lett. 123, 131804 (2019).
  * Brevik  _et al._ [2023] I. Brevik, A. M. Favitta, and M. Chaichian, Phys. Rev. D 107, 043522 (2023).
  * Thompson and Hanrahan [1977] W. Thompson and S. Hanrahan, Journal of Vacuum Science and Technology 14, 643 (1977).
  * Redhead [1999] P. Redhead,  _Extreme high vacuum_ , Tech. Rep. (Cern, 1999).
  * Yan  _et al._ [2021] Q. Yan, H. Li, J. Zeng, Q.-F. Sun, and X. Xie, Communications Physics 4, 1 (2021).
  * Prati [2003] E. Prati, Journal of electromagnetic waves and applications 17, 1177 (2003).
  * Mansuripur [1999] M. Mansuripur, Optics and Photonics News 10, 32 (1999).
  * Jiang and Wilczek [2019] Q.-D. Jiang and F. Wilczek, Phys. Rev. B 99, 125403 (2019).
  * Høye and Brevik [2020] J. S. Høye and I. Brevik, The European Physical Journal Plus 135, 1 (2020).
  * Carroll  _et al._ [1990] S. M. Carroll, G. B. Field, and R. Jackiw, Phys. Rev. D 41, 1231 (1990).
  * Ema  _et al._ [lack] Y. Ema, M. Hazumi, H. Iizuka, K. Mukaida, and K. Nakayama, Zero casimir force in axion electrodynamics and new force search (2023), arXiv:2302.14676 [hep-ph] .