# Asymmetries in triboelectric charging: generalizing mosaic models to
different-material samples and sliding contacts

Galien Grosjean<EMAIL_ADDRESS>Scott Waitukaitis  Institute of
Science and Technology Austria (ISTA), Lab Building West, Am Campus 1, 3400
Klosterneuburg, Austria

(August 30, 2024)

###### Abstract

Nominally identical materials exchange net electric charge during contact
through a mechanism that is still debated. ‘Mosaic models’, in which surfaces
are presumed to consist of a random patchwork of microscopic donor/acceptor
sites, offer an appealing explanation for this phenomenon. However, recent
experiments have shown that global differences persist even between same-
material samples, which the standard mosaic framework does not account for.
Here, we expand the mosaic framework by incorporating global differences in
the densities of donor/acceptor sites. We develop an analytical model, backed
by numerical simulations, that smoothly connects the global and deterministic
charge transfer of different materials to the local and stochastic mosaic
picture normally associated with identical materials. Going further, we extend
our model to explain the effect of contact asymmetries during sliding,
providing a plausible explanation for reversal of charging sign that has been
observed experimentally.

##  I Introduction

Contact electrification (CE), also known as triboelectrification or
tribocharging, is the exchange of electric charge between objects after
touching or rubbing against each other. It occurs in settings as varied as
industrial powder flows [1], volcanic plumes [2], dust storms [3], or the
early stages of planet formation [4]. Fundamental issues, including the
identity of the charge carrier(s) and the property (or properties) driving
charge exchange, are still vigorously debated [5]. In addition to these
issues, one particularly confounding aspect of CE is that samples made of the
same material systematically exchange charge [6]. This is particularly
prevalent in granular systems, where large amount of charge can build up from
repeated collisions, leading for example to the electrification of dust storms
and thunderclouds [1, 2, 3, 4].

Figure 1:  Mosaic models. (a) In the standard mosaic formulation, two
surfaces, A and B, of size $L\times L$ are thought to consist of donor sites
with probability $p$, and acceptors with probability $1-p$. Sites have a size
$l_{0}$ and can be spatially correlated over a length $l$. During contact,
donors give a charge $e$ to facing acceptors, with a certain probability
$\alpha$. In this paper, we adapt this model in two ways. (b) First, we study
contact between samples of different donor probability, e.g. grid A vs. C with
$p_{\mathrm{A}}\neq p_{\mathrm{C}}$. (c) Second, we study situations were one
smaller grid slides a distance $d$ over another longer one, e.g. grid A vs. D.

Models for different-material CE and same-material CE are usually considered
separately, in large part due to a geometric distinction in their underlying
frameworks. Models for different-material CE often consider transfer of charge
due to an underlying global parameter, i.e. one that is spatially uniform over
the surface of a material. For example, one can point to the work function in
the case of different metals [7], or Lewis basicity as one candidate (of very
many) to explain the case of insulators [8]. On the other hand, models for
same-material CE often consider transfer due to an underlying parameter that
is local, i.e. allowed to vary over the surface of a material. Surfaces can be
thought of as consisting of a very large number of microscopic donor and
acceptor sites, dividing the surface of the material into a grid or mosaic.
During contact, if a donor on one surface touches an acceptor on the other,
charge exchange occurs locally, and net charge exchange becomes a stochastic
result of all such ‘micro-exchanges’ [9]. In Fig. 1, we illustrate this
process with two square grids, A and B, randomly populated with donor and
acceptor sites. The exact nature of the hypothesized, locally-varying charge-
driving parameter is unknown, but proposed candidates include adsorbed water
islands [10], material from the sample [11], patterns of strain in the
material [12], anions and cations produced by heterolytic cleavage of polymers
[13, 14], just to name a few. Regardless of what the parameter is, mosaic
models are appealing and widely favored because they permit identical
materials to be ‘identical’ in a global sense, but nonetheless still exchange
charge [9, 15, 16, 17, 18].

Mosaic models are qualitatively consistent with experimental data obtained
with Kelvin Probe Force Microscopy (KPFM) and the Scanning Kelvin Probe method
(SKP), which indicate regions of alternating polarity on CE-charged surfaces
with lengthscales ranging from nanometers to centimeters [19, 20, 21, 10, 11,
13, 22]. At the same time, there is ample room for interpretation. For
example, converting KPFM and SKP voltage data to charge data is non-trivial
[23], and experiments to date have not been able to connect the fluctuating
features seen at smaller local scales to the net charge transferred at the
global scale. Recent results by Sobolev et al. suggest that much of the mosaic
patterns observed could be explained by discharge that occurs when samples are
separated after CE [24], and experiments that look at charge transfer from
surface-exploring contacts on granular particles show features reminiscent of
both global and local mechanisms [25].

In this paper, we build a bridge between the local models of same-material CE
and the global models of different material CE by extending mosaic models to
account for asymmetries in the densities of donors/acceptors. For instance, we
can illustrate this by considering the grids A and C in Fig. 1, where C has a
much higher density of donors than A. We reveal a continuous transition from a
‘same-material’ regime where charge exchange is stochastic to a ‘different
material’ regime where it is deterministic. We validate our analytical
treatment of the problem with numerical simulations that mimic the essential
ingredients of the hypothesized charge transfer process. We illustrate the
value of these results by showing they can qualitatively explain real world
situations, e.g. our own experiments with surface-exploring granular contacts
[25].

Going further, we also show that introducing donor/acceptor asymmetries in
mosaic models allows us to make sense of older experimental results regarding
polarity reversal during sliding [26]. Sliding contacts [6, 26, 27, 28] and,
to some extent, rolling contacts [29, 30, 31] are situations where the contact
itself introduces an asymmetry. Contact is spread over a much larger area on
the ‘rubbed’ sample than on the ‘rubber’, e.g. grids A and D in Fig. 1. Using
our generalization of mosaic models, we show that charging can be non-
monotonous with sliding distance, granted the initial donor/acceptor
probabilities differ. In some cases, this can lead the net charge of the
‘rubber’ to flip sign. Both non-monotonous charging and sign flips have been
observed experimentally, but were previously attributed to contaminants on the
surface [26, 28]. The value of our results here, therefore, is to offer an
alternative explanation to this phenomenon which naturally emerges from
donor/acceptor asymmetries.

##  II Accounting for donor/acceptor asymmetry

Outcome | 1\. Strict probabilities | 2\. Relaxed probabilities | 3\. Approximate relaxed probabilities
---|---|---|---
transfer A$\rightarrow$B | $p_{\mathrm{A}}(1-p_{\mathrm{B}})\,\alpha$ | $p_{\mathrm{A}}(1-p_{\mathrm{B}})\,\alpha\left(1-p_{\mathrm{B}}(1-p_{\mathrm{A}})\,\alpha\right)$ | $\sim p_{\mathrm{A}}(1-p_{\mathrm{B}})\,\alpha$
transfer A$\leftarrow$B | $p_{\mathrm{B}}(1-p_{\mathrm{A}})\,\alpha$ | $p_{\mathrm{B}}(1-p_{\mathrm{A}})\,\alpha\left(1-p_{\mathrm{A}}(1-p_{\mathrm{B}})\,\alpha\right)$ | $\sim p_{\mathrm{B}}(1-p_{\mathrm{A}})\,\alpha$
transfer A$\rightleftarrows$B | 0 | $p_{\mathrm{A}}p_{\mathrm{B}}(1-p_{\mathrm{B}})(1-p_{\mathrm{B}})\,\alpha^{2}$ | $\sim 0$
no transfer | $1-(p_{\mathrm{A}}+p_{\mathrm{B}}-2p_{\mathrm{A}}p_{\mathrm{B}})\,\alpha$ | $\left(1-p_{\mathrm{A}}(1-p_{\mathrm{B}})\right)\left(1-p_{\mathrm{B}}(1-p_{\mathrm{A}})\right)\alpha^{2}$ | $\sim 1-(p_{\mathrm{A}}+p_{\mathrm{B}}-2p_{\mathrm{A}}p_{\mathrm{B}})\,\alpha$
Table 1: Transfer probabilities can be considered in three different ways.
Strictly, there are three possible outcomes for each contacting pair of sites:
A$\rightarrow$B transfer, A$\leftarrow$B transfer, or no transfer. The
probabilities of each outcome are presented in column 1. Column 2 presents the
relaxed probabilities, where we assume that A$\rightarrow$B and A$\leftarrow$B
transfers are independent. This means that the probability of simultaneous
transfer in both directions is nonzero. However, as shown in column 3, if we
ignore the $\alpha^{2}$ terms, we recover the strict probabilities.

We begin by considering donor and acceptor sites that are randomly assigned on
a square lattice on each surface and defined by three lengthscales. The first
one is the side length of the whole surface $L$. The second one is the size of
a microscopic site $l_{0}$, which can donate/accept one elementary charge $e$.
Finally, a third, intermediate lengthscale must be included to obtain
realistic values for charge density [18]. It is the correlation length between
sites $l$, which corresponds to the typical size of features on the mosaic.
Once a donor site touches an acceptor site, there is a probability $\alpha$
that the transfer occurs. We will routinely express $L$ and $l$ in units of
$l_{0}$, and charge in units of $e$. Before any contacts, the donor
probability on a given grid M is $p_{\mathrm{M}}^{0}$, and the corresponding
acceptor probability is $s_{\mathrm{M}}^{0}=1-p_{\mathrm{M}}^{0}$.

Let us first consider a single contact between two surfaces, meaning that
$p_{\mathrm{M}}=p_{\mathrm{M}}^{0}$. Denoting the surfaces $A$ and $B$, we
write the donor probabilities $p_{\mathrm{A}}$ and $p_{\mathrm{B}}$, and
acceptor probabilities $s_{\mathrm{A}}=1-p_{\mathrm{A}}$ and
$s_{\mathrm{B}}=1-p_{\mathrm{B}}$. We first look at the case where $l=l_{0}$.
Focusing on an individual pair of opposing sites, three situations can occur
during contact. First, there can be transfer from A to B, which occurs with
probability $p_{\mathrm{A}}(1-p_{\mathrm{B}})\,\alpha$. Second, there can be
transfer from B to A, which occurs with probability
$p_{\mathrm{B}}(1-p_{\mathrm{A}})\,\alpha$. Finally, there can be no transfer
at all, which has probability
$1-(p_{\mathrm{A}}+p_{\mathrm{B}}-2p_{\mathrm{A}}p_{\mathrm{B}})\,\alpha$.
These values are summarized in the first column of Table 1. For large
surfaces, i.e. when $L\gg l_{0}$, the probability of the net charge transfer
from A$\rightarrow$B is Gaussian with mean
$eL^{2}/l_{0}^{2}\,p_{\mathrm{A}}\left(1-p_{\mathrm{B}}\right)\alpha$ and
width
$eL/l_{0}\sqrt{p_{\mathrm{A}}(1-p_{\mathrm{B}})\,\alpha\,(1-p_{\mathrm{A}}(1-p_{\mathrm{B}})\,\alpha)}$.
The expression for A$\leftarrow$B transfer is the same with $p_{\mathrm{A}}$
and $p_{\mathrm{B}}$ permuted.

Determining the distribution for the total charge exchange is more complicated
as it involves the difference of A$\rightarrow$B and A$\leftarrow$B, and at a
particular site transfer both ways cannot technically happen simultaneously. A
helpful simplification comes if we approximate the A$\rightarrow$B and
A$\leftarrow$B events as independent. The probabilities for this scenario are
shown in column 2 of Table 1, where there is now a possibility of simultaneous
transfer A$\rightleftarrows$B. Considering experiments suggest $\alpha\ll 1$
[9], and knowing that $p_{\mathrm{A}}\leq 1$ and $p_{\mathrm{B}}\leq 1$, we
make the assumption that $\alpha^{2}$ terms can be ignored. As shown in column
3 of Table 1, the relaxed case with independent events is then approximately
equal to the strict case of column 1.

Now treating the distributions of A$\rightarrow$B and A$\leftarrow$B transfer
as independent Gaussians, the distribution for the total charge exchanged
after one contact $Q_{1}$ is also a Gaussian, with mean

$\overline{Q_{1}}_{\,l=l_{0}}=e\alpha\frac{L^{2}}{l_{0}^{2}}\left(p_{\mathrm{A}}-p_{\mathrm{B}}\right)$
(1)

and standard deviation

$\sigma_{l=l_{0}}\approx
e\frac{L}{l_{0}}\sqrt{\alpha\left(p_{\mathrm{A}}+p_{\mathrm{B}}-2p_{\mathrm{A}}p_{\mathrm{B}}\right)}.$
(2)

Again, these expressions assume we can ignore $\alpha^{2}$ terms in $\sigma$
and take the convention that the positive direction for charge exchange is
A$\rightarrow$B.

The case where $l>l_{0}$ is treated similarly, but considering a rescaled
matrix containing large patches of $n=(l/l_{0})^{2}$ sites. The distance $l$
can be understood as the typical size of a patch, the average distance at
which sites become uncorrelated or, equivalently, the first zero of the
spatial correlation function [18]. To calculate the distribution, we will not
immediately consider charge transfer, and will first consider instead the
number of patches of $n$ donors on A facing patches of $n$ acceptors on B,
which we denote $N_{\rightarrow}$. There is probability
$p_{\mathrm{A}}(1-p_{\mathrm{B}})$ that a donor patch on A would face an
acceptor patch on B. Assuming $L\gg l$, the values of $N_{\rightarrow}$ follow
a normal distribution with mean
$\mu_{\rightarrow}=L^{2}/l^{2}p_{\mathrm{A}}(1-p_{\mathrm{B}})$ and width
$\sigma_{\rightarrow}=L/l\sqrt{p_{\mathrm{A}}(1-p_{\mathrm{B}})(1-p_{\mathrm{A}}(1-p_{\mathrm{B}}))}$.
Denoting $N_{\leftarrow}$ the number donor patches on B facing acceptor
patches on A, we find the same expressions for $\mu_{\leftarrow}$ and
$\sigma_{\leftarrow}$ with A and B permuted. Treating the two distributions as
independent, the distribution for $N_{\rightarrow}-N_{\leftarrow}$ is also a
Gaussian with mean $\mu_{\rightarrow}-\mu_{\leftarrow}$ and standard deviation
$\sqrt{\sigma_{\rightarrow}^{2}+\sigma_{\leftarrow}^{2}}$.

We note that each patch can transfer or receive a maximum of $n$ elementary
charges. If $n\gg 1$, then the actual number of charges transferred by a donor
patch is close to $\alpha n$, with fluctuations
$\sqrt{n\alpha(1-\alpha)}\ll\alpha n$. This means that $\alpha$ is treated
slightly differently here, determining the average transfer rate between
patches. We find the same expression for $\overline{Q_{1}}_{\,l>l_{0}}$ as the
one from Eq. (1), whereas the standard deviation can be written

$\begin{split}\sigma_{l>l_{0}}\approx
e\frac{Ll}{l_{0}^{2}}\alpha\sqrt{\left(p_{\mathrm{A}}+p_{\mathrm{B}}-2p_{\mathrm{A}}p_{\mathrm{B}}\right)}\end{split}$
(3)

assuming that $p_{\mathrm{A}}\approx p_{\mathrm{B}}$. In this case, the
differences with Eq. (2) lay in the exponent $1/2$ of $\alpha$, and a factor
$l/l_{0}$.

Finally, we express the donor probabilities in terms of their difference
$\delta=p_{\mathrm{B}}-p_{\mathrm{A}}$ and average
$p=(p_{\mathrm{A}}+p_{\mathrm{B}})/2$. It follows that the charge exchanged
after one contact follows a Gaussian distribution with mean

$\overline{Q_{1}}=-e\alpha\frac{L^{2}}{l_{0}^{2}}\delta\equiv\mu$ (4)

and width

$\sigma\approx
e\frac{Ll}{l_{0}^{2}}\alpha^{m}\sqrt{2\left(p\left(1-p\right)+\frac{\delta^{2}}{4}\right)}$
(5)

where $m=0.5$ for $l=l_{0}$ and $1$ for $l>l_{0}$. While $\delta$ appears in
$\sigma$, its most noticeable effect compared to the known case of
$p_{\mathrm{A}}=p_{\mathrm{B}}$ is that the distribution is no longer centered
around zero. If $\delta=0$, charge transfer is only powered by the statistical
differences between two surfaces, which lead to a nonzero charge exchange
magnitude $\overline{|Q_{1}|}=\sigma\sqrt{2/\pi}$. In the general case
however, charge exchange magnitude follows a folded normal distribution which
can be written

$\overline{\left|Q_{1}\right|}=\sigma\sqrt{\frac{2}{\pi}}\exp{\frac{-\mu^{2}}{2\sigma^{2}}}+\mu\operatorname{erf}{\frac{\mu}{\sqrt{2\sigma^{2}}}}$
(6)

where $\operatorname{erf}{(x)}$ denotes the Gauss error function. One can
easily show that Eqs. (4) through (6) recover the results from Ref. [18] when
$\delta=0$.

Figure 2: Model for donor probability asymmetry. (a) The charge exchanged
during a single contact $Q_{1}$ follows a Gaussian distribution, no longer
centered around $0$ when $\delta\neq 0$. Median and widths are given by Eqs.
(4) and (5), respectively. Statistics are shown for 5000 pairs of surfaces
with $p=0.5$, $\alpha=0.1$, $l=5\,l_{0}$ and $L=1000\,l_{0}$. (b) Charge
exchange magnitude $|Q_{1}|$ follows a folded normal distribution. (c) Average
magnitude $\overline{|Q_{1}|}$ is described by Eq. (6). The intercept
$\sigma\sqrt{2/\pi}$ corresponds to the same-material case. The ratio
$\mu/\sigma$, proportional to $l/l_{0}$, determines whether charging is
powered by fluctuations or not. (d) Charging can be linear or quadratic in
$L$, depending on whether charging is average-driven or variance-driven. The
same system can go from linear to quadratic as $L$ increases.

To validate this analysis, we perform numerical simulations, where surfaces
are generated and ‘contacts’ performed using the method described in Ref.
[18]. Features of size $l$ can be included by interpolating a coarse
$L/l\times L/l$ matrix filled with random numbers between 0 and 1, then
calculating the threshold value that leads to a binary field with the desired
donor probability $p$ [18]. Figure 2(a) shows the distribution of charge
transferred after contact between two square surfaces of side lengths
$L=1000\,l_{0}$, with different colors corresponding to different values of
the donor density difference, $\delta$. Each distribution is determined from
simulating 2000 pairs of such surfaces. As can be seen, introducing $\delta>0$
shifts the Gaussian profile, which is no longer centered around zero. The
solid lines correspond to normal distributions whose means and widths are
given by Eqs. (4) and (5). Figure 2(b) shows the magnitude of charge exchange
$|Q_{1}|$ for the same set of data. For small values of $\delta$, the
distributions are folded, their average value $\overline{|Q_{1}|}$ being given
by Eq. (6). As $\delta$ increases and the distributions move away from zero,
they become essentially Gaussian.

The folded Gaussian and Gaussian distributions can be understood as two
regimes, one where the standard deviation, $\sigma$, most significantly
contributes to $\overline{|Q_{1}|}$, and the other where the average $\mu$,
and therefore $\delta$, is primarily responsible. These variance-driven and
average-driven regimes, at their extremes, can be understood as same-material
CE with macroscopically identical samples, and CE between two different
materials, respectively. One possible way to identify where the regime change
occurs is to compare the two terms of Eq. (6). Posing $\chi=\mu/\sigma$, the
terms are equal if

$\sqrt{\frac{2}{\pi}}\exp{\frac{-\chi}{2}}=\chi\operatorname{erf}{\frac{\chi}{\sqrt{2}}}$
(7)

which can be solved numerically to find $\chi\approx\pm 0.8769$. Note that
$\chi$ is proportional to $l/l_{0}$ and inversely proportional to $L/l_{0}$,
meaning that systems where correlations lengths are small compared to their
overall size tend to be more average-dominated. Figure 2(c) compares Eq. (6)
with simulations. Each point corresponds to 1000 first contacts between two
$1000\times 1000$ grids, with $\alpha=1$ and $p=0.5$, averaged over 30 pairs
of grids. Five values of $l/l_{0}$ are shown, namely 1, 20, 50, 70 and 100.
The intercepts in $\delta=0$ correspond to the known case of
$\overline{|Q_{1}|}=\sigma\sqrt{2/\pi}$. For large $\delta$, the effect of
correlation length $l/l_{0}$ disappears as charge exchange is purely driven by
average differences between surfaces.

We can point to experiments whose outcomes are consistent with the scaling
laws obtained for the variance-dominated and average-dominated regimes. For
instance, we can note that $\mu\sim L^{2}$ and $\sigma\sim L$. In experiments
with polymers, Apodaca _et al._ measured $\overline{|Q_{1}|}$ as a function of
$L$ [9]. When contacting Polydimethylsiloxane (PDMS) with either Polyvinyl
chloride (PVC), stainless steel or copper, they found that
$\overline{|Q_{1}|}\sim L^{\eta}$ with values of $\eta$ between 1.6 and 1.9,
slightly below 2. With PDMS-PDMS contacts, however, they observed that $\eta$
dropped to approximately 1. Equation (6) predicts such linear and quadratic
behaviors of $\overline{|Q_{1}|}$ with $L$, with a continuous transition
between them when the contributions of $\mu$ and $\sigma$ are comparable. The
fact that $\mu\sim L^{2}$ and $\sigma\sim L$ also means that the same
materials could behave differently at different sizes, i.e. charging could be
variance-dominated for smaller samples, but average-dominated for larger ones.
Figure 2(d) illustrates this, comparing simulations with Eqs. (6) and (7).

Figure 3: Comparison with experiments. (a) Single collision statistics
measured experimentally, taken from Ref. [25], for three pairs of samples. The
dashed line is a Gaussian fit. (b) Simulations based on donor/acceptor
asymmetry proposed in this paper, performed with parameters inferred from the
fits using Eqs. (4) and (5). The model is able to explain the experimental
data, with donor/acceptor donor density differences $\delta$ on the order of
$<1$%.

We can also use this framework to make sense of recent experiments [25] where
single-contact charge-transfer distributions were measured for same-material
samples. Figure 3 shows three such distributions taken from Ref. [25],
measured using acoustic levitation to perform collisions between a sphere and
a plate, both made of fused silica. These distributions are not always
centered around zero, meaning that there are global differences between
samples which can potentially dominate over local fluctuations. Incorporating
a difference in donor probability $\delta\neq 0$ in the model makes off-center
distributions possible. We know from Ref. [25] that $L\approx 20$ µm. Using
Eqs. (4) and (5), we can extract $\mu$ and $\sigma$ from the distributions. We
then translate those values in terms of $l/l_{0}$ and $\delta$ by making
reasonable assumptions for the remaining parameters. Inspired by measurements
from Refs. [9] and [11], we pose $\alpha=0.1$, $p=0.5$ and $l_{0}=1$ nm. From
Eq. (5), we find that

$l=\frac{\sigma l_{0}^{2}}{\alpha
L\sqrt{2\left(p(1-p)+\frac{\delta^{2}}{4}\right)}},$ (8)

and from equ. (4), we get

$\delta=\frac{\mu l_{0}^{2}}{e\alpha L^{2}}.$ (9)

Assuming by convention that the distributions were measured for surface A and
that the charge carriers are positive, we obtain the values for $l/l_{0}$ and
$\delta$ shown in Fig. 3(b). We can then use these parameters as inputs for
numerical simulations which we contrast against the experiments, performing
the same number of measurements $N$.

Thanks to the introduction of $\delta$, the simulations are able to reproduce
the experimental results well, as illustrated by plotting the Gaussian fits
from the experiment over the simulated data. The width of the distributions is
relatively constant, leading to little variation of $l/l_{0}\approx 73\pm 15$,
which is not too far from the size of the features observed in Ref. [11]. We
also note that $|\delta|=0.003\pm 0.002$, which is small compared to $p$. This
means that very slight, but nonetheless global differences in donor/acceptor
densities can significantly affect the outcome of experiments, in a way that
is distinct from purely local fluctuations. The simulations from Fig. 3(b)
show that the generalized mosaic framework can produce realistic results,
which could for example be incorporated into simulations of granular media.
Nonetheless, we had to resort to assumptions for parameters that are difficult
to access experimentally. Dedicated measurements could solve this issue. For
instance, $\alpha$ can be obtained by repeating contacts if the location on
the samples is fixed. Atomic force microscopy could be used to measure $l$ and
perhaps estimate $p$, assuming such features can be observed on fused silica
samples. The regime change from $|Q|\sim L$ to $|Q|\sim L^{2}$ from Fig. 2(d),
which depends on the ratio $\mu/\sigma$, could, if observed, be contrasted
with independent measurements of $\mu$ and $\sigma$. Regime change or not,
reproducing the $|Q|\sim L$ scaling from Ref. [9] would be highly valuable.
Considering that very small ($<$1%) differences in donor/acceptor densities
are sufficient to leave the variance-driven regime, it is possible that this
regime would only exist in carefully controlled conditions with strictly
identical samples. It is, however, difficult to estimate a plausible range for
$\delta$ without knowing more about the nature of the donors/acceptors.

##  III Asymmetry during sliding

We now show how thinking about donor/acceptor density differences can be used
to treat another issue in CE, namely when a smaller area on one surface is
slid over a larger area on another. We consider two samples of different size,
divided into microscopic square sites of side length $l_{0}$. Sample A of size
$w\times L$ is slid along a larger sample B of size $W\times L$, with $W>w$.
We denote $d$ the distance over which A has slid on B. If we consider that A
and B are fully in contact at the start and end of the sliding, we have $d\leq
W-w$. The first contact occurs at sliding distance $d=0$. Subsequently, every
time $d$ increases by $l_{0}$, a new contact occurs. Figure 4 illustrates
sliding using a simulation where $w=L=50\,l_{0}$ and $W=250\,l_{0}$.

Figure 4: Asymmetry during sliding. (a) A square grid A with $50\times 50$
sites slides on a rectangular grid B with $50\times 250$ sites. Parameters are
$\alpha=0.1$, $p_{\mathrm{A}}^{0}=0.45$, $p_{\mathrm{B}}^{0}=0.25$, and
$l/l_{0}=4$. Every time sliding distance $d$ is increased by one site length
$l_{0}$, a new contact occurs. After one contact, some donor and acceptor
sites become charged and no longer contribute. (b) As grid A slides along B,
it is exposed to new regions with unused sites. (c) Around $d=80\,l_{0}$, most
donors on A have been used. Only the acceptor sites contribute, causing
charging direction to reverse. (d) At $d=200\,l_{0}$, most acceptors on A have
been used, and charging stops. (e) Densities $p_{\mathrm{A}}$ and
$s_{\mathrm{A}}$ both decrease monotonously with $d$, but at different rates.
On the other hand, the section of B facing A is continually renewed. At
$d=200\,l_{0}$, we see that $p_{\mathrm{A}}\approx 0$ and
$s_{\mathrm{A}}\approx 0$, but $s_{\mathrm{B}}\approx s_{\mathrm{B}}^{0}$ and
$p_{\mathrm{B}}\approx p_{\mathrm{B}}^{0}$. Vertical lines indicate the values
of $d$ which correspond to panels (a) through (d).

During sliding, certain donors and acceptors will contribute and become
inactive, leading donor and acceptor densities $p_{\mathrm{A}}$ and
$s_{\mathrm{A}}$ to decrease with $d$. On B, however, the sliding continually
brings new, unused sites. Note that we only consider the section of B which
faces A at any given time when calculating $p_{\mathrm{B}}$ and
$s_{\mathrm{B}}$. This means that, even in a situation where the
donor/acceptor densities on A and B are initially equal, they will eventually
become unequal—hence recreating the donor/acceptor asymmetry discussed in the
previous section. In general, we have $p_{\mathrm{M}}(d,x)$ and
$s_{\mathrm{M}}(d,x)$, where $x$ is the horizontal position on M along the
sliding direction. We denote $p_{\mathrm{M}}^{0}$ and
$s_{\mathrm{M}}^{0}=1-p_{\mathrm{M}}^{0}$ the initial donor and acceptor
densities. Note that for $d>0$, we can no longer assume that
$s_{\mathrm{M}}(d,x)=1-p_{\mathrm{M}}(d,x)$. In Fig. 4(a), $d=0$ and only one
contact has occurred. Most sites are still in their unused state (shown in
white for donors and gray for acceptors).

Figure 4(e) shows the evolution with $d$ of donor/acceptor densities on A and
B. First, one can see that $p_{\mathrm{A}}$ and $s_{\mathrm{A}}$ decrease
monotonously, but at different rates. This is also visible in Fig. 4(c), where
all donor sites on A have already transferred their charge (in dark purple),
while there are still many acceptor sites available (in white). This can cause
a reversal of the charging direction. Sample A is not uniform at this stage,
with most of the remaining acceptor sites (in white) located toward the
trailing edge, i.e. toward smaller $x$. The evolution of donor/acceptor
densities on B is qualitatively different from A, as can be seen in Fig. 4(e).
After an initial decrease of $p_{\mathrm{B}}$ and $s_{\mathrm{B}}$, they
slowly come back to approximately their initial values $p_{\mathrm{B}}^{0}$
and $s_{\mathrm{B}}^{0}$. This is because the part of B contributing to charge
exchange is constantly renewed. Looking at B on Fig. 4(d), we can see that the
acceptors which have acquired charge (light yellow) are concentrated toward
the beginning of the track, below $x\approx 50\,l_{0}$, while the used up
acceptors (dark purple) are spread over a much larger portion of the track, up
to $x\approx 200\,l_{0}$.

In normal, forward facing contacts, the number of sites contributing to charge
exchange at every contact is the number of donors facing acceptors (and vice-
versa) multiplied by the transfer probability $\alpha$. Writing the number of
contacts $n$ as a continuous variable, one can show than the number of charges
being exchanged at every step decreases like $\exp\left(-\alpha n\right)$ [9,
18]. However, when sliding, the number of contacts $n$ experienced by a given
site is now a function of the location on A or B, as discussed above. Every
time A moves by $l/l_{0}$, the first column of sites in the direction of the
sliding faces ‘new’ sites on B. As a consequence, sites on the leading edge of
A will initially experience more charge exchange than sites on the trailing
edge.

The whole surface of A contributes to every contact, as every site on A
experiences $n_{\mathrm{A}}=d/l_{0}$ contacts. On B, however, the number of
contacts is a function of the position $x$ along the sliding direction, which
can be written as the convolution product
$n_{\mathrm{B}}\left(x\right)=\int_{0}^{x}H(w-\tau)H(d-(x-\tau))\,\mathrm{d}\tau$
where $H$ is the Heaviside step function. This causes, in turn,
$p_{\mathrm{B}}$ and $s_{\mathrm{B}}$ to depend on $x$, as well as
$p_{\mathrm{A}}$ and $s_{\mathrm{A}}$ through the contacts with B, which makes
the analytical calculation of $Q$ significantly more difficult in the general
case.

The case $w=l_{0}$, however, is comparatively straightforward. Sample A
becomes as thin as possible and the dependency on $x$ disappears. In this
case, $p_{\mathrm{B}}$ and $s_{\mathrm{B}}$ remain constant, and every
additional contact causes $p_{\mathrm{A}}$ to decrease by a factor
$(1-s_{\mathrm{B}}^{0})\,\alpha$, and $s_{\mathrm{A}}$ by a factor
$(1-p_{\mathrm{B}}^{0})\,\alpha$. To illustrate this, we write the evolution
of $p_{\mathrm{A}}$ for $d={0,1,2\dots n}$, keeping in mind that
$s_{\mathrm{B}}=s_{\mathrm{B}}^{0}$ as the surface is continually renewed. At
every contact, we subtract the number of transfers going from A to B. This
leads to

$\begin{split}p_{\mathrm{A}}(0)&=p_{\mathrm{A}}^{0}\\\
p_{\mathrm{A}}(1)&=p_{\mathrm{A}}(0)-p_{\mathrm{A}}(0)\,s_{\mathrm{B}}^{0}\,\alpha=p_{\mathrm{A}}^{0}\left(1-s_{\mathrm{B}}^{0}\,\alpha\right)\\\
p_{\mathrm{A}}(2)&=p_{\mathrm{A}}(1)-p_{\mathrm{A}}(1)\,s_{\mathrm{B}}^{0}\,\alpha=p_{\mathrm{A}}^{0}\left(1-2\,s_{\mathrm{B}}^{0}\,\alpha+\left(s_{\mathrm{B}}^{0}\,\alpha\right)^{2}\right)\\\
&\leavevmode\nobreak\ \leavevmode\nobreak\ \vdots\\\
p_{\mathrm{A}}(n)&=p_{\mathrm{A}}^{0}\left(1-n\,s_{\mathrm{B}}^{0}\,\alpha\right)+\mathcal{O}(\alpha^{2}).\end{split}$

The reasoning is the same for $s_{\mathrm{A}}$. We ignore the $\alpha^{2}$
terms and consider $d$ as a continuous variable. When $d$ increases by a small
increment $\delta d$, $p_{\mathrm{A}}$ increases by

$\delta
p_{\mathrm{A}}=-p_{\mathrm{A}}\left(1-s_{\mathrm{B}}^{0}\right)\alpha\,\delta
d$

which we integrate to obtain

$p_{\mathrm{A}}\left(d\right)=p_{\mathrm{A}}^{0}\exp{\left(-\alpha\left(1-p_{\mathrm{B}}^{0}\right)\frac{d}{l_{0}}\right)}=p_{\mathrm{A}}^{0}\exp{-\frac{d}{\lambda_{p}}}.$
(10)

Similarly, we find

$s_{\mathrm{A}}\left(d\right)=\left(1-p_{\mathrm{A}}^{0}\right)\exp{\left(-\alpha\,p_{\mathrm{B}}^{0}\frac{d}{l_{0}}\right)}=s_{\mathrm{A}}^{0}\exp{-\frac{d}{\lambda_{s}}}$
(11)

where $\lambda_{p}=l_{0}/\alpha\,(1-p_{\mathrm{B}}^{0})$ and
$\lambda_{s}=l_{0}/\alpha\,p_{\mathrm{B}}^{0}$ are the characteristic sliding
lengths over which donors and acceptors on A get depleted.

Figure 5 compares Eq. (10) with numerical simulations, with
$L=W=10000\,l_{0}$,
$p_{\mathrm{A}}^{0}=s_{\mathrm{A}}^{0}=p_{\mathrm{B}}^{0}=p_{\mathrm{B}}^{0}=0.5$
and $l=l_{0}$. The width $w$ is increased from $l_{0}$ to $L$. One can see
that when $w>l_{0}$ and $d>l_{0}$, we depart from the exponential behavior.
The decrease of $p_{\mathrm{A}}$ is slowed down by the presence of used up
sites on B. In other words, we locally have
$p_{\mathrm{B}}<p_{\mathrm{B}}^{0}$ and $s_{\mathrm{B}}<s_{\mathrm{B}}^{0}$,
leading $\lambda_{p}$ and $\lambda_{s}$ to effectively increase. As can be
seen on Fig. 5(b), this causes the exponential to plateau until about $d\geq
w$, where the remaining active sites on A become sufficiently sparse so that
we once more have $p_{\mathrm{B}}\approx p_{\mathrm{B}}^{0}$ and
$s_{\mathrm{B}}\approx s_{\mathrm{B}}^{0}$. Figure 5(c) compares Eq. (10) for
a wider range of parameters. We see a good agreement for any combination of
$\delta$ and $p$, granted that $\alpha\ll 1$ and $l\ll L$. The curves are
plotted until $p_{\mathrm{A}}$ has decreased by 99%, as fluctuations become
large when only a few sites remain.

Figure 5: Simple sliding model. (a) The decrease of $p_{\mathrm{A}}$ follows
Eq. (10) (dashed line) only if grid A is infinitely narrow, i.e. when
$w=l_{0}$. We start with a $1\times 10000$ grid and increase the width
progressively up to $10000\times 10000$. (b) On a semi-log scale, one can see
that Eq. (10) works for sliding distances $d\ll L$, as A is still
approximately uniform, and for $d\geq w$. (c) We verify Eq. (10) for multiple
values of $\alpha$, $p$ and $\delta$, but keeping $L=10000\,l_{0}$, $w=l_{0}$
and $l=l_{0}$.

From Eqs. (10) and (11), we see that if $p_{\mathrm{B}}=0.5$, then donors and
acceptors on A are depleted at exactly the same rate. If $p_{\mathrm{B}}\neq
0.5$, there can exist a sliding distance $d^{*}$ such that
$p_{\mathrm{A}}\left(d^{*}\right)=s_{\mathrm{A}}\left(d^{*}\right)$. From Eqs.
(10) and (11), we get

$d^{*}=\frac{\alpha\,l_{0}}{2p_{\mathrm{B}}-1}\ln{\left(\frac{1-p_{\mathrm{A}}^{0}}{p_{\mathrm{A}}^{0}}\right)}.$
(12)

This means that past a certain sliding distance, the average direction of
charge exchange can reverse. Note that $d^{*}$ is positive only if either
$p_{\mathrm{A}}>0.5$ and $p_{\mathrm{B}}<0.5$, or $p_{\mathrm{A}}<0.5$ and
$p_{\mathrm{B}}>0.5$. Such reversals were observed experimentally when sliding
various materials on glass [26] or with same-material samples [28], but were
generally attributed to the presence of contamination on the samples’
surfaces. However, Eq. (12) shows that reversals can naturally arise as a
consequence of different donor/acceptor probabilities on the two samples.

Figure 6: Reversals of charging direction and sign during sliding. We measure
the charge of the square grid A as a function of sliding distance $d/l_{0}$,
with $\alpha=0.1$, $L=500\,l_{0}$, and either $l=5\,l_{0}$ (solid red lines),
$l=5\,l_{0}$ (dashed orange lines) or $l=20\,l_{0}$ (dotted purple lines). We
show six trajectories for each set of parameters to illustrate inherent
variability. (a, b) When $p_{\mathrm{A}}^{0}=p_{\mathrm{B}}^{0}=p$, we observe
that (a) $Q_{\mathrm{A}}=0$ on average when $p=0.5$, (b) $Q_{\mathrm{A}}$
increases (resp. decreases) monotonously for $p>0.5$ (resp. $p<0.5$). (c-f)
When $p_{\mathrm{A}}^{0}\neq p_{\mathrm{B}}^{0}$, a much wider range of charge
behaviors can occur. This includes (c) monotonous charging, (d) non-monotonous
charging with no zero-crossing, (e) zero-crossings or (f) combinations of
those depending on random fluctuations, which increase with $l/l_{0}$.

By running simulations, we can show that the necessary condition that
$p_{\mathrm{A}}^{0}\lessgtr 0.5$ and $p_{\mathrm{B}}^{0}\gtrless 0.5$ can be
somewhat relaxed for wider samples $w>l_{0}$. Indeed, as $p_{\mathrm{A}}$ and
$p_{\mathrm{B}}$ become functions of the position on A and B, the condition
for the reversal of charging direction can be met locally. In Fig. 6, we
measure the charge of A during sliding for various values of
$p_{\mathrm{A}}^{0}$ and $p_{\mathrm{B}}^{0}$. The different behaviors
observed in Refs. [26, 28] are reproduced, including reversal of charging
direction (non-monotonous charging) and reversal of charge (zero-crossings).
In general, non-monotonous charging only occurs if $p_{\mathrm{A}}^{0}\neq
p_{\mathrm{B}}^{0}$. For larger values of lengthscale $l/l_{0}$, the larger
fluctuations between surfaces can also cause that condition to be met.

##  IV Discussion and conclusions

Looking at the literature, it is still unclear whether the patterns of
positive and negative voltage sometimes observed after contact are a sign of
actual charge donor/acceptor mosaics [11], a consequence of discharge during
separation [24], or something else entirely. Regardless, mosaic models
naturally follow from the assumption that charging is not a perfectly uniform
process, and provide a useful conceptual tool for analytical and numerical
studies. By accounting for global differences between samples, we extend this
framework from strictly identical to different-material samples, and anything
in between. Our model shows that whether charging is driven by global
differences or local fluctuations depends not only on material properties, but
on scale. Local fluctuations tend to matter more for smaller samples, and the
transition might be within reach of experimental observations.

Accounting for differences in donor/acceptor densities allows us to not only
reproduce same-material CE measurements [25], but also make sense of contact
sliding experiments [26, 28]. Combined with the asymmetric geometry due to the
sliding, donor/acceptor asymmetries are sufficient to explain the sign flips
that were observed experimentally, as an alternative explanation to impurities
or contamination.

One might argue that the flexibility of mosaic models is also their weakness.
Some parameters are difficult to access experimentally, which so far has
prevented the field from drawing conclusions about the underlying mechanism.
However, we suggest experiments based on our results which might allow to
build a fuller picture of CE. Repeatedly contacting two samples, resetting
between contacts by discharging, provides the mean $\mu$ and standard
deviation $\sigma$ from Eqs. (4) and (5), as established here with the data
from Ref. [25]. If multiple contacts are performed on the same location, one
can determine transfer probability $\alpha$. Varying contact area over a wide
enough range would allow to verify Eq. (6) by providing an independent
measurement of the ratio $\mu/\sigma$. Deviations from the law could reveal
information about $l$ and $l_{0}$, or at the very least provide an upper
bound. Finally, sliding experiments provide an underutilized tool to test CE
theories, as the charging profiles depend on the donor/acceptor probabilities
of the two surfaces in non-trivial ways. Performing same-material and
different-material sliding experiments in a variety of conditions, including
changing the aspect ratio of the contact area, and noting if and where the
sign flips occur could be a new way to estimate donor/acceptor probabilities
to probe the underlying principles governing tribocharging.

## Acknowledgements

This project has received funding from the European Research Council Grant
Agreement No. 949120 and from the European Union’s Horizon 2020 research and
innovation program under the Marie Skłodowska-Curie Grant Agreement No.
754411.

## References

  * Rescaglio  _et al._ [2017] A. Rescaglio, J. Schockmel, N. Vandewalle, and G. Lumay, Combined effect of moisture and electrostatic charges on powder flow, EPJ Web Conf. 140, 13009 (2017).
  * Cimarelli  _et al._ [2022] C. Cimarelli, S. Behnke, K. Genareau, J. M. Harper, and A. R. Van Eaton, Volcanic electrification: recent advances and future perspectives, Bull Volcanol 84, 78 (2022).
  * Zhang and Zhou [2020] H. Zhang and Y.-H. Zhou, Reconstructing the electrical structure of dust storms from locally observed electric field data, Nat Commun 11, 5072 (2020).
  * Steinpilz  _et al._ [2020] T. Steinpilz, K. Joeris, F. Jungmann, D. Wolf, L. Brendel, J. Teiser, T. Shinbrot, and G. Wurm, Electrical charging overcomes the bouncing barrier in planet formation, Nat. Phys. 16, 225 (2020).
  * Lacks and Shinbrot [2019] D. J. Lacks and T. Shinbrot, Long-standing and unresolved issues in triboelectric charging, Nat. Rev. Chem. 3, 465 (2019).
  * Shaw [1926] P. E. Shaw, The Electrical Charges from Like Solids, Nature 118, 659 (1926).
  * Lowell [1975] J. Lowell, Contact electrification of metals, J. Phys. D: Appl. Phys. 8, 53 (1975).
  * Zhang  _et al._ [2019] X. Zhang, L. Chen, Y. Jiang, W. Lim, and S. Soh, Rationalizing the Triboelectric Series of Polymers, Chem. Mater. 31, 1473 (2019).
  * Apodaca  _et al._ [2010] M. Apodaca, P. Wesson, K. Bishop, M. Ratner, and B. Grzybowski, Contact Electrification between Identical Materials, Angew. Chem. Int. Ed. 49, 946 (2010).
  * Knorr [2011] N. Knorr, Squeezing out hydrated protons: low-frictional-energy triboelectric insulator charging on a microscopic scale, AIP Advances 1, 022119 (2011).
  * Baytekin  _et al._ [2011] H. T. Baytekin, A. Z. Patashinski, M. Branicki, B. Baytekin, S. Soh, and B. A. Grzybowski, The Mosaic of Surface Charge in Contact Electrification, Science 333, 308 (2011).
  * Sow  _et al._ [2012] M. Sow, D. J. Lacks, and R. Mohan Sankaran, Dependence of contact electrification on the magnitude of strain in polymeric materials, J. Appl. Phys. 112, 084909 (2012).
  * Burgo  _et al._ [2012] T. A. L. Burgo, T. R. D. Ducati, K. R. Francisco, K. J. Clinckspoor, F. Galembeck, and S. E. Galembeck, Triboelectricity: Macroscopic Charge Patterns Formed by Self-Arraying Ions on Polymer Surfaces, Langmuir 28, 7407 (2012).
  * Balestrin  _et al._ [2014] L. B. d. S. Balestrin, D. D. Duque, D. S. d. Silva, and F. Galembeck, Triboelectricity in insulating polymers: evidence for a mechanochemical mechanism, Faraday Discuss. 170, 369 (2014).
  * Xie  _et al._ [2013] L. Xie, G. Li, N. Bao, and J. Zhou, Contact electrification by collision of homogenous particles, Journal of Applied Physics 113, 184908 (2013).
  * Yu  _et al._ [2017] H. Yu, L. Mu, and L. Xie, Numerical simulation of particle size effects on contact electrification in granular systems, J. Electrostat. 90, 113 (2017).
  * Harris  _et al._ [2019] I. A. Harris, M. X. Lim, and H. M. Jaeger, Temperature dependence of nylon and PTFE triboelectrification, Phys. Rev. Mater. 3, 085603 (2019).
  * Grosjean  _et al._ [2020] G. Grosjean, S. Wald, J. C. Sobarzo, and S. Waitukaitis, Quantitatively consistent scale-spanning model for same-material tribocharging, Phys. Rev. Materials 4, 082602(R) (2020).
  * Hull [1949] H. H. Hull, A Method for Studying the Distribution and Sign of Static Charges on Solid Materials, Journal of Applied Physics 20, 1157 (1949).
  * Terris  _et al._ [1989] B. D. Terris, J. E. Stern, D. Rugar, and H. J. Mamin, Contact electrification using force microscopy, Phys. Rev. Lett. 63, 2669 (1989).
  * Shinbrot  _et al._ [2008] T. Shinbrot, T. S. Komatsu, and Q. Zhao, Spontaneous tribocharging of similar materials, EPL 83, 24004 (2008).
  * Barnes and Dinsmore [2016] A. M. Barnes and A. D. Dinsmore, Heterogeneity of surface potential in contact electrification under ambient conditions: A comparison of pre- and post-contact states, Journal of Electrostatics 81, 76 (2016).
  * Pertl  _et al._ [2022] F. Pertl, J. C. Sobarzo, L. Shafeek, T. Cramer, and S. Waitukaitis, Quantifying nanoscale charge density features of contact-charged surfaces with an FEM/KPFM-hybrid approach, Phys. Rev. Mater. 6, 125605 (2022).
  * Sobolev  _et al._ [2022] Y. I. Sobolev, W. Adamkiewicz, M. Siek, and B. A. Grzybowski, Charge mosaics on contact-electrified dielectrics result from polarity-inverting discharges, Nat. Phys. 18, 1347 (2022).
  * Grosjean and Waitukaitis [2023] G. Grosjean and S. Waitukaitis, Single-Collision Statistics Reveal a Global Mechanism Driven by Sample History for Contact Electrification in Granular Media, Phys. Rev. Lett. 130, 098202 (2023).
  * Shaw  _et al._ [1928] P. E. Shaw, C. S. Jex, and W. B. Hardy, Tribo-electricity and friction. II.—Glass and solid elements, Proc. R. Soc. A 118, 97 (1928).
  * Lowell and Truscott [1986a] J. Lowell and W. S. Truscott, Triboelectrification of identical insulators. I. An experimental investigation, J. Phys. D: Appl. Phys. 19, 1273 (1986a).
  * Lowell and Truscott [1986b] J. Lowell and W. S. Truscott, Triboelectrification of identical insulators. II. Theory and further experiments, J. Phys. D: Appl. Phys. 19, 1281 (1986b).
  * Grzybowski  _et al._ [2003] B. A. Grzybowski, J. A. Wiles, and G. M. Whitesides, Dynamic Self-Assembly of Rings of Charged Metallic Spheres, Phys. Rev. Lett. 90, 083903 (2003).
  * Wiles  _et al._ [2004] J. A. Wiles, M. Fialkowski, M. R. Radowski, G. M. Whitesides, and B. A. Grzybowski, Effects of Surface Modification and Moisture on the Rates of Charge Transfer between Metals and Organic Materials, J. Phys. Chem. B 108, 20296 (2004).
  * Thomas III  _et al._ [2008] S. W. Thomas III, S. J. Vella, G. K. Kaufman, and G. M. Whitesides, Patterns of Electrostatic Charge and Discharge in Contact Electrification, Angewandte Chemie 120, 6756 (2008).


# A quick algorithm to compute an approximated
power spectral density from an arbitrary Allan deviation

Fabrizio De Marchi<EMAIL_ADDRESS>Department of Mechanical and
Aerospace Engineering, Sapienza University of Rome, Via Eudossiana, 18, 00184
Rome, Italy  Michael K. Plumaris  Department of Mechanical and Aerospace
Engineering, Sapienza University of Rome, Via Eudossiana, 18, 00184 Rome,
Italy  Eric A. Burt  Jet Propulsion Laboratory, California Institute of
Technology, Pasadena, CA, USA  Luciano Iess  Department of Mechanical and
Aerospace Engineering, Sapienza University of Rome, Via Eudossiana, 18, 00184
Rome, Italy

###### Abstract

Complex architectures for wireless communications, digital electronics and
space-based navigation interlink several oscillator-based devices such as
clocks, transponders and synthesizers. Estimators characterizing their
stability are critical for addressing the impact of random fluctuations
(noise) on the overall system performance. Manufacturers typically specify
this as an Allan/Hadamard Variance (AVAR/HVAR) profile in the -integration -
_time_ domain, yet, stochastic processes governing the noise take place in the
-Fourier - _frequency_ domain in the shape of a Power Spectral Density (PSD)
function. Both are second-moment measures of the time series, however, it is
only possible to translate unambiguously _from_ the PSD _to_ the AVAR/HVAR,
not vice versa, except in the case of a single noise type, which is severely
limiting in real-life applications.

This note elaborates an analytical method to generate an approximated PSD
expressed as a set of power-laws defined in specific intervals in the
frequency domain, starting from an AVAR/HVAR expressed a set of power-laws in
the time domain. The proposed algorithm is straightforward to implement,
applicable to all noise types (and combinations thereof) and can be self-
validated by reconstructing the corresponding AVAR/HVAR by direct calculus. We
also report on its limitations of and analytical expressions of the continuous
version of this algorithm. Coupling with well-established algorithms relying
on the PSD for power-law noise generation, the ensuing method encompasses the
capability for generating multi-colored noise in end-to-end simulations, as
demonstrated hereby for NASA’s Deep Space Atomic Clock.

##  I Introduction

Stochastic processes governing the noise of oscillator-based devices are of
great concern to engineers and scientists relying on systems distributing
high-quality timing and frequency information. An illustrative example lies in
the monitoring of frequency standards onboard GNSS, the main providers of
navigation and timing signals. Aside from the clock’s inherent instability,
periodic effects such as temperature and radiation pressure induce
fluctuations on the clock signal. Moreover, artificial fluctuations can arise
as a consequence of the the radiometric/optical link used in the estimation
process Sesia  _et al._ (2011), especially when using the same data type for
orbit determination and time transfer Dirkx, D.  _et al._ (2016). These
effects result in a convoluted noise spectrum which is difficult to model in
numerical simulations.

Considering the breadth and diversity of manufacturers and users of frequency
and timing signals, an IEEE special issue iee (2022) standardizes the
techniques and physical quantities used to measure and characterize
instabilities in instruments, across time and frequency domains. The Allan
Variance (AVAR), introduced in 1966 Allan (1966), is the most widespread
statistical measure for clock and oscillator instabilities in the _time_
domain. It intuitively displays the uncertainty one may expect after an
integration time $\tau$ following the last synchronization, and helps to
isolate the effect of spurs in the signal Rubiola (2008). Furthermore, the
AVAR and its derivatives (see Riley and Howe (2008) for a concise yet complete
description) are straightforwardly measured for long averaging times using a
time interval counter. Algorithms have been proposed Zucca and Tavella (2005);
Galleani (2008) that can generate clock noise directly from the ADEV, however,
these do not guarantee a high degree of flexibility when dealing with a
multitude of noise types. In fact, these variances cannot completely encompass
the covariance properties of the noise: the AVAR, for instance, is impervious
to the even-symmetry components of time signals (with respect to $t=\tau$) and
the Modified VAR is insensitive to $t=3/2\tau$ Rubiola (2008), which is why
some noise types have the same, indistinguishable signature.

The Power Spectral Density (PSD) must be identified for a comprehensive
characterization of the dominant noise type(s), setting confidence intervals
of the device under test, or making bias corrections for specific variance
estimators Riley and Howe (2008). This is conveniently assumed to follow (a
combination) of power-laws

$S_{y}(f)=hf^{\alpha}$ (1)

where $f$ is the Fourier (or sideband) frequency in Hertz, $h$ the intensity
coefficient and $\alpha$ is the power law exponent, which distinguishes the
noise processes for integer values. The main advantage of the PSD is its
deterministic nature -in the limit of an infinite sample sequence- yet its
measurement is limited by the bandwidth of the signal and/or the measurement
system, which is why a cutoff frequency $f_{H}$ is defined as the minimum of
the two iee (2022).

In light of these arguments, one may prefer to characterize clock and
oscillator behaviour in both time and frequency domain, hence the need to
translate between the two. Although the conversion from PSD to ADEV is
uniquely defined (assuming a $f_{H}$) the opposite is not true. Essentially,
this is an ill-posed problem Greenhall (1997). A direct conversion may only
take place in the case of a single slope of the time variance in log-log
plane, corresponding to a single noise exponent $\alpha$, which is is a
constricting simplification for systems affected by a multitude of noise
types, most notably GNSS.

In efforts to circumvent this limitation, this technical note elaborates an
analytical method to generate an approximated PSD expressed as a set of power-
laws defined in specific frequency intervals, starting from an Allan/Hadamard
Variance (AVAR/HVAR) expressed a set of power-laws in the time domain.

The manuscript is sectioned as follows. Sec. II outlines the theoretical basis
of the method for a set of cases: 1) the AVAR $\iff$ PSD for a single slope 2)
the exact PSD $\Rightarrow$ AVAR when the former is a continuous function
formed by a set of power laws valid in discrete intervals 3) the novel
proposed algorithm for ADEV/HDEV $\Rightarrow$ PSD when the former is composed
by a set of power laws valid in discrete intervals. In Sec. III we show a set
of numerical tests based on Allan variances of real clocks and we discuss the
limits of our method: asymptotic behavior (Sec. III.2) and the passage to the
continuous case (Sec. IV). Finally, in Sec. V we draw the possible
applications and conclusions of this work.

##  II Mathematical framework

The mathematical model of timing signals and the applicability of the AVAR is
briefly recalled hereby.
For a given clock, represented as an oscillator with nominal frequency
$\nu_{0}$, we define the phase as

$\Phi(t)=2\pi\nu_{0}t+\varphi(t)=2\pi\nu_{0}\left[t+x(t)\right]$ (2)

where $\varphi(t)$ is the random component and $x(t)$ is the time-error
function (i.e. the difference between the time of the clock and a reference
”real” time $t$). The instantaneous frequency is defined as

$\nu(t)=\frac{1}{2\pi}\frac{\mathrm{d}\Phi}{\mathrm{d}t}.$ (3)

For specific applications, the fractional-frequency error

$y(t)=\dot{x}(t)=\frac{\nu(t)}{\nu_{0}}-1$ (4)

is more convenient.
From the above definitions, the relations among the (single-sided) phase,
time, and frequency fluctuation power spectra ($S_{\varphi}(f),S_{x}(f)$ and
$S_{y}(f)$, respectively) are

$S_{\varphi}(f)=(2\pi\nu_{0})^{2}S_{x}(f);\qquad S_{y}(f)=(2\pi
f)^{2}S_{x}(f).$ (5)

The autocorrelation function of the fractional-frequency error is (for real
data)

$R_{y}(t)=\lim_{T\to+\infty}\frac{1}{2T}\int_{-T}^{+T}y(t^{\prime})y(t^{\prime}+t)\mathrm{d}t^{\prime}$
(6)

where $T\to+\infty$ indicates the passage from periodic to nonperiodic
signals. For the Wiener-Khinchin theorem $R_{y}(t)$ and $S_{y}(f)$ are Fourier
Transform/Inverse Fourier Transform couple

$S_{y}(f)=\int_{0}^{+\infty}R_{y}(t)\cos(2\pi ft)\mathrm{d}t;\qquad
R_{y}(t)=\int_{0}^{+\infty}S_{y}(f)\cos(2\pi ft)\mathrm{d}f$ (7)

(single-sided). Since measurements are spaced by an interval of time $\tau$
(i.e. integration time), we define the following estimators

$y_{n+1}(t,\tau)=\frac{y_{n}(t+\tau)-y_{n}(t)}{c_{n+1}/c_{n}};\qquad n=1,...$
(8)

where

$y_{1}(t,\tau)=\frac{x(t+\tau)/\tau-x(t)/\tau}{c_{1}}$ (9)

and $c_{n}$ are numerical coefficients (see below).
The autocorrelation of $y_{n}(t,\tau)$ is

$R_{y,n}(t,\tau)=\lim_{T\to+\infty}\frac{1}{2T}\int_{-T}^{+T}y_{n}(t^{\prime},\tau)y_{n}(t^{\prime}+t,\tau)\mathrm{d}t^{\prime};\qquad
n=1,...$ (10)

By inserting Eq. (8) into Eq. (10) it can be demonstrated that (see e.g. Van
Vliet and Handel (1982); Greenhall (1997))

$R_{y,n}(t,\tau)=\frac{2^{2n-2}}{c_{n}^{2}}\int_{0}^{+\infty}S_{y}(f)\cos(2\pi
ft)\frac{\sin^{2n}(\pi f\tau)}{(\pi f\tau)^{2}}\mathrm{d}f;\qquad n=1,...$
(11)

and the coefficients $c_{n}$ are obtained by conventionally imposing that

$R_{y,n}(0,\tau)=R_{y,1}(0,\tau)\quad\forall\,n$ (12)

in the case of white noise (i.e. $S_{y}(f)$=constant). This corresponds to
Makdissi  _et al._ (2010)

$c_{n}^{2}=c_{1}^{2}\frac{2^{2n-2}\Gamma(n-1/2)}{\sqrt{\pi}\,\Gamma(n)}\qquad
n=2,...$ (13)

assuming $c_{1}=1$ one finds $c_{2}^{2}=2$, $c_{3}^{2}=6$, $c_{4}^{2}=20$,
$c_{5}^{2}=70$, etc.
By inserting Eq. (13) into Eq. (11) one obtains

$R_{y,n}(t,\tau)=\frac{\sqrt{\pi}\,\Gamma(n)}{\Gamma(n-1/2)}\int_{0}^{+\infty}S_{y}(f)\cos(2\pi
ft)\frac{\sin^{2n}(\pi f\tau)}{(\pi f\tau)^{2}}\mathrm{d}f.$ (14)

Finally, we define $R_{y,2}(0,\tau)$ as the Allan variance (AVAR) and
$R_{y,3}(0,\tau)$ as the Hadamard variance (HVAR) (Cutler and Searle (1966);
Van Vliet and Handel (1982)). This latter, is able to deal with more divergent
noise sources, characterized by power spectra $\propto f^{\alpha}$ with
$-5<\alpha\leq-2$ that cannot be handled by the AVAR Riley and Howe (2008).
From Eq. (8), the corresponding estimators are

$\displaystyle y_{2}(t,\tau)$
$\displaystyle=\frac{y_{1}(t+\tau)-y_{1}(t)}{\sqrt{2}}=\frac{x(t+2\tau)-2x(t+\tau)+x(t)}{\tau\sqrt{2}}$
(15) $\displaystyle y_{3}(t,\tau)$
$\displaystyle=\frac{y_{2}(t+\tau)-y_{2}(t)}{\sqrt{3}}=\frac{x(t+3\tau)-3x(t+2\tau)+3x(t+\tau)-x(t)}{\tau\sqrt{6}}.$
(16)

Basing on Eq. (14) with $n=(2,3)$ the single-sided AVAR/HVAR
$\sigma_{y}^{2}(\tau$) can be expressed as

$\sigma_{y}^{2}(\tau)=2\int_{0}^{\infty}S_{y}(f)q(\pi\tau
f)\frac{\sin^{4}(\pi\tau f)}{(\pi\tau f)^{2}}\,\mathrm{d}f;$ (17)

where $q(\pi\tau f)=1$ for the AVAR and $q(\pi\tau f)=4/3\sin^{2}(\pi\tau f)$
for the HVAR.
The conditions for the convergence of the integral into Eq. (17) are discussed
in the following sections.

###  II.1 Single slope case

A power-law in the frequency domain, describing a power/amplitude spectral
density (PSD/ASD), has a direct correspondence to a power-law in the
integration time domain, describing an Allan variance/deviation (AVAR/ADEV),
and vice versa if and only if there is a single slope in the $]0,+\infty[$ $f$
and $\tau$ domains Burgoon and Fischer (1978) (the lower limit is included if
$\alpha\geq 0$). This univocal correspondence is easily demonstrated by
substituting Eq. (1) into Eq. (17).
By making the change $z=\pi\tau f$, Eq. (17) becomes

$\sigma_{y}^{2}(\tau)=2h\frac{I^{\infty}(\alpha)}{\pi^{\alpha+1}}\tau^{-\alpha-1}$
(18)

where the integral

$I^{\infty}(\alpha)=\int_{0}^{\infty}q(z)\frac{\sin^{4}z}{z^{2-\alpha}}\,\mathrm{d}z$
(19)

is is convergent if $-3<\alpha<+1$ (AVAR case) and $-5<\alpha<+1$ (HVAR case).
Analogously, if we define a single-sloped AVAR/HVAR

$\sigma^{2}_{y}(\tau)=B\tau^{\mu}$ (20)

(where $B$ is constant), $h$ and $\alpha$ in Eq. (1) are given by

$h=\frac{B}{2\pi^{\mu}J^{\infty}(\mu)};\qquad\alpha=-\mu-1;$ (21)

where we defined

$J^{\infty}(\mu)=I^{\infty}(-\mu-1)=\int_{0}^{\infty}q(z)\frac{\sin^{4}z}{z^{\mu+3}}\,\mathrm{d}z.$
(22)

The integral $J^{\infty}(\mu)$ is convergent if $-2<\mu<+2$ (AVAR case) and
$-2<\mu<+4$ (HVAR case).
Therefore, there is a biunivocal correspondance between $(h,\alpha)$
coefficients of the PSD and the $(B,\mu)$ coefficients of the AVAR/HVAR given
by Eq. (21) QED.
Details about the functions $I^{\infty}(\alpha)=J^{\infty}(-\alpha-1)$ are
reported in Appendix A.
To deal with $\alpha\geq 1$ (or $\mu\leq-2$), an upper cutoff frequency
$f_{H}$ to assure the convergence is introduced. In this case the integral
contains oscillating terms which are usually neglected if we limit to $\tau\gg
1/(2\pi f_{H})$. As a consequence of this approximation, different slopes in
frequency domain can correspond to the same slope in time domain, making
impossible to obtain the PSD from the AVAR. This is the case of the flicker
and white phase modulation noises ($\alpha=$ +1 and +2, respectively) both of
which correspond, neglecting oscillating terms, to an AVAR $\propto\tau^{-2}$.
This problem can be avoided by introducing the modified Allan variance Riley
and Howe (2008).
In the following we will limit our analysis to the cases for which the
integral of Eq. (19) is convergent.

###  II.2 Several power laws case: exact formula to obtain AVAR from PSD

Defining an arbitrary (i.e. not necessarily equally-spaced) set of $n$
coefficients $\alpha_{i}$ and $n-1$ frequencies $f_{i}$ with
$f_{1}<f_{2}<...<f_{n-1}$, we may write the PSD as a continuous function
$S_{y}(f)$ expressed as a set of power laws

$S_{y}(f)=\begin{cases}h_{1}f^{\alpha_{1}}&\mbox{if $f<f_{1}$}\\\
h_{i}f^{\alpha_{i}}&\mbox{if $f_{i-1}<f<f_{i}$, with $2\leq i\leq n-1$}\\\
h_{n}f^{\alpha_{n}}&\mbox{if $f>f_{n-1}$}.\end{cases}$ (23)

In the log-log plane (here and below $\log=\log_{10}$ and $\ln=\log_{e}$)

$\log S_{y}(f)=\log h_{i}+\alpha_{i}\log f;\qquad f_{i-1}<f<f_{i}.$ (24)

Defining $Q_{asd}=\sqrt{S_{y}(f_{1})}$ (i.e. the value of the ASD at
$f=f_{1}$) we obtain $h_{1}=Q_{asd}^{2}f_{1}^{-\alpha_{1}}$.
Since $S_{y}(f)$ is a continuous function we have

$h_{i+1}=h_{i}f_{i}^{\alpha_{i}-\alpha_{i+1}};\qquad i=1,...,n-1.$ (25)

To pass from PSD to AVAR we insert Eq. (23) into Eq. (17)

$\sigma_{y}^{2}(\tau)=2\sum_{i=1}^{n}\int_{f_{i-1}}^{f_{i}}h_{i}f^{\alpha_{i}}q(\pi\tau
f)\frac{\sin^{4}(\pi\tau f)}{(\pi\tau f)^{2}}\mathrm{d}f$ (26)

where we defined $f_{0}=0$ and $f_{n}=\infty$. With the variable change
$z=\pi\tau f$ we obtain

$\sigma_{y}^{2}(\tau)=2\sum_{i=1}^{n}\frac{h_{i}}{(\pi\tau)^{\alpha_{i}+1}}\int_{\pi\tau
f_{i-1}}^{\pi\tau f_{i}}q(z)\frac{\sin^{4}z}{z^{2-\alpha_{i}}}\,\mathrm{d}z.$
(27)

Therefore, a PSD in the form of Eq. (23) exactly corresponds to an AVAR in the
form

$\sigma^{2}_{y}(\tau)=\sum_{i=1}^{n}C_{i}(\tau)\tau^{-\alpha_{i}-1}$ (28)

where

$C_{i}(\tau)=\frac{2h_{i}}{{\pi}^{\alpha_{i}+1}}\int_{\pi\tau
f_{i-1}}^{\pi\tau
f_{i}}q(z)\frac{\sin^{4}z}{z^{2-\alpha_{i}}}\,\mathrm{d}z\qquad i=1,...,n$
(29)

act as a weight function. In particular

$C_{1}(\tau)=\frac{2h_{1}}{{\pi}^{\alpha_{1}+1}}\int_{0}^{\pi\tau
f_{1}}q(z)\frac{\sin^{4}z}{z^{2-\alpha_{1}}}\,\mathrm{d}z;\qquad
C_{n}(\tau)=\frac{2h_{n}}{{\pi}^{\alpha_{n}+1}}\int_{\pi\tau
f_{n-1}}^{\infty}q(z)\frac{\sin^{4}z}{z^{2-\alpha_{n}}}\,\mathrm{d}z.$ (30)

For $\tau\rightarrow 0$ the $n$-th term of Eq. (28) becomes the dominant one
while the integral tends to $I^{\infty}(\alpha_{n})$ while, for
$\tau\rightarrow\infty$, the first term becomes dominant and the integral
tends to $I^{\infty}(\alpha_{1})$. Therefore,

$\displaystyle\sigma_{y}^{2}(\tau\rightarrow 0)$
$\displaystyle\approx\tau^{-\alpha_{n}-1}\frac{2h_{n}}{\pi^{\alpha_{n}+1}}I^{\infty}(\alpha_{n});$
(31) $\displaystyle\sigma_{y}^{2}(\tau\rightarrow\infty)$
$\displaystyle\approx\tau^{-\alpha_{1}-1}\frac{2h_{1}}{\pi^{\alpha_{1}+1}}I^{\infty}(\alpha_{1}).$
(32)

Independently of $n$, the asymptotic behaviors at $\tau\rightarrow 0$ and
$\tau\rightarrow\infty$ are given by adopting the single-slope formula.
Finally, we want to express the AVAR given by Eq. (28) and Eq. (29) in the
form

$\sigma_{y}^{2}(\tau)=\sum_{i=1}^{n}B_{i}(\tau)\tau^{\mu_{i}}$ (33)

with $B_{1}\tau^{\mu_{1}}$ at $\tau\approx 0$ and $B_{n}\tau^{\mu_{n}}$ for
large $\tau$. Therefore, we need to reverse the indices by passing from $i$ to
$n-i+1$

$\mu_{i}=-\alpha_{n-i+1}-1;\qquad i=1,...,n$ (34)

and, after substituting Eq. (34) into Eq. (29), we obtain

$B_{i}(\tau)=C_{n-i+1}(\tau)=\frac{2h_{n-i+1}}{{\pi}^{\alpha_{n-i+1}+1}}\int_{\pi\tau
f_{{n-i}}}^{\pi\tau
f_{n-i+1}}q(z)\frac{\sin^{4}z}{z^{2-\alpha_{n-i+1}}}\,\mathrm{d}z.$ (35)

It is easy to verify that for $\tau\to 0$ the AVAR reduces to $\approx
B_{1}\tau^{\mu_{1}}$, and for $\tau\to\infty$ the AVAR reduces to $\approx
B_{n}\tau^{\mu_{n}}$.

###  II.3 Several power laws case: approximated formula to obtain PSD from
AVAR

We hereby propose the novel method to obtain a reasonable approximation of the
PSD corresponding to a given AVAR/HVAR. Although this is not an exact
conversion, one may always apply the exact formula Eq. (35) on the yielded PSD
to validate the output against the original input AVAR/HVAR.
For a given set of $n$ values $\\{\mu_{1},...,\mu_{n}\\}$ and $n-1$ nodes
$\\{\tau_{1},...,\tau_{n-1}\\}$ with $\tau_{1}<\tau_{2}<...<\tau_{n-1}$ we
define the AVAR as

$\sigma^{2}_{y}(\tau)=\begin{cases}B_{1}\tau^{\mu_{1}}&\mbox{if
$\tau<\tau_{1}$}\\\ B_{i}\tau^{\mu_{i}}&\mbox{if $\tau_{i-1}<\tau<\tau_{i}$,
with $2\leq i\leq n-1$}\\\ B_{n}\tau^{\mu_{n}}&\mbox{if
$\tau>\tau_{n-1}$.}\end{cases}$ (36)

As was for $S_{y}(f)$, in the log-log plane we have

$\log[\sigma_{y}^{2}(\tau)]_{i}=\log
B_{i}+\mu_{i}\log\tau;\qquad\tau_{i-1}<\tau<\tau_{i}.$ (37)

By defining $Q_{adev}=\sigma_{y}(\tau_{1})$ we obtain
$B_{1}=Q_{adev}^{2}\tau_{1}^{-\mu_{1}}$ and the continuity constraint of
$\sigma^{2}_{y}(\tau)$ leads to

$B_{i}\tau_{i}^{\mu_{i}}=B_{i+1}\tau_{i}^{\mu_{i+1}};\qquad i=1,...,n-1.$ (38)

From Eq. (31) and Eq. (32) we obtain $h_{n}$ and $h_{1}$, respectively, as
functions of $B_{1}$ and $B_{n}$. These are the asymptotic behaviors at
$f\rightarrow\infty$ and $f\rightarrow 0$, respectively.
The fundamental step is to extend this approach as to include the intermediate
$h_{i}$. Therefore, the approximated PSD will resemble the form of Eq. (23).
From Eq. (34) we have

$\alpha_{i}=-\mu_{n-i+1}-1;\qquad i=1,...,n.$ (39)

By defining

$I_{i}^{\infty}=\int_{0}^{\infty}q(z)\frac{\sin^{4}z}{z^{2-\alpha_{i}}}\,\mathrm{d}z;\qquad
J_{i}^{\infty}=\int_{0}^{\infty}q(z)\frac{\sin^{4}z}{z^{3+\mu_{i}}}\,\mathrm{d}z\qquad
i=1,...,n;$ (40)

we obtain

$I_{i}^{\infty}=J_{n-i+1}^{\infty};\qquad J_{i}^{\infty}=I_{n-i+1}^{\infty}.$
(41)

From Eq. (35), and integrating between zero and $+\infty$ we obtain 111The
choice of integrating between zero and infinity has been done for the sake of
simplicity. However, this imply that an exact convergence to the real PSD is
impossible, also in the continuous case (see Sec. IV). The discrepancy being
proportional to the degree of convexity/concavity of the AVAR.

$h_{i}=\frac{B_{n-i+1}}{2J_{n-i+1}^{\infty}\pi^{\mu_{n-i+1}}}=\frac{B_{n-i+1}\pi^{\alpha_{i}+1}}{2I_{i}^{\infty}};\qquad
i=1,...,n.$ (42)

Or, equivalently

$h_{n-i+1}=\frac{B_{i}}{2J_{i}^{\infty}\pi^{\mu_{i}}};\qquad i=1,...,n.$ (43)

The final step is to obtain the frequency nodes $f_{i}$ ($1\leq i\leq n-1$)
where the $S_{y}(f)$ change-of-slopes occur. Using Eq. (25) (i.e. the
continuity constraint for the PSD) we obtain

$f_{i}=\left(\frac{h_{i}}{h_{i+1}}\right)^{1/(\alpha_{i+1}-\alpha_{i})}=\left[\frac{B_{n-i+1}}{B_{n-i}}\frac{\pi^{\alpha_{i}+1}}{\pi^{\alpha_{i+1}+1}}\frac{I_{i+1}^{\infty}}{I_{i}^{\infty}}\right]^{1/(\alpha_{i+1}-\alpha_{i})};\qquad
i=1,...,n-1.$ (44)

Using also Eq. (38) (i.e. the continuity constraint of the AVAR/HVAR), we
obtain

$\frac{B_{n-i+1}}{B_{n-i}}=\tau_{n-i}^{\mu_{n-i}-\mu_{n-i+1}}=\tau_{n-i}^{\alpha_{i}-\alpha_{i+1}};\qquad
i=1,...,n-1.$ (45)

Therefore by introducing Eq. (45) into Eq. (44) we obtain the frequency nodes

$f_{i}=\frac{1}{\pi\tau_{n-i}}\left[\frac{I_{i+1}^{\infty}}{I_{i}^{\infty}}\right]^{1/(\alpha_{i+1}-\alpha_{i})}=\frac{1}{\pi\tau_{n-i}}\left[\frac{J_{n-i}^{\infty}}{J_{n-i+1}^{\infty}}\right]^{1/(\mu_{n-i+1}-\mu_{n-i})};\qquad
i=1,...,n-1$ (46)

or

$f_{n-i}=\frac{1}{\pi\tau_{i}}\left[\frac{J_{i}^{\infty}}{J_{i+1}^{\infty}}\right]^{1/(\mu_{i+1}-\mu_{i})};\quad
i=1,...,n-1.$ (47)

Hereinabove, the PSD is expressed as a set of power-law functions
$h_{i}f^{\alpha_{i}}$ in the intervals $[f_{i-1},f_{i}]$ (i.e. as Eq. (23))
where $h_{i}$ coefficients are given by Eq. (42), $f_{i}$ nodes by Eq. (46)
and $\alpha_{i}$ coefficients by Eq. (39). This represents a good
approximation of the input AVAR/HVAR($B_{i},\tau_{i},\mu_{i}$) defined in Eq.
(36), as will be demonstrated in the following section.
It must be said that the problem of converting the intersection nodes
$f\iff\tau$ was tackled by Burgoon and Fischer (1978), which was limited to
the (single-node) 2-sloped case. Table 2 of Burgoon and Fischer (1978) reports
a set of coefficients to convert the nodes for integer values of $\mu$
(including the cases $\mu=-3,-2$ that require the introduction of an upper
cutoff frequency). The other reported coefficients can be easily obtained by
Eq. (46) with $n=2$ and $\mu=-1,0,+1$. The conversion of the ”nodes” is not
trivial as it strongly depends on the values of the slopes. As pointed out by
Burgoon and Fischer (1978), it can lead to large errors if, for example, one
adopts $f_{i}=1/\tau_{n-1}$.
Furthermore, it should be noted that Vernotte  _et al._ (1993) describe the
”multivariance method” to infer the noise coefficient set which best fits a
set of time variance measurements in a weighted-least-squares sense. However,
for high-fidelity results one must implement a large set of variances, and the
frequency intervals where the coefficient sets are valid is not specified.

##  III Numerical tests

Having defined the recipe for converting the AVAR to PSD, the efficacy of this
conversion is hereby tested, considering two schemes:

  1. 1.

AVAR$\to$PSD$\to$AVAR: given an arbitrary AVAR/HVAR in the form of Eq. (36),
by means of our algorithm in Eq. (43) we generate the approximate PSD, and
then compare the corresponding AVAR/HVAR with the original one.

  2. 2.

PSD$\to$AVAR$\to$PSD: given an input PSD spectrum, by Eq. (17) we calculate
the corresponding AVAR and, using the method described in Sec. II.3, we re-
obtain the PSD and we compare it with the input one.

###  III.1 Real Allan deviations (Algorithm Fidelity in the domain?)

In Fig. 1 we report some applications of the scheme AVAR$\to$PSD$\to$AVAR. We
consider three AVARs relative to the following cases:

  * •

AccuBeat Ultra Stable Oscillator (USO)222https://www.accubeat.com/uso, a high
stability crystal quartz oscillator tailored for deep space exploration, and
recently selected to fly onboard the ESA JUICE mission to the Galilean moons
Shapira  _et al._ (2016);

  * •

Orolia?s space-qualifed RAFS333Orolia datasheet: https://safran-navigation-
timing.com/product/rafs/, currently state-of-the-art rubidium clock selected
by ESA to provide the frequency reference onboard the Galileo Second
Generation navigation satellites (in this case we considered the HDEV).

  * •

A fictitious and willingly complex ADEV with $\mu=\\{-3/2,$ $-1,$ $0,$ $+1,$
$+3/2\\}$, $\tau=\\{2^{-2},2,2^{2},2^{4}\\}$ and $Q_{adev}=10^{-14.5}$. This
corresponds to $\alpha=\\{-5/2,$ $-2,$ $-1,$ $0,$ $+1/2\\}$, $f=\\{4.473\times
10^{-6},$ $1.053\times 10^{-3},$ $0.180,$ $33.87\\}$ and
$Q_{asd}=10^{-11.71}$.

In all cases we found a good agreement between the input and calculated Allan
variances.

Figure 1: Left: ASD obtained by the input ADEV (right, red lines). Right:
input ADEV vs. ADEV (black line) corresponding to the ASD on the left. Top to
bottom: AccuBeat, RAFS clocks (Hadamard deviation in this case). The last one
been invented (see text for details).

Moreover, we tested the scheme PSD$\to$AVAR$\to$PSD in the case of a
Lorentzian frequency noise (i.e. $S_{y}(f)=(1+f^{2})^{-1}$). Using Eq. (17),
we numerically calculated the AVAR at 100 values of the integration time. We
used this set of nodes and slopes to re-obtain the PSD using the method
described above. In Fig. 2 (left panel) we compare the Lorentzian spectrum and
the calculated one (red and black lines, respectively). The corresponding ADEV
is reported in the right panel.

Figure 2: Left: Lorentzian noise spectrum (red line) and the profile
calculated basing on the ADEV reported on right panel (black line).

Finally, we perform a more complete test. Its purpose is to show that the PSD
obtained by the algorithm can be used to generate random noise whose ADEV
agrees with the original one.
We consider three datasets, all relative to the trapped-ion Deep Space Atomic
Clock (DSAC, that completed its two-year mission in low-Earth orbit on Sept.
18, 2021 Burt  _et al._ (2021)): 1) expected time fluctuations of DSAC-1 in
flight (without other noise sources) 2) the same for DSAC-2 next generation
clock and 3) ?raw? offsets, which includes measurement system noise, GPS
noise, orbit determination noise, and DSAC-1 clock noise (this represents an
upper bound for the clock itself).
For each dataset, starting from $N$ time-fluctuations data $x_{i}$ spaced by
$\tau$, we directly obtain the AVAR as

$\sigma^{2}_{y}(\tau)={1\over
2(N-2)\tau^{2}}\sum_{n=1}^{N-2}({x}_{n+2}-2x_{n+1}+x_{n})^{2}$ (48)

for a set of integration times $\tau_{j}$.
By applying the algorithm described in Sec. II.3 we obtain the set of
parameters $(h_{j},\alpha_{j},f_{j})$ relative to the PSDs $S_{y}(f)$ and
$S_{x}(f)=S_{y}(f)/(2\pi f)^{2}$.
By generalising the noise-generating algorithms of Timmer and Koenig (1995) to
to an arbitrary PSD (as opposed to a single slope) we interpolate the
$S_{x}(f)$ function with $n$ (even) frequencies between 0 and $f_{n}=1/\Delta
t$, where $\Delta t$ is the desired spacing of the output data. We therefore
generate the following complex vector $\mathbf{F}$

$X_{j}=\frac{\sqrt{S_{x}(f_{j})}}{2}\left[N_{R}^{i}(0,1)+iN_{j}^{i}(0,1)\right]\quad\mbox{for
}j=1,n/2+1;\qquad X_{n/2+1+j}=X^{*}_{n/2-j+1}\quad\mbox{for }j=1,...,n/2$ (49)

where used the symmetry property of the discrete Fourier transform of real
data. For this reason (Nyquist-Shannon theorem) all information is contained
below $f_{sup}=f_{n/2+1}=1/(2\Delta t)$. Coefficients
$N_{R}^{j}(0,1),N_{I}^{j}(0,1)$ are two normally distributed random numbers
and the factor 2 instead of $\sqrt{2}$ at the denominator is necessary since
$S_{x}(f)$ is a single-sided spectrum. It can be demonstrated that
$<X_{j}^{*}X_{j}>=S_{x}(f_{j})$.
Finally, we obtain the $n$-elements ($n$ is similar to the number of elements
of the input datasets) vectors $\mathbf{x}$ and $\mathbf{y}$ of the simulated
phase/time noises

$\mathbf{x}=IDFT(\mathbf{X})\sqrt{\frac{n-1}{\Delta t}};\qquad
y_{i}=\frac{x_{i}-x_{i-1}}{\Delta t}\qquad i=2,...,n$ (50)

where $IDFT$ is the inverse discrete Fourier transform.
Finally, by applying Eq. (48) to the simulated $\mathbf{y}$ vector, we obtain
the stability corresponding to the noise. In the left panel of Fig. 3, the
ADEVs from the input data (colored lines) are compared to the ADEVs of the
noise. In all cases, a good agreement is found, indicating that the algorithm
is able to deduce, from the input ADEVs, a faithful representation of their
spectral contents (the deduced ASDs are reported in right panel). For an
improved confidence at higher integration times, one may choose to compute the
OADEV, albeit this entails a computational expense.

Figure 3: Left: Allan deviations of the input data (black: DSAC-1, blue:
DSAC-2, red: ”raw” data) compared with the simulated time fluctuations noises
(grey lines). Right: ASDs estimated by the algorithm starting from the ADEVs
of the input data. Input data spacings are 1s (DSAC-1 and 2) and 30s (”raw”
data).

###  III.2 Asymptotic behaviour at zero and infinity

By inserting Eq. (43) into Eq. (35) we obtain

$\hat{\sigma}_{y}^{2}(\tau)=\sum_{i=1}^{n}\frac{B_{i}\tau^{\mu_{i}}}{J_{i}^{\infty}}\int_{\pi\tau
f_{n-i}}^{\pi\tau
f_{n-i+1}}q(z)\frac{\sin^{4}z}{z^{3+\mu_{i}}}\,\mathrm{d}z;\qquad(f_{0}=0;\quad
f_{n}=+\infty).$ (51)

Therefore, the approximated AVAR/HVAR is in the form

$\hat{\sigma}_{y}^{2}(\tau)=\sum_{i=1}^{n}W_{i}(\tau)B_{i}\tau^{\mu_{i}}$ (52)

where

$W_{i}(\tau)=\frac{1}{J_{i}^{\infty}}\int_{\pi\tau f_{n-i}}^{\pi\tau
f_{n-i+1}}q(z)\frac{\sin^{4}z}{z^{3+\mu_{i}}}\,\mathrm{d}z.$ (53)

By isolating the $i=1$ and $i=n$ terms into Eq. (51) we obtain

$\hat{\sigma}_{y}^{2}(\tau)=\frac{B_{1}\tau^{\mu_{1}}}{J_{1}^{\infty}}\int_{\pi\tau
f_{n-1}}^{\infty}q(z)\frac{\sin^{4}z}{z^{3+\mu_{1}}}\,\mathrm{d}z+\sum_{i=2}^{n-1}\frac{B_{i}\tau^{\mu_{i}}}{J_{i}^{\infty}}\int_{\pi\tau
f_{n-i}}^{\pi\tau
f_{n-i+1}}q(z)\frac{\sin^{4}z}{z^{3+\mu_{i}}}\,\mathrm{d}z+\frac{B_{n}\tau^{\mu_{n}}}{J_{n}^{\infty}}\int_{0}^{\pi\tau
f_{1}}q(z)\frac{\sin^{4}z}{z^{3+\mu_{n}}}\,\mathrm{d}z.$ (54)

Independently of the number of frequency nodes (i.e. the size of $\Delta
f_{n-i}$) when $\tau\rightarrow 0$, the calculated
$\hat{\sigma}_{y}^{2}(\tau)$ cannot converge to the input AVAR/HVAR since
$f_{1}$ and $f_{n-1}$ are not zero and infinity, respectively.
However, the integral into the first term tends to $J_{1}^{\infty}$, while
others tend to zero in all cases. Therefore, $\hat{\sigma}^{2}(\tau)$ is
asymptotic to $B_{1}\tau^{\mu_{1}}$ for $\tau\to 0$. On the contrary, when
$\tau\rightarrow+\infty$ the last integral tends to $J_{n}^{\infty}$ and the
others tend to zero (because we are integrating in a part of the domain where
the function is small). Therefore, $\hat{\sigma}^{2}(\tau)\rightarrow
B_{n}\tau^{\mu_{n}}$.
In Fig. 4 we report the results for an input ADEV with uniformly randomly
chosen $\mu$ parameters between -1.5 and 1.5 (red line). The asymptotic
behavior at zero and infinity is always assured.
The calculated  PSD (black line, left panel) corresponds to the ADEV in the
right panel (black line). Note that the larger is the degree of
concavity/convexity in the input ADEV, the larger are the discrepances.

Figure 4: Left: ASD obtained by the input ADEV (right, red lines). Right:
input ADEV vs. reconstructed ADEV (black line) from ASD on the left. The ADEV
profile has been invented (coefficients $\mu$ uniformly randomly chosen
between -1.5 and +1.5).

In Appendix B we demonstrate that if $n=2$ (i.e. two slopes only), the
reconstructed AVAR/HVAR, in an interval around $\tau_{1}$, is above the input
one if this latter is convex. It is below if the input function is concave.

##  IV Continuous case

We will calculate the continuous version of the method described above (i.e.:
when the number of sampling points tends to infinity but also the interval
$[f_{1},f_{n-1}]$ is expanded to $[0,+\infty[$).
From Eq. (28) and Eq. (29), the AVAR/HVAR corresponding to a PSD defined in
terms of power-laws is

$\hat{\sigma}^{2}_{y}(\tau)=\sum_{i=1}^{n}\tau^{-\alpha_{i}-1}\frac{2h_{i}}{{\pi}^{\alpha_{i}+1}}\int_{\pi\tau
f_{i-1}}^{\pi\tau
f_{i}}q(z)\frac{\sin^{4}z}{z^{2-\alpha_{i}}}\,\mathrm{d}z\,\qquad i=1,...,n;$
(55)

where $f_{0}=0$ and $f_{n}=+\infty$. Since (for $x$ near to $a$)

$\int_{a}^{x}f(t)\,\mathrm{d}t\approx f(a)(x-a)+\dots,$ (56)

by assuming an infinite number of frequency nodes, we get

$\sigma^{2}_{y}(\tau)=2\sum_{i=1}^{+\infty}h_{i}f_{i-1}^{\alpha_{i}}q(\pi\tau
f_{i-1})\frac{\sin^{4}(\pi\tau f_{i-1})}{(\pi\tau f_{i-1})^{2}}\,\Delta
f_{i-1}$ (57)

where $\Delta f_{i-1}=f_{i}-f_{i-1}$. In the continuous case $f_{i}\to f$
(e.g.: a variable defined in the $[0,+\infty[$ domain),
$\alpha_{i}\to\alpha(f)$ and $\Delta f_{i-1}\to\mathrm{d}f$. The summation
becomes an integral and we re-obtain Eq. (17) where $S_{y}(f)=h(f)$.
In the case of AVAR/HVAR, the passage to the continuous corresponds to
$\tau_{i}\to\tau$, $\mu_{i}\to\mu(\tau)$, $J_{i}^{\infty}\to
J^{\infty}[\mu(\tau)]$, $B_{i}\to B(\tau)$ and
$\sigma_{y}^{2}(\tau)=B(\tau)\tau^{\mu(\tau)}$.
By inserting Eq. (42) into Eq. (55) and using Eq. (39) we get

$\hat{\sigma}^{2}_{y}(\tau)=\sum_{i=1}^{n}\tau^{\mu_{n-i+1}}\frac{B_{n-i+1}}{J_{n-i+1}^{\infty}}\int_{\pi\tau
f_{i-1}}^{\pi\tau
f_{i}}q(z)\frac{\sin^{4}z}{z^{3+\mu_{n-i+1}}}\,\mathrm{d}z\,\qquad i=1,...,n.$
(58)

The subscript $n-i+1$ indicates that $B,\mu$ and $J^{\infty}$ are relative to
$\tau=\tau_{n-i+1}$ (e.g.: $B_{n-i+1}=B(\tau_{n-i+1})$). From Eq. (46) we
obtain

$\tau_{n-i+1}=\frac{1}{\pi
f_{i-1}}\left[\frac{J_{i-1}^{\infty}}{J_{i}^{\infty}}\right]^{1/(\mu_{i}-\mu_{i-1})};\qquad
i=1,...,n-1.$ (59)

Therefore, we can express $B,J^{\infty}$ and $\mu$ as functions of $f_{n-i}$.
Assuming infinite frequency nodes we get

$\hat{\sigma}^{2}_{y}(\tau)=\sum_{i=1}^{+\infty}\frac{B(f_{i-1})}{J^{\infty}(f_{i-1})\pi^{\mu(f_{i-1})}}f_{i-1}^{-1-\mu(f_{i-1})}q(\pi\tau
f_{i-1})\frac{\sin^{4}(\pi\tau f_{i-1})}{(\pi\tau f_{i-1})^{2}}\,\Delta
f_{i-1}.$ (60)

When passing to the continuous case, we define

$Z(\mu_{i})=\lim_{\Delta\mu_{i}\to
0}\left[\frac{J_{i}^{\infty}}{J_{i+1}^{\infty}}\right]^{1/\Delta\mu_{i}};\qquad\mbox{where
}\Delta\mu_{i}=\mu_{i+1}-\mu_{i}.$ (61)

The limit is an indeterminate form ($1^{\infty}$) and it is easy to
demonstrate that

$Z(\mu)=\exp{\left[-\frac{\mathrm{d}\ln
J^{\infty}(\mu)}{\mathrm{d}\mu}\right]}.$ (62)

Note that the limit $\Delta\mu\to 0$ means that the input function must be not
only continue but also ”smooth” (i.e. also its first derivative must be
continue).
Therefore, in the continuous case, Eq. (46) and Eq. (47) become

$f(\tau^{\prime})=\frac{Z[\mu(\tau^{\prime})]}{\pi\tau^{\prime}}.$ (63)

If allowed (i.e. if and only if $f(\tau^{\prime})$ is bijective), this
relation can be inverted to calculate $\tau^{\prime}(f)$.
Details about the function $Z(\mu)$ are reported in Appendix C.
Passing to the continuous case, Eq. (60) becomes

$\sigma^{2}_{y}(\tau)\approx\int_{0}^{+\infty}\frac{B(f)}{J^{\infty}[\mu(f)]\pi^{\mu(f)}}f^{-1-\mu(f)}q(\pi\tau
f)\frac{\sin^{4}(\pi\tau f)}{(\pi\tau f)^{2}}\,\mathrm{d}f$ (64)

the corresponding PSD is therefore

$S_{y}(f)=\frac{B[\tau^{\prime}(f)]}{2J^{\infty}[\mu(\tau^{\prime}(f))]\pi^{\mu[\tau^{\prime}(f)]}}f^{-1-\mu[\tau^{\prime}(f)]}$
(65)

where $\tau^{\prime}(f)$ is obtained by inverting Eq. (63).
To invert Eq. (63) and to apply Eq. (65), we need to calculate $\mu(\tau)$ and
$B(\tau)$ from a given input AVAR/HVAR $\sigma_{y}^{2}(\tau)$.
To this aim, we interpolate it with a set of power laws as in Eq. (36) with
the same continuity constraint. We define a set of nodes
$\\{\tau_{1},...,\tau_{n-1}\\}$ and we pass to the continuous case

$\mu_{i}=\frac{\log\sigma_{y}^{2}(\tau_{i})-\log\sigma_{y}^{2}(\tau_{i-1})}{\log\tau_{i}-\log\tau_{i-1}}\quad\to\quad\mu(\tau)=\frac{\mathrm{d}\log\sigma_{y}^{2}(\tau)}{\mathrm{d}\log\tau}\\\
$ (66)

and

$\log B_{i}=\log\sigma_{y}^{2}(\tau_{i})-\mu_{i}\log\tau_{i}\quad\to\quad\log
B(\tau)=\log\sigma_{y}^{2}(\tau)-\mu(\tau)\log\tau.$ (67)

If the AVAR/HVAR is given as an analytical function, Eq. (65) provides the
corresponding analytical PSD. It represents the limit of the (approximated)
PSD described in Sec. II.3 when the input AVAR/HVAR is interpolated by a
number of nodes that tends to infinity.
In Appendix D we report two checks of Eq. (65).

##  V Conclusions

Depending on the scenario, most notably due to equipment availability or
capability, it may be preferable to characterize clock and oscillator
instabilities in time or frequency domains, hence the need to translate
between the two. In this work we describe a simple algorithm to numerically
compute an approximated power spectral density (PSD) corresponding to an Allan
(or Hadamard) variance (AVAR/HVAR) given as an input. This may be applied
directly if the Allan deviation is given in terms of an arbitrary set of
joined power-laws defined in contiguous intervals of time, but we also report
the formula to be used if the Allan deviation is expressed in terms of an
analytical function.

The algorithm is summarized as follows:

  1. 1.

Express the input ADEV/HDEV in power-law form by selecting a set
($\tau_{i},\mu_{i})$ and computing $B_{i}$ by means of Eq. (38);

  2. 2.

Use Eq. (42) to compute the set $h_{i}$ from $B_{i}$;

  3. 3.

Translate the slopes to the frequency domain with Eq. (39);

  4. 4.

Compute frequency nodes from Eq. (46).

Having obtained the PSD, one may apply the exact formulation Eq. (29) to
reconstruct the input AVAR/HVAR. This validation exercise has demonstrated the
applicability of the proposed algorithm for ”spectralizing” a series of clocks
affected by combinations of noises, is illustrated in Eq. (1). Furthermore,
its limitations are investigated by applying it to ”extreme” (and irrealistic)
Allan deviations. The discrepancy between inverse-calculated ADEV/HDEV and the
input ADEV/HDEV is in general proportional to the local degree of
concavity/convexity, nonetheless, in all cases, the calculated PSD tends to
converge to the real one at the extremes of the frequency domain (zero and
infinity).

Foreseeable uses of this algorithm are twofold. A coupling with well-
established algorithms relying on PSD Timmer and Koenig (1995); Kasdin (1995)
would constitute a versatile and effective tool for generating multi-colored
noise series in time domain, for addressing the deleterious effects that
oscillator instabilities (standalone or embedded in a system with external
disturbances) produce on the overall performance. In turn, this may help to
optimize the operation of reference time scales for applications such as GNSS,
where predictions of the time deviation for free-running clocks are required.

A second foreseeable use is the computation of the autocorrelation matrix for
estimation filters which process non-gaussian observables, as in the context
of spacecraft navigation. When using a single datatype for the estimation of
orbit and clock parameters, as is prospected for the one-way navigation of the
next-generation of deep-space probes Ely  _et al._ (2018), errors can inject
in the clock estimates some residuals of a non-perfect orbit estimation. The
”true” autocorrelation may therefore help to effectively disentangle orbit and
clock parameters in the estimation filter, for an effective decoupling of
transnational and time dynamics.

##  Appendix A Coefficients of $I^{\infty}(\alpha)$

The integrai $I^{\infty}(\alpha)$ can be analytically expressed as

$I^{\infty}_{ADEV}(\alpha)=\Gamma[\alpha-1]\sin\left(\frac{\alpha\pi}{2}\right)\times
2^{-1-2\alpha}(1-2^{1+\alpha})$ (68)

and

$I^{\infty}_{HVAR}(\alpha)=\Gamma[\alpha-1]\sin\left(\frac{\alpha\pi}{2}\right)\times
2^{-1-2\alpha}\times 3^{1-\alpha}(4\times 3^{\alpha}-2^{\alpha}-5\times
6^{\alpha})$ (69)

where $\Gamma$ is the gamma function, which is undefined for non-positive
integers. In the range of values of our interest ($-3<\alpha<1$, AVAR;
$-5<\alpha<1$, HVAR), the functions reported into Eq. (68) and Eq. (69) are
undefined at $\alpha=\\{-2,-1,0\\}$ and $\alpha=\\{-4,-3,-2,-1,0\\}$ (HDEV),
respectively. However, it is easy to demonstrate that their domains can be
extended to these points by analytical continuation.
Plots of $J^{\infty}(\mu)=I^{\infty}(-\alpha-1)$ functions for Allan (left)
and Hadamard (right) deviations are reported in Fig. 5.
Coefficients of $I^{\infty}(\alpha)=J^{\infty}(-\mu-1)$ for some values of
$\alpha$ (or $\mu=-\alpha-1$) for Allan and Hadamard deviations are reported
in Table 1.

Figure 5: Allan (left) and Hadamard (right) functions $J^{\infty}(\mu)$. The
minimum is at: $\bar{\mu}=-0.25677$, $J(\bar{\mu})=0.682881$ (Allan) and
$\bar{\mu}=0.83161$, $J^{\infty}(\bar{\mu})=3.13152$ (Hadamard).
$\begin{array}[]{|l |l |l l |l l|
}\hline\cr\hline\cr\alpha&\mu&I^{\infty}_{AVAR}(\alpha)&&I^{\infty}_{HVAR}(\alpha)&\\\\[8.0pt]
\hline\cr-5&4&-&-&+\infty&+\infty\\\ -4&3&-&-&\dfrac{11\pi}{5}&6.9115\\\
-3&2&+\infty&+\infty&\log\left(\dfrac{1594323\sqrt{3}}{65536}\right)&3.74091\\\
-2&1&\dfrac{\pi}{3}&1.0472&\pi&3.14159\\\
-1&0&\log(2)&0.693147&\dfrac{3}{2}\log\left(\dfrac{256}{27}\right)&3.37401\\\
0&-1&\dfrac{\pi}{4}&0.785398&\dfrac{3\pi}{2}&4.71239\\\
1&-2&+\infty&+\infty&+\infty&+\infty\\\ \hline\cr\end{array}$ (70) Table 1:
Coefficients of $I^{\infty}(\alpha)=J^{\infty}(-\mu-1)$ for integer values of
$\alpha=-\mu-1$ for Allan and Hadamard deviations.

##  Appendix B Two slopes case ($n=2$)

Here we demonstrate that, in the two-slopes case,
$\hat{\sigma}_{y}^{2}(\tau_{1})>B_{1}\tau_{1}^{\mu_{1}}$ if and only if
$\mu_{2}>\mu_{1}$ and vice versa.
The reconstructed AVAR/HVAR is given by

$\hat{\sigma}_{y}^{2}(\tau)=\frac{B_{1}\tau^{\mu_{1}}}{J_{1}^{\infty}}\int_{\pi\tau
f_{1}}^{\infty}q(z)\frac{\sin^{4}z}{z^{3+\mu_{1}}}\,\mathrm{d}z+\frac{B_{2}\tau^{\mu_{2}}}{J_{2}^{\infty}}\int_{0}^{\pi\tau
f_{1}}q(z)\frac{\sin^{4}z}{z^{3+\mu_{2}}}\,\mathrm{d}z$ (71)

where

$\pi\tau
f_{1}=\frac{\tau}{\tau_{1}}\left(\frac{J_{1}^{\infty}}{J_{2}^{\infty}}\right)^{1/(\mu_{2}-\mu_{1})}.$
(72)

Since $B_{2}=B_{1}\tau_{1}^{\mu_{1}-\mu_{2}}$, we have

$\hat{\sigma}_{y}^{2}(\tau)=B_{1}\tau^{\mu_{1}}\left[\frac{1}{J_{1}^{\infty}}\int_{\pi\tau
f_{1}}^{\infty}q(z)\frac{\sin^{4}z}{z^{3+\mu_{1}}}\,\mathrm{d}z+\frac{(\tau/\tau_{1})^{\mu_{2}-\mu_{1}}}{J_{2}^{\infty}}\left(1-\int_{\pi\tau
f_{1}}^{\infty}q(z)\frac{\sin^{4}z}{z^{3+\mu_{2}}}\,\mathrm{d}z\right)\right].$
(73)

At $\tau=\tau_{1}$

$\hat{\sigma}_{y}^{2}(\tau_{1})-B_{1}\tau_{1}^{\mu_{1}}=B_{1}\tau_{1}^{\mu_{1}}\int_{0}^{\pi\tau_{1}f_{1}}q(z)\frac{\sin^{4}z}{z^{3}}\left(\frac{1}{J_{2}^{\infty}z^{\mu_{2}}}-\frac{1}{J_{1}^{\infty}z^{\mu_{1}}}\right)\,\mathrm{d}z.$
(74)

The integrand is equal to zero at

$z=\left(\frac{J_{1}^{\infty}}{J_{2}^{\infty}}\right)^{1/(\mu_{2}-\mu_{1})}=\pi\tau_{1}f_{1}$
(75)

(i.e.: the upper limit of integration). If $\mu_{2}>\mu_{1}$ the integrand is
always positive, so the value of the reconstructed AVAR/HVAR is larger than
the one of the input function (at the node $\tau_{1}$ but, for continuity,
also in an interval that contains $\tau_{1}$). On the contrary, if
$\mu_{2}<\mu_{1}$ the trial AVAR/HVAR is locally smaller than the input one.
The discrepancy being proportional to the difference between $\mu_{1}$ and
$\mu_{2}$ (in Fig. 6 we report two examples with a big change in slope
$|\mu_{2}-\mu_{1}|=3$). This is not, in general, true if $n>2$.

Figure 6: Top ASD (left) and ADEV (right) for a 2-slopes case. The ADEV
corresponding to the ASD on the left is above the input ADEV (red line) if
$\mu_{2}>\mu_{1}$ (here $\mu_{1}=-1.5;\mu_{2}=+1.5$) and below otherwise
(bottom panels, where $\mu_{2}=-1.5;\mu_{1}=+1.5$).

##  Appendix C The function $Z(\mu)$

From Eq. (62), the natural logarithms of $Z(\mu)$ are

$\ln[Z^{\infty}_{AVAR}(\mu)]=\frac{\pi}{2}\tan\left(\frac{\pi\mu}{2}\right)+\psi(-\mu-2)+\left(\frac{1}{1-2^{\mu}}-2\right)\ln(2)$
(76)

and

$\ln[Z^{\infty}_{HVAR}(\mu)]=\frac{\pi}{2}\tan\left(\frac{\pi\mu}{2}\right)+\psi(-\mu-2)+\left(\frac{2^{\mu+3}}{-2^{\mu+3}+3^{\mu+1}+5}-1\right)\ln(2)+\frac{\ln(3)}{3^{-\mu-1}\left(2^{\mu+3}-5\right)-1}$
(77)

where $\psi(x)=\Gamma^{\prime}(x)/\Gamma(x)$ is the digamma function.
The function $Z(\mu)$ is plotted in Fig. 7 for the Allan and Hadamard cases.
As was for $J^{\infty}(\mu)$, also $Z(\mu)$ must be extended to some points by
analytical continuation.
For example, in the ADEV case the points are $\mu=\\{-1,0,1\\}$ where the
limits are

$Z(-1)=e^{1-\gamma};\qquad Z(0)=\frac{e^{3/2-\gamma}}{2\sqrt{2}};\qquad
Z(1)=\frac{1}{8}e^{11/6-\gamma}$ (78)

where $\gamma=0.577215664...$ is the Euler?Mascheroni constant.

Figure 7: Allan (left) and Hadamard (right) functions $Z(\mu)$.

##  Appendix D Check of the continuous case formula

To check Eq. (65), we assume an AVAR in the form

$\sigma_{y}^{2}(\tau)=\sum_{i=1}^{N}b_{i}\tau^{\mu_{i}}$ (79)

and the corresponding PSD can be analytically calculated as

$S_{y}(f)=\sum_{i=1}^{N}\frac{b_{i}}{2\pi^{\mu_{i}}J^{\infty}(\mu_{i})}f^{-\mu_{i}-1}.$
(80)

We compare the PSD generated by Eq. (65) with Eq. (80). We consider a
”realistic” case: $N=3$, $b_{i}=1\,\forall i$ and $\mu_{i}=\\{-1,0,1\\}$. The
corresponding PSD is

$S_{y}(f)=2+\frac{1}{2f\ln 2}+\frac{3}{2\pi^{2}f^{2}}.$ (81)

By Eq. (66) and Eq. (67) we obtain

$\mu(\tau)=\frac{\tau^{2}-1}{\tau^{2}+\tau+1};\qquad
B(\tau)=\tau^{\frac{\tau+2}{\tau^{2}+\tau+1}-2}\left(\tau^{2}+\tau+1\right).$
(82)

In Fig. 8 (left panel) we compare the analytical ASD ($\sqrt{S_{y}(f)}$ from
Eq. (81), black line) with the one calculated by using Eq. (65) where
$\mu(\tau)$ and $B(\tau)$ are given by Eq. (82) (red line). In the right panel
we report the input ADEV ($\sigma_{y}(\tau)=\sqrt{\tau+1+1/\tau}$).
Finally, we consider the case $N=2,b_{i}=1\,\forall i$ and $\mu_{1}=-1.9$,
$\mu_{2}=+1.9$. In Fig. 9 we report the result: the discrepancy is due to the
strong difference between $\mu_{1}$ and $\mu_{2}$ (since $\mu$ must be between
-2 and 2, we are considering an extremely large degree of convexity).

Figure 8: Left: comparison between the analytical ASD (black line, i.e.: the
square root of Eq. (81)) and the calculated one (red line). Right: the input
ADEV (i.e.: $\sigma_{y}(\tau)=\sqrt{\tau+1+1/\tau}$).

Figure 9: The same as Fig. 8 but in the case $N=2$ and $\mu_{1}=-1.9$ and
$\mu_{2}=+1.9$ (see text for details).

## References

  * Sesia  _et al._ (2011) I. Sesia, L. Galleani, and P. Tavella, IEEE Transactions on Aerospace and Electronic Systems 47, 884 (2011).
  * Dirkx, D.  _et al._ (2016) Dirkx, D., Noomen, R., Visser, P. N. A. M., Gurvits, L. I., and Vermeersen, L. L. A., A&A 587, A156 (2016).
  * iee (2022) IEEE P1139/D17, November 2022 , 1 (2022).
  * Allan (1966) D. W. Allan, Proceedings of the IEEE 54, 221 (1966).
  * Rubiola (2008) E. Rubiola,  _Phase noise and frequency stability in oscillators_ (Cambridge University Press, 2008).
  * Riley and Howe (2008) W. Riley and D. Howe, “Handbook of frequency stability analysis,”  (2008).
  * Zucca and Tavella (2005) C. Zucca and P. Tavella, IEEE transactions on ultrasonics, ferroelectrics, and frequency control 52, 289 (2005).
  * Galleani (2008) L. Galleani, Metrologia 45, S175 (2008).
  * Greenhall (1997) C. Greenhall, in  _Proceedings of International Frequency Control Symposium_ (1997) pp. 358–365.
  * Van Vliet and Handel (1982) C. M. Van Vliet and P. H. Handel, Physica A: Statistical Mechanics and its Applications 113, 261 (1982).
  * Makdissi  _et al._ (2010) A. Makdissi, F. Vernotte, and E. D. Clercq, IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control 57, 1011 (2010).
  * Cutler and Searle (1966) L. Cutler and C. Searle, Proceedings of the IEEE 54, 136 (1966).
  * Burgoon and Fischer (1978) R. Burgoon and M. Fischer, in  _32nd Annual Symposium on Frequency Control_ (1978) pp. 514–519.
  * Note (1) The choice of integrating between zero and infinity has been done for the sake of simplicity. However, this imply that an exact convergence to the real PSD is impossible, also in the continuous case (see Sec. IV). The discrepancy being proportional to the degree of convexity/concavity of the AVAR.
  * Vernotte  _et al._ (1993) F. Vernotte, E. Lantz, J. Groslambert, and J. Gagnepain, IEEE Transactions on Instrumentation and Measurement 42, 342 (1993).
  * Note (2) https://www.accubeat.com/uso.
  * Shapira  _et al._ (2016) A. Shapira, A. Stern, S. Prazot, R. Mann, Y. Barash, E. Detoma, and B. Levy, in  _2016 European Frequency and Time Forum (EFTF)_ (2016) pp. 1–5.
  * Note (3) Orolia datasheet: https://safran-navigation-timing.com/product/rafs/.
  * Burt  _et al._ (2021) E. A. Burt, J. D. Prestage, R. L. Tjoelker, D. G. Enzer, D. Kuang, D. W. Murphy, D. E. Robison, J. M. Seubert, R. T. Wang, and T. A. Ely, Nature 595, 43 (2021).
  * Timmer and Koenig (1995) J. Timmer and M. Koenig, Astronomy and Astrophysics 300, 707 (1995).
  * Kasdin (1995) N. J. Kasdin, Proc. IEEE 83, 802 (1995).
  * Ely  _et al._ (2018) T. A. Ely, E. A. Burt, J. D. Prestage, J. M. Seubert, and R. L. Tjoelker, IEEE transactions on ultrasonics, ferroelectrics, and frequency control 65, 950 (2018).