††thanks: email<EMAIL_ADDRESS>
# Nonlinear optics from hybrid dispersive orbits

Yongjun Li ID Brookhaven National Laboratory, Upton 11973, New York, USA
Derong Xu  Brookhaven National Laboratory, Upton 11973, New York, USA  Victor
Smaluk  Brookhaven National Laboratory, Upton 11973, New York, USA  Robert
Rainer ID Brookhaven National Laboratory, Upton 11973, New York, USA

###### Abstract

In this paper we present an expansion of the technique of characterizing
nonlinear optics from off-energy orbits (NOECO) Olsson  _et al._ (2020) to
cover harmonic sextupoles in storage rings. The existing NOECO technique has
been successfully used to correct the chromatic sextupole errors on the MAX-IV
machine, however, it doesn’t account for harmonic sextupoles, which are widely
used on many other machines. Through generating vertical dispersion with
chromatic skew quadrupoles, a measurable dependency of nonlinear optics on
harmonic sextupoles can be observed from hybrid horizontal and vertical
dispersive orbits. Proof of concept of our expanded technique was accomplished
by simulations and beam measurements on the National Synchrotron Light Source
II (NSLS-II) storage ring.

##  I introduction

Characterizing the nonlinear optics of storage rings is becoming more
essential with the introduction of higher order multipole magnets in
accelerator design. Errors from the higher order multipoles have been observed
to degrade machine performance, such as reduction of dynamic aperture, energy
acceptance, etc. Some efforts have been made to identify the nonlinear
multipole errors by measuring distorted resonance driving terms Franchi  _et
al._ (2014), which requires a complicated Hamiltonian dynamics analysis. A
more practical technique for measuring the nonlinear optics from off-energy
closed orbits (NOECO) was reported and demonstrated on the MAX-IV ring Olsson
_et al._ (2020). Significant improvements on its dynamic aperture and beam
lifetime were observed after correcting sextupole errors. Desired results were
obtained while testing the NOECO technique on the ESRF-EBS ring as well Liuzzo
_et al._ (2022). However, the dependency of nonlinear optics on off-energy
orbits is only measurable for chromatic sextupoles. The horizontal dispersion
seen by chromatic sextupoles are usually quite large, as to effectively
correct the chromaticity. This technique, however, doesn’t apply to harmonic
sextupoles, which do not see the first order linear dispersion. Harmonic
sextupoles are used in almost every third-generation light source ring, and
some fourth-generation diffraction-limited machines, such as the ALS-U ring
Steier  _et al._ (2019). They are even being used in the design of a future
electron-ion collider ring Cai  _et al._ (2022). As such, an expansion of the
existing NOECO technique to correct for the harmonic sextupoles would be
useful due to their common, integral use in current and future accelerator
design. In the National Synchrotron Light Source II (NSLS-II) ring Dierker
(2007), the number of harmonic sextupoles are greater than the number of
chromatic sextupoles (180:90). Therefore, correcting harmonic sextupole errors
is important for improving machine performance due to their greater influence.
In this paper, we outline our expansion on the capabilities of existing
sextupole correction techniques to accommodate for the harmonic sextupoles.

A straightforward method for calibrating harmonic sextupoles for correction
would be to temporarily convert them to chromatic ones. This could be achieved
by tuning the quadrupoles inside achromats to generate a commensurate amount
of dispersion at the locations of harmonic sextupoles Olsson . However, this
method would require a significant modification of the original linear
lattice. Implementations during online measurements, such as updating the
nonlinear optics dependency for different leaked dispersion bumps would also
be complicated. Another method would be to generate local orbit bumps,
calibrated through the sextupoles, and then measuring the optics distortion
with different bump parameters. This method would not only require sufficient
beam position monitors (BPMs) that neighbor the sextupoles, but would also be
complicated to implement. In real-world applications, it is time-consuming to
form perfectly closed local bumps with orbit correctors, and then to update
the optics dependence on these bump settings Choi . While the above methods
would be capable of achieving the desired outcome, they are not practical when
considering the limitations of routine operations of user facilities.

When a sextupole sees vertical dispersion, the nonlinear optics of the off-
energy orbits will also depend on its gradient
$K_{2}=\frac{1}{(B\rho)_{0}}\frac{\partial^{2}B_{y}}{\partial x^{2}}$,
normalized with the beam rigidity $(B\rho)_{0}$. A vertical dispersive wave
can be generated through chromatic skew quadrupoles. In most light source
rings, skew quadrupoles are widely equipped to control the residual vertical
dispersion and linear coupling. Usually, a considerable amount of vertical
dispersion can be generated, but only introduces weak coupling when the
Betatron tune has sufficiently diverged from the linear difference/sum
resonance. Thus, the nonlinear off-energy optics depends on not only chromatic
sextupoles, but also on the original harmonic ones. In other words, horizontal
harmonic sextupoles are converted into vertical chromatic ones, which makes
their calibration and correction possible on hybrid dispersive orbits. In our
studies, the NSLS-II ring double-bend achromat lattice was used to demonstrate
these expanded capabilities.

The remainder of this paper is outlined as follows: Sect. II introduces the
principle of the technique in conjunction with the NSLS-II lattice. In Sect.
III we demonstrate our technique with some simulations. Some beam measurements
to calibrate both the chromatic and harmonic sextupole errors are given in
Sect. IV, with the caveat that no real sextupole correction can be implemented
at this time due to their in-series power supplies. Sect. V discusses the
hardware requirements necessary to apply this technique. A brief summary is
given in Sect. VI.

##  II Nonlinear optics on hybrid dispersive orbit

Chromatic skew quadrupoles (located at horizontally dispersive sections) can
couple the dispersion function between the horizontal and vertical planes.
This property is widely used to minimize the vertical beam size in most light
source rings. At the NSLS-II ring, each odd-numbered cell is equipped with one
$0.2\;m$ long chromatic skew quadrupole (see Fig. 1). Their maximum gradients
are $g_{1}=0.35\;T\cdot m^{-1}$, which is limited by the capacity of their
power supplies. Assuming we can double their gradients to $g_{1}=0.70\;T\cdot
m^{-1}$, a vertical dispersion wave with a $\sim 0.1\;m$ amplitude can be
generated. The necessity for a double gradient will be discussed in Sect. V.
Although these gradients are twice as large as the maximum output of their
power supplies, they are still quite weak compared to other operational
quadrupoles with a maximum gradient of $g_{1,max}=22\;T\cdot m^{-1}$. Under
these conditions, the exact coupled optics computed with the Ripken
parameterization Borchardt  _et al._ (1988); Willeke and Ripken (1989)
indicates that the linear optics remain weakly coupled. In Fig. 1, the non-
dominated functions $\beta_{1,y}$ and $\beta_{2,x}$ (dashed lines) are
observed as very close to zero, while the dominated $\beta_{1,x}$ and
$\beta_{2,y}$ (solid lines) are almost the same as in the uncoupled case. The
skew quadrupoles also cause a small amount of horizontal dispersion to be
leaked into the straight sections. Although such small residual dispersion
could not be solely used to measure the off-energy nonlinear optics, its
effect is accounted for in our method because the exact parameterization has
been used. In short, when the machine tune is configured to avoid linear
coupling resonances, chromatic skew quadrupoles can generate a considerable
amount of vertical dispersion, but only introduce relatively small linear
coupling. The newly generated vertical dispersion seen by the original
harmonic sextupoles can make the nonlinear optics on off-energy orbits (as
illustrated in Fig. 2) dependent on their gradients. Therefore, this
dependence can be utilized for their calibration and correction.

Figure 1:  Exact Ripken Twiss functions for a supercell at NSLS-II when a
vertical dispersion wave is generated with chromatic skew quadrupoles (with a
normalized gradient $K_{1,s}=\frac{1}{(B\rho)_{0}}\frac{\partial
B_{y}}{\partial x}=0.070\;m^{-2}$). The location of the skew quadrupole is
marked with a red vertical dash-dot line. Two non-dominated $\beta$-functions
(dashed lines) indicate that this optics configuration remains weakly coupled.
Figure 2:  Three-dimensional view of the hybrid dispersive orbit with
$\delta=0.5\%$ for one NSLS-II supercell.

When horizontal dispersion is seen by sextupoles in an uncoupled linear optics
configuration, the dependence of the $\beta$-function on the beam energy
deviation, $\delta=\frac{\Delta P}{P}$, and the normalized sextupole gradient,
$K_{2}$, can be formulated Wiedemann (2015). When both horizontal and vertical
dispersion can be seen by sextupoles in a weakly coupled optics configuration,
no such simple formulae are available. However, it can be numerically computed
with the two following methods. Method 1: First, for a given energy offset
$\delta$ and skew quadrupole settings $K_{1,s}$, a hybrid dispersive closed
orbit can be obtained with iterative tracking. This hybrid dispersive orbit
now has both the horizontal and vertical offsets. Then a one-turn matrix $R$
can be obtained along the dispersive orbit with a truncated power series
algorithm technique Berz (1988). From the linear components, four coupled
Ripken Twiss functions can then be extracted and propagated around the whole
ring. By slightly tweaking the settings of an arbitrary sextupole with a
$dK_{2}$, and repeating the above procedure, the dependence of
$\frac{d\beta}{d\delta}$ on $K_{2}$ can be determined. Method 2: A direct
particle tracking can be implemented with the same lattice setting as
described in the previous method. When the initial particle coordinates are
confined within the linear regime, the linear one-turn matrix $R$ can also be
fitted from turn-by-turn trajectories. Then the Ripken Twiss functions can be
parameterized. After comparing the results for one sextupole using these two
methods, which were yielded respectively with the mad-x PTC analysis Schmidt
(2005) and the particle trajectory tracking with the code elegant Borland
(2000), the results of both methods were consistent 3.

Figure 3:  Comparison of the optics dependence on one harmonic sextupole
computed with different methods. Top: the response vector in the horizontal
plane. Bottom: the response vector in the vertical plane

For demonstration purposes, we choose a harmonic sextupole “SH3” and a
chromatic sextupole “SM1” to compute their linear dependence on the off-energy
optics, i.e., the so-called response vectors, as seen below. Only two
dominated optics functions $\frac{d\beta_{1,x}}{d\delta}$ and
$\frac{d\beta_{2,y}}{d\delta}$ observed at their corresponding BPMs were
computed. If no skew quadrupoles are used to excite the beam, the optics
functions degenerate to the uncoupled $\beta_{x}$ and $\beta_{y}$. The
response vectors computed with and without the vertical dispersion are
compared in Fig. 4. With horizontal-only dispersive orbits, the dependence of
off-energy optics on “SH3(N)” is not measurable in both the horizontal and
vertical planes. On the hybrid dispersive orbits, a measurable dependence on
“SH3(Y)” can be observed. Note that, for both cases, the dependence of the
chromatic “SM1(Y/N)” is always measurable because it sees a large horizontal
dispersion. In the meantime, the dependencies are quite similar since the
optics are only slightly altered.

Figure 4:  Comparison of the off-energy optics dependence on a chromatic
sextupole “SM1” and a harmonic sextupole “SH3”. On horizontal-only dispersive
orbits, the optics dependence on the harmonic sextupole (labeled as “SH3(N)”)
is too small to measure. On hybrid dispersive orbits, a measurable dependence
(labeled as “SH3(Y)”) is observed.

In principle, chromatic and harmonic sextupole errors can be calibrated
simultaneously with a sufficiently large vertical dispersion. However,
chromatic sextupoles usually have stronger responses than harmonic ones,
particularly when the dispersion in the vertical plane is coupled from the
horizontal plane. Therefore, for practical purposes, we can uncouple the
chromatic and harmonic sextupole correction via a two-stage correction.
Stage-1: correcting chromatic sextupoles first with the existing technique
Olsson  _et al._ (2020). Stage-2: generating a vertical dispersion wave, then
calibrating harmonic sextupoles from hybrid dispersive orbits. In this paper,
we only focus on the $2^{nd}$ stage because the $1^{st}$ stage has already
been well studied with the linear optics from closed orbits (LOCO) algorithm
Safranek (1997) in the ref. Olsson  _et al._ (2020).

The $\beta$-functions and phase advances can be measured directly from turn-
by-turn data using the harmonic analysis Borer  _et al._ (1992), or the
numerical analysis of fundamental frequencies (NAFF) algorithm Laskar  _et
al._ (1992); Zisopoulos  _et al._ (2013). Therefore, instead of the LOCO
algorithm, the dependence of $\frac{d\beta}{d\delta}$ on the harmonic
sextupole settings $dK_{2}$, which were deliberately mis-set, was used in our
simulations. Given a vertical dispersion wave pattern as shown in Fig. 1, the
response matrices of $\frac{d\beta_{(1,x),(2,y)}}{d\delta}$ dependence on 180
harmonic sextupoles were computed with the lattice model and illustrated in
Fig. 5. Here, we only used two dominated $\beta_{(1,x),(2,y)}$-functions, and
the other two non-dominated ones, $\beta_{(1,y),(2,x)}$, were ignored because
they are too small to measure accurately.

Figure 5:  Given a hybrid dispersion wave in Fig. 1, the harmonic sextupole
response matrices (RM) observed at the BPMs are obtained by computing the
deviation of $\frac{d\beta}{d\delta}$ while slightly changing each individual
harmonic sextupole from its nominal setting. Left: the dominated $\beta_{1,x}$
function’s RM in the horizontal plane. Right: the dominated $\beta_{2,y}$
function’s RM in the vertical plane.

##  III Simulations

Below we simulated two specific cases to study the performance of our expanded
sextupole correction technique.

###  III.1 Case 1: two individual isolated errors

First, we studied a case in which two isolated (far apart from each other)
sextupoles errors $\Delta K_{2}=+1.5,\;-1m^{-3}$ were added onto the
$32^{th},\;154^{th}$ harmonic sextupoles. The distortions of
$\Delta\frac{d\beta}{d\delta}$ observed at the BPMs are shown with the dashed
lines in Fig. 6. The needed corrections $\Delta K_{2}$ were obtained by
solving the following linear regression problem with the response matrices
computed in the previous section,

$\Delta\frac{d\beta}{d\delta}\approx M_{x,y}\Delta K_{2},$ (1)

here, $M_{x,y}$ represents a vertically stacked matrix with the horizontal and
vertical response matrices.

Figure 6:  Comparison of nonlinear optics distortions before and after
correction in the presence of two isolated sextupole errors. The dashed lines
are the uncorrected distortions, and the solid lines represent the corrected
ones.

The correction scheme obtained with Eq. (1) for 180 harmonic sextupoles is
shown in Fig. 7. Due to the high degeneracy among the neighboring sextupoles,
the scheme doesn’t reproduce the original error distributions. Instead, they
spread to their neighbors. Nevertheless, the errors were localized, and after
applying the correction scheme, the nonlinear optics were recovered as
illustrated in Fig. 6.

Figure 7:  Correction scheme for two isolated sextupole errors obtained by the
linear regression algorithm. Despite high degeneracy among sextupoles, two
error sources are precisely localized.

###  III.2 Case 2: normally distributed errors

In this case, random distributed errors on all 180 harmonic sextupoles are
introduced and the distortion of off-energy optics are computed. Then the same
correction procedure is employed. For comparison, the optics distortions
before and after correction, and the real error distributions and computed
correction scheme are illustrated in Fig. 8 and Fig. 9), respectively.

Figure 8:  Distortions of the nonlinear optics with 180 randomly added errors.
The dashed lines are uncorrected distortions, and the solid lines represent
the distortions after correction. Figure 9:  Comparison of the added errors
(labeled as “real”) on 180 harmonic sextupoles and the obtained correction
scheme (shown with red lines labeled as “fitting”).

In both cases, as seen in Fig. 7 and 9, the obtained correction schemes only
approximately follow the real errors that were added in advance. This is due
to the strong degeneracy that exists among sextupoles in the NSLS-II lattice.
Minor imperfections of the BPMs, and other errors can also result in some
degeneracy. However, the distortion of nonlinear optics can still be well
corrected in Fig. 6 and 8. It is also worth mentioning that the dependence of
nonlinear optics on sextupoles is not purely linear, therefore, an iterative
correction might be necessary in online measurements.

###  III.3 Improvement on dynamic aperture degradation

As observed in the previous simulations, strong degeneracy among sextupoles
prevents reproducing the real error distributions accurately. It is because,
on the NSLS-II ring, every three harmonic sextupoles on the same girder are
closely assembled. In our case, what is actually corrected is the distorted
nonlinear optics dependence on beam energy deviations seen by the BPMs. The
correction scheme based on the BPM observations, therefore, might only be able
to recover the optics distortion rather than the dynamic aperture. To
illustrate this, the dynamic apertures of the ideal machine, and
uncorrected/corrected nonlinear lattices for the $2^{nd}$ simulation were
computed for comparison (Fig. 10). Although the degraded dynamic aperture due
to sextupole errors could not be fully recovered through the correction
scheme, a significant improvement was achieved. Such improvement is the main
purpose of calibrating and correcting the distorted nonlinear optics. If we
could distinguish between the degeneracy among the sextupoles, further
improvement could be made. This topic is slightly beyond the scope of this
paper, however, but worth more study.

Figure 10:  Comparison of the dynamic apertures for the ideal machine, and
uncorrected/corrected nonlinear lattices. Although the degraded dynamic
aperture (yellow line) could not be well recovered to that of the ideal
machine (blue line), significant improvement (green line) was achieved,
through correcting the distorted nonlinear optics.

##  IV Measurements

###  IV.1 Two-stage measurements

A two-stage proof-of-principle through online calibration of sextupole errors
was implemented at the NSLS-II storage ring. As the sextupoles are powered in
series, the sextupoles lack individual configurability. Therefore, no actual
nonlinear optics correction can be implemented with these limitations. For
stage-1, we calibrated 90 chromatic sextupoles with existing techniques.
First, spurious vertical dispersion was minimized using 15 chromatic skew
quadrupoles, and the global linear coupling was well corrected with another 15
non-dispersive skew quadrupoles. The $\frac{d\beta}{d\delta}$ seen by the BPMs
were measured from horizontal dispersive orbits through varying the beam
energies. By comparing the measured nonlinear optics against the design model,
the chromatic sextupole errors (red bars in Fig. 11) were obtained using the
model response matrices, and then incorporated into the lattice model. The
updated model would be used as the reference for the stage-2 calibration.

Figure 11:  Measured off-energy optics at two stages. Top: off-energy optics
based on the ideal model and online measurements from horizontal dispersive
orbits which can be used to calibrate chromatic sextupoles. Bottom: online
measurements of off-energy optics from horizontal and hybrid dispersive
orbits, which can be used to calibrate harmonic sextupoles.

For stage-2, a vertical dispersion wave was generated with 15 dispersive skew
quadrupoles to their maximum capacity. Based on the measured dispersion, the
15 skew quadrupole settings and the vertical dispersion at the BPMs were
reproduced with the lattice model as illustrated in Fig. 12. To achieve
greater accuracy, a large amplitude vertical dispersion wave is preferred.
However, it is limited by the capacity of the skew quadrupole power supply.
Under the current configuration, $\sim 0.05\;m$ is the maximum amplitude that
can be generated. The $\frac{d\beta}{d\delta}$ seen by the BPMs were re-
measured, but from hybrid dispersive orbits this time. With the updated
lattice model (incorporated with skew quadrupoles and chromatic sextupole
errors) as the new reference, 180 harmonic sextupole errors were obtained
(blue bars in Fig. 11.

Figure 12:  Measured vertical dispersion seen by the BPMs, and its
reproduction with the lattice model.

The off-energy optics ($\beta$-functions) were measured with different RF
frequencies $f_{0}+\Delta f$, i.e., different energy
$\delta=-\frac{1}{\alpha_{c}}\frac{\Delta f}{f_{0}}$, with $\alpha_{c}$ the
momentum compaction factor, and $f_{0}$ the nominal RF oscillator frequency
for on-momentum electrons. Using Eq. (1), the sextupole errors were calibrated
as illustrated in Fig. 13.

Figure 13:  Calibrated sextupole errors from the two-stage measurements. The
chromatic sextupoles marked with red bars are from the first stage, while the
harmonic ones (blue bars) are from the second stage.

###  IV.2 Validation of stage-2 measurement

Skew quadrupoles used as correctors are usually operated with dual polarity
power supplies. By flipping skew quadrupole polarities, hybrid dispersive
orbits are also flipped in the vertical plane. In Appendix Appendix: flipped
dispersive orbit in the vertical plane, we prove that the off-energy optics
dependence on sextupoles in the flipped orbit remains unchanged when the
Ripken parameterization is used. Therefore, we can repeat the stage-2
measurement on the flipped vertical dispersive orbits as a validation. In Fig.
14, two out-of-phase dispersion waves were obtained by flipping all skew
quadrupole outputs $K_{1}$ from $-0.035\;m^{-2}$ to $0.035\;m^{-2}$. On the
flipped vertical dispersive orbits, stage-2 measurements of off-energy optics
changed with respect to stage-1, were repeated as illustrated in Fig. 15. A
similar pattern can be recognized in two independently measured optics
distortions. This indicates that the dependence of off-energy optics on
harmonic sextupoles is measurable on the hybrid dispersive orbit, although the
precision is limited by the low capacity of skew quadrupole power supplies.

Figure 14:  Out-of-phase vertical dispersion waves were generated with flipped
skew quadrupole polarities. Figure 15:  Similar off-energy optics distortions
were observed when the hybrid dispersive orbit was flipped in the vertical
plane.

##  V Requirements on power supplies of magnets

In this section, we estimate the requirements of the power supplies for the
skew quadrupoles and sextupoles to make this calibration and correction
practical, specifically on the NSLS-II ring. First, sextupoles need to be
powered independently, or equipped with back-leg windings to allow individual
corrections. Second, the skew quadrupoles that generate the vertical
dispersion wave should be sufficiently strong for better resolution. Next, we
use the NSLS-II lattice to estimate the required skew quadrupole strengths.

The measurement accuracy of $\beta$-functions using turn-by-turn data was
found at the level of about $0.03\;m$ on the NSLS-II storage ring. If we
adjust the RF frequency by $\pm 1,000\;Hz$, the corresponding energy deviation
is $\Delta\delta=-\frac{1}{\alpha_{c}}\frac{\Delta f}{f_{0}}\approx 1.0\%$
with $\alpha_{c}$ momentum compaction factor. To resolve a single harmonic
sextupole error to the level of 1 unit $\Delta K_{2}=1\;m^{-3}$, the magnitude
of $\frac{d}{dK_{2}}\frac{d\beta}{d\delta}$ needs to be greater than
$3\;m^{4}$. To generate such a strong dependence, the amplitude of the
vertical dispersion wave is required to be greater than $0.1\;m$. In order to
couple $0.1\;m$ dispersion to the vertical plane, $K_{1,sm}\geq 0.07\;m^{-2}$
must be applied to all 15 skew quadrupoles. Therefore, strong skew quadrupoles
were chosen to identify two isolated errors in our simulations. Although a
$\frac{\partial B_{x}}{\partial x}=0.70\;T\cdot m^{-1}$ quadrupolar gradient
is still quite weak, it already exceeds the capacity of our skew quadrupole
power supplies $g_{2,max}=0.35\;T\cdot m^{-1}$. In other words, stronger skew
quadrupoles are needed to generate larger vertical dispersion for better
resolution. As all 270 sextupoles have some associated errors, the accumulated
magnitude of $\Delta\frac{d\beta}{d\delta}$ are at the level of 100-200
meters, which allows us to calibrate the approximate error distribution in
Sect. IV.

The above estimation doesn’t consider even higher orders of nonlinear optics
from off-energy orbits. Once the higher order nonlinearities
$\frac{d^{n}\beta}{d\delta^{n}}$ with $n\geq 2$ appear, we couldn’t improve
the measurement accuracy by increasing beam energy off-sets. On the other
hand, increasing skew quadrupole strengths to generate higher vertical
dispersive orbits can significantly improve the sensitivity of off-energy
optics to sextupole settings without introducing too much nonlinearity.
Therefore, having sufficiently strong skew quadrupoles should be considered a
necessary condition for this technique.

##  VI Summary

We expanded the capability of the technique for measuring nonlinear optics
distortions from off-energy orbits to account for the harmonic sextupole
contribution. Using hybrid dispersive off-energy nonlinear optics, the errors
of the harmonic sextupoles can be measured. The corresponding correction can
be more effectively implemented if they are independently configurable. A
practical benefit of our expanded method is that a considerable amount of
vertical dispersion can be generated with weak skew quadrupoles. Meanwhile,
because the original lattice is already weakly coupled, its optics properties
can still be well maintained. Thus far, only sextupole calibration was
considered in our studies, and higher order nonlinear magnets, such as
octupoles, we have not yet investigated. This technique might be applicable if
their contributions were sufficiently strong.

###### Acknowledgements.

We would like to thank the collaborative and productive discussion with Dr. D.
Olsson (MAX IV, Lund Uni.), Prof. Y. Hao (MSU), and some NSLS-II colleagues,
Dr. J. Choi, Dr. Y. Hidaka, Dr. M. Song, Dr. G. Tiwari, Dr. X. Yang, et al.
This research used resources of the National Synchrotron Light Source II, a
U.S. Department of Energy (DOE) Office of Science User Facility operated for
the DOE Office of Science by Brookhaven National Laboratory under Contract No.
DE-SC0012704.

## Appendix: flipped dispersive orbit in the vertical plane

The hybrid dispersive orbit can be flipped only in the vertical plane by
changing the skew quadrupole polarities. Usually skew quadrupoles used for the
correction purposes are operated with dual polarity power supplies. In this
appendix, we prove that an exact flipping of the vertical dispersive orbit
doesn’t change the dependence of off-energy optics on sextupoles which can be
used to validate online measurements.

###  VI.1 Transfer matrix with single skew quadrupole

The $4\times 4$ transfer matrix of a normal quadrupole reads as

$T_{n}=\left[\begin{matrix}T_{x}&\mathbf{0}\\\
\mathbf{0}&T_{y}\end{matrix}\right]$ (2)

with two $2\times 2$ zero blocks $\mathbf{0}$ as its off-diagonal matrices.

After rotating it by $\frac{\pi}{4}$ around the longitudinal axis, it becomes
a skew with a transfer matrix

$T_{s}=R_{z}\left(\frac{\pi}{4}\right)T_{n}R_{z}\left(-\frac{\pi}{4}\right)=\left[\begin{matrix}\frac{T_{x}+T_{y}}{2}&\frac{-T_{x}+T_{y}}{2}\\\
\frac{-T_{x}+T_{y}}{2}&\frac{T_{x}+T_{y}}{2}\end{matrix}\right]$ (3)

where $R_{z}\left(\theta\right)$ represents the transfer matrix of a pure
rotation with an angle of $\theta$.

Assuming there is only one skew quadrupole inside a periodic lattice cell with
the layout “normal section 1 – skew quad – normal section 2”, the transfer
matrix of the whole section is

$M_{t}=\left[\begin{matrix}M_{x,2}&\mathbf{0}\\\
\mathbf{0}&M_{y,2}\end{matrix}\right]T_{s}\left[\begin{matrix}M_{x,1}&\mathbf{0}\\\
\mathbf{0}&M_{y,1}\end{matrix}\right]=\left[\begin{matrix}A&B\\\
C&D\end{matrix}\right]$ (4)

When the skew quadrupole polarity is flipped from $K_{1,s}$ to $-K_{1,s}$, the
new transfer matrix can be obtained by swapping $T_{x}$ and $T_{y}$ in Eq.
(3),

$T_{s}^{\prime}=\left[\begin{matrix}\frac{T_{y}+T_{x}}{2}&\frac{-T_{y}+T_{x}}{2}\\\
\frac{-T_{y}+T_{x}}{2}&\frac{T_{y}+T_{x}}{2}\end{matrix}\right]=ST_{s}S^{-1}$
(5)

with

$S=\left[\begin{matrix}\mathbf{1}&\mathbf{0}\\\
\mathbf{0}&-\mathbf{1}\end{matrix}\right]$ (6)

where $\mathbf{1}$ is the $2\times 2$ identity matrix. The transfer matrix of
the whole section then becomes

$\displaystyle M_{t}^{\prime}$
$\displaystyle=\left[\begin{matrix}M_{x,2}&\mathbf{0}\\\
\mathbf{0}&M_{y,2}\end{matrix}\right]T_{s}^{\prime}\left[\begin{matrix}M_{x,1}&\mathbf{0}\\\
\mathbf{0}&M_{y,1}\end{matrix}\right]=S\left(S^{-1}\left[\begin{matrix}M_{x,2}&\mathbf{0}\\\
\mathbf{0}&M_{y,2}\end{matrix}\right]S\right)T_{s}\left(S^{-1}\left[\begin{matrix}M_{x,1}&\mathbf{0}\\\
\mathbf{0}&M_{y,1}\end{matrix}\right]S\right)S^{-1}$ (7)
$\displaystyle=S\left[\begin{matrix}A&B\\\
C&D\end{matrix}\right]S^{-1}=\left[\begin{matrix}A&-B\\\
-C&D\end{matrix}\right]$

By comparing Eq. (4) and Eq. (7), we can conclude that when a periodic cell
contains only one skew quadrupole, only the signs of the off-diagonal blocks
$B$ and $C$ flip when the skew quadrupole polarity is flipped.

###  VI.2 Ripken Twiss functions with single flipped skew

Using the Ripken parameterizationWilleke and Ripken (1989); Lebedev and Bogacz
(2010) on $M_{t}$, the transfer matrix is diagonalized respectively as

$M_{t}=UVU^{-1},$ (8)

with

$V=\left[\begin{matrix}\cos\mu_{1}&\sin\mu_{1}&0&0\\\
-\sin\mu_{1}&\cos\mu_{1}&0&0\\\ 0&0&\cos\mu_{2}&\sin\mu_{2}\\\
0&0&-\sin\mu_{2}&\cos\mu_{2}\end{matrix}\right],U=\left[\begin{matrix}\sqrt{\beta_{1x}}&0&\sqrt{\beta_{2x}}\cos
v_{2}&-\sqrt{\beta_{2x}}\sin v_{2}\\\
-\frac{\alpha_{1x}}{\sqrt{\beta_{1x}}}&\frac{1-u}{\beta_{1x}}&\frac{-\alpha_{2x}\cos
v_{2}+u\sin v_{2}}{\sqrt{\beta_{2x}}}&\frac{\alpha_{2x}\sin v_{2}+u\cos
v_{2}}{\sqrt{\beta_{2x}}}\\\ \sqrt{\beta_{1y}}\cos
v_{1}&-\sqrt{\beta_{1y}}\sin v_{1}&\sqrt{\beta_{2y}}&0\\\
\frac{-\alpha_{1y}\cos v_{1}+u\sin
v_{1}}{\sqrt{\beta_{1y}}}&\frac{\alpha_{1y}\sin v_{1}+u\cos
v_{1}}{\sqrt{\beta_{1y}}}&-\frac{\alpha_{2y}}{\sqrt{\beta_{2y}}}&\frac{1-u}{\sqrt{\beta_{2y}}}\end{matrix}\right]$
(9)

where $\beta,\alpha,\mu,\nu,u$ are defined in Lebedev and Bogacz (2010).

Now we implement the same parameterization on $M_{t}^{\prime}$ with a flipped
skew quadrupole polarity,

$\displaystyle M_{t}^{\prime}$ $\displaystyle=SM_{t}S^{-1}=SUVU^{-1}S^{-1}$
(10)
$\displaystyle=\left(SUS^{-1}\right)\left(SVS^{-1}\right)\left(SUS^{-1}\right)^{-1}$
$\displaystyle=U^{\prime}V^{\prime}U^{\prime-1}$

with

$V^{\prime}=V$ (11)

and

$U^{\prime}=SUS^{-1}=\left[\begin{matrix}\sqrt{\beta_{1x}}&0&-\sqrt{\beta_{2x}}\cos
v_{2}&\sqrt{\beta_{2x}}\sin v_{2}\\\
-\frac{\alpha_{1x}}{\sqrt{\beta_{1x}}}&\frac{1-u}{\beta_{1x}}&\frac{\alpha_{2x}\cos
v_{2}-u\sin v_{2}}{\sqrt{\beta_{2x}}}&-\frac{\alpha_{2x}\sin v_{2}+u\cos
v_{2}}{\sqrt{\beta_{2x}}}\\\ -\sqrt{\beta_{1y}}\cos
v_{1}&\sqrt{\beta_{1y}}\sin v_{1}&\sqrt{\beta_{2y}}&0\\\ \frac{\alpha_{1y}\cos
v_{1}-u\sin v_{1}}{\sqrt{\beta_{1y}}}&-\frac{\alpha_{1y}\sin v_{1}+u\cos
v_{1}}{\sqrt{\beta_{1y}}}&-\frac{\alpha_{2y}}{\sqrt{\beta_{2y}}}&\frac{1-u}{\sqrt{\beta_{2y}}}\end{matrix}\right]$
(12)

By comparing Eq. (8) and Eq. (12), we get

$\begin{gathered}\mu_{1}^{\prime}=\mu_{1},\qquad\mu_{2}^{\prime}=\mu_{2},\qquad
u^{\prime}=u,\qquad v_{1}^{\prime}=v_{1}+\pi,\qquad
v_{2}^{\prime}=v_{2}+\pi\\\
\beta_{1x}^{\prime}=\beta_{1x},\qquad\beta_{1y}^{\prime}=\beta_{1y},\qquad\beta_{2x}^{\prime}=\beta_{2x},\qquad\beta_{2y}^{\prime}=\beta_{2y},\\\
\alpha_{1x}^{\prime}=\alpha_{1x},\qquad\alpha_{1y}^{\prime}=\alpha_{1y},\qquad\alpha_{2x}^{\prime}=\alpha_{2x},\qquad\alpha_{2y}^{\prime}=\alpha_{2y}\end{gathered}$
(13)

Therefore, flipping the skew quadrupole polarity doesn’t change the signs or
values of $\beta$ and $\alpha$ functions in the Ripken parameterization.

###  VI.3 Off-energy optics dependence on sextupoles with flipped dispersive
orbit

Now we consider the off-energy optics dependence on sextupoles with the
flipped dispersive orbit. When an off-energy particle passes through a
sextupole with a vertical offset $y_{0}$, it sees a skew quadrupolar
component, which is proportional to $K_{2}Ly_{0}=K_{2}L\eta_{y}\delta$ with
$K_{2}L$ as the sextupole’s effective field integral, and $\delta$ as the
particle momentum offset. On the flipped vertical dispersion orbit, the one-
turn transfer matrix, which is composed of a sequence of matrices (each of
them has only one sextupole included)

$M=\left[\begin{matrix}A&B\\\
C&D\end{matrix}\right]=\left[\begin{matrix}A_{1}&B_{1}\\\
C_{1}&D_{1}\end{matrix}\right]\left[\begin{matrix}A_{2}&B_{2}\\\
C_{2}&D_{2}\end{matrix}\right]\cdots\left[\begin{matrix}A_{n}&B_{n}\\\
C_{n}&D_{n}\end{matrix}\right]$ (14)

Here each section’s transfer matrix $\left[\begin{matrix}A_{i}&B_{i}\\\
C_{i}&D_{i}\end{matrix}\right]$ includes only one skew quadrupole or
sextupole. After flipping vertical dispersion with skew quadrupoles, the one-
turn matrix on the closed dispersive orbit with $-y_{0}$ can be obtained by
following the same rule as in the presence of one skew quadrupole,

$\displaystyle M^{\prime}$ $\displaystyle=S\left[\begin{matrix}A_{1}&B_{1}\\\
C_{1}&D_{1}\end{matrix}\right]S^{-1}S\left[\begin{matrix}A_{2}&B_{2}\\\
C_{2}&D_{2}\end{matrix}\right]S^{-1}\cdots S\left[\begin{matrix}A_{n}&B_{n}\\\
C_{n}&D_{n}\end{matrix}\right]S^{-1}$ (15)
$\displaystyle=S\left[\begin{matrix}A&B\\\
C&D\end{matrix}\right]S^{-1}=\left[\begin{matrix}A&-B\\\
-C&D\end{matrix}\right]$

Therefore, flipping the vertical dispersion doesn’t change the sign or value
of coupled Ripken Twiss functions, and neither does the off-energy optics
dependence on sextupole strength
$\frac{1}{\partial{K_{2}}}\frac{\partial{\beta}}{\partial\delta}$. This
property has been numerically confirmed with the mad-x and elegant
computations. It can also be used to validate online measurements.

## References

  * Olsson  _et al._ (2020) D. Olsson, Å. Andersson, and M. Sjöström, “Nonlinear optics from off-energy closed orbits,” Physical Review Accelerators and Beams 23, 102803 (2020).
  * Franchi  _et al._ (2014) A. Franchi, L. Farvacque, F. Ewald, G. Le Bec, and K.B. Scheidt, “First simultaneous measurement of sextupolar and octupolar resonance driving terms in a circular accelerator from turn-by-turn beam position monitor data,” Physical Review Special Topics-Accelerators and Beams 17, 074001 (2014).
  * Liuzzo  _et al._ (2022) S.M. Liuzzo, N. Carmignani, L.R. Carver, L. Hoummi, T. Perron, B. Roche, and S. White, “Lifetime correction using fast-off-energy response matrix measurements,” in  _13th International Particle Accelerator Conference_ (2022) p. TUPOMS008.
  * Steier  _et al._ (2019) C. Steier, P. Amstutz, K. Baptiste, P. Bong, E. Buice, P. Casey, K. Chow, R. Donahue, M. Ehrlichman, J. Harkins, _et al._ , “Design progress of ALS-U, the soft x-ray diffraction limited upgrade of the advanced light source,” in  _10th Int. Particle Accelerator Conf.(IPAC’19), Melbourne, Australia_ (2019).
  * Cai  _et al._ (2022) Y. Cai, Y. Nosochkov, S. Berg, J. Kewisch, Y. Li, D. Marx, C. Montag, S. Tepikian, F. Willeke, G. Hoffstaetter, _et al._ , “Optimization of chromatic optics in the electron storage ring of the electron-ion collider,” Physical Review Accelerators and Beams 25, 071001 (2022).
  * Dierker (2007) S. Dierker,  _NSLS-II preliminary design report_ , Tech. Rep. (Brookhaven National Lab, Upton, NY, USA, 2007).
  * (7) D. Olsson, personal communication.
  * (8) J. Choi, personal communication.
  * Borchardt  _et al._ (1988) I. Borchardt, E. Karantzoulis, H. Mais, and G. Ripken, “Calculation of beam envelopes in storage rings and transport systems in the presence of transverse space charge effects and coupling,” Zeitschrift für Physik C Particles and Fields 39, 339–349 (1988).
  * Willeke and Ripken (1989) F. Willeke and G. Ripken, “Methods of beam optics,” in  _AIP Conference Proceedings_ , Vol. 184 (American Institute of Physics, 1989) pp. 758–819.
  * Wiedemann (2015) H. Wiedemann,  _Particle accelerator physics_ (Springer Nature, 2015).
  * Berz (1988) A. Berz, “Differential algebraic description of beam dynamics to very high orders,” Part. Accel. 24, 109–124 (1988).
  * Schmidt (2005) F. Schmidt, “MAD-X PTC integration,” in  _Proceedings of the 2005 Particle Accelerator Conference_ (IEEE, 2005) pp. 1272–1274.
  * Borland (2000) M. Borland, “elegant: A flexible SDDS-compliant code for accelerator simulation,” Advanced Photon Source, LS-287 (2000).
  * Safranek (1997) J. Safranek, “Experimental determination of storage ring optics using orbit response measurements,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 388, 27–36 (1997).
  * Borer  _et al._ (1992) J. Borer, C. Bovet, A. Burns, and G. Morpurgo, “Harmonic analysis of coherent bunch oscillations in lep,” Conf. Proc. C, 3rd European Particle Accelerator Conference (EPAC 92) 920324, 1082–1084 (1992).
  * Laskar  _et al._ (1992) J. Laskar, C. Froeschlé, and A. Celletti, “The measure of chaos by the numerical analysis of the fundamental frequencies. application to the standard mapping,” Physica D: Nonlinear Phenomena 56, 253–269 (1992).
  * Zisopoulos  _et al._ (2013) P. Zisopoulos, Y. Papaphilippou, A. Streun, and V. Ziemann, “Beam Optics Measurements through Turn by turn Beam Position Data in the SLS,” in  _4th International Particle Accelerator Conference_ (2013) p. WEPEA067.
  * Lebedev and Bogacz (2010) V. Lebedev and S. Bogacz, “Betatron motion with coupling of horizontal and vertical degrees of freedom,” Journal of Instrumentation 5, P10010 (2010).


# On Image Registration and Subpixel Estimation

Serap A. Savari
Texas A&M University
College Station, TX 77843-3128, USA

###### Abstract

Image registration is a classical problem in machine vision which seeks
methods to align discrete images of the same scene to subpixel accuracy in
general situations. As with all estimation problems, the underlying difficulty
is the partial information available about the ground truth. We consider a
basic and idealized one-dimensional image registration problem motivated by
questions about measurement and about quantization, and we demonstrate that
the extent to which subinterval/subpixel inferences can be made in this
setting depends on a type of complexity associated with the function of
interest, the relationship between the function and the pixel size, and the
number of distinct sampling count observations available.

##  1 Introduction

The correlation of different discrete images of a setting is an extensively
studied problem for fifty years [1, 2] with diverse applications in numerous
fields. The review article Ref. 2 highlights three main objectives for these
investigations. First, there is an interest in similarity measures between
images. Second, one may wish to match up a pair of images to achieve global
registration. Third, there is sometimes an interest in displacement or
translation estimation. The main focus of Ref. 2 is on applications of the
Fourier-Mellin transform. However, other techniques like wavelets have also
been considered [3].

This paper is motivated in part by literature in different areas. There is a
precise and well-known theorem from digital communications stating that about
$2\tau W$ sample points are needed to provide a good approximation to a very
long segment of duration $\tau$ of a waveform bandlimited to frequencies with
magnitude at most $W$. [4] Therefore, for that problem the complexity of the
function of interest is important in planning how much information is needed
to understand it. Returning to image registration, the authors of [2, p. 4074]
acknowledge that “different tasks and scenarios have different requirements.
In addition to the methods, further improvements should also concentrate on
developing advanced task-oriented matching frameworks.” In particular, data
collection and computational considerations are critical in most image
registration applications. However, there does not appear to be any theory
analogous to Ref. 4 even for the simplest examples of image registration. The
goal of this paper is to begin to develop such a theory in an idealized but
practically motivated situation. Like Ref. 4, the focus here will be on a one-
dimensional problem, but the signal will not be bandlimited. As we will see,
there are constraints on learning from observations even in the absence of
computational considerations and complications like noise.

The need for measurement dates at least as far back as the early Phoenicians
trading around the Mediterranean Sea [5, p. 1]. Suppose our measuring
instrument has graduated marks separated by a standard length $T$, and we wish
to determine the length of a line segment longer than $T$. Further suppose
that we cannot precisely control where the ends of the line segment fall
relative to the graduated marks of the measuring device; for example, the
images captured from the scanning electron microscope are subject to drift.
[6]. Let us also assume that there is no ambiguity in the binary decision of
whether or not a mark meets the line segment. If we have a single observation
that $n$ marks of the device meet the feature, then we can only infer that the
length of the line segment is at least $(n-1)T$ and less than $(n+1)T$.
However, assuming that repeated measurements do not alter the length of the
line segment or its remaining in its entirety within the extent of the
measuring device, then a subsequent measurement could have either one more or
one less mark meeting the feature. In this case, the binary images with
different numbers of marks could be aligned up to a partial order, and the
uncertainty about the length of the line segment would decrease from $2T$ to
$T$. While we may be interested in also measuring objects that are shorter
than $T$, we would need a more accurate instrument to reliably do so.
Furthermore, additional observations could not decrease this level of
uncertainty in the absence of other prior information about the line segment
and/or the measurement process. However, next suppose that within the extent
of our measuring device we have two stationary line segments each of length at
least $T$ separated by a distance which is also at least $T$ long. Now a
single observation provides not only an uncertainty of $2T$ for the lengths of
each line segment and the space between them but also a $2T$ uncertainty about
the length of a line segment together with the space separating it from the
other line segment and a $2T$ uncertainty about the total length from the
beginning of one line segment to the end of the other one. As we will describe
in Section 2, for this example we can capture at most four different
observations of sampling patterns for the line segment and the space between
them which can be aligned to within a partial ordering of marks. In this case
the four distinct patterns would reveal relationships about subinterval
information up to one of six partitions of a unit cube. A need for greater
accuracy would require different types of prior information.

Quantization has been part of signal processing practice for decades [7, 8]
and is a different motivation to study the results of ideal sampling at
intervals of length $T$ on piecewise constant functions. Digital images are
generally quantized [9, §2.4] and have a quantization error relative to the
original analog signal even in the absence of other noise or distortions.
Because information about the original analog signal is lost during
quantization, it may be more realistic to view the quantized signal as the
ground truth, and that is the approach of this work. Unlike the interpolation
schemes discussed in Ref. 2, we assume that the pieces of our piecewise
constant function each have length at least $T$. For ease of alignment we also
focus here on spatially-limited signals. Since our primary focus is on the
effect of displacement, we will assume ideal sampling with no noise or other
distortions; our goal is to infer as much as we can about the lengths of each
piece of the signal.

The plan for the rest of the paper is as follows. In Section 2, we introduce
most of our notation and build intuition by studying the simplest examples in
detail. In Section 3, we descibe the possible observations corresponding to a
general spatially-limited piecewise constant function (with piece lengths at
least $T$ each) and discuss how the inverse problem relates to partitions of a
unit hypercube. In Section 4, we make some prefatory remarks about
generalizing to noisy observations and conclude.

##  2 Preliminaries

Assume a spatially-limited piecewise constant function $g(t)$ with $m$
contiguous pieces or regions of lengths $R_{1},\ R_{2},\ \dots,\ R_{m}$ with
respective values $g_{1},\ g_{2},\ \dots,g_{m}$ defined by

$g(t)\;=\;\left\\{\begin{array}[]{ll}g_{1},&0\leq t<R_{1}\\\
g_{i},&\sum_{j=1}^{i-1}R_{j}\leq t<\sum_{j=1}^{i}R_{j},\;i\geq 2\\\
0,&\mbox{otherwise};\end{array}\right.$

here $g_{1}\neq 0,\;g_{m}\neq 0,$ and $g_{i}\neq g_{i+1},i\in\\{1,\ \dots,\
m-1\\}$. The sampling interval is known to be of length $T$, and we represent
region length $R_{i}$ as

$R_{i}=(n_{i}-f_{i})T,$

where for each $i,\;n_{i}$ is an integer that is at least two and $0<f_{i}\leq
1$; at least one sample is taken from each region. The exact number $\eta_{i}$
of samples from region $i$ depends on the placement of the sampling grid.
Define the offset $\Delta_{i}\in[0,\ T)$ as the distance from the left
endpoint of region $i$ to the first sampling point within the region. We have
the following basic result.

Proposition 1: The number $\eta_{i}$ of samples from region $i,\;i\in\\{1,\
\dots,\ m\\}$, satisfies

$\eta_{i}\;=\;\left\\{\begin{array}[]{ll}n_{i},&0\leq\Delta_{i}<(1-f_{i})T\\\
n_{i}-1,&(1-f_{i})T\leq\Delta_{i}<T.\end{array}\right.$

Proof: If we start counting samples from a sample point at distance
$\Delta_{i}$ from the left endpoint of region $i$, then the location of sample
$n_{i}$ occurs at $\Delta_{i}+(n_{i}-1)T$ while the length of $R_{i}$ is
$(n_{i}-f_{i})T$ with $0<f_{i}\leq 1.\;\Box$

Proposition 1 indicates that if within our set ${\cal S}$ of observations of
the sample counts
$(\eta_{1},\eta_{2},\ \dots,\eta_{m})$ within the support of $g(t)$ we have at
least one with $\eta_{i}=n_{i}$ and one with $\eta_{i}=n_{i}-1$, then we know
the value of $n_{i}$. When $R_{i}$ is an integer multiple of $T,\;\eta_{i}$
can only take on one value and $n_{i}$ would be unclear without other prior
information. Proposition 1 also demonstrates that studying
$\eta_{i},\;i\in\\{1,\ \dots,\ m\\}$, in isolation involves a subdivision of
offset values $\Delta_{i}$ in $[0,\ T)$ into at most two parts with a boundary
depending on $f_{i}$. We will see further subdivisions of offset intervals
when we study multiple regions together.

When $m\geq 2$, the set ${\cal S}$ may offer certain types of subinterval
information based on Propostion 1 and the following relationships among
successive offsets.

Lemma 2: Suppose $m\geq 2$, The offsets $\Delta_{i},\;i\in\\{2,\ \dots,\
m\\}$, satisfy

$\Delta_{i}=(\Delta_{i-1}+f_{i-1}T)\bmod T.$

Proof: For convenience, define $R_{0}=0$. For $i\geq 2$, the location of the
final sample point of region $i-1$ is at

$\sum_{j=0}^{i-2}R_{j}+\Delta_{i-1}+(\eta_{i-1}-1)T,$

so the location of the first sample point of region $i$ is at

$\sum_{j=0}^{i-2}R_{j}+\Delta_{i-1}+\eta_{i-1}T.$

Therefore,

$\displaystyle\Delta_{i}$ $\displaystyle=$
$\displaystyle\left(\sum_{j=0}^{i-2}R_{j}+\Delta_{i-1}+\eta_{i-1}T\right)-\sum_{j=0}^{i-1}R_{j}$
$\displaystyle=$ $\displaystyle(\Delta_{i-1}+(\eta_{i-1}-n_{i-1})T+f_{i-1}T),$

which implies the result. $\Box$

We next show

Proposition 3: When $m\geq 2$, the pair $(\eta_{i},\ \eta_{i+1}),\;i\in\\{1,\
\dots,\ m-1\\}$, of sample counts from regions $i$ and $i+1$ satisfies one of
three cases:
Case 1: If $0<f_{i}<1$ and $f_{i}+f_{i+1}\leq 1$, then

$(\eta_{i},\ \eta_{i+1})\;=\;\left\\{\begin{array}[]{ll}(n_{i},\
n_{i+1}),&0\leq\Delta_{i}<(1-f_{i}-f_{i+1})T\\\ (n_{i},\
n_{i+1}-1),&(1-f_{i}-f_{i+1})T\leq\Delta_{i}<(1-f_{i})T\\\ (n_{i}-1,\
n_{i+1}),&(1-f_{i})T\leq\Delta_{i}<T;\end{array}\right.$

Case 2: If $0<f_{i}<1$ and $f_{i}+f_{i+1}\geq 1$, then

$(\eta_{i},\ \eta_{i+1})\;=\;\left\\{\begin{array}[]{ll}(n_{i},\
n_{i+1}-1),&0\leq\Delta_{i}<(1-f_{i})T\\\ (n_{i}-1,\
n_{i+1}),&(1-f_{i})T\leq\Delta_{i}<(2-f_{i}-f_{i+1})T\\\ (n_{i}-1,\
n_{i+1}-1),&(2-f_{i}-f_{i+1})T\leq\Delta_{i}<T;\end{array}\right.$

Case 3: If $f_{i}=1$, then

$(\eta_{i},\ \eta_{i+1})\;=\;\left\\{\begin{array}[]{ll}(n_{i}-1,\
n_{i+1}),&0\leq\Delta_{i}<(1-f_{i+1})T\\\ (n_{i}-1,\
n_{i+1}-1),&(1-f_{i+1})T\leq\Delta_{i}<T.\end{array}\right.$

Proof: For all three cases the values of $\eta_{i}$ immediately follow from
Proposition 1, For Case 1, Lemma 2 implies that when
$(1-f_{i}-f_{i+1})T\leq\Delta_{i}<(1-f_{i})T$ we have
$(1-f_{i+1})T\leq\Delta_{i}<T.$ Similarly, for Case 2, Lemma 2 implies that
when $(1-f_{i})T\leq\Delta_{i}<(2-f_{i}-f_{i+1})T$, we have
$0\leq\Delta_{i+1}<(1-f_{i+1})T.$ For Case 3, $\Delta_{i+1}=\Delta_{i}$. The
values of $\eta_{i+1}$ for all three cases follow from Proposition 1. $\Box$

Proposition 3 offers a subinterval relationship between region lengths $R_{i}$
and $R_{i+1}$ when there are three distinct sample count pairs $(\eta_{i},\
\eta_{i+1})$. If this is not available, then there will be ambiguity about
$n_{i},\ n_{i+1}$ and/or $f_{i}+f_{i+1}$ that can only be resolved with other
prior information.

The next result is key to understanding the possible sample counts for a set
of $K$ contiguous regions. Let $\lceil x\rceil$ denote the least integer
greater than or equal to $x$. We have

Proposition 4: The projection $(\eta_{i},\ \eta_{i+1},\ \dots,\
\eta_{i+k-1}),\;i\geq 1,\ i+k-1\leq m$, of $(\eta_{1},\ \dots,\ \eta_{m})$
will be defined by at most $k+1$ $\Delta_{i}$ offset intervals $[0,\
\Delta_{i}^{(1)}),\;[\Delta_{i}^{(1)},\ \Delta_{i}^{(2)}),\
\dots,[\Delta_{i}^{(k)},\ T),$ where the values of $\Delta_{i}^{(1)},\
\Delta_{i}^{(2)},\ \dots,\\\ \Delta_{i}^{(k)}$ are

$(1-f_{i})T,\;(\lceil f_{i}+f_{i+1}\rceil-(f_{i}+f_{i+1})T,\ \dots,\
\left(\lceil\sum_{j=0}^{k}f_{i+j}\rceil-\sum_{j=0}^{k}f_{i+j}\right)T$

sorted in ascending order.

Proof (by induction): Propositions 1 and 3 establish the result for $k=1$ and
$k=2$. Suppose it is true at $k=K$ and we wish to show that it continues to
hold at $k=K+1.$ By Lemma 2, the offset $\tilde{\Delta}_{i+K}=(1-f_{K})T$ for
region $i+K$ corresponds to the offset
$\tilde{\Delta}_{i}=\left(\lceil\sum_{j=0}^{K}f_{i+j}\rceil-\sum_{j=0}^{K}f_{i+j}\right)T$
for region $i$. For $k=K$, there is one offset interval, say
$[\Delta_{i}^{(h)},\ \Delta_{i}^{(h+1)})$ which contains $\tilde{\Delta}_{i}$.
For $k=K+1$ we will divide $[\Delta_{i}^{(h)},\ \Delta_{i}^{(h+1))}$ into
$[\Delta_{i}^{(h)},\ \tilde{\Delta}_{i})$ and $[\tilde{\Delta}_{i},\
\Delta_{i}^{(h+1)})$ and leave the other intervals as they are. This
collection of intervals meets the requirements on endpoints and cardinality.
If for $k=K$ the interval $[\Delta_{i}^{(h)},\ \Delta_{i}^{(h+1)})$
corresponds to the sample counts $(\eta_{i}^{(h)},\ \eta_{i+1}^{(h)},\ \dots,\
\eta_{i+K-1}^{(h)}),$ then for $k=K+1$ the intervals for $[\Delta_{i}^{(h)},\
\tilde{\Delta}_{i})$ and $[\tilde{\Delta}_{i},\ \Delta_{i}^{(h+1)})$
correspond to the sample counts $(\eta_{i}^{(h)},\ \eta_{i+1}^{(h)},\ \dots,\
\eta_{i+K-1}^{(h)},\ n_{i+K}),$ and
$(\eta_{i}^{(h)},\ \eta_{i+1}^{(h)},\ \dots,\ \eta_{i+K-1}^{(h)},\
n_{i+K}-1)$. To show that the remaining offset intervals for $k=K$ correspond
to projections with a unique extension at $k=K+1$, it suffices to show that
the offset $\tilde{\Delta}_{i+K}=0$ for region $i+K$ corresponds to one of the
region $i$ offset endpoints $\Delta_{i}^{(1)},\ \Delta_{i}^{(2)},\
\dots,\Delta_{i}^{(K)}$. By Lemma 2,

$\Delta_{i+K}=(\Delta_{i}+\left(\sum_{j=0}^{K-1}f_{i+j}\right)T)\bmod T.$

By the inductive hypothesis, one of $\\{\Delta_{i}^{(1)},\ \Delta_{i}^{(2)},\
\dots,\Delta_{i}^{(K)}\\}$ is

$\left(\lceil\sum_{j=0}^{K-1}f_{i+j}\rceil-\sum_{j=0}^{K-1}f_{i+j}\right)T,$

which corresponds to ${\Delta}_{i+K}=0\bmod T.\;\Box$

If we omit the special cases of regions or clusters of contiguous regions with
length that is an integer multiple of $T$, then for a spatially-limited,
piecewise constant function with $m$ regions there will be $m+1$ possible
distinct observation vectors of sample counts of the support of the function.
If all of them are available, then the subinterval relationships can be
determined up to one of the $m!$ partitions of a unit hypercube. When $m=3$,
the six sets of four observation vectors and the corresponding portion of the
unit cube are as follows:

  * •

$\\{(n_{1},\ n_{2},\ n_{3}),\ (n_{1},\ n_{2},\ n_{3}-1),\ (n_{1},\ n_{2}-1,\
n_{3}),\ (n_{1}-1,\ n_{2},\ n_{3})\\}$
$f_{1}+f_{2}+f_{3}<1.$

  * •

$\\{(n_{1},\ n_{2},\ n_{3}-1),\ (n_{1},\ n_{2}-1,\ n_{3}),\ (n_{1}-1,\ n_{2},\
n_{3}),(n_{1}-1,\ n_{2},\ n_{3}-1)\\}$
$f_{1}+f_{2}<1,\;f_{2}+f_{3}<1,\;f_{1}+f_{2}+f_{3}>1.$

  * •

$\\{(n_{1},\ n_{2},\ n_{3}-1),\ (n_{1},\ n_{2}-1,\ n_{3}),\ (n_{1},\ n_{2}-1,\
n_{3}-1),\\\ (n_{1}-1,\ n_{2},\ n_{3}-1)\\}$
$f_{1}+f_{2}<1,\;f_{2}+f_{3}>1.$

  * •

$\\{(n_{1},\ n_{2}-1,\ n_{3}),\ (n_{1}-1,\ n_{2},\ n_{3}),\ (n_{1}-1,\ n_{2},\
n_{3}-1),\\\ (n_{1}-1,\ n_{2}-1,\ n_{3})\\}$
$f_{1}+f_{2}>1,\;f_{2}+f_{3}<1.$

  * •

$\\{(n_{1},\ n_{2}-1,\ n_{3}),\ (n_{1},\ n_{2}-1,\ n_{3}-1),\ (n_{1}-1,\
n_{2},\ n_{3}-1),\\\ (n_{1}-1,\ n_{2}-1,\ n_{3})\\}$
$f_{1}+f_{2}>1,\;f_{2}+f_{3}>1,\;f_{1}+f_{2}+f_{3}<2.$

  * •

$\\{(n_{1},\ n_{2}-1,\ n_{3}-1),\ (n_{1}-1,\ n_{2},\ n_{3}-1),\ (n_{1}-1,\
n_{2}-1,\ n_{3}),\\\ (n_{1}-1,\ n_{2}-1,\ n_{3}-1)\\}$
$f_{1}+f_{2}+f_{3}>2.$

In the next section we will apply and extend the preceding results to generate
the collection of observable vectors of sample counts corresponding to
$R_{1},\ \dots,\ R_{m}$.

##  3 Main Results

To simplify the discussion we will henceforth assume that there are no
integers $i$ and $h$ for which $\sum_{j=0}^{h}f_{i+j}$ is an integer, but the
set of observable vectors of sample counts in such cases is an intersection of
sets of observable vectors in the cases we consider. The main idea in the
construction of observable vectors of sample counts is to use Propositions 1
and 3 as base cases and to progressively extend the $K+1$ projections of
$(\eta_{1},\ \dots,\ \eta_{m})$ to $(\eta_{i},\ \eta_{i+1},\ \dots,\
\eta_{i+K-1}),\;i\geq 1,\ i+K\leq m$ to the set of $K+2$ projections
$(\eta_{i},\ \eta_{i+1},\ \dots,\ \eta_{i+K}).$ We introduce the following
notation. Let $\lfloor x\rfloor,\;x\geq 0,$ denote the integer part of $x$.
For non-negative integers $K$, define

$\kappa(i,\ K)\;=\;\lfloor\sum_{j=0}^{K}f_{i+j}\rfloor.$

If $\kappa(i,\ K)\geq 1,$ then for $h\in\\{1,\ \dots,\kappa(i,\ K)\\}$ let
$\beta_{h}(i,\ K)$ be the non-negative integer for which

$\sum_{j=\beta_{h}(i,\ K)+1}^{K}f_{i+j}<h\;\mbox{and}\;\sum_{j=\beta_{h}(i,\
K)}^{K}f_{i+j}>h.$

To simplify notation, we will assume in the following that $i$ is fixed and
write $\kappa(K)=\kappa(i,\ K)$ and $\beta_{h}(K)=\beta_{h}(i,\ K)$. We first
identify vectors of sample counts which have two extensions.

Theorem 5: Suppose we have all $K+1$ vectors of sample counts
$(\eta_{i},\ \eta_{i+1},\ \dots,\ \eta_{i+K-1})$ and $i+K\leq m$. If
$\kappa(K)=0$, then the vector with both extensions $\eta_{i+K}=n_{i+K}$ and
$\eta_{i+K}=n_{i+K}-1$ satisfies $\eta_{i+j}=n_{i+j},\;j\in\\{0,\ \dots,\
K-1\\}$ . If $\kappa(K)\geq 1$, then the vector with both extensions
$\eta_{i+K}=n_{i+K}$ and $\eta_{i+K}=n_{i+K}-1$ satisfies

$\displaystyle\eta_{i+j}=n_{i+j}-1,$ $\displaystyle
j\in\\{\beta_{\kappa(K)}(K),\ \dots,\beta_{1}(K)\\}$
$\displaystyle\eta_{i+j}=n_{i+j},$ $\displaystyle
j\not\in\\{\beta_{\kappa(K)}(K),\ \dots,\beta_{1}(K)\\}$

Proof: We begin by showing

$\sum_{j=0}^{K-1}\eta_{i+j}=\left(\sum_{j=0}^{K-1}n_{i+j}\right)-{\kappa(K)}.$

Since
$\sum_{j=0}^{K}R_{i+j}=\left(\sum_{j=0}^{K}n_{i+j}\right)-\left(\sum_{j=0}^{K}f_{i+j}\right)$,
by an extension of Proposition 1, every observation $(\eta_{i}^{{}^{\prime}},\
\eta_{i+1}^{{}^{\prime}},\ \dots,\ \eta_{i+K}^{{}^{\prime}})$ of sample counts
from regions $i$ to $i+K$ of the support of $g(t)$ satisfies

$\left(\sum_{j=0}^{K}n_{i+j}\right)-({\kappa(K)}+1)\leq\sum_{j=0}^{K}\eta_{i+j}^{{}^{\prime}}\leq\left(\sum_{j=0}^{K}n_{i+j}\right)-{\kappa(K)}.$
(1)

Suppose we consider an observation with

$\sum_{j=0}^{K-1}\eta_{i+j}^{{}^{\prime}}\leq\left(\sum_{j=0}^{K-1}n_{i+j}\right)-({\kappa(K)}+1).$

Then the extension $(\eta_{i}^{{}^{\prime}},\ \eta_{i+1}^{{}^{\prime}},\
\dots,\ \eta_{i+K-1}^{{}^{\prime}},\ n_{i+K}-1)$ satisfies

$\sum_{j=0}^{K}\eta_{i+j}^{{}^{\prime}}\leq\left(\sum_{j=0}^{K}n_{i+j}\right)-({\kappa(K)}+2),$

which violates (1). Next consider an observation with

$\sum_{j=0}^{K-1}\eta_{i+j}^{{}^{\prime}}\geq\left(\sum_{j=0}^{K-1}n_{i+j}\right)-({\kappa(K)}-1).$

Then the extension $(\eta_{i}^{{}^{\prime}},\ \eta_{i+1}^{{}^{\prime}},\
\dots,\ \eta_{i+K-1}^{{}^{\prime}},\ n_{i+K})$ satisfies

$\sum_{j=0}^{K}\eta_{i+j}^{{}^{\prime}}\geq\left(\sum_{j=0}^{K}n_{i+j}\right)-({\kappa(K)}-1),$

which also violates (1). Hence, for the observation of interest,

$\sum_{j=0}^{K-1}\eta_{i+j}=\left(\sum_{j=0}^{K-1}n_{i+j}\right)-{\kappa(K)}.$

If ${\kappa(K)}=0$, then by Proposition 1, $\eta_{i+j}=n_{i+j},\;j\in\\{0,\
\dots,\ K-1\\}$. To complete the proof, assume ${\kappa(K)}\geq 1$. It remains
to show that

$\eta_{i+j}=n_{i+j}-1,j\in\\{\beta_{\kappa(K)}(K),\ \dots,\beta_{1}(K)\\}.$

We proceed with an argument by contradiction. Suppose instead that there are
$b_{\kappa(K)},\ \dots,b_{1}$ with $b_{\kappa(K)}<b_{\kappa(K)-1}<\dots<b_{1}$
for which

$\eta_{i+j}=n_{i+j}-1,j\in\\{b_{\kappa(K)},\ \dots,b_{1}\\},$

and $(b_{\kappa(K)},\ \dots,b_{1})\neq(\beta_{\kappa(K)}(K),\
\dots,\beta_{1}(K))$. Then there is a smallest index $k\in\\{1,\ \dots,\
{\kappa(K)}\\}$ for which $b_{k}\neq\beta_{k}(K)$. If $b_{k}>\beta_{k}(K)$,
then

$k-1<\sum_{j=b_{k}}^{K}f_{i+j}<k.$

By an extension of Proposition 1,

$\left(\sum_{j=b_{k}}^{K}n_{i+j}\right)-k\leq\sum_{j=b_{k}}^{K}\eta_{i+j}\leq\left(\sum_{j=b_{k}}^{K}n_{i+j}\right)-(k-1).$
(2)

If we choose the extension with $\eta_{i+K}=n_{i+K}-1$, then

$\sum_{j=b_{k}}^{K}\eta_{i+j}=\left(\sum_{j=b_{k}}^{K}n_{i+j}\right)-(k+1),$

contradicting (2). Therefore, $b_{k}<\beta_{k}(K)$. Hence, by an extension of
Proposition 1,

$\left(\sum_{j=\beta_{k}(K)}^{K}n_{i+j}\right)-(k+1)\leq\sum_{j=\beta_{k}(K)}^{K}\eta_{i+j}\leq\left(\sum_{j=\beta_{k}(K)}^{K}n_{i+j}\right)-k.$
(3)

If we choose the extension with $\eta_{i+K}=n_{i+K}$, then

$\sum_{j=\beta_{k}(K)}^{K}\eta_{i+j}=\left(\sum_{j=\beta_{k}(K)}^{K}n_{i+j}\right)-(k-1),$

contradicting (3). $\Box$

The final step in generating all possible sample counts given region lengths
$R_{i},\ \dots,\ R_{i+K}$ is determining how to uniquely extend an observation
$(\eta_{i},\ \eta_{i+1},\ \dots,\ \eta_{i+K-1})$ of sample counts which is not
covered by Theorem 5. We have the following result.

Theorem 6: Suppose we have all $K+1$ vectors of sample counts
$(\eta_{i},\ \eta_{i+1},\ \dots,\ \eta_{i+K-1})$.

  * •

If $\sum_{j=0}^{K-1}\eta_{i+j}=(\sum_{j=0}^{K-1}n_{i+j})-(\kappa(K)+1)$, then
$\eta_{i+K}=n_{i+K}$.

  * •

If $\sum_{j=0}^{K-1}\eta_{i+j}=(\sum_{j=0}^{K-1}n_{i+j})-(\kappa(K)-1)$, then
$\eta_{i+K}=n_{i+K}-1$.

  * •

If there are $b_{\kappa(K)},\ \dots,b_{1}$ with
$b_{\kappa(K)}<b_{\kappa(K)-1}<\dots<b_{1}$ for which

$\eta_{i+j}=n_{i+j}-1,\;j\in\\{b_{\kappa(K)},\ \dots,b_{1}\\},$

and $(b_{\kappa(K)},\ \dots,b_{1})\neq(\beta_{\kappa(K)}(K),\
\dots,\beta_{1}(K))$, then let
$k\in\\{1,\ \dots,\ {\kappa(K)}\\}$ be the smallest index for which
$b_{k}\neq\beta_{k}(K)$. If $b_{k}>\beta_{k}(K)$, then $\eta_{i+K}=n_{i+K}$.
If $b_{k}<\beta_{k}(K)$, then $\eta_{i+K}=n_{i+K}-1$.

Proof: Since $\kappa(K)<\sum_{j=0}^{K}f_{i+j}<\kappa(K)+1,$ it follows from an
extension of Proposition 1 that

$\left(\sum_{j=0}^{K}n_{i+j}\right)-({\kappa(K)}+1)\leq\sum_{j=0}^{K}\eta_{i+j}\leq\left(\sum_{j=0}^{K}n_{i+j}\right)-{\kappa(K)}.$

Furthermore, either $\kappa(K-1)=\kappa(K)$ or $\kappa(K-1)=\kappa(K)-1.$
These facts imply the first two cases. For the final case, if
$b_{k}>\beta_{k}(K)$, then $k-1<\sum_{j=b_{k}}^{K}f_{i+j}<k.$ By (2),
$\eta_{i+K}=n_{i+K}$. Finally, if $b_{k}<\beta_{k}(K)$, then (3) implies that
$\eta_{i+K}=n_{i+K}-1.\;\Box$

Observe that if we have the full set of $m+1$ possible observation vectors of
sample counts of the regions in the support of $g(t)$, then Propositions 1 and
3 and Theorems 5 and 6 define $n_{i},\ i\in\\{1,\ \dots,\ m\\},$ and in which
of the $m!$ partitions of the unit hypercube $(f_{1},\ \dots,\ f_{m})$
resides. However, we may not have the full set, and then the best we can do is
to learn from projections corresponding to single regions or to contiguous
sets of regions. Continuing the discussion following Proposition 3, if our
projections for regions $1$ and $2$ only yield $(n_{1},n_{2}-1)$ and
$(n_{1}-1,\ n_{2})$, then we learn the values of $n_{1}$ and $n_{2}$, but we
do not obtain the relationship between $f_{1}$ and $f_{2}$. If instead we only
see $(n_{1},n_{2})$ and $(n_{1},n_{2}-1)$, then we learn the value of $n_{2}$,
but if there is a third vector to complete the set, then it could be either
$(n_{1}-1,\ n_{2})$ or $(n_{1}+1,n_{2}-1)$, so we cannot infer $n_{1}$ or the
relationship between $f_{1}$ and $f_{2}$. However, sometimes two observations
can provide some subinterval information. For example, if $m=3$ and we observe
the pair of sample counts $(n_{1},n_{2}-1,n_{3})$ and
$(n_{1}-1,n_{2},n_{3}-1)$, then we know $n_{1},\ n_{2}$ and $n_{3}$;
furthermore, since in one case
$\eta_{1}+\eta_{2}+\eta_{3}=n_{1}+n_{2}+n_{3}-1$ and in the other case
$\eta_{1}+\eta_{2}+\eta_{3}=n_{1}+n_{2}+n_{3}-2$, then we also know that
$1<f_{1}+f_{2}+f_{3}<2$. To look at a simlar example for $m=5$, if we receive
the pair of sample counts $(n_{1},n_{2}-1,n_{3},n_{4}-1,n_{5})$ and
$(n_{1}-1,n_{2},n_{3}-1,n_{4},n_{5}-1)$, then we can infer $n_{1},\ n_{2},\
n_{3},\ n_{4},\ n_{5},1<f_{1}+f_{2}+f_{3}<2,\ 1<f_{2}+f_{3}+f_{4}<2,\\\
1<f_{3}+f_{4}+f_{5}<2,$ and $2<f_{1}+f_{2}+f_{3}+f_{4}+f_{5}<3$.

##  4 Conclusions and a Few Remarks on Noise

We have initiated a new theoretical approach to the study of discrete image
registration. For the basic and idealized problem that we considered, we
established limits on subinterval/subpixel estimation in the absence of prior
information about the function of interest or the displacement process as well
as the absence of distortions or considerations about the procurement of
sample count data. However, a brief commentary on noise may be helpful. First
consider the simplest case where all but one image has no noise and we are
given a single image with noise. Then we know the number of regions of the
support and the value of the function in each region. Therefore, the goal is
is determine the sample counts for the regions in the support of the noisy
image; there are finitely many possibilities with constraints on sample
counts, $\sum_{j=0}^{K}\eta_{i+j},i\in\\{1,\ \dots,\ m\\},\;K\in\\{0,\ \dots,\
m-i\\}$ that are consistent with the other data. Part of the solution involves
the segmentation problem of establishing where the support begins and ends.
The overall problem is a detection/hypothesis testing/classification problem,
and there is a large literature on the detection of signals in noise [10]. The
more realistic problem is global image registration when all of the images are
noisy[6]. For example, in semiconductor metrology the critical dimension
scanning electron microscope can interact with and alter the measurement
sample [11, p. 11], and so there is a desire to minimize these effects; this
may compromise image fidelity, which is also a source of concern [11, p. 120].
Returning to our one-dimensional problem where the support of $g(t)$ does not
change over multiple observations apart from translation, it would be ideal to
jointly process the images to obtain a consistent composite belief about the
number of regions $m$ in the support of $g(t)$, the value of the function in
each region, and the bounds on $\sum_{j=0}^{K}\eta_{i+j},i\in\\{1,\ \dots,\
m\\},\;K\in\\{0,\ \dots,\ m-i\\}$. Deep learning [12] approaches may be
natural candidates for this objective. Finally, we end this paper on alignment
by mentioning that the term “noise” has also been applied to problems
involving simpler human judgment[13]; there are questions in Ref. 13 about the
practice of working with a single ground truth.

## References

  * [1] P. E. Anuta, “Spatial registration of multispectral and multitemporal digital imagery using fast Fourier transform techniques,” IEEE Trans. Geosci. Electron. 8(4), 353-368, Oct. 1970.
  * [2] X. Tong, Z. Ye, Y. Xu, S. Gao, H. Xie, Q. Du, S. Liu, X. Xu, S. Liu, K. Luan, and U. Stilla, “Image registration with Fourier-based image correlation: a comprehensive review of developments and applications,” IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 12(10), 4062-4081, Oct. 2019.
  * [3] R. DeVore, W. Shao, J. Pierce, E. Kaymaz, B. Lerner, and W. Campbell, “Using nonlinear wavelet compression to enhance image registration,” Proc. SPIE, 3028, 539-551, April 1997.
  * [4] H. J. Landau and H. O. Pollak, “Prolate spheroidal wave functions, Fourier analysis, and uncertainty – III,” Bell Sys. Tech. J. 41, 1295-1336, 1962.
  * [5] C. L. Dotson, Fundamentals of Dimensional Metrology, Fifth Edition, Thomson Delmar Learning, New York, 2006.
  * [6] M. T. Postek, A. E. Vládar and P. Cizmar, “Nanomanufacturing concerns about measurements made in the SEM Part III: Vibration and drift,” Proc. SPIE, 9173, 917306, 2014.
  * [7] R. M. Oliver, J. R. Pierce, and C. E. Shannon, “The philosophy of PCM,” Proc. IRE, 36, pp. 1324-1331, Nov. 1948.
  * [8] N. S. Jayant, ed., Waveform Quantization and Coding, IEEE Press, New York, 1976.
  * [9] R. C. Gonzalez and R. E. Woods, Digital Image Processing, Fourth Edition, Pearson, New York, NY, 2018.
  * [10] R. N. McDonough and A. D. Whalen, Detection of Signals in Noise, Second Edition, Academic Press, San Diego, CA, 1995.
  * [11] B. Su, E. Solecky, and A. Vaid, Introduction to Metrology Applications in IC Manufacturing, SPIE Press, Bellingham, WA, 2015.
  * [12] Y. LeCun, Y. Bengio, and G. Hinton, “Deep learning,” Nature, 521(7553), pp. 436-444, 2015.
  * [13] M. Kejriwal, H. Santos, K. Shen, A. M. Mulvehill, and D. L. McGuiness, “A noise audit of human-labeled benchmarks for machine commensense reasoning,” Nature Portfolio Scientific Reports, 14:8609, 2024.


# On the Analysis of AoI-Reliability Tradeoff in Heterogeneous IIoT Networks

Hossam Farag, Syed Muhammad Ali and Čedomir Stefanović  Department of
Electronic Systems, Aalborg University, Denmark
Email: {hmf, syedma<EMAIL_ADDRESS>
###### Abstract

Age of information (AoI) and reliability are two critical metrics to support
real-time applications in Industrial Internet of Things (IIoT). These metrics
reflect different concepts of timely delivery of sensor information.
Monitoring traffic serves to maintain fresh status updates, expressed in a low
AoI, which is important for proper control and actuation actions. On the other
hand, safety-critical information, e.g., emergency alarms, is generated
sporadically and must be delivered with high reliability within a predefined
deadline. In this work, we investigate the AoI-reliability trade-off in a
real-time monitoring scenario that supports two traffic flows, namely AoI-
oriented traffic and deadline-oriented traffic. Both traffic flows are
transmitted to a central controller over an unreliable shared channel. We
derive expressions of the average AoI for the AoI-oriented traffic and
reliability, represented by Packet Loss Probability (PLP), for the deadline-
oriented traffic using Discrete-Time Markov Chain (DTMC). We also conduct
discrete-event simulations in MATLAB to validate the analytical results and
evaluate the interaction between the two types of traffic flows. The results
clearly demonstrate the tradeoff between the AoI and PLP in such heterogeneous
IIoT networks and give insights on how to configure the network to achieve a
target pair of AoI and PLP.

##  I Introduction

Industrial Internet of Things (IIoT) is a key pillar for the Industry 4.0
paradigm to enable smart manufacturing [1]. In IIoT applications, the field
network, realized by Industrial Wireless Sensor Network (IWSN), allows remote
monitoring of wide-ranging industrial processes, thereby enabling efficient
and sustainable production. IWSNs typically serve the communication of
different traffic flows that are characterized by different generation
patterns (event-triggered and time-triggered) and communication requirements
in terms of latency, throughput and timeliness [2]. In process monitoring and
control scenarios, emergency alarms and safety-critical information represent
the event-triggered flows that must be transmitted within stringent deadline
constraints to maintain system stability and avoid dangerous consequences. On
the other hand, regular, periodic monitoring traffic represents the time-
triggered traffic of sensor readings to be transmitted to the central
controller. The freshness of such traffic is crucial to keep the central
controller updated with the status of the industrial process to drive
decisions or feedback loops. The information freshness however, cannot be
captured by traditional performance metrics, such as delay and throughput. Age
of Information (AoI) [3] was introduced as a relevant metric for quantifying
information freshness from the perspective of the central controller. AoI is
defined as the time elapsed since the latest received packet was generated,
thereby it is intrinsically different from classical performance metrics like
throughput and delay, which focuses on a single packet and only captures its
network time. It has been proven that the adopted strategies to minimize delay
or to maximize throughput are not necessarily optimal for minimizing AoI [4].

Several research works were devoted to the study and improvement of the
performance of IWSNs with heterogeneous traffic [2, 5, 6, 7]. The main goal of
these approaches is to give critical traffic the highest transmission priority
while sacrificing the delivery and reliability of less-critical traffic.
However, none of these works consider the AoI performance, where the freshness
of less critical traffic might be subject to AoI constraints. In [8], the
authors propose improved multiple access schemes to improve the AoI
performance of energy harvesting IWSNs. The authors in [9] study the AoI
performance in UAV-aided IWSNs where they formulate two optimization problems
to minimize the average and maximal AoI. The work in [10] presents two greedy
scheduling policies to minimize the AoI and jitter in industrial cyber-
physical systems. All these works focus on analyzing or optimizing the
performance of AoI in IWSNs where the network supports only a single traffic
flow of AoI-oriented data, which is not the case in most IIoT applications. In
heterogeneous IIoT networks, different traffic flows are associated with
different performance metrics and their corresponding sources contend over
unreliable multi-access channel. The authors in [11] propose an optimal
generation policy for the status updates in a heterogeneous IoT network that
serves AoI-sensitive and AoI-insensitive traffic flows, however, this work
focuses only on the optimal AoI performance of the AoI-sensitive traffic and
disregards the performance of the AoI-insensitive one. In typical IIoT
applications, the AoI-insensitive flow could be a critical traffic that is
characterized by stringent deadline constraint, and its reliability (i.e., the
ratio of packets delivered within the deadline bound) is crucial for the
stability and functionality of the system.

In this work, we investigate the AoI-reliability tradeoff in heterogeneous
IIoT network that supports two traffic flows, namely AoI-oriented traffic and
deadline-oriented traffic. The packets from both traffic flows are transmitted
to a common central controller via unreliable multi-access channel. We derive
the average AoI for the AoI-oriented traffic and the reliability, represented
by Packet Loss Probability (PLP), for the deadline-oriented traffic using
Discrete-Time Markov Chain (DTMC). We also validate our analysis through
discrete-event simulations via MATLAB. The obtained results demonstrate the
interaction between the two traffic flows and show that the overall network
performance is mainly influenced by the access probabilities. The analysis and
the results in this work give insights on how to configure the system to
achieve a target pair of AoI and PLP, which could be through adopting enhanced
channel access and/or queue management strategies.

The remainder of the text is organized as follows. Section II describes the
network model and the basic system parameters. Section III presents the DTMC
analysis of the AoI, followed by the DTMC analysis of the queue size and PLP
in Section IV. The results are given in Section V, and finally the paper is
concluded in Section VI.

##  II Network Model

We consider a static deployment of an IWSN in a process monitoring scenario
where a set of sensor nodes are randomly distributed to monitor an industrial
process. All sensory information is transmitted to a central controller using
a time-slotted random access via a fading channel. Specifically, we consider a
slow and flat Rayleigh fading channel with additive white Gaussian noise
within a time slot. The network supports two types of traffic flows, namely
deadline-oriented traffic ($\mathrm{T_{D}}$) and AoI-oriented traffic
($\mathrm{T_{AoI}}$). The $\mathrm{T_{D}}$ traffic represents alarms and
safety-critical data that are generated at emergency events and must be
delivered to the central controller within a predefined deadline. The
$\mathrm{T_{AoI}}$ traffic represents time-triggered sensory data where the
goal is to keep the central controller’s received information as fresh as
possible. Without loss of generality, we consider one node that transmits
$\mathrm{T_{D}}$ traffic and $N$ nodes transmitting the $\mathrm{T_{AoI}}$
traffic. At each time slot, a packet arrives to the infinite buffer of the
$\mathrm{T_{D}}$ with probability $\lambda$. Further, each $\mathrm{T_{D}}$
packet is attached with a constant deadline of $D$, such that it should be
delivered within $D$ time slots since its arrival, otherwise it is dropped. At
the beginning of a time slot, the $\mathrm{T_{D}}$ node (if its queue is non-
empty) attempts to transmit the packet at the head of the queue with
probability $p_{1}$. When the transmission of the $\mathrm{T_{D}}$ fails, the
packet is retransmitted until it is either successfully received or its
deadline has expired. The $\mathrm{T_{AoI}}$ traffic follows the generate-at-
will model [12], in which a $\mathrm{T_{AoI}}$ node generates a fresh sample
when it decides to transmit, i.e., access the channel with probability
$p_{2}$. After the transmission attempt, the $\mathrm{T_{AoI}}$ node discards
the packet, i.e., there is no retransmissions. All such nodes are assumed to
be synchronized, and the packet arrivals align with the boundary of the time
slot. Acknowledgments of successful transmission are received within the same
time slot via an error-free channel.

We consider the capture effect [13], where the central controller can
successfully decode a packet if the received Signal-to-Interference-plus-Noise
Ratio (SINR) exceeds a certain threshold $\gamma$ (capture ratio). The
threshold $\gamma$ is determined according to a certain packet error
probability as a function of packet length, modulation, channel coding,
diversity and receiver design [14]. Let $S$ denotes the set of nodes
concurrently transmitting within the same time slot, the $\text{SINR}_{i}$ at
the central controller corresponding to an arbitrary transmitting node $i$ is

$\text{SINR}_{i}=\frac{P_{i}|h_{i,c}|^{2}d_{i,c}^{-\alpha}}{\sigma^{2}+\sum_{j\in
S\setminus\\{i\\}}P_{j}|h_{j,c}|^{2}d_{j,c}^{-\alpha}},$ (1)

where $P_{i}$ is the transmitting power of node $i$, $h_{i,c}$ is the Rayleigh
random variable of the channel between node $i$ and the central controller $c$
($|h_{i,c}|^{2}$ is exponentially distributed [15]), $d_{i,c}$ denotes the
distance between the node $i$ and the central controller, $\sigma^{2}$ is the
noise power and $\alpha$ is the path loss exponent. In order for the central
controller to successfully decode a received packet, it should be
$\text{SINR}_{i}>\gamma$ and the probability of this event is given as

$\begin{split}\mathds{P}(\text{SINR}_{i}>\gamma)&=\mathrm{exp}\left(-\frac{\gamma\sigma^{2}}{P_{i}|h_{i,c}|^{2}d_{i,c}^{-\alpha}}\right)\times\\\
&\prod_{j\in
S\setminus\\{i\\}}\left(1+\gamma\frac{P_{j}|h_{j,c}|^{2}d_{j,c}^{-\alpha}}{P_{i}|h_{i,c}|^{2}d_{i,c}^{-\alpha}}\right)^{-1}.\end{split}$
(2)

Based on (2), the successful update probability per time slot $q_{D}$ of the
$T_{D}$ node, i.e, the probability of successfully receiving a $T_{D}$ packet
within a time slot, is given as

$q_{D}=\sum_{k=0}^{N}\binom{N}{k}p_{1}(1-p_{2})^{N-k}{p_{2}}^{k}P_{S_{k}},$
(3)

where $P_{S_{k}}=\mathds{P}(\text{SINR}_{i}>\gamma\,\,|\,\,|S|=k)$, which can
be obtained via (2).

The successful update probability $q_{{AoI}}$ of a $\mathrm{T_{AoI}}$ node,
which reflects a successful delivery of an update at the end of a time slot,
depends on the queue status of the $T_{D}$ node. Let $Q$ be a random variable
that represents the status of the output queue of the $T_{D}$ node (empty or
non empty), then $q_{{AoI}}$ can be expressed as

$\begin{split}q_{{AoI}}&=p_{2}[\underbrace{q_{S_{0}}\mathds{P}(Q=0)}_{L_{1}}+\underbrace{q_{S_{0}}(1-p_{1})\mathds{P}(Q>0)}_{L_{2}}\\\
&+\underbrace{q_{S_{1}}p_{1}\mathds{P}(Q>0)}_{L_{3}}],\end{split}$ (4)

where $q_{S_{0}}$ and $q_{S_{1}}$ represent the successful probabilities when
the $\mathrm{T_{D}}$ node is idle (not transmitting) and active
(transmitting), respectively. The successful update probability $q_{{AoI}}$ in
(4) is calculated by three terms corresponding to three different cases. The
term $L_{1}$ denotes an empty queue of the $\mathrm{T_{D}}$ node, thus not
transmitting in the current time slots, and only the $\mathrm{T_{AoI}}$ nodes
attempting to transmit. The term $L_{2}$ refers to non-empty queue of the
$\mathrm{T_{D}}$ node, but it decides not to transmit in the current time
slot. Finally, the term $L_{3}$ represents the case when the $\mathrm{T_{D}}$
node decides to transmit the packet at the head of the queue and contend with
the attempting $\mathrm{T_{AoI}}$ nodes.

$q_{S_{0}}=\sum_{k=0}^{N-1}\binom{N-1}{k}(1-p_{2})^{N-k-1}{p_{2}}^{k}P_{S_{k}},$
(5)
$q_{S_{1}}=\sum_{k=0}^{N-1}\binom{N-1}{k}(1-p_{2})^{N-k-1}{p_{2}}^{k}P_{S_{k+1}}.$
(6)

##  III Analysis of the AoI

In this section, we analyze the average AoI at the central controller
corresponding to the $\mathrm{T_{AoI}}$ traffic. We assume that the status
updates from all $\mathrm{T_{AoI}}$ nodes are equally important. Therefore,
without loss of generality, we consider an arbitrary $\mathrm{T_{AoI}}$ node
and evaluate the AoI in discrete time. The AoI represents the number of time
slots elapsed since the last received packet was generated. If $\Delta(t)$
denotes the AoI at the end of time slot $t$, then we have

$\Delta(t+1)=\begin{cases}\Delta(t)+1&\text{unsuccessful transmission}\\\
1&\text{successful transmission}.\end{cases}$ (7)

Figure 1: DTMC model of the $\mathrm{T_{AoI}}$ traffic.

Fig. 1 shows the Discrete-Time Markov Chain (DTMC) model of the AoI
corresponds to an arbitrary $\mathrm{T_{AoI}}$ node where each state
represents the AoI at the central controller. As illustrated in Fig. 1, the
AoI-based DTMC transits from any state $n$ to 1 only upon a successful
reception of a $\mathrm{T_{AoI}}$ packet, otherwise it transits to state
$n+1$. Let $X_{t}$ represents the value of $\Delta(t)$ at time slot $t$, then
the transition probability from state $n$ to state $m$ is
$P_{nm}=\mathds{P}(X_{t+1}=m\,\,|X_{t}=n)$, and the transition matrix
$\mathbf{P_{AoI}}$ is written as

$\mathbf{P_{AoI}}=\begin{bmatrix}q_{{AoI}}&1-q_{{AoI}}&0&0&\dots\\\
q_{{AoI}}&0&1-q_{{AoI}}&0&\dots\\\
\vdots&\vdots&\vdots&\ddots&\ddots\end{bmatrix}.$ (8)

The row vector
$\boldsymbol{\pi^{AoI}}=[\pi_{1},\pi_{2},......,\pi_{n-1},\pi_{n},...]$
represents the steady-state probability vector of the DTMC in Fig. 1, where
$\pi_{n}=\lim_{t\to+\infty}\mathds{P}(X_{t})=n$ denotes the probability that
the AoI is equal to $n$ at the steady state. Using the set of equations
$\boldsymbol{\pi^{AoI}}\mathbf{P_{AoI}}=\boldsymbol{\pi^{AoI}}$ and
$\sum_{n}\pi_{n}=1$, $\pi_{n}$ is obtained as

$\pi_{n}=q_{AoI}(1-q_{{AoI}})^{(n-1)}$ (9)

Accordingly, the average AoI $\overline{\Delta}$ is calculated as

$\begin{split}\overline{\Delta}&=\sum_{n=1}^{\infty}\pi_{n}n=\sum_{n=1}^{\infty}nq_{{AoI}}(1-q_{{AoI}})^{(n-1)}\\\
&=\frac{q_{{AoI}}}{1-q_{{AoI}}}\sum_{n=1}^{\infty}n(1-q_{{AoI}})^{n}.\end{split}$
(10)

Since we have $q_{{AoI}}<1$, (10) can be rewritten as

$\overline{\Delta}=\frac{q_{{AoI}}}{1-q_{{AoI}}}\frac{1-q_{{AoI}}}{q_{{AoI}}^{2}}=\frac{1}{q_{{AoI}}}.$
(11)

However, the average AoI cannot account for extreme AoI events occurring with
very low probabilities at IWSNs. As mentioned in Section I, the received
updates of the monitoring traffic ($\mathrm{T_{AoI}}$) are used for control
and actuation actions, which implies certain requirements on tolerated values
of AoI. Here we analyze the AoI violation probability, which is the
probability that the AoI exceeds a certain constraint, which can be expressed
as

$\begin{split}&\mathds{P}(\Delta>c)=1-\mathds{P}(\Delta\leq
c)=1-\sum_{n=1}^{c}\pi_{n}\\\
&=\frac{q_{AoI}}{1-q_{AoI}}\sum_{n=1}^{c}(1-q_{AoI})^{n}=(1-q_{AoI})^{c},\end{split}$
(12)

where $c$ is the target AoI constraint, which is specified according to a
considered application scenario.

Figure 2: DTMC model of the $\mathrm{T_{D}}$ traffic.

$\displaystyle\mathbf{P_{(D=3)}}=\begin{bmatrix}(1-\lambda)&\lambda&0&0\\\
(1-\lambda)q_{D}&\lambda q_{D}&(1-q_{D})&0\\\
q_{D}(1-\lambda)^{2}&\lambda(1-\lambda)q_{D}&\lambda q_{{}_{D}}&(1-q_{D})\\\
(1-\lambda)^{3}&\lambda(1-\lambda)^{2}&\lambda(1-\lambda)&\lambda\end{bmatrix}.$
(13)

$\displaystyle\mathbf{P_{D}}=\begin{bmatrix}(1-\lambda)&\lambda\\\
(1-\lambda)q_{D}&\lambda q_{D}&(1-q_{D})\\\
q_{D}(1-\lambda)^{2}&\lambda(1-\lambda)q_{D}&\lambda q_{D}&(1-q_{D})\\\
\vdots&\vdots&\vdots&\ddots&\ddots\\\
q_{D}(1-\lambda)^{D-1}&q_{D}\lambda(1-\lambda)^{D-2}&q_{D}\lambda(1-\lambda)^{D-3}&\dots&\lambda
q_{D}&(1-q_{D})\\\
(1-\lambda)^{D}&\lambda(1-\lambda)^{D-1}&\lambda(1-\lambda)^{D-2}&\dots&\lambda(1-\lambda)&\lambda\end{bmatrix}.$
(14)

##  IV Queue Size and Packet Loss Probability of $\mathrm{T_{D}}$ Traffic

The successful update probability $q_{AoI}$ in (4) and (5) is based on the
queue status of the $\mathrm{T_{D}}$ user $(\mathds{P}(Q>0))$. In this
section, we use a DTMC model to evaluate the distribution of the queue size of
the $\mathrm{T_{D}}$ node. The DTMC model is shown in Fig. 2, where each state
represents the waiting time (in number of time slots) of the Head-of-Line
(HoL) $\mathrm{T_{D}}$ packet since its arrival. The HoL packet is kept
waiting in the queue until it is successfully received, i.e., its
acknowledgment is received, hence the waiting time accounts for the time slot
in which the packet is transmitted. Since the $\mathrm{T_{D}}$ packet is
dropped if its waiting time exceeds the predefined constraint $D$, we have a
finite DTMC with $D+1$ states, and the state transitions can be illustrated as
follows.

Assuming that $D=3$ in the example in Fig. 2, there are 4 states for the DTMC.
The state 0 represents empty queue as there is no packet waiting, and remains
in the same state as long as there is no arrivals. The system transits from
state 0 to state one when a packet arrives (the HoL packet has the chance to
be delivered within one time slot). The system remains in state 1 when the HoL
packet is successfully transmitted and a new packet arrives ($\lambda q_{D}$).
A transition from state 1 to state 0 occurs when the HoL packet is
successfully transmitted and no packet arrives ($(1-\lambda)q_{D}$). The
system transits from state 1 to state 2 when the HoL is not successfully
transmitted. The system remains in state 2 when the HoL packet is transmitted
successfully and a new packet arrived in the previous slot ($\lambda q_{D}$).
A transition from state 2 to state 1 occurs when the HoL is successfully
transmitted and one packet arrives in the current slot while no packet arrived
in the previous slot ($\lambda(1-\lambda)q_{D}$). The system transits from
state 2 to state 0 when the HoL packet is successfully transmitted and there
were no arrivals within the previous two slots ($q_{D}(1-\lambda)^{2}$). A
transition from state 2 to state 3 occurs when the HoL packet is not
successfully transmitted, and the packet is now dropped. The system remains in
state 3 when the HoL packet arrived 3 slots before the current slot. The
system transits from state 3 to state 0 when no packet arrived within the
previous three slots. A transition from state 3 to state 1 occurs when a
packet arrived in the current slot and no packets arrived in the previous two
slots. Finally, the system transits from state 3 to state 2 when a packet
arrived in the previous slot and no packets arrived in the previous two slots.
The transition matrix for the DTMC with $D=3$ can be given as (13). For an
arbitrary deadline constraint $D$, the general transition matrix is expressed
as (14).

The vector
$\boldsymbol{\pi^{D}}=[\grave{\pi_{0}},\grave{\pi_{1}},......,\grave{\pi_{D-1}},\grave{\pi_{D}}]$
represents the steady-state probability vector of the general DTMC in Fig. 2.
$\boldsymbol{\pi^{D}}$ can be derived from the set of linear equations
$\boldsymbol{\pi^{D}}\mathbf{P_{D}}~{}=~{}\boldsymbol{\pi^{D}}$ and
$\sum_{n=0}^{D}\grave{\pi_{n}}=1$. Then, we have
$\mathds{P}(Q>0)=1-\grave{\pi_{0}}$. Moreover, the reliability of the
$\mathrm{T_{D}}$ traffic can be represented by $PLP$ with
$PLP=\grave{\pi_{D}}(1-q_{D})$, which also reflects the deadline violation
probability.

##  V Results and Discussion

In this section, we evaluate the $\Delta$-$PLP$ trade-off of the considered
heterogeneous network based on the presented analysis in Section III and
Section IV. We also validate our analysis by comparing the numerical results
with simulation results obtained via discrete-event simulations in MATLAB. In
the following results, we consider one $\mathrm{T_{D}}$ node located at
distance 30 m from the central controller and 5 $\mathrm{T_{AoI}}$ nodes
distributed in an isotropic directions around the central controller with
equal distance of 40 m. We set the transmission power for all nodes to 10 dbm
and the receiver noise power to -80 dbm.

Figure 3: The effect of $p_{1}$ on $\Delta$-$PLP$ tradeoff under varying
$\gamma$ with $0.1\leq p_{1}\leq 1$, $p_{2}=0.6$ and $\lambda=0.7$.

Fig. 3 and Fig. 4 show the evaluation of $\Delta$ and $PLP$ for different
values of capture threshold $\gamma$ with varying $p_{1}$ and $p_{2}$,
respectively. From these figures, we can see that the analytical results match
well with the simulation results, which validates our analysis in Section III
and Section IV. In Fig. 3 we plot the average AoI ($\Delta$) and $PLP$ by
varying $0.1\leq p_{1}\leq 1$ with 0.1 increment. From this figure, we can
observe that with high capture capability of the central controller (i.e.,
$\gamma=-2$ dbm and $\gamma=-1$ dbm), the transmission probability of the
$\mathrm{T_{D}}$ node has insignificant effect on the AoI of the
$\mathrm{T_{AoI}}$ nodes, while for low capture capability (i.e., $\gamma=0$
dbm and $\gamma=2$ dbm), the AoI of the $\mathrm{T_{AoI}}$ nodes is highly
affected by higher values of $p_{1}$. For instance, with $\gamma=-2\,$dbm, the
AoI increases by 38% when $p_{1}$ changes from 0.1 to 0.8, while it increases
by 120% at $\gamma=2\,$dbm.

Figure 4: The effect of $p_{2}$ on the AoI and PLP trade-off under varying
$\gamma$ with $p_{1}=0.6$, $0.1\leq p_{2}\leq 1$ and $\lambda=0.5$.

In Fig. 4, we show the $\Delta$-$PLP$ tradeoff by fixing $p_{1}$ to 0.6 and
varying $p_{2}$ under different values of $\gamma$. We can observe that for
low values of $\gamma$, we can decrease $\Delta$ while maintaining low $PLP$
(below 0.2) for the $\mathrm{T_{D}}$ node. On the other side, at high values
of $\gamma$, the increase in $p_{2}$ would decrease $\Delta$ but at the cost
of higher $PLP$ of the $\mathrm{T_{D}}$ node. The figure can be used as a
reference to select the best transmission strategy to achieve a target pair of
$\Delta$ and $PLP$. For instance, with $\gamma=-1\,$dbm, a target AoI below 10
can be achieved at $p_{2}=0.3$ and $PLP$ = 0.133. Therefore, the transmission
strategy in this case would be to allow the $\mathrm{T_{AoI}}$ nodes and the
$\mathrm{T_{D}}$ node transmit together. With $\gamma=2\,$dbm, to keep
$\Delta$ below 10, the transmission probability $p_{2}$ increases to 0.7 while
$PLP$ increases to 0.36, meaning that more $\mathrm{T_{D}}$ packets missing
the deadline and dropped. In this case, it would be better to adopt a
scheduled access strategy (e.g., round-robin) to achieve the target AoI while
maintaining low PLP.

Figure 5: The AoI violation probability under varying $p_{2}$ with $p_{1}=0.5$
and $\lambda=0.5$.

Fig. 5 shows the AoI violation probability over time obtained from (5) with
the target AoI constraints $c=3\,\textrm{and}\,5$ and
$\gamma=-2\,\textrm{dbm}\,\textrm{and}\,0\,\textrm{dbm}$. As we can see, the
AoI violation probability decreases as $p_{2}$ increases as the
$\mathrm{T_{AoI}}$ node attempts to transmit its packet more often. Moreover,
the AoI violation probability increases significantly for high values of
$\gamma$ as it is likely that a $\mathrm{T_{AoI}}$ packet is lost due to the
low capture capability. The figure gives a more detailed view of the
performance when the target is to keep the AoI below a certain constraint.

##  VI Conclusion

In this paper, we investigated the AoI-reliability tradeoff in heterogeneous
IIoT networks. Specifically, we considered two traffic flows transmitted via
unreliable multi-access channel to a central controller, where one flow
represents deadline-oriented traffic and the other flow represents AoI-
oriented traffic. We derived the average AoI for the AoI-oriented traffic and
the PLP for the deadline-oriented traffic using DTMC and validated our
analysis through simulations. The obtained results showed that the AoI-PLP
tradeoff is mainly influenced by the access probabilities and give insights on
how to configure the heterogeneous network to achieve a target AoI-PLP
performance. As future work, the analysis in this work can be extended and
utilized to propose channel access and /or queue management strategies that
improve the AoI performance while maintaining low PLP for more general network
setups.

## Acknowledgement

This paper has received funding from the European Union’s Horizon 2020
research and innovation programme under grant agreement No. 883315.

## References

  * [1] E. Sisinni, A. Saifullah, S. Han, U. Jennehag, and M. Gidlund, “Industrial internet of things: Challenges, opportunities, and directions,” _IEEE Transactions on Industrial Informatics_ , vol. 14, no. 11, pp. 4724–4734, 2018\.
  * [2] H. Farag, E. Sisinni, M. Gidlund, and P. Österberg, “Priority-aware wireless fieldbus protocol for mixed-criticality industrial wireless sensor networks,” _IEEE Sensors Journal_ , vol. 19, no. 7, pp. 2767–2780, 2019\.
  * [3] Y. Sun, E. Uysal-Biyikoglu, R. D. Yates, C. E. Koksal, and N. B. Shroff, “Update or wait: How to keep your data fresh,” _IEEE Transactions on Information Theory_ , vol. 63, no. 11, pp. 7492–7508, 2017.
  * [4] A. Maatouk, S. Kriouile, M. Assaad, and A. Ephremides, “The age of incorrect information: A new performance metric for status updates,” _IEEE/ACM Transactions on Networking_ , vol. 28, no. 5, pp. 2215–2228, 2020.
  * [5] H. Farag, S. Grimaldi, M. Gidlund, and P. Österberg, “Rea-6tisch: Reliable emergency-aware communication scheme for 6tisch networks,” _IEEE Internet of Things Journal_ , vol. 8, no. 3, pp. 1871–1882, 2021.
  * [6] H. Farag, M. Gidlund, and P. Österberg, “A delay-bounded mac protocol for mission- and time-critical applications in industrial wireless sensor networks,” _IEEE Sensors Journal_ , vol. 18, no. 6, pp. 2607–2616, 2018\.
  * [7] T. Qiu, K. Zheng, M. Han, C. L. P. Chen, and M. Xu, “A data-emergency-aware scheduling scheme for internet of things in smart cities,” _IEEE Transactions on Industrial Informatics_ , vol. 14, no. 5, pp. 2042–2051, 2018\.
  * [8] N. Hirosawa, H. Iimori, K. Ishibashi, and G. T. F. D. Abreu, “Minimizing age of information in energy harvesting wireless sensor networks,” _IEEE Access_ , vol. 8, pp. 219 934–219 945, 2020.
  * [9] J. Liu, P. Tong, X. Wang, B. Bai, and H. Dai, “Uav-aided data collection for information freshness in wireless sensor networks,” _IEEE Transactions on Wireless Communications_ , vol. 20, no. 4, pp. 2368–2382, 2021.
  * [10] D. Sinha and R. Roy, “Scheduling status update for optimizing age of information in the context of industrial cyber-physical system,” _IEEE Access_ , vol. 7, pp. 95 677–95 695, 2019.
  * [11] G. Stamatakis, N. Pappas, and A. Traganitis, “Optimal policies for status update generation in an iot device with heterogeneous traffic,” _IEEE Internet of Things Journal_ , vol. 7, no. 6, pp. 5315–5328, 2020.
  * [12] Y. Sun, E. Uysal-Biyikoglu, R. D. Yates, C. E. Koksal, and N. B. Shroff, “Update or wait: How to keep your data fresh,” _IEEE Transactions on Information Theory_ , vol. 63, no. 11, pp. 7492–7508, 2017.
  * [13] M. Zorzi and R. Rao, “Capture and retransmission control in mobile radio,” _IEEE Journal on Selected Areas in Communications_ , vol. 12, no. 8, pp. 1289–1298, 1994.
  * [14] G. Proakis and M. Salehi, _Digital Communications_. New York, USA: McGraw-Hill, 2000.
  * [15] D. Tse and P. Viswanath, _Fundamentals of wireless communication_. UK: Cambridge university press, 2005.