# DDPEN: Trajectory Optimisation With Sub Goal Generation Model

Aleksander Gamayunov1,2 Aleksey Postnikov1,2 Gonzalo Ferrer2 1 The authors are
with the Sberbank Robotics Laboratory, Moscow, Russia. {gamayunov.a.r,
<EMAIL_ADDRESS>2Skolkovo Institute of Science and Technology,
Moscow, Russia<EMAIL_ADDRESS>
###### Abstract

Differential dynamic programming (DDP) is a widely used and powerful
trajectory optimization technique, however, due to its internal structure, it
is not exempt from local minima. In this paper, we present Differential
Dynamic Programming with Escape Network (DDPEN) - a novel approach to avoid
DDP local minima by utilising an additional term used in the optimization
criteria pointing towards the direction where robot should move in order to
escape local minima.

In order to produce the aforementioned directions, we propose to utilize a
deep model that takes as an input the map of the environment in the form of a
costmap together with the desired goal position. The Model produces possible
future directions that will lead to the goal, avoiding local minima which is
possible to run in real time conditions. The model is trained on a synthetic
dataset and overall the system is evaluated at the Gazebo simulator.

In this work we show that our proposed method allows avoiding local minima of
trajectory optimization algorithm and successfully execute a trajectory 278 m
long with various convex and nonconvex obstacles.

##  I Introduction

The need for service professions - couriers, cleaners, etc. is growing along
with the population of the planet. At the same time, mobile robotics continue
to develop and can replace at least part of the monotonous work, while
providing people with new highly skilled jobs.

We are witnessing news about novel robot couriers, and robot vacuum cleaners
that have already become commonplace while autonomous tractors are starting to
drive through the agricultural fields.

Navigation is one of the core subsystems for robot’s autonomy which is largely
dependent on path planning and trajectory planning algorithms.

In this work, we consider optimizing algorithms and propose a solution to the
important but common problem of nonconvex trajectory optimization.

Optimization trajectory algorithms iteratively modify the initial assumption
of the trajectory according to the specified criteria, and gradually reduce
the error value. One of the important issues is that the free state space is
usually nonconvex and presents multiple local minima. An optimization
algorithm can reach one of the local minima and stay in it without ever
reaching the global minima.

For a trajectory optimization task, reaching local minima can lead to a
situation where the algorithm will not lead the mobile robot to the goal, even
if the optimal trajectory exists.

Figure 1: Comparison of DDPEN (left side) with DDP (right side) on trajectory
optimization with nonconvex obstacle. The robot is in the center of gray
costmap. Black part of the costmap represents obstacles. Green line is the
optimized trajectory. Green rectangles is predicted sub goals $x_{s}$. Blue
circle is desired goal position $x_{goal}$.

In that case, the mobile robot will not be able to complete its task.

A number of works are devoted to this problem. Zhaoming Xie et al. in his work
of extension the differential dynamic programming trajectory optimization
approach (DDP[1, 2]) with nonlinear constraints (CDDP)[3], mention that “like
all nonconvex optimization, we need a good initial trajectory to ensure that
the algorithm would not get stuck in a bad local minima”.

In that case sampling trajectory planning algorithm like $RRT^{*}$[4] or some
sort of fast neural path planning $A^{*}$[5] algorithms for generating a safe
initial trajectory for optimization can be used. Another approach of local
minima avoidance is based on DDP maximum entropy formulation (MEDDP)[6], with
multimodal policy exploration.

Escaping from local minima is a trivial task for human in case of 2d
environment but the overall process of escaping from the local minimum is ill-
posed with the current gradient-based optimisation tools in DDP.

In this work, we have proposed an approach for escaping from local minima that
provides additional criteria of optimisation, calculated by a neural network,
that “suggests” the best direction for the optimization process.

The Neural network learns to imitate optimal planning algorithm and provides
intermediate goal as additional optimisation criteria that leads the
optimisation process out from local minima.

We have evaluated vanilla DDP optimisation algorithm and our proposed approach
(DDPEN) that takes additional optimization criteria from neural network at
simulated scenarios in Gazebo[7] and have shown that our proposed approach is
being able to complete complex trajectories even with nonconvex obstacles.

##  II Method

###  II-A Trajectory optimisation problem

The state of an agent at time $t$ is represented by $x_{t}$, where $x$ is the
two dimensional vector in the Cartesian coordinate system. Transitions between
states are provided by the dynamic function $x_{t+1}=f(x_{t},u_{t})$ where
$u_{t}$ is a control action at time $t$. Sequence of control action
$U=[u_{t_{0}},u_{t_{1}},...,u_{t_{N-1}}]$ applied to dynamic function forms
the trajectory $X=[x_{t_{0}},x_{t_{1}},...,x_{t_{N}}]$ from initial time
$t_{0}$ to prediction horizon $t_{N}$. The total cost of trajectory $J_{0}$ is
the sum of running costs $l$ and the final cost $l_{f}$. The solution of the
optimal control problem is to find the minimum cost control sequence.

$J_{0}(X,U)=\sum_{i=0}^{N-1}l(x_{i},u_{i})+l_{f}(x_{N})$ (1)
$U^{*}(X)\equiv\textit{arg}\min_{{}_{U}}J_{0}(X,U)$ (2)

The current work focuses on the cost function $l(x,u)$ that forms the
optimisation direction and can move the optimisation process into a local
minimum.

The cost function, in our vanilla realisation of DDP optimisation method,
includes $c_{g}$ \- cost of distance to the goal state $x_{goal}$, $c_{c}$ \-
control cost and $c_{o}$ \- cost representing the distance to obstacles $obs$
at distance less than $d$ meters.

$l(x,u)=c_{c}(u)+c_{g}(x,x_{goal})+c_{o}(x,obs)$ (3)
$l_{f}(x)=c_{g}(x,x_{goal})+c_{o}(x,obs)$ (4)

$c_{c}(u)=\sum u^{2}$ (5) $c_{g}(x,x_{goal})=\left\lVert
x-x_{goal}\right\rVert_{2}$ (6) $c_{o}(x,obs,d)=\sum ReLU(d-\left\lVert
x-obs\right\rVert_{2})$ (7)

###  II-B DDPEN

The proposed method DDPEN takes into account an additional cost $c_{s}$ \-
distance to sub goal $x_{s}$, generated by the neural network Eq 10, which
represents the movement direction of the mobile robot that allows avoiding
local minima.

$l(x,u)=c_{c}(u)+c_{g}(x,x_{goal})+c_{o}(x,obs)+c_{s}(x,x_{s})$ (8)
$c_{s}(x,x_{s})=\left\lVert x-x_{s}\right\rVert_{2}$ (9)

Figure 2: Sub Goal Generation Model architecture.

###  II-C Sub Goal Generation Model

In order to generate future positions that will force the trajectory
optimization algorithm to avoid local minima we propose to utilise deep model
Eq. 10, that is trained to imitate goal proposals generated by $A^{*}$
algorithm, that would satisfy real-time restrictions and whose performance
would not depend on the position of obstacles.

$x_{s}=SG_{\theta}(map,x_{goal}),$ (10)

where $SG_{\theta}$ \- deep Sub Goal Generation Model, map \- local costmap of
environment, $x_{goal}$ \- desired goal position of robot.

Architecture of Sub Goal Generation Model, shown in Figure 2, consist of
convolution based encoder of obstacle map, namely EfficientNet B0[8], linear
layer that embeds goal positions, and transformer based decoder that generates
sub goal positions.

###  II-D Dataset

We have created a synthetic dataset that captures randomly generated obstacle
maps and optimal path to randomly sampled goal on that map. As a base for
obstacle map generation it is taken a costmap with 200x200 cells with the cell
side length of 5cm. Predefined obstacles are created at random positions with
non-blocking center position constraint. Obstacles are positioned through a
random affine transformation for augmentation purposes. After generation of
obstacles, their costs are dilated from occupied cost to free moving space.

We have used $A^{*}$ path planning algorithm to generate examples of paths
from starting position, which is in our case always the grid center, to
randomly sampled goal position near borders of the costmap. After processing
the optimal path from randomly sampled goal positions on a randomly generated
costmap, intermediate positions of sub goals at 30, 50, 70 steps of the
optimal path are saved together with goals and costmaps, forming the overall
dataset. Synthetically generated samples from the dataset is shown in Figure
3.

Figure 3: Synthetically generated samples from dataset. Left: Example of randomly generated costmap Right: Example of paths obtained by A* algorithm to randomly sampled goals. TABLE I: Comparison of DDP and DDPEN time execution of path. The values is the mean$\pm\sigma$ in seconds Type of obstacles | DDP | DDPEN
---|---|---
| Forward pass | Backward pass | Forward pass | Backward pass
Free | $\textbf{191}\pm 1$ | $\textbf{188}\pm 1$ | $\textbf{191}\pm 1$ | $198\pm 1$
Cylinders |  243$\pm$1 | $313\pm 27$ | 255$\pm 8$ | 281$\pm 32$
Cubes without minima | $284\pm 54$ | $\textbf{262}\pm 31$ | $\textbf{243}\pm 15$ | $266\pm 6$
Cubes with minima | fail | fail | $\textbf{245}\pm 7$ | $\textbf{275}\pm 20$

###  II-E Simulator

For conducting experiments and evaluation, we use Gazebo [8] simulator with an
environment map reconstructed from lidar scans and four wheels sidewalk robot
model with a lidar sensor. Environment map and model of robot is shown in
Figure 4. Point cloud data from simulated lidar sensor after processing by
ROS2[9] node provides costmaps that uses as an input to trajectory
optimization algorithm and sub goal generation model. We use a global path
passing through 8 points on the reconstructed map for experiments, as shown in
Figure. 4. $x_{goal}$ is located on the intersection of the local costmap
border and the line between next and previous global path points, and is used
as input to the trajectory optimization algorithm and sub goal generation
model.

Figure 4: Simulated street on the top on figure, simulated sidewalk robot in
environment on the right down and the real robot at the same place on the left
down. Green circles on the simulated street show the path points for
experiments. Figure 5: Robot trajectories in a simulated environment with non-
convex obstacles. The green circle is the beginning of the paths, the red
circle is the end of the paths. The solid and dotted lines are paths created
by DDPEN and DDP respectively. The cross marker is the end of the paths in
case of hitting the local minimum.

###  II-F Running cycle

The implementation of a particular optimization algorithm takes 75ms for the
optimization epoch, where generation of intermediate goal by neural network
takes 9.03 ms ± 139 µs. We use 4 optimisation epochs before executing
calculated trajectory. For the next cycle of 4 optimisations we use previously
generated trajectory with cut past trajectory predictions and added new
assumption of trajectory by applying last control, as shown in Figure 6.
Executor uses generated trajectory for sending controls to the robot at stable
10Hz and approximate controls between predictions while next optimisation
cycle will end.

Figure 6: Running cycle - execution trajectory while optimisation process in
progress

##  III Evaluation

We evaluate vanilla DDP and proposed DDPEN trajectory optimization algorithms
by executing path of 278 meters in both directions 5 times. The path that the
robot must follow passes through the green dots as shown in the Fig. 4.

We manually place different types of obstacles along the path which results in
4 different environments:

  * •

Free \- original obstacle free environment.

  * •

Cylinders \- environment with additionally placed cylinders with the radius 1m
along the path.

  * •

Cubes without minima  \- environment with additionally placed cubes with a
side 2m forming convex obstacles.

  * •

Cubes with minima \- environment with additionally placed cubes with a side 2m
forming nonconvex obstacles.

Resulting trajectories of DDP and DDPEN in an environment with Cubes with
minima is shown in Figure 5.

Both DDP and DDPEN have equal configuration and dynamic model with linear max
velocity limit of 1.5 m/s and angular velocity limit of 0.5 rad/s.

The main goal that we have focused on in this work is the ability to complete
the entire path trajectory even with non-convex obstacles, and we also
evaluated the path execution time as an indicator of system efficiency.

###  III-A Results

Table I shows the average time $\pm$ standard deviation in seconds for passing
4 different scenarios, where longer execution times mean more iterations were
required to converge, failures happens due to weak gradient signals in
configurations that are close a local minima. Table I shows expected result
for DDP which can’t complete trajectory on map with local minima, while our
proposed method DDPEN is being able to complete trajectory with nonconvex
obstacles.

Time of path execution in forward and backward directions are not equal,
because trajectories in different directions have minor differences that
affect the optimization process.

We can see time execution deviation in the same scenarios that shows how
simulation errors affects trajectory optimisation process. A small error in
the the simulation process leads to differences in the execution time of the
trajectory, especially in the case of trajectories located near obstacles.

The sum of trajectory execution mean time have differences on cylinders and
cubes maps - DDP $1102\pm 68$s comparing to DDPEN $1045\pm 37$s, because DDPEN
starts avoiding obstacles based on additional sub goal provided by Sub goal
generation model earlier than obstacle cost $c_{o}$ starts to affect the
vanilla DDP optimization process as shown in Figure 1.

##  IV Conclusion

In this paper we have shown a novel approach DDPEN that allows trajectory
optimization algorithms to escape from local minima. Our approach utilizes an
additional optimization term representing the direction towards the robot
should move in order to escape local minima. This additional optimization term
based at deep sub goal generation model that predicts future directions that
the robot should follow in order to avoid local minima. The model is trained
with synthetic dataset and overall system is evaluated on a simulated scene at
Gazebo simulator.

The evaluation shows that the proposed method DDPEN is more reliable against
the local minima problem, while vanilla DDP fails.

## References

  * [1] D. Mayne, “A second-order gradient method for determining optimal trajectories of non-linear discrete-time systems,” _International Journal of Control_ , vol. 3, no. 1, pp. 85–95, 1966.
  * [2] D. H. Jacobson and D. Q. Mayne, _Differential dynamic programming_. Elsevier Publishing Company, 1970, no. 24.
  * [3] Z. Xie, C. K. Liu, and K. Hauser, “Differential dynamic programming with nonlinear constraints,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2017, pp. 695–702.
  * [4] I.-B. Jeong, S.-J. Lee, and J.-H. Kim, “Rrt*-quick: A motion planning algorithm with faster convergence rate,” in _Robot Intelligence Technology and Applications 3_. Springer, 2015, pp. 67–76.
  * [5] R. Yonetani, T. Taniai, M. Barekatain, M. Nishimura, and A. Kanezaki, “Path planning using neural a* search,” in _International conference on machine learning_. PMLR, 2021, pp. 12 029–12 039.
  * [6] O. So, Z. Wang, and E. A. Theodorou, “Maximum entropy differential dynamic programming,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, 2022, pp. 3422–3428.
  * [7] N. Koenig and A. Howard, “Design and use paradigms for gazebo, an open-source multi-robot simulator,” in _2004 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)(IEEE Cat. No. 04CH37566)_ , vol. 3. IEEE, 2004, pp. 2149–2154.
  * [8] M. Tan and Q. Le, “Efficientnet: Rethinking model scaling for convolutional neural networks,” in _International conference on machine learning_. PMLR, 2019, pp. 6105–6114.
  * [9] S. Macenski, F. Martín, R. White, and J. G. Clavero, “The marathon 2: A navigation system,” in _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2020, pp. 2718–2725.


# Understanding Deep Learning using Topological Dynamical Systems, Index
Theory, and Homology

Bill Basener  Professor of Data Science
University of Virginia, School of Data Science
Charlottesville, VA
<EMAIL_ADDRESS>
###### Abstract

In this paper we investigate Deep Learning Models using topological dynamical
systems, index theory, and computational homology. These mathematical
machinery was invented initially by Henri Poincaré around 1900 and developed
over time to understand shapes and dynamical systems whose structure and
behavior is too complicated to solve for analytically but can be understood
via global relationships. In particular, we show how individual neurons in a
neural network can correspond to simplexes in a simplicial complex manifold
approximation to the decision surface learned by the NN, and how these
simplexes can be used to compute topological invariants from algebraic
topology for the decision manifold with an explicit computation of homology
groups by hand in a simple case. We also show how the gradient of the
probability density function learned by the NN creates a dynamical system,
which can be analyzed by a myriad of topological tools such as Conley Index
Theory, Morse Theory, and Stable Manifolds. We solve analytically for
associated the differential equation for a trained NN with a single hidden
layer of 256 Neurons applied to the MINST digit dataset, and approximately
numerically that it a sink and basin of attraction for each of the 10 classes,
but the sinks and strong attracting manifolds lie in regions not corresponding
to images of actual digits. Index theory implies the existence of saddles.
Level sets of the probability functions are 783-dimensional manifolds which
can only change topology at critical points of the dynamical system, and these
changes in topology can be investigated with Morse Theory.

###### Index Terms:

Machine Learning, Deep Learning, Topology, Index Theory, Conley Index,
Homology, Computational Homology, algebraic Topology, Stable Manifold,
Manifold, Gradient Flow

##  I Introduction

Topology is the study of shape apart from angles and distance. For example, a
triangle is topologically equivalent to a circle, but distinct from a figure-
eight. Combinatorial topology is the study of combinations of simple pieces;
for example collections of line segments connected at endpoints, surfaces
constructed by connecting triangles connected along their edges, or higher
dimensional manifold created be connecting tetrahedra or higher dimensional
analogs called simplexes along their faces. In a similar way, neural networks
are constructed of simple components (neurons), and it is the global
properties between the neurons from the weights that define the profound
learned information. In this paper we provide computation of the set
separating classes learned by a neural network as a set of line segments
connected at their endpoints, and explicitly compute a topological invariant
for this set called the homology groups. The example is simple for the sake of
exposition, but the general principle that weights and architecture of a
neural network can be used to compute topological invariants is very
interesting.

An important use of topology that motivated its development is the ability to
understand behavior in dynamical systems (for example, flows or solutions to
differential equations of a vector field) even when the computation of
analytic solutions unreasonably complicated or impossible. There is a rich
variety of mathematical machinery that has been developed along this generally
calculus-based approach. Morse theory describes the topology of a manifold
based on a gradient function, and index theory gives relationships between the
topology of manifolds and vector fields. For example the ’hairy ball’ theorem
says that any vector field on a 2-dimensional sphere with have locations where
the vector length is zero, and in addition the sum of the indexes for these
locations will equal the Euler characteristic for the sphere, which is 2. The
hairy ball theorem is an example of the Poincaré–Hopf index theorem (the sum
of the indexes of a vector field on a manifold equals the Euler characteristic
of the manifold), and more generally Conley Index Theory and The Fundamental
Theorem of Dynamical Systems [Norton(1995)].

In this paper we consider the vector field on the input space $X$ for a neural
network defined by the gradient of the probability function for a class $G$,
specifically $\nabla P_{G}(x)$ for $x\in X$ and where $P_{G}(x)$ is the
probability for class $G$ assigned to input $x$. One can also study the
gradient field $\nabla P(x)$ where $P(x)=\textrm{max}\\{P_{G}(x)\\}$ where the
max is over all classes $G$; this vectorfield is defined everywhere except the
separating manifolds. We construct a neural network trained on the standard
MINST dataset, explicitly compute the gradient vectorfield from the weights in
the network, and follow the orbits of the MINST data under the flow of this
gradient vector field. It appears there are ten sinks (maxima for the
gradient), and the basic of attraction for each is the classification region
for each class. Interestingly, digit image for every sink has no resemblance
to the digit it is representing, an interesting connection to adversarial
networks. We suggest that, in general, a neural network learns weights, but at
the same time what is being learned is probability density functions on the
(often high-dimensional) input space, and the topology of these pdfs can be
studied by looking at the gradient functions, with a rich mathematical theory
being available for connecting topology with vector fields.

Topological methods have been used to study neural networks recently. For
example, it has been shown that viewing the data passung through successive
layers, the topology of the classes tends to become simpler from each layer to
the next through decreasing Betti numbers, computational homology, and visual
inspection [Naitzat et al.(2020)Naitzat, Zhitnikov, and Lim]. Neural networks
have been used to learn topological signatures or features in images, with
layers to exploit these features for learning tasks [hofer2017deep]. Manifold
geometry and topology has been used to study manifolds that approximate data
clouds [hofer2017deep]. In this paper, our emphasis is on using topology and
dynamical systems theory to understand the classification structures,
separating manifolds, and probability density functions learned by a network
and used for further prediction.

While beyond the scope of this paper, it is worth noting that our two
approaches, combinatorial topology from discrete pieces of a network/manidolf
and vectorfields with calculus, are connected in a deeply mathematical sense.
The homology groups we compute are computed from combinations of simplicies.
There are theories of cohomology groups defined using integrals of vectors
fields and more generally differential forms on manifolds. De Rham’s Theorem
says that under certain definitions for these groups and assumptions, the
discrete-based homology groups are isomorphic to the calculus-based chomology
groups. Hopefully this gives some sense of coherence between our two
topological approaches to study ’what is learned’ by a neural network.

##  II Homology of Separating Manifold: A Simple Example

We constructed the neural netowork with two classes, a single hidden layer of
3 ReLU neurons, and a final classification layer with a single neuron for
classification. The network with trained weights on a simple 2-class problem
is shown in Figure 1.

Figure 1: A simple neural network with trained weights.

It is not hard then to show that the boundary between the classes is the set
$f_{1}^{2}(x_{1},x_{2})=0$ consists of the line segments labeled with their
formulas as shown in Figure 2, where $f_{i}^{1}(x_{1},x_{2})$ is the output of
the ReLU function with the inputs and biases shown in figure 1.

Figure 2: The separating surface where the class probabilities are equal is a
set of line segments which make up the hexegon. The vertices of the hexagon
are labeled $a,b,c,d,e,f$. The output of the neural network is positive inside
the hexagon and zero outside. For each $f_{i}^{1}$, the line where
$f_{i}^{1}=0$ is shown as a bold dashed line, and the shaded side indicates
the side where the output of the ReLU $F_{i}^{1}$ is increasing. The arrows
show the direction of increase in output of the neural network. Inside the
center triangle, $f_{1}=f_{2}=f_{3}=0$ and this the output of the network is
$b_{2}*W^{2}=3\times 3.4=11.4$.

For a definition of homology and background material, see [Hatcher(2000)], but
some definitions are as follows. An $n-$dimensional simplex (or just
$n-$simplex) $(n\geq 0)$ is the convex hull of $n+1$ points (called vertices)
such that no three of the vertices are coplanar; for example a point, a line
segment, a triangle, tetrahedron, etc. A simplicial complex is a set of
simplexes such that if any two simplexes intersect, they do so along one faces
(the lower dimensional sub-simplexes). The separating set
$S=\\{(x_{1},x_{2})|f_{1}^{2}(x_{1},x_{2})=0\\}$ is a simplicial complex
consisting of six 1-dimensional simplexes and six 0-dimensional simplexes
(points). It is a continuous (but not smooth) 1-dimensional manifold.

For a simplicial complex, the group of $n-$chains is

$C_{n}(X)=\\{\sum_{i}m_{i}\sigma_{i}|m_{i}\in\mathbb{Z}\\}$

and each $\sigma_{i}$ is an $n-$dimensional simplex in the complex. The
boundary $\partial$ of an $n-$simplex is the set of $(n-1)-$dimensional faces
(any convex hull of $(n-1)$ of the vertices, which is of course itself an
$(n-1)$-dimensional simplex) with the orientation inherited from the simplex.
Then $\partial_{k}:C_{k}(X)\to C_{k-1}(X)$, and the $k^{\textrm{th}}$ homology
group of $X$ is defined to be

$H_{k}(X)=\text{ker}(\partial_{k})/\textrm{Im}(\partial_{k+1}).$

Clearly, for the separating manifold $S$,

$H_{n}(S)=0,\textrm{ for }n>1$

since there are no $n-$chains with $n>1$. The group $C_{1}$ has the basis
$\\{\overline{ab},bc,cd,de,ef,fa\\}$ and $C_{0}$ has the basis
$\\{a,b,c,d,e,f\\}$. Clearly $\partial_{2}(mmab+mbc+mcd+mde+mef+mfa)=0$ for
any integer $m$ and this defines $\textrm{ker}(\partial_{1})$. Then

$H_{1}(S)=\text{ker}(\partial_{1})/\textrm{Im}(\partial_{2}\cong\mathbb{Z}/0\cong\mathbb{Z}$

and

$H_{1}(S)=C_{0}/\textrm{Im}(\partial_{1})\cong\mathbb{Z}.$

This shows that the separating manifold (curve) has the homology groups of a
circle, and thus is (isomorphic to) a circle. While this example is trivial,
it suggests the tools of computational topology (an excellent text is
[Edelsbrunner & Harer(2010)Edelsbrunner and Harer] while introductory topics
are covered in [Basener(2006)] and also see [Kaczynski & Mrozek(2004)Kaczynski
and Mrozek]) may be useful in determining homology for more complex
situations. At minimum, it shows that the homology groups can be computed from
just knowing the weights and architecture.

##  III Dynamical Systems and Index Theory

A central theorem in dynamical systems and index theory is the Poincaré-Hopf
Index theorem, which says that for a vectorfield $v:M\to\mathbb{R}^{n}$ on an
$n-$dimensional manifold $M$, the sum of the index of the fixed points (aka
zeros, equilibria, or singularities) is equal to the Euler characteristic of
the manifold,

$\sum_{i}I(x_{i})=\chi(M),$

where the index $I$ of a fixed point (or zero) of the vector field is defined
by taking an isolating sphere $S$ around the fixed point and the index is
equal to the degree of the map $S\to S$ defined by

$x\mapsto v(x)/||V(x)||.$

Poincaré proved this in two dimensions, and Hopf extended it to arbitrary
dimensions. Examples are shown in Figure 3, which is in two dimensions where
the index is equal to the number of counterclockwise rotations of the vector
as you travel counterclockwise around a circle isolating the fixed point.

Figure 3: The index for a few different vector fields in 2-dimensions, in
which case this index is equal to the number of counterclockwise rotations of
the vector as you travel counterclockwise around a circle isolating the fixed
point. (Figure is from [Basener(2006)] used with permission.)

There are numerous good sources on the Index Theory for vectorfields and
applications such as [Zhang et al.(2006)Zhang, Mischaikow, and Turk,
Josevich(), Reininghaus & Hotz(2011)Reininghaus and Hotz, Libgober(2012)] and
work on computing and visualizing the topology of vector fields [Scheuermann
et al.(1998)Scheuermann, Kruger, Menzel, and Rockwood, Reininghaus et
al.(2011)Reininghaus, Lowen, and Hotz]. This theory was extended by Cayley and
James Clerk Maxwell to study topography, and generalized to gradients of
functions on manifolds by Marsten Morse in what is now called Morse Thoery
[Milnor(2016)]. This line of theory culminated in Conley Index theory
published by Conley in [Conley(1978)] focusing on isolated invariant sets (not
just fixed points), and also see [Mischaikow & Mrozek(2002)Mischaikow and
Mrozek], and is connected to Floer homology [Salamon(1990)]. One of Conley’s
theorems proving all dynamical systems are composed of (chain recurrent)
invariant sets and gradient-like sets known as the The Fundamental Theorem of
Dynamical Systems [Norton(1995)]. The triangle in the center of Figure 2 is an
isolated invariance set under the flow moving in the direction of increasing
probability. A deep and profound concept throughout these topics, evident even
in the 2-dimensional vector fields, is that discrete objects like the index of
a fixed point or Euler characteristic of a manifold, which are stable under
continuous perturbations, are connected to integrals and differential
constructions.

Figure 4: Computations of the degrees of isolated invariant sets by counting
rotations of the vector field around a isolating curve (topological circle),
where the vectorfield is the gradient of the probability output of the NN for
the class shown in blue. Observe that the index around the outer cueve equals
the sum of the indices of the smaller circles. This is the Poincaré Hopf
Theorem applied to the probability gradient of a NN. (This is a notional
depiction of a neural network for the sake of illustrating the theorem created
by taking the output of a NN and modifying so the probabilities close to zero
around the edges of the region $[-6,6]\times[-6,6]$ and arrows are added
visually to illustrate the theorem.)

We suggest that the concepts of index theory and related topological tools can
be used to study the structure of ’what is learned’ by a neural network.
Specifically, we can apply index theory and general theory for dynamical
systems to the gradient vector field of probability functions on the state
space that are learned by the network. We do this by training a neural network
on the MINST dataset and solving explicitly for the gradient of the
probability functions for the classes. The input images are $28$ by $28$
pixels (so our vector field is in $28\times 28=784$-dimensional space) and we
have a single hidden layer of 256 ReLU neurons and an output layer of 10
softmax classification neurons (one for each of the 0 through 9 digit
classes). The network has a $97\%$ validation accuracy which is reasonable but
relatively well below what can be achieved by deep CNNS which can achieve
$99.77\%$ accuracy or better [Ciregan et al.(2012)Ciregan, Meier, and
Schmidhuber].

Our input space is $\mathbb{R}^{7}84$ (restricted the unti cube with all
values in $[0,1]$), and we denote the weight from input $i$ to neuron $j$ in
the first layer by $W_{ij}^{1}$, with these weights making the $784\times 256$
matrix $W^{1}$. The weight from hidden layer neuron $i$ to softmax neuron $k$
is denoted $W_{jk}^{2}$, and these weights make up the $256\times 10$ matrix
$W^{2}$. The biases are denoted by $b_{j}^{1}$ for the $j$-th neuron in the
first layer and $b_{k}^{2}$ for the $k$-th neuron in the classification layer.
We denote the $j$-th neuron in the hidden layer by $f_{j}$ and the $k$-th
neuron function in the softmax classification layer by $g_{k}$. The gradient
is a simple application of the chain rule, with

$\nabla g_{k}(x_{1},..,x_{7}84)=\left(\frac{\partial g_{k}}{\partial
x_{1}},...,\frac{\partial g_{k}}{\partial x_{784}}\right)$

where

$\frac{\partial g_{k}}{\partial x_{i}}=\sum_{j=1}^{256}\frac{\partial
g_{k}}{\partial f_{j}}\frac{\partial f_{j}}{\partial x_{i}}.$

In Figure 5 we show the projection of the MINST data on the first two PCAs
colored by class. Figure 6 shows this data on the same PCA projection after 9
iterations of Eulers method approximation to the gradient flow
$x^{\prime}=\nabla g_{k}(x)$ stepping by $x\mapsto x+0.05\nabla g_{k}(x)$.
Figure 7 show the data after 99 iterations, projected on the first two PCs
computed from that data.

Figure 5: The MIST data shown with points colored by class membership. This
data is shown projected onto its first two PCs.

Figure 6: The MIST data after 9 iterations of the gradient flow Euler’s
approximation. This data is shown projected onto the first two PCs computed
from the original data.

Figure 7: The MIST data after 99 iterations of the gradient flow Euler’s
approximation. This data is shown projected onto its first to PCs.

In Figure 8, we shown the digit images for 10 iterations each of a 4 (top two
rows) and a 0 (bottom two rows). In each iteration the probability is
increasing. However, observe that the visual look of the image is becoming
less like a digit. This same phenomena is consistent across all digits
observed. Clearly, this is not desirable behaviour as the non-digit-like
images are being classified with very high probability to a single class. It
seems there is a string attracting stable manifold for the probability max for
each class; that is, there is a submanifold in the data space that is far from
the actual data, and under the gradient flow orbits are attracted to this
manifold and approach the max along this manifold. The behavior is depicted
notionally in Figre 9 using the linear differential equation
$x^{\prime}=-4x+6y,y^{\prime}=x-y2$. While our network architecture is known
to be sub-optimal, observed phenomena like this could be used to improve
optimization strategies and cost functions.

Figure 8: The first 10 iterations of a handwritten 4 digit is shown in the top
2 rows, and the first 10 iterations of a handwritten 0 digit is shown in the
bottom 2 rows. Note that with each successive iteration the probability
increases but the image looks less like a number.

Figure 9: A vector field and phase plane for a 2D linear differential equation
that has a strong attracting direction. This is a notional depiction of what
seems to be happening in the probability gradient flow for the NN on the MINST
data, where the circles in the figure notionally represent a distribution of
data that is off the strong attracting manifold, and the 1d strong attracting
direction notionally represents a higher dimensional strong stable manifold.

The function $g_{k}(f_{1},...,f_{256})$ is defined by

$g_{k}(f_{1},...,f_{256})=\frac{e^{\sum_{j}W_{jk}^{2}f_{j}+b_{k}}}{\sum_{k^{\prime}}e^{\sum_{j}W_{jk^{\prime}}^{2}f_{j}+b_{k^{\prime}}}}$

which leads to

$\frac{\partial g_{k}}{\partial
f_{j}}=g_{k}\left[W_{jk}^{2}-\sum_{k^{\prime}\neq
k}W{jk^{\prime}}^{2}g_{k^{\prime}}\right]$

Define $\delta_{j}(x)$ to be the function defined by

$\delta_{j}(x)=\begin{cases}1&\text{if }\sum_{i}W_{ij}^{2}x_{i}+b_{j}^{1}>0\\\
\textrm{undefined}&\text{if }\sum_{i}W_{ij}^{2}x_{i}+b_{j}^{1}=0\\\
0&\text{otherwise}\end{cases}$

(For numerical stability we use $\delta_{j}(x)=0$ when
$\sum_{i}W_{ij}^{2}x_{i}+b_{j}^{1}>0$ in our simulations.) The the partial
derivatives of $f_{j}$ are defined by

$\nabla f_{j}=\left(\frac{\partial f_{j}}{\partial x_{1}},...,\frac{\partial
f_{j}}{\partial x_{256}}\right)=\delta_{j}(x)(W_{1j}^{1},...,W_{256j}^{1})$

So

$\frac{\partial g_{k}}{\partial
x_{i}}=\sum_{j=1}^{256}\delta_{j}(x)g_{k}\left[W_{jk}^{2}-\sum_{k^{\prime}\neq
k}W{jk^{\prime}}^{2}g_{k^{\prime}}\right]W_{jk}^{1}$

More concisely

$\nabla
g_{k}(x)=g_{k}(x)\widetilde{W}^{1}(x)\left[\Phi_{k}(x)+\Gamma(x)\right]$

where $\Phi_{k}$ is the $256\times 1$ column vector whose $j$-th value is
$(1+g_{k}(x))W_{jk}^{2}$; $\Gamma$ is the $256\times 1$ column vector
$W^{2}g(x)$; and $\widetilde{W}^{1}(x)$ is a $784\times 256$ matrix

$\widetilde{W}^{1}_{ij}(x)=\delta_{j}(x)W_{ij}^{1}.$

We now consider the zeros and/or singularities for $\nabla g_{k}(x)$. It is
clear that $g_{k}(x)>0$ for all $x$. The vector $\Phi_{k}$ is nonzero (except
in the degenerate case where $W_{jk}^{2}=0$ for all $j=1,...,256$). The vector
$\Gamma$ is also nonzero. Thus, the attracting sets (possibly simplexes) of
the gradient flow are the points where $\widetilde{W}^{1}_{ij}(x)$ is
undefined or zero.

##  IV Conclusions

In this paper we presented how the homology of a separating manifold can be
computed from the weights in a simple neural network. We suggest that this may
be helpful in understanding the complexity of the separating
surfaces/manifolds and complexity of learned information. Our example used
ReLU activation functions so the results were combinatorial and readily
understood with combinatorial homology, but other activation functions would
lead to smooth gradients and differential tolopogy tools that mirror the
combinatorial ones (cohomology, Jacobians of equilibria points, etc.).

We then showed how dynamical systems and index theory can be applied to give
insight to what is leaned by a neural network. We showed an expository example
in 2 dimensions where the NN and theorem can be understood and computed
visually. We then gave an example in 784 dimensions where numerically
numerical approximations suggest that there are 10 attracting invariant sets.
(Notably there are 10 classes, so the gradient at each point is the gradient
of the class probability for the most probable class.) We showed that the
analytic formula for the gradient can be computed explicitly for this example.
We also showed that this dynamical systems approach uncovers an apparent
weakness in the neural network; although the accuracy is over 97% on
validation data, the attracting sets where each class probiablity is maximized
correspond to images that look nothing line the original data of handwritten
digits.

This work suggests a number of lines for future research. Can topology be used
to better understand the probability functions - for example are the basins of
attraction for the attractors simply connected, in our case for MINST or under
certain circumstances? Examples shown in Figures 3 and 4 show the obvious fact
that basins of attraction are neither always connected nor always comprised of
simply connected components. The differential equations for the gradient flow
are exact and the probability is a Lyapunov function, constraining the type of
dynamics (no periodic orbits or chaos, for example). The most general question
seems to be what topological information can generally be computed for neural
networks, and what are the connections between this topology and the weights
and architecture of the network.

Computing topology is interesting, but is not only an end in itself. It is
likely that some topological / dynamical system / geometric information is
desirable, and thus can be used to either understand current optimization
methods or heuristically observed improvements. For example, rapidly
fluctuating topology durring optimization is an indicator of instability and
need for regularization. In general, perhaps simple topology is preferable for
robustness. Our observation that data flowing under the gradient system is
attract to regions far from the original data as shown in Figure 8 may inspire
new methods. It appears that there the attracting set has a strongly
attracting direction/manifold that is far from the data; data moves quickly to
this manifold and then follows this manifold slowly to the attracting set. It
may be preferable to have this attracting manifold approximate the original
data, and it may be that the strong manifold is connected to the weights - for
example, a submanifold where some of the neuron outputs are zero. Perhaps
these are useful for deriving improved optimization methods or give improved
understanding to currently observed useful methods.

## References

  * [Basener(2006)] Basener, W.  _Toplogy and Its Applications_.  John Wiley & Sons, Inc., 2006.  doi: 10.1002/9780470067949.
  * [Ciregan et al.(2012)Ciregan, Meier, and Schmidhuber] Ciregan, D., Meier, U., and Schmidhuber, J.  Multi-column deep neural networks for image classification.  In _2012 IEEE conference on computer vision and pattern recognition_ , pp. 3642–3649. IEEE, 2012.
  * [Conley(1978)] Conley, C. C.  _Isolated invariant sets and the Morse index_.  Number 38. American Mathematical Soc., 1978.
  * [Edelsbrunner & Harer(2010)Edelsbrunner and Harer] Edelsbrunner, H. and Harer, J.  _Computational Topology - an Introduction._ American Mathematical Society, 2010.  ISBN 978-0-8218-4925-5.
  * [Hatcher(2000)] Hatcher, A.  _Algebraic topology_.  Cambridge Univ. Press, Cambridge, 2000.  URL https://pi.math.cornell.edu/ hatcher/AT/ATpage.html.
  * [Josevich()] Josevich, P.  The poincaré-hopf index theorem and the fundamental theorem of algebra.  _Brown University, articolo disponibile in rete_.  URL https://www2.math.upenn.edu/ pjmcgrat/research/pumj.pdf.
  * [Kaczynski & Mrozek(2004)Kaczynski and Mrozek] Kaczynski, T., M. K. and Mrozek, M.  _Computational Homology_.  Springer-Verlag, New York, 2004.
  * [Libgober(2012)] Libgober, J.  The euler characteristic, poincaré-hopf theorem and applications.  _Uchicago math REU_ , 17(6), 2012.  URL http://www.math.uchicago.edu/ may/VIGRE/VIGRE2010/REUPapers/Libgober.pdf.
  * [Milnor(2016)] Milnor, J.  _Morse Theory.(AM-51), Volume 51_.  Princeton university press, 2016.
  * [Mischaikow & Mrozek(2002)Mischaikow and Mrozek] Mischaikow, K. and Mrozek, M.  Conley index.  _Handbook of dynamical systems_ , 2:393–460, 2002.
  * [Naitzat et al.(2020)Naitzat, Zhitnikov, and Lim] Naitzat, G., Zhitnikov, A., and Lim, L.-H.  Topology of deep neural networks.  _J. Mach. Learn. Res._ , 21(184):1–40, 2020.
  * [Norton(1995)] Norton, D. E.  The fundamental theorem of dynamical systems.  _Commentationes Mathematicae Universitatis Carolinae_ , 36(3):585–597, 1995.
  * [Reininghaus & Hotz(2011)Reininghaus and Hotz] Reininghaus, J. and Hotz, I.  Combinatorial 2d vector field topology extraction and simplification.  In _Topological Methods in Data Analysis and Visualization_ , pp. 103–114. Springer, 2011.
  * [Reininghaus et al.(2011)Reininghaus, Lowen, and Hotz] Reininghaus, J., Lowen, C., and Hotz, I.  Fast combinatorial vector field topology.  _IEEE Transactions on Visualization and Computer Graphics_ , 17(10):1433–1443, 2011.  doi: 10.1109/TVCG.2010.235.
  * [Salamon(1990)] Salamon, D.  Morse theory, the conley index and floer homology.  _Bulletin of the London Mathematical Society_ , 22(2):113–140, 1990.
  * [Scheuermann et al.(1998)Scheuermann, Kruger, Menzel, and Rockwood] Scheuermann, G., Kruger, H., Menzel, M., and Rockwood, A.  Visualizing nonlinear vector field topology.  _IEEE Transactions on Visualization and Computer Graphics_ , 4(2):109–116, 1998.  doi: 10.1109/2945.694953.
  * [Zhang et al.(2006)Zhang, Mischaikow, and Turk] Zhang, E., Mischaikow, K., and Turk, G.  Vector field design on surfaces.  _ACM Transactions on Graphics (ToG)_ , 25(4):1294–1326, 2006.


Send correspondence to:
N.M<EMAIL_ADDRESS>// W.G<EMAIL_ADDRESS>
# Intensity Interferometry at Calern and beyond:
progress report

Nolan Matthews  Université Côte d’Azur, CNRS, Institut de Physique de France,
France  Jean-Pierre Rivet  Université Côte d’Azur, Observatoire de la Côte
d’Azur, CNRS, Laboratoire Lagrange, France  Mathilde Hugbart  Université Côte
d’Azur, CNRS, Institut de Physique de France, France  Guillaume Labeyrie
Université Côte d’Azur, CNRS, Institut de Physique de France, France  Robin
Kaiser  Université Côte d’Azur, CNRS, Institut de Physique de France, France
Olivier Lai  Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS,
Laboratoire Lagrange, France  Farrokh Vakili  Université Côte d’Azur,
Observatoire de la Côte d’Azur, CNRS, Laboratoire Lagrange, France  David
Vernet  Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, UMS
Galilée, France  Julien Chabé  Université Côte d’Azur, Observatoire de la Côte
d’Azur, CNRS, Laboratoire Géoazur, France  Clément Courde  Université Côte
d’Azur, Observatoire de la Côte d’Azur, CNRS, Laboratoire Géoazur, France
Nicolas Schuhler  European Southern Observatory, Vitacura, Chile  Pierre
Bourget  European Southern Observatory, Vitacura, Chile  William Guerin
Université Côte d’Azur, CNRS, Institut de Physique de France, France

###### Abstract

We present the current status of the I2C stellar intensity interferometer used
towards high angular resolution observations of stars in visible wavelengths.
In these proceedings, we present recent technical improvements to the
instrument, and share results from ongoing campaigns using arrays of small
diameter optical telescopes. A tip-tilt adaptive optics unit was integrated
into the optical system to stabilize light injection into an optical fiber.
The setup was successfully tested with several facilities on the Calern
Plateau site of the Observatoire de la Côte d’Azur. These include one of the 1
m diameter telescopes of the C2PU observatory, a portable 1 m diameter
telescope, and also the 1.5 m MéO telescope. To better constrain on-sky
measurements, the spectral transmission of instrument was characterized in the
laboratory using a high resolution spectrograph. The system was also tested
with two of the auxiliary telescopes of the VLTI resulting in successful
temporal and spatial correlation measurements of three stars.

###### keywords:

Intensity Interferometry, Optical Interferometry, I2C

##  1 INTRODUCTION

Intensity interferometry was originally pioneered by Robert Hanbury Brown and
Richard Q. Twiss (HBT) in the 1950s, leading to the successful implementation
and operations of the Narrabri Stellar Intensity Interferometer that measured
the angular diameters of 32 stars [1]. The technique was subsequently
abandoned in astronomy in favor of the far more sensitive techniques that had
emerged in stellar amplitude interferometry with separated telescopes[2] of
which has steadily evolved into the current exceptional interferometric
facilities present worldwide such as the CHARA111https://www.chara.gsu.edu/
array and the
VLTI222https://www.eso.org/sci/facilities/paranal/telescopes/vlti.html.
However, over the last decade and a half there has been strong experimental
efforts in astronomical HBT interferometry, that can be broadly divided into
two complementary approaches. One is to use imaging atmospheric Cherenkov
telescopes (IACTs), principally designed for ground based gamma-ray astronomy
for intensity interferometry by taking advantage of their large diameters ($>$
10 m), and array based layouts[3]. Successful on-sky measurements have been
demonstrated with both the MAGIC[4] and VERITAS[5] observatories with key
implications for a possible implementation on the future Cherenkov Telescope
Array that holds significant promise if operated additionally as an HBT
observatory [6]. An alternative approach to implementing HBT interferometry is
to use optical telescopes whose significantly better optical quality allows
the usage of advanced technologies such as state of the art single photon
counting detectors. Our I2C consortium333 https://inphyni.univ-
cotedazur.eu/sites/cold-atoms/research/i2c, primarily consisting of a
collaboration of physicists and astronomers in the network of the Université
Côte d’Azur, is focused on this latter approach, demonstrating leadership
through observations with the facilities on the Calern Plateau site of the
Observatoire de la Côte d’Azur.

Our group has made steady progress in the development of an on-sky HBT
instrument demonstrating both technical and scientific achievements. Single
telescope temporal intensity interferometry measurements were first
successfully performed with a 1 m diameter telescope[7], and soon thereafter
extended to spatial intensity interferometry, measuring the HBT correlations
between a pair of 1 m telescopes with a 15 m baseline at the C2PU
observatory[8]. Subsequently, the angular diameter of the bright luminous blue
variable (LBV) star P Cyg was measured in its H-alpha emission, which provided
a distance measurement when combined with a linear diameter obtained from
spectroscopic modeling[9]. Using an improved version of our instrument, this
technique was reapplied to P Cyg and additionally to $\beta$ Ori.[10] The long
term goal of our group is to continue the progress made thus far towards an
instrument which can perform interferometric measurements in visible
wavelengths (B, V, R, I photometric bands) with sensitivity to stars brighter
than at least magnitude 8 if installed onto major facilities.

The focus of this proceeding is to provide an update on recent technical
improvements made to extend the I2C instrument to additional optical
telescopes and to increase the sensitivity to fainter stars, or equivalently
with improved precision. A first order tip-tilt adaptive optics unit was
integrated into the optical system to stabilize and improve the overall light
collection efficiency of the system. The setup was sequentially tested on
several of the telescopes available on the Calern Plateau and allows for
coordinated campaigns between the different telescopes, potentially providing
up to 4 telescope intensity interferometry observations with semi-configurable
baselines ranging approximately from 15 to 150 meters at visible wavelengths.
Furthermore, the temporal and spectral characteristics of the instrument were
extensively measured in the laboratory. Laboratory HBT correlation peaks were
measured at the level of 1%, thus demonstrating the potential of squared
visibility measurements with a precision of 10-2. Finally, the system was
successfully tested at the VLTI using two of the auxiliary telescopes in
observations of three stars.

###  1.1 Brief Theory of Intensity Interferometry

An intensity interferometer performs the correlation of starlight intensity
fluctuations between two telescopes in order to measure the squared
visibility. For two telescopes with a projected baseline $r$ between them, the
second order coherence function is

$g^{(2)}(r,\tau)=\frac{\langle I_{1}(t)I_{2}(r,t+\tau)\rangle}{\langle
I_{1}\rangle\langle I_{2}\rangle}$ (1)

where $I_{1}$ and $I_{2}$ are time-dependent intensities recorded at the two
telescopes with a time-lag $\tau$ between them, and the brackets indicate
averaging over time $t$. The second-order coherence function “ $g^{(2)}$ ” is
related to the first order coherence function $g^{(1)}$ through the Siegert
relation[11, 12] by

$g^{(2)}(r,\tau)=1+|g^{(1)}(r,\tau)|^{2}$ (2)

where the first-order coherence function can be separated into spatial and
temporal components

$g^{(1)}(r,\tau)=V(r)g^{(1)}(\tau)$ (3)

where $V(r)$ is the interferometric visibility of the source, given by the
Fourier transform of the source sky brightness distribution. In the case of an
unresolved source $V(r)=1$, and the measured second order coherence function
will depend only on the temporal component $g^{(1)}(\tau)$ given by the
Fourier transform of the measured light spectral density by the Wiener-
Khintchine theorem. For linearly polarized thermal light at zero optical path
delay $g^{(1)}(\tau=0)=1$ whereas at large time-lags much greater than the
coherence or correlation time $g^{(1)}=0$ such that there is a “bunching peak”
centered about zero optical path delay with an effective temporal width given
by the coherence time. The amplitude of this peak at zero time-lag measures
the squared visibility under the condition that the instrumental resolving
time is shorter than the coherence time. The coherence time as defined by
Mandel & Wolf (1995)[13] can be described by the integral of the squared
first-order coherence function

$T_{c}=\int|g^{(1)}(\tau)|^{2}d\tau=\int|s(\nu)|^{2}d\nu$ (4)

which by Parseval’s theorem is equal to the integral of the squared normalized
spectral density $s(\nu)$. For visible light with a bandpass of
$\Delta\lambda\sim 1\,$nm the corresponding coherence time is of order 1 ps,
much shorter than what can be achieved with conventional detectors. In this
case, a measurement averages over many coherence times and reduces the value
of the $g^{(2)}$ peak amplitude at $\tau=0$ by a factor of $\sim T_{c}/T_{d}$
where $T_{d}$ is the effective time-resolution of the detector. The result is
that amplitude of the $g^{(2)}$ peak does not directly measure the squared
visibility, but instead the squared visibility times this loss of contrast.
The squared visibility can be extracted by dividing the value of the
$g^{(2)}(r)-1$ peak measured between telescopes to the $g^{(2)}-1$ peak
measured at zero-baseline

$|V(r)|^{2}=\frac{g^{(2)}(r,\tau=0)-1}{g^{(2)}(r=0,\tau=0)-1}$ (5)

under the assumption that the profile of the $g^{(2)}(\tau)$ peak is constant.
In practice, we measure the ratio of the area of the $g^{(2)}(r)-1$ peak, to
the area of the $g^{(2)}(r=0)-1$ peak for the squared visibility, where at
zero-baseline this area should correspond to the coherence time that can be
calculated from the detected spectrum as given by Equation 4.

---
Figure 1:  Photographs of facilities available for intensity interferometry
observations on the Calern Plateau site of the Observatoire de la Côte d’Azur.
These include the MéO optical metrology station (top-left), the T1M portable
telescope (top-right), and the C2PU observatory (bottom right). A satellite
photograph (Google) of the facilities are shown in the bottom-left with C2PU
in the Southeast, and MéO at the Northwest. The yellow crosses indicate the
locations at which differential GPS positions were obtained, and serve at the
current possible positions of the T1M telescope.

##  2 Intensity Interferometry on the Calern Plateau

Observations with the I2C instrument are primarily conducted using the
facilities of the Observatoire de la Côte d’Azur available on the Calern
Plateau site located in southern France. Previously reported observations were
performed with the Centre Pédagogique Planète Univers (C2PU) facility
consisting of a pair of effectively identical 1 m diameter optical telescopes
that are separated by a fixed 15 m separation on a nearly East-West
baseline444https://www.oca.eu/fr/c2pu-accueil. Recently, technical intensity
interferometry tests were successfully performed with the 1.5 m diameter MéO
(Métrologie Optique) telescope that is located approximately 150 m to the
North-West of the Omicron (westernmost) telescope of C2PU. Similar tests were
performed with the T1M telescope, a portable 1 m diameter Newtonian telescope
on a Dobson-type fully motorized azimuthal mount that enables configurable
baselines to extend the accessible sampling of the uv-plane. These facilities
are shown in Figure 1. To help aid in the planning of future observations we
have been able to incorporate a custom array configuration file within the
interferometric observational planning tool ASPRO[14] to coordinate the
optimal baseline configurations for a given target.

###  2.1 Experimental setup

Figure 2:  The left figure sketches the general operation of the coupling
assemblies (CAs). The tip-tilt correction is applied to the incident beam to
provide corrections in the XY plane. The diverging lens ($L_{1}$) produces a
collimated beam required to properly filter the light with the narrowband
filter ($NF$), typically of width $\Delta\lambda\sim 1\,$nm. A dichroic
beamsplitter ($D$) reflects blue light towards a guiding camera used for the
tip-tilt control. The transmitted red light is passed through the narrowband
filter, and then is separated into two orthogonal linear polarization states
by a polarizing beam splitter ($PBS$). On the reflected beam, an additional
linear polarizer is included ($P$) to further reject the complementary
polarization state. The light is then coupled to a graded-index multi-mode
fiber (GRIN-MMF) with a 100$\,\mu$m core diameter using a converging lens
($L_{2}$). The right image displays a photograph of the CA for the T1M
telescope that adapts to the output port of the telescope. Note that the
incident beam does not always come directly from the telescope: The light is
first modified by either a Barlow system (shown by the leftmost cylindrical
element in the right photograph) or focusing lens, as described in Section
2.1.

To enable each telescope for intensity interferometry observations we install
a coupling assembly (CA) at the focus of each telescope. These CAs differ
slightly for each telescope to account for differences in the primary focal
length. The general structure of the coupling assembly is pictured in Figure
2. The main purpose of the CA is to perform spectral and polarization
filtering, and subsequent injection of the filtered starlight into a multimode
optical fiber. Each CA separates the light into two orthogonal polarization
states which are each injected into an APC optical fiber to allow for either
polarization resolved measurements or improve sensitivity by co-adding the
$g^{(2)}$ measurements made in each polarization mode. With exception to the
C2PU telescopes, the input converging beam shown in Figure 2 does not come
directly from the telescopes. On the T1M telescope ($F=3\,$m), a corrected
Barlow device (Powermate x2) is inserted upstream from the CA, to enhance the
focal length. On the MéO telescope ($F=31\,$m), the focal length is first
reduced by a converging lens placed upstream from the CA. Furthermore, the
focal length of the collimating lens (optic $L_{1}$ in Figure 2) is different
on both telescopes ($f=-50$ mm for the T1M and C2PU, and $f=-75\,$mm on MéO).

The fiber for both polarization states are fed to 50/50 fiber splitters that
allows a zero-baseline $g^{(2)}$ measurement to be compared against
expectation from the spectrum. The outputs of the fiber splitters are
connected to single photon avalanche diode (SPAD) detectors (Excelitas SPCM-
ARQH) with a characteristic timing resolution of $\sim 500\,$ps. The detection
of a photon generates an electronic pulse that is then fed to a time to
digital converter (Swabian Instruments) that records the pulse arrival time
and performs the cross-correlation between all pairs of detectors measuring
parallel polarization states. The $g^{(2)}$ correlations are recorded to disk
typically every 10 seconds over a time-lag window of $\pm 300\,$ns with 10 ps
binning. A software analysis then shifts each correlation function in time by
the instrumental and geometrical optical path delays such that they are
centered about zero optical path delay. The time-corrected correlations are
then averaged together in order to resolve the $g^{(2)}$ bunching peak.

###  2.2 Inclusion of tip-tilt adaptive optics

In our setup, the size of the multimode fiber core corresponds to an on-sky
angle of $\sim 7\,"$. Any imperfection in the telescope tracking results in a
loss of collected flux, degrading the overall sensitivity. The tracking
precision of each telescope varies, but is especially pronounced in the case
of the T1M due to a non-permanent mounting that allows for portability but
makes it especially subject to environmental conditions, particularly wind. In
past observations, any small imperfections in the tracking were corrected
using manual commands sent to the telescope that were typically required every
few minutes.

To remove the need for manual corrections, and to improve the overall light
collection efficiency, a first-order tip-tilt adaptive optics unit (Starlight
Xpress AO-USB) was integrated into the CA and can be seen in Figure 2. The
tip-tilt optic consists of a 60 mm diameter glass plate with a thickness of 13
mm that can be tilted by small stepper motors in order to displace the optical
beam in the XY plane. The range of displacements in the XY plane is
$\pm\,0.15$ mm or equivalently 1.5 times larger than the MMF core diameter. If
a larger correction is needed, a command can be sent from the ST4 port output
of the tip-tilt module to adjust the telescope position in order to bring the
image of the star back within the range of accessible displacements provided
by the tip-tilt.

The tip-tilt is controlled via a program written in Python that applies
commands based on the position of the star image in the CA guiding camera
relative to a desired setpoint. Both the guiding camera and the tip-tilt are
connected via USB to a laptop that runs the tip-tilt control software. The
tip-tilt can receive commands up to 5 ms per increment limiting the response
to frequencies less than 200 Hz. In practice, we are able to achieve a rate
correction of $\sim 40\,$Hz which takes into account the image processing
time. This rate is more than sufficient to correct for the slow imperfections
in the telescope tracking, but cannot account for faster atmospheric effects.

Initial tests of the tip-tilt were performed in late Summer 2021 / Fall 2021
with several of the telescopes available at Calern. To quantify the
improvement in light collection efficiency we recorded data, taken using the
Omicron telescope of the C2PU, where the tip-tilt control was alternated on
and off approximately every 5 minutes for a duration of one hour. The main
results of these tests are shown in Figure 3. We find a modest improvement of
8% in the overall efficiency when comparing the mean count rate during the on
cycles to the off cycles. In addition, the average displacement of the star
image centroid on the guiding camera was measured over a duration of
approximately 1 minute. The tip-tilt reduced the root mean square radial
fluctuations about the mean from 15.0 $\mu$m to 8.4 $\mu$m, or equivalently
from 15 % to 8.4 % of the fiber core diameter. However, we note that the
amount of improvement in fiber injection efficiency is strongly dependent on
additional factors, and in particular wind and seeing conditions which affect
each telescope differently to varying degrees. Nevertheless, the inclusion of
the tip-tilt provided stable light injection into the fiber for periods of
several hours without any manual correction, thus achieving the primary goal.

During January 2022 an observational campaign was performed on $\gamma\,$Cas
in it’s H$\alpha$ emission line using the MéO and T1M telescopes. The source
was observed in multiple baseline configurations over a duration of
approximately 78 hours of observation. A clear $g^{(2)}$ temporal bunching
peak was observed with both telescopes. Here, we show the result for the T1M
telescope in Figure 4. A full analysis of the data, including reduction of the
spatial correlations between the two telescopes with subsequent modeling, is
still pending and will be reported in the future. However, the zero-baseline
$g^{(2)}$ peak reported here demonstrates the usefulness of the tip-tilt optic
to enable intensity interferometry observations onto additional telescopes.

Figure 3:  The left figure shows the detected count rate in units of photons
per second as a function of time. The regions shaded in blue/white indicate
when the tip-tilt control was active/deactive, respectively. The right figure
shows the displacement of the star centroid about the mean along the X and Y
axes over a duration of 1 minute for when the tip-tilt control was applied
(orange) and deactivated (blue). Figure 4: Temporal $g^{(2)}(\tau)$
measurement of Gam Cas observed by the portable T1M telescope. The data (solid
blue line) are fit by a Gaussian (dotted black line).

##  3 Laboratory Developments

In parallel to the observational campaigns, there have been several
developments made in the laboratory to provide a robust characterization of
the instrument and examine any limiting sensitivity factors. In particular, a
high resolution spectrograph was introduced into our laboratory setup to
experimentally characterize the spectral transmission of the narrowband
interferometric filters used in our CA. In addition, we perform high precision
measurements of the $g^{(2)}$ bunching peak at the level of 1-2%.

###  3.1 Spectral characterization of narrow-band filters

Squared visibilities are extracted by the ratio of the area of the $g^{(2)}$
peak measured between telescopes to that of a zero-baseline correlation. The
area of the $g^{(2)}$ peak at zero-baseline is in principle given by the
coherence time of the light that is dependent on the detected optical
spectrum. To validate the on-sky zero baseline measurement it is therefore
useful to properly characterize the spectral transmission of the system. In
general, this necessarily includes all instrumental and physical contributions
such as the source spectrum and instrumental throughput. The instrumental
throughput is primarily given by the transmission of our narrowband
interferometric filters, typically of width $\Delta\lambda\sim 1$ nm. Other
instrumental factors such as the mirror reflectivity, detector response, and
fiber transmission display weak wavelength dependence over the transmission
band of the filter and are assumed as constant.

To characterize the filter transmission, we constructed a high resolution
spectrograph by modifying a monochromator (THR 1500) in which the output beam,
normally focused onto the output slit, is instead focused onto a camera. The
spectrograph was calibrated using a wavelength tunable laser diode in
conjunction with an optical spectrum analyzer. The attained spectral
resolution was less than 0.01 nm at 780 nm ($R\sim 80000$) which is more than
sufficient to characterize our filters. To test the spectral transmission of
the filters, we first inject thermal light from a halogen lamp into a
multimode fiber. The output beam of the fiber is collimated and directed into
our CA. The CA performs the spectral and polarization filtering, and then re-
injects the light into a different multi-mode fiber. This fiber is fed to an
optics bench at the input of the spectrograph, which conjugates the image of
the output fiber and refocuses it onto the input slit of the spectrograph.

To measure the spectra, we take a series of three images consisting of raw
spectra, a flat-field correction, and a dark frame. The flat-field image is
performed by recording an image with the narrowband filter removed from the
beam path, and the dark image is formed by recording an image with the input
beam blocked. The corrected spectrum is then obtained by

$T(\lambda)=\frac{\sum_{i}R(i,\lambda)-\sum_{i}D(i,\lambda)}{\sum_{i}F(i,\lambda)-\sum_{i}D(i,\lambda)}$
(6)

where $R,D,F$ are the raw, dark, and flat-field images with spatial coordinate
$i$ and wavelength coordinate $\lambda$. An example of the spectral
transmission measured in the laboratory are shown in Figure 5. The data
correspond to the measurement of two spectral filtering setups, each one used
in separate coupling assemblies. Each setup consists of two sequential
narrowband filters where both are specified with a central wavelength of 780
nm, but with different bandpasses of 1 nm and 10 nm. The wider bandpass filter
is included in order to block out of band transmission, which is pronounced
for the 1 nm bandpass filter. We find that the two filtering setups have a
comparable peak transmission of $\sim 82\,$%. This value is compatible with
the vendor specified transmission of $>90\,$% for a single filter, such that
for both filters in a given setup the total transmission should be $>81\,$%.
In Figure 5 the theoretical spectra shown assumes 95 % transmission in both
filters. In addition, the measured spectral transmission shows that the
central transmitted wavelength are shifted with respect to the theoretical
expectation, and furthermore to one another by approximately 0.1 nm. We note
these measurements are consistent with similar shifts seen in lower spectral
resolution measurements from the vendor. If these non overlapping transmission
bands are unaccounted for, this will lead to a underestimation of the squared
visibility. Similar tests were performed on other filters, including the
H$\alpha$ filters used in on-sky observations[9, 10], and were found to be in
close agreement with each other.

Figure 5:  Measured spectral transmission of two $\lambda_{0}\sim$ 780 nm
filters (in blue and orange) compared with the theoretical spectra (dotted
black line).

###  3.2 Precision measurements of HBT peaks

To fully characterize the experimental setup we perform measurements of
artificial stars in the laboratory. The artificial star simulates a point-like
thermal source and thus can be used to compare the measured $g^{(2)}$ bunching
area to the coherence time expected from the spectral transmission, and
against on-sky observations of the $g^{(2)}$ zero-baseline measurement of
stellar targets. Furthermore, it allows testing for any systematic effects
such as spurious correlations, or detector variability, that may be difficult
to characterize during on-sky observations. The artificial star is generated
by first injecting light from a halogen lamp into a single mode fiber. The
output beam of the fiber is collimated and then directed into our CA in order
to measure the $g^{(2)}$ function. We record the $g^{(2)}$ correlation over a
wide range of time lags from $\pm$300 ns with 10 ps resolution. Correlations
are typically recorded every minute over duration a few to a hundred hours.
Once the acquisition is complete, we analyze the data by summing all
correlations together, and then normalizing in order to obtain $g^{(2)}$. The
reduced data is fitted to a Gaussian in order to extract a bunching peak area
to measure the coherence time.

Figure 6:  The left figure shows the measured $g^{(2)}$ bunching peak in the
blue points as a function of the time-lag, binned into 100 ps bins. A data is
fit by a Gaussian which is shown by the dotted orange line. The right figure
shows the measured uncertainty in $g^{(2)}$ as a function of the number of
coincidences shown in blue where the dotted black line is the expected
uncertainty due to shot-noise fluctuations. The uncertainty is measured by
taking the root mean square fluctuations of $g^{(2)}$ over time-lags in which
the bunching peak is negligible.

Figure 6 shows the measured $g^{(2)}$ bunching peak for one of the
$\lambda_{0}=780\,$nm filters described in Section 3.1. The data was taken
continuously over a duration of 73.9 hours. Figure 6 also shows the measured
uncertainty as function of the total number of coincidences. This uncertainty
is obtained by taking the root mean square fluctuations of the accumulated
$g^{(2)}$ value over time-lags in which the bunching peak is negligible, taken
here as all time-lags between $3<|\tau|<20$ ns, where $\tau$ is the time-lag
difference from the $g^{(2)}$ peak center. The accumulated $g^{(2)}$ peak is
fitted by a Gaussian function in which we find a standard deviation of width
$392\pm 6\,$ps, and an amplitude of $\left(1.504\pm 0.017\right)\times
10^{-3}$ corresponding to a bunching area or coherence time of $1.477\pm
0.019$ ps. This value is in very good agreement with the expectation from the
measured spectra of 1.474 ps but notably lower than what is expected from the
theoretical spectra of 1.69 ps.

The uncertainty in the measured bunching area can be directly corresponded to
an uncertainty in squared visibility. Here, we demonstrate in the laboratory a
measurement of the bunching area at a level of 1.3 %, suggesting that a
similar uncertainty can be obtained on stellar sources, provided similar
number of photon statistics. While such precision with 1-meter class telescope
and single wavelength channel observations under reasonable observations times
(e.g. $<100\,$hours) are so-far limited to bright stars, the measurements are
nevertheless robust. The data indicates we are free from systematic effects in
our correlation to the level of our measurement precision since the
uncertainty is at the theoretical shot-noise limit.

##  4 Technical tests at the VLTI

Figure 7:  The left photograph displays a CAD drawing of the I2C-VLTI
interface + coupling module. The incident beam arrives from the top of the
image where it is reflected off a dichroic onto a modified coupling assembly.
The coupling assembly performs spectral and polarization filtering of the
light, and focuses the light into a multimode fiber that is connected to the
detectors. The right photograph shows a picture of the entire interface as
mounted in the ROS of one of the ATs.

During March 2022 we were given time to install and perform technical
observations with our I2C instrument using two of the 1.8 m diameter auxilary
telescopes (ATs) of the Very Large Telescope Interferometer while they were
located in their maintenance stations. The main purpose of this mission was to
establish the technical feasibility of performing intensity interferometry
observations with the ATs. These tests occurred during time in which the ATs
were not scheduled for other science programs due to the delay lines being in
use with the unit telescopes (UTs). To adapt our setup to the ATs, a custom
mechanical interface was designed to direct starlight into a modified version
of our fiber injection and spectral filtering coupling assemblies. For these
observations the coupling assemblies were set up with H $\alpha$ filters with
a bandpass of $\Delta\lambda\sim 1$ nm. The fiber injection setup for both
telescopes were installed inside the Relay Optics Structure (ROS) of the ATs.
The interface is easily removable such that the optical configuration of the
ROS can be quickly reset to its original state. The interface consists
primarily of a dichroic that is mounted in the beam path that redirects light
onto our coupling assembly. A CAD drawing and photograph of the interface as
installed in the ROS are pictured in Figure 7.

Figure 8:  Temporal and spatial $g^{(2)}$ measurements of three different
stars using two of the VLTI-AT telescopes. The top-left, top-right, and
bottom-left plots show the zero-baseline correlations measured for Sirius,
Rigel, and Spica, each fit with a Gaussian function shown in the dotted black
line. The bottom-right plot shows the averaged spatial correlation between the
ATs of Spica that reveals a loss of contrast in the $g^{(2)}$ peak associated
with a reduction in the visibility.

The entire system was installed during a single afternoon using an artificial
light beacon to test and align the modules. The first night was used to ensure
proper on-sky alignment and to assess the count rate on several sources. The
following 5 nights were devoted to performing HBT measurements of stellar
targets. In total we obtained 2.3 hours on Alf CMa (Sirius), 12.9 hours on Bet
Ori (Rigel), and 34.3 hours on Alf Vir (Spica). For each source, we were able
to resolve the zero-baseline $g^{(2)}(\tau)$ bunching peaks shown in Figure 8.
To evaluate the bunching area we performed a Gaussian fit to the averaged
correlation function obtained for each star. The bunching area measured was
then compared against the expected value calculated from the spectrum. The
spectrum was formed by combining the narrow-band filter transmission and
target spectra that was contemporaneously measured with the interferometric
observations. For the bunching areas we find for Sirius $\tau_{m}=1.11\pm
0.09\,$ps ($\tau_{s}=1.14\,$ps), for Rigel $\tau_{m}=1.07\pm 0.07$ ps
($\tau_{s}=1.17\,$ps), and for Spica $\tau_{m}=1.22\pm 0.09$ ps
($\tau_{s}=1.12\,$ps), where $\tau_{m}$ is the measured area, and $\tau_{s}$
is the expected value from the spectrum. All measured values are in reasonable
agreement with the expectation within 1.4$\,\sigma$.

In addition, for both Rigel and Spica, $g^{(2)}(r)$ peaks for the correlations
between the two ATs were measured, each showing a contrast less than for the
zero-baseline calibration. For Rigel, a marginal peak was observed consistent
with the low expected squared visibilities ($V^{2}\sim$ 0.15) for the obtained
projected baselines. In the case of Spica, a clear $g^{(2)}$ peak was observed
with a significance of $7.5$ $\sigma$ and shown in Figure 8. The bunching area
corresponds to a measured squared visibility of $V^{2}=0.42\pm 0.06$ that
includes all data sampled over the baseline range from $\sim$ 23 to 48 m. The
expectation of the visibility of Spica is highly non-trivial as the system
consists of a short period interacting binary of two hot B-type stars [15, 16]
and a full analysis is beyond the scope of these proceedings.

##  5 Conclusion and Outlook

The I2C system has been successfully tested with several of the facilities
available on the Calern Plateau and additionally with two of the ATs of the
VLTI. These results have demonstrated the versatility of the intensity
interferometry technique to be utilized on telescopes of different optical
designs. Several technical improvements are being pursued in order to improve
sensitivity to allow for the observation of fainter systems or equivalently
increase the precision of the squared visibility on stars of similar
magnitude. Our current detectors exhibit a high quantum efficiency $>60\%$ but
with a fairly modest temporal resolution of $\sim 500$ ps. In the near future
we aim towards on-sky tests using super conducting nanowire single photon
detectors which exhibit similar, if not better, quantum efficiency but with
drastically better timing $<50\,$ps[17]. Furthermore, in order to have a
sensitivity competitive with other HBT observatories using Cherenkov
telescopes and ultimately with direct interferometry, it is necessary to
implement spectral multiplexing. The sensitivity of an intensity
interferometer is to first-order independent of the width of the spectral
passband. Therefore, the sensitivity can be improved by performing $g^{(2)}$
measurements in many spectral passbands and co-adding the results [18]. Some
initial tests towards multiplexed observations are being explored in the
laboratory using the spectrograph described earlier, as well as with
comercially available monochromators. In addition, we will work towards
scaling our system to work with 3 or 4 telescopes in order to provide improved
simultaneous coverage of the uv-plane during observations.

###### Acknowledgements.

We acknowledge the financial support of the Région PACA (project I2C), the
French ANR (project I2C, ANR-20-CE31-0003), OCA, Doeblin federation and UCA
science councils grants. The authors would like to thank Andreas Kaufer for
the opportunity to test our equipment with the VLTI-AT telescopes. We
acknowledge the use of spectral data obtained from the Southern Spectroscopic
Project Observatory team555https://2spot.org/FR/. This research has made use
of the Jean-Marie Mariotti Center Aspro service666Available at
http://www.jmmc.fr/aspro

. The authors would also like to thank Jason Aufdenberg for his assistance in
our analysis of Spica data.

## References

  * [1] Hanbury Brown, R., Davis, J., and Allen, L. R., “The angular diameters of 32 stars,” MNRAS 167(1), 121 (1974).
  * [2] Labeyrie, A., “Interference fringes obtained on Vega with two optical telescopes.,” ApJ 196, L71–L75 (Mar 1975).
  * [3] LeBohec, S. and Holder, J., “Optical intensity interferometry with atmospheric cerenkov telescope arrays,” ApJ 649, 399–405 (Sep 2006).
  * [4] Acciari, V. A., Bernardos, M. I., Colombo, E., Contreras, J. L., Cortina, J., De Angelis, A., Delgado, C., Díaz, C., Fink, D., Mariotti, M., Mangano, S., Mirzoyan, R., Polo, M., Schweizer, T., and Will, M., “Optical intensity interferometry observations using the MAGIC Imaging Atmospheric Cherenkov Telescopes,” MNRAS 491, 1540–1547 (Jan. 2020).
  * [5] Abeysekara, A. U., Benbow, W., Brill, A., Buckley, J. H., Christiansen, J. L., Chromey, A. J., Daniel, M. K., Davis, J., Falcone, A., Feng, Q., Finley, J. P., Fortson, L., Furniss, A., Gent, A., Giuri, C., Gueta, O., Hanna, D., Hassan, T., Hervet, O., Holder, J., Hughes, G., Humensky, T. B., Kaaret, P., Kertzman, M., Kieda, D., Krennrich, F., Kumar, S., LeBohec, T., Lin, T. T. Y., Lundy, M., Maier, G., Matthews, N., Moriarty, P., Mukherjee, R., Nievas-Rosillo, M., O’Brien, S., Ong, R. A., Otte, A. N., Pfrang, K., Pohl, M., Prado, R. R., Pueschel, E., Quinn, J., Ragan, K., Reynolds, P. T., Ribeiro, D., Richards, G. T., Roache, E., Ryan, J. L., Santander, M., Sembroski, G. H., Wakely, S. P., Weinstein, A., Wilcox, P., Williams, D. A., and Williamson, T. J., “Demonstration of stellar intensity interferometry with the four VERITAS telescopes,” Nature Astronomy 4, 1164–1169 (July 2020).
  * [6] Dravins, D., LeBohec, S., Jensen, H., and Nuñez, P. D., “Optical intensity interferometry with the cherenkov telescope array,” APh 43, 331 – 347 (2013).  Seeing the High-Energy Universe with the Cherenkov Telescope Array - The Science Explored with the CTA.
  * [7] Guerin, W., Dussaux, A., Fouché, M., Labeyrie, G., Rivet, J.-P., Vernet, D., Vakili, F., and Kaiser, R., “Temporal intensity interferometry: photon bunching in three bright stars,” MNRAS 472(4), 4126–4132 (2017).
  * [8] Guerin, W., Rivet, J.-P., Fouché, M., Labeyrie, G., Vernet, D., Vakili, F., and Kaiser, R., “Spatial intensity interferometry on three bright stars,” MNRAS 480, 245–250 (07 2018).
  * [9] Rivet, J. P., Siciak, A., de Almeida, E. S. G., Vakili, F., Domiciano de Souza, A., Fouché, M., Lai, O., Vernet, D., Kaiser, R., and Guerin, W., “Intensity interferometry of P Cygni in the H $\alpha$ emission line: towards distance calibration of LBV supergiant stars,” MNRAS 494, 218–227 (Feb. 2020).
  * [10] de Almeida, E. S. G., Hugbart, M., Domiciano de Souza, A., Rivet, J.-P., Vakili, F., Siciak, A., Labeyrie, G., Garde, O., Matthews, N., Lai, O., Vernet, D., Kaiser, R., and Guerin, W., “Combined spectroscopy and intensity interferometry to determine the distances of the blue supergiants P Cygni and Rigel,” Monthly Notices of the Royal Astronomical Society (June 2022).  https://doi.org/10.1093/mnras/stac1617.
  * [11] Siegert, A. J. F., “On the fluctuations in signals returned by many independently moving scatterers, report / radiation laboratory, massachusetts institute of technology,” (1943).
  * [12] Ferreira, D., Bachelard, R., Guerin, W., Kaiser, R., and Fouché, M., “Connecting field and intensity correlations: The siegert relation and how to test it,” American Journal of Physics 88(10), 831–837 (2020).
  * [13] Mandel, L. and Wolf, E., [Optical Coherence and Quantum Optics ], Cambridge University Press (1995).
  * [14] Bourgès, L. and Duvert, G., “ASPRO2: get ready for VLTI’s instruments GRAVITY and MATISSE,” in [Optical and Infrared Interferometry and Imaging V ], Malbet, F., Creech-Eakman, M. J., and Tuthill, P. G., eds., 9907, 249 – 259, International Society for Optics and Photonics, SPIE (2016).
  * [15] Herbison-Evans, D., Hanbury Brown, R., Davis, J., and Allen, L. R., “A Study of $\alpha$ Virginis with an Intensity Interferometer,” MNRAS 151, 161–176 (01 1971).
  * [16] Aufdenberg, J. P. and Hammill, J. M., “Modeling the h$\alpha$ emission surrounding spica using the lyman continuum from a gravity-darkened central star,” The Astrophysical Journal 923, 10 (dec 2021).
  * [17] Chang, J., Zadeh, I. E., Los, J. W. N., Zichi, J., Fognini, A., Gevers, M., Dorenbos, S., Pereira, S. F., Urbach, P., and Zwiller, V., “Multimode-fiber-coupled superconducting nanowire single-photon detectors with high detection efficiency and time resolution,” Appl. Opt. 58, 9803–9807 (Dec 2019).
  * [18] Trippe, S., Kim, J.-Y., Lee, B., Choi, C., Oh, J., Lee, T., Yoon, S.-C., Im, M., and Park, Y.-s., “Optical multi-channel intensity interferometry - or: how to resolve O-stars in the Magellanic Clouds,” Journal of The Korean Astronomical Society 47, 235–253 (2014).