# Fermionic dynamical Casimir effect: Magnus expansion

C. D. Fosco and G. Hansen
Centro Atómico Bariloche and Instituto Balseiro
Comisión Nacional de Energía Atómica
R8402AGP S. C. de Bariloche, Argentina.

###### Abstract

We study pair creation out of the vacuum, for a system consisting of a massive
Dirac field in $1+1$ dimensions, contained between a pair of perfectly
reflecting boundaries, one of them oscillating. After analyzing some general
properties of the vacuum-decay process, we evaluate the corresponding
transition amplitude in a Magnus expansion of the $S$-matrix. We show how this
yields, besides the single-pair creation amplitude, multipair ones, as well as
corrections to the single pair amplitude.

We also apply it to obtain an approximate, yet explicitly unitary expression
for the Bogoliubov transformation between the in and out Fock spaces.

##  1 Introduction

Startling macroscopic phenomena arise when quantum fields are subjected to
non-trivial boundary conditions. Among them, the Casimir effect [1] stands out
as a prime example. In its static version, it manifests itself in the
existence of forces which reflect the dependence of the vacuum energy on the
boundary conditions. This subject of study is the so called static Casimir
effect (SCE) [2]. On the other hand, the dynamical Casimir effect [3] (DCE),
delves with time-dependent changes in the boundary conditions. Within the
general framework of the DCE, we consider here a Dirac field, in a situation
where the geometry of the boundaries changes, while the (bag model) conditions
they impose is unaltered. Under these assumptions, the DCE observables
manifest features of the space-time geometry swept by the boundaries. A
particularly interesting effect is, under the appropriate circumstances, the
creation of particles out of the vacuum state. Those particles are real quanta
of the field, their nature being inherited from that of the field upon which
the boundary conditions are imposed. For instance, photons are emitted when
the boundaries are perfect electromagnetic conductors, and fermion pairs when
bag conditions [4] are imposed upon a Dirac field. The latter is the kind of
system we are interested in, among other reasons because it may be relevant in
Condensed Matter Physics, since Dirac fields, specially in lower dimensions,
appear as ingredients of many continuum models [5].

In [6] the fermionic DCE for a massless Dirac field in $1+1$ dimensions,
satisfying bag conditions on two moving boundaries, has been been studied. In
[7] we used a perturbative approach, applicable to either massive or massless
Dirac fermions, to study fermion pair creation in an oscillating cavity, also
in $1+1$ dimensions, to the lowest non-trivial order in the amplitude of
oscillation. It is the aim of this paper to improve upon that calculation, by
performing a different, yet related expansion: firstly, we just assume the
part of the action taken as perturbation to be small (in comparison to the
unperturbed action). This approach keeps reparametrization invariance for the
curves swept by the boundaries, even when working at the lowest non-trivial
order in the new expansion. By the same token, relativistic corrections are
incorporated, since the perturbation involves a Lorentz factor. Secondly, we
use a Magnus expansion to derive transition probabilities, and interpret its
implications from the point of view of the vacuum decay amplitudes, and for
the mapping between in and out Fock spaces.

The paper is structured as follows: in Sect. 2 we define the model to be
considered, briefly reviewing its main aspects. Then, in Sect. 3, we discuss
some general aspects of the pair creation phenomenon, and in Sect. 4 we deal
with the perturbative evaluation of the decay amplitude, in the context of the
Magnus expansion. The latter is applied to obtain an explicitly unitary
approximation for the Bogoliubov transformation, in Sect. 5. Finally, in Sect.
6 we convey our conclusions.

##  2 The model

The model we consider here has already been studied, to the lowest non-trivial
order, in [7] (we refer to that work for the properties which here are applied
without demonstration). The system consists of a Dirac field
$(\psi,\bar{\psi})$ in $1+1$ dimensions which is subjected to bag boundary
conditions by means of its coupling to a ‘potential’ $V$, the role of which is
to impose the appropriate boundary conditions of the two borders (see below).
The real-time action $\mathcal{S}$ for the system is given by

${\mathcal{S}}(\bar{\psi},\psi;V)\;=\;\int
d^{2}x\,\bar{\psi}\,\big{(}i\not\\!\partial-m-V\big{)}\,\psi\;,$ (1)

where $m$ is the mass of the fermions. We adopt conventions such that Planck’s
constant $\hbar$ and the speed of light $c$, are both set equal to $1$. Space-
time coordinates are denoted by $x^{\mu}$ (with $\mu=0,1$ and $x^{0}=t$), and
the Minkowski metric is $(g_{\mu\nu})\equiv\text{diag}(1,-1)$. For Dirac’s
$\gamma$ matrices, we choose the representation $\gamma^{0}\equiv\sigma_{1}$,
$\gamma^{1}\equiv i\sigma_{3}$, where $\sigma_{i}\,(i=1,2,3)$ are the standard
Pauli’s matrices.

Finally, we deal with the form of $V$ which is required in order to impose the
desired boundary conditions, in the approach outlined in [8, 9]. The field is
assumed to be in the presence of two walls, labeled by $L$ and $R$, but,
contrary to what we did in [7], we now assume that $L$ moves while $R$ remains
static. The first moves along the trajectory $x^{1}=\eta(x^{0})$, while $R$
remains at the position $x^{1}=a$. This motion is constrained so that the
walls do not overlap, i.e., $\eta(x^{0})<a$.

Recalling [7, 10], the potential $V$ for this scenario assumes the form:

$V(x)=g_{L}\,\gamma^{-1}(u)\,\delta[x^{1}-\eta(x^{0})]+g_{R}\,\delta(x^{1}-a)\,,$
(2)

where $\gamma(u)=1/\sqrt{1-u^{2}}$ denotes the instantaneous Lorentz factor as
a function of $u=u(x^{0})\equiv\dot{\eta}(x^{0})$. The dimensionless constants
$g_{L}$ and $g_{R}$ determine the nature of the resulting boundary conditions:
values differing from $2$, result in imperfect boundary conditions, allowing
for the possibility of some current escaping from, or coming into, the cavity.
Conversely, setting $g_{L}=g_{R}=2$ establishes ‘bag’ boundary condition as
shown in [8]. For the sake of completeness, we summarize those conditions, for
the case at hand. To that end, we introduce the projectors:
$P_{L}(x^{0})=\frac{1+i\not\hat{n}_{L}(x^{0})}{2}$, for the moving boundary at
$L$, and $P_{R}=\frac{1+i\not\hat{n}_{R}}{2}$, for the static one at $R$.
Here, the unit normals are given by:

$\hat{n}_{L}(x^{0})\,=\,\gamma(\dot{\eta}(x^{0}))\,\big{(}\dot{\eta}(x^{0}),1\big{)}\;,\;\;\;\hat{n}_{R}\,=\,(0,-1)\;,$
(3)

and the boundary conditions become:

$\displaystyle
P_{L}(x^{0})\,\psi\big{(}x^{0},\eta{(x^{0})}^{+}\big{)}\,=\,0\,=\,P_{R}\,\psi(x^{0},a^{-})\;,$
(4)

where the $+$ and $-$ superscripts indicate ‘lateral’ limits, approaching the
respective point from the right and from the left, respectively.

##  3 Pair creation

We want to compute the probability of vacuum decay, due to the motion of the
$L$ wall. We make the usual assumptions of a finite temporal extent for the
external perturbation of a system, so that both the initial state (vacuum) and
the final state can be characterized in terms of the states of the Fock space
of a free Hamiltonian: the in and out spaces. In our case, that Hamiltonian
shall be the one corresponding to a situation where both boundaries are
static, and at the same distance before and after the motion occurs.

To be more explicit, we split the action as follows:

$\displaystyle\mathcal{S}=\mathcal{S}_{0}+\mathcal{S}_{I}\;,$ (5)

with

$\displaystyle\mathcal{S}_{0}(\bar{\psi},\psi)\equiv\int
d^{2}x\,\bar{\psi}(x)\,\big{(}\,i\not{\partial}-m-V_{0}(x)\,\big{)}\,\psi(x)\,,$
(6) $\displaystyle\mathcal{S}_{I}(\bar{\psi},\psi,\eta)\equiv-\int
d^{2}x\,\bar{\psi}(x)\,\varphi(x)\,\psi(x)\,,$ (7)

and:

$\displaystyle V_{0}(x)$ $\displaystyle\equiv
2\,\big{[}\delta(x^{1})\,+\,\delta(x^{1}-a)\big{]}\;,$ (8)
$\displaystyle\varphi(x)$ $\displaystyle\equiv
2\,\big{[}\gamma^{-1}(\dot{\eta}(x^{0}))\,\delta(x^{1}-\eta(x^{0}))\,-\,\delta(x^{1})\big{]}\;.$

Note that, contrary to our previous work with similar models, we are not
expanding the perturbation in powers of the departure $\eta$, keeping just the
lowest (linear) order term. On the contrary, we develop the expansion keeping
the structure for $\varphi$, and in principle also the $\gamma$ factor which
accounts for relativistic corrections.

The free piece, $\mathcal{S}_{0}$, determines the basis of eigenstates of the
Dirac operator with the boundary conditions which result from (4) when
$\eta\equiv 0$. In the case at hand, this implies finding the solutions to:

$(i\not{\partial}-m)\,\psi(x)=0\;,$ (9)

with the boundary conditions:

$\displaystyle\frac{1+i\gamma_{1}}{2}\,\psi(x^{0},0^{+})\,=\,0\,=\,\frac{1-i\gamma_{1}}{2}\,\psi(x^{0},a^{-})\;.$
(10)

A convenient basis of solutions to this system is the following:

$\displaystyle\psi_{p,\pm}(x^{1})=\frac{1}{\sqrt{m+aE_{p}^{2}}}\begin{pmatrix}\pm
E_{p}\sin(px^{1})\\\ p\cos(px^{1})+m\,\sin(px^{1})\end{pmatrix}\,,$ (11)

where the $+$ sign corresponds to positive ($E_{p}$) energy solutions and $-$
to negative ($-E_{p}$) energy ones. In both cases
$E_{p}\equiv\sqrt{p^{2}+m^{2}}$ and the values of $p$ are determined by the
transcendental equation:

$ma\,\text{sinc}(pa)+\cos(pa)=0\,,$ (12)

where the sinc function is, in our conventions, defined as:
$\text{sinc}\,x=\frac{\sin x}{x}$.

The solutions are orthogonal and normalized as follows:

$\int_{0}^{a}\,dx^{1}\,\psi_{p,\sigma}^{\dagger}(x^{1})\,\psi_{p^{\prime},\sigma^{\prime}}^{\vphantom{\dagger}}(x^{1})=\delta_{p,p^{\prime}}\,\delta_{\sigma,\sigma^{\prime}}\,.$
(13)

The spectrum of the system is discrete, and the interchange $p\to-p$ does not
yield an independent solution. Thus, it is sufficient to consider (as we do)
$p$ assuming non-negative values only.

When $m=0$, the spectrum becomes:
$|E_{p_{n}}|=p_{n}=(n+\frac{1}{2})\frac{\pi}{a}$ with $n=0,1,\ldots$ Solutions
are labeled by the index $p_{n}$, which is related to the particles’ spatial
momenta when the size of the cavity tends to infinity, i.e., in the limit
$a\to\infty$. For practical purposes, we will use the term momenta to refer
specifically to the index $p_{n}$.

Both the in and out spaces are isomorphic, and the associated annihilation and
creation operators correspond to a conveniently chosen basis of solutions to
the above equation. For instance, the mode expansion for the in field operator
(we omit the ‘in’ label in both the field and the operators),

$\psi(x)=\sum_{k}\big{[}\,b_{k}\,e^{-iE_{k}x^{0}}\,u_{k}(x^{1})+\,d_{k}^{\dagger}\,e^{iE_{k}x^{0}}\,v_{k}(x^{1})\,\big{]}\,,$
(14)

where $u_{k}(x^{1})\equiv\psi_{k,+}(x^{1})$ and
$v_{k}(x^{1})\equiv\psi_{k,-}(x^{1})$, with $b_{k}^{\dagger}$ and
$d_{k}^{\dagger}$ being creation operators of fermions and anti-fermions,
respectively. With the normalization we have chosen for the solutions to the
Dirac equation, the only non-vanishing anti-commutators among these operators
are:

$\\{b_{k}^{\vphantom{\dagger}}\,,\,b^{\dagger}_{p}\\}\;=\;\delta_{kp}\;,\quad\\{d_{k}^{\vphantom{\dagger}}\,,\,d^{\dagger}_{p}\\}\;=\;\delta_{kp}\;.$
(15)

Also note that a single non-negative discrete label is sufficient to
characterize either a fermion or anti-fermion state.

Before actually calculating the decay probability amplitudes, we see that,
since fermion number is conserved, the vacuum can only decay into one or more
fermion anti-fermion pairs. Let us first see how, for the action we are
considering, the knowing of just the single pair creation amplitude, all the
multiple pair ones can be obtained by Wick’s theorem combinatorics and
knowedge of the single pair decay amplitude. Indeed, assume that we want to
consider a rather general vacuum decay amplitude, from the vacuum to a certain
number, $n$ say, of fermion pairs. A rather straightforward calculation shows
that the reduction formulae [11] can be adapted to this system: indeed, the
connected transition amplitudes are determined by the corresponding Green’s
function. In this case ($n$ pairs), it then involves a number $n$ of $\psi$
and $\bar{\psi}$ legs. To obtain the amplitude from such a Green’s function,
one acts with the inverse of the free propagator (in our case: the one in the
presence of static walls) on the legs of the function, and then attaches the
respective free spinor. In other words, the procedure is rather similar to the
one used in order to extract scattering amplitudes in free space, the
difference being in fact that the free evolution and wave functions correspond
to a static cavity.

Since the full action of the system is quadratic in the Dirac field, Wick’s
theorem holds true for that Green’s function: it can be written in terms of
combinations invoving pairings, albeit with a contraction $S_{F}$, the exact
Dirac propagator in the presence of a moving boundary:

$S_{F}(x,y)\;\equiv\;\bra{0}T[\psi(x)\bar{\psi}(y)]\ket{0}\;,$ (16)

which satisfies:

$(i\not\\!\partial_{x}-m-V(x))\,S_{F}(x,y)=i\,\delta^{2}(x-y)\,.$ (17)

All the contributions to the $n$-pair amplitude can then be written in terms
of that propagator (which, in general, one cannot calculate exactly in closed
form). On the other hand, using the reduction formalism, the basic amplitude,
corresponding to the creation of a fermion anti-fermion pair, with labels $k$
and $p$, respectively, may be obtained as follows:

$\displaystyle\langle k,p;{\rm out}|0;{\rm in}\rangle\,=\,\int d^{2}x\int
d^{2}y\,\bar{u}_{k}(x^{1})\,e^{iE_{k}x^{0}}\times\,$
$\displaystyle(i\not\\!\partial_{x}-m-V_{0}(x^{1}))\,S_{F}(x,y)\,(i\,\overleftarrow{\not\\!\partial_{y}}+m+V_{0}(y^{1}))\,e^{iE_{p}y^{0}}\,v_{p}(y^{1})\,.$
(18)

This allows us to write in a more explicit form the corresponding exact (yet
formal) matrix elements of the transition amplitude $T$, ($S=I+iT$, $S$:
S-matrix):

$\displaystyle T_{fi}=\int dx^{1}\int
dy^{1}\;\bar{u}_{k}(x^{1})\,M(x^{1},y^{1})\,v_{p}(y^{1})$ $\displaystyle
M(x^{1},y^{1})=\int dx^{0}\int
dy^{0}\,e^{i(E_{k}x^{0}+E_{p}y^{0})}\,\Big{[}\varphi(x)\delta^{2}(x-y)\,+\,\varphi(x)S_{F}^{(0)}(x,y)\varphi(y)$
$\displaystyle+\,\int
d^{2}z\varphi(x)S_{F}^{(0)}(x,z)\varphi(z)S_{F}^{(0)}(z,y)\varphi(y)\,+\ldots\Big{]}\,.$
(19)

Here, we have denoted the initial and final states by $\ket{i}$ and $\ket{f}$
respectively. We shall have: $\ket{i}=\ket{0}$ and
$\ket{f}=b_{k\vphantom{p}}^{\dagger\vphantom{p}}d_{p\vphantom{k}}^{\dagger\vphantom{k}}\ket{0}$,
where $b_{p}^{\dagger}$ and $d_{k}^{\dagger}$ create fermions and anti-
fermions, respectively. $S_{F}^{(0)}$ is the Dirac propagator for the fermions
in the presence of the static boundaries, which is determined by the equation:

$(i\not\\!\partial_{x}-m-V_{0}(x))\,S_{F}^{(0)}(x,y)=i\,\delta^{2}(x-y)\;.$
(20)

Knowing the result for a single pair, it is possible to calculate the
probability amplitude for $n$-pairs. For example, for the creation of two
particle pairs, we consider the states $\ket{i}=\ket{0}$ and
$\ket{f}=b_{k_{1}\vphantom{p_{1}}}^{\dagger\vphantom{p}}d_{p_{1}\vphantom{k_{1}}}^{\dagger\vphantom{k}}b_{k_{2}\vphantom{p_{2}}}^{\dagger\vphantom{p}}d_{p_{2}\vphantom{k_{2}}}^{\dagger\vphantom{k}}\ket{0}$.
Using the reduction formula, the probability amplitude for the creation of two
fermion pairs can be expressed in terms of the single pair amplitude:

$\displaystyle\langle k_{1},p_{1},k_{2},p_{2};{\rm out}|0;{\rm in}\rangle\,=$
$\displaystyle\langle k_{1},p_{1};{\rm out}|0;{\rm in}\rangle\times\langle
k_{2},p_{2};{\rm out}|0;{\rm in}\rangle$ (21) $\displaystyle\,+$
$\displaystyle\langle k_{1},p_{2};{\rm out}|0;{\rm in}\rangle\times\langle
k_{2},p_{1};{\rm out}|0;{\rm in}\rangle\,.$

Thus, the transition amplitude for multiple pairs can be systematically
constructed from the single pair amplitude by considering all possible
combinations of pairings. In the next section, we evaluate the transition
probability within the framework of the interaction representation.

##  4 Perturbation theory

The $S$-matrix may be written as follows:

$S=1+i\,T=U_{I}(\infty,-\infty)\,,$ (22)

with $U_{I}(t_{f},t_{i})$ denoting the evolution operator, from $t_{i}$ to
$t_{f}$, in the interaction representation, while $T$ is the transition
matrix. The probability of fermion and anti-fermion pair production, will be
denoted as $P_{kp}$, in consonance with the indices characterizing the
respective momenta $p$ and $k$.

In the interaction picture, the Hamiltonian for the perturbation becomes:

$\displaystyle\mathcal{H}_{I}^{\prime}(x^{0})=2\,\gamma^{-1}(x^{0})\,\sum_{k,p}\big{[}b_{k\vphantom{p}}^{\dagger}b_{p\vphantom{k}}^{\vphantom{\dagger}}\,e^{i(E_{k}-E_{p})x^{0}}\,g_{kp}(\eta(x^{0}))+b_{k\vphantom{p}}^{\dagger}d_{p\vphantom{k}}^{\dagger}\,e^{i(E_{k}+E_{p})x^{0}}\,h_{kp}(\eta(x^{0}))$
(23) $\displaystyle-
d_{k\vphantom{p}}b_{p\vphantom{k}}\,e^{-i(E_{k}+E_{p})x^{0}}\,h_{kp}(\eta(x^{0}))-d_{k\vphantom{p}}^{\vphantom{\dagger}}d_{p\vphantom{k}}^{\dagger}\,e^{-i(E_{k}-E_{p})x^{0}}\,g_{kp}(\eta(x^{0}))\,\big{]}\;,$

where, for an arbitrary argument $x^{1}$, the functions appearing above are
given by:

$\displaystyle g_{kp}(x^{1})$
$\displaystyle\equiv\bar{u}_{k}(x^{1})u_{p}(x^{1})=-\bar{v}_{k}(x^{1})v_{p}(x^{1})\;,$
(24) $\displaystyle h_{kp}(x^{1})$
$\displaystyle\equiv\bar{u}_{k}(x^{1})v_{p}(x^{1})=-\bar{v}_{k}(x^{1})u_{p}(x^{1})\,.$
(25)

One can check that under the transposition of suffixes they satisfy:

$g_{kp}(x^{1})\,=\,g_{pk}(x^{1})\;,\;\;\;h_{kp}(x^{1})\,=\,-h_{pk}(x^{1})\,,$
(26)

and that, when $m\equiv 0$, they adopt a rather simple form:

$\displaystyle g_{kp}(x^{1})$
$\displaystyle=\frac{1}{a}\,\sin((k+p)\,x^{1})\,,$ (27) $\displaystyle
h_{kp}(x^{1})$ $\displaystyle=\frac{1}{a}\,\sin((k-p)\,x^{1})\,.$ (28)

For later use, we render also their approximate form for small amplitudes of
motion, namely, for $|\eta(t)|<<a$:

$g_{kp}(\eta(t))\,\simeq\frac{1}{a}\,(k+p)\,\eta(t)\;,\;\;h_{kp}(\eta(t))\,\simeq\frac{1}{a}\,(k-p)\,\eta(t)\,.$
(29)

To recall the relation between the order in which is calculated the
probability, and the one of the Dyson expansion for the $S$-matrix, given in
terms of the evolution operator,

$S\;=\;\sum_{n=0}^{\infty}S^{(n)}\;,\quad
S^{(n)}\;=\;U_{I}^{(n)}(+\infty,-\infty)\,,$ (30)

we see that:

$P_{kp}\,=\,\sum_{n=1}^{\infty}P_{kp}^{(n)}\;,\;\;P_{kp}^{(n)}\,=\,\sum_{l=1}^{n-1}{S_{kp}^{(n-l)}}^{*}\,S_{kp}^{(l)}\;,$
(31)

with $S_{kp}^{(l)}$ denoting the $l^{\rm th}$-order term in the matrix element
of the $S$-matrix.

In a previous work [7], we dealt with the lowest (second) order contribution
to the probability, which involves just the first order term in the Dyson
expansion.

Let us now see how the Dyson expansion above is related to the Magnus
expansion [12], which in this case allows for the evolution operator to be
expressed as an exponential:

$U_{I}(t,t_{0})\,=\,\exp\\{-i\,\Omega(t,t_{0})\\}\,,\quad\Omega(t_{0},t_{0})=0\;,$
(32)

where the operator $\Omega(t,t_{0})$ can be written as a series:

$\Omega(t,t_{0})=\sum_{n=1}^{\infty}\Omega_{n}(t,t_{0})\;,$ (33)

where the first few terms are:

$\displaystyle\Omega_{1}(t,t_{0})$
$\displaystyle=\int_{t_{0}}^{t}dt_{1}\,\mathcal{H}_{I}^{\prime}(t_{1})\,,$
(34) $\displaystyle\Omega_{2}(t,t_{0})$
$\displaystyle=\frac{1}{2}\,\int_{t_{0}}^{t}dt_{1}\,\int_{t_{0}}^{t_{1}}dt_{2}\,\left[\mathcal{H}_{I}^{\prime}(t_{1}),\mathcal{H}_{I}^{\prime}(t_{2})\right]\,,$
$\displaystyle\Omega_{3}(t,t_{0})$
$\displaystyle=\frac{1}{6}\,\int_{t_{0}}^{t}dt_{1}\,\int_{t_{0}}^{t_{1}}dt_{2}\,\int_{t_{0}}^{t_{2}}dt_{3}\,\big{(}\,[\mathcal{H}_{I}^{\prime}(t_{1}),[\mathcal{H}_{I}^{\prime}(t_{2}),\mathcal{H}_{I}^{\prime}(t_{3})]]$
$\displaystyle+[\mathcal{H}_{I}^{\prime}(t_{3}),[\mathcal{H}_{I}^{\prime}(t_{2}),\mathcal{H}_{I}^{\prime}(t_{1})]]\,\big{)}\,.$

An important aspect of this expansion for $\Omega(t,t_{0})$ is its hermiticity
at all orders; namely, unitarity of the evolution operator at every order. In
spite of the previous assertion, of course one can connect the terms in both
expansions, albeit, as we shall see, there is more content to be found than
making this connection. We see that equating the Dyson perturbative series for
$U_{I}(\infty,-\infty)$ (abbreviated as $U_{I}$) and the Magnus series,

$U_{I}=1+(-i)\,U_{I}^{(1)}+{(-i)}^{2}\,U_{I}^{(2)}+\ldots\,=\,e^{-i\,(\Omega_{1}+\Omega_{2}+\ldots)}\;,$
(35)

where $\Omega_{n}\equiv\Omega_{n}(\infty,-\infty)$, with the first terms of
the series given by:

$\displaystyle U_{I}^{(1)}$ $\displaystyle=\Omega_{1}\,,$ (36) $\displaystyle
U_{I}^{(2)}$ $\displaystyle=\Omega_{2}+\frac{1}{2!}\,{(\Omega_{1})}^{2}\,,$
(37) $\displaystyle U_{I}^{(3)}$
$\displaystyle=\Omega_{3}+\frac{1}{2!}\,(\Omega_{1}\Omega_{2}+\Omega_{2}\Omega_{1})+\frac{1}{3!}\,{(\Omega_{1})}^{3}\;.$
(38)

Let us evaluate now $\Omega_{1}$ and the corresponding contribution to the
transition amplitudes to the first non-trivial order.

###  4.1 First order

The first term of the Magnus series (33) becomes

$\displaystyle\Omega_{1}=\sum_{k,p}\big{(}\,G_{kp}^{\vphantom{\dagger}}\,b_{k}^{\dagger}b_{p}^{\vphantom{\dagger}}-G_{kp}^{*}\,d_{k}^{\vphantom{\dagger}}d_{p}^{\dagger}+H_{kp}^{\vphantom{\dagger}}\,b_{k\vphantom{p}}^{\dagger\vphantom{p}}d_{p\vphantom{k}}^{\dagger\vphantom{k}}-H_{kp}^{*}\,d_{k}^{\vphantom{p\dagger}}b_{p}^{\vphantom{k\dagger}}\,\big{)}\,,$
(39)

where the coefficients $G_{kp}$ and $H_{kp}$ are defined as:

$\displaystyle G_{kp}$ $\displaystyle\equiv
2\int_{-\infty}^{\infty}dt\,e^{i(E_{k}-E_{p})t}\,\gamma^{-1}(\dot{\eta}(t))\,g_{kp}(\eta(t))\,,$
(40) $\displaystyle H_{kp}$ $\displaystyle\equiv
2\int_{-\infty}^{\infty}dt\,e^{i(E_{k}+E_{p})t}\,\gamma^{-1}(\dot{\eta}(t))\,h_{kp}(\eta(t))\,.$
(41)

Note that $\Omega_{1}$ inherits the quadratic structure, in creation and
annihilation operators, of the interaction Hamiltonian in (23). In other
words, the same mononials in those operators appear. It is rather
straightforward to see that that kind of quadratic form also appears in all
the $\Omega_{n}$: they are also bilinear, and the same combinations of
operators appear (the only ones which preserve the fermion number). Of course
the kernels defining the quadratic form will in general be different for each
order.

The operator structure of (39) implies that, at this order, the only
possibility is the creation of a single particle pair, with the correspoding
probability given by:

$P_{kp}^{(1)}=|\braket{f}{\Omega_{1}}{i}|^{2}=|H_{kp}|^{2}\,.$ (42)

We note that this result coincides with what we would obtain by considering
the modulus squared of the first-order term in the expansion of the transition
amplitude $T_{fi}$ in (3).

Given the symmetry properties (26), the coefficients $G_{kp}$ are hermitian
and $H_{kp}$ are antisymmetric with respect to index swapping:

$G_{kp}={(G_{pk})}^{*}\,,\quad H_{kp}=-H_{pk}\,.$ (43)

In scenarios involving small amplitudes and non-relativistic speeds, i.e.,
$|\eta(t)|<<a$ and $|\dot{\eta}(t)|<<1$, these coefficients simplify to:

$\displaystyle G_{pk}$
$\displaystyle\simeq\frac{2}{a}\,(p+k)\,\tilde{\eta}(p-k)\,,$ (44)
$\displaystyle H_{pk}$
$\displaystyle\simeq\frac{2}{a}\,(p-k)\,\tilde{\eta}(p+k)\,,$ (45)

where $\tilde{\eta}(\nu)$ denotes the Fourier transform of $\eta(t)$, with the
convention: $\tilde{\eta}(\nu)\equiv\int dt\,e^{i\nu t}\,\eta(t)$.

As a particular example, for a periodic motion such as:

$\eta=\eta_{0}\,\sin(\omega t),\quad\eta_{0}>0,\,\omega>0\,,$ (46)

under conditions of small amplitude $\eta_{0}<<a$ and non-relativistic motion
$\omega\eta_{0}<<1$, the coefficients take the form:

$\displaystyle G_{pk}$
$\displaystyle\simeq\frac{2\pi}{i}\,\frac{\eta_{0}}{a}\,(p+k)\,\big{(}\delta(k-(p+\omega))-\delta(k-(p-\omega))\big{)}\,,$
(47) $\displaystyle H_{pk}$
$\displaystyle\simeq\frac{2\pi}{i}\,\frac{\eta_{0}}{a}\,(p-k)\,\delta(k-(\omega-p))\,.$
(48)

Thus, the first-order probability per unit time for pair creation is:

$\displaystyle\frac{P_{kp}^{(1)}}{T}\,=$
$\displaystyle\frac{2\pi}{a^{2}}\,\eta_{0}^{2}\,{(k-p)}^{2}\,\delta(k+p-\omega)\,.$
(49)

This result indicates that particle production is restricted to specific
frequencies of the boundary motion, and it only occurs when the oscillation
frequency $\omega$ matches the sum of the momenta of the particle pair.
Furthermore, due to the discrete nature of these values, particle production
is constrained to frequencies that satisfy the relation
$\omega=\frac{n\pi}{a}$, with $n=2,3,\ldots$.

Now, let’s consider scenarios with small amplitude oscillations where
velocities are not necessarily small. The Lorentz factor, can be expanded into
a Fourier series:

$\gamma^{-1}(\dot{\eta}(t))=\sqrt{1-\dot{\eta}^{2}(t)}=\sum_{l=0}^{\infty}a_{l}\,\cos(2l\,\omega
t)\,.$ (50)

The coefficients $a_{l}$ are defined as:

$a_{l}=\begin{cases}{}_{2}F_{1}\left(-\frac{1}{2},\frac{1}{2};1;{(\omega\eta_{0})}^{2}\right)&\text{for}\quad
l=0\,,\\\\[5.0pt]
\frac{2}{l!}\,{(-\frac{1}{2})}_{l}\,{}_{2}F_{1}\left(l-\frac{1}{2},l+\frac{1}{2};2l+1;{(\omega\eta_{0})}^{2}\right)&\text{for}\quad
l\geq 1\,,\end{cases}$ (51)

where ${(\cdot)}_{l}$ is the Pochhammer symbol, and ${}_{2}F_{1}$ is the
Gaussian hypergeometric function. For this setup, the pair-probability per
unit of time to the first order in the Magnus expansion, results:

$\displaystyle\frac{P_{kp}^{(1)}}{T}\,=$
$\displaystyle\frac{2\pi}{a^{2}}\,\eta_{0}^{2}\,{(k-p)}^{2}\,\big{[}\,a_{0}^{2}\,\delta(k+p-\omega)$
(52)
$\displaystyle+\frac{1}{4}\sum_{l=1}^{\infty}a_{l}^{2}\,\bigl{(}\delta(k+p-(2l+1)\,\omega)+\delta(k+p-(2l-1)\,\omega)\bigr{)}\,\big{]}\,.$

The coefficients $a_{l}$ represent the weight of the different contributions
to the creation of particles with higher momenta $k$ and $p$. Figure 1
displays these coefficients, showing $|a_{l}|$ as a function of the peak
velocity $\omega\eta_{0}$. The plot shows that $a_{0}\simeq 1$ and
$a_{l}\simeq 0$ for $l\geq 1$, under small velocities, aligning with the non-
relativistic outcome (49). As velocities increase, higher $a_{l}$ terms become
more significant, indicating that particles with $k+p\neq\omega$ can be
created.

Figure 1: Modulus of the Fourier expansion coefficients, $a_{l}$, as a
function of the mirror’s peak velocity, $\omega\eta_{0}$.

###  4.2 Second order

As it is evident in (37), the second order of the Dyson series, $U_{I}^{(2)}$,
receives contributions from both the first order term in the Magnus expansion,
$\Omega_{1}$, and from $\Omega_{2}$. Let us consider the first one, which
appears squared, and its action on the vacuum:

$\displaystyle{(\Omega_{1})}^{2}\,|0\rangle\;=$
$\displaystyle\Big{\\{}\sum_{k,l}\,{[(G-{\rm
tr}\,G)H]}_{kl}\,(b_{k}^{\dagger}\,d_{l}^{\dagger}+d_{k}^{\dagger}\,b_{l}^{\dagger})\,+\,{\big{(}\sum_{k,l}b_{k}^{\dagger}H_{kl}d_{l}^{\dagger}\big{)}}^{2}$
(53) $\displaystyle-{\rm tr}(H^{*}H)+{({\rm
tr}\,G)}^{2}\,\Big{\\}}\,|0\rangle\,,$

where ${\rm tr}$ denotes the trace.

Let us see how keeping the exponential form for the $S$ matrix resulting from
the Magnus expansion yields the multipair transition amplitudes, a property
which, in fact, holds true as a consequence of the quadratic structure of the
$\Omega_{n}$ operators-

The key is to identify each contribution to the multipair transition
amplitudes by analyzing the operator structure applied to the vacuum in (53).
The first term, with its quadratic structure, contributes to the creation of
one pair, serving as a higher-order correction to the result in (42). The
second term, characterized by a quartic structure, leads to the creation of
two pairs. Lastly, the term proportional to the identity accounts for the
vacuum-to-vacuum transition amplitude.

On the other and, the second-order term in the Magnus expansion, $\Omega_{2}$,
contributes just to the probability amplitude of creating one pair. Its
contribution is given by:

$\displaystyle\braket{f}{\Omega_{2}}{i}=$ $\displaystyle
2\sum_{l}\int_{0}^{\infty}dt_{1}\,\big{[}e^{i(E_{k}+E_{l})t_{1}}\,\left(g_{lp}*h_{kl}\right)(t_{1})-e^{i(E_{k}-E_{l})t_{1}}\,\left(h_{lp}*g_{kl}\right)(t_{1})$
(54)
$\displaystyle+e^{i(E_{p}+E_{l})t_{1}}\,\left(g_{kl}*h_{lp}\right)(t_{1})-e^{i(E_{p}-E_{l})t_{1}}\,\left(h_{kl}*g_{lp}\right)(t_{1})\big{]}\,,$
(55)

with the definitions:

$\displaystyle\begin{aligned}
\left(g_{lp}*h_{kl}\right)(t_{1})\,\equiv\,\int_{-\infty}^{\infty}dt_{2}\,e^{i(E_{k}+E_{p})t_{2}}\,&\gamma^{-1}(\dot{\eta}(t_{2}))\,\gamma^{-1}(\dot{\eta}(t_{1}+t_{2}))\\\
\,&g_{lp}(t_{2})\,h_{kl}(t_{1}+t_{2})\,,\end{aligned}$ (56)
$\displaystyle\begin{aligned}
\left(g_{kl}*h_{lp}\right)(t_{1})\,\equiv\,\int_{-\infty}^{\infty}dt_{2}\,e^{i(E_{k}+E_{p})t_{2}}\,&\gamma^{-1}(\dot{\eta}(t_{2}))\,\,\gamma^{-1}(\dot{\eta}(t_{1}+t_{2}))\\\
\,&g_{kl}(t_{2})\,h_{lp}(t_{1}+t_{2})\,.\end{aligned}$ (57)

As we did with the first order term, we can trace the origin of these
contributions to the one pair creation amplitude back to the terms in the
expansion of the transition amplitude, $T_{fi}$. In this case, these
contributions originate from the second term of the expansion of
$M(x^{1},y^{1})$ in (3).

##  5 Bogoliubov transformation

In this section, we examine how the dynamic boundary conditions affect vacuum
states and the creation and annihilation operators of the fermionic field
within the cavity. We specifically focus on how the in and out operators are
related, emphasizing the transformations induced by changes in boundary
conditions.

The in vacuum, $|0\rangle_{\text{in}}$, is the lowest energy state defined
before the boundaries begin to oscillate. It serves as the foundation for the
in Fock space, which is constructed by the successive application of in
particle and antiparticle creation operators to this vacuum state. Similarly,
the out vacuum, $|0\rangle_{\text{out}}$, is the lowest energy state after the
boundary motion ceases, and serves as the basis for constructing the out Fock
space through similar applications of out creation operators. These two vacuum
states are not equivalent, which is manifested in the possibility of the in
vacuum evolving into a state containing particles in the out vacuum.

The in and out operators are connected through the Bogoliubov transformation
$U$:

$\begin{pmatrix}b\\\
d^{\dagger}\end{pmatrix}_{\text{out}}=\,U\,{\begin{pmatrix}b\\\
d^{\dagger}\end{pmatrix}}_{\text{in}},$ (58)

where $U$ is a unitary matrix that preserves the anticommutation relations
between operators in (15), for both in and out operators.

To explicitly express the Bogoliubov transformation, we organize the particle
creation and annihilation operators into vectors. The components of these
vectors are the operators arranged in ascending order of their momenta, as
follows:

${(b,d^{\dagger})}^{T}={(b_{p_{0}},b_{p_{1}},\ldots,d_{p_{0}}^{\dagger},d_{p_{1}}^{\dagger},\ldots)}^{T}\,.$
(59)

Due to the confined nature of the system, these operators do not correspond to
well-defined momentum states. Instead, they are linked to quantized energy
states that arise from the perfect reflective ‘bag’ boundary conditions
specified in equation (4). Moreover, we emphasize that the fermions are
strictly confined within the cavity due to these boundary conditions, and
there are no fermions outside of it.

To derive the transformation $U$, we examine a general operator $A$, for which
we know that the $A_{\text{in}}$ and $A_{\text{out}}$ operators are related by
a canonical transformation represented by the unitary operator $U_{I}$:

$A_{\text{out}}=U_{I}^{\dagger}\,A_{\text{in}}\,U_{I}\,.$ (60)

At first order in the Magnus expansion, we have:

$A_{\text{out}}\,\simeq\,e^{i\Omega_{1}}\,A_{\text{in}}\,e^{-i\Omega_{1}}\,.$
(61)

Employing Hadamard’s formula, we express the right-hand side of (61) as an
infinite series of nested commutators of operators $\Omega_{1}$ and
$A_{\text{in}}$:

$A_{\text{out}}\;\simeq\;\sum_{n=0}^{\infty}\frac{i^{n}}{n!}\underbrace{[\Omega_{1},\ldots[\Omega_{1},[\Omega_{1}}_{\text{n
times}},A_{\text{in}}]]\ldots]\,.$ (62)

We then introduce a parameter-dependent operator $A(s)$, where
$A(1)=A_{\text{out}}$ and $A(0)=A_{\text{in}}$. Expanding $A(s)$ as a Taylor
series around $s=0$, yields:

$A(s)=\sum_{n=0}^{\infty}\frac{s^{n}}{n!}\,A^{(n)}(0)\,,$ (63)

where $A^{(n)}$ represents the nth derivative with respect to $s$.
Substituting in (63) the successive derivatives of the differential equation:

$i\frac{d}{ds}A(s)=[A(s),\Omega_{1}]\,,$ (64)

and evaluating at $s=1$, recovers (62).

Extending this approach to our basis operators,
${(b_{k}^{\vphantom{\dagger}}(s),d_{k}^{\dagger}(s))}^{T}$, leads to a
differential equation for the operator vector:

$i\frac{d}{ds}\begin{pmatrix}b_{k}(s)\\\
d_{k}^{\dagger}(s)\end{pmatrix}\;\simeq\;\sum_{p}M_{kp}\begin{pmatrix}b_{p}(s)\\\
d_{p}^{\dagger}(s)\end{pmatrix},$ (65)

where $p$ is summed over the discrete spectrum that results from solving (12),
and the operator mixing matrix $M$ is:

$M_{kp}\equiv\begin{pmatrix}G_{kp}&H_{kp}\\\
-H_{kp}^{*}&-G_{kp}^{*}\end{pmatrix}.$ (66)

Finally, proceding as we did with $A$ in (63), we obtain the Bogoliubov
transformation to the first order in the Magnus expansion for the basis of
operators:

$\begin{pmatrix}b\\\
d^{\dagger}\end{pmatrix}_{\text{out}}\,\simeq\,e^{-iM}\,{\begin{pmatrix}b\\\
d^{\dagger}\end{pmatrix}}_{\text{in}}.$ (67)

We observe that, to this order in the Magnus expansion, particle pair
production requieres the mixing of the $b$ and $d^{\dagger}$ operators;
therefore, the mixing matrix $M$ must be non-diagonal to enable pair creation.
Consequently, to produce a particle and anti-particle, the matrix element
$H_{kp}$ must be non-zero for some values of $k$ and $p$.

Now, we will calculate the Bogoliubov transformation $U$ for the non-
relativistic oscillatory movement previously discussed, defined by equation
(46). For this motion, the values of $G_{kp}$ and $H_{kp}$ are provided in
equations (47) and (48). To compute the transformation, we must evaluate the
exponential of the $M$ matrix, (66). Given that the momenta spectrum is
unbounded, the $M$ matrix is infinite dimensional. To derive an explicit
expression for $U$, we introduce a cutoff $p_{N}$ for the state index. This
results in the vectors ${(b,d^{\dagger})}^{T}$ in (67) becoming finite
dimensional:

${(b,d^{\dagger})}^{T}={(b_{p_{0}}^{\vphantom{\dagger}},b_{p_{1}}^{\vphantom{\dagger}},\ldots,b_{p_{N}}^{\vphantom{\dagger}},d_{p_{0}}^{\dagger},d_{p_{1}}^{\dagger},\ldots,d_{p_{N}}^{\dagger})}^{T}\,,$
(68)

and consequently, leading to a finite dimensional Bogoliubov transformation
matrix.

We note that applying a momenta cutoff is equivalent to considering a theory
with a finite set of fermionic degrees of freedom. This is reflected in the
Fourier representation of the interaction Hamiltonian, (23), where the sum
over the mixing of particle and anti-particle creation and destruction
operators becomes finite.

We conclude by presenting results for the Bogoliubov transformation $U$ at the
threshold frequency $\omega=\frac{2\pi}{a}$, and for the two lowest non-
trivial cutoffs.

For a cutoff $p_{N}=p_{1}=\frac{\pi}{2a}$, we have:

$U\,=\,\left(\begin{array}[]{cccc}\cos(\mathcal{N})&0&0&\sin(\mathcal{N})\\\
0&\cos(\mathcal{N})&-\sin(\mathcal{N})&0\\\
0&\sin(\mathcal{N})&\cos(\mathcal{N})&0\\\
-\sin(\mathcal{N})&0&0&\cos(\mathcal{N})\\\ \end{array}\right)\,,$ (69)

where $\mathcal{N}\equiv\frac{\pi\eta_{0}T}{a^{2}}$, and $T$ is the total
time.

For $p_{N}=p_{2}=\frac{3\pi}{2a}$, we have:

$U\,=\,\left(\begin{array}[]{cccccc}c&0&\frac{3s}{\sqrt{10}}&0&\frac{s}{\sqrt{10}}&0\\\
0&\frac{c+9}{10}&0&-\frac{s}{\sqrt{10}}&0&\frac{3(c-1)}{10}\\\
-\frac{3s}{\sqrt{10}}&0&\frac{9c+1}{10}&0&\frac{3(c-1)}{10}&0\\\
0&\frac{s}{\sqrt{10}}&0&c&0&\frac{3s}{\sqrt{10}}\\\
-\frac{s}{\sqrt{10}}&0&\frac{3(c-1)}{10}&0&\frac{c+9}{10}&0\\\
0&\frac{3(c-1)}{10}&0&-\frac{3s}{\sqrt{10}}&0&\frac{9c+1}{10}\\\
\end{array}\right)\,,$ (70)

where $s\equiv\sin\left(\sqrt{10}\,\mathcal{N}\right)$, and
$c\equiv\cos\left(\sqrt{10}\,\mathcal{N}\right)$.

##  6 Conclusions

In this work, we analyzed the dynamical Casimir effect for fermions confined
by moving boundaries in a $1+1$ dimensional setup, with a particular focus on
the role of higher-order perturbative terms and the Bogoliubov transformation
between the creation and annihilation operators. Our methodology incorporated
both the Dyson series and the Magnus expansion to compute the probability of
fermion pair creation due to boundary oscillations. A notable advantage of the
Magnus expansion is its capability to ensure a unitary evolution of the system
at all orders.

We derived formulas for the first and second orders in the Magnus expansion,
specifically calculating the pair creation probabilities at the first order
for both relativistic and non-relativistic periodic movements. These
calculations revealed the dependence of pair creation probabilities on the
speed of the boundary’s movements. This analysis demonstrates how the velocity
of boundary movement influences the various contributions to the probability
of pair creation. Furthermore, we examined how each order contributes not only
to the single-pair creation amplitude but also to multipair creation
processes, along with corrections to the single-pair amplitude.

Additionally, we explicitly computed the Bogoliubov transformation for the
case of non-relativistic periodic movement. This computation further
elucidated the mechanism by which the in operators are transformed into a
mixture of out operators due to the dynamic boundary conditions. This shows
how quantum states evolve in response to changes in the geometry and dynamics
of the boundary.

## Acknowledgements

The authors thank ANPCyT, CONICET and UNCuyo for financial support.

## References

  * [1] H. B. G. Casimir, Indagat. Math. 10, 261 (1948); Kon. Ned. Akad. Wetensch. Proc. 51, 793 (1948); Front. Phys. 65, 342 (1987); Kon. Ned. Akad. Wetensch. Proc. 100, 61 (1997).
  * [2] For a review, see, for example: M. Bordag, U. Mohideen and V. M. Mostepanenko, Phys. Rept. 353, 1 (2001).
  * [3] For a review, see, for example: V. Dodonov, MDPI Phys. 2, 67 (2020).
  * [4] A. Chodos, R. L. Jaffe, K. Johnson and C. B. Thorn, Phys. Rev. D 10, 2599 (1974).
  * [5] E. Fradkin, Field Theories of Condensed Matter Systems (Cambridge University Press, New York, 2013).
  * [6] F. D. Mazzitelli, J. P. Paz and M. A. Castagnino, Phys. Lett. B 189, 132 (1987).
  * [7] C. D. Fosco and G. Hansen, Phys. Rev. D 105, 016004 (2022).
  * [8] C. D. Fosco and E. Losada, Phys. Rev. D 78, 025017 (2008).
  * [9] C. C. Ttira, C. D. Fosco, and E. Losada, Phys. Rev. D 82, 085008 (2010).
  * [10] C. D. Fosco and G. Hansen, Phys. Rev. D 108, 056005 (2023).
  * [11] C. Itzykson and J. B. Zuber, Quantum Field Theory (McGraw-Hill, New York, 1980).
  * [12] S. Blanes, F. Casas, J. A. Oteo, and J. Ros, Eur. J. Phys. 31, 907 (2010).