††thanks: Corresponding Author, Email<EMAIL_ADDRESS>
# New Spin Structure Constraints on Hyperfine Splitting and Proton Size

D. Ruth  University of New Hampshire, Durham, New Hampshire 03824, USA  K.
Slifer  University of New Hampshire, Durham, New Hampshire 03824, USA  J.-P.
Chen  Thomas Jefferson National Accelerator Facility, Newport News, Virginia
23606, USA  C. E. Carlson  William & Mary, Williamsburg, Virginia 23187, USA
F. Hagelstein  Institute of Nuclear Physics, Johannes Gutenberg
Universit$\ddot{a}$t Mainz, 55099 Mainz, Germany  PRISMA+ Cluster of
Excellence, Johannes Gutenberg Universit$\ddot{a}$t Mainz, 55099 Mainz,
Germany  Laboratory for Particle Physics, Paul Scherrer Institute, 5232
Villigen PSI, Switzerland  V. Pascalutsa  Institute of Nuclear Physics,
Johannes Gutenberg Universit$\ddot{a}$t Mainz, 55099 Mainz, Germany  A. Deur
Thomas Jefferson National Accelerator Facility, Newport News, Virginia 23606,
USA  S. Kuhn  Old Dominion University, Norfolk, Virginia 23529, USA  M. Ripani
Istituto Nazionale di Fisica Nucleare, Sezione di Genova, I-16146 Genova,
Italy  X. Zheng  University of Virginia, Charlottesville, Virginia 22903, USA
R. Zielinski  University of New Hampshire, Durham, New Hampshire 03824, USA
C. Gu  Duke University, Durham, NC 27708, USA

###### Abstract

The 1S hyperfine splitting in hydrogen is measured to an impressive ppt
precision and will soon be measured to ppm precision in muonic hydrogen. The
latter measurement will rely on theoretical predictions, which are limited by
knowledge of the proton polarizability effect $\Delta_{\text{pol}}$. Data-
driven evaluations of $\Delta_{\text{pol}}$ have long been in significant
tension with baryon chiral perturbation theory. Here we present improved
results for $\Delta_{\text{pol}}$ driven by new spin structure data, reducing
the long-standing tension between theory and experiment and halving the
dominating uncertainty in hyperfine splitting calculations.

###### pacs:

11.55.Hx,25.30.Bf,29.25.Pj,29.27.Hj

The hyperfine splitting in hydrogen, the renowned 21 cm line arising from the
magnetic dipole interaction of electron and proton, stands out as one of the
best-measured quantities in physics, currently known to an impressive 12
digits. This level of precision is challenging for theory to match,
particularly in accounting for the effects of proton structure [1]. These
structure effects are amplified in muonic hydrogen, where the hydrogen’s
electron is replaced by a muon. Presently, several high-profile experiments
[2, 3] are aiming at a first-ever measurement of the ground-state hyperfine
splitting in muonic hydrogen. Their success in finding this forbiddingly
narrow transition crucially depends on an accurate assessment of proton
structure effects. However, previous determinations [4, 5, 6] of the leading
uncertainty among these proton structure effects, the proton polarizability
effect $\Delta_{\text{pol}}$, have large error bars and are in significant
tension with corresponding theoretical calculations [4, 7, 8]. In this letter,
we improve upon the evaluation of the proton polarizability contribution with
new experimental proton spin structure data from Thomas Jefferson National
Laboratory, significantly reducing this long-standing tension and halving this
quantity’s uncertainty.

The classical picture of the hydrogen atom is fairly simple: a pointlike
spinless electron bound by a pointlike spinless proton via the Coulomb force.
The Schrödinger equation gives the energy spectrum in natural units as
$E_{n}=-\alpha^{2}m_{r}/(2n^{2})$, with $n$ the principal quantum number and
$m_{r}$ the reduced mass. The more sophisticated picture of the atom is
modeled as a correction to this simple picture, including the effects due to
spin and the structure of the proton, to be discussed here. The proton
structure effects are small, but clearly seen in the hydrogen spectrum at the
current level of precision. They are more prominent in muonic hydrogen,
because of a much smaller Bohr radius, $a_{B}=1/(\alpha m_{r})$, given that
$m_{r}$ goes roughly as the lepton mass, which is 200 times heavier for the
muon. The muon has thus 2003 greater probability [as given by the wave-
function squared at the origin, $\left|\Psi_{n}(0)\right|^{2}=1/(\pi
a_{B}^{3}n^{3})$] to be probing the proton substructure.

Because of this heightened sensitivity to nuclear structure details in muonic
atoms, the recent breakthrough in the laser spectroscopy of muonic hydrogen
($\mu$H) by the CREMA Collaboration led to an order-of-magnitude improvement
in the measurement of the proton charge radius [9]. Surprisingly to many, it
appeared to be $7\,\sigma$ smaller than the most-recent CODATA recommended
value of the charge radius at the time [10]. This spectacular discrepancy,
dubbed as the “proton radius puzzle” (see [11, 12] for an early review) has
since been largely resolved and, from 2018 onwards, CODATA recommends the
smaller, and more precise, $\mu$H value. This chapter is not yet closed, with
many new measurements of the proton charge radius underway using the
conventional methods of normal hydrogen (H) spectroscopy, elastic electron-
proton ($ep$) scattering and even muon-proton ($\mu p$) scattering, see [7,
13, 14] for recent reviews. Here we concern ourselves with the calculation of
the $\mu$H hyperfine splitting (HFS), and its implications for the next
milestone of $\mu$H spectroscopy: the upcoming experimental measurement of the
$\mu$H ground-state HFS.

Two different collaborations are competing to provide this first-ever HFS
measurement: CREMA [2] at the Paul Scherrer Institut (PSI) and FAMU [3] at the
RIKEN-RAL Muon Facility. Given the extreme narrowness of this transition,
their success depends in part on how well the proton-structure corrections are
understood, since the searches can only be done over a very limited range of
frequencies.

To leading order, $\mathcal{O}(\alpha^{4})$, the HFS is given by the Fermi
energy

$E_{\mathrm{F}}=\frac{8\alpha}{3a_{B}^{3}}\frac{1+\kappa_{p}}{m_{l}M_{p}},$
(1)

where $m_{l}$ is the lepton mass (either $m_{e}$ in H or $m_{\mu}$ in $\mu$H),
and $M_{p}$ is the proton mass. Here the proton structure is only represented
through the anomalous magnetic moment.

At the next order, $\mathcal{O}(\alpha^{5})$, the proton structure effects can
all be computed via the two-photon (2$\gamma$) exchange diagram of Fig. 1,
which usually is split into three contributions:

$\displaystyle
E_{nS\text{-HFS}}^{2\gamma}=\frac{E_{\mathrm{F}}}{n^{3}}\big{(}\Delta_{\mathrm{Z}}+\Delta_{\mathrm{recoil}}+\Delta_{\mathrm{pol}}\big{)}$
(2)

The largest, $\Delta_{Z}$, comes from the proton Zemach radius $R_{Z}$, a
measure of how far the electric and magnetic distributions of the proton are
correlated with each other. $R_{Z}$ is expressed in terms of the elastic
electric and magnetic form factors $G_{E}(Q^{2})$ and $G_{M}(Q^{2})$:

$R_{\mathrm{Z}}=-\frac{4}{\pi}\int_{0}^{\infty}\frac{dQ}{Q^{2}}\left[\frac{G_{E}(Q^{2})G_{M}(Q^{2})}{g_{p}}-1\right],$
(3)

with $g_{p}$ as the proton gyromagnetic g-factor. The recoil contribution
$\Delta_{\mathrm{recoil}}$ can likewise be expressed in terms of the form
factors, but the final contribution from the polarizability,
$\Delta_{\mathrm{pol}}$, is more complicated.

The polarizability effect is caused by the proton’s moments induced by the
electromagnetic fields of the bound leptons. Unlike the radius, the
polarizability contribution is not given by the form factors, but rather by
the inelastic structure functions $g_{1,2}(x,Q^{2})$, which are functions of
the Bjorken $x$, a variable which tracks the fraction of the interaction’s
momentum carried by one of the proton’s quark constituents. This contribution
is more difficult to obtain, due to the necessity to cover a 2-dimensional
phase space, while the required spin structure function data are relatively
sparse, especially at low $Q^{2}$ which dominates the determination.
Previously, there had been only limited $g_{1}$ data and a complete lack of
$g_{2}$ data in the kinematic region most relevant to the HFS. Nonetheless, a
data-driven evaluation of this contribution has been attempted in the past [4,
5, 6].

The present status is that the existing data-driven evaluations, while
consistent with each other, are in disagreement with chiral perturbation
theory ($\chi$PT), which predicts a significantly smaller contribution of this
effect [15, 8].

Under the general assumptions of unitarity (optical theorem) and analyticity
(dispersion relations) of the forward Compton scattering, the contribution of
the spin structure functions has the following form [4, 7]:

$\displaystyle\Delta_{\mathrm{pol}}$ $\displaystyle=$
$\displaystyle\frac{\alpha
m_{l}}{2\pi(1+\kappa_{p})M_{p}}\big{(}\Delta_{1}+\Delta_{2}\big{)},$ (4)
$\displaystyle\Delta_{1}$ $\displaystyle=$
$\displaystyle\int\limits_{0}^{\infty}\frac{dQ^{2}}{Q^{2}}\bigg{[}\beta_{1}(\tau_{l})F_{2}^{2}(Q^{2})$
(5) $\displaystyle+$
$\displaystyle\frac{8M_{p}^{2}}{Q^{2}}\int\limits_{0}^{x_{\mathrm{th}}}dx\,\tilde{\beta}_{1}(\tau,\tau_{l})\,g_{1}(x,Q^{2})\bigg{]},$
$\displaystyle\Delta_{2}$ $\displaystyle=$
$\displaystyle-24M_{p}^{2}\int\limits_{0}^{\infty}\frac{dQ^{2}}{Q^{4}}\int\limits_{0}^{x_{\mathrm{th}}}dx\,\tilde{\beta}_{2}(\tau,\tau_{l})\,g_{2}(x,Q^{2})$
(6)

where $\beta$’s are elementary kinematic functions, $\tau$’s are kinematic
variables, and $F_{2}$ is the Pauli form factor, the explicit definitions of
which are in the Supplemental Materials. $x_{\mathrm{th}}$ corresponds to the
minimum energy necessary to generate a pion, at an invariant mass $W$ of
1073.2 MeV.

tpe (150,60) i2,K,i1 o2,Kbar,o1
fermion,label.side=left,label=$e^{-}$/$\mu^{-}$i1,t1
fermion,label.side=left,label=pi2,t4 10 leftK input vovalblob; vovalblobz
((.5w), (.4), K); l.d=0,l.a=180,l=H/$\mu$HK phantom,tension=1t1,t2
phantom,tension=0.5t4,t3 phantom,tension=1t4,ti,t3
fermion,tension=1,label.side=left,label=$e^{-}$/$\mu^{-}$t2,o1
fermion,tension=1,label.side=left,label=pt3,o2 10 rightKbar input vovalblob;
vovalblobz ((.5w), (.4), Kbar); l.d=0,l.a=3,l=H/$\mu$HKbar
fermion,tension=0,label.side=leftt1,t2 fermion,tension=0,label.side=leftt4,t3
input vovalblob; vovalblob ((.2w), (1.8), ti);
boson,tension=0,label=$\gamma$,label.side=rightt1,t4
boson,tension=0,label=$\gamma$t3,t2

Figure 1: Two-photon exchange diagram showing the interaction between electron
($e^{-}$) and proton (p) in hydrogen, or muon ($\mu^{-}$) and proton in muonic
hydrogen. Figure 2: The hyperfine contribution integrands for $\Delta_{1}$ and
$\Delta_{2}$ in Eqs. (5) and (6), weighted by $Q^{2}$, for muonic hydrogen.
Results from the g2p experiment [16] are shown in blue squares. The results of
the EG4 experiment [17] are shown in orange triangles. The inner error bars
represent the statistical uncertainty, while the outer error bars represent
the total uncertainty including systematic error. The green dash-dot and cyan
dashed lines represent the phenomenological Hall B and MAID models [18, 19]
respectively. The form factor term of the integrand for $\Delta_{1}$ is
constructed using the Arrington form factor fit [20]. The red line indicates a
new phenomenological fit to the data and extrapolation to low $Q^{2}$=0 and
high $Q^{2}$, with the red band representing the uncertainty of the
calculation. The results are similar but have different mass scaling in
electronic hydrogen.

$\Delta_{\mathrm{pol}}$ currently dominates the theoretical uncertainty of HFS
calculations [7, 1] and it is evident that to calculate this contribution
accurately, we must examine experimental measurements of the spin structure
functions $g_{1}$ and $g_{2}$. Notably, the low-$Q^{2}$ regime dominates the
integrals of Eqs. (5) and (6) due to the $\frac{1}{Q^{4}}$ factors, so it is
especially vital to determine the spin structure functions at low $Q^{2}$ if
we wish to fully understand the hydrogen atom, and by extension the HFS effect
in general.

Figure 3: The polarizability contribution to the hyperfine splitting for
muonic hydrogen. The analysis of this work is shown in a red circle, and is
compared to previous data-driven dispersion relation calculations [4, 5, 6]
shown in orange squares, and the baryon chiral perturbation theory calculation
[21] shown in blue triangles.

In the following, we focus on a new empirical input for the proton spin
structure functions. Our evaluation details new results from Jefferson Lab
Experiments E03-006 (EG4) and E08-027 (g2p), two complementary experiments
both aimed at collecting low $Q^{2}$ data with longitudinally and transversely
polarized proton (NH3) targets, respectively [16, 17]. Here, we present these
experiments’ contribution to the hyperfine integrals above, with the g2p data
providing the first data-driven extraction of $\Delta_{2}$, and the EG4 data
providing new $\Delta_{1}$ data with unprecedented coverage in the low-$Q^{2}$
region.

In the EG4 experiment, 1-3.5 nA of longitudinally polarized electron beam was
incident on a longitudinally polarized NH3 target. The scattered electrons
were detected using the CEBAF Large Acceptance Spectrometer (CLAS). The
longitudinal polarized cross section difference $\Delta\sigma_{\parallel}$ was
directly extracted from the yield difference between left- and right-handed
beam electrons, such that contributions from the unpolarized material cancel.
Combined with an estimation of the (small) transverse contribution based on a
parameterization of world data, the proton structure function $g_{1}$ was
extracted. Four different beam energies, along with the very small scattering
angle down to 6∘, allowed the minimum $Q^{2}$ to reach a very low 0.012 GeV2.
These $g_{1}$ results were used to form the bulk of the low $Q^{2}$
$\Delta_{1}$ data presented in this letter, see Fig. 2 (top panel).

In the g2p experiment, the parallel and perpendicular double spin asymmetries
$A_{\parallel}$ and $A_{\perp}$ were measured for the scattering of 50 nA
polarized electrons on longitudinally and transversely polarized NH3 targets,
respectively. Scattered electrons were detected at an angle of $\approx$6.5∘
using the Hall A High Resolution Spectrometers and a Septa Magnet. Measured
asymmetries were combined with unpolarized cross section models from the
Bosted-Christy phenomenological fit [22] to form polarized cross section
differences, which were used to extract the spin structure functions. By
varying the polarized target magnetic field and electron beam energy, five
different kinematic settings were measured ranging from $Q^{2}$ of 0.02 GeV2
to 0.12 GeV2. Four of these settings were measured with a transverse polarized
target field, giving rise to a perpendicular polarized cross section
difference and a $g_{2}$ result, and one setting with a longitudinally
polarized target field, which provides a parallel polarized cross section
difference and a $g_{1}$ result. The results from the g2p experiment are the
first data in a range relevant to the HFS, and so are used to form the
$\Delta_{2}$ results in this letter.

Results for the $\Delta_{1}$ integrand are shown in the top of Fig. 2. The
unmeasured part of the integral, largely at low Bjorken-$x$, is estimated
using the CLAS Hall B model [18]. This is the best available model, containing
significantly more modern $g_{1}$ data than the Simula parametrization [23]
used in previous analyses [4]. A new phenomenological fit, shown in red, is
generated to extrapolate to the low $Q^{2}$ region. Details on the fitting
procedure can be found in the Supplemental Materials. Numerical results for
these contributions are obtained by integrating over the data where they
exist, primarily the EG4 data shown [17], as well as data from the EG1b
experiment in the $Q^{2}$ = 1.0-5.0 GeV2 region [18]. The contribution from
the low-$Q^{2}$ regime is calculated by integrating the displayed
extrapolation fit, while the high-$Q^{2}$ contribution above $Q^{2}$ = 5.0
GeV2 is calculated using the Hall B Model [18].

Results for the $\Delta_{2}$ integrand in Eq. (6) are shown in the bottom of
Fig. 2. The unmeasured part of the $dx$ integral is again estimated using the
Hall B model [18]. The results of g2p shown are the first ever direct
experimental extractions of this quantity. The low and high $Q^{2}$ regions
are calculated using the displayed fit, which is described in detail in the
Supplemental Materials. Due to the comparative lack of $g_{2}$ data, the
extrapolation has a somewhat larger error than for the $\Delta_{1}$ results.

This historical lack of $g_{2}$ data makes it difficult to conclude if the
Hall B model [18] is a good estimation of the low-$x$ region or not. To
account for this, we compare the result using the older Simula parametrization
[23], which contains a significantly different prediction for the low-$x$
behaviour of $g_{2}$, and include the difference in our extrapolation error by
comparing the upper and lower error bands of our extrapolating fit to the data
in each case. Despite the very different models, this error contribution is
relatively small, because the low-$x$ region is suppressed for $\Delta_{2}$.

The integrated results for $\Delta_{1}$, $\Delta_{2}$, and
$\Delta_{\mathrm{pol}}$ are as follows:

$\displaystyle\Delta_{1}^{e}$ $\displaystyle=$ $\displaystyle 6.78\pm
1.02(\text{data})\pm 0.24(\text{extrapolation})$ (7)
$\displaystyle\Delta_{1}^{\mu}$ $\displaystyle=$ $\displaystyle 5.69\pm
0.84(\text{data})\pm 0.20(\text{extrapolation})$ (8)
$\displaystyle\Delta_{2}^{e}$ $\displaystyle=$ $\displaystyle-1.98\pm
0.16(\text{data})\pm 0.38(\text{extrapolation})$ (9)
$\displaystyle\Delta_{2}^{\mu}$ $\displaystyle=$ $\displaystyle-1.40\pm
0.11(\text{data})\pm 0.31(\text{extrapolation})\qquad$ (10)
$\displaystyle\Delta_{\text{pol}}^{e}$ $\displaystyle=$ $\displaystyle
1.09\text{ ppm}\pm 0.31\text{ ppm}$ (11)
$\displaystyle\Delta_{\text{pol}}^{\mu}$ $\displaystyle=$ $\displaystyle
200.6\text{ ppm}\pm 52.4\text{ ppm}\qquad$ (12)

The total polarizability contribution to the hyperfine splitting is provided
in parts per million (ppm) of the Fermi energy $E_{F}$. The uncertainties for
$\Delta_{1}$ and $\Delta_{2}$ are divided into uncertainty coming directly
from the data, and a combined systematic uncertainty coming from the
extrapolations into high and low $Q^{2}$ regions and into the low-$x$ regime.
The extrapolation error is calculated by generating pseudo-data within the
data’s error bars, and calculating a new fit to this pseudo-data. This
procedure is repeated 1000 times, and the standard deviation in the resulting
fits is taken as the extrapolation error band. The total extrapolation error
in the table also includes a contribution from the choice of low-$x$ fill-in
model, where the highest upper band and lowest lower band achievable with
different choices of fill-in model are taken as the absolute limits of the
error band. The data error is a combination of statistical and systematic
error from the respective experiments contributing [16, 17, 18].

As can be seen in Fig. 3, the new data from EG4 and g2p dramatically reduce
the long standing discrepancy between the leading-order (LO) $\chi$PT
prediction [15, 8] and earlier data-driven dispersive evaluations [4, 5, 6] of
the polarizability contribution. The large difference from the earlier
dispersive results is illustrative of the importance of low-$Q^{2}$ data for
$\Delta_{\mathrm{pol}}$, and the improvement in the available phenomenological
models [18], which are constrained by a larger amount of data compared to
earlier parametrizations used in previous analyses [23, 24].

Our new data-driven evaluations of $\Delta_{\mathrm{pol}}$ put us in the
unique position to update the theoretical predictions of the HFS in (muonic)
hydrogen, as well as the extractions of the proton Zemach radius from
measurements of the HFS. The $1S$ HFS in H is extraordinarily well-measured
[25, 26]:

$E^{\,\text{exp.}}_{1S\text{-HFS}}(\mathrm{H})=1\,420.405\,751\,768(2)\,\text{MHz}.$
(13)

Therefore, the presently most precise extraction of $R_{\mathrm{Z}}$ from
spectroscopy is achieved when comparing the measured $1S$ HFS in H to the full
theory prediction including QED, electroweak and strong interaction effects:

$E_{nS\text{-HFS}}=\frac{E_{\mathrm{F}}}{n^{3}}\left(1+\Delta_{\mathrm{QED}}+\Delta_{\mathrm{weak}}+\Delta_{\mathrm{strong}}\right).$
(14)

For details on the numerical factors entering Eq. (14), we refer to the
compilations in [7, Eq. (40) and (42)]. Here, $\Delta_{\mathrm{strong}}$
contains the 2$\gamma$-exchange contributions introduced in Eq. (2), as well
as other hadronic corrections such as hadronic vacuum polarization. The
$\mathcal{O}(\alpha^{5})$ recoil corrections $\Delta_{\mathrm{recoil}}$ are
taken from [27, Eq. (14) and (15)] and are consistent with our choice of
$F_{2}$. Since the $F_{2}$ term in $\Delta_{\mathrm{pol}}$ cancels exactly
with a corresponding term in $\Delta_{\mathrm{recoil}}$, we do not need to
take into account uncertainties of the $F_{2}$ parametrization in our
$\Delta_{\mathrm{pol}}$ evaluation. Note that radiative corrections (e.g.,
through electronic vacuum polarization) to the $2\gamma$-exchange diagram are
taken into account as well. The proton Zemach radius is then extracted from
the $1S$ HFS in H as:

$R_{\mathrm{Z}}=1.036(8)\,\text{fm. }$ (15)

This result is more precise than previous extractions [1, 28], as can be seen
from the top panel of Fig. 4, where extractions of the Zemach radius from the
measured $1S$ HFS in H are shown, assuming the same theoretical prediction of
the HFS, but different values of $\Delta_{\mathrm{pol}}$ as shown in Fig. 3.
Our evaluation of $\Delta_{\mathrm{pol}}$ based on new data for the proton
spin structure functions (red), can be considered an update of the previous
dispersive analysis [1] (orange open circle). It has moved closer to the
extraction based on LO $\chi$PT (blue), but still does not agree. In the
bottom panel of Fig. 4, recent precise results from lattice QCD suggest a
small $R_{\mathrm{Z}}$ [29], in perfect agreement with the LO $\chi$PT
extraction from the 1S H HFS, but in tension with evaluations based on proton
form factors measured in scattering [30, 31]. Our work is compatible with both
[29] and [31].

Given the limited beam time and required tunability of the laser setup,
precise theory guidance is crucial for the experiments planned by the CREMA
[2] and FAMU [3] collaborations. Here we present an updated theory prediction
for the $1S$ HFS in $\mu$H, based on the theory compilation in [7, Eq. (40)],
and substituting our $\Delta_{\mathrm{pol}}$, as well as our extraction of the
Zemach radius from the $1S$ HFS in H, cf. Eq. (15):

$E^{\,\text{th.}}_{1S\text{-HFS}}(\mu\mathrm{H})=182.636(16)\,\text{meV}.$
(16)

This result is in perfect agreement with the presently most precise prediction
presented in [7, Eq. (48)] that is based on scaling the result from the very
accurate electronic hydrogen experimental HFS measurement, cf. Eq. (13), and
calculating the small (including small in the uncertainty limits) nonscaling
corrections, as suggested in [32, 33].

Figure 4: Comparison of various extractions of the proton Zemach radius
$R_{\mathrm{Z}}$ [4, 8, 29, 30, 31].

We present results for the first ever experimental data in a regime which
contributes significantly to the integrals of the $\Delta_{1}$ and
$\Delta_{2}$ Hyperfine Splitting contributions. This new data provides
previously lacking guidance on how to constrain theoretical calculations of
the Hyperfine Splitting effect. Previous data-driven work to determine these
quantities [4] has been limited by older results lacking inelastic proton spin
structure function data in the low $Q^{2}$ regime. These new results are much
closer to agreement with $\chi$PT calculations of $\Delta_{\mathrm{pol}}$,
strongly reducing the long-standing tension between methods [7] and reducing
the overall error on the polarizability contribution by a factor of two. The
reduction in uncertainty provided by our data is crucial in order to
facilitate the search for the narrow $\mu$H $1S$ HFS in the planned
experiments by the CREMA [2] and FAMU [3] Collaborations, as well as to
interpret these future measurements. With these experiments aiming at up to 1
ppm relative precision, they have the potential to provide novel insights into
the magnetic structure of the proton.

###### Acknowledgements.

We thank Aldo Antognini and Vadim Lensky for their comments and discussions.
This work was supported by the United States Department of Energy (DOE) under
award numbers DE-FG02-88ER40410, DE-SC0024665, and DE–SC0014434, the National
Science Foundation under grant PHY-181236, and by the Schweizerische
Nationalfonds. This work is also supported by the Deutsche
Forschungsgemeinschaft (DFG) through the Emmy Noether Programme Grant
449369623, the Research Unit FOR5327 Grant 458854507, and by the Swiss
National Science Foundation (SNSF) through the Ambizione Grant PZ00P2_193383.
This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Office of Nuclear Physics under contract DE-AC05-06OR23177.

## References

  * Volotka  _et al._ [2005] A. V. Volotka, V. M. Shabaev, G. Plunien, and G. Soff, The European Physical Journal D - Atomic, Molecular, Optical and Plasma Physics, 33, 23 (2005), ISSN 1434-6079.
  * Amaro  _et al._ [2022] P. Amaro _et al._ , SciPost Phys., 13, 020 (2022), arXiv:2112.00138 [physics.atom-ph] .
  * Pizzolotto  _et al._ [2020] C. Pizzolotto _et al._ , Eur. Phys. J. A, 56, 185 (2020).
  * Carlson  _et al._ [2008] C. E. Carlson, V. Nazaryan, and K. Griffioen, Phys. Rev. A, 78, 022517 (2008).
  * Tomalak, Oleksandr [2019] Tomalak, Oleksandr, Eur. Phys. J. A, 55, 64 (2019).
  * Faustov, R and Gorbacheva I and Martynenko A [2006] Faustov, R and Gorbacheva I and Martynenko A, Bellingham, WA (2006).
  * Antognini  _et al._ [2022] A. Antognini, F. Hagelstein, and V. Pascalutsa, Annual Review of Nuclear and Particle Science, 72, 389 (2022).
  * Hagelstein  _et al._ [2023] F. Hagelstein, V. Lensky, and V. Pascalutsa, The European Physical Journal C, 83, 762 (2023a), ISSN 1434-6052.
  * Pohl  _et al._ [2010] R. Pohl _et al._ , Nature, 466, 213 (2010).
  * Mohr  _et al._ [2012] P. J. Mohr, B. N. Taylor, and D. B. Newell, Rev. Mod. Phys., 84, 1527 (2012).
  * Carlson [2015] C. E. Carlson, Prog. Part. Nucl. Phys., 82, 59 (2015), arXiv:1502.05314 [hep-ph] .
  * Pohl  _et al._ [2013] R. Pohl, R. Gilman, G. A. Miller, and K. Pachucki, Ann. Rev. Nucl. Part. Sci., 63, 175 (2013), arXiv:1301.0905 [physics.atom-ph] .
  * Karr  _et al._ [2020] J.-P. Karr, D. Marchand, and E. Voutier, Nature Rev. Phys., 2, 601 (2020).
  * Gao and Vanderhaeghen [2022] H. Gao and M. Vanderhaeghen, Rev. Mod. Phys., 94, 015002 (2022), arXiv:2105.00571 [hep-ph] .
  * Hagelstein and Pascalutsa [2016] F. Hagelstein and V. Pascalutsa, PoS, CD15, 077 (2016).
  * Ruth  _et al._ [2022] D. Ruth, R. Zielinski, C. Gu, _et al._ , Nature Physics, 18, 1441 (2022).
  * Zheng  _et al._ [2021] X. Zheng, A. Deur, H. Kang, _et al._ , Nature Physics, 17, 736 (2021).
  * Fersch  _et al._ [2017] R. Fersch _et al._ (CLAS Collaboration), Phys. Rev. C, 96, 065208 (2017).
  * Drechsel  _et al._ [2007] D. Drechsel, S. S. Kamalov, and L. Tiator, The European Physical Journal A, 34, 69 (2007), ISSN 1434-601X.
  * Arrington [2004] J. Arrington, Phys. Rev. C, 69, 022201 (2004).
  * Hagelstein  _et al._ [2023] F. Hagelstein, V. Lensky, and V. Pascalutsa, Eur. Phys. J. C, 83, 762 (2023b), arXiv:2305.09633 [nucl-th] .
  * Christy and Bosted [2010] M. E. Christy and P. E. Bosted, Phys. Rev. C, 81, 055213 (2010).
  * Simula  _et al._ [2002] S. Simula, M. Osipenko, G. Ricco, and M. Taiuti, Phys. Rev., D65, 034017 (2002a).
  * Simula  _et al._ [2002] S. Simula _et al._ , in  _Proceedings, 2nd International Symposium GDH 2002, Genova, Italy, July 3-6, 2002INSPIRE_ (2002) arXiv:nucl-th/0212031 [nucl-th] .
  * Hellwig  _et al._ [1970] H. Hellwig, R. F. C. Vessot, M. W. Levine, P. W. Zitzewitz, D. W. Allan, and D. J. Glaze, IEEE Transactions on Instrumentation and Measurement, 19, 200 (1970), ISSN 0018-9456.
  * Karshenboim [2000] S. G. Karshenboim, Can. J. Phys., 78, 639 (2000), arXiv:physics/0008051 .
  * Antognini  _et al._ [2022] A. Antognini, Y.-H. Lin, and U.-G. Meißner, Phys. Lett. B, 835, 137575 (2022), arXiv:2208.04025 [nucl-th] .
  * Antognini  _et al._ [2013] A. Antognini, F. Kottmann, F. Biraben, P. Indelicato, F. Nez, and R. Pohl, Annals Phys., 331, 127 (2013), arXiv:1208.2637 [physics.atom-ph] .
  * Djukanovic  _et al._ [2023] D. Djukanovic, G. von Hippel, H. B. Meyer, K. Ottnad, M. Salg, and H. Wittig, (2023), arXiv:2309.17232 [hep-lat] .
  * Lin  _et al._ [2022] Y.-H. Lin, H.-W. Hammer, and U.-G. Meißner, Phys. Rev. Lett., 128, 052002 (2022), arXiv:2109.12961 [hep-ph] .
  * Distler  _et al._ [2011] M. O. Distler, J. C. Bernauer, and T. Walcher, Phys. Lett. B, 696, 343 (2011), arXiv:1011.1861 [nucl-th] .
  * Tomalak [2018] O. Tomalak, Eur. Phys. J. A, 54, 3 (2018), arXiv:1709.06544 [hep-ph] .
  * Peset and Pineda [2017] C. Peset and A. Pineda, JHEP, 04, 060 (2017), arXiv:1612.05206 [nucl-th] .


are uniformly bounded, we can argue as in the case $p=4$ to obtain the desired
result.

##### Covariances

In this part we verify that (58) holds when $p\neq q$, for $p=2,\ldots,4$.
Below we will use the notation $v_{\\_}{3}=\int x^{3}\nu(dx)$.

$p=1$:

The case $q=2$ can be analysed similarly to the case $p=q=2$. Furthermore,
since

$\displaystyle\frac{1}{\Delta_{\\_}{n}n}\sum_{\\_}{j=l_{\\_}{n}+1}^{n-1}\mathbb{E}(\xi_{\\_}{j,n}^{(1)}\xi_{\\_}{j,n}^{(3)}\mid\mathscr{G}_{\\_}{j-1}^{n})=$
$\displaystyle
v_{\\_}{3}\frac{1}{\Delta_{\\_}{n}n}\sum_{\\_}{k=l_{\\_}{n}}^{n-2}\chi_{\\_}{k}\theta_{\\_}{k-1,k}^{(4)},$
$\displaystyle\frac{1}{\Delta_{\\_}{n}n}\sum_{\\_}{j=l_{\\_}{n}+1}^{n-1}\mathbb{E}(\xi_{\\_}{j,n}^{(1)}\xi_{\\_}{j,n}^{(4)}\mid\mathscr{G}_{\\_}{j-1}^{n})=$
$\displaystyle-\frac{v_{\\_}{3}}{n\Delta_{\\_}{n}}\sum_{\\_}{i=0}^{n-1-l_{\\_}{n}}\sum_{\\_}{m=i}^{n-2}\alpha_{\\_}{i,m}\theta_{\\_}{i,m,m^{\prime}}^{(1)},$
$\displaystyle\frac{1}{\Delta_{\\_}{n}n}\sum_{\\_}{j=l_{\\_}{n}+1}^{n-1}\mathbb{E}(\xi_{\\_}{j,n}^{(1)}\xi_{\\_}{j,n}^{(5)}\mid\mathscr{G}_{\\_}{j-1}^{n})=$
$\displaystyle-\frac{v_{\\_}{3}}{n\Delta_{\\_}{n}}\sum_{\\_}{k=l_{\\_}{n}+1}^{n-2}\sum_{\\_}{m=k}^{n-2}\alpha_{\\_}{k,m}\theta_{\\_}{k,m}^{(4)}$

in which $\theta_{\\_}{i,m,m^{\prime}}^{(1)}$ and $\theta_{\\_}{k,m}^{(4)}$
are is as in (63) and (61), respectively, we can reproduce the argument used
above for $p=q$, $p=3,4,5$, to conclude that (58) is also satisfied in this
situation.

$p=2$:

By the Cauchy-Schwartz inequality and Assumption (3)in [25], for all
$l_{\\_}{n}+1\leq j<2l_{\\_}{n}$

$\mathbb{E}\left[\left|\mathbb{E}(\xi_{\\_}{j,n}^{(2)}\xi_{\\_}{j,n}^{(3)}\mid\mathscr{G}_{\\_}{j-1}^{n})\right|\right]\leq
C\mathbb{E}\left[\left(\sum_{\\_}{k=l_{\\_}{n}}^{j-1}\chi_{\\_}{k}\mathrm{Leb}(\cup_{\\_}{i=0}^{k-l_{\\_}{n}}\mathcal{\mathcal{P}}_{\\_}{A}^{n}(i,j))\right)^{2}\right]^{1/2}\leq
C\Delta_{\\_}{n},$

which easily implies that

$\frac{1}{\Delta_{\\_}{n}n}\sum_{\\_}{j=l_{\\_}{n}+1}^{2l_{\\_}{n}}\mathbb{E}(\xi_{\\_}{j,n}^{(2)}\xi_{\\_}{j,n}^{(3)}\mid\mathscr{G}_{\\_}{j-1}^{n})\overset{\mathbb{P}}{\rightarrow}0.$

For $j\geq 2l_{\\_}{n}+1$ we have that

$\displaystyle\mathbb{E}(\xi_{\\_}{j,n}^{(2)}\xi_{\\_}{j,n}^{(3)}\mid\mathscr{G}_{\\_}{j-1}^{n})=$
$\displaystyle\beta_{\\_}{j,j}^{(2)}\vartheta_{\\_}{j}^{(1)}+\beta_{\\_}{j,j}^{(2)}\vartheta_{\\_}{j}^{(2)},$

where

$\vartheta_{\\_}{j}^{(1)}:=\sum_{\\_}{k=l_{\\_}{n}}^{j-l_{\\_}{n}-1}\chi_{\\_}{k}\mathbb{E}[(\beta_{\\_}{k,j}^{(1)})^{2}];\,\,\,\vartheta_{\\_}{j}^{(2)}:=\sum_{\\_}{k=j-l_{\\_}{n}}^{j-1}\chi_{\\_}{k}\mathbb{E}[(\beta_{\\_}{k,j}^{(1)})^{2}]$

Note that $\beta_{\\_}{j,j}^{(2)}\vartheta_{\\_}{j}^{(1)}$ and
$\beta_{\\_}{j,j}^{(2)}\vartheta_{\\_}{j}^{(2)}-\mathbb{E}(\beta_{\\_}{j,j}^{(2)}\vartheta_{\\_}{j}^{(2)})$
are $l_{\\_}{n}$-uncorrelated. Furthermore, by Rosenthal’s inequality and
relation (44), we get that

$\mathbb{E}[\mid\vartheta_{\\_}{j}^{(v)}\mid{}^{4}]\leq
C\Delta_{\\_}{n}^{4},\,\,\,v=1,2,$

which in turn implies that

$\displaystyle\mathbb{E}[(\beta_{\\_}{j,j}^{(2)}\vartheta_{\\_}{j}^{(v)})^{2}]$
$\displaystyle\leq C\Delta_{\\_}{n}^{2},\,\,\,v=1,2.$

Using this bound, we deduce as in (59) that

$\frac{1}{n\Delta_{\\_}{n}}\sum_{\\_}{j=2l_{\\_}{n}+1}^{n-1}\beta_{\\_}{j,j}^{(2)}\vartheta_{\\_}{j}^{(1)}+\frac{1}{n\Delta_{\\_}{n}}\sum_{\\_}{j=2l_{\\_}{n}+1}^{n-1}\left[\beta_{\\_}{j,j}^{(2)}\vartheta_{\\_}{j}^{(2)}-\mathbb{E}(\beta_{\\_}{j,j}^{(2)}\vartheta_{\\_}{j}^{(2)})\right]=\mathrm{o}_{\\_}{\mathbb{P}}(1),$

which is clearly enough for (58) for the case $q=3$. Furthermore, since for
$j\geq 2l_{\\_}{n}+1$

$\displaystyle\mathbb{E}(\xi_{\\_}{j,n}^{(2)}\xi_{\\_}{j,n}^{(4)}\mid\mathscr{G}_{\\_}{j-1}^{n})=$
$\displaystyle-\beta_{\\_}{j,j}^{(2)}(\vartheta_{\\_}{j}^{(3)}+\vartheta_{\\_}{j}^{(4)}),$
$\displaystyle\mathbb{E}(\xi_{\\_}{j,n}^{(2)}\xi_{\\_}{j,n}^{(5)}\mid\mathscr{G}_{\\_}{j-1}^{n})=$
$\displaystyle-\beta_{\\_}{j,j}^{(2)}(\vartheta_{\\_}{j}^{(5)}+\vartheta_{\\_}{j}^{(6)}+\vartheta_{\\_}{j}^{(7)});$

where

$\displaystyle\vartheta_{\\_}{j}^{(3)}:=\sum_{\\_}{i=0}^{j-l_{\\_}{n}-1}\sum_{\\_}{m=i}^{j-l_{\\_}{n}-1}\alpha_{\\_}{i,m}\theta_{\\_}{i,m,j}^{(6)};$
$\displaystyle\,\,\vartheta_{\\_}{j}^{(4)}:=\sum_{\\_}{i=0}^{j-l_{\\_}{n}-1}\sum_{\\_}{m=j-l_{\\_}{n}}^{j-1}\alpha_{\\_}{i,m}\theta_{\\_}{i,m,j}^{(6)};$
$\displaystyle\vartheta_{\\_}{j}^{(5)}:=\sum_{\\_}{k=l_{\\_}{n}+1}^{j-l_{\\_}{n}}\sum_{\\_}{m=k}^{j-l_{\\_}{n}-1}\alpha_{\\_}{k,m}\theta_{\\_}{k,j}^{(7)};$
$\displaystyle\,\,\vartheta_{\\_}{j}^{(6)}:=\sum_{\\_}{k=l_{\\_}{n}+1}^{j-l_{\\_}{n}}\sum_{\\_}{m=j-l_{\\_}{n}}^{j-1}\alpha_{\\_}{k,m}\theta_{\\_}{k,j}^{(7)},$

and
$\vartheta_{\\_}{j}^{(7)}=\sum_{\\_}{k=j-l_{\\_}{n}+1}^{j-1}\sum_{\\_}{m=k}^{j-1}\alpha_{\\_}{k,m}\theta_{\\_}{k,j}^{(7)}$,
in which

$\displaystyle\theta_{\\_}{i,m,j}^{(6)}\lor\theta_{\\_}{i,m,j}^{(7)}$
$\displaystyle\leq
C\Delta_{\\_}{n}a(t_{\\_}{j}-t_{\\_}{m}-(l_{\\_}{n}+1)\Delta_{\\_}{n}),$

and
$\beta_{\\_}{j,j}^{(2)}\vartheta_{\\_}{j}^{(l)}-\mathbb{E}(\beta_{\\_}{j,j}^{(2)}\vartheta_{\\_}{j}^{(v)})$
are $l_{\\_}{n}$-uncorrelated for $v=3,4,5,6$, we can repeat the preceding
arguments to conclude that (58) is also valid for these cases.

$p=3$:

By decomposing

$\beta_{\\_}{k^{\prime},j}^{(2)}=\beta_{\\_}{k^{\prime},k}^{(1)}+\beta_{\\_}{k^{\prime},k}^{(2)}+\sum_{\\_}{m=k+1}^{j-1}\beta_{\\_}{k^{\prime},m}^{(1)},$
(64)

for $k^{\prime}\leq k\leq j-2$, we obtain, due to (57), that

$\displaystyle\sum_{\\_}{j=l_{\\_}{n}+1}^{n-1}\mathbb{E}(\xi_{\\_}{j,n}^{(3)}\xi_{\\_}{j,n}^{(4)}\mid\mathscr{G}_{\\_}{j-1}^{n})=$
$\displaystyle\sum_{\\_}{j=l_{\\_}{n}+1}^{n-1}\mathbb{E}(\xi_{\\_}{j,n}^{(3)}\xi_{\\_}{j,n}^{(4)})$
(65)
$\displaystyle+\sum_{\\_}{m=2l_{\\_}{n}+2}^{n-1}\sum_{\\_}{k=2l_{\\_}{n}+1}^{m-1}\sum_{\\_}{k^{\prime}=l_{\\_}{n}}^{k-l_{\\_}{n}-1}\chi_{\\_}{k}\beta_{\\_}{k^{\prime},m}^{(1)}\theta_{\\_}{m,k^{\prime}}^{(8)}$
$\displaystyle+\sum_{\\_}{k=2l_{\\_}{n}+1}^{n-2}\sum_{\\_}{k^{\prime}=l_{\\_}{n}}^{k-l_{\\_}{n}-1}\chi_{\\_}{k}\beta_{\\_}{k^{\prime},k}^{(2)}\theta_{\\_}{k,k^{\prime}}^{(8)}$
$\displaystyle+\sum_{\\_}{k=2l_{\\_}{n}+1}^{n-2}\sum_{\\_}{k^{\prime}=l_{\\_}{n}}^{k-l_{\\_}{n}-1}[\chi_{\\_}{k}\beta_{\\_}{k^{\prime},k}^{(1)}-\mathbb{E}(\chi_{\\_}{k}\beta_{\\_}{k^{\prime},k}^{(1)})]\theta_{\\_}{k,k^{\prime}}^{(8)},$

where
$\theta_{\\_}{k,k^{\prime}}^{(8)}:=\mathrm{Leb}(\cup_{\\_}{j=k+1}^{n}\mathcal{\mathcal{P}}_{\\_}{A}^{n}(k^{\prime}+1,j)).$
Note that once again all the error terms are sums of
$(\mathscr{G}_{\\_}{k})$-increment martingales. Due to Lemma 2, we only need
to verify that the second moment of each of these summands is of order
$\mathrm{O}(\Delta_{\\_}{n})$. For the last two sums, this follows easily by
using that $\mathbb{E}[(\beta_{\\_}{k^{\prime},k}^{(2)})^{2}]\leq C,$ and

$\mathbb{E}(\chi_{\\_}{k}^{2}\beta_{\\_}{k^{\prime},k}^{(1)}\beta_{\\_}{m,k}^{(1)})\leq
C\Delta_{\\_}{n}a(t_{\\_}{k}-t_{\\_}{k^{\prime}}\land
t_{\\_}{m}),\,\,\,\left|\theta_{\\_}{k,k^{\prime}}^{(8)}\right|\leq
C\Delta_{\\_}{n}a(t_{\\_}{k}-t_{\\_}{k^{\prime}}).$

The required bound for the second sum in (65) is a bit more involved. Let
$\vartheta_{\\_}{m}^{(8)}=\sum_{\\_}{k=2l_{\\_}{n}+1}^{m-1}\sum_{\\_}{k^{\prime}=l_{\\_}{n}}^{k-l_{\\_}{n}-1}\chi_{\\_}{k}\beta_{\\_}{k^{\prime},m}^{(1)}\theta_{\\_}{m,k^{\prime}}^{(8)}$.
In view that $\mathbb{E}[(\beta_{\\_}{q,m}^{(1)})^{2}]\leq
C\Delta_{\\_}{n}a(t_{\\_}{m}-t_{\\_}{q})$, the sought-after bound is obtained
by the following estimates

$\displaystyle\mathbb{E}[(\vartheta_{\\_}{m}^{(8)})^{2}]\leq$ $\displaystyle
C\Delta_{\\_}{n}\sum_{\\_}{k=2l_{\\_}{n}+1}^{m-1}\sum_{\\_}{k^{\prime}=l_{\\_}{n}}^{k-l_{\\_}{n}-1}\sum_{\\_}{q=l_{\\_}{n}}^{k^{\prime}}\mathbb{E}[(\beta_{\\_}{q,m}^{(1)})^{2}]\theta_{\\_}{m,k^{\prime}}^{(8)}\theta_{\\_}{m,q}^{(8)}$
$\displaystyle\leq$ $\displaystyle
C\Delta_{\\_}{n}^{2}\sum_{\\_}{k=2l_{\\_}{n}+1}^{m-1}\int_{\\_}{l_{\\_}{n}\Delta_{\\_}{n}}^{t_{\\_}{k+1}-l_{\\_}{n}\Delta_{\\_}{n}}a(t_{\\_}{m}-r)\left(\int_{\\_}{l_{\\_}{n}\Delta_{\\_}{n}}^{r}a(t_{\\_}{m}-s)^{2}ds\right)dr$
$\displaystyle=$ $\displaystyle
C\Delta_{\\_}{n}^{2}\sum_{\\_}{k=2l_{\\_}{n}+1}^{m-1}\int_{\\_}{t_{\\_}{m}-t_{\\_}{k+1}+l_{\\_}{n}\Delta_{\\_}{n}}^{t_{\\_}{m}-l_{\\_}{n}\Delta_{\\_}{n}}a(x)\left(\int_{\\_}{x}^{t_{\\_}{m}-l_{\\_}{n}\Delta_{\\_}{n}}a(y)^{2}dy\right)dx$
$\displaystyle\leq$ $\displaystyle
C\Delta_{\\_}{n}^{2}\sum_{\\_}{k=2l_{\\_}{n}+1}^{m-1}\int_{\\_}{t_{\\_}{m}-t_{\\_}{k+1}+l_{\\_}{n}\Delta_{\\_}{n}}^{t_{\\_}{m}-l_{\\_}{n}\Delta_{\\_}{n}}a(y)^{2}dy$
$\displaystyle\leq$ $\displaystyle
C\Delta_{\\_}{n}\int_{\\_}{(2l_{\\_}{n}+1)\Delta_{\\_}{n}}^{t_{\\_}{m}}\int_{\\_}{t_{\\_}{m}-u+(l_{\\_}{n}-1)\Delta_{\\_}{n}}^{t_{\\_}{m}-l_{\\_}{n}\Delta_{\\_}{n}}a(y)^{2}dydu$
$\displaystyle\leq$ $\displaystyle
C\Delta_{\\_}{n}\int_{\\_}{(2l_{\\_}{n}+1)\Delta_{\\_}{n}}^{t_{\\_}{m}-l_{\\_}{n}\Delta_{\\_}{n}}\int_{\\_}{w}^{t_{\\_}{m}-l_{\\_}{n}\Delta_{\\_}{n}}a(y)^{2}dydw$
$\displaystyle\leq$ $\displaystyle
C\Delta_{\\_}{n}\int_{\\_}{0}^{t_{\\_}{m}}a(y)^{2}ydy\leq C\Delta_{\\_}{n},$

where in the last step we used the fact that $a(s)=\mathrm{O}(s^{-\alpha})$,
as $s\rightarrow+\infty$ for some $\alpha>1$. For $q=5$, we use the
decomposition

$\sum_{\\_}{q=k+1}^{j-2}\chi_{\\_}{q}=L(\cup_{\\_}{i=0}^{k}\cup_{\\_}{q=k+1}^{j-2}\mathcal{\mathcal{P}}_{\\_}{A}^{n}(i,q))+\sum_{\\_}{i=k}^{j-3}\delta_{\\_}{i,j-1},$

to deduce that

$\displaystyle\sum_{\\_}{j=l_{\\_}{n}+1}^{n-1}\mathbb{E}(\xi_{\\_}{j,n}^{(3)}\xi_{\\_}{j,n}^{(5)}\mid\mathscr{G}_{\\_}{j-1}^{n})=$
$\displaystyle\sum_{\\_}{j=l_{\\_}{n}+1}^{n-1}\mathbb{E}(\xi_{\\_}{j,n}^{(3)}\xi_{\\_}{j,n}^{(4)})+\sum_{\\_}{j=l_{\\_}{n}+1}^{n-1}\vartheta_{\\_}{j}^{(9)}$
$\displaystyle+\sum_{\\_}{k=l_{\\_}{n}}^{n-3}\sum_{\\_}{m=k+1}^{n-2}\alpha_{\\_}{k+1,m}\vartheta_{\\_}{k,m}^{(10)}$

where
$\vartheta_{\\_}{k,m}^{(10)}=\sum_{\\_}{q=l_{\\_}{n}}^{k}\chi_{\\_}{q}\theta_{\\_}{q+1,m,m}^{(5)}$,
and

$\displaystyle\vartheta_{\\_}{j}^{(9)}=$
$\displaystyle\sum_{\\_}{k=l_{\\_}{n}}^{j-3}\sum_{\\_}{i=k}^{j-3}\mathbb{E}[(\beta_{\\_}{k,j}^{(1)})^{2}][\delta_{\\_}{k,j}\delta_{\\_}{k,j-1}-\mathbb{E}(\delta_{\\_}{k,j}\delta_{\\_}{k,j-1})]$
$\displaystyle+\sum_{\\_}{k=l_{\\_}{n}}^{j-3}\delta_{\\_}{k,j}L(\cup_{\\_}{i=0}^{k}\cup_{\\_}{q=k+1}^{j-2}\mathcal{\mathcal{P}}_{\\_}{A}^{n}(i,q)).$

Note that the sum associated with $\vartheta_{\\_}{j}^{(9)}$ can be analysed
in the same way as in the case $p=q=5$. Moreover, the last error term are sums
of $(\mathscr{F}_{\\_}{k})$-increment martingale. After some simple changes of
variables and order of integration we further have that

$\sum_{\\_}{k=l_{\\_}{n}}^{n-3}\sum_{\\_}{m=k+1}^{n-2}\mathbb{E}[(\alpha_{\\_}{k+1,m}\vartheta_{\\_}{k,m}^{(10)})^{2}]\leq
C\int_{\\_}{0}^{n\Delta_{\\_}{n}}\int_{\\_}{0}^{r}\varGamma(s)^{2}dsdr=\mathrm{o}[(n\Delta_{\\_}{n})^{2}],$

which according to Lemma 3, let us conclude that (58) is also valid in this
case.

$p=4$

Lastly, from (64), we obtain easily that

$\sum_{\\_}{j=l_{\\_}{n}+1}^{n-1}\mathbb{E}(\xi_{\\_}{j,n}^{(4)}\xi_{\\_}{j,n}^{(5)}\mid\mathscr{G}_{\\_}{j-1}^{n})=\sum_{\\_}{k=2l_{\\_}{n}+1}^{n-3}\sum_{\\_}{m=k+1}^{n-2}\alpha_{\\_}{k+1,m}\vartheta_{\\_}{k,m}^{(11)},$

with

$\vartheta_{\\_}{k,m}^{(11)}=\sum_{\\_}{q=l_{\\_}{n}}^{k-l_{\\_}{n}-1}\sum_{\\_}{p=q-l_{\\_}{n}}^{n-2}\beta_{\\_}{q,p}^{(1)}\mathrm{Leb}(\cup_{\\_}{j=m\lor
p+1}^{n-2}\mathcal{\mathcal{P}}_{\\_}{A}^{n}(q+1,j)).$

Since $\alpha_{\\_}{k+1,m}\vartheta_{\\_}{k,m}^{(11)}$ is an
$(\mathscr{F}_{\\_}{k})$-increment martingale and for $m\geq k+1$,
$\mathbb{E}[(\vartheta_{\\_}{k,m}^{(11)})^{2}]$ is bounded up to a constant by

$\displaystyle\sum_{\\_}{q=l_{\\_}{n}}^{k-l_{\\_}{n}-1}\sum_{\\_}{q^{\prime}=l_{\\_}{n}}^{q}$
$\displaystyle\mathrm{Leb}(\cup_{\\_}{i=0}^{q^{\prime}-l_{\\_}{n}}\cup_{\\_}{p=q-l_{\\_}{n}}^{n-2}\mathcal{\mathcal{P}}_{\\_}{A}^{n}(q+1,p))\mathrm{Leb}(\cup_{\\_}{j=m+1}^{n-2}\mathcal{\mathcal{P}}_{\\_}{A}^{n}(q+1,j))$
$\displaystyle\times\mathrm{Leb}(\cup_{\\_}{j=m+1}^{n-2}\mathcal{\mathcal{P}}_{\\_}{A}^{n}(q^{\prime}+1,j))$

which is in turn uniformly bounded, we deduce as before that (58) holds in
this situation.

∎

###  Appendix C Simulation study: Set-up

####  C.1 Choice of the marginal distributions and trawl functions

In the simulation study, we consider three different marginal distributions
for the Lévy seed $L^{\prime}$:

  * •

Negative binomial distribution: $L^{\prime}\sim\mathrm{NegBin}(m,\theta)$,
where the corresponding probability mass function is given by
$\mathrm{P}(L^{\prime}=x)=\frac{1}{x!}\frac{\Gamma\left(m+x\right)}{\Gamma\left(m\right)}\left(1-\theta\right)^{m}\theta^{x}$
for $x\in\\{0,1,\ldots\\}$. We choose $m=3.2$ and $\theta=0.2$. Then

$\displaystyle\mathbb{E}(L^{\prime})$ $\displaystyle=m\theta/(1-\theta)=0.8,$
$\displaystyle\mathrm{Var}(L^{\prime})$ $\displaystyle=1,$ $\displaystyle
c_{\\_}4(L^{\prime})$
$\displaystyle=m\theta(\theta^{2}+4\theta+1)/(\theta-1)^{4}=2.875.$

  * •

Gamma distribution: $L^{\prime}\sim\Gamma(\alpha_{\\_}g,\sigma_{\\_}g)$, where
$\alpha_{\\_}g>0$ is the shape parameter and $\sigma_{\\_}g>0$ the scale
parameter. Then the corresponding density is given by

$f(x)=\frac{1}{\sigma_{\\_}g^{\alpha_{\\_}g}\Gamma(\alpha_{\\_}{g})}x^{\alpha_{\\_}g-1}e^{-x/\sigma_{\\_}g},$

for $x>0$. We choose $\alpha_{\\_}g=0.64$ and $\sigma_{\\_}g=1.25$. We note
that

$\displaystyle\mathbb{E}(L^{\prime})$
$\displaystyle=\alpha_{\\_}g\sigma_{\\_}g=0.8,$
$\displaystyle\mathrm{Var}(L^{\prime})$
$\displaystyle=\alpha_{\\_}g\sigma_{\\_}g^{2}=1,$ $\displaystyle
c_{\\_}4(L^{\prime})$ $\displaystyle=6\alpha\sigma^{4}=9.375.$

  * •

Gaussian distribution: $L^{\prime}\sim\mathrm{N}(\mu,\sigma^{2})$. We choose
$\mu=0.8$ and $\sigma=1$. We note that

$\displaystyle\mathbb{E}(L^{\prime})$ $\displaystyle=\mu=0.8,$
$\displaystyle\mathrm{Var}(L^{\prime})$ $\displaystyle=\sigma^{2}=1,$
$\displaystyle c_{\\_}4(L^{\prime})$ $\displaystyle=0.$

For the trawl function, we consider two parametric settings:

  * •

Exponential trawl function: $a(x)=\exp(-\lambda x)$, for $\lambda>0$ and
$x\geq 0$. We choose $\lambda=1$.

  * •

SupGamma trawl function: $a(x)=(1+x/\overline{\alpha})^{-H}$, for
$\overline{\alpha}>0,H>1$ and $x\geq 0$. Note that for $H\in(1,2]$, the
resulting trawl process exhibits long memory and for $H>2$ short memory. We
choose $\overline{\alpha}=0.1,H=1.5$ to study the long-memory setting.

Using the R package ambit, we simulate trawl processes on a grid of width
$\Delta_{\\_}n\in\\{0.1,0.01\\}$ for $n=10,000$ and consider relevant subsets
of the trawl processes when $n\in\\{1000,2000,5000,10000\\}$.

We will consider the observations on the grid $i\Delta_{\\_}n$ for
$i\in\\{0,1,\ldots,n-1\\}$ for the trawl function estimation in the following.

All results are reported based on 200 Monte Carlo runs.

####  C.2 Fixed-$t$ versus fixed-$i$ reporting

Recall that our asymptotic theory relies on a double-asymptotic scheme where
we need simultaneously that $\Delta_{\\_}n\to 0$, $n\to\infty$ and
$n\Delta_{\\_}n\to\infty$.

Moreover, recall that we set $n\Delta_{\\_}n^{3}\to\mu$ and
$n\Delta_{\\_}n^{2}\to\mu_{\\_}0$. Throughout our simulation experiments, we
shall focus on the case when $\mu=0$ (and $\mu_{\\_}0=0$) which suggests a
relation of $\Delta_{\\_}n=O(n^{-\beta})$ for $\beta\in(1/3,1)$ (or
$\beta\in(1/2,1)$). In finite samples, we will not be able to determine the
particular rate $\beta$ of the double-asymptotic regime, hence we focus on two
choices of $\Delta_{\\_}n\in\\{0.1,0.01\\}$ and three choices of $n$ (either
$n\in\\{1000,5000,10000\\}$ or $n\in\\{2000,5000,10000\\}$) to cover different
in-fill and long-span asymptotic settings.

We will be comparing results using a ”fixed-$t$” and a ”fixed-$i$” notion,
which we describe next.

We are interested in estimating the trawl function $a(t)$ at time $t$. Since
the estimator of the trawl function is piecewise constant between
observations, we will focus on estimating $a$ on the grid points only, i.e.
when $t=i\Delta_{\\_}n$ for $i\in\\{0,1,\ldots,n-1\\}$.

We will first consider the case when we fix $t$ and compare the performance of
the estimators for two choices of $\Delta_{\\_}n$ and 3 choices of $n$.

We note that the observations are on the grid $i\Delta_{\\_}n$ for
$i=0,1,\ldots,n-1$. We estimate $a(t)$ for suitable $t=i\Delta_{\\_}n$.

#####  C.2.1 Fixed-$t$ reporting

First of all, we fix the time points $t$ and consider two different values for
the width of the observation grid, i.e. $\Delta_{\\_}n\in\\{0.1,0.01\\}$. We
then choose $i$ such that the fixed-$t$ values are obtained, please see Table
2 for details.

$t$ | 0 | 0.1 | 0.5 | 1 | 5 | 10
---|---|---|---|---|---|---
If $\Delta_{\\_}n=0.1$, then $i$ | 0 | 1 | 5 | 10 | 50 | 100
If $\Delta_{\\_}n=0.01$, then $i$ | 0 | 10 | 50 | 100 | 500 | 1000
Table 2:  Recall that $t=i\Delta_{\\_}n$. n the simulation study in this
subsection, we fix the values of $t$.

When $t$ is fixed, we check the impact of the sample size $n$ and consider
three settings: $n\in\\{2000,5000,10000\\}$333Note that we are not choosing
$n=1000$ here, as in the case when we fix $i$, since that $i=1000$ is required
to reach the time point $t=10$ when using the finer grid $\Delta_{\\_}n=0.01$.

The fixed-$t$ reporting scheme has the clear advantage that we can compare
estimates of the same object, $a(t)$ for fixed-$t$, when varying both
$\Delta_{\\_}n$ and $n$.

Let us now consider the two examples of the trawl functions used in the
simulation study.

For an exponential trawl function with parameter $\lambda=1$, we present the
true values for which we estimate the trawl function in Table 3.

$t$ | 0 | 0.1 | 0.5 | 1 | 5 | 10
---|---|---|---|---|---|---
$a(t)=\exp(-\lambda t)$ | 1 | 0.905 | 0.607 | 0.368 | 0.007 | 0.000
for $\lambda=1$ |  |  |  |  |  |
Table 3:  Fixed-$t$ reporting: True values for which we estimate the
exponential trawl function in the following simulation study.

For a supGamma trawl function with parameters $\overline{\alpha}=0.1,H=1.5$,
we present the true values for which we estimate the trawl function in Table
4.

$t$ | 0 | 0.1 | 0.5 | 1 | 5 | 10
---|---|---|---|---|---|---
$a(t)=(1+t/\alpha)^{-H}$, | 1 | 0.354 | 0.068 | 0.027 | 0.003 | 0.001
for $\alpha=0.1,H=1.5$ |  |  |  |  |  |
Table 4:  Fixed-$t$ reporting: True values for which we estimate the supGamma
trawl function in the following simulation study.

#####  C.2.2 Fixed-$i$ reporting

In the case of a fixed-$i$ reporting scheme, we again focus on estimating
$a(t)=a(i\Delta_{\\_}n)$. However, now we will fix $i$ rather than $t$ and
hence we estimated the trawl function at different time points when varying
$\Delta_{\\_}n$. To see this, as before, consider two different values for the
width of the observation grid, i.e. $\Delta_{\\_}n\in\\{0.1,0.01\\}$. We then
choose fixed values of $i$ and report in Table 5 the corresponding values of
$t$.

$i$ | 0 | 1 | 5 | 10 | 20 | 30
---|---|---|---|---|---|---
If $\Delta_{\\_}n=0.1$, then $t$ | 0 | 0.1 | 0.5 | 1 | 2 | 3
If $\Delta_{\\_}n=0.01$, then $t$ | 0 | 0.01 | 0.05 | 0.1 | 0.2 | 0.3
Table 5:  Recall that $t=i\Delta_{\\_}n$. In the simulation study in this
subsection, we fix the values of $i$.

Moreover, we check the impact of the sample size $n$ and consider three
settings: $n\in\\{1000,5000,10000\\}$.

In the fixed-$i$ reporting scheme, we cannot directly compare
$a(i\Delta_{\\_}n)$ with $a(i\widetilde{\Delta}_{\\_}n)$ for
$\Delta_{\\_}n\not=\widetilde{\Delta}_{\\_}n$, since we are effectively
evaluating (and estimating) the trawl function at different points in time.
However, here we rather ask the question of how well we can estimate the trawl
function for different lags given fixed width $\Delta_{\\_}n$.

Let us now consider the two examples of the trawl functions used in the
simulation study.

For an exponential trawl function with parameter $\lambda=1$, we present the
true values for which we estimate the trawl function in Table 6.

$i$ | 0 | 1 | 5 | 10 | 20 | 30
---|---|---|---|---|---|---
If $\Delta_{\\_}n=0.1$, then $a(i\Delta_{\\_}n)=\exp(-i\Delta_{\\_}n\lambda)$ | 1 | 0.905 | 0.607 | 0.368 | 0.135 | 0.050
If $\Delta_{\\_}n=0.01$, then $a(i\Delta_{\\_}n)=\exp(-i\Delta_{\\_}n\lambda)$ | 1 | 0.990 | 0.951 | 0.905 | 0.819 | 0.741
for $\lambda=1$ |  |  |  |  |  |
Table 6:  Fixed-$i$ reporting: True values for which we estimate the
exponential trawl function in the following simulation study.

For a supGamma trawl function with parameters $\overline{\alpha}=0.1,H=1.5$,
we present the true values for which we estimate the trawl function in Table
7.

$i$ | 0 | 1 | 5 | 10 | 20 | 30
---|---|---|---|---|---|---
If $\Delta_{\\_}n=0.1$, then $a(i\Delta_{\\_}n)=(1+i\Delta_{\\_}n/\alpha)^{-H}$ | 1 | 0.354 | 0.068 | 0.027 | 0.010 | 0.006
If $\Delta_{\\_}n=0.01$, then $a(i\Delta_{\\_}n)=(1+i\Delta_{\\_}n/\alpha)^{-H}$ | 1 | 0.867 | 0.544 | 0.354 | 0.192 | 0.125
for $\alpha=0.1,H=1.5$ |  |  |  |  |  |
Table 7:  Fixed-$i$ reporting: True values for which we estimate the supGamma
trawl function in the following simulation study.

###  Appendix D Simulation study: Consistency results

First of all, we study the consistency and bias of the proposed estimators in
finite samples. We consider the trawl function estimator $\widehat{a}(t)$ and
the bias-corrected version
$\hat{a}(t)-\frac{1}{2}\Delta_{\\_}n\hat{a}^{\prime}(t)$.

For both estimators, we report their empirical mean, bias and standard
deviation based on 200 Monte Carlo runs.

We report the results for an exponential trawl function, for both fixed-$t$,
see Tables LABEL:tab:Con-NB-Exp-fixedt LABEL:tab:Con-Gamma-Exp-fixedt
LABEL:tab:Con-Gauss-Exp-fixedt, and fixed-$i$, see Tables LABEL:tab:Con-NB-Exp
LABEL:tab:Con-Gamma-Exp LABEL:tab:Con-Gauss-Exp. We observe that for an
exponential trawl, the estimation bias is generally small for both estimators,
with the one for the bias-corrected version, as expected, being even smaller.

We report the results for a supGamma trawl function with a parameter choice in
the long-memory regime in Tables LABEL:tab:Con-NB-LM-fixedt LABEL:tab:Con-
Gamma-LM-fixedt LABEL:tab:Con-Gauss-LM-fixedt, for fixed-$t$, and, for
fixed-$i$, in Tables LABEL:tab:Con-NB-LM LABEL:tab:Con-Gamma-LM LABEL:tab:Con-
Gauss-LM. Again, the results are generally very good, apart from one
particular setting: When $\Delta_{\\_}n=0.1$ and $t=0$ (which is equivalent to
$i=0$), we note that $\hat{a}(0)$ is strongly biased. This effect appears both
when using the recommended realised variance-based estimator and the
estimator-type used for time points $t>0$. The bias appears for all three
marginal distributions considered. That said, the bias-corrected estimator
$\hat{a}(0)-\frac{1}{2}\Delta_{\\_}n\hat{a}(0)$ performs noticeably better in
this particular setting. This rather poor performance goes away when
$\Delta_{\\_}n=0.01$ is considered instead.

We note that the choice of $n$ has almost no impact on the size of the bias,
but an increase in $n$ decreases the standard deviation of the estimators.

Overall, we conclude that the proposed estimators work well in finite samples
across different marginal distributions, short- and long-memory settings and
different in-fill and long-span asymptotic settings. In the long-memory case,
we advise the use of the bias-corrected estimator, in particular for small
times $t$.

####  D.1 Fixed-$t$

In the following, we will be presenting simulation results where we estimate
the trawl function $a$ at fixed time points $t$.

#####  D.1.1 Exponential trawl function

#####  D.1.2 SupGamma trawl function

####  D.2 Fixed-$i$

In the following, we will be presenting simulation results where we estimate
the trawl function $a$ at time points $t=i\Delta_{\\_}n$ for fixed values of
$i$.

#####  D.2.1 Exponential trawl function

#####  D.2.2 SupGamma trawl function

###  Appendix E Simulation study: Asymptotic normality results

Next, we investigate the finite sample behaviour of the central limit
theorems. We consider three statistics and report their mean, standard
deviation and coverage probabilities at levels 75%, 80%, 85%, 90%, 95% and 99%
based on 200 Monte Carlo samples:

Infeasible statistic:

The infeasible statistic is given by

$T^{IF}(t):=\sqrt{n\Delta_{\\_}{n}}\frac{\left(\hat{a}(t)-a(t)\right)}{\sqrt{\sigma_{\\_}a^{2}(t)}},$

where the true (theoretical) variance $\sigma_{\\_}a^{2}(t)$ is used.

Feasible statistic without bias correction:

The feasible statistic without bias correction is given by

$T^{F}:=\sqrt{n\Delta_{\\_}{n}}\frac{\left(\hat{a}(t)-a(t)\right)}{\sqrt{\hat{\sigma}_{\\_}a^{2}(t)}},$

where $\hat{\sigma}_{\\_}a^{2}(t)$ is the estimator proposed in the main
article.

Feasible statistic with bias correction:

The feasible statistic with bias correction is given by

$T^{F}_{\\_}{\mathrm{bias\,corrected}}:=\sqrt{n\Delta_{\\_}{n}}\frac{\left(\hat{a}(t)-a(t)-\frac{1}{2}\Delta_{\\_}n\hat{a}^{\prime}(t)\right)}{\sqrt{\hat{\sigma}_{\\_}a^{2}(t)}},$

where the estimator $\hat{a}$ proposed in the main article is used.

Based on our asymptotic theory, we expect the three statistics to follow a
standard normal distribution for large enough $n$ and small enough
$\Delta_{\\_}n$.

There is one case, however, which needs to be treated separately: In the
purely Gaussian case with $t=0$, we multiply the test statistics by $\sqrt{3}$
to reflect the fact that they would otherwise follow an $N(0,1/3)$
distribution. Note also that in that particular case, there is no infeasible
statistic corresponding to the central limit theorems we have presented in the
main article and hence the cell in the table is left empty.

######  Remark E.1.

In the purely Gaussian case, when $t=0$, we could consider the following
infeasible statistic

$\displaystyle\sqrt{n}\frac{(\hat{a}(0)-a(0))}{\sqrt{2(a(0))^{2}}},$

which asymptotically follows a standard normal distribution under the
assumptions of Theorem 3 (Part 1.) when $\mu=0$ and $\mu_{\\_}0=0$.

####  E.1 Fixed-$t$

#####  E.1.1 Exponential trawl function

We start off with the negative binomial trawl process with an exponential
trawl function, see Table LABEL:tab:NB-Exp-fixedt. For the larger value of
$\Delta_{\\_}n=0.1$, we observe a noticeable finite sample bias for all three
statistics, in particular, when $t$ is small, but at the same time, the
coverage probabilities are reasonable. The best results are obtained for the
bias-corrected statistic. We note, however, that the estimate of the
asymptotic variance gets less reliable when $t$ increases. When moving to the
case of the smaller $\Delta_{\\_}n=0.01$, the finite sample bias is reduced.
The problem that the estimate of the asymptotic variance is slightly less
reliable for increasing $t$ persists. We have to keep in mind that, when
decreasing $\Delta_{\\_}n$, while keeping $n$ fixed as in the set-up here,
$\Delta_{\\_}nn$ decreases, so we might not necessarily get an improvement in
the finite sample performance since the performance will depend on the balance
between the infill and long-span asymptotic regime.

When moving to the case of a Gamma marginal distribution with Exponential
trawl, see Table LABEL:tab:Gamma-Exp-fixedt, the results are similar, but
generally slightly worse than in the negative binomial case. Here we note that
the bias correction seems to be more useful for the case when
$\Delta_{\\_}n=0.1$ than for $\Delta_{\\_}n=0.01$. Also, when $\Delta_{\\_}nn$
is relatively small (e.g. for $\Delta_{\\_}n=0.1,n=2000$ or
$\Delta_{\\_}n=0.01,n\in\\{2000,5000\\}$) the finite sample performance of the
test statistic for large $t\in\\{5,10\\}$ is significantly worse than for
small $t$.

Moving on to the Gaussian case, see Table LABEL:tab:Gaussian-Exp-fixedt, we
find again that the bias-corrected statistic performs better, in particular
for small values of $t$. We also note that the case $t=0$ appears generally
tricky, where the best performance out of the cases considered are the cases
$\Delta_{\\_}n=0.1,n=2000$ and $\Delta_{\\_}n=0.01,n=2000$. It is important to
remember that this finding is not surprising given the additional (stronger)
assumptions we have in the central limit theorem for this particular case.

#####  E.1.2 SupGamma trawl function

Let us next study the finite sample performance in the long-memory regime, see
Tables LABEL:tab:NB-LM-fixedt, LABEL:tab:Gamma-LM-fixedt and
LABEL:tab:Gaussian-LM-fixedt. We recall that we found that, in this case,
there was a noticeable estimation bias for all three marginal distributions
when $t=0$. These results are reflected in our findings for the three
statistics. When $t=0$ and $\Delta_{\\_}n=0.1$, the finite sample results are
very poor. However, when $\Delta_{\\_}n=0.01$, the results are generally much
better and improve significantly as $n$ increases. We also recall the
observation from the exponential case that, when $t$ is increasing, the
estimates of the mean and standard deviation of the test statistic deviate
greatly from 0 and 1, respectively, when $\Delta_{\\_}nn$ is comparably small.
We observe the same pattern in the supGamma setting, where the findings are
even more pronounced. We conclude here that for larger values of $t$, the most
reliable performance of the test statistics is in the setting when
$\Delta_{\\_}n=0.1$ and $n\in\\{5000,10000\\}$.

####  E.2 Fixed-$i$

The finite sample performance of the three test statistics when reported for
fixed-$i$ is similar to as described above.

#####  E.2.1 Exponential trawl function

#####  E.2.2 SupGamma trawl function

###  Appendix F Simulation study: Consistency of the slice (ratio) estimation

As we have seen in the main article, for forecasting purposes it is of
interest to estimate slices of the trawl set. Hence we investigate in our
simulation study the finite sample performance of our proposed slice
estimators. We concentrate on the case which is relevant for our empirical
study and choose $\Delta_{\\_}n=0.1$ and $n=5000$. Recall that, in the main
article, we proposed the following estimators:

$\displaystyle\widehat{\mathrm{Leb}(A)}$
$\displaystyle=\sum_{\\_}{l=0}^{n-1}\hat{a}(l\Delta_{\\_}n)\Delta_{\\_}n,$
$\displaystyle\widehat{\mathrm{Leb}(A\cap A_{\\_}h)}$
$\displaystyle=\sum_{\\_}{\begin{subarray}{c}l\in\\{0,\ldots,n-1\\}:\\\
l\Delta_{\\_}n\geq h\end{subarray}}\hat{a}(l\Delta_{\\_}n)\Delta_{\\_}n,$
$\displaystyle\widehat{\mathrm{Leb}(A\setminus A_{\\_}h)}$
$\displaystyle=\sum_{\\_}{\begin{subarray}{c}l\in\\{0,\ldots,n-1\\}:\\\
l\Delta_{\\_}n<h\end{subarray}}\hat{a}(l\Delta_{\\_}n)\Delta_{\\_}n.$

We note that the estimates could become negative, so in our simulation
experiment, we made the following adjustment to the estimators:

$\displaystyle T_{\\_}h$
$\displaystyle:=\sum_{\\_}{\begin{subarray}{c}l\in\\{0,\ldots,n-1\\}:\\\
l\Delta_{\\_}n\geq h\end{subarray}}\hat{a}(l\Delta_{\\_}n)\Delta_{\\_}n$
$\displaystyle\widehat{\mathrm{Leb}(A\cap A_{\\_}h)}$
$\displaystyle=\left\\{\begin{array}[]{ll}T_{\\_}h,&\mathrm{if}\;T_{\\_}h<\widehat{\mathrm{Leb}(A)}\\\
\widehat{\mathrm{Leb}(A)},&\mathrm{otherwise}.\end{array}\right.$
$\displaystyle\widehat{\mathrm{Leb}(A\setminus A_{\\_}h)}$
$\displaystyle=\widehat{\mathrm{Leb}(A)}-T_{\\_}h.$

Based on these estimators, we estimate $\mathrm{\mathrm{Leb}}(A)$,
$\mathrm{\mathrm{Leb}}(A\cap A_{\\_}h)$, $\mathrm{\mathrm{Leb}}(A\setminus
A_{\\_}h)$, $\mathrm{\mathrm{Leb}}(A\cap A_{\\_}h)/\mathrm{\mathrm{Leb}}(A)$
and $\mathrm{\mathrm{Leb}}(A\setminus A_{\\_}h)/\mathrm{\mathrm{Leb}}(A)$. For
each of the five estimators, we report their empirical mean, bias and standard
deviation based on 200 Monte Carlo samples. We observe that the slices and
corresponding ratios can be estimated with high precision for all three
marginal distributions and for both the short-memory and the long-memory
settings, see Tables LABEL:tab:slices-exp and LABEL:tab:slices-lm,
respectively.

When estimating the slices, we also implemented the estimator using a bias
correction, but find that the bias-corrected estimator performs worse when
estimating the slices and ratios.

As an alternative estimator, we use

$\displaystyle\widehat{\mathrm{Leb}(A)}$
$\displaystyle=\widehat{\mathrm{Var}(X)}$ $\displaystyle S_{\\_}h$
$\displaystyle:=\widehat{\mathrm{Cov}(X_{\\_}0,X_{\\_}h)},$
$\displaystyle\widehat{\mathrm{Leb}(A\cap A_{\\_}h)}$
$\displaystyle=\left\\{\begin{array}[]{ll}S_{\\_}h,&\mathrm{if}\;S_{\\_}h<\widehat{\mathrm{Leb}(A)}\\\
\widehat{\mathrm{Leb}(A)},&\mathrm{otherwise}.\end{array}\right.$
$\displaystyle\widehat{\mathrm{Leb}(A\setminus A_{\\_}h)}$
$\displaystyle=\widehat{\mathrm{Leb}(A)}-S_{\\_}h,$

where $\widehat{\mathrm{Var}(X)}$ and
$\widehat{\mathrm{Cov}(X_{\\_}0,X_{\\_}h)}$ denote the empirical variance and
covariance, respectively.

Based on these alternative estimators, we estimate $\mathrm{\mathrm{Leb}}(A)$,
$\mathrm{\mathrm{Leb}}(A\cap A_{\\_}h)$, $\mathrm{\mathrm{Leb}}(A\setminus
A_{\\_}h)$, $\mathrm{\mathrm{Leb}}(A\cap A_{\\_}h)/\mathrm{\mathrm{Leb}}(A)$
and $\mathrm{\mathrm{Leb}}(A\setminus A_{\\_}h)/\mathrm{\mathrm{Leb}}(A)$. For
each of the five estimators, we report their empirical mean, bias and standard
deviation based on 200 Monte Carlo samples. We observe that the slices and
corresponding ratios can be estimated with high precision for all three
marginal distributions and for both the short-memory and the long-memory
setting, see Tables LABEL:tab:slices-exp-acfbased and LABEL:tab:slices-lm-
acfbased, respectively. In particular, it appears that these empirical acf-
based estimators outperform the estimator based on summing up the individual
trawl function estimates.

### References

  * [1] [author] Bacro, J.J., Gaetan, C.C., Opitz, T.T. and Toulemonde, G.G. (2020). Hierarchical Space-Time Modeling of Asymptotically Independent Exceedances With an Application to Precipitation Data. Journal of the American Statistical Association 115 555-569. 10.1080/01621459.2019.1617152
  * [2] [author] Barndorff-Nielsen, O. E.O. E. (2011). Stationary infinitely divisible processes. Brazilian Journal of Probability and Statistics 25 294–322. 10.1214/11-BJPS140
  * [3] [author] Barndorff-Nielsen, Ole E.O. E., Benth, Fred EspenF. E. and Veraart, Almut E. D.A. E. D. (2018). Ambit stochastics. Probability Theory and Stochastic Modelling 88\. Springer, Cham. 10.1007/978-3-319-94129-5 3839270
  * [4] [author] Barndorff-Nielsen, Ole E.O. E., Lunde, AsgerA., Shephard, NeilN. and Veraart, AlmutA. (2014). Integer-valued Trawl Processes: A Class of Stationary Infinitely Divisible Processes. Scandinavian Journal of Statistics 41 693-724. https://doi.org/10.1111/sjos.12056
  * [5] [author] Basse-O’Connor, AndreasA. and Rosiński, JanJ. (2016). On infinitely divisible semimartingales. Probability Theory and Related Fields 164 133–163.
  * [6] [author] Bennedsen, MikkelM., Lunde, AsgerA., Shephard, NeilN. and Veraart, AlmutA. (2022). Inference and forecasting for continuous-time integer-valued trawl processes and their use in financial economics. Preprint. 10.48550/ARXIV.2107.03674
  * [7] [author] Bennedsen, MikkelM., Lunde, AsgerA., Shephard, NeilN. and Veraart, AlmutA. (2022). Inference and forecasting for continuous-time integer-valued trawl processes and their use in financial economics: Supplementary material. Preprint. 10.48550/ARXIV.2107.03674
  * [8] [author] Bollen, Nicolas P. B.N. P. B., Smith, TomT. and Whaley, Robert E.R. E. (2004). Modeling the bid/ask spread: measuring the inventory-holding premium. Journal of Financial Economics 72 97–141. https://doi.org/10.1016/S0304-405X(03)00169-7
  * [9] [author] Courgeau, ValentinV. and Veraart, Almut E. D.A. E. D. (2022). Asymptotic theory for the inference of the latent trawl model for extreme values. Scandinavian Journal of Statistics n/a. To appear. https://doi.org/10.1111/sjos.12563
  * [10] [author] Diebold, Francis X.F. X. and Mariano, Roberto S.R. S. (1995). Comparing Predictive Accuracy. Journal of Business & Economic Statistics 13 253-263. 10.1080/07350015.1995.10524599
  * [11] [author] Doukhan, P.P., Jakubowski, A.A., Lopes, S. R. C.S. R. C. and Surgailis, D.D. (2019). Discrete-time trawl processes. Stochastic Processes and their Applications 129 1326-1348. https://doi.org/10.1016/j.spa.2018.05.004
  * [12] [author] Doukhan, P.P., Roueff, F.F. and Rynkiewicz, J.J. (2020). Spectral estimation for non-linear long range dependent discrete time trawl processes. Electronic Journal of Statistics 14 3157-3191. 10.1214/20-EJS1742
  * [13] [author] Grahovac, D.D., Leonenko, N. N.N. N. and Taqqu, M. S.M. S. (2018). Intermittency of trawl processes. Statistics and Probability Letters 137 235 - 242. https://doi.org/10.1016/j.spl.2018.01.030
  * [14] [author] Hall, P.P. and Heyde, C. C.C. C. (1980). Martingale Limit Theory and Its Application. Probability and mathematical statistics. Academic Press.
  * [15] [author] Huang, Roger D.R. D. and Stoll, Hans R.H. R. (1997). The Components of the Bid-Ask Spread: A General Approach. The Review of Financial Studies 10 995–1034.
  * [16] [author] Jacod, J.J. and Protter, P.P. (2011). Discretization of Processes. Stochastic Modelling and Applied Probability. Springer, Berlin, Heidelberg.
  * [17] [author] Leonte, DanD. (2022). Simulation and estimation of trawl processes (working title), PhD thesis, Imperial College London Chapter: Forecasting trawl processes, in preparation.
  * [18] [author] Lin, Z.Z. and Bai, Z.Z. (2010). Probability Inequalities. Springer, Berlin, Heidelberg.
  * [19] [author] Luschgy, HaraldH. and Pagès, GillesG. (2008). Moment estimates for Lévy Processes. Electronic Communications in Probability 13 422 – 434. 10.1214/ECP.v13-1397
  * [20] [author] Noven, R. C.R. C., Veraart, A. E. D.A. E. D. and Gandy, A.A. (2018). A latent trawl process model for extreme values. Journal of Energy Marktes 11 1-24.
  * [21] [author] Pakkanen, Mikko S.M. S., Passeggeri, RiccardoR., Sauri, OrimarO. and Veraart, Almut E. D.A. E. D. (2021). Limit theorems for trawl processes. Electronic Journal of Probability 26 1 – 36. 10.1214/21-EJP652
  * [22] [author] Paulauskas, V.V. (2020). A note on linear processes with tapered innovations. Lithuanian Mathematical Journal 60 64-79.
  * [23] [author] Rosiński, JanJ. and Żak, TomaszT. (1996). Simple conditions for mixing of infinitely divisible processes. Stochastic Processes and their Applications 61 277–288. https://doi.org/10.1016/0304-4149(95)00083-6
  * [24] [author] Sauri, O.O. and Veraart, A. E. D.A. E. D. (2022). Supplementary Appendix to Nonparametric estimation of trawl processes: Theory and Applications. Supplementary material.
  * [25] [author] Sauri, OrimarO. and Veraart, Almut E. D.A. E. D. (2022). Nonparametric estimation of trawl processes: Theory and Applications. Preprint.
  * [26] [author] Shephard, NeilN. and Yang, Justin J.J. J. (2016). Likelihood Inference for Exponential-Trawl Processes. In The Fascination of Probability, Statistics and their Applications: In Honour of Ole E. Barndorff-Nielsen (MarkM. Podolskij, RobertR. Stelzer, SteenS. Thorbjørnsen and Almut E. D.A. E. D. Veraart, eds.) 251–281. Springer International Publishing, Cham. 10.1007/978-3-319-25826-312
  * [27] [author] Shephard, N.N. and Yang, J. J.J. J. (2017). Continuous Time Analysis of Fleeting Discrete Price Moves. Journal of the American Statistical Association 112 1090-1106.
  * [28] [author] Talarczyk, AA. and Treszczotko, L.L. (2020). Limit theorems for integrated trawl processes with symmetric Lévy bases. Electronic Journal of Probability 25 1 – 24. 10.1214/20-EJP509
  * [29] [author] Turner, Matthew DM. D. (2011). Explicit $L_{\\_}p$-norm estimates of infinitely divisible random vectors in Hilbert spaces with applications, PhD thesis, University of Tennessee.
  * [30] [author] Veraart, A. E. D.A. E. D. (2019). Modeling, simulation and inference for multivariate time series of counts using trawl processes. Journal of Multivariate Analysis 169 110 - 129.
  * [31] [author] Veraart, Almut E. D.A. E. D. (2022). ambit: Simulation and Estimation of Ambit Processes R package version 0.1.2.
  * [32] [author] Veraart, Almut E. D.A. E. D. (2022). Nonparametric-trawl-estimation R code, release v1.0.0.
  * [33] [author] Wolpert, Robert L.R. L. and Brown, Lawrence D.L. D. (2011). Stationary Infinitely-Divisible Markov Processes with Non-negative Integer Values. Working paper, April 2011.
  * [34] [author] Wolpert, Robert L.R. L. and Taqqu, Murad S.M. S. (2005). Fractional Ornstein–Uhlenbeck Lévy processes and the Telecom process: Upstairs and downstairs. Signal Processing 85 1523–1545.


# Persistence of spectral projections for stochastic operators on large tensor
products

R.S.MacKay  Mathematics Institute, University of Warwick, Coventry CV4 7AL,
U.K<EMAIL_ADDRESS>
###### Abstract.

In this paper it is proved that for families of stochastic operators on a
countable tensor product, depending smoothly on parameters, any spectral
projection persists smoothly, where smoothness is defined using norms based on
ideas of Dobrushin. A rigorous perturbation theory for families of stochastic
operators with spectral gap is thereby created. It is illustrated by deriving
an effective slow 2-state dynamics for a 3-state probabilistic cellular
automaton. Some further potential applications are discussed.

###### Key words and phrases:

Stochastic operator, Probabilistic cellular automaton, Spectral projection

###### 2010 Mathematics Subject Classification:

60J05, 60G60

##  1\. Introduction

The problem of persistence of spectral projections for large tensor products
is crucial in many domains. Perhaps the most significant domain is many-
particle quantum systems. In addition to condensed matter physics, this has
taken on enhanced interest because of the problem of designing quantum
registers for quantum computing.

Yet a parallel problem arises for stochastic systems with many components.
This paper specialises to Markov processes and mostly to discrete time, such
as probabilistic cellular automata (PCA). Then the transition operator $T$
acts on the space $F$ of real-valued continuous functions of the state of the
whole system. The space $F$ can be considered to be the tensor product of the
spaces of real-valued continuous functions of the state of the individual
units.

Even if the system is geometrically ergodic (meaning there is a unique
stationary probability and it attracts every probability exponentially in an
appropriate metric) and the update of each unit is independent of the state
outside a bounded neighbourhood, when one changes parameters in a reasonable
way the stationary probability may move at a speed going to infinity with the
number of units if distances between probability distributions are measured in
any of the standard ways (e.g. total variation, Jensen-Shannon, Hellinger,
Kantorovich, Fisher information [M1], and Prokhorov [M2]). The solution
proposed in [M1] was to introduce a new111After publication of [M1], I found
that Steif had proposed a solution 20 years earlier [Sf]. Although his metric
is defined differently and in a restricted context, I was able to prove that
it is in fact equal to mine in the finite case (Appendix to [M2]), and with
Armstrong-Goodall we have now proved they are equal under the general
conditions for definition of mine [AM]. metric for probabilities on large
product systems, christened “Dobrushin metric” as most of the ingredients were
already in Dobrushin’s work, but credit should also be given to Vasershtein
[Vn] (more commonly transliterated now as Wasserstein). With respect to this
metric, smooth variation of the stationary probability was proved for the
class of geometrically ergodic PCA, uniformly in the size of the system [M1].
A slightly more sophisticated way of viewing this result is as persistence of
the (rank-1) projection $P$ onto the space of stationary measures, for which
complementary projection $Q=I-P$ sends constant functions to zero. Thus $P$ is
a spectral projection, a projection operator onto a subspace corresponding to
a closed subset of the spectrum of $T$, whose complementary projection $Q=I-P$
is onto a complementary subspace corresponding to the disjoint closed
complement of the spectrum.

A question that the work of [M1] prompted is whether other spectral
projections for stochastic operators might also persist uniformly smoothly in
the size of the system. Here suitable conditions are formulated and a proof of
persistence is given.

Borel measures $p$ are dual to continuous functions $f$, in that one can take
$p(f)$ to be the integral of $f$ with respect to $p$. We think of functions as
column vectors and measures as row vectors. A transition operator $T$ acts to
the right on functions and to the left on measures. Transition operators
preserve total probability, which can be written as preservation of the
function ${\bf 1}$, defined to take the value $1$ everywhere.

The outline of the paper is that firstly, Dobrushin metric is reviewed
(sec.2). Then the continuation problem for spectral projections of a class of
stochastic operators is formulated and solved (sec.3). An illustration of the
result is given (sec.4), followed by a general development of second-order
perturbation theory for stochastic operators (sec.5) and then a discussion of
further potential applications, including to metastability (sec.6). The paper
ends with a short discussion (sec.7).

##  2\. Dobrushin metric

Consider transition operators $T$ for functions and probability distributions
on the product $X$ of a set of metric spaces $(X_{s},d_{s})$, for sites $s$ in
a countable set $S$. The spaces $(X_{s},d_{s})$ are assumed to be Polish
(complete separable metric spaces) with bounded diameter, $\sup_{s\in
S}\mbox{diam}_{s}(X_{s})<\infty$. The product $X$ is endowed with product
topology and with Borel measures. The set of Borel probabilities on $X$
(measures $p$ satisfying $p(X)=1$ and $p(Y)\geq 0$ for all Borel subsets $Y$)
is denoted by ${\mathcal{P}}$.

For a function $f:X\to\mathbb{R}$, define its Lipschitz constant with respect
to variations on site $s\in S$ by

$\Delta_{s}(f)=\sup\frac{f(x)-f(x^{\prime})}{d_{s}(x_{s},x^{\prime}_{s})}$

over $x,x^{\prime}\in X$ differing at site $s$ and agreeing elsewhere. Define
the set $F$ of Dobrushin smooth functions to be those for which the semi-norm

$|f|_{F}=\sum_{s\in S}\Delta_{s}(f)$

is finite. Define the space $Z$ of zero-charge measures on $X$ to be the
signed Borel measures $\mu$ for which $\mu(X)=0$. Define a norm

$|\mu|_{Z}=\sup_{f\in F\setminus C}\frac{\mu(f)}{|f|_{F}}$

on $Z$, where $C$ is the set of constant functions. Then define the Dobrushin
distance between any two $\rho,\sigma\in{\mathcal{P}}$ by

$D(\rho,\sigma)=|\rho-\sigma|_{Z}.$

This makes ${\mathcal{P}}$ into a complete metric space, with diameter
$\sup_{s\in S}\mbox{diam}_{s}(X_{s})$ [M1].

The point of this metric is that if $(I-T)$ is invertible on $Z$ then $T$ has
a unique stationary probability $p\in\mathcal{P}$ and it varies smoothly with
respect to changes in $T$:

$p^{\prime}=pT^{\prime}(I-T)_{Z}^{-1},$

where ′ denotes derivative with respect to parameters and $(I-T)_{Z}$ means
the restriction of $I-T$ to $Z$ (which is necessary to take its inverse) [M1].
Also there are conditions for invertibility of $I-T$ in terms of Dobrushin’s
dependency matrix, which are verifiable in relevant classes of system.

##  3\. Spectral projections

Given a bounded transition operator $T$ on a space $F$ of continuous
functions, a spectral projection for $T$ is a bounded operator $P$ on $F$ such
that $P^{2}=P$, $PT=TP$ and the parts of the spectrum of $T$ corresponding to
the restrictions of $T$ to the image222commonly called the “range”, but
“range” is also used for the whole space into which $P$ maps. of $P$ and the
image of $Q=I-P$ (which are invariant under $T$) are disjoint. Because they
are closed and bounded, it follows that the distance between them is positive,
called a spectral gap.

An example of a spectral projection is $P={\bf 1}p$ for the stationary
probability $p$ for a geometrically ergodic operator $T$. This is because
$p{\bf 1}=1$ for any probability $p$, $T{\bf 1}={\bf 1}$ by definition of a
stochastic operator, $pT=p$ for a stationary measure, and geometric ergodicity
implies the eigenspace for eigenvalue 1 is one-dimensional and the rest of the
spectrum is in a disk of radius less than 1.

To extend the persistence theory to more general spectral projections, one
needs to define norms on changes to transition operators and on tangent
vectors to the manifold $M$ of projections.

Note that the set $M$ of projections on $F$ is indeed a manifold (a variant of
a Grassmann manifold). Here is an outline proof, because ingredients are
useful later. Let $B(F)$ be the space of bounded linear operators $P$ on the
space $F$ and define $\Phi$ on $B(F)$ by

(1) $\Phi(P)=P^{2}-P.$

So $M=\Phi^{-1}(0)$. Given $P\in M$, let $R=\mbox{im}\ P$ and $K=\ker
P=\mbox{im}\ Q$. Both are closed. Then $F=R\oplus K$. Relative to this direct-
sum decomposition,

(2) $P=\left[\begin{array}[]{cc}I&0\\\ 0&0\end{array}\right].$

For an arbitrary operator $\pi=\left[\begin{array}[]{cc}\pi_{1}&\pi_{2}\\\
\pi_{3}&\pi_{4}\end{array}\right]$ on $F$,

(3) $D\Phi_{P}(\pi)=\left[\begin{array}[]{cc}\pi_{1}&0\\\
0&-\pi_{4}\end{array}\right].$

With respect to the same direct-sum decomposition, let $T_{P}M$ be333Note the
re-use of the symbol $T$, which will signify tangent space as well as the
transition operator. the space of operators on $F$ of the form

(4) $\left[\begin{array}[]{cc}0&\pi_{2}\\\ \pi_{3}&0\end{array}\right],$

and $N_{P}M$ be those of the form $\left[\begin{array}[]{cc}\pi_{1}&0\\\
0&\pi_{4}\end{array}\right]$. Consequently, $D\Phi_{P}$ maps $N_{P}M$ to
itself by $\left[\begin{array}[]{cc}I&0\\\ 0&-I\end{array}\right]$, which has
bounded inverse, namely itself. So the implicit function theorem shows that
$M$ is locally the graph of a $C^{1}$ function $\psi:T_{P}M\to N_{P}M$. This
completes the outline proof. Furthermore, $D\psi=0$, so $T_{P}M$ is the
tangent space to $M$ at $P$ (hence the notation). Note, however, that the
manifold $M$ of projections has many components of different dimensions,
corresponding to the rank of the projection.

Now the paper moves to the promised definition of norms on changes to
transition operators and on tangents to $M$.

Firstly, a change $T^{\prime}$ to a transition operator satisfies
$T^{\prime}{\bf 1}={\bf 0}$, where ${\bf 0}$ is the function taking value $0$
everywhere, so $T^{\prime}$ takes any measure into $Z$. One can quantify the
size of its effect on $Z$ by the operator norm of its restriction to $Z$:

$|T^{\prime}_{Z}|_{Z}=\sup_{\mu\in Z\setminus 0}\frac{|\mu
T^{\prime}|_{Z}}{|\mu|_{Z}}.$

One needs also, however, to measure the size of $\rho T^{\prime}$ for measures
$\rho$ outside $Z$, for example the Borel probabilities $\mathcal{P}$, which
can be non-zero even if $\mu T^{\prime}=0$ for all $\mu\in Z$. A suitable
quantification of this is

$L(T^{\prime})=\sup_{p\in\mathcal{P}}|pT^{\prime}|_{Z}.$

These two quantities were used in [M1] to define continuous change of a
transition operator and smooth change of transition operator. So take the norm

(5) $|T^{\prime}|_{*}=|T^{\prime}_{Z}|_{Z}+L(T^{\prime}).$

It is a norm because each part is non-negative, homogeneous and satisfies the
triangle inequality, and $|T^{\prime}|_{*}=0$ implies $T^{\prime}=0$. To see
the latter, $|T^{\prime}|_{*}=0$ implies both $|T^{\prime}_{Z}|_{Z}$ and
$L(T^{\prime})$ are zero. $T^{\prime}$ maps any measure into $Z$, so for any
$\mu\in Z$ then $\mu T^{\prime}=0$ and for any $p\in\mathcal{P}$ then
$pT^{\prime}=0$. Now any measure $\nu$ can be written as $kp+\mu$ for some
$k\in\mathbb{R}$, $p\in\mathcal{P}$ and $\mu\in Z$, so $\nu T^{\prime}=0$,
thus $T^{\prime}=0$.

Secondly, the question of defining a norm on tangents to the manifold $M$ of
projections is addressed. Infinitesimal changes $P^{\prime}$ to a projection
$P$ do not necessarily map measures to $Z$. They are characterised by
$P^{\prime}P+PP^{\prime}=P^{\prime}$, which can equivalently be written as
$P^{\prime}P=QP^{\prime}$ or as $PP^{\prime}=P^{\prime}Q$ and hence in the
form (4) with respect to the direct sum decomposition $F=R\oplus K$.

Nevertheless, for spectral projections one can restrict attention to the
submanifolds $M_{1},M_{0}$ of projections $P$ that fix ${\bf 1}$ or send it to
${\bf 0}$, respectively, according as eigenvalue $1$ is or is not in the
spectrum of the restriction of $T$ to the image of $P$. In either case,
tangents satisfy the additional condition $P^{\prime}{\bf 1}={\bf 0}$. Thus
$P^{\prime}$ satisfies the same conditions as a change to a transition
operator above. So, use the norm $|P^{\prime}|_{*}$ on the tangent spaces to
$M_{0}$ and $M_{1}$. For the special case of $P^{\prime}={\bf 1}p^{\prime}$,
where $p^{\prime}$ is a change in a probability, $|P^{\prime}|_{*}$ is
precisely $|p^{\prime}|_{Z}$.

So now the question of persistence of a spectral projection reduces to
application of the implicit function theorem to the equation $PT=TP$ for $P$
in the manifold $M_{0}$ or $M_{1}$. The tangent space $TM_{i}$ to $M_{i}$
($i=0,1$) at $P$ consists of the operators $P^{\prime}$ such that
$PP^{\prime}+P^{\prime}P=P^{\prime}$ and $P^{\prime}{\bf 1}={\bf 0}$. Consider
the functions $G_{i}:M_{i}\times N\to TM_{i}$ where $N$ is the space (affine
with boundary) of transition operators $T$ taking non-negative functions to
non-negative ones and satisfying $T{\bf 1}={\bf 1}$, defined by the commutator

$G_{i}(P,T)=[P,T]=PT-TP.$

It is easily checked that the right hand side is in $TM_{i}$. $G_{i}$ is
jointly $C^{1}$ with derivatives

(6) $\displaystyle\frac{\partial G_{i}}{\partial P}(P^{\prime})$
$\displaystyle=$ $\displaystyle[P^{\prime},T],$ $\displaystyle\frac{\partial
G_{i}}{\partial T}(T^{\prime})$ $\displaystyle=$
$\displaystyle[P,T^{\prime}].$

Tangents $T^{\prime}$ to $N$ are characterised by $T^{\prime}{\bf 1}={\bf 0}$.
If $\frac{\partial G_{i}}{\partial P}:TM_{i}\to TN$ is invertible with bounded
inverse then a solution $P$ to $G_{i}(P,T)=0$ has locally unique continuation
$P(T)$ for nearby $T$ and depends $C^{1}$ on $T$, with

(7) $P^{\prime}=-\left(\frac{\partial G_{i}}{\partial
P}\right)^{-1}\frac{\partial G_{i}}{\partial T}T^{\prime}.$

The condition of bounded inverse boils down to a spectral gap. To see this,
take a direct-sum decomposition in which

(8) $P=\left[\begin{array}[]{cc}I&0\\\ 0&0\end{array}\right],\
T=\left[\begin{array}[]{cc}T_{P}&0\\\ 0&T_{Q}\end{array}\right].$

Then as in (4), $P^{\prime}$ has the form

(9) $P^{\prime}=\left[\begin{array}[]{cc}0&U\\\ V&0\end{array}\right].$

So

(10) $\frac{\partial G}{\partial
P}(P^{\prime})=\left[\begin{array}[]{cc}0&UT_{Q}-T_{P}U\\\
VT_{P}-T_{Q}V&0\end{array}\right],$

and for (7) one wants to be able to solve $\frac{\partial G}{\partial
P}(P^{\prime})$ equal to an arbitrary element of $TM_{i}$ for $P^{\prime}$.
The elements of $TM_{i}$ have the same zero-diagonal block form.

Thus this reduces to solving two “Sylvester equations”

$AX-XB=C$

for operators $X=U,V$ respectively, with $A,B$ equal to $T_{P},T_{Q}$ or vice
versa, and arbitrary bounded $C$. There is a unique bounded solution $X$ to
each of the Sylvester equations iff the spectra of $T_{P}$ and $T_{Q}$ are
disjoint [BR]. In the disjoint case, it is automatic that the resulting
inverse operator $C\mapsto X$ is bounded [K].

For practical purposes, one needs a bound on the inverse. Define

$\mbox{sep}(A,B)=\inf_{X\neq 0}\frac{\|AX-XB\|}{\|X\|}.$

This generalises the definition of the “separation” of two matrices $A,B$,
which is usually defined using Frobenius norm [Sw, Vh]. Note that
$\mbox{sep}(B,A)$ is not necessarily equal to $\mbox{sep}(A,B)$. From the
above discussion, $\mbox{sep}(A,B)>0$ iff the spectra of $A$ and $B$ are
disjoint.

In some cases one can be explicit about its size. For example, if the spectral
radii $\rho_{A},\rho_{B^{-1}}$ of $A$ and $B^{-1}$ satisfy
$\rho_{A}\rho_{B^{-1}}<1$ then for $\lambda>1$ there exist $C_{A},C_{B^{-1}}$
(dependent on $\lambda$) such that $\|A^{n}\|\leq C_{A}(\lambda\rho_{A})^{n}$
and $\|B^{-n}\|\leq C_{B^{-1}}(\lambda\rho_{B^{-1}})^{n}$. The explicit
(convergent) solution

$X=-\sum_{n=0}^{\infty}A^{n}CB^{-n-1}$

of the Sylvester equation shows that

$\mbox{sep}(A,B)\geq\frac{1-\lambda^{2}\rho_{A}\rho_{B^{-1}}}{C_{A}C_{B^{-1}}\lambda\rho_{B^{-1}}}.$

A method to obtain similar bounds for some more general forms of separation of
the spectra is described in [N].

The important point is that for a family of examples with growing system size
$N$, if the separation is bounded away from zero then the continuation above
is uniform in $N$. Furthermore, it can apply to infinite systems.

The same can be done in continuous time, with the transition operator replaced
by a transition generator, but I do not spell it out here, save to mention
that if the spectra of $A$ and $B$ lie in $\Re z\geq r_{A}$, $\Re z\leq r_{B}$
respectively, with $r_{A}>r_{B}$, the explicit (convergent) solution

$X=\int_{0}^{\infty}e^{(-A+r)t}Ce^{(B-r)t}dt$

for $r\in(r_{B},r_{A})$ provides a bound of the form

$\mbox{sep}(A,B)\geq\frac{r_{A}-r_{B}-2\varepsilon}{c_{A}c_{B}}$

for any $\varepsilon>0$, with $c_{A},c_{B}$ depending on $\varepsilon$.

##  4\. An illustration

As an example, consider a $3$-state PCA, with local state space $\\{+,0,-\\}$
at each site of a finite undirected graph with $N$ nodes and bounded degree,
say by $m$. The transition probabilities for the state at a given site are
taken to be

(11) $\left[\begin{array}[]{ccc}1-(1+\alpha n_{-})\varepsilon&(1+\alpha
n_{-})\varepsilon&0\\\ 1/2&0&1/2\\\ 0&(1+\alpha n_{+})\varepsilon&1-(1+\alpha
n_{+})\varepsilon\end{array}\right],$

in basis $(+,0,-)$, where $n_{\pm}$ are the numbers of neighbours in states
$\pm$ respectively, $\alpha\geq 0$ and $\varepsilon\in[0,(1+\alpha m)^{-1}]$.
For $\varepsilon=0$ the system has eigenvalue $1$ with multiplicity $2^{N}$
and eigenvalue $0$ with multiplicity $3^{N}-2^{N}$.

From the above theory, the spectral projection to the subspace for eigenvalue
$1$ has a continuation for $\varepsilon$ small. It is uniform in $N$ because
the separation of the relevant operators can be bounded away from $0$
uniformly in $N$. The dynamics on the image of the projection is slow because
the spectrum moves from $\\{1\\}$ by at most $O(\varepsilon)$. It is still
Markovian. It might loosely be considered as an effective PCA with two states
$\\{+,-\\}$ on each site but $R=\mbox{im}\ P$ is not spanned by
$\delta$-functions on states, so such a description requires interpretation
analogous to quasiparticles in quantum mechanics.

If $\alpha<1/m$ and $\varepsilon>0$, it is geometrically ergodic, because the
whole system is. This can be proved using Dobrushin’s dependency matrix, as
follows. Take discrete metric on the local state spaces. Then the update
probability distributions $p(\sigma)$ for the state at a site given its
current state $\sigma$ and the numbers $n_{\pm}$ are the rows of the matrix
(11) and so the variation distance for a change of state on the given site is
at most $(1-\varepsilon)$ and for a change on a neighbouring site is at most
$\alpha\varepsilon$. Thus the dependency matrix has $\ell_{\infty}$-norm at
most $1-(1-m\alpha)\varepsilon$. If $\alpha<1/m$ and $\varepsilon>0$, this is
less than $1$ and so the system is geometrically ergodic.

It is interesting to compute approximate dynamics on the image of the
projection corresponding to the spectrum near $1$. This is analogous to a
second-order perturbation theory computation in quantum mechanics, e.g. the
derivation of the $t-J$ model from the Hubbard model [Sp]. Here, the $0$ state
mediates interactions between the $\pm$ states. A general treatment of second-
order perturbation theory for families of stochastic operators is given in the
next section.

##  5\. Second-order perturbation theory

Under the spectral gap condition at $\varepsilon=0$ and with the norm (5), it
has been proved above that for a smooth family of stochastic operators
$T_{\varepsilon}$ and a spectral projection $P_{0}$ for $T_{0}$, there is a
range of $\varepsilon$ for which $P_{0}$ continues smoothly to a spectral
projection $P_{\varepsilon}$ for $T_{\varepsilon}$. It is convenient to write

$P_{\varepsilon}=\psi_{\varepsilon}P_{0}{\psi_{\varepsilon}}^{-1}$

for an $\varepsilon$-dependent invertible bounded linear map $\psi$, with
$\psi_{0}$ the identity. There is a lot of freedom in the choice of $\psi$,
the only constraints being that
$\psi_{\varepsilon}(\mbox{im}P_{\varepsilon})=\mbox{im}P_{0}$ and
$\psi_{\varepsilon}(\ker P_{\varepsilon})=\ker P_{0}$, but it is good also to
take $\psi$ to be as smooth in $\varepsilon$ as is $T$. Then one can define

$S=\psi^{\prime}\psi^{-1},$

where ′ denotes $\frac{d}{d\varepsilon}$, and the above conditions reduce to

(12) $[S,P]=P^{\prime},$

where $[,]$ again denotes the commutator. A convenient solution is

$S=P^{\prime}(P-Q),$

which can be checked to satisfy the condition (12).

Then it is desired to compute $\hat{T}=\psi^{-1}T\psi$ on $\mbox{im}\ P_{0}$.
This is the stochastic operator that represents the dynamics of
$T_{\varepsilon}$ on $\mbox{im}\ P_{\varepsilon}$, using the coordinate system
$\psi_{\varepsilon}$.

Start from $\hat{T}_{0}=T_{0}$. The first derivative
$\hat{T}^{\prime}=\psi^{-1}(T^{\prime}+[T,S])\psi$. Using the above choice of
$S$, one obtains $\hat{T}^{\prime}=\psi^{-1}J\psi$, where
$J=PT^{\prime}P+QT^{\prime}Q$. The second derivative can be evaluated to
$\hat{T}^{\prime\prime}=\psi^{-1}(PT^{\prime\prime}P+QT^{\prime\prime}Q+[[T,P^{\prime}],P^{\prime}])\psi$.

Evaluating these at $\varepsilon=0$, one obtains $\hat{T}$ to second order in
$\varepsilon$ as

(13)
$\hat{T}_{\varepsilon}=P_{0}T_{\varepsilon}P_{0}+Q_{0}T_{\varepsilon}Q_{0}+\frac{\varepsilon^{2}}{2}[[T_{0},P_{0}^{\prime}],P_{0}^{\prime}].$

It is perhaps more useful to substitute $[T,P^{\prime}]=[T^{\prime},P]$, since
the righthand side of this is readily computable, but the second occurrence of
$P^{\prime}$ above means that solution of this equation for $P^{\prime}$ is
unavoidable.

For the example of section 4, there is already an effect at first order in
$\varepsilon$, so second-order perturbation theory is perhaps unnecessary, but
it still serves as an illustration of the procedure.

$T_{0}=P_{0}=\otimes_{s\in S}\left[\begin{array}[]{ccc}1&0&0\\\
\frac{1}{2}&0&\frac{1}{2}\\\
0&0&1\end{array}\right],T_{0}^{\prime}=\otimes_{s\in
S}\left[\begin{array}[]{ccc}-\beta_{-}&\beta_{-}&0\\\ 0&0&0\\\
0&\beta_{+}&-\beta_{+}\end{array}\right],$

where $\beta_{\pm}=1+\alpha n_{\pm}$. Thus

$P_{0}T_{\varepsilon}P_{0}+Q_{0}T_{\varepsilon}Q_{0}=\otimes_{s\in
S}\left[\begin{array}[]{ccc}1-\beta_{-}\varepsilon/2&0&\beta_{-}\varepsilon/2\\\
\tfrac{1}{2}+\beta_{+}\varepsilon/2&(\beta_{+}+\beta_{-})\varepsilon/2&\tfrac{1}{2}+\beta_{-}\varepsilon/2\\\
\beta_{+}\varepsilon/2&0&1-\beta_{+}\varepsilon/2\end{array}\right].$

To compute the second-order term, begin with

$[T_{0}^{\prime},P_{0}]=\otimes_{s\in
S}\left[\begin{array}[]{ccc}\frac{1}{2}\beta_{-}&-\beta_{-}&\frac{1}{2}\beta_{-}\\\
\frac{1}{2}\beta_{-}&-\frac{1}{2}(\beta_{-}+\beta_{+})&\frac{1}{2}\beta_{+}\\\
\frac{1}{2}\beta_{+}&-\beta_{+}&\frac{1}{2}\beta_{+}\end{array}\right].$

Using $P_{0}P_{0}^{\prime}=P_{0}^{\prime}Q_{0}$, $P_{0}^{\prime}$ has just
four independent parameters on each site, and solving
$[T,P^{\prime}]=[T^{\prime},P]$ for them yields

$P_{0}^{\prime}=\otimes_{s\in
S}\left[\begin{array}[]{ccc}\frac{1}{2}\beta_{-}&-\beta_{-}&\frac{1}{2}\beta_{-}\\\
\frac{1}{2}\beta_{+}&-\frac{1}{2}(\beta_{-}+\beta_{+})&\frac{1}{2}\beta_{-}\\\
\frac{1}{2}\beta_{+}&-\beta_{+}&\frac{1}{2}\beta_{+}\end{array}\right].$

To compute the $\varepsilon^{2}$ term in (13), we can first subtract
$[T_{0}^{\prime},P_{0}]$ from $P_{0}^{\prime}$, leaving site-wise

$[[T_{0},P_{0}^{\prime}],P_{0}^{\prime}]=\frac{\beta_{+}-\beta_{-}}{2}[[T_{0}^{\prime},P_{0}],\left[\begin{array}[]{ccc}0&0&0\\\
1&0&-1\\\
0&0&0\end{array}\right]]=\frac{\beta_{+}-\beta_{-}}{2}\left[\begin{array}[]{ccc}-\beta_{-}&0&\beta_{-}\\\
-\beta_{-}&\beta_{-}-\beta_{+}&\beta_{+}\\\
-\beta_{+}&0&\beta_{+}\end{array}\right].$

Thus the final result to second order, using the $(+,-)$ part of the
$\psi_{\varepsilon}$ basis, is the effective PCA

$\hat{T}_{\varepsilon}=\otimes_{s\in
S}\left[\begin{array}[]{cc}1-\varepsilon\beta_{-}/2+\varepsilon^{2}\beta_{-}(\beta_{-}-\beta_{+})/4&\varepsilon\beta_{-}/2-\varepsilon^{2}\beta_{-}(\beta_{-}-\beta_{+})/4\\\
\varepsilon\beta_{+}/2-\varepsilon^{2}\beta_{+}(\beta_{+}-\beta_{-})/4&1-\varepsilon\beta_{+}/2+\varepsilon^{2}\beta_{+}(\beta_{+}-\beta_{-})/4\end{array}\right].$

The outcome is a two-state model in which to leading order there is a small
probability $\varepsilon\beta_{\mp}/2=\varepsilon(1+\alpha n_{\mp})/2$ per
timestep for transition from state $+$ to $-$, respectively $-$ to $+$.

Perhaps a reader can suggest (and apply the method to) a more significant
example.

##  6\. Further potential applications

One area to which the above results might be usefully applied is
metastability. Metastability is the phenomenon that an ergodic process may
spend long times exploring restricted subsets of the support of the
probability distribution, switching between them infrequently but sufficiently
to achieve ergodicity in the long run. One reference is [BdH].

As an application of the result, one can start from the paper [D] in which the
hypothesis is a reversible Markov process in continuous time with spectrum in
$[-\varepsilon,0]\cup[-\infty,-1]$ such that each function in the image of the
spectral projection $P$ for $[-\varepsilon,0]$ is bounded, and it is proved
that there is a partition of the state spaces into “metastable regions”. It is
not clear to me, however, whether there are relevant examples satisfying the
hypotheses. One might think that Glauber dynamics of the 2D Ising model below
the critical temperature would qualify, but I’m not aware that it is proved to
have a spectral gap. Nevertheless, if there are examples on product spaces
then the result of the present paper proves robustness of spectral gap and
hence of the phenomenon of metastability.

One could envisage the result also being useful to treat perturbation of
Markov dynamics with more than one stationary distribution, for example with
more than one communicating component. Addition of some interaction between
the communicating components typically makes the system have a single
communicating component but the result of this paper shows there is a
continuation of the spectral projection to a spectral projection with spectrum
contained near $1$ and the same rank (equal to the original number of
communicating components), and it gives strong control over the resulting
continuation. More substantially, one would like to deduce something about
metastability in ergodic finite versions of infinite PCA with non-unique
stationary distribution.

Another possible application is to perturbations of product systems in which
the units all have simple eigenvalue +1 and isolated spectrum near some
$\lambda$ in the open unit disk. The conclusion is that there persists an
invariant subspace with decay constant near $\lambda$.

##  7\. Discussion

It has been proved here that every spectral projection of a stochastic
operator on a product space persists $C^{r}$-smoothly with respect to the norm
(5) for $C^{r}$-smooth changes in the transition operator, again measured
using (5). This generalises the case of the rank-one projection onto a
stationary distribution, treated in [M1].

Second-order perturbation theory has been developed for families of such
operators and an example treated. Potential applications have also been
suggested to robustness of metastability and some other uses for multi-
component stochastic processes.

## Acknowledgements

I am grateful to Nick Higham for pointers to the literature on Sylvester
equations.

## References

  * [AM] Armstrong-Goodall J, MacKay RS, Dobrushin and Steif metrics are equal, arXiv REF.
  * [BR] Bhatia R, Rosenthal P, How and why to solve the operator equation $AX-XB=Y$? Bull Lond Math Soc 29 (1997) 1–21.
  * [BdH] Bovier A, den Hollander F, Metastability (Springer, 2015)
  * [D] Davies EB, Metastable states of symmetric Markov semigroups II, J London Math Soc 26 (1982) 541–50.
  * [K] Kato T, Perturbation theory for linear operators (Springer, 1980).
  * [M0] MacKay RS, Discrete breathers: classical and quantum, Physica A 288 (2000) 174–198.
  * [M1] MacKay RS, Robustness of Markov processes on large networks, J Difference Eqns & Applns 17 (2011) 1155–67.
  * [M2] MacKay RS, Management of complex systems, Nonlinearity 31 (2018) R52–65.
  * [N] Nevanlinna O, Sylvester equations and polynomial separation of spectra, arxiv:1904.07549 (2019)
  * [Sp] Spalek J, $t-J$ model Then and Now: a personal perspective from the pioneering times, Acta Phys Polonica A 111 (2007) 409–24.
  * [Sf] Steif JE, Convergence to equilibrium and space-time Bernoullicity for spin systems in the $M<\varepsilon$ case, Ergod Th Dyn Sys 11 (1991) 547–575.
  * [Sw] Stewart GW, Error and perturbation bounds for subspaces associated with certain eigenvalue problems, SIAM Rev 15 (1973) 727–64.
  * [Vh] Varah JM, On the separation of two matrices, SIAM J Num Anal 16 (1979) 216–22.
  * [Vn] Vasershtein LN, Markov processes over denumerable products of spaces describing large systems of automata, Problemy Peredachi Informatskii 5 (1969) 64–72.