# Inertial spin dynamics in epitaxial cobalt films

Vivek Unikandanunni  Department of Physics, Stockholm University, 10691
Stockholm, Sweden  Rajasekhar Medapalli  Center for Memory and Recording
Research, University of California San Diego, San Diego, CA 92093, USA
Department of Physics, Lancaster University, Bailrigg, Lancaster LA1 4YW,
United Kingdom  Marco Asa  Department of Physics, Politecnico di Milano
Technical University, Milano, Italy  Edoardo Albisetti  Department of Physics,
Politecnico di Milano Technical University, Milano, Italy  Daniela Petti
Department of Physics, Politecnico di Milano Technical University, Milano,
Italy  Riccardo Bertacco  Department of Physics, Politecnico di Milano
Technical University, Milano, Italy  Eric E. Fullerton  Center for Memory and
Recording Research, University of California San Diego, San Diego, CA 92093,
USA  Stefano Bonetti<EMAIL_ADDRESS>Department of Physics,
Stockholm University, 10691 Stockholm, Sweden  Department of Molecular
Sciences and Nanosystems, Ca’ Foscari University of Venice, 30172 Venice,
Italy

###### Abstract

We investigate the spin dynamics driven by terahertz magnetic fields in
epitaxial thin films of cobalt in its three crystalline phases. The terahertz
magnetic field generates a torque on the magnetization which causes it to
precess for about 1 ps, with a sub-picosecond temporal lag from the driving
force. Then, the magnetization undergoes natural damped THz oscillations at a
frequency characteristic of the crystalline phase. We describe the
experimental observations solving the inertial Landau-Lifshitz-Gilbert
equation. Using the results from the relativistic theory of magnetic inertia,
we find that the angular momentum relaxation time $\eta$ is the only material
parameter needed to describe all the experimental evidence. Our experiments
suggest a proportionality between $\eta$ and the strength of the magneto-
crystalline anisotropy.

††preprint: APS/123-QED

The fundamental understanding of magnetism has much improved since the first
experiments on magnetic materials at femto- and picosecond time scales, the
so-called ultrafast regime. The pioneering work of Beaurepaire et al.
Beaurepaire  _et al._ (1996) demonstrated that sub-picosecond magnetization
dynamics is possible, against the prediction of the textbook Landau-Lifshitz-
Gilbert (LLG) equation. Those results have since then been confirmed in
several experiments and using both optical Koopmans  _et al._ (2010); Kirilyuk
_et al._ (2010); Carva  _et al._ (2013); Bigot  _et al._ (2009); Kim  _et al._
(2012); Eschenlohr  _et al._ (2013); Bigot  _et al._ (2005); Koopmans  _et
al._ (2000); Hohlfeld  _et al._ (1997); Cinchetti  _et al._ (2006); Guidoni
_et al._ (2002) and X-ray techniques Stamm  _et al._ (2007); Boeglin  _et al._
(2010); Graves  _et al._ (2013); Iacocca  _et al._ (2019), making the field of
ultrafast magnetism, and of magnetization dynamics beyond the LLG equation, an
active field of research with possible implications for novel data storage
technologies.

Another recent analysis of the LLG equation showed that from a classical
mechanics point of view, the LLG equation assumes an unphysical inertial
tensor Ciornei  _et al._ (2011), also suggested by Gilbert himself in a
footnote almost two decades ago Gilbert (2004). However, Ciornei _et al._
Ciornei  _et al._ (2011) re-derived the equation with a realistic inertial
tensor which resulted in a slightly revised equation known as the inertial LLG
(iLLG) equation, which predicts the appearance of a nutation resonance at a
frequency much higher than the ferromagnetic resonance (FMR) one. After that
pioneering work, much theoretical effort has been performed trying to identify
the frequency regime where such nutation resonance was to be expected, with
predictions varying over a few orders of magnitude Fähnle  _et al._ (2011);
Böttcher and Henk (2012); Li  _et al._ (2015); Kikuchi and Tatara (2015);
Bastardis  _et al._ (2018). The experimental confirmation has been achieved
only recently and detected nutation dynamics in the $\sim 1$ THz range Neeraj
_et al._ (2021). Despite this novel experimental evidence, and the significant
theoretical progress to understand the microscopic origin of inertia
Bhattacharjee  _et al._ (2012); Olive  _et al._ (2012, 2015); Fransson  _et
al._ (2017); Mondal  _et al._ (2017); Cherkasskii  _et al._ (2020, 2021);
Thibaudeau and Nicolis (2021); Mondal and Oppeneer (2021); Mondal (2021);
Mondal  _et al._ (2021); Lomonosov  _et al._ (2021); Anders  _et al._ (2020);
Ruggeri (2021); Titov  _et al._ (2021a, b); Giordano and Déjardin (2020);
Rahman and Bandyopadhyay (2021); Gupta  _et al._ (2021); Jhuria  _et al._
(2020), a complete picture of how different material parameters affect the
nutation dynamics is still missing.

In this Letter, we provide the first experimental data on the dependence of
inertial spin dynamics on a key magnetic property, i.e. the magneto-
crystalline anisotropy. It has been suggested that spin-orbit coupling is the
fundamental interaction needed to derive a correct inertial tensor from first
principles. Hence, experiments where the magneto-crystalline anisotropy is
well defined and controlled, are expected to return important insights on this
open question. Similar to Ref. Neeraj  _et al._ (2021), we perform THz pump /
optical probe time-resolved magneto-optical measurements to trigger and detect
inertial spin dynamics. In contrast to the approach of Ref. Neeraj  _et al._
(2021), where a narrowband ($\Delta f/f\approx 0.1$) THz source was
implemented, we use intense single-cycle terahertz radiation Hoffmann and
Fülöp (2011) whose broad band ($\Delta f/f>1$) allows to cover the frequency
range where the nutation resonance is expected to appear. In addition, due to
the impulsive character of the driving pulse, our measurements are expected to
detect not only the forced response of the system, but also its natural one.
We model the magnetization dynamics solving the inertial LLG equation
numerically and we contextualize our results with the existing microscopic
theory of magnetic inertia.

Figure 1: (a) Geometry of THz pump-MOKE probe setup. (b) Frequency spectrum of
terahertz pump pulse (c) Magnetization loops for fcc, bcc, and hcp cobalt
measured using the longitudinal MOKE.

We choose to investigate three epitaxial cobalt thin films grown on MgO
substrate with face-centered cubic (fcc), body-centered cubic (bcc), and
hexagonal close packed (hcp) crystal structures. The fabrication details for
the three samples are given in the Supplemental Material SM . The hcp sample
has a strong in-plane magneto-crystalline anisotropy characterized by an easy
magnetization axis along the $c$-direction of the hcp structure (which lies in
the film plane in our films), and a hard axis orthogonal to it. For the two
cubic crystal structures, the anisotropy is still in-plane, but its strength
is much reduced, and a hard magnetization direction is not clearly identified
SM . The hcp and fcc samples were grown under similar deposition conditions
and are respectively 10 nm and 15 nm in thickness, while the bcc sample was
grown in a different laboratory and it has a thickness of 8 nm. Fig. 1(a)
shows the geometry of the single cycle THz pump - optical probe experiment.
The magnetization $\mathbf{M}$ of the sample is aligned along the
$x$-direction by means of an external bias field $|H_{ext}|=100$ mT, kept
constant during the experiment. The single-cycle THz pump pulses are generated
in the organic crystal OH1 by the optical rectification of 1300 nm radiation
from an optical parametric amplifier Jazbinsek  _et al._ (2019). The pump
pulse has a peak magnetic field of 0.3 T parallel to the $y$-direction, which
maximizes the torque on the magnetization, and impinges on the sample at an
angle of incidence $\theta_{inc}=45$ degrees. Fig. 1(b) is the Fourier
transform of the electro-optical sampling measurement in a a 50 $\mu$m-thick
GaP crystal Nahata  _et al._ (1996), used to characterize the THz pulse. It
shows that the pump field is peaked at around 2 THz, and has a bandwidth
exceeding 1 THz. The magnetization dynamics is probed using the time-resolved
magneto-optical Kerr effect (MOKE), in the specific measuring the Kerr
rotation angle $\theta_{K}$ of a nominally 40 fs, 800 nm probe beam, using a
balanced detection scheme. All radiation is derived from the same amplified
laser system ensuring intrinsic synchronization, with the relative delay
between the beams controlled by a mechanical translation stage, and with the
pump modulated at a frequency equal to half the laser repetition rate. Fig.
1(c) shows the easy axis magnetization loops for the three samples
investigated in this work. The coercive field for the hcp sample is about 50
mT, whereas is approximately 30 mT for both fcc and bcc samples SM .

Figure 2: Solid symbols: time-resolved Kerr rotation measurements on fcc, bcc
and hcp cobalt thin films. Dashed line: integral of the pump THz magnetic
field $H_{THz}$. Inset: zoomed-in main panel data for $t>1.7$ ps. The data is
shifted vertically for clarity. The continuous lines the best fits obtained
using Eq. (1).

Fig. 2 shows the time-resolved MOKE measurements of the terahertz-field
induced dynamics in all three samples. The plotted traces represent the
difference of the data recorded with magnetic fields of equal magnitude but
opposite polarity, ensuring the magnetic character of the signal Hudl  _et
al._ (2019). For all samples, the MOKE response is dominated by the coherent
precession of the magnetization around the applied THz magnetic field $H_{\rm
THz}$, which in fcc and bcc films is larger in amplitude than in the hcp one.
A very small demagnetization, showing up as a lingering non-zero average MOKE
signal, is also present. The presence of both coherent (precession) and
incoherent (demagnetization) effects in the observed THz-driven dynamics is
consistent with Ref. Bonetti  _et al._ (2016), where it was also observed that
in epitaxial films the demagnetization signal was negligible. Similar to Refs.
Bonetti  _et al._ (2016) and Neeraj  _et al._ (2021), no coherent precession
is observed when $H_{\rm THz}\parallel\mathbf{M}$, since in that case the
torque acting on the magnetization is zero. We present this measurement in the
Supplemental Material SM .

The dashed grey line in the same plot is the integral of $H_{\rm THz}$ over
time, obtained numerically from the electro-optic sampling measurement used to
characterize the THz field. Ref. Bonetti  _et al._ (2016) demonstrated that
the coherent response of the magnetization to an off-resonant THz field can be
obtained simply by integrating $H_{\rm THz}$, which is the solution of the LLG
equation for small and off-resonant excitations. However, while in Ref.
Bonetti  _et al._ (2016) the temporal overlap between the MOKE data and the
integral of $H_{\rm THz}$ was exact within the experimental error, here we
notice a substantial lag between them, approximately 200 fs for the fcc and
bcc samples, and 400 fs for the hcp one, highlighted by the vertical lines. In
other words, the shape of the MOKE response is still consistent with the
integral of the THz field if properly scaled, however its phase is not. This
phase shift is particularly dramatic for the hcp sample, where it looks as if
the magnetization precesses in the opposite direction as compared to the fcc
and bcc samples. We repeated the experiment on all samples, and the evidence
is robust. The sign of the magneto-optical coefficient does not change either
between the different samples, as demonstrated by the magneto-optical
hysteresis loops in Fig. 1(c). Hence, the observed phase shift is real and it
appears to be strongly dependent on the crystalline structure of the sample.

Another intriguing observation from the data in Fig. 2 is found in the inset,
where we zoom in on the main panel data at temporal delays $t>1.7$ ps. When
the pump field has left the sample, a comparatively tiny, yet detectable,
damped ringing of the magnetization can be observed. We can fit such behavior
with the phenomenological formula

$\Delta\theta_{K}(t)=Ae^{-t/\tau_{1}}+Be^{-t/\tau_{2}}\sin(2\pi ft)$ (1)

where $\tau_{1}$ is the recovery time of the incoherent demagnetization
dynamics, $\tau_{2}$ is the decay time of the sinusoidal oscillation, $f$ is
the frequency of the oscillations, and $A$, $B$ are the constants describing
the the amplitude of the demagnetization and, respectively, of the sinusoidal
oscillations. The fit returns $f_{fcc}\approx 1.3$ THz, $f_{bcc}\approx 1.4$
THz and $f_{hcp}\approx 2.1$ THz. $\tau_{2}$ is found to be approximately
0.82, 0.70 and 0.72 ps for the fcc, bcc and, respectively, hcp samples,
corresponding to damping coefficients $\alpha=1/\omega\tau_{2}$ which are
$\alpha_{fcc}\approx 0.15$, $\alpha_{bcc}\approx 0.16$ and
$\alpha_{hcp}\approx 0.10$.

Figure 3:  (a) Symbols: time-resolved Kerr signal at $t>1.7$ ps for THz
magnetic field values of different maximum amplitude. The data is vertically
shifted for clarity. Solid lines: best fit obtained using Eq. (1). (b)
Symbols: extracted oscillation amplitude $B$ as a function of THz magnetic
field and corresponding standard deviation. Dashed line: linear fit to the
data with imposed zero offset.

Before discussing these results, we present in Fig. 3 the THz field dependent
measurement on the fcc sample, which showed the largest signal in Fig. 2. The
terahertz field strength is controlled through the relative orientation of a
pair of wire-grid polarizers in the THz pump path. The second polarizer was
kept fixed in order to preserve the polarization of the THz field impinging on
the sample. The time-resolved MOKE signal is shown in Fig. 3(a) for the
maximum field strength, 75% and 50% of it, below which we were at the noise
level of our setup. We used again Eq. (1) to fit the oscillations and to
extract the amplitudes and recovery times as a function of THz field strength.
Fig. 3(b) shows the extracted oscillations amplitude $B$ as a function of THz
field strength, which can be fitted with a linear function with no offset.

The evidence presented so far is consistent with the presence of a sizeable
magnetic inertia in crystalline cobalt films, manifesting itself with a
lagging response to an external field and to the appearance of nutation
oscillations. In order to investigate this hypothesis thoroughly, we performed
numerical simulations using the inertial LLG equation, written in a slightly
different form than the one given in Ref. Neeraj  _et al._ (2021)

$\frac{d\mathbf{M}}{dt}=-\absolutevalue{\gamma}\mathbf{M}\times\mathbf{H}_{\rm
eff}+\mathbf{M}\times\left(\alpha\frac{d\mathbf{M}}{dt}-\eta\frac{d^{2}\mathbf{M}}{dt^{2}}\right),$
(2)

where $\absolutevalue{\gamma}$/$2\pi$ = 28 GHz/T is the gyromagnetic ratio,
$\mathbf{H}_{\rm eff}=(H_{\rm bias}+H_{K})\mathbf{x}+H_{\rm
THz}(t)\mathbf{y}+H_{d}\mathbf{z}$ is the effective magnetic field which
comprises of the external bias field $H_{\rm bias}$, the anisotropy field
$H_{K}$, the applied THz field $H_{\rm THz}(t)$, and the demagnetizing field
$H_{d}$; $M_{s}$ is the saturation magnetization of the sample, $\alpha$ is
the Gilbert damping parameter, and $\eta$ is the angular momentum relaxation
time defined as in Ref. Mondal  _et al._ (2021), i.e. $\eta=\alpha\tau$. Since
$\alpha\ll 1$, the absolute values $\eta$ are much smaller than $\tau$
reported in Ref. Neeraj  _et al._ (2021). The last term on the right hand side
of Eq. (2) is the nutation term that is present only when $\eta\neq 0$. In the
following, we solve this equation in the macrospin approximation (i.e. the
sample is considered as a homogeneous ferromagnet) and using a conventional
fourth-order Runge-Kutta method.

Figure 4:  Experimental (solid) and simulated (semi-transparent) Fourier
transform of the magnetization dynamics in (a) fcc (blue) (b) bcc (orange) and
(c) hcp (green) cobalt thin films. The full experimental trace is used for the
experimental data in Fig. 2, and the $M_{z}$ component for the simulations.
(d) Solid lines: simulated response of $M_{z}$ to a single cycle terahertz
field $H_{\rm THz}$ in the time-domain using the same parameters. Dashed line:
integral over time of $H_{\rm THz}$.

In Fig. 4(a)-(c), we compare the results from these numerical simulations to
the experimental results in frequency domain. For both simulations and
experiments, we Fourier transform the temporal traces obtained at time delays
$t>1.7$ ps, when the THz pump field has left the sample. Including the full
temporal trace would hide the small features below the broad single-cycle
response. In the numerical simulations, we calculate $\mathbf{H}_{\rm eff}$
solely from experimentally measured quantities found in the Supplemental
Material SM  or in previous references Unikandanunni  _et al._ (2021). This
allow us to estimate $H_{d}\approx 1.6$ T for all samples, and $H_{K}\approx
0.8$ T for the hcp sample and one order of magnitude smaller for the other two
samples. We used the nominal values for $H_{\rm bias}$ and $H_{\rm THz}$. The
only two free parameters are then $\eta$ and $\alpha$, which can be
independently tuned to match the peak frequency and, respectively, linewidth.
Using $\eta_{fcc}=120$ fs, $\eta_{bcc}=110$ fs and $\eta_{hcp}=75$ fs we can
reproduce the main experimental peak frequency, and assuming
$\alpha_{fcc}=0.15$, $\alpha_{bcc}=0.16$ and $\alpha_{hcp}=0.10$ from the fits
using Eq. (1), we can also match the linewidth of the main peak. No observable
difference was found within 5-10 fs for $\eta$ and within 0.01 for $\alpha$,
giving an approximate 10% relative accuracy. We have performed additional
simulations (not shown) and we also observe that the effective field does not
affect the nutation frequency and linewidth in a noticeable way unless it
reaches values of the order of a few Tesla.

We discuss below the plausibility of these values; assuming for the time being
that they are reasonable, and looking at Fig. 4(d), we obtain the remarkable
result that the inertial LLG equation is able to reproduce all the
experimental evidence of Fig. 2: the presence of a damped nutation oscillation
and the temporal shift of the coherent magnetization precession. In this small
amplitude limit, the inertial LLG also predicts a linear scaling of the
coherent precession and of the nutation amplitude with terahertz field
strength, as shown experimentally in Fig. 3(b). None of these experimental
evidences can be reproduced solving the standard LLG equation, proving that
the additional inertial term is necessary.

The only experimental evidence which is not reproduced by the inertial LLG
equation, in the currently known form and in the macrospin approximation, is
the presence of higher order harmonic peaks in the frequency response seen in
Fig. 4(a)-(c). We can clearly identify the second and third harmonics for the
fcc and bcc samples, and the second harmonic for the hcp one. However, Kikuchi
et al. Kikuchi and Tatara (2015) predicted such possibility if the third and
other higher order time derivatives of the magnetization, not included in the
standard framework of the inertial LLG model, are considered. We leave this
question open to future theoretical and experimental investigations, here we
simply note that the presence of harmonics at integer multiples $n$ of the
fundamental frequency could be consistent with nutation dynamics. It is not
consistent with the presence of standing waves across the film thickness,
which show instead a $n^{2}$ dependence due to confinement Razdolski  _et al._
(2017). We also do not observe any apparent inverse thickness dependence,
which is instead expected in the case of standing waves.

As a final control to test the general validity of our experimental results
and of the inertial LLG equation, we performed additional measurements using a
different THz single-cycle pump field with a bandwidth extending from 2 to 4
THz instead, i.e. with negligible overlap with the nutation resonances. This
is achieved by replacing the nonlinear crystal generating the THz radiation
and by adjusting the corresponding pump wavelength, leaving the rest of the
setup unchanged. The results are reported in the Supplemental Material SM ,
and they show that neither THz oscillations nor phase shift of the coherent
precession is observed when the pump field does not match the nutation
resonance. This is also in agreement with previous measurements done in fcc
cobalt driven by a THz field with similar bandwidth Shalaby  _et al._ (2018).
The inertial LLG equation with the same parameters reproduces even these
experimental data to an excellent degree, with no nutation oscillations nor
phase shift observed in this case.

We now turn the discussion to the two free parameters in the inertial LLG
equation, namely the damping $\alpha$ and the angular relaxation time $\eta$.
In order to match the experimental linewidth, we used for all three films a
damping parameter which is an order of magnitude larger than the typical FMR
Gilbert damping of these materials Unikandanunni  _et al._ (2021). While we do
not have a microscopic explanation for these large values, we notice that the
same issue was found in the first experimental report of nutation in
ferromagnets and left as an open question Neeraj  _et al._ (2021). Our
experiments, which are able to observe the natural nutation oscillations,
allow us to extract the damping factor directly from the data. In order to
estimate the magnitude of $\alpha$ for the inertial dynamics, a microscopic
theoretical investigation is needed, which is beyond the scope of our work.
Our experiments suggest though that a complete inertial LLG equation may
contain either distinct Gilbert and inertial damping coefficients, or a time-
dependent one, in order to fully describe the magnetization dynamics. A time-
dependent $\alpha$ can be qualitatively linked to a damping mechanism
dominated by comparatively stronger electron-phonon scattering at sub-
picosecond time scales, and weaker spin-lattice relaxation at longer time
scales Suhl (1998); Gilbert (2004); Ritzmann  _et al._ (2020). A time-
dependent $\alpha$ has also been recently suggested to include all damping
mechanisms, including time-retardation effects Bajpai and Nikolić (2019).

The most important experimental observation of this work, which is expected to
contribute to a microscopic understanding of inertial dynamics, is the strong
dependence of the nutation frequency on the different cobalt samples with
different magneto-crystalline anisotropy, which in turn is dependent on the
strength of the spin-orbit coupling. The relativistic theory of magnetic
inertia at ultrafast time scales Mondal  _et al._ (2017) demonstrated that the
presence of a finite angular momentum relaxation time is due to a spin-orbit
coupling effect of order $1/c^{4}$, whereas the Gilbert damping is of order
$1/c^{2}$. These two quantities are therefore dependent on each other, and it
was suggested in Ref. Mondal  _et al._ (2017) that the ratio $\eta/\alpha$
should be a constant. From our data, we calculate $\eta/\alpha=746\pm 46$ fs
from the three cobalt films, which is constant within the accuracy of our
estimates of the two parameters ($\sim 10\%$). This reduces the number of free
parameters in the inertial LLG equation to only one, at least within the same
$3d$ element, and it further strengthens the interpretation of our results in
terms of relativistic spin dynamics. We note that the magneto-crystalline
anisotropy energy is about one order of magnitude larger in hcp cobalt than in
the two cubic phases, while the nutation frequency differs by less than a
factor of two among them. This preliminary observation suggests that a
relation of proportionality may exist between the frequency of nutation and
the strength of the magneto-crystalline anisotropy, and that it is sub-linear.
We also note that fcc and bcc phases are energetically close; a small change
in lattice parameter can induce a so-called Bain transformation between them
Cuenya  _et al._ (2001). Hence, it is not too surprising that also their
magneto-crystalline anisotropy and nutation frequencies are similar. We
anticipate that future works will shed light on how to derive the nutation
frequency from first principles or from other magnetic properties of the
material.

In summary, we measured the temporal evolution of terahertz-field driven spin
dynamics in three epitaxial cobalt samples with fcc, bcc and hcp crystal
structures. We observed the appearance of THz oscillations with distinct
frequencies for the three samples and of a delayed coherent magnetization
response, which could be naturally described in the framework of the inertial
LLG equation assuming a magnetic damping one order of magnitude larger than
the conventional Gilbert damping at FMR frequencies. While surprising, this
evidence may be consistent with recent theoretical works suggesting a time-
dependent damping coefficient. We could also estimate a constant ratio between
the angular momentum relaxation time and the measured damping, in agreement
with the prediction of the full relativistic theory of magnetic inertia.
Finally, we could observe higher harmonics of the nutation oscillations, not
described by the currently accepted inertial LLG equation with temporal
derivatives up to the second order, but possibly consistent with a higher
order extension of the same equation. Our work provides the strongest evidence
for inertial spin dynamics so far, where all the experimental results can be
reproduced with a single free parameter. We envisage that our results will
trigger future experimental and theoretical investigations towards a deeper
microscopic understanding of magnetic inertia at ultrafast time scales.

V.U. and S.B. acknowledge support from the European Research Council, Starting
Grant 715452 “MAGNETIC-SPEED-LIMIT”. R.M. and E.E.F. were supported by U.S.
Department of Energy, Office of Science, Office of Basic Energy Sciences,
under Contract No. DE-SC0003678.

## References

  * Beaurepaire  _et al._ (1996) E. Beaurepaire, J.-C. Merle, A. Daunois, and J.-Y. Bigot, Physical review letters 76, 4250 (1996).
  * Koopmans  _et al._ (2010) B. Koopmans, G. Malinowski, F. Dalla Longa, D. Steiauf, M. Fähnle, T. Roth, M. Cinchetti, and M. Aeschlimann, Nature materials 9, 259 (2010).
  * Kirilyuk  _et al._ (2010) A. Kirilyuk, A. V. Kimel, and T. Rasing, Reviews of Modern Physics 82, 2731 (2010).
  * Carva  _et al._ (2013) K. Carva, M. Battiato, D. Legut, and P. M. Oppeneer, Physical Review B 87, 184425 (2013).
  * Bigot  _et al._ (2009) J.-Y. Bigot, M. Vomir, and E. Beaurepaire, Nature Physics 5, 515 (2009).
  * Kim  _et al._ (2012) J.-W. Kim, M. Vomir, and J.-Y. Bigot, Physical review letters 109, 166601 (2012).
  * Eschenlohr  _et al._ (2013) A. Eschenlohr, M. Battiato, P. Maldonado, N. Pontius, T. Kachel, K. Holldack, R. Mitzner, A. Föhlisch, P. M. Oppeneer, and C. Stamm, Nature materials 12, 332 (2013).
  * Bigot  _et al._ (2005) J.-Y. Bigot, M. Vomir, L. Andrade, and E. Beaurepaire, Chemical physics 318, 137 (2005).
  * Koopmans  _et al._ (2000) B. Koopmans, M. Van Kampen, J. Kohlhepp, and W. De Jonge, Physical Review Letters 85, 844 (2000).
  * Hohlfeld  _et al._ (1997) J. Hohlfeld, E. Matthias, R. Knorren, and K. Bennemann, Physical review letters 78, 4861 (1997).
  * Cinchetti  _et al._ (2006) M. Cinchetti, M. S. Albaneda, D. Hoffmann, T. Roth, J.-P. Wüstenberg, M. Krauß, O. Andreyev, H. Schneider, M. Bauer, and M. Aeschlimann, Physical review letters 97, 177201 (2006).
  * Guidoni  _et al._ (2002) L. Guidoni, E. Beaurepaire, and J.-Y. Bigot, Physical review letters 89, 017401 (2002).
  * Stamm  _et al._ (2007) C. Stamm, T. Kachel, N. Pontius, R. Mitzner, T. Quast, K. Holldack, S. Khan, C. Lupulescu, E. Aziz, M. Wietstruk, _et al._ , Nature materials 6, 740 (2007).
  * Boeglin  _et al._ (2010) C. Boeglin, E. Beaurepaire, V. Halté, V. López-Flores, C. Stamm, N. Pontius, H. Dürr, and J.-Y. Bigot, Nature 465, 458 (2010).
  * Graves  _et al._ (2013) C. Graves, A. Reid, T. Wang, B. Wu, S. De Jong, K. Vahaplar, I. Radu, D. Bernstein, M. Messerschmidt, L. Müller, _et al._ , Nature materials 12, 293 (2013).
  * Iacocca  _et al._ (2019) E. Iacocca, T.-M. Liu, A. H. Reid, Z. Fu, S. Ruta, P. Granitzka, E. Jal, S. Bonetti, A. Gray, C. Graves, _et al._ , Nature communications 10, 1 (2019).
  * Ciornei  _et al._ (2011) M.-C. Ciornei, J. Rubí, and J.-E. Wegrowe, Physical Review B 83, 020410 (2011).
  * Gilbert (2004) T. L. Gilbert, IEEE transactions on magnetics 40, 3443 (2004).
  * Fähnle  _et al._ (2011) M. Fähnle, D. Steiauf, and C. Illg, Physical Review B 84, 172403 (2011).
  * Böttcher and Henk (2012) D. Böttcher and J. Henk, Physical Review B 86, 020404 (2012).
  * Li  _et al._ (2015) Y. Li, A.-L. Barra, S. Auffret, U. Ebels, and W. E. Bailey, Physical Review B 92, 140413 (2015).
  * Kikuchi and Tatara (2015) T. Kikuchi and G. Tatara, Physical Review B 92, 184410 (2015).
  * Bastardis  _et al._ (2018) R. Bastardis, F. Vernay, and H. Kachkachi, Physical Review B 98, 165444 (2018).
  * Neeraj  _et al._ (2021) K. Neeraj, N. Awari, S. Kovalev, D. Polley, N. Z. Hagström, S. S. P. K. Arekapudi, A. Semisalova, K. Lenz, B. Green, J.-C. Deinert, _et al._ , Nature Physics 17, 245 (2021).
  * Bhattacharjee  _et al._ (2012) S. Bhattacharjee, L. Nordström, and J. Fransson, Physical review letters 108, 057204 (2012).
  * Olive  _et al._ (2012) E. Olive, Y. Lansac, and J.-E. Wegrowe, Applied Physics Letters 100, 192407 (2012).
  * Olive  _et al._ (2015) E. Olive, Y. Lansac, M. Meyer, M. Hayoun, and J.-E. Wegrowe, Journal of Applied Physics 117, 213904 (2015).
  * Fransson  _et al._ (2017) J. Fransson, D. Thonig, P. Bessarab, S. Bhattacharjee, J. Hellsvik, and L. Nordström, Physical Review Materials 1, 074404 (2017).
  * Mondal  _et al._ (2017) R. Mondal, M. Berritta, A. K. Nandy, and P. M. Oppeneer, Physical Review B 96, 024425 (2017).
  * Cherkasskii  _et al._ (2020) M. Cherkasskii, M. Farle, and A. Semisalova, Physical Review B 102, 184432 (2020).
  * Cherkasskii  _et al._ (2021) M. Cherkasskii, M. Farle, and A. Semisalova, arXiv preprint arXiv:2103.07864 (2021).
  * Thibaudeau and Nicolis (2021) P. Thibaudeau and S. Nicolis, arXiv preprint arXiv:2103.04787 (2021).
  * Mondal and Oppeneer (2021) R. Mondal and P. M. Oppeneer, arXiv preprint arXiv:2107.07939 (2021).
  * Mondal (2021) R. Mondal, Journal of Physics: Condensed Matter 33, 275804 (2021).
  * Mondal  _et al._ (2021) R. Mondal, S. Großenbach, L. Rózsa, and U. Nowak, Physical Review B 103, 104404 (2021).
  * Lomonosov  _et al._ (2021) A. M. Lomonosov, V. V. Temnov, and J.-E. Wegrowe, arXiv preprint arXiv:2105.07376 (2021).
  * Anders  _et al._ (2020) J. Anders, C. Sait, and S. Horsley, arXiv preprint arXiv:2009.00600 (2020).
  * Ruggeri (2021) M. Ruggeri, arXiv preprint arXiv:2103.09888 (2021).
  * Titov  _et al._ (2021a) S. Titov, W. Coffey, Y. P. Kalmykov, and M. Zarifakis, Physical Review B 103, 214444 (2021a).
  * Titov  _et al._ (2021b) S. Titov, W. Coffey, Y. P. Kalmykov, M. Zarifakis, and A. Titov, Physical Review B 103, 144433 (2021b).
  * Giordano and Déjardin (2020) S. Giordano and P.-M. Déjardin, Physical Review B 102, 214406 (2020).
  * Rahman and Bandyopadhyay (2021) R. Rahman and S. Bandyopadhyay, Journal of Physics: Condensed Matter (2021).
  * Gupta  _et al._ (2021) R. Gupta, S. Husain, A. Kumar, R. Brucas, A. Rydberg, and P. Svedlindh, Advanced Optical Materials 9, 2001987 (2021).
  * Jhuria  _et al._ (2020) K. Jhuria, J. Hohlfeld, A. Pattabi, E. Martin, A. Y. A. Córdova, X. Shi, R. L. Conte, S. Petit-Watelot, J. C. Rojas-Sanchez, G. Malinowski, _et al._ , Nature Electronics 3, 680 (2020).
  * Hoffmann and Fülöp (2011) M. C. Hoffmann and J. A. Fülöp, Journal of Physics D: Applied Physics 44, 083001 (2011).
  * (46) See Supplemental Material for (i) magnetic and structural sample characterization and (ii) supplemental pump-probe measurements, which include Refs. Unikandanunni  _et al._ (2021); Bonetti  _et al._ (2016).
  * Jazbinsek  _et al._ (2019) M. Jazbinsek, U. Puc, A. Abina, and A. Zidansek, Applied Sciences 9, 882 (2019).
  * Nahata  _et al._ (1996) A. Nahata, D. H. Auston, T. F. Heinz, and C. Wu, Applied physics letters 68, 150 (1996).
  * Hudl  _et al._ (2019) M. Hudl, M. d’Aquino, M. Pancaldi, S.-H. Yang, M. G. Samant, S. S. Parkin, H. A. Dürr, C. Serpico, M. C. Hoffmann, and S. Bonetti, Physical review letters 123, 197204 (2019).
  * Bonetti  _et al._ (2016) S. Bonetti, M. Hoffmann, M.-J. Sher, Z. Chen, S.-H. Yang, M. Samant, S. Parkin, and H. Dürr, Physical review letters 117, 087205 (2016).
  * Unikandanunni  _et al._ (2021) V. Unikandanunni, R. Medapalli, E. E. Fullerton, K. Carva, P. M. Oppeneer, and S. Bonetti, Applied Physics Letters 118, 232404 (2021).
  * Razdolski  _et al._ (2017) I. Razdolski, A. Alekhin, N. Ilin, J. P. Meyburg, V. Roddatis, D. Diesing, U. Bovensiepen, and A. Melnikov, Nature communications 8, 1 (2017).
  * Shalaby  _et al._ (2018) M. Shalaby, A. Donges, K. Carva, R. Allenspach, P. M. Oppeneer, U. Nowak, and C. P. Hauri, Physical Review B 98, 014405 (2018).
  * Suhl (1998) H. Suhl, IEEE Transactions on Magnetics 34, 1834 (1998).
  * Ritzmann  _et al._ (2020) U. Ritzmann, P. M. Oppeneer, and P. Maldonado, Physical Review B 102, 214305 (2020).
  * Bajpai and Nikolić (2019) U. Bajpai and B. K. Nikolić, Physical Review B 99, 134409 (2019).
  * Cuenya  _et al._ (2001) B. R. Cuenya, M. Doi, S. Löbus, R. Courths, and W. Keune, Surface science 493, 338 (2001).


used $\delta t/t_{R}=0.1$, where $t_{R}$ is the duration of the Ramsey
measurement. As everywhere else, we used the relative length of the cycle
$t_{\mathrm{cyc}}/t_{R}=3$

At each time step a TLS can switch between its states $\ket{0}^{(n)}$ and
$\ket{1}^{(n)}$. The switching probabilities are

$\displaystyle p_{01}^{(n)}=W_{01}^{(n)}\delta t,\qquad
p_{10}^{(n)}=W_{10}^{(n)}\delta t\ .$ (95)

where $W_{01}^{(n)}$ and $W_{10}^{(n)}$ are the switching rates. We determined
numerically whether the TLS switches or not in a standard way by comparing
$W_{ij}^{(n)}\delta t$ with a random number from the uniform distribution
$U(0,1)$.

The generated states of the TLS produce an array $d^{(n)}(m)$ of random
numbers that take values $\pm 1$. Here $m$ enumerates the time steps. We
collect $N=10^{5}$ outcomes of the simulated Ramsey measurements , which means
that we use $\tilde{N}=(t_{\mathrm{cyc}}/\delta t)\times 10^{5}$ samples, that
is, $1\leq m\leq\tilde{N}$. A $k$th Ramsey measurement, which is done in the
time interval $kt_{\mathrm{cyc}}\leq t\leq kt_{\mathrm{cyc}}+t_{R}$,
corresponds to the range of steps $k(t_{\mathrm{cyc}}/\delta t)<m\leq
k(t_{\mathrm{cyc}}/\delta t)+(t_{R}/\delta t)$. For our discretized time
sequence, the random phase accumulated by the qubit in this time interval is

$\displaystyle\theta{}_{k}=\sum_{m=1}^{\lceil t_{R}/\delta
t\rceil}\sum_{n=1}^{N_{\mathrm{TLS}}}V^{(n)}d^{(n)}(\lfloor
kt_{\mathrm{cyc}}/\delta t\rfloor+m).$ (96)

We used the probability of obtaining “1” in a $k$th measurement
$p(\theta_{k})$ as given by Eq. (7). We compared $p(\theta_{k})$ with a random
numbers $o_{k}$ from $U(0,1)$. If $p(\theta_{k})>o_{k}$, we set the outcome of
the $k$th Ramsey measurement to $x_{k}=1$, otherwise we set $x_{k}=0$. The
whole procedure was independently repeated 300 times for statistical
averaging. This allowed us to numerically analyze and compare the parameters
$r_{1},r_{2}(k),r_{3}(k,l)$ with the theory, as well as to investigate other
parameters of interest, as discussed in the main text.

Of primary interest to us was the analysis of the TLSs that produce noise of
the $1/f$ type in a reasonably broad frequency range. There are many ways to
obtain such noise. The results presented in the main text refer to the noise
in which the coupling of the TLSs to the qubit was the same for all TLSs,
$V^{(n)}=V$, but the distribution of the switching rates
$W^{(n)}=W_{01}^{(n)}+W_{10}^{(n)}$ was log-normal. We also simulated other
types of the TLSs. They lead to qualitatively similar results.

###  B.2 Simulating Gaussian noise

The effect of a Gaussian noise on the outcomes of Ramsey measurements is fully
characterized by the correlation function
$f_{k}=\langle\theta_{n}\theta_{n+k}\rangle$ of the phases $\theta{}_{k}$
acquired by the qubit in the measurements. Since the measurements are
periodically repeated, one has to sample $\theta{}_{k}$ for successive $k$.
The probability to have a given phase $\theta_{k}$ depends on the entire
history of the previously “observed” phases $\theta_{k^{\prime}}$ with
$k^{\prime}<k$. This means that the quantity of interest is the conditional
probability

$\displaystyle
P(\theta{}_{k}|\theta{}_{0},...,\theta{}_{k-1})=P(\theta{}_{0},...,\theta{}_{k})$
$\displaystyle\times[P(\theta{}_{0},...,\theta{}_{k-1})]^{-1}$ (97)

where $\theta_{0}$ is the outcome of the first measurement. The probability
(B.2) has to be evaluated recursively starting with the probability of
$\theta_{0}$. We will use that the distribution of the phases is stationary,

$P(\\{\theta\\})=Z^{-1}\exp\left[-\frac{1}{2}\sum_{k,k^{\prime}}(\hat{f}{}^{-1})_{k-k^{\prime}}\theta_{k}\theta_{k^{\prime}}\right],$

where $Z$ is the normalization constant and $\hat{f}{}^{-1}$ is the matrix
reciprocal to the matrix $f_{kk^{\prime}}$; we note that the matrix elements
of the latter matrix are $f_{kk^{\prime}}=f_{|k-k^{\prime}|}$, cf. Eq. (X).

In evaluating the conditional probability given by Eq. (B.2) one should keep
in mind that the values of $\theta{}_{k}$ are correlated at a finite distance
$k_{\mathrm{corr}}$, which is determined by the relation between the
correlation time of the underlying noise and the period of the sequence
$t_{\mathrm{cyc}}$. In other words, it means that, to a good approximation
(which needs to be checked) $|f_{k}|$ can be set equal to zero for
$k>k_{\mathrm{corr}}$. Then for $k>k_{\mathrm{corr}}$ one can approximate the
conditional probability $P(\theta{}_{k}|\theta{}_{0},...,\theta{}_{k-1})$ with
$P(\theta{}_{k}|\theta{}_{k-k_{\mathrm{corr}}},\theta{}_{k-k_{\mathrm{corr}}+1},...,\theta{}_{k-1})$,
i.e., instead of the probability (B.2) we have to calculate

$\displaystyle
P(\theta{}_{k}|\theta{}_{k-k_{\mathrm{corr}}},...,\theta{}_{k-1})=P(\theta{}_{k-k_{\mathrm{corr}}},...,\theta{}_{k})$
$\displaystyle\times[P(\theta{}_{k-k_{\mathrm{corr}}},...,\theta{}_{k-1})]^{-1}.$
(98)

It is important that all conditional probabilities in Eq. (B.2) have a
Gaussian form, albeit they are not zero-mean, because of the correlations. As
we now show, we can sample each $\theta{}_{k}$ from
$\mathcal{N}(\mu_{k},\sigma_{k})$, where $\mu_{k}$ and $\sigma_{k}$ are
respectively the mean value and the standard deviation for $\theta_{k}$. They
depend on the values of $\theta_{k^{\prime}}$ with $k^{\prime}<k$.

The first phase to be sampled, $\theta{}_{0}$, is sampled with $\mu_{0}=0$ and
$\sigma_{0}=1/\sqrt{f_{0}}$. To find the distribution of the phases
$\theta{}_{k}$ with $k>0$ we note that, when calculating
$P(\theta{}_{k-k_{\mathrm{corr}}},...,\theta{}_{k})$, rather than using the
full reciprocal matrix $(\hat{f}{}^{-1})_{|k-k^{\prime}|}$ we should use a
$(k_{\mathrm{corr}}+1)\times(k_{\mathrm{corr}}+1)$ matrix reciprocal to the
$(k_{\mathrm{corr}}+1)\times(k_{\mathrm{corr}}+1)$ part of the matrix
$f_{kk^{\prime}}$. This matrix $\Psi_{kk^{\prime}}$ is defined by the equation

$\sum_{m=k-k_{\mathrm{corr}}}^{k}\Psi_{km}f_{mk^{\prime}}=\delta_{kk^{\prime}},$

Along with $\Psi_{kk^{\prime}}$ we need the matrix $\psi_{kk^{\prime}}$, which
is the reciprocal of the $k_{\mathrm{corr}}\times k_{\mathrm{corr}}$ part of
$f_{kk^{\prime}}$,

$\sum_{m=k-k_{\mathrm{corr}}}^{k-1}\psi_{km}f_{mk^{\prime}}=\delta_{kk^{\prime}}.$

The matrices $\Psi_{kk^{\prime}}$ and $\psi_{kk^{\prime}}$ are symmetric.
However, even though $f_{kk}=f_{0}$ is independent of $k$, the diagonal matrix
elements of the matrices $\Psi_{kk}$ and $\psi_{kk}$ depend on $k$.

There is an important relation between the matrices $\hat{\psi}$ and
$\hat{\Psi}$,

$\displaystyle\psi_{kk^{\prime}}=\Psi_{kk^{\prime}}-(\Psi_{mk}\Psi_{mk^{\prime}}/\Psi_{mm})$
$\displaystyle k,k^{\prime}=m-k_{\mathrm{corr}},\ldots,m-1;$ (99)

This relation can be checked by multiplying from the left by $f_{k_{1}k}$ and
summing over $k=m-k_{\mathrm{corr}},\ldots,m-1$.

Taking the relation (B.2) into account, we can write the exponential in
$P(\theta{}_{k-k_{\mathrm{corr}}},...,\theta{}_{k})$ as

$\displaystyle\exp\left[-\frac{1}{2}\Psi_{kk}(\theta{}_{k}-\mu_{k})^{2}\right.$
$\displaystyle\left.-\frac{1}{2}\sum_{m,m^{\prime}=k-k_{\mathrm{corr}}}^{k-1}\psi_{mm^{\prime}}\theta{}_{m}\theta{}_{m^{\prime}}\right],$

where

$\displaystyle\mu_{k}=-\Psi_{kk}^{-1}\sum_{m=k-k_{\mathrm{corr}}}^{k-1}\Psi_{km}\theta{}_{m}$
(100)

is the phase accumulated over $k_{\mathrm{corr}}$ steps that preceded the
$k$th step.

Ultimately, we have for the conditional probability of $\theta{}_{k}$ the
distribution

$\displaystyle
P(\theta{}_{k}|\theta{}_{k-k_{\mathrm{corr}}},\theta{}_{k-k_{\mathrm{corr}}+1},...,\theta{}_{k-1})$
$\displaystyle=(\Psi_{kk}/2\pi)^{1/2}\exp\left[-\frac{1}{2}\Psi_{kk}(\theta{}_{k}-\mu_{k})^{2}\right]$
(101)

We have also used here the Cramer rule that relates the matrix element
$\Psi_{kk}$ to the ratio of the determinants of the matrices $\hat{\Psi}$ and
$\hat{\psi}$.

The above prescription allowed us to sample a sequence of random phases
$\theta{}_{k}$. Each obtained $\theta{}_{k}$ was used to determine whether the
outcome of the simulated Ramsey measurement gives “0” or “1” based on the
probability (7). From the observed outcomes, we could calculate
$r_{1},r_{2}(k)$, and $r_{3}(k,l)$ as well as other statistical
characteristics of the simulated sequence of periodic Ramsey measurements.

The value of $k_{\mathrm{corr}}$ depends on a particular type of noise. We
chose it in such a way that the results became virtually independent of the
choice. For exponentially correlated frequency noise with the correlation time
$\tau_{\mathrm{corr}}$ one can choose
$k_{\mathrm{corr}}=a_{\mathrm{corr}}\tau_{\mathrm{corr}}/t_{\mathrm{cyc}}$
with a sufficiently large $a_{\mathrm{corr}}$. For the $1/f$-type noise we
studied, with the characteristic minimal frequency $\omega_{\min}$, one can
choose $k_{\mathrm{corr}}=a_{1/f}(\omega_{\min}t_{\mathrm{cyc}})^{-1}$ with a
sufficiently large $a_{1/f}$. However, for the parameters used in the
simulations it was checked that the results become independent of
$k_{\mathrm{corr}}$ once one sets $k_{\mathrm{corr}}\approx 10$ for the
exponentially correlated noise and $k_{\mathrm{corr}}=400$ for the $1/f$
noise.

##  Appendix C Power spectrum of the TLS-induced frequency fluctuations

The power spectrum of the qubit frequency fluctuations due to the coupling to
TLSs becomes particularly simple if the TLSs are symmetric, the coupling is
the same for all TLSs, and the distribution of the switching rates is log-
normal, i.e., the probability density of having a given switching rate $W$ is
$\propto 1/W$. Such a distribution is natural if one thinks of a uniform
distribution of the barriers which the TLSs have to overcome in switching. For
a finite set of the TLSs the log-normal distribution is mimicked by the
switching probabilities of the form

$W^{{}^{(n)}}t_{R}=\exp\left(-\alpha(n+n_{0})\right),$ (102)

where $\alpha$ and $n_{0}$ are control parameters. In the simulations we used
$\alpha=0.75$ and considered $n_{0}=0,5$ or $7$.

The selected switching rates correspond to the power spectrum $S_{q}(\omega)$,
which is a superposition of Lorentzian peaks given by Eq. (V.1). This
superposition provides a $1/f$-type behavior. It is clearly seen in Fig. 7
already for 10 TLSs, in which case the $1/f$ spectrum spreads over two decades
in frequency. The range of $1/f$-type behavior can be extended by altering
$\alpha$ and $n_{0}$ or by including more TLSs. The inset in Fig. 7 shows
that, for 20 TLSs, $1/f$-behavior spreads over five decades.

Figure 7: Simulated power spectrum of 10 symmetric TLSs with
$V^{(n)}=V=0.2/t_{R}$ and $W_{01}^{{}^{(n)}}=W_{10}^{{}^{(n)}}=\exp(-\alpha
n)/t_{R}$ for $n=1,2,\ldots,10$ and $\alpha=0.75$ (green data points). The
spectrum displays $1/f$ character over more than two decades. The results
agree with the analytical result (V.1). Inset: the power spectrum for 20 TLSs,
$n=1,...,20$ with the same $V$ and $\alpha$ as given by Eq. (V.1). In this
case the spectrum displays more than five decades of $1/f$ behavior. The
dashed lines in the main plot and in the inset show $1/f$ noise.

The spectrum of the main plot in Fig. 7 is obtained via sampling the TLSs
states $d^{{}^{(n)}}(m)=\pm 1$ with a time step $\delta t$. At each time step
$m$ the outcome is

$D(m)=\sum_{n=1}^{N_{\mathrm{TLS}}}V^{{}^{(n)}}t_{R}d^{{}^{(n)}}(m).$ (103)

The fast Fourier transform of $D(m)$ is calculated in the standard way as

$\displaystyle\tilde{D}(\omega_{k})=\left|\sum_{m}\exp(-2\pi
imk/\tilde{N})D(m)\right|^{2}/\tilde{N}$ (104)

where $\omega_{k}=2\pi k/\tilde{N}$ and, as before,
$\tilde{N}=Nt_{\mathrm{cyc}}/\delta t$ with $N=10^{5}$.

To obtain the power spectrum
$S_{q}(\omega_{k})=\langle\tilde{D}(\omega_{k})\rangle$, this procedure is
repeated 300 times to get sufficient statistics for averaging. The result is
shown in the main plot in Fig. 7. For weak couplng to the TLSs the spectrum
$S_{q}(\omega)$ is immediately related to the centered correlator
$\tilde{r}_{2}(k)$.

##  Appendix D Broadening of the distribution by the noise correlations

Noise correlations lead to broadening of the distribution of the measurement
outcomes compared to the standard binomial distribution. Here we illustrate
this well-known effect for our setting. For a random variable $x_{n}$ that
takes on the values 0 and 1, the probability of observing $m$ ones in $M$
measurements is

$\langle m/M\rangle\equiv
M^{-1}\mathbb{E}\left[\sum_{n=1}^{M}x_{n}\right]=r_{1},$ (105)

whereas the variance is

$\displaystyle\sigma_{M}^{2}=\langle(m/M)^{2}\rangle-r_{1}^{2}\equiv
M^{-2}\mathbb{E}\left[\left(\sum_{n=1}^{M}x_{n}\right)^{2}\right]-r_{1}^{2}$
$\displaystyle=M^{-1}r_{1}(1-r_{1})+2M^{-2}\sum_{k=1}^{M-1}(M-k)\tilde{r}_{2}(k).$
(106)

For the case where $Mt_{\mathrm{cyc}}$ is much larger than the correlation
time of the noise one should replace in the second line $M-k$ with $M$.

The term $r_{1}(1-r_{1})/M$ in Eq. (D) is the standard variance of the
binomial distribution. It falls off as $M^{-1}$ with the increasing number of
measurements $M$. Equation (D) shows that the correction to the binomial
distribution due to the noise correlations, which is determined by the
centered correlator $\tilde{r}_{2}(k)$, also scales as $M^{-1}$. This
correction can be significant. The change of the standard deviation due to the
noise correlation for several types of noise that we studied is seen from
Table 1.

Noise type | Setup | $M^{1/2}\sigma_{\mathrm{simul}}$ | $M^{1/2}\sigma_{\mathrm{binom}}$
---|---|---|---
TLS | $n_{0}=5$ | 1.283 | 0.379
TLS | $n_{0}=7$ | 1.696 | 0.379
Gaussian EC | $\tau_{\mathrm{corr}}/t_{R}=20$ | 0.605 | 0.379
Gaussian EC | $\tau_{\mathrm{corr}}/t_{R}=100$ | 1.120 | 0.379
Gaussian $1/f$ | $\omega_{\mathrm{min}}t_{R}=1\times 10^{-3}$ | 0.956 | 0.379
Gaussian $1/f$ | $\omega_{\mathrm{min}}t_{R}=1\times 10^{-5}$ | 1.378 | 0.379
Table 1: Standard deviations of the simulated probabilities of the measurement
outcomes scaled by the number of measurements $M^{1/2}\sigma_{\mathrm{simul}}$
and of the binomial distributions $M^{1/2}\sigma_{\mathrm{binom}}$ for several
types of noise. The results for the noise from the TLS refer to 4 symmetric
TLSs with scaled coupling $V^{(n)}t_{R}=Vt_{R}=0.2$ and the switching rates
$W^{{}^{(n)}}t_{R}=\exp[-0.75(n+n_{0})]$ ($n=1,2,3,4$). The results for the
Gaussian exponentially correlated noise (Gaussian EC) and the $1/f$ noise with
the spectrum (IV.4) refer to the noise intensities that give the phase
variance $\langle\theta_{0}^{2}\rangle\equiv f_{0}=0.16$.

Figure 8: Broadening of the distribution of the outcomes of $M=10^{4}$
measurements. Panels (a) and (b) refer to the exponentially correlated and
$1/f$ noise, respectively; the noise intensity is chosen such that the phase
variance is $\langle\theta_{0}^{2}\rangle\equiv f_{0}=0.16$. (c) noise from 4
TLSs with the same parameters as in Table 1. Solid lines show the Gaussian
distribution, which approximates the binomial distribution with $r_{1}=0.826$
and $M=10^{4}$. Each histogram is aggregated into 100 bins, where the lower
and upper ends are determined by the minimal and maximal sampled values,
respectively. The theoretical curves use the same binning.

The effect of the correlations on the variance can be studied by evaluating
the probability distribution $\rho(m|M)$ of having $m$ ones in $M$
measurements. This distribution for $M=10^{4}$ is shown in Fig. 8 for several
types of noise. It was obtained by repeating $10^{4}$ measurements $10^{5}$
times to have a good statistics. The distribution is compared with the
binomial distribution for the same $M=10^{4}$, which is very well approximated
by a Gaussian distribution $\propto\exp[-(m-Mr_{1})^{2}/2Mr_{1}(1-r_{1})]$.

##  Appendix E TLS dynamics during the Ramsey measurement

We consider here the dynamics of the TLS-dependent components of the density
matrix during the Ramsey measurement, i.e., in the time intervals
$nt_{\mathrm{cyc}}<t<nt_{\mathrm{cyc}}+t_{R}$. These components,
$\rho_{\lambda}$ with $\lambda=I,z,\pm$, are defined in Eqs. (VIII.2) and
(VIII.2). As indicated in the main text, the equations for $\rho_{\lambda}$
with different $\lambda$ are obtained by substituting Eq. (VIII.2) into the
full master equation (VIII.2), multiplying the left- and right-hand sides by
$\hat{I},\hat{I}_{q}-\sigma_{z},\sigma_{\pm}$, and taking trace over the qubit
states. The left-hand side of the resulting equation for $\rho_{\lambda}$ is

$\displaystyle\partial_{t}\rho_{\lambda}=\sum_{n}\sum_{\varkappa=0,z}\dot{C}_{\lambda\varkappa}^{(n)}{}\hat{\tau}_{\varkappa}^{(n)}\prod_{m\neq
n}\sum_{\varkappa^{\prime}=0,z}C_{\lambda\varkappa^{\prime}}^{(m)}\hat{\tau}_{\varkappa^{\prime}}^{(m)}.$
(107)

As seen from Eq. (VIII.2, the term $\sum_{n}\mathcal{L}^{(n)}\rho_{\lambda}$
in the right-hand side of the equation for $\partial_{t}\rho_{\lambda}$ has
the same structure as Eq. (107): a sum over $n$ multiplied by the product over
$m\neq n$ of
$\sum_{\varkappa^{\prime}=0,z}C_{\lambda\varkappa^{\prime}}^{(m)}\hat{\tau}_{\varkappa^{\prime}}^{(m)}$.
The commutator $[\rho,H_{\mathrm{q-TLS}}]$ also has the same structure. One
can divide the both sides of the equation for $\partial_{t}\rho_{\lambda}$ by
$\rho_{\lambda}$, reminiscent of the standard trick of separation of variables
in a differential equation. This gives equations (VIII.2).

The equations for $C_{\lambda\varkappa}^{(n)}$ are split into sets of pairs of
coupled equations. For $\lambda=0,z$ we have

$\displaystyle\dot{C}_{\lambda 0}^{(n)}=0,\quad\dot{C}_{\lambda
z}^{(n)}=\Delta W^{(n)}{}C_{\lambda 0}^{(n)}{}-W^{(n)}{}C_{\lambda z}^{(n)}{}\
,$ $\displaystyle\Delta W^{(n)}{}=W_{10}^{(n)}{}-W_{01}^{(n)}{}\ .$ (108)

The parameter $\Delta W^{(n)}{}$ characterizes the asymmetry of the TLS,
whereas $W^{(n)}{}=W_{10}^{(n)}{}+W_{01}^{(n)}{}$ is the TLS relaxation rate.

If $\lambda=+$, Eq. (VIII.2) for $C_{\lambda\varkappa}^{(n)}$ reads

$\displaystyle\dot{C}_{+0}^{(n)}{}=iV^{(n)}{}C_{+z}^{(n)}{},$
$\displaystyle\dot{C}_{+z}^{(n)}{}=(\Delta
W^{(n)}{}+iV^{(n)}{})C_{+0}^{(n)}{}-W^{(n)}{}C_{+z}^{(n)}{}\ ,$ (109)

and by construction, $C_{-\varkappa}^{(n)}=C_{+\varkappa}^{(n)}{}^{*}$.

Solving these equations, we obtain for $\lambda=I,z$

$\displaystyle C_{\lambda z}^{(n)}{}(t)=\left[C_{\lambda
z}^{(n)}{}(0^{+})-\frac{\Delta W^{(n)}{}}{W^{(n)}{}}C_{\lambda
0}^{(n)}{}(0^{+})\right]e^{-W^{(n)}{}t}$ $\displaystyle+\frac{\Delta
W^{(n)}{}}{W^{(n)}{}}C_{\lambda 0}^{(n)}{}(0^{+}),\quad C_{\lambda
0}^{(n)}{}(t)=C_{\lambda 0}^{(n)}{}(0^{+})\ ,$ (110)

whereas

$\displaystyle
C_{+\varkappa}^{(n)}{}(t)=\sum_{k=1,2}B_{\varkappa}^{(k;n)}\exp(\nu_{k}^{(n)}{}t),$
$\displaystyle\nu_{1,2}^{(n)}{}=-\frac{1}{2}W^{(n)}{}\pm\gamma^{(n)},$
$\displaystyle\gamma^{(n)}=\frac{1}{2}\left[(W^{(n)}{})^{2}+4iV^{(n)}{}(\Delta
W^{(n)}{}+iV^{(n)}{})\right]^{1/2}\ ,$ (111)

here

$\displaystyle
B_{0}^{(k;n)}=[\nu_{3-k}^{(n)}C_{+0}^{(n)}{}(0^{+})-iV^{(n)}{}C_{+z}^{(n)}{}(0^{+})]$
$\displaystyle\times(\nu_{3-k}^{(n)}{}-\nu_{k}^{(n)}{})^{-1},$ $\displaystyle
B_{z}^{(k;n)}=-i(\nu_{k}^{(n)}{}/V^{(n)}{})B_{0}^{(k;n)}\quad(k=1,2).$ (112)

These expressions are used in the main text to obtain
$C_{\lambda\varkappa}^{(n)}(t)$ in the explicit form.

To find the initial conditions for the equations for
$C_{\lambda\varkappa}^{(n)}$ in the above expressions we take into account
that at the instant $t=0^{+}$ the qubit is in the state
$(\ket{0}+\ket{1})/\sqrt{2}$, whereas the stationary TLS populations are given
by Eq. (57). Therefore

$\displaystyle C_{Iz}^{(n)}{}(0^{+})=C_{\pm
z}(0^{+})=-C_{zz}^{(n)}{}(0^{+})=\Delta W^{(n)}/2W^{(n)}\ ,$ $\displaystyle
C_{I0}^{(n)}{}(0^{+})=C_{\pm 0}(0^{+})=-C_{z0}^{(n)}{}(0^{+})=\frac{1}{2},$
(113)

With these initial conditions, we have in particular

$\displaystyle C_{Iz}^{(n)}{}(t_{R})=-C_{zz}^{(n)}{}(t_{R})=\Delta
W^{(n)}{}/2W^{(n)}{}\ ,$ $\displaystyle
C_{I0}^{(n)}{}(t_{R})=-C_{z0}^{(n)}{}(t_{R})=\frac{1}{2}\ ,$ (114)

and

$\displaystyle
C_{+z}^{(n)}{}(t_{R})=\frac{1}{2}e^{-W^{(n)}{}t_{R}/2}\left[\frac{\Delta
W^{(n)}{}}{W^{(n)}{}}\cosh(\gamma^{(n)}t_{R})\right.$
$\displaystyle\left.+\left(i\frac{V^{(n)}{}}{\gamma^{(n)}}+\frac{\Delta
W^{(n)}{}}{2\gamma^{(n)}}\right)\sinh(\gamma^{(n)}t_{R})\right];$ (115)

the expression for $C_{+0}^{(n)}(t_{R})$ is given in the main text, Eq.
(III.1).

##  Appendix F Time evolution of the pair correlator

The operator $\mathfrak{R}_{\mathrm{corr}}$ defined in Eq. (IX.1) describes
the decay of the pair correlations function of Ramsey measurement. The
parameters $\mathbb{K}_{\\{m\\}_{s}}$ in the expression for
$\mathfrak{R}_{\mathrm{corr}}$ read

$\displaystyle\mathbb{K}_{\\{m\\}_{s}}=\frac{1}{2}e^{-t_{R}/T_{2}}\mathrm{Re}\Bigl{[}e^{i\tilde{\phi}_{R}}\mathcal{J}^{(m_{1})}...\mathcal{J}^{(m_{s})}$
$\displaystyle\times\prod_{n\neq
m_{1},...,m_{s}}e^{-W^{(n)}t_{R}/2}\Xi^{(n)}(t_{R})\Bigr{]}$ (116)

where

$\displaystyle\mathcal{J}^{(n)}{}=C_{+z}^{(n)}{}(t_{R})-\frac{\Delta
W^{(n)}{}}{W^{(n)}{}}C_{+0}^{(n)}{}(t_{R}),$
$\displaystyle=2i\frac{V^{(n)}{}}{\gamma^{(n)}}\,\frac{W_{01}^{(n)}{}W_{10}^{(n)}{}}{W^{(n)}{}{}^{2}}e^{-W^{(n)}{}t_{R}/2}\sinh\gamma^{(n)}t_{R}$
(117)

We use here the Hamiltonian $H_{R}$ in which $\phi_{R}$ was replaced by the
auxiliary phase $\tilde{\phi}_{R}$ to compensate the coupling-induced shift of
the average qubit frequency.

###  F.1 Time evolution of the terms $\propto\sigma_{\pm}$ in
$\mathfrak{R}_{\mathrm{corr}}$

As indicated in the main text, the evolution of the terms
$\propto\sigma_{\pm}$ in $\mathfrak{R}_{\mathrm{corr}}(t)$ in the time
interval $kt_{\mathrm{cyc}}<t<kt_{\mathrm{cyc}}+t_{R}$ have to be studied
separately for each term in the sum over $m_{1},...,m_{s}$ in
$\mathfrak{R}_{\mathrm{corr}}(t)$. We denote the corresponding terms as
$\mathfrak{R}_{\mathrm{corr}}(t|\\{m\\}_{s})$,

$\mathfrak{R}_{\mathrm{corr}}(t)=\sum_{s\geq
1}\sum_{\\{m\\}_{s}}\mathfrak{R}_{\mathrm{corr}}(t|\\{m\\}_{s}).$

Similar to Eqs. (VIII.2) and (60), we seek
$\mathfrak{R}_{\mathrm{corr}}(t|\\{m\\}_{s})$ in the form

$\displaystyle\mathfrak{R}_{\mathrm{corr}}(t|\\{m\\}_{s})=\frac{1}{4}e^{-(t-kt_{\mathrm{cyc}})/T_{2}}\mathbb{K}_{\\{m\\}_{s}}\sum_{\alpha=\pm}\sigma_{\alpha}$
$\displaystyle\times\prod_{n}\sum_{\varkappa=0,z}\tilde{C}_{\alpha\varkappa}^{(n)}(t|\\{m\\}_{s})\hat{\tau}_{\varkappa}^{(n)}\
.$ (118)

The equations for the coefficients
$\tilde{C}_{\alpha\varkappa}^{(n)}{}(t|\\{m\\}_{s})$ have the same form as
Eqs. (E) for $C_{\alpha\varkappa}^{(n)}(t)$. The initial conditions for these
equations are set at $t=t_{\mathrm{cyc}}^{+}$. They follow from the expression
(IX.2.2) for $\mathfrak{R}_{\mathrm{corr}}(kt_{\mathrm{cyc}})$ and Eq. (68).
They are different for the values of $n$ that coincide with one of the
components $m_{i}$ of the vector $\\{m\\}_{s}$ and for those values of $n$
that differ from the components of $\\{m\\}_{s}$. For the first group we have

$\displaystyle\tilde{C}_{+0}^{(m_{i})}(kt_{\mathrm{cyc}}^{+}|\\{m\\}_{s})=0,\qquad\tilde{C}_{+z}^{(m_{i})}(kt_{\mathrm{cyc}}^{+}|\\{m\\}_{s})$
$\displaystyle=\exp[-W^{(m_{i})}(kt_{\mathrm{cyc}}-t_{R})],\quad
m_{i}\in\\{m\\}_{s}\ .$ (119)

In contrast, for the second group we have

$\displaystyle\tilde{C}_{+0}^{(n)}{}(kt_{\mathrm{cyc}}^{+}|\\{m\\}_{s})=1/2,\quad\tilde{C}_{+z}^{(n)}{}(kt_{\mathrm{cyc}}^{+}|\\{m\\}_{s})$
$\displaystyle=\Delta W^{(n)}{}/2W^{(n)}{}\quad n\neq m_{1},...,m_{s}\ .$
(120)

It then follows from Eqs. (E) - (E) that

$\displaystyle\tilde{C}_{+0}^{(m_{j})}(kt_{\mathrm{cyc}}+t_{R}|\\{m\\}_{s})=i(V_{m_{j}}/\mu^{(m_{j})})\sinh\mu^{(m_{j})}t_{R}$
$\displaystyle\times\exp\left[-W^{(m_{j})}(kt_{\mathrm{cyc}}-t_{R}/2)\right],\quad
m_{j}\in\\{m\\}_{s}\ ,$ (121)

whereas

$\displaystyle\tilde{C}_{+0}^{(n)}{}(kt_{\mathrm{cyc}}+t_{R}|\\{m\\}_{s})=\frac{1}{2}\Xi^{(n)}(t_{R})\
,$ $\displaystyle n\notin\\{m\\}_{s}\ ,$ (122)

where $\Xi^{(n)}(t_{R})$ is given by Eq. (III.1).

At time $kt_{\mathrm{cyc}}+t_{R}$ the qubit undergoes a rotation about the
$z$-axis by the angle $\tilde{\phi}_{R}$ followed by a $\pi/2$ rotation about
the $y$-axis, as seen from Eqs. (VIII.1) and (53). As a result the matrices
$\sigma_{\alpha}$ in the operator
$\mathfrak{R}_{\mathrm{corr}}(kt_{\mathrm{cyc}}+t_{R}|\\{m\\}_{s})$, Eq.
(F.1), transform as

$\sigma_{\alpha}\to e^{i\alpha\tilde{\phi}_{R}}(i\alpha\sigma_{y}-\sigma_{z})\
.$

This transformation corresponds to the transformation of
$\mathfrak{R}_{\mathrm{corr}}(kt_{\mathrm{cyc}}+t_{R}|\\{m\\}_{s})$ into
$\mathfrak{R}_{\mathrm{corr}}(kt_{\mathrm{cyc}}+t_{R}^{+}|\\{m\\}_{s})$. As
indicated in the main text, if the qubit is reset at
$kt_{\mathrm{cyc}}+t_{R}^{+}$, then
$\mathfrak{R}_{\mathrm{corr}}(kt_{\mathrm{cyc}}+t_{R}^{+}|\\{m\\}_{s})$
decays. However, if the measurement is performed at time
$kt_{\mathrm{cyc}}+t_{R}^{+}$, i.e.,
$\mathfrak{R}_{\mathrm{corr}}(kt_{\mathrm{cyc}}+t_{R}^{+}|\\{m\\}_{s})$ is
multiplied by $\hat{\pi}=(\hat{I}_{q}-\sigma_{z})/2$ and the trace is taken
over the qubit and the TLSs, the contribution of
$\mathfrak{R}_{\mathrm{corr}}(kt_{\mathrm{cyc}}+t_{R}^{+}|\\{m\\}_{s})$ to the
result of the measurement is

$\displaystyle\tilde{r}_{2}(k)=\frac{1}{2}e^{-t_{R}/T_{2}}\sum_{s}\sum_{\\{m\\}_{s}}\mathbb{K}_{\\{m\\}_{s}}$
$\displaystyle\times\mathrm{Re}\left[e^{i\tilde{\phi}_{R}}\prod_{n}2\tilde{C}_{+0}^{(n)}(kt_{\mathrm{cyc}}+t_{R}|\\{m\\}_{s})\right]\
.$ (123)

This result depends only on the coefficients $\tilde{C}_{+\varkappa}^{(n)}$
with $\varkappa=0$, which describe the contribution of the unit operators of
the TLSs; the trace of the terms proportional to $\hat{\tau}_{z}^{(n)}$ is
zero.

##  Appendix G Explicit expressions for the phase correlators for $1/f$ noise

For completeness, here we provide the explicit expressions for the correlators
$f_{k}=\langle\theta_{0}\theta_{k}\rangle$ of the phases accumulated by the
qubit during Ramsey measurements separated by $k$ cycles, i.e., separated by
the time $kt_{\mathrm{cyc}}$. The expressions are related to Gaussian noise
$\delta\omega_{\mathrm{q}}(t)$ produced by a large number of TLSs with the
same coupling to the qubit and with the log-normal distribution of the
switching rates. The power spectrum $S_{q}(\omega)$ of
$\delta\omega_{\mathrm{q}}(t)$ in this case is given by Eq. (IV.4).

The integral over $\omega$, Eq. (IV.1), that relates $f_{k}$ to
$S_{q}(\omega)$, can be expressed in terms of the exponential integral and
hyperbolic sine and cosine integral functions
$\mathrm{Ei}(z),\mathrm{shi}(z)$, and $\mathrm{chi}(z)$ as follows

$\displaystyle f_{k}=\frac{D_{\rm
fl}t_{R}^{2}}{2\pi}\\{2e^{-a_{k}b_{\min}}[(1-a_{k}b_{\min})[\cosh
b_{\min}-1]+b_{\min}\sinh b_{\min}]/b_{\min}^{2}$
$\displaystyle+2a_{k}^{2}\mathrm{Ei}(-a_{k}b_{\min})-(a_{k}+1)^{2}\mathrm{Ei}(-a_{k}b_{\min}-b_{\min})-(a_{k}-1)^{2}\mathrm{Ei}(b_{\min}-a_{k}b_{\min})\\},\quad
k>0\ ,$ (124)

and

$\displaystyle f_{0}=\frac{D_{\rm
fl}t_{R}^{2}}{\pi}[-\frac{1}{b_{\min}^{2}}+\frac{2}{b_{\min}}-(b_{\min}-1)\exp(-b_{\min})/b_{\min}^{2}-\mathrm{chi}(b_{\min})+\mathrm{shi}(b_{\min})]\
.$ (125)

These expressions depend on two dimensionless parameters

$a_{k}=\frac{kt_{\mathrm{cyc}}}{t_{R}}\ ,\qquad b_{\min}=\omega_{\min}t_{R}\
,$

that are convenient for a numerical evaluation of the correlators $f_{k}$.

## References

  * Álvarez and Suter [2011] G. A. Álvarez and D. Suter, Measuring the Spectrum of Colored Noise by Dynamical Decoupling, Phys. Rev. Lett. 107, 230501 (2011).
  * Bylander  _et al._ [2011] J. Bylander, S. Gustavsson, F. Yan, F. Yoshihara, K. Harrabi, G. Fitch, D. G. Cory, Y. Nakamura, J.-S. Tsai, and W. D. Oliver, Noise spectroscopy through dynamical decoupling with a superconducting flux qubit, Nat. Phys. 7, 565 (2011).
  * Sank  _et al._ [2012] D. Sank, R. Barends, R. C. Bialczak, Y. Chen, J. Kelly, M. Lenander, E. Lucero, M. Mariantoni, A. Megrant, M. Neeley, P. J. J. O’Malley, A. Vainsencher, H. Wang, J. Wenner, T. C. White, T. Yamamoto, Y. Yin, A. N. Cleland, and J. M. Martinis, Flux Noise Probed with Real Time Qubit Tomography in a Josephson Phase Qubit, Phys. Rev. Lett. 109, 067001 (2012).
  * Yan  _et al._ [2012] F. Yan, J. Bylander, S. Gustavsson, F. Yoshihara, K. Harrabi, D. G. Cory, T. P. Orlando, Y. Nakamura, J.-S. Tsai, and W. D. Oliver, Spectroscopy of low-frequency noise and its temperature dependence in a superconducting qubit, Phys. Rev. B 85, 174521 (2012).
  * Paz-Silva and Viola [2014] G. A. Paz-Silva and L. Viola, General Transfer-Function Approach to Noise Filtering in Open-Loop Quantum Control, Phys. Rev. Lett. 113, 250501 (2014).
  * Yoshihara  _et al._ [2014] F. Yoshihara, Y. Nakamura, F. Yan, S. Gustavsson, J. Bylander, W. D. Oliver, and J.-S. Tsai, Flux qubit noise spectroscopy using Rabi oscillations under strong driving conditions, Phys. Rev. B 89, 020503 (2014).
  * Kim  _et al._ [2015] M. Kim, H. J. Mamin, M. H. Sherwood, K. Ohno, D. D. Awschalom, and D. Rugar, Decoherence of Near-Surface Nitrogen-Vacancy Centers Due to Electric Field Noise, Phys. Rev. Lett. 115, 087602 (2015).
  * Brownnutt  _et al._ [2015] M. Brownnutt, M. Kumph, P. Rabl, and R. Blatt, Ion-trap measurements of electric-field noise near surfaces, Rev. Mod. Phys. 87, 1419 (2015).
  * O’Malley  _et al._ [2015] P. J. J. O’Malley, J. Kelly, R. Barends, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. G. Fowler, I.-C. Hoi, E. Jeffrey, A. Megrant, J. Mutus, C. Neill, C. Quintana, P. Roushan, D. Sank, A. Vainsencher, J. Wenner, T. C. White, A. N. Korotkov, A. N. Cleland, and J. M. Martinis, Qubit Metrology of Ultralow Phase Noise Using Randomized Benchmarking, Phys. Rev. Applied 3, 044009 (2015).
  * Szańkowski  _et al._ [2016] P. Szańkowski, M. Trippenbach, and Ł. Cywiński, Spectroscopy of cross correlations of environmental noises with two qubits, Phys. Rev. A 94, 012109 (2016).
  * Yan  _et al._ [2016] F. Yan, S. Gustavsson, A. Kamal, J. Birenbaum, A. P. Sears, D. Hover, T. J. Gudmundsen, D. Rosenberg, G. Samach, S. Weber, J. L. Yoder, T. P. Orlando, J. Clarke, A. J. Kerman, and W. D. Oliver, The flux qubit revisited to enhance coherence and reproducibility, Nat Commun 7, 1 (2016).
  * Myers  _et al._ [2017] B. A. Myers, A. Ariyaratne, and A. C. B. Jayich, Double-Quantum Spin-Relaxation Limits to Coherence of Near-Surface Nitrogen-Vacancy Centers, Phys. Rev. Lett. 118, 197201 (2017).
  * Quintana  _et al._ [2017] C. M. Quintana, Y. Chen, D. Sank, A. G. Petukhov, T. C. White, D. Kafri, B. Chiaro, A. Megrant, R. Barends, B. Campbell, Z. Chen, A. Dunsworth, A. G. Fowler, R. Graff, E. Jeffrey, J. Kelly, E. Lucero, J. Y. Mutus, M. Neeley, C. Neill, P. J. J. O’Malley, P. Roushan, A. Shabani, V. N. Smelyanskiy, A. Vainsencher, J. Wenner, H. Neven, and J. M. Martinis, Observation of Classical-Quantum Crossover of $1/f$ Flux Noise and Its Paramagnetic Temperature Dependence, Phys. Rev. Lett. 118, 057702 (2017).
  * Ferrie  _et al._ [2018] C. Ferrie, C. Granade, G. Paz-Silva, and H. M. Wiseman, Bayesian quantum noise spectroscopy, New J. Phys. 20, 123005 (2018).
  * Noel  _et al._ [2019] C. Noel, M. Berlin-Udi, C. Matthiesen, J. Yu, Y. Zhou, V. Lordi, and H. Häffner, Electric-field noise from thermally activated fluctuators in a surface ion trap, Phys. Rev. A 99, 063427 (2019).
  * von Lüpke  _et al._ [2020] U. von Lüpke, F. Beaudoin, L. M. Norris, Y. Sung, R. Winik, J. Y. Qiu, M. Kjaergaard, D. Kim, J. Yoder, S. Gustavsson, L. Viola, and W. D. Oliver, Two-Qubit Spectroscopy of Spatiotemporally Correlated Quantum Noise in Superconducting Qubits, PRX Quantum 1, 010305 (2020).
  * Wolfowicz  _et al._ [2021] G. Wolfowicz, F. J. Heremans, C. P. Anderson, S. Kanai, H. Seo, A. Gali, G. Galli, and D. D. Awschalom, Qubit guidelines for solid-state spin defects, Nat. Rev. Mater. 6, 906 (2021), comment: 40 pages, 7 figures, 259 references.
  * Wang and Clerk [2021] Y.-X. Wang and A. A. Clerk, Intrinsic and induced quantum quenches for enhancing qubit-based quantum noise spectroscopy, Nat. Commun. 12, 6528 (2021).
  * Ristè  _et al._ [2013] D. Ristè, C. C. Bultink, M. J. Tiggelman, R. N. Schouten, K. W. Lehnert, and L. DiCarlo, Millisecond charge-parity fluctuations and induced decoherence in a superconducting transmon qubit, Nat Commun 4, 1913 (2013).
  * Serniak  _et al._ [2018] K. Serniak, M. Hays, G. de Lange, S. Diamond, S. Shankar, L. D. Burkhart, L. Frunzio, M. Houzet, and M. H. Devoret, Hot Nonequilibrium Quasiparticles in Transmon Qubits, Phys. Rev. Lett. 121, 157701 (2018).
  * Christensen  _et al._ [2019] B. G. Christensen, C. D. Wilen, A. Opremcak, J. Nelson, F. Schlenker, C. H. Zimonick, L. Faoro, L. B. Ioffe, Y. J. Rosen, J. L. DuBois, B. L. T. Plourde, and R. McDermott, Anomalous Charge Noise in Superconducting Qubits, Phys. Rev. B 100, 140503 (2019), comment: 10 pages, 7 figures, arXiv:1905.13712 .
  * Schlör  _et al._ [2019] S. Schlör, J. Lisenfeld, C. Müller, A. Bilmes, A. Schneider, D. P. Pappas, A. V. Ustinov, and M. Weides, Correlating Decoherence in Transmon Qubits: Low Frequency Noise by Single Fluctuators, Phys. Rev. Lett. 123, 190502 (2019).
  * Paladino  _et al._ [2002] E. Paladino, L. Faoro, G. Falci, and R. Fazio, Decoherence and $1/f$ Noise in Josephson Qubits, Phys. Rev. Lett. 88, 228304 (2002).
  * Galperin  _et al._ [2004] Y. M. Galperin, B. L. Altshuler, and D. V. Shantsev, Low-frequency noise as a source of dephasing of a qubit, in  _Fundamental Problems of Mesoscopic Physics_ , edited by I. V. Lerner (Kluwer Academic Publishing, The Netherlands, 2004) pp. 141–165, comment: 18 pages, 8 figures, Proc. of NATO/Euresco Conf. ”Fundamental Problems of Mesoscopic Physics: Interactions and Decoherence”, Granada, Spain, Sept.2003.
  * Galperin  _et al._ [2006] Y. M. Galperin, B. L. Altshuler, J. Bergli, and D. V. Shantsev, Non-Gaussian Low-Frequency Noise as a Source of Qubit Decoherence, Phys. Rev. Lett. 96, 097009 (2006).
  * Paladino  _et al._ [2014] E. Paladino, Y. M. Galperin, G. Falci, and B. L. Altshuler, $1/f$ noise: Implications for solid-state quantum information, Rev. Mod. Phys. 86, 361 (2014).
  * Müller  _et al._ [2019] C. Müller, J. H. Cole, and J. Lisenfeld, Towards understanding two-level-systems in amorphous solids: Insights from quantum circuits, Rep. Prog. Phys. 82, 124501 (2019).
  * Herzog and Hahn [1956] B. Herzog and E. L. Hahn, Transient Nuclear Induction and Double Nuclear Resonance in Solids, Phys. Rev. 103, 148 (1956).
  * Klauder and Anderson [1962] J. R. Klauder and P. W. Anderson, Spectral Diffusion Decay in Spin Resonance Experiments, Phys. Rev. 125, 912 (1962).
  * Li  _et al._ [2013] F. Li, A. Saxena, D. Smith, and N. A. Sinitsyn, Higher-order spin noise statistics, New J. Phys. 15, 113038 (2013).
  * Norris  _et al._ [2016] L. M. Norris, G. A. Paz-Silva, and L. Viola, Qubit Noise Spectroscopy for Non-Gaussian Dephasing Environments, Phys. Rev. Lett. 116, 150503 (2016).
  * Szańkowski  _et al._ [2017] P. Szańkowski, G. Ramon, J. Krzywda, D. Kwiatkowski, and L. Cywiński, Environmental noise spectroscopy with qubits subjected to dynamical decoupling, J. Phys.: Condens. Matter 29, 333001 (2017).
  * Sung  _et al._ [2019] Y. Sung, F. Beaudoin, L. M. Norris, F. Yan, D. K. Kim, J. Y. Qiu, U. von Lüpke, J. L. Yoder, T. P. Orlando, S. Gustavsson, L. Viola, and W. D. Oliver, Non-Gaussian noise spectroscopy with a superconducting qubit sensor, Nat. Commun. 10, 3715 (2019).
  * Maizelis  _et al._ [2011] Z. A. Maizelis, M. L. Roukes, and M. I. Dykman, Detecting and Characterizing Frequency Fluctuations of Vibrational Modes, Phys. Rev. B 84, 144301 (2011).
  * Sun  _et al._ [2015] F. Sun, J. Zou, Z. A. Maizelis, and H. B. Chan, Telegraph Frequency Noise in Electromechanical Resonators, Phys. Rev. B 91, 174102 (2015).
  * Wilen  _et al._ [2021] C. D. Wilen, S. Abdullah, N. A. Kurinsky, C. Stanford, L. Cardani, G. D’Imperio, C. Tomei, L. Faoro, L. B. Ioffe, C. H. Liu, A. Opremcak, B. G. Christensen, J. L. DuBois, and R. McDermott, Correlated charge noise and relaxation errors in superconducting qubits, Nature 594, 369 (2021).
  * Liu  _et al._ [2022] C.-H. Liu, D. C. Harrison, S. Patel, C. D. Wilen, O. Rafferty, A. Shearrow, A. Ballard, V. Iaia, J. Ku, B. L. T. Plourde, and R. McDermott, Quasiparticle Poisoning of Superconducting Qubits from Resonant Absorption of Pair-breaking Photons, ArXiv220306577 Quant-Ph (2022), comment: 13 pages, 11 figures, 3 tables, arXiv:2203.06577 [quant-ph] .
  * Anderson  _et al._ [1972] P. W. Anderson, B. I. Halperin, and C. M. Varma, Anomalous low-temperature thermal properties of glasses and spin glasses, Philos. Mag. J. Theor. Exp. Appl. Phys. 25, 1 (1972).
  * Phillips [1972] W. A. Phillips, Tunneling states in amorphous solids, J Low Temp Phys 7, 351 (1972).
  * Phillips [1987] W. A. Phillips, Two-Level States in Glasses, Rep. Prog. Phys. 50, 1657 (1987).
  * Nielsen and Chuang [2011] M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information: 10th Anniversary Edition_ , 1st ed. (Cambridge University Press, Cambridge ; New York, 2011).
  * You  _et al._ [2021] X. You, A. A. Clerk, and J. Koch, Positive- and negative-frequency noise from an ensemble of two-level fluctuators, Phys Rev Res. 3, 013045 (2021), comment: 15 pages, 8 figures, arXiv:2005.03591 .
  * Ithier  _et al._ [2005] G. Ithier, E. Collin, P. Joyez, P. J. Meeson, D. Vion, D. Esteve, F. Chiarello, A. Shnirman, Y. Makhlin, J. Schriefl, and G. Schön, Decoherence in a superconducting quantum bit circuit, Phys. Rev. B 72, 134519 (2005).
  * Gurvitz [2019] S. Gurvitz, Generalized Landauer formula for time-dependent potentials and noise-induced zero-bias dc current, J. Phys. A: Math. Theor. 52, 175301 (2019).
  * Shapiro and Loginov [1978] V. E. Shapiro and V. M. Loginov, “Formulae of differentiation” and their use for solving stochastic equations, Physica A: Statistical Mechanics and its Applications 91, 563 (1978).


# Technical Report for Argoverse Challenges on 4D Occupancy Forecasting

Pengfei Zheng 1, Kanokphan Lertniphonphan11footnotemark: 1 3, Feng Chen
11footnotemark: 1 3, Siwei Chen ††footnotemark:  2,
Bingchuan Sun3, Jun Xie3, Zhepeng Wang3
1University of Science and Technology Beijing, 2Tsinghua University
3Lenovo Research
<EMAIL_ADDRESS>Equal contributionWork done as an intern at Lenovo
Research.

###### Abstract

This report presents our Le3DE2E_Occ solution for 4D Occupancy Forecasting in
Argoverse Challenges at CVPR 2023 Workshop on Autonomous Driving (WAD). Our
solution consists of a strong LiDAR-based Bird’s Eye View (BEV) encoder with
temporal fusion and a two-stage decoder, which combines a DETR head and a UNet
decoder. The solution was tested on the Argoverse 2 sensor dataset to evaluate
the occupancy state 3 seconds in the future. Our solution achieved 18% lower
L1 Error (3.57) than the baseline on the 4D Occupancy Forecasting task in
Argoverse Challenges at CVPR 2023.

##  1 Introduction

4D occupancy forecasting is a new challenge introduced in CVPR 2023 WAD. The
task is to understand how an environment evolves with time which is crucial
for motion planning in autonomous driving. However, conventional methods
require costly human annotations, such as detection bounding boxes, tracking
ids, or semantic segmentation, which make it difficult to scale up to a large
labeled dataset. This challenge aims to learn future occupancy forecasting
from unlabeled datasets.

In this challenge, given a particular agent’s observation of the world in the
past n seconds, we need to predict the space-time evolution of the world in
the future n seconds. Specifically, the occupancy state of 5 frames in the
next 3s is predicted by observing the point cloud of 5 frames in the past and
present timestamp within 3s (at a frequency of 5/3Hz). The future frame point
cloud is obtained by rendering the occupancy state from a given query ray.
Then, all point clouds and occupancies are aligned to the current frame under
the LIDAR coordinate system.

Our solution adopts a voxel feature encoder to transform LiDAR point clouds
into a 2D Bird’s Eye View (BEV) feature map and uses a DETR[1] structure head
to predict voxel-wise future occupancy. To refine the results, a UNet[4] head
is employed as a second-stage decoder to produce more accurate forecasting
results.

Figure 1: System overview. Firstly the LiDAR point clouds of the current frame are voxelized and encoded to the BEV feature map. Secondly, the historical frame BEV feature maps are fused with the current frame by a temporal encoder. Thirdly, Spatial-temporal fused BEV is fed into OccFormer which generates voxel-wise occupancy forecasting results. Fourthly, UNet works as a second-stage decoder to refine the forecasts. Last, voxel rendering generates point-wise depth along the given ray direction as post-processing.  Team | L1 ($\downarrow$) | AbsRel | NFCD | CD
---|---|---|---|---
Host_34336_Team (Raytracing) | 6.68 | 0.48 | 8.79 | 49.88
Host_34336_Team (Point Cloud Forecasting as a P) | 4.34 | 0.26 | 5.23 | 92.08
Le3DE2E_Occ (Ours) | 3.57 | 0.22 | 3.35 | 91.61
Table 1: 4D Occupancy Forecasting on the official leaderboard

##  2 Method

We employ the bev feature as a unified representation and adopted the
OccFormer module as an occupancy head for 4d occupancy forecasting based on
UniAD[2]. OccFormer module consists of a transformer decoder with T sequential
blocks. The pipeline of our method is shown in Figure 1.

###  2.1 BEV feature

We employ the LIDAR BEV encoder based on SECOND[6] to generate the LIDAR BEV
feature $B_{l}$. The LIDAR encoder takes as input the current and past $T=5$
frames of the point cloud with a time step of 0.6 s. All point clouds have
been aligned to the current frame. Afterward, we fuse the past T frames with
the current frame to aggregate temporal features by a 2D convolutional block.
The spatial-temporal fused BEV feature map is fed to the occupancy decoder.

###  2.2 OccFormer

Given the BEV feature from upstream modules, we feed it to OccFormer to get
the multi-time-step occupancy results. OccFormer consists of $T$ sequential
blocks, where $T=5$ is a number of future frames. Each block is responsible
for generating the occupancy of a particular frame. Unlike the instance
occupancy output by OccFormer in UniAD, we need dense voxel-wise occupancy.

In each sequential block, the BEV features first perform self-queries through
the self-attentive layer to compute the similarity metric within the feature.
Then, we randomly initialize the instance query $Q_{I}$, which are track query
$Q_{A}$ , agent position $P_{A}$ and motion query $Q_{X}$. The instance query
$Q_{I}$ and BEV feature interact using cross-attention so that the dense scene
feature and sparse agent feature benefit from each other. The interacted BEV
feature is sent to the next sequential block, which is cycled in RNN fashion.

All $T$ frame BEV features are upsampled to obtain future occupancy $O_{T}$.
Dimensions of all dense features and instance features are 256 in OccFormer. A
UNet head is used as a second-stage decoder to enhance the multiscale
forecasting results. After generating voxel-wise occupancy forecasting, voxel
rendering is then performed by the query rays to get the point-wise estimated
depth.

##  3 Experiments

###  3.1 Dataset

The competition used the Argoverse 2 Sensor Dataset[5], which consisted of
1000 scenes (750 for training, 150 for validation, and 150 for testing) with a
total of 4.2 hours of driving data. The total dataset is extracted in the form
of 1 TB of data. Each vehicle log has a duration of approximately 15 seconds
and includes an average of approximately 150 LiDAR scans with 10 FPS LiDAR
frames.

###  3.2 Evaluation Metrics

The metric is performed in the range of [-70m, -70m, -4.5m, 70m, 70m, 4.5m]
with a voxel resolution of 0.2m. The absolute L1 distance between the true
expected depth along a given ray direction and the predicted expected depth
obtained by rendering along the same ray direction is used as the main metric.
Absolute relative error (AbsRel), near-field chamfer distance (NFCD), and
vanilla chamfer distance (CD) are measured together as other metrics. More
details of the evaluation can be found in the paper [3].

###  3.3 Implementation details

We followed the baseline [3] for data preparation for training and evaluation.
Point shuffle is used for data augmentation in the training stage. The voxel
size is set to (0.075m, 0.075m, 0.2m) and the resulting BEV feature $B_{l}$
has a shape of 240*240. We train the model with a cosine annealing policy with
a 0.001 learning rate and use L1 loss and Adam optimizer. We train the model
from scratch for 20 epochs on 8 V100 GPUs with a total batch size of 8.

We submitted our results to the testing server and got a 3.57 L1 Error. The
results are shown in Table 1. Our results outperform the baseline with $18\%$
and $15\%$ improvements in L1 Error and Absolute Relative L1 Error,
respectively.

##  4 Conclusion

In our model, we employ a LIDAR encoder to encode spatial features and then do
a temporal fusion of historical BEV features. We use the BEV feature as a
unified intermediate representation. We employ an OccFormer as 4D occupancy
prediction head to loop out future occupancy in RNN style. Following the work
of [3], we ”render” point cloud data from 4D occupancy predictions and
estimate depth for supervision. The experimental results indicate that our
model achieved better scores than the baseline with $18\%$ and $15\%$
improvements in L1 Error and Absolute Relative L1 Error, respectively.

## References

  * [1] Nicolas Carion, Francisco Massa, Gabriel Synnaeve, Nicolas Usunier, Alexander Kirillov, and Sergey Zagoruyko.  End-to-end object detection with transformers.  In Andrea Vedaldi, Horst Bischof, Thomas Brox, and Jan-Michael Frahm, editors, Computer Vision – ECCV 2020, pages 213–229, Cham, 2020. Springer International Publishing.
  * [2] Yihan Hu, Jiazhi Yang, Li Chen, Keyu Li, Chonghao Sima, Xizhou Zhu, Siqi Chai, Senyao Du, Tianwei Lin, Wenhai Wang, Lewei Lu, Xiaosong Jia, Qiang Liu, Jifeng Dai, Yu Qiao, and Hongyang Li.  Planning-oriented autonomous driving.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2023.
  * [3] Tarasha Khurana, Peiyun Hu, David Held, and Deva Ramanan.  Point cloud forecasting as a proxy for 4d occupancy forecasting.  In IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2023.
  * [4] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  ArXiv, abs/1505.04597, 2015.
  * [5] Benjamin Wilson, William Qi, Tanmay Agarwal, John Lambert, Jagjeet Singh, Siddhesh Khandelwal, Bowen Pan, Ratnesh Kumar, Andrew Hartnett, Jhony Kaesemodel Pontes, Deva Ramanan, Peter Carr, and James Hays.  Argoverse 2: Next generation datasets for self-driving perception and forecasting.  In Proceedings of the Neural Information Processing Systems Track on Datasets and Benchmarks (NeurIPS Datasets and Benchmarks 2021), 2021\.
  * [6] Yan Yan, Yuxing Mao, and Bo Li.  Second: Sparsely embedded convolutional detection.  Sensors, 18(10):3337, 2018.