# Eigenfunction and eigenmode-spacing statistics in chaotic photonic crystal
graphs

Shukai Ma<EMAIL_ADDRESS>Quantum Materials Center, University of Maryland,
College Park, Maryland 20742, USA  Department of Physics, University of
Maryland, College Park, Maryland 20742, USA  Thomas M. Antonsen  Department of
Physics, University of Maryland, College Park, Maryland 20742, USA  Department
of Electrical and Computer Engineering, University of Maryland, College Park,
Maryland 20742-3285, USA  Steven M. Anlage  Quantum Materials Center,
University of Maryland, College Park, Maryland 20742, USA  Department of
Physics, University of Maryland, College Park, Maryland 20742, USA  Department
of Electrical and Computer Engineering, University of Maryland, College Park,
Maryland 20742-3285, USA

###### Abstract

The statistical properties of wave chaotic systems of varying dimensionalities
and realizations have been studied extensively. These systems are commonly
characterized by the statistics of the eigenmode-spacings and the statistics
of the eigenfunctions. Here, we propose photonic crystal (PC) defect waveguide
graphs as a new physical setting for chaotic graph studies. Photonic crystal
waveguides possess a dispersion relation for the propagating modes which is
engineerable. Graphs constructed by joining these waveguides possess junctions
and bends with distinct scattering properties. We present numerically
determined statistical properties of an ensemble of such PC-graphs including
both eigenfunction amplitude and eigenmode-spacing studies. Our proposed
system is compatible with silicon nanophotonic technology and opens chaotic
graph studies to a new community of researchers.

##  I I. Introduction

Wave-chaotic phenomena have been studied in various complex scattering
systems, ranging from 1D graphs Hul  _et al._ (2004); Martínez-Argüello  _et
al._ (2018); Rehemanjiang  _et al._ (2020); Lu  _et al._ (2020); Chen  _et
al._ (2020), 2D billiards Wu  _et al._ (1998); Dembowski  _et al._ (2003,
2005); Hemmady  _et al._ (2005); Kuhl  _et al._ (2008); Shinohara  _et al._
(2010); Dietz and Richter (2019) to 3D enclosures Hemmady  _et al._ (2012);
Xiao  _et al._ (2018); Ma  _et al._ (2019a, 2020a, 2020b); Frazier  _et al._
(2020). The statistical properties of many system quantities, such as the
closed system eigenvalues and the open system scattering/impedance matrices,
exhibit universal characteristics, which only depend on general symmetries
(e.g., time-reversal, symplectic) and the degree of system loss. Random Matrix
Theory (RMT) has enjoyed great success in describing the energy level spacing
statistics of large nuclei Mehta (2004); Dietz  _et al._ (2017a). The
statistics of complex systems that show time-reversal invariance (TRI) are
described by the Gaussian orthogonal ensemble (GOE) of random matrices, and
the statistics of systems showing broken time-reversal invariance are
described by the Gaussian unitary ensemble (GUE), and the statistics of
systems with TRI and antiunitary symmetry show Gaussian symplectic ensemble
(GSE) statistics So  _et al._ (1995); Hemmady  _et al._ (2005); Ma  _et al._
(2020a). It was later conjectured by Bohigas, Giannoni, and Schmit that any
system with chaotic dynamics in the classical limit would also have semi-
classical wave properties whose statistics are governed by RMT associated with
one of these three ensembles Bohigas  _et al._ (1984). In practice,
identifying RMT statistical properties in experimental data is difficult due
to the presence of non-universal features that obscure the universal
fluctuations. The Random Coupling Model (RCM) has found great success in
characterizing the statistical properties of a variety of experimental systems
by removing the non-universal effects induced by port coupling and short-orbit
effects Hemmady  _et al._ (2005); Zheng  _et al._ (2006); Hart  _et al._
(2009); Yeh  _et al._ (2010); Hemmady  _et al._ (2012); Li  _et al._ (2015);
Aurégan and Pagneux (2016); Xiao  _et al._ (2018); Zhou  _et al._ (2019); Ma
_et al._ (2020a, b).

Chaotic microwave graphs support complex scattering phenomena despite their
relatively simple structure Hul  _et al._ (2004); Martínez-Argüello  _et al._
(2018); Rehemanjiang  _et al._ (2020); Lu  _et al._ (2020); Chen  _et al._
(2020); Zhang  _et al._ (2022). We refer loosely to a chaotic graph as one in
which there are multiple paths of differing lengths from one node to another,
such that waves traveling along these paths interfere when arriving at nodes.
A microwave graph structure can be realized with coaxial cables as graph bonds
connected at T-junctions. The graph architecture allows for various useful
circuit components (such as variable phase shifters and attenuators, as well
as microwave circulators) to be incorporated into the structure Rehemanjiang
_et al._ (2020); Lu  _et al._ (2020); Chen  _et al._ (2020). Recent studies
show that non-universal statistical features exist in chaotic graph systems,
and these are hypothesized to be caused by the non-zero reflection at the
graph vertices. These reflections create trapped modes that impact the
spectral statistics of the graph Białous  _et al._ (2016); Dietz  _et al._
(2017b); Lu  _et al._ (2020); Che  _et al._ (2021).

Figure 1:  a. The open-plate view of the photonic crystal lattice with an
L-shaped defect waveguide region. b. The side view of the PC unit cell. A
dielectric rod is sandwiched between two PEC (perfect electric conductor)
surfaces. The quantities $h_{0}$ and $d_{0}$ are the height and diameter of
the rod. c. The photonic band structure from a supercell defect waveguide
simulation. The waveguide modes appear in the bulk bandgap region (delineated
by the blue dashed lines). d. The E-field profile of one of the waveguide mode
solutions (the red circle in c). The black dashed line marks the center of the
waveguide.

Here, we introduce an alternative type of chaotic graph built with photonic
crystal (PC) waveguides. Photonic crystals find extensive use in
semiconductor-based nanophotonic technology, where they act as low-loss
waveguides for visible and infrared light Johnson  _et al._ (2001). The PC
graph bonds are realized with defect waveguides, and the nodes are formed by
the waveguide junctions. As TRI is preserved and spin-1/2 behavior is absent,
chaotic PC-graph systems are expected to fall into the GOE universality class
of RMT. With numerical simulation tools, we conduct a series of statistical
tests of chaotic photonic crystal graph systems including both eigenvalue and
eigenfunction studies of closed graphs.

The unique properties of PC structures for quantum chaos studies are as
follows. First, the bond lengths of PC-graphs can be altered by means of
lithographic fabrication Johnson  _et al._ (2001). Second, the scattering
properties of PC-graph nodes can be engineered by changing the rod properties
and geometry. Moreover, PC defect waveguides and cavities are technologies
that are widely used in the field of integrated nanophotonics. To our
knowledge, PC waveguide graphs have not been previously utilized for quantum
graph studies.

The paper is organized as follows. In Section II, we introduce the design
details of the proposed PC graph structure as well as its numerical
implementation methods. We present the chaotic closed-graph mode-spacing study
in section III, and focus on the discussion of closed-graph eigenfunction
studies in section IV. Different methods of conducting eigenfunction
statistical studies are also discussed. We summarize the paper and discuss the
future applications of the proposed PC graphs in section V.

##  II II. Photonic crystal graph

A photonic crystal system consists of a regular lattice of artificial atoms
(or scatterers) whose spacing is comparable to the operating wavelength
Joannopoulos  _et al._ (1997, 2008). The material properties and the
geometrical details of the atoms are carefully designed in order to achieve a
specific functionality. A PC system is usually constructed as a 2D planar
structure, which makes it especially good for lithographic fabrication and
photonics applications. Importantly, a 2D PC can show a complete bulk bandgap
in its electromagnetic excitation spectrum Joannopoulos  _et al._ (1997,
2008). In PC-based devices, waveguides and cavities can be constructed for
example by making air defects (removing a certain number of the atoms) in the
original lattice. This creates guided propagating modes in the bulk bandgap,
which ensures that the modes are confined to the defect region. A variety of
defect waveguides can be realized by changing the atom properties Joannopoulos
_et al._ (1997). Recent photonic topological insulator (PTI) studies present
an alternative form of PC waveguide using the interface between two different
topological domains Xiao  _et al._ (2016); Lai  _et al._ (2016); Ma  _et al._
(2019b); Ma and Anlage (2020); Liu  _et al._ (2021). Such PTI-based waveguides
can also be utilized as the bonds for realizing graphs. The interesting
properties of the topologically protected modes may further enhance the
quantum graph studies. Here, we will utilize the defect waveguide modes to
build chaotic graph structures. Note that these PC defect waveguide graphs are
distinctly different from photonic crystal slabs and billiards Dietz and
Richter (2019); Maimaiti  _et al._ (2020); Zhang and Dietz (2021), and from
coupled resonant dielectric cylinders Rehemanjiang  _et al._ (2020); Reisner
_et al._ (2021) studied previously.

Figure 2:  a. The simulated S-parameters of the PC-graph right-angle junction
(solid), as well as the S-parameter from a $3.5mm$ uniform coaxial cable
2-port right-angle connector (dashed), which serves as a reference. Inset: the
schematics of the PC right-angle junction and the 2-port coaxial cable
junction. b. The simulated S-parameters of the PC-graph T-junction (solid) and
the S-parameters from a coaxial cable 3-port ‘T’ connector (dashed). Inset:
schematics of the PC T-junction and the 3-port coaxial cable junction. All PC-
junction S-parameters are frequency averaged with a 500MHz window to remove
the effect of spurious reflections.

The construction of the chaotic PC graph starts by building a square lattice
with identical dielectric rods (the blue lattice in Fig. 1(a)). The lattice
constant is $a_{0}=36.8mm$. The dielectric rod lattice is located between two
metallic surfaces and installed in a vacuum background. The defect waveguide
is created by simply removing one row or column of the dielectric rods. The
detailed shape of the dielectric rod is shown in Fig. 1(b). The diameter of
the dielectric rod is $d_{0}=0.36a_{0}=13.2mm$ and the height is
$h_{0}=0.1a_{0}=3.68mm$. Because the PC lattice is thin in the vertical
direction ($z$-direction), the waveguide modes considered here are polarized
such that $E_{z}\neq 0$, $E_{x,y}=0$, $H_{z}=0$ and $H_{x,y}\neq 0$ (TM with
respect to $z$), and these are the only modes that propagate for frequencies
$f<12$ GHz, well above the band gap of interest here. The relative permitivity
and permeability of the dielectric rods are $\epsilon_{r}=11.56$ and
$\mu_{r}=1$. We realize the proposed PC structure numerically with COMSOL
Multiphysics Software. The presence of the waveguide mode is clearly
demonstrated by the supercell photonic band structure (PBS) simulation (Fig.
1(c)). The supercell simulation model consists of a single column of PC
lattice with Floquet periodic boundaries on the two long sides. That is,
fields on opposite long sides differ by a phase factor $\theta=k_{x}a_{0}$.
Fields on the opposite short sides are assigned totally absorbing boundary
conditions. We remove the center rod to create the defect waveguide region.
The PBS simulation is conducted by computing the system eigenmodes while
varying the wavenumber $k_{x}$ in the range $[-\pi/a_{0},\pi/a_{0}]$. As shown
in Fig. 1(c), the defect waveguide modes have emerged inside of the bulk
bandgap region from $2.5\sim 3.6$GHz. The eigenfrequencies in the band-gap
region are essentially real, as the absorbing boundaries are well separated
from the defect creating the waveguide.

We see from Fig. 1(c) that the waveguide modes in the PC have a distinct
dispersion relation (frequency versus $k_{x}$) that differs from that of a
transmission line, $\omega=k_{x}c$. The dispersion relation is parabolic at
the lower band-edge, as in a metallic waveguide. However, the group velocity
vanishes at the upper band-edge, in contrast with that of a metallic waveguide
for which there is no upper band-edge, rather just a boundary where new
propagating modes are possible Zhang  _et al._ (2022). In Fig. 1(d), the
$|E_{z}|$ profile of the entire simulation domain shows clearly that the mode
solution is indeed a guided wave because its amplitude is highly localized in
the defect region.

In addition to the propagation on the graph bonds, it is also important to
characterize the scattering properties of graph nodes Bittner  _et al._
(2013). The PC graph structure is realized by connecting multiple straight
defect waveguides with both right-angle and T-shaped junctions. We have
characterized the scattering matrix of both types of junctions with COMSOL and
CST Microwave STUDIO. For the right-angle junctions, non-zero transmission is
found only in the bulk bandgap region. As shown in Fig. 2(a), the transmission
$|S_{21}|$ and the reflection $|S_{11}|$ for right angle bends vary
systematically as functions of frequency. Thus, right-angle bends are vertices
with a non-trivial scattering matrix in a graph. The PC-graph T-junction also
presents a complex scattering profile over the entire bandgap region. As shown
in Fig. 2(b), the magnitude of the reflection (transmission) coefficient
deviates systematically from 1/3 (2/3) as a function of frequency.

It should be noted that for the PC graph junctions the sum of the squares of
the displayed scattering coefficients is close to unity. This implies that at
the frequencies of interest the photonic waveguides are operating in a single
transverse mode. For an ideal transmission line there is no reflection at a
bend. Also, at an ideal T-junction the reflection coefficient
$|S_{11}|=|-1/3|$, and the transmission coefficient is $|S_{12}|=|2/3|$. The
ideal scattering parameters apply to transmission lines whose transverse
dimensions are much smaller than an axial wavelength (equivalently the
frequency is well below the first cut-off frequency of non-TEM modes.). The
calculated scattering properties of an air-filled coaxial ($r_{o}=2$ mm,
$r_{i}=0.86$ mm) cable bend and T-junction are shown as well in Fig. 2, and,
by contrast, they show minimal variation in the frequency range of interest.
The cut-off frequency of the $TE_{11}$ mode in this cable is 33.34 GHz, well
above the operating frequency. Consequently, the scattering parameters are
close to the ideal values. Graphs constructed of metallic waveguides, as
opposed to transmission lines, can be expected to have scattering parameters
with characteristics similar to PC graphs for frequencies between the first
and second lowest cut-off frequencies where single mode operation is possible
Zhang  _et al._ (2022). These nontrivial junction scattering parameters are a
second unique feature of PC graphs, and can be expected to influence the
properties of the graph eigenmodes. We note that alternative methods of making
waveguide bends, or connecting waveguides, can be applied, for example by
removing or adding dielectric rods at the right-angle turn, in order to tune
the transmission property of the waveguide joints Mekis  _et al._ (1996). One
may adjust the degree of wave localization of graph nodes by engineering the
scattering properties of the PC waveguide junctions. The engineering of the
node reflection and transmission properties is beyond the scope of this paper.

We simulate the closed PC graphs with the COMSOL eigenvalue solver. The graph
topology is that of a flattened tetrahedral graph having $14$ straight
segments and $13$ junctions (including both 3-way junctions and right-angle
junctions). The four exterior faces of the parallel-plate PC structure (see
Fig. 3(a)) are assigned totally absorbing boundary conditions. The total
length of the simulated graph is on the scale of $\sim 9$ m, which hosts about
80 eigenmodes (within the bulk PC bandgap) in a typical realization. We note
that one is able to decrease the mode-spacing by enlarging the size of the PC
graph, and we choose the current system scale due to limited computational
power. We have created a statistical ensemble of chaotic PC graphs by changing
the length of the bonds for a given graph topology.

A particular eigenmode solution of a graph (shown in Fig. 3(a)) is shown in
Fig. 3(b). It is clear that the graph mode is localized to the defect
waveguide region and displays a longitudinal sinusoidal standing wave pattern.
The mode amplitude on a single bond is uniform but varies between different
bonds (shown by the color differences), where bonds are defined as straight
waveguide regions between 3-way junctions and right-angle bends.

Figure 3:  a. Schematic diagram of the PC waveguide graph system. The graph
bonds are shown as clear gray channels, and the rectangular lattice is shown
as a lattice of dielectric rods (black circles). b. The simulated $|E_{z}|$
profile of the graph system eigenmode at 2.8GHz. c shows the statistics of the
normalized mode spacing from graph simulation (histogram). The theoretical
predictions for the GOE and GUE systems are shown in red and yellow curves. d
shows the spectral rigidity $\Delta_{3}(L)$ of the normalized spectrum from
graph simulation (blue circles). The theoretical predictions for the GOE, GUE,
and the Poisson systems are shown in red, yellow, and black curves. e shows
the statistics of the consecutive mode spacing ratio $r$ (histogram) and the
theoretical predictions for the GOE and GUE systems.

##  III III. Eigenmode Spacing Statistical Analysis

In order to characterize the statistical properties of PC defect waveguide
graphs, we start by conducting the nearest neighbor mode spacing analysis of
the proposed PC graphs. An ensemble of 10 different graph realizations is
studied numerically, and we obtained $\sim 800$ eigenfrequency values from the
ensemble. The graph topology ranges from $13\sim 16$ straight segments and
$11\sim 14$ junctions. We note that this ensemble of graphs is conducted with
the eigenvalue simulation in COMSOL. The simulation results are used for mode-
spacing studies here, as well as the eigenfunction statistics in Sec. IV. The
model is composed of lossless dielectric structures. However, the presence of
radiating boundary conditions on the perimeter of the finite-size structure
shown in Fig. 3(a) will give rise to a finite loss. The quality factors of the
simulated graph modes are on the order of $10^{5}$ to $10^{7}$, allowing for
clear identification of all eigenfrequencies. The graph topology is kept fixed
in the eigenfunction studies discussed below. Because the topology and the
total length $L_{tot}$ of each graph realization are different, we normalized
the system eigenmode solutions using the following method. We considered 10
different graphs. For each graph, we calculated eigenfrequencies in the range
2.8 GHz to 3.6 GHz and placed the eigenfrequencies in ascending order. The
eigenfrequencies were then converted to eigen-wavenumbers ($k_{i}$, $i$ is the
mode index) using the dispersion relation for the PC-waveguide PBS shown in
Fig. 1(c). This results in a tabulation of consecutive wave-numbers. We then
fit the mode-number index to the wave-number using a quadratic fitting
function $n(k)=c_{2}k^{2}+c_{1}k+c_{0}$, where $c_{0,1,2}$ are fitting
parameters and $n(k)$ is the total number of modes below $k$ Rehemanjiang  _et
al._ (2016). For a typical graph we found $c_{2}=0.0054m^{2},c_{1}=2.21m$ and
$c_{0}=-55.28$. The small value of $c_{2}$ indicates the $n(k)$ is nearly
linear in $k$ as expected for quasi-one-dimensional structures. We then
evaluate the fitting functions at the calculated $k$-values for each graph,
creating a list of normalized eigenvalues $e_{i}=n_{fit}(k_{i})$. The nearest
neighbor spacings are finally computed as $s_{i}=e_{i+1}-e_{i}$.

The distribution of the normalized nearest-neighbor mode spacing values of the
entire ensemble is shown as the histogram in Fig. 3(c). We have included the
approximate theoretical predictions of the mode-spacing statistics for the GOE
and GUE systems in the figure. These theoretical predictions are based on RMT
Berry (1981); Kottos and Smilansky (1997, 1999), which are methods of
understanding the universal statistical properties of wave-chaotic systems
Casati  _et al._ (1980); Bohigas  _et al._ (1984). As shown in Fig. 3(c), the
distribution of the normalized mode-spacing matches reasonably well with the
theoretical prediction for the GOE universality class. Good agreement between
the graph nearest neighbor mode-spacing statistics and the RMT theoretical
prediction is also reported in various works on graphs, although long-range
statistical quantities tend to show non-universal behavior due to the trapped-
mode issue mentioned above Białous  _et al._ (2016); Dietz  _et al._ (2017b);
Lu  _et al._ (2020); Che  _et al._ (2021).

To test the long-range statistical properties of the eigenmodes, we compute
the spectral rigidity of the graph spectrum $\Delta_{3}(L)$, where $L$
measures the length of a segment in the normalized spectrum. Our procedure is
as follows. The quantity $\Delta_{3}(e,L)$ representing the deviation of the
spectral staircase from a straight line, is evaluated for each value of $L$ by

$\Delta_{3}(e,L)=\frac{1}{L}min_{a,b}\int_{e}^{e+L}[N(x)-ax-b]^{2}dx,$ (1)

where $N(e)=n(k)$ and $e$ is chosen at random. We use 100 values of e for each
graph. Then we average the values of e over 10 graph realizations Haq  _et
al._ (1982). The resulting values of $\Delta_{3}(e,L)$ are then plotted in
Fig. 3(d) along with theoretical predictions based on GOE, and GUE random
matrices Che  _et al._ (2021), and the predictions of Poisson statistics Berry
(1985) (appropriate for classically integrable systems).

The $\Delta_{3}(L)$ curve in Fig. 3(d) is linear in $L$ for small values of
$L$. As shown, systems of all symmetries (Poisson, GOE, GUE) show a linear-
in-$L$ dependence for $L<1$ Berry (1985), and the data is consistent with this
trend. Chaotic systems are then expected to show a logarithmic increase in
$\Delta_{3}(L)$ Berry (1985), and this is evident in the data up to
approximately $L=5$. It is also predicted that $\Delta_{3}(L)$ will saturate
due to the presence of short orbits in any experimental realization of a wave
chaotic system Berry (1985), and this is also seen in the data. Previous
studies find that the long-range statistics of quantum graphs usually show
some degree of agreement with the GOE prediction in the range of roughly
$L\in[2,5]$, along with deviations beyond this range Hul  _et al._ (2004);
Białous  _et al._ (2016); Dietz  _et al._ (2017b). We observe in Fig. 3(d) a
deviation from the GOE RMT predictions, especially for larger $L$, similar to
cable graph studies in Ref. Hul  _et al._ (2004). Because the operating
bandwidth of the proposed PC-graph system is relatively narrow (Fig. 1 (c)),
we believe the long-range statistics may also be affected by finite-size
effects because the number of available eigenmodes is bounded above and below
by the edge of the PC bandgap. The solution is to increase the graph size
substantially, but this is beyond our computational resources. In summary, the
$\Delta_{3}(L)$ behavior is very similar to that observed in wave chaotic
systems of higher dimensionality.

| PC-Graph | GOE | GUE | GSE
---|---|---|---|---
$\left<r\right>$ | 1.89 | 1.75 | 1.37 | 1.18
$\left<\tilde{r}\right>$ | 0.51 | 0.54 | 0.60 | 0.68
Table 1:  The summarized consecutive mode spacing ratios of the photonic
crystal graph (PC-Graph) system, and the theoretical predictions for the GOE,
GUE, and GSE systems Atas  _et al._ (2013).

In addition to the mode-spacing distribution test, we note that the method of
consecutive mode spacing ratios $r_{i}=\frac{s_{i}}{s_{i-1}}$ and
$\tilde{r}_{i}=min\left(r_{i},\frac{1}{r_{i}}\right)$ can also be adopted Atas
_et al._ (2013). Here the ratio statistic has the feature that it is a measure
of the correlation between adjacent spacings. One advantage of this type of
statistical study is that it does not require unfolding the spectrum Atas  _et
al._ (2013). However, here we use the unfolded spectrum $s_{i}$ to compute the
spacing ratios. For the PC graphs, the averaged values of
$\left<r\right>=1.89$ and $\left<\tilde{r}\right>=0.51$, are closer to the GOE
theoretical predictions than the GUE prediction (Table. 1). We note that the
high value of $\left<r\right>$ in the PC-graph is related to the correlation
between consecutive level-spacings. We further present the distribution of
mode-spacing ratios $r$ of the PC graphs and corresponding theoretical
predictions of the GOE and GUE systems in Fig. 3(e). The statistics of $r$
match reasonably well with the GOE prediction.

##  IV IV. Eigenfunction Analysis

We next study the statistics of the closed PC graph eigenfunctions. The
eigenfunction statistics have been studied experimentally in 2D chaotic
systems by probing the electromagnetic (EM) standing wave field inside a
microwave cavity Kudrolli  _et al._ (1995); Wu  _et al._ (1998); Gokirmak  _et
al._ (1998); Chung  _et al._ (2000); Kuhl  _et al._ (2000); Hlushchuk  _et
al._ (2001). In those studies, the experimental probability amplitude
distribution and two-point correlation function agree well with the random
plane wave conjecture. Here, the wave properties of the entire closed PC graph
(nodes and bonds) can be faithfully simulated. We employ the same eigenvalue
simulation model used in the mode-spacing studies above. For a graph system,
the telegrapher’s equation is formally equivalent to the 1D Schrödinger
equation, where the wavefunction $\psi$ is represented by the wave voltage.
For thin parallel plate waveguides, the wave voltage difference between the
top and bottom metallic plates is represented by the $E_{z}$ value at the
middle cutting plane (at the height of $z=h_{0}/2$). Because the PC graph
bonds have a finite width, the bond eigenfunction will be evaluated along a 1D
line at the center of the waveguide (shown as the black dashed line in Fig.
1(d)). We next examine two PC graph eigenfunction characterization methods.

Method I: grid-wise representation. For each eigenmode of each graph
realization, we will use the entire set of the graph bond $|E_{z}|^{2}$ vs.
longitudinal position along the center-line of the waveguide data points to
represent the eigenfunction. The name ‘grid-wise’ comes from the granular
nature of this method, where the total number of eigenfunction data points is
inversely proportional to the computational grid size. We study the
wavefunction statistics of the PC graph by computing the distribution of the
normalized probability density $v$, which is defined as the square of the
eigenfunction values

$v_{j}=|\psi(r_{j})|^{2}=\frac{|E_{z}(r_{j})|^{2}\cdot
L_{tot}}{\sum_{j}|E_{z}(r_{j})|^{2}\cdot\Delta L_{j}},$ (2)

where $L_{tot}$ is the total length of the graph, $r_{j}$ and $\Delta L_{j}$
are the location and the grid size of the $j^{\prime}th$ grid point,
$|E_{z}(r_{j})|$ is the z-directed electric field magnitude at the
$j^{\prime}th$ grid point, and the summation in the denominator runs over
every graph grid point. In our simulation, the grid size $\Delta L_{j}\sim
0.05\lambda_{op}$ where $\lambda_{op}=10cm$ is the operating wavelength at
$3GHz$. The distribution of the probability density values is computed using
the data from all simulated eigenmode solutions (48 to be specific) from a
single realization, and the results are shown in Fig. 4(a), and discussed
below.

Figure 4:  a. Probability amplitude distribution of photonic crystal defect
waveguide graph eigenmodes values ($v$) obtained from the grid-wise method
(symbols connected by blue line). The theoretical prediction for systems with
GOE statistics is in red. b. shows the statistics of the normalized bond
amplitude value $|a|$ (symbols connected by blue line) and the Gaussian
distribution (red dashed). The inset shows the distribution of the quantity
$|a|^{2}$ and the Porter-Thomas theoretical prediction for the GOE system.

Method II: bond-value representation. For each graph realization, the
eigenfunction of a mode is represented by a set of ‘bond-values’
$E_{z}(b_{m})$ which is defined as the amplitude of the standing-wave
wavefunction on the graph bond $b_{m}$. The quantity $m$ is the index of the
bond and runs from 1 to 14. The standing wave on the bond is made up of two
counter propagation waves $\psi(x)=a_{m}\,e^{ikx}+a_{m}^{*}\,e^{-ikx}$, where
$a_{m}$ is the wavefunction amplitude at bond $b_{m}$ and $x$ is the distance
from a vertex along the bond. We first conduct a sine-fit of the raw
$E_{z}(x)$ values on each bond, which yields the amplitude, and the value of
$|a_{m}|$ is obtained as 1/2 of the amplitude value. The normalization process
follows the same method as in Ref. Kaplan (2001) which ensures that
$\sum_{m}L(b_{m})|a_{m}|^{2}=L_{tot}$ where $L(b_{m})$ is the length of bond
$b_{m}$. Here the distribution of the probability density values is computed
over $14\times 78$ data points, where 14 is the number of the bonds and 78 is
the number of eigenmode solutions from one graph configuration, and the
results for $P(|a|)$ are shown in Fig. 4(b).

The random plane wave hypothesis underlying chaotic eigenfunction treatments
asserts that the eigenfunctions in a 2D or 3D cavity can be thought of as
superpositions of plane waves with random directions and random phases.
Consequently, the value of an eigenfunction at any point is a sum of many
random waves and is thus a zero mean Gaussian random variable Berry (1981).
For systems with TRS the plane waves come in counter-propagating pairs, and
the wave function is real. For systems with TRS broken, wave directions are
completely random, and the wave function is complex with independent real and
imaginary parts. For the PC graph the situation is changed in that at any
point in the graph there are only two counter-propagating waves of equal
amplitude. The amplitudes of these waves, $|a_{m}|$ will be different on the
different bonds leading to statistical variations.

Figure 5:  a. The inverse participation ratio (IPR) of graph eigenmodes
obtained from the bond-value method, computed based on method IIa. The
simulated $|E_{z}|$ profile of the red (black) circled mode in a is shown in c
(d). b. The inverse participation ratio (IPR) of graph eigenmodes is computed
based on method IIb. The simulated $|E_{z}|$ profile of the red (black)
circled mode in b is shown in e (f).

We first discuss the eigenfunction statistics obtained using the grid-wise
method in Fig. 4(a). For waves in a 2D or 3D cavity the distribution of values
of $v$ based on the random plane wave hypothesis is $P(v)=(2\pi
v)^{-1/2}e^{-v/2}$. This is shown as a solid red line in Fig. 4(a). We find
that the simulated PC graph result matches to some degree the GOE prediction
in that for $v>1$ the curve approaches a straight line on a log-linear plot,
and for $v$ small there is an indication of square root singularity. We notice
that the mismatch between the graph data and the GOE prediction exists at the
small and large probability density values, indicating an absence of both very
small and very large $v$ values. One possibility is that this method tends to
over-count the appearance of medium-sized eigenfunction values (similar to the
systematic errors experienced in Ref. Wu  _et al._ (1998)). It may also
indicate that the data simply does not match the GOE prediction. In addition
to the discrete eigenfunction imaging method, we note that the eigenfunction
statistics may also be tested through resonance width distributions in
transmission measurements (Porter-Thomas statistics) Alt  _et al._ (1995);
Dembowski  _et al._ (2005).

Here we now use Method II to study the distribution of normalized bond values,
$|a|$ that are displayed in Fig. 4(b). Also plotted in Fig. 4(b) is a Gaussian
distribution, which would be characteristic of an eigenfunction of a GOE
random matrix. However, a clear deviation from Gaussian statistics is seen for
low amplitudes, similar to the deviation observed with Method I. The inset to
Fig. 4(b) shows a plot of $P(|a|^{2})$ compared to the Porter-Thomas
distribution expected for GOE systems. There is an agreement between the two
except for deviations between data and the model distribution at small values
of $|a|^{2}$, similar to that seen in Fig. 4(a). Together, these results
suggest that the wavefunction statistics of this simple graph are not fully
consistent with the random plane wave hypothesis.

We next present the inverse participation ratio (IPR) computation based on the
bond-value method (Method II) in Fig. 5. IPR is a measure of the degree of
localization for a wave function Kaplan (2001); Hul  _et al._ (2009), and can
be used to quantify the degree of deviation of wavefunction statistics from
the random plane wave hypothesis. Based on its IPR value, the eigenfunction
behavior varies between two limits, namely the maximum ergodic limit where the
wave function occupies each graph bond with equal chance, and the maximum
localization limit where the eigenmode is confined to only one bond. RMT
predicts an IPR value by assuming Gaussian random fluctuation of the
eigenfunctions. Two different IPR definitions are tested here. Method IIa
follows the definition in Ref. Kaplan (2001), where the IPR value for each
graph mode is evaluated using the formula
$IPR_{IIa}=\langle\left|a_{m}\right|^{4}\rangle$. As shown in Fig. 5(a), the
IPR values of the photonic crystal graph modes vary erratically but lie
between the maximum ergodic limit ($IPR_{IIa}$ = 1) and the RMT prediction
limit ($IPR_{IIa}$ = 2) Kaplan (2001). In the maximum localization limit the
$IPR_{IIa}=B$, where $B=14$ is the number of graph bonds, and the results are
far from this limit. Method IIb follows the definition in Ref. Hul  _et al._
(2009) where
$IPR_{IIb}=\sum_{m}\left|\tilde{a}_{m}\right|^{4}/\left[\sum_{m}\left|\tilde{a}_{m}\right|^{2}\right]^{2}$.
The quantity $\tilde{a}_{m}$ is the un-normalized bond wavefunction amplitude.
Here, the graph IPR values also vary erratically from mode to mode (Fig.
5(b)), but lie well below the maximum localization limit ($IPR_{IIb}$ = 1) and
closer to the RMT prediction limit ($IPR_{IIb}=1/B=0.07$) Hul  _et al._
(2009).

The conclusions we draw from the above two methods are not exactly the same.
For Method IIa, two exemplary eigenfunction profiles are shown in Figs. 5(c)
and d, which correspond to the RMT and ergodic limits, respectively. One may
directly spot the different nature of these two modes based on their
eigenfunction patterns, which is a convenient feature of the photonic crystal
graph system. For Method IIb, we present the eigenfunction profile of two
neighboring graph modes in Figs. 5(e) and (f). The eigenmode in Fig. 5(e) has
a larger value of $IPR_{IIb}$ and shows a strongly localized distribution.
That in Fig. 5(f) has a small value of $IPR_{IIb}$ and is more evenly
distributed over the bonds. The average value of IPR over all the modes is
more meaningful in this case, and we have $\langle IPR_{IIa}\rangle=1.39$ and
$\langle IPR_{IIb}\rangle=0.10$ which are closer to the ergodic and RMT
limits, respectively, than to the localization limit. Previous work Kaplan
(2001); Hul  _et al._ (2004) on IPR on quantum graphs shows the surprising
results that larger graphs, with the number of vertices $>$ 10, tend to show
strong deviations from RMT predictions, while smaller graphs show better
agreement. We note that by computing the $IPR_{IIa}$ on single loops rather
than the entire graph, one may identify the modes trapped inside these loops
from $IPR_{IIa}\sim B_{sl}$, where $B_{sl}$ is the total bond number of the
single loop. This study is discussed in the Appendix.

Taking into account the eigenfunction and IPR statistics, it would appear that
PC graphs are close to the random plane wave condition for tetrahedral-like
graphs, but clear systematic differences from RMT predictions remain. The
complex vertex scattering properties of right-angle bends and T-junctions,
along with the $\omega(k)$ dispersion relation of PC defect waveguide modes,
suggest that PC defect graphs may be very effective for future wave chaotic
studies.

##  V V. Conclusion

To conclude, we have designed and simulated an alternative chaotic graph
system with photonic crystal defect waveguides that show a dispersion relation
that is different from that utilized in coaxial cable microwave graphs. We
show that a series of statistical studies can be carried out on an ensemble of
closed graphs, including nearest-neighbor spacing statistics and eigenfunction
statistics studies. Because both the graph bonds and nodes can be probed, one
may better analyze the non-universal features of chaotic graphs using a PC
system. It is possible to adjust the degree of wave localization by
engineering the scattering properties of the PC waveguide junctions. This
property of PC waveguides may facilitate further studies of localization
phenomena in graphs, for example, the emergence or suppression of trapped
modes.

##  VI Acknowledgements

This work was supported by ONR under Grant No. N000141912481, AFOSR COE Grant
FA9550-15-1-0171, DARPA WARDEN Grant HR00112120021, and the Maryland Quantum
Materials Center.

*

##  Appendix A Appendix A: Sub-graph IPR Analysis

Figure 6:  a. The inverse participation ratio (IPR) of the left side of the
graph, computed based on method IIa. The simulated $|E_{z}|$ profile of the
circled modes in a is shown in c (d). b. The IPR of the right side of the
graph is computed based on method IIa. The solid and dashed lines in a. and b.
refer to the RMT and the maximum ergodic limits, respectively.

Here we give a brief discussion on applying the inverse participation ratio
(IPR) analysis to only a sub-region of the graph structure. Although the
original definition of IPR is for the entire graph, we may apply the same
computation in the context of subgraphs. Here we conduct subgraph IPR studies
with bonds that only belong to the left half (one single loop) or the right
half (two loops) of the PC graph shown in Fig. 6. The left half loop IPR is
shown in Fig. 6 (a), computed based on method IIa. The IPRs computed with
method IIb are not shown here. The majority of the eigenmodes have IPRs
between the RMT (IPR = 2) and the maximum ergodic limits (IPR = 1). In the
left half loop IPR test, we find three modes that have IPR values that reach
the maximum localization limit (IPR = B, B = 7 for the left single loop). We
show the eigenfunction profile of the two circled modes in Fig. 6 (c) and (d).
These two modes indeed show a high degree of localization behavior, where the
eigenfunction amplitude is mainly concentrated on one short bond. Fig. 6 (b)
shows the IPR analysis of the right half graph. Here the modes with a higher
degree of localization do not stand out in this subgraph IPR analysis. Thus we
note that applying the IPR analysis to single loops can bring out the modes
with a higher degree of localization.

## References

  * Hul  _et al._ (2004) O. Hul, S. Bauch, P. Pakoński, N. Savytskyy, K. Życzkowski, and L. Sirko, Physical Review E 69, 056205 (2004).
  * Martínez-Argüello  _et al._ (2018) A. M. Martínez-Argüello, A. Rehemanjiang, M. Martínez-Mares, J. A. Méndez-Bermúdez, H.-J. Stöckmann, and U. Kuhl, Physical Review B 98, 075311 (2018).
  * Rehemanjiang  _et al._ (2020) A. Rehemanjiang, M. Richter, U. Kuhl, and H.-J. Stöckmann, Physical Review Letters 124, 116801 (2020).
  * Lu  _et al._ (2020) J. Lu, J. Che, X. Zhang, and B. Dietz, Physical Review E 102, 022309 (2020).
  * Chen  _et al._ (2020) L. Chen, T. Kottos, and S. M. Anlage, Nature Communications 11, 5826 (2020).
  * Wu  _et al._ (1998) D. H. Wu, J. S. A. Bridgewater, A. Gokirmak, and S. M. Anlage, Physical Review Letters 81, 2890 (1998).
  * Dembowski  _et al._ (2003) C. Dembowski, B. Dietz, H.-D. Gräf, A. Heine, F. Leyvraz, M. Miski-Oglu, A. Richter, and T. H. Seligman, Physical Review Letters 90, 014102 (2003).
  * Dembowski  _et al._ (2005) C. Dembowski, B. Dietz, T. Friedrich, H.-D. Gräf, H. L. Harney, A. Heine, M. Miski-Oglu, and A. Richter, Physical Review E 71, 046202 (2005).
  * Hemmady  _et al._ (2005) S. Hemmady, X. Zheng, E. Ott, T. M. Antonsen, and S. M. Anlage, Physical Review Letters 94, 014102 (2005).
  * Kuhl  _et al._ (2008) U. Kuhl, R. Höhmann, J. Main, and H.-J. Stöckmann, Physical Review Letters 100, 254101 (2008).
  * Shinohara  _et al._ (2010) S. Shinohara, T. Harayama, T. Fukushima, M. Hentschel, T. Sasaki, and E. E. Narimanov, Physical Review Letters 104, 163902 (2010).
  * Dietz and Richter (2019) B. Dietz and A. Richter, Physica Scripta 94, 014002 (2019).
  * Hemmady  _et al._ (2012) S. Hemmady, T. M. Antonsen, E. Ott, and S. M. Anlage, IEEE Transactions on Electromagnetic Compatibility 54, 758 (2012).
  * Xiao  _et al._ (2018) B. Xiao, T. M. Antonsen, E. Ott, Z. B. Drikas, J. G. Gil, and S. M. Anlage, Physical Review E 97, 062220 (2018).
  * Ma  _et al._ (2019a) S. Ma, B. Xiao, R. Hong, B. Addissie, Z. Drikas, T. Antonsen, E. Ott, and S. Anlage, Acta Physica Polonica A 136, 757 (2019a).
  * Ma  _et al._ (2020a) S. Ma, B. Xiao, Z. Drikas, B. Addissie, R. Hong, T. M. Antonsen, E. Ott, and S. M. Anlage, Physical Review E 101, 022201 (2020a).
  * Ma  _et al._ (2020b) S. Ma, S. Phang, Z. Drikas, B. Addissie, R. Hong, V. Blakaj, G. Gradoni, G. Tanner, T. M. Antonsen, E. Ott, and S. M. Anlage, Physical Review Applied 14, 14022 (2020b).
  * Frazier  _et al._ (2020) B. W. Frazier, T. M. Antonsen, S. M. Anlage, and E. Ott, Physical Review Research 2, 043422 (2020).
  * Mehta (2004) M. L. Mehta,  _Random matrices_ (Elsevier, 2004).
  * Dietz  _et al._ (2017a) B. Dietz, A. Heusler, K. H. Maier, A. Richter, and B. A. Brown, Physical Review Letters 118, 012501 (2017a).
  * So  _et al._ (1995) P. So, S. M. Anlage, E. Ott, and R. N. Oerter, Physical Review Letters 74, 2662 (1995).
  * Bohigas  _et al._ (1984) O. Bohigas, M. J. Giannoni, and C. Schmit, Physical Review Letters 52, 1 (1984).
  * Zheng  _et al._ (2006) X. Zheng, S. Hemmady, T. M. Antonsen, S. M. Anlage, and E. Ott, Physical Review E 73, 046208 (2006).
  * Hart  _et al._ (2009) J. A. Hart, T. M. Antonsen, and E. Ott, Physical Review E 80, 041109 (2009).
  * Yeh  _et al._ (2010) J.-H. Yeh, J. A. Hart, E. Bradshaw, T. M. Antonsen, E. Ott, and S. M. Anlage, Physical Review E 82, 041114 (2010).
  * Li  _et al._ (2015) X. Li, C. Meng, Y. Liu, E. Schamiloglu, and S. D. Hemmady, IEEE Transactions on Electromagnetic Compatibility 57, 448 (2015).
  * Aurégan and Pagneux (2016) Y. Aurégan and V. Pagneux, Acta Acustica united with Acustica 102, 869 (2016).
  * Zhou  _et al._ (2019) M. Zhou, E. Ott, T. M. Antonsen, and S. M. Anlage, Chaos: An Interdisciplinary Journal of Nonlinear Science 29, 033113 (2019).
  * Zhang  _et al._ (2022) W. Zhang, X. Zhang, J. Che, J. Lu, M. Miski-Oglu, and B. Dietz, Physical Review E 106, 044209 (2022).
  * Białous  _et al._ (2016) M. Białous, V. Yunko, S. Bauch, M. Ławniczak, B. Dietz, and L. Sirko, Physical Review Letters 117, 144101 (2016).
  * Dietz  _et al._ (2017b) B. Dietz, V. Yunko, M. Białous, S. Bauch, M. Ławniczak, and L. Sirko, Physical Review E 95, 052202 (2017b).
  * Che  _et al._ (2021) J. Che, J. Lu, X. Zhang, B. Dietz, and G. Chai, Physical Review E 103, 042212 (2021).
  * Johnson  _et al._ (2001) S. Johnson, A. Mekis, S. Fan, and J. Joannopoulos, Computing in Science and Engineering 3, 38 (2001).
  * Joannopoulos  _et al._ (1997) J. D. Joannopoulos, P. R. Villeneuve, and S. Fan, Nature 386, 143 (1997).
  * Joannopoulos  _et al._ (2008) J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade,  _Photonic Crystals_ (Princeton University Press, 2008).
  * Xiao  _et al._ (2016) B. Xiao, K. Lai, Y. Yu, T. Ma, G. Shvets, and S. M. Anlage, Physical Review B 94, 195427 (2016).
  * Lai  _et al._ (2016) K. Lai, T. Ma, X. Bo, S. Anlage, and G. Shvets, Scientific Reports 6, 28453 (2016), 1601.01311 .
  * Ma  _et al._ (2019b) S. Ma, B. Xiao, Y. Yu, K. Lai, G. Shvets, and S. M. Anlage, Physical Review B 100, 085118 (2019b).
  * Ma and Anlage (2020) S. Ma and S. M. Anlage, Applied Physics Letters 116, 250502 (2020).
  * Liu  _et al._ (2021) Y. G. N. Liu, P. S. Jung, M. Parto, D. N. Christodoulides, and M. Khajavikhan, Nature Physics 17, 704 (2021).
  * Maimaiti  _et al._ (2020) W. Maimaiti, B. Dietz, and A. Andreanov, Physical Review B 102, 214301 (2020).
  * Zhang and Dietz (2021) W. Zhang and B. Dietz, Physical Review B 104, 064310 (2021).
  * Reisner  _et al._ (2021) M. Reisner, M. Bellec, U. Kuhl, and F. Mortessagne, Optical Materials Express 11, 629 (2021).
  * Bittner  _et al._ (2013) S. Bittner, B. Dietz, M. Miski-Oglu, A. Richter, C. Ripp, E. Sadurní, and W. P. Schleich, Physical Review E 87, 042912 (2013).
  * Mekis  _et al._ (1996) A. Mekis, J. C. Chen, I. Kurland, S. Fan, P. R. Villeneuve, and J. D. Joannopoulos, Physical Review Letters 77, 3787 (1996).
  * Rehemanjiang  _et al._ (2016) A. Rehemanjiang, M. Allgaier, C. H. Joyner, S. Müller, M. Sieber, U. Kuhl, and H.-J. Stöckmann, Physical Review Letters 117, 064101 (2016).
  * Berry (1981) M. V. Berry, European Journal of Physics 2, 91 (1981).
  * Kottos and Smilansky (1997) T. Kottos and U. Smilansky, Physical Review Letters 79, 4794 (1997).
  * Kottos and Smilansky (1999) T. Kottos and U. Smilansky, Annals of Physics 274, 76 (1999).
  * Casati  _et al._ (1980) G. Casati, F. Valz-Gris, and I. Guarnieri, Lettere al Nuovo Cimento 28, 279 (1980).
  * Haq  _et al._ (1982) R. Haq, A. Pandey, and O. Bohigas, Physical Review Letters 48, 1086 (1982).
  * Berry (1985) M. V. Berry, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 400, 229 (1985).
  * Atas  _et al._ (2013) Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux, Physical Review Letters 110, 084101 (2013).
  * Kudrolli  _et al._ (1995) A. Kudrolli, V. Kidambi, and S. Sridhar, Physical Review Letters 75, 822 (1995).
  * Gokirmak  _et al._ (1998) A. Gokirmak, D.-H. Wu, J. S. A. Bridgewater, and S. M. Anlage, Review of Scientific Instruments 69, 3410 (1998).
  * Chung  _et al._ (2000) S.-H. Chung, A. Gokirmak, D.-H. Wu, J. S. A. Bridgewater, E. Ott, T. M. Antonsen, and S. M. Anlage, Physical Review Letters 85, 2482 (2000).
  * Kuhl  _et al._ (2000) U. Kuhl, E. Persson, M. Barth, and H.-J. Stöckmann, The European Physical Journal B 17, 253 (2000).
  * Hlushchuk  _et al._ (2001) Y. Hlushchuk, L. Sirko, U. Kuhl, M. Barth, and H.-J. Stöckmann, Physical Review E 63, 046208 (2001).
  * Kaplan (2001) L. Kaplan, Physical Review E 64, 036225 (2001).
  * Alt  _et al._ (1995) H. Alt, H. D. Gräf, H. L. Harney, R. Hofferbert, H. Lengeler, A. Richter, P. Schardt, and H. A. Weidenmüller, Physical Review Letters 74, 62 (1995).
  * Hul  _et al._ (2009) O. Hul, P. Šeba, and L. Sirko, Physical Review E 79, 066204 (2009).


# Losses, Dissonances, and Distortions

Pablo Samuel Castro
Google Research, Brain Team
<EMAIL_ADDRESS>

##  1 Introduction

In recent years, there has been a growing interest in using machine learning
models for creative purposes. In most cases, this is with the use of large
generative models which, as their name implies, can generate high-quality and
realistic outputs in music (Huang et al., 2019), images (Esser et al., 2021),
text (Brown et al., 2020), and others. The standard approach for artistic
creation using these models is to take a pre-trained model (or set of models)
and use them for producing output. The artist directs the model’s generation
by “navigating” the latent space (Castro, 2020), fine-tuning the trained
parameters (Dinculescu et al., 2019), or using the model’s output to steer
another generative process (White, 2019; Castro, 2019).

At a high-level what all these approaches are doing is converting the
numerical signal of a machine learning model’s output into art, whether
implicitly or explicitly. However, in most (if not all) cases they only do so
after the initial model has been trained. This is somewhat unfortunate, as
there are plenty of numerical signals available during the training process,
such as the loss and gradient values, that can be used for creative purposes.

In this paper I present a study in using the losses and gradients obtained
during the training of a simple function approximator as a mechanism for
creating musical dissonance and visual distortion in a solo piano performance
setting. These dissonances and distortions become part of an artistic
performance not just by affecting the visualizations, but also by affecting
the artistic musical performance. The system is designed such that the
performer can in turn affect the training process itself, thereby creating a
closed feedback loop between two processes: the training of a machine learning
model and the performance of an improvised piano piece.

##  2 Components

#### Losses and gradients

Let $f_{\theta}:X\rightarrow Y$ denote a function parameterized by a
$d$-dimensional vector of weights $\theta\in\mathbb{R}^{d}$ that aims to
approximate a “true” function $f:X\rightarrow Y$. We improve the approximation
by updating the parameters $\theta$ so as to minimize a loss function
$\mathcal{L}(f_{\theta},f)\rightarrow\mathbb{R}$. This is typically done using
gradient descent, where we use the derivative (or gradient) of the loss
function to update the parameters:
$\theta\leftarrow\theta-\alpha\nabla\mathcal{L}(f_{\theta},f)$, where $\alpha$
is a learning rate. If set properly, this process will result in
$\mathcal{L}\rightarrow 0$.

Thus, at every iteration of the learning process we have $d+1$ values at our
disposal: the $d$ partial gradients from $\nabla\mathcal{L}$ and the loss
itself. In the next sections I will describe how I use these values as part of
a performance, but of course there are an infinitude of ways that artists can
incorporate these as part of their work.

#### Cubics and Lissajous knots

In order for the learning dynamics to be clearly observed during the
performance, it is preferable that the learning process is able to converge
relatively quickly. For this reason I chose two relatively simple functions to
learn: Cubic polynomials and Lissajous knots. The polynomials are single-
valued functions $f_{a,b,c,d}:\mathbb{R}\rightarrow\mathbb{R}$, while
Lissajous knots are multi-valued functions
$g_{n_{x},n_{y},n_{z},a,b,c}:\mathbb{R}\rightarrow\mathbb{R}^{3}$, where
$n_{x}$, $n_{y}$, and $n_{z}$ are integers (see the top row of Figure 1 for an
example of each):

$\displaystyle f_{{\bf a,b,c,d}}(x)$ $\displaystyle=ax^{3}+bx^{2}+cx+d$
$\displaystyle g_{n_{x},n_{y},n_{z},{\bf a,b,c}}(t)$ $\displaystyle=\langle
cos(n_{x}t+a),cos(n_{y}t+b),cos(n_{z}t+c)\rangle$

In both cases, the parameters $\theta$ of the learned function aim to approach
the true values of $a$, $b$, $c$, and $d$ (e.g. the integer-parameters of the
Lissajous knots are not learned). We use the mean-squared error loss for both:
$\mathbb{E}_{x}\left[\sqrt{\left(f_{a,b,c,d}(x)-f_{\theta}\right(x))^{2}}\right]$
and
$\mathbb{E}_{t}\left[\sqrt{\left(g_{n_{x},n_{y},n_{z},a,b,c}(t)-g_{\theta}\right(t))^{2}}\right]$.

#### Dissonances

Music is made of the combination of individual notes played on a variety of
instruments. Each note is actually a combination of a number of pitches or
frequencies: the fundamental frequency111This is typically what is referred to
as ”the pitch” of a played note.; and a series of overtone frequencies, that
are pitches at higher frequencies than the fundamental. A well-tuned
instrument will have overtones that are multiples of the fundamental frequency
(and in this case, these are called harmonics). For example, a well-tuned A
note may have the following frequencies (one fundamental and three overtones):
$\\{440,880,1320,1760\\}$. If we detune the overtones by an amount
proportional to the loss:
$\\{440,880(1+\mathcal{L}),1320(1+\mathcal{L}),1760(1+\mathcal{L})\\}$, then
what we will hear throughout the learning process is the sound “converging” to
its well-tuned state, starting from a detuned state.

#### Distortions

In addition to creating dissonance, we can create visual distortions using the
partial gradients of $\nabla\mathcal{L}$, and two instances of this are
explored:
1) The video input is split into its RGB components and each is translated by
an amount proportional to the first three partial gradients of
$\nabla\mathcal{L}$. Thus, when fully converged, each of these gradients will
be zero, and each of the three RGB frames will be exactly superimposed,
resulting in an unaltered image.
2) The previous distortion distorted the placement of the RGB components but
kept the aspect ratios of each unaltered. In this distortion the RGB
components are unaltered, but the $(x,y)$ positions of each pixel are
distorted by an amount proportional to
$(cos(\nabla\mathcal{L}_{1}),cos(\nabla\mathcal{L}_{2}))$, where
$\nabla\mathcal{L}_{i}$ denotes the $i$-th partial derivative of
$\nabla\mathcal{L}$.

##  3 Performance and Conclusion

The above ideas are combined into a musical performance, played on a
Disklavier piano, which is a regular acoustic piano that can also send MIDI
signal to the computer. The performance is organized into 4 parts (see Figure
1 in the appendix for screenshots of each):
Part 1: Every time a bass note is played, a new polynomial is generated by
sampling the coefficients $a,b,c,d$, and a new approximant is generated by
sampling $\theta$. Every note played on the upper half of the piano induces a
gradient step, and the loss of each step is used to detune the played note’s
overtones. The target and learned polynomials are displayed on a black
background.
Part 2: Every time a chord is played on the left hand, a new target Lissajous
knot is generated by samplinig $n_{x},n_{y},n_{z},a,b,c$, and a new
approximant is generated by sampling $\theta$. Gradient steps are continuously
performed as long as the chord is held, with the loss detuning the overtones
of the notes being played.
Part 3: Same as part 2, but we also display a video of the performer in the
background and use Distortion (1) to affect the RGB channels.
Part 4: Same as part 1, but with a video of the performer in the background.
Additionally, each note played triggers a “bubble” superimposed on the video
which is distorted using Distortion (2).

This is meant to be an improvised “process piece” that is different evey time
it is performed. An example performance is provided at the following link:
https://youtu.be/Qjg0bt5hgi4.

#### Conclusion

In this work I have explored using the training process (as opposed to just
the output) of a machine learning model as a key component of an artistic
performance. Note that the models being learned were unrelated to musical
performance, yet their numerical signals can be exploited in creative ways. As
machine learning models become more ubiquituous, I hope artists of all forms
can make use of all the available signals in artistic and creative ways.

##  4 Ethical considerations

This work uses only simulated/mathematical signals as training data, video
footage of myself, and all music performed is original. As such, there are no
real ethical concerns with this particular piece. However, as this is also a
call for artists to consider training dynamics as a key component of their
craft, we urge them to be thoughtful and respectful in their use of machine
learning models, data, and in the spirit of their expressiveness.

## References

  * Brown et al. [2020] Tom Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared D Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, Sandhini Agarwal, Ariel Herbert-Voss, Gretchen Krueger, Tom Henighan, Rewon Child, Aditya Ramesh, Daniel Ziegler, Jeffrey Wu, Clemens Winter, Chris Hesse, Mark Chen, Eric Sigler, Mateusz Litwin, Scott Gray, Benjamin Chess, Jack Clark, Christopher Berner, Sam McCandlish, Alec Radford, Ilya Sutskever, and Dario Amodei.  Language models are few-shot learners.  In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, _Advances in Neural Information Processing Systems_ , volume 33, pages 1877–1901. Curran Associates, Inc., 2020.  URL https://proceedings.neurips.cc/paper/2020/file/1457c0d6bfcb4967418bfb8ac142f64a-Paper.pdf.
  * Castro [2019] Pablo Samuel Castro.  Performing structured improvisations with pre-trained deep learning models.  In _Proceedings of the International Conference on Computational Creativity (ICCC)_ , 2019.
  * Castro [2020] Pablo Samuel Castro.  Ganterpretations.  In _4th Workshop on Machine Learning for Creativity and Design at NeurIPS 2020_ , 2020.
  * Dinculescu et al. [2019] Monica Dinculescu, Jesse Engel, and Adam Roberts, editors.  _MidiMe: Personalizing a MusicVAE model with user data_ , 2019.
  * Esser et al. [2021] Patrick Esser, Robin Rombach, and Björn Ommer.  Taming transformers for high-resolution image synthesis.  In _Proceedings of CVPR’21_ , 2021.
  * Huang et al. [2019] Cheng-Zhi Anna Huang, Ashish Vaswani, Jakob Uszkoreit, Ian Simon, Curtis Hawthorne, Noam Shazeer, Andrew M. Dai, Matthew D. Hoffman, Monica Dinculescu, and Douglas Eck.  Music transformer.  In _International Conference on Learning Representations_ , 2019.
  * White [2019] Tom White.  Shared visual abstractions.  _CoRR_ , abs/1912.04217, 2019.  URL http://arxiv.org/abs/1912.04217.

##  Appendix A Screenshots

Screenshots of a performance.

Figure 1: Screenshots of parts 1-4 (from left to right, top to bottom) of a
performance. True functions are in white, learned functions are in red.

##  Appendix B Code details

All the code was custom-built using SuperCollider and OpenFrameworks. Both
pieces of software communicate with each other via OSC. The code will be made
available on GitHub.

###  B.1 SuperCollider

The MIDI I/O between the computer and the Disklavier was done with code
written in SuperCollider. It is generally difficult to detune an acoustic
piano in real-time, so in addition to the regular acoustic sound of the piano,
I superimpose synthesized sounds generated in SuperCollider, which I can
detune using the loss provided by the OpenFrameworks software.

###  B.2 OpenFrameworks

The machine learning and visualization code is written in OpenFrameworks.
Because of the simplicity of the functions being approximated, I calculate the
gradients and losses “by hand” instead of using a machine learning framework.

Distortion 2 is achieved by projecting the “texture” of the webcam onto each
bubble, and passing each of these bubbles through custom shaders (written in
GLSL) that can alter the color of each pixel and placement of each vertex.


# Causal State Estimation and Heisenberg Uncertainty Principle

Junxin Chen  Benjamin B. Lane  LIGO, Massachusetts Institute of Technology,
Cambridge, MA 02139, USA  Su Direkci  The Division of Physics, Mathematics and
Astronomy, California Institute of Technology, CA 91125, USA  Dhruva Ganapathy
Xinghui Yin  Nergis Mavalvala  LIGO, Massachusetts Institute of Technology,
Cambridge, MA 02139, USA  Yanbei Chen  The Division of Physics, Mathematics
and Astronomy, California Institute of Technology, CA 91125, USA  Vivishek
Sudhir<EMAIL_ADDRESS>LIGO, Massachusetts Institute of Technology,
Cambridge, MA 02139, USA  Department of Mechanical Engineering, Massachusetts
Institute of Technology, Cambridge, MA 02139, USA

###### Abstract

The observables of a noisy quantum system can be estimated by appropriately
filtering the records of their continuous measurement. Such filtering is
relevant for state estimation and measurement-based quantum feedback control.
It is therefore imperative that the observables estimated through a causal
filter satisfy the Heisenberg uncertainty principle. In the Markovian setting,
prior work implicitly guarantees this requirement. We show that any causal
estimate of linear observables of a linear, but not necessarily Markovian,
system will satisfy the uncertainty principle. In particular, this is true
irrespective of any feedback control of the system and of where in the
feedback loop — inside or outside — the measurement record is accessed.
Indeed, causal estimators using the in-loop measurement record can be as
precise as those using the out-of-loop record. These results clarify the role
of causal estimators to a large class of quantum systems, restores the
equanimity of in-loop and out-of-loop measurements in their estimation and
control, and simplifies future experiments on measurement-based quantum
feedback control.

_Introduction._ Estimating the state of a system from a noisy measurement
record is ubiquitous in engineering and science. For a classical linear system
driven by stationary Gaussian noise, the optimal estimator — in the sense of
minimizing the mean-square error — is given by the Wiener filter [1, 2]. For
prediction problems where only the past measurement record is available, the
Wiener filter is causal and is given in terms of the spectrum of the
measurement record and a model of its cross-spectrum with the state.
Additionally, if the dynamics are Markovian, a state-space model of the system
and its observation is available and the Kalman filter uses that description
to produce an equivalent, more tractable estimate of the state [3]. State
estimation is also a crucial element in optimal control: the separation
principle [4, 5, 6] asserts that the feedback controller that minimizes a
quadratic cost function of a linear system driven by Gaussian noise can be
split into an optimal state estimator and a linear regulator.

These ideas have been transposed to quantum systems: the theory of quantum
state estimation [7, 8, 9, 10, 11, 12], feedback control [13, 14, 15, 16, 17],
and the separation principle [18] have all been developed. A central
distinction between classical and quantum systems is that the state of the
latter cannot be determined with arbitrary precision, even in principle. The
Heisenberg uncertainty principle imposes a limit on the precision of
accessible knowledge. Consequently, it is believed that any causal estimator
of a quantum state will not violate the uncertainty principle.

The purpose of this paper is to prove that this intuitive expectation is true.
Recent experiments in quantum state estimation and feedback control [19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34] call for a simple yet
general guarantee that causal state estimation and control of quantum systems
will not violate the basic tenet of the Heisenberg uncertainty principle. We
show that causal quantum state estimation using an open-loop measurement
record satisfies the uncertainty principle irrespective of whether the system
dynamics are Markovian or not. If the measurement is used for continuous
feedback, the in-loop measurement record can have features like noise
squashing [35]; thus, an additional out-of-loop signal is often used as an
independent check of the prepared state in measurement-based quantum control
experiments [29, 34, 33]. We show that the in-loop record can be used for
state estimation, just as the out-of-loop record, without sacrificing
precision. Thus, our work guarantees that causal quantum state estimation —
with or without feedback, whether the system is Markovian or not — will
produce a physical state; further, it demonstrates the underestimated value of
the in-loop record in quantum feedback experiments.

_Causal estimation without feedback._ Let us first recall the premise of the
Heisenberg uncertainty principle. The positivity of the quantum state
$\hat{\rho}$ implies that $\Tr[\hat{M}(t)^{\dagger}\hat{M}(t)\hat{\rho}]\geq
0$ for any operator $\hat{M}(t)$ (thought of as being in the Heisenberg
picture). Choosing $\hat{M}=\hat{A}+\lambda\hat{B}$, the inequality must still
hold for the minimum value of its left-hand side as a function of $\lambda$.
This gives the uncertainty principle [36]

$\begin{split}\sigma_{A(t)}^{2}\sigma_{B(t)}^{2}\geq&\frac{1}{4}\absolutevalue{\expectationvalue{\anticommutator{\hat{A}(t)}{\hat{B}(t)}}}^{2}\\\
&+\frac{1}{4}\absolutevalue{\expectationvalue{\commutator{\hat{A}(t)}{\hat{B}(t)}}}^{2}.\end{split}$
(1)

where $\sigma_{A}^{2}\equiv\langle(\hat{A}-\langle\hat{A}\rangle)^{2}\rangle$
is the variance of the operator $\hat{A}$, and $\commutator*{\cdot}{\cdot}$ is
the commutator and $\anticommutator*{\cdot}{\cdot}$ the anti-commutator. That
is, (i) as long as the quantum state is positive, the uncertainty principle
holds, and (ii) if either the anti-/commutator is proportional to the identity
operator, the lower bound is state-independent.

Positivity of the state is implicit in the Markovian setting [7, 8, 9, 11,
37]. The state is driven by Markovian noises that describe the measurement
(and any other environment the system interacts with); when averaged over the
realizations of these noises, the state satisfies a Lindblad-type master
equation, which preserves positivity. However, experiments are susceptible to
non-Markovian noises, so the results of prior work need not apply.

We consider the problem without assuming Markovianity at the expense of
assuming linearity. By linearity, we mean that the observables of interest
have commutators proportional to their identity at all times [38, 39] and that
the estimator is linear in the measurement. In this setting, the uncertainty
principle written in the form
$\sigma_{A(t)}^{2}\sigma_{B(t)}^{2}\geq\tfrac{1}{4}|\langle[\hat{A}(t),\hat{B}(t)]\rangle|^{2}$
is state independent, and its sanctity boils down to the equal-time
commutators of the relevant observables.

Concretely, we focus attention on a bosonic system specified by a pair of
canonically conjugate variables $(\hat{x}(t),\hat{p}(t))$. The linearity
condition is then equivalent to asserting that the observables of interest are
linear combinations of the canonical pair and that the system Hamiltonian is
at most quadratic in the pair. The system is monitored continuously, the
result of which is described by the measurement record $\hat{y}$ — itself a
linear combination of the canonical pair. The continuous monitoring condition
is that [38]

$[\hat{y}(t),\hat{y}(t^{\prime})]=0\quad\mathrm{for\;all}\quad t,t^{\prime}.$
(2)

The observables are estimated by filters $H_{x,p}$ acting linearly on the past
record:

$\hat{x}_{e}(t)=(H_{x}*\hat{y})(t),\quad\hat{p}_{e}(t)=(H_{p}*\hat{y})(t),$
(3)

where $*$ stands for convolution. The estimation errors, i.e. the uncertainty
in the obervables of the state conditioned on the measurement record [11], are

$\displaystyle\Delta\hat{x}(t)$
$\displaystyle=\hat{x}(t)-\hat{x}_{e}(t)=\hat{x}(t)-(H_{x}*\hat{y})(t)$ (4)
$\displaystyle\Delta\hat{p}(t)$
$\displaystyle=\hat{p}(t)-\hat{p}_{e}(t)=\hat{p}(t)-(H_{p}*\hat{y})(t).$ (5)

If the uncertainty principle is to be upheld, the equal-time commutator of the
errors should be

$[\Delta\hat{x}(t),\Delta\hat{p}(t)]\stackrel{{\scriptstyle?}}{{=}}[\hat{x}(t),\hat{p}(t)]=i\hbar.$
(6)

Using eq. 4 together with eqs. 2 and 6

$\displaystyle[\Delta\hat{x}(t),\Delta\hat{p}(t)]$
$\displaystyle=i\hbar-\commutator{\hat{x}(t)}{\hat{p}_{e}(t)}-\commutator{\hat{x}_{e}(t)}{\hat{p}(t)}.$
(7)

The key in determining the remaining equal-time commutators is an
understanding of the _unequal-time_ commutators between the system observable
and the measurement record. For this, we resort to linear response theory [40,
38, 39, 41]. The key result is Kubo’s formula, which gives the linear response
of a system observable (not necessarily linear) $\hat{A}(t)$ to an external
perturbation $\hat{F}(t)$, whose effect is described by the Hamiltonian
$\hat{H}(t)=\hat{H}_{0}-\hat{B}\hat{F}(t)$, where $\hat{H}_{0}$ is the free
Hamiltonian and $\hat{B}$ is a system observable. Kubo’s formula states

$\hat{A}(t)=\hat{A}^{(0)}(t)+\frac{i}{\hbar}\int_{-\infty}^{t}\chi_{AB}(t,t^{\prime})\hat{F}(t^{\prime})\differential{t^{\prime}}$
(8)

where $\hat{A}^{(0)}(t)$ is the evolution under $\hat{H}_{0}$ and

$\chi_{AB}(t,t^{\prime})=\expectationvalue{\commutator{\hat{A}^{(0)}(t)}{\hat{B}^{(0)}(t^{\prime})}}$
(9)

is the linear response susceptibility, evaluated on the state the system is in
before the perturbation $\hat{F}$ acts. For a linear system with linear
observables $\hat{A}$ and $\hat{B}$, it can be shown (see Appendix) that the
susceptibility is

$\chi_{AB}(t,t^{\prime})=\commutator{\hat{A}^{(0)}(t)}{\hat{B}^{(0)}(t^{\prime})}.$
(10)

Figure 1: (a) Schematic diagram for a general open-loop measurement. $\hat{F}$
is a general force coupling to the system through system parameter $\hat{x}$.
The effect of this interaction is encoded in the outgoing field
$\hat{y}_{\mathrm{out}}$ together with other system information. A general
force $\hat{G}$ couples to the system through $\hat{y}_{\mathrm{out}}$ does
not influence the system due to causality. (b) Schematic diagram for a
measurement-based feedback controlled system. $H$ and $K$ stand for the system
and controller transfer functions respectively. $r$ and $y_{\mathrm{c}}$ are
the set point and the in-loop measurement record respectively.

Consider now the schematic diagram of a general measurement as shown in fig. 1
(a) where $\hat{y}_{\mathrm{out}}=\hat{y}$ is the outgoing measurement record
and $\hat{G}$ is the generalized force coupling to the system through
$\hat{y}_{\mathrm{out}}$ by the interaction Hamiltonian
$-\hat{y}_{\mathrm{out}}\hat{G}$. Focusing on $\hat{x}$ and using eq. 8, we
have

$\hat{x}(t)=\hat{x}^{(0)}(t)+\frac{i}{\hbar}\int_{-\infty}^{t}\chi_{xy}(t,t^{\prime})\hat{G}(t^{\prime})\differential{t^{\prime}},$
(11)

with
$\chi_{xy}(t,t^{\prime})=[\hat{x}^{(0)}(t),\hat{y}_{\mathrm{out}}^{(0)}(t^{\prime})]$.
Now we bring in the crucial ingredient of causality: for an open-loop system,
the measurement record $\hat{y}$ cannot influence the system dynamics at a
later time [42], i.e.

$\chi_{xy}(t,t^{\prime})=[\hat{x}^{(0)}(t),\hat{y}^{(0)}(t^{\prime})]=0\mbox{\quad
for\quad}t\geq t^{\prime}.$ (12)

Since the susceptibility is zero for $t\geq t^{\prime}$, the force $\hat{G}$
does not influence the system for such times. Thus

$\chi_{xy}(t,t^{\prime})=[\hat{x}(t),\hat{y}(t^{\prime})]=0\mbox{\quad
for\quad}t\geq t^{\prime}.$ (13)

Let us now consider $[\hat{x}(t),\hat{p}_{e}(t)]$: using eq. 3,

$[\hat{x}(t),\hat{p}_{e}(t)]=\int_{-\infty}^{+\infty}H_{p}(\tau)[\hat{x}(t),\hat{y}(t-\tau)]\differential{\tau}.$
(14)

By causality of the estimation filter, $H_{p}(\tau)=0$ for $\tau<0$. For
$\tau\geq 0$, eq. 13 implies that
$[\hat{x}(t),\hat{y}(t-\tau)]=\chi_{xy}(t,t-\tau)=0$. Thus causality — of the
measurement interaction and of the estimation filter — implies that
$[\hat{x}(t),\hat{p}_{e}(t)]=0$. Similarly, it can be shown that
$[\hat{x}_{e}(t),\hat{p}(t)]=0$.

In summary, for causal estimation using open-loop measurements,

$[\Delta\hat{x}(t),\Delta\hat{p}(t)]=i\hbar.$ (15)

The Heisenberg uncertainty principle holds for linear causal estimates based
on linear open-loop measurements of a linear system.

_Case with feedback._ We now consider the case where the system is feedback-
controlled using the in-loop measurement record $\hat{y}_{c}$ (see fig. 1(b)).
In this case the in-loop record does affect the system at later times after
the measurement interaction. However the out-of-loop record $\hat{y}_{o}$ does
not. In order to use this fact, note that the in-loop and out-of-loop
measurement records, are related to each other as

$\hat{y}_{c}(t)=(K_{c}*\hat{y}_{\mathrm{o}})(t),$ (16)

where $K_{c}(t)$ is the inverse Fourier transform of
$(1+H[\omega]K[\omega])^{-1}$ ($H[\omega]$ and $K[\omega]$ are the transfer
functions of the system and the controller respectively). Importantly, if
$H,K$ are causal, then so is $K_{c}$.

Imagine now estimation based on the in-loop measurement record $\hat{y}_{c}$:
$\hat{x}_{\mathrm{ec}}(t)=(H_{x}^{\mathrm{c}}*\hat{y}_{\mathrm{c}})(t)$ and
$\hat{p}_{\mathrm{ec}}(t)=(H_{p}^{\mathrm{c}}*\hat{y}_{\mathrm{c}})(t)$.
Following the same line of reasoning as in the open-loop case, the question of
whether the estimated observables respect the Heisenberg uncertainty principle
boils down to whether $[\hat{x}(t),\hat{y}_{c}(t^{\prime})]=0$ for $t\geq
t^{\prime}$. Clearly,

$\displaystyle[\hat{x}(t),\hat{y}_{c}(t^{\prime})]=\int_{-\infty}^{\infty}K_{c}(t^{\prime}-\tau)[\hat{x}(t),\hat{y}_{\mathrm{o}}(\tau)]\differential{\tau}.$

We consider the integrand in two complementary intervals, $t^{\prime}<\tau$
and $t^{\prime}\geq\tau$. In the former region,
$K_{\mathrm{c}}(t^{\prime}-\tau)=0$ since $K_{\mathrm{c}}$ is causal. In the
latter region, since $t\geq t^{\prime}\geq\tau$, we have that
$[\hat{x}(t),\hat{y}_{\mathrm{o}}(\tau)]=0$. In sum, the integral is zero in
both intervals due to causality. Thus $[\hat{x}(t),\hat{y}_{c}(t^{\prime})]$,
vanishes as long as $t\geq t^{\prime}$. A similar argument holds for
$[\hat{p}(t),\hat{y}_{c}(t^{\prime})]$. Thus
$[\hat{x}(t),\hat{y}_{c}(t^{\prime})]=0=[\hat{p}(t),\hat{y}_{c}(t^{\prime})]$
for $t\geq t^{\prime}$. In summary, the estimation errors
$\Delta\hat{x}_{c}=\hat{x}-\hat{x}_{\mathrm{ec}},\Delta\hat{p}_{c}=\hat{p}-\hat{p}_{\mathrm{ec}}$
satisfy

$[\Delta\hat{x}_{\mathrm{c}}(t),\Delta\hat{p}_{\mathrm{c}}(t)]=i\hbar,$ (17)

just as in the open-loop case.

_Effect of feedback on estimation error._ The in-loop measurement record
$\hat{y}_{\mathrm{c}}$ is often deemed untrustworthy for state estimation.
Historically, this view originated from the apparent violation of the
Heisenberg uncertainty principle by the in-loop field in linear measurement-
based feedback control of optical fields [43]. This behavior, called “noise
squashing” [35], arises because fields inside a feedback loop are not freely
propagating and so need not satisfy canonical commutation relations [44, 45].

In fact, if a model of the feedback loop is available, then state estimation
using the in-loop measurement record is as accurate as an estimate based on
the out-of-loop record. We can always write the system observable in the
presence of feedback, $\hat{x}_{\mathrm{c}}$, as the sum of the open-loop one
$\hat{x}_{\mathrm{o}}$ plus a feedback term:

$\hat{x}_{\mathrm{c}}(t)=\hat{x}_{\mathrm{o}}(t)-\hat{x}_{\mathrm{fb}}(t)=\hat{x}_{\mathrm{o}}(t)-(K*\hat{y}_{\mathrm{c}})(t),$
(18)

where $K$ is a causal feedback control filter. Now consider two filters
$H_{x}^{i}[\omega]$ $(i=\mathrm{c,o})$, one estimating the observable
$\hat{x}$ from the in-loop record and the other from the out-of-loop record.
The respective errors are
$\Delta\hat{x}_{i}(t)=\hat{x}_{i}(t)-(H_{x}^{i}*\hat{y}_{i})(t)$. It is
straightforward to show that if the in-loop filter is chosen to be

$H_{x}^{\mathrm{c}}[\omega]=\frac{H_{x}^{\mathrm{o}}[\omega]}{K_{\mathrm{c}}[\omega]}-K[\omega]=(H_{x}^{\mathrm{o}}[\omega]H[\omega]-1)K[\omega]+H_{x}^{\mathrm{o}}[\omega],$

then $\Delta\hat{x}_{\mathrm{c}}=\Delta\hat{x}_{\mathrm{o}}$. Note that if the
feedback loop and $H_{x}^{\mathrm{o}}$ are stable and causal, so is
$H_{x}^{\mathrm{c}}$. Therefore by proper filter design, state estimation can
be performed using the in-loop measurement record without loss of fidelity.

_Conclusion._ We have produced a straightforward argument that shows that
causal quantum state estimation respects the Heisenberg uncertainty principle
for the estimated observables. This is done without assuming the system or
measurement dynamics are Markovian and applies equally to systems with and
without feedback. We assume the quantum system to be linear and utilize linear
response theory, both widely applicable for current experimental measurement-
based quantum control. Furthermore in the scenario of feedback control, we
clarify that despite “noise squashing”, the in-loop measurement record can
provide as faithful an estimate of an observable as the out-of-loop record. As
a matter of practice, this vastly simplifies experiments and extends the reach
of quantum state estimation to non-Markovian scenarios [46]. Even more
important, in principle, this insight eliminates the compromise in measurement
efficiency that is required in having two simultaneous measurements on the
same system.

_Acknowledgements._ We wish to thank C. Meng and W. Bowen for fruitful
discussions. We acknowledge the generous support from the National Science
Foundation (NSF). Especially, JC and BBL are supported through the NSF grant
PHY-20122088, DG by the NSF grant PHY-0823459, and SD and YC by NSF grants
PHY-2011961, PHY-2011968, and PHY–1836809. In addition SD and YC acknowledge
the support by the Simons Foundation (Award Number 568762).

##  Appendix A Susceptibility of linear observables of linear systems

The general expression for the susceptibility $\chi_{AB}(t)$ is given in
linear response theory as [40, Theorem 1]

$\begin{split}\chi_{AB}(t)&=\mathrm{Tr}\\{[\hat{\rho},\hat{A}(0)]\hat{B}(t)\\}=\expectationvalue{[\hat{A}(0),\hat{B}(t)]}.\end{split}$
(19)

where $\hat{\rho}$ is the density operator and $\hat{A},\hat{B}$ are arbitrary
operators. The second equality follows from using the cyclic property of the
trace.

In the following, we will show that for linear observables of a linear system,
the susceptibility $\chi_{AB}(t)$ is equal to the second expression above
_without_ the expectation value bracket.

First, let us clarify the definitions of linear systems and linear observables
[39]. A linear quantum system is any quantum system whose Hamiltonian is at
most quadratic in its canonical coordinates and momenta. By a linear
observable of such a system, we mean any linear combination, not necessarily
time-independent, of the system’s canonical coordinates, momenta, and identity
operator.

That is, if the canonical coordinates and momenta are collected as
$\\{\hat{C}_{i}\\}_{i=1,2,\cdots}$, with

$[\hat{C}_{i},\hat{C}_{j}]=c_{ij}\mathds{1},$ (20)

where $\\{c_{ij}\\}$ are scalars, a linear system is described by the
Hamiltonian

$\hat{H}(t)=h_{ij}^{(2)}(t)\hat{C}_{i}\hat{C}_{j}+h_{i}^{(1)}(t)\hat{C}_{i}+h^{(0)}(t)\mathds{1},$
(21)

where $h_{ij}^{(k)}(t)$ are some functions. A generic linear observable is
then

$\hat{A}(t)=\alpha_{i}^{(A)}(t)\hat{C}_{i}+\beta^{(A)}(t)\mathds{1};$ (22)

here $\\{\alpha_{i}^{(A)}(t),\beta^{(A)}(t)\\}$ are functions.

Given a linear observable $\hat{B}(0)$ in terms of its constant coefficients
$\\{\alpha_{i}^{(B)}(0),\beta^{(B)}(0)\\}$, its time evolution is given by the
Heisenberg equation $\dot{\hat{B}}(t)=(i/\hbar)[\hat{H}(t),\hat{B}(t)]$. On
the one hand,
$\dot{\hat{B}}(t)=\dot{\alpha}_{i}^{(B)}(t)\hat{C}_{i}+\dot{\beta}^{(B)}\mathds{1}$,
while on the other (using eqs. 21 and 20),

$\displaystyle[\hat{H}(t),\hat{B}(t)]=$
$\displaystyle\left(h_{ij}^{(2)}(t)+h_{ji}^{(2)}(t)\right)c_{jk}\alpha_{k}^{(B)}(t)\,\hat{C}_{i}$
$\displaystyle\qquad+h_{i}^{(1)}(t)c_{ij}\alpha_{k}^{(B)}(t)\,\mathds{1}.$

Comparing the two expressions shows that: (i) the ansatz eq. 22 is sufficient
for evolution under $\hat{H}$ for time $t>0$ if it holds at $t=0$, and (ii)
the coefficients $\\{\alpha_{i}^{(B)}(t),\beta^{(B)}(t)\\}$ are given by
solutions of the differential equations:

$\displaystyle\dot{\alpha}_{i}^{(B)}$
$\displaystyle=\left(h_{ij}^{(2)}(t)+h_{ji}^{(2)}(t)\right)c_{jk}\alpha_{k}^{(B)}\equiv
M_{ik}(t)\alpha_{k}^{(B)}$ $\displaystyle\dot{\beta}^{(B)}$
$\displaystyle=h_{i}^{(1)}(t)c_{ij}\alpha_{k}^{(B)}(t).$

These have solutions as long as the entries of the matrix $M_{ik}(t)$ are
continuous functions of time[47, Chap. 4] — which is assumed.

Having secured that the ansatz in eq. 22 is consistent, a simple calculation
shows that

$[\hat{A}(0),\hat{B}(t)]=\alpha_{i}^{(A)}(0)c_{ij}\alpha_{j}^{(B)}(t)\,\mathds{1}.$
(23)

Since this is proportional to the identity operator, the expectation value
$\expectationvalue*{[\hat{A}(0),\hat{B}(t)]}=\alpha_{i}^{(A)}(0)c_{ij}\alpha_{j}^{(B)}(t)$
is state-independent. So we may define $\chi_{AB}(t)=[\hat{A}(0),\hat{B}(t)]$
.

## References

  * Wiener [1942] N. Wiener,  _Extrapolation, Interpolation, and Smoothing of Stationary Time Series_ (MIT Press, 1942).
  * Kailath [1981] T. Kailath,  _Lectures on Wiener and Kalman Filtering_ (Springer, 1981).
  * Kalman [1960] R. E. Kalman, A New Approach to Linear Filtering and Prediction Problems, Journal of Basic Engineering 82, 35 (1960).
  * Kushner [1967] H. J. Kushner,  _Stochastic Stability and Control_ (Academic Press, 1967).
  * Wonham [1968] W. M. Wonham, On the Separation Theorem of Stochastic Control, SIAM Journal on Control 6, 312 (1968).
  * Lindquist [1973] A. Lindquist, On Feedback Control of Linear Stochastic Systems, SIAM Journal on Control 11, 323 (1973).
  * Belavkin [1989] V. P. Belavkin, A new wave equation for a continuous nondemolition measurement, Physics Letters A 140, 355 (1989).
  * Barchielli and Belavkin [1991] A. Barchielli and V. Belavkin, Measurements continuous in time and a posteriori states in quantum mechanics, Journal of Physics A 24, 1495 (1991).
  * Belavkin [1992] V. P. Belavkin, Quantum Continual Measurements and a Posteriori Collapse on CCR, Communications in Mathematical Physics 146, 611 (1992).
  * Wiseman [1996] H. M. Wiseman, Quantum trajectories and quantum measurement theory, Quantum and Semiclassical Optics: Journal of the European Optical Society Part B 8, 205 (1996).
  * Doherty  _et al._ [1999] A. C. Doherty, S. M. Tan, A. S. Parkins, and D. F. Walls, State determination in continuous measurement, Phys. Rev. A 60, 2380 (1999).
  * Bouten  _et al._ [2004] L. Bouten, M. Guta, and H. Maassen, Stochastic Schrödinger equations, Journal of Physics A: Mathematical and General 37, 3189 (2004).
  * Wiseman [1994] H. M. Wiseman, Quantum theory of continuous feedback, Physical Review A 49, 2133 (1994).
  * Wiseman and Doherty [2005] H. M. Wiseman and A. C. Doherty, Optimal Unravellings for Feedback Control in Linear Quantum Systems, Physical Review Letters 94, 070405 (2005).
  * Genoni  _et al._ [2013] M. G. Genoni, S. Mancini, and A. Serafini, Optimal feedback control of linear quantum systems in the presence of thermal noise, Physical Review A 87, 042333 (2013).
  * Wiseman and Milburn [2010] H. M. Wiseman and G. J. Milburn,  _Quantum Measurement and Control_ (Cambridge University Press, 2010).
  * Jacobs [2014] K. Jacobs,  _Quantum Measurement Theory and its Applications_ (Cambridge University Press, 2014).
  * Bouten and Handel [2008] L. Bouten and R. V. Handel, On the separation principle in quantum control, in  _Quantum Stochastics and Information_ (World Scientific, 2008) pp. 206–238.
  * Hood  _et al._ [2000] C. J. Hood, T. W. Lynn, A. C. Doherty, A. S. Parkins, and H. J. Kimble, The Atom-Cavity Microscope: Single Atoms Bound in Orbit by Single Photons, Science 287, 1447 (2000).
  * Kubanek  _et al._ [2009] A. Kubanek, M. Koch, C. Sames, A. Ourjoumtsev, P. W. H. Pinkse, K. Murr, and G. Rempe, Photon-by-photon feedback control of a single-atom trajectory, Nature 462, 898 (2009).
  * Sayrin  _et al._ [2011] C. Sayrin, I. Dotsenko, X. Zhou, B. Peaudecerf, T. Rybarczyk, S. Gleyzes, P. Rouchon, M. Mirrahimi, H. Amini, M. Brune, J.-M. Raimond, and S. Haroche, Real-time quantum feedback prepares and stabilizes photon number states, Nature 477, 73 (2011).
  * Zhou  _et al._ [2012] X. Zhou, I. Dotsenko, B. Peaudecerf, T. Rybarczyk, C. Sayrin, S. Gleyzes, J. M. Raimond, M. Brune, and S. Haroche, Field Locked to a Fock State by Quantum Feedback with Single Photon Corrections, Physical Review Letters 108, 243602 (2012).
  * Vijay  _et al._ [2012] R. Vijay, C. Macklin, D. H. Slichter, S. J. Weber, K. W. Murch, R. Naik, A. N. Korotkov, and I. Siddiqi, Stabilizing Rabi oscillations in a superconducting qubit using quantum feedback, Nature 490, 77 (2012).
  * Hatridge  _et al._ [2013] M. Hatridge, S. Shankar, M. Mirrahimi, F. Schackert, K. Geerlings, T. Brecht, K. M. Sliwa, B. Abdo, L. Frunzio, S. M. Girvin, R. J. Schoelkopf, and M. H. Devoret, Quantum Back-Action of an Individual Variable-Strength Measurement, Science 339, 178 (2013).
  * Murch  _et al._ [2013] K. Murch, S. Weber, C. Macklin, and I. Siddiqi, Observing single quantum trajectories of a superconducting quantum bit, Nature 502, 211 (2013).
  * Weber  _et al._ [2014] S. Weber, A. Chantasri, J. Dressel, A. N. Jordan, K. Murch, and I. Siddiqi, Mapping the optimal route between two quantum states, Nature 511, 570 (2014).
  * Wilson  _et al._ [2015] D. J. Wilson, V. Sudhir, N. Piro, R. Schilling, A. Ghadimi, and T. J. Kippenberg, Measurement-based control of a mechanical oscillator at its thermal decoherence rate, Nature 524, 325 (2015).
  * Hacohen-Gourgy  _et al._ [2016] S. Hacohen-Gourgy, L. S. Martin, E. Flurin, V. V. Ramasesh, K. B. Whaley, and I. Siddiqi, Quantum dynamics of simultaneously measured non-commuting observables, Nature 538, 491 (2016).
  * Sudhir  _et al._ [2017] V. Sudhir, D. Wilson, R. Schilling, H. Schütz, S. Fedorov, A. Ghadimi, A. Nunnenkamp, and T. Kippenberg, Appearance and Disappearance of Quantum Correlations in Measurement-Based Feedback Control of a Mechanical Oscillator, Physical Review X 7, 011001 (2017).
  * Rossi  _et al._ [2018] M. Rossi, D. Mason, J. Chen, Y. Tsaturyan, and A. Schliesser, Measurement-based quantum control of mechanical motion, Nature 563, 53 (2018).
  * Rossi  _et al._ [2019] M. Rossi, D. Mason, J. Chen, and A. Schliesser, Observing and verifying the quantum trajectory of a mechanical resonator, Phys. Rev. Lett. 123, 163601 (2019).
  * Minev  _et al._ [2019] Z. K. Minev, S. O. Mundhada, S. Shankar, P. Reinhold, R. Gutiérrez-Jáuregui, R. J. Schoelkopf, M. Mirrahimi, H. J. Carmichael, and M. H. Devoret, To catch and reverse a quantum jump mid-flight, Nature 570, 200 (2019).
  * Tebbenjohanns  _et al._ [2021] F. Tebbenjohanns, M. L. Mattana, M. Rossi, M. Frimmer, and L. Novotny, Quantum control of a nanoparticle optically levitated in cryogenic free space, Nature 595, 378 (2021).
  * Magrini  _et al._ [2021] L. Magrini, P. Rosenzweig, C. Bach, A. Deutschmann-Olek, S. G. Hofer, S. Hong, N. Kiesel, A. Kugi, and M. Aspelmeyer, Real-time optimal quantum control of mechanical motion at room temperature, Nature 595, 373 (2021).
  * Wiseman [1999] H. Wiseman, Squashed states of light: theory and applications to quantum spectroscopy, Journal of Optics B: Quantum and Semiclassical Optics 1, 459 (1999).
  * Robertson [1929] H. P. Robertson, The uncertainty principle, Phys. Rev. 34, 163 (1929).
  * Tsang  _et al._ [2009] M. Tsang, J. H. Shapiro, and S. Lloyd, Quantum theory of optical temporal phase and instantaneous frequency. ii. continuous-time limit and state-variable approach to phase-locked loop design, Phys. Rev. A 79, 053843 (2009).
  * Braginsky and Khalili [1992] V. B. Braginsky and F. Khalili,  _Quantum Measurement_ (Cambridge University Press, 1992).
  * Buonanno and Chen [2002] A. Buonanno and Y. Chen, Signal recycled laser-interferometer gravitational-wave detectors as optical springs, Physical Review D 65, 042001 (2002).
  * Kubo [1957] R. Kubo, Statistical-mechanical theory of irreversible processes. i. general theory and simple applications to magnetic and conduction problems, Journal of the Physical Society of Japan 12, 570 (1957).
  * Clerk [2004] A. A. Clerk, Quantum-limited position detection and amplification: A linear response perspective, Physical Review B 70, 245306 (2004).
  * Gardiner and Collett [1985] C. W. Gardiner and M. J. Collett, Input and output in damped quantum systems: Quantum stochastic differential equations and the master equation, Physical Review A 31, 3761 (1985).
  * Yamamoto  _et al._ [1986] Y. Yamamoto, N. Imoto, and S. Machida, Amplitude squeezing in a semiconductor laser using quantum nondemolition measurement and negative feedback, Physical Review A 33, 3243 (1986).
  * Shapiro  _et al._ [1987] J. Shapiro, G. Saplakoglu, S.-T. Ho, P. Kumar, B. Saleh, and M. Teich, Theory of light detection in the presence of feedback, JOSA B 4, 1604 (1987).
  * Taubman  _et al._ [1995] M. S. Taubman, H. Wiseman, D. E. McClelland, and H.-A. Bachor, Intensity feedback effects on quantum-limited noise, Journal of the Optical Society of America B 12, 1792 (1995).
  * Whittle  _et al._ [2021] C. Whittle, E. D. Hall, S. Dwyer, N. Mavalvala, V. Sudhir, _et al._ , Approaching the motional ground state of a 10-kg object, Science 372, 1333 (2021).
  * Arnold [1992] V. I. Arnold,  _Ordinary Differential Equations_ , 3rd ed. (Springer, 1992).