# Acceleration of a polarized neutron by its weak nuclear self-interaction

M. Donaire<EMAIL_ADDRESS>Departamento de Física Teórica, Atómica y
Óptica, Universidad de Valladolid, 47011 Valladolid, Spain  CNRS, Université
Grenoble Alpes, Institut Néel, 38000, Grenoble, France

###### Abstract

It is proven that the rotation of the spin of a polarized neutron is
accompanied by a net nuclear force upon it. This force arises from the weak
nuclear self-interaction of its constituent quarks, whose chiral nature
induces the transfer of a net momentum to the fields of Z and W-bosons. This
effect is linear in Fermi’s constant. As a result, it is estimated that along
the spin-flip of a polarized neutron its velocity undergoes a variation of the
order of meters per second.

##  I Introduction

The weak nuclear interaction of the Standard Model is chiral [1, 2, 3]. This
implies that its Hamiltonian is not parity invariant, which manifests in the
preferential coupling of left-handed particles to right-handed antiparticles.
In turn, that results in a net alignment between the spins and the momenta of
the hadrons and the leptons which participate in the interaction. It explains,
for instance, the celebrated experiment of Wu _et al._ [4] on the beta decay
of polarized Co60 nuclei, where left-handed electrons are preferably emitted
in the direction opposite to the spins of the nuclei [5]. Parity violation in
the electroweak (EW) interaction can be also tested through kinematical
effects in the scattering of leptons with hadrons. This is the case of
polarized electron scattering [6], where the EW parameters and the hadron
structure are probed through the small asymmetry in the scattering cross-
section of electrons with opposite helicities. Along these lines, the aims of
this article are twofold. From a fundamental perspective we will prove that
the weak nuclear self-interaction of a polarized neutron may generate a net
force upon it, transferring an equivalent linear momentum to virtual bosons.
From a practical perspective we will show that this force has macroscopic
kinematical effects, since it is linear in $G_{F}$, causing a variation on the
neutron velocity of the order of meters per second. Incidentally, the
dependence of this estimate on the hadron structure makes it suitable for
probing hadronic features, like the spin structure of nucleons, with low-
energy experimental techniques.

Homogeneity of space and translation invariance of the vacuum state of any
quantum field theory imply that the net momentum of vacuum field fluctuations
vanishes identically. However, when field fluctuations couple to current
fluctuations through a Hamiltonian which violates both parity and time-
reversal, it is symmetry allowed for the field fluctuations to carry a net
linear momentum. In virtue of total momentum conservation, the currents must
gain a kinetic momentum of equal strength in the opposite direction. The
possibility that a net momentum could be carried by vacuum field fluctuations
was first envisaged by A. Feigel in a magnetoelectric medium [7, 8]. Using a
semiclassical approach he showed that, in crossed electric and magnetic
fields, $\mathbf{E}\perp\mathbf{B}$, the medium generates a magnetoelectric
birefringence along $\mathbf{E}\times\mathbf{B}$. As a result, the momenta
carried by the modes of the quantum electromagnetic (EM) field propagating
parallel and antiparallel to $\mathbf{E}\times\mathbf{B}$ do not cancel each
other and, in virtue of total momentum conservation, a momentum equivalent to
their difference is transferred to the magnetoelectric medium in the opposite
direction. Later, B. van-Tiggelen and collaborators, using a fully QED
approach, have computed the momentum gained by a hydrogen atom in crossed
electric and magnetic fields [9], and the momentum of a chiral molecule in a
magnetic field [10]. In the former case parity and time-reversal are broken by
the electric and the magnetic fields, respectively, while in the latter parity
is broken by the chiral distribution of the ligands around the active ion of
the chiral molecule. It was found in Ref.[11] that, in its major part, the
acceleration of the chiral molecule is accompanied by a net EM momentum which
arises from the crossed electric and magnetic fields sourced by the charges
and the helical currents within the chiral molecule. In either case, the
kinetic momenta gained by atoms and molecules were found to be too small to be
experimentally observed.

Further, it has been proved in Ref.[12] that an analogous phenomenon may take
place in nuclear physics, where the EW Hamiltonian itself is not parity
invariant. Using the Fermi four-fermion interaction, it was shown that an
unpolarized proton in a magnetic field accelerates along the magnetic field,
transferring an equivalent momentum to virtual positrons and neutrinos. The
explanation given in Ref.[12] is analogous to that which explains the
directionality of the beta decay, but for the fact that the leptons in
Ref.[12] are virtual instead of actual particles. The velocity of the proton
computed that way happens to be of order $G_{F}^{2}$, thus too small to be
detected. The reason for this was the inappropriate use of the Fermi
interaction in a process that, in contrast to the beta decay, is not resonant
and does not produce actual leptons. Hence, the estimate in Ref.[12] is not
the leading order one but a second order correction.

In this article we use the gauge theory of EW interactions and prove that,
throughout the spin-flip of a polarised neutron, a net internal force
accelerates it. This force arises from the weak nuclear self-interaction of
the constituent or valence quarks. In virtue of total momentum conservation,
an equivalent linear momentum is transferred to virtual Z and W-bosons. Using
two simple models of hadron confinement, namely, the MIT-bag model and the
harmonic oscillator model, we estimate a velocity variation of the order of
meters per second along the spin direction. If confirmed experimentally, this
finding would prove the existence of net self-interaction forces and the
transfer of momentum to virtual particles, and would open up the possibility
of probing the hadron structure with low-energy techniques.

##  II Fundamentals

For simplicity and without much loss of generality, let us consider the
Lagrangian of a neutron as the sum of the Lagrangians of its three valence or
constituent quarks, two down and one up, which are bound by some effective QCD
confinement potential, $\mathcal{V}_{QCD}$, and are minimally coupled to the
EW bosons,

$\displaystyle L=$ $\displaystyle\sum_{j=1}^{3}\int
d^{3}r_{j}\Bigl{[}\bar{\Psi}_{j}\bigl{(}i\gamma^{\mu}D^{j}_{\mu}-m_{j}\mathbb{I}\bigr{)}\Psi_{j}-\sum_{k<j}^{3}\int
d^{3}r_{k}$ $\displaystyle\times$
$\displaystyle\Psi^{\dagger}_{j}\Psi_{k}\mathcal{V}_{QCD}(\mathbf{r}_{j}-\mathbf{r}_{k})\Psi^{\dagger}_{k}\Psi_{j}\Bigr{]}.$
(1)

In this equation, given in natural units, $\Psi$ denotes a generic isospin
doublet, $\Psi_{j}=\bigl{(}\psi^{u}_{j}(t,\mathbf{r}_{j})$
$\psi^{d}_{j}(t,\mathbf{r}_{j})\bigr{)}$, and $D^{j}_{\mu}$ is the covariant
derivative operator upon quark $j$. In terms of the EW gauge fields it reads,

$\displaystyle D^{j}_{\mu}$
$\displaystyle=\partial_{\mu}^{j}\mathbb{I}_{\textrm{iso}}-ie\>A_{\mu}(\mathbf{r}_{j})\mathbb{Q}-\frac{igZ_{\mu}(\mathbf{r}_{j})}{2\cos{\theta_{W}}}\bigl{[}\mathbb{I}\bigl{(}\tau_{3}-2\sin^{2}{\theta_{W}}\mathbb{Q}\bigr{)}$
$\displaystyle-\gamma_{5}\tau_{3}\bigr{]}-\frac{ig}{\sqrt{2}}\bigl{[}W^{+}_{\mu}(\mathbf{r}_{j})\frac{\mathbb{I}-\gamma_{5}}{2}(\tau_{1}+i\tau_{2})+\textrm{h.c.}\bigr{]},$
(2)

where $\mathbb{I}_{\textrm{iso}}$, $\mathbb{Q}$ and $\tau_{i}$ are two-by-two
matrices acting on the isospin indices, namely, the identity, the electric
charge operator and the isospin generators, respectively; $\mathbb{I}$ and
$\gamma_{5}$ act on the four-spinor components of $\psi_{j}^{u,d}$; and $e$,
$g/(2\cos{\theta_{W}})$, $g/\sqrt{2}$ are the coupling constants of A, Z and W
gauge fields, respectively, with $\theta_{W}$ the weak mixing angle.

As for the confinement potential, without much loss of generality and for
numerical estimations we will consider two simple models. These are, the MIT-
bag model for valence and ultrarelativistic quarks [13], and the harmonic-
oscillator model for constituent and nonrelativistic quarks [14] –see appendix
A for further details. In the MIT-bag the three quarks are massless and are
confined within a static spherical cavity of radius $R_{n}$ by means of
boundary conditions at the surface. These conditions prevent the current flow
through the surface and impose the balance between the outwards pressure of
the quarks and the inward QCD vacuum pressure. As for the nonrelativistic
harmonic oscillator, all the quarks have an equivalent mass $m$ and interact
pairwise through harmonic potentials of equivalent natural frequencies.

In the remainder of this article we will work in the Hamiltonian formalism in
Schrödinger’s picture. Thus, we identify the non-perturbative Hamiltonians
with those of the neutron quarks and the free boson fields. As for the former,
it can be derived from Eq.(1),

$\displaystyle H_{n}=$ $\displaystyle\sum_{j=1}^{3}\int
d^{3}r_{j}\Bigl{[}\sum_{f=u,d}\bar{\psi}^{f}_{j}\bigl{(}-i\gamma^{m}\partial^{j}_{m}+m_{j}\bigr{)}\psi^{f}_{j}+\sum_{k<j}^{3}\int
d^{3}r_{k}$ $\displaystyle\times$
$\displaystyle\sum_{\tilde{f}=u,d}\psi^{f\dagger}_{j}\psi^{\tilde{f}}_{k}\mathcal{V}_{QCD}(\mathbf{r}_{j}-\mathbf{r}_{k})\psi^{\tilde{f}\dagger}_{k}\psi^{f}_{j}\Bigr{]}.$
(3)

As for the free EW bosons, the gauge fields $A_{\mu}$, $Z_{\mu}$ and
$W^{\pm}_{\mu}$ that we will denote generically by $V_{\mu}$ can be decomposed
in sums over normal modes of wave vector $\mathbf{q}$, polarization $\lambda$,
creation/annihilation operators $a_{V,\mathbf{q}}^{\lambda(,\dagger)}$, and
energy $\mathcal{E}^{V}_{\mathbf{q}}=\sqrt{q^{2}+M^{2}_{V}}$, with $M_{V}$
being the mass [3]. Hence, the Hamiltonian of free bosons reads

$H_{F}^{V}=\sum_{\mathbf{q},\lambda}\mathcal{E}^{V}_{\mathbf{q}}a_{V,\mathbf{q}}^{\lambda,\dagger}a_{V,\mathbf{q}}^{\lambda},\quad
V=\\{A,Z,W^{\pm}\\}.$ (4)

The Hamiltonian densities of the interactions between quarks and EW bosons can
be read from the minimal coupling term of Eq.(1) together with Eq.(2),

$\displaystyle\mathcal{W}_{EM}$
$\displaystyle=e\sum_{i=1}^{3}j_{EM}^{i,\mu}A_{\mu}(\mathbf{r}_{i}),\>\>\mathcal{W}_{nc}=\frac{g}{2\cos{\theta_{W}}}\sum_{i=1}^{3}j_{nc}^{i,\mu}Z_{\mu}(\mathbf{r}_{i}),$
$\displaystyle\mathcal{W}_{cc}$
$\displaystyle=\frac{g}{\sqrt{2}}\sum_{i=1}^{3}j_{cc}^{i,\mu}W^{+}_{\mu}(\mathbf{r}_{i})+\textrm{h.c.}$
(5)

In these equations $j_{EM}^{i,\mu}$ is the neutral and parity-invariant EM
current density of quark $i$,
$j_{EM}^{i,\mu}=q_{u}\bar{\psi}^{u}_{i}\gamma^{\mu}\psi^{u}_{i}+q_{d}\bar{\psi}_{i}^{d}\gamma^{\mu}\psi^{d}_{i}$,
where $q_{u,d}$ are the electric charges;
$j_{nc}^{i,\mu}=\bar{\psi}^{d}_{i}\gamma^{\mu}(g_{v}^{d}\mathbb{I}-g_{a}^{d}\gamma_{5})\psi^{d}_{i}+\bar{\psi}^{u}_{i}\gamma^{\mu}(g_{v}^{u}\mathbb{I}-g_{a}^{u}\gamma_{5})\psi^{u}_{i}$
is the neutral and chiral current, with $g_{v,a}^{u,d}$ being the vector $(v)$
and axial $(a)$ couplings; and
$j_{cc}^{i,\mu}=\frac{1}{2}\bar{\psi}^{d}_{i}\gamma^{\mu}(\mathbb{I}-\gamma_{5})\psi^{u}_{i}$
is the charged and chiral current.

The system of hadrons and EW bosons governed by the Hamiltonians of
Eqs.(3)-(5) possesses a conserved linear momentum $\mathbf{K}$ which comprises
the canonical conjugate momentum of the three quarks bound in the neutron and
the _transverse_ momenta of the EW bosons,
$\mathbf{K}=\mathbf{P}+\mathbf{P}^{A}_{\perp}+\mathbf{P}^{Z}_{\perp}+\mathbf{P}^{W}_{\perp}$.
The transverse momentum operator of each boson V reads

$\mathbf{P}_{\perp}^{V}=\sum_{\mathbf{q},\lambda}\mathbf{q}a_{V,\mathbf{q}}^{\lambda,\dagger}a_{V,\mathbf{q}}^{\lambda},\quad
V=\\{A,Z,W^{\pm}\\},$ (6)

whereas the total conjugate momentum $\mathbf{P}$ is made of a kinetic
momentum $\mathbf{P}_{\textrm{kin}}$ and the _longitudinal_ momenta associated
to the EW bosons,
$\mathbf{P}=\mathbf{P}_{\textrm{kin}}+\mathbf{P}^{A}_{\parallel}+\mathbf{P}^{Z}_{\parallel}+\mathbf{P}^{W}_{\parallel}$
[15]. The expressions for the longitudinal momenta in terms of boson fields
can be deduced from the spatial components of the covariant derivative in
Eq.(2). It suffices to realize that, in the minimal coupling formalism, the
total conjugate and kinetic momentum operators read
$\mathbf{P}=-i\sum_{j=1}^{3}\boldsymbol{\nabla}^{j}$,
$\mathbf{P}_{\textrm{kin}}=-i\sum_{j=1}^{3}\mathbf{D}^{j}$, from which,

$\mathbf{P}^{A}_{\parallel}+\mathbf{P}^{Z}_{\parallel}+\mathbf{P}^{W}_{\parallel}=-i\sum_{j=1}^{3}(\boldsymbol{\nabla}^{j}-\mathbf{D}^{j}).$
(7)

Therefore, conservation of total momentum $\mathbf{K}$ implies

$\Delta\langle\mathbf{P}_{\textrm{kin}}\rangle=-\Delta\langle\mathbf{P}^{A}_{\parallel}+\mathbf{P}^{Z}_{\parallel}+\mathbf{P}^{W}_{\parallel}+\mathbf{P}^{A}_{\perp}+\mathbf{P}^{Z}_{\perp}+\mathbf{P}^{W}_{\perp}\rangle.$
(8)

For future purposes we write the densities of Z and W longitudinal momenta in
terms of the quarks’ wavefunctions [16],

$\boldsymbol{\mathcal{P}}_{\parallel}^{Z}=\sum_{i=1}^{3}\frac{g}{2\cos{\theta_{W}}}j_{nc}^{i0}\mathbf{Z}(\mathbf{r}_{i}),\>\>\boldsymbol{\mathcal{P}}_{\parallel}^{W}=\sum_{i=1}^{3}\frac{g}{\sqrt{2}}j_{cc}^{i0}\mathbf{W}^{+}(\mathbf{r}_{i})+\textrm{h.c.}$
(9)

##  III Calculation

Let us consider a free neutron initially polarized along the
$-\hat{\mathbf{z}}$ direction in the presence of an intense constant magnetic
field $\mathbf{B}_{0}=B_{0}\hat{\mathbf{z}}$. The wavefunction of the neutron,
made of three quarks, must be totally antisymmetric. For simplicity, let us
consider that the neutron spin results entirely from the addition of the spins
of its quarks. The color function is an antisymmetric SU(3) singlet,
$\xi_{color}$. In its ground state, the spatial wavefunction of each quark is
an $l=0$ function, and the spin-flavour wavefunction must be symmetrized
($\mathcal{S}$). All together, the ground state of a spin-down neutron can be
written as [1]

$|n^{0}_{\downarrow}\rangle=\frac{\xi_{color}}{\sqrt{3}}\left[|d\downarrow
d\uparrow u\downarrow\rangle^{0}_{\mathcal{S}}-2|d\downarrow d\downarrow
u\uparrow\rangle^{0}_{\mathcal{S}}\right],$ (10)

and an equivalent expression holds for $|n^{0}_{\uparrow}\rangle$ but for the
inversion of all the quarks’ spins –cf. appendix B.

We aim at proving that, as the neutron spin is inverted, a variation in the
neutron kinetic momentum is induced by its weak nuclear self-interaction. A
simple manner to reverse the spin consists of applying an additional magnetic
field $\mathbf{B}_{1}$ much weaker than $\mathbf{B}_{0}$ which rotates at the
Larmor frequency $\omega_{0}=\gamma_{n}B_{0}$ in the $xy$ plane, with
$\gamma_{n}$ being the neutron gyromagnetic ratio and $\omega_{0}$ being
typically a radiofrequency. This field makes the neutron spin oscillate at
frequency $\omega_{1}=\gamma_{n}B_{1}$, inverting its direction in a time
interval $\pi/\omega_{1}$ and absorbing one photon of frequency $\omega_{0}$
–see appendix B.

Next, we compute the kinetic momentum gained by the neutron throughout its
spin-flip (sf), $\Delta\langle\mathbf{P}^{\textrm{sf}}_{\textrm{kin}}\rangle$.
Since the EM interaction is parity invariant, the transfer of EM momentum is
only caused by the absorption of an actual photon of frequency $\omega_{0}$,
which is negligibly small for radiofrequencies. In contrast, the weak nuclear
interaction is chiral and its mediating Z and W-bosons are only virtual, so in
order for them to carry a net momentum it is the weak self-interaction of the
quarks that must provide it. We will find that this is indeed the case along
the direction of the quarks’ spins, and that it is the dominant contribution
to Eq.(8). Hence, we may write

$\displaystyle\Delta\langle\mathbf{P}^{\textrm{sf}}_{\textrm{kin}}\rangle$
$\displaystyle\simeq\langle
n^{0}_{\downarrow}|\mathbf{P}^{Z}_{\parallel}+\mathbf{P}^{Z}_{\perp}+\mathbf{P}^{W}_{\parallel}+\mathbf{P}^{W}_{\perp}|n^{0}_{\downarrow}\rangle$
$\displaystyle-\langle
n^{0}_{\uparrow}|\mathbf{P}^{Z}_{\parallel}+\mathbf{P}^{Z}_{\perp}+\mathbf{P}^{W}_{\parallel}+\mathbf{P}^{W}_{\perp}|n^{0}_{\uparrow}\rangle.$
(11)

Finally, symmetry considerations imply that the two terms on the right-hand-
side of Eq.(11) have equivalent strengths and opposite signs. Hence, it
suffices to compute any of them.

Figure 1: Diagrammatic representation of two of the processes contributing to
Eq.(12) for $\langle
n^{0}_{\uparrow}|\mathbf{P}^{\textrm{W}}_{\parallel}|n^{0}_{\uparrow}\rangle$
[(a)] and to Eq.(13) for $\langle
n^{0}_{\uparrow}|\mathbf{P}^{\textrm{Z}}_{\perp}|n^{0}_{\uparrow}\rangle$
[(b)]. Time runs along the vertical from $t_{0}$ up to the time of observation
$t$ at which the momentum operators apply. Interaction vertices enter at the
intermediate time $\tau$.

The calculation of Eq.(11) can be performed adiabatically, since the spin
oscillation frequency, $\omega_{1}$, is much smaller than the transition
frequencies and decay rates of the quantum intermediate processes. These
frequencies are the masses $M_{Z,W}$ and decay rates $\Gamma_{W,Z}$ of the
bosons, satisfying $M_{Z,W}\gg\Gamma_{W,Z}\gg\omega_{1}$ –cf. appendix C for
further explanations. Using time-dependent perturbation theory, the
longitudinal and transverse momenta of Eq.(11) are computed at first and
second order in $\mathcal{W}_{nc,cc}$, respectively. The diagrams in Fig.1
depict two of the processes contributing to $\langle
n^{0}_{\uparrow}|\mathbf{P}^{\textrm{W}}_{\parallel}|n^{0}_{\uparrow}\rangle$
[(a)] and $\langle
n^{0}_{\uparrow}|\mathbf{P}^{\textrm{Z}}_{\perp}|n^{0}_{\uparrow}\rangle$
[(b)]. In all the contributing intermediate processes the spins of all the
quarks are preserved by the interaction with Z and W-bosons. Hence, the
intermediate states $\\{|I\rangle\\}$ are spin-up protons
$|p^{I}_{\uparrow}\rangle$ for
$\langle\mathbf{P}^{W}_{\parallel,\perp}\rangle$ and spin-up neutrons
$|n^{I}_{\uparrow}\rangle$ for
$\langle\mathbf{P}^{Z}_{\parallel,\perp}\rangle$.

Writing generically $V_{\mu}$ in the place of $Z_{\mu}$ and $W_{\mu}$ fields,
and keeping nonvanishing terms only, we obtain for the longitudinal momenta,

$\displaystyle\langle
n^{0}_{\uparrow}|\mathbf{P}^{\textrm{V}}_{\parallel}|n^{0}_{\uparrow}\rangle=-2\textrm{Re}\sum_{I}\int
d^{3}x\int
d^{3}y\left\langle\boldsymbol{\mathcal{P}}_{\parallel,I}^{V\dagger}(\mathbf{x})\frac{\mathcal{W}_{I}(\mathbf{y})}{H^{V}_{F}}\right\rangle$
$\displaystyle=-2g^{2}_{V}\textrm{Re}\sum_{I}\int d^{3}x\int
d^{3}y\left\langle
j_{I}^{0}(\mathbf{x})\mathbf{V}^{\dagger}(\mathbf{x})\frac{\mathbf{V}(\mathbf{y})\cdot\mathbf{j}_{I}(\mathbf{y})}{H^{V}_{F}}\right\rangle,$
(12)

where $g_{V}$ is the coupling constant of each field, and
$\boldsymbol{\mathcal{P}}_{\parallel,I}^{V}(\mathbf{x})$ and
$\mathcal{W}_{I}(\mathbf{y})$ are the longitudinal momentum and interaction
potential densities, respectively, evaluated between
$|n^{0}_{\uparrow}\rangle$ and the intermediate nucleon state $|I\rangle$.
They have been rewritten after the second equality in terms of the
$0^{th}$-component of the chiral current $j_{I}^{0}$–i.e., the weak charge
density, and the transient chiral current $\mathbf{j}_{I}(\mathbf{y})$ which
goes along the neutron spin. The big brackets stand for the expectation value
of the boson field operators in the vacuum. As for the transverse momenta,
their net contribution is

$\displaystyle\langle
n^{0}_{\uparrow}|\mathbf{P}^{\textrm{V}}_{\perp}|n^{0}_{\uparrow}\rangle=\sum_{I}\int
d^{3}x\int
d^{3}y\left\langle\frac{\mathcal{W}^{\dagger}_{I}(\mathbf{x})}{H^{V}_{F}}\mathbf{P}_{\perp}^{V}\frac{\mathcal{W}_{I}(\mathbf{y})}{H^{V}_{F}}\right\rangle$
$\displaystyle=2g^{2}_{V}\textrm{Re}\sum_{I}\int d^{3}x\int
d^{3}y\left\langle\frac{j_{I}^{0}(\mathbf{x})V^{\dagger}_{0}(\mathbf{x})}{H^{V}_{F}}\mathbf{P}_{\perp}^{V}\frac{\mathbf{V}(\mathbf{y})\cdot\mathbf{j}_{I}(\mathbf{y})}{H^{V}_{F}}\right\rangle.$
(13)

Expanding the quadratic vacuum fluctuations of the bosons in terms of normal
modes and adding up Eqs.(12) and (13), we arrive at (cf. appendix D)

$\displaystyle\langle
n^{0}_{\uparrow}|\mathbf{P}^{\textrm{V}}_{\parallel}+\mathbf{P}^{\textrm{V}}_{\perp}|n^{0}_{\uparrow}\rangle$
$\displaystyle=-g^{2}_{V}\textrm{Re}\sum_{I}\int d^{3}x\int
d^{3}y\int\frac{d^{3}q}{(2\pi)^{3}}$ $\displaystyle\times
j_{I}^{0}(\mathbf{x})\frac{e^{i\mathbf{q}\cdot(\mathbf{x}-\mathbf{y})}}{q^{2}+M^{2}_{V}}\>\mathbf{j}_{I}(\mathbf{y}).$
(14)

From Eq.(14) together with Eq.(11) we interpret that the variation of the
kinetic momentum of a polarized neutron results from the variation of the
(self)interaction between the transient weak charges $j_{I}^{0}$ and chiral
currents $\mathbf{j}_{I}$ of its constituent quarks. The Z and W-bosons which
mediate the interaction carry off an equivalent momentum in the opposite
direction. Since the chiral currents reverse direction when the quarks’ spins
are inverted, the spin-flip of a neutron results in a net variation of its
kinetic momentum, and explains the equivalent contribution of the two terms in
Eq.(11). The finiteness of the calculation is guaranteed by the short-range of
the weak nuclear forces. Hence, the momentum integral of Eq.(14) scales as
$\sim e^{-|x-y|M_{V}}/(xyM_{V})$ –see appendix D for further explanations.
Also worth noting is that Eq.(14) is the EW analog of the longitudinal
momentum of the EM field sourced by a chiral molecule,
$\mathbf{E}\times\mathbf{B}$, where the electrostatic field is sourced by the
electric charges and the magnetic field is created by their helical currents
[11].

Finally, considering that the major contribution from intermediate states
comes from the ground states, of neutrons for
$\mathbf{P}_{\parallel,\perp}^{Z}$ and protons for
$\mathbf{P}_{\parallel,\perp}^{W}$, the sum of the momenta in Eq.(11) yields

$\displaystyle\Delta\langle\mathbf{P}^{\textrm{sf}}_{\textrm{kin}}\rangle\simeq\hat{\mathbf{z}}\frac{-20\pi
G_{F}}{\sqrt{2}R_{n}^{3}}\mathcal{F}_{n},\textrm{ with }$ (15)
$\displaystyle\mathcal{F}_{n}=\int
d\tilde{r}\>\tilde{r}^{2}\left[F^{2}(\tilde{r})+G^{2}(\tilde{r})\right]\left[F^{2}(\tilde{r})-\frac{G^{2}(\tilde{r})}{3}\right].$

The factor $\mathcal{F}_{n}$ is a dimensionless integral which depends on the
upper and lower components of the quarks’ ground state wavefunction,
$\psi^{0,\uparrow}_{u,d}\approx\bigl{(}F(\tilde{r})_{\uparrow}\>\>G(\tilde{r})_{\uparrow}\bigr{)}$,
with $\tilde{r}=r/R_{n}$. In particular, for the MIT-bag relativistic model,
the total variation of the neutron velocity along the spin-flip is $13\>$m/s,
while for the nonrelativistic harmonic oscillator it is $7\>$m/s
approximately.

##  IV Conclusions and discussion

From a fundamental perspective, if this effect is verified experimentally, it
would demonstrate the existence of net internal forces in a system governed by
chiral interactions. The apparent contradiction with the action-reaction
principle is resolved by considering the momentum transferred to the virtual
bosons which mediate the interaction. In this respect, net internal forces
have also been found theoretically in binary systems of excited atoms [17].
However, the virtual character of the mediating photons is there not clear
since they are resonant. In the present case, Z and W-bosons are undoubtedly
virtual.

Concerning the accuracy of our estimate and its dependence on the hadron
model, such a dependence enters Eq.(15) through the spatial wavefunctions in
$\mathcal{F}_{n}$ only. That means that the energy spectrum of the nucleons is
irrelevant, as the energy intervals are negligible in Eqs.(12), (13) in
comparison to the boson energies in the denominators. Therefore, any
refinement upon our simplified models of the sort of one-gluon-exchange
potentials, spin-spin couplings, etc. [18], would be relevant in considering
their impact in the quarks’ wavefunctions only. More relevant is the spin
structure of the neutron, which has an impact in the spin-flavour functions of
$|n^{0}_{\uparrow,\downarrow}\rangle$ as well as in the orbital angular
momentum of the spatial wavefunctions in $\mathcal{F}_{n}$. We have considered
the simplest wavefunction in Eq.(10), but from the advent of the _spin crisis_
[19] it is recognized that the total spin of nucleons is not equal to the sum
of the spins of the valence quarks. In this respect, our finding may be used
to probe the hadronic structure and the angular momentum content of the
nucleons, and as a test for hadron models. Hence, the validity of our estimate
for the velocity reduces to its orders of magnitude. More realistic models for
nucleons are needed to improve its accuracy [20]. Ultracold neutron traps and
neutron spectroscopy are the experimental techniques proposed for its
observation [21].

###### Acknowledgements.

We thank B. van Tiggelen and G. Rikken for fruitful discussions on the
fundamentals of the phenomenon, R. Ballou for enlightening discussions on the
experimental techniques for its observation, and A. Cano for his hospitality
at Institut Néel. Financial support from the NextGenerationEU funds, through
the Spanish Ministerio de Universidades, is acknowledged.

##  Appendix A Hadron confinement models

For simplicity, two effective models of hadron confinement have been
considered in the article for the numerical estimations, namely, that of the
MIT-bag for valence and relativistic quarks [13], and the harmonic-oscillator
model for constituent and nonrelativistic quarks [14]. Formally, the only
remarkable difference is that only the latter allows for an explicit
separation between internal and external degrees of freedom.

###  A.1 Relativistic MIT-bag model

In the MIT-bag model (B), the Lagrangian density of Ref.[13] for the three
valence quarks bound in the neutron is, in the massless approximation,

$\mathcal{L}_{0}^{\textrm{B}}=\sum_{j=1}^{3}\Bigl{[}\frac{i}{2}\left(\bar{\psi}_{j}\gamma^{\mu}\partial_{\mu}\psi_{j}-\left(\partial_{\mu}\bar{\psi}_{j}\right)\gamma^{\mu}\partial_{\mu}\psi_{j}\right)-B\Bigr{]}\Theta(|\mathbf{r}_{j}-\mathbf{R}|-R_{n})-\frac{1}{2}\bar{\psi}_{j}\psi_{j}\delta(|\mathbf{r}_{j}-\mathbf{R}|-R_{n}).$
(16)

From this Lagrangian, one derives the Hamiltonian,

$H_{n}^{\textrm{B}}=-i\sum_{j=1}^{3}\int
d^{3}r_{j}\bar{\psi}_{j}\gamma^{k}\partial_{k}\psi_{j}\Theta(|\mathbf{r}_{j}-\mathbf{R}|-R_{n}),$
(17)

together with the boundary conditions,

$-i\hat{\mathbf{R}}_{n}\cdot\boldsymbol{\gamma}\psi_{j}|_{\partial\Omega_{n}}=\psi_{j}|_{\partial\Omega_{n}},\quad\hat{\mathbf{R}}_{n}\cdot\boldsymbol{\nabla}(\bar{\psi}_{j}\psi_{j})|_{\partial\Omega_{n}}=-2B.$
(18)

In the above equations the subscript $j$ runs over the three valence quarks,
with position vectors $\mathbf{r}_{j}$, $\mathbf{R}$ being the center of the
bag; $\\{\psi_{j}=\psi(\mathbf{r}_{j})\\}$ are Dirac spinors; $R_{n}$ is the
static radius of the spherical cavity with outwards unitary vector
$\hat{\mathbf{R}}_{n}$ at the surface $\partial\Omega_{n}$. The linear
boundary condition in Eq.(18) implies no current flow through the surface,
whereas the nonlinear one implies the balance between the outwards pressure of
the quarks and the inward QCD vacuum pressure $B$ at the surface.

The components of the spinor eigenfunctions $\\{\psi(\mathbf{r}_{j})\\}$ can
be split in upper and lower components, each one labeled with four common
numbers, $n$, $\kappa$, $J$ and $J_{3}$; and different orbital angular
momentum numbers, $l$ and $l^{\prime}$.

$\psi(\mathbf{r})=\left(\begin{array}[]{c}F_{n\kappa}(r)\mathcal{Y}_{l,1/2}^{J,J_{3}}(\hat{\mathbf{r}})\\\
\\\
i\>G_{n\kappa}(r)\mathcal{Y}_{l^{\prime},1/2}^{J,J_{3}}(\hat{\mathbf{r}})\end{array}\right)=\mathcal{N}_{n,\kappa}\left(\begin{array}[]{c}j_{l}(\mathcal{E}_{n,\kappa}r)\mathcal{Y}_{l,1/2}^{J,J_{3}}(\hat{\mathbf{r}})\\\
\\\
(-1)^{\eta}i\>j_{l^{\prime}}(\mathcal{E}_{n,\kappa}r)\mathcal{Y}_{l^{\prime},1/2}^{J,J_{3}}(\hat{\mathbf{r}})\end{array}\right).$

In these equations $\mathcal{N}_{n,\kappa}$ is a normalization constant,
$j_{l^{(^{\prime})}}$ are spherical Bessel functions,
$\mathcal{Y}_{l^{(^{\prime})},1/2}^{J,J_{3}}$ are spinor spherical harmonics,
and we use a representation in which Dirac’s matrices are given by

$\gamma_{0}=\left({\begin{array}[]{cc}\mathbb{I}_{2}&0\\\
0&-\mathbb{I}_{2}\end{array}}\right),\quad\gamma_{5}=\left({\begin{array}[]{cc}0&\mathbb{I}_{2}\\\
\mathbb{I}_{2}&0\end{array}}\right),\quad\gamma^{i}=\left({\begin{array}[]{cc}0&\sigma_{i}\\\
-\sigma_{i}&0\end{array}}\right).$

$\mathcal{E}_{n,\kappa}$ is the $n$-eigenenergy value, i.e., the $n^{th}$
solution of a transcendental equation for certain eigenvalue $\kappa$ of the
operator
$\mathcal{K}=\gamma_{0}(\boldsymbol{\Sigma}\cdot\mathbf{L}+\mathbb{I})=\gamma_{0}(\boldsymbol{\Sigma}\cdot\mathbf{J}-\mathbb{I}/2)$,
which derives from the equations of the boundary coniditions. For masseless
quarks,

$j_{l}(\mathcal{E}_{n,\kappa}R_{n})=(-1)^{\eta}j_{l+1}(\mathcal{E}_{n,\kappa}R_{n}),\textrm{
with }\eta=0\textrm{ for }\kappa<0,\eta=1\textrm{ for }\kappa>0.$ (19)

The number $\kappa$ may take positive and negative values which depend on the
total angular momentum $J$, $\kappa=\pm(J+1/2)$. The orbital angular momenta
are $l=-(\kappa+1)$ and $l^{\prime}=-\kappa$ for $\kappa<0$; and $l=\kappa$,
$l^{\prime}=\kappa-1$ for $\kappa>0$. In particular, the wavefunctions of the
first two levels with $J=1/2$ and opposite parity (i.e., $1S_{1/2}$ and
$1P_{1/2}$ states) are

$\psi_{1,-1}^{s}(\mathbf{r})=\mathcal{N}_{1,-1}\left(\begin{array}[]{c}j_{0}(\mathcal{E}_{1,-1}r)\>\chi_{s}\\\
\\\
i\>j_{1}(\mathcal{E}_{1,-1}r)\hat{\mathbf{r}}\cdot\boldsymbol{\sigma}\>\chi_{s}\end{array}\right),\quad\psi_{1,+1}^{s}(\mathbf{r})=\mathcal{N}_{1,+1}\left(\begin{array}[]{c}j_{1}(\mathcal{E}_{1,+1}r)\hat{\mathbf{r}}\cdot\boldsymbol{\sigma}\>\chi_{s}\\\
\\\ -ij_{0}(\mathcal{E}_{1,+1}r)\>\chi_{s}\end{array}\right),$

where $\chi_{s}$ are unit two-component spinors, $s=+$ for spin-up $\uparrow$,
and $s=-$ for-spin down $\downarrow$; and the eigenenergies are
$\mathcal{E}_{1,-1}=2.04/R_{n}$, $\mathcal{E}_{1,+1}=3.81/R_{n}$. It is
$\psi^{s}_{1,-1}(\mathbf{r}_{j})$ the one to be identified with
$\psi_{f_{j}}^{0,s_{j}}$ for any valence quark $j$ of flavour $f_{j}$ and spin
$s_{j}$, in the wavefunction of the neutron ground state –see below.

Finally, in order to fit the nucleon radius for any combination of
wavefunctions of the three valence quarks, one requires the balance between
the internal pressure exerted by the quarks and the external vacuum pressure,
$-B$, at the surface of the bag. This implies the boundary condition

$\frac{d\mathcal{E}_{tot}}{dr}|_{\partial\Omega_{n}}=\frac{d}{dr}\left[\frac{4\pi}{3}Br^{3}+\sum_{j=1}^{3}\mathcal{E}^{j}R_{n}/r-z_{0}/r\right]\Bigr{|}_{r=R_{n}}=0,$
(20)

where $\mathcal{E}^{j}$ is the energy of the wavefunction of quark $j$ and
$z_{0}$ is a numerical parameter of value 1.5 approximately which weights the
contribution of an additional QCD vacuum energy needed to fit the mass
spectrum together with the value of $B^{1/4}\approx 145\>$MeV [2, 22]. Note
that the normalization constants $\mathcal{N}_{n,\kappa}$ of each quark’s
wavefunction depend on the radius of the whole nucleus computed this way.

###  A.2 Nonrelativistic harmonic oscillator model

As for the nonrelativistic harmonic oscillator (HO) model, its Hamiltonian is

$\displaystyle H_{n}^{\textrm{HO}}$ $\displaystyle=\sum_{j=1}^{3}\int
d^{3}r_{j}\psi^{\dagger}_{j}\Bigl{[}-\frac{\nabla^{2}}{2m}+\frac{m\omega^{2}}{4}\sum_{k<j}^{3}\int
d^{3}r_{k}$
$\displaystyle\times\psi^{\dagger}_{k}|\mathbf{r}_{j}-\mathbf{r}_{k}|^{2}\psi_{k}\Bigr{]}\psi_{j},$
(21)

where the masses $m$ and the oscillator frequencies $\omega$ are equivalent
for the three constituent quarks [24]. In the nonrelativistic limit, the lower
component of the Dirac spinors is negligible, and the spin-space
eigenfunctions of the quarks can be approximated by tensor products of complex
scalar functions and two-component spinors. Further, the above Hamiltonian is
separable in terms of Jacobi’s coordinates. The EW interactions take place at
the location of one of the constituent quarks at a time. We refer to it as
_active quark_ and denote its position vector by $\mathbf{r}_{A}$. The
positions of the remaining quarks, i.e., the _spectators_ , will be denoted by
$\mathbf{r}_{1}$, $\mathbf{r}_{2}$ in the laboratory frame. Their
corresponding conjugate momenta are $\mathbf{p}_{A}$, $\mathbf{p}_{1}$ and
$\mathbf{p}_{2}$, respectively. The Jacobi position vectors and corresponding
conjugate momenta read [12, 23],

$\displaystyle\mathbf{R}$
$\displaystyle=\frac{\mathbf{r}_{1}+\mathbf{r}_{2}+\mathbf{r}_{A}}{3},\quad\mathbf{r}_{\rho}=\mathbf{r}_{2}-\mathbf{r}_{1},\quad\mathbf{r}_{\lambda}=\mathbf{r}_{A}-\frac{\mathbf{r}_{1}+\mathbf{r}_{2}}{2},$
$\displaystyle\mathbf{p}$
$\displaystyle=\mathbf{p}_{1}+\mathbf{p}_{2}+\mathbf{p}_{A},\>\>\mathbf{p}_{\rho}=\frac{\mathbf{p}_{2}-\mathbf{p}_{1}}{2},\>\>\mathbf{p}_{\lambda}=\frac{2\mathbf{p}_{A}-\mathbf{p}_{1}-\mathbf{p}_{2}}{3},$

where $\mathbf{R}$ is the centre of mass of the nucleon and $\mathbf{p}$ its
conjugate momentum. As for the active quark, its position vector
$\mathbf{r}_{A}$ reads, in terms of Jacobi’s coordinates,
$\mathbf{r}_{A}=\mathbf{R}+2\mathbf{r}_{\lambda}/3$. In terms of Jacobi’s
coordinates and momenta, the Hamiltonian is

$H_{n}^{\textrm{HO}}=p^{2}/6m+p_{\rho}^{2}/m+3p_{\lambda}^{2}/4m+m\omega^{2}(r_{\rho}^{2}/4+r_{\lambda}^{2}/3).$
(22)

The eigenfunctions in terms of Jacobi’s coordinates reduce to the product of
the wavefunctions of two independent and isotropic harmonic oscillators times
plane wavefunctions for the free motion of the center of mass. In particular,
the wavefunction for the ground state of the neutron with total conjugate
momentum $\mathbf{p}$ is

$\Psi^{0,s}_{n}(\mathbf{R},\mathbf{r}_{\lambda},\mathbf{r}_{\rho})=\frac{e^{-i\mathbf{p}\cdot\mathbf{R}}}{\sqrt{\mathcal{V}}}\frac{\beta^{3}}{3^{3/4}\pi^{3/2}}e^{-\beta^{2}r_{\lambda}^{2}/3}e^{-\beta^{2}r_{\rho}^{2}/4}\chi_{s},$

with $\beta=\sqrt{m\omega}$, $\mathcal{V}$ being an infinite volume, and
$\chi_{s}$ a unit two-component spinor. Using the above formula for the
calculation of $\Delta\langle\mathbf{P}_{\textrm{kin}}^{\textrm{sf}}\rangle$,
integrating in $\mathbf{R}$ and $\mathbf{r}_{\rho}$, getting back to the
coordinates of the laboratory frame with $\tilde{x}\equiv
2r_{\lambda}/3R_{n}$, and taking $R_{n}=2\beta^{-1}\sqrt{2/3\pi}$ as the mean
distance of the active quark to the center of the neutron, we find out
$F(\tilde{x})=(2/\pi)^{3/2}e^{-2\tilde{x}^{2}/\pi}$ in the expression for
$\mathcal{F}_{n}$ in Eq.(15). Note also that, for consistency with the actual
values of the mass and the radius of the neutron, $\omega\approx 0.4\>m$,
which implies that excited states would be inconsistent with the
nonrelativistic assumption and thus are to be discarded.

##  Appendix B Spin-flip procedure and EM momentum

It is mentioned in the article that a simple and realistic manner to invert
the spin of the neutron in a coherent manner is by means of a $\pi$-pulse with
a radiofrequency (RF) field. Here, we detail the features of this field as
well as its action upon the neutron.

Let us consider a free neutron initially polarized along the
$-\hat{\mathbf{z}}$ direction in the presence of an intense constant magnetic
field $\mathbf{B}_{0}=B_{0}\hat{\mathbf{z}}$. To invert the spin, an
additional magnetic field $\mathbf{B}_{1}$ is introduced which rotates at the
Larmor frequency $\omega_{0}=-\gamma_{n}B_{0}$ in the $xy$ plane, with
$\gamma_{n}$ being the neutron gyromagnetic ratio,
$\mathbf{B}_{1}(t)=B_{1}[\cos{(\omega_{0}t)}\hat{\mathbf{x}}+\sin{(\omega_{0}t)}\hat{\mathbf{y}}]$.
This field is circularly polarized around $\hat{\mathbf{z}}$, and its strength
$B_{1}$ is much weaker than $B_{0}$. The resultant external EM vector
potential can be written as

$\mathbf{A}_{\textrm{ext}}=\frac{B_{0}}{2}(-y\hat{\mathbf{x}}+x\hat{\mathbf{y}})+\frac{B_{1}}{\omega_{0}}[\cos{(k_{0}z-\omega_{0}t)}\hat{\mathbf{x}}-\sin{(k_{0}z-\omega_{0}t)}\hat{\mathbf{y}}],\quad
k_{0}z\rightarrow 0,$ (23)

with photon density $\mathcal{N}_{\gamma}=B_{1}^{2}/\omega_{0}$. The coupling
of these photons to the constituent/valence quarks is through the Hamiltonian
$\mathcal{W}_{EM}$. All in all that interaction simplifies to
$W_{EM}\approx-\gamma_{n}\boldsymbol{\sigma}\cdot\mathbf{B}/2$, where
$\boldsymbol{\sigma}$ applies to the two-component spinor wavefunction of the
neutron, $\alpha_{1}|n\uparrow\rangle+\alpha_{2}|n\downarrow\rangle$,
$|\alpha_{1}|^{2}+|\alpha_{2}|^{2}=1$, where the spin-flavour wavefunction of
Eq.(10) is [1]

$\displaystyle|n\downarrow\rangle=\frac{1}{\sqrt{18}}\Bigl{[}$
$\displaystyle|u\downarrow d\uparrow d\downarrow\rangle+|d\downarrow d\uparrow
u\downarrow\rangle+|d\uparrow d\downarrow u\downarrow\rangle+|d\uparrow
u\downarrow d\downarrow\rangle+|u\downarrow d\downarrow
d\uparrow\rangle+|d\downarrow u\downarrow d\uparrow\rangle$
$\displaystyle-2(|u\uparrow d\downarrow d\downarrow\rangle+|d\downarrow
u\uparrow d\downarrow\rangle+|d\downarrow d\downarrow
u\uparrow\rangle)\Bigr{]},$ (24)

and likewise for $|n\uparrow\rangle$ with all the spins reversed –the
superscript $0$ in the expression for $|n^{0}_{\uparrow}\rangle$ in the
article refers to the $l=0$ ground state of the spatial wavefunctions.
Further, replacing $\mathbf{B}$ in $W_{EM}$ with the classical expression of
$\mathbf{B}_{0}+\mathbf{B}_{1}$ or using the quantum expression for
$\mathbf{B}$ and the EM states containing the photons of the field
$\mathbf{B}_{1}$ otherwise, one can apply time-dependent perturbation theory
and compute an effective Rabi (R) Hamiltonian for the interaction between the
RF field and the neutron. That is,

$W^{\textrm{R}}_{EM}\approx-\frac{\omega_{0}}{2}|n\downarrow\rangle\langle
n\downarrow|+\frac{\omega_{0}}{2}|n\uparrow\rangle\langle
n\uparrow|+-\frac{\omega_{1}}{2}(|n\downarrow\rangle\langle
n\uparrow|+|n\uparrow\rangle\langle n\downarrow|),$ (25)

where $\omega_{1}=-\gamma_{n}B_{1}$ is the Rabi frequency. This Hamiltonian
accounts for an infinite number of interaction processes involving multiple
and alternate absorptions and re-emissions of photons of the field
$\mathbf{B}_{1}$. Starting with the neutron ground state
$|n^{0}_{\downarrow}\rangle$ at the initial time $t_{0}=0$, the neutron spin
wavefunction reads, at a later time $t$,

$|n^{0}(t)\rangle=\cos{(\omega_{1}t/2)}|n^{0}_{\downarrow}\rangle-i\sin{(\omega_{1}t/2)}e^{-i\omega_{0}t}|n^{0}_{\uparrow}\rangle,$
(26)

In turn, this makes the total neutron spin oscillate at frequency
$\omega_{1}$,

$\langle\mathbf{S}_{n}(t)\rangle=\frac{-1}{2}\left[\sin{(\omega_{1}t)}\sin{(\omega_{0}t)}\hat{\mathbf{x}}-\sin{(\omega_{1}t)}\cos{(\omega_{0}t)}\hat{\mathbf{y}}+\cos{(\omega_{1}t)}\hat{\mathbf{z}}\right],$
(27)

inverting its direction in a time interval $\pi/\omega_{1}$ and absorbing one
of the photons of frequency $\omega_{0}$ and momentum $k_{0}\hat{\mathbf{z}}$.
Since the EM interaction is parity invariant, in the computation of
$\Delta\langle\mathbf{P}_{\textrm{kin}}^{\textrm{sf}}\rangle$ along the spin-
flip the transfer of EM momentum takes place only through the absorption of an
actual photon of frequency $\omega_{0}$,

$\langle
n^{0}_{\downarrow},N_{\gamma}|\mathbf{P}_{\perp}^{A}|n^{0}_{\downarrow},N_{\gamma}\rangle-\langle
n^{0}_{\uparrow},N_{\gamma}-1|\mathbf{P}_{\perp}^{A}|n^{0}_{\uparrow},N_{\gamma}-1\rangle=k_{0}\hat{\mathbf{z}}.$
(28)

Typical values for $\omega_{0}$ in experiments range between $10^{5}$ and
$10^{8}\>$Hz, which corresponds to a variation of the neutron velocity of the
order of $10^{-10}$-$10^{-7}\>$m/s, which is negligibly small in comparison to
that provided by the virtual Z and W-bosons

##  Appendix C Quantum perturbative calculation

As indicated in the article, the calculation of the right-hand-side of Eq.(11)
for $\Delta\langle\mathbf{P}_{\textrm{kin}}^{\textrm{sf}}\rangle$ is performed
in the framework of time-dependent quantum perturbation theory in
Schrödinger’s picture. Also, since the frequency of the spin oscillations,
$\omega_{1}$, is much smaller than any transition frequency and decay rate of
the intermediate processes involved, the calculation is adiabatic with respect
to the spin dynamics. In good approximation, the transition frequencies
correspond to the bosons’ masses, $M_{Z}=91.2\>$GeV, $M_{W}=80.4\>$GeV, and
their decay rates are $\Gamma_{Z}\simeq 2.5\>GeV$, $\Gamma_{W}\simeq
2.1\>GeV$, satisfying $M_{Z,W}\gg\Gamma_{W,Z}\gg\omega_{1}$. Defining the non-
perturbative Hamiltonian $H_{0}$ as that made of the neutron Hamiltonian, the
Hamiltonians of the free boson fields and the effective Rabi Hamiltonian which
causes the oscillations of the neutron spin,
$H_{0}=H_{F}^{A}+H_{F}^{Z}+H_{F}^{W^{\pm}}+H_{n}+W^{R}_{EM}$, the
nonperturbative time-evolution propagator is
$\mathbb{U}_{0}(\tau)=\exp{(-i\tau\>H_{0})}$. At leading order in the
perturbative potentials, $\mathcal{W}_{nc}$ and $\mathcal{W}_{cc}$, the
expectation values of the longitudinal momenta of Z and W-bosons are computed
at first order, while those of the transverse momenta are computed at second
order. These are, respectively,

$\displaystyle\langle\mathbf{P}^{\textrm{Z,W}}_{\parallel}(t)\rangle$
$\displaystyle=2\textrm{Re}\langle
n^{0}_{\downarrow}|\mathbb{U}_{0}^{\dagger}(t)\mathbf{P}^{\textrm{Z,W}}_{\parallel}\int_{0}^{t}-i\>d\tau
e^{-\Gamma_{Z,W}(t-\tau)/2}\mathbb{U}_{0}(t-\tau)W_{nc,cc}\mathbb{U}_{0}(\tau)|n^{0}_{\downarrow}\rangle,$
(29) $\displaystyle\langle\mathbf{P}^{\textrm{Z,W}}_{\perp}(t)\rangle$
$\displaystyle=\langle
n^{0}_{\downarrow}|\int_{0}^{t}d\tau^{\prime}e^{-\Gamma_{Z,W}(t-\tau^{\prime})/2}\mathbb{U}_{0}^{\dagger}(\tau^{\prime})W^{\dagger}_{nc,cc}\mathbb{U}_{0}^{\dagger}(t-\tau^{\prime})\>\mathbf{P}^{\textrm{Z,W}}_{\perp}$
$\displaystyle\times\int_{0}^{t}d\tau
e^{-\Gamma_{Z,W}(t-\tau)/2}\mathbb{U}_{0}(t-\tau)W_{nc,cc}\mathbb{U}_{0}(\tau)|n^{0}_{\downarrow}\rangle.$
(30)

The variation of the kinetic momentum all along the spin-flip,
$\Delta\langle\mathbf{P}_{\textrm{kin}}^{\textrm{sf}}\rangle$, comprises the
difference between the evaluation of the above equations for
$t\rightarrow\pi/\omega_{1}$ (i.e., for
$|n(\pi/\omega_{1})\rangle=|n^{0}_{\uparrow}\rangle$) and $t\rightarrow 0^{+}$
with $t\Gamma_{Z,W}\gg 1$ (i.e., for
$|n(0)\rangle=|n^{0}_{\downarrow}\rangle$). The diagrams of Fig.1 depict two
of the processes contributing to $\langle
n^{0}_{\uparrow}|\mathbf{P}^{\textrm{W}}_{\parallel}|n^{0}_{\uparrow}\rangle$
[(a)] and $\langle
n^{0}_{\uparrow}|\mathbf{P}^{\textrm{Z}}_{\perp}|n^{0}_{\uparrow}\rangle$
[(b)]. Note that a physically equivalent calculation consists of implementing
the spin-flip process adiabatically. To do so it suffices to replace in the
above formulas the attenuating exponential factors
$e^{-\Gamma_{Z,W}(t-\tau)/2}$ with $e^{-\eta(t-\tau)/2}$, with
$\eta\rightarrow 0^{+}$, and to replace the initial time $0$ with
$t_{0}\rightarrow-\infty$ such that $-\eta t_{0}\gg 1$.

##  Appendix D Calculation of
$\Delta\langle\mathbf{P}_{\textrm{kin}}^{\textrm{sf}}\rangle$

The evaluation of Eq.(29) at $t=\pi/\omega_{1}$ yields, in each case,

$\langle
n^{0}_{\uparrow}|\mathbf{P}^{\textrm{Z}}_{\parallel}|n^{0}_{\uparrow}\rangle=-2\textrm{Re}\sum_{I}\left\langle\frac{\langle
n^{0}_{\uparrow}|\mathbf{P}^{\textrm{Z}}_{\parallel}|n_{I}\rangle\langle
n_{I}|W_{nc}|n^{0}_{\uparrow}\rangle}{H^{Z}_{F}+\mathcal{E}_{I}-m_{n}-i\Gamma_{Z}/2}\right\rangle,\quad\langle
n^{0}_{\uparrow}|\mathbf{P}^{\textrm{W}}_{\parallel}|n^{0}_{\uparrow}\rangle=-2\textrm{Re}\sum_{I}\left\langle\frac{\langle
n^{0}_{\uparrow}|\mathbf{P}^{\textrm{W}}_{\parallel}|p_{I}\rangle\langle
p_{I}|W_{cc}|n^{0}_{\uparrow}\rangle}{H^{W}_{F}+\mathcal{E}_{I}-m_{n}-i\Gamma_{W}/2}\right\rangle,$
(31)

where the big brackets denote the expectation value of the boson field
operators in the vacuum state. As for the intermediate nucleon states, labeled
with the script $I$, they are neutron states $\\{|n_{I}\rangle\\}$ for
$\langle\mathbf{P}^{\textrm{Z}}_{\parallel}\rangle$ and protons
$\\{|p_{I}\rangle\\}$ for $\langle\mathbf{P}^{\textrm{W}}_{\parallel}\rangle$,
with energies $\mathcal{E}_{I}$. The terms
$\mathcal{E}_{I}-m_{n}-i\Gamma_{Z,W}/2$ in the denominators of Eq.(31) are
negligible, since they are much less than the boson masses $M_{Z,W}$. Next, it
can be verified that, for the MIT-bag model, the overlap between the spatial
wavefunction of the ground state and those of excited states is small. On the
other hand, in the HO model the excited sates exceed the nonrelativistic
limit. Hence, it suffices hereafter to consider only the spin-up ground states
of neutrons and protons in the sums over intermediate states,
$|n^{0}_{\uparrow}\rangle$ and $|p^{0}_{\uparrow}\rangle$. Disregarding the
irrelevant color function, the wavefunction of the neutron is

$\displaystyle\Psi^{0}_{n\uparrow}(\mathbf{r}_{1},\mathbf{r}_{2},\mathbf{r}_{3})$
$\displaystyle=\frac{1}{\sqrt{18}}\Bigl{[}\psi^{u,+}_{1,-1}(\mathbf{r}_{1})\psi^{d,-}_{1,-1}(\mathbf{r}_{2})\psi^{d,+}_{1,-1}(\mathbf{r}_{3})+\psi^{d,+}_{1,-1}(\mathbf{r}_{1})\psi^{d,-}_{1,-1}(\mathbf{r}_{2})\psi^{u,+}_{1,-1}(\mathbf{r}_{3})+\psi^{d,-}_{1,-1}(\mathbf{r}_{1})\psi^{d,+}_{1,-1}(\mathbf{r}_{2})\psi^{u,+}_{1,-1}(\mathbf{r}_{3})$
$\displaystyle+\psi^{d,-}_{1,-1}(\mathbf{r}_{1})\psi^{u,+}_{1,-1}(\mathbf{r}_{2})\psi^{d,+}_{1,-1}(\mathbf{r}_{3})+\psi^{u,+}_{1,-1}(\mathbf{r}_{1})\psi^{d,+}_{1,-1}(\mathbf{r}_{2})\psi^{d,-}_{1,-1}(\mathbf{r}_{3})+\psi^{d,+}_{1,-1}(\mathbf{r}_{1})\psi^{u,+}_{1,-1}(\mathbf{r}_{2})\psi^{d,-}_{1,-1}(\mathbf{r}_{3})$
(32)
$\displaystyle-2\bigl{(}\psi^{u,-}_{1,-1}(\mathbf{r}_{1})\psi^{d,+}_{1,-1}(\mathbf{r}_{2})\psi^{d,+}_{1,-1}(\mathbf{r}_{3})+\psi^{d,+}_{1,-1}(\mathbf{r}_{1})\psi^{u,-}_{1,-1}(\mathbf{r}_{2})\psi^{d,+}_{1,-1}(\mathbf{r}_{3})+\psi^{d,+}_{1,-1}(\mathbf{r}_{1})\psi^{d,+}_{1,-1}(\mathbf{r}_{2})\psi^{u,-}_{1,-1}(\mathbf{r}_{3})\bigr{)}\Bigr{]},$

and likewise for $\Psi^{0}_{p\uparrow}$ but for the exchange $d\leftrightarrow
u$. Taking all this into account, Eq.(31) can be written as

$\displaystyle\langle
n^{0}_{\uparrow}|\mathbf{P}^{\textrm{Z}}_{\parallel}|n^{0}_{\uparrow}\rangle$
$\displaystyle=\frac{-g}{\cos{\theta_{W}}}\int d^{3}x\int
d^{3}y\left\langle\boldsymbol{\mathcal{P}}_{\parallel,I}^{Z\dagger}(\mathbf{x})\frac{\mathbf{Z}(\mathbf{y})\cdot\mathbf{j}_{nc}^{I}(\mathbf{y})}{H^{Z}_{F}}\right\rangle,$
(33) $\displaystyle\langle
n^{0}_{\uparrow}|\mathbf{P}^{\textrm{W}}_{\parallel}|n^{0}_{\uparrow}\rangle$
$\displaystyle=-\sqrt{2}g\int d^{3}x\int
d^{3}y\left\langle\boldsymbol{\mathcal{P}}_{\parallel,I}^{W^{\pm}dagger}(\mathbf{x})\frac{\mathbf{W}^{\mp}(\mathbf{y})\cdot\mathbf{j}_{cc}^{I(\dagger)}(\mathbf{y})}{H^{W}_{F}}\right\rangle,$
(34)

where the transient momentum density vectors and currents are evaluated
between the initial state $|n^{0}_{\uparrow}\rangle$ and
$|I\rangle=\\{|n^{0}_{\uparrow}\rangle,|p^{0}_{\uparrow}\rangle\\}$ in each
case,

$\displaystyle\boldsymbol{\mathcal{P}}_{\parallel,I}^{Zdagger}(\mathbf{x})$
$\displaystyle=\frac{g}{2\cos{\theta_{W}}}j_{nc,v}^{0,I}(\mathbf{x})\>\mathbf{Z}(\mathbf{x})=\frac{g}{2\cos{\theta_{W}}}(g_{v}^{u}+2g_{v}^{d})\psi^{s\dagger}_{1,-1}(\mathbf{x})\psi^{s^{\prime}}_{1,-1}(\mathbf{x})\mathbf{Z}(\mathbf{x})\delta_{ss^{\prime}},\textrm{
and}$ (35)
$\displaystyle\boldsymbol{\mathcal{P}}_{\parallel,I}^{W^{\pm}dagger}(\mathbf{x})$
$\displaystyle=\frac{g}{\sqrt{2}}j_{cc,v}^{0,I}(\mathbf{x})\mathbf{W}^{\pm}(\mathbf{x})=\frac{g}{2\sqrt{2}}\psi^{s\dagger}_{1,-1}(\mathbf{x})\psi^{s^{\prime}}_{1,-1}(\mathbf{x})\mathbf{W}^{\pm}(\mathbf{x})\delta_{ss^{\prime}}.$
(36)

Note that the weak charge density factors are only vector-like ($v$),
$j_{nc/cc,v}^{0,I}$, which are accompanied by delta functions of spin
conservation, and whose signs do not depend on the spin orientation. Using the
nomenclature outlined for the MIT-bag wavefunctions,
$|\psi^{s}_{1,-1}(\mathbf{x})|^{2}=F_{1-1}^{2}(x)+G_{1-1}^{2}(x)$. As for the
transient currents, they are axial currents which also preserve the spin.
Their only nonvanishing contributions come from their components along the
neutron spin, i.e., along $\hat{\mathbf{z}}$,

$\displaystyle\hat{\mathbf{z}}\cdot\mathbf{j}_{nc}^{I}(\mathbf{y})$
$\displaystyle=-\frac{4g_{a}^{d}-g_{a}^{u}}{3}\bar{\psi}^{s^{\prime}}_{1,-1}(\mathbf{y})\gamma^{3}\gamma_{5}\psi^{s}_{1,-1}(\mathbf{y})\delta_{ss^{\prime}},$
$\displaystyle\hat{\mathbf{z}}\cdot\mathbf{j}_{cc}^{I(\dagger)}(\mathbf{y})$
$\displaystyle=-\frac{5}{6}\bar{\psi}^{s^{\prime}}_{1,-1}(\mathbf{y})\gamma^{3}\gamma_{5}\psi^{s}_{1,-1}(\mathbf{y})\delta_{ss^{\prime}}.$
(37)

These currents are the result of the addition of the contributions of each
active quark with position vector $\mathbf{y}$, whose signs do depend on the
spin orientation, $+$ for $\uparrow$ and $-$ for $\downarrow$,
$\bar{\psi}^{s}_{1,-1}(\mathbf{y})\gamma^{3}\gamma_{5}\psi^{s}_{1,-1}(\mathbf{y})=F_{1-1}^{2}(y)+\cos{2\theta_{y}}\>G_{1-1}^{2}(y)(\delta_{s+}-\delta_{s-})$
[25]. It is the prefactors of the above equations that account for the
combination of the different quark spins and weights of each of the components
of $|n^{0}_{\uparrow}\rangle$ and $|p^{0}_{\uparrow}\rangle$. It is worth
noting that the combination of transient vector-like charges and transient
axial-like currents is a consequence of the fact that the quarks’
wavefunctions have all the same parity in $\Psi_{n\uparrow}^{0}$ and
$\Psi_{p\uparrow}^{0}$. If they had different parity the combination would be
axial-like charges with vector-like currents instead –eg., for a transition
between states with wavefunctions $\psi_{1,-1}^{s}$ and $\psi_{1,+1}^{s}$.

Writing now the quadratic vacuum fluctuations of the boson fields in terms of
normal modes and replacing the field Hamiltonian by their corresponding
eigenvalues [3],

$\left\langle
V^{\dagger}_{\mu}(\mathbf{x})\frac{1}{H_{F}^{V}}V_{\nu}(\mathbf{y})\right\rangle=\int\frac{d^{3}q}{(2\pi)^{3}}\frac{e^{\pm
i\mathbf{q}\cdot(\mathbf{x}-\mathbf{y})}}{-2(\mathcal{E}^{V}_{q})^{2}}\left(\eta_{\mu\nu}-\frac{q_{\mu}q_{\nu}}{M_{V}^{2}}\right),\quad
V=\\{Z,W^{\pm}\\},\textrm{and }\eta_{\mu\nu}\textrm{ being the Minkowski
metric,}$ (38)

and inserting the values for $g_{v,a}^{u,d}$ [1] we arrive at

$\displaystyle\langle
n^{0}_{\uparrow}|\mathbf{P}^{\textrm{Z}}_{\parallel}|n^{0}_{\uparrow}\rangle$
$\displaystyle=\frac{\hat{\mathbf{z}}5g^{2}}{24\cos^{2}{\theta_{W}}}\int\frac{d^{3}q}{(2\pi)^{3}}\int
d^{3}x\int
d^{3}y\>\frac{e^{i\mathbf{q}\cdot(\mathbf{x}-\mathbf{y})}\left(1+\frac{q_{z}^{2}}{M_{Z}^{2}}\right)}{2(q^{2}+M_{Z}^{2})}$
$\displaystyle\times[F_{1-1}^{2}(x)+G_{1-1}^{2}(x)][F_{1-1}^{2}(y)+\cos{2\theta_{y}}G_{1-1}^{2}(y)],$
(39) $\displaystyle\langle
n^{0}_{\uparrow}|\mathbf{P}^{\textrm{W}}_{\parallel}|n^{0}_{\uparrow}\rangle$
$\displaystyle=\frac{\hat{\mathbf{z}}5g^{2}}{12}\int\frac{d^{3}q}{(2\pi)^{3}}\int
d^{3}x\int
d^{3}y\>\frac{e^{i\mathbf{q}\cdot(\mathbf{x}-\mathbf{y})}\left(1+\frac{q_{z}^{2}}{M_{W}^{2}}\right)}{2(q^{2}+M_{W}^{2})}$
$\displaystyle\times[F_{1-1}^{2}(x)+G_{1-1}^{2}(x)][F_{1-1}^{2}(y)+\cos{2\theta_{y}}G_{1-1}^{2}(y)].$
(40)

Likewise, the evaluation of Eq.(30) at $t=\pi/\omega_{1}$ yields

$\langle
n^{0}_{\uparrow}|\mathbf{P}^{\textrm{Z}}_{\perp}|n^{0}_{\uparrow}\rangle=\sum_{I}\left\langle\frac{\langle
n^{0}_{\uparrow}|W_{nc}|n_{I}\rangle\mathbf{P}^{\textrm{Z}}_{\perp}\langle
n_{I}|W_{nc}|n^{0}_{\uparrow}\rangle}{(H^{Z}_{F}+\mathcal{E}_{I}-m_{n})^{2}+\Gamma_{Z}^{2}}\right\rangle,\quad\langle
n^{0}_{\uparrow}|\mathbf{P}^{\textrm{W}}_{\perp}|n^{0}_{\uparrow}\rangle=\sum_{I}\left\langle\frac{\langle
n^{0}_{\uparrow}|W_{cc}|p_{I}\rangle\mathbf{P}^{\textrm{W}}_{\perp}\langle
p_{I}|W_{cc}|n^{0}_{\uparrow}\rangle}{(H^{W}_{F}+\mathcal{E}_{I}-m_{n})^{2}+\Gamma_{W}^{2}}\right\rangle,$
(41)

which, restricting the sums over intermediate nucleon states to
$|n_{I}\rangle=|n^{0}_{\uparrow}\rangle$,
$|p_{I}\rangle=|p^{0}_{\uparrow}\rangle$ in either case, give

$\displaystyle\langle
n^{0}_{\uparrow}|\mathbf{P}^{\textrm{Z}}_{\perp}|n^{0}_{\uparrow}\rangle$
$\displaystyle=\frac{g^{2}}{2\cos^{2}{\theta_{W}}}\int d^{3}x\int
d^{3}y\left\langle\frac{j_{nc,v}^{0,I}(\mathbf{x})Z_{0}(\mathbf{x})}{H^{Z}_{F}}\mathbf{P}_{\perp}^{Z}\frac{\mathbf{Z}(\mathbf{y})\cdot\mathbf{j}_{nc}^{I}(\mathbf{y})}{H^{Z}_{F}}\right\rangle,$
(42) $\displaystyle\langle
n^{0}_{\uparrow}|\mathbf{P}^{\textrm{W}}_{\perp}|n^{0}_{\uparrow}\rangle$
$\displaystyle=g^{2}\int d^{3}x\int
d^{3}y\left\langle\frac{j_{cc,v}^{0,I}(\mathbf{x})W^{\pm}_{0}(\mathbf{x})}{H^{W}_{F}}\mathbf{P}_{\perp}^{W}\frac{\mathbf{W}^{\mp}(\mathbf{y})\cdot\mathbf{j}_{cc}^{I(\dagger)}(\mathbf{y})}{H^{W}_{F}}\right\rangle.$
(43)

Analogous steps to the ones taken for the longitudinal momenta lead to

$\displaystyle\langle
n^{0}_{\uparrow}|\mathbf{P}^{\textrm{Z}}_{\perp}|n^{0}_{\uparrow}\rangle$
$\displaystyle=\frac{-\hat{\mathbf{z}}5g^{2}}{24\cos^{2}{\theta_{W}}}\int\frac{d^{3}q}{(2\pi)^{3}}\int
d^{3}x\int
d^{3}y\>\frac{q_{z}^{2}\>e^{i\mathbf{q}\cdot(\mathbf{x}-\mathbf{y})}}{2M_{Z}^{2}(q^{2}+M_{Z}^{2})}$
$\displaystyle\times[F_{1-1}^{2}(x)+G_{1-1}^{2}(x)][F_{1-1}^{2}(y)+\cos{2\theta_{y}}G_{1-1}^{2}(y)],$
(44) $\displaystyle\langle
n^{0}_{\uparrow}|\mathbf{P}^{\textrm{W}}_{\perp}|n^{0}_{\uparrow}\rangle$
$\displaystyle=\frac{-\hat{\mathbf{z}}5g^{2}}{12}\int\frac{d^{3}q}{(2\pi)^{3}}\int
d^{3}x\int
d^{3}y\>\frac{q_{z}^{2}\>e^{i\mathbf{q}\cdot(\mathbf{x}-\mathbf{y})}}{2M_{W}^{2}(q^{2}+M_{W}^{2})}$
$\displaystyle\times[F_{1-1}^{2}(x)+G_{1-1}^{2}(x)][F_{1-1}^{2}(y)+\cos{2\theta_{y}}G_{1-1}^{2}(y)],$
(45)

By comparing Eqs.(39) and (40) with Eqs.(44) and (45) we note that the
transverse momenta cancel out identically with the second $q$-dependent terms
of the longitudinal momenta. Finally, we calculate the remaining integrals
that account for the (self)interaction of the transient weak charges with the
chiral currents, mediated by Z and W-bosons. That is,

$\displaystyle\int\frac{d^{3}q}{(2\pi)^{3}}\int d^{3}x\int
d^{3}y\>\frac{e^{i\mathbf{q}\cdot(\mathbf{x}-\mathbf{y})}}{q^{2}+M_{V}^{2}}\mathcal{M}(x)\mathcal{N}(y)$
$\displaystyle=\int dx\>x^{2}\int
dy\>y^{2}\frac{2\pi}{M_{V}xy}\left[e^{-M_{V}|x-y|}-e^{-M_{V}(x+y)}\right]\mathcal{M}(x)\mathcal{N}(y)$
$\displaystyle\simeq\frac{4\pi}{M_{V}}\int
dx\>x^{2}\mathcal{M}(x)\mathcal{N}(x),\quad V=\\{Z,W^{\pm}\\},$ (46)

where in the last approximation we have made use of the fact that the
integrals extend over a distance of the order of the neutron radius, and
$R_{n}\gg 1/M_{V}$. The insertion of this formula into Eqs.(39), (40), (44)
and (45) leads to the formula for $\mathcal{F}_{n}$ of Eq.(15). The numerical
estimations were made there, for the MIT-bag model, with
$F_{1-1}(x)=j_{0}(2.04\>x/R_{n}),G_{1-1}(x)=j_{1}(2.04\>x/R_{n})$ and
$R_{n}=1.1$ fm.

## References

  * [1] J.F. Donoghe, E. Golowich and B.R. Holstein, Dynamics of the Standard Model, Cambridge University Press, Cambridge (2022); J.C. Burgess and G. Moore, The Standard Model: A Primer, Cambridge University Press, Cambridge (2012).
  * [2] W. Greiner and B. Müller, Gauge Theory of Weak Interactions, Springer-Verlag, Berlin (19966).
  * [3] M.E. Peskin and D.V. Schroeder, An Introduction to Quantum Field Theory, Westview Press, Chicago (1995); C. Itzykson and J.B. Zuber, Quantum Field Theory, McGraw-Hill Inc. (1980).
  * [4] C.S. Wu, E. Ambler, R.W. Hayward, D.D. Hoppes, and R.P. Hudson, Experimental Test of Parity Conservation in Beta Decay, Phys. Rev. 105, 1413 (1957).
  * [5] J.D. Jackson, S.B. Treiman and H.W. Wyld. Jr., Possible Test of Time Reversal Invariance in Beta Decay, Phys. Rev. 106, 517 (1957); J.S. Nico, Neutron Beta Decay, J. Phys. G: Nucl. Part. Phys. 36, 104001 (2009); M.S. Sozzi, Discrete Symmetries and CP Violation. From Experiment to Theory, Oxford University Press Inc., New York (2008).
  * [6] P. Souder and K.D. Paschke, Parity Violation in Electron Scattering, Front. Phys. 11, 111301 (2016); J. Erler, C.J. Horowitz, S. Mantry and P. A. Souder, Weak Polarized Electron Scattering, Annu. Rev. Nucl. Part. Sci. 64, 269 (2014).
  * [7] A. Feigel, Quantum Vacuum Contribution to the Momentum of Dielectric Media, Phys. Rev. Lett. 92, 020404 (2004).
  * [8] B.A. van Tiggelen, G.L.J.A. Rikken and V. Krstić, Momentum Transfer from Quantum Vacuum to Magnetoelectric Matter, Phys. Rev. Lett. 96, 130402 (2006); O.A. Croze, Alternative derivation of the Feigel effect and call for its experimental verification, Proc. R. Soc. A 468, 429 (2012).
  * [9] S. Kawka and B.A. van Tiggelen, Quantum Electrodynamics of Casimir Momentum: Momentum of the Quantum Vacuum?, EPL 89, 11002 (2010). B.A. van Tiggelen, S. Kawka and G.L.J.A. Rikken, QED Corrections to the Electromagnetic Abraham Force. Casimir Momentum of the Hydrogen Atom?, Eur. Phys. J. D 66, 272 (2012).
  * [10] M. Donaire, B.A. van Tiggelen and G.L.J.A. Rikken, Casimir Momentum of a Chiral Molecule in a Magnetic Field, Phys. Rev. Lett. 111, 143602 (2013).
  * [11] M. Donaire, B.A. van Tiggelen and G.L.J.A. Rikken, Transfer of linear momentum from the quantum vacuum to a magnetochiral molecule, _J. Phys.: Condens. Matter_ 27, 214002 (2015).
  * [12] M. Donaire, Acceleration of an Unpolarized Proton along a Uniform Magnetic Field: Casimir Momentum of Leptons, JHEP 10, 041 (2019).
  * [13] A. Chodos, R.L. Jaffe, K. Johnson, C.B Thorn and V.F. Weisskopf, New Extended Model of Hadrons, Phys. Rev. D 9, 3471 (1974); A. Chodos, R.L. Jaffe, K. Johnson and C.B Thorn, Baryon Structure in the Bag Theory, Phys. Rev. D 10, 2599 (1974); P. Hasenfratz and J. Kuti, The Quark Bag Model, Phys. Rep. 40, 75 (1978).
  * [14] N. Isgur and G. Karl, P-wave Baryons in the Quark Model, Phys. Rev. D 18, 4187 (1978); N. Isgur and G. Karl, Positive-Parity Excited Baryons in a Quark Model with Hyperfine Interactions, Phys. Rev. D 19, 2653 (1979); N. Isgur, Mesonlike baryons and the spin-orbit puzzle, Phys. Rev. D62, 014025 (2000).
  * [15] C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg, _Photons and Atoms. Introduction to Quantum Electrodynamics_ (WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim, 2004), Ch. I.B.5.
  * [16] Despite being written in terms of gauge fields, the longitudinal momenta are gauge-invariant vectors. For instance, for the EM field $\mathbf{P}_{\parallel}^{A}=\mathbf{E}\times\mathbf{B}$ [15].
  * [17] M. Donaire, Net Force on an Asymmetrically Excited Two-Atom System from Vacuum Fluctuations, Phys. Rev. A 94, 062701 (2016).; J. Sánchez-Cánovas and M. Donaire, Two-Atom van der Waals Forces with One Atom Excited: The Identical-Atom Limit, Phys. Rev. A 104, 052814 (2021).
  * [18] A. de Rújula, H. Georgi and S.L. Glashow, Hadron Masses in a Gauge Theory, Phys. Rev. D 12, 147 (1975); E. Eichten, K. Gottfried,T. Kinoshita, K.D.L ane and T.M. Yan, Charmonium: Comparison with Experiment, Phys. Rev. D 21, 203 (1980); E. Eichten, S. Godfrey, H. Mahlke, J.L. Rosner, Quarkonia and theirTransitions, Rev. Mod. Phys. 80, 1161 (2008); G. Yang, J. Ping and J. Segovia, Tetra- and Penta-Quark Structures in the Constituent Quark Model, Symmetry 12, 1869 (2020).
  * [19] J. Ashman, _et al._ (The European Muon Collaboration), A Measurement of the Spin Asymmetry and Determination of the Structure Function g1 in Deep Inelastic Muon-Proton Scattering, Phys. Lett. B 206, 364 (1988); J. Ashman, _et al._ (The European Muon Collaboration), An Investigation of the Spin Structure of the Proton in Deep Inelastic Scattering of Polarised Muons on Polarised Protons, Nucl. Phys. B 328, 1 (1989); N.J. Abdulameer _et al._ (PHENIX Collaboration), Measurement of Direct-Photon Cross Section and Double-Helicity Asymmetry at $\sqrt{s}=510$ GeV in $\vec{p}+\vec{p}$ Collisions, Phys. Rev. Lett. 130, 251901 (2023); L.E. Gordon and W. Vogelsang, Polarized and Unpolarized Isolated Prompt Photon Production Beyond the Leading Order, Phys. Rev. D 50, 1901 (1994).
  * [20] A. W. Thomas and W. Weise, The Structure of the Nucleon, Wiley-VCH, Berlin, Germany (2001); A.W. Thomas, S. Theberge, and G.A. Miller, The Cloudy Bag Model of the Nucleon, Phys. Rev. D 24, 216 (1981); F. Myhrer and A.W. Thomas, A possible resolution of the proton spin problem, Phys. Lett. B 663, 302 (2008); A. Deur, S.J. Brodsky, G.F. De Téramond, The Spin Structure of the Nucleon, Rept. Prog. Phys. 82, 076201 (2019); R.L. Jaffe and A. Manohar, The g1 Problem: Deep Inelastic Electron Scattering and the Spin of the Proton, Nucl. Phys. B 337, 506 (1990); R.L. Jaffe and H.J. Lipkin, The Spin on the Sea in a Simple Quark Model, Phys. Lett. B 266, 458 (1991); R.L. Jaffe, Gluon Spin in the Nucleon, Phys. Lett. B 365, 359 (1996).
  * [21] R. Ballou, private communication.
  * [22] K. Johnson, _The MIT Bag Model_ , Acta Phys. Pol. B6, 865 (1975).
  * [23] L. Conci and M. Traini, Quark Momentum Distribution in Nucleons, Few-Body Systems , 8, 123 (1990).
  * [24] For consistency, the interaction Hamiltonians with the gauge fields are also to be modified in the nonrelativistic limit. However, since we are only interested in the spatial wavefunction of the nucleon, it suffices to consider this limit upon $H_{n}^{\textrm{HO}}$ only.
  * [25] In Eqs. (35)-(37), the spatial wavefunctions $\psi^{s(^{\prime})}_{1,-1}(\mathbf{x,y})$ refer to those of the individual valence quarks. They are common to all the quarks in the states $|n^{0}_{\uparrow\downarrow}\rangle$, $|p^{0}_{\uparrow\downarrow}\rangle$, but for their spin. The net result of the spin difference on the evaluation of the axial currents is the multiplicative prefactors which account for the combination of up and down quarks in each of the three-quark components of those states. Thus, the remaining factor depends only on the net spin of the initial state, i.e., $|n^{0}_{\uparrow}\rangle$ in our case. Hence, it must be understood that $s=+$ in all those equations.


# The value of the early-time Lieb-Robinson correlation function for qubit
arrays

Brendan J. Mahoney  Craig S. Lent  Department of Electrical Engineering
University of Notre Dame
Notre Dame, IN 46556, USA

(August 28, 2024)

###### Abstract

The Lieb-Robinson correlation function captures propagation of quantum
correlations in many-body system. We calculate the value of the leading order
of the correlation function, not its bound, for a system of interacting qubits
at early times. The general analytical result is compared to numerical
calculations and is applied to regular qubit lattices in one, two, and three
dimensions. We find an exact expression value of the Lieb-Robinson velocity
for near-neighbor coupled systems and see the approximately exponential
leading edge of correlations emerge for large arrays.

Lieb-Robinson, quantum correlations, entanglement

##  I Introduction

Entanglement in many-body quantum systems and the spread of quantum
correlations has been of considerable interest for both fundamental reasons
and for the possible applications in quantum computing. Focus on quantum
information naturally employs entropic measures derived from the von Neumann
entropy of a state. Multiple measures have been proposed to capture the
quantum entanglement of spatially separated systems and this remains an area
of great activity Tsomokos  _et al._ (2007); Horodecki  _et al._ (2009); Amico
_et al._ (2008). Information so conceived is a property of the state of the
system, characterized by either a state vector for pure states or a density
operator for pure or mixed states.

Another avenue, not tied to the state of the system, considers the commutation
relations between local operators that act on spatially separated subsystems.
Preeminent here is investigation of the Lieb-Robinson operator Lieb and
Robinson (1972), which quantifies the quantum correlation between a Heisenberg
operator $\hat{A}_{j}(t)$ which acts on subsystem $j$ of the composite many-
body system, and an operator $\hat{B}_{k}$ which acts on a distinct subsystem
$k$. The operator is given by

$\hat{C}_{A_{j},B_{k}}(t)=[\hat{A}_{j}(t),\hat{B}_{k}]$ (1)

The Lieb-Robinson correlation function is given by the norm of the operator.

${C}_{A_{j},B_{k}}(t)=\left\|\hat{C}_{A_{j},B_{k}}(t)\right\|$ (2)

At $t=0$, the operators $\hat{A}$ and $\hat{B}$, having different support,
commute and the correlation function is zero. We can think of the operator
$\hat{A}_{j}(t)$ as spreading out into the system and as its support comes to
include that of $\hat{B}$ the correlation increases. This captures the quantum
correlation, between the two systems which is the source of entanglement. The
Lieb-Robinson correlation function depends only on the choice of the operators
and the Hamiltonian of the system, independent of any initial state.

Lieb and Robinson established, for quite general conditions, that the
correlation so defined propagates outward with a finite speed and with an
exponentially bounded leading edge in space:

$C_{AB}(t)\leq c\,e^{-a(d(j,k)-v_{\text{\tiny LR}}t)}$ (3)

where $d(j,k)$ is an appropriate distance measure between the two subsystems
and $v_{\text{\tiny LR}}$ is the Lieb-Robinson velocity.

Because of the importance of this result, there has been a substantial effort
to refine the bound in particular circumstances Wang and Hazzard (2020);
Edwards  _et al._ (2010); Else  _et al._ (2020); Prémont-Schwarz  _et al._
(2010); Nachtergaele  _et al._ (2006); Hastings (2004) and to understand the
connection between entropic measures and correlation measures Bravyi  _et al._
(2006); Kastner (2015). Spin models and equivalent qubit arrays have played a
prominent role in studying bounds for specific systems Schneider  _et al._
(2021) and can also be explored experimentally King  _et al._ (2021).

Following the work of Luitz et al. Luitz and Bar Lev (2017); Hémery  _et al._
(2019); Luitz and Bar Lev (2019); Colmenarez and Luitz (2020), we focus here
on a specific class of Hamiltonians and compute the correlation function
itself, rather than a bound on the correlation function. This comes at the
cost of generality, of course, but yields insight into the way correlations
propagate. In Section II we consider a linear array of near-neighbor coupled
qubits—the 1D tranverse field Ising model Schneider  _et al._ (2021), and
calculate the Lieb-Robinson correlation function numerically. This is only
practical for small systems. We then derive an analytical expression for the
correlation function which we can compare with the numerical results. The
analysis is extended to an arbitrary array of qubits in Section III, and then
applied to regular square lattices in two and three dimensions in Section IV.

For the purposes of stating a bound, the norm in Eq. (3) need not be
specified, but to calculate the value of the correlation function a particular
norm must be chosen. Here we use the normalized Frobenius norm

$\left\|\hat{Q}\right\|=\sqrt{\frac{\operatorname{Tr}\left(\hat{Q}^{\dagger}\hat{Q}\right)}{\mathcal{N}}},$
(4)

where $\mathcal{N}$ is the dimension of the Hilbert space. The normalization
ensures that the value is independent of the size of the space.

##  II Lieb Robinson correlation function for a regular linear array

###  II.1 Direct numerical calculation

We consider first a linear array of $N_{q}$ qubits with an energy $\Delta$
coupling adjacent z-components, as described by the following Hamiltonian:

$\hat{H}=-\gamma\sum\limits_{k}^{{N_{q}}}{\hat{\sigma}^{(k)}_{x}-\frac{\Delta}{2}}\sum\limits_{k}^{{N_{q}}-1}{\hat{\sigma}^{(k)}_{z}}\,\hat{\sigma}^{(k+1)}_{z}.$
(5)

Here $\hat{\sigma}^{(k)}_{\alpha}$ is the Pauli spin operator with
$\alpha=(x,y,z)$ operating on the $k^{th}$ qubit. Each operator is understood
to be embedded in the full direct-product Hilbert space of the array with
dimension $\mathcal{N}=2^{N_{q}}$. This system is the 1D transverse field
Ising model. The first term represent the internal dynamics which allow each
qubit to flip its state. The energy $\gamma$ sets the characteristic time for
the dynamics

$\tau\equiv\pi\hbar/\gamma.$ (6)

We have the usual commutation relations for Pauli operators.

$\displaystyle\left[\hat{\sigma}_{x}^{(j)},\hat{\sigma}_{y}^{(k)}\right]$
$\displaystyle=$ $\displaystyle 2\,i\,\hat{\sigma}_{z}^{(k)}\delta_{j,k}$
$\displaystyle\left[\hat{\sigma}_{y}^{(j)},\hat{\sigma}_{z}^{(k)}\right]$
$\displaystyle=$ $\displaystyle 2\,i\,\hat{\sigma}_{x}^{(k)}\delta_{j,k}$
$\displaystyle\left[\hat{\sigma}_{z}^{(j)},\hat{\sigma}_{x}^{(k)}\right]$
$\displaystyle=$ $\displaystyle 2\,i\,\hat{\sigma}_{y}^{(k)}\delta_{j,k}$ (7)

We consider the Lieb-Robinson operator between the $z$ components of the first
qubit and qubit $k$.

$\hat{C}_{k}(t)={[\hat{\sigma}_{1}^{z}(t),\hat{\sigma}_{k}^{z}(0)]}$ (8)

The Lieb-Robinson correlation function is given by the norm of the operator.

$C_{k}(t)=\left\|\hat{C}_{k}(t)\right\|$ (9)

There is no loss of generality here by choosing the first qubit as the
reference because the coupling is the same between any adjacent qubits. The
correlation function $C_{k}$ captures the quantum correlation between two
qubits that are $(k-1)$ apart.

The time-development of the Heisenberg operator $\hat{\sigma}_{1}^{(z)}$ is
given by

$\hat{\sigma}_{z}^{(1)}(t)=e^{i\frac{\hat{H}}{\hbar}t}\hat{\sigma}_{z}^{(1)}e^{-i\frac{\hat{H}}{\hbar}t}$
(10)

The correlation function can then be written

$C_{k}(t)=\left\|\left[e^{i\frac{\hat{H}}{\hbar}t}\hat{\sigma}_{z}^{(1)}e^{-i\frac{\hat{H}}{\hbar}t},\hat{\sigma}_{z}^{(k)}(0)\right]\right\|.$
(11)

Figure 1 shows the Lieb Robinson correlation function for a line of nine
qubits by a direct numerical evaluation of Eq. (11) for the Hamiltonian in Eq.
(5).

Figure 1: The Lieb-Robinson correlation function between the reference qubit 1
and qubit $k$ for a linear array of 9 qubits with near-neighbor interactions.
The calculation is done by numerical evaluation of $C_{q}$ using Eq. (11) and
the Hamiltonian of Eq. (5) with $\Delta/\gamma=1.$

We are interested in the early-time behavior as correlations initially spread
down the line. By “early-time” behavior we mean times such that the
correlation with the reference qubit 1 is small. Our focus is on the leading
order term in time as the correlations between the first site and the $k^{th}$
site initially grow. Figure 2 shows a view of the correlations zoomed in to
the lower left-hand corner of Fig. 1.

Figure 2: Detail of the early-time behavior Lieb-Robinson correlation function
as the quantum correlations move down the line of nine qubits. The calculation
and parameters are as in Figure (1).

The early time behavior of the Lieb-Robinson correlation function has a power-
law dependence, as is evident in Fig. 2 and in the log-log plot shown in Fig.
3. The dots in the figure represent the numerical solution of equation (11)
for $\Delta/\gamma=1$. The lines are the results of the analytical expression
derived below. Figure 4 similarly shows the numerical result for the case when
$\Delta/\gamma=5$. The larger interaction strengths makes the correlations
spread more rapidly down the line.

Figure 3: The Lieb-Robinson early-time correlation function for the nine-qubit
line demonstrates the power-law dependence on time. The qubit index is $k$ and
the dots represent numerical evaluation of Eq. (11) for the Hamiltonian in Eq.
(5) with $\Delta/\gamma=1$. The lines represent the analytic result of Eq.
(30).  Figure 4: The Lieb-Robinson early-time correlation function for the
nine-qubit line with qubit index $k$. The dots represent the direct numerical
evaluation of Eq. (11) with $\Delta/\gamma=5$ and the lines represent the
analytic result of Eq. (30). The stronger coupling between neighboring qubits
yields a faster propagation of correlations down the chain.

###  II.2 Analytic solution for early-time behavior

We adopt the following notation for an $n$-nested commutator.

$\left[(\hat{A})^{n},\hat{B}\right]\equiv\left[\hat{A},\left[\hat{A},\ldots\left[\hat{A}[\hat{A},\hat{B}]\right]\right]\ldots\right]$
(12)

It will also be useful to define a notation for a nested commutator with a
sequence of operators.

$\left[(\hat{\mathcal{O}}_{1},\hat{\mathcal{O}}_{2},\ldots,\hat{\mathcal{O}}_{n}),\hat{B}\right]\equiv\left[\hat{\mathcal{O}}_{1},\left[\hat{\mathcal{O}}_{2},,\ldots\left[\hat{\mathcal{O}}_{n},\hat{B}]\right]\right]...\right]$
(13)

Using the notation of Eq. (12), we can write the time dependence of any
Heisenberg operator as

$\hat{A}(t)=\hat{A}(0)+\sum_{n=1}^{\infty}\frac{1}{n!}\left(\frac{it}{\hbar}\right)^{n}\left[(\hat{H})^{n},\hat{A}\right].$
(14)

Applying this to $\hat{\sigma}_{z}^{(1)}(t)$ we have

$\displaystyle\hat{C}_{k}(t)$ $\displaystyle=$
$\displaystyle\left[\left(\sum_{n=1}^{\infty}\frac{1}{n!}\left(\frac{it}{\hbar}\right)^{n}\left[(\hat{H})^{n},\hat{\sigma}^{(1)}_{z}\right]\right),\hat{\sigma}^{(k)}_{z}\right]$
(15) $\displaystyle=$
$\displaystyle\sum_{n=1}^{\infty}\frac{1}{n!}\left(\frac{it}{\hbar}\right)^{n}\left[\left[(\hat{H})^{n},\hat{\sigma}^{(1)}_{z}\right],\hat{\sigma}^{(k)}_{z}\right]$
(16) $\displaystyle=$
$\displaystyle\sum_{n=1}^{\infty}\hat{C}_{k}^{(n)}\left(\frac{t}{\tau}\right)^{n}$
(17)

The early time Lieb-Robinson operator for qubit $k$ is the lowest-order
$\hat{C}_{k}^{(n)}$ which is not zero. All other terms in Eq. (17) will be
vanishingly small for early enough times.

For convenience, we define

$\hat{G}^{(n)}_{k}=\left[\left[(\hat{H})^{(n)},\hat{\sigma}_{z}^{(1)}\right],\hat{\sigma}_{z}^{(k)}\right],$
(18)

so that,

$\hat{C}^{\,(n)}_{k}=\frac{1}{n!}\frac{\pi^{n}i^{n}}{\gamma^{n}}\hat{G}^{(n)}_{k}.$
(19)

To calculate the self-correlation (bit 1 with itself), we first evaluate

$\displaystyle\left[(\hat{H})^{(1)},\hat{\sigma}_{z}^{(1)}\right]$
$\displaystyle=$
$\displaystyle(-\gamma)[\hat{\sigma}_{x}^{(1)},\hat{\sigma}_{z}^{(1)}]$ (20)
$\displaystyle\;-\frac{\Delta}{2}[\hat{\sigma}_{z}^{(1)}\hat{\sigma}_{z}^{(2)},\hat{\sigma}_{z}^{(1)}]$
$\displaystyle=$ $\displaystyle-(-2i)\gamma\hat{\sigma}_{y}^{(1)}$

so

$\displaystyle\hat{G}^{(1)}_{1}$ $\displaystyle=$
$\displaystyle\left[\left[(\hat{H})^{(1)},\hat{\sigma}_{z}^{(1)}\right],\hat{\sigma}_{z}^{(1)}\right]$
(21) $\displaystyle=$ $\displaystyle(-2i)^{2}\gamma\hat{\sigma}_{x}^{(1)}$

and

$\displaystyle\hat{C}^{\,(1)}_{1}$ $\displaystyle=$
$\displaystyle-i4\pi\hat{\sigma}_{x}^{(1)}$ $\displaystyle C^{\,(1)}_{1}$
$\displaystyle=$ $\displaystyle\|\hat{C}^{\,(n)}_{k}\|=4\pi.$ (22)

The early-time Lieb-Robinson self-correlation is therefore

$C_{1}(t)\approx 4\pi\left(\frac{t}{\tau}\right)$ (23)

where in this case the $n=1$ term is the first nonzero term which dominates at
early times. The linear increase in $C_{1}(t)$ is evident in Figs. (1), (3),
and (4).

For the near-neighbor correlation $C_{2}(t)$, the first non-zero term in Eq.
(18) occurs for $n=3$.

$\hat{G}^{(3)}_{2}=\left[\left[(\hat{H})^{(3)},\hat{\sigma}_{z}^{(1)}\right],\hat{\sigma}_{z}^{(2)}\right],$
(24)

The term $\left[(\hat{H})^{(3)},\hat{\sigma}_{z}^{(1)}\right]$ contains many
cross terms from the three iterated commutators with $\hat{H}$ in Eq. (5). But
only one term survives the commutator with $\hat{\sigma}_{z}^{(2)}$ in Eq.
(24), with the result that

$\hat{G}^{(3)}_{2}=(-2i)^{4}(-\gamma)^{2}\left(\frac{\Delta}{2}\right)\hat{\sigma}_{x}^{(1)}\hat{\sigma}_{x}^{(2)}$
(25)

and the early-time result for the correlation function is

$C_{2}(t)\approx\frac{2^{3}}{3!}\pi^{2}\left(\frac{\Delta}{\gamma}\right)\left(\frac{t}{\tau}\right)^{3}.$
(26)

For $C_{3}$, the $n=5$ term is the leading term. The structure of the only
surviving term for $\hat{G}_{3}^{(5)}$ (Eq. (18)) in this case is
illustrative. We write it using the notation of Eq. (13).

$\hat{G}^{(5)}_{3}\propto\left[\left[\left(\hat{\sigma}_{x}^{(3)},\hat{\sigma}_{z}^{(3)}\hat{\sigma}_{z}^{(2)},\hat{\sigma}_{x}^{(2)},\hat{\sigma}_{z}^{(2)}\hat{\sigma}_{z}^{(1)},\hat{\sigma}_{x}^{(1)}\right),\hat{\sigma}_{z}^{(1)}\right],\hat{\sigma}_{z}^{(3)}\right]$
(27)

The nested commutator has the structure
$(X_{3},\,\,\,Z_{3}Z_{2},X_{2},\,\,\,Z_{2}Z_{1},X_{1})$ traversing the
shortest path (directly) between qubit $3$ and qubit $1$. The alternation
between terms requires 5 nestings of commutators with the the Hamiltonian. A
nesting of this form is necessary to avoid the whole expression becoming zero.
The general feature is that to connect each pair of qubits along the path
requires two nested commutators: a $\hat{\sigma}_{x}^{(k)}$ term followed by a
$\hat{\sigma}_{z}^{(k)}\hat{\sigma}_{z}^{(k+1)}$ term. Therefore, only odd-
order terms contribute.

Each commutator in Eq. (27) with $\hat{\sigma}_{x}$ generates a factor of
$(-2i)\gamma$ and each commutator with the product
$\hat{\sigma}_{z}\hat{\sigma}_{z}$ generates a factor of $(-2i)(-\Delta/2)$ so

$\hat{G}^{(5)}_{3}=(-2i)^{6}\gamma^{4}(-\Delta/2)^{2}\hat{\sigma}_{x}^{(3)}\hat{\sigma}_{x}^{(2)}\hat{\sigma}_{x}^{(1)}.$
(28)

The norm of the Pauli string is 1. Using (19), we then obtain

${C_{3}}(t)\approx\frac{{{2^{4}}{\pi^{5}}}}{{(5)!}}{\left({\frac{\Delta}{\gamma}}\right)^{2}}{\left({\frac{t}{\tau}}\right)^{5}}.$
(29)

In general, the lowest-order nonzero $\hat{G}^{(n)}_{k}$ is proportional to a
Pauli string of $k$ $\hat{\sigma}_{x}$ terms, one for each qubit between 1 and
k, and $n=2k-1$. The result for the early-time Lieb-Robinson correlation
function between the reference qubit 1 and qubit $k$ is given by:

${C_{k}}(t)\approx\frac{{{2^{k+1}}{\pi^{2k-1}}}}{{(2k-1)!}}{\left({\frac{\Delta}{\gamma}}\right)^{k-1}}{\left({\frac{t}{\tau}}\right)^{2k-1}}.$
(30)

Correlations calculated with Eq. (30) yield the solid lines in Figs. 3 and 4
and the dots show the values from numerical solution of Eq. (11). The
numerical calculations include all orders in time, and thus are correct at all
times, not just early times, as illustrated in Fig. 1. Equation (30) shows the
leading order term, which is arbitrarily close to the exact result for early
enough times. The excellent agreement between the two seen in Figs. 3 and 4
shows that the higher order terms are indeed negligible.

###  II.3 Asymptotic behavior

Figures 1 to 4 show the Lieb Robinson correlation $C_{k}(t)$ for individual
qubits with index $k$ in the linear array as a function of time. Figure 5
shows snapshots of the correlation function, computed from Eq. (30), as a
function of $k$ at particular times. This lets us see the “correlation wave”
propagating from the reference bit down the line. Values of $C_{k}$ larger
than $10^{-2}$ are clipped and not shown here because our focus is on early-
time behavior for each qubit. At the times shown in the figure, the
correlation front decays rapidly in space–faster than an exponential.

Figure 5: Snapshots of the early-time Lieb-Robinson correlation at different
times for a line of qubits as a function of qubit index $k$. The correlation
front propagates down the line and, for the times shown, the shape of the
front is attenuated faster than exponentially. The solid lines are the results
of Eq.(30) with coupling $\Delta/\gamma=1$. The dotted line shows the
correlation front at one time evaluated using Eq. (31).

We examine the asymptotic shape of the correlation front propagating down the
line quantitatively by taking the large $k$ limit of Eq. (30). The factorial
can then be replaced by Stirlings approximation to obtain

$C_{k}(t)\approx\sqrt{\frac{2}{\pi}}\sqrt{\frac{\gamma}{\Delta}}\frac{1}{\sqrt{k}}\left(\frac{v_{\text{\tiny
LR}}}{(k-1/2)}\frac{t}{\tau}\right)^{2k-1},$ (31)

where $v_{\text{\tiny LR}}$ is the Lieb-Robinson velocity given by

$v_{\text{\tiny LR}}=e\pi\sqrt{\frac{\Delta}{2\gamma}}.$ (32)

This velocity is expressed in terms of the natural dimensionless time
$t/\tau$. In dimensional form it would be $v_{\text{\tiny
LR}}=(e/\sqrt{2})\sqrt{\Delta\gamma}/\hbar$ qubits/s.

Figure 6: Velocity saturation to the Lieb-Robinson velocity for the one-
dimensional qubit line. The finite-difference velocity defined by Eq. (33) is
calculated from Eq. (30) for for $\Delta/\gamma=1$ and
$C_{\text{thresh}}=10^{-25}$. The values marked by an $\times$ are the results
of numerical solution of Eq. (11)

We can examine how the propagation of the correlation front converges to a
constant speed. Let $C_{\text{thresh}}$ be a threshold value of the
correlation function and let $t_{k}/\tau$ be the time at which qubit $k$
crosses that threshold of correlation, so
$C_{k}(t_{k}/\tau)=C_{\text{thresh}}$. We define the backwards finite-
difference velocity of the propagating front by

$v_{k}=\frac{1}{t_{k}/\tau-t_{k-1}/\tau}.$ (33)

Figure 6 shows this velocity $v_{k}$ as a function of qubit index down the
line. The velocity saturates to the value $v_{\text{\tiny LR}}$ given by Eq.
(32). The line shows values calculated from Eq. (30) for $\Delta/\gamma=1$ and
$C_{\text{thresh}}=10^{-25}$. The saturated Lieb-Robinson velocity is then
$e\pi/\sqrt{2}\approx 6.04$ in terms of the dimensionless time $t/\tau$.
Figure 6 also shows the values of $v_{k}$ calculated numerically directly from
the defining Eq. (11).

Figure 7: Snapshots of the Lieb-Robinson correlation as in Fig. 5, but now for
much later times and much farther down the qubit line. Snapshots are shown for
even values of $t/\tau$ between $1360$ and $1400$. Note that the qubit index
$k$ is offset by $10^{4}$. By this point the correlation front, calculated
from Eq. (30), is well-approximated by the exponential dependence of Eq. (34)
and propagates at the saturated Lieb-Robinson velocity given by Eq. (32). The
dotted line shows the correlation front at one time evaluated in the limit of
very large $k$ using Eq. (34).

Figure 7 shows the correlation front much farther down the line (qubits with
index more than $10,250$) and at considerably longer times ($t/\tau$ more than
$1360$). By this point, the shape of the leading edge of the front is
approximately exponential and moving at a steady speed. Taken together, Figs.
5 and 7 also illustrate that “early-time” does not refer to small values of
$t/\tau$, but rather to the span of time before the quantum correlations
between qubit $k$ and the reference qubit 1 become large.

In the region of the advancing correlation front for very large values of $k$,
Eq. (31) can be written:

$C_{k}(t)\approx
e\sqrt{\frac{2}{\pi}}\sqrt{\frac{\gamma}{\Delta}}\frac{1}{\sqrt{k}}e^{-2\left(k-v_{\text{\tiny
LR}}\left(\frac{t}{\tau}\right)\right)}$ (34)

This matches the slope of the lines and the spacing in Fig. 7.

We note that the magnitude of early-time correlations become smaller as the
front propagates down the chain because of the $k^{-1/2}$ dependence in
Eq.(34). This also means the shape of the front differs slightly from a simple
$e^{-2kx}$ in each snapshot shown in Fig. 7. The attenuation of the Lieb-
Robinson correlation function was also noted in Ref. Colmenarez and Luitz
(2020) for the Heisenberg XXX model with short-range interactions.

##  III Correlation function for an arbitrary qubit array

The result in Eq. (30) can be generalize for any network with arbitrary
interaction strengths $\Delta_{j,k}$ between qubits $j$ and $k$. The
Hamiltonian for the network is

$\hat{H}=-\gamma\sum\limits_{k}^{{N_{q}}}{\hat{\sigma}_{{}_{k}}^{x}-}\frac{1}{2}\sum\limits_{j,k>j}^{{N_{q}}}{{\Delta_{j,k}}\;\hat{\sigma}_{z}^{(j)}}\,\hat{\sigma}_{z}^{(k)}.$
(35)

Figure 8 shows a network of nine qubits with the strength of the coupling
$\Delta_{j,k}$ between each indicated. The dots in Figure 9 show the result of
numerical calculation of the Lieb-Robinson correlation function between the
first qubit and the $k^{th}$ qubit in the array using

$C_{1,k}(t)=\left\|\left[e^{i\frac{\hat{H}}{\hbar}t}\hat{\sigma}_{z}^{(1)}e^{-i\frac{\hat{H}}{\hbar}t},\hat{\sigma}_{z}^{(k)}(0)\right]\right\|.$
(36)

The solid lines are the result of the analytic expression derived below.

Figure 8: An arbitrary network of qubits with different coupling parameters
$\Delta$.  Figure 9: The Lieb-Robinson correlation function for the qubit
network shown in Fig. 8. $C_{1,k}$ is the correlation between qubit 1 and
qubit $k$. The dots are the result of numerical evaluation of Eq. (36) and the
lines represent the general analytic result of Eq. (38). Results for some
values of $k$ nearly overlap: $k=$2, 3; $k=$ 4, 5, 6; $k=$7, 8.

For the linear array studied in the previous section, there was one shortest
path between qubit $1$ and qubit $k$, the direct path along the line, that led
to the dominant term in $\hat{G}$ (Eq. 18). This was seen, for example, in the
calculation for $G^{(5)}_{3}$ in Eq. (27), for which the path was simply
$\\{1,2,3\\}$. We must now generalize this analysis to include more than one
path through the network. A path is then a list of indices that begin at qubit
$j$ and end at qubit $k$. The length of a path is the number of edges that
connect successive qubits.

For an arbitrary array of interacting qubits, there can be several paths with
the same minimum length $L$ that contribute to the leading order early-time
correlation function. We denote each minimum-length path by $P_{\alpha}(j,k)$
where $\alpha$ is the index of the particular path. The index of the $m^{th}$
qubit along path $P_{\alpha}(j,k)$ is denoted $q^{\alpha}_{m}$. The lowest
order nonzero operator

$\hat{G}_{j,k}^{(2L+1)}=\left[\left[(\hat{H})^{(2L+1)},\hat{\sigma}_{z}^{(j)}\right],\hat{\sigma}_{z}^{(k)}\right],$
(37)

then becomes a weighted sum of Pauli strings of $\hat{\sigma}_{x}$ operators
for each qubit traversed in the path. The weights are given by the products of
the couplings $\Delta_{q^{\alpha}_{m},q^{\alpha}_{m+1}}$ along the path. We
then obtain for the early-time Lieb-Robinson correlation function between
qubits $j$ and $k$

$C_{j,k}(t)\approx\frac{2^{L+2}\pi^{2L+1}}{(2L+1)!}\left(\frac{t}{\tau}\right)^{2L+1}\\\
\times\sqrt{\sum_{P_{\alpha}(j,k)}\left(\prod_{m=1}^{L}\frac{\Delta_{q^{\alpha}_{m},q^{\alpha}_{m+1}}}{\gamma}\right)^{2}}$
(38)

where the sum is over all paths with the same minimum length $L$. Equation
(38) makes no assumption about the locality of the coupling $\Delta$ between
qubits, which could extend far beyond near-neighbors. Equation (38) reduces to
Eq. (30) when all the interaction strengths $\Delta$ are the same and connect
only near neighbors. In that case, there is just the single minimum path
between qubit $1$ and qubit $k$ which has length $L=k-1$. The results shown in
Fig. 9 are in excellent agreement with the numerical calculations.

##  IV Correlation function for 2D and 3D regular lattices

Figure 10: Isocontours of the Lieb-Robinson correlation function for a square
two-dimensional array of qubits with near-neighbor coupling. Contours are
shown for correlation values of $10^{-60}$ at several values of time.

We consider a two-dimensional square lattice of qubits with uniform near-
neighbor coupling.

$\hat{H}=-\gamma{\sum}\limits_{k_{x},k_{y}}\hat{\sigma}^{(k_{x},k_{y})}_{x}\\\
-\left(\frac{\Delta}{2}\right){\sum}\limits_{\overset{\text{near
neighbors}}{\overset{k_{x},k_{y}}{k^{\prime}_{x},k^{\prime}_{y}}}}\hat{\sigma}^{(k_{x},k_{y})}_{z}\,\hat{\sigma}^{(k^{\prime}_{x},k^{\prime}_{y})}_{z}$
(39)

The qubit indices $k_{x}$ and $k_{y}$ go from from $-N$ to $+N$. The reference
qubit is at the origin $(k_{x},k_{y})=(0,0)$. Clearly there is no distinction
between positive and negative directions, so the value of the Lieb-Robinson
correlation function at a particular qubit depends only on
$(n,m)=(|k_{x}|,|k_{y}|)$. The minimum path length from the origin to $(n,m)$
is simply $m+n$, and the number of paths with that length is $(n+m)!/(n!m!)$.
Applying the general result of Eq. (38), we obtain for the correlation between
the reference qubit at the origin and the $(n,m)$ qubit (in any quadrant):

$C_{n,m}(t)\approx\frac{2^{n+m+2}\;\pi^{2(n+m)+1}}{(2(n+m)+1)!}\sqrt{\frac{(n+m)!}{n!m!}}\\\
\times\left(\frac{\Delta}{\gamma}\right)^{n+m}\left(\frac{t}{\tau}\right)^{2(n+m)+1}.$
(40)

Figure 10 shown snapshots of isocontours of $C$ evaluated using Eq. (40) for
several values of $t/\tau$. The Lieb-Robinson correlation front expands
outward from the origin. All the minimum-length paths are Manhattan paths,
with the most direct route for expanding correlations along the coordinate
axes. The correlation function along coordinate axes reduces to the one-
dimensional result of Eq. (30). (Note that the reference qubit is here indexed
$(0,0)$ rather than 1.)

As with the one-dimensional case, the speed of the correlation front motion
saturates to the Lieb-Robinson velocity, but for for the square 2D lattice the
velocity depends on direction. We can express this in terms of the angle
$\theta$ with the $x$-axis and find:

$v_{\text{\tiny
LR}}(\theta)=e\pi\left(\sqrt{\frac{\Delta}{2\gamma}}\right)\left(\frac{1+\tan(\theta)}{\tan(\theta)^{\frac{\tan(\theta)}{\tan(\theta)+1}}}\right)^{\frac{1}{4}}\\\
\times\left(\frac{\sqrt{1+\tan(\theta)^{2}}}{1+\tan(\theta)}\right).$ (41)

This expression is valid for $0\leq\theta\leq\pi/4$; other angles can be
deduced by symmetry.

Figure 11: Isosurfaces of the early-time Lieb-Robinson correlation function
for a three-dimensional regular square lattice of qubits. The surfaces are
shown for near-neighbor coupling $\Delta/\gamma=1$ and $C=1\times 10^{-170}$
at times (a) $t/\tau=1$, (b) $t/\tau=5$, and (c) $t/\tau=11$.

The Lieb-Robinson correlation for a three dimensional square lattice with
uniform near-neighbor coupling can similarly obtained from the general result
of Eq. (38). The correlation function between a qubit at the origin and one
with coordinates $(n,m,p)=(|k_{x}|,|k_{y}|,|k_{z}|)$ (in any octant) can be
written

$C_{n,m,p}(t)\approx\frac{2^{n+m+p+2}\;\pi^{2(n+m+p)+1}}{(2(n+m+p)+1)!}\\\
\times\sqrt{\frac{(n+m+p)!}{m!n!p!}}\left(\frac{\Delta}{\gamma}\right)^{(n+m+p)}\left(\frac{t}{\tau}\right)^{2(n+m+p)+1}.$
(42)

Figure 11 shows snapshots of the isosurfaces of the Lieb-Robinson correlation
function $C$ for three successive times. The influence of the Manhattan paths
is clear in the rounded octahedral shape.

The direction-dependent Lieb-Robinson velocity for the 3D square lattice is
expressed in terms of the polar angle $\phi$ and azimuthal angle $\theta$:

$v_{\text{\tiny
LR}}(\theta,\phi)=e\pi\left(\sqrt{\frac{\Delta}{2\gamma}}\right)\left(\frac{P(\theta,\phi)}{Q(\theta,\phi)}\right)^{\frac{1}{4}}R(\theta,\phi)$
(43)

where

$\displaystyle
P(\theta,\phi)=1+\tan(\theta)+\sqrt{\frac{1+\tan(\theta)^{2}}{\tan(\phi)^{2}}}$
(44)
$Q(\theta,\phi)=\tan(\theta)^{\left(\frac{\tan(\theta)}{1+\tan(\theta)+\sqrt{\frac{1+\tan(\theta)^{2}}{\tan(\phi)^{2}}}}\right)}\\\
\times\left(\sqrt{\frac{1+\tan(\theta)^{2}}{\tan(\phi)^{2}}}\right)^{\frac{\sqrt{\frac{1+\tan(\theta)^{2}}{\tan(\phi)^{2}}}}{1+\tan(\theta)+\sqrt{\frac{1+\tan(\theta)^{2}}{\tan(\phi)^{2}}}}}$
(45)

and

$R(\theta,\phi)=\frac{\sqrt{1+\tan^{2}(\theta)+\tan^{2}(\phi)+\tan^{2}(\theta)\tan^{2}(\phi)}}{\tan(\phi)+\tan(\phi)\tan(\theta)+\sqrt{1+\tan(\theta)^{2}}}$
(46)

##  V Discussion

Calculating the value of the Lieb-Robinson correlation function itself, rather
than a bound, for a specific Hamiltonian illuminates some aspects of the way
quantum information propagates. We are able to see how the Lieb-Robinson
velocity emerges as the saturation velocity of the correlation fronts
propagation. A specific value for this velocity, and its angular dependence,
is also obtained for regular lattices. We also see that the leading edge of
the correlation front is super-exponential as the propagation starts, and
acquires the exponential dependence in the bound of Eq. (3) only later.

Equation (38) is the main result of this work and applies to a large class of
ZZ-coupled Hamiltionians given by Eq. (35). We have specialized this result to
address the case of regular square qubit arrays in 1D, 2D and 3D, with near-
neighbor couplings.

There is some competition between the strength of coupling between qubits and
the number of paths through intermediate qubits. But for early times the cost
of adding path length is very high because of the $2L+1$ exponent in the time
exponent of Eq. (38), and because of the factorial in the denominator. This
latter factor also results in the decay of correlations in space as the front
propagates outward. This is reflected the $k^{-1/2}$ in Eq. (31) for the
Ising-coupled qubit line.

## References

  * Tsomokos  _et al._ (2007) D. I. Tsomokos, M. J. Hartmann, S. F. Huelga, and M. B. Plenio, Entanglement dynamics in chains of qubits withnoise and disorder, New J. of Phys. 9, 79 (2007).
  * Horodecki  _et al._ (2009) R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Quantum entanglement, Rev. Mod. Phys. 81, 865 (2009).
  * Amico  _et al._ (2008) L. Amico, R. Fazio, A. Osterloh, and V. Vedral, Entanglement in many-body systems, Rev. Mod. Phys. 80, 517 (2008).
  * Lieb and Robinson (1972) E. Lieb and D. Robinson, The finite group velocity of quantum spin systems, Commun. Math. Phys . 28, 251 (1972).
  * Wang and Hazzard (2020) Z. Wang and K. R. Hazzard, Tightening the lieb-robinson bound in locally interacting systems, PRX Quantum 1, 010303 (2020).
  * Edwards  _et al._ (2010) E. E. Edwards, S. Korenblit, K. Kim, R. Islam, M.-S. Chang, J. K. Freericks, G.-D. Lin, L.-M. Duan, and C. Monroe, Quantum simulation and phase diagram of the transverse-field ising model with three atomic spins, Phys. Rev. B 82, 060412 (2010).
  * Else  _et al._ (2020) D. V. Else, F. Machado, C. Nayak, and N. Y. Yao, Improved lieb-robinson bound for many-body hamiltonians with power-law interactions, Phys. Rev. A 101, 022333 (2020).
  * Prémont-Schwarz  _et al._ (2010) I. Prémont-Schwarz, A. Hamma, I. Klich, and F. Markopoulou-Kalamara, Lieb-robinson bounds for commutator-bounded operators, Phys. Rev. A 81, 040102 (2010).
  * Nachtergaele  _et al._ (2006) B. Nachtergaele, Y. Ogata, and R. Sims, Propagation of correlations in quantum lattice systems, J. Stat. Phys. 124, 1 (2006).
  * Hastings (2004) M. B. Hastings, Locality in quantum and markov dynamics on lattices and networks, Phys. Rev. Lett. 93, 140402 (2004).
  * Bravyi  _et al._ (2006) S. Bravyi, M. B. Hastings, and F. Verstraete, Lieb-robinson bounds and the generation of correlations and topological quantum order, Phys. Rev. Lett. 97, 050401 (2006).
  * Kastner (2015) M. Kastner, Entanglement-enhanced spreading of correlations, New. J. Phys. 17, 123024 (2015).
  * Schneider  _et al._ (2021) J. T. Schneider, J. Despres, S. J. Thomson, L. Tagliacozzo, and L. Sanchez-Palencia, Spreading of correlations and entanglement in the long-range transverse ising chain, Phys. Rev. Research 3, L012022 (2021).
  * King  _et al._ (2021) A. D. King, C. Nisoli, E. D. Dahl, G. Poulin-Lamarre, and A. Lopez-Bezanilla, Qubit spin ice, Science 373, 576 (2021), https://www.science.org/doi/pdf/10.1126/science.abe2824 .
  * Luitz and Bar Lev (2017) D. J. Luitz and Y. Bar Lev, Information propagation in isolated quantum systems, Phys. Rev. B 96, 020406 (2017).
  * Hémery  _et al._ (2019) K. Hémery, F. Pollmann, and D. J. Luitz, Matrix product states approaches to operator spreading in ergodic quantum systems, Phys. Rev. B 100, 104303 (2019).
  * Luitz and Bar Lev (2019) D. J. Luitz and Y. Bar Lev, Emergent locality in systems with power-law interactions, Phys. Rev. A 99, 010105 (2019).
  * Colmenarez and Luitz (2020) L. Colmenarez and D. J. Luitz, Lieb-robinson bounds and out-of-time order correlators in a long-range spin chain, Phys. Rev. Research 2, 043047 (2020).