# Multi-ion frequency reference using dynamical decoupling

Lennart Pelzer  Physikalisch-Technische Bundesanstalt, Bundesallee 100, 38116
Braunschweig, Germany  Kai Dietze  Physikalisch-Technische Bundesanstalt,
Bundesallee 100, 38116 Braunschweig, Germany  Institut fÃ¼r Quantenoptik,
Leibniz UniversitÃ¤t Hannover, Welfengarten 1, 30167 Hannover, Germany  VÃctor
JosÃ© MartÃnez-Lahuerta  Institut fÃ¼r Theoretische Physik, Leibniz
UniversitÃ¤t Hannover, AppelstraÃe 2, 30167 Hannover  Institut fÃ¼r
Quantenoptik, Leibniz UniversitÃ¤t Hannover, Welfengarten 1, 30167 Hannover,
Germany  Ludwig Krinner  Physikalisch-Technische Bundesanstalt, Bundesallee
100, 38116 Braunschweig, Germany  Institut fÃ¼r Quantenoptik, Leibniz
UniversitÃ¤t Hannover, Welfengarten 1, 30167 Hannover, Germany  Johannes
Kramer  Physikalisch-Technische Bundesanstalt, Bundesallee 100, 38116
Braunschweig, Germany  Institut fÃ¼r Quantenoptik, Leibniz UniversitÃ¤t
Hannover, Welfengarten 1, 30167 Hannover, Germany  Fabian Dawel  Physikalisch-
Technische Bundesanstalt, Bundesallee 100, 38116 Braunschweig, Germany
Institut fÃ¼r Quantenoptik, Leibniz UniversitÃ¤t Hannover, Welfengarten 1,
30167 Hannover, Germany  Nicolas C. H. Spethmann  Physikalisch-Technische
Bundesanstalt, Bundesallee 100, 38116 Braunschweig, Germany  Klemens Hammerer
Institut fÃ¼r Theoretische Physik, Leibniz UniversitÃ¤t Hannover, AppelstraÃe
2, 30167 Hannover  Piet O. Schmidt  Physikalisch-Technische Bundesanstalt,
Bundesallee 100, 38116 Braunschweig, Germany  Institut fÃ¼r Quantenoptik,
Leibniz UniversitÃ¤t Hannover, Welfengarten 1, 30167 Hannover, Germany

(August 28, 2024)

###### Abstract

We present the experimental realization of a continuous dynamical decoupling
scheme which suppresses leading frequency shifts in a multi-ion frequency
reference based on ${}^{40}\mathrm{Ca}^{+}$. By near-resonant magnetic
coupling of the ${}^{2}\mathrm{S}_{1/2}$ and ${}^{2}\mathrm{D}_{5/2}$ Zeeman
sub-levels using radio-frequency dressing fields, engineered transitions with
reduced sensitivity to magnetic-field fluctuations are obtained. A second
stage detuned dressing field reduces the influence of amplitude noise in the
first stage driving fields and decreases 2nd-rank tensor shifts, such as the
electric quadrupole shift. Suppression of the quadratic dependence of the
quadrupole shift to
$3(2)\text{\,}\mathrm{m}\mathrm{Hz}\mathrm{/}\mu\mathrm{m}^{2}$ and coherence
times of $290(20)\text{\,}\mathrm{ms}$ on the optical transition are
demonstrated even within a laboratory environment with significant magnetic
field noise. Besides removing inhomogeneous line shifts in multi-ion clocks,
the demonstrated dynamical decoupling technique may find applications in
quantum computing and simulation with trapped ions by a tailored design of
decoherence-free subspaces.

Decoherence and dephasing of quantum systems through interaction with the
environment pose a major challenge in the implementation of quantum sensors,
quantum computers and quantum simulators [1]. In particular, magnetic and
electric fields can cause significant shifts that induce decoherence and
dephasing due to their spectral noise properties and inhomogeneity across the
quantum system. For example, in optical clocks based on trapped ions electric
field gradients of the trapping potential couple to the electric quadrupole
moment of the clock states causing a quadrupole shift (QPS) [2]. Several
mitigation strategies have been developed to overcome these limitations.
Technological solutions in the form of active stabilization and passive
shielding have been demonstrated for magnetic fields in a variety of quantum
systems [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]. In
combination with permanent magnets, second long coherence times between
magnetically sensitive levels could be observed, e.g. in trapped calcium ions
[20]. However, most of these technical measures increase size and complexity
of a setup, or can not be eliminated by shielding, such as the QPS. In atomic
clocks, mitigation strategies based on atomic properties include choosing
magnetic field insensitive Zeeman levels, or averaging over several Zeeman
components, which at the same time eliminates the QPS [21, 22, 23, 24]. The
latter can also be eliminated by averaging the transition frequencies of a
single transition with three mutually orthogonal magnetic fields applied to
the system [2, 22, 25, 26], or by applying the magnetic field at an angle of
$54.7^{\circ}$ to the electric field gradient [27, 28]. However, many of these
strategies fail for multi-ion clocks that aim at improving the signal-to-noise
ratio of optical ion clocks by employing more than one ion [29, 30, 31, 32,
33], since the shifts are inhomogeneous across the ion crystal, resulting in
an effective broadening of the observed clock transition. A solution is
provided by dynamical decoupling (DD), a technique that dates back to Hahn
spin echoes in nuclear magnetic resonance spectroscopy [34].

Since then it has been further developed using pulsed [35, 36, 35] and
continuous schemes [36, 37, 38, 29, 39, 40, 41, 42, 43, 44, 45, 46], and
employed experimentally in various systems to eliminate noise sources in
quantum simulations [47], quantum computing [48, 49, 50, 51, 52, 53, 54, 55,
56, 57, 58, 59], quantum memories [60, 61, 62, 63, 64, 65, 66, 67, 68, 69,
70], and quantum sensing [71, 72, 73, 74, 75, 71, 76]. In the context of
precision spectroscopy and optical clocks, DD is implemented, e.g. by changing
the magnetic field within a Ramsey sequence [77] or through precisely timed
radio-frequency (rf) pulses that populate different Zeeman states to cancel
undesired shifts within each individual interrogation. The latter has been
demonstrated for QPS and magnetic field shift in ${}^{88}\mathrm{Sr}^{+}$ [78]
and ${}^{176}\mathrm{Lu}^{+}$ [79], as well as employed to eliminate magnetic
field shifts in a precision measurement of local Lorentz violation [80, 81].
These pulsed schemes can suppress noise sources that vary on time scales
slower than the time between pulses.

Here, we demonstrate the suppression of linear Zeeman shifts as well as the
QPS on an optical transition in a trapped ${}^{40}\mathrm{Ca}^{+}$ frequency
reference using a cascaded continuous dynamical decoupling (CDD) scheme with
up to five ions. In contrast to similar pulsed schemes [78, 79], here the
levels are continuously coupling with rf-fields to engineer an ensemble of
robust artificial transitions by choosing adequate driving parameters [29,
82]. This suppresses inhomogeneous shifts across the ion crystal and mitigates
noise up to half the Fourier frequencies corresponding to the rf Rabi
frequency coupling of the Zeeman states. These artificial transitions can be
treated similar to bare transitions, so e.g., Rabi and Ramsey spectroscopy
sequences as well as entangling operations on motional sidebands are possible
[82].

In the following we briefly summarize the CDD technique employed here [29, 82]
and which are shown schematically in Fig 1. The Zeeman manifolds of ground
${}^{2}\mathrm{S}_{1/2}$ and excited ${}^{2}\mathrm{D}_{5/2}$ levels are
energetically split by a static external magnetic field with magnitude $B_{0}$
in $z$-direction

$H_{0}=\sum_{i=S,D}g_{i}\mu_{B}B_{0}J_{z}^{i}=\sum_{i=S,D}\omega_{0}^{i}J_{z}^{i}$
(1)

Here, the gyromagnetic factor $g_{i}$ and spin operator $z$-component
$J_{z}^{i}$ for each ensemble
$i\in[\text{S}_{\frac{1}{2}},\text{D}_{\frac{5}{2}}]$ as well as the Bohr
magneton $\mu_{B}$ are used. These so-called bare atomic levels are coupled by
near-resonant rf-fields of the form

$H_{\text{rf}}=\sum_{i=S,D}g_{i}\Omega_{1}^{i}\cos(\omega_{1}^{i}t)J_{x}^{i}\\\
-g_{i}\Omega_{2}^{i}\sin(\omega_{1}^{i}t)\cos(\omega_{2}^{i}t)J_{x}^{i}\;.$
(2)

The first drive (first line in Eq. 2)) with frequency
$\omega_{1}^{i}=\omega_{0}^{i}-\Delta_{1}^{i}$ and amplitude $\Omega_{1}^{i}$
produces a ladder of dressed states with frequency splitting

$\bar{\omega}_{0}^{i}=\sqrt{(g_{i}\Omega_{1}^{i}/2)^{2}+(\Delta_{1}^{i})^{2}}\;.$
(3)

The second stage drive (2nd line in Eq. 2) with frequency
$\omega_{2}^{i}=\bar{\omega}_{0}^{i}-\Delta_{2}^{i}$ and amplitude
$\Omega_{2}^{i}$ couples the first stage dressed states analogously to the
first drive if the hierarchy of coupling frequencies
$\omega_{2}^{i}\ll\omega_{1}^{i}$ is fulfilled. Here, colinear alignment of
all rf-drives is assumed for simplicity. For a detailed derivation of the
general geometry we refer the reader to [82]. The two-stage design serves two
purposes. Firstly, temporal or spatial drive strength variations
$\delta\Omega_{1}^{i}(\vec{r},t)$ will result in undesired energy shifts of
the resulting states. Choosing a weaker second drive can mitigate these
shifts, since assuming a constant fractional variation results in smaller
absolute frequency shifts. Secondly, the additional drive adds control
parameters to suppress Zeeman and 2nd-rank tensor shifts simultaneously as
will be shown below.

Figure 1: a) Continuous dynamical decoupling scheme. Near resonant rf-coupling
with driving frequency $\omega_{j}^{i}$ and strength $\Omega_{j}^{i}$ results
in an ensemble of $2(J_{j}^{i}+1)$ dressed states with spacing
$\bar{\omega}_{j}^{i}$ and numbering $m_{j}^{i}$ where $j=1,2$ marks the 1st
or 2nd layer of dressing and $i=S,D$ the Zeeman manifold. Here, the quantum
numbers $J_{1}^{S}=1/2=J_{2}^{S}$ and $J_{1}^{D}=5/2=J_{2}^{D}$ match these of
the spin of the underlying bare atomic states. Cross-coupling and Bloch-
Siegert shifts result in a shift of the bare states’ transition frequency
$\omega_{0}^{j}\rightarrow\omega_{0}^{j\\#}$. The doubly-dressed transitions
are probed with $729\text{\,}\mathrm{nm}$ light at frequency $\nu_{L}$ and
coupling strength $\Omega_{L}$. Adapted from [83] b) Sketch of the
spectroscopy sequence with typical durations. The lower part shows a sketch of
the employed CDD waveforms with signal amplitudes $y_{i}$. Two successive
frequency sweeps transfer the population adiabatically from the
${}^{2}\mathrm{S}_{1/2}$ via the first dressed to the doubly-dressed target
state.

The so-called mixing angle
$\cos(\theta_{j}^{i})=\Delta_{j}^{i}/\bar{\omega}_{j}^{i}$ with $j\in[1,2]$
and $i\in[S_{\frac{1}{2}},D_{\frac{5}{2}}]$ determines the suppression of
these shifts for each stage and manifold. In a doubly-dressed basis, the
response of the system to sufficiently slow magnetic field variations $\delta
B(\vec{r},t)$ can be described by

$V_{\delta
B}=\sum_{i=S,D}\cos(\theta_{1}^{i})\cos(\theta_{2}^{i})g_{i}\mu_{B}\delta
B(\vec{r},t)J_{z}^{i}$ (4)

and the quadrupole shift of the excited level by [82]

$V_{Q}=\left(1-3\cos^{2}(\theta_{1}^{D})\right)\left(1-3\cos^{2}(\theta_{2}^{D})\right)\\\
\times\Theta\frac{3}{8}\frac{\partial E_{z}}{\partial
z}\left(J^{D}(J^{D}+1)-3(J_{z}^{D})^{2}\right)\,,$ (5)

using the spin of the D-level $J^{D}$. The magnitude of the shift depends on
the quadrupole moment $\Theta$ and the electric field gradient $\partial
E/\partial z$ [2], which is given by the joint potential of the trap and
neighbouring ions. It is worth noting that the suppression of both shifts is
independent for both stages. This allows e.g. suppression of the Zeeman shift
using a resonant drive with $\theta_{1}^{S}=\theta_{2}^{S}=\theta_{1}^{D}=0$
for all drives except the 2nd-stage $D$-level drive. QPS suppression can be
attained using the so-called magic mixing angle on this last drive with
$\cos^{2}(\theta_{2}^{D})=\frac{1}{3}$ at the expense of a reduced Zeeman
shift suppression. The achievable degree of suppression depends on the
deviation of the angles $\theta_{j}^{i}$ from their ideal values. The
detunings $\Delta_{j}^{i}$ are affected by magnetic field fluctuations acting
on the bare atomic states, the QPS of the D states, frequency shifts from off-
resonant cross-coupling of the S drive on the D level and vice versa, as well
as the correction of counter-rotating terms in the rotating wave approximation
(Bloch-Siegert shift) [29, 82]. Since $\theta_{j}^{i}$ depends inversely
proportional on $\Omega_{j}^{i}$, the scheme becomes more robust against
fluctuations in magnitude and noise frequency of these shifts for higher
coupling strength at the expense of larger frequency shifts caused by
amplitude noise of the rf-drives, assuming a fixed fractional stability. The
parameter set for optimum suppression is thus a trade-off between detuning
fluctuations and driving strength, requiring a detailed knowledge of the noise
properties of both.

In the setup presented here, the magnetic field noise causes line broadening
on the order of $\nu_{\delta B}\lesssim$1\text{\,}\mathrm{kHz}$$ with a
timescale on the order of milliseconds given by the mains supply and its
harmonics. In contrast, the QPS is near static with an inhomogeneous shift of
$\nu_{\text{QPS}}\lesssim$100\text{\,}\mathrm{Hz}$$ for typical trapping
potentials. Assuming a realistic mixing angle mismatch of $\delta\theta=1\%$,
a Zeeman shift suppression on the order of
$\cos(\theta_{1}^{i})\cos(\theta_{2}^{i})\approx 10^{-4}$ (compare Eq. 4) for
a doubly resonant drive and $\left(1-3\cos^{2}(\theta_{2}^{D})\right)\approx
3\times 10^{-2}$ QPS suppression (Eq. 5) can be expected for small deviations.
The sensitivities of the dressed states to magnetic field changes can be fine-
tuned separately for S and D-levels, such that a residual linear sensitivity
due to the magic-angle detuned second D-drive can be partially compensated
with a detuned second S-drive.

The majority of the experimental setup is described in detail elswhere [84,
83], therefore only key components are summarized here.
${}^{40}\mathrm{Ca}^{+}$ ions are confined in a segmented Paul trap [85] with
low excess micromotion even for axially extended crystals [84, 30] as well as
low motional heating due to on-chip filtering and a low noise voltage supply
[86]. All lasers needed for cooling, detection and state preparation are
locked to a wavemeter111High Finesse WS U10 and most beams are supplied via
fibre-coupled components to obtain a robust setup with small footprint. The
individual beams are switched and frequency steered by acousto-optic
modulators (AOM) that are controlled by a pulse sequencer [84]. Light shifts
during probing of the clock transition are minimized by using shutters in all
beam paths.

The interrogation laser for the
${}^{2}\mathrm{S}_{1/2}\leftrightarrow^{2}\mathrm{D}_{5/2}$ transition at
$729\text{\,}\mathrm{n}\mathrm{m}$ is based on an amplified extended cavity
diode laser222TA pro Toptica, which is first pre-stabilized to a
$10\text{\,}\mathrm{c}\mathrm{m}$ cavity resulting in a fractional linewidth
of $170(4)\text{\,}\mathrm{Hz}$. Further reduction of the laser linewidth is
established by transfer locking [87, 88] to a highly stable reference cavity
[89] via a fibre-based frequency comb. An independent measurement revealed an
upper limit for the fractional frequency instability of the
${}^{40}\mathrm{Ca}^{+}$ clock laser against the reference laser of
$\sigma_{y}<2\times 10^{-15}$ for $\tau>\qty{1}{}$. The reference laser
reaches typically an instability of $5\times 10^{-17}$ for up to
$10\text{\,}\mathrm{s}$ and a drift of $50\text{\,}\mathrm{\SIUnitSymbolMicro
Hz}\text{\,}{\mathrm{s}}^{-1}$. Therefore, the employed laser is not limiting
for the following measurements.

A single aspheric lens is used for collecting fluorescence light of the ions.
It is focused simultaneously onto a photomultiplier tube (PMT) and a
scientific complementary metal-oxide-semiconductor (sCMOS) camera using a beam
splitter. The PMT allows for discrimination of bright and dark states within
$100\text{\,}\mathrm{\SIUnitSymbolMicro s}$ detection time [84]. Using the
camera, spatially resolved detection of single ion excitation in a linear ion
crystal is possible at the cost of a longer detection time of
$10\text{\,}\mathrm{ms}$. For this, the counts in a region of interest around
each ion are analysed.

Three pairs of orthogonal magnetic field coils generate a static offset field
of 357 along the axial trap direction resulting in 10 splitting of the
${}^{2}\mathrm{S}_{1/2}$ levels and an average quadratic Zeeman shift of
$1.37\text{\,}\mathrm{Hz}$ [23]. The axial magnetic field component is
actively stabilized against slow drifts using a fluxgate sensor and feedback
to the current of the coils. Mains power synchronous magnetic field noise
limits the linewidth of the
${}^{2}\mathrm{S}_{1/2}\leftrightarrow^{2}\mathrm{D}_{5/2}$ clock transition
typically to $\Delta\nu\approx$100\text{\,}\mathrm{Hz}$$. In addition, on
timescales of ten seconds fluctuations on the order of $|\Delta
B_{0}|\approx$5\text{\,}\mathrm{nT}$$ are observed.

Resonant magnetic field coils are employed for the CDD to shape and enhance
the generated coupling waveforms. These rf-coils consist of two separate LCR
circuits with tunable capacitors to match the resonance frequency of the
Zeeman ensembles. The quality factors $Q_{S}=8.59(3)$, $Q_{D}=15.95(4)$ of the
coils are chosen as a compromise between required field strength and
minimizing signal distortion by its transfer function. Duty cycle and
laboratory temperature change the temperature-dependent transfer function of
the coils, resulting in a variation in coupling strength. Therefore, the
capacitive and resistive parts are mounted on a temperature-controlled
baseplate and the coils are designed to minimize warm-up during operation.

The mounts are placed on translation stages and positioned in a distance of
approximately $5\text{\,}\mathrm{cm}$ to the ion crystal in an inverted
viewport outside the vacuum chamber. The signal for each coil is supplied via
an inductively-coupled, impedance-matched in-coupling coil, which is driven by
an amplifier. The signal source is a two-channel arbitrary voltage
generator333Keysight 33622A (AWG) allowing flexible waveforms, e.g. for state
preparation using a rapid adiabatic passage (see below). To avoid distortion
of the two-tone signal by the phase and amplitude response of the coils, their
transfer function is measured and pre-compensated in the AWG waveform by
applying the inverse transfer function to the waveform (compare Eq. S1).

Figure 2: $729\text{\,}\mathrm{nm}$ pulse time scan of a doubly-dressed
artificial transition. The strong Zeeman shift suppression with resonant
driving parameters allows a long coherence time. The data is in good agreement
with a combined Gaussian and exponential decay model
$P_{e}(t_{r})=\frac{1}{2}\left(1-\cos(2\pi\Omega_{L}t_{r})e^{-t_{r}/\Gamma}e^{-t_{r}^{2}/(2\gamma^{2})}\right)$
with given natural decay $\Gamma=$1.168\text{\,}\mathrm{s}$$ [90] and
$\gamma=$0.29(2)\text{\,}\mathrm{s}$$. The Gaussian decay shape indicates long
noise correlation times with respect to the duration of an individual data
point [91].

The CDD sequence begins with the preparation of the target state in the
dressed S-levels. This is achieved through an adiabatic ramp from the bare
ground-state to the dressed level by applying a tailored waveform to the
S-coils (see supplementary material Eqs. S2 \- S5). A continuous transition
without phase jumps from the preparation to the final second-stage waveform is
crucial to maximize the contrast. The D-level waveform can be switched on
instantaneously, since the D state is not populated and the cross-driving
effect of the D-level drive onto the population in the S-level is negligible.
Spectroscopy within the CDD spectrum with the ${}^{2}\mathrm{S}_{1/2}$ to
${}^{2}\mathrm{D}_{5/2}$ laser at $729\text{\,}\mathrm{nm}$ is performed after
the applied rf-drives have reached their steady-state values. As long as the
optical coupling is weak compared to the rf-coupling $\Omega_{L}\ll\Omega_{i}$
the fast dynamics of the manifold averages out and the artificial transitions
can be treated analogously to a natural ${}^{40}\mathrm{Ca}^{+}$ clock
transition. In the following the
$\ket{S_{1/2},m_{0}^{S}=-\frac{1}{2},m_{1}^{S}=\frac{1}{2},m_{2}^{S}=\frac{1}{2}}\leftrightarrow\ket{D_{5/2},m_{0}^{D}=-\frac{3}{2},m_{1}^{D}=\frac{1}{2},m_{2}^{D}=\frac{1}{2}}$
transition (compare Fig. 1) is probed with the clock laser directed along the
static magnetic field and trap axis. With the parameters given in Table S1, a
magnetic field shift sensitivity of
$-1.2(1)\text{\times}{10}^{-4}\text{\,}\mathrm{Hz}\mathrm{/}\mathrm{n}\mathrm{T}^{2}$
was determined in a separate measurement. This represents an almost four
orders of magnitude improvement over the sensitivity of the bare states of
$5.6\text{\,}\mathrm{Hz}\text{\,}{\mathrm{nT}}^{-1}$ for typical field
variations of $5\text{\,}\mathrm{nT}$ and is comparable to clock species with
magnetic-field sensitivity only from nuclear spin, such as
${}^{27}\mathrm{Al}^{+}$ and ${}^{115}\mathrm{In}^{+}$ [92, 93]. Therefore, a
coherence time of $0.29(2)\text{\,}\mathrm{s}$ can be observed on one of the
artificial transitions (see Fig. 2). The optical coupling strength of each
transition between dressed states depends on the selection rules of the
underlying Zeeman transition as well as on mixing angles and quantum numbers
of the dressed state [82].

Figure 3: Comparison of the relative frequency shifts of the artificial clock
transition along a linear five-ion crystal for different CDD-parameters. The
individual line centers are obtained from a scan of the clock laser with
fitting of individual peaks using ion-resolved camera data. a) Calculated QPS
for five-ion crystal for bare Zeeman levels. Ion positions are simulated in a
pure harmonic potential with given trapping frequencies to extract field
gradients at the ion positions. b) Doubly resonant parameter set
($\Delta_{j}^{i}=0\quad\forall i,j$). Tensorial shifts are of equal size
compared to the calculated $|m_{J}|=\frac{3}{2}$ bare level shift (green
crosses). c) Magic angle detuned second stage D-level CDD drive. The visible
linear dependency arises from magnetic field inhomogeneity in combination with
lower Zeeman shift suppression compared to (b). d) Same CDD drive as for (c),
but optimized offset B-field for field inhomogeneity compensation.

In order to demonstrate the capabilities of this scheme we probe the line
shifts of a linear five-ion crystal. Such a system is a stringent test for a
multi-ion frequency reference. In Fig. 3 the comparison of inhomogeneous line
shifts along the crystal for two different CDD-parameter sets is shown. The
frequencies of the first stage rf-fields are held close to the respective
splitting of the Zeeman states $\Delta_{1}^{S}=0=\Delta_{1}^{D}$. The detuning
of the 2nd D-level stage is set for the two cases $\Delta_{2}^{D}=0$ and
$\Delta_{2}^{D}=1/\sqrt{8}g_{D}\Omega_{2}^{D}$, corresponding to the resonant
and magic mixing angle, respectively. In the first case only Zeeman shifts are
suppressed, the second also suppresses 2nd-rank tensor shifts at the cost of
reduced Zeeman shift suppression. A relative QPS suppression factor of 12 was
measured between resonant and magic detuned parameter set (Fig. 3 b, c),
measured by the quadratic fitting parameter of
$3(2)\text{\,}\mathrm{m}\mathrm{Hz}\mathrm{/}\mu\mathrm{m}^{2}$.
Simultaneously, the linear term increases by a factor of 3, caused by lower
suppression of the magnetic field inhomogeneity. In a third setting, the CDD
parameters are fine-tuned by an offset of the static B-field such that the
imperfections in residual inhomogeneous static and dynamic B-field partially
compensate, which allows to suppress the overall frequency spread across the
$20\text{\,}\mathrm{\SIUnitSymbolMicro m}$ long crystal to
$1.8(7)\text{\,}\mathrm{Hz}$ (Fig. 3 d). Further optimization by variation of
all CDD parameters is likely to improve the results. However, technical
measures for enhanced homogeneity and stability of the fields will loosen the
requirements to find the optimal parameters.

Figure 4: Overlapping Allan deviation of the relative frequency stability of a
linear five-ion crystal. Two-point sampling Rabi interrogation with
$100\text{\,}\mathrm{ms}$ interrogation time and duty cycle of $33\%$ was
used. The excitation probability of the individual ions was measured using an
sCMOS camera. The sensitivity to small dynamic or offset magnetic field
variations depends on the absolute field value. Therefore, inhomogeneous
fields cause stability differences for the individual ions.

The slow temporal variations of the leading frequency shifts, namely residual
amplitude fluctuations of the CDD-drive and of the magnetic field, are
measured to be on the level of $|\frac{\delta\Omega_{i}}{\Omega_{i}}|<6\times
10^{-5}$ and $|\frac{\delta B_{0}}{B_{0}}|<2\times 10^{-5}$, respectively.
These were determined by observing the splitting of bare and first-stage
dressed states over several $100\text{\,}\mathrm{s}$.

The frequency instability of the doubly-dressed transition of a five-ion
crystal is displayed in Fig. 4. With the chosen parameters, a gain in
statistical uncertainty due to increased ion number is given up to
$10\text{\,}\mathrm{s}$ averaging time before residual magnetic field drift
compromise the instability. The Zeeman shift sensitivity depends on the
absolute field at the position of each ion, therefore the stability varies for
the individual ions. However, the short-term instability of $3\times
10^{-15}/\sqrt{\tau}$ is comparable to ${}^{40}\mathrm{Ca}^{+}$ clocks with
multi-layer magnetic shielding [94]. Increasing the interrogation duty cycle
from $T/T_{c}=33\%$ to $\approx 90\%$ will be possible with a faster
experimental control system together with faster camera detection and readout.

In conclusion, we have demonstrated the suppression of linear Zeeman shifts
and QPS on an artificial clock transition through CDD. The presented CDD
scheme is especially useful for transitions with strong electronic linear
Zeeman shift. This lifts the need for elaborate shielding and reduces the
statistical uncertainty in clocks by enabling longer interrogation times. By
reducing magnetic field noise and drive field inhomogeneity, as well as
implementing faster camera detection, a larger number of ions would enable
this frequency reference to be used as a flywheel in composite, multi-species
clocks [95, 96, 97, 98]. To the best of our knowledge, this is the first time
a 2-stage CDD scheme has been successfully implemented, which might prove
useful for other quantum sensors or qubit systems that are limited by magnetic
field sensitivity.

##  I acknowledgments

We thank PTBâs unit-of-length working group for providing the stable silicon
referenced laser source. Fruitful discussions with Nati Aharon, Alex Retzker
and the group of Roee Ozeri are acknowledged. This joint research project was
financially supported by the State of Lower Saxony, Hanover, Germany through
NiedersÃ¤chsisches Vorab and by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) â Project-ID 274200144 â SFB 1227 (DQ-mat,
projects A06 and B03). This project also received funding from the European
Metrology Programme for Innovation and Research (EMPIR) cofinanced by the
Participating 5 States and from the European Unionâs Horizon 2020 research
and innovation programme (Project No. 20FUN01 TSCAC).

## References

  * Acín  _et al._ [2018] A. Acín, I. Bloch, H. Buhrman, T. Calarco, C. Eichler, J. Eisert, Daniel Esteve, N. Gisin, S. J. Glaser, F. Jelezko, S. Kuhr, M. Lewenstein, M. F. Riedel, P. O. Schmidt, R. Thew, A. Wallraff, I. Walmsley, and F. K. Wilhelm, New J. Phys. 20, 080201 (2018).
  * Itano [2000] W. M. Itano, J. Res. NIST 105, 829 (2000).
  * Merkel  _et al._ [2019] B. Merkel, K. Thirumalai, J. E. Tarlton, V. M. Schäfer, C. J. Ballance, T. P. Harty, and D. M. Lucas, Rev. Sci. Instrum. 90, 044702 (2019).
  * Zeng  _et al._ [2023] M. Zeng, Z. Ma, R. Hu, B. Zhang, Y. Hao, H. Zhang, Y. Huang, H. Guan, and K. Gao, Chinese Phys. B 10.1088/1674-1056/acf5d5 (2023).
  * Borkowski  _et al._ [2023] M. Borkowski, L. Reichsöllner, P. Thekkeppatt, V. Barbé, T. van Roon, K. van Druten, and F. Schreck, Rev. Sci. Instrum. 94, 073202 (2023).
  * Duan  _et al._ [2022] Z.-X. Duan, W.-T. Wu, Y.-T. Lin, and S.-J. Yang, Rev. Sci. Instrum. 93, 123201 (2022).
  * Altarev  _et al._ [2014] I. Altarev, E. Babcock, D. Beck, M. Burghoff, S. Chesnevskaya, T. Chupp, S. Degenkolb, I. Fan, P. Fierlinger, A. Frei, E. Gutsmiedl, S. Knappe-Grüneberg, F. Kuchler, T. Lauer, P. Link, T. Lins, M. Marino, J. McAndrew, B. Niessen, S. Paul, G. Petzoldt, U. Schläpfer, A. Schnabel, S. Sharma, J. Singh, R. Stoepler, S. Stuiber, M. Sturm, B. Taubenheim, L. Trahms, J. Voigt, and T. Zechlau, Rev. Sci. Instrum. 85, 075106 (2014).
  * Dijck  _et al._ [2023] E. A. Dijck, C. Warnecke, M. Wehrheim, R. B. Henninger, J. Eff, K. Georgiou, A. Graf, S. Kokh, L. P. Kozhiparambil Sajith, C. Mayo, V. M. Schäfer, C. Volk, P. O. Schmidt, T. Pfeifer, and J. R. Crespo López-Urrutia, Rev. Sci. Instrum. 94, 083203 (2023).
  * Ji  _et al._ [2021] Y.-H. Ji, L. Zhou, S.-T. Yan, C. He, C. Zhou, S. Barthwal, F. Yang, W.-T. Duan, W.-D. Zhang, R.-D. Xu, Q. Wang, D.-X. Li, J.-H. Gao, X. Chen, J. Wang, and M.-S. Zhan, Rev. Sci. Instrum. 92, 083201 (2021).
  * Leopold  _et al._ [2019] T. Leopold, S. A. King, P. Micke, A. Bautista-Salvador, J. C. Heip, C. Ospelkaus, J. R. Crespo López-Urrutia, and P. O. Schmidt, Rev. Sci. Instrum. 90, 073201 (2019).
  * Farolfi  _et al._ [2019] A. Farolfi, D. Trypogeorgos, G. Colzi, E. Fava, G. Lamporesi, and G. Ferrari, Rev. Sci. Instrum. 90, 115114 (2019).
  * Devlin  _et al._ [2019] J. A. Devlin, E. Wursten, J. A. Harrington, T. Higuchi, P. E. Blessing, M. J. Borchert, S. Erlewein, J. J. Hansen, J. Morgner, M. A. Bohman, A. H. Mooser, C. Smorra, M. Wiesinger, K. Blaum, Y. Matsuda, C. Ospelkaus, W. Quint, J. Walz, Y. Yamazaki, and S. Ulmer, Phys. Rev. Applied 12, 044012 (2019).
  * Wodey  _et al._ [2020] E. Wodey, D. Tell, E. M. Rasel, D. Schlippert, R. Baur, U. Kissling, B. Kölliker, M. Lorenz, M. Marrer, U. Schläpfer, M. Widmer, C. Ufrecht, S. Stuiber, and P. Fierlinger, Rev. Sci. Instrum. 91, 035117 (2020).
  * Vovrosh  _et al._ [2018] J. Vovrosh, G. Voulazeris, P. G. Petrov, J. Zou, Y. Gaber, L. Benn, D. Woolger, M. M. Attallah, V. Boyer, K. Bongs, and M. Holynski, Sci. Rep. 8, 2023 (2018).
  * Brandl  _et al._ [2016] M. F. Brandl, M. W. van Mourik, L. Postler, A. Nolf, K. Lakhmanskiy, R. R. Paiva, S. Möller, N. Daniilidis, H. Häffner, V. Kaushal, T. Ruster, C. Warschburger, H. Kaufmann, U. G. Poschinger, F. Schmidt-Kaler, P. Schindler, T. Monz, and R. Blatt, Rev. Sci. Instrum. 87, 113103 (2016).
  * Altarev  _et al._ [2015] I. Altarev, M. Bales, D. H. Beck, T. Chupp, K. Fierlinger, P. Fierlinger, F. Kuchler, T. Lins, M. G. Marino, B. Niessen, G. Petzoldt, U. Schläpfer, A. Schnabel, J. T. Singh, R. Stoepler, S. Stuiber, M. Sturm, B. Taubenheim, and J. Voigt, J. Appl. Phys. 117, 183903 (2015).
  * Lee and Romalis [2008] S.-K. Lee and M. V. Romalis, J. Appl. Phys. 103, 084904 (2008).
  * Mager [1970] A. Mager, IEEE Trans. Magn. 6, 67 (1970).
  * Hanson and Pipkin [1965] R. J. Hanson and F. M. Pipkin, Rev. Sci. Instrum. 36, 179 (1965).
  * Ruster  _et al._ [2016] T. Ruster, C. T. Schmiegelow, H. Kaufmann, C. Warschburger, F. Schmidt-Kaler, and U. G. Poschinger, Appl. Phys. B 122, 254 (2016).
  * Dubé  _et al._ [2005] P. Dubé, A. Madej, J. Bernard, L. Marmet, J.-S. Boulanger, and S. Cundy, Phys. Rev. Lett. 95, 033001 (2005).
  * Margolis  _et al._ [2004] H. S. Margolis, G. P. Barwood, G. Huang, H. A. Klein, S. N. Lea, K. Szymaniec, and P. Gill, Science 306, 1355 (2004).
  * Chwalla  _et al._ [2009] M. Chwalla, J. Benhelm, K. Kim, G. Kirchmair, T. Monz, M. Riebe, P. Schindler, A. Villar, W. Hänsel, C. Roos, R. Blatt, M. Abgrall, G. Santarelli, G. Rovera, and Ph. Laurent, Phys. Rev. Lett. 102, 023002 (2009).
  * Dubé  _et al._ [2013] P. Dubé, A. A. Madej, Z. Zhou, and J. E. Bernard, Phys. Rev. A 87, 023806 (2013).
  * Oskay  _et al._ [2005] W. Oskay, W. Itano, and J. Bergquist, Phys. Rev. Lett. 94, 163001 (2005).
  * Schneider  _et al._ [2005] T. Schneider, E. Peik, and Chr. Tamm, Phys. Rev. Lett. 94, 230801 (2005).
  * Tan  _et al._ [2019] T. R. Tan, R. Kaewuam, K. J. Arnold, S. R. Chanu, Z. Zhang, M. S. Safronova, and M. D. Barrett, Phys. Rev. Lett. 123, 063201 (2019).
  * Steinel  _et al._ [2023] M. Steinel, H. Shao, M. Filzinger, B. Lipphardt, M. Brinkmann, A. Didier, T. E. Mehlstäubler, T. Lindvall, E. Peik, and N. Huntemann, Phys. Rev. Lett. 131, 083002 (2023).
  * Aharon  _et al._ [2019] N. Aharon, N. Spethmann, I. D. Leroux, P. O. Schmidt, and A. Retzker, New J. Phys. 21, 083040 (2019).
  * Keller  _et al._ [2019] J. Keller, T. Burgermeister, D. Kalincev, A. Didier, A. P. Kulosa, T. Nordmann, J. Kiethe, and T. E. Mehlstäubler, Phys. Rev. A 99, 013405 (2019).
  * Arnold  _et al._ [2015] K. Arnold, E. Hajiyev, E. Paez, C. H. Lee, M. D. Barrett, and J. Bollinger, Phys. Rev. A 92, 032108 (2015).
  * Herschbach  _et al._ [2012] N. Herschbach, K. Pyka, J. Keller, and T. E. Mehlstäubler, Appl. Phys. B 107, 891 (2012).
  * Champenois  _et al._ [2010] C. Champenois, M. Marciante, J. Pedregosa-Gutierrez, M. Houssin, and M. Knoop, Phys. Rev. A 81, 10.1103/PhysRevA.81.043410 (2010).
  * Hahn [1950] E. L. Hahn, Phys. Rev. 80, 580 (1950).
  * Viola  _et al._ [1999] L. Viola, E. Knill, and S. Lloyd, Phys. Rev. Lett. 82, 2417 (1999).
  * Valahu  _et al._ [2022] C. H. Valahu, I. Apostolatos, S. Weidt, and W. K. Hensinger, J. Phys. B: At. Mol. Opt. Phys. 55, 204003 (2022).
  * Cai and Xia [2022] M. Cai and K. Xia, Phys. Rev. A 106, 042434 (2022).
  * Yalçınkaya  _et al._ [2019] İ. Yalçınkaya, B. Çakmak, G. Karpat, and F. F. Fanchini, Quantum Inf. Process. 18, 156 (2019).
  * Cohen  _et al._ [2017] I. Cohen, N. Aharon, and A. Retzker, Fortschritte Phys. 65, 1600071 (2017).
  * Mikelsons  _et al._ [2015] G. Mikelsons, I. Cohen, A. Retzker, and M. B. Plenio, New J. Phys. 17, 053032 (2015).
  * Aharon  _et al._ [2013] N. Aharon, M. Drewsen, and A. Retzker, Phys. Rev. Lett. 111, 230507 (2013).
  * Bermudez  _et al._ [2012] A. Bermudez, P. O. Schmidt, M. B. Plenio, and A. Retzker, Phys. Rev. A 85, 040302 (2012).
  * Cai  _et al._ [2012] J.-M. Cai, B. Naydenov, R. Pfeiffer, L. P. McGuinness, K. D. Jahnke, F. Jelezko, M. B. Plenio, and A. Retzker, New J. Phys. 14, 113023 (2012).
  * Chaudhry and Gong [2012] A. Z. Chaudhry and J. Gong, Phys. Rev. A 85, 012315 (2012).
  * Facchi  _et al._ [2004] P. Facchi, D. A. Lidar, and S. Pascazio, Phys. Rev. A 69, 032314 (2004).
  * Viola and Lloyd [1998] L. Viola and S. Lloyd, Phys. Rev. A 58, 2733 (1998).
  * Morong  _et al._ [2023] W. Morong, K. Collins, A. De, E. Stavropoulos, T. You, and C. Monroe, PRX Quantum 4, 010334 (2023).
  * Barthel  _et al._ [2023] P. Barthel, P. H. Huber, J. Casanova, I. Arrazola, D. Niroomand, T. Sriarunothai, M. B. Plenio, and C. Wunderlich, New J. Phys. 25, 063023 (2023).
  * Leu  _et al._ [2023] A. D. Leu, M. F. Gely, M. A. Weber, M. C. Smith, D. P. Nadlinger, and D. M. Lucas, Phys. Rev. Lett. 131, 120601 (2023).
  * Fang  _et al._ [2022] C. Fang, Y. Wang, S. Huang, K. R. Brown, and J. Kim, Phys. Rev. Lett. 129, 240504 (2022).
  * Pokharel  _et al._ [2018] B. Pokharel, N. Anand, B. Fortman, and D. A. Lidar, Phys. Rev. Lett. 121, 220502 (2018).
  * Stark  _et al._ [2018] A. Stark, N. Aharon, A. Huck, H. A. R. El-Ella, A. Retzker, F. Jelezko, and U. L. Andersen, Sci. Rep. 8, 14807 (2018).
  * Randall  _et al._ [2015] J. Randall, S. Weidt, E. D. Standing, K. Lake, S. C. Webster, D. F. Murgia, T. Navickas, K. Roth, and W. K. Hensinger, Phys. Rev. A 91, 012322 (2015).
  * Tan  _et al._ [2013] T. R. Tan, J. P. Gaebler, R. Bowler, Y. Lin, J. D. Jost, D. Leibfried, and D. J. Wineland, Phys. Rev. Lett. 110, 263002 (2013).
  * Webster  _et al._ [2013] S. C. Webster, S. Weidt, K. Lake, J. J. McLoughlin, and W. K. Hensinger, Phys. Rev. Lett. 111, 140501 (2013).
  * Souza  _et al._ [2012] A. M. Souza, G. A. Álvarez, and D. Suter, Phys. Rev. A 86, 050301 (2012).
  * Noguchi  _et al._ [2012] A. Noguchi, S. Haze, K. Toyoda, and S. Urabe, Phys. Rev. Lett. 108, 060503 (2012).
  * Xu  _et al._ [2012] X. Xu, Z. Wang, C. Duan, P. Huang, P. Wang, Y. Wang, N. Xu, X. Kong, F. Shi, X. Rong, and J. Du, Phys. Rev. Lett. 109, 070502 (2012).
  * Zalivako  _et al._ [2023] I. V. Zalivako, A. S. Borisenko, I. A. Semerikov, A. Korolkov, P. L. Sidorov, K. Galstyan, N. V. Semenin, V. Smirnov, M. A. Aksenov, A. K. Fedorov, K. Y. Khabarova, and N. N. Kolachevsky, Frontiers in Quantum Science and Technology 2, 1228208 (2023), arxiv:2305.06071 [quant-ph] .
  * Wang  _et al._ [2021] P. Wang, C.-Y. Luan, M. Qiao, M. Um, J. Zhang, Y. Wang, X. Yuan, M. Gu, J. Zhang, and K. Kim, Nat. Commun. 12, 233 (2021).
  * Sinuco-Leon  _et al._ [2021] G. A. Sinuco-Leon, H. Mas, S. Pandey, G. Vasilakis, B. M. Garraway, and W. von Klitzing, Phys. Rev. A 104, 033307 (2021).
  * Sinuco-Leon  _et al._ [2019] G. A. Sinuco-Leon, B. M. Garraway, H. Mas, S. Pandey, G. Vasilakis, V. Bolpasi, W. von Klitzing, B. Foxon, S. Jammi, K. Poulios, and T. Fernholz, Phys. Rev. A 100, 053416 (2019).
  * Trypogeorgos  _et al._ [2018] D. Trypogeorgos, A. Valdés-Curiel, N. Lundblad, and I. B. Spielman, Phys. Rev. A 97, 013407 (2018).
  * Laucht  _et al._ [2017] A. Laucht, R. Kalra, S. Simmons, J. P. Dehollain, J. T. Muhonen, F. A. Mohiyaddin, S. Freer, F. E. Hudson, K. M. Itoh, D. N. Jamieson, J. C. McCallum, A. S. Dzurak, and A. Morello, Nat. Nanotechnol. 12, 61 (2017).
  * Zhang  _et al._ [2014] J. Zhang, A. M. Souza, F. D. Brandao, and D. Suter, Phys. Rev. Lett. 112, 050502 (2014).
  * Golter  _et al._ [2014] D. A. Golter, T. K. Baldwin, and H. Wang, Phys. Rev. Lett. 113, 237601 (2014).
  * Liu  _et al._ [2013] G.-Q. Liu, H. C. Po, J. Du, R.-B. Liu, and X.-Y. Pan, Nat. Commun. 4, 2254 (2013).
  * Biercuk  _et al._ [2009] M. J. Biercuk, H. Uys, A. P. VanDevender, N. Shiga, W. M. Itano, and J. J. Bollinger, Nature 458, 996 (2009).
  * Timoney  _et al._ [2011] N. Timoney, I. Baumgart, M. Johanning, A. F. Varon, M. B. Plenio, A. Retzker, and Ch. Wunderlich, Nature 476, 185 (2011).
  * de Lange  _et al._ [2010] G. de Lange, Z. H. Wang, D. Ristè, V. V. Dobrovitski, and R. Hanson, Science 330, 60 (2010).
  * Zhang  _et al._ [2021] M. Zhang, Y. Xie, J. Zhang, W. Wang, C. Wu, T. Chen, W. Wu, and P. Chen, Phys. Rev. Applied 15, 014033 (2021).
  * Baumgart  _et al._ [2016] I. Baumgart, J.-M. Cai, A. Retzker, M. B. Plenio, and Ch. Wunderlich, Phys. Rev. Lett. 116, 240801 (2016).
  * Kotler  _et al._ [2011] S. Kotler, N. Akerman, Y. Glickman, A. Keselman, and R. Ozeri, Nature 473, 61 (2011).
  * Hall  _et al._ [2010] L. T. Hall, C. D. Hill, J. H. Cole, and L. C. L. Hollenberg, Phys. Rev. B 82, 045208 (2010).
  * Pham  _et al._ [2012] L. M. Pham, N. Bar-Gill, C. Belthangady, D. Le Sage, P. Cappellaro, M. D. Lukin, A. Yacoby, and R. L. Walsworth, Phys. Rev. B 86, 045214 (2012).
  * Bishof  _et al._ [2013] M. Bishof, X. Zhang, M. J. Martin, and J. Ye, Phys. Rev. Lett. 111, 093604 (2013).
  * Lange  _et al._ [2020] R. Lange, N. Huntemann, C. Sanner, H. Shao, B. Lipphardt, Chr. Tamm, and E. Peik, Phys. Rev. Lett. 125, 143201 (2020).
  * Shaniv  _et al._ [2019] R. Shaniv, N. Akerman, T. Manovitz, Y. Shapira, and R. Ozeri, Phys. Rev. Lett. 122, 223204 (2019).
  * Kaewuam  _et al._ [2020] R. Kaewuam, T. R. Tan, K. J. Arnold, S. R. Chanu, Z. Zhang, and M. D. Barrett, Phys. Rev. Lett. 124, 083202 (2020).
  * Yeh  _et al._ [2023] C.-H. Yeh, K. C. Grensemann, L. S. Dreissen, H. A. Fürst, and T. E. Mehlstäubler, New J. Phys. 25, 093054 (2023).
  * Dreissen  _et al._ [2022] L. S. Dreissen, C.-H. Yeh, H. A. Fürst, K. C. Grensemann, and T. E. Mehlstäubler, Nat. Commun. 13, 7314 (2022).
  * Martínez-Lahuerta  _et al._ [2024] V. J. Martínez-Lahuerta, L. Pelzer, K. Dietze, L. Krinner, P. O. Schmidt, and K. Hammerer, Quantum Sci. Technol. 9, 015013 (2024).
  * Pelzer [2023] L. Pelzer,  _Robust Artificial Clock Transition by Continuous Dynamical Decoupling in Multi-Ion Calcium Crystals_ , Ph.D. thesis, Leibniz Unversität Hannover (2023).
  * Hannig  _et al._ [2019] S. Hannig, L. Pelzer, N. Scharnhorst, J. Kramer, M. Stepanova, Z. T. Xu, N. Spethmann, I. D. Leroux, T. E. Mehlstäubler, and P. O. Schmidt, Rev. Sci. Instrum. 90, 053204 (2019).
  * Pyka  _et al._ [2014] K. Pyka, N. Herschbach, J. Keller, and T. E. Mehlstäubler, Appl. Phys. B 114, 231 (2014).
  * Beev  _et al._ [2017] N. Beev, J.-A. Fenske, S. Hannig, and P. O. Schmidt, Rev. Sci. Instrum. 88, 054704 (2017).
  * Stenger  _et al._ [2002] J. Stenger, H. Schnatz, C. Tamm, and H. Telle, Phys. Rev. Lett. 88, 073601 (2002).
  * Scharnhorst  _et al._ [2015] N. Scharnhorst, J. B. Wübbena, S. Hannig, K. Jakobsen, J. Kramer, I. D. Leroux, and P. O. Schmidt, Opt. Express 23, 19771 (2015).
  * Matei  _et al._ [2017] D. G. Matei, T. Legero, S. Häfner, C. Grebing, R. Weyrich, W. Zhang, L. Sonderhouse, J. M. Robinson, J. Ye, F. Riehle, and U. Sterr, Phys. Rev. Lett. 118, 263202 (2017).
  * Kreuter  _et al._ [2005] A. Kreuter, C. Becher, G. P. T. Lancaster, A. B. Mundt, C. Russo, H. Häffner, C. Roos, W. Hänsel, F. Schmidt-Kaler, R. Blatt, and M. S. Safronova, Phys. Rev. A 71, 032504 (2005).
  * Monz [2011] T. Monz,  _Quantum Information Processing beyond Ten Ion-Qubits_ , Ph.D. thesis, University of Innsbruck (2011).
  * Rosenband  _et al._ [2007] T. Rosenband, P. O. Schmidt, D. B. Hume, W. M. Itano, T. M. Fortier, J. E. Stalnaker, K. Kim, S. A. Diddams, J. C. J. Koelemeij, J. C. Bergquist, and D. J. Wineland, Phys. Rev. Lett. 98, 220801 (2007).
  * Peik  _et al._ [1994] E. Peik, G. Hollemann, and H. Walther, Phys. Rev. A 49, 402 (1994).
  * Huang  _et al._ [2021] Y. Huang, B. Zhang, M. Zeng, H. Zhang, Y. Hao, Z. Chen, M. Wang, H. Guan, and K. Gao, preprint 10.21203/rs.3.rs-120082/v1 (2021).
  * Hume and Leibrandt [2016] D. B. Hume and D. R. Leibrandt, Phys. Rev. A 93, 032138 (2016).
  * Dörscher  _et al._ [2020] S. Dörscher, A. Al-Masoudi, M. Bober, R. Schwarz, R. Hobson, U. Sterr, and C. Lisdat, Commun. Phys. 3, 1 (2020).
  * Rosenband and Leibrandt [2013] T. Rosenband and D. R. Leibrandt, arXiv:1303.6357 (2013), arxiv:1303.6357 .
  * Borregaard and Sørensen [2013] J. Borregaard and A. S. Sørensen, Phys. Rev. Lett. 111, 090802 (2013).
  * Keller  _et al._ [2003] M. Keller, B. Lange, K. Hayasaka, W. Lange, and H. Walther, Appl. Phys. B: Lasers Opt. 76, 125 (2003).

Supplemental Material for: “Multi-ion frequency reference using dynamical
decoupling”

##  II RF-Coils

The design and placement of the rf-magnetic field coils for continuous
dynamical decoupling (CDD) followed considerations as summarized in the
following. The segmented Paul trap used for the experiments was designed to
obtain low micromotion along the axial direction [99, 84]. We therefore
refrained from placing an rf-coil or wire close to the trapping zone, to avoid
changing the trapping potential. Instead, the signal coils were placed outside
the vacuum vessel inside an inverted viewport in a distance of approximately
$5\text{\,}\mathrm{cm}$ to the trap center (compare Fig. S1a). A resonant LCR
circuit for signal amplification is necessary to achieve the required field
strength. We choose quality factors of $Q_{S}=8.59(3)$ and $Q_{D}=15.95(4)$.
The signal distortion by the transfer function of the coils was taken into
account in the applied waveform. For this, the transfer function $H(s)$ was
measured with a receiver coil and modeled using

$H(s)=g\frac{(s-z_{r}+iz_{i})(s-z_{r}-iz_{i})}{(s-p_{r}+ip_{i})(s-p_{r}-ip_{i})}\,,$
(S1)

where the gain $g$, and a second order zero with real and imaginary values
$z_{r},z_{i}$ and a second order pole $p_{r},p_{i}$ are chosen to fit the
measured coil response. The inverse of the transfer function was applied to
the desired waveforms (Eqs. S2 \- S5) to pre-compensate the signal distortion.
A two-channel arbitrary waveform generator44433622A, Keysight (AWG) is used to
generate the calculated waveform. An rf-switch, controlled by the experimental
control system, ensures strict synchronization of the rf-pulses within the
sequence. The signal is further amplified using an rf-amplifier555ZHL-3A,
Mini-Circuits (compare Fig. S1b).

Figure S1: a) Cross-section of the ion trap setup. The alignment of the rf-
coils to the ion trap inside the vacuum vessel is displayed. b) Electric
schematic of the rf-coil signal path. The alignment of the incoupling and
signal coils is optimized for impedance matching between the
$50\text{\,}\mathrm{\SIUnitSymbolOhm}$ input and the low-impedance signal
coil.

With a typical input power ranging from
$20\text{\,}\mathrm{dB}25\text{\,}\mathrm{dB}$, the coils heat up by up to
$10\text{\,}\mathrm{K}$, depending on the duty cycle of the rf-pulse. The
temperature dependence of the resonance frequency
$f_{0}(T)=\frac{1}{\sqrt{L(T)C(T)}}$ of the circuit parts causes a variation
of the coupling strength with altering duty cycle and changes in laboratory
temperature. Therefore, the circuit design includes a temperature-controlled
baseplate for the capacitive and resistive parts. The inductive part of the
circuits consists of a copper coil held by a 3D printed polylactide (PLA)
mount. For minimal heat build-up, the mount is constructed to provide
efficient air cooling by convection.

##  III Waveforms

In the following, the employed waveforms for the CDD sequence are displayed.
The required signal to obtain an interaction of the form given in Eq. 2 is
displayed in Eq. S5. If only the first stage CDD is used it simplifies to S3.
Here, the Zeeman manifold index has been omitted for clarity. Eqs. S2 and S4
are employed as preparation steps for first and second stage, respectively.
Furthermore, the adiabatic transfer of the population is used for preparation
of the desired dressed state. The waveforms are combined continuously by
adapting the phases $\phi_{i}$ according to the preceding waveform. Note, that
the signal amplitudes $A_{1},A_{2}$ must be calibrated in situ to match the
desired coupling amplitudes $\Omega_{1},\Omega_{2}$. The frequencies
$\omega_{i}$, sweep durations $t_{\text{sw},i},\sigma$ and phases $\phi_{i}$
can be chosen with high precision. The frequency range of the frequency sweeps
$\Delta\omega_{\text{sw},i}=\omega_{i}-\omega_{i,\text{init}}$ must be small
compared to the energetic gap of the levels to avoid unwanted excitation. In
addition, amplitude and frequency changes
$\delta\Delta\omega_{\text{sw},i}/\delta t$ must be small enough to ensure
adiabaticity. The preparation durations $t_{\text{sw},i}$ have been optimized
in situ for maximal peak excitation of the target state. The sign of the
frequency sweep range $\Delta\omega_{\text{sw},i}$ determines which dressed
ground state is populated. Using the parameters given in Table S1, a target
state population $>97\%$ is obtained.

$\displaystyle y(t)=$
$\displaystyle\frac{1}{2}e^{-\frac{(t-t_{\text{sw},1})^{2}}{\sigma^{2}}}A_{1}\sin\left(\left[\omega_{1,\text{init}}t+\frac{(\omega_{1}-\omega_{1,\text{init}})}{2t_{\text{sw},1}}t^{2}\right]\right)$
(S2) $\displaystyle y(t)=$
$\displaystyle\frac{1}{2}A_{1}\sin(t\omega_{1}+\phi_{1})$ (S3) $\displaystyle
y(t)=$
$\displaystyle\frac{1}{2}A_{1}\sin(\omega_{1}t+\phi_{1})+\frac{1}{2}e^{-\frac{(t-t_{\text{sw},2})^{2}}{\sigma^{2}}}A_{2}\sin(\omega_{1}t+\frac{\pi}{2}+\phi_{1})\times\sin\left(\left[\omega_{2,\text{init}}t+\frac{\omega_{2}-\omega_{2,\text{init}}}{2t_{\text{sw},2}}t^{2}\right]+\phi_{2}\right)$
(S4) $\displaystyle y(t)=$
$\displaystyle\frac{1}{2}A_{1}\sin(\omega_{1}t+\phi_{\text{fast}})+\frac{1}{2}A_{2}\sin(\omega_{1}t+\frac{\pi}{2}+\phi_{\text{fast}})\times\sin(\omega_{2}t+\phi_{\text{slow}})$
(S5)

##  IV Continuous dynamical decoupling parameter

Parameter Set | Resonant | Magic mixing angle
---|---|---
$\omega_{S1}$ | 10002089 | 10002089
$\omega_{D1}$ | 5994834 | 5994834
$\Omega_{S1}$ | 46862 | 46862
$\Omega_{D1}$ | 115446 | 115446
$\Omega_{S2}$ | 3469 | 3469
$\Omega_{D2}$ | 6809 | 6809
$\omega_{S2}$ | 46915 | 46951
$\omega_{D2}$ | 69287 | 70731
$\Delta\omega_{\text{sw1}}$ | -150000 | -1500000
$t_{\text{sw},1}$ | 500 | 500
$\Delta\omega_{\text{sw2}}$ | 80000 | 30000
$t_{\text{sw},2}$ | 7000 | 7000
Table S1: Parameter sets for CDD waveforms. All units of frequencies
$\omega_{i}$, coupling strength $\Omega_{i}$ and sweep range
$\Delta\omega_{\text{sw},i}$ are in , all units for sweep times
$t_{\text{sw},i}$ are given in $\unit{\micro}$.


# Ultraheavy Ultrahigh-Energy Cosmic Rays

B. Theodore Zhang1, Kohta Murase2,3,4,1, Nick Ekanger5, Mukul Bhattacharya2,3,
Shunsaku Horiuchi5,6 1Center for Gravitational Physics and Quantum
Information, Yukawa Institute for Theoretical Physics, Kyoto University,
Kyoto, Kyoto 606-8502, Japan  2Department of Physics, The Pennsylvania State
University, University Park, PA 16802, USA  3Department of Astronomy &
Astrophysics, The Pennsylvania State University, University Park, PA 16802,
USA  4Center for Multimessenger Astrophysics, Institute for Gravitation and
the Cosmos, The Pennsylvania State University, University Park, PA 16802, USA
5Center for Neutrino Physics, Department of Physics, Virginia Tech,
Blacksburg, VA 24061, USA  6Kavli IPMU (WPI), UTIAS, The University of Tokyo,
Kashiwa, Chiba 277-8583, Japan

###### Abstract

We investigate the propagation of ultraheavy (UH) nuclei as ultrahigh-energy
cosmic rays (UHECRs). We show that their energy loss lengths at $\lesssim 300$
EeV are significantly longer than those of protons and intermediate nuclei,
and that the highest-energy cosmic rays with energies beyond $\sim 100$ EeV,
including the Amaterasu particle, may originate from such UH-UHECRs. We derive
constraints on the contribution of UH-UHECR sources, and find that they are
consistent with energy generation rate densities of UHECRs from collapsars and
neutron star mergers.

††preprint: APS/123-QED

##  I Introduction

The origin of ultrahigh-energy cosmic rays (UHECRs) has been a long-standing
mystery for over 50 years since the first detection of $\sim 100\rm~{}EeV$
cosmic rays [1, 2, 3]. The observed spectrum of UHECRs shows a hardening at
the ankle around $4\rm~{}EeV$ and a cutoff at $\sim 50\rm~{}EeV$ [4, 5, 6].
There are possible discrepancies in the measured energy spectra [7] between
the Pierre Auger Observatory (Auger) [8] and Telescope Array (TA) [9, 10]
observations. Even though the differences at low energies may be explained by
the systematic effect on the energy scale, an excess remains in the TA data at
the highest energies [11]. Recently, TA reported the detection of an extremely
energetic UHECR event with an energy of $244\pm 29{}_{-76}^{+51}\rm~{}EeV$,
dubbed the “Amaterasu” particle [12].

The composition of UHECRs is important to unveiling their origin (see e.g.,
Ref. [13]). The depth of the cosmic-ray shower maximum, $X_{\rm max}$, is a
measurable quantity to infer the particle composition [14]. The Auger data
favor a mixed composition of UHECRs, and intermediate-mass (like carbon and
oxygen) and/or heavy (like iron) nuclei make a dominant contribution beyond
$10\rm~{}EeV$ [15, 16, 17]. In particular, the fraction of protons gradually
decreases above the ankle, while intermediate-mass nuclei may become dominant
at higher energies and the contribution of heavy nuclei seems negligible
within the energy range of $10^{18.4}-10^{19.4}\rm~{}eV$ [17]. However, these
results are largely affected by hadronic interaction models. The distribution
of the depth of the cosmic-ray shower maximum, measured by the TA is
consistent with the Auger data, but the interpretation is still under debate
[18].

Not all UHECR source candidates can generate dominantly heavy nuclei UHECRs.
Acceleration of heavy nuclei has been considered in the context of collapsars
[19, 20, 21, 22, 23, 24, 25, 26, 27] and compact binary mergers involving a
neutron star [28, 29, 30, 31]. The energy generation rate density of UHECRs at
$10^{19.5}$ eV is $EQ_{E}^{19.5}\approx(0.2-2)\times
10^{43}\rm~{}erg~{}Mpc^{-3}~{}yr^{-1}$ (e.g., Refs. [32, 33], which is
consistent with energy budgets of collapsars including gamma-ray bursts (GRBs)
and hypernovae, as well as BNS mergers (see Table II of Ref. [34]). With the
luminosity requirement obtained from the Hillas condition [Blandford:1999hi],
the possible maximum energy of the accelerated cosmic rays is estimated to be

$\displaystyle E_{A,\rm max}$ $\displaystyle\approx 9.8\times{10}^{20}~{}{\rm
eV}~{}\left(\frac{Z}{40}\right){\left(\frac{\epsilon_{B}}{0.01}\right)}^{1/2}$
$\displaystyle\times{\left(\frac{L}{{10}^{51}~{}{\rm
erg/s}}\right)}^{1/2}{\left(\frac{\Gamma}{{10}^{2.5}}\right)}^{-1}\beta^{1/2}$
(1)

where $Z$ is the nuclei charge, $\epsilon_{B}$ is the energy fraction by the
magnetic field, $L$ is the total luminosity, $\Gamma$ is the Lorentz factor,
and $\beta$ is the characteristic velocity normalized by the speed of light
[1, 35].

Ultraheavy (UH) nuclei, which are defined as nuclei heavier than iron-group
nuclei throughout this work, are believed to be synthesized due to the
$r$-process inside neutron-rich environments [36, 37, 38, 39, 40, 41]. The
sources of UH nuclei can be compact binary mergers including binary neutron
star (BNS) and neutron-star–black-hole (NSBH) mergers (e.g., Refs. [42, 43,
44, 45, 46]), as well as collapsars including GRBs and magnetorotational
supernovae (e.g., Refs. [22, 47, 48, 49, 50, 27, 51, 52, 53]). UH nuclei have
an advantage of being accelerated to energies beyond $100\rm~{}EeV$, compared
to conventional light- and intermediate-mass group nuclei, which could provide
an additional contribution to the highest-energy cosmic rays.

In this work, we study the fate of UH-UHECRs during their intergalactic
propagation and constrain their contributions to the highest-energy cosmic
rays beyond $10^{20}$ eV. We also discuss the Amaterasu particle as a UH-UHECR
event and its implications for future observations.

##  II Propagation of UH-UHECRs

UH nuclei can be photodisintegrated and spalled into lighter nuclei due to
their interaction with background target photons and matter, respectively. For
the intergalactic propagation of UHECRs from the source to Earth, interactions
with the cosmic microwave background (CMB) and extragalactic background light
(EBL) are dominant, and we utilize the public code CRPropa 3.2 to propagate
UHECRs through intergalactic space [54]. However, because CRPropa 3.2 does not
provide a module for nuclei with mass numbers of $A>56$, we newly generate the
photodisintegration cross section table of UH-UHECRs using the nuclear
reaction network Talys 1.96 [55, 56] as an extension to CRPropa 3.2. In our
photodisintegration reaction network, the maximum atomic number and mass
number of nuclei are $Z=92$ and $A=146$, respectively, with a total of 2434
isotopes. In addition, the decay of unstable UH-UHECRs is implemented with the
data table taken from NuDat 3 [57].

Figure 1: Total energy loss lengths for various nuclei: p, O, Fe, and Pt. The
black dashed line is the energy loss length due to the adiabatic expansion of
the universe. CMB and EBL [58] are considered as target photons.

The dominant energy loss processes during the intergalactic propagation of
UHECR nuclei are photodisintegration, photomeson production, Bethe-Heitler
pair production, and adiabatic losses due to the expansion of the universe. In
Fig. 1, we show total energy loss lengths for light nuclei represented by
proton (p) and helium (He), oxygen (O) as intermediate-mass nuclei, iron (Fe)
as heavy nuclei, and platinum (Pt) as UH nuclei (see SM for other UH nuclei).
The inverse of the total energy loss length is written as $\lambda_{\rm
loss}^{-1}=\lambda_{\rm loss,phdis}^{-1}+\lambda_{\rm
loss,phmes}^{-1}+\lambda_{\rm loss,BH}^{-1}+\lambda_{\rm loss,ad}^{-1}$, where
$\lambda_{\rm loss,ad}=c/H_{0}\sim 4000\rm~{}Mpc$ is the adiabatic energy loss
length. The photodisintegration energy loss length is $\lambda_{\rm loss}^{\rm
phdis}\approx\left(n_{\rm CMB}\hat{\sigma}_{\rm phdis}\right)^{-1}\simeq
1.3(A/195)^{-0.21}\rm~{}Mpc$, where $\hat{\sigma}_{\rm
dis}=\sigma_{A\gamma}\kappa_{A\gamma}$ is the effective photodisintegration
cross section and $\kappa_{A\gamma}=1/A$ is the inelasticity at the giant
dipole resonance (GDR). The effective photodisintegration cross section can be
analytically approximated as $\sigma_{A\gamma}\approx\sigma_{\rm
GDR}\Delta\bar{\varepsilon}_{\rm GDR}/\bar{\varepsilon}_{\rm GDR}\approx
120(A/195)^{1.21}\,{\rm mb}$ for $A\gtrsim 10$, where $\sigma_{\rm GDR}\approx
4.3\times 10^{-28}A^{1.35}\rm\ {cm}^{2}$ is the GDR cross section,
$\bar{\varepsilon}_{\rm GDR}\approx 42.65A^{-0.21}\rm~{}MeV$ is the resonance
energy in the nuclear rest frame, and $\Delta\bar{\varepsilon}_{\rm
GDR}\approx 21.05A^{-0.35}\rm~{}MeV$ is the width [59]. The typical resonance
energy is $E_{A}^{\rm phdis}\approx 0.5Am_{p}c^{2}\bar{\varepsilon}_{\rm
GDR}/\varepsilon_{t}\simeq 1.9\times 10^{21}~{}{\rm
eV}~{}(A/195)^{0.79}(\varepsilon_{t}/6.6\times 10^{-4}~{}{\rm eV})^{-1}$,
where $\varepsilon_{t}$ is the target photon energy. Energy losses due to the
Bethe-Heitler pair production process are also important for UH-UHECRs because
of their large atomic numbers. The energy loss length is estimated to be
$\lambda_{\rm loss}^{\rm BH}\approx\left(n_{\rm CMB}\hat{\sigma}_{\rm
BH}\right)^{-1}\simeq 32\ (Z/78)^{-2}(A/195)\rm~{}Mpc$, where
$\hat{\sigma}_{\rm BH}\sim 8\times 10^{-31}(Z^{2}/A)\rm~{}cm^{2}$ is the
effective cross section [34], $E_{A}^{\rm BH}\approx
0.5Am_{p}c^{2}\bar{\varepsilon}_{\rm BH}/\varepsilon_{\rm CMB}\simeq 1.4\times
10^{21}~{}{\rm eV}~{}(A/195)(\varepsilon_{t}/6.6\times 10^{-4}\rm~{}eV)^{-1}$
is the typical energy, and $\bar{\varepsilon}_{\rm BH}\sim 10\rm~{}MeV$.
Although the photodisintegration process is expected to be dominant, which is
indeed the case at $E_{A}\sim{\rm a~{}few}\times 10^{21}$ eV, we find that the
Bethe-Heitler pair production process is the most important at
$E_{A}\sim(1-5)\times 10^{20}\rm~{}eV$ [see Fig. 6 in Supplementary Material
(SM)]. The photomeson production process is irrelevant when considering the
propagation of UH-UHECRs due to the high-energy threshold at $\sim 4\times
10^{22}\rm~{}(A/195)(\varepsilon_{t}/6.6\times 10^{-4}\rm~{}eV)^{-1}\rm~{}eV$.

As seen in Fig 1, UH-UHECRs can travel longer distances than the GZK (Greisen-
Zatsepin-Kuzmin) distance for UHECR protons [60, 61] and the energy loss
length of iron-group-mass nuclei [62, 63]. Thanks to their large mass numbers,
it takes longer for them to be completely disintegrated, by which rarer
sources can contribute or UH-UHECRs could appear beyond the cutoff in spectra
of UHECR protons and conventional nuclei. To calculate the UHECR spectrum and
composition on Earth, we take into account the change in the composition of
UH-UHECRs during propagation, as the UH nuclei become progressively lighter
with the ejection of one or more nucleons. See the Supplementary Material for
details.

##  III UH-UHECR Spectrum and Composition

Figure 2: Energy spectrum and the first/second moments of $X_{\rm max}$
distribution considering both conventional and UH nuclei. The energy spectrum
of UHECRs is obtained from Ref. [64], and the first/second moments of $X_{\rm
max}$ distribution is taken from Ref. [65].

Using the available Auger and TA data on the energy spectrum and composition
of UHECRs, we examine possible contributions of the sources of UH-UHECRs. In
this work, we consider two representative composition models of conventional
nuclei, depending on whether there are heavy iron-group nuclei or not:

  * •

Model A: conventional nuclei (p, He, O, Si) + UH nuclei (Se, Te, Pt)

  * •

Model B: conventional nuclei (p, He, O, Si, Fe) + UH nuclei (Se, Te, Pt)

The origin of intermediate-mass nuclei can be associated with material from
the inner core of massive stars collapsing into black holes [66, 25]. Iron-
group nuclei are naturally expected in the supernova/hypernova ejecta [19, 20,
23, 25, 27] or the surface of neutron stars [67]. Alternatively, UHECR nuclei
can be provided by the reacceleration of Galactic cosmic rays [68, 69]. We
treat UH-UHECRs as an additional component, assuming that they are produced by
either the same or different populations. We consider Se, Te, and Pt to
represent the first, second, and third peaks of nuclei synthesized via the
$r$-process.

Figure 3: Same as Fig. 2, but for the TA data. Note that the data in the red
vertical band corresponds to the Amaterasu particle [12].

The spectrum of UHECRs escaping from the sources is assumed to have a power-
law spectrum with an exponential cutoff,

$\frac{d\dot{N_{A}}}{d\mathcal{R}}=f_{A}\dot{N}_{0}\left(\frac{\mathcal{R}}{\mathcal{R}_{0}}\right)^{-s_{\rm
CR}}{\rm exp}\left(-\frac{\mathcal{R}}{\mathcal{R}_{\rm max}}\right),$ (2)

where $\mathcal{R}=E/Z$ is particle rigidity, $E$ is particle energy,
$\mathcal{R}_{0}$ is a reference rigidity, $f_{A}$ is the number fraction of
nuclei, $s_{\rm CR}$ is the spectral index, and $\dot{N}_{0}$ is a
normalization constant. To reduce the number of free parameters, we assume all
of the accelerated nuclei have the same rigidity $\mathcal{R}_{\rm max}$ and
spectral index $s_{\rm CR}$. Large-scale magnetic fields in structured regions
of the Universe could further modify the observed spectrum and composition,
where CRs with low rigidity cannot arrive on Earth, known as the magnetic
horizon effect [70, 71]. However, this effect is important if the average
source separation is $d_{s}\gtrsim 20\sqrt{l_{c}/0.1\rm~{}Mpc}\rm~{}Mpc$,
where $l_{c}$ is the coherence length [72]. Because transient sources like
collapsars and BNS mergers may lead to sufficient source densities [73, 74],
it is reasonable to adopt the one-dimensional propagation scenario assuming
UHECR sources have an average distance less than $20\rm~{}Mpc$. We neglect the
effect of redshift evolution. In this work, we adopt the hadronic interaction
model Epos-LHC [75, 13].

To constrain the energy generation density of UH-UHECRs, $Q_{\rm UH-
UHECR}\equiv\int_{E_{\rm min}}^{E_{\rm
max}}d\mathcal{R}E\frac{d\dot{N_{A}}}{d\mathcal{R}}$, with $E_{\rm
min}=10^{18}\rm~{}eV$ and $E_{\rm max}=10^{21}\rm~{}eV$. We perform a combined
fit with both conventional and UH nuclei components. For simplicity, we only
consider the combination of conventional nuclei with one of the UH nuclei,
among Se, Te, and Pt, for both Model A and Model B, which is helpful to
explore the effect of the mass and charge of UH nuclei on the fitting. For the
TA data, we show the best-fit values of the energy injection rate of UH-
UHECRs, where the uncertainty is estimated under the condition that the total
chi-square value $\chi_{\rm tot}^{2}\leqslant\chi_{\rm tot,min}^{2}+1$. For
the Auger data, we use $\chi_{\rm tot}^{2}\leqslant\chi_{\rm tot,min}^{2}+10$
to estimate the uncertainty range because of much smaller statistic errors in
the measured flux and composition. Our results are summarized in Table 1.

In Fig. 2, we provide the results of our fitting for the energy spectrum and
composition of UHECRs measured by Auger [64, 65] (see SM for details). With
Model A and Model B, we see that the energy generation rate density of the
three UH nuclear species is constrained to be

$Q_{\rm UH-UHECR}^{\rm Auger}\lesssim(0.1-8)\times
10^{42}\rm~{}erg~{}Mpc^{-3}~{}yr^{-1}.$ (3)

Note that the best-fit values derived from Model B can be about one order
magnitude smaller than those from Model A, because the presence of heavy iron-
group nuclei gives tighter constraints on the fraction of UH-UHECRs. We note
that the values have significant uncertainties, and more stringent constraints
on the energy budget of UH-UHECRs are obtained when we use the best-fit models
only with conventional nuclei (see SM for details).

Nuclei | $Q_{\rm UH-UHECR}^{\rm Auger}$ | $Q_{\rm UH-UHECR}^{\rm TA}$
---|---|---
| $[\rm erg~{}Mpc^{-3}~{}yr^{-1}]$ | $[\rm erg~{}Mpc^{-3}~{}yr^{-1}]$
Model A |  |
Se | $3.5_{-2.0}^{+0.0}\times 10^{42}$ | $9.9_{-7.6}^{+10}\times 10^{42}$
Te | $3.8_{-1.3}^{+4.0}\times 10^{41}$ | $3.8_{-2.9}^{+2.9}\times 10^{42}$
Pt | $1.7_{-1.0}^{+4.1}\times 10^{42}$ | $1.9_{-1.4}^{+1.7}\times 10^{42}$
Model B |  |
Se | $2.0_{-1.6}^{+0.0}\times 10^{42}$ | $7.3_{-6.8}^{+4.2}\times 10^{42}$
Te | $0.6_{-0.0}^{+4.1}\times 10^{41}$ | $1.9_{-1.1}^{+5.0}\times 10^{42}$
Pt | $1.9_{-1.7}^{+0.0}\times 10^{41}$ | $1.5_{-1.1}^{+1.3}\times 10^{42}$
Table 1: Energy injection rate densities of UH-UHECRs, allowed by the two
composition models considered in this work.

The results for the TA data are shown in Fig. 3, where we find that the energy
generation rate densities of UH-UHECRs are constrained to be

$Q_{\rm UH-UHECR}^{\rm TA}\lesssim(0.4-20)\times
10^{42}\rm~{}erg~{}Mpc^{-3}~{}yr^{-1},$ (4)

for both composition models, which is about 3 times larger than that derived
based on the Auger data for both composition models. In particular, we find
that including UH-UHECRs as a second population can provide better fits to the
observed energy spectrum and composition measured by TA with $\chi_{\rm
min}^{2}/{\rm d.o.f.}\sim(0.2-3)$, compared to the fits only with conventional
nuclei for both composition models, which give $\chi_{\rm min}^{2}/{\rm
d.o.f.}\sim(4-6)$. The inclusion of additional UH-UHECRs does not show
improvements when fitting to Auger data with $\chi_{\rm min}^{2}/{\rm
d.o.f.}\sim(3-4)$ (cf. $\chi_{\rm min}^{2}/{\rm d.o.f.}\sim(3-3.5)$ for the
fits with only conventional nuclei.). These results may shed light on the
detected number of the highest-energy events of UHECRs between Auger and TA,
where UH-UHECRs may have a larger effect on the latter.

Recently, the TA Collaboration reported the detection of an extremely
energetic UHECR event by the surface detector, with an energy of $244\pm
29{}_{-76}^{+51}\rm~{}EeV$ [12]. The nature of the particle is unclear at
present. Even though the primary photon is excluded at the 99.985% confidence
level, it is difficult to distinguish whether it is a proton or heavy nucleus
[12]. Fig. 3 shows that the Amaterasu particle may be explained as a UH-UHECR
event. To explore this possibility, we also show the back-tracked direction of
the Amaterasu particle for different nuclear species, as shown in Fig. 4,
where we adopt the Galactic magnetic field model provided in Ref. [76]. Thanks
to larger atomic numbers, it may come from the supergalactic plane rather than
the local void region, which could also help us to explain the discrepancy
between the TA and Augur spectra.

If UH-UHECRs contribute to the observed highest-energy cosmic rays, the
required energy budget is $Q_{\rm UH-UHECR}\sim
10^{42}\rm~{}erg~{}Mpc^{-3}~{}yr^{-1}$, which can be satisfied by both
collapsars and BNS/NSBH mergers [34]. For example, the cosmic-ray luminosity
density of BNS mergers is $\sim 10^{43.5}\rm~{}erg~{}Mpc^{-3}~{}yr^{-1}$,
which is consistent with the kinetic energy $\mathcal{E}_{\rm
ej}\approx(1/2)M_{\rm ej}c^{2}\beta_{\rm ej}^{2}\sim 2\times 10^{51}~{}{\rm
erg}~{}(M_{\rm ej}/0.05~{}M_{\odot}){(\beta_{\rm ej}/0.2)}^{2}$, the rate
density $\rho\sim 300\rm~{}Gpc^{-3}~{}yr^{-1}$, and the energy fraction by
cosmic rays, $\epsilon_{\rm CR}\sim 5$%. Short GRBs have isotropic-equivalent
energies of ${\mathcal{E}}_{\gamma}^{\rm iso}\sim 10^{51}-10^{52}~{}{\rm erg}$
and rate densities of $\rho\sim 10~{}{\rm Gpc}^{-3}~{}{\rm yr}^{-1}$ [77, 78],
inferring that the gamma-ray luminosity density is
$\sim{10}^{43}-{10}^{44}~{}{\rm erg}~{}{\rm Mpc}^{-3}~{}{\rm yr}^{-1}$. Thus,
it is possible for UH nuclei from neutron star mergers and short GRBs to
contribute to UHECRs above $\sim 10^{20}$ eV.

##  IV Summary and discussion.

We numerically studied the propagation of UH-UHECRs and derived constraints on
their contribution to the observed UHECR flux. Thanks to their energy loss
lengths at $\lesssim 10^{21}$ eV, which are longer than those of protons and
intermediate-mass nuclei, UH-UHECRs could potentially contribute to the
highest-energy cosmic rays beyond $\sim 10^{20}$ eV, including the Amaterasu
event. The allowed energy generation rate densities are consistent with those
of collapsars and BNS mergers.

Establishing the highest-energy cosmic rays as UH nuclei may infer that UHECRs
are produced by transients rather than steady sources such as active galactic
nuclei. The candidate sources include BNS mergers [29, 30, 34] and collapsars
[27]. UH nuclei presumably originate from outflows and ejecta rather than the
ambient interstellar material. Thus, promising acceleration sites would be an
external reverse shock [19, 25] or internal dissipation [66, 26], rather than
external forward shock (e.g., Ref. [29]). More detailed implications for the
sources, including the nucleus survival problem, will be discussed in future
work [59].

The UH nuclear origin of UHECRs can be tested with future measurements of the
composition at the highest-energies, as shown in Figs. 2 and 3. This behavior
is contrast to a high-luminosity GRB model that predicts a lighter composition
at the highest energies [24] and a reacceleration model that changes toward an
iron-group composition [29]. Multimessenger observations with neutrinos and
gamma rays would also be useful. The UH-UHECR scenario requires the survival
of nuclei, in which TeV gamma rays are likely to escape from the sources. The
gamma-ray signal from nuclear deexcitation and electromagnetic cascades
induced by the Bethe-Heitler process could be interesting targets for future
gamma-ray detectors such as the Cherenkov Telescope Array [79].

Figure 4: The skymap of backtracking particles with mean energy
$E=243\rm~{}EeV$ and variation $E=70\rm~{}EeV$ for p ($Z=1$, black), Fe
($Z=56$, red), Zr ($Z=40$, green) and Pt ($Z=78$, blue) in equatorial
coordinates. The black circle indicates the arrival direction of the Amaterasu
particle. The supergalactic plane is shown as dotted curves, and the Galactic
plane is shown as solid curves. This figure is generated with the analysis
code provided by Ref. [12].

###### Acknowledgements.

We thank Toshihiro Fujii, Xilu Wang, and Kunihito Ioka for the useful
information and discussions. This work is supported by NSF Grants Nos.
AST-1908689 (K.M., N.E., M.B and S.H.), AST-2108466 (K.M.), AST-2108467
(K.M.), and AST-1908960 (K.M., N.E., M.B. and S.H.), PHY-2209420 (S.H.), and
U.S. DOE Office of Science award number DE-SC0020262 (S.H.), and KAKENHI Nos.
20H01901 (K.M. and B.T.Z.), 20H05852 (K.M. and B.T.Z.), 22K03630 (S.H.), and
23H04899 (S.H.). M.B. acknowledges support from the Eberly Postdoctoral
Fellowship at the Pennsylvania State University. This work was supported by
World Premier International Research Center Initiative (WPI Initiative), MEXT,
Japan. Note added. — While the paper was being finalized, the related work
[80] appeared. Our work is independent and focuses on the general aspects of
UH-UHECRs that can be produced by not only BNS mergers but also collapsars.

## References

  * Hillas [1984] A. M. Hillas, The Origin of Ultrahigh-Energy Cosmic Rays, Ann. Rev. Astron. Astrophys. 22, 425 (1984).
  * Nagano and Watson [2000] M. Nagano and A. A. Watson, Observations and implications of the ultrahigh-energy cosmic rays, Rev. Mod. Phys. 72, 689 (2000).
  * Kotera and Olinto [2011] K. Kotera and A. V. Olinto, The Astrophysics of Ultrahigh Energy Cosmic Rays, Ann. Rev. Astron. Astrophys. 49, 119 (2011), arXiv:1101.4256 [astro-ph.HE] .
  * Abbasi  _et al._ [2008] R. U. Abbasi _et al._ (HiRes Collaboration), First observation of the Greisen-Zatsepin-Kuzmin suppression, Phys. Rev. Lett. 100, 101101 (2008), arXiv:astro-ph/0703099 [astro-ph] .
  * Abraham  _et al._ [2008] J. Abraham _et al._ (Pierre Auger Collaboration), Observation of the suppression of the flux of cosmic rays above $4\times 10^{19}$eV, Phys. Rev. Lett. 101, 061101 (2008), arXiv:0806.4302 [astro-ph] .
  * Abraham  _et al._ [2010a] J. Abraham _et al._ (Pierre Auger Collaboration), Measurement of the energy spectrum of cosmic rays above $10^{18}$ eV using the Pierre Auger Observatory, Phys. Lett. B685, 239 (2010a), arXiv:1002.1975 [astro-ph.HE] .
  * Ivanov [2018] D. Ivanov (Telescope-Array, Pierre Auger), Report of the Telescope Array - Pierre Auger Observatory Working Group on Energy Spectrum, PoS ICRC2017, 498 (2018).
  * Aab  _et al._ [2015] A. Aab _et al._ (Pierre Auger), The Pierre Auger Cosmic Ray Observatory, Nucl. Instrum. Meth. A 798, 172 (2015), arXiv:1502.01323 [astro-ph.IM] .
  * Tokuno  _et al._ [2012] H. Tokuno _et al._ , New air fluorescence detectors employed in the Telescope Array experiment, Nucl. Instrum. Meth. A 676, 54 (2012), arXiv:1201.0002 [astro-ph.IM] .
  * Abu-Zayyad  _et al._ [2013] T. Abu-Zayyad _et al._ (Telescope Array), The surface detector array of the Telescope Array experiment, Nucl. Instrum. Meth. A 689, 87 (2013), arXiv:1201.4964 [astro-ph.IM] .
  * Bergman  _et al._ [2024] D. R. Bergman _et al._ (Pierre Auger, Telescope Array), Measurement of UHECR energy spectrum with the Pierre Auger Observatory and the Telescope Array, PoS ICRC2023, 406 (2024).
  * Abbasi  _et al._ [2023] R. U. Abbasi _et al._ (Telescope Array), An extremely energetic cosmic ray observed by a surface detector array, Science 382, abo5095 (2023), arXiv:2311.14231 [astro-ph.HE] .
  * Halim  _et al._ [2023] A. A. Halim _et al._ (Pierre Auger), Constraining the sources of ultra-high-energy cosmic rays across and above the ankle with the spectrum and composition data measured at the Pierre Auger Observatory, JCAP 05, 024, arXiv:2211.02857 [astro-ph.HE] .
  * Greisen [1960] K. Greisen, Cosmic ray showers, Ann. Rev. Nucl. Part. Sci. 10, 63 (1960).
  * Abraham  _et al._ [2010b] J. Abraham _et al._ (Pierre Auger), Measurement of the Depth of Maximum of Extensive Air Showers above $10^{18}$ eV, Phys. Rev. Lett. 104, 091101 (2010b), arXiv:1002.0699 [astro-ph.HE] .
  * Aab  _et al._ [2014] A. Aab _et al._ (Pierre Auger), Depth of Maximum of Air-Shower Profiles at the Pierre Auger Observatory: Measurements at Energies above $10^{17.8}$ eV, Phys. Rev. D 90, 122005 (2014), arXiv:1409.4809 [astro-ph.HE] .
  * Abdul Halim  _et al._ [2023a] A. Abdul Halim _et al._ (Pierre Auger), Studies of the mass composition of cosmic rays and proton-proton interaction cross-sections at ultra-high energies with the Pierre Auger Observatory, PoS ICRC2023, 438 (2023a).
  * Abdul Halim  _et al._ [2023b] A. Abdul Halim _et al._ (Pierre Auger), Depth of maximum of air-shower profiles: testing the compatibility of the measurements at the Pierre Auger Observatory and the Telescope Array, PoS ICRC2023, 249 (2023b).
  * Murase  _et al._ [2008] K. Murase, K. Ioka, S. Nagataki, and T. Nakamura, High-energy cosmic-ray nuclei from high- and low-luminosity gamma-ray bursts and implications for multi-messenger astronomy, Phys. Rev. D78, 023005 (2008), arXiv:0801.2861 [astro-ph] .
  * Wang  _et al._ [2008] X.-Y. Wang, S. Razzaque, and P. Meszaros, On the Origin and Survival of UHE Cosmic-Ray Nuclei in GRBs and Hypernovae, Astrophys. J. 677, 432 (2008), arXiv:0711.2065 [astro-ph] .
  * Liu  _et al._ [2011] R.-Y. Liu, X.-Y. Wang, and Z.-G. Dai, Nearby low-luminosity GRBs as the sources of ultra-high energy cosmic rays revisited, Mon. Not. Roy. Astron. Soc. 418, 1382 (2011), arXiv:1108.1551 [astro-ph.HE] .
  * Metzger  _et al._ [2011] B. D. Metzger, D. Giannios, and S. Horiuchi, Heavy Nuclei Synthesized in Gamma-Ray Burst Outflows as the Source of UHECRs, Mon. Not. Roy. Astron. Soc. 415, 2495 (2011), arXiv:1101.4019 [astro-ph.HE] .
  * Horiuchi  _et al._ [2012] S. Horiuchi, K. Murase, K. Ioka, and P. Meszaros, The survival of nuclei in jets associated with core-collapse supernovae and gamma-ray bursts, Astrophys. J. 753, 69 (2012), arXiv:1203.0296 [astro-ph.HE] .
  * Zhang  _et al._ [2018] B. T. Zhang, K. Murase, S. S. Kimura, S. Horiuchi, and P. Mészáros, Low-luminosity gamma-ray bursts as the sources of ultrahigh-energy cosmic ray nuclei, Phys. Rev. D97, 083010 (2018), arXiv:1712.09984 [astro-ph.HE] .
  * Zhang and Murase [2019] B. T. Zhang and K. Murase, Ultrahigh-energy cosmic-ray nuclei and neutrinos from engine-driven supernovae, Phys. Rev. D 100, 103004 (2019), arXiv:1812.10289 [astro-ph.HE] .
  * Boncioli  _et al._ [2019] D. Boncioli, D. Biehl, and W. Winter, On the common origin of cosmic rays across the ankle and diffuse neutrinos at the highest energies from low-luminosity Gamma-Ray Bursts, Astrophys. J. 872, 110 (2019), arXiv:1808.07481 [astro-ph.HE] .
  * Bhattacharya  _et al._ [2022] M. Bhattacharya, S. Horiuchi, and K. Murase, On the synthesis of heavy nuclei in protomagnetar outflows and implications for ultra-high energy cosmic rays, Mon. Not. Roy. Astron. Soc. 514, 6011 (2022), arXiv:2111.05863 [astro-ph.HE] .
  * Kyutoku and Ioka [2016] K. Kyutoku and K. Ioka, The unreasonable weakness of r-process cosmic rays in the neutron-star-merger nucleosynthesis scenario, Astrophys. J. 827, 83 (2016), arXiv:1603.00467 [astro-ph.HE] .
  * Kimura  _et al._ [2018a] S. S. Kimura, K. Murase, and P. Mészáros, Super-Knee Cosmic Rays from Galactic Neutron Star Merger Remnants, Astrophys. J. 866, 51 (2018a), arXiv:1807.03290 [astro-ph.HE] .
  * Rodrigues  _et al._ [2019] X. Rodrigues, D. Biehl, D. Boncioli, and A. M. Taylor, Binary neutron star merger remnants as sources of cosmic rays below the “Ankle”, Astropart. Phys. 106, 10 (2019), arXiv:1806.01624 [astro-ph.HE] .
  * Murase and Fukugita [2019] K. Murase and M. Fukugita, Energetics of High-Energy Cosmic Radiations, Phys. Rev. D 99, 063012 (2019), arXiv:1806.04194 [astro-ph.HE] .
  * Katz  _et al._ [2009] B. Katz, R. Budnik, and E. Waxman, The energy production rate & the generation spectrum of UHECRs, JCAP 0903, 020, arXiv:0811.3759 [astro-ph] .
  * Jiang  _et al._ [2021] Y. Jiang, B. T. Zhang, and K. Murase, Energetics of ultrahigh-energy cosmic-ray nuclei, Phys. Rev. D 104, 043017 (2021), arXiv:2012.03122 [astro-ph.HE] .
  * Murase  _et al._ [2018] K. Murase, F. Oikonomou, and M. Petropoulou, Blazar Flares as an Origin of High-Energy Cosmic Neutrinos?, Astrophys. J. 865, 124 (2018), arXiv:1807.04748 [astro-ph.HE] .
  * Kotera  _et al._ [2010] K. Kotera, D. Allard, and A. V. Olinto, Cosmogenic Neutrinos: parameter space and detectabilty from PeV to ZeV, JCAP 1010, 013, arXiv:1009.1382 [astro-ph.HE] .
  * Burbidge  _et al._ [1957] E. M. Burbidge, G. R. Burbidge, W. A. Fowler, and F. Hoyle, Synthesis of the elements in stars, Rev. Mod. Phys. 29, 547 (1957).
  * Thielemann  _et al._ [2017] F. K. Thielemann, M. Eichler, I. V. Panov, and B. Wehmeyer, Neutron Star Mergers and Nucleosynthesis of Heavy Elements, Ann. Rev. Nucl. Part. Sci. 67, 253 (2017), arXiv:1710.02142 [astro-ph.HE] .
  * Horowitz  _et al._ [2019] C. J. Horowitz _et al._ , $r$-Process Nucleosynthesis: Connecting Rare-Isotope Beam Facilities with the Cosmos, J. Phys. G 46, 083001 (2019), arXiv:1805.04637 [astro-ph.SR] .
  * Kajino  _et al._ [2019] T. Kajino, W. Aoki, A. B. Balantekin, R. Diehl, M. A. Famiano, and G. J. Mathews, Current status of r -process nucleosynthesis, Prog. Part. Nucl. Phys. 107, 109 (2019), arXiv:1906.05002 [astro-ph.HE] .
  * Cowan  _et al._ [2021] J. J. Cowan, C. Sneden, J. E. Lawler, A. Aprahamian, M. Wiescher, K. Langanke, G. Martínez-Pinedo, and F.-K. Thielemann, Origin of the heaviest elements: The rapid neutron-capture process, Rev. Mod. Phys. 93, 15002 (2021), arXiv:1901.01410 [astro-ph.HE] .
  * Arcones and Thielemann [2023] A. Arcones and F.-K. Thielemann, Origin of the elements, Astron. Astrophys. Rev. 31, 1 (2023).
  * Freiburghaus  _et al._ [1999] C. Freiburghaus, J.-F. Rembges, T. Rauscher, E. Kolbe, F.-K. Thielemann, K.-L. Kratz, B. Pfeiffer, and J. J. Cowan, The astrophysical r-process: A comparison of calculations following adiabatic expansion with classical calculations based on neutron densities and temperatures, The Astrophysical Journal 516, 381 (1999).
  * Wanajo  _et al._ [2014] S. Wanajo, Y. Sekiguchi, N. Nishimura, K. Kiuchi, K. Kyutoku, and M. Shibata, Production of all the $r$-process nuclides in the dynamical ejecta of neutron star mergers, Astrophys. J. Lett. 789, L39 (2014), arXiv:1402.7317 [astro-ph.SR] .
  * Bartos and Marka [2019] I. Bartos and S. Marka, A nearby neutron-star merger explains the actinide abundances in the early Solar System, Nature 569, 85 (2019).
  * Ekanger  _et al._ [2023] N. Ekanger, M. Bhattacharya, and S. Horiuchi, Nucleosynthesis in outflows of compact objects and detection prospects of associated kilonovae, Mon. Not. Roy. Astron. Soc. 525, 2040 (2023), arXiv:2303.00765 [astro-ph.HE] .
  * Fujibayashi  _et al._ [2023] S. Fujibayashi, K. Kiuchi, S. Wanajo, K. Kyutoku, Y. Sekiguchi, and M. Shibata, Comprehensive Study of Mass Ejection and Nucleosynthesis in Binary Neutron Star Mergers Leaving Short-lived Massive Neutron Stars, Astrophys. J. 942, 39 (2023), arXiv:2205.05557 [astro-ph.HE] .
  * Siegel  _et al._ [2019] D. M. Siegel, J. Barnes, and B. D. Metzger, Collapsars as a major source of r-process elements, Nature 569, 241 (2019), arXiv:1810.00098 [astro-ph.HE] .
  * Barnes and Metzger [2022] J. Barnes and B. D. Metzger, Signatures of r-process Enrichment in Supernovae from Collapsars, Astrophys. J. Lett. 939, L29 (2022), arXiv:2205.10421 [astro-ph.HE] .
  * Nishimura  _et al._ [2015] N. Nishimura, T. Takiwaki, and F. K. Thielemann, The r-process nucleosynthesis in the various jet-like explosions of magnetorotational core-collapse supernovae, Astrophys. J. 810, 109 (2015), arXiv:1501.06567 [astro-ph.SR] .
  * Yong  _et al._ [2021] D. Yong _et al._ , r-Process elements from magnetorotational hypernovae, Nature 595, 223 (2021), arXiv:2107.03010 [astro-ph.SR] .
  * Ekanger  _et al._ [2022] N. Ekanger, M. Bhattacharya, and S. Horiuchi, Systematic exploration of heavy element nucleosynthesis in protomagnetar outflows, Mon. Not. Roy. Astron. Soc. 513, 405 (2022), arXiv:2201.03576 [astro-ph.HE] .
  * Zha  _et al._ [2024] S. Zha, B. Müller, and J. Powell, Nucleosynthesis in the Innermost Ejecta of Magnetorotational Supernova Explosions in 3-dimensions (2024), arXiv:2403.02072 [astro-ph.HE] .
  * Reichert  _et al._ [2024] M. Reichert, M. Bugli, J. Guilet, M. Obergaulinger, M. A. Aloy, and A. Arcones, Nucleosynthesis in magnetorotational supernovae: impact of the magnetic field configuration, Mon. Not. Roy. Astron. Soc. 529, 3197 (2024), arXiv:2401.14402 [astro-ph.HE] .
  * Alves Batista  _et al._ [2022] R. Alves Batista _et al._ , CRPropa 3.2 — an advanced framework for high-energy particle propagation in extragalactic and galactic spaces, JCAP 09, 035, arXiv:2208.00107 [astro-ph.HE] .
  * Koning  _et al._ [2005] A. J. Koning, S. Hilaire, and M. C. Duijvestijn, TALYS: Comprehensive Nuclear Reaction Modeling, AIP Conf. Proc. 769, 1154 (2005).
  * Koning, A. J.  _et al._ [2007] Koning, A. J., Hilaire, S., and Duijvestijn, M. C., TALYS-1.0, International Conference on Nuclear Data for Science and Technology , 211 (2007).
  * [57] NuDat 3 Data Base, National Nuclear Data Center (NNDC), Brookhaven National Laboratory, Upton, NY, USA, https://www.nndc.bnl.gov/nudat3/.
  * Gilmore  _et al._ [2012] R. C. Gilmore, R. S. Somerville, J. R. Primack, and A. Dominguez, Semi-analytic modeling of the EBL and consequences for extragalactic gamma-ray spectra, Mon. Not. Roy. Astron. Soc. 422, 3189 (2012), arXiv:1104.0671 [astro-ph.CO] .
  * [59] N. Ekanger _et al._ , in preparation .
  * Greisen [1966] K. Greisen, End to the cosmic ray spectrum?, Phys. Rev. Lett. 16, 748 (1966).
  * Zatsepin and Kuzmin [1966] G. T. Zatsepin and V. A. Kuzmin, Upper limit of the spectrum of cosmic rays, JETP Lett. 4, 78 (1966), [Pisma Zh. Eksp. Teor. Fiz.4,114(1966)].
  * Stecker [1998] F. W. Stecker, On the origin of the highest energy cosmic rays, Phys. Rev. Lett. 80, 1816 (1998), arXiv:astro-ph/9710353 .
  * Stecker and Salamon [1999] F. W. Stecker and M. H. Salamon, Photodisintegration of ultrahigh-energy cosmic rays: A New determination, Astrophys. J. 512, 521 (1999), arXiv:astro-ph/9808110 .
  * Aab  _et al._ [2020] A. Aab _et al._ (Pierre Auger), Measurement of the cosmic-ray energy spectrum above $2.5{\times}10^{18}$ eV using the Pierre Auger Observatory, Phys. Rev. D 102, 062005 (2020), arXiv:2008.06486 [astro-ph.HE] .
  * Yushkov [2020] A. Yushkov (Auger), Mass Composition of Cosmic Rays with Energies above 1017.2 eV from the Hybrid Data of the Pierre Auger Observatory, PoS ICRC2019, 482 (2020).
  * Zhang  _et al._ [2017] B. T. Zhang, K. Murase, F. Oikonomou, and Z. Li, High-energy cosmic ray nuclei from tidal disruption events: Origin, survival, and implications, Phys. Rev. D96, 063007 (2017), arXiv:1706.00391 [astro-ph.HE] .
  * Fang  _et al._ [2012] K. Fang, K. Kotera, and A. V. Olinto, Newly-born pulsars as sources of ultrahigh energy cosmic rays, Astrophys. J. 750, 118 (2012), arXiv:1201.5197 [astro-ph.HE] .
  * Caprioli [2015] D. Caprioli, ”Espresso” Acceleration of Ultra-high-energy Cosmic Rays, Astrophys. J. 811, L38 (2015), arXiv:1505.06739 [astro-ph.HE] .
  * Kimura  _et al._ [2018b] S. S. Kimura, K. Murase, and B. T. Zhang, Ultrahigh-energy Cosmic-ray Nuclei from Black Hole Jets: Recycling Galactic Cosmic Rays through Shear Acceleration, Phys. Rev. D 97, 023026 (2018b), arXiv:1705.05027 [astro-ph.HE] .
  * Lemoine [2005] M. Lemoine, Extra-galactic magnetic fields and the second knee in the cosmic-ray spectrum, Phys. Rev. D 71, 083007 (2005), arXiv:astro-ph/0411173 .
  * Berezinsky and Gazizov [2007] V. Berezinsky and A. Z. Gazizov, Diffusion of Cosmic Rays in the Expanding Universe. 2. Energy Spectra of Ultra-High Energy Cosmic Rays, Astrophys. J. 669, 684 (2007), arXiv:astro-ph/0702102 .
  * Halim  _et al._ [2024] A. A. Halim _et al._ (Pierre Auger), Impact of the Magnetic Horizon on the Interpretation of the Pierre Auger Observatory Spectrum and Composition Data (2024), arXiv:2404.03533 [astro-ph.HE] .
  * Murase and Takami [2009] K. Murase and H. Takami, Implications of Ultra-High-Energy Cosmic Rays for Transient Sources in the Auger Era, Astrophys. J. 690, L14 (2009), arXiv:0810.1813 [astro-ph] .
  * Takami and Murase [2012] H. Takami and K. Murase, The Role of Structured Magnetic Fields on Constraining Properties of Transient Sources of Ultra-high-energy Cosmic Rays, Astrophys. J. 748, 9 (2012), arXiv:1110.3245 [astro-ph.HE] .
  * Pierog  _et al._ [2015] T. Pierog, I. Karpenko, J. M. Katzy, E. Yatsenko, and K. Werner, EPOS LHC: Test of collective hadronization with data measured at the CERN Large Hadron Collider, Phys. Rev. C 92, 034906 (2015), arXiv:1306.0121 [hep-ph] .
  * Jansson and Farrar [2012] R. Jansson and G. R. Farrar, The Galactic Magnetic Field, Astrophys. J. Lett. 761, L11 (2012), arXiv:1210.7820 [astro-ph.GA] .
  * Wanderman and Piran [2015] D. Wanderman and T. Piran, The rate, luminosity function and time delay of non-Collapsar short GRBs, Mon. Not. Roy. Astron. Soc. 448, 3026 (2015), arXiv:1405.5878 [astro-ph.HE] .
  * Escorial  _et al._ [2023] A. R. Escorial _et al._ , The Jet Opening Angle and Event Rate Distributions of Short Gamma-Ray Bursts from Late-time X-Ray Afterglows, Astrophys. J. 959, 13 (2023), arXiv:2210.05695 [astro-ph.HE] .
  * Murase and Beacom [2010] K. Murase and J. F. Beacom, Very-High-Energy Gamma-Ray Signal from Nuclear Photodisintegration as a Probe of Extragalactic Sources of Ultrahigh-Energy Nuclei, Phys. Rev. D82, 043008 (2010), arXiv:1002.3980 [astro-ph.HE] .
  * Farrar [2024] G. R. Farrar, Binary neutron star mergers as the source of the highest energy cosmic rays (2024), arXiv:2405.12004 [astro-ph.HE] .
  * Globus  _et al._ [2023] N. Globus, A. Fedynitch, and R. D. Blandford, Treasure Maps for Detections of Extreme Energy Cosmic Rays, Astrophys. J. 945, 12 (2023), arXiv:2210.15885 [astro-ph.HE] .
  * Olive  _et al._ [2014] K. A. Olive _et al._ (Particle Data Group), Review of Particle Physics, Chin. Phys. C38, 090001 (2014).
  * Kim  _et al._ [2024] J. Kim _et al._ (Telescope Array), Highlights from the Telescope Array Experiment, PoS ICRC2023, 008 (2024).
  * Abbasi  _et al._ [2018] R. U. Abbasi _et al._ (Telescope Array), Depth of Ultra High Energy Cosmic Ray Induced Air Shower Maxima Measured by the Telescope Array Black Rock and Long Ridge FADC Fluorescence Detectors and Surface Array in Hybrid Mode, Astrophys. J. 858, 76 (2018), arXiv:1801.09784 [astro-ph.HE] .

##  V Supplemental Material

###  V.1 UH nuclei

The photodisintegration cross sections for UH nuclei derived from Talys 1.96
are shown in Fig. 5. We also added the photodisintegration cross section of
iron nuclei, for comparison purposes. We can see a reasonably good match
between the analytical estimate and the numerical results. In Fig. 6, we show
the energy loss length for three typical UH nuclei, namely Selenium (Se),
Tellurium (Te), and Platinum (Pt), respectively.

Figure 5: The photodisintegration cross section for UH nuclei is shown. Solid
lines are from the output of Talys 1.96, while the dashed lines are estimated
using the Delta-function approximation [59].

Figure 6: Energy loss length for photodisintegration (thick dashed lines),
photomeson production (dotted lines), Bethe-Heitler pair production (dotted-
dashed lines), and adiabatic expansion (black-dashed line) of UH nuclei of Se,
Te, and Pt as a function of energy. We also indicate the interaction length of
the photodisintegration process (thin dashed lines).

To study the change of composition of UH-UHECRs during propagation, we
estimate the fraction of survived nuclei $f_{\rm survival}=N_{\rm
obs}(E>E_{\rm obs},A>A_{\rm obs})/N_{\rm inj}(E_{\rm inj})$ numerically with
CRPropa 3 [e.g. 81], where $N_{\rm inj}$ is the number of nuclei injected at a
given distance, $N_{\rm obs}$ is the observed number of particles, $A_{\rm
obs}$ and $E_{\rm obs}$ are the observed nuclei mass number and energy,
respectively. The propagation distance $d_{95\%}$ can be inferred when the
value of $f_{\rm survival}=5\%$, which implies that $95\%$ of the nuclei are
lost during propagation. In Fig. 7, we show the propagation distance
$d_{95\%}$ for three typical UH nuclei species, namely Selenium (Se),
Tellurium (Te), and Platinum (Pt). For UH nuclei with primary energy $E_{\rm
inj}=300\rm~{}EeV$, the propagation distance is $d_{95\%}\sim 1.5\rm~{}Mpc$
for Se($A_{\rm obs}\geqslant 80$), $d_{95\%}\sim 6.5\rm~{}Mpc$ for Te($A_{\rm
obs}\geqslant 130$) and $d_{95\%}\sim 7.5\rm~{}Mpc$ for Pt($A_{\rm
obs}\geqslant 195$). The propagation distance increased to $d_{95\%}\sim
65\rm~{}Mpc$ for Se($A_{\rm obs}\geqslant 56$), $d_{95\%}\sim 320\rm~{}Mpc$
for Te($A_{\rm obs}\geqslant 56$) and $d_{95\%}\sim 480\rm~{}Mpc$ for
Pt($A_{\rm obs}\geqslant 56$).

Figure 7: Propagation distance $d_{95\%}$, corresponding to $f_{\rm
survival}=5\%$, is shown for Se, Te and Pt nuclei. The black-dashed lines
correspond to a value of 0.05 indicating that 95% of the nuclei are lost
during propagation.

###  V.2 Fitting procedure

Here, we describe the details of the fitting procedure to the observed UHECR
spectrum and composition. The observed flux is estimated with the following
formula [e.g., 66],

$\Phi_{A}=\sum_{A^{\prime}}\Phi_{AA^{\prime}}=\sum_{A^{\prime}}\frac{c}{4\pi}\frac{dn_{AA^{\prime}}}{dE},$
(5)

and

$\displaystyle dn_{AA^{\prime}}(E)$ $\displaystyle=$
$\displaystyle\int_{z_{\rm min}}^{z_{\rm
max}}dz\left|\frac{dt}{dz}\right|\xi(z)\rho_{0}$ (6) $\displaystyle\times$
$\displaystyle\int_{E^{\prime}_{\rm min}}^{E^{\prime}_{\rm
max}}dE^{\prime}\frac{dN_{A^{\prime}}}{dE^{\prime}}\eta_{AA^{\prime}}(E,E^{\prime},z),$

where $z$ is redshift, $\xi(z)$ is the redshift evolution factor, $\rho_{0}$
is the local event rate is the local event rate, and
$dt/dz=1/[H_{0}(1+z)(\sqrt{\Omega_{\Lambda}+\Omega_{k}(1+z)^{2}+\Omega_{m}(1+z)^{3}})]$.
We assume continuous source distribution with minimum source distance $d_{\rm
min}=1.5\rm~{}Mpc$ and maximum source distance $d_{\rm max}=3854\rm~{}Mpc$,
corresponding to redshift $z_{\rm min}=0.0003$ and $z_{\rm max}=4.5$. The
cosmological parameters adopted in this work are $H_{0}=67.3\ \rm km\ s^{-1}\
Mpc^{-1}$, $\Omega_{m}=0.315$, and $\Omega_{\Lambda}=0.685$ [82].

We adopt the chi-square method to fit the observed energy spectrum and
composition [33],

$\displaystyle\chi^{2}_{\text{spec}}=$ $\displaystyle\sum_{i}\frac{(f_{\rm
norm}\Phi^{\text{model}}(\hat{E}_{i};s,\mathcal{R}_{\text{max}})-\Phi^{\text{Auger/TA}}(E_{i}))^{2}}{\sigma_{i}^{2}}$
(7) $\displaystyle+\left(\frac{\delta_{E}}{\sigma_{E}}\right)^{2},$

where $f_{\rm norm}$ is a free normalization parameter, $\Phi^{\text{model}}$
is the simulation results, $\hat{E}_{i}\equiv(1+\delta_{E})E_{i}$, the range
of $\delta_{E}$ is from $-14\%$ to $14\%$, $\sigma_{E}$ is the systematic
uncertainty of the measured energy scale. The variation of the observed flux
$E^{3}\Phi$ has an uncertainty of $\sim 30\%$ because of the systematic
uncertainty $\sigma_{E}$ on the energy scale of measured particles [13].

Similarly, we estimate the $\chi^{2}$ value when fitting to the composition,
specifically, the first and second moments of $X_{\rm max}$ distribution,
$\langle X_{\text{max}}\rangle$ and $\sigma(X_{\text{max}})$.

$\displaystyle\chi^{2}_{\langle X_{\text{max}}\rangle}=$
$\displaystyle\sum_{i}\frac{(\langle
X_{\text{max}}^{\text{model}}\rangle(\hat{E}_{i};s,\mathcal{R}_{\text{max}})-\langle
X_{\text{max}}^{\text{Auger/TA}}\rangle(E_{i}))^{2}}{\sigma_{i}^{2}}.$ (8)

The total $\chi^{2}$ value value is

$\chi^{2}=\chi^{2}_{\text{spec}}+\chi^{2}_{\langle
X_{\text{max}}\rangle}+\chi^{2}_{\sigma(X_{\text{max}})},$ (9)

and the best-fit results are obtained when $\chi^{2}$ reaches its minimum
value $\chi^{2}_{\rm min}$. We fit the observed spectrum and composition with
energy above $10^{18.75}\rm~{}eV$ for Auger data, there are $16$ data points
on the spectrum and $10$ data points on the moments of $X_{\rm max}$
distribution. For TA data, we consider energy above $10^{18.45}\rm~{}eV$ with
$21$ data points on the spectrum and $10$ data points on the moments of
$X_{\rm max}$ distribution. The degree of freedom $\rm d.o.f.$ is defined as
the total number of data points minus the fitting parameter. In this work, we
treat the maximum particle energy $\mathcal{R}_{\rm max}$, spectral index
$s_{\rm CR}$ and nuclei number fraction $f_{A}$ as the free parameters. Note
$\sum f_{A}=1$.

The fitting parameters with Auger data shown in Fig. 2 is $\mathcal{R}_{\rm
max}=10^{18.3}\rm~{}eV$, $s_{\rm CR}=-0.3$, $f_{\rm P}=0.03$, $f_{\rm
He}=0.08$, $f_{\rm O}=0.87$, $f_{\rm Si}=0.08$, $f_{\rm Se}=0.007$ for Model A
with UH-UHECR nuclei Se, $\mathcal{R}_{\rm max}=10^{18.3}\rm~{}eV$, $s_{\rm
CR}=-0.4$, $f_{\rm P}=0.02$, $f_{\rm He}=0.009$, $f_{\rm O}=0.79$, $f_{\rm
Si}=0.17$, $f_{\rm Te}=0.007$ for Model A with UH-UHECR nuclei Te,
$\mathcal{R}_{\rm max}=10^{18.2}\rm~{}eV$, $s_{\rm CR}=-1.1$, $f_{\rm
P}=0.009$, $f_{\rm He}=0.1$, $f_{\rm O}=0.72$, $f_{\rm Si}=0.16$, $f_{\rm
Pt}=0.005$ for Model A with UH-UHECR nuclei Pt, $\mathcal{R}_{\rm
max}=10^{18.3}\rm~{}eV$, $s_{\rm CR}=-0.3$, $f_{\rm P}=0.02$, $f_{\rm
He}=0.25$, $f_{\rm O}=0.68$, $f_{\rm Si}=0.03$, $f_{\rm Fe}=0.007$, $f_{\rm
Se}=0.004$ for Model B with UH-UHECR nuclei Se, $\mathcal{R}_{\rm
max}=10^{18.3}\rm~{}eV$, $s_{\rm CR}=-0.7$, $f_{\rm P}=0.15$, $f_{\rm
He}=0.3$, $f_{\rm O}=0.48$, $f_{\rm Si}=0.057$, $f_{\rm Fe}=0.012$, $f_{\rm
Te}=0.0001$ for Model B with UH-UHECR nuclei Te, $\mathcal{R}_{\rm
max}=10^{18.3}\rm~{}eV$, $s_{\rm CR}=-0.3$, $f_{\rm P}=0.18$, $f_{\rm
He}=0.0009$, $f_{\rm O}=0.74$, $f_{\rm Si}=0.06$, $f_{\rm Fe}=0.012$, $f_{\rm
Pt}=0.0004$ for Model B with UH-UHECR nuclei Pt.

The fitting parameters with TA data shown in Fig. 3 is $\mathcal{R}_{\rm
max}=10^{18.5}\rm~{}eV$, $s_{\rm CR}=0.1$, $f_{\rm P}=0.8$, $f_{\rm He}=0.1$,
$f_{\rm O}=0.03$, $f_{\rm Si}=0.06$, $f_{\rm Se}=0.004$ for Model A with UH-
UHECR nuclei Se, $\mathcal{R}_{\rm max}=10^{18.7}\rm~{}eV$, $s_{\rm CR}=0.4$,
$f_{\rm P}=0.28$, $f_{\rm He}=0.001$, $f_{\rm O}=0.2$, $f_{\rm Si}=0.47$,
$f_{\rm Te}=0.005$ for Model A with UH-UHECR nuclei Te, $\mathcal{R}_{\rm
max}=10^{18.7}\rm~{}eV$, $s_{\rm CR}=0.4$, $f_{\rm P}=0.39$, $f_{\rm
He}=0.06$, $f_{\rm O}=0.04$, $f_{\rm Si}=0.5$, $f_{\rm Pt}=0.003$ for Model A
with UH-UHECR nuclei Pt, $\mathcal{R}_{\rm max}=10^{18.7}\rm~{}eV$, $s_{\rm
CR}=0.4$, $f_{\rm P}=0.4$, $f_{\rm He}=0.15$, $f_{\rm O}=0.2$, $f_{\rm
Si}=0.24$, $f_{\rm Fe}=0.002$, $f_{\rm Se}=0.01$ for Model B with UH-UHECR
nuclei Se, $\mathcal{R}_{\rm max}=10^{18.7}\rm~{}eV$, $s_{\rm CR}=0.5$,
$f_{\rm P}=0.35$, $f_{\rm He}=0.32$, $f_{\rm O}=0.01$, $f_{\rm Si}=0.3$,
$f_{\rm Fe}=0.001$, $f_{\rm Te}=0.0002$ for Model B with UH-UHECR nuclei Te,
$\mathcal{R}_{\rm max}=10^{18.7}\rm~{}eV$, $s_{\rm CR}=0.5$, $f_{\rm P}=0.4$,
$f_{\rm He}=0.05$, $f_{\rm O}=0.19$, $f_{\rm Si}=0.3$, $f_{\rm Fe}=0.00003$,
$f_{\rm Pt}=0.002$ for Model B with UH-UHECR nuclei Pt.

###  V.3 Constraints from the best-fit results with conventional nuclei

We try to constrain the energy generation density of UH-UHECRs, $Q_{\rm UH-
UHECR}$ based on the best-fit results with conventional nuclei. Specifically,
we begin with the best-fit results considering conventional nuclei only. Then,
we inject UH-UHECRs as an additional source population, but with the same
rigidity and spectral index, see Figure 2-3. We can see the obtained results
are consistent with the combined fit method, as shown in Table 2. With Model
A, we can see that the energy injection rate of the three ultra-heavy nuclei
species is constrained to be located in a narrow range with $Q_{\rm UH-
UHECR}\sim(4-8)\times 10^{42}\rm~{}erg~{}Mpc^{-3}~{}yr^{-1}$. However, we only
obtain the upper limits of the energy injection rate of UH-UHECRs when
adopting Model B, $Q_{\rm UH-UHECR}\lesssim(1.2-1.4)\times
10^{42}\rm~{}erg~{}Mpc^{-3}~{}yr^{-1}$. The reason is that including heavy
iron-group nuclei could provide a better fit for the observed energy spectrum,
with $\chi_{\rm tot,min}=64/{\rm d.o.f.}$ for conventional nuclei in Model B
and $\chi_{\rm tot,min}=72/{\rm d.o.f.}$ in Model A, where $\rm d.o.f.$ is the
degree of freedom.

Nuclei | $Q_{\rm UH-UHECR}^{\rm Auger}$ | $Q_{\rm UH-UHECR}^{\rm TA}$
---|---|---
| $[\rm erg~{}Mpc^{-3}~{}yr^{-1}]$ | $[\rm erg~{}Mpc^{-3}~{}yr^{-1}]$
Model A, fixed |  |
Se | $7.1_{-2.1}^{+1.6}\times 10^{42}$ | $2.2_{-1.1}^{+1.3}\times 10^{43}$
Te | $5.7_{-1.1}^{+0.0}\times 10^{42}$ | $1.7_{-0.8}^{+0.4}\times 10^{43}$
Pt | $5.0_{-0.6}^{+0.0}\times 10^{42}$ | $1.5_{-0.6}^{+0.4}\times 10^{43}$
Model B, fixed |  |
Se | $\lesssim 1.6\times 10^{42}$ | $\lesssim(3.8-13)\times 10^{42}$
Te | $\lesssim 1.4\times 10^{42}$ | $\lesssim(4.8-11)\times 10^{42}$
Pt | $\lesssim 1.2\times 10^{42}$ | $\lesssim(4.7-9)\times 10^{42}$
Table 2: Energy injection rate of UH-UHECRs.

Figure 8: Energy spectrum and the first/second moments of $X_{\rm max}$
distribution considering both conventional and UH-UHECRs. Left panel: Model A
with $E_{p,\rm max}=10^{18.1}\rm~{}eV$ and $s_{\rm CR}=-2$. Right panel: Model
B with $E_{p,\rm max}=10^{18.2}\rm~{}eV$ and ${}_{\rm CR}=-1.3$. The energy
spectrum of UHECRs is obtained from [64] and the first/second moments of
$X_{\rm max}$ distribution is taken from [65].

Figure 9: Same as Fig. 8, but for TA data. Left panel: Model A with $E_{p,\rm
max}=10^{18.7}\rm~{}eV$ and $s_{\rm CR}=0.2$. Right panel: Model B with
$E_{p,\rm max}=10^{18.6}\rm~{}eV$ and ${}_{\rm CR}=0$. The energy spectrum of
UHECRs is obtained from [83] and the first/second moments of $X_{\rm max}$
distribution is taken from [84].