# Combined Effects of $f(R)$ Gravity and Massive Neutrinos on the Turn-Around
Radii of Dark Matter Halos

Jounghun Lee11affiliation: Astronomy Program, Department of Physics and
Astronomy, FPRD, Seoul National University, Seoul 08826, Korea  , Marco
Baldi22affiliation: Dipartimento di Fisica e Astronomia, Alma Mater Studiorum
Università di Bologna, viale Berti Pichat, 6/2, I-40127 Bologna, Italy
33affiliation: INAF - Osservatorio Astronomico di Bologna, via Ranzani 1,
I-40127 Bologna, Italy 44affiliation: INFN - Sezione di Bologna, viale Berti
Pichat 6/2, I-40127 Bologna, Italy<EMAIL_ADDRESS>
###### Abstract

We present a new statistics based on the turn-around radii of cluster halos to
break the dark sector degeneracy between the $\Lambda{\rm CDM}$ model and the
alternative ones with $f(R)$ gravity and massive neutrinos ($\nu$)
characterized by the strength of the fifth force, $|f_{R0}|$, and the total
neutrino mass, $M_{\nu}$. Analyzing the rockstar halo catalogs at the present
epoch from the DUSTGRAIN-pathfinder $N$-body simulations performed for four
different cosmologies, namely, $\Lambda{\rm CDM}$ ($|f_{R0}|=0$, $\sum
m_{\nu}=0.0\,{\rm eV}$), fR6 ($|f_{R0}|=10^{-6}$, $\sum m_{\nu}=0.0\,{\rm
eV}$), fR6+$0.06\,{\rm eV}$ ($|f_{R0}|=10^{-6}$, $\sum m_{\nu}=0.06\,{\rm
eV}$) and fR5+$0.15\,{\rm eV}$ ($|f_{R0}|=10^{-5}$, $\sum m_{\nu}=0.15\,{\rm
eV}$), which are known to yield very similar conventional statistics to one
another. For each model, we select those cluster halos which do not neighbor
any other larger halos in their bound zones and construct their bound-zone
peculiar velocity profiles at $z=0$. Then, we determine the radial distance of
each selected halo at which the bound-zone velocity becomes equal to the
recession speed of the Hubble flow as its turn around radius, and evaluate the
cumulative probability distribution of the ratios of the turn-around radii to
the virial counterparts, $P(r_{t}/r_{v}\geq\alpha)$. The degeneracy between
the fR6 and fR5+$0.15\,{\rm eV}$ models is found to be readily broken by the
$10\sigma_{\Delta P}$ difference in the value of $P(\alpha=4)$, while the
$3.2\sigma_{\Delta P}$ difference between the $\Lambda{\rm CDM}$ and
fR6+$0.06\,{\rm eV}$ models is detected in the value of $P(\alpha=8.5)$. It is
also found that the four models yield smaller differences in $P(\alpha)$ at
higher redshifts.

Unified Astronomy Thesaurus concepts: Large-scale structure of the universe
(902); Cosmological models (337)

##  1 Introduction

The turn-around radius of a dark matter (DM) halo is a characteristic distance
scale at which the velocity field around the halo has a vanishingly small
value in the radial direction due to the complete counter-balance between its
inward gravity and the outward repulsion of the Hubble flow. Even though the
turn-around radius is a property of a highly nonlinear structure, its value
can in principle be theoretically predictable from the first principles as far
as the halo forms through the spherically symmetric gravitational collapse
process (Pavlidou & Tomaras, 2014; Pavlidou et al., 2014). This advantageous
aspect of the turn-around radius has motivated many authors to examine its
potential as a probe of cosmology. For example, Pavlidou & Tomaras (2014)
analytically evaluated the upper limit on the turn-around radii for the
standard $\Lambda{\rm CDM}$ cosmology, where the gravitational law is
described by Einstein’s general relativity (GR), the present acceleration of
the universe is driven by the cosmological constant ($\Lambda$) with equation
of state $w=-1$, and the most dominant matter content is the collisionless
cold DM (CDM) particles having negligibly low speed at the moment of their
decoupling.

What Pavlidou & Tomaras (2014) proved was that the spherical upper limit on
the turn-around radii sensitively depends on the amount of $\Lambda$ (see also
Pavlidou et al., 2014; Bhattacharya & Tomaras, 2017) and thus that a bound
violation, if observed to occur, could in principle challenge the $\Lambda{\rm
CDM}$ cosmology. Here, a bound violation is a term coined by Pavlidou &
Tomaras (2014) to describe an event of observing a cosmic structure whose
turn-around radius exceeding the analytically found spherical upper limit of
the $\Lambda{\rm CDM}$ cosmology. Later, Lopes et al. (2018) theoretically
proved that the upper limit of the turn-around radii can be used to detect the
presence of modified gravity (MG, Clifton et al., 2012, for a review) which
has an effect of significantly increasing the turn-around radii (see also
Lopes et al., 2019).

The aforementioned theoretical works were based on the simple top-hat
spherical dynamics, from which the real gravitational dynamics is in fact well
known to depart (e.g., Bond & Myers, 1996). To take into account the non-
spherical nature of gravitational collapse for the determination of the turn-
around radii and their upper limit, the numerical experiments had to be
employed (Pavlidou & Tomaras, 2014). For instance, Lee & Yepes (2016) used a
high-resolution N-body simulation to measure the turn-around radii of DM halos
located in the cosmic web (Bond et al., 1996) and demonstrated that the
anisotropic merging along the filamentary structures has an effect of
enlarging the turn-around radii. Their result implied that the occurrence of a
bound violation is not unconditionally prohibited but occasionally possible
even in the $\Lambda{\rm CDM}$ cosmology since the numerically determined non-
spherical upper limit on the turn-around radii turned out to be higher than
the analytically predicted spherical limit (see also Bhattacharya & Tomaras,
2021; Faraoni, 2021; Giusti & Faraoni, 2021).

Nevertheless, the usefulness of the turn-around radii as a cosmological probe
is not necessarily undermined by the fact that their upper limit cannot be
treated in a purely analytical way. Lee & Li (2017) numerically found that it
becomes significantly more probable for a bound violation to occur in the
presence of MG and thus that the frequency of the occurrence of the bound
violations should be a powerful test of GR. Their claim was supported by
several follow-up works which theoretically proved that the alternative
cosmologies including quintessence dark energy (DE), scalar-tensor theory, and
phantom brane world induce more frequent occurrence of the bound violations
(Bhattacharya & Tomaras, 2017; Bhattacharya & Kousvos, 2017; Nojiri et al.,
2018; Lopes et al., 2018, 2019).

In light of the aforementioned works which disclosed the sensitivity of the
turn-around radii especially to the nature of gravity, we attempt here to
numerically explore if the turn-around radii is capable of discriminating the
alternative MG models that have been known to be degenerate with the
$\Lambda{\rm CDM}$ cosmology by the conventional statistics such as the linear
and nonlinear density power spectra, cluster mass function, halo bias factor
and redshift distortion effect (Baldi et al., 2014; Hagstotz et al., 2019a).
For this exploration, our analysis will focus on a particular class of MG
models, namely, the $\nu{\rm CDM}$+$f(R)$ gravity model, where the massive
neutrinos ($\nu$) with non-zero total mass $\sum m_{\nu}$ is present along
with CDM and the apparent acceleration of space time at the present epoch is
caused by the failure of GR on the cosmological scales (see De Felice &
Tsujikawa, 2010, for a review).

The gravitational dynamics of this alternative model is dictated by the
modified Einstein-Hilbert action in which some specified function, $f(R)$,
substitutes for the Ricci scalar $R$. An additional fifth force is generated
by its extra degree of freedom, $f_{R}\equiv df/dR$, dubbed the scalaron,
whose absolute value at the present epoch, $|f_{R0}|$, quantifies how strong
fifth force the $f(R)$ gravity can exert (see Hu & Sawicki, 2007, and
references therein). Although the Chameleon shielding mechanism turns off the
fifth force in the high-density regions, the overall effect of $f(R)$ gravity
alone is to enhance the density growth via its fifth force compared with the
$\Lambda{\rm CDM}$ case (see Khoury & Weltman, 2004, and references therein).
However, in the presence of massive neutrinos which has an effect of
suppressing the density growth (e.g., see Lesgourgues & Pastor, 2014), this
effect of $f(R)$ gravity can be severely attenuated. It was indeed numerically
shown that a proper combination of $\sum m_{\nu}$ with $|f_{R0}|$ can make a
$\nu{\rm CDM}$+$f(R)$ gravity model to yield very similar conventional
statistics to the $\Lambda{\rm CDM}$ case (Baldi et al., 2014; Hagstotz et
al., 2019a).

In this Paper, we will provide a numerical evidence supporting that the turn-
around radii may break this degeneracy between the $\Lambda{\rm CDM}$ and
$\nu{\rm CDM}$+$f(R)$ gravity models. The organization of this paper is as
follows. In Section 2, we briefly review the previously developed algorithm
for the estimation of the turn-around radii of DM halos. In Section 3, we
describe the numerical data used for our analysis and explain how well the
cumulative probability of the turn-around radii of DM halos differentiate
between the $\Lambda{\rm CDM}$ and $\nu{\rm CDM}$+$f(R)$ gravity models. In
Section 4, we summarize the results and discuss the caveats and limitations of
our statistics as a cosmological discriminator.

##  2 A Review of the TRE algorithm

The neighborhood around a DM halo is often divided into three distinct sectors
called the infall, bound and Hubble zones, depending on which effect is more
dominant between the gravity and the cosmic expansion. The infall (Hubble)
zone corresponds to the radial distance range, $r\leq 2r_{v}$ ($r\geq
10r_{v}$), in which the effect of the gravitational attraction of the halo on
the radial components of the peculiar velociteis, $v_{r}$, completely
surpasses (surrenders) that of the receding Hubble flow, where $r_{v}$ is the
halo virial radius. Meanwhile, the in-between bound zone corresponds to the
region where the two competing forces are so well balanced that $v_{r}$ can be
tractable in the linear perturbation theory.

Falco et al. (2014) showed that the profile of the radial components of the
peculiar velocity field in the bound zone, $v_{r}(r)$, around the cluster
halos with virial mass $M_{v}\gtrapprox 0.5\times 10^{14}\,h^{-1}\\!M_{\odot}$
has a universal shape, well approximated by the following formula,

$\frac{v_{r}(r)}{V_{c}}=-A\left(\frac{r_{v}}{r}\right)^{n}\,,$ (1)

where two adjustable parameters, $A$ and $n$, quantify the amplitude and slope
of the profile, respectively, and $V_{c}$ is the circular velocity equivalent
to $(GM_{v}/r_{v})^{1/2}$. The negative sign in the right-hand side of
Equation (1) indicates that the bound-zone neighbors still feel the net
gravitational force of the halo. From here on, we let $v_{r}(r)$ exclusively
denote the profile of the radial components of the peculiar velocity field in
the bound zone around the cluster halos and call it the bound-zone velocity
profile.

Falco et al. (2014) claimed the universality of Equation (1) based on their
numerical finding that the stacked bound-zone velocity profiles over the
cluster halos has a constant slope and amplitude, being almost independent of
the cluster masses and redshifts. It was also found by Falco et al. (2014)
that not only the stacked ones but also the bound-zone velocity profiles,
$v_{r}(r)$, around individual cluster halos follow well the above formula,
although the best-fit values of $A$ and $n$ exhibited substantial scatters
around the mean values. The pioneering work of Falco et al. (2014) motivated
further numerical investigations of the bound-zone velocity profiles around
the cluster halos, which all confirmed the validity of Equation (1) for the
description of $v_{r}(r)$ (Lee et al., 2015; Lee, 2016; Lee & Yepes, 2016;
Lee, 2018; Hansen et al., 2020).

Lee & Yepes (2016) demonstrated with the help of a N-body simulations that the
best-fit values of $A$ and $n$ in Equation (1) depend on the halo environments
and that the best-agreements between the numerically obtained $v_{r}(r)$ and
Equation (1) can be achieved for the case that the cluster halos are located
in the relatively low-density environments, having no larger neighbor halos in
their bound zones. They also showed that even when $v_{r}(r)$ is constructed
not directly from the DM particles but only from the distinct neighbor halos
in the bound zones, it is well described by Equation (1), proving the
feasibility of the observational application of Equation (1) to real data. Lee
(2016) confirmed the universality of Equation (1), showing that the best-fit
values of $A$ and $n$ in Equation (1) are quite insensitive to the variation
of the key cosmological parameters, $\sigma_{8}$ and $\Omega_{m}$. It was also
found by Lee (2016) that Equation (1) is valid even on the lower mass scales
corresponding to the group-size halos, $5\times 10^{12}\leq
M_{v}/(\,h^{-1}\\!M_{\odot})\leq 10^{13}$, as far as the halos are located in
the isolated regions.

Lee et al. (2015) proposed an algorithm based on Equation (1) to estimate the
turn-around radii, $r_{t}$, of the cluster halos, calling it the turn-around
radius estimator (TRE). By definition, the magnitude of the bound-zone
velocity, $v_{r}$, at $r_{t}$ becomes equal to the speed of the Hubble flow,
$H_{0}r_{t}$. By equation (1), however, $|v_{r}(r_{t})|$ is nothing but
$A\,\left[r_{v}/r_{t}\right]^{n}$. Given that $A$ and $n$ have constant
values, the turn-around radius of a cluster halo can be estimated simply by
solving the following equation.

$A\,\left[\frac{r_{v}}{r_{t}}\right]^{n}\ =\frac{H_{\rm 0}r_{t}}{V_{c}}\,.$
(2)

Lee & Li (2017) applied the TRE to the numerical data from a high resolution
N-body simulation and found that the TRE worked better when they used the
best-fit values of $A$ and $n$ obtained separately for each cluster halo
rather than using their constant mean values.

For the application of the TRE to the real clusters from observations, the
critical issue to address was whether or not the values of $A$ and $n$ for the
individual clusters could be obtained without measuring the bound-zone
peculiar velocities. According to Falco et al. (2014), for those cluster halos
embedded in cosmic filaments, it is possible to construct $v_{r}(r)$ from
limited information only on the redshift space positions of the bound-zone
galaxies. Once $v_{r}(r)$ is constructed for the filament clusters whose
virial mass and radius are known, then the values of $A$ and $n$ can be
readily obtained through fitting of $v_{r}(r)$ to Equation (1). Lee (2018)
applied this TRE to the local galaxy clusters located in the straight
filamentary structures and successfully estimated their turn-around radii,
validating its practical usefulness. In Section 3, we will apply the TRE to
the numerical data for the investigation of the combined effects of $f(R)$
gravity and massive neutrinos on the turn-around radii of the cluster halos.

##  3 Physical Analysis

The DUSTGRAIN-pathfinder simulation project aimed at keeping track of
$768^{3}$ CDM particles of individual mass $8.1\times
10^{10}\,h^{-1}\\!M_{\odot}$ under the influence of $f(R)$ gravity in the
presence of massive neutrinos (Giocoli et al., 2019) on the periodic box of
volume $750^{3}\,h^{-3}$Mpc3. The MG-GADGET encoded by Puchwein et al. (2013)
was implemented for the computation of the Hu-Sawicki $f(R)$ gravity (Hu &
Sawicki, 2007), while the incorporation of massive neutrinos was achieved via
the particle-based routine programmed by Viel et al. (2010). See Giocoli et
al. (2019) and Puchwein et al. (2013) for detailed information on the
DUSTGRAIN-pathfinder project and MG-GADGET code, respectively. Among many
$\nu$CDM+$f(R)$ gravity models simulated by the DUSTGRAIN-pathfinder, we
consider the following three, fR6 ($|f_{R0}|=10^{-6}$, $\sum m_{\nu}=0.0\,{\rm
eV}$), ${\rm fR6+0.06}\,{\rm eV}$ ($|f_{R0}|=10^{-6}$, $\sum
m_{\nu}=0.06\,{\rm eV}$), and ${\rm fR5+0.15}\,{\rm eV}$ with
$|f_{R0}|=10^{-5}$ and $\sum m_{\nu}=0.15\,{\rm eV}$) models. As mentioned in
Section 1, it was shown by the previous works of Baldi et al. (2014) and
Hagstotz et al. (2019a) that the conventional statistics can hardly
discriminate these three models from the $\Lambda$CDM cosmology, which was
also simulated by the DUSTGRAIN-pathfinder project, setting the initial
conditions at the Planck values (Planck Collaboration et al., 2016).

For each of the four models, i.e., the $\Lambda{\rm CDM}$ and the three
$\nu{\rm CDM}$+$f(R)$ gravity cosmologies, Lee et al. (2022) identified the DM
halos by applying the Rockstar algorithm (Behroozi et al., 2013) to the
snapshot data of the DUSTGRAIN-pathfinder simulations. From the Rockstar halo
catalogs at $z=0$, we extract the cluster-size distinct halos with $M_{v}\geq
4.05\times 10^{13}\,h^{-1}\\!M_{\odot}$ enclosing $500$ or more DM particles
within their virial radii $r_{v}$. Among them, we select only those cluster
halos which do not neighbor any higher-mass halos in their bound zones. From
here on, the selected cluster halos will be referred to as the host halos, for
which we find the neighbor halos located in their bound-zones and containing
$20$ or more DM particles. Then, we compute the radial components of the
relative peculiar velocities of the bound-zone neighbors around each host as

$v_{r}=\hat{\bf r}\cdot\left({\bf V}_{h}-{\bf V}_{b}\right)\,,$ (3)

where ${\bf V}_{h}$ and ${\bf V}_{b}$ denote the comoving peculiar velocities
of a host halo and its bound-zone neighbor halo, respectively, and $\hat{\bf
r}$ is the unit vector in the direction from the host halo center to the
bound-zone neighbor separated by a distance $r$. The rescaled bound-zone
velocities and separation distances, $\tilde{v}_{r}$ and $\tilde{r}$, are
defined as $\tilde{v}_{r}\equiv v_{r}/V_{c}$ and $\tilde{r}\equiv r/r_{v}$,
respectively. From here on, we will call $\tilde{v}_{r}$ the bound-zone
velocity profile, dropping the term, ”rescaled”, unless otherwise stated.

Breaking up the bound-zone distance range of $2<\tilde{r}<10$ into several
intervals of equal length, $\Delta\tilde{r}=1$, we take the average of
$\tilde{v}_{r}$ over those neighbors with $\tilde{r}$ falling in each interval
to determine the bond-zone velocity profile, $\tilde{v}_{r}(r)$, for each
host. Then, we take its ensemble average over all of the hosts to determine
the stacked bound-zone velocity profile, $\langle\tilde{v}_{r}\rangle$, the
result of which is shown in Figure 1. As can be seen, the average bound-zone
velocity profiles conspicuously differ between the ${\rm fR5+0.15}\,{\rm eV}$
and other three models. The former case yields significantly higher values of
$|\langle\tilde{v}_{r}\rangle|$ in the whole range of $\tilde{r}$ than the
latter case, which must be caused by the difference in the strength of the
fifth force between the two cases. This result implies that the free streaming
of more massive neutrinos present in the ${\rm fR5+0.15}\,{\rm eV}$ model do
not severely attenuate the effect of it stronger fifth force on the bound-zone
velocity profiles.

We fit $\langle\tilde{v}_{r}\rangle$ to Equation (1) by simultaneously
adjusting $A$ and $n$ with the help of the $\chi^{2}$-minimization method for
each model. Table 1 lists their best-fit values in the third and fourth
columns for the four cosmologies. Figure 2 compares the numerically obtained
$\langle\tilde{v}_{r}\rangle$ (black filled circles) to the analytic formula
with the best-fit values of $A$ and $n$ (red solid lines), revealing the
excellent agreements between the numerical and analytical results. This result
confirms the validity and usefulness of Equation (1) not only for the
$\Lambda{\rm CDM}$ model but also for the $\nu$CDM+$f(R)$ gravity models.
Figure 3 plots the $68\%$, $95\%$ and $99\%$ contours of $\chi^{2}(A,n)$,
quantifying how significant the differences in the best-values of $A$ and $n$
are between the ${\rm fR5+0.15}\,{\rm eV}$ and the other three models.

Recall that the ${\rm fR5+0.15}\,{\rm eV}$ model is very similar to the other
three ones in the conventional statistics (Baldi et al., 2014). Especially the
${\rm fR6}$ model has been shown to be almost indistinguishable from the ${\rm
fR5+0.15}\,{\rm eV}$, as the effect of the stronger fifth force,
$|f_{R0}|=10^{-5}$ is so severely attenuated by the free streaming massive
neutrinos with $\sum m_{\nu}=0.15\,{\rm eV}$ that the palpable effect amounts
only to that of $|f_{R0}|=10^{-6}$. However, the average bound-zone velocity
profile, unlike the aforementioned conventional statistics, is capable of
disentangling the effect of the fifth force from that of the massive
neutrinos, sensitively varying with the former, but not with the latter. Yet,
given that the other three models still remain mutually indistinguishable even
by $\langle v_{r}(r)\rangle$, we now investigate if any other statistics based
on $v_{r}(r)$ beyond its ensemble average can break the degeneracy among the
other three models. Basically, we consider the turn-around radii of the hosts
as such a statistics, and estimate them by applying the TRE reviewed in
Section 2 to the individual bound-zone velocity profiles.

For each host, we fit the individual bound-zone velocity profile to Equation
(1) and separately determine the best-fit values of $A$ and $n$. While
performing this fitting procedure, we exclude a small fraction of the hosts
whose bound-zone velocity profiles fail to be fitted by Equation (1) due to
the low number of their bound-zone neighbors. Table 1 lists the number of the
hosts, $N_{h}$, whose $v_{r}(r)$ matches Equation (1) for the four cosmologies
in the second column. Plugging the best-fit values of $A$ and $n$ into
Equation (2) and solving it, we estimate the turn-around radius, $r_{t}$, of
each of the included hosts. Counting the host halos whose turn-around radii
exceed $\alpha r_{v}$ as a function of a dimensionless variable $\alpha$, we
obtain the cumulative probability $P(r_{t}\geq\alpha r_{v})$. To assess the
errors, $\sigma_{P}$, in $P(r_{t}\geq\alpha r_{v})$, we create $10,000$
bootstrap resamples composed of equal number of the hosts and obtain the
cumulative probabilities from each resample. The one standard scatter among
the resamples from the average is taken as $\sigma_{P}$.

Figure 4 plots the cumulative probabilities, $P(r_{t}\geq\alpha r_{v})$, with
the bootstrap errors, $\sigma_{P}$, for the four models. As expected, the
${\rm fR5+0.15}\,{\rm eV}$ model yields the most conspicuously different
$P(r_{t}\geq\alpha r_{v})$ from the other three models. The host halos in the
${\rm fR5+0.15}\,{\rm eV}$ model seem to have much larger turn-around radii
than in the other three models, producing $22\%$ higher value of
$P(\alpha=4)$, corresponding to a $10\sigma_{\Delta P}$ signal, where the
uncertainties, $\sigma_{\Delta p}$, associated with the measurement of the
difference between the cumulative probabilities, are computed through the
propagation of the bootstrap errors.

Regarding the other three models, they yield almost identical values of
$P(r_{t}\geq\alpha r_{v})$ up to $\alpha=8$. However, the two models,
$\Lambda{\rm CDM}$ and ${\rm fR6+0.06}\,{\rm eV}$, which are in fact
indistinguishable by the conventional statistics (Baldi et al., 2014; Hagstotz
et al., 2019a), show different behaviors in the limit of $\alpha>8$. The
difference in the value of $P(\alpha=8.5)$ between the two models is found to
be as high as $3.2\sigma_{\Delta P}$ in spite of the strongest degeneracy
between them, yielding the same value of $\sigma_{8}$ (Baldi et al., 2014;
Hagstotz et al., 2019a). Table 1 lists $P(\alpha=4)$ and $P(\alpha=8.5)$ for
the four cosmologies in the fifth and sixth columns. Noting that the ${\rm
fR6+0.06}\,{\rm eV}$ model yields a smaller value of $P(\alpha=8.5)$ than the
other models and that the ${\rm fR6}$ model does not show any significant
difference in the whole range of $\alpha$ from the $\Lambda{\rm CDM}$ case, we
suspect the following: Even though the bound-zone velocity profiles are less
sensitive to the presence of massive neutrinos, the effect of $f(R)$ gravity
combined with massive neutrinos on the turn-around radii is different from
that of $f(R)$ gravity alone especially for the case that the $f(R)$ gravity
is not strong enough to prostrate the effect of free streaming massive
neutrinos.

We make the same analysis but of the cluster halos identified at two different
redshifts, $z=0.2$ and $0.4$ to see how $P(\alpha)$ changes with redshifts,
the results of which are shown in Figures 5 and 6, respectively. As can be
seen, $P(\alpha)$, diminishes more rapidly with $\alpha$ at higher redshifts.
While the statistical significance of the difference in $P(\alpha=4)$ between
the ${\rm fR5+0.15}\,{\rm eV}$ and the other three models are quite robust
against the redshift variation, that in $P(\alpha\geq 8)$ between the
$\Lambda{\rm CDM}$ and ${\rm fR6+0.06}\,{\rm eV}$ models drops to
$2.3\sigma_{\Delta P}$ at $z=0.2$ and to negligible level at $z=0.4$. The low
abundance of the cluster halos at higher redshifts contribute to the large
uncertainties in $P(\alpha)$, rendering it inconclusive whether or not the
strongest degeneracy between $\Lambda{\rm CDM}$ and ${\rm fR6+0.06}\,{\rm eV}$
can be broken by the turn-around radii of the cluster halos at $z\geq 0.2$.

##  4 Summary and Conclusion

We have numerically demonstrated that the turn-around radii of cluster halos
can in principle be useful to detect the effect of $f(R)$ gravity attenuated
by the presence of massive neutrinos. The samples of the cluster halos with
$M_{v}\geq 4.0\times 10^{14}\,h^{-1}\\!M_{\odot}$ at $z=0$ were obtained from
the DUSTGRAIN-pathfinder simulations (Giocoli et al., 2019) performed for four
different cosmologies: the Planck $\Lambda{\rm CDM}$ and three $\nu{\rm
CDM}$+$f(R)$ gravity models having different strength of fifth force and total
neutrino mass: ${\rm fR6}$, ${\rm fR6+0.06}\,{\rm eV}$, ${\rm fR5+0.15}\,{\rm
eV}$, which were known to be degenerate with the $\Lambda{\rm CDM}$ model and
with one another, yielding very similar conventional statistics (Baldi et al.,
2014; Hagstotz et al., 2019a).

For the determination of the turn-around radii, $r_{t}$, of the cluster halos
at which the value of the peculiar velocity field becomes equal to the
recession speed of the Hubble flow, we have employed the TRE (turn-around
radius estimator) developed by Lee et al. (2015), which in turn utilizes the
universal analytic formula for the peculiar velocity profile in the bound zone
around the cluster halos, put forth by Falco et al. (2014) for the
$\Lambda{\rm CDM}$ case. Our comparison of the analytic formula with the
average bound-zone velocity profile through the $\chi^{2}$ statistics has
confirmed its validity not only for the Planck $\Lambda{\rm CDM}$ but also for
the three $\nu{\rm CDM}$+$f(R)$ gravity models. It has also revealed that the
amplitudes and slopes of the bound-zone velocity profiles, quantified by its
two adjustable parameters, significantly differ between the ${\rm
fR5+0.15}\,{\rm eV}$ and the other three models (Figures 1-3). This result
implies that the bound-zone velocities of the cluster halos must be much more
susceptible to the presence of the strong fifth force than to that of massive
neutrinos and thus that some statistics based on them may be useful to
disentangle the former from the latter.

With the turn-around radii of the cluster halos estimated by the application
of the TRE to the bound-zone velocity profiles of individual cluster halos, we
have determined the cumulative probability distributions of the turn-around to
virial radius ratios, $P(r_{t}\geq\alpha r_{v})$. With the help of the
bootstrap statistics, we have shown that the ${\rm fR5+0.15}\,{\rm eV}$ model
can be plainly differentiated by $P(r_{t}/r_{v}\geq 4)$ from the other three
cases, with statistical significance as high as $10\sigma_{\Delta P}$ (Figure
4). We have also detected a $3.2\sigma_{\Delta P}$ difference in
$P(r_{t}/r_{v}\geq 8.5)$ between the $\Lambda{\rm CDM}$ and ${\rm
fR6+0.06}\,{\rm eV}$ models, in spite of the strongest degeneracy between the
two cases. Yet, given the low value of $P(r_{t}/r_{v}\geq 8.5)\sim{\cal
O}(10^{-4})$, it is not conclusive whether the $3.2\sigma_{\Delta P}$
difference in $P(r_{t}/r_{v}\geq 8.5)$ between the $\Lambda{\rm CDM}$ and
${\rm fR6+0.06}\,{\rm eV}$ models is a real signal or just a spurious one
produced by the shot noise. A larger sample of the massive cluster halos with
$r_{t}\geq 8r_{v}$) will be required to confirm the statistical significance
of the difference in $P(\alpha>8)$ between the two models. It has been also
shown that the significance of the difference between the ${\rm
fR5+0.15}\,{\rm eV}$ and the other three models is robust against the
variation of the redshifts from $z=0.0$ to $0.4$ (Figures 5-6).

The advantage of using the turn-around radii estimated by the TRE to
distinguish among the $\Lambda{\rm CDM}$ and $\nu{\rm CDM}$+$f(R)$ gravity
models comes from the universality of the bound-zone velocity profile on which
the TRE is based. As revealed by the previous works (Falco et al., 2014; Lee,
2016), the shape of the bound-zone velocity profile is insensitive to the key
cosmological parameters of the $\Lambda{\rm CDM}$ model. Thus, the variation
of $\sigma_{8}$ and $\Omega_{m}$ in the $\Lambda{\rm CDM}$ model cannot
produce the same effect on the turn-around radii as the $\nu{\rm CDM}$+$f(R)$
gravity. Furthermore, it does not require to track down the redshift evolution
of $r_{t}$ unlike the previously suggested statistics as a possible
discriminator of the $\nu{\rm CDM}$+$f(R)$ gravity model such as the size
evolutions of galaxy voids, nonlinear growth rates and redshift distortions,
evolution of the drifting average coefficient of the field cluster mass
function, and high order weak lensing statistics (Peel et al., 2018; Hagstotz
et al., 2019b; Giocoli et al., 2019; Wright et al., 2019; Ryu et al., 2020;
Contarini et al., 2021).

It is, however, worth discussing the practical difficulties of our statistics.
To detect a signal of the difference in $P(\alpha\geq 8)$ strong enough to
distinguish among the degenerate models, what is required is to measure the
turn-around radii of as many galaxy clusters as possible in the local
universe. However, as shown in Lee (2018) and Hansen et al. (2020), the TRE is
applicable only to those isolated galaxy clusters whose bound zone neighbor
galaxies exhibit very high degree of anisotropy in their spatial
distributions. Due to this limitation of the TRE, it would be difficult to
obtain a large sample of the galaxy clusters with their turn-around radii
measurable without information on the peculiar velocity field.
Notwithstanding, we expect that the large peculiar velocity dataset available
from the future galaxy surveys like The Large Synoptic Survey Telescope (LSST)
survey (Tyson, 2002) should allow us to directly measure the turn-around radii
of almost all of the galaxy clusters at low-redshifts $z\leq 0.2$, making our
statistics based on the turn-around radii to be practically useful as a
powerful discriminator of the $\nu{\rm CDM}$+$f(R)$ gravity models.

JL acknowledges the support by Basic Science Research Program through the
National Research Foundation (NRF) of Korea funded by the Ministry of
Education (No.2019R1A2C1083855). JL thanks S.Ryu for having been helpful in
obtaining the rockstar halo catalogs from the Dustgrain-Pathfinder
simulations. MB acknowledges support by the project ”Combining Cosmic
Microwave Background and Large Scale Structure data: an Integrated Approach
for Addressing Fundamental Questions in Cosmology”, funded by the MIUR
Progetti di Ricerca di Rilevante Interesse Nazionale (PRIN) Bando 2017 - grant
2017YJYZAH. MB also acknowledges the use of computational resources from the
parallel computing cluster of the Open Physics Hub at the Physics and
Astronomy Department in Bologna (https://site.unibo.it/openphysicshub/en) .

## References

  * Amendola (2000) Amendola, L. 2000, Phys. Rev. D, 62, 043511
  * Amendola (2004) Amendola, L. 2004, Phys. Rev. D, 69, 103524
  * Amendola et al. (2013) Amendola, L., Appleby, S., Bacon, D., et al. 2013, Living Reviews in Relativity, 16, 6
  * Baldi et al. (2010) Baldi, M., Pettorino, V., Robbers, G., et al. 2010, MNRAS, 403, 1684
  * Baldi (2012a) Baldi, M. 2012a, MNRAS, 420, 430
  * Baldi (2012b) Baldi, M. 2012b, MNRAS, 422, 1028
  * Baldi et al. (2014) Baldi, M., Villaescusa-Navarro, F., Viel, M., et al. 2014, MNRAS, 440, 75
  * Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013, ApJ, 762, 109
  * Bhattacharya & Kousvos (2017) Bhattacharya, S. & Kousvos, S. R. 2017, Phys. Rev. D, 96, 104006
  * Bhattacharya & Tomaras (2017) Bhattacharya, S. & Tomaras, T. N. 2017, European Physical Journal C, 77, 526
  * Bhattacharya & Tomaras (2021) Bhattacharya, S. & Tomaras, T. N. 2021, Annals of Physics, 427, 168427
  * Blanchard et al. (2022) Blanchard, A., Héloret, J.-Y., Ilić, S., et al. 2022, arXiv:2205.05017
  * Bond et al. (1996) Bond, J. R., Kofman, L., & Pogosyan, D. 1996, Nature, 380, 603
  * Bond & Myers (1996) Bond, J. R. & Myers, S. T. 1996, ApJS, 103, 1
  * Boubekeur et al. (2014) Boubekeur, L., Giusarma, E., Mena, O., et al. 2014, Phys. Rev. D, 90, 103512
  * Brax et al. (2004) Brax, P., van de Bruck, C., & Davis, A.-C. 2004, Reports on Progress in Physics, 67, 2183
  * Buchdahl (1970) Buchdahl, H. A. 1970, MNRAS, 150, 1
  * Capozziello et al. (2019) Capozziello, S., Dialektopoulos, K. F., & Luongo, O. 2019, International Journal of Modern Physics D, 28, 1950058
  * Carroll (2001) Carroll, S. M. 2001, Living Reviews in Relativity, 4, 1
  * Chuang et al. (2021) Chuang, Y.-T., Okumura, T., & Shirasaki, M. 2021, arXiv:2111.01417
  * Clifton et al. (2012) Clifton, T., Ferreira, P. G., Padilla, A., et al. 2012, Phys. Rep., 513, 1
  * Contarini et al. (2021) Contarini, S., Marulli, F., Moscardini, L., et al. 2021, MNRAS, 504, 5021
  * De Felice & Tsujikawa (2010) De Felice, A. & Tsujikawa, S. 2010, Living Reviews in Relativity, 13, 3. doi:10.12942/lrr-2010-3
  * Dvali et al. (2000) Dvali, G., Gabadadze, G., & Porrati, M. 2000, Physics Letters B, 485, 208
  * Falco et al. (2014) Falco, M., Hansen, S. H., Wojtak, R., et al. 2014, MNRAS, 442, 1887
  * Faraoni (2021) Faraoni, V. 2021, Journal of Physics Conference Series, 2156, 012017
  * Giocoli et al. (2019) Giocoli, C., M., Baldi, M., & Moscardini, L. 2018, MNRAS, 481, 2813
  * Giusti & Faraoni (2019) Giusti, A. & Faraoni, V. 2019, Physics of the Dark Universe, 26, 100353
  * Giusti & Faraoni (2021) Giusti, A. & Faraoni, V. 2021, Phys. Rev. D, 103, 044049
  * Gunn & Gott (1972) Gunn, J. E. & Gott, J. R. 1972, ApJ, 176, 1
  * Hagstotz et al. (2019a) Hagstotz, S., Costanzi, M., Baldi, M., et al. 2019, MNRAS, 486, 3927
  * Hagstotz et al. (2019b) Hagstotz, S., Gronke, M., Mota, D. F., et al. 2019, A&A, 629, A46
  * Hansen et al. (2020) Hansen, S. H., Hassani, F., Lombriser, L., et al. 2020, JCAP, 2020, 048. doi:10.1088/1475-7516/2020/01/048
  * Hojjati et al. (2011) Hojjati, A., Pogosian, L. & Zhao, G. 2011, JCAP, 1108, 005
  * Hu & Sawicki (2007) Hu, W., & Sawicki, I. 2007, Phys. Rev. D, 76, 064004
  * Jibrail et al. (2020) Jibrail, A. W., Elahi, P. J., & Lewis, G. F. 2020, MNRAS, 492, 2369
  * Khoury & Weltman (2004) Khoury, J. & Weltman, A. 2004, Phys. Rev. D, 69, 044026
  * Komatsu et al. (2011) Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18
  * Lee et al. (2015) Lee, J., Kim, S., & Rey, S.-C. 2015, ApJ, 815, 43
  * Lee (2016) Lee, J. 2016, ApJ, 832, 123
  * Lee & Yepes (2016) Lee, J. & Yepes, G. 2016, ApJ, 832, 185
  * Lee & Li (2017) Lee, J. & Li, B. 2017, ApJ, 842, 2
  * Lee (2018) Lee, J. 2018, ApJ, 856, 5
  * Lee et al. (2022) Lee, J., Ryu, S., & Baldi, M. 2022, arXiv:2206.03406
  * Lesgourgues & Pastor (2014) Lesgourgues, J. & Pastor, S. 2014, New Journal of Physics, 16, 065002. doi:10.1088/1367-2630/16/6/065002
  * Li & Barrow (2007) Li, B., & Barrow, J. D. 2007, Phys. Rev. D, 75, 084010
  * Lopes et al. (2018) Lopes, R. C. C., Voivodic, R., Abramo, L. R., et al. 2018, JCAP, 2018, 010
  * Lopes et al. (2019) Lopes, R. C. C., Voivodic, R., Abramo, L. R., et al. 2019, JCAP, 2019, 026
  * Macciò et al. (2004) Macciò, A. V., Quercellini, C., Mainini, R., et al. 2004, Phys. Rev. D, 69, 123516
  * Nojiri et al. (2018) Nojiri, S., Odintsov, S. D., & Faraoni, V. 2018, Phys. Rev. D, 98, 024005
  * Onuora & Thomas (2000) Onuora, L. I. & Thomas, P. A. 2000, MNRAS, 319, 614
  * Oyaizu (2008) Oyaizu, H. 2008, Phys. Rev. D, 78, 123523
  * Pavlidou & Tomaras (2014) Pavlidou, V. & Tomaras, T. N. 2014, JCAP, 2014, 020
  * Pavlidou et al. (2014) Pavlidou, V., Tetradis, N., & Tomaras, T. N. 2014, JCAP, 2014, 017
  * Peel et al. (2018) Peel, A., Pettorino, V., Giocoli, C., et al. 2018, A&A, 619, A38
  * Peel et al. (2019) Peel, A., Lalande, F., Stark, J.-L., et al. 2018, Phys. Rev. D, 100, 023508
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A24
  * Puchwein et al. (2013) Puchwein, E., Baldi, M., & Springel, V. 2013, MNRAS, 436, 348
  * Ryu et al. (2020) Ryu, S., Lee, J., & Baldi, M. 2020, ApJ, 904, 93
  * Sotiriou & Faraoni (2010) Sotiriou, T. P. & Faraoni, V. 2010, Reviews of Modern Physics, 82, 451
  * Springel (2005) Springel, V. 2005, MNRAS, 364, 1105
  * Springel et al. (2001) Springel, V., White, S. D. M., Tormen, G., et al. 2001, MNRAS, 328, 726
  * Starobinsky (1980) Starobinsky, A. A. 1980, Phys. Rev. B, 91, 99
  * Tanoglidis et al. (2015) Tanoglidis, D., Pavlidou, V., & Tomaras, T. N. 2015, JCAP, 2015, 060
  * Tyson (2002) Tyson, J. A. 2002, Proc. SPIE, 4836, 10. doi:10.1117/12.456772
  * Viel et al. (2010) Viel, M., Haehnelt, M. G., & Springel, V. 2010, JCAP, 2010, 015
  * Weinberg (1989) Weinberg, S. 1989, Reviews of Modern Physics, 61, 1
  * Wetterich (1995) Wetterich, C. 1995, A&A, 301, 321
  * Wright et al. (2019) Wright, B. S., Koyama, K., Winther, H. A., et al. 2019, JCAP, 06, 040
  * Zucca et al. (2019) Zucca, A., Pogosian, L., Silvestri, A., et al. 2019, JCAP, 1905, 001

Figure 1: Bound-zone velocity profiles around the cluster halos with
$M_{v}\geq 4.05\times 10^{13}\,h^{-1}\,M_{\odot}$ at $z=0$ for four different
cosmological models.

Figure 2: Comparison of the numerically obtained bound-zone velocity profiles
(black filled circles) with the best-fit analytic formula (red solid lines) at
$z=0$.

Figure 3: $68\%$, $95\%$ and $99\%$ contours from the $\chi^{2}$-statistics
for two parameters, $A$ and $n$, that characterize the analytic formula for
the bound-zone velocity profiles, Equation (1), at $z=0$.

Figure 4: Cumulative probability function of the ratios of the turn-around
radii to the virial counterparts at $z=0$.

Figure 5: Same as Figure 4 but at $z=0.2$.

Figure 6: Same as Figure 4 but at $z=0.4$.

Table 1: Best Parameters and Probabilities of $r_{t}\geq\alpha r_{v}$ at $z=0$ Cosmology | $N_{h}$ | $A$ | $n$ | $P(\alpha=4)$ | $P(\alpha=8.5)$
---|---|---|---|---|---
|  |  |  | ($10^{-2})$ | ($10^{-3}$)
$\Lambda{\rm CDM}$ | $18188$ | $1.061\pm 0.045$ | $0.580\pm 0.023$ | $43.50\pm 0.37$ | $0.87\pm 0.22$
fR6 | $19317$ | $1.086\pm 0.041$ | $0.590\pm 0.021$ | $43.69\pm 0.36$ | $0.50\pm 0.16$
fR6+$0.06\,{\rm eV}$ | $18811$ | $1.096\pm 0.043$ | $0.595\pm 0.021$ | $43.61\pm 0.37$ | $0.11\pm 0.08$
fR5+$0.15\,{\rm eV}$ | $19291$ | $1.245\pm 0.053$ | $0.625\pm 0.035$ | $53.74\pm 0.36$ | $1.37\pm 0.28$


# Gauge dependence of spontaneous radiation spectrum of relativistic atomic
beam under non-uniform electrostatic field

Xue-Nan Chen  For correspondence<EMAIL_ADDRESS>Hubei Key Laboratory of
Gravitational and Quantum Physics, School of Physics, Huazhong University of
Science and Technology, Wuhan 430074, China  Yu-Hang Luo  Hubei Key Laboratory
of Gravitational and Quantum Physics, School of Physics, Huazhong University
of Science and Technology, Wuhan 430074, China  Xiang-Song Chen  For
correspondence<EMAIL_ADDRESS>Hubei Key Laboratory of Gravitational and
Quantum Physics, School of Physics, Huazhong University of Science and
Technology, Wuhan 430074, China

###### Abstract

Gauge theory requires physical observables to be gauge-independent. However,
ever since Lamb noticed the problem of gauge selection in calculating atomic
spontaneous radiation spectrum, the problem of gauge dependence was
encountered in many fields of physics research. Therefore, it is important to
test the self-consistency of gauge symmetry for various physical systems. In
this paper, we calculate the transient spontaneous radiation spectrum of a
relativistic hydrogen atom in the non-uniform electrostatic field under the
atomic self-reference frame. The physical system studied in this paper is a
frame-transformed version of our recent work [Euro. J. Phys. D 76, 84(2022)]
where the radiating object is static while the charge is moving
relativistically. The obtained peak frequency can differ by about $413$
$\mathrm{KHz}$ or larger for the commonly used Coulomb, Lorentz, and
multipolar gauges. This observation can be significant not only for studying
how the gauge field interacts with the quantum system in theory, but also for
practical experimental applications, such as the timing accuracy of atomic
clocks in the external electromagnetic field.

##  I Introduction

The gauge theory constructs the interaction between matter and gauge field,
and the self-consistency of gauge theory requires physical observables to be
gauge-independent. For classical physics and scattering problems of elementary
particles, gauge symmetry can be guaranteed. However, when studying the
processes involving quantum bound states, which was first noticed by Lamb in
1952 Lamb52 , the issue of gauge symmetry seems to be tricky. Lamb found that
the experimental results matched well with the results of length gauge
$-q\vec{E}\cdot\vec{r}$ rather than velocity gauge $-(q/m)\vec{A}\cdot\vec{p}$
when studying atomic spectral lines. The problem of gauge dependence also
appears in the research fields of the multi-photon process Bassani77 ; Kobe78
; Olariu79 , the multi-electron system Starace71 ; Anderson74 ; Lin78 ; Kobe79
, strong field physics Bauer06 ; Faisal07a ; Faisal07b ; Vanne09 , cavity QED
Dicke54 ; Rzazewski75 ; Keeling07 ; Bernardis18 , and relativistic Coulomb
excitation Bayman05 ; Dasso06 . Different gauge conditions give different
predictions, which contradicts the fact that the experimental results are
unique.So people tended to think that there is a kind of “preferred gauge
condition” in these physical research fields Forney77 . At the same time, many
scholars constructed the so-called “gauge-independent” Hamiltonian operator to
solve this conflict Yang76 ; Kobe82 ; Lamb87 ; Funai19 .

However, the “gauge-invariant” Hamiltonian operators can currently only
applies to the study of quantum systems in weak radiation fields under the
dipole approximation. Once the electromagnetic scalar potential is time-
varying,the “gauge-invariant” Hamiltonian operators cannot be applicable. In
this case, the Hamiltonian of the system is

${}H(t)=\frac{1}{2m_{q}}(\vec{p}-q\vec{A}(\vec{x},t))^{2}+q\phi(\vec{x},t)+qV(\vec{x}),$
(1)

Certainly, the “total” $H(t)$ is gauge dependent. When using $H(t)$ to define
eigenstates and calculate the transition rates, we found they both are gauge
dependent. In our recent work Chen22 , we constructed a quantum system
interacting with a relativistic strong Coulomb field to show the
inapplicability of the “gauge-invariant” Hamiltonian. Although the Schrödinger
equation under the external field approximation is gauge covariant, there is a
problem that the definition of quantum states depends on gauge. Using the
external field approximation method, we calculated the transient spontaneous
emission spectrum of a one-dimensional quantum harmonic oscillator in the
field produced by relativistic moving charge clusters under multipolar gauge,
Coulomb gauge, and Lorentz gauge, and found that different gauge conditions
give different predictions. We need to use various physical systems to obtain
the effective external field while an appropriate and full quantum-field
method is still lacking Gross82  and compare the transient radiation spectrum
with the result under different gauge.

To test the self-consistency of gauge symmetry, we require that there are
relativistic particles in the physical system, and the time-varying scalar
potential affects the system adiabatically. Under this experimental demand,
the results calculated under different gauge conditions can be well
distinguished within the experimental error range. A realistic system is a
relativistic atomic beam. In the 1970s, LAMPF could prepare relativistic
$\mathrm{H}^{-}$ ions, and relativistic atomic beams could be obtained through
photodetachment Bryant71 . Moreover, with the implementation of ion storage
ring projects such as ASTRID ASTRID , CRYRING CRYRING , ESR ESR , TARN II TARN
II , TSR TSR , IUCF IUCF , CSR CSR , etc., we will get relativistic ion beams
more conveniently. Using these experimental facilities, we only need to
consider letting relativistic ions or atoms pass through a non-uniform
electromagnetic field area to realize our experimental proposal.

In this paper, we adopt the experimental scheme of relativistic atoms to test
the self-consistency of gauge symmetry. The paper is organized as follows: In
Sect. II, we introduce our physical system: a non-uniform electrostatic field
acting adiabatically on a relativistic atomic beam. In the atomic self-
reference, parameters can be adjusted to make $\phi_{\text{ext}}(\vec{x},t)$
of the charged ring to be adiabatic for the moving atom. In Sect. III, we
solve the instantaneous eigen-states of our system by using external field
approximation. The solutions indeed differ significantly for the commonly used
Coulomb, Lorentz, and multipolar gauges. Then in Sect. IV, we compute the
transient spontaneous radiation spectrum explicitly and find again a
significant gauge dependence. In the last section, we summarize our results
and discuss their physical implications.

##  II A non-uniform electrostatic field acting adiabatically on a
relativistic atomic beam

Figure 1: Schematic design of the physical system. The blue solid circle
represents a relativistic hydrogen atom with the velocity $v$. The blue coil
represents a charged ring with the charge $+Ne$ and the radius $R$.

Figure 1 shows a schematic design of our physical system. The non-uniform
electrostatic field is produced by a charged ring with charge $+Ne$ and radius
$R$. A relativistic hydrogen atom with velocity $v$ passes through the charged
ring. In order to meet the needs of observing transient spontaneous radiation
spectrum, we need to meet the adiabatic evolution of excited states.
Degenerate quantum states are difficult to meet the requirements of adiabatic
conditions, so we need to open the degeneracy of quantum states. To ensure
that the energy levels are not degenerate when relativistic atoms pass through
the charged ring, we provide a uniform electric and magnetic field along the
direction of relativistic atomic motion. In principle, only in the case of
non-relativsitic motion, we can use centroid separation method to study the
quantum system under the laboratory reference frame Kaiser17 . In the case of
relativistic motion, we can only study the problem in atomic self-reference
frame of the quantum system Bryant83 . So in this paper, we work in the atomic
self-reference frame. In the actual data processing, we can transform the
experimental results into the atomic self-reference frame. In the atomic self-
reference frame, the coordinate system is established with the center of the
hydrogen atom as the coordinate origin and the moving direction as the
positive direction of the $z$-axis. Because the applied uniform electric field
and magnetic fields are along the $z$-axis, the magnitude of the electric and
magnetic fields will not change when transforming from the laboratory
reference system to the atomic self-reference frame. After the coordinate
system is established, the specific expressions of the electric and magnetic
fields are $\vec{E}=(0,0,E)$ and $\vec{B}=(0,0,B)$ respectively.

If excited, the relativistic atom can emit a photon by spontaneous radiation.
It can be expected that by adjusting the parameter $N$ and $R$ and selecting
the appropriate excited state, in the atomic self-reference frame, the scalar
potential $\phi_{\text{ext}}(\vec{x},t)$ of the charged ring can be adiabatic
for the relativistic atom for duration which is sufficiently long for the
excited relativistic atom to emit a transient photon. The velocity of the
relativistic atom $v$ should be chosen large enough to produce a considerable
difference for $\phi_{\text{ext}}(\vec{x},t)$ among various gauges, and at the
same time small enough to allow for a rough external-field approximation.
Under the external-field approximation method, the Hamiltonian describing the
physical system of our schematic design is

${}H_{G}(t)=\frac{1}{2m_{e}}\left[\vec{p}+e\frac{1}{2}(\vec{B}\times\vec{x})+e\vec{A}_{\text{ext}}(\vec{x},t)\right]^{2}+V(\vec{x})+eEz-e\phi_{\text{ext
}}(\vec{x},t),$ (2)

where, $V(\vec{x})=-e^{2}/r$ is the binding potential energy of hydrogen atom.
$(\phi_{\text{ext}}(\vec{x},t),\vec{A}_{\text{ext}}(\vec{x},t))$ is the gauge
potential generated by a charged ring in the atomic self-reference frame.

##  III The instantaneous eigen-states of the relativistic atom in atomic
self-reference frame

To compute the radiation spectrum of the relativistic atom , one must first
define its eigen-states. In our case, this is done by solving

${}H_{G}(t)|\alpha_{G}(t)\rangle=W_{\alpha G}|\alpha_{G}(t)\rangle,$ (3)

According to Chen22 , We may always apply to Eq. (3) the following unitary
transformation:

${}|\alpha_{G}(t)\rangle=U(t)|\alpha^{\prime}_{G}(t)\rangle,\quad
H_{G}(t)=U(t)H^{\prime}_{G}(t)U^{-1}(t),$ (4)

then, Eq. (3) becomes

${}H^{\prime}_{G}(t)|\alpha^{\prime}_{G}(t)\rangle=W_{\alpha
G}|\alpha^{\prime}_{G}(t)\rangle,$ (5)

where the new operator is expressed as

${}H^{\prime}_{G}(t)=\frac{1}{2m_{e}}\left[\vec{p}+\frac{e}{2}\left(\left(\vec{B}+\vec{B}_{\mathrm{ext}}(0,t)\right)\times\vec{x}\right)\right]^{2}+V(\vec{x})+eEz-e\phi_{\mathrm{ext}}(\vec{x},t),$
(6)

when choosing appropriate unitary transformation $U(t)$. Notice that unitary
transformation and operator $H^{\prime}_{G}(t)$ are merely a mathematical
technique, after fixing a gauge, for the convenience of solving the eigen-
equation. Hence, when solving Eq. (3) of the relativistic atom, the vector
potential $\vec{A}_{\text{ext}}(\vec{x},t)$ is trivial because after unitary
transformation, vector potential $\vec{A}_{\text{ext}}$ is transformed into
the gauge-invariant form $(1/2)\vec{B}_{\text{ext}}\times\vec{x}$. And only
the scalar potential $\phi_{\text{ext}}(\vec{x},t)$ is essential.

In our schematic physical system, the relativistic atomic beam moves along the
symmetry axis of the charged ring, Hence, $\vec{B}_{\text{ext }}(0,t)=0$. When
the uniform magnetic field $\vec{B}$ is weak, we can ignore the effect of the
$\vec{B}^{2}$ term. For the effect of the external electric field, we only
consider the effect of the electric dipole. Hence, we obtain

${}H_{G}^{\prime}(t)\approx\frac{\vec{p}^{2}}{2m_{e}}+V(\vec{x})+\frac{e}{2m_{e}}B_{z}l_{z}+eEz-e(\phi_{\text{ext
}}(0,t)+\partial_{z}\phi_{\text{ext}}(\vec{x},t)|_{\vec{x}=0}z),$ (7)

where $l_{z}$ is the angular momentum operator in $z$ direction of the
relativistic atom. In this paper, we mainly study the definition of quantum
state under Coulomb, Lorentz, and multipolar gauge. The specific forms and
Taylor expansion forms of scalar potential $\phi_{\text{ext}}(\vec{x},t)$ of
three gauge conditions are shown in Appendix A.

Different from Chen22  in solving the quantum state of a one-dimensional
harmonic oscillator, the solution of the relativistic atom’s state affected by
non-uniform electrostatic field is difficult to solve by the analytic method.
Hence, we use the method of numerical calculation to solve it. We use the
eigenstates $|nlm\rangle$ of free Hamiltonian
$H_{0}=\vec{p}^{2}/(2m_{e})+V(\vec{x})$ to expand quantum states

${}\left|\alpha_{G}^{\prime}(t)\right\rangle=\sum_{nlm}c_{\alpha_{G}^{\prime},nlm}(t)|nlm\rangle.$
(8)

Here, operator $H_{G}^{\prime}(t)$ is expressed as a matrix by using the
eigenstate of free Hamiltonian, in which one matrix element is expressed as

$\displaystyle\left\langle
n^{\prime}l^{\prime}m^{\prime}\left|H_{G}\left(t^{\prime}\right)\right|nlm\right\rangle=$
$\displaystyle\left(W_{n}+m\hbar\omega_{L}-e\phi_{\mathrm{ext}}(0,t)\right)\delta_{n^{\prime}n}\delta_{l^{\prime}l}\delta_{m^{\prime}m}$
(9)
$\displaystyle+e\left(E-\left.\partial_{z}\phi_{\mathrm{ext}}(\vec{x},t)\right|_{\vec{x}=0}\right)\left\langle
n^{\prime}l^{\prime}m^{\prime}|z|nlm\right\rangle,$

where, $\omega_{L}=eB/(2m_{e})$.

By diagonalizing the matrix form of the operator $H_{G}^{\prime}(t)$, we can
obtain the eigenvalues $W_{\alpha G}(t)$ and corresponding eigenvectors
$|\alpha_{G}(t)\rangle$ under the external field approximation. Here, the
expansion coefficient $c_{\alpha_{G}^{\prime},nlm}(t)$ of the eigenstate is
given by each component of the eigenvector obtained by diagonalization method.
After some rough estimation, we find that the following values suffice our
study: $N=10^{12}$, $\beta=v/c=0.1$, and $R=2.8$ $\mathrm{cm}$. If the
diameter of the charged ring is considered to be $1$ $\mathrm{mm}$, the
voltage on the charged ring is in the order of $10^{5}$ $\mathrm{V}$. The
intensity of uniform electric and magnetic fields applied in our physical
system are $E=1.03\times 10^{6}$ $\mathrm{V}/\mathrm{m}$ and $B=2.12$
$\mathrm{T}$ respectively. The relativistic atom moves from $z_{i}=-0.33$
$\mathrm{m}$ to $z_{f}=0.33$ $\mathrm{m}$. In the following calculation and
data expression, we will use the atomic unit system. Hence,
$W_{n}=-1/(2n^{2})$ and the operator $H_{G}^{\prime}(t)$ can be expressed as

${}H_{G}^{\prime}(t)=H_{0}+E_{G}(t)H_{z}+4.5\times 10^{-6}m,$ (10)

where,
$E_{G}(t)=E-\left.\partial_{z}\phi_{\mathrm{ext}}(\vec{x},t)\right|_{\vec{x}=0}$.
$n$, $l$ and $m$ are the principal quantum number, the orbital quantum number,
and the magnetic quantum number respectively. And
$H_{z,n^{\prime}l^{\prime}m^{\prime},nlm}=\left\langle
n^{\prime}l^{\prime}m^{\prime}|z|nlm\right\rangle$.

We first show the results of $E_{G}(t)$ under the multipolar gauge condition
with Fig. 2. Figure 3 shows the results of $E_{G}(t)$ under the Coulomb and
Lorentz gauge conditions as differences from the results under the multipolar
gauge conditions.

Figure 2: The results of $E_{G}(t)$ under the multipolar gauge condition.
Figure 3: The results of $E_{G}(t)$ under the Coulomb and Lorentz gauge
conditions as differences from the results under the multipolar gauge
conditions.

The numerical diagonalization method can only deal with finite-dimensional
matrices, so we need to truncate the principal quantum number when calculating
quantum states and eigenvalues. In this paper, we truncate the principal
quantum number at $n=15$. Because only the electric dipole effect of the non-
uniform electrostatic field is considered, the matrix element $E_{G}(t)H_{z}$
contributed by the non-uniform electrostatic field is non-zero only when
$l^{\prime}-l=\pm 1$ and $m^{\prime}=m$. Hence, we can use the magnetic
quantum number $m$ to divide the matrix $H_{z}(t)$ and diagonalize it in
block. We can obtain the definition of the eigenvector after diagonalization,
that is, the quantum state of relativistic atoms in a non-uniform
electrostatic field. These quantum states can rank the eigenvalues from small
to large under the same magnetic quantum number $m$, and will be given the
serial number $N$. In this way, the quantum state can be defined by two
parameters, serial number $N$ and magnetic quantum number $m$, namely
$|\alpha_{G}(t)\rangle=|(N,m)_{G}(t)\rangle$.

In order to maximize the difference of transient spontaneous emission spectrum
under different gauge conditions under the parameters of our physical system,
and to satisfy the adiabatic conditions in the following, we choose $(53,0)$
as the excited state to be studied. Figure 4 shows the result of the
eigenvalue of $(53,0)$ under the multipolar gauge condition. Figure 5 shows
the results of the eigenvalue of $(53,0)$ under the Coulomb and Lorentz gauge
conditions as differences from the results under the multipolar gauge
condition.

Figure 4: The result of the eigenvalue of $(53,0)$ under the multipolar gauge
condition. Figure 5: The results of the eigenvalue of $(53,0)$ under the
Coulomb and Lorentz gauge conditions as differences from the results under the
multipolar gauge condition.

Before solving the eigen-equations, we must check that with our chosen
parameters, $\phi_{\text{ext}}(\vec{x},t)$ is indeed an adiabatic potential
for the relativistic atom, in all three gauges we use. The criterion for the
adiabatic approximation is the parameter

${}r_{\beta\alpha}=\left|\frac{\hbar\langle\beta_{G}(t)\mid\dot{\alpha}_{G}(t)\rangle}{W_{\beta
G}(t)-W_{\alpha G}(t)}\right|.$ (11)

where $\alpha,\beta$ label two different states. We evaluate the adiabatic
parameters of the quantum state $|\alpha_{G}(t)\rangle=|(53,0)_{G}(t)\rangle$.
The maximum adiabatic parameters are of the order of $10^{-4}$ under multipole
gauge and Coulomb gauge, and of $10^{-3}$ under Lorentz gauge. Thus, if
prepared in the state $|(53,0)\rangle$ in the initial moment, the relativistic
atom will largely stay in the state $|(53,0)_{G}(t)\rangle$ during the whole
process.

With the adiabatic conditions verified, the gauge dependence of the
instantaneous eigenvalues, as we just aforementioned for our designed
relativistic atom, may in principle be tested experimentally. A possible
experimental observable is the transient spectrum of spontaneous radiation, as
we will compute it in the next section.

##  IV The transient spontaneous radiation spectrum and its gauge dependence
in atomic self-reference frame

We now add into our system the coupling to the background vacuum
electromagnetic field. Because in Chen22 , we find that “Coulomb-type” gauge
dependence is more obvious than “Lamb-type” gauge dependence, so in this
paper, the gauge condition of background electromagnetic field is taken as
multipolar gauge. The total Hamiltonian of the system coupled with the vacuum
background field is

${}H_{\mathrm{tot}}=H_{G}(t)-q\vec{E}_{\text{vac}}\cdot\vec{r}+H_{B}$ (12)

where, $H_{G}(t)$ is the Hamiltonian of the relativistic atoms under the non-
uniform electrostatic field. $E_{\text{vac}}$ is vacuum background field
strength

${}\vec{E}_{\text{vac}}=i\sum_{\vec{k}}\sum_{\lambda}\sqrt{\frac{\hbar\omega_{\vec{k}}}{2\varepsilon_{0}V}}\left(a_{\vec{k},\lambda}\vec{e}_{\vec{k}}e^{i\vec{k}\cdot\vec{r}}-a_{\vec{k},\lambda}^{\dagger}\vec{e}_{\vec{k}}^{*\lambda}e^{-i\vec{k}\cdot\vec{r}}\right).$
(13)

Here,

${}\sum_{\vec{k}}\rightarrow\frac{V}{(2\pi)^{3}}\int d^{3}k,$ (14)

and physically detectable photons should satisfy
$\vec{k}\cdot\vec{e}_{\vec{k}}^{\lambda}=0$. The Hamiltonian of the vacuum
background field is

${}H_{B}=\sum_{\vec{k}}\sum_{\lambda}\hbar\omega_{\vec{k}}\left(a_{\vec{k},\lambda}^{\dagger}a_{\vec{k},\lambda}+\frac{1}{2}\right).$
(15)

The second term is often omitted in quantum optical calculations. Since the
scale of atoms is much smaller than the wavelength scale of the vacuum
background field, only dipole interactions

${}\vec{E}_{\mathrm{vac}}\approx
i\sum_{\vec{k}}\sum_{\lambda}\sqrt{\frac{\hbar\omega_{\vec{k}}}{2\varepsilon_{0}V}}\left(a_{\vec{k},\lambda}\vec{e}_{\vec{k}}-a_{\vec{k},\lambda}^{\dagger}\vec{e}_{\vec{k}}^{*\lambda}\right)$
(16)

need to be considered. By using the eigenstates of relativistic atoms in the
non-uniform electrostatic field solved in the previous section, we can rewrite
the total Hamiltonian as

${}H_{\text{tot}}=\sum_{\alpha}E_{G,\alpha}(t)\sigma_{\alpha\alpha}^{G}(t)+\sum_{\vec{k}}\sum_{\lambda}\hbar\omega_{\vec{k}}a_{\vec{k},\lambda}^{\dagger}a_{\vec{k},\lambda}+\sum_{\vec{k}}\sum_{\lambda}\sum_{\alpha\beta}\left(g_{\vec{k},\alpha\beta}^{\lambda,G}(t)\sigma_{\alpha\beta}^{G}(t)a_{\vec{k},\lambda}+h.c.\right),$
(17)

where,
$\vec{r}^{G}_{\alpha\beta}(t)\equiv\langle\alpha_{G}(t)|\vec{r}|\beta_{G}(t)\rangle$,
$\sigma^{G}_{\alpha\beta}(t)\equiv|\alpha_{G}(t)\rangle\langle\beta_{G}(t)|$
and
$g_{\vec{k},\alpha\beta}^{\lambda,G}(t)=ie\left[\left(\hbar\omega_{\vec{k}}\right)/\left(2\varepsilon_{0}V\right)\right]^{1/2}$
$\left(e_{\vec{k}}^{\lambda}\cdot r_{\alpha\beta}^{G}(t)\right)$. For the
interaction Hamiltonian, we take the rotating wave approximation

${}H_{I}\approx\sum_{\vec{k}}\sum_{\lambda}\sum_{\alpha\beta}\left(g_{\vec{k},\alpha\beta}^{\lambda,G}(t)\sigma_{\alpha\beta}^{G+}(t)a_{\vec{k},\lambda}+g_{\vec{k},\beta\alpha}^{*\lambda,G}(t)\sigma_{\alpha\beta}^{G-}(t)a_{\vec{k},\lambda}^{\dagger}\right),$
(18)

where

${}\sigma_{\alpha\beta}^{G+}(t)=\left\\{\begin{array}[]{rl}\left|\alpha_{G}(t)\right\rangle\left\langle\beta_{G}(t)\right|,&\text{
when }E_{G,\alpha}(t)>E_{G,\beta}(t)\\\ 0,&\text{ when
}E_{G,\alpha}(t)<E_{G,\beta}(t)\end{array}\right.,$ (19)

and

${}\sigma_{\alpha\beta}^{G-}(t)=\left\\{\begin{array}[]{rl}0,&\text{ when
}E_{G,\alpha}(t)>E_{G,\beta}(t)\\\
\left|\alpha_{G}(t)\right\rangle\left\langle\beta_{G}(t)\right|,&\text{ when
}E_{G,\alpha}(t)<E_{G,\beta}(t)\end{array}\right..$ (20)

Considering only the single-photon process, the whole state can be
approximated as

${}|\Psi(t)\rangle\approx\sum_{\alpha}b_{\alpha}(t)|\alpha_{G}(t),0\rangle+\sum_{\alpha,\vec{k}}\sum_{\lambda}b_{\alpha_{G},\vec{k}}^{\lambda}(t)\left|\alpha(t),\gamma_{\vec{k}}^{\lambda}\right\rangle,$
(21)

where, $|\alpha_{G}(t),0\rangle=|\alpha_{G}(t)\rangle\otimes|0\rangle$ and
$|\alpha(t),\gamma_{\vec{k}}^{\lambda}\rangle=|\alpha(t)\rangle\otimes|\gamma^{\lambda}_{\vec{k}}\rangle$.
$|n_{G}(t)\rangle$ is the instantaneous eigen-state of the relativistic atom
as we constructed in the previous section, and $|0\rangle$,
$|\gamma_{\vec{k}}^{\lambda}\rangle$ are the Fock states of the photon. In
order to derive the coefficients $b_{\alpha}$ and
$b^{\lambda}_{\alpha,\vec{k}}$, we introduce the dynamic phase
$\theta_{\alpha}=-\int_{t_{i}}^{t_{f}}E_{\alpha}(s)ds$ and adiabatic phase
$\gamma_{\alpha}=i\int_{t_{i}}^{t_{f}}\langle\alpha(s)\mid\dot{\alpha}(s)\rangle$
and write

$\displaystyle b_{\alpha}(t)$
$\displaystyle=\exp\left[i\left(\gamma_{\alpha}(t)+\theta_{\alpha}(t)\right)\right]c_{\alpha}(t),$
(22) $\displaystyle b_{\alpha,\vec{k}}^{\lambda}(t)$
$\displaystyle=\exp\left[i\left(\gamma_{\alpha}(t)+\theta_{\alpha}(t)\right)-i\omega_{\vec{k}}\left(t-t_{i}\right)\right]c_{\alpha,\vec{k}}^{\lambda}(t),$

Using the state-evolution equation
$i\partial_{t}|\psi(t)\rangle=H_{\text{tot}}(t)|\psi(t)\rangle$ and neglecting
the transition terms which conserve the photon number because the non-uniform
electrostatic field adiabatically affects relativistic atoms in atomic self-
reference frame, we obtain differential equations for the new coefficients
$c_{\alpha}$ and $c^{\lambda}_{\alpha,\vec{k}}$ given by

$\displaystyle\dot{c}_{\alpha}(t)$
$\displaystyle=-i\frac{1}{\hbar}\sum_{\vec{k},\lambda}\overline{\sum}_{\beta}g_{\vec{k},\alpha\beta}^{\lambda,G}(t)\exp\left(i\left(\gamma_{\beta}(t)+\theta_{\beta}(t)-\left(\gamma_{\alpha}(t)+\theta_{\alpha}(t)\right)\right)-i\omega_{\vec{k}}\left(t-t_{i}\right)\right)c_{\beta,\vec{k}}^{\lambda}(t),$
(23) $\displaystyle\dot{c}_{\beta,\vec{k}}^{\lambda}(t)$
$\displaystyle=-i\frac{1}{\hbar}g_{\vec{k},\beta\alpha}^{*\lambda,G}(t)\exp\left(-i\left(\gamma_{\beta}(t)+\theta_{\beta}(t)-\left(\gamma_{\alpha}(t)+\theta_{\alpha}(t)\right)\right)+i\omega_{\vec{k}}\left(t-t_{i}\right)\right)c_{\alpha}(t),$

where, $\overline{\sum}_{\beta}$ stands for the sum of all terms satisfying
$E_{\beta}(t)<E_{\alpha}(t)$. To handle the decay of the excited state, we
follow a method similar to the Weisskopf-Wigner approximation Weisskopf30
viz.

$\displaystyle\dot{c}_{\alpha}(t)$
$\displaystyle=-\frac{1}{\hbar^{2}}\sum_{\vec{k},\lambda}\overline{\sum}_{\beta}g_{\vec{k},\alpha\beta}^{\lambda,G}(t)\exp\left(i\left(\gamma_{\beta}(t)+\theta_{\beta}(t)-\left(\gamma_{\alpha}(t)+\theta_{\alpha}(t)\right)\right)-i\omega_{\vec{k}}\left(t-t_{i}\right)\right)$
(24)
$\displaystyle\left.\times\int_{t_{i}}^{t}g_{\vec{k},\beta\alpha}^{*\lambda,G}\left(t^{\prime}\right)\exp\left(-i\left(\gamma_{\beta}\left(t^{\prime}\right)+\theta_{\beta}\left(t^{\prime}\right)-\left(\gamma_{\alpha}\left(t^{\prime}\right)+\theta_{\alpha}\left(t^{\prime}\right)\right)\right)+i\omega_{\vec{k}}\left(t^{\prime}-t_{i}\right)\right)c_{\alpha}\left(t^{\prime}\right)\right).$

To proceed with the computation in a clearer form, we introduce
$\Delta_{G,\alpha\beta}(t)=[E_{\alpha,G}(t)/\hbar-i\langle\alpha_{G}(t)|\dot{\alpha}_{G}(t)\rangle]-[E_{\beta,G}(t)/\hbar-i\langle\beta_{G}(t)|\dot{\beta}_{G}(t)\rangle]$.
Since the radiation spectrum typically has a peak frequency, we can
effectively perform the frequency integration in the range
$(-\infty,+\infty)$, yielding

${}\int_{-\infty}^{\infty}d\omega_{k}e^{i\left(\omega-\omega_{\vec{k}}\right)\left(t-t^{\prime}\right)}=2\pi\delta\left(t-t^{\prime}\right).$
(25)

Using Eqs. (24) and (25), we thereby obtain

${}\dot{c}_{\alpha}(t)=-\sum_{\lambda}\overline{\sum}_{\beta}\int\frac{\Delta_{G,\alpha\beta}(t)^{3}\mathrm{~{}d}\cos\theta\mathrm{d}\phi}{16\hbar\pi^{2}c^{3}\varepsilon_{0}}e^{2}\left|\vec{e}_{\vec{k}}^{\lambda}\cdot\vec{r}_{\alpha\beta}^{G}(t)\right|^{2}c_{\alpha}(t).$
(26)

See Appendix B for the detail of the calculation process of
$\left|\vec{e}_{\vec{k}}^{\lambda}\cdot\vec{r}_{\alpha\beta}^{G}(t)\right|^{2}$.
Hence

${}\dot{c}_{\alpha}(t)=-\overline{\sum}_{\beta}\frac{e^{2}\Delta_{G,\alpha\beta}(t)^{3}\left|\vec{r}^{G}_{\alpha\beta}(t)\right|^{2}}{6\hbar\pi
c^{3}\varepsilon_{0}}c_{\alpha}(t).$ (27)

Then, to obtain the solutions for $c_{\alpha}(t)$ and
$c^{\lambda}_{\beta,\vec{k}}(t)$, we combine Eq. (27) with the second line of
Eqs. (23), compute the relevant matrix elements, and perform the numerical
integration step by step. Ultimately, the to-be-observed radiation spectrum,
accumulated until a time $t_{f}$, is computed via

$\displaystyle\int
d\omega_{\vec{k}}S_{\beta}\left(\omega_{\vec{k}},t_{f}\right)$
$\displaystyle=\sum_{\lambda}\int\frac{d^{3}k}{(2\pi)^{3}}\left|c^{\lambda}_{\beta,\vec{k}}\left(t_{f}\right)\right|^{2}$
(28)
$\displaystyle=\sum_{\lambda}\int\frac{\omega_{\vec{k}}^{2}\mathrm{d}\omega_{\vec{k}}\mathrm{d}\cos\theta\mathrm{d}\phi}{(2\pi)^{3}c^{3}}\left|c_{\beta,\vec{k}}^{\lambda}(t_{f})\right|^{2},$

where, $S_{\beta}$ is the spontaneous radiation spectrum of the transition
process $(53,0)\rightarrow\beta$. The selection rules of the transition
process can be obtained: $\Delta m=0,\pm 1$. Here, the process of $\Delta m=0$
is called “$\pi$-type” spontaneous radiation. The process of $\Delta m=\pm 1$
is called “$\sigma$-type” spontaneous radiation. According to Eq. (28), the
results of the spontaneous radiation spectrum of the maximum branching ratio
process in the two types can be obtained

${}S_{\beta}\left(\omega_{\vec{k}},t_{f}\right)=\sum_{\lambda}\int\frac{\omega_{\vec{k}}^{2}\mathrm{~{}d}\cos\theta\mathrm{d}\phi}{(2\pi)^{3}c^{3}}\left|c_{\beta,\vec{k}}^{\lambda}(t_{f})\right|^{2},$
(29)

where, $\beta=(43,0)$ or $(34,\pm 1)$.

Figure 6 shows the result of the “$\pi$-type” spontaneous radiation spectrum
of the maximum branching ratio process $(53,0)\rightarrow(43,0)$ under
multipolar gauge. Figure 7 shows the result of the “$\pi$-type” spontaneous
radiation spectrum of the process $(53,0)\rightarrow(43,0)$ under three gauge
conditions near the peak frequency. Figure 8 and Fig. 9 show the result of the
“$\sigma$-type” spontaneous radiation spectrum of the maximum branching ratio
process $(53,0)\rightarrow(34,\pm 1)$ under multipolar gauge. Figure 10 and
Fig. 11 show the result of the “$\sigma$-type” spontaneous radiation spectrum
of the process $(53,0)\rightarrow(34,\pm 1)$ under three gauge conditions near
the peak frequency.

Figure 6: The result of the spontaneous radiation spectrum of the process
$(53,0)\rightarrow(43,0)$ under multipolar gauge. Figure 7: The result of the
spontaneous radiation spectrum of the process $(53,0)\rightarrow(43,0)$ under
three gauge conditions near the frequency $\omega_{0}=1.18825\times 10^{-3}$
$\mathrm{a.u.}$. Figure 8: The result of the spontaneous radiation spectrum of
the process $(53,0)\rightarrow(34,-1)$ under multipolar gauge. Figure 9: The
result of the spontaneous radiation spectrum of the process
$(53,0)\rightarrow(34,1)$ under multipolar gauge. Figure 10: The result of the
spontaneous radiation spectrum of the process $(53,0)\rightarrow(34,-1)$ under
three gauge conditions near the frequency $\omega_{0}=1.2207\times 10^{-3}$
$\mathrm{a.u.}$. Figure 11: The result of the spontaneous radiation spectrum
of the process $(53,0)\rightarrow(34,1)$ under three gauge conditions near the
frequency $\omega_{0}=1.21170\times 10^{-3}$ $\mathrm{a.u.}$.

In Fig. 6, Fig. 8, and Fig. 9, we plot the unperturbed spectrum in the absence
of non-uniform electrostatic field generated by charged ring. We find that the
spontaneous radiation spectrum is obviously different from the spontaneous
radiation spectrum of the stationary state. The peak frequency corresponding
to the maximum radiation intensity under multipolar gauge is about $1.4\times
10^{-7}$ $\mathrm{a.u.}$ for “$\pi$-type” and about $5.4\times 10^{-7}$
$\mathrm{a.u.}$ for “$\sigma$-type” higher than that of unperturbed spectrum.
We also find that under the influence of non-uniform electrostatic field, the
radiation intensity of “$\pi$-type” is enhanced, while that of “$\sigma$-type”
is significantly weakened.

We compare the peak frequency corresponding to the maximum radiation intensity
under three gauge conditions, and find that the peak frequency corresponding
to the maximum radiation under the Coulomb gauge is about $9\times 10^{-12}$
$\mathrm{a.u.}$ for “$\pi$-type” and about $8\times 10^{-12}$ $\mathrm{a.u.}$
for “$\sigma$-type” higher than that under the multipolar gauge and about
$3.9\times 10^{-11}$ $\mathrm{a.u.}$ for “$\pi$-type” and about $3.8\times
10^{-11}$ $\mathrm{a.u.}$ for “$\sigma$-type” lower than that under Lorentz
gauge.

##  V Summary and discussion

In this paper, we proposed a physical system in which a relativistic atomic
beam passes through the non-uniform electrostatic field produced by a charged
ring. We calculated the spontaneous radiation spectrum of our physical system
and found gauge dependence which cannot be cured by existing methods. The
obtained peak frequency can differ by about $413$ $\mathrm{KHz}$ or larger for
the commonly used Coulomb, Lorentz, and multipolar gauges by using SI units.
Compared with the physical system we studied before Chen22 , there are mature
experimental methods to prepare the physical system we need. This provides a
feasible experimental scheme for us to test the self-consistency of gauge
theory under the external field approximation method.

Since we calculate the spontaneous emission in the atomic self-reference
frame, one has to transform measurement results from the laboratory reference
frame to this frame when comparing the measurement with our calculation.
Although there are some risks in defining quantum states by using semi-
classical Hamiltonian (2) obtained by the external field approximation method
in the atomic self-reference frame, We can still get the definition of
effective external field potential by fitting experimental data in the absence
of complete quantum field theory. This effective external potential will help
us to understand the physical significance of the gauge potential.

Meanwhile, the design we proposed also makes a tunable system to study the
relativistic bound-state problems. Because the quantum state of the
relativistic atom is defined in the atomic self-reference frame, in our
design, the acting charged ring is relativistic. Moreover, the design makes
both the relativity parameter $\beta=v/c$ and the interaction strength
continuously adjustable.

Our design also has value in practical application. For example, atomic clocks
in satellites are often influenced by the flux of charged particles from the
universe. In this case, it has become a realistic problem to determine which
gauge should be applied to ensure the accuracy of timing.

###### Acknowledgements.

The work had been, and was partly supported by the China NSF via Grants No.
11275077 and No. 11535005.

## Author Contribution Statement

X.-N. Chen and Y.-H. Luo are the co-first authors of this paper. X.-N. Chen
and X.-S. Chen are co-corresponding authors of this paper.

##  Appendix A The scalar potential in different gauges

For Lorentz gauge,

$\phi_{L}=\frac{Q}{4\pi\varepsilon_{0}\sqrt{(z+vt)^{2}+(1-\beta^{2})R^{2}}}.$
(30)

For Coulomb guage,

$\phi_{C}=\frac{Q}{4\pi\varepsilon_{0}\sqrt{(z+vt)^{2}+R^{2}}}.$ (31)

At $t$, the expanded forms of these two scalar potentials are

$\phi_{L}\approx\frac{Q}{4\pi\varepsilon_{0}}[\frac{1}{\sqrt{(vt)^{2}+(1-\beta^{2})R^{2}}}-\frac{(vt)z}{\sqrt{(vt)^{2}+(1-\beta^{2})R^{2}}^{3}}+\cdots],$
(32)

and

$\phi_{C}\approx\frac{Q}{4\pi\varepsilon_{0}}[\frac{1}{\sqrt{(vt)^{2}+R^{2}}}-\frac{(vt)z}{\sqrt{(\beta
t)^{2}+R^{2}}^{3}}+\cdots].$ (33)

By using PZW transformation Power59 ; Wolley71 , we obtain the result of the
multipolar gauge

$\phi_{M}\approx\frac{Q}{4\pi\varepsilon_{0}}[-\frac{(1-\beta^{2})(vt)z}{\sqrt{(vt)^{2}+(1-\beta^{2})R^{2}}^{3}}+\cdots].$
(34)

##  Appendix B Calculation of
$\left|\vec{e}_{\vec{k}}^{\lambda}\cdot\vec{r}_{\alpha\beta}^{G}(t)\right|^{2}$

The constraint relationship between the momentum and the polarization vector
of the physical photon state is $\vec{k}\cdot\vec{e}^{\lambda}_{\vec{k}}=0$.
For example, for a photon with momentum $(0,0,k)$, its polarization vector is
$\vec{e}^{\pm}_{(0,0,k)}=\frac{1}{\sqrt{2}}(1,\pm i,0)$. For a photon with
arbitrary momentum $k(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$, its
polarization vector can be determined by the three-dimensional represatation
of $\mathrm{SO}(3)$ $R$ is obtained. The matrix representation in $R$ is

$R=\begin{pmatrix}\sin\phi&\cos\theta\cos\phi&\sin\theta\cos\phi\\\
-\cos\phi&\cos\theta\sin\phi&\sin\theta\sin\phi\\\
0&-\sin\theta&\cos\theta\end{pmatrix}.$ (35)

Hence,

$\vec{e}^{\pm}_{\vec{k}}=\vec{e}^{\pm}_{(0,0,k)}R^{T}=\frac{1}{\sqrt{2}}(\sin\phi\pm
i\cos\theta\cos\phi,-\cos\phi\pm i\cos\theta\sin\phi,\mp i\sin\theta).$ (36)

Then, we obtain

$\vec{e}^{\pm}_{\vec{k}}\cdot\vec{r}^{G}_{\alpha\beta}(t)=\frac{1}{\sqrt{2}}[(\sin\phi\pm
i\cos\theta\cos\phi)x^{G}_{\alpha\beta}(t)+(-\cos\phi\pm
i\cos\theta\sin\phi)y^{G}_{\alpha\beta}(t)\mp i\sin\theta
z^{G}_{\alpha\beta}(t)].$ (37)

Hence,

$\displaystyle|\vec{e}^{\pm}_{\vec{k}}\cdot\vec{r}^{G}_{\alpha\beta}(t)|^{2}$
$\displaystyle=\frac{1}{\sqrt{2}}[(\sin\phi\pm
i\cos\theta\cos\phi)x^{G}_{\alpha\beta}(t)+(-\cos\phi\pm
i\cos\theta\sin\phi)y^{G}_{\alpha\beta}(t)\mp i\sin\theta
z^{G}_{\alpha\beta}(t)]$ (38)
$\displaystyle\cdot\frac{1}{\sqrt{2}}[(\sin\phi\mp
i\cos\theta\cos\phi)x^{G}_{\alpha\beta}(t)+(-\cos\phi\mp
i\cos\theta\sin\phi)y^{G}_{\alpha\beta}(t)\pm i\sin\theta
z^{G}_{\alpha\beta}(t)]$
$\displaystyle=\frac{1}{2}[(\sin^{2}\phi+\cos^{2}\theta\cos^{2}\phi)(x^{G}_{\alpha\beta}(t))^{2}-2\sin^{2}\theta\cos\phi\sin\phi(x^{G}_{\alpha(t)})(y^{G}_{\alpha\beta}(t))$
$\displaystyle+(\cos^{2}\phi+\cos^{2}\theta\sin^{2}\phi)(y^{G}_{\alpha\beta}(t))^{2}-2\cos\theta\sin\theta\sin\phi(y^{G}_{\alpha\beta}(t))(z^{G}_{\alpha\beta}(t))$
$\displaystyle+\sin^{2}\theta(z^{G}_{\alpha\beta}(t))^{2}-2\cos\theta\sin\theta\cos\phi(z^{G}_{\alpha\beta}(t))(x^{G}_{\alpha\beta}(t))].$

The integral of the cross term coefficient $f(\theta,\phi)$ to the solid angle
is $0$, while the integarl of the square term coefficient $f(\theta,\phi)$ to
the solid angle is $4\pi/3$.

## References

  * (1) W. E. Lamb, Phys. Rev. 85, 259 (1952)
  * (2) F. Bassani, J. J. Forney, A. Quattropani, Phys. Rev. Lett. 39, 1070(1977)
  * (3) D. H. Kobe, Phys.Rev. Lett., 40, 538(1978)
  * (4) S. Olariu, I. Popescu, C. B. Collins, Phys. Rev. D 20, 3095(1979)
  * (5) A. F. Starace, Phys. Rev. A 3, 1242(1971)
  * (6) M. T. Anderson, F. Weinhold, Phys. Rev. A 10, 1457(1974)
  * (7) D. L. Lin, Phys. Rev. A 17, 1939(1978)
  * (8) D. H. Kobe, Phys. Rev. A 19, 205(1979)
  * (9) J. Bauer, Phys. Rev. A 73, 023421(2006)
  * (10) F. Faisa,.J. Phys. B: At., Mol. and Opt.40, F145(2007)
  * (11) F. H. M. Faisal, Phys. Rev. A 75, 063412(2007)
  * (12) Y. V. Vanne, A. Saenz, Phys. Rev. A 79, 023421(2009)
  * (13) R. H. Dicke, Phys. Rev. 93, 99(1954)
  * (14) K. Rzażewski, K. Wódkiewiczand and W. Żakowicz, Phys. Rev. Lett. 35, 432(1975)
  * (15) J. Keeling, J. Phys.: Cond. Mat., 2007, 19, 295213(2007)
  * (16) D. De Bernardis, T. Jaako, P. Rabl, Phys. Rev. A 97, 043820(2018)
  * (17) B. F. Bayman and F. Zardi, Phys. Rev. C 71, 014904 (2005)
  * (18) C. H. Dasso, M.I. Gallardo, H. M. Sofia and A. Vitturi, Phys. Rev. C 73, 034612 (2006)
  * (19) J. Forney, A. Quattropani, F. Bassani, Il Nuovo Cimento B 37, 78(1977)
  * (20) K. H. Yang, Ann. Phys. 101, 62 (1976)
  * (21) D. H. Kobe, E. C. T. Wen, K. H. Yang, Phys. Rev. D, 26, 1927(1982)
  * (22) W. E. Lamb, R. R. Schlicher, and M. O. Scully, Phys. Rev. A 36, 2763 (1987)
  * (23) N. Funai, J. Louko, and E. Martín-Martínez, Phys. Rev. D 99 065014 (2019)
  * (24) X.-N. Chen, Y.-H. Luo and X.-S. Chen, Euro. J. Phys. D 76, 84(2022)
  * (25) H. C. Bryant, P. A. Lovoi, G. G. Ohlsen, Phys. Rev. Lett. 27, 1628(1971)
  * (26) F. Gross, Phys. Rev. C 26, 2203 (1982); ibid. 26, 2226 (1982) .
  * (27) S. Möller, in: Proc. First Eur. Part. Accel. Conf., ed. S. Tazzari (World Scientific, Singapore 1989), p. 112-114
  * (28) C. J. Herlander, in: Cooler Rings and their Applications, ed. T. Katayama and A. Noda (World Scientific, Singapore 1991), p. 15-19
  * (29) H. Eickhoff, K. Beckert, B. Franczak, B. Franzke, F. Nolden, H. Poth, et al., in: Cooler Rings and their Applications, ed. T. Katayama and A. Noda (World Scientific, Singapore 1991), p. 11 - 14
  * (30) T. Tanabe, K. Noda T. Honma M. Kodaira K. Chida , T. Watanabe, et al., Nucl. Instr. Meth., Phys. Res. A 307, 7 (1991)
  * (31) E. Jaeschke, M. Blum, R. E. Pollock, F. Schmitt, M. Steck, B. Holzer et al., Part. Accel. 32, 97 (1990)
  * (32) R.E. Pollock, in: Cooler Rings and their Applications, ed. T. Katayama and A. Noda (World Scientific, Singapore 1991), p. 4-10
  * (33) X. Ma, X. H. Cai, X. L. Zhu, S. F. Zhang, L. J. Meng, D. C. Zhang et. al., Int. J. Mod. Phys. E 18, 373(2009)
  * (34) Manuel Kaiser, Jens Grimmel, Lara Torralbo-Campo, Markus Mack, Florian Karlewski, Florian Jessen, et. al., Phys. Rev. A, 96, 043401 (2017).
  * (35) H. C. Bryant, D. A. Clark, K. B. Butterfield, C. A. Frost, H. Sharifian, H. Tootoonchi, J. B. et.al., Phys. Rev. A, 27, 2889 (1983).
  * (36) V. F. Weisskopf and E. P. Wigner, Z. Phys. 63, 54 (1930)
  * (37) E. A. Power and S. Zienau. Phil. Trans. R . Soc. Lond. A 251, 427 (1959)
  * (38) R. G. Woolley, Proc. R . Soc. Lond. A 321, 557 (1971)