# Anisotropic magnetotransport realized in doped Hematite

E.F. Galindez-Ruales1, L. Šmejkal1,2, S. Das1, E. Baek1,3,  C. Schmitt1,  F.
Fuhrmann1,  A. Ross1,
R. González-Hernández1,4, A. Rothschild5, J. Sinova1, C.-Y. You3,  G. Jakob1,
M. Kläui1,6*

1 Institute of Physics, Johannes Gutenberg University Mainz, Staudingerweg 7,
55128 Mainz, Germany.
2 Institute of Physics, Academy of Sciences of the Czech Republic, 182 00
Prague 8, Czech Republic.
3 Department of Physics and Chemistry, DGIST, Daegu 42988, Republic of Korea.
4 Grupo de Investigacíon en Física Aplicada, Departamento de Física,
Universidad del Norte, Barranquilla, Colombia.
5 Department of Materials Science and Engineering, Technion-Israel Institute
of Technology, Haifa 32000, Israel.
6 Center for Quantum Spintronics, Norwegian University of Science and
Technology, Trondheim 7491, Norway.

*<EMAIL_ADDRESS>
###### Abstract

Conventional antiferromagnetic materials have long been recognized for their
time-reversal symmetry, resulting in a zero anomalous Hall coefficient.
However, a paradigm shift occurs when examining easy-axis antiferromagnets and
their spin-flop transition. This transition introduces a magnetic canted
moment, leading to the emergence of a non-zero anomalous Hall signal and the
generation of a non-dissipative transversal current. While high symmetry
systems typically manifest an isotropic Hall effect, our study unveils the
extraordinary behavior exhibited by hematite that becomes conductive due to
small Ti doping. We investigate the magnetotransport in Titanium-doped
hematite, uncovering a highly pronounced and unconventional symmetry. Notably,
this effect displays a remarkable dependence on the crystal orientation of the
material. We establish a compelling correlation between our experimental
observations and the predicted anomalous Hall effect in altermagnets through
symmetry analysis. This study expands our understanding of the Hall effect in
antiferromagnetic materials and sheds light on the intricate interplay between
crystal orientation and unconventional Hall phenomena.

_Keywords_ Magneto transport $\cdot$ Hall effect $\cdot$ anisotropy $\cdot$
hematite $\cdot$ antiferromagnets $\cdot$ altermagnets $\cdot$ spin-flop
transition

## Introduction

Hematite ($\alpha-$Fe2O3) is a widely studied material in the field of solid-
state physics due to its rich electronic and magnetic properties. It has
garnered significant attention in recent years, particularly in the realm of
spintronics and other emerging applications [1, 2, 3, 4, 5, 6, 7, 8, 9]. The
antiferromagnetic structure of hematite has been extensively investigated,
primarily focusing on its relationship to spin transport based on magnons [5].
At temperatures between 265 and 960 K, hematite exhibits weak ferromagnetism,
with the antiferromagnetically ordered spins slightly tilted in the basal
plane. However, below the Morin transition (mt) temperature of 265 K, the
spins reorient to align parallel to the c-axis, forming a collinear
antiferromagnetic magnetic order. This transition is attributed to the
delicate balance between magnetic-dipole and single-ion anisotropy energies
[10]. These two distinct states, referred to as the weak ferromagnetic (wfm)
and collinear antiferromagnetic (afm) states, respectively, have been
extensively investigated [11]. To tune the properties of hematite, doping has
been used [12, 13]. In particular, it was found that the introduction of low
percentages of certain dopants can induce changes in anisotropies and shift
the Morin transition temperature while having minimal impact on the overall
magnetic properties, such as the magnitude of the magnetic moment [13] and the
presence of the Morin transition (mt) [5, 14]. Additionally, hematite films
may exhibit altered Morin transition behavior due to growth-induced strain,
which can modify lattice constants and anisotropies [5].

In addition to these magnetic properties, reports in the literature have
discussed the observation of the Hall effect in doped hematite that can become
electrically conducting. Anomalous-like Hall behavior, similar to that
typically seen in ferromagnets when a magnetic field is applied along the
c-axis, has been observed in hematite [15, 12]. However, in many experimental
investigations, the resistivity of the samples has been too high, resulting in
a zero Hall coefficient when the field was oriented along the c-axis [16, 17].
Nonetheless, heavily doped samples have exhibited more conventional Hall
effect behavior, with the Hall voltage being directly proportional to the
applied field within a temperature range of 190 to 290 K [18]. Several
materials have shown new Hall effects based on magnetic properties, such as
the Spin Hall effect, inverse Spin Hall effect, or quantum spin Hall effect,
among others [19].

Recently, novel transport effects have been predicted in systems with
particular symmetry. Hematite is considered to belong to the class of
altermagnets [20, 21, 22] however, work is still required to ascertain the
altermagnetic properties. The Anomalous Hall Effect (ahe), initially studied
in ferromagnetic materials, has gained importance in the investigation of
antiferromagnets (afms) [23, 24]. It is known that the AHE can occur in
noncoplanar antiferromagnets even without the presence of spin-orbit
interaction. Altermagnets have properties similar to AFMs, so similar methods
could be exploited. In collinear antiferromagnets, the ahe can arise due to
the induced canted moment during spin reorientation, breaking the symmetry and
resulting in a non-zero ahe. Furthermore, the structural order of non-magnetic
atoms and their chirality can also contribute to the origin of the induced
moment and the ahe [22, 25, 26].

Given these intriguing observations, a key open question is the origin and
symmetry of the magnetotransport in hematite and the understanding of the
mechanisms. Specifically, elucidating the contribution from the Anomalous Hall
Effect (ahe) and its origin. So far, either the induced magnetization on the
spin reorientation or from some hidden mechanism from its altermagnetic
nature. Experimentally, so far, the origin could not be determined as the ahe
has only been probed in simple geometries without the necessary angular
dependence to allow for the required symmetry analysis. So, one could not
relate the experimental results to the predicted ahe in altermagnets is
attributed exclusively to the crystal structure [22]. The non-cero ahe
coefficient there comes from the presence of non-magnetic atoms at
noncentrosymmetric positions, despite the inversion center at the magnetic
atoms. However, measuring a Hall voltage on hematite presents challenges due
to its insulating nature. Doping the material becomes necessary to overcome
this limitation and access the transport properties. It is crucial to
carefully select magnetic and non-magnetic dopants to avoid inducing stress
and altering the properties of hematite [13, 12, 27]. Based on previous
reports [15, 28], it has been concluded that when the doping level is below
2%, the dopants added to the material will be fully ionized at all our
temperatures. Consequently, the number of carriers present in the material is
expected to remain relatively constant with temperature, changing by no more
than one order of magnitude, and the low doping level does not significantly
affect the magnetism.

In this article, we ascertain the anisotropic anomalous Hall effect in
hematite with low Ti doping. This material possesses the necessary crystal
structure and magnetic ordering symmetries that can lead to generating a
transverse voltage when particular symmetries are broken. By introducing low
doping levels, we effectively close the bandgap by putting defect states into
the energy gap, enabling access to the electronic transport properties of the
material. Specifically, we study 1% Ti-doped hematite thin films with a
thickness of 150 nm, highlighting the intriguing characteristics of the
anomalous Hall effect in this system. Our measurements reveal a strong
anisotropy in the crystal orientation dependence of the Hall effect,
demonstrating a pronounced sensitivity to the crystallographic axis. We
observe a sign inversion in the antisymmetric part of the signal, which defies
intuitive expectations. This investigation reveals novel symmetry-related
properties of hematite’s electronic and magnetic properties. In particular, we
discover an unconventional crystal dependence of the transverse voltage and an
unintuitive sign inversion in the antisymmetric part that is elucidated by
high-precision angular-dependent Hall measurements.

## Results and discussion

The discovery of altermagnets suggests a new classification of magnetic
materials that exhibit fascinating characteristics beyond those of traditional
antiferromagnets [29]. These altermagnetic materials, such as ruthenium oxide,
hematite, and various other compounds, display intriguing traits such as spin
splitting and non-zero anomalous Hall coefficients [29]. Hematite, as an
example of an altermagnetic material, potentially exhibits these unique
properties. Its hexagonal crystalline structure (Fig. 1.a) consists of two
sub-lattices of iron atoms with oxygen atoms positioned at non-centrosymmetric
sites. At different temperatures, hematite exhibits three distinct non-
paramagnetic phases, as illustrated in Fig. 1.b, which can be described using
magnetic space and point groups. By controlling the percentage of material
doping, in our case, we can shift the Fermi level by some tenths of eV [30,
31], providing access to the electronic transport properties of the material
while keeping the antiferromagnetic properties of hematite.

Figure 1:  (a) Visual representation of hematite’s hexagonal crystal
structure, showcasing the two iron sublattices (red and blue) and the oxygen
atoms (black) in non-centrosymmetric positions. (b) The different magnetic
phases of hematite, with corresponding magnetic space and point groups. The
Morin temperature defines the altermagnetic phases that have and lack weak
ferromagnetism, resulting in a non-zero magnetization at high temperatures.

High-quality 150 nm films of 1% Ti-doped hematite were fabricated (refer to
the Materials and Methods section for detailed fabrication procedures). Hall
bars were patterned at various in-plane orientations relative to the a-axis.
Magnetization measurements (Fig. 2.b) revealed a decrease in the magnetic
moment below 338 K, which was identified as the Morin transition (MT). Unlike
a single crystal sample [10], this transition exhibits a broader nature,
attributed to the existence of a multidomain state in our samples [5]. The
easy-axis phase (shown in Figure 2.a) does not exhibit a significant magnetic
moment as expected for a system with magnetic properties of hematite.
Furthermore, the magnitude of the canted moment agrees with the bulk value, as
reported in references [4] and [11]. The absence of secondary phases (see
Supporting Information), highly epitaxial thin films, well-defined in-plane
orientation, the presence of the Morin transition, and the lack of a higher
canted magnetic moment suggest that the doping level does not significantly
impact the magnetic properties of hematite.

Figure 2: Measurements of SQUID magnetometry for a magnetic field applied
perpendicular to the c-axis. (a) Magnetization per formula unit as a function
of the field for the Ti-doped samples in both the collinear phase (300 K) and
the weak ferromagnetic phase (350 K). At high temperatures, the hysteresis
loop for the canted moment saturates at 0.5 T. The signal observed at lower
temperatures is attributed to the nucleation of the $\gamma$-Fe2O3 phase at
the substrate interface. (b) The magnetization as a function of temperature
under an external magnetic field of 20 mT reveals a decrease in magnetization
as temperature increases due to the Morin transition. This transition from the
collinear phase (left panel) to the weak ferromagnetic phase (right panel)
occurs over a temperature range of 20 Kelvin. (c) The critical field ($H_{c}$)
required to reach the spin flop state varies depending on the temperature, as
well as the longitudinal resistivity decreases proportionally. $H_{c}$ is well
described by $H_{cr}(T)\propto\sqrt{T_{M}-T}$ following the Landau theory,
while the resistivity behaves as transport due to hopping in the impurity band
of the donors.

Determining the conducting regime is essential in selecting the appropriate
analysis approach. In the case of our poorly conducting samples, this is
achieved by studying the resistivity behavior as a function of temperature
(blue line in Fig. 2.c). While a semiconductor-like model fits reasonably
well, the slope yields a gap of 1600 $k_{B}$, which does not align with the
gap of hematite. Moreover, the thermal activation factor value ($5\times
10^{-3}$) suggests that the transport mechanism is likely variable-range
hopping [32].

Given the percentage of dopants in our samples (1%), we expect a relatively
large number of charge carriers (approximately $10^{15}\ \mathrm{cm}^{-3}$ in
our case). Thus, conduction via hopping in the impurity "band" of the donors
would be most probable. The relatively low doping level allows for the
measurement of some Hall conductivity in the poorly conducting regime [33],
without significantly altering the intriguing symmetry properties of hematite.
This is evident in the very low conventional Hall coefficient (approximately
$5\ \mathrm{cm}^{3}\mathrm{C}^{-1}$) observed over a wide range of
temperatures (Figure 3) [18, 28]. However, the anomalous Hall signal increases
around the spin flop transition (Fig. 3.a).

Figure 3: Hall resistivity in relation to the oop field for temperatures below
(a) and above (b) the Morin transition. The shifting of the critical field
with temperature in the collinear phase indicates the equivalence between the
spin-flop and the Morin transition as expected for hematite. On the right side
of the graph, the magnetic configurations of the Fe atoms in both phases are
displayed while maintaining crystal symmetries.

The Hall signal in the collinear phase is almost zero at low fields, followed
by a substantial anomalous behavior around the spin-flop transition. At the
spin-flop transition, the Néel vector undergoes a rotation from the out-of-
plane (OOP) direction into the basal plane, resulting in a canted magnetic
moment. The critical field is determined where the transverse resistivity
increases, and it is energy-equivalent to the Morin transition. The
temperature dependence of the critical field is well described by
$H_{cr}\propto\sqrt{T_{M}-T}$ (red dotted line in Fig 2.c), with a critical
temperature ($T_{M}$) of $337\pm 8$ K. This behavior aligns with a typical
second-order phase transition within the Landau theory, as previously reported
[34, 35]. However, the presence of a multidomain state in our sample leads to
a different definition of the Morin transition compared to this critical
temperature. In our case, the critical point is determined as the point where
the slope starts changing, either when decreasing the temperature (MT) or
increasing the field ($H_{c}$).

In the collinear phase (Fig 3.a), when an out-of-plane (OOP) field is applied,
a fully antisymmetric signal is observed only when the current direction is
parallel to one of the three a-axes. A small hysteretic behavior in the Hall
signal may be attributed to possible magnetic domain configurations [5, 13].
The size of the domains present in each structure can differ due to strain-
induced anisotropy introduced during the patterning process [36].

In the weak ferromagnetic phase (Fig 3.b), the Hall signal agrees with the
expected behavior from the Anomalous Hall Effect (AHE), as it saturates in
accordance with the canted moment of the samples.

Interesting behavior arises when the Hall bar is patterned in an orientation
that deviates from the high symmetry a-axis (Fig. 4).

Figure 4: (a) The Hall conductivity at 300 K of different devices measured at
various angles relative to the in-plane crystallographic axis, denoted as a
$(10\bar{1}0)$. A diagram illustrating the different device orientations is
provided below. Following symmetrization (as stated in Eq. 4), the even (b)
and odd (c) contributions are distinguished. The even contribution exhibits
significantly different weight in the Hall signal across the various relative
orientations, whereas the odd portion’s sign inversion cannot be attributed to
weak moment-induced anomalous Hall effect contribution.

When the current channel is along a lower symmetry axis, it allows access to
the symmetric components of the conductivity tensor (Fig. 4.b). This effect is
only observed when the conductivity components are anisotropic and the
orientation of the Hall bar is appropriate. This effect is zero in highly
symmetric orientations (black and gray lines), while in less symmetrical
orientations (purple and blue lines), the even part of the conductivity has a
larger contribution than the odd part.

The odd part of the Hall conductivity for Hall bars with different
orientations relative to the a-axis exhibits a sign inversion at 30, 45, and
90 degrees but preserves the sign at rotations of 120 degrees (Figure 4.c).
This behavior is unexpected within the framework of conventional hopping
transport and suggests that it may be attributed to spin-splitting from the
altermagnetic nature. However, we note that in our poorly conducting transport
regime, we cannot quantify altermagnetic spin-splitting effects.

In terms of symmetry, the shift from the a-axis oriented Hall bar to the $\pm
45^{\circ}$ orientation results in a change from a predominantly antisymmetric
signal to a symmetric behavior, accompanied by a sign inversion. The odd
conductivity components display a decrease in amplitude when the Hall bar is
oriented along the a-axis compared to other orientations. However, the crystal
equivalent orientations of the Hall bars exhibit similar symmetry behavior
(-30 and 90, 45 and -45, 0 and 120 degrees). An inversion of the sign is
observed when the current is at 45 and -45 degrees relative to the a-axis, as
well as from 0 and -30, or 0 and 90 degrees, but not between Hall bars
oriented 120 degrees relative to each other.

In the poorly conducting regime of our samples, the dominant transport
mechanism is variable-range hopping, accompanied by a low Hall mobility of
approximately $0.1\ \mathrm{cm}^{2}\mathrm{V}^{-1}\mathrm{s}^{-1}$. Despite
the limited conductivity, we are able to access both the symmetric and
antisymmetric components of the conductivity tensor. The symmetries observed
in the Hall effect are directly correlated with the Hall bars’ relative
orientation and the material’s crystallographic structure. Significantly, this
dependence on orientation is not influenced by the canted moments induced
during spin reorientation, as the induced canted moment remains consistent for
all device orientations at a fixed field strength. It can be concluded that
the conductivity tensor is non-symmetric, and these non-symmetric terms can
only be accessed when the symmetry of the system is lowered, i.e., the current
vector is not along a high symmetry axis.

In conclusion, the inversion of the sign in the odd part of the Hall
conductivity in hematite may originate from the altermagnetic magnetic
properties present in the inaccessible but still representative bands of the
material. This sign inversion signifies the presence of unique underlying
mechanisms governing the anomalous Hall effect (AHE) in this particular
material.

These findings highlight the intricate interplay between crystal structure,
magnetic properties, and conductivity in our samples. Investigating the
anisotropic anomalous Hall effect in hematite provides significant insights
into this material’s electronic and magnetic properties. By unraveling the
complex relationship between symmetry breaking and the resulting Hall signals,
we advance our understanding of this phenomenon and its potential applications
in various fields. The intricate symmetry of the Hall signal observed is
intriguing and will spark further investigations, particularly for different
doping levels that go beyond the current work, to fully understand and
elucidate the specific mechanisms responsible for the observed symmetries and
sign inversion in the Hall conductivity. Our findings contribute to a deeper
knowledge of hematite’s behavior, laying the groundwork for its utilization in
areas such as electronics, energy storage, and data storage.

Furthermore, this study serves as a stepping stone for exploring similar
materials with analogous properties and provides valuable insights that
contribute to our broader understanding of electronic and magnetic phenomena.
Our research establishes a foundation for future investigations into the
anomalous Hall effect in altermagnets and its implications for emerging
technologies.

## Acknowledgments

This work was supported by the Max Planck Graduate Center with the Johannes
Gutenberg-Universität Mainz (MPGC). The authors in Mainz acknowledge support
from the DFG project number 423441604. All authors from Mainz also acknowledge
support from SPIN + X (DFG SFB TRR 173 No. 268565370, projects A01, A03, and
B02). M.K. acknowledges support by the Research Council of Norway through its
Centers of Excellence funding scheme, project number 262633 "QuSpin." C.-Y.Y.
and E.B. acknowledge support from the National Research Foundation of Korea
(Grant No. 2021M3F3A2A01037525).

## Supporting Information

## Materials and Methods

Thin (0001)-oriented hematite films of 150 nm thickness doped with 1% Ti
[$\alpha$-Fe1.99Ti0.01O3:(Ti)$\alpha$-Fe2O3] were deposited on c-plane (0001)
sapphire substrates using pulsed laser deposition (pld). The deposition was
carried out at a set point temperature of 800 ∘C under an oxygen partial
pressure of 10 mTorr. The films were grown from a stoichiometric doped target,
and a total of 20,000 pulses were used to achieve the desired thickness of
approximately 150 nm, as described in [37]. The crystallographic purity of the
thin films was confirmed using X-ray diffraction (xrd). Also, the in-plane
crystallographic orientation of the sample was determined (See Crystalline and
epitaxial-quality determination), and Hall bars were patterned in different
orientations with respect to the crystallographic direction a.

The temperature range for the colinear phase of the hematite films was
determined using superconducting quantum interference device (squid)
magnetometry measurements (See Magnetometry measurements). All structures were
fabricated using two electron beam lithography steps. First, a protective
resist mask was defined to allow ion-beam etching until the substrate. Then,
contacts were implemented using a lift-off technique with a Cr/Au bilayer. Six
different orientations were patterned in the samples, with angles of
$\phi=0,30,45,60,90,120^{\circ}$ in the plane, where $\phi=0$ points along the
a axis.

The longitudinal resistance of the Hall bar was measured using a high-
impedance electrometer while applying a current of $I=2\times 10^{6}\
\mathrm{A}/\mathrm{cm}^{2}$. A highly sensitive nanovoltmeter measured the
transverse voltage ($V_{xy}$). Various voltage components, including the
voltmeter offset ($V_{0}$), thermoelectric voltages ($V_{S}$), Nernst effect
voltage, and Righi-Leduc voltage, were subtracted using a current
symmetrization technique with the current $I$ of the Hall bar. The Hall
resistivity ($\rho_{xy}^{Hall}$) as a function of the external magnetic field
$H$ was calculated by multiplying by the geometrical factor of the Hall bar
[38]:

$\rho_{xy}^{Hall}=R_{xy}^{Hall}\cdot t$ (1)

The longitudinal resistivity is given by:

$\rho_{xx}=\frac{R_{xx}\cdot w\cdot t}{a}$ (2)

In our case, $t=150$ nm represents the thickness of the sample, $w=10\ \mu$m
is the width of the Hall bar, and $a=80\ \mu$m is the distance between the two
channels of the 6-point Hall bar.

The Hall conductivity can be calculated using the values of $\rho_{xx}$ and
$\rho_{xy}$:

$\sigma_{xy}=\frac{\rho_{xy}}{\rho^{2}_{xx}+\rho^{2}_{xy}}$ (3)

To analyze the results under symmetry arguments, the transversal conductivity
is symmetrized/anti-symmetrized with respect to the external magnetic field
using the equation:

$\sigma^{A(S)}{xy}=\sigma{xy}(H)-(+)\sigma_{xy}(-H)$ (4)

The electrical measurements were performed in a ${}^{4}\mathrm{He}$ cryostat,
capable of temperatures between $1.6$ K $-360$ K and magnetic fields up to 12
Teslas. Thermal stabilization was ensured during all measurements, monitored
by a Cernox sensor element with a variation of less than 0.03% throughout the
measurements. The tendency lines were calculated based on b-splines, and the
dispersion across five different devices was used to calculate the errors.

### Crystalline and epitaxial-quality determination

X-ray measurements were carried out on the thin film samples before and after
the etching process using a four-circle Bruker D8 Discover thin film
diffractometer. $\theta/2\theta$ measurements (Figure S1.a) were used to
evaluate the highly oriented Ti-doped hematite films without any impurity
phase present. Phi-scans (Figure S1.c) on the ($10\bar{1}4$) peak allowed us
to identify the orientations of the [$10\bar{1}0$] planes. Figure S1.a shows
the xrd pattern, where a diffraction peak at $2\theta=39^{\circ}$ is observed,
slightly smaller than the bulk value [39]. This confirms previous reports on
homogeneous doping of hematite grown by pld [40]. The absence of impurity
phases and the (0001)-orientation of the films were further confirmed by
2$\theta=\omega$ measurements (inset in Figure S1.a). The lattice parameters
were calculated from the peak positions and found to reproduce the literature
values of values of a = 5.05(1) Å and c = 13.72(2) Å[41]. Additionally, a
$\phi$ scan (Figure S1.c) allows us to determine the in-plane orientation of
the samples, which was used to pattern Hall bars in different orientations
with respect to the crystallographic direction a.

Figure S1: (a) The X-ray diffraction analysis of Ti-doped $\alpha$-Fe2O3 with
(0001) orientation did not reveal any secondary phases. An inset of the film’s
(0006) peak showed that the films were epitaxial in nature. (b) The diagram
illustrates the crystallographic plane’s 3-fold symmetry, where $\alpha_{2}$
and $\alpha_{4}$ angles were equivalent with a field inversion, as were
$\alpha_{1}$ and $\alpha_{3}$. (c) The in-plane xrd scan further confirmed the
crystal’s symmetry, with a 120-degree separation and similar heights
indicating a well-defined in-plane symmetry.

### Magnetometry measurements

The squid measurements were performed using a commercial superconducting
quantum interference device (squid) from Quantum Design. The films were
saturated with 5 T at the starting temperature of 355 K, and then were cooled
at a constant rate of 2 K/min under a static applied field of $\mu_{0}H=20$
mT. Due to the diamagnetic contribution to the signal from the sapphire
substrates used to grow the films and the mounting of the sample in the SQUID,
a temperature-independent subtraction to the signal was made from the data
from the average value of the magnetization below the observed transition for
the temperature scans, and as a linear subtraction calculated at higher fields
($\pm$3 T to $\pm$ 5 T) for the field scans. The Morin transition of the films
was observed, characterized by a drop in magnetic moment to effectively zero
above the background substrate contributions. The transition occurred at
around 330 K, and at 320 K, the transition was complete, indicating that the
films were in the purely easy-axis phase. The squid measurements agreed with
other reports, where dopants were found to increase the Morin transition
temperature even above room temperature [13]. The presence of Ti atoms did not
significantly affect the magnetic moments in the samples, as the magnitude of
the canted moment in the wfm phase fell within the range reported for pure
single crystals [39, 4]. Both the temperature increment and the broadening of
the Morin transition can be attributed to the smaller unit cell of the film as
compared to the bulk values, leading to a change in the temperature dependence
of the competing anisotropies that underlie the physics of the Morin
transition[42]. We also note a soft loop at low magnetic fields in the easy-
axis phase, possibly due to a thin ferrimagnetic seed layer that nucleates at
the hematite/substrate interface [13, 43]

### Transport properties on the weak ferromagnetic phase

The Hall conductivity also shows a dependence on the relative crystal
orientation in the weak ferromagnetic phase, also showing a change of sign
between different orientations of the Hall bar. Those aspects do not appear at
lower fields on the EA phase and are not only contributions from the canted
moment.

Figure S2: (a) Hall conductivity of devices along the easy-axis phase (a-axis)
was measured and compared with the conductivity in the WFM phase and at a
45-degree relative orientation. The even (b) and odd (c) parts of the signal
reveal a sign inversion in the odd contribution with the two relative
orientations.

## References

  * [1] V. Baltz et al. “Antiferromagnetic spintronics” In _Rev. Mod. Phys._ 90 American Physical Society, 2018, pp. 015005 DOI: 10.1103/RevModPhys.90.015005
  * [2] W. R. W. Ahmad et al. “A review on hematite $\alpha$-Fe2O3 focusing on nanostructures, synthesis methods and applications” In _2016 IEEE Student Conference on Research and Development (SCOReD)_ , 2016, pp. 1–6 DOI: 10.1109/SCORED.2016.7810090
  * [3] A. Mirzaei, B. Hashemi and K. Janghorban “$\alpha$-Fe2O3 based nanomaterials as gas sensors” In _J. Mater. Sci. Mater. Electron._ 27.4, 2015, pp. 3109–3144 DOI: 10.1007/s10854-015-4200-z
  * [4] R. Lebrun et al. “Tunable long-distance spin transport in a crystalline antiferromagnetic iron oxide” In _Nature_ 561.7722, 2018, pp. 222–225 DOI: 10.1038/s41586-018-0490-7
  * [5] A. Ross et al. “Propagation Length of Antiferromagnetic Magnons Governed by Domain Configurations” PMID: 31809058 In _Nano Lett._ 20.1, 2020, pp. 306–313 DOI: 10.1021/acs.nanolett.9b03837
  * [6] A. Ross et al. “Structural sensitivity of the spin Hall magnetoresistance in antiferromagnetic thin films” In _Phys. Rev. B_ 102 American Physical Society, 2020, pp. 094415 DOI: 10.1103/PhysRevB.102.094415
  * [7] H. Meer, O. Gomonay, A. Wittmann and M. Kläui “Antiferromagnetic insulatronics: Spintronics in insulating 3d metal oxides with antiferromagnetic coupling” 080502 In _Appl. Phys. Lett._ 122.8, 2023 DOI: 10.1063/5.0135079
  * [8] Y. Cheng et al. “Electrical Switching of Tristate Antiferromagnetic Néel Order in $\alpha$-Fe2O3 Epitaxial Films” In _Phys. Rev. Lett._ 124 American Physical Society, 2020, pp. 027202 DOI: 10.1103/PhysRevLett.124.027202
  * [9] E. Cogulu et al. “Direct imaging of electrical switching of antiferromagnetic Néel order in $\alpha$-Fe2O3 epitaxial films” In _Phys. Rev. B_ 103 American Physical Society, 2021, pp. L100405 DOI: 10.1103/PhysRevB.103.L100405
  * [10] F. J. Morin “Magnetic Susceptibility of $\alpha$-Fe2O3 and $\alpha$-Fe2O3 with Added Titanium” In _Phys. Rev._ 78 American Physical Society, 1950, pp. 819–820 DOI: 10.1103/PhysRev.78.819.2
  * [11] T. Dannegger et al. “Magnetic properties of hematite revealed by an ab initio parameterized spin model” In _Phys. Rev. B_ 107 American Physical Society, 2023, pp. 184426 DOI: 10.1103/PhysRevB.107.184426
  * [12] G. A. Acket and J. Volger “Electric transport in N-type Fe2O3” In _Physica_ 32.9, 1966, pp. 1543–1550 DOI: https://doi.org/10.1016/0031-8914(66)90034-6
  * [13] A. Ross et al. “An insulating doped antiferromagnet with low magnetic symmetry as a room temperature spin conduit” In _Appl. Phys. Lett._ 117.24, 2020, pp. 242405 DOI: 10.1063/5.0032940
  * [14] J. Han et al. “Birefringence-like spin transport via linearly polarized antiferromagnetic magnons” In _Nat. Nanotechnol._ 15.7, 2020, pp. 563–568 DOI: 10.1038/s41565-020-0703-8
  * [15] A. J. E. Rettie et al. “Synthesis, electronic transport and optical properties of Si:$\alpha$-Fe2O3 single crystals” In _J. Mater. Chem. C_ 4 The Royal Society of Chemistry, 2016, pp. 559–567 DOI: 10.1039/C5TC03368C
  * [16] K.B. Vlasov, E. A. Rosenberg, A. G. Titova and Yu. M. Yakovlev “Hall effect in hematite single crystals in the region of the Morin transition” In _Sov. Phys.-Solid State_ 22 American Institute of Physics, 1980, pp. 967 —969
  * [17] Sergey V. Ovsyannikov, Natalia V. Morozova, Alexander E. Karkin and Vladimir V. Shchennikov “High-pressure cycling of hematite Fe2O3: Nanostructuring, in situ electronic transport, and possible charge disproportionation” In _Phys. Rev. B_ 86 American Physical Society, 2012, pp. 205131 DOI: 10.1103/PhysRevB.86.205131
  * [18] B. Zhao et al. “Electrical transport properties of Ti-doped Fe2O3(0001) epitaxial films” In _Phys. Rev. B_ 84 American Physical Society, 2011, pp. 245325 DOI: 10.1103/PhysRevB.84.245325
  * [19] J. Sinova et al. “Spin Hall effects” In _Rev. Mod. Phys._ 87 American Physical Society, 2015, pp. 1213–1260 DOI: 10.1103/RevModPhys.87.1213
  * [20] L. Wang et al. “Crystal-induced transverse current in collinear antiferromagnetic $\gamma$-FeMn” In _Appl. Phys. Lett._ 120.1, 2022, pp. 012403 DOI: 10.1063/5.0069504
  * [21] J. Kipp et al. “The chiral Hall effect in canted ferromagnets and antiferromagnets” In _Commun. Phys._ 4.1, 2021 DOI: 10.1038/s42005-021-00587-3
  * [22] L. Šmejkal, R. González-Hernández, T. Jungwirth and J. Sinova “Crystal time-reversal symmetry breaking and spontaneous Hall effect in collinear antiferromagnets” In _Sci. Adv._ 6.23, 2020, pp. eaaz8809 DOI: 10.1126/sciadv.aaz8809
  * [23] R. Shindou and N. Nagaosa “Orbital Ferromagnetism and Anomalous Hall Effect in Antiferromagnets on the Distorted fcc Lattice” In _Phys. Rev. Lett._ 87 American Physical Society, 2001, pp. 116801 DOI: 10.1103/PhysRevLett.87.116801
  * [24] L. Šmejkal, Y. Mokrousov, B. Yan and A. H. MacDonald “Topological antiferromagnetic spintronics” In _Nat. Phys._ 14.3, 2018, pp. 242–251 DOI: 10.1038/s41567-018-0064-5
  * [25] K. Samanta et al. “Crystal Hall and crystal magneto-optical effect in thin films of SrRuO3” In _J. Appl. Phys._ 127.21, 2020, pp. 213904 DOI: 10.1063/5.0005017
  * [26] D.-F. Shao et al. “Interfacial Crystal Hall Effect Reversible by Ferroelectric Polarization” In _Phys. Rev. Appl._ 15 American Physical Society, 2021, pp. 024057 DOI: 10.1103/PhysRevApplied.15.024057
  * [27] R.N. Bhowmik and A. G. Lone “Electric field controlled magnetic exchange bias and magnetic state switching at room temperature in Ga-doped $\alpha$-Fe2O3 oxide” In _J. Magn. Magn. Mater._ 462, 2018, pp. 105–118 DOI: https://doi.org/10.1016/j.jmmm.2018.05.007
  * [28] E. Gharibi, A. Hbika, B. Dupre and C. Gleitzer “ChemInform Abstract: Electrical Properties of Pure and Titanium-Doped Hematite Single Crystals, in the Basal Plane, at Low Oxygen Pressure.” In _Chem. Inform._ 21.43, 1990 DOI: https://doi.org/10.1002/chin.199043022
  * [29] L. Šmejkal, J. Sinova and T. Jungwirth “Beyond Conventional Ferromagnetism and Antiferromagnetism: A Phase with Nonrelativistic Spin and Crystal Rotation Symmetry” In _Phys. Rev. X_ 12 American Physical Society, 2022, pp. 031042 DOI: 10.1103/PhysRevX.12.031042
  * [30] T. R. Paudel et al. “Doping Rules and Doping Prototypes in A2BO4 Spinel Oxides” In _Advanced Functional Materials_ 21.23, 2011, pp. 4493–4501 DOI: https://doi.org/10.1002/adfm.201101469
  * [31] J. Robertson and S. J. Clark “Limits to doping in oxides” In _Phys. Rev. B_ 83 American Physical Society, 2011, pp. 075205 DOI: 10.1103/PhysRevB.83.075205
  * [32] Y. L. Zhao et al. “Variable range hopping in TiO2 insulating layers for oxide electronic devices” In _AIP Advances_ 2.1, 2012, pp. 012129 DOI: 10.1063/1.3682346
  * [33] S. Onoda, N. Sugimoto and N. Nagaosa “Quantum transport theory of anomalous electric, thermoelectric, and thermal Hall effects in ferromagnets” In _Phys. Rev. B_ 77 American Physical Society, 2008, pp. 165103 DOI: 10.1103/PhysRevB.77.165103
  * [34] P. J. Besser, A. H. Morrish and C. W. Searle “Magnetocrystalline Anisotropy of Pure and Doped Hematite” In _Phys. Rev._ 153 American Physical Society, 1967, pp. 632–640 DOI: 10.1103/PhysRev.153.632
  * [35] S. Foner “High-Field Antiferromagnetic Resonance in Cr2O3” In _Phys. Rev._ 130 American Physical Society, 1963, pp. 183–197 DOI: 10.1103/PhysRev.130.183
  * [36] H. Meer et al. “Strain-induced shape anisotropy in antiferromagnetic structures” In _Phys. Rev. B_ 106 American Physical Society, 2022, pp. 094430 DOI: 10.1103/PhysRevB.106.094430
  * [37] D. A. Grave et al. “Heteroepitaxial hematite photoanodes as a model system for solar water splitting” In _J. Mater. Chem. A_ 4 The Royal Society of Chemistry, 2016, pp. 3052–3060 DOI: 10.1039/C5TA07094E
  * [38] D. S. Perloff “Four-point sheet resistance correction factors for thin rectangular samples” In _Solid State Electron._ 20, 1977, pp. 681–687
  * [39] A.H. Morrish “Canted Antiferromagnetism: Hematite” World Scientific, 1994 URL: https://books.google.de/books?id=E2X8uW8UxmEC
  * [40] K. D. Malviya et al. “Systematic comparison of different dopants in thin film hematite ($\alpha$-Fe2O3) photoanodes for solar water splitting” In _J. Mater. Chem. A_ 4 The Royal Society of Chemistry, 2016, pp. 3091–3099 DOI: 10.1039/C5TA07095C
  * [41] O. M. Lemine et al. “Rietveld analysis and Mössbauer spectroscopy studies of nanocrystalline hematite $\alpha$-Fe2O3” In _J. Alloys Compd._ 502.2, 2010, pp. 279–282 DOI: https://doi.org/10.1016/j.jallcom.2010.04.175
  * [42] J. O. Artman, J. C. Murphy and S. Foner “Magnetic Anisotropy in Antiferromagnetic Corundum-Type Sesquioxides” In _Phys. Rev._ 138 American Physical Society, 1965, pp. A912–A917 DOI: 10.1103/PhysRev.138.A912
  * [43] I. Dzyaloshinsky “A thermodynamic theory of “weak” ferromagnetism of antiferromagnetics” In _J. Phys. Chem. Solids_ 4.4, 1958, pp. 241–255 DOI: https://doi.org/10.1016/0022-3697(58)90076-3


[a,b]S. Kováčik

# The Fuzzy Onion: An Initial Study

J. Tekel  M. Hrmo

###### Abstract

In our previous contribution, we introduced a matrix formulation of a three-
dimensional quantum space named the fuzzy onion. The novel part of the
construction is the radial derivative term, which has been defined to recover
the correct continuum limit. Here, we describe a numerical simulation of the
scalar field theory in this space and test some physical properties of the
model with emphasis on the interaction between neighbouring layers.

##  1 Introduction

A long-term goal of physics is to understand the phenomenon of quantum space
to see if consistent physical theories can be built on it and then to draw
verifiable predictions from that theory. One of the most studied examples is
the fuzzy sphere, which has been studied in various contexts [2, 1, 3, 4, 5,
6, 11, 12, 10, 7, 8, 9, 13, 14, 15, 16].

The upside of the fuzzy sphere model is its simplicity and transparency, which
is the reason behind the numerous studies performed on it. The downside is
obvious ––– the corresponding geometry differs from the space we live in. Had
we been living in a two-dimensional space, the fuzzy sphere would be a good
candidate for the model of quantum space. Instead, it is either an interesting
mathematical model to study basic properties of quantum space or a useful
regularisation tool, [17, 18]. In our previous contribution, we proposed a way
of connecting fuzzy spherical layers to form the fuzzy onion space. This model
has two advantages: it is still nearly as simple as the original fuzzy sphere
model, and it describes a three-dimensional space.

The construction of the fuzzy onion as presented there, had two nontrivial
steps. The first was connecting various layers. If the fundamental length
scale is fixed along increasing sphere radius, one has to use matrices of
increasing size. Comparing them, for example, to calculate the radial
derivative is naively impossible. This has been overcome by using the Fourier
transform and discarding or adding the unmappable degrees of freedom. The
second step was to define the full Laplace operator that captures the geometry
of the space. There is some ambiguity in this definition; ours differs from
that of [19, 20]. Our goal was to reproduce the results of the model of three-
dimensional quantum space studied in [3, 4] that has been studied thoroughly
and has well-established commutative limit; some preliminary tests of the
construction have been passed, see [23].

##  2 The fuzzy sphere and the fuzzy onion

We begin with a brief reminder of the fuzzy sphere and the fuzzy onion
constructions. Fields on a sphere can be expanded into spherical harmonics,
which form an infinite-dimensional representation of the $SO(3)$ generators
$\hat{L}_{i}$:

$\displaystyle f(\theta,\varphi)$ $\displaystyle=$
$\displaystyle\sum_{lm}c_{lm}Y_{lm}(\theta,\varphi),$
$\displaystyle\hat{{\cal{L}}}\ Y_{lm}$ $\displaystyle=$
$\displaystyle\hat{L}_{i}\hat{L}_{i}\ Y_{lm}=l(l+1)Y_{lm},$ (1)
$\displaystyle\hat{L}_{3}\ Y_{lm}$ $\displaystyle=$ $\displaystyle m\ Y_{lm},$

Finite-dimensional representations with momentum cut-off exist; that is, they
have a maximum value of $l$. In this case, the fields can be represented as
$N\times N$ matrices $Y_{lm}$, where $l=N-1$ is the maximal angular momentum.
Fields are encoded into $N\times N$ matrices, and the angular momentum
eigenstates now form the basis for Hermitian matrices of this size.

The matrix algebra generated by $Y_{lm}$ is closed under multiplication and
noncommutative. Both fields on the sphere and Hermitian matrices can form an
$so(3)$ representation space, but they differ in the momentum cut-off. The
expansion coefficients serve as a map between functions on the sphere and
Hermitian matrices. If one imposes a momentum cut-off on the sphere, this map
can be bijective. However, this comes with a loss — with a limited number of
expansion terms, one can only approximate the Dirac $\delta$-function, or in
other words, the spatial resolution is limited. In this way, Hermitian
matrices can be used to encode fields on a sphere with limited spatial
resolution. The minimal distinguishable scale is given by the parameter
$\lambda$, where
$[\hat{x}_{i},\hat{x}_{j}]=2i\lambda\varepsilon_{ijk}\hat{x}_{k}$ and
$\hat{x}_{i}$ are position operators.

The radius of the fuzzy sphere is $r\sim N\lambda$, which means that if the
minimal length scale is fixed, to describe a larger sphere, one needs a larger
matrix. Layers of fuzzy spheres with increasing radius, denoted by the
superscript (i), can then be put together into a larger matrix of the form

$\Psi=\begin{pmatrix}\Phi^{(1)}&&&\\\ &\Phi^{(2)}&&\\\ &&\ddots&\\\
&&&\Phi^{(N_{m})}\end{pmatrix},$ (2)

where the angular part of the Laplace operator acts on each of the layers
separately

${\cal{L}}\Psi=\begin{pmatrix}\hat{{\cal{L}}}^{(1)}\Phi^{(1)}&&&\\\
&\hat{{\cal{L}}}^{(2)}\Phi^{(2)}&&\\\ &&\ddots&\\\
&&&\hat{{\cal{L}}}^{(N_{m})}\Phi^{(N_{m})}.\end{pmatrix}.$ (3)

To define the radial part of the Laplace operator, first, we need to define an
operator that connects submatrices $\Phi^{(i)}$ on consecutive layers:

$\displaystyle{\cal{U}}\Phi^{(i)}$ $\displaystyle=$
$\displaystyle\Phi^{(i+1)}\in{\cal{H}}(N+1),$ (4)
$\displaystyle{\cal{D}}\Phi^{(i)}$ $\displaystyle=$
$\displaystyle\Phi^{(i-1)}\in{\cal{H}}(N-1),$ (5)

where ${\cal{H}}(N)$ is the set of Hermitian matrices of size $N$. This can be
done by first doing the Fourier expansion and then discarding the unmappable
coefficients when going one layer down or adding zeroes when going one layer
up:

$\displaystyle\Phi^{(N)}$ $\displaystyle=$
$\displaystyle\sum\limits_{l=1}^{N-1}\sum\limits_{m=-l}^{l}c^{(N)}_{lm}Y_{lm}^{(N)},$
$\displaystyle{\cal{D}}\uparrow\hskip 14.22636pt{\cal{U}}\downarrow$
$\displaystyle\Phi^{(N+1)}$ $\displaystyle=$
$\displaystyle\sum\limits_{l=1}^{N-1}\sum\limits_{m=-l}^{l}c^{(N+1)}_{lm}Y_{lm}^{(N+1)},$
where $\displaystyle c_{l,m}^{(N)}$ $\displaystyle=$ $\displaystyle
c_{l,m}^{(N+1)}\text{ for: }l\leq N-1,$ $\displaystyle c_{N,m}^{(N+1)}$
$\displaystyle=$ $\displaystyle 0.$

Using this, we are in a position to define the radial derivative and the
radial part of the Laplace operator. This is the source of some ambiguity as
infinitely many definitions of the derivative coincide in the
$\lambda\rightarrow 0$ limit. We chose a simple option that reproduces the
ordinary energy spectrum, but a more detailed study is needed here:

$\partial_{r}\Phi^{(N)}=\frac{{\cal{D}}\phi^{(N+1)}-{{\cal{U}}}\phi^{(N-1)}}{2\lambda},$
(6)

and

$\partial^{2}_{r}\Phi^{(N)}=\frac{{\cal{D}}\phi^{(N+1)}-2\phi^{(N)}+{{\cal{U}}}\phi^{(N-1)}}{\lambda^{2}}.$
(7)

Using those, we can define the full Laplace operator that acts on the $\Psi$
that describe fields that exist on a series of concentric fuzzy spheres of
increasing radius, that is, on the fuzzy onion space:

$\Delta\Psi=\left(\frac{1}{r^{2}}\frac{\partial}{\partial_{r}}\left(r^{2}\frac{\partial}{\partial
r}\right)-\frac{{\cal{L}}^{2}}{r^{2}}\right)\Psi$ (8)

It is trivial to define a potential on the fuzzy sphere using polynomials

$V(\phi)=\mbox{Tr}\ P(\Phi),$ (9)

where the trace operation replaces the integration — the same can be done on
the fuzzy onion. Still, we must include $4\pi r$ as an integration measure as
it differs across various fuzzy spheres. With this definition, there is
something disconcerting. While the fields are smeared on each layer due to
nonlocality, they are not smeared across various layers. This can be changed
artificially by a smearing procedure

${\cal{S}}\phi^{(n)}=\frac{\phi^{(n)}+\sum\limits_{i}\alpha_{i}\left({\cal{U}}^{i}\phi^{(n-i)}+{\cal{D}}^{i}\phi^{(n+i)}\right)}{1+\sum\limits_{i}\alpha_{i}},$
(10)

for example with $\alpha_{1}=\frac{1}{2},\alpha_{2+}=0$:

${\cal{S}}\phi^{(n)}=\frac{\phi^{(n)}+\frac{1}{2}{\cal{D}}\phi^{(n+1)}+\frac{1}{2}{\cal{U}}\phi^{(n-1)}}{2},$
(11)

and with it, we can take a potential that is nonlocal in each direction as

$V(\Psi)=\sum\limits_{j=1}^{N_{m}}V\left({\cal{S}}\phi^{(j)}\right).$ (12)

The other option is to consider the model for strings stretching between
various layers as in [21].

##  3 Scalar field theory

The scalar field theory is defined by the action

$S[\Psi]=4\pi\lambda^{2}\mbox{Tr }r\left(a\ \Psi{\cal{K}}\Psi+b\ \Psi^{2}+c\
\Psi^{4}\right),\ \ {\cal{K}}=-\frac{1}{2}\Delta,$ (13)

where $r$ is a diagonal matrix that multiplies each layer with the
corresponding radius. The mean values of observables are defined as:

$\displaystyle\left\langle{\cal{O}}(\Psi)\right\rangle=\frac{1}{Z}\int d\Psi
e^{-S(\Psi)}{\cal{O}}(\Psi),\ d\Psi=\prod_{N=1}^{M}d\Phi^{(N)}\ .$ (14)

Basically, this is a Hermitian matrix model in which we consider only a subset
of matrices of a certain type and have a specific action that serves to define
the integration measure.

Hermitian matrices are partially described by their eigenvalues, and for a
pure potential matrix model, the action depends only on them. The
diagonalisation procedure and integrating the rest of the degrees of freedom
produces a new term in the integration measure that can be interpreted as a
logarithmic repulsion between the eigenvalues [22]. In the case of the field
theories on fuzzy spaces, the action involves a kinetic term that prevents
diagonalisation, but analysing the behaviour of the eigenvalues still proves
to be useful for understanding the model.

Usually, one takes $a=1$ and then produces the phase diagram of the model with
respect to $b,c$. While $c$ is positive, $b$ can be positive, zero or
negative. For sufficiently large values of $b$, all eigenvalues reside close
to the origin; this is called the disordered phase. There are two options for
sufficiently negative values of $b$. There are two potential minima and either
all eigenvalues are in one of them, this is called the uniformly ordered
phase, or are split between them due to repulsion; this is the nonuniformly
ordered phase.

As was shown in great detail in [13], the phase diagram for the scalar field
theory on the fuzzy sphere converges when expressed in rescaled parameters

$\tilde{b}=\frac{b}{aN^{3/2}},\ \tilde{c}=\frac{c}{aN^{2}}.$ (15)

This is interesting because, in our case, the value of $N$ varies across
layers. That means that for fixed values of $b$ and $c$, different layers
might prefer different phases. The assumption in this regard is that the
outermost layer would set the behaviour as it has the largest number of
degrees of freedom.

To understand the behaviour of the model, we have performed Hamiltonian Monte
Carlo (HMC) studies of the model. While the fuzzy sphere model has been
studied with $N>100$, the fuzzy onion model is more computationally demanding
as one has to perform the Fourier transformation at each step. Since this is
an initial study, we chose $N=10$ for the largest matrix.

There were two questions we wanted to answer. Does the radial part of the
Laplace operator work as intended and cause interaction between consecutive
layers? Does the outermost layer dictate the behaviour when there is a
mismatch between preferred phases across various layers?

To answer the first question, we have set a simulation in a regime where we
expected the field theory to be in the nonuniformly ordered phase and saved a
single configuration of the fields on the fuzzy onion after a sufficient
number of steps after the thermalisation had been achieved. In the first
simulation, the radial part of the kinetic term has been removed, and we can
observe a range of wild patterns on each layer, completely independent from
each other, see Figure 1. This is an obvious outcome as, without the radial
term, the simulation is basically a simulation of distinct fuzzy spheres. This
serves as a benchmark of what noninteraction between the layers looks like.

Then, we repeated the same simulation but with the radial term included. The
alignment is now obvious; see Figure 1. In principle, we could have calculated
correlation coefficients to be more precise, but in this case, we believe the
plots are sufficient to prove the point.

(a)

(b)

Figure 1: Plots of $\Phi(\theta,\varphi)$ in a nonuniformly ordered phase with
$a=1,b=-3,c=5$ and $N=10$ and the radial term missing (a) or present (b). We
can see that including the radial term caused alignment across consecutive
layers. Figure taken from [23].

To answer the second question, we ran a simulation with the largest matrix of
size $N=8$, the radial term included and values $b=-3.07,c=12.8$. The
eigenvalues were measured after each $10^{3}$ steps taken. The simulation was
initiated in a uniformly ordered phase; that is, for each of the layers, all
of the eigenvalues were in the same potential well. The simulations have been
pushed strongly, meaning the acceptance rate was close to $50\%$, and the
phase seemed stable. However, it turned out to be actually meta-stable only as
around the step $220\times 10^{3}$, first one and then shortly after, half of
the eigenvalues on the outermost layer moved to the other potential well, and
the field moved to the nonuniformly ordered phase. This effect then cascaded
down and shifted all but two innermost layers to the same phase. We expect
that in a properly ergodic simulation, actually, even the second innermost
layers would shift and then the innermost layer would move between two layers.

Figure 2: HMC trajectories for eigenvalues of $\Phi^{(N)}$ on each of the
layers. We can see that while all of them started in the same phase, the
change first happened at the outermost layer $\Phi^{(8)}$ and moved down to
the $\Phi^{(3)}$ layer.

##  4 Conclusion

We have performed an initial study of the model of the scalar field theory on
the fuzzy onion space. So far, the model behaves as expected. As we have shown
in [23], the model reproduces the results of the hydrogen atom problem in
noncommutative space and can be used to study some classical problems, such as
the heat transfer.

Here, we have discussed in more detail the behaviour of the scalar field
theory. Again, the model seems to be working as intended. The radial term
allows the interaction of neighbouring layers, and it causes not only the
alignment of the fields but also a rather abrupt alignment of the phases.

Obviously, more detailed studies are warranted. In the first regard, we are
interested in the large $N$ limit of the model. Also, it would be interesting
to see if fine-tuning the coefficients $(b,c)$ would allow for the coexistence
of multiple phases. Another line of research is to verify various ways of
connecting the potential terms across various layers, either in the smearing
or string formalism mentioned earlier in this report.

## Acknowledgments

This research was supported by the VEGA 1/0703/20 grant and the MUNI Award for
Science and Humanities funded by the Grant Agency of Masaryk University.

## References

  * [1] J. Madore, Class. Quant. Grav. 9 (1992), 69-88
  * [2] J. Hoppe, Ph.D. thesis, MIT, 1982
  * [3] S. Kováčik and P. Prešnajder, J. Math. Phys. 54 (2013), 102103
  * [4] V. Gáliková, S. Kováčik and P. Prešnajder, J. Math. Phys. 54 (2013), 122106
  * [5] P. Vitale and J. C. Wallet, JHEP 04 (2013), 115
  * [6] A. B. Hammou, M. Lagraa and M. M. Sheikh-Jabbari, Phys. Rev. D 66 (2002), 025025 doi:10.1103/PhysRevD.66.025025 [arXiv:hep-th/0110291 [hep-th]].
  * [7] J. Bagger and N. Lambert, JHEP 02 (2008), 105
  * [8] K. Skenderis and M. Taylor, Phys. Rept. 467 (2008), 117-171
  * [9] S. Baez, A. P. Balachandran, B. Ydri and S. Vaidya, Commun. Math. Phys. 208 (2000), 787-798
  * [10] R. J. Szabo, Phys. Rept. 378 (2003), 207-299
  * [11] H. Steinacker, PoS QGQGS2011 (2011), 004
  * [12] S. S. Gubser and S. L. Sondhi, Nucl. Phys. B 605 (2001), 395-424
  * [13] S. Kováčik and D. O’Connor, JHEP 10 (2018), 010
  * [14] B. Ydri, JHEP 03 (2014), 065
  * [15] M. Panero, JHEP 05 (2007), 082
  * [16] F. Garcia Flores, X. Martin and D. O’Connor, Int. J. Mod. Phys. A 24 (2009), 3917-3944
  * [17] T. Kuroki, Nucl. Phys. B 543 (1999), 466-484 doi:10.1016/S0550-3213(98)00815-3 [arXiv:hep-th/9804041 [hep-th]].
  * [18] W. Zhu, C. Han, E. Huffman, J. S. Hofmann and Y. C. He, Phys. Rev. X 13 (2023) no.2, 021009 doi:10.1103/PhysRevX.13.021009 [arXiv:2210.13482 [cond-mat.stat-mech]].
  * [19] A. B. Hammou, M. Lagraa and M. M. Sheikh-Jabbari, _Coherent state induced star product on R**3(lambda) and the fuzzy sphere,_ Phys. Rev. D 66 (2002), 025025, [arXiv:hep-th/0110291 [hep-th]].
  * [20] P. Vitale and J. C. Wallet, _Noncommutative field theories on $R^{3}_{\lambda}$: Toward UV/IR mixing freedom,_ JHEP 04 (2013), 115, [arXiv:1212.5131 [hep-th]].
  * [21] H. C. Steinacker and J. Tekel, JHEP 06 (2022), 136 doi:10.1007/JHEP06(2022)136 [arXiv:2203.02376 [hep-th]].
  * [22] J. Tekel, Acta Phys. Slov. 65 (2015) no.5, 369-468 [arXiv:1512.00689 [hep-th]].
  * [23] S. Kováčik and J. Tekel, [arXiv:2309.00576 [hep-th]].


# Higher-Order Allan Variance for Atomic Clocks of Arbitrary Order:
Mathematical Foundation

Takayuki Ishizaki1, Taichi Ichimura1, Takahiro Kawaguchi2, Yuichiro Yano3,
Yuko Hanado3 1 Tokyo Institute of Technology, 2-12-1, Ookayama, Meguro, Tokyo,
152-8552, Japan2 Gunma University, 1-5-1, Tenjincho, Kiryu, Gunma, 376-8515,
Japan3 National Institute of Information and Communications Technology, 4-2-1,
Nukui-Kitamachi, Koganei, Tokyo, 184-8795 , JapanManuscript received April 19,
2005; revised August 26, 2015.

###### Abstract

In this paper, we perform a time-domain analysis of the higher-order Allan
variance for atomic clock models of arbitrary order. Adopting a standard
atomic clock model where the time series of the clock reading deviation is
expressed as a Wiener or integrated Wiener process, we define the higher-order
Allan variance as the mean squared higher-order difference of the clock
reading deviation. The main results of this paper are threefold. First, we
prove that the higher-order difference operation of the clock reading
deviation, which can be interpreted as a linear aggregation with binomial
coefficients, is not only sufficient but also necessary for a resulting
aggregated time series to be an independent and identically distributed
Gaussian process. Second, we derive a complete analytical expression of the
higher-order Allan variance, which consists of both time-dependent and time-
independent terms. Third, we prove that the higher-order Allan variance is
time-independent if and only if the order of difference operation is greater
than or equal to the order of the atomic clock model.

###### Index Terms:

Higher-order Allan variance, higher-order difference, atomic clock models of
arbitrary order, analytical expression.

##  1 Introduction

Atomic clocks use the quantum transitions of atoms and ions as their frequency
reference to output extremely precise and accurate waves at a constant
frequency [1, 2, 3, 4]. This characteristic of atomic clocks supports modern
society in a wide range of fields where time and frequency are important[5].
Examples include communications, positioning, broadcasting, undersea
exploration, finance, space exploration, and geophysics[6, 7, 8].

Frequency stability quantifies the performance of a clock, which must be
evaluated statistically due to the stochastic nature of atomic clocks. In
order to assess frequency stability from time series data, a clock statistical
parameter is estimated, assuming that the corresponding stochastic process has
“ergodicity,” a property such that a time average is equal to an ensemble
average [9]. Ergodicity also requires that the stochastic process be at least
weakly stationary. Without such stationarity, frequency domain analysis based
on the Wiener-Khinchin theorem, and existing measures of frequency stability
such as the Allan variance and its variants [10], the total variance [11], and
the parabolic variance [12], are invalid. Therefore, a rational estimation of
the clock statistical parameter requires that the time series data of the
clock reading deviation be made into a stationary process.

It is well known that the clock reading deviation data of atomic clocks is not
stationary. In particular, it was recently shown that small atomic clocks have
not only linear frequency drift, but also non-negligible nonlinear frequency
drift [13]. These drift components confirm that it is not rational to estimate
a clock statistical parameter directly from non-stationary time series data,
because the estimated value does not converge to the actual parameter even in
the limit of increasing the number of data. Therefore, a rational estimate
includes generating stationary data by appropriately transforming the obtained
non-stationary data while preserving equivalent information about the
estimated clock statistical parameter.

The first contribution of this paper proves the “uniqueness” of the linear
aggregation that makes the clock reading deviation data stationary. One
possible transformation that can make non-stationary data stationary is to
take weighted sums of the data divided into fixed intervals, i.e., to apply
linear aggregation to the time series data. In this paper, we mathematically
prove that the aggregation weights, such that the resulting aggregated time
series is an independent and identically distributed Gaussian process
(implying ergodicity), are essentially unique. Interestingly, such weights can
be found as binomial coefficients. This result is stated in Theorem 1.

Theorem 1 shows that the unique aggregation operation consists of binomial
coefficients and can be interpreted as a higher-order difference operation,
which supports the rationality of existing frequency stability measures in
terms of “necessity.” For example, the Allan variance applies a second-order
difference operation to define a clock statistical parameter, and the Hadamard
variance applies a third-order difference operation. From these facts, it can
be deduced from mathematical necessity that the existing measures are rational
in the sense of making the clock reading deviation data stationary. However,
one should keep in mind that it is necessary to apply a higher-order
difference operation according to the order of the atomic clock, i.e., the
number of its differential equations.

The second contribution of this paper is the derivation of a complete
analytical expression of the higher-order Allan variance, which is a clock
statistical parameter defined as the mean squared higher-order difference of
the clock reading deviation. The derived expression is a generalization of the
Allan and Hadamard variances in [14]. The definition of the higher-order Allan
variance is consistent with the Allan and Hadamard variances in that the
dependence of clock noise parameters on the averaging time is identical, and
the coefficient for the white-frequency noise is normalized. This result is
stated in Theorem 2.

The higher-order Allan variance is expected to be effective for analyzing the
order of atomic clocks based on the separation of noise components in long-
term data. In addition, the analytical expression derived in this paper can be
used to verify the reliability of numerical analysis, since it is independent
of numerical errors in difference operations. We note that especially for
higher-order difference operations, numerical errors due to digit drop
generally become pronounced leading to unreliabile numerical calculations.

The higher-order Allan variance, deduced from the necessity for
stationarization, is essentially equivalent to the higher-order structure
function of the clock reading deviation. The structure function is a classical
concept introduced by Kolmogorov in 1941 for the analysis of turbulence [15,
16], and is still used as a powerful tool in this field [17]. It has also been
shown that the structure function corresponds to a generalization of the Allan
variance [9, 18]. Theorem 2 in this paper derives a complete analytical
expression of the higher-order structure function of the clock reading
deviation as a function of clock noise parameters from the time-domain state-
space model. We note that the evaluation of oscillator stability by structure
functions in previous studies assumes that the power spectral density of the
noise components is given in the frequency domain. This is clearly different
from the modeling and analysis in the time domain of this paper.

The third and final contribution of this paper is to derive a necessary and
sufficient condition for the higher-order Allan variance to be independent of
time. In particular, we prove that the higher-order Allan variance is time-
independent if and only if the order of the difference operation is equal to
or larger than the order of the atomic clock. This result is stated in Theorem
3. Although structure functions for time series data of clock reading
deviations can be made stationary by a higher-order difference operation, to
the best of the authors’ knowledge, there is no result that mathematically
proves that a higher-order difference operation is “necessary” for the
variance to be time-independent.

The remainder of this paper is organized as follows. In Section 2, we
introduce a standard atomic clock model of arbitrary order, where the time
series of the clock reading deviation is expressed as a Wiener or integrated
Wiener process. An equivalent discrete-time model is also introduced. In
Section 3, we mathematically prove that the higher-order difference operation
of the clock reading deviation is necessary and sufficient for a discrete-time
aggregated process to be an independent and identically distributed Gaussian
process. In Section 4, we derive the complete higher-order Allan variance
analytical expression, based on the fact proven in Section 3, for atomic clock
models of arbitrary order. Furthermore, we analyze the relationship between
the difference operation and atomic clock model orders with application
examples provided at the end of the section. Section 5 concludes this paper
with a summary of the results and brief description of future work.

Notation We denote the set of real numbers by $\mathbb{R}$, the $i$th
canonical unit vector, i.e., the $i$th column of the identity matrix, by
$e_{i}$, the transpose of a matrix $A$ by $A^{\sf T}$, the subspace spanned by
vectors $v_{1},\ldots,v_{n}$ by
$\operatorname{\mathsf{span}}\\{v_{1},\ldots,v_{n}\\}$, the ensemble mean of a
random variable $X$ by $\mathbb{E}[X]$, the Dirac delta function by
$\delta(\cdot)$, and the binomial coefficient by

$\left(\;\begin{matrix}n\\\
k\end{matrix}\>\right):=\frac{n(n-1)\cdots(n-k+1)}{k!}.$

The equality “$:=$” is used in the sense that a new symbol on the left-hand
side is defined by known symbols on the right-hand side.

##  2 Description of Atomic Clock Models

###  2.1 Continuous-Time Models

First, we introduce a standard atomic clock model. From experimental
observations, it is known [19, 20] that the time series of clock reading
deviation relative to an ideal time coordinate can be mathematically modeled
as a continuous-time stochastic process

$\begin{split}y(t)&=\sum_{i=1}^{n}\frac{c_{i}}{(i-1)!}t^{i-1}\\\
&+\sum_{i=1}^{n}\int_{0}^{t}\int_{0}^{t_{1}}\cdots\int_{0}^{t_{i-1}}v_{i}(t_{i})dt_{i}\cdots
dt_{2}dt_{1}\end{split}$ (1)

where $c_{i}$ is a scalar constant corresponding to an initial value, $n$
denotes the order of the atomic clock model, and

$v_{i}=\\{v_{i}(t):t\geq 0\\}$

denotes a continuous-time white Gaussian process such that

$\begin{split}&\mathbb{E}\bigl{[}v_{i}(t)\bigr{]}=0\\\
&\mathbb{E}\bigl{[}v_{i}(t)v_{i}(t+\tau)\bigr{]}=q_{i}^{2}\delta(\tau)\\\
&\mathbb{E}\bigl{[}v_{i}(t)v_{j}(t+\tau)\bigr{]}=0,\quad\forall\tau\in\mathbb{R},\quad
i\neq j.\end{split}$ (2)

The integral equation model in (1) can be represented by a system of
stochastic differential equations. In particular, we express its state
equation as

$\left\\{\>\begin{aligned} \dot{x}_{1}(t)&=x_{2}(t)+v_{1}(t)\\\
\dot{x}_{2}(t)&=x_{3}(t)+v_{2}(t)\\\ &~{}~{}\vdots\\\
\dot{x}_{n-1}(t)&=x_{n}(t)+v_{n-1}(t)\\\
\dot{x}_{n}(t)&=v_{n}(t).\end{aligned}\right.$ (3)

If the initial values are given such that

$x_{i}(0)=c_{i},\quad i=1,2,\ldots,n,$

then the first state variable $x_{1}$ satisfies

$x_{1}(t)=y(t),\quad\forall t\geq 0.$

Furthermore, we denote the stacked vectors by

$x(t):=\left[\>\begin{matrix}x_{1}(t)\\\ \vdots\\\
x_{n}(t)\end{matrix}\>\right],\quad v(t):=\left[\>\begin{matrix}v_{1}(t)\\\
\vdots\\\ v_{n}(t)\end{matrix}\>\right].$

Then, the stochastic differential equation system model in (3) can be written
in the matrix form as

$\Sigma:\left\\{\>\begin{aligned} \dot{x}(t)&=Ax(t)+v(t)\\\
y(t)&=Cx(t)\end{aligned}\right.$ (4)

where $A\in\mathbb{R}^{n\times n}$ and $C\in\mathbb{R}^{1\times n}$ are
defined as

$A:=\left[\>\begin{matrix}0&1&0&\cdots&0\\\ 0&0&1&\cdots&0\\\
\vdots&&\ddots&\ddots&\vdots\\\ \vdots&&&0&1\\\
0&0&\cdots&\cdots&0\end{matrix}\>\right],\quad
C:=\left[\>\begin{matrix}1&0&\cdots&0\end{matrix}\>\right].$ (5)

The system noise $v$ is a white Gaussian process such that

$\mathbb{E}\bigl{[}v(t)\bigr{]}=0,\quad\mathbb{E}\bigl{[}v(t)v^{\sf
T}(t+\tau)\bigr{]}=Q\delta(\tau),\quad\forall t\geq 0.$

where the covariance matrix is diagonal and given as

$Q:=\mathsf{diag}(q_{1}^{2},\ldots,q_{n}^{2}).$

In the rest of this paper, we assume that all noise variances
$q_{1}^{2},\ldots,q_{n}^{2}$ are nonzero without loss of generality.

###  2.2 Equivalent Discrete-Time Models

For the discussion below, we introduce the following standard notion of a
discrete-time stochastic process [21].

###### Definition 1

A discrete-time stochastic process

$f=\\{f[k]:k=0,1,\ldots\\}$

is said to be an independent and identically distributed Gaussian process (IID
Gaussian process) if

$\mathbb{E}\bigl{[}f[k]f^{\sf T}[k^{\prime}]\bigr{]}=0,\quad\forall k\neq
k^{\prime},$

and every $f[k]$ follows an identical Gaussian distribution.

An important property of the IID Gaussian process is that the statistics about
ensemble can be characterized by the statistics about time, i.e., the
ergodicity is satisfied. In particular, it follows from the the law of large
numbers that

$\frac{1}{K}\sum_{k=0}^{K-1}f[k]\xrightarrow{{\rm
a.s.}}\mathbb{E}\bigl{[}f[k]\bigr{]}$

as $K$ goes to infinity for every time instant $k$, where “a.s.” indicates
“almost surely.” This property is fundamental to estimate statistical
parameters, such as mean and variance, from an observed time series data, as
pointed out in [9].

Consider a discrete-time sequence $\\{t_{0},t_{1},\ldots\\}$ in an ideal time
coordinate. The interval between two adjacent times is assumed to be constant,
or equivalently

$t_{k}:=\tau k,\quad k=0,1,\ldots,$ (6)

where $\tau$ denotes a sampling period. The relation between $x(t_{k})$ and
$x(t_{k+1})$ in (4) can be formally expressed as

$x(t_{k+1})=e^{A\tau}x(t_{k})+\int_{t_{k}}^{t_{k+1}}e^{A(t_{k+1}-t)}v(t)dt.$

Based on this expression, we can derive a discrete-time model that is
equivalent to the continuous-time model in (4) over the discrete time sequence
as

$\left\\{\>\begin{aligned} x[k+1]&=A_{\tau}x[k]+v[k]\\\
y_{\tau}[k]&=Cx[k]\end{aligned}\right.$ (7)

where the discrete-time transition matrix $A_{\tau}$ is defined as

$A_{\tau}:=e^{A\tau},$ (8)

the discrete-time state $x$ and clock reading deviation $y_{\tau}$ as

$x[k]:=x(t_{k}),\quad y_{\tau}[k]:=y(t_{k}).$

Furthermore, the system noise

$v=\\{v[k]:k=0,1,\ldots\\}$

is an IID Gaussian process such that every $v[k]$ follows the Gaussian
distribution such that

$\mathbb{E}\bigl{[}v[k]\bigr{]}=0,\quad\mathbb{E}\bigl{[}v[k]v^{\sf
T}[k]\bigr{]}=Q_{\tau},\quad\forall k=0,1,\ldots$

where the covariance matrix is given as

$\displaystyle Q_{\tau}:=\int_{0}^{\tau}e^{At}Qe^{A^{\sf T}t}dt.$ (9)

It should be noted that, because of the structure of $A$ in (5), the power of
$A_{\tau}$ can be represented as

$A_{\tau}^{k}=\sum_{m=0}^{n-1}\frac{(\tau k)^{m}}{m!}A^{m},$ (10)

which is proven by the fact that

$A^{m}=0,\quad\forall m\geq n.$

Furthermore, the characteristic polynomial of $A_{\tau}$ is

$\operatorname{\mathsf{det}}(\lambda I-A_{\tau})=(\lambda-1)^{n}.$

This can be proven by the fact that $A_{\tau}$ is an upper triangular matrix
whose diagonal elements are all one. Thus, the lowest order annihilator
polynomial

$(A_{\tau}-I)^{n}=0$ (11)

is obtained by Cayley-Hamilton theorem [22].

##  3 Stationarizing Operation

###  3.1 Formulation

Consider the discrete-time stochastic process

$y_{\tau}:=\\{y_{\tau}[k]:k=1,2,\ldots\\}.$

It is clear from the integral equation in (1) that this $y_{\tau}$ is not a
stationary process because both mean and covariance are time dependent.
Therefore, in the following, we consider “stationarizing” it by linear
aggregation. To this end, we introduce an aggregated discrete-time process as

$Y_{N}:=\Biggl{\\{}\underbrace{\sum_{i=0}^{N}\alpha_{i}y_{\tau}\bigl{[}N(K+1)-i\bigr{]}}_{Y_{N}[K]}:K=0,1,\ldots\Biggr{\\}}$
(12)

where $\\{\alpha_{0},\ldots,\alpha_{N}\\}$ is the set of nonzero weights for
the linear aggregation, and $N$ specifies the number of signals to be
aggregated. This aggregation process is depicted in Fig. 1, where $N$ is
supposed to be 3.

Figure 1: Construction of aggregated discrete-time process.

By the definition of $Y_{N}[K]$ and the state equation of $\Sigma$ in (4), we
see that

$\begin{split}Y_{N}[K]&=\underbrace{\left(\sum_{i=0}^{N}\alpha_{i}CA_{\tau}^{N-i}\right)x[NK]}_{X_{N}[K]}\\\
&+\underbrace{\sum_{i=0}^{N-1}\sum_{j=0}^{i}\alpha_{j}CA_{\tau}^{i-j}v[N(K+1)-1-i]}_{V_{N}[K]}.\end{split}$
(13)

This expression will be used for the statistic analysis of the aggregated
series $Y_{N}$. It should be noted that the first and second terms are
uncorrelated, namely

$\mathbb{E}\bigl{[}X_{N}[K]V_{N}[K]\bigr{]}=0,\quad\forall K=0,1,\ldots$ (14)

because those are, respectively, linear combinations of

$\\{v[0],\ldots,v[NK-1]\\},\quad\\{v[NK],\ldots,v[N(K+1)-1]\\},$

which are independent each other.

###  3.2 Analysis

In this subsection, we analyze a condition that should be satisfied by the
weight set for the aggregated process $Y_{N}$ in (12) to be an IID Gaussian
process. Such a weight set is essentially unique as shown in the following
theorem.

###### Theorem 1

For a an atomic clock model $\Sigma$ in (4), consider the discrete-time
aggregated process $Y_{N}$ in (12). Then, there exists a weight set
$\\{\alpha_{0},\ldots,\alpha_{N}\\}$ such that $Y_{N}$ is an IID Gaussian
process for any initial values of $\Sigma$ if and only if $N$ is greater than
or equal to $n$. In particular, if $N$ is equal to $n$, then such a weight set
is uniquely determined except for constant multiples, and is found to be

$\alpha_{i}=\left(\;\begin{matrix}N\\\ i\end{matrix}\>\right)(-1)^{i},\quad
i=0,\ldots,N.$ (15)

###### Proof:

We first prove that, if $N\geq n$, then $Y_{N}$ with the weight set in (15) is
found to be an IID Gaussian process. For $X_{N}[K]$ in (13), it follows that

$\sum_{i=0}^{N}\left(\;\begin{matrix}N\\\
i\end{matrix}\>\right)(-1)^{i}A_{\tau}^{N-i}=\underbrace{(A_{\tau}-I)^{n}}_{0}(A_{\tau}-I)^{N-n}=0,$

where (11) is used. Furthermore, $V_{N}[K]$ in (13) is a linear combination of
white Gaussian noises with zero means and constant variances. Hence, $Y_{N}$
is proven to be an IID Gaussian process for any initial conditions.

Next, we prove that, if $Y_{N}$ for $N=n$ is supposed to be an IID Gaussian
process, then the weight set is uniquely determined. To this end, we first
prove that, if $Y_{n}$ is an IID Gaussian process, then it follows that

$\beta:=\sum_{i=0}^{n}\alpha_{i}CA_{\tau}^{n-i}$

is zero. Denoting the $i$th element of $\beta$ by $\beta_{i}$, we have

$\mathbb{E}\bigl{[}Y_{n}[K]\bigr{]}=\sum_{i=1}^{n}\beta_{i}\mathbb{E}\bigl{[}x_{i}[nK]\bigr{]},$
(16)

where the mean of $x_{i}[nK]$ is found to be

$\mathbb{E}\bigl{[}x_{i}[nK]\bigr{]}=\sum_{m=0}^{n-i}\frac{(\tau
nK)^{m}}{m!}c_{i+m}.$

Note that only the mean of $x_{n}[nK]$ is constant. This implies that, if any
of $\beta_{1},\ldots,\beta_{n-1}$ is not zero, then there exists some initial
value such that (16) depends on $K$. Thus, all $\beta_{1},\ldots,\beta_{n-1}$
must be zero if $Y_{n}$ is an IID Gaussian process. Furthermore, in the case
of only $\beta_{n}$ is nonzero, we have

$\mathbb{E}\left[\left(Y_{n}[K]-\mathbb{E}\bigl{[}Y_{n}[K]\bigr{]}\right)^{2}\right]=\beta_{n}^{2}\left(\mathbb{E}\bigl{[}x_{n}^{2}[nK]\bigr{]}-c_{n}^{2}\right)+c$

where $c$ is a constant representing the variance of $V_{n}[K]$ in (13). This
depends on $K$ because $x_{n}$ is a random walk process having a time-
dependent variance. This proves that $\beta$ must be zero if $Y_{n}$ is an IID
Gaussian process.

In addition, we prove that, if $\beta$ is zero, then (15) is the unique weight
set, except for constant multiples. Using the expression in (10), we can
rewrite $\beta$ as

$\beta=\sum_{i=0}^{n}\sum_{m=0}^{n-1}\alpha_{i}\frac{\\{(n-i)\tau\\}^{m}}{m!}CA^{m}.$

Thus, the fact that $\beta$ is zero can be written as

$\underbrace{\left[\>\begin{matrix}1&1&\cdots&1\\\ 0^{1}&1^{1}&\cdots&n^{1}\\\
0^{2}&1^{2}&\cdots&n^{2}\\\ \vdots&\vdots&\vdots&\vdots\\\
0^{n-1}&1^{n-1}&\cdots&n^{n-1}\end{matrix}\>\right]}_{W_{n}}\left[\>\begin{matrix}\alpha_{n}\\\
\alpha_{n-1}\\\ \alpha_{n-2}\\\ \vdots\\\ \alpha_{0}\end{matrix}\>\right]=0.$

We notice that $W_{n}\in\mathbb{R}^{n\times(n+1)}$ is a Vandermonde matrix
[22] being of row full rank. Therefore, the kernel of $W_{n}$ is necessarily
one-dimensional. This proves that the weight set such that $\beta$ is zero is
unique, except for constant multiples.

Finally, we prove that there exists a weight set such that $Y_{N}$ is an IID
Gaussian process only if $N\geq n$. This is equivalent to prove that, if
$N<n$, then, for every nonzero weight set, $Y_{N}$ is not an IID Gaussian
process for some initial value. Analyzing a resultant Vandermonde matrix, we
can verify that there does not exist a weight set such that

$\sum_{i=0}^{N}\alpha_{i}CA_{\tau}^{N-i}=0.$

Therefore, at least either of the mean or variance of $Y_{N}[K]$ depends on
$K$. This proves the claim. ∎

Theorem 1 shows that there exists an essentially unique weight set for the
aggregated process $Y_{N}$ in (12) to be an IID Gaussian process. Furthermore,
the number $N$ for the aggregation must be greater than or equal to the order
$n$ of the atomic clock model. The major novelty of this theorem is to show
that, starting from the discussion of “necessity” for stationarizing the time
series of the sampled clock reading deviation, the linear aggregation to apply
is uniquely found as the binomial coefficient in (15).

##  4 Higher-Order Allan Variance

###  4.1 Definition

In this subsection, we introduce a generalized notion of the Allan variance
for higher-order atomic clocks, based on the fact that the linear aggregation
with the weight set in (15) can be interpreted in terms of a “higher-order
difference operation.” To this end, let us introduce the operation of
difference for a continuous time series $f(\cdot)$. In particular, the first
order difference of $f(\cdot)$ is denoted as

$\triangle_{\tau}f(t):=f(t+\tau)-f(t).$

Then, the $N$th order difference is recursively defined as

$\triangle_{\tau}^{N}f(t):=\triangle_{\tau}\left\\{\triangle_{\tau}^{N-1}f(t)\right\\}.$
(17)

Using this difference operation, $Y_{N}[K]$ in (12) with the weight set in
(15) can be represented as

$Y_{N}[K]=\triangle_{\tau}^{N}y(t_{NK}).$

From this viewpoint, Theorem 1 states that the higher-order difference of the
sampled clock reading deviation is found to be an IID Gaussian process, and
conversely, the higher-order difference must be used for the aggregation of
the sampled clock reading deviation to be an IID Gaussian process.

Based on the notion of this higher-order difference, we introduce the
following generalization of the Allan variance.

###### Definition 2

Consider the atomic clock model $\Sigma$ in (4). The higher-order Allan
variance of $\Sigma$ is defined by

$\sigma^{2}_{N}(\tau;t):=\frac{1}{r_{0}^{(N)}\tau^{2}}\mathbb{E}\left[\left\\{\triangle^{N}_{\tau}y(t)\right\\}^{2}\right],$
(18)

where $\triangle^{N}_{\tau}$ is the $N$th order difference operator in (17),
and

$r_{0}^{(N)}:=\sum_{i=0}^{N-1}\sum_{j_{1}=0}^{i}\sum_{j_{2}=0}^{i}(-1)^{j_{1}+j_{2}}\left(\;\begin{matrix}N\\\
j_{1}\end{matrix}\>\right)\left(\;\begin{matrix}N\\\
j_{2}\end{matrix}\>\right).$

We remark that $y(t)$ represents the clock reading deviation relative to an
ideal time coordinate in Definition 2. The constant $r_{0}^{(N)}$ is
introduced for normalization. The higher-order Allan variance coincides with
the Allan and Hadamard variances for $N=2$ and $N=3$, respectively. In
particular, the Allan variance is obtained as

$\sigma_{2}^{2}(\tau;t)=\frac{1}{2\tau^{2}}\mathbb{E}\left[\left\\{y(t+2\tau)-2y(t+\tau)+y(t)\right\\}^{2}\right],$

and the Hadamard variance is obtained as

$\begin{split}\sigma_{3}^{2}(\tau;t)=\frac{1}{6\tau^{2}}\mathbb{E}\Bigl{[}\bigl{\\{}y(t+3\tau)&-3y(t+2\tau)\\\
&+3y(t+\tau)-y(t)\bigr{\\}}^{2}\Bigr{]}.\end{split}$

Interestingly, the higher-order Allan variance, mathematically deduced from
the necessity for stationarization, also coincides with the higher-order
structure function of the clock reading deviation, except for the difference
of constant multiples [9]. This coincidence affirms the rationality of
existing oscillator stability measures. We further remark that the
normalization constant $r_{0}^{(N)}$ in (18), which appears in the Allan and
Hadamard variances, is first found by our time-domain analysis in Theorem 2
below.

It has been pointed out in [14] that the Allan variance of third order atomic
clock models depends not only on the sampling period $\tau$, but also on the
time $t$, while the Hadamard variance depends only on $\tau$. This can be
explained by the fact that the mean and variance of the second-order
difference of the clock reading deviation $y(t)$ is time dependent while those
of the third-order difference is not. In the following, we aim at generalizing
such a fact for atomic clock models of arbitrary order.

###  4.2 Analytical Expression

An analytical expression of the higher-order Allan variance can be derived as
in the following theorem.

###### Theorem 2

Consider the atomic clock model $\Sigma$ in (4). The higher-order Allan
variance of $\Sigma$ is found to be

$\begin{split}\sigma^{2}_{N}(\tau;t)=&\sum_{i=1}^{n}\sum_{j=1}^{n}\frac{\gamma_{i}^{(N)}\gamma_{j}^{(N)}}{r_{0}^{(N)}}\mathbb{E}\bigl{[}x_{i}(t)x_{j}(t)\bigr{]}\\\
+&\sum_{m=0}^{n-1}\frac{r_{m}^{(N)}}{r_{0}^{(N)}}q_{1+m}^{2}\tau^{2m-1}\end{split}$
(19)

where $\gamma_{m}^{(N)}$ is a constant defined as

$\gamma_{m}^{(N)}:=\frac{\tau^{m-2}}{(m-1)!}\sum_{i=0}^{N}\left(\;\begin{matrix}N\\\
i\end{matrix}\>\right)(-1)^{N-i}i^{m-1},$

and $r_{m}^{(N)}$ is a constant defined as

$\begin{split}r_{m}^{(N)}:=&\left(\frac{1}{m!}\right)^{2}\sum_{i=0}^{N-1}\sum_{j_{1}=0}^{i}\sum_{j_{2}=0}^{i}(-1)^{j_{1}+j_{2}}\left(\;\begin{matrix}N\\\
j_{1}\end{matrix}\>\right)\left(\;\begin{matrix}N\\\
j_{2}\end{matrix}\>\right)\\\ \times&\sum_{l=0}^{m}\left(\;\begin{matrix}m\\\
l\end{matrix}\>\right)\bigl{\\{}(i-j_{1})(i-j_{2})\bigr{\\}}^{m-l}\\\
\times&\sum_{k=0}^{l}\left(\;\begin{matrix}l\\\
k\end{matrix}\>\right)\frac{(2i-j_{1}-j_{2})^{k}}{2l-k+1}.\end{split}$

###### Proof:

We consider the case of $t=t_{NK}$. This supposition does not lose generality
because we can arbitrarily set the origin of the time variable such that $t$
is equal to $\tau NK$ with some $K$. From (13) and (14), we can see that

$r_{0}^{(N)}\tau^{2}\sigma^{2}_{N}(\tau;t_{NK})=\mathbb{E}\bigl{[}X_{N}^{2}[K]\bigr{]}+\mathbb{E}\bigl{[}V_{N}^{2}[K]\bigr{]}.$

In the following, we evaluate these terms. Using the expression in (10) and
the weight set in (15), we have

$\begin{split}X_{N}[K]&=\sum_{i=0}^{N}\alpha_{i}CA_{\tau}^{N-i}x(t_{NK})\\\
&=\sum_{i=0}^{N}\alpha_{N-i}CA_{\tau}^{i}x(t_{NK})\\\
&=\tau\sum_{m=1}^{n}\gamma_{m}^{(N)}x_{m}(t_{NK}).\end{split}$

Therefore, we obtain

$\mathbb{E}\bigl{[}X_{N}^{2}[K]\bigr{]}=\tau^{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\gamma_{i}^{(N)}\gamma_{j}^{(N)}\mathbb{E}\bigl{[}x_{i}(t_{NK})x_{j}(t_{NK})\bigr{]},$

which corresponds to the first term of (19).

We evaluete the mean of $V_{N}^{2}[K]$. From (13), we have

$\begin{split}\mathbb{E}\bigl{[}V_{N}^{2}[K]\bigr{]}&=\sum_{i=0}^{N-1}\underbrace{\Biggl{(}\sum_{j=0}^{i}\alpha_{j}CA_{\tau}^{i-j}\Biggr{)}}_{X_{i}}Q_{\tau}\Biggl{(}\sum_{j=0}^{i}\alpha_{j}CA_{\tau}^{i-j}\Biggr{)}^{\sf
T}\\\
&=\sum_{i=0}^{N-1}\int^{\tau}_{0}X_{i}e^{At}Q\left(X_{i}e^{At}\right)^{\sf
T}dt,\end{split}$

where we have used the definition of $Q_{\tau}$ in (9). Note that

$\begin{split}X_{i}e^{At}&=\sum_{j=0}^{i}\alpha_{j}e_{1}^{\sf
T}e^{A\left\\{(i-j)\tau+t\right\\}}\\\
&=\sum_{j=0}^{i}\alpha_{j}\sum_{m=0}^{n-1}\frac{\left\\{(i-j)\tau+t\right\\}^{m}}{m!}e_{1+m}^{\sf
T}.\end{split}$

Thus, using the fact that

$e_{i}^{\sf T}Qe_{j}=\left\\{\begin{array}[]{cl}q_{i}^{2},&i=j\\\ 0,&{\rm
otherwise},\end{array}\right.$

we obtain

$\begin{split}\mathbb{E}\bigl{[}V_{N}^{2}[K]\bigr{]}&=\sum_{i=0}^{N-1}\sum_{m=0}^{n-1}\int^{\tau}_{0}\Biggl{(}\sum_{j=0}^{i}\alpha_{j}\tfrac{\left\\{(i-j)\tau+t\right\\}^{m}}{m!}q_{1+m}\Biggr{)}^{2}dt\\\
&=\sum_{i=0}^{N-1}\sum_{m=0}^{n-1}\left(\frac{q_{1+m}}{m!}\right)^{2}I_{im}\end{split}$

where $I_{im}$ denotes the integral term of

$I_{im}:=\int^{\tau}_{0}\Biggl{(}\sum_{j=0}^{i}\alpha_{j}\left\\{(i-j)\tau+t\right\\}^{m}\Biggr{)}^{2}dt.$

To evaluate the integral $I_{im}$, we use the expansions of

$\begin{split}\Biggl{(}\sum_{j=0}^{i}x_{j}\Biggr{)}^{2}&=\sum_{j_{1}=0}^{i}\sum_{j_{2}=0}^{i}x_{j_{1}}x_{j_{2}},\\\
(\alpha+\beta+\gamma)^{m}&=\sum_{l=0}^{m}\left(\;\begin{matrix}m\\\
l\end{matrix}\>\right)\alpha^{m-l}\sum_{k=0}^{l}\left(\;\begin{matrix}l\\\
k\end{matrix}\>\right)\beta^{l-k}\gamma^{k}.\end{split}$

Applying these equalities, we can rewrite the integrand in $I_{im}$ just as a
polynomial of $t$. Thus, calculating the integral, we have the second term of
(19). ∎

The analytical expression in Theorem 2 is composed of two terms. The first
term is a time-dependent term that is related to the variance and covariance
of atomic clock states. The second term is a time-independent term that is
represented as a polynomial of the sampling period $\tau$.

It is clear that whether the higher-order Allan variance is time dependent or
not depends on whether the first term is nonzero or not. In fact, an explicit
expression of the time-dependent term can also be derived as follows.

$N$ | $q_{1}^{2}\tau^{-1}$ | $q_{2}^{2}\tau$ | $q_{3}^{2}\tau^{3}$ | $q_{4}^{2}\tau^{5}$ | $q_{5}^{2}\tau^{7}$ | $q_{6}^{2}\tau^{9}$ | $q_{7}^{2}\tau^{11}$ | $q_{8}^{2}\tau^{13}$ | $q_{9}^{2}\tau^{15}$ | $q_{10}^{2}\tau^{17}$
---|---|---|---|---|---|---|---|---|---|---
2 | 1.0000e+00 | 3.3333e-01 |  |  |  |  |  |  |  |
3 | 1.0000e+00 | 1.6667e-01 | 9.1667e-02 |  |  |  |  |  |  |
4 | 1.0000e+00 | 1.3333e-01 | 3.3333e-02 | 2.3968e-02 |  |  |  |  |  |
5 | 1.0000e+00 | 1.1905e-01 | 2.2619e-02 | 6.9444e-03 | 6.1488e-03 |  |  |  |  |
6 | 1.0000e+00 | 1.1111e-01 | 1.8254e-02 | 4.1005e-03 | 1.4863e-03 | 1.5632e-03 |  |  |  |
7 | 1.0000e+00 | 1.0606e-01 | 1.5909e-02 | 3.0123e-03 | 7.7687e-04 | 3.2460e-04 | 3.9542e-04 |  |  |
8 | 1.0000e+00 | 1.0256e-01 | 1.4452e-02 | 2.4531e-03 | 5.2278e-04 | 1.5218e-04 | 7.2018e-05 | 9.9720e-05 |  |
9 | 1.0000e+00 | 1.0000e-01 | 1.3462e-02 | 2.1170e-03 | 3.9850e-04 | 9.4365e-05 | 3.0604e-05 | 1.6180e-05 | 2.5098e-05 |
10 | 1.0000e+00 | 9.8039e-02 | 1.2745e-02 | 1.8943e-03 | 3.2660e-04 | 6.7492e-05 | 1.7582e-05 | 6.2864e-06 | 3.6723e-06 | 6.3080e-06
TABLE I: Coefficients of higher-order Allan variance.

###### Lemma 1

Consider the atomic clock model $\Sigma$ in (4). Then

$\begin{split}\mathbb{E}\bigl{[}x_{i}(t)x_{j}(t)\bigr{]}&=\Biggl{(}\sum_{m=0}^{n-i}\frac{t^{m}}{m!}c_{i+m}\Biggr{)}\Biggl{(}\sum_{m=0}^{n-j}\frac{t^{m}}{m!}c_{j+m}\Biggr{)}\\\
&+\sum_{m=\operatorname*{\mathsf{max}}\\{i,j\\}}^{n}\rho_{ij}^{(m)}q_{m}^{2}t^{2m-i-j+1},\end{split}$
(20)

where $\rho^{(m)}_{ij}$ is a positive constant defined as

$\rho^{(m)}_{ij}:=\frac{1}{(m-i)!(m-j)!(2m-i-j+1)}.$

###### Proof:

The state $x_{i}(t)$ can be formally written as

$x_{i}(t)=\underbrace{e_{i}^{\sf
T}e^{At}x(0)}_{\overline{x}_{i}(t)}+\underbrace{\int_{0}^{t}e_{i}^{\sf
T}e^{A(t-\tau)}v(\tau)d\tau}_{\xi_{i}(t)}.$

Because the mean of $\xi_{i}(t)$ is zero, we have

$\mathbb{E}\bigl{[}x_{i}(t)x_{j}(t)\bigr{]}=\overline{x}_{i}(t)\overline{x}_{j}(t)+\mathbb{E}\bigl{[}\xi_{i}(t)\xi_{j}(t)\bigr{]}.$

The first term here corresponds to the first term of (20). On the other hand,
the second term can be written as

$\mathbb{E}\bigl{[}\xi_{i}(t)\xi_{j}(t)\bigr{]}=\int_{0}^{t}\underbrace{e_{i}^{\sf
T}e^{A(t-\tau)}}_{X_{i}(\tau)}Q\left(e_{j}^{\sf T}e^{A(t-\tau)}\right)^{\sf
T}d\tau.$

Note that $X_{i}(\tau)$ can be represented as

$X_{i}(\tau)=\sum_{m=i}^{n}\frac{(t-\tau)^{m-i}}{(m-i)!}e_{m}^{\sf T}.$

Therefore, we have

$\begin{split}\int_{0}^{t}X_{i}(\tau)QX_{j}^{\sf
T}(\tau)d\tau&=\int_{0}^{t}X_{\operatorname*{\mathsf{max}}\\{i,j\\}}(\tau)QX_{\operatorname*{\mathsf{max}}\\{i,j\\}}^{\sf
T}(\tau)d\tau\\\
&=\int_{0}^{t}\textstyle\sum_{m=\operatorname*{\mathsf{max}}\\{i,j\\}}^{n}\limits\tfrac{(t-\tau)^{2m-i-j}}{(m-i)!(m-j)!}q_{m}^{2}d\tau.\end{split}$

This integral corresponds to the second term of (20). ∎

From Lemma 1, we see that the variance and covariance of atomic clock states
are nonzero polynomials of $t$. In particular, the first term of (20)
corresponds to the deterministic trend of atomic clock states, which is
dependent on a part of the initial values $c_{1},\ldots,c_{n}$, while the
second term corresponds to their stochastic behavior, which is dependent on a
part of the noise variances $q_{1}^{2},\ldots,q_{n}^{2}$.

Furthermore, the following lemma proves an essential fact to determine the
time dependence of (19).

###### Lemma 2

The constant $\gamma_{m}^{(N)}$ in (19) is zero if and only if $m$ is less
than or equal to $N$.

###### Proof:

Let us first prove that

$\gamma_{m}^{(N)}=0,\quad\forall m=1,\ldots,N.$ (21)

From the binomial theorem, we have

$(z-1)^{N}=\sum_{i=0}^{N}\left(\;\begin{matrix}N\\\
i\end{matrix}\>\right)(-1)^{N-i}z^{i}.$ (22)

Considering $z=1$, we see that $\gamma_{1}^{(N)}$ is zero. Furthermore, the
differentiation of (22) with respect to $z$ leads to

$N(z-1)^{N-1}=\sum_{i=0}^{N}\left(\;\begin{matrix}N\\\
i\end{matrix}\>\right)(-1)^{N-i}iz^{i-1}.$

Considering $z=1$ again, we see that $\gamma_{2}^{(N)}$ is zero. Repeating
such a procedure proves (21).

Next, we prove that $\gamma_{m}^{(N)}$ is not zero for $m>N$. Using the
notation of (15), we define

$w_{m}:=\left[\>\begin{matrix}0^{m-1}\\\ 1^{m-1}\\\ \vdots\\\
N^{m-1}\end{matrix}\>\right],\quad\alpha:=\left[\>\begin{matrix}\alpha_{N}\\\
\alpha_{N-1}\\\ \vdots\\\ \alpha_{0}\end{matrix}\>\right].$

Note that $\gamma_{m}^{(N)}$ is a constant multiple of $\alpha^{\sf T}w_{m}$.
Then, we see that (21) is equivalent to

$\operatorname{\mathsf{span}}\\{\alpha\\}^{\perp}=\operatorname{\mathsf{span}}\\{w_{1},\ldots,w_{N}\\},$

where $\perp$ denotes the orthogonal complement. Furthermore, $w_{i}$ and
$w_{j}$ are linearly independent if $i\neq j$. Thus

$w_{m}\notin\operatorname{\mathsf{span}}\\{\alpha\\}^{\perp},\quad\forall
m=N+1,N+2,\ldots,$

or equivalently

$\alpha^{\sf T}w_{m}\neq 0,\quad\forall m=N+1,N+2,\ldots.$

Hence, the claim is proven. ∎

Lemma 2 in conjunction with Lemma 1 implies that the time-dependent term in
(19) disappears if and only if $N$ is greater than or equal to $n$, proving
the following theorem.

###### Theorem 3

Consider the atomic clock model $\Sigma$ in (4). The higher-order Allan
variance of $\Sigma$ is time independent if and only if $N$ is greater than or
equal to $n$.

###### Proof:

The sufficiency is directly proven by Theorem 2 and Lemma 2. This is because,
if $N\geq n$, then

$\gamma^{(N)}_{m}=0,\quad\forall m=1,\ldots,n.$

To prove the necessity, we suppose that $N<n$. From the analytical expression
in Lemma 1, we see that

$\sigma^{2}_{N}(\tau;t)=\tfrac{\left(\gamma_{N+1}^{(N)}\right)^{2}}{r_{0}^{(N)}}\\!\rho_{(N+1)(N+1)}^{(n)}q_{n}^{2}t^{2n-2N-1}\\!+o(t^{2n-2N-2}),$

which is time dependent because $\gamma_{N+1}^{(N)}$ is not zero as shown in
Lemma 2. Hence, the claim follows. ∎

###  4.3 Example

We demonstrate application examples of Theorem 2. Consider the third-order
atomic clock model. For the second-order difference, we have

$\sigma^{2}_{2}(\tau;t)=\frac{1}{2}(c_{3}^{2}+q_{3}^{2}t)\tau^{2}+q_{1}^{2}\tau^{-1}+\frac{1}{3}q_{2}^{2}\tau+\frac{23}{60}q_{3}^{2}\tau^{3}.$
(23a) Furthermore, for the third-order difference, we have
$\sigma^{2}_{3}(\tau;t)=q_{1}^{2}\tau^{-1}+\frac{1}{6}q_{2}^{2}\tau+\frac{11}{120}q_{3}^{2}\tau^{3}.$
(23b)

These two values coincide with the Allan variance and the Hadamard variance
shown in [14].

Figure 2: Logarithmic plot of coefficients of higher-order Allan variance.

For reference, we numerically calculate the coefficients of the time-
independent term of (19) varying the difference order $N$. We suppose that the
order of atomic clock models is also $N$. In this case, the higher-order Allan
variance is not time dependent, and is simplified as

$\sigma_{N}^{2}(\tau)=\sum_{m=0}^{N-1}\frac{r_{m}^{(N)}}{r_{0}^{(N)}}q_{1+m}^{2}\tau^{2m-1}.$

The result is shown in Tab. I, whose logarithmic plot is given in Fig. 2. This
result shows that coefficients with respect to noise variances for higher-
order atomic clock models are exponentially smaller. This trend is consistent
with the empirical fact that longer term measurements are required to assess
the accuracy of higher-order atomic clock models.

##  5 Concluding Remarks

In this paper, we have mathematically deduced the definition of the higher-
order Allan variance in terms of “necessity” to make the stochastic process of
clock reading deviation stationary. In particular, we have proven that the
aggregation weights are essentially unique for the resulting clock reading
deviation to be an independent and identically distributed Gaussian process.
Furthermore, we have derived a complete analytical expression of the higher-
order Allan variance for a standard atomic clock model of arbitrary order.

The main focus of this paper is to analyze a rational statistical parameter
for the atomic clock model of arbitrary order. Building an efficient algorithm
to estimate the statistical parameter and extending the results to the
overlapping version of the higher-order Allan variance is part of future work.

## Acknowledgments

This paper includes the results of research and development conducted by the
Ministry of Internal Affairs and Communications (MIC) under its ”Research and
Development for Expansion of Radio Resources (JPJ000254)” program.

## References

  * [1] Masao Takamoto, Feng-Lei Hong, Ryoichi Higashi, and Hidetoshi Katori.  An optical lattice clock.  Nature, 435(7040):321–324, 2005.
  * [2] N Huntemann, C Sanner, B Lipphardt, Chr Tamm, and E Peik.  Single-ion atomic clock with 3$\times$ 10- 18 systematic uncertainty.  Physical review letters, 116(6):063001, 2016.
  * [3] Peter DD Schwindt, Yuan-Yu Jau, Heather L Partner, Darwin K Serkland, Aaron Ison, Andrew McCants, Edward Winrow, John Prestage, James Kellogg, Nan Yu, et al.  Miniature trapped-ion frequency standard with 171 yb+.  In 2015 Joint Conference of the IEEE International Frequency Control Symposium & the European Frequency and Time Forum, pages 752–757. IEEE, 2015.
  * [4] John Kitching.  Chip-scale atomic devices.  Applied Physics Reviews, 5(3), 2018.
  * [5] James Camparo.  The rubidium atomic clock and basic research.  Physics today, 60(11):33–39, 2007.
  * [6] S Rentsch, J Thompson, R Adams, D Moore, and R Raborn.  100 ocean bottom nodes with their csac clock drift analysed for 91 days: three controlled experiments.  In 84th EAGE Annual Conference & Exhibition, volume 2023, pages 1–5. European Association of Geoscientists & Engineers, 2023.
  * [7] Tyler G Reid, Andrew M Neish, Todd F Walter, and Per K Enge.  Leveraging commercial broadband leo constellations for navigating.  In Proceedings of the 29th International Technical Meeting of the Satellite Division of the Institute of Navigation (Ion Gnss+ 2016). 29th International Technical Meeting of the Satellite Division of the Institute of Navigation (Ion Gnss+ 2016), Portland, Oregon, volume 12, pages 2016–16, 2016\.
  * [8] Ruxandra Bondarescu, Mihai Bondarescu, György Hetényi, Lapo Boschi, Philippe Jetzer, and Jayashree Balakrishna.  Geophysical applicability of atomic clocks: direct continental geoid mapping.  Geophysical Journal International, 191(1):78–82, 2012.
  * [9] William C Lindsey and Chak Ming Chie.  Theory of oscillator instability based upon structure functions.  Proceedings of the IEEE, 64(12):1652–1666, 1976.
  * [10] David W Allan and Judah Levine.  A historical perspective on the development of the allan variances and their strengths and weaknesses.  IEEE transactions on ultrasonics, ferroelectrics, and frequency control, 63(4):513–519, 2016.
  * [11] Charles A. Greenhall, Dave A. Howe, and Donald B. Percival.  Total variance, an estimator of long-term frequency stability.  IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 46:1183–1191, 1999.
  * [12] François Vernotte, Michel Lenczner, Pierre-Yves Bourgeois, and Enrico Rubiola.  The parabolic variance (pvar): A wavelet variance based on the least-square fit.  IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 63(4):611–623, 2015.
  * [13] Alan T Gardner and John A Collins.  A second look at chip scale atomic clocks for long term precision timing.  In OCEANS 2016 MTS/IEEE Monterey, pages 1–9. IEEE, 2016.
  * [14] Cristina Zucca and Patricia Tavella.  The clock model and its relationship with the allan and related variances.  IEEE transactions on ultrasonics, ferroelectrics, and frequency control, 52(2):289–296, 2005.
  * [15] A N Kolmogorov.  The local structure of turbulence in incompressible viscous fluid for very large reynolds numbers.  Doklady Akademiia Nauk SSSR, 30:301–305, 1941.
  * [16] A N Kolmogorov.  Turbulence and stochastic process: Kolmogorov’s ideas 50 years on.  Doklady Akademiia Nauk SSSR, 32:16–18, 1941.
  * [17] Alexandre Arenas and Alexandre J Chorin.  On the existence and scaling of structure functions in turbulence according to the data, 2006.
  * [18] Charles A. Greenhall.  The fundamental structure function of oscillator noise models, 12 1982\.
  * [19] Lorenzo Galleani.  A tutorial on the two-state model of the atomic clock noise.  Metrologia, 45(6):S175, 2008.
  * [20] Lorenzo Galleani and Patrizia Tavella.  Time and the kalman filter.  IEEE Control Systems Magazine, 30(2):44–65, 2010.
  * [21] Geoffrey Grimmett and David Stirzaker.  Probability and random processes.  Oxford university press, 2020.
  * [22] Dennis S Bernstein.  Matrix mathematics: theory, facts, and formulas.  Princeton University Press, 2009.