# Exploring Strange Origin of Dirac Neutrino Masses at Hadron Colliders

Hooman Davoudiasl111email<EMAIL_ADDRESS>High Energy Theory Group, Physics
Department, Brookhaven National Laboratory, Upton, New York 11973, USA  Ian M.
Lewis222email<EMAIL_ADDRESS>Department of Physics and Astronomy,
University of Kansas, Lawrence, Kansas, 66045 USA  Matthew Sullivan333email:
<EMAIL_ADDRESS>High Energy Theory Group, Physics Department, Brookhaven
National Laboratory, Upton, New York 11973, USA

###### Abstract

We consider the possibility that Dirac neutrino masses may be a manifestation
of chiral symmetry breaking via non-perturbative QCD dynamics. The key role
played by light quarks in this mechanism can naturally lead to signals that
are accessible to hadron colliders. Bounds from charged meson decays imply a
dominant effect from the strange quark condensate. We propose a model for
Dirac neutrino mass generation with an extra Higgs doublet at the TeV scale
and significant coupling to strange quarks and leptons. Current data on
$D-\bar{D}$ mixing constrain the allowed parameter space of the model, and a
100 TeV $pp$ collider would either discover or largely exclude it. A distinct
feature of this scenario is that measurements of the of charged Higgs leptonic
branching ratios can distinguish between “normal” and “inverted” neutrino mass
hierarchies, complementing future determinations at neutrino oscillation
experiments.

##  I Introduction

The origin of non-zero but tiny neutrino masses $m_{\nu}\lesssim 1$ eV remains
an open fundamental question in particle physics. The smallness of $m_{\nu}$
implies very suppressed interactions with the Higgs and its vacuum expectation
value $\langle{H}\rangle=v/\sqrt{2}\approx 174$ GeV. Hence, typically,
neutrino mass models either involve very massive fermions which are
inaccessible to direct observations, or else negligible Yukawa couplings.
These considerations generally leave only indirect rare processes, such as
neutrinoless double $\beta$ decay, as the possible signals of the underlying
neutrino mass generation mechanism.

Given that $m_{\nu}\neq 0$ requires a connection with a source of electroweak
symmetry breaking (EWSB), the involvement of $\langle{H}\rangle$ seems like a
necessary ingredient in any plausible model. However, there is another source
of EWSB in the Standard Model (SM), the light quark condensate
$\langle{\bar{q}\,q}\rangle\sim-(300~{}\text{MeV})^{3}$, with
$q\in\\{u,d,s\\}$. This effect – generated by non-perturbative QCD
interactions – has been largely ignored as a possible origin of neutrino
masses. The much smaller contribution from $\langle{\bar{q}\,q}\rangle$ to
EWSB suggests that it could have far less suppressed interactions with
neutrinos in an underlying mass generation mechanism. Therefore, one may
generically expect that key states in such a mechanism would be accessible to
hadron colliders, where light quarks are important initial states for
producing new particles.

In this work, we will consider generation of Dirac neutrino masses through QCD
quark condensation. First, we will examine the general aspects of this setup
using effective field theory (EFT). We will find that constraints from meson
decays lead us to consider the strange quark condensate
$\langle{\bar{s}\,s}\rangle$ as the dominant source of neutrino masses in this
scenario. In order to realize the mechanism, we will introduce a model that
includes a TeV scale Higgs doublet $H_{2}$ with $\mathcal{O}{(1)}$ couplings
to strange quarks as well as right-handed partners assumed for each SM
neutrino flavor. In our convention, the SM-like Higgs doublet $H$ is denoted
by $H_{1}$. We will show that, depending on the specific values of couplings,
heavy scalars can be resonantly produced and detected in di-jet or “lepton +
missing energy” final states, under minimal assumptions regarding the model.
Prospects for flavor experiments and potential consequences for supernova
observations will also be briefly discussed. Prior works that have considered
neutrino mass generation from QCD dynamics include Refs. Thomas and Xu (1992);
McDonald (1996); Ibanez et al. (2001); Davoudiasl and Everett (2006); Babič et
al. (2021).

##  II EFT Framework

Let us begin with an effective theory description and consider the following
dimension-6 interaction

$O_{D}=\zeta\frac{[\bar{Q}\,s]\epsilon[\bar{L}\,\nu_{R}]}{M_{D}^{2}}+{\small\rm
H.C.}\,,$ (1)

where $Q$ and $s$ are the second generation quark doublet and the strange
quark singlet, respectively, $L$ is a lepton doublet and $\epsilon$ is the
2-dimensional Levi-Civita symbol contracting the $SU(2)_{L}$ indices. The
right-handed field $\nu_{R}$ is the Dirac partner of the left-handed $\nu_{L}$
component of $L$. In Eq. (1), $\zeta$ is a $3\times 3$ matrix corresponding to
three generations of $L$ and $\nu_{R}$, assuming three massive neutrinos. Note
that $[\bar{Q}\,s]$ has the electroweak quantum numbers of an anti-Higgs
doublet and condensation of its iso-spin $-1/2$ component can generate a
neutrino mass term. We have suppressed all generational indices in Eq. (1).
The effective scale $M_{D}$ is related to mass scales in an ultraviolet (UV)
completion and generally includes the effect of couplings in such a theory. In
principle, one could also have a similar operator with $s$ replaced by the
right-handed up or down quark. However, as will be discussed below, those
choices will lead to significant constraints. To obtain a neutrino mass
$m_{\nu}\sim 0.1$ eV upon quark condensation, we then need $M_{D}\sim 16$ TeV,
where we have used $\langle{\bar{s}\,s}\rangle\approx-(300~{}\text{MeV})^{3}$
Davies et al. (2019).

Constraints from charged meson $P^{+}$ decays of the type $P^{+}\to
e+\text{``missing energy"}$ were considered in Ref. Davoudiasl and Everett
(2006). We adapt the results of Ref. Davoudiasl and Everett (2006) for the
partial width for meson decay $P^{+}\to e^{+}\,\nu_{R}$

$\Gamma(P^{+}\to
e^{+}\,\bar{\nu}_{R})=\frac{\sum_{i}|\zeta_{ei}|^{2}}{64\pi\,M_{D}^{4}}f_{P}^{2}\mu_{P}^{2}m_{P}\,,$
(2)

where the sum over $i$ is over right-handed neutrinos, $f_{P}$ and $m_{P}$ are
the decay constant and mass of $P$. The parameter $\mu_{P}$ is given by

$\frac{m_{\pi}^{2}}{2\bar{m}}\quad;\quad\frac{m_{K}^{2}}{m_{s}+\bar{m}}\quad;\quad\frac{m_{D_{s}}^{2}}{m_{c}+m_{s}}\,,$
(3)

for $P=\pi,K$, or $D_{s}$. The mass of the $D_{s}^{+}$ is $m_{D_{s}}=2.0$ GeV.
In the above, $\bar{m}=(m_{u}+m_{d})/2$; $m_{u}$, $m_{d}$, $m_{s}$, and
$m_{c}$ are the up, down, strange, and charm quark masses, respectively. The
current relevant branching ratio limits are Zyla et al. (2020)

$\displaystyle{\rm Br}(\pi^{+}\rightarrow e^{+}\,\nu)$
$\displaystyle=\,(1.23\pm 0.004)\times 10^{-4}$ (4) $\displaystyle{\rm
Br}(K^{+}\rightarrow e^{+}\,\nu)$ $\displaystyle=\,(1.582\pm 0.007)\times
10^{-5}$ (5) $\displaystyle{\rm Br}(D_{s}^{+}\rightarrow e^{+}\,\nu)$
$\displaystyle<\,8.3\times 10^{-5}\quad\text{(90\% C.L.).}$ (6)

The decay constants are $f_{\pi}=130$ MeV, $f_{K}=156$ MeV, and
$f_{D_{s}}=250$ MeV; and the lifetimes are $\tau_{\pi^{+}}=2.6\times 10^{-8}$
s, $\tau_{K^{+}}=1.2\times 10^{-8}$ s, and $\tau_{D_{s}^{+}}=5.0\times
10^{-13}$ s Zyla et al. (2020). The total widths of the mesons can be
calculated by $\Gamma_{P}=1/\tau_{P}$.

To find bounds, we calculate when the partial width in Eq. (2) saturate the
uncertainty bands of the branching ratios444black While the quoted
uncertainties on the branching ratios are purely experimental, we have also
checked our calculations against the $R$-ratios Zyla et al. (2020) where the
theory uncertainties Cirigliano and Rosell (2007); Bryman et al. (2011, 2021)
are subdominant. The conclusions are the same. and assume
$\sum_{i}|\zeta_{ei}|^{2}\sim 1$. The results then suggest when $P=\pi^{+}$
the current limits on these decay require $M_{D}\gtrsim 60$ TeV Davoudiasl and
Everett (2006) and for $P=K^{+}$ we have $M_{D}\gtrsim 100$ TeV, which would
not allow large enough neutrino masses, based on the preceding discussion.
Note, however, that if $Q$ is the second generation quark doublet $(c,s)$,
then $P^{+}=D_{s}^{+}$ and as we will show later, the current limits on
$D_{s}^{+}\to e^{+}+\text{``missing energy"}$ would allow $M_{D}\gtrsim 3$
TeV, which can lead to the observed neutrino masses. Hence, we will assume the
dominance of the second generation interactions in Eq. (1), which will then
require a singlet $s$ to get a quark condensate, since the charm quark $c$ is
too heavy to condense.

Another aspect of phenomenology in our framework is to make sure that right-
handed neutrinos do not get populated in the early Universe McDonald (1996)
and become extra degrees of freedom during Big Bang Nucleosynthesis (BBN) and
cosmic microwave background eras. This can be arranged for $M_{D}\gtrsim 10$
TeV Davoudiasl et al. (2020, 2021), by assuming a low reheat temperature
$T_{\rm rh}\lesssim 0.1$ GeV which nonetheless is high enough to allow a
consistent cosmology, given that BBN only requires $T_{\rm rh}\gtrsim 4$ MeV
Hannestad (2004) and there is currently no known cosmological input requiring
a higher reheat temperature.

We also point out that supernova constraints from SN1987A can be satisfied if
$M_{D}\gtrsim 13$ TeV Thomas and Xu (1992), which does not exclude the range
of parameters relevant to our model and can be accommodated. Note that the
dominant contribution in our model is from the strange quark, which is
presumably not as abundant in a supernova as the up and down quarks making up
the nucleons. Hence, the above supernova constraint – which was derived
assuming neutrino masses from $\langle{\bar{u}u}\rangle$ and the coupling of
first generation quarks to right-handed neutrinos – could be somewhat relaxed.

Before going further, we would like to make a comment about the possibility of
ruling out our model, solely by improving the bound (6) on ${\rm
Br}(D_{s}^{+}\rightarrow e^{+}\,\nu)$. This bound, which from the above
discussion yields $M_{D}\gtrsim 3$ TeV, is based on the Belle experiment data
set of $913$ fb-1 Zupanc et al. (2013). Given that our model requires
$M_{D}\sim 16$ TeV, and optimistically assuming that future measurements are
only statistically limited, to probe our model one needs a data set that is
larger by $(16/3)^{4}$, i.e. $\sim 700$ ab-1. The Belle II experiment is
expected to accumulate $50$ ab-1 over its entire run Altmannshofer et al.
(2019), and is currently the only foreseeable machine that could improve Belle
results. Hence, we do not expect that the model proposed here be ruled out in
the foreseeable future, as the constraints on ${\rm Br}(D_{s}^{+}\rightarrow
e^{+}\,\nu)$ get more stringent. This suggests that a 100 TeV collider could
be the only envisioned facility that can potentially discover or largely rule
out our scenario555We would like to note that since our model consists of
Dirac neutrinos, discovery of neutrinoless double beta decay would rule out
the scenario proposed here., in the context of a UV model, as will be
described later in this work.

###  II.1 Flavor Model

Introduction of a new Higgs doublet which couples to light fermions could
induce unwanted flavor violating effects. To mitigate such effects, we will
adopt the Spontaneous Flavor Violation (SFV) Egana-Ugrinovic et al. (2019a, b)
framework. This would allow us to have interactions of light fermions with
$H_{2}$ with $\mathcal{O}{(1)}$ strength, necessary to obtain sufficiently
large neutrino masses, without causing large deviations from flavor
constraints. These couplings, in particular interaction with light quarks,
would then allow significant resonant production of $H_{2}$ at high energy
hadron colliders, which can lead to a testable neutrino mass mechanism. We use
SFV of the up-type, since the down-type would require large $b$ quark Yukawa
coupling $\gtrsim 1$.

The general 2HDM Yukawa sector with neutrinos would look like

$\displaystyle\displaystyle\sum_{a=1}^{2}$
$\displaystyle-\lambda^{a}_{u}\bar{Q}\,\epsilon\,H_{a}^{*}\,u-\lambda^{a}_{d}\bar{Q}\,H_{a}\,d-\lambda^{a}_{\nu}\bar{L}\,\epsilon\,H_{a}^{*}\,\nu_{R}-\lambda^{a}_{\ell}\bar{L}\,H_{a}\,\ell\,$
(7) $\displaystyle+{\rm H.C.}$

In the up-type SFV framework, these Yukawa coupling matrices will take the
form

$\begin{array}[]{ccc}\lambda^{1}_{u}=V_{CKM}^{\dagger}Y_{u},&\quad\lambda^{1}_{d}=Y_{d},&\quad\lambda^{1}_{\ell}=Y_{\ell}\\\
\lambda^{2}_{u}=\xi
V_{CKM}^{\dagger}Y_{u},&\quad\lambda^{2}_{d}=K_{d},&\quad\lambda^{2}_{\ell}=\xi_{\ell}Y_{\ell},\end{array}$
(8)

where $Y_{u}$, $Y_{d}$, and $Y_{\ell}$ are the diagonal SM Yukawa couplings in
the mass basis for up-type quarks, down-type quarks, and charged leptons,
respectively. Here, $V_{CKM}$ is the CKM matrix,
$K_{d}=diag(\kappa_{d},\kappa_{s},\kappa_{b})$ is flavor diagonal with real
entries $\kappa_{d,s,b}$; $\xi$ and $\xi_{\ell}$ are real constants. For the
neutrino Yukawa couplings, we will use the scheme

$\displaystyle\lambda^{1}_{\nu}$ $\displaystyle=$ $\displaystyle 0$
$\displaystyle\lambda^{2}_{\nu}$ $\displaystyle=$ $\displaystyle
V_{PMNS}K_{\nu},$ (9)

where $K_{\nu}=diag(\kappa_{\nu,1},\kappa_{\nu,2},\kappa_{\nu,3})$ with
$\kappa_{\nu,i}$ real and, in our model, proportional to the neutrino masses
$m_{i}$. Assuming the above structure of Yukawa couplings, we then find that
the effective operator in Eq. (1) can be obtained by integrating out $H_{2}$,
with the identifications $\zeta\to\kappa_{s}\lambda^{2}_{\nu}$ and $M_{D}\to
M_{H_{2}}$, where $M_{H_{2}}$ is the heavy Higgs doublet mass. Inserting
$\langle s\bar{s}\rangle$ and rotating to the mass basis, Eq. (1) becomes a
neutrino mass term and we can make the identification

$\displaystyle\kappa_{\nu,i}=\frac{M_{H_{2}}^{2}}{\kappa_{s}\,\langle\bar{s}s\rangle}m_{i}.$
(10)

Interestingly, perturbative unitarity gives an upper and lower bound on
$\kappa_{s},\kappa_{\nu,i}$ from Eq. (10). Through the usual perturbative
unitarity arguments Chanowitz et al. (1979), the scattering
$f\bar{f}\rightarrow H_{2}\rightarrow f\bar{f}$ gives an upper bound
$|\kappa_{s}|,|\kappa_{\nu,i}|\leq\sqrt{8\pi}$. Using Eq. (10), this
translates into a lower bound on $\kappa_{s}$:

$\displaystyle|\kappa_{s}|=\frac{M_{H_{2}}^{2}}{|\kappa_{\nu,i}|\,|\langle\bar{s}s\rangle|}m_{i}\geq\frac{M_{H_{2}}^{2}}{\sqrt{8\pi}\,|\langle\bar{s}s\rangle|}m_{i},$
(11)

and similarly for $\kappa_{\nu,i}$. The most constraining limit comes from the
heaviest neutrino, which depends on if the neutrino masses have a normal or
inverted hierarchy. As discussed below, the cosmological bound of $\sum
m_{i}<0.12$ eV Ivanov et al. (2020) can place an upper bound of $0.03$ eV and
$0.016$ eV on the heaviest neutrino in the normal and inverted hierarchies,
respectively. Hence, normalizing the neutrino mass to $0.03$ eV, we find

$\displaystyle\sqrt{8\pi}>|\kappa_{s}|,\,|\kappa_{\nu,i}|>5.5\times
10^{-3}\left(\frac{M_{H_{2}}}{5~{}{\rm
TeV}}\right)^{2}\frac{m_{i}}{0.03~{}{\rm eV}}.$ (12)

Although the lower bound is small, in the scenario presented here, it is not
possible to decouple $H_{2}$ by setting its couplings with fermions to zero.

For the up-type parameters, we focus on $\xi\ll 1$ and
$\kappa_{d},\kappa_{b}\ll\kappa_{s}$, with $\kappa_{s}\lesssim 1$. We will
justify these parameter values, as well as $\lambda^{1}_{\nu}=0$, within a UV
framework later. Although flavor-changing neutral currents are suppressed in
the SFV framework, many processes still constrain the allowed parameters. For
our strange-dominant scheme, though, the only additional relevant constraint
comes from $D-\bar{D}$ mixing. Reference Egana-Ugrinovic et al. (2019b)
calculates the contribution from charged Higgs diagrams in this SFV model to
$D-\bar{D}$ mixing. Since the SM contribution is not currently known, we bound
the new physics contribution to the mass difference by the 95.5% upper limit
on the observed value from Ref. Aaij et al. (2021). The bounds from this
procedure are shown in Fig. 1.

Figure 1: Limits on the strange quark coupling $\kappa_{s}$ as a function of
the mass $M_{H_{2}}$ of the heavy Higgs doublet from $D-\bar{D}$ mixing
constraints.

##  III UV Framework

Next, we will introduce a UV framework for how the assumed Yukawa structures
in Eqs. (7), (8), and (9) may emerge at the TeV scale. We have implicitly
assumed that the two Higgs doublets $H_{1}$ and $H_{2}$ have negligible mixing
and hence are, to a good approximation, mass eigenstates. Here,
$\langle{H_{1}}\rangle\neq 0$ is responsible for masses of all fermions except
for neutrinos.

We note that, so far, it is not forbidden to have the usual Dirac neutrino
Yukawa coupling to the SM Higgs doublet $H_{1}$, containing the real scalar
observed at $m_{H}\approx 125$ GeV. In order to eliminate this interaction –
which requires an exponentially small Yukawa coupling $\sim 10^{-12}$ – in
favor of the QCD-mediated mechanism, we assume that both $H_{2}$ and $\nu_{R}$
are charged under a global symmetry, which we will take to be a $Z_{2}$
parity:

$Z_{2}(H_{2})=Z_{2}(\nu_{R})=-1.$ (13)

Other SM fields are $Z_{2}$-even.

The above charge assignment, however, does not allow the coupling to the $s$
quark assumed in Eq. (7), unless it has a chiral $Z_{2}$-odd parity. Since
$H_{1}$ is $Z_{2}$-even, a $Z_{2}$-odd strange quark would be massless, which
is not acceptable. We note that Refs. Thomas and Xu (1992); Davoudiasl and
Everett (2006) had entertained the possibility that the condensate was formed
by an up quark, which would remain massless. This possibility is disfavored by
hadron phenomenology and lattice investigations Bazavov et al. (2018);
Alexandrou et al. (2020). In this work, we will describe possible model
building solutions that will not lead to a massless light quark. Hence, the
above $s$ quark coupling to $H_{2}$ must be induced through spontaneous
$Z_{2}$ breaking. To accomplish this, we postulate a new scalar $\phi$, with
$Z_{2}(\phi)=-1$, interacting through the dim-5 operator,

$c_{s}\frac{\phi\,H_{2}\bar{Q}s}{\Lambda}\,.$ (14)

Once $\langle{\phi}\rangle\neq 0$, it leads to the assumed coupling in Eq.
(7), given by $\lambda^{2}_{s}=c_{s}\langle{\phi}\rangle/\Lambda$, with
$c_{s}$ being a complex constant.

Our UV framework also needs to explain why a similar dim-5 operator
$\sim\phi\bar{L}\epsilon H_{1}^{*}\nu_{R}$, leading us back to a Yukawa
coupling to the SM Higgs, is not generated. This could be simply a result of
mass hierarchies in the underlying model. That is, one could generate the
interaction in Eq. (14), for example, through the exchange of a $Z_{2}$-odd
vector-like quark $\chi$, of mass $m_{\chi}\gtrsim 1$ TeV, that has the SM
quantum numbers of a right-handed strange quark $s$. In the far UV, this
requires the following interactions

$y_{s}\phi\bar{\chi}\frac{(1+\gamma_{5})}{2}s+y_{\chi}H_{2}\bar{Q}\frac{(1+\gamma_{5})}{2}\chi+{\small\rm
H.C.}\,,$ (15)

which would yield $c_{s}/\Lambda=y_{s}y_{\chi}/m_{\chi}$. If there is no
analogue $Z_{2}$-odd vector-like lepton, or if it is very heavy, then we would
not obtain a dim-5 operator $\sim\phi\bar{L}\epsilon H_{1}^{*}\nu_{R}$, or
else it would have negligible effect.

One may also worry that non-perturbative $Z_{2}$-violating gravitational
effects of the type

$\frac{\phi\bar{L}\,\epsilon\,H_{1}^{*}\,\nu_{R}}{M_{\rm P}}\,,$ (16)

suppressed by the Planck mass $M_{\rm P}\sim 10^{19}$ GeV, could induce large
competing effects. However, for $\langle{\phi}\rangle/M_{\rm P}\ll 10^{-12}$
these effects are expected to be negligible. In what follows, we will assume
$\langle{\phi}\rangle\sim$ (1-1000) TeV, which satisfies the above condition.
In Eq. (14), it is implicitly assumed that $\chi$ can be integrated out to
yield an effective theory description valid for $\langle{\phi}\rangle\lesssim
m_{\chi}$. Therefore, if $\langle{\phi}\rangle$ is near the lower end of this
interval, we could also expect to have TeV-scale $\chi$ particles, which are
$SU(3)_{c}$ triplets and can provide another potential signal at hadron
colliders. Given the above model building considerations, we can account for
the features of the assumed Yukawa sector in Eq. (7).

##  IV Hadron collider phenomenology

Now we discuss the hadron collider phenomenology of our model. This model
consists of four heavy scalars from $H_{2}$: charged Higgs bosons $H^{\pm}$, a
scalar $H$, and a pseudoscalar $A$. The major production modes of these scalar
are through the $\kappa_{s}$ coupling, introduced in the flavor model:

$\displaystyle s\bar{s}\rightarrow H/A,c\bar{s}\rightarrow H^{+},\,{\rm
and}\,s\bar{c}\rightarrow H^{-}.$ (17)

The distinctive feature of this model, is that in a wide range of the
parameter space the charged Higgs decay into leptons
$H^{-}\rightarrow\ell_{i}\bar{\nu}_{j}$ is substantial. At hadron colliders,
the neutrino flavors of the final state cannot be determined and must be
summed over. The partial widths of the charged Higgs bosons into leptons are
then

$\displaystyle\Gamma(H^{\pm}\rightarrow\ell_{i}\nu)=\frac{1}{16\,\pi}\left(\lambda^{2}_{\nu}{\lambda^{2}_{\nu}}^{\dagger}\right)_{ii}\,M_{H_{2}},$
(18)

and the partial width into jets is

$\displaystyle\Gamma(H^{\pm}\rightarrow c\,s)$ $\displaystyle=$
$\displaystyle\frac{3}{16\,\pi}\kappa_{s}^{2}\,M_{H_{2}}$ (19)
$\displaystyle\approx$ $\displaystyle 2.4\times
10^{-3}\,\left(\frac{\kappa_{s}}{0.2}\right)^{2}\,M_{H_{2}}.$

Using the values for the neutrino mixing matrices and mass differences found
in Ref. Esteban et al. (2020), as well as the $\lambda_{\nu}^{2}$ parameter
values above, we find for the normal hierarchy the partial widths into
different lepton final states are:

$\displaystyle\Gamma(H^{\pm}\rightarrow e\,\nu)$ $\displaystyle\approx$
$\displaystyle 3.3\times 10^{-5}\left[1+12\left(\frac{m_{1}}{0.03~{}{\rm
eV}}\right)^{2}\right]\left(\frac{M_{H_{2}}}{5~{}{\rm
TeV}}\right)^{4}\left(\frac{\kappa_{s}}{0.2}\right)^{-2}\,M_{H_{2}}$ (20)
$\displaystyle\Gamma(H^{\pm}\rightarrow\mu\,\nu)$ $\displaystyle\approx$
$\displaystyle 6.1\times 10^{-4}\left[1+0.63\,\left(\frac{m_{1}}{0.03~{}{\rm
eV}}\right)^{2}-3.5\times
10^{-3}\cos\delta_{CP}\right]\left(\frac{M_{H_{2}}}{5~{}{\rm
TeV}}\right)^{4}\left(\frac{\kappa_{s}}{0.2}\right)^{-2}\,M_{H_{2}}$ (21)
$\displaystyle\Gamma(H^{\pm}\rightarrow\tau\,\nu)$ $\displaystyle\approx$
$\displaystyle 4.6\times 10^{-4}\left[1+0.83\,\left(\frac{m_{1}}{0.03~{}{\rm
eV}}\right)^{2}+4.7\times
10^{-3}\cos\delta_{CP}\right]\left(\frac{M_{H_{2}}}{5~{}{\rm
TeV}}\right)^{4}\left(\frac{\kappa_{s}}{0.2}\right)^{-2}\,M_{H_{2}},$ (22)

where $m_{1}$ is the lightest neutrino mass. The cosmological bound on the sum
of neutrino masses is $\sum m_{i}\lesssim 0.12$ eV Ivanov et al. (2020). For
the normal hierarchy, this translates to a bound on the lightest neutrino mass
of $m_{1}\lesssim 0.030$ eV, to which we have normalized the lightest neutrino
mass. As can be seen immediately, the dependence of $\delta_{CP}$ is minimal.
Hence, we set $\delta_{CP}=0$ for simplicity.

Figure 2:  (a,b) Branching ratios of $H^{\pm}$ as a function of the lightest
neutrino mass into (solid) electrons, (dashed) muons, and (dot-dashed) taus.
(c,d) Ratios of branching ratios between different lepton final states as a
function of the lightest neutrino mass: (solid) electron to muon, (dashed)
muon to tau, and (dot-dashed) electron to tau. The normal hierarchy is shown
in (a,c) and the inverted hierarchy in (b,d). The $H_{2}$ mass is set to 5
TeV, and gray regions correspond to when $\sum m_{i}>0.12$ eV.

Figures 2(a,c) shows the (a) branching ratios of $H^{\pm}$ into leptons and
(c) ratios of branching ratios into leptons as a function of hte lightest
neutrino mass in the normal hierarchy. As is clear in these figures, as well
as Eqs. (20-22), the branching ratio into electrons is by far the smallest
and, while we expect more $\mu$s than $\tau$s, the branching ratios into muons
and taus are comparable. This can be understood by noting that in the normal
hierarchy $m_{3}$ is the heaviest neutrino and couples most strongly to heavy
Higgses. Since $\theta_{13}$ is the smallest mixing angle, the electron
couples least strongly. Since $\theta_{12}$ and $\theta_{23}$ are comparable,
the coupling to muons and taus are comparable. However, as is also clear in
the figures and Eqs. (20-22), the branching ratio into electrons is most
sensitive to the value of the lightest neutrino mass.

In the inverted hierarchy, the situation is very different. Now $m_{2}$ is the
heaviest neutrino with $m_{1}$ closest in mass. In this case, the partial
width into electrons is larger than muons and taus:

$\displaystyle\Gamma(H^{\pm}\rightarrow e\,\nu)$ $\displaystyle\approx$
$\displaystyle 1.0\times 10^{-3}\left[1+0.11\left(\frac{m_{3}}{0.016~{}{\rm
eV}}\right)^{2}\right]\left(\frac{M_{H_{2}}}{5~{}{\rm
TeV}}\right)^{4}\left(\frac{\kappa_{s}}{0.2}\right)^{-2}\,M_{H_{2}}$ (23)
$\displaystyle\Gamma(H^{\pm}\rightarrow\mu\,\nu)$ $\displaystyle\approx$
$\displaystyle 4.6\times 10^{-4}\left[1+0.24\,\left(\frac{m_{3}}{0.016~{}{\rm
eV}}\right)^{2}-4.7\times
10^{-3}\cos\delta_{CP}\right]\left(\frac{M_{H_{2}}}{5~{}{\rm
TeV}}\right)^{4}\left(\frac{\kappa_{s}}{0.2}\right)^{-2}\,M_{H_{2}}$ (24)
$\displaystyle\Gamma(H^{\pm}\rightarrow\tau\,\nu)$ $\displaystyle\approx$
$\displaystyle 6.2\times 10^{-4}\left[1+0.18\,\left(\frac{m_{3}}{0.016~{}{\rm
eV}}\right)^{2}+3.5\times
10^{-3}\cos\delta_{CP}\right]\left(\frac{M_{H_{2}}}{5~{}{\rm
TeV}}\right)^{4}\left(\frac{\kappa_{s}}{0.2}\right)^{-2}\,M_{H_{2}},$ (25)

now $m_{3}$ is the lightest neutrino mass and the cosmological bound $\sum
m_{i}\lesssim 0.12$ eV requires $m_{3}\lesssim 0.016$ eV. Figures 2(b,d) show
the branching ratio results for the inverted hierarchy. As is clearly seen, in
the inverted hierarchy charged Higgs decays into electrons always dominates
the lepton modes. Additionally, the dependence on $\delta_{CP}$ is still very
weak, and all modes depend on the exact value of the lightest neutrino mass at
a comparable level.

Our model of neutrino masses has distinct predictions for charged Higgs decays
that separate the normal and inverted hierarchies. In the normal hierarchy the
charged Higgs decays into muons dominate, followed by taus, and electrons. In
the inverted hierarchy, things are, well, inverted: decays into electrons
dominate, then taus, then muons. Hence, our model has the remarkable feature,
that if true, the neutrino hierarchy could be determined by counting the
number of electrons and muons coming from charged Higgs decays.

The neutral scalars, $H$ and $A$, can only decay into neutrinos and
$s\bar{s}$. Again, since neutrinos are not detectable, we must sum over them
in the final state. The partial widths are then

$\displaystyle\Gamma(A/H\rightarrow s\bar{s})$ $\displaystyle=$
$\displaystyle\frac{3\,\kappa_{s}^{2}}{32\,\pi}M_{H_{2}}\approx 1.1\times
10^{-3}\left(\frac{\kappa_{s}}{0.2}\right)^{2}M_{H_{2}}$ (26)
$\displaystyle\Gamma(A/H\rightarrow\nu\bar{\nu})$ $\displaystyle=$
$\displaystyle\frac{\kappa_{\nu,1}^{2}+\kappa_{\nu,2}^{2}+\kappa_{\nu,3}^{2}}{32\pi}M_{H_{2}}$
(27) $\displaystyle\approx$ $\displaystyle\begin{cases}\displaystyle 5.5\times
10^{-4}\left[1+1.0\left(\frac{m_{1}}{0.03~{}{\rm
eV}}\right)\right]\left(\frac{M_{H_{2}}}{5~{}{\rm
TeV}}\right)^{4}\left(\frac{\kappa_{s}}{0.2}\right)^{-2}\,M_{H_{2}}\quad&{\rm
Normal~{}hierarchy}\\\ \displaystyle 1.0\times
10^{-3}\left[1+0.16\left(\frac{m_{3}}{0.016~{}{\rm
eV}}\right)\right]\left(\frac{M_{H_{2}}}{5~{}{\rm
TeV}}\right)^{4}\left(\frac{\kappa_{s}}{0.2}\right)^{-2}\,M_{H_{2}}&{\rm
Inverted~{}hierarchy}\end{cases}$

Figure 3:  (a) black$5\sigma$ discovery reaches of $\kappa_{s}$ at a 100 TeV
pp collider as a function of the $H_{2}$ mass with (black) 3 and (blue) 10
ab-1 of data. The normal hierarchy is used, $\delta_{CP}=0$, and the lightest
neutrino mass is set to its maximum value. Dijet discovery is possible for
parameter values above the solid lines, and within the dashed (dot-dashed)
curve 10 (25) events can be observed in $pp\rightarrow H^{\pm}\rightarrow
e\nu+\mu\nu$. The shaded region is ruled out by current flavor constraints.
(b) Regions of parameters space needed to observe 10 events with 10 ab-1 in
the $e\nu+\mu\nu$ channel for (black) the normal hierarchy and (blue) the
inverted hierarchy. The region inside the solid lines are for when neutrino
masses satisfy $\sum m_{i}=0.12$ eV and dashed lines for when the lightest
neutrino is massless. We have set $\delta_{CP}=0$. Dijet projections are the
same as in panel (a).

As discussed previously, the major production mode of the heavy Higgs bosons
is via $s$-channel production. Pair production through EW gauge bosons only
depends on gauge couplings and is therefore more model independent than
$s$-channel production. However, the pair production rates are exceedingly
small. At a 14 TeV LHC, pair production for two 2 TeV scalars is $\sim 0.05$
ab while at 100 TeV the pair production rate is $\sim 0.1$ fb Davoudiasl et
al. (2020). For $\kappa_{s}=0.2$, the single production of a 2 TeV heavy
scalar boson is $\sim 4-8$ fb at the 14 TeV LHC and black$\sim 0.9$ pb at 100
TeV666blackThe range of cross sections denote the difference between charged
and neutral scalar production. Due to the high $x$ pdfs, at the LHC this
difference is more pronounced. At 100 TeV, these rates are within $10\%$ of
each other.. Hence, we focus on single production. Additionally, a $4$ fb
cross section can create $\sim 10,000$ heavy Higgs states with 3 ab-1 of data
at the LHC. blackThere are two major signal modes: di-jet with
$H/A/H^{\pm}\rightarrow jj$ and lepton plus missing energy with
$H^{\pm}\rightarrow e\nu/\mu\nu$. While this is a relatively large rate, the
QCD di-jet backgrounds are large and branching ratios into lepton plus
neutrino are percent level. In the di-jet channel, extrapolating to 3 ab-1 of
data from current ATLAS Aad et al. (2020) and CMS Sirunyan et al. (2020)
results, this translates to black ruling out $\kappa_{s}\gtrsim 0.3-0.6$ at
the $2\sigma$-level.777blackThe range of $\kappa_{s}$ is due to ATLAS and CMS
having different exclusion regions. For details on the extrapolation method,
see Ref. Davoudiasl et al. (2021). However, as seen in Fig. 1, these values
are already ruled out by flavor constraints. As the mass increases, the HL-LHC
di-jet exclusion become significantly worse than the flavor constraints. black
In the lepton plus neutrino channels, for a 2 TeV charged scalar with
$\kappa_{s}=0.2$, the $pp\rightarrow H^{\pm}\rightarrow e\nu/\mu\nu$ cross
sections are $\sim(3-7)\times 10^{-5}$ pb, where this range encompasses both
hierarchies and electrons/muons. Using a simple root luminosity extrapolation
of current ATLAS Aad et al. (2019) and CMS CMS (2021) searches for
$W^{\prime}\rightarrow e\nu/\mu\nu$, we find with 3 ab-1 of data the
experiments could be sensitive to cross sections of $(8-17)\times 10^{-5}$ pb.
This is tantalizingly close to our signal cross section, with the inverted
hierarchy in the $H^{\pm}\rightarrow e\nu$ mode being the most promising
signal to observe. Although this may bear more detailed studies, we have
assumed that both statistical and systematic uncertainties scale with
luminosity. The true sensitivity will depend on the details of improvements in
systematic uncertainties in these channels. Hence, we focus on signals at a
100 TeV pp collider.

In Fig. 3 we show the discovery reach as a function of the heavy Higgs mass at
a 100 TeV pp collider with 3 and 10 ab-1 of data. Both (a) and (b) consider
single production for $H^{\pm},H,A$. Single production cross sections are
generated via MadGraph5_aMC@NLO Alwall et al. (2014) with the model
implemented using FeynRules Christensen and Duhr (2009); Alloul et al. (2014).
We consider all four heavy states, $H^{\pm},H,A$, to have degenerate masses
and add their cross sections incoherently. Their collective mass is labelled
as $M_{H_{2}}$. Indeed, this assumption greatly alleviates the bounds from the
oblique parameters Peskin and Takeuchi (1990, 1992); Barbieri et al. (2006);
Haber and O’Neil (2011); Ahriche et al. (2015).

The discovery reach in the normal hierarchy is shown in Fig. 3. We set
$\delta_{CP}=0$ and $m_{1}$ to the maximum that saturates the cosmological
bound $\sum m_{i}\lesssim 0.12$ eV. Areas above the solid curves could be
discovered in the dijet channel, where we have adapted the projected dijet
search reach Golling et al. (2017)888See Ref. Davoudiasl et al. (2021) for
details on how we project dijet searches.. Areas within the dotted and dash-
dot curves can be discovered with 10 or 25 total events, respectively, in the
$pp\rightarrow H^{\pm}\rightarrow e\nu$ and $pp\rightarrow
H^{\pm}\rightarrow\mu\nu$ channels. blackThe number of events for discovery is
chosen considering zero background, we need 25 events to discover with
Gaussian statistics and $\sim 10$ with Poisson (assuming 1 background event).
A more complete study with a full background simulation would be needed to
determine the precise discovery region. However, we expect our estimates to
encompass the true range since background will be small.

blackWe can understand the curves in 3 by noting that $\kappa_{s}$ governs the
production cross section. Hence, at low $\kappa_{s}$ not enough $H^{\pm}$ are
produced to create a signal. At high $\kappa_{s}$ there is enough cross
section, but decays into jets dominate, and the branching ratios ${\rm
BR}(H^{\pm}\rightarrow e\nu)$ and ${\rm BR}(H^{\pm}\rightarrow\mu\nu)$ are too
small to generate enough events. Hence, the leptonic searches enclose an area
in the $\kappa_{s}-M_{H_{2}}$ instead of providing the simple lower bound of
the dijet searches. The regions ruled out by flavor searches discussed
earlier.

In Fig. 3, the regions enclosed by the curves are the parameter regions that
will provide 10 events in the $e\nu$ and $\mu\nu$ channels with 10 ab-1 of
data. We show the regions for the normal and inverted hierarchies, with the
respective lightest neutrino masses set to zero or saturating the cosmological
bound $\sum m_{i}\lesssim 0.12$ eV. As can be seen, the 100 TeV collider is
always more sensitive to the inverted hierarchy. This can be understood by
noting that Fig. 2 and Eqs. (20-22,23-25) clearly indicate that in both
hierarchies we obtain similar amounts of muons, but the inverted hierarchy
produces many more electrons.

The dijet bounds depend very weakly on the hierarchy. For large $\kappa_{s}$
this is unsurprising since dijet decays dominate the total width, and the
branching ratio into dijets will depend very little on the details of the
neutrino hierarchy. However, for small $\kappa_{s}$ the partial widths into
neutrinos can dominate the total width. As Eqs. (20-22,23-25,27) indicate, the
widths into neutrinos depend strongly on the neutrino hierarchy. In this case,
the production cross section times branching ratio into jets is schematically

$\displaystyle\sigma(pp\rightarrow H^{\pm}/H/A){\rm BR}(H^{\pm}/H/A\rightarrow
jj)\sim\frac{\kappa_{s}^{4}}{\Gamma_{\nu}},$ (28)

where $\Gamma_{\nu}$ is a generic partial width into neutrino final states. At
a fixed mass, there is an upper bound on this cross section times branching
ratio. However, even a factor of two change in $\Gamma_{\nu}$ can be
compensated by a $2^{1/4}\approx 1.2$ change in $\kappa_{s}$. That is, a large
change in the branching fraction into neutrinos can be compensated by a small
$\mathcal{O}(10\%)$ change in $\kappa_{s}$. In the region where this
important, $\kappa_{s}\ll 1$, this change is not visible on the plots.

##  V Summary and Conclusions

In this work, we considered the possibility that the small neutrino masses
deduced from flavor oscillation experiments originate from the QCD quark
condensate. We focused on the case of Dirac neutrinos, where the requisite
effective interaction is mediated by a dimension-6 operator. The required
coupling to light quarks makes hadron colliders natural settings for
investigation of the underlying model. We showed that charged meson decay
constraints can be avoided if strange quarks provide the dominant effect. This
scenario can be realized with one additional Higgs doublet at the TeV scale,
giving rise to the dimension-6 interactions of quarks, leptons, and right-
handed neutrinos. We described how the structure of the heavy doublet Yukawa
couplings may be realized by appropriate UV model building.

Flavor data from $D-\bar{D}$ meson mixing constrain the allowed parameter
space of the model. However, much of the remaining parameter space can be
tested through resonant production of heavy Higgs states from $s\bar{s}$ and
$c\bar{s}/s\bar{c}$ initial states at a future 100 TeV hadron collider. Both
dijet and lepton + missing energy final states can in general be accessible,
for generic choices of parameters. An interesting and uncommon aspect of our
model is that measurements of the heavy charged Higgs branching ratios into
the lepton flavor final states can distinguish between the underlying “normal”
and “inverted” neutrino mass hierarchies. Specifically, in our model leptonic
charged Higgs decays are dominated by muons in the case of normal hierarchy,
followed by taus and electrons; for inverted hierarchy the charged lepton
ordering is reversed. Thus, collider probes of our scenario can be
complementary to corresponding tests of the mass ordering at future neutrino
oscillation experiments, and shed light on the microscopic neutrino mass
generation mechanism.

###### Acknowledgements.

Acknowledgments: The work of H.D. and M.S. is supported by the United States
Department of Energy under Grant Contract DE-SC0012704. I.M.L. is supported in
part by the United States Department of Energy grant number DE-SC0017988. The
data to reproduce the plots are available upon request.

## References

  * Thomas and Xu (1992) S. D. Thomas and R.-M. Xu, Phys. Lett. B 284, 341 (1992).
  * McDonald (1996) J. McDonald (1996), eprint hep-ph/9610324.
  * Ibanez et al. (2001) L. E. Ibanez, F. Marchesano, and R. Rabadan, JHEP 11, 002 (2001), eprint hep-th/0105155.
  * Davoudiasl and Everett (2006) H. Davoudiasl and L. L. Everett, Phys. Lett. B 634, 55 (2006), eprint hep-ph/0512188.
  * Babič et al. (2021) A. Babič, S. Kovalenko, M. I. Krivoruchenko, and F. Šimkovic, Phys. Rev. D 103, 015007 (2021), eprint 1911.12189.
  * Davies et al. (2019) C. T. H. Davies, K. Hornbostel, J. Komijani, J. Koponen, G. P. Lepage, A. T. Lytle, and C. McNeile (HPQCD), Phys. Rev. D 100, 034506 (2019), eprint 1811.04305.
  * Zyla et al. (2020) P. Zyla et al. (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Cirigliano and Rosell (2007) V. Cirigliano and I. Rosell, Phys. Rev. Lett. 99, 231801 (2007), eprint 0707.3439.
  * Bryman et al. (2011) D. Bryman, W. J. Marciano, R. Tschirhart, and T. Yamanaka, Ann. Rev. Nucl. Part. Sci. 61, 331 (2011).
  * Bryman et al. (2021) D. Bryman, V. Cirigliano, A. Crivellin, and G. Inguglia (2021), eprint 2111.05338.
  * Davoudiasl et al. (2020) H. Davoudiasl, I. M. Lewis, and M. Sullivan, Phys. Rev. D 101, 055010 (2020), eprint 1909.02044.
  * Davoudiasl et al. (2021) H. Davoudiasl, I. M. Lewis, and M. Sullivan, Phys. Rev. D 104, 015024 (2021), eprint 2103.12089.
  * Hannestad (2004) S. Hannestad, Phys. Rev. D 70, 043506 (2004), eprint astro-ph/0403291.
  * Zupanc et al. (2013) A. Zupanc et al. (Belle), JHEP 09, 139 (2013), eprint 1307.6240.
  * Altmannshofer et al. (2019) W. Altmannshofer et al. (Belle-II), PTEP 2019, 123C01 (2019), [Erratum: PTEP 2020, 029201 (2020)], eprint 1808.10567.
  * Egana-Ugrinovic et al. (2019a) D. Egana-Ugrinovic, S. Homiller, and P. Meade, Phys. Rev. Lett. 123, 031802 (2019a), eprint 1811.00017.
  * Egana-Ugrinovic et al. (2019b) D. Egana-Ugrinovic, S. Homiller, and P. R. Meade, Phys. Rev. D 100, 115041 (2019b), eprint 1908.11376.
  * Chanowitz et al. (1979) M. S. Chanowitz, M. A. Furman, and I. Hinchliffe, Nucl. Phys. B 153, 402 (1979).
  * Ivanov et al. (2020) M. M. Ivanov, M. Simonović, and M. Zaldarriaga, Phys. Rev. D 101, 083504 (2020), eprint 1912.08208.
  * Aaij et al. (2021) R. Aaij et al. (LHCb), Phys. Rev. Lett. 127, 111801 (2021), eprint 2106.03744.
  * Bazavov et al. (2018) A. Bazavov et al. (Fermilab Lattice, MILC, TUMQCD), Phys. Rev. D 98, 054517 (2018), eprint 1802.04248.
  * Alexandrou et al. (2020) C. Alexandrou, J. Finkenrath, L. Funcke, K. Jansen, B. Kostrzewa, F. Pittler, and C. Urbach, Phys. Rev. Lett. 125, 232001 (2020), eprint 2002.07802.
  * Esteban et al. (2020) I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz, and A. Zhou, JHEP 09, 178 (2020), eprint 2007.14792.
  * Aad et al. (2020) G. Aad et al. (ATLAS), JHEP 03, 145 (2020), eprint 1910.08447.
  * Sirunyan et al. (2020) A. M. Sirunyan et al. (CMS), JHEP 05, 033 (2020), eprint 1911.03947.
  * Aad et al. (2019) G. Aad et al. (ATLAS), Phys. Rev. D 100, 052013 (2019), eprint 1906.05609.
  * CMS (2021) Tech. Rep., CERN, Geneva (2021), URL https://cds.cern.ch/record/2758647.
  * Alwall et al. (2014) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro, JHEP 07, 079 (2014), eprint 1405.0301.
  * Christensen and Duhr (2009) N. D. Christensen and C. Duhr, Comput. Phys. Commun. 180, 1614 (2009), eprint 0806.4194.
  * Alloul et al. (2014) A. Alloul, N. D. Christensen, C. Degrande, C. Duhr, and B. Fuks, Comput. Phys. Commun. 185, 2250 (2014), eprint 1310.1921.
  * Peskin and Takeuchi (1990) M. E. Peskin and T. Takeuchi, Phys. Rev. Lett. 65, 964 (1990).
  * Peskin and Takeuchi (1992) M. E. Peskin and T. Takeuchi, Phys. Rev. D 46, 381 (1992).
  * Barbieri et al. (2006) R. Barbieri, L. J. Hall, and V. S. Rychkov, Phys. Rev. D 74, 015007 (2006), eprint hep-ph/0603188.
  * Haber and O’Neil (2011) H. E. Haber and D. O’Neil, Phys. Rev. D 83, 055017 (2011), eprint 1011.6188.
  * Ahriche et al. (2015) A. Ahriche, G. Faisel, S.-Y. Ho, S. Nasri, and J. Tandean, Phys. Rev. D 92, 035020 (2015), eprint 1501.06605.
  * Golling et al. (2017) T. Golling et al., CERN Yellow Rep. pp. 441–634 (2017), eprint 1606.00947.


# Topological State Reconstruction For Wireless Stabilization of Distant
Atomic Clocks

Adam Walton and Benjamin T. H. Varcoe<EMAIL_ADDRESS>University of Leeds

###### Abstract

High-precision frequency alignment with classical communication channels is
difficult due to noise, propagation delays, and signal degradation. Current
optical methods [1], commonly involving frequency combs, are capable of
synchronising clocks with exceptional precision up to the region of a part in
$10^{20}$. Alternatively, wireless methods see use where this is not
practical, with achievable precision within the nanosecond region [2]. This
leaves few options for achieving high-precision clock synchronisation without
requiring specialised equipment, a fibre connection, or a line of sight
communication channel. Here we present a novel approach combining quantum
state reconstruction with feedback controls to stabilize the frequency of two
atomic clocks separated by a 900 MHz free space radio link. Quantum state
reconstruction enables tracking of phase and frequency fluctuations during
transmission. We see that a part in $10^{16}$ precision in frequency alignment
of the clocks can be achieved using commonly-available radio equipment,
allowing precise timekeeping and synchronization over long distances provided
a radio communications channel can be established, with potential applications
in a wide variety of timekeeping applications.

††preprint: APS/123-QED

Accurate time measurement and synchronization are essential in various
scientific and technological applications [3, 4]. Hence, atomic clocks have
proven to be extremely valuable in maintaining precise timekeeping across a
wide range of applications, from fundamental physics to navigation.
Synchronizing atomic clocks has therefore also become a crucial problem in
various physics-related fields, such as geodesy [5, 6, 7] and fundamental
physics measurements [8, 9, 10]. Additionally, it plays a vital role in
national timekeeping and providing accurate timing information for GPS. In the
UK, the National Physical Laboratory (NPL) offers a service with clocks
capable of detecting frequency differences on the order of one part in
$10^{15}$ over a day [11].

Precise frequency distribution can be achieved through optical methods
provided that placement is practical for fibre cables [12] or line-of-sight
[13, 14] is available, but for longer distances, wireless methods like Two-Way
Satellite Time and Frequency Transfer (TWSTFT) [15, 16] and Global Navigation
Satellite System (GNSS) [17] are employed. However, these wireless methods are
less precise, providing time measurements accurate to the nanosecond [11].

In this paper we introduce an innovative approach that combines tomographic
state reconstruction with a phase-locked loop (PLL) to stabilize the frequency
of two atomic clocks separated by a 900 MHz free-space radio link. Tomographic
reconstruction allows the retrieval of information about the phase and
amplitude of electromagnetic states, even in the presence of inherent noise
and measurement limitations. The aim is to detect phase errors in the
broadcast channel between two parties using imperfect atomic clocks as
frequency references. These phase errors stem from frequency differences
between the clocks. By employing these errors as feedback for a controller, we
can stabilize the clocks and achieve clock synchronization with impressive
precision, achieving drifts below the nano-Hertz level. Such corrections
already happen in ordinary communications systems and form a normal part of
device operation, though these corrections are a part of signal processing and
do not make changes to the clock source. The novel feature introduced here is
to correct the frequency of an atomic clock, effectively treating the clock as
a flywheel oscillator and thereby allowing the frequency reference to be used
by other elements of a system. In particular this would allow for the creation
of a wireless ”White rabbit” type protocol,

The conventional synchronization of atomic clocks using classical
communication channels encounters issues related to noise, propagation delays,
and signal degradation. In contrast, combining a topological state
reconstruction with a PLL presents a promising solution to overcome these
challenges. The PLL is a widely used technique for generating a stable output
frequency by comparing it with a reference signal. By incorporating quantum
state reconstruction into the PLL, we harness the advantages of quantum
measurement techniques to enhance stability and accuracy during the frequency-
locking process. This integration opens up new possibilities for achieving
precise synchronization between atomic clocks over long distances with
unprecedented accuracy.

Topological state reconstruction combined with PLL stabilization of atomic
clocks through a free space radio link operating at higher frequencies, such
as 2.4 GHz, holds tremendous potential for time synchronization across a
network. This approach allows compensation for phase and frequency
fluctuations induced by the transmission medium and environmental factors,
resulting in improved accuracy and stability. Highly precise synchronization
over long distances is crucial for various applications, including global
positioning systems, telecommunications, and fundamental scientific research.
This paper offers a comprehensive analysis of the proposed system,
highlighting its advantages, potential challenges, and experimental results,
all of which demonstrate the feasibility and practicality of this approach for
future implementations.

Figure 1: A comparison of an interferometer used for topological state
reconstruction (top) to IQ demodulation used for radio signals (bottom). At
$\phi=\frac{\pi}{2}$, topological state reconstruction is equivalent to IQ
demodulation in radio broadcasting.

We present a measurement of the phase of a 900 MHz microwave link using a
homodyne POVM measurement that simultaneously extracts the phase and amplitude
of a microwave channel. We use this to align the phase of separated atomic
clocks. We demonstrate this technique over a free space broadcast channel,
achieving a fractional stability of $10^{-16}$ range.

In Homodyne detection, the signal is split into two parts. Each component is
separately detected in a balanced detector measuring the input electric fields
$E_{s}=A_{s}(t)\exp(i\omega_{s}t)$ and $E_{l}=A_{l}\exp(i\omega_{l}t)$ via the
combinations

$E_{1}=\frac{1}{2}(E_{s}+E_{l}),~{}E_{2}=\frac{1}{2}(E_{s}-E_{l})$
$E_{3}=\frac{1}{2}(E_{s}+iE_{l}),~{}E_{4}=\frac{1}{2}(E_{s}-iE_{l})$

The $n$th receiver output intensity is therefore given by

$\displaystyle I_{n}(t)=$ $\displaystyle R~{}Re$
$\displaystyle\left\\{\frac{A_{s}(t)exp(i\omega_{s}t)+(\pm
1)^{n}A_{l}(t)exp(i\omega_{l}t)}{\sqrt{2}}\right\\}$

$\displaystyle I_{1}(t)$ $\displaystyle=R[P_{s}(t)+P_{l}$
$\displaystyle+2\sqrt{P_{s}(t)P_{l}}\cos(\omega_{IF}+\theta_{s}(t)-\theta_{l})]$
$\displaystyle I_{2}(t)$ $\displaystyle=R[P_{s}(t)+P_{l}$
$\displaystyle-2\sqrt{P_{s}(t)P_{l}}\cos(\omega_{IF}+\theta_{s}(t)-\theta_{l})]$

where $\omega_{IF}=\omega_{s}-\omega_{l}$, $R$ is a normalisation constant,
$P_{s}$ and $P_{l}$ are respectively the power of the signal and local
oscillator and $\theta_{s}$ and $\theta_{l}$ are the phase of the signal and
local oscillator. Typically we make the assumptions that the local oscillator
phase and power are both stable and unchanging, removing time dependence.

The output of the balanced detector is:

$I_{I}(t)=I_{1}(t)-I_{2}(t)$ $\therefore
I_{I}(t)=R\sqrt{P_{s}(t)P_{l}}\cos(\omega_{IF}+\theta_{s}(t)-\theta_{l})]$

Inserting a $\pi/2$ phase shift into the local oscillator delivers the output
of the other quadrature (Fields 3 and 4):

$I_{Q}(t)=I_{3}(t)-I_{4}(t)$ $\therefore
I_{Q}(t)=R\sqrt{P_{s}(t)P_{l}}\sin(\omega_{IF}+\theta_{s}(t)-\theta_{l})]$

These two intensities, $I_{I}(t)$ and $I_{Q}(t)$, can be viewed as coordinates
on the unit circle, allowing us to uniquely determine the phase of the local
oscillator.

In two-bit digital transmission, the phase $\theta_{s}$ is restricted to one
of two values $\theta_{s}\in{0,\pi}$. This is known as binary phase-shift
keying (BPSK), the pattern of shifts is called a ‘cluster’. The benefit of
this keying method is that we are able to accurately determine the location of
the origin as the center of the cluster and thereby accurately determine the
relative phase between the sender and receiver. Under this circumstance, the
frequency $\omega_{IF}$ will be observed as a rotation of the cluster via the
relationship,

$\frac{d\theta_{s}}{dt}=\omega_{IF}.$

The experimental method that we have used to synchronise the atomic clocks is
to draw a feedback signal from the cluster rotation. The phase-locked loop
detects this rotation and feeds this information back to the atomic clock,
adjusting the frequency (and phase) to match that of the transmitter.

Our approach draws parallels with PLLs that are already widely used in digital
broadcast channels to ensure consistent frequencies between transmitters and
receivers, facilitating smooth radio broadcasting operations. However, a
crucial distinction lies in two areas: firstly, our use of frequency scaling,
where, rather than relying on the 10 MHz reference provided by an atomic
clock, we scale up the frequency to 900MHz. Secondly the aim of the PLL is to
stabilise an external clock. This allows us to extend the reach of a single
master oscillator to a network of local oscillators.

The rationale behind this frequency scaling stems from the assumption that any
error in the 10 MHz reference is a result of changes in the time measurement
of the clocks. As the frequency increases, the phase shift for a given time
error also increases. As we are sensitive to phase, multiplying the frequency
amplifies the time error, granting us the ability to achieve proportionally
higher precision for frequency distribution. This enhancement in precision
holds great promise for various applications that require accurate
synchronization over long distances.

The broadcasts to be performed consist of a pair of thermal distributions. For
a plot in phase space, characterised by:

$\Psi\left(t\right)=X\cos\left(2\pi ft\right)+Pi\sin\left(2\pi ft\right),$ (1)

the measurement distribution is a pair of clusters centered on two points:

$\Psi\left(t\right)=A\cos\left(2\pi
ft+n\pi\right)~{}~{}~{}~{}~{}~{}n\in\\{0,1\\},$ (2)

as shown in Figure 2. A phase error is visible in this Figure, which shows the
clusters rotated from an ideal position on the X-axis, caused by differences
of frequencies between clocks used as references for the broadcasting and
receiving parties. We omit typical forms of phase error correction used in
broadcasting, as these are methods of signal processing which have no effect
on the reference clocks. Using a feedback controller instead to adjust the
frequency of a secondary clock in response to phase errors allows for clock
synchronisation.

Figure 2: A sample measurement distribution for a broadcast provided by
GNURadio [18]. A phase difference between the measurement cluster centers and
the X-axis due to the rotations that this experiment is based on is clearly
visible.

Taking repeat measurements over the duration of a broadcast between a
transmitter and receiver, using a primary and secondary clock respectively, we
are able to reconstruct the state as frequency errors cause the measurement
clusters to rotate about a centre point between them. For a difference in
frequency $\Delta f$ between the two reference clocks, we see a rate of change
of phase $\frac{\Delta\phi}{\Delta t}$ given by:

$\frac{\Delta\phi}{\Delta t}=2\pi\Delta f.$ (3)

For an ideal, perfectly synchronised pair of clocks, we expect no change in
phase over time, therefore we are limited by the minimum detectable phase
difference. Using the atomic clocks as references for a higher frequency
broadcast allows measurement of changes in broadcast frequencies instead,
which magnifies the phase errors, reducing the minimum size of phase errors
that can be observed. For some broadcast with a frequency scaled up to
$f_{Broadcast}=a\times f_{clock}$, where $f_{clock}=10$ MHz, we instead see a
rate of phase drift of:

$\frac{\Delta\phi_{a}}{\Delta t}=2\pi\Delta f_{Broadcast}=2\pi a~{}\Delta
f_{clock}.$ (4)

Given that radio equipment can be used to broadcast signals far in excess of
10 MHz, this provides an easily adjustable variable which can be increased to
improve the minimum detectable frequency errors between a pair of clocks. The
limiting factor of this method then becomes the inherent stability of the
clocks, and the responsiveness of the feedback controller used to correct the
phase errors.

## The Wired Channel

Figure 3: Tracking the average phase of one cluster in a 70 MHz broadcast. The
feedback controller was disabled for a period to demonstrate it’s effects.

A test with a wired 70 MHz broadcast is shown in Figure 3. The Figure shows
the effect of the feedback controller on the phase of the reconstructed state,
in which the broadcast is performed with the feedback controller for 120
seconds before the controller is disabled for an additional 120 seconds.
During this period, phase errors accumulated due to lingering frequency
differences. Enabling the controller again readjusts the secondary clock
frequency, rotating the measurement cluster back towards a mean of zero phase,
where it is maintained. At this point the clocks are once again stabilised.

## The Wireless Channel

We then replace the wired channel between the pair of radios with a 2.5 meter
wireless channel, and increase the broadcast frequency to 900 MHz. This is an
increase in frequency of a factor of 12.9 compared to the 70 MHz wired test,
and the drift in phase caused by differences in the reference frequencies
increases by the same factor and the gain of the feedback loop therefore also
increases. The results of this are shown in Figure 4, accompanied by an Allan
deviation plot derived from the same data. These results are shown alongside
the measurements from a wired broadcast at 70 MHz for comparison. In this
figure you can see the increased phase noise arising from the increase in
frequency, however for a longer averaging duration the Allan variance was not
noticeably affected by the change from wired to short range wireless
broadcasting. In both cases, a small peak is seen in the Allan variance plot,
this is likely displaying the delay before feedback from the controller
affects the broadcast. A linear fit shows that the total average rate of
change of phase for the wireless broadcast is estimated to be $3.609\times
10^{-06}\pm 1.192\times 10^{-7}$, which corresponds to a frequency error of
574.3 nHz. Given the source frequency of 900 MHz this corresponds to a free
space fractional frequency drift $\frac{\Delta f}{f}=6.382\times 10^{-16}$
accumulated over a 5 hour scan time.

Using this wireless data we can also plot the modified Allan variance showing
$1.701\times 10^{-15}$ at 4136 seconds. These measurements are an improvement
on the Allan variance provided in the datasheet for the PRS10 clock used in
this test for most averaging times, with the Allan variance never exceeding a
part in $10^{13}$ even at the clock’s ideal averaging period [19]. A longer
broadcasting time would be required in order to locate the point at which the
Allan variance no longer improves, though this will also be dependent on clock
performance and the chosen frequency. While the wired broadcast did show an
improved Allan deviation compared to the wireless broadcast for a longer
averaging duration, demonstrating wireless synchronisation is of greater
relevance.

Synchronisation was also achieved at 2.4 GHz, a frequency commonly used in
communications equipment, though the increased sensitivity to errors at higher
frequencies meant that the instability of the clocks used caused less reliable
broadcasts. Clocks with improved short-term stability would make higher
frequency testing more reliable.

Figure 4: Tracking the average phase of one cluster in a 900 MHz wireless
broadcast, and a 70 MHz wired broadcast. Due to the higher frequency, greater
instability in the phase drift is visible in the wireless broadcast than
wired. The Allan variance for both broadcasts are also shown.

This method only relies on radio broadcasts, which are already employed over
long distances and at a wide array of frequencies due to their ubiquitous use
in communications, and a suitably responsive feedback controller, which are
widely used in a variety of applications and are well understood. This
provides clear avenues though which the experimental setup employed here may
be improved, through use of broadcasting equipment suitable for longer ranges,
and through increasing the broadcast frequency to magnify the effects of clock
frequency errors. A clear next step would also be to perform tests with more
stable clocks, ideally focusing on a superior primary clock. This is because
instability in the primary and secondary clocks are not differentiated between
with this method, and are both ascribed to errors in the secondary clock. This
may result in errors in stability estimates.

Fractional error estimates for frequency measurements in the order of
$10^{-16}$ were observed, an improvement of approximately a factor of $10^{3}$
over the minimum advertised Allan variance for the clocks used. This was
performed across 2.5m of free space using easily accessible radio equipment,
and offers obvious avenues for improvement through longer ranges, higher
frequency broadcasts, as well as replacement of the PRS10C clocks with time
references that have superior short-term stability.

## Author Information

Adam Walton and Benjamin T. H. Varcoe: These authors contributed equally to
the work.

## References

  * Tian  _et al._ [2017] H. Tian, Y. Song, J. Yu, H. Shi, and M. Hu, IEEE Photonics Journal 9, 1 (2017).
  * GPS [2022] Gps accuracy, https://www.gps.gov/systems/gps/performance/accuracy/ (2022), accessed: 15/12/23.
  * Gozzard [2023] D. Gozzard, Nature 618, 680 (2023).
  * Ludlow  _et al._ [2015] A. D. Ludlow, M. M. Boyd, J. Ye, E. Peik, and P. O. Schmidt, Rev. Mod. Phys. 87, 637 (2015).
  * McGrew  _et al._ [2018] W. F. McGrew, X. Zhang, R. J. Fasano, S. A. Schäffer, K. Beloy, D. Nicolodi, R. C. Brown, N. Hinkley, G. Milani, M. Schioppo, T. H. Yoon, and A. D. Ludlow, Nature 564, 87 (2018).
  * Mehlstäubler  _et al._ [2018] T. E. Mehlstäubler, G. Grosche, C. Lisdat, P. O. Schmidt, and H. Denker, Reports on Progress in Physics 81, 64401 (2018).
  * Riehle [2017] F. Riehle, Nature Photonics 11, 25 (2017).
  * G. Barontini [2022] V. B. G. Barontini, L. Blackburn, EPJ Quantum Technology 9, 10.1140/epjqt/s40507-022-00130-5 (2022).
  * Kolkowitz  _et al._ [2016] S. Kolkowitz, I. Pikovski, N. Langellier, M. D. Lukin, R. L. Walsworth, and J. Ye, Phys. Rev. D 94, 124043 (2016).
  * Wcisło  _et al._ [2016] P. Wcisło, P. Morzyński, M. Bober, A. Cygan, D. Lisak, R. Ciuryło, and M. Zawada, Nature Astronomy 1, 0009 (2016).
  * NPL [2023] Time scales, https://www.npl.co.uk/time-frequency/time-scales (2023), accessed: 04/07/23.
  * Cizek  _et al._ [2022] M. Cizek, L. Pravdova, T. M. Pham, A. Lesundak, J. Hrabina, J. Lazar, T. Pronebner, E. Aeikens, J. Premper, O. Havlis, R. Velc, V. Smotlacha, L. Altmannova, T. Schumm, J. Vojtech, A. Niessner, and O. Cip, Opt. Express 30, 5450 (2022).
  * Caldwell  _et al._ [2023] E. Caldwell, J.-D. Deschênes, J. Ellis, W. Swann, B. Stuhl, H. Bergeron, N. Newbury, and L. Sinclair, Nature 618, 721 (2023).
  * Deschênes  _et al._ [2016] J.-D. Deschênes, L. C. Sinclair, F. R. Giorgetta, W. C. Swann, E. Baumann, H. Bergeron, M. Cermak, I. Coddington, and N. R. Newbury, Phys. Rev. X 6, 021016 (2016).
  * Guo  _et al._ [2022] Y. Guo, S. Gao, Y. Bai, Z. Pan, Y. Liu, X. Lu, and S. Zhang, Remote Sensing 14, 528 (2022).
  * Piester  _et al._ [2011] D. Piester, M. Rost, M. Fujieda, T. Feldmann, and A. Bauch, Advances in Radio Science 9, 1 (2011).
  * GNS [2023] Gnss and satellite time transfer, https://www.nist.gov/pml/time-and-frequency-division/time-services/gnss-and-satellite-time-transfer (2023), accessed: 05/07/23.
  * GNU [2023] Gnuradio, https://www.gnuradio.org (2023), accessed: 09/06/23.
  * PRS [2023] Prs10 — rubidium frequency standard with low phase noise, https://www.thinksrs.com/downloads/pdfs/catalog/PRS10c.pdf (2023), accessed: 18/07/23.


# The appearance function for paper-folding words

Rob Burns

###### Abstract

We provide a complete characterisation of the appearance function for paper-
folding sequences for factors of any length. We make use of the software
package Walnut to establish these results.

##  1 Introduction

The regular paper-folding sequence begins  1, 1, -1, 1, 1, -1, -1, 1, …. It is
derived from the hills and valleys created when a piece of paper is folded
length-wise multiple times. In the limit, it consists of an infinite sequence
of $1$’s and $-1$’s, in which a $1$ corresponds to a hill and a $-1$ to a
valley in the unfolded paper. This sequence, with $-1$’s replaced by $0$’s,
appears as sequence A014577 in the On-Line Encyclopedia of Integer Sequences
(OEIS).[1] In a more general form, introduced by Davis and Knuth,[2] a paper-
folding sequence is derived from an infinite folding instruction set. This set
of instructions also consists of an infinite sequence of $1$’s and $-1$’s and
determines the way in which the paper is folded. The regular paper-folding
sequence is produced by an instruction set consisting of an infinite sequence
of $1$’s.

In this paper, we study the appearance function of the paper-folding
sequences. We show that, when $n\geq 7$, the appearance function is determined
in a simple way by the folding instruction set. The connection between the
folding instructions and the appearance function for smaller values of $n$ is
not quite as simple, but can still be described.

We make use of the software package Walnut. Hamoon Mousavi, who wrote the
program, has provided an introductory article [7]. Papers that have used
Walnut include [8], [11], [9], [6], [4]. Further resources related to Walnut
can be found at Jeffrey Shallit’s page.

The free open-source mathematics software system SageMath [12] was used to
check some of the calculations.

##  2 Background and notation

Let $w$ be an infinite word. The first element of $w$ is denoted by $w[1]$,
the second element by $w[2]$ etc.. A sub-word of $w$ is a continuous set of
elements contained within $w$. The sub-word $w[i],w[i+1],w[i+2],\dots,w[j]$
will be abbreviated to $w[i:j]$. A finite word is called a factor of $w$ if it
appears as a sub-word of $w$. The length $k$ prefix of $w$ is the length $k$
sub-word of $w$ starting with the element $w[1]$, i.e. the subword $w[1:k]$.

The appearance function of $w$, denoted $A_{w}(n)$, is defined to be the least
integer k such that a copy of each length $n$ factor of $w$ is contained in
the prefix $w[1:k]$. For convenience we will use a related function which we
call $S_{w}(n)$. $S_{w}(n)$ is defined to be the least integer k such that a
copy of each length $n$ factor of $w$ starts somewhere within the prefix
$w[1:k]$. The two functions are connected by the equation
$A_{w}(n)=S_{w}(n)+n-1$.

The set of folding instructions associated with a paper-folding sequence will
be denoted by $f=(f_{0},f_{1},f_{2},\dots)$. The paper-folding sequence
associated to the folding instructions $f$ will be denoted by
$P_{f}=P_{f}[1],P_{f}[2],\dots$. The appearance function of $P_{f}$ will be
abbreviated to $A_{f}$ and $S_{P_{f}}$ will be abbreviated to $S_{f}$.

Dekking, Mendés France and van der Poorten [3] showed that the value of
$P_{f}[k]$ can be written in terms of $f$ in a fairly simple way. If
$k=2^{s}\cdot r$ where $r$ is odd, then

$P_{f}[k]=\begin{cases}f_{s},&\text{if }\,\,r\equiv 1\pmod{4}\\\
-f_{s},&\text{if }\,\,r\equiv 3\pmod{4}\end{cases}$ (1)

Schaeffer [10] observed that equation (1) leads to a 5-state deterministic
finite automaton that takes, as input the base-2 expansion of an integer $k$
in parallel with the folding instructions $f$ and outputs $P_{f}[k]$. The
automaton outputs the correct value of $P_{f}[k]$ provided that enough folding
instructions have been read in. In particular, more than $\log_{2}(k)$ folding
instructions must be included in the input. The Walnut software package
includes this automaton and can be used to investigate its behaviour.

For integers $n$, define the function $\phi(n)$ to be the least integer $r$
such that $r\geq n$ and $r$ is a power of $2$. So, if $2^{k-1}<n\leq 2^{k}$,
then $\phi(n)=2^{k}$. An alternative definition is that

$\phi(n)=2^{k},\text{ where }k=\lceil\log_{2}(n)\rceil.$

Schaeffer [10] showed that, when $n\geq 3$:

$\max_{f}S_{f}(n)=6\cdot\phi(n).$ (2)

Goč et al. [5] showed that when $n\geq 7$,

$\min_{f}S_{f}(n)=4\cdot\phi(n).$ (3)

##  3 Formula for $S_{f}$ and $A_{f}$

We begin this section with some Walnut commands. Walnut includes an automaton
which takes as parallel input a folding instruction set $f$ and the base-2
representation of an integer $k$, written in least significant digit (lsd)
first order, and outputs the value of $P_{f}[k]$. The Walnut representation of
$P_{f}[k]$ is $PF[f][k]$. We now introduce some Walnut formulae which will be
useful. Firstly we define the automaton pffaceq which takes as parallel input
the folding instruction set $f$ and three integers $i$, $j$ and $n$, written
in base-2 lsd format. The resulting automaton accepts the input if and only if
the length $n$ subwords of $P_{f}$ starting at indices $i$ and $j$ are
identical. Since it has 153 states, it cannot be displayed here.

    def pffaceq  "?lsd_2 Ak (k < n) => PF[f][i+k] = PF[f][j+k]":

The following code creates an automaton related to the function $\phi$ which
was defined in section 2. The automaton takes two integers $x$ and $y$ as
input, written in base-2 lsd form. It accepts the input if $x$ is a power of
$2$ and $\phi(y)=x$. We use the name pfphi for this automaton. It is pictured
in figure 1.

    reg power2 lsd_2 "0*10*":
    def pfphi "?lsd_2 $power2(x) & (x >= y) & x < 2*y":

Figure 1: Automaton pfphi.

We now create an automaton pfapp which takes as parallel input a folding
instruction set $f$ and two integers $i$ and $n$, written in base-2 lsd
format. The resulting automaton accepts the input if and only if the length
$n$ factor of $P_{f}$ starting at index $i$ does not appear earlier within
$P_{f}$. This automaton has $121$ states.

    def pfapp "?lsd_2 (Aj (j<i) => (Et t<n & PF[f][i+t] != PF[f][j+t]))":

We next establish some preliminary results.

######  Lemma 3.1.

When $n\geq 7$, the length $n$ factor $P_{f}[6\cdot\phi(n):6\cdot\phi(n)+n-1]$
first appears in $P_{f}$ starting at either index $4\cdot\phi(n)$ or index
$6\cdot\phi(n)$. This factor appears nowhere else in the prefix
$P_{f}[1:6\cdot\phi(n)+n-1]$ of $P_{f}$.

###### Proof.

We create an automaton which takes as input the folding instructions $f$ and
an integer $n$ written in base-2 lsd form. It accepts the input if there is an
index $k<6\cdot\phi(n)$ such that $k\neq 4\cdot\phi(n)$ and the two factors
$P_{f}[6\cdot\phi(n):6\cdot\phi(n)+n-1]$ and $P_{f}[k:k+n-1]$ are identical.

    eval pftemp "?lsd_2 (n >= 7) & (Ex, k x >= 1 & $pfphi(x,n) &
       (k != 4*x) & (k<6*x) & (Ai (i<n) => PF[f][k+i]=PF[f][6*x+i]))":

The automaton accepts no input showing that the factor
$P_{f}[6\cdot\phi(n):6\cdot\phi(n)+n-1]$ first appears either at index
$4\cdot\phi(n)$ or $6\cdot\phi(n)$ and appears nowhere else in the prefix
$P_{f}[1:6\cdot\phi(n)+n-1]$. ∎

######  Lemma 3.2.

When $n\geq 7$, the factor $P_{f}[6\cdot\phi(n):6\cdot\phi(n)+n-1]$ is always
the last length $n$ factor to appear in $P_{f}$.

###### Proof.

We know from (2) that no length $n$ factor of $P_{f}$ can first begin later
than index $6\cdot\phi(n)$. We create an automaton which takes as parallel
input a folding instruction set $f$ and two integers $i$ and $n$, written in
base-2 lsd format. The automaton accepts the input if and only if the length
$n$ factors beginning at $4\cdot\phi(n)$ and $6\cdot\phi(n)$ are identical,
$4\cdot\phi(n)<i<6\cdot\phi(n)$ and the factor $P_{f}[i:i+n-1]$ appears no
earlier within $P_{f}$.

    eval pftemp "?lsd_2 (Ex (x >= 1) & $pfphi(x,n) &
       $pffaceq(f,4*x,6*x,n) & (Ei (i>4*x) & (i<6*x) & $pfapp(f,i,n)))":

The automaton accepts no input. So, if the factor
$P_{f}[6\cdot\phi(n):6\cdot\phi(n)+n-1]$ first appears starting at
$6\cdot\phi(n)$ then it is the last factor to appear in $P_{f}$ because of
(2). If, on the other hand, it first appears starting at index $4\cdot\phi(n)$
(the only other possibility due to lemma 3.1), it is again the last factor to
appear, otherwise pftemp would accept some input. ∎

######  Lemma 3.3.

Let $n\geq 7$ with $2^{k-1}<n\leq 2^{k}$, so $\phi(n)=2^{k}$. Then the factors
$P_{f}[6\cdot 2^{k}:6\cdot 2^{k}+n-1]$ and $P_{f}[6\cdot 2^{k}:6\cdot
2^{k}+2^{k}-1]$ first appear in $P_{f}$ at the same starting index.

###### Proof.

If the factor $P_{f}[6\cdot 2^{k}:6\cdot 2^{k}+n-1]$ first appears at index
$6\cdot 2^{k}$ then the factor $P_{f}[6\cdot 2^{k}:6\cdot 2^{k}+2^{k}-1]$ must
also first appear at index $6\cdot 2^{k}$ since $n\leq 2^{k}$. So, assume the
factor $P_{f}[6\cdot 2^{k}:6\cdot 2^{k}+n-1]$ first appears at index $4\cdot
2^{k}$. By lemma 3.1, this is the only other possible starting index. We
create an automaton which accepts the pair $f$ and $n$ when $P_{f}[6\cdot
2^{k}:6\cdot 2^{k}+n-1]$ first appears at index $4\cdot 2^{k}$ and
$P_{f}[6\cdot 2^{k}:6\cdot 2^{k}+2^{k}-1]$ does not first appear at index
$4\cdot 2^{k}$.

    eval pftemp "?lsd_2 (n >= 7) & (Ex x >= 1 & $pfphi(x,n) &
       (Ak (k<n) & PF[f][4*x+k] = PF[f][6*x+k]) &
       (Er (r<x) & (PF[f][4*x+k] != PF[f][6*x+k])))":

The automaton accepts no input showing that, when $P_{f}[6\cdot 2^{k}:6\cdot
2^{k}+n-1]$ first appears at index $4\cdot 2^{k}$, then so does the factor
$P_{f}[6\cdot 2^{k}:6\cdot 2^{k}+2^{k}-1]$. ∎

Our main result is an exact formula for $S_{f}$ (and therefore $A_{f}$).

######  Theorem 3.4.

For $n\geq 7$,

$S_{f}(n)=\begin{cases}4\cdot\phi(n),&\text{if }\,\,\phi(n)=2^{k}\text{ and
}f_{k+1}\neq f_{k+2}\\\ 6\cdot\phi(n),&\text{if }\,\,\phi(n)=2^{k}\text{ and
}f_{k+1}=f_{k+2}.\end{cases}$

###### Proof.

We start by showing that the theorem holds when $n$ is a power of $2$, so that
$\phi(n)=n$. By lemmas 3.1 and 3.2, the factor
$P_{f}[6\cdot\phi(n):6\cdot\phi(n)+n-1]$ is always the last length $n$ factor
to appear in $P_{f}$ and appears first at either index $4\cdot\phi(n)$ or
$6\cdot\phi(n)$. The following automaton takes as parallel input the folding
instructions $f$ and an integer $n$ and accepts the input if $n$ is a power of
$2$ and the last length $n$ factor to appear in $P_{f}$ (i.e. the factor
$P_{f}[6\cdot\phi(n):6\cdot\phi(n)+n-1]$) starts at index $4\cdot\phi(n)$.

    eval pfpow24 "?lsd_2 (n >= 7) & $power2(n) &
       (Ex x >= 1 & $pfphi(x,n) & $pffaceq(f,6*x, 4*x, n))":

Figure 2: Automaton pfpow24.

The automaton pfpow24 is pictured in figure 2. Remembering that the first
element of the folding instructions $f$ has index $0$, it is clear from the
picture that the pair $f$ and $n$ is accepted if and only if $n=2^{k}$ for
some $k\geq 3$ and $f_{k+1}\neq f_{k+2}$.

The next automaton takes as parallel input the folding instructions $f$ and an
integer $n$ and accepts the input if $n$ is a power of $2$ and the last length
$n$ factor to appear in $P_{f}$ starts at index $6\cdot\phi(n)$.

    eval pfpow26 "?lsd_2 (n >= 7) & $power2(n) &
       (Ex, k ((x >=1) & (k<n)) & $pfphi(x,n) &
       (PF[f][4*x+k] != PF[f][6*x+k]))":

Figure 3: Automaton pfpow26.

The automaton is pictured at figure 3. It is clear from the picture that the
pair $f$ and $n$ is accepted if and only if $n=2^{k}$ for some $k\geq 3$ and
$f_{k+1}=f_{k+2}$.

This completes the proof when $n$ is a power of $2$. The general case follows
from lemma 3.3. If $n\geq 7$ and $\phi(n)=2^{k}$, then the last length $n$
factor to appear is $P_{f}[6\cdot\phi(n):6\cdot\phi(n)+n-1]$. If
$\phi(n)=2^{k}$, lemma 3.3 says that this factor first appears at the same
index as the factor $P_{f}[6\cdot 2^{k}:6\cdot\phi(n)+2^{k}-1]$. From above,
this starting index is $4\cdot\phi(n)$ when $f_{k+1}\neq f_{k+2}$ and is
$6\cdot\phi(n)$ when $f_{k+1}=f_{k+2}$. ∎

######  Corollary 3.5.

For $n\geq 7$,

$A_{f}(n)=\begin{cases}4\cdot 2^{k}+n-1,&\text{if }\,\,2^{k-1}<n\leq
2^{k}\text{ and }f_{k+1}\neq f_{k+2}\\\ 6\cdot 2^{k}+n-1,&\text{if
}\,\,2^{k-1}<n\leq 2^{k}\text{ and }f_{k+1}=f_{k+2}.\end{cases}$

######  Corollary 3.6.

$A_{f}(n)=4\cdot\phi(n)+n-1$ for all $n\geq 7$ if and only if

$f=(f_{0},f_{1},f_{2},f_{3},1,-1,1,-1,\dots)\text{ or
}f=(f_{0},f_{1},f_{2},f_{3},-1,1,-1,1,\dots)$

where $f_{0},f_{1},f_{2},f_{3}\in\\{-1,1\\}$.

$A_{f}(n)=6\cdot\phi(n)+n-1$ for all $n\geq 7$ if and only if

$f=(f_{0},f_{1},f_{2},f_{3},1,1,1,\dots)\text{ or
}f=(f_{0},f_{1},f_{2},f_{3},-1,-1,-1,\dots)$

where $f_{0},f_{1},f_{2},f_{3}\in\\{-1,1\\}$.

##  4 What happens when $n<7$?

When the word length is less than $7$, the appearance function displays a
number of different behaviours. A calculation shows that, for folding
instructions $f$,

$\displaystyle S_{f}(1)\in\\{2,3\\}\,\,\,:$
$\displaystyle\,\,\,A_{f}(1)\in\\{2,3\\}$ $\displaystyle
S_{f}(2)\in\\{4,5,6\\}\,\,\,:$ $\displaystyle\,\,\,A_{f}(2)\in\\{5,6,7\\}$
$\displaystyle S_{f}(3)\in\\{14,16,22,24\\}\,\,\,:$
$\displaystyle\,\,\,A_{f}(3)\in\\{16,18,24,26\\}$ $\displaystyle
S_{f}(4)\in\\{14,16,22,24\\}\,\,\,:$
$\displaystyle\,\,\,A_{f}(4)\in\\{17,19,25,27\\}$ $\displaystyle
S_{f}(5)\in\\{28,32,44,48\\}\,\,\,:$
$\displaystyle\,\,\,A_{f}(5)\in\\{32,36,48,52\\}$ $\displaystyle
S_{f}(6)\in\\{31,32,47,48\\}\,\,\,:$
$\displaystyle\,\,\,A_{f}(6)\in\\{36,37,52,53\\}.$

To see how the choice of folding sequence $f$ determines $S_{f}(n)$ and
$A_{f}(n)$, we use the automaton

    eval pftemp# "?lsd_2 Ak (k < 50) => ($pfapp(f,r,#) &
       (Es (s <=r) & (At (t<#) => PF[f][k+t] = PF[f][s+t])))":,

replacing the symbol $\\#$ with the integers $\\{1,2,3,4,5,6\\}$ as
appropriate. The automaton accepts the pair $(f,r)$ when $S_{f}(\\#)=r$. We
can restrict the search to $k<50$ because we know that, for each $f$,
$S_{f}(n)$ is an increasing sequence and $S_{f}(7)=48$.

We start with the case $n=1$. The automaton

    eval pftemp1 "?lsd_2 Ak (k < 50) => ($pfapp(f,r,1) &
       (Es (s <=r) & (At (t<1) => PF[f][k+t] = PF[f][s+t])))":,

is pictured in figure 4.

Figure 4: Automaton for $S_{f}(1)$.

State $5$ is the only accepting state because of the requirement that more
than $\log_{2}(k)$ folding instructions are read into the automaton. The
automaton shows that

$\displaystyle S_{f}(1)=2\,\,\,$ $\displaystyle\text{when}\,\,\,f_{0}\neq
f_{1}$ $\displaystyle S_{f}(1)=3\,\,\,$
$\displaystyle\text{when}\,\,\,f_{0}=f_{1}.$

The automaton for $S_{f}(2)$ is displayed in figure 5. It shows that

$\displaystyle S_{f}(2)=4\,\,\,$
$\displaystyle\text{when}\,\,\,(f_{0},f_{1},f_{2})\in\\{(-1,-1,1),(-1,1,-1),(1,-1,1),(1,1,-1)\\}$
$\displaystyle S_{f}(2)=5\,\,\,$
$\displaystyle\text{when}\,\,\,(f_{0},f_{1},f_{2})\in\\{(-1,1,1),(1,-1,-1)\\}$
$\displaystyle S_{f}(2)=6\,\,\,$
$\displaystyle\text{when}\,\,\,(f_{0},f_{1},f_{2})\in\\{(-1,-1,-1),(1,1,1)\\}.$

Figure 5: Automaton for $S_{f}(2)$.

Automata for $S_{f}(n)$ when $n\in\\{3,4,5,6\\}$ are pictured in figures 6, 7,
8 and 9. In summary, formulae can be derived for $S_{f}(3)$ and $S_{f}(4)$ in
terms of $(f_{1},f_{2},f_{3},f_{4})$. A formula can be derived for $S_{f}(5)$
in terms of $(f_{1},f_{2},f_{3},f_{4},f_{5})$. Finally, a formula can be
derived for $S_{f}(6)$ in terms of $(f_{0},f_{1},f_{2},f_{3},f_{4},f_{5})$.

Figure 6: Automaton for $S_{f}(3)$. Figure 7: Automaton for $S_{f}(4)$. Figure
8: Automaton for $S_{f}(5)$. Figure 9: Automaton for $S_{f}(6)$.

Observe that the corresponding automaton for $S_{f}(7)$, which is shown in
figure 10, looks simple compared to that of smaller values of the factor
length $n$.

Figure 10: Automaton for $S_{f}(7)$.

## References

  * [1] OEIS Foundation Inc. (2022).  The on-line encyclopedia of integer sequences.
  * [2] C. Davis and D. E. Knuth.  Number representations and dragon curves – 1.  J. Recreat. Math., 3:66–81, 1970.
  * [3] F.M. Dekking, M. Mendès France, and A.J.v.d Poorten.  Folds!  Math. Intell., 4:30–138, 173–181, 190–195, 1982.
  * [4] Chen Fei Du, Hamoon Mousavi, Eric Rowland, Luke Schaeffer, and Jeffrey Shallit.  Decision algorithms for Fibonacci-automatic words, II: Related sequences and avoidability.  Theoretical Computer Science, 657:146–162, 2017.
  * [5] Daniel Goč, Hamoon Mousavi, Luke Schaeffer, and Jeffrey Shallit.  A new approach to the paperfolding sequences.  In Evolving Computability, pages 34–43. Springer International Publishing, 2015.
  * [6] Daniel Goč, Hamoon Mousavi, and Jeffrey Shallit.  On the number of unbordered factors.  In Language and Automata Theory and Applications, pages 299–310. Springer Berlin Heidelberg, 2013.
  * [7] Hamoon Mousavi.  Automatic theorem proving in Walnut.  arXiv, 2016.
  * [8] Hamoon Mousavi, Luke Schaeffer, and Jeffrey Shallit.  Decision algorithms for Fibonacci-automatic words, I: Basic results.  RAIRO Theor. Informatics Appl., 50:39–66, 2016.
  * [9] Narad Rampersad and Jeffrey Shallit.  Congruence properties of combinatorial sequences via Walnut and the Rowland-Yassawi-Zeilberger automaton.  arXiv, 2021.
  * [10] Luke Schaeffer.  Deciding properties of automatic sequences.  Master’s thesis, University of Waterloo, School of Computer Science, 2013\.
  * [11] Jeffrey Shallit.  Frobenius numbers and automatic sequences.  arXiv, 2021.
  * [12] The Sage Developers.  SageMath, the Sage Mathematics Software System (Version 8.2), 2018.  https://www.sagemath.org.