# Polarization entanglement with highly non-degenerate photon pairs
enhanced by effective walk-off compensation method

Sungeun Oh<EMAIL_ADDRESS>Thomas Jennewein
<EMAIL_ADDRESS>Department of Physics and Astronomy, and
Institute for Quantum Computing, University of Waterloo, Waterloo, Ontario,
Canada N2L3G1

###### Abstract

We demonstrate a novel entangled photon source (EPS) that produces
polarization entanglement in highly non-degenerate photon pairs, generated
through Type-0 spontaneous parametric down conversion (SPDC) using bulk
periodically poled Lithium Niobate (PPLN). This source is specifically
designed for entanglement-based quantum key distribution (EBQKD) between
ground and satellite, as part of the Quantum Encryption and Science Satellite
(QEYSSat) mission funded by the Canadian Space Agency (CSA). Through the
utilization of birefringent beam displacers, a hybrid Mach-Zehnder a Sagnac
interferometer, we ensure high polarisation contrast and stable interference
of the highly non-degenerate photon pairs, which however causes substantial
spatial and temporal walk-offs of the photon paths which poses a formidable
challenge. We introduce an effective compensation method using birefringent
crystal wedges to eliminate spatial and temporal walk-offs simultaneously.
From only $1.0\text{\,}\mathrm{mW}$ of pump power we observed a coincidence
rate of $N=(33.33\pm 0.05)$kHz, a significant improvement compared to the
absence of the spatial compensation. We also observed an estimated pair
generation rate of $(2.92\pm 0.12)$MHz and an entanglement visibility of
$(96.6\pm 0.3)\%$ , making it a promising source for long-distance quantum
communication for ground-to-satellite and fiber optic links.

††preprint: APS/123-QED

##  I Introduction

While numerous source for polarization entanglement with degenerate or near-
degenerate photon pairs have been demonstrated[1, 2, 3, 4], polarization
entanglement with highly non-degenerate photon pairs has received less
attentions. However, sources that produce highly non-degenerate photon pairs
is relevant[5, 6, 7, 8] in the context of quantum communication scenarios that
require distinct wavelengths. Of our particular interest is the establishment
of ground to space entanglement-based quantum key distribution (EBQKD)[9, 10],
which necessitates a high pair production rate of around 100 MHz[11], as well
as the utilization of connectivity to the satellite free-space link as well as
to the ground-based optical fiber link. The photon to be transmitted through
the ground-based optical fiber will preferably have a wavelength around
$1550\text{\,}\mathrm{nm}$, to make it compatible with existing optical fiber
infrastructure and make use of the low-loss transmission window. Photons with
such a wavelength are most effectively received by superconducting nanowire
single photon detectors (SNSPDs). However, because deploying SNSPDs on
satellites is challenging due cryogenic cooling requirements [12] , typically
Silicon based single-photon avalanche diodes (SPADs) are chosen for a space
based quantum receiver, such as in the QEYSsat payload. The Silicon-based
photon detectors operate at a high performance and have been proven to
operating in space[13, 14, 15, 16] . As a result, the ground to space link
will operate with a wavelength in the visible (VIS) or near-infrared (NIR)
range for efficient photodetection. We present an EPS specifically designed
for implementing EBQKD between ground-based and satellite-based systems, as
part of the Quantum Encryption and Science Satellite (QEYSSat) mission funded
by the Canadian Space Agency (CSA)[17, 18, 19]. We implement a novel EPS
system based on birefringent beam displacers and pairs of birefringent
compensation wedges that generates photon pairs at the two distinct photon
wavelengths. In particular, our novel compensation method effectively
addresses the spatial and temporal walk-off of the photons caused by the
birefringent beam displacer material and provides the degrees of freedom
required to attain a high photon pair coupling as well as high entanglement
quality in our EPS.

##  II Scheme to Generate Polarization Entangled non-degenerate Photon Pairs

Our EPS generates photon pairs through spontaneous parametric down conversion
(SPDC)[20] in optically nonlinear crystals. A continuous wave (CW) laser
(Frankfurt laser, Germany) provides a polarized pump beam at a wavelength of
$\lambda_{p}=$523.6\text{\,}\mathrm{nm}$$, with a bandwidth of
$\Delta\lambda_{p}=$0.1\text{\,}\mathrm{nm}$$. A $10\text{\,}\mathrm{mm}$ long
Type-0 PPLN crystal (Covesion, UK) with a poling period of
$\Lambda=$7.1\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ (at room temperature),
converts pump(p) photons into signal(s)-idler(i) photon pairs with the same
polarization under the following quasi-phase matching (QPM) conditions:

$\displaystyle\frac{2\pi}{\Lambda(T)}=k_{p}-k_{s}-k_{i}$ (1)

$k=2\pi n(\lambda,T)/\lambda$ are the wave vectors, where $n(\lambda,T)$ are
the wavelength- and temperature-dependent refractive indices of the crystals
for the corresponding fields. The reason for the crystal choice is that PPLNs
can achieve a higher photon pair generate rate compared to other conventional
nonlinear crystals, due to its high nonlinear coefficient[21]. With the
presented values, we obtain signal photons at a wavelength of
$\lambda_{s}=790.8nm$ suitable for free-space transmission, and idler photons
at a wavelength of $\lambda_{i}=1550nm$ suitable for ground transmission
through optical fiber. The associated spectral bandwidths $\Delta\omega$ (at
1/e) are estimated by the first-order Taylor expansion of equation 1.

$\displaystyle\Delta k\approx\left.\dfrac{\partial
k_{p}}{\partial\omega}\right|_{\omega_{p}}\Delta\omega_{p}-\left.\dfrac{\partial
k_{s}}{\partial\omega}\right|_{\omega_{s}}\Delta\omega_{s}-\left.\dfrac{\partial
k_{i}}{\partial\omega}\right|_{\omega_{i}}\Delta\omega_{s}=\dfrac{2\pi}{L}$
(2)

Using equation 2, the temporal widths $\tau=2\pi/\Delta\omega$ (at 1/e) are
$\tau_{s}=$2.68\text{\,}\mathrm{ps}$$ for the signal photon, and
$\tau_{i}=$2.61\text{\,}\mathrm{ps}$$ for the idler photon. We consider
entanglement of signal and idler photons with horizontal ($H$) and vertical
($V$) polarizations using an optical interferometer. The resulting Bell state
can be written as

$\displaystyle\ket{\text{Bell}}=\dfrac{1}{\sqrt{2}}\left(\ket{V_{s}V_{i}}+e^{i\left(\phi_{r}+\Delta\phi_{r}\right)}\ket{H_{s}H_{i}}\right)$
(3)

The overall relative phase between the $\ket{VV}$ and the $\ket{HH}$ states
can be represented as a combination of the relative phase resulting from the
path difference (denoted as $\phi_{r}$) and the relative phase variation
induced by temperature fluctuations in optical components, contributing to the
additional difference in paths (denoted as $\Delta\phi_{r}$). For instance,
maximal entanglement in the linear diagonal polarization is achieved when the
overall phase $\Delta\phi_{r}=0$. The phase variation from each beam path in
the interferometer due to thermal and other fluctuations becomes
$\Delta\phi=\delta\phi_{p}+\delta\phi_{s}+\delta\phi_{i}$, where each term is

$\displaystyle\delta\phi=\sum_{j}\dfrac{2\pi
L_{j}}{\lambda}\left(\dfrac{\partial n_{j}}{\partial
T}+n_{j}\alpha\right)\Delta T,$ (4)

representing the phase shift introduced by all transmissive optical
components, each indexed with $j$, within the interferometer, each with a
length of $L_{j}$. The thermo-optic coefficients $\partial n_{j}/\partial T$
and the thermal expansion coefficient $\alpha$ inside the bracket in equation
4 correspond to the changes in the refractive index and the material length,
respectively. One interferometer design widely used for quantum entanglement
generation is the Sagnac interferometer, which incorporates self-compensation
for relative phase stabilization [8]. However, a Sagnac interferometer with a
polarization beamsplitter, which is commonly used, has technical difficulties
in achieving entanglement in non-degenerate photon pairs with high-efficiency
transmission for both signal and idler photons, while simultaneously
accommodating their two orthogonal polarizations. The reason lies in the
wavelength dependency of the transmittivity and reflectivity ratio of the
polarization beamsplitter [22].

Figure 1: Schematic diagram of the interferometer implementing two beam
displacers and a half waveplate and including Sagnac loop. Two PPLN crystals
are placed in the loop, each responsible for the SPDC of one of the polarized
pump beams.

To overcome this challenge, we chose to build an optical interferometer
employing two birefringent beam displacers capable of splitting and
recombining two orthogonally polarized beams, and inserting a half waveplate
(HWP) between the two beam displacers in order to have two identical path
lengths [23, 24, 25]. Thereby the beam that was originally an ordinary ray
(non-deflected, denoted with subscript ’$o$’) and the extraordinary ray
(deflected ray, denoted with subscript ’$e$’) switch roles as their
polarizations are interchanged. Hence, this scheme allows for the effective
transmissions of all three waves through the interferometer. The perfect
alignment of the two beam displacers can be seen as an equivalent
configuration to the Sagnac interferometer, but this is always an ideal case.
To ensure self-phase compensation, we add an actual Sagnac loop into the
interferometer, by inserting a dichroic mirror along the beam path between the
half waveplate and the second beam displacer. This way, the dichroic mirror
reflects the two orthogonally polarized pump beams, allowing them to enter the
Sagnac loop, and transmits the signal and idler beams with both polarizations
effectively. Another advantage of including the Sagnac loop is the prevention
of back reflection from the interferometer to the pump laser. If we solely
apply the beam-displacing interferometer scheme, the returning pump beam,
which has not undergone down-conversion, will be recombined at the first beam
displacer through the same beam paths. However, with the presence of the
additional Sagnac loop, the beam disperses upon returning to the first beam
displacer, eliminating the risk of laser instability or even damage.

The phase variation of each of the separated beam paths along the Sagnac loop
(denoted with subscripts, CW (clock-wise) and CCW (counter clock-wise)) is
expressed as follows.

$\displaystyle\Delta\phi_{CW}=\delta\phi_{e}(\lambda_{p},T)+\delta\phi_{o}(\lambda_{s},T)+\delta\phi_{o}(\lambda_{i},T)$
(5)
$\displaystyle\Delta\phi_{CCW}=\delta\phi_{o}(\lambda_{p},T)+\delta\phi_{e}(\lambda_{s},T)+\delta\phi_{e}(\lambda_{i},T)$

The relative phase variation is then

$\displaystyle\Delta\phi_{r}$ $\displaystyle=\Delta\phi_{CW}-\Delta\phi_{CCW}$
(6)
$\displaystyle=\delta\phi_{r}(\lambda_{p},T)-\delta\phi_{r}(\lambda_{s},T)-\delta\phi_{r}(\lambda_{i},T)$

where $\delta\phi_{r}=\delta\phi_{e}-\delta\phi_{o}$. We make the assumption
that the phases introduced by the half waveplate and the dichroic mirror are
negligible in comparison to those induced by the beam displacers, due to their
much smaller thicknesses. Phase variations are then calculated to be
$\delta\phi_{r}(\lambda_{p},T)=6.22\pi$,
$\delta\phi_{r}(\lambda_{s},T)=4.09\pi$ and
$\delta\phi_{r}(\lambda_{i},T)=2.08\pi$. Hence, the relative phase is
$\Delta\phi_{r}=0.052\pi$, ensuring sufficient phase stability. No self-
compensation effect would introduce a phase shift of
$\Delta\phi_{r}=12.39\pi$. As a result, while self-compensation may not be
entirely perfect, it certainly helps quantum states remain usable for a longer
period of time.

##  III Spatiotemporal Walk-off Compensation

Achieving a high entangled photon pair generation rate requires high spatial
and temporal (phase) overlap of both the signal beams and the idler beams.
Even if the relative phase ($\phi_{r}$ in equation 3) is minimized to some
extent in the Sagnac loop of our source, there is still phase introduced by
the two beam displacers. The amount of displacement of the extraordinary ray
in the beam displacers can be determined from its angle of refraction, denoted
as $\theta_{e}$. This angle can be obtained from the refractive indices of the
ordinary and extraordinary rays (denoted as $n_{o}$ and $n_{e}$,
respectively), and the optic angle (angle between the optic axis and the axis
normal to the surface of the medium for the incident rays, denoted as
$\theta$) of the material.

$\displaystyle\tan{\theta_{e}}=\left(1-\dfrac{n_{o}^{2}}{n_{e}^{2}}\right)\dfrac{\tan{\theta}}{1+\dfrac{n_{o}^{2}}{n_{e}^{2}}\tan^{2}{\theta}}$
(7)

It is worth noting that equation 7 assumes the incident beam entering the
medium perpendicular to the surface of the medium. The separation between the
two rays created by the beam displacer with a length of $L$ is then
$D=\tan{\theta_{e}}$.

Figure 2: Diagrams illustrating how the signal and idler beams undergo spatial
$\Delta D$ and temporal $\Delta T$ walk-offs in the second beam displacer. If
the beams had maintained the pump wavelength, then no walk-off would have
occurred at the second beam displacer. The degree of displacement decreases as
the wavelength increases.

When the down-converted signal and idler photons pass through the second beam
displacer, they do not return to their original position, given that their
wavelengths are no longer the pump wavelength. This deviation from the
original position is referred to as spatial walk-off, denoted as $\Delta D$.

$\displaystyle\Delta
D_{s,i}=L\tan{\theta_{e,2}^{p}}-L\tan{\theta_{e,1}^{s,i}}$ (8)

$\theta_{e,1}$ and $\theta_{e,2}$ are the refraction angles of the first and
the second beam displacers, respectively, as shown in Figure 2. We use two
identical $\alpha$BBO beam displacers (Newlight Photonics inc, Canada) with
the lengths of $L=$39.4\text{\,}\mathrm{mm}$$ for the interferometer.
According to equation 8, the spatial walk-offs of the signal and the idler
beams are $0.10\text{\,}\mathrm{mm}$ and $0.17\text{\,}\mathrm{mm}$,
respectively. The time taken for the beam to pass through the beam displacer
can be easily conceptualized through a geometric approach.

$\displaystyle T_{o\rightarrow e}^{s,i}=$
$\displaystyle\dfrac{\mathcal{N}_{o}^{p}L}{c}+\dfrac{\mathcal{N}_{e}^{s,i}L}{c\cos{\theta_{e,2}^{s,i}}}$
(9) $\displaystyle T_{e\rightarrow o}^{s,i}=$
$\displaystyle\dfrac{\mathcal{N}_{e}^{p}L}{c\cos{\theta_{e,1}^{p}}}+\dfrac{\mathcal{N}_{o}^{s,i}L}{c}$

$T_{o\rightarrow e}$ is the time taken by the ray, which initially travels
straight and then becomes displaced. Similarly, $T_{e\rightarrow o}$ is the
time taken by the ray, which initially becomes displaced then travels
straight. $\mathcal{N}_{o}$ and $\mathcal{N}_{e}$ are group indices of the
o-ray and e-ray, respectively. The temporal walk-off is then the time
difference between the two travel time.

$\displaystyle\Delta T_{s,i}=T_{o\rightarrow e}^{s,i}-T_{e\rightarrow
o}^{s,i}$ (10)

Using equation 10, we obtain the temporal walk-offs of
$0.65\text{\,}\mathrm{ps}$ for the signal beam and $1.06\text{\,}\mathrm{ps}$
for the idler beam. Both spatial and temporal walk-offs not only reduce the
entangled photon pair rate but also degrade the entanglement visibility, and
this worsens with more extreme degenerate photon pairs. Employing birefringent
materials to compensate for the temporal walk-offs is a commonly used method.
Here, we additionally employ a separation between two birefringent crystal
wedges, allowing simultaneous compensations for both spatial and temporal
walk-offs.

Figure 3: Schematic diagram of two Calcite wedges used to compensate for both
spatial and temporal walk-offs. The optic axis of our Calcite wedges is
aligned perpendicular to the direction of ray propagation such that both the
ordinary (purple) and extraordinary (red) rays undergo no beam displacement
when entering the Calcite perpendicularly. However, upon exiting the first
wedge, the e-ray experiences more refraction than the o-ray due to negative
birefringence ($n_{e}<n_{o}$). The path lengths of the e-ray and the o-ray in
between the two wedges are shown as $P_{2,e}$ and $P_{2,o}$, respectively.

The choice of a wedge shape offers several advantages. As illustrated in
Figure 3 one side is flat, allowing the ray to enter perpendicularly. This
simplifies calculations for group indices appearing in equation 9. Moreover,
even if the ray exits the first wedge at a non-perpendicular angle, when it
enters the second wedge that is identical in shape, it regains its original
propagation direction. For our source, we use Calcite as birefringent wedges
to slow down the o-rays. A wedges angle
$\phi=$15\text{\,}\mathrm{\SIUnitSymbolDegree}$$ is chosen arbitrarily. Since
Calcite is a negative uniaxial material, o-rays experience more refraction
than the e-rays, in both the signal and the idler path. It is worth noting
that in the first Calcite wedge (on the left), the e-ray travels a slightly
longer distance than the o-ray by $\Delta D\tan{\phi}$. Also, when a gap
between the two crystal wedges is introduced, the distance the two beams
travel within the gap is slightly different, due to the different bendings of
the paths. Consequently, the o-ray travels a longer distance than the e-ray,
resulting in a reduction of the temporal walk-off. This leads to either an
increase or a decrease in the original temporal walk-off, and all of these
should be taken into account when determining the thickness of the crystal. By
identifying the correct lateral separation (shown as $d$ in Figure 3) between
the two wedges, we effectively merge the two rays at the second wedge,
compensating for the spatial walk-offs. At the same time, temporal walk-offs
can easily be compensated by the thickness of the wedges.

We determine the lateral separation $d$ needed between the wedges for a
perfect spatial overlap. The vertical displacement of the o-ray is
$b=a\tan{\left(\theta_{o}-\phi\right)}$ such that the vertical shift of the
e-ray is $b-\Delta D$. In doing so, a relationship between the two vertical
displacements can be established as follows.

$\displaystyle a\tan{\left(\theta_{o}-\phi\right)}-\Delta D=\left(a-\Delta
D\tan{\phi}\right)\tan{\left(\theta_{e}-\phi\right)}$ (11)

$\theta_{o,e}$, $\phi$ and $\Delta D$ are known quantities, so we can
rearrange equation 11 to determine $a$.

$\displaystyle
a=\dfrac{1-\tan{\phi}\tan{\left(\theta_{e}-\phi\right)}}{\tan{\left(\theta_{o}-\phi\right)}-\tan{\left(\theta_{e}-\phi\right)}}\Delta
D$ (12)

Finally, we calculate the lateral separation $d$ as follows.

$\displaystyle d=a-b\tan{\phi}$ (13)

The total path lengths for each ray can be expressed as follows.

$\displaystyle P(\text{o-ray})$ $\displaystyle=P_{1}+P_{2,o}+P_{3}$ (14)
$\displaystyle P(\text{e-ray})$ $\displaystyle=P_{1}+\Delta
D\tan{\phi}+P_{2,e}+P_{3}$

The difference in travel times between the two paths is then

$\displaystyle\Delta T$
$\displaystyle=\left[\dfrac{n_{o}}{c}\left(P_{1}+P_{3}\right)+\dfrac{1}{c}P_{2,o}\right]$
(15) $\displaystyle-\left[\dfrac{n_{e}}{c}\left(P_{1}+\Delta
D\tan{\phi}+P_{3}\right)+\dfrac{1}{c}P_{2,e}\right]$

By setting the overall thickness of the wedges to

$\displaystyle L$ $\displaystyle=\dfrac{c\Delta T}{\left|n_{e}-n_{o}\right|}$
(16)
$\displaystyle=\dfrac{\left(\mathcal{N}_{o}-\mathcal{N}_{e}\right)(P_{1}+P_{3})+P_{2,o}-P_{2,e}-\mathcal{N}_{e}\Delta
D\tan{\phi}}{\left|n_{e}-n_{o}\right|}$

we can effectively cancel out the temporal walk-off. To mention again, all
derivations take into account compensations involving negative birefringent
wedges (one can derive for positive birefringent wedges in a similar manner).

The less overlap between the two beams, the higher the relative phase in
equation 3. If we already know the ideal photon pair rate (representing
perfect beam overlap), we can intuitively analyze the impact of spatial and
temporal walk-offs by multiplying the ideal photon rate with a beam overlap
factor. We consider Gaussian fields and characterize the spatial modes of the
ordinary and extraordinary rays.

$\displaystyle
g(x,y)=\dfrac{1}{\sqrt{2\pi\sigma^{2}}}e^{-(x^{2}+y^{2})/2\sigma^{2}}$ (17)

The standard deviation can be obtained from the full-width-half-maximum,
FWHM$=2\sqrt{2\ln{2}}\sigma$. The sptial overlap of the two beam separated by
$\Delta D$ is then the product of the two Gaussians integrated over the
transverse plane (x,y) as follows.

$\displaystyle\int\limits_{-\infty}^{\infty}g_{o}(x,y)g_{e}(x-\Delta
D_{x},y-\Delta D_{y})dxdy$ (18)

The temporal overlap can be determined in a similar manner, using equation 17
for the temporal modes of the ordinary and extraordinary rays. The standard
deviation can be derived from the temporal width $\tau=2\sqrt{2\ln{2}}\sigma$.

##  IV Experiments and Results

We demonstrate polarization-entanglement of photons with our entangled photon
source, with the hybrid interferometer implemented (as shown in Figure 4).

Figure 4: Design of the EPS. The optical components are : Filter (F), Optical
Isolator (OA), Linear Polarizer (LP), Lens (L), Beam Displacer (BD), Half
waveplate (HWP), Dichroic Mirror (DM), Polarizing Beam splitter (PBS)

The pump beam collimated by an aspheric lens ($f_{p}=$8\text{\,}\mathrm{mm}$$)
passes through several optical components to be diagonally polarized. It is
then focused by the magnifying lens ($f_{p}=$40\text{\,}\mathrm{cm}$$) to the
PPLN crystals in the Sagnac loop. The first beam displacer splits the
diagonally polarized pump beam into two beams with horizontal and vertical
polarization by displacing one of them that is extraordinary with respect to
the displacer’s optic axis. The polarization of the two beams is switched by
passing through a half waveplate before they enter the Sagnac loop. The first
dichroic mirror reflects the pump beam into the Sagnac loop, while allowing
the signal and idler beams that are produced from down conversion by the PPLN
crystals to pass through and reach the second beam displacer. We place two
identical PPLN crystals perpendicular to each other such that one is
responsible for $\ket{o}_{p}\rightarrow\ket{o}_{s}\ket{o}_{i}$ while the other
is responsible for $\ket{e}_{p}\rightarrow\ket{e}_{s}\ket{e}_{i}$. The two
signal and idler beams combined by the second beam displacer are separated by
the second dichroic mirror (signal beam is transmitted and idler beam is
reflected). Along each beam path, there is another lens
($f_{s}=$60\text{\,}\mathrm{cm}$$ and $f_{i}=$50\text{\,}\mathrm{cm}$$) to re-
collimate the beam before entering the Calcite wedges for the walk-off
compensation. It is worth noting that the lenses are placed before the wedges
in order to prevent from the beam entering non-perpendicular to the wedges’
surfaces. At this point, we would expect the signal and the idler photons to
be completely entangled. A long-pass filter and a bandpass filter are placed
on each beam path to ensure that any noises are reduced. The polarization-
entanglement test is proceeded by measuring the polarization correlation of
the Bell states. A half-wave plate and a polarizing beam splitter are placed
on the end of each beam path (signal and idler). The arrived signal and idler
photons are focused by the aspheric lenses ($f_{s}=$11\text{\,}\mathrm{cm}$$
and $f_{i}=$15\text{\,}\mathrm{cm}$$) and collected through the SMFs. SPAD
(Exceltilas, USA) and SNSPD (Quantum Opus, USA) are used to collect signal and
idler photons, respectively. The detection efficiencies from the signal and
the idler beam paths are estimated to be $\eta_{s}=0.38$ and $\eta_{i}=0.217$
considering all the optical components along the ways including the detectors
(These values are estimated based on the specifications of the optical
components provided by the manufacturers). The coincidence counts within a
$1.5\text{\,}\mathrm{ns}$ detection window are measured using a time-tagging
unit.

Figure 5: (a) Spatial walk-off compensation. While fixing the signal side
(fully compensated), we shift the position of one of the Calcite wedges on the
idler side, as shown in the diagram, varying the gap. The optical fiber
collecting idler photons is realigned, and the thickness is simultaneously
corrected for each gap variation. (b) Theoretical curve of the spatial overlap
for idler photons (solid line) as a function of lateral separation, along with
the calculated coincidences (dot) for each lateral separation. Note that both
the vertical and horizontal error bars are too short to be visible. (c)
Temporal walk-off compensation. While fixing the idler side (fully
compensated), we shift the position of one of the Calcite wedges on the signal
side, as shown in the diagram, varying only the thickness. (d) Theoretical
curve of the temporal overlap for signal photons (solid line) as a function of
compensation length, along with the calculated visibility (dot) for each
compensation length. Note that the vertical error bars are too short to be
visible.

In actual experiments, higher values of $\Delta
D_{s}=$0.145\text{\,}\mathrm{mm}$$ and $\Delta
D_{i}=$0.325\text{\,}\mathrm{mm}$$ are observed due to additional divergence
at the two dichroic mirrors. Nevertheless, compensating for the spatial walk-
off should pose no issues with the additional divergence and should simply be
taken into account. By scanning the beams in the transverse directions with
our detectors, we identify the FWHM of the signal and idler beams at the
receiving fibers to be $0.6\text{\,}\mathrm{mm}$ and
$0.8\text{\,}\mathrm{mm}$, respectively. With the actual walk-offs obtained
from the experiment, the spatial overlap factors calculated from equation 18
are 52.4% for the signal photon, and 17.1% for the idler photon. Using
equation 13, the lateral separation needed for the signal and idler paths are
calculated to be $d_{s}=$2.75\text{\,}\mathrm{mm}$$ and
$d_{i}=$6.6\text{\,}\mathrm{mm}$$, respectively. To validate our compensation
for the spatial walk-off, the lateral separation between the wedges on the
idler side is varied as $0\text{\,}\mathrm{mm}$, $1.6\text{\,}\mathrm{mm}$,
$3.6\text{\,}\mathrm{mm}$, $6.6\text{\,}\mathrm{mm}$,
$9.6\text{\,}\mathrm{mm}$ and $11.6\text{\,}\mathrm{mm}$, while the spatial
walk-off on the signal side is fully compensated (see Figure 5a). With the
pump power of $1.0\text{\,}\mathrm{mW}$, background subtracted coincidence
counts at each variation is measured and plotted on Figure 5b. The peak
experimental value plotted on the idler graph corresponds to when spatial
walk-off is fully compensated for both the signal and idler beams, with the
coincidence counts measuring $N=(33.33\pm 0.05)$kHz. This marks almost eight
times enhancement in the coincidence counts, as without compensation on the
idler side (while the signal side remained fully compensated), the coincidence
count was $N=(4.20\pm 0.01)$kHz. The additional beam divergence (resulting in
spatial walk-off) did introduce some ambiguities when determining temporal
walk-off compensation. If we neglect any additional temporal walk-offs other
than those caused by the beam displacers, the temporal overlaps of two
Gaussian intensity profiles, as calculated using temporal version of equation
18, are 88.8% for the signal photon, and 72.1% for the idler photon. To
validate our compensation for the temporal walk-off, the compensation length
(the overall thickness of the wedges) on the signal side is varied as
$0.68\text{\,}\mathrm{mm}$, $1.18\text{\,}\mathrm{mm}$,
$1.68\text{\,}\mathrm{mm}$, $2.18\text{\,}\mathrm{mm}$ and
$2.68\text{\,}\mathrm{mm}$ (see Figure 5c). The values are chosen simply to
plot the entanglement visibility with its optimal value centered on the graph
shown in Figure 5d. The correlations are measured by rotating the half
waveplate on the signal arm while the half waveplate on the idler arm are
fixed to $0^{o}$ (V basis), $45^{o}$ (A basis), $90^{o}$ (H basis) and
$135^{o}$ (D basis). Then the visibility
$V=(N_{c,max}-N_{c,min})/(N_{c,max}+N_{c,min})$ of each of the four bases is
obtained from the minimum and maximum coincidence counts of the correlation
measurements. The result shows the peak average value of $V=(96.6\pm 0.3)\%$,
clearly showing the enhancement in the entanglement visibility by the Calcite
wedges. Furthermore, based on the compensation results, the additional beam
divergence appears to have had minimal, if any, impact on the trend.

Figure 6: Polarization entanglement measurement results. The correlation
measurement was conducted using four different polarization bases. The
coincidence counts, measured at $1.0\text{\,}\mathrm{mW}$ of pump power while
changing the polarization angle from $0^{o}$ to $180^{o}$, are plotted on the
graph.

Figure 6 shows the optimal entanglement test results. The measured
visibilities with the full walkoff compensation are $V=(96.8\pm 0.5)\%$,
$V=(95.5\pm 0.5)\%$, $V=(98.6\pm 0.5)\%$ and $V=(95.4\pm 0.5)\%$. In order to
confirm the entanglement, we perform the CHSH test [26]. The CHSH inequality
is expressed as

$\displaystyle
S=\left|E_{\alpha,\beta}-E_{\alpha,\delta}+E_{\gamma,\delta}+E_{\gamma,\beta}\right|\leq
2$ (19)

where we have the correlation estimates expressed as

$\displaystyle
E_{\alpha,\beta}=\dfrac{N\left(\alpha,\beta\right)-N\left(\alpha,\beta^{\perp}\right)-N\left(\alpha^{\perp},\beta\right)+N\left(\alpha^{\perp},\beta^{\perp}\right)}{N\left(\alpha,\beta\right)+N\left(\alpha,\beta^{\perp}\right)+N\left(\alpha^{\perp},\beta\right)+N\left(\alpha^{\perp},\beta^{\perp}\right)}$
(20)

Here, we define additional angles, $\alpha^{\perp}=\alpha+90^{o}$ and
$\beta^{\perp}=\beta+90^{o}$. We chose $\alpha=0^{o}$, $\beta=22.5^{o}$,
$\gamma=45^{o}$ and $\delta=67,5^{o}$ and confirmed that the inequality is
strongly violated at $S=2.75\pm 0.03$ by 25 standard deviations, which also
agrees nicely with the expected value from the average visibility
$S_{exp}=2\sqrt{2}V_{avg}\approx 2.75$, indicating the quantum entanglement
[27, 26].

Figure 7: The single counts of both signal (red) and idler (blue) photons, as
well as the coincidence (green) counts are measured from four different
polarization basis while varying the pump power. After the necessary
background subtractions, the data are plotted on these graphs. Graphs showing
the photon pair emission rates $R_{pair}$ (green), along with the heralding
efficiencies of the signal photon $\eta_{s}$ (blue) and the idler photon
$\eta_{i}$ (red) for different pump powers, calculated from the data shown in
(a).

We also examine the consistency of experimental data quality with increasing
pump power levels. The experimental data includes the single counts of
signal($N_{s}$) and idler($R_{i}$) photons, as well as coincidence counts($N$)
representing simultaneous detection of the signal and the idler photons. The
results are plotted in Figure 7a. At a pump power of
$1.0\text{\,}\mathrm{mW}$, the average background subtracted values for four
different bases were $N_{s}=(460.7\pm 3.5)$kHz, $N_{i}=(210.7\pm 1.0)$kHz and
$N=(33.33\pm 0.05)$kHz. These are the values that also corresponds to the
maximum counts shown in Figure 6. The calculated pair rates and the heralding
efficiencies are plotted in Figure 7b. The pair rate from the V-basis at
$1.0\text{\,}\mathrm{mW}$ was $R_{pair}=(2.92\pm 0.12)$MHz. Also, signal and
idler heralding efficiencies of $\eta_{s}=15.8\%$ and $\eta_{i}=7.2\%$,
respectively. The pair rate from all four bases showed a linear increase with
increasing pump power, while the heralding efficiencies remained fairly
constant as one would expect.

##  V Discussion and Outlook

The most significant advantage of our compensation method is its ability to
effectively and easily compensate for both the spatial and temporal walk-off
of photon fields at any wavelength. With our novel compensation method, the
coincidence rate of signal and idler photons at $1.0\text{\,}\mathrm{mW}$ of
pump power was $N=(33.33\pm 0.05)$kHz, showing a significant improvement from
the coincidence rate observed without the compensation which was $N=(4.20\pm
0.01)$kHz. At the same power, the average entanglement visibility of $(96.6\pm
0.3)\%$, and the average pair generation rate of $(2.92\pm 0.12)$MHz from all
four bases were observed. The trends of all these bases in Figure 7b indicate
that $100\text{\,}\mathrm{MHz}$ of pair rate could possibly be achieved at the
power level of only around $50\text{\,}\mathrm{mW}$, making it a promising
source for long distance QKD [11]. With the validation from our experiment, we
anticipate this compensation technique contributing not only to advancements
in long-distance quantum networks but also to other practical quantum
entanglement applications that require non-degenerate entangled photon pairs.

##  VI Acknowledgement

The authors acknowledge funding from the Natural Sciences and Engineering
Research Council of Canada (NSERC), Canada Foundation for Innovation (CFI),
Ontario Research Fund (ORF) and Canadian Space Agency (CSA). We thank Younseok
Lee for assistance with the purchases of PPLN crystal and the CW laser. We
thank Ramy Tannous and Henri Morin for assistance with handling equipments in
the lab.

## References

  * Tang  _et al._ [2015] Z. Tang, R. Chandrasekara, Y. C. Tan, C. Cheng, K. Durak, and A. Ling, The photon pair source that survived a rocket explosion (2015), arXiv:1512.08834 [physics.ins-det] .
  * Kim  _et al._ [2006] T. Kim, M. Fiorentino, and F. N. C. Wong, Phase-stable source of polarization-entangled photons using a polarization sagnac interferometer, Phys. Rev. A 73, 012316 (2006).
  * Shimizu  _et al._ [2008] R. Shimizu, T. Yamaguchi, Y. Mitsumori, H. Kosaka, and K. Edamatsu, Generation of polarization entanglement from spatially correlated photons in spontaneous parametric down-conversion, Phys. Rev. A 77, 032338 (2008).
  * Chapman  _et al._ [2022] J. C. Chapman, C. C. Lim, and P. G. Kwiat, Hyperentangled time-bin and polarization quantum key distribution, Phys. Rev. Appl. 18, 044027 (2022).
  * Shen  _et al._ [2022] L. Shen, C. H. Chow, J. Y. X. Peh, X. J. Yeo, P. K. Tan, and C. Kurtsiefer, Distributing polarization-entangled photon pairs with high rate over long distances through standard telecommunication fiber, Physical Review Applied 18, 10.1103/physrevapplied.18.044075 (2022).
  * Sidhu  _et al._ [2021] J. S. Sidhu, S. K. Joshi, M. Gündoğan, T. Brougham, D. Lowndes, L. Mazzarella, M. Krutzik, S. Mohapatra, D. Dequal, G. Vallone, P. Villoresi, A. Ling, T. Jennewein, M. Mohageg, J. G. Rarity, I. Fuentes, S. Pirandola, and D. K. L. Oi, Advances in space quantum communications, IET Quantum Communication 2, 182–217 (2021).
  * Ursin  _et al._ [2007] R. Ursin, F. Tiefenbacher, T. Schmitt-Manderbach, H. Weier, T. Scheidl, M. Lindenthal, B. Blauensteiner, T. Jennewein, J. Perdigues, P. Trojek, B. Ömer, M. Fürst, M. Meyenburg, J. Rarity, Z. Sodnik, C. Barbieri, H. Weinfurter, and A. Zeilinger, Entanglement-based quantum communication over 144km, Nature Physics 3, 481–486 (2007).
  * Hentschel  _et al._ [2009] M. Hentschel, H. Hübel, A. Poppe, and A. Zeilinger, Three-color sagnac source of polarization-entangled photon pairs, Optics Express 17, 23153 (2009).
  * Li  _et al._ [2023] Z. Li, K. Xue, J. Li, L. Chen, R. Li, Z. Wang, N. Yu, D. S. L. Wei, Q. Sun, and J. Lu, Entanglement-assisted quantum networks: Mechanics, enabling technologies, challenges, and research directions, IEEE Communications Surveys amp; Tutorials 25, 2133–2189 (2023).
  * Ma  _et al._ [2007] X. Ma, C.-H. F. Fung, and H.-K. Lo, Quantum key distribution with entangled photon sources, Physical Review A 76, 10.1103/physreva.76.012307 (2007).
  * Bourgoin  _et al._ [2013] J.-P. Bourgoin, E. Meyer-Scott, B. L. Higgins, B. Helou, C. Erven, H. Hübel, B. Kumar, D. Hudson, I. DSouza, R. Girard, R. Laflamme, and T. Jennewein, A comprehensive design and performance analysis of low earth orbit satellite quantum communication, New Journal of Physics 15, 023006 (2013).
  * Esmaeil Zadeh  _et al._ [2017] I. Esmaeil Zadeh, J. W. N. Los, R. B. M. Gourgues, V. Steinmetz, G. Bulgarini, S. M. Dobrovolskiy, V. Zwiller, and S. N. Dorenbos, Single-photon detectors combining high efficiency, high detection rates, and ultra-high timing resolution, APL Photonics 2, 111301 (2017), https://pubs.aip.org/aip/app/article-pdf/doi/10.1063/1.5000001/14567918/111301_1_online.pdf .
  * Yang  _et al._ [2024] M. Yang, S.-K. Liao, W.-S. Tang, W.-Q. Cai, Y. Li, J. Yin, and C.-Z. Peng, In-orbit dark count rate performance and radiation damage high-temperature annealing of silicon avalanche photodiode single-photon detectors of the micius satellite, Opt. Express 32, 12601 (2024).
  * [14] J. B. Abshire, X. Sun, H. Riris, J. M. Sirota, J. F. McGarry, S. Palm, D. Yi, and P. Liiva, Geoscience laser altimeter system (glas) on the icesat mission: On-orbit measurement performance, Geophysical Research Letters 32, https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/2005GL024028 .
  * DSouza  _et al._ [2021] I. DSouza, J.-P. Bourgoin, B. L. Higgins, J. G. Lim, R. Tannous, S. Agne, B. Moffat, V. Makarov, and T. Jennewein, Repeated radiation damage and thermal annealing of avalanche photodiodes, EPJ Quantum Technology 8, 10.1140/epjqt/s40507-021-00103-0 (2021).
  * Ribezzo  _et al._ [2023] D. Ribezzo, M. Zahidy, G. Lemmi, A. Petitjean, C. De Lazzari, I. Vagniluca, E. Conca, A. Tosi, T. Occhipinti, L. K. Oxenløwe, A. Xuereb, D. Bacco, and A. Zavatta, Quantum key distribution over 100 km of underwater optical fiber assisted by a fast-gated single-photon detector, Phys. Rev. Appl. 20, 044052 (2023).
  * Jennewein  _et al._ [2014] T. Jennewein, J. Bourgoin, B. Higgins, C. Holloway, E. Meyer-Scott, C. Erven, B. Heim, Z. Yan, H. Hübel, G. Weihs, E. Choi, I. D’Souza, D. Hudson, and R. Laflamme, Qeyssat: A mission proposal for a quantum receiver in space, Proceedings of SPIE - The International Society for Optical Engineering 8997 (2014).
  * Pugh  _et al._ [2017] C. J. Pugh, S. Kaiser, J.-P. Bourgoin, J. Jin, N. Sultana, S. Agne, E. Anisimova, V. Makarov, E. Choi, B. L. Higgins, and T. Jennewein, Airborne demonstration of a quantum key distribution receiver payload, Quantum Science and Technology 2, 024009 (2017).
  * [19] Qeyssat mission - canadian space agency, online source (last accessed in april 2024) : $www.asc-csa.gc.ca/eng/satellites/qeyssat.asp$.
  * Steinlechner  _et al._ [2012] F. Steinlechner, P. Trojek, M. Jofre, H. Weier, D. Perez, T. Jennewein, R. Ursin, J. Rarity, M. W. Mitchell, J. P. Torres, H. Weinfurter, and V. Pruneri, A high-brightness source of polarization-entangled photons optimized for applications in free space, Optics Express 20, 9640 (2012).
  * Gayer  _et al._ [2008] O. Gayer, Z. Sacks, E. Galun, and A. Arie, Temperature and wavelength dependent refractive index equations for MgO-doped congruent and stoichiometric LiNbO3, Applied Physics B: Lasers and Optics 91, 343 (2008).
  * Makarov [2022] D. N. Makarov, Theory for the beam splitter in quantum optics: quantum entanglement of photons and their statistics, hom effect (2022), arXiv:2211.03359 [quant-ph] .
  * Shi  _et al._ [2020] Y. Shi, S. Moe Thar, H. S. Poh, J. A. Grieve, C. Kurtsiefer, and A. Ling, Stable polarization entanglement based quantum key distribution over a deployed metropolitan fiber, Applied Physics Letters 117, 10.1063/5.0021755 (2020).
  * Lee  _et al._ [2020] Y. S. Lee, M. Xie, R. Tannous, and T. Jennewein, Sagnac-type entangled photon source using only conventional polarization optics (2020), arXiv:2007.05095 [quant-ph] .
  * Lohrmann  _et al._ [2020] A. Lohrmann, C. Perumangatt, A. Villar, and A. Ling, Broadband pumped polarization entangled photon-pair source in a linear beam displacement interferometer, Applied Physics Letters 116, 10.1063/1.5124416 (2020).
  * Clauser  _et al._ [1969] J. F. Clauser, M. A. Horne, A. Shimony, and R. A. Holt, Proposed experiment to test local hidden-variable theories, Phys. Rev. Lett. 23, 880 (1969).
  * Peruzzo and Sorella [2023] G. Peruzzo and S. Sorella, Entanglement and maximal violation of the chsh inequality in a system of two spins j: A novel construction and further observations, Physics Letters A 474, 128847 (2023).


Consensus-Based Distributed Estimation in the presence of Heterogeneous, Time-Invariant Delays

Mohammadreza  Doostmohammadian, Member, IEEE,    Usman A. Khan, Senior Member, IEEE, Mohammad Pirani, Member, IEEE, and Themistoklis Charalambous, Senior Member, IEEE

M.  Doostmohammadian and T. Charalambous are with the School of Electrical Engineering at Aalto University, Finland (name.surname@aalto.fi). M.  Doostmohammadian is also with the Faculty of Mechanical Engineering at Semnan University, Iran (doost@semnan.ac.ir).
U. A. Khan is with the Electrical and Computer Engineering Department at Tufts University, MA, USA (khan@ece.tufts.edu).
M. Pirani is with the Department of Electrical and Computer Engineering, University of Waterloo, Canada (mpirani@uwaterloo.ca).

Classical distributed estimation scenarios typically assume timely and reliable exchanges of information over the sensor network.
This paper, in contrast, considers single time-scale distributed estimation via a sensor network subject to transmission time-delays. The proposed discrete-time networked estimator consists of two steps: (i) consensus on (delayed) a-priori estimates, and (ii) measurement update.  The sensors only share their a-priori estimates with their  out-neighbors  over (possibly) time-delayed transmission links. The delays are assumed to be fixed over time, heterogeneous, and known. We assume distributed observability instead of local observability, which significantly reduces the communication/sensing loads on sensors. Using the notions of augmented matrices and Kronecker product, the convergence of the proposed estimator over strongly-connected networks is proved for a specific upper-bound on the time-delay.

istributed estimation, consensus, Kronecker product, communication time-delays

§ INTRODUCTION

Latency in data transmission networks may significantly affect the performance of decision-making over sensor networks and multi-agent systems <cit.>. In particular, time-delays may cause instability in networked control systems which are originally stable in the corresponding delay-free case. For example, the consequence of communication delays on the consensus stability  are discussed in [1, 2, 3] among others,

and centralized observer design

are  discussed in [4, 5]. This work extends to distributed estimation over a sensor network with random communication time-delays.

The literature on distributed estimation spans from multi time-scale scenarios to single time-scale methods.

The former case  requires many iterations of averaging/data-sharing (consensus/communication time-scale) between two consecutive system time-steps (system time-scale) [6, 7], where the estimation performance tightly depends on the number of consensus iterations. This  is less efficient in terms of computational and communication loads on sensors and, further, requires much faster data sharing/processing rate which might be inaccessible over large networks. In terms of observability, in the multi time-scale method, number of communication/consensus iterations is greater than the network diameter, and therefore, all sensors eventually gain all state information (and system observability)  between every two system time-steps. In the single time-scale, however, every sensor performs only one iteration of consensus, and therefore, many works require the system to be locally observable in the neighborhood of the sensors [8, 9, 10, 11, 12, 13, 14, 15]; in this work, we assume global observability as in [16, 17, 18, 19]. Recall that local observability mandates: (i) more network connectivity, and/or (ii)  access to more system outputs at each sensor, and may considerably increase the communication/sensing-related costs [20, 21, 22]. This work, however, considers least connectivity requirement   (strong-connectivity) and least outputs at each sensor (one output), while addressing transmission delays.

The networked estimator in this paper is single time-scale, where sensors perform one consensus iteration on (possibly) delayed a-priori estimates in their in-neighborhood, and then,  measurement-update using their own outputs.
As in[1], we consider  arbitrary time-delays  at every communication link, but the delays are time-invariant and known. The delays are bounded so that no information is lost over the network and the data would eventually reach the recipient sensor.
To avoid considering a trivial case, this work makes no assumption on the stability of the  linear system. Further, similar to [9, 11], we assume that the system is full-rank as in structurally-cyclic/self-damped systems [23, 24]. We adopt the notions of augmented representation [1] and the Kronecker network product [25] to simplify the convergence analysis. We show that feedback gain design in the absence of delays via the Linear-Matrix-Inequality (LMI) in [16, 26] also results in  stable  estimation  for  some upper-bounded delayed cases. Further, we provide a solution to design delay-tolerant networked estimators for a given bound on the delays. Therefore, the gain design requires no information other than the bound  on the delays and, the LMI complexity is determined by the original low-order system and not the high-order augmented one. Note that in this work, we assume no measurement delays; this is because the sensors take direct measurements, while spatially distributed in large-scale with (possibly) delayed communications. Further, as proved in [5], stability  depends only on the measurement packet loss, not the packet delays.

Paper organization: Section <ref> provides some preliminaries and problem statement. Section <ref> states our main results on delay-tolerant distributed estimation. Section <ref> provides the simulations, and Section <ref> concludes the paper.

§ THE FRAMEWORK

§.§ System-Output Model

We consider discrete-time system and measurements as,

\begin{align}\label{eq_sys}
\mb{x}_{k} &= A \mb{x}_{k-1} + \mb{\nu}_{k},~k\geq 1\\
\mb{y}_{k} &= C \mb{x}_{k} + \mb{\zeta}_k, \label{eq_y}
\end{align}

with $\mb{x}_k \in \mathbb{R}^n$ as system states, $\mb{y}_k \in \mathbb{R}^N$ as system outputs, $\mb{\zeta}_k\sim\mc{N}(0,R)$ and $ \mb{\nu}_k \sim \mc{N}(0,Q)$ as independent  noise variables, all at time-step $k$. It is not assumed that $\rho(A)<1$ (potentially unstable system), while $\det (A)\neq 0$ (full-rank),
where $\rho(\cdot)$ and $\det(\cdot)$  are the spectral radius and determinant, respectively. Examples of such full-rank systems are structurally cyclic [24] and self-damped[A linear system  is self-damped if its matrix $A$ has non-zero diagonal entries [24, 23]. Similarly, network $\mc{G}=\{\mc{V},\mc{E}\}$ is self-damped if for every node $i \in \mc{V}$ we have $(i,i) \in \mc{E}$, i.e., there is a self-link at every node.] systems [24, 23] among others, prevalent in, e.g., social opinion dynamics [27].
In this work, without loss of generality, we assume $N$ sensors each with one output ${y}_{k}^i$.

Every sensor knows the system matrix $A$. The pair $(A,C)$ is observable (similar arguments for detectable case), implying global observability. However, in general, the pair $(A,  C_i)$ is not necessarily observable at any sensor $i$.

Note that the output matrix $C$ can be defined via the graph-theoretic methods  to ensure (structural)  $(A,C)$-observability. From [16, 28, 22], having one output from (at least) one state associated with every irreducible block of the adjacency matrix $A$ ensures structural observability. The optimal output selection strategies are also of interest as in [29, 24, 30].

§.§ Preliminaries on Consensus Algorithms

Consider discrete-time consensus algorithms over a network of sensors $\mc{G}=\{\mc{V},\mc{E}\}$ ($\mc{V}$ as the node set and $\mc{E}$ as the link set) with $\mb{z}_{k}$ as the state of  sensors at time $k$, which evolves as $\mb{z}_{k} = P \mb{z}_{k-1}$. Matrix $P$ (as the consensus weight) represents the communication between the sensors via graph $\mathcal{G}$.

The sensor network is in general directed. For notation simplicity denote $P(i,j)$ by $p_{ij}$, where $0<p_{ij}<1$ if $(j,i) \in \mc{E}$ and $0$ otherwise. $P$  is row-stochastic, i.e., $\sum_{j=1}^N p_{ij} = 1$, and $p_{ii} \neq 0$ for all $i$. Further,  the network $\mc{G}$ needs to be (at least) strongly-connected (SC), i.e., there is a path from every node $i$ to every node $j$ in $\mathcal{V}$, implying  that fusion matrix $P$ is irreducible, also called stochastic, indecomposable, and aperiodic (SIA) and  $\lim_{k \rightarrow \infty} P^k = \mb{d}\mb{1}_N^\top
$ [1], with $\mathbf{1}_N$ as all-ones vector of size $N$. For such SIA matrix,  $\rho(P)=1$.

§.§ Delay Model

In this work, it is assumed that the data-transmission over the link $(j,i)$ from sensor $j$ to sensor $i$  has a-priori unknown bounded (integer)  time-delay, $\tau_{ij}$, where $0\leq \tau_{ij} \leq \overline{\tau} <\infty$, and $\overline{\tau}$ is an upper bound to the delays in all links. The messages are time-stamped, so the recipient knows the time-step the data was sent.
Further, $\tau_{ii}=0$, i.e., every sensor $i$ knows its own state with no delay. To model the delayed state vectors we adopt the notations in [1]. In a network of $N$ sensors, define an augmented state vector $\underline{\mb{x}}_k = \left(\mb{x}_k; \mb{x}_{k-1}; \dots; \mb{x}_{k-\overline{\tau}} \right)$  with '$;$' as column concatenation, and $\mb{x}_{k-r} = \left({x}^1_{k-r};  \dots; {x}^n_{k-r} \right)$ for $0\leq r \leq \overline{\tau}$. Then, for a given $N$-by-$N$ matrix $P$ and maximum delay $\overline{\tau}$, define the augmented matrix $\overline{P}$ as,

\begin{align} \label{eq_aug_W} \small
\overline{P}= \left(
\begin{array}{cccccc}
P_0 & P_1 & P_{2} & \hdots & P_{\overline{\tau}-1} & P_{\overline{\tau}} \\
I_N &   0_N & 0_N &\hdots  & 0_N& 0_N\\
0_N & I_N & 0_N &  \hdots  & 0_N & 0_N  \\
0_N &  0_N & I_N  &  \hdots  & 0_N & 0_N  \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
0_N & 0_N & 0_N & \hdots & I_N & 0_N
\end{array}
\right), \normalsize
\end{align}

with $I_{N}$ and $0_{N}$ as the identity and zero matrix of size $N$. The non-negative matrices $P_r$ are defined based on the time delay $0\leq r \leq \overline{\tau} $ on the network links as follows,

\begin{align}
P_r(i,j) = \left\{
\begin{array}{ll}
p_{ij}, & \text{If}~ \tau_{ij}=r  \\
0, & \text{Otherwise}.
\end{array}\right.
\end{align}

Assuming fixed delays, for any $(j,i) \in \mc{E}$ in the given communication network $\mc{G}$, only one of $P_0(i,j),P_1(i,j), \hdots, P_{\overline{\tau}}(i,j)$ is  equal to  $p_{ij}$  and the rest are zero. This implies that the row-sum of each of the first $N$ rows of $\overline{P}$ and $P$ are equal, i.e., $\sum_{j=1}^{N(\overline{\tau}+1)} \overline{P}(i,j) =  \sum_{j=1}^{N} {P}(i,j)$ for $1\leq i \leq N$ and $P = \sum_{r=0}^{\overline{\tau}} {P}_r$ for $k \geq 0$. Therefore, in case of having a row-stochastic matrix $P$, the augmented matrix $\overline{P}$ is also row-stochastic[Note the subtle difference between our notation vs. Ref. [1]. In [1] column augmented matrix is introduced, while we consider row augmented matrices.]. In the proposed estimator, we do not need the matrix $\overline{P}$ and it is only defined to simplify the mathematical analysis. Similar to [1], the knowledge of the probability distribution of $\tau_{ij}$s is not needed for the analysis in this paper.

For the time-delay $\tau_{ij}$ on  link $(j,i)$:

The delay is known and bounded $\tau_{ij} \leq \overline{\tau}$.
The upperbound $\overline{\tau}$ guarantees no lost information, i.e., the data sent from sensor $j$ at time
$k$ would eventually reach the recipient sensor $i$ (at most) at time $k+\overline{\tau}$ (a $\overline{\tau}+1$-slot transmission buffer).

* Delay $\tau_{ij}$ is arbitrary, fixed, and  may/may-not differ for  $(j,i)$ links (heterogeneous/homogeneous delays).

§.§ Problem Statement

The problem in this work is to design a networked estimator for the system-output model (<ref>)-(<ref>) satisfying Assumption <ref>,

where every sensor relies only on its partial system output (partial observability) and the received  (possibly delayed satisfying Assumption <ref>) information from its in-neighbors. This work particularly differs from  [8, 9, 10, 11, 12, 13, 14, 15] via the following remark.

Let $\mc{N}_i=\{j|(j,i) \in \mc{E}\}$ denote the set of in-neighborhood of  sensor $i$ over the network $\mc{G}$. The pair $(A, \sum_{j \in \mc{N}_i} C_j)$ is not necessarily observable at any sensor $i$, implying no local observability assumption.

§ DISTRIBUTED ESTIMATION IN PRESENCE OF DELAYS

Every sensor $i$ performs the following two steps for distributed state estimation in the presence of time-delays,

\begin{align}\label{eq_p}
\widehat{\mb{x}}^i_{k|k-1} =& p_{ii}A\widehat{\mb{x}}^i_{k-1|k-1} + \sum_{j\in\mathcal{N}_i} \sum_{r=0}^{\overline{\tau}} p_{ij}A^{r+1}\widehat{\mb{x}}^j_{k-r|k-r} \mb{I}_{k-r,ij}(r),
\\\label{eq_m}
\widehat{\mb{x}}^i_{k|k} =& \widehat{\mb{x}}^i_{k|k-1} + K_i C_i^\top \left({y}^i_k-C_i\widehat{\mb{x}}^i_{k|k-1}\right),
\end{align}
where  $\mb{I}_{k,ij}(r)$ is the indicator function  defined as [1],

\begin{align}
\mb{I}_{k,ij}(r) = \left\{
\begin{array}{ll}
1, & \text{if}~ \tau_{ij}=r  \\
0, & \text{otherwise}.
\end{array}\right.
\end{align}

In (<ref>), $\widehat{\mb{x}}^i_{k|k-1}$ denotes the sensor $i$'s a-priori state estimate at time $k$ given all the (possibly delayed) information up to time $k-1$ from its in-neighbors $\mathcal{N}_i$. Step (<ref>) represents one iteration of consensus-based information-fusion on all the received information, where sensor $i$ sums the weighted estimates
of sensors $j\in \mc{N}_i$ as they arrive  knowing the delays. The performance analysis of (<ref>)-(<ref>) in terms of mean-square stability (for the delay-free case) is given in [16]. Recall that, for every link $(i,j)$ the indicator $\mb{I}_{k-r,ij}(r)$ is only non-zero for one $r$ between $0$ and $\overline{\tau}$ (due to fixed delay assumption). The second step (<ref>) is a measurement-update (also known as innovation) to modify the a-priori estimate based on the new measurement of sensor $i$. Clearly, the  protocol (<ref>)-(<ref>) is single time-scale with one step of information-sharing/consensus-update between every $k-1$ and $k$.
Using the notion of augmented vector, define $\underline{\widehat{\mb{x}}}_{k|k-1} = \left(\widehat{\mb{x}}_{k|k-1}; \widehat{\mb{x}}_{k-1|k-2}; \dots; \widehat{\mb{x}}_{k-\overline{\tau}|k-\overline{\tau}-1} \right)$ and similarly $\underline{\widehat{\mb{x}}}_{k|k}$. Then, the augmented version of (<ref>)-(<ref>) is,

\begin{align}\label{eq_p_aug}
\underline{\widehat{\mb{x}}}_{k|k-1} =&\overline{PA} \underline{\widehat{\mb{x}}}_{k-1|k-1},
\\\label{eq_m_aug}
\underline{\widehat{\mb{x}}}_{k|k} =&  \underline{\widehat{\mb{x}}}_{k|k-1} + \mb{b}^{\overline{\tau}+1}_1 \otimes K D_C^\top \left(\mb{y}_k-D_C\Xi^{Nn}_{1,\overline{\tau}}\underline{\widehat{\mb{x}}}_{k|k-1}\right),
\end{align}

$\overline{PA}$ is the modified augmented version of $P \otimes A$ as,

\begin{align} \label{eq_aug_WA}
\overline{PA}= \left(
\begin{array}{cccccc}
P_0 \otimes A & P_1\otimes A^2  &  \hdots & P_{\overline{\tau}-1}\otimes A^{\overline{\tau}}  & P_{\overline{\tau}} \otimes A^{\overline{\tau}+1}  \\
I_{Nn} &   0_{Nn}  &\hdots  & 0_{Nn}& 0_{Nn}\\
0_{Nn} & I_{Nn} &   \hdots  & 0_{Nn} & 0_{Nn}  \\
0_{Nn} &  0_{Nn} &  \ddots  & 0_{Nn} & 0_{Nn}  \\
\vdots & \vdots &  \hdots & \vdots & \vdots \\
0_{Nn} & 0_{Nn} &  \hdots & I_{Nn} & 0_{Nn}
\end{array}
\right),
\end{align}
${D}_C=\mbox{blockdiag}(C_i)$, $K = \mbox{blockdiag}(K_i)$, and the  auxiliary matrix $\Xi^m_{i,\overline{\tau}}$ is an $m \times (\overline{\tau}+1)m $ matrix defined as
$\Xi^m_{i,\overline{\tau}}= (\mb{b}^{\overline{\tau}+1}_i \otimes I_m)^\top $
with $\mb{b}^{\overline{\tau}+1}_i$ as the unit column-vector of the $i$'th coordinate ($1\leq i \leq {\overline{\tau}+1}$).

Define the augmented state  $\underline{\mb{x}}_{k} = \left(1_{N} \otimes   \mb{x}_k; 1_{N} \otimes \mb{x}_{k-1}; \dots; 1_{N} \otimes \mb{x}_{k-\overline{\tau}} \right) $ and the augmented error $\underline{\mb{e}}_{k}$ at time $k$ as follows,

\begin{align}
\underline{\mb{e}}_{k}  &= \underline{\mb{x}}_{k} - \underline{\widehat{\mb{x}}}_{k|k}  =\underline{\mb{x}}_{k} - \Bigl(\underline{\widehat{\mb{x}}}_{k|k-1} + \mb{b}^{\overline{\tau}+1}_1 \otimes K D_C^\top \mb{y}_k-D_C\Xi^{Nn}_{1,\overline{\tau}}\underline{\widehat{\mb{x}}}_{k|k-1}\Bigr) \nonumber
\\
&= \underline{\mb{x}}_{k}  - \overline{PA}\underline{\widehat{\mb{x}}}_{k-1|k-1} - \mb{b}^{\overline{\tau}+1}_1 \otimes K D_C^\top \left(\mb{y}_k-D_C\Xi^{Nn}_{1,\overline{\tau}}\overline{PA} \underline{\widehat{\mb{x}}}_{k-1|k-1}\right). \label{eq_love}
\end{align}

Define $\widetilde{\mb{\nu}}_{k} =  1_{N} \otimes \mb{\nu}_{k} $

and $\underline{\mb{\nu}}_{k} =\mb{b}^{\overline{\tau}+1}_i \otimes\widetilde{\mb{\nu}}_{k}$.
We have $D_C^\top \mb{y} =D_C^\top D_C(1_{N} \otimes\mb{x}_k)+  D_C^\top \mb{\zeta} $.
Recall that the row-stochasticity of $P$ and $\overline{P}$ along with the system dynamics (<ref>) implies that $\underline{\mb{x}}_{k}   = \overline{PA} \underline{\mb{x}}_{k-1} +  \underline{\mb{\nu}}_{k}  $; substituting this  along with (<ref>)-(<ref>) into (<ref>),

\begin{align}
\underline{\mb{e}}_{k} &= \overline{PA} \underline{\mb{x}}_{k-1} + \underline{\mb{\nu}}_{k} -  \overline{PA}\underline{\widehat{\mb{x}}}_{k-1|k-1} -  \mb{b}^{\overline{\tau}+1}_1 \otimes K D_C^\top \left(D_C(1_{N} \otimes\mb{x}_k)+  \mb{\zeta}_{k}-D_C\Xi^{Nn}_{1,\overline{\tau}}\overline{PA}\underline{\widehat{\mb{x}}}_{k-1|k-1}\right) \nonumber \\
%	&= \overline{PA} \underline{\mb{x}}_{k-1} + \underline{\mb{\nu}}_{k}
%	-  \overline{PA}\underline{\widehat{\mb{x}}}_{k-1|k-1} - \mb{b}^{\overline{\tau}+1}_1 \otimes K D_C^\top \Bigl(D_C(1_{N} \otimes A\mb{x}_{k-1}) + D_C\widetilde{\mb{\nu}}_{k}+  \mb{\zeta}_k
%	- D_C\Xi^{Nn}_{1,\overline{\tau}}\overline{PA}\underline{\widehat{\mb{x}}}_{k-1|k-1}\Bigr) \nonumber \\
&= \overline{PA} \underline{\mb{e}}_{k-1}
- \mb{b}^{\overline{\tau}+1}_1 \otimes K \overline{D}_C \Bigl(1_{N} \otimes\mb{x}_{k}
- \Xi^{Nn}_{1,\overline{\tau}}\overline{PA}\underline{\widehat{\mb{x}}}_{k-1|k-1}\Bigr) + \underline{\mb{\eta}}_k \nonumber  \\
&= \overline{PA} \underline{\mb{e}}_{k-1}
- \mb{b}^{\overline{\tau}+1}_1 \otimes K \overline{D}_C \Bigl(\Xi^{Nn}_{1,\overline{\tau}}\overline{PA}\underline{\mb{x}}_{k-1}
- \Xi^{Nn}_{1,\overline{\tau}}\overline{PA}\underline{\widehat{\mb{x}}}_{k-1|k-1}\Bigr) + \underline{\mb{\eta}}_k \nonumber  \\
&= \overline{PA} \underline{\mb{e}}_{k-1}
- \mb{b}^{\overline{\tau}+1}_1 \otimes K \overline{D}_C \Xi^{Nn}_{1,\overline{\tau}}\overline{PA}\underline{\mb{e}}_{k-1}  + \underline{\mb{\eta}}_k   = \underline{\widehat{A}} \underline{\mb{e}}_{k-1}  + \underline{\mb{\eta}}_k, \label{eq_err1}
\end{align}

where $\underline{\widehat{A}}= \overline{PA}
- \mb{b}^{\overline{\tau}+1}_1 \otimes K \overline{D}_C \Xi^{Nn}_{1,\overline{\tau}}\overline{PA}$ is the closed-loop matrix, $\overline{D}_C= D_C^\top D_C$, and $\underline{\mb{\eta}}_k$ collects the noise terms as,

\begin{align}
\underline{\mb{\eta}}_k &= \underline{\mb{\nu}}_{k} - \mb{b}^{\overline{\tau}+1}_1 \otimes K \overline{D}_C\widetilde{\mb{\nu}}_{k} -\mb{b}^{\overline{\tau}+1}_1\otimes  K {D}_C^\top\mb{\zeta}_{k} =  \mb{b}^{\overline{\tau}+1}_1 \otimes (\widetilde{\mb{\nu}}_{k} - K \overline{D}_C\widetilde{\mb{\nu}}_{k} - K {D}_C^\top\mb{\zeta}_{k}).
\label{eq_eta}
\end{align}

The error dynamics in the absence of any delay is as follows,

\begin{align} \nonumber
\mb{e}_{k} &= (P\otimes A - K \overline{D}_C (P\otimes A))\mb{e}_{k-1} +
\mb{\eta}_k \\
&= \widehat{A} \mb{e}_{k-1} + 	\mb{\eta}_k
\label{eq_err_nodelay}
\end{align}

where  $\mb{\eta}_k$ follows  the formulation (<ref>) with $\overline{\tau}=0$ and $\widehat{A}=P\otimes A - K \overline{D}_C (P\otimes A)$ is the  delay-free closed-loop matrix.
For the Schur stability of the error dynamics (<ref>) and (<ref>), we need $\rho(\underline{\widehat{A}})<1$ and $\rho(\widehat{A})<1$, respectively. We first discuss the condition for Schur stability of $\widehat{A}$ and then extend the results to Schur stability of $\underline{\widehat{A}}$. Following Kalman theorem  and  justification in [16, 26], for stability of (<ref>) the pair $(P\otimes A, \overline{D}_C)$ needs to be observable (or detectable); this is known as distributed observability [16], discussed next.

Given a full-rank matrix $A$ and output matrix $C$, following Assumption <ref> and Remark <ref>, the pair $(P\otimes A, \overline{D}_C) $ is (structurally) observable if the matrix $P$ is irreducible.

The proof follows the results in [25] on the (structural) observability of  composite Kronecker-product networks. Given a system digraph $\mc{G}_1$ associated with full-rank system $A$ and measurement matrix $C$ satisfying $(A,C)$-observability, the minimum sufficient condition for observability of the Kronecker-product network (denoted by $\mc{G} \times \mc{G}_1$) is that $\mc{G}$ be strongly-connected and self-damped (see Theorem 4 in [25]). Following the definition of consensus matrix $P$, we have $p_{ii} \neq 0$ (satisfying the self-damped condition). The strong-connectivity of the  sensor network $\mc{G}$ is equivalent with irreducibility of matrix $P$, which completes the proof.

For observable $(P\otimes A, \overline{D}_C) $, the gain matrix $K$ can be designed such that  $\rho(\widehat{A})<1$.

§.§ Constrained Feedback Gain Design

It is known that for an observable pair $(P \otimes A, \overline{D}_C)$, the feedback gain ${K}$ can be designed to ensure Schur stability of the error dynamics (<ref>) (Corollary <ref>), i.e., $\rho(\widehat{A})<1$. Typically, such $K$ is designed via solving the following LMI,

\begin{eqnarray} \label{eq_succ}
X - \widehat{A}^\top X  \widehat{A} \succ 0,
\end{eqnarray}

for some $X\succ 0 $ with $\succ$ implying positive-definiteness. The solution of (<ref>) is, in general, a full matrix. However, for distributed estimation, we need the state feedback to be further localized, i.e., the gain matrix ${K}$ needs to be block-diagonal so every sensor uses its own state-feedback.  Such a constrained feedback gain design is proposed in [26, 16] based on cone-complementarity LMI algorithms, which are known to be of polynomial-order complexity for application in large scale.

§.§ Stability of the Delayed Estimator Dynamics

Following the Schur stability of the delay-free error dynamics (<ref>)  (via LMI design of  $K$), we extend the results to stability of the delayed dynamics (<ref>), i.e., to get  $\rho(\underline{\widehat{A}})<1$ in the presence of  delays.

Let conditions in Lemma <ref> hold and the feedback gain $K$ is designed such that $\rho({\widehat{A}})<1$ from Corollary <ref>. The networked estimator (<ref>)-(<ref>)  successfully tracks the system (<ref>) (subject to delays satisfying Assumption <ref> and possibly with $\rho({A})>1$) with stable error for any $\overline{\tau}\leq \overline{\tau}^*$, where $\overline{\tau}^* = \argmax_{\overline{\tau}} %\underset{\overline{\tau}}{\argmax}
\{\rho(P\otimes A^{\overline{\tau}+1} - K \overline{D}_C (P\otimes A^{\overline{\tau}+1}) )< 1\}$.

For the proof, following from Lemma <ref> (in the Appendix), we show that $\rho(\underline{\widehat{A}})<1$ for $\overline{\tau}\leq \overline{\tau}^*$, implying Schur stable error dynamics (<ref>). Recall that $\rho({\widehat{A}})<1$ implies that the networked estimator (<ref>)-(<ref>)  successfully tracks  system  (<ref>) for the delay-free case. From (<ref>) and Lemma <ref>, for the closed-loop matrix $\underline{\widehat{A}}$ (as modified augmented version of ${\widehat{A}}$) and $\rho({A})>1$,

\begin{align} \label{eq_tau*}
\rho(\underline{\widehat{A}})\leq \rho(P\otimes A^{\overline{\tau}+1} - K \overline{D}_C (P\otimes A^{\overline{\tau}+1}) )^{\frac{1}{\overline{\tau}+1}},
\end{align}

which implies that $\rho(\underline{\widehat{A}})<1$ for any $\overline{\tau}\leq \overline{\tau}^*$. In case $\rho({A})<1$, since  $\rho({A})^{\overline{\tau}+1}<\rho({A})$, Schur stability of $\widehat{A}$ also ensures the stability of $\underline{\widehat{A}}$ (for all $\overline{\tau}$). This completes the proof.

This theorem gives a sufficient condition for  stable tracking in the presence of  heterogeneous, time-invariant delays $\tau_{ij} = \overline{\tau}$.

§.§ Convergence Rate

Note that, in general, the exact characterization of the convergence rate/time of the linear systems is difficult. The following lemma gives the order of convergence time.

The convergence time of the error dynamics (<ref>) and (<ref>) are  of order[Given functions $f(\cdot)$ and $g(\cdot)$, say $f = O(g)$ if
$\sup_n |\frac{f(n)}{g(n)}| < \infty$ and $f = \Omega(g)$ if $g = O(f)$. If both $f = O(g)$ and
$f = \Omega(g)$ holds, then $f = \Theta(g)$ [27].]  $\Theta(\frac{1}{1-\rho(\underline{\widehat{A}})})$ and $\Theta(\frac{1}{1-\rho(\widehat{A})})$, respectively.

The proof follows from Lemma 3 in [27].

Following Lemma <ref> and <ref>, the geometric decay rate of  (<ref>) (for $\tau_{ij}=\overline{\tau}$) is  proportional to,

\begin{align} \label{eq_rate}
1-\rho(P\otimes A^{\overline{\tau}+1} - K \overline{D}_C (P\otimes A^{\overline{\tau}+1}) )^{\frac{1}{\overline{\tau}+1}}
\end{align}

For longer delays (i.e., greater $\overline{\tau},\overline{\tau}^*$ while $\overline{\tau}\leq \overline{\tau}^*$), the consensus rate in (<ref>)  and, in turn,  the decay rate of overall error (<ref>) is slower. In fact, the convergence rate  is lower-bounded by (<ref>) (for time-invariant delays).
For linear feedback systems, one can easily adjust the closed-loop eigenvalues and the convergence rate by  design of the desired feedback gain $K$. However, in the decentralized case,  $K$ is constrained to be block-diagonal. Such LMI-design with additional bound-constraint on closed-loop eigenvalues is a complex problem.

§.§ Discussions

* Each sensor processes the a-priori estimates of its in-neighbors as they arrive. The messages are time-stamped and the sensor knows the time-step (and hence the delay) of the received information. The  proposed solution works for both heterogeneous and homogeneous delays.

* The observability results are independent of the specific choice of the consensus weights (e.g., lazy Metropolis [31] or simply $\frac{1}{|\mc{N}_i|+1}$ [1])  to satisfy stochasticity of $P$. This is due to  generic/structural observability results which hold for almost all choices of numerical entries of $P$ as long as its structure (the sensor network) is fixed/time-invariant [32, 33]. Recall from structural observability that the weights for which the system is unobservable are of zero Lebesgue measure, i.e., if we choose the weights randomly, then the system will be almost surely observable. However, the structure of the consensus fusion matrix $P$ may affect  LMI-based gain design, the bound in Eq. (<ref>)$, \rho(\underline{\widehat{A}})$, and the convergence rate (<ref>). Note from Lemma <ref> (in the Appendix) that $\rho(P \otimes A)= \rho(A)$ since $\rho(P)=1$ for any choice of row-stochastic $P$.

* The cost-optimal design  of the sensor network structure and sensor placement [21, 20] can be considered to reduce the communication-related and/or sensing-related costs. In general, the cost-optimal design  subject to strong-connectivity is NP-hard. However, considering bidirectional links among sensors, it has a solution with polynomial-order complexity $\mc{O}(N^2)$ [21].

* In case of sensor failure, the concept of observational equivalence in both centralized [32] and distributed [33] scenarios can recover the loss of observability.

* To design a distributed estimator to tolerate time-delays bounded by $\overline{\tau}_1$, one can redesign the LMI gain matrix $K$ by replacing $\widehat{A}= P\otimes A^{\overline{\tau}_1+1} - K \overline{D}_C (P\otimes A^{\overline{\tau}_1+1})$ in (<ref>). Clearly, from Theorem <ref> and (<ref>), such $K$ results in $\rho(\underline{\widehat{A}})<1$ for $\overline{\tau}\leq \overline{\tau}_1$ (simply replace $\overline{\tau}^*=\overline{\tau}_1$). Such LMI gain design for the delay-free closed-loop matrix $\widehat{A}$ of size $nN$ instead of the delayed  matrix $\underline{\widehat{A}}$ of size $nN(\overline{\tau}+1)$, significantly reduces the  complexity order  with no need of using the augmented matrix $\overline{P}$.

* In <cit.>,  it is claimed that if a consensus matrix $P$ is weakly diagonally dominant, then the off-diagonal delays in $P$ are harmless for stability. However, in error dynamics (<ref>), the entries of $\underline{\widehat{A}}$ (and the weak/strong diagonal dominance of the closed-loop system) depends on the feedback gain $K$ and cannot be evaluated only based on the open-loop matrices $A$ and $P$.

§ SIMULATION

For MATLAB simulation, we consider a linear structurally-cyclic system of $n=6$ states with $4$ irreducible sub-systems and a group of $N=4$ sensors each taking one system output (from each irreducible block) satisfying Assumption <ref>. The system is full-rank and unstable with $\rho(A)=1.04$. System and output noise are considered as $\mc{N}(0,0.004)$. The network of sensors $\mc{G}$ is considered as a simple directed self-damped cycle $1 \rightarrow 2 \rightarrow 3 \rightarrow 4 \rightarrow 1$, where  the system is not observable in the neighborhood of any sensor  (Remark <ref>). With this structure, the consensus weights are considered random while satisfying row-stochasticity of $P$. It is clear that such irreducible $P$ matrix satisfies Lemma <ref> and, therefore, using the LMI strategy in [26, 16], the block-diagonal gain matrix $K$ is designed such that $\rho(\widehat{A})=0.64<1$ and $\overline{\tau}^* =10$, implying stable error dynamics for any $\overline{\tau} \leq 10$ (sufficiency from Theorem <ref>). Next, considering both  heterogeneous delays (uniformly distributed between $0$ and $\overline{\tau}$ for different links [1]) and  homogeneous delays (equal to $\overline{\tau}$ at all links) following Assumption <ref>,  the performance of the distributed estimation is analyzed. Fig. <ref> shows the Monte-Carlo simulation ($100$ trials) of mean-squared error (MSE) over the network with, (i) no time-delay, (ii) homogeneous/heterogeneous delays with $\overline{\tau}=3,8,19$  and (sufficiency) bounds in (<ref>) as  ${\rho(\underline{\widehat{A}})\leq 0.74, 0.94, 1.45}$. For fixed homogeneous delays $\overline{\tau}=19$, we have ${\rho(\underline{\widehat{A}})=1}$. From Fig. <ref> and Eq. (<ref>), longer delays decrease the MSE decay rate for $\overline{\tau}\leq \overline{\tau}^*$ (with $\overline{\tau}^*=10$ for this example), while for $\overline{\tau}> \overline{\tau}^*$ the error  may not necessarily converge.

MSEs are bounded steady-state stable for $\overline{\tau}=0,3,8$. For $\overline{\tau}=19> \overline{\tau}^*=10$ the MSE is marginally unstable for the homogeneous case, while it is stable for the heterogeneous case. Note that  Theorem <ref> gives the sufficiency condition (on $\overline{\tau}$) for convergence not necessity.

§ CONCLUSIONS AND FUTURE DIRECTIONS

This paper extends the recent literature on distributed estimation over  linear networks to time-delayed ones. We assume heterogeneous time-invariant  delays for the communication links. For a given bound on the delays, we provide a solution to design distributed estimators enabling all sensors to successfully track the system state over delayed networks.

Part of ongoing research focuses on rank-deficient systems, which adds more complexity to the problem in terms of system outputs, network connectivity, and data-sharing  [19, 28].
Other promising research directions are (i) detecting sensor faults/attacks [34, 35] along with considering latency on the distributed estimation networks, (ii) extension to time-varying delays as in [1], and (iii)

pruning the network to improve the observability properties and convergence [22, 36].

§ APPENDIX

Some of the following lemmas can be found in standard matrix theory books, e.g., in [37].

Consider two square matrices $P$ and $A$ of size $N$ and $n$, respectively, with the set of eigenvalues $\{\lambda_1,\hdots,\lambda_N\}$ and $\{\mu_1,\hdots,\mu_n\}$. Then, the set of eigenvalues of $P \otimes A$ is $\{\lambda_i \mu_j| i=1,\hdots,N,~j=1,\hdots,n \}$.

Define the following $nN$-by-${nN}$ block matrix,

\begin{align}\label{eq_aug_A}
\overline{A}_{n,i}= \left(
\begin{array}{ccccc}
0_N & \hdots & A_{i} & \hdots & 0_N \\
I_N &   0_N & \hdots &\hdots & 0_N\\
0_N &  I_N & \ddots   &  \hdots   & 0_N  \\
\vdots & \vdots & \ddots &  \ddots & \vdots \\
0_N & 0_N &  \hdots & I_N & 0_N
\end{array}
\right),
\end{align}
where $N$-by-${N}$ matrix $A_{i}$ is located at the $i$th block (and the only non-zero block) in the first block-row of $\overline{A}_{n,i}$. Let $p(\lambda)$ and $q(\overline{\lambda})$  represent the characteristic polynomials of $A_{i}$ and $\overline{A}_{n,i}$, respectively. Then, $q(\overline{\lambda}) = \overline{\lambda}^{N(n-i)}p(\overline{\lambda}^i)$.

\begin{align}
\overline{\lambda} I_{nN} - \overline{A}_{n,i}=
\left(\begin{array}{cc}
E  & F\\
G & H
\end{array}\right),
\end{align}

where block-matrix $E$ is $N(i-1)$-by-$N(i-1)$,  $F$ is $N(i-1)$-by-$N(n-i+1)$,  $G$ is $N(n-i+1)$-by-$N(i-1)$, and $H$ is $N(n-i+1)$-by-$N(n-i+1)$ defined as,

\begin{align}
E &= \left(
\begin{array}{ccccc}
\overline{\lambda} I_N & 0_N & \hdots & \hdots &  0_N \\
-I_N &  \overline{\lambda} I_N  & \hdots &\hdots & 0_N\\
0_N &  -I_N & \ddots   &  \hdots   & 0_N  \\
\vdots & \vdots & \ddots &  \ddots & \vdots \\
0_N & 0_N &  \hdots & -I_N & \overline{\lambda} I_N
\end{array}
\right)\\
F &= \left(
\begin{array}{cccc}
-A_i & 0_N & \hdots  & 0_N \\
0_N &   0_N & \hdots & 0_N\\
\vdots & \vdots & \vdots &   \vdots \\
0_N & 0_N &  \hdots  & 0_N
\end{array}
\right) G =\left(
\begin{array}{cccc}
0_N &  \hdots & 0_N & -I_N   \\
0_N & \hdots & 0_N & 0_N\\
\vdots & \vdots & \vdots & \vdots \\
0_N & 0_N &  \hdots &  0_N
\end{array}
\right) \\
H &= \left(
\begin{array}{ccccc}
\overline{\lambda} I_N & 0_N & \hdots & \hdots & 0_N \\
-I_N &  \overline{\lambda} I_N & \hdots &\hdots & 0_N\\
0_N &  -I_N & \ddots   &  \hdots   & 0_N  \\
\vdots & \vdots & \ddots &  \ddots & \vdots \\
0_N & 0_N &  \hdots & -I_N & \overline{\lambda} I_N
\end{array}
\right).
\end{align}
Recall that $p(\lambda)= |\lambda I_N - A_i|$, and

\begin{align} \label{eq_gfeh}
q(\overline{\lambda}) = |\overline{\lambda} I_{nN} - \overline{A}_{n,i}|= |E||H-GE^{-1}F|.
\end{align}

We have,

\begin{align}
E^{-1} &=& \left(
\begin{array}{ccccc}
\frac{I_N}{\overline{\lambda}} & 0_N & \hdots & \hdots & 0_N \\
\frac{I_N}{\overline{\lambda}^2} &   \frac{I_N}{\overline{\lambda}} & \hdots &\hdots & 0_N\\
\vdots & \vdots & \ddots &  \ddots & \vdots \\
\frac{I_N}{\overline{\lambda}^{i-1}} & \frac{I_N}{\overline{\lambda}^{i-2}} &  \hdots & \frac{I_N}{\overline{\lambda}^{2}} & \frac{I_N}{\overline{\lambda}}
\end{array}
\right),
\end{align}
and $H-GE^{-1}F$ is equal to,

\begin{align}
\left(
\begin{array}{ccccc}
\overline{\lambda} I_N-\frac{A_i}{\overline{\lambda}^{i-1}} & 0_N & \hdots & \hdots & 0_N \\
-I_N &  \overline{\lambda} I_N & \hdots &\hdots & 0_N\\
0_N &  -I_N & \ddots   &  \hdots   & 0_N  \\
\vdots & \vdots & \ddots &  \ddots & \vdots \\
0_N & 0_N &  \hdots & -I_N & \overline{\lambda} I_N
\end{array}
\right).
\end{align}

Then, $|E|= \overline{\lambda}^{N(i-1)}$, $ |H-GE^{-1}F|=\overline{\lambda}^{N(n-i)}|\overline{\lambda} I_N-\frac{A_i}{\overline{\lambda}^{i-1}}|$, and substituting these in (<ref>),

\begin{align}
q(\overline{\lambda}) =  \overline{\lambda}^{N(n-i)}|\overline{\lambda}^i I_N - A_i| = \overline{\lambda}^{N(n-i)} p(\overline{\lambda}^i).
\end{align}

The proof is complete.

Given matrix $A$ with $\rho(A)<1$, we have $\rho(\overline{A})\leq \rho(A)^{\frac{1}{\overline{\tau}+1}}<1$ with  $\overline{A}$ as the augmented form of $A$ via Eq. (<ref>).

The characteristic polynomial of $\overline{A}$ can be defined based on Lemma <ref>. Let  $p(\lambda)$ and $q(\overline{\lambda})$ respectively  represent the characteristic polynomial of $A$ and $\overline{A}$.

For $\tau_{ij} = \overline{\tau}$ for all $i,j$ and $A_{\overline{\tau}} = A$. Therefore, $ q(\overline{\lambda}) = p(\overline{\lambda}^{\overline{\tau}+1})$ and $	\rho(\overline{A}) = \rho(A)^{\frac{1}{\overline{\tau}+1}}<1.$
We know that the function $\rho(A)^{\frac{1}{{\tau}+1}}$ is an increasing function of $\tau$ (given $\rho(A)<1$), then,

for $\tau_{ij} =r<\overline{\tau}$ for all $i,j$ and $A_{r} = A$, we have $\rho(\overline{A}) = \rho(A)^{\frac{1}{r+1}}< \rho(A)^{\frac{1}{\overline{\tau}+1}} <1.$
This can be generalized for any choice of bounded time-delay and associated augmented matrix
in the form of (<ref>). This completes the proof.

[1]
C. N. Hadjicostis and T. Charalambous,
“Average consensus in the presence of delays in directed graph
IEEE Transactions on Automatic Control, vol. 59, no. 3, pp.
763–768, 2013.

[2]
M. Sharifi and M. J. Yazdanpanah,
“Finite time consensus of nonlinear multi-agent systems in the
presence of communication time delays,”
European Journal of Control, vol. 53, pp. 10–19, 2020.

[3]
Y. Ghaedsharaf, M. Siami, C. Somarakis, and N. Motee,
“Centrality in time-delay consensus networks with structured
Automatica, vol. 125, pp. 109378, 2021.

[4]
S. Sundaram and C. N. Hadjicostis,
“Delayed observers for linear systems with unknown inputs,”
IEEE Transactions on Automatic Control, vol. 52, no. 2, pp.
334–339, 2007.

[5]
L. Schenato,
“Optimal estimation in networked control systems subject to random
delay and packet loss,”
in 45th IEEE Conference on Decision and Control, 2006, pp.

[6]
X. He, X. Ren, H. Sandberg, and K. H. Johansson,
“Secure distributed filtering for unstable dynamics under
compromised observations,”
in 58th IEEE Conference on Decision and Control, 2019, pp.

[7]
S. Battilotti, F. Cacace, and M. d’Angelo,
“A stability with optimality analysis of consensus-based distributed
filters for discrete-time linear systems,”
Automatica, vol. 129, pp. 109589, 2021.

[8]
S. Kar, G. Hug, J. Mohammadi, and J. M. F. Moura,
“Distributed state estimation and energy management in smart grids:
A consensus + innovations approach,”
IEEE Journal of selected topics in signal processing, vol. 8,
no. 6, pp. 1022–1038, 2014.

[9]
F. S. Cattivelli and A. H. Sayed,
“Diffusion strategies for distributed kalman filtering and
IEEE Transactions on Automatic Control, vol. 55, no. 9, pp.
2069–2084, Sep. 2010.

[10]
S. Das and J. M. F. Moura,
“Consensus+ innovations distributed kalman filter with optimized
IEEE Transactions on Signal Processing, vol. 65, no. 2, pp.
467–481, 2017.

[11]
A. Mohammadi and A. Asif,
“Distributed consensus $+ $ innovation particle filtering for
bearing/range tracking with communication constraints,”
IEEE Trans. on Signal Processing, vol. 63, no. 3, pp. 620–635,

[12]
A. Abdelmawgoud, M. Jamshidi, and P. Benavidez,
“Distributed estimation in multimissile cyber-physical systems with
time delay,”
IEEE Systems Journal, vol. 14, no. 1, pp. 1491–1502, 2020.

[13]
A. Jenabzadeh and B. Safarinejadian,
“Distributed estimation and control for nonlinear multi-agent
systems in the presence of input delay or external disturbances,”
ISA transactions, vol. 98, pp. 198–206, 2020.

[14]
L. Zou, Z. Wang, Q. Han, and D. Zhou,
“Moving horizon estimation for networked time-delay systems under
round-robin protocol,”
IEEE Trans. on Automatic Control, vol. 64, no. 12, pp.
5191–5198, 2019.

[15]
L. Liu, A. Yang, X. Tu, M. Fei, and W. Naeem,
“Distributed weighted fusion estimation for uncertain networked
systems with transmission time-delay and cross-correlated noises,”
Neurocomputing, vol. 270, pp. 54–65, 2017.

[16]
M. Doostmohammadian and U. Khan,
“On the genericity properties in distributed estimation: Topology
design and sensor placement,”
IEEE Journal of Selected Topics in Signal Processing, vol. 7,
no. 2, pp. 195–204, 2013.

[17]
H. Silm, D. Efimov, W. Michiels, R. Ushirobira, and J. Richard,
“A simple finite-time distributed observer design for linear
time-invariant systems,”
Systems & Control Letters, vol. 141, pp. 104707, 2020.

[18]
A. Mitra and S. Sundaram,
“Distributed observers for LTI systems,”
IEEE Trans. on Automatic Control, vol. 63, no. 11, pp.
3689–3704, 2018.

[19]
M. Doostmohammadian and U. A. Khan,
“On the distributed estimation of rank-deficient dynamical systems:
A generic approach,”
in IEEE Conf. on Acoustics, Speech, and Signal Processing,
2013, pp. 4618–4622.

[20]
S. Pequito, F. Rego, S. Kar, A. P. Aguiar, A. Pascoal, and C. Jones,
“Optimal design of observable multi-agent networks: A structural
system approach,”
in IEEE European Control Conference, 2014, pp. 1536–1541.

[21]
M. Doostmohammadian, H. R. Rabiee, and U. A. Khan,
“Structural cost-optimal design of sensor networks for distributed
IEEE Signal Processing Letters, vol. 25, no. 6, pp. 793–797,

[22]
M. Doostmohammadian and H. R. Rabiee,
“On the observability and controllability of large-scale iot
networks: Reducing number of unmatched nodes via link addition,”
IEEE Control Systems Letters, vol. 5, no. 5, pp. 1747–1752,

[23]
A. Chapman and M. Mesbahi,
“On strong structural controllability of networked systems: A
constrained matching approach,”
in American Control Conference. IEEE, 2013, pp. 6126–6131.

[24]
M. Doostmohammadian and U. A. Khan,
“On the complexity of minimum-cost networked estimation of
self-damped dynamical systems,”
IEEE Transactions on Network Science and Engineering, vol. 7,
no. 3, pp. 1891–1900, 2019.

[25]
M. Doostmohammadian and U. A. Khan,
“Minimal sufficient conditions for structural
observability/controllability of composite networks via Kronecker
IEEE Transactions on Signal and Information processing over
Networks, vol. 6, pp. 78–87, 2019.

[26]
U. A. Khan and A. Jadbabaie,
“Coordinated networked estimation strategies using structured
systems theory,”
in 49th IEEE Conference on Decision and Control, Orlando, FL,
Dec. 2011, pp. 2112–2117.

[27]
J. Ghaderi and R. Srikant,
“Opinion dynamics in social networks with stubborn agents:
Equilibrium and convergence rate,”
Automatica, vol. 50, no. 12, pp. 3209–3215, 2014.

[28]
S. Pequito, S. Kar, and A. P. Aguiar,
“A framework for structural input/output and control configuration
selection in large-scale systems,”
IEEE Trans. on Automatic Control, vol. 61, no. 2, pp. 303–318,

[29]
Y. Ru and C. N. Hadjicostis,
“Sensor selection for structural observability in discrete event
systems modeled by Petri nets,”
IEEE Transactions on Automatic Control, vol. 55, no. 8, pp.
1751–1764, 2010.

[30]
M. Pirani, J. A. Taylor, and B. Sinopoli,
“Strategic sensor placement on graphs,”
Systems & Control Letters, vol. 148, pp. 104855, 2021.

[31]
J Bu, M Fazel, and M Mesbahi,
“Accelerated consensus with linear rate of convergence,”
in American control conference. IEEE, 2018, pp. 4931–4936.

[32]
C. Commault and J. Van Der Woude,
“A classification of nodes for structural controllability,”
IEEE Transactions on Automatic Control, vol. 64, no. 9, pp.
3877–3882, 2018.

[33]
M. Doostmohammadian and U. A. Khan,
“Measurement partitioning and observational equivalence in state
in IEEE Conference on Acoustics, Speech and Signal Processing,
2016, pp. 4855–4859.

[34]
M. Deghat, V. Ugrinovskii, I. Shames, and C. Langbort,
“Detection and mitigation of biasing attacks on distributed
estimation networks,”
Automatica, vol. 99, pp. 369–381, 2019.

[35]
M. Doostmohammadian and N. Meskin,
“Sensor fault detection and isolation via networked estimation:
Full-rank dynamical systems,”
IEEE Transactions on Control of Network Systems, 2020,
Early access.

[36]
M. Pirani and S. Sundaram,
“On the smallest eigenvalue of grounded laplacian matrices,”
IEEE Trans. on Automatic Control, vol. 61, no. 2, pp. 509–514,

[37]
R. Bhatia,
Matrix analysis,
Springer Science & Business Media, 2013.