# Lotka-Volterra population dynamics in coherent and tunable oscillators of
trapped polariton condensates

J. D. Töpfer<EMAIL_ADDRESS>School of Physics and Astronomy,
University of Southampton, Southampton, SO171BJ, United Kingdom  H. Sigurdsson
School of Physics and Astronomy, University of Southampton, Southampton,
SO171BJ, United Kingdom  Skolkovo Institute of Science and Technology, Moscow,
Bolshoy Boulevard 30, bld. 1, 121205, Russia  S. Alyatkin  Skolkovo Institute
of Science and Technology, Moscow, Bolshoy Boulevard 30, bld. 1, 121205,
Russia  P. G. Lagoudakis<EMAIL_ADDRESS>School of Physics and
Astronomy, University of Southampton, Southampton, SO171BJ, United Kingdom
Skolkovo Institute of Science and Technology, Moscow, Bolshoy Boulevard 30,
bld. 1, 121205, Russia

(August 27, 2024)

We demonstrate a regime in which matter-wave condensates of exciton-polaritons
trapped in an elliptically shaped two-dimensional potential appear as a
coherent mixture of ground and first-excited state of the quantum harmonic
oscillator. This system resembles an optically controllable two-level system
and produces near terahertz harmonic oscillations of the condensate’s center
of mass along the major axis of the elliptical trapping potential. The
population ratio between the two trap levels is tunable through the excitation
laser power and is shown to follow Lotka-Volterra dynamics. We demonstrate
coherence formation between two spatially displaced trapped condensate
oscillators - the polaritonic analogue of Huygen’s clock synchronization for
coupled condensate oscillators.

###### pacs:

##  I Introduction

The competitive Lotka–Volterra model describes coupled autonomous systems (or
species in terms of predator-prey settings) which fight over a common
resource. Such competition can be found in brain activity Rabinovich  _et al._
(2008), economics Goodwin (1982), and laser systems described by Maxwell-Bloch
equations where different modes compete over the laser gain Hacinliyan  _et
al._ (2010). However, dissipative Bose-Einstein condensates can also display
similar competative population dynamics given the right setting. The purest
form of condensation in the quantum sense refers to an equilibrium bosonic gas
macroscopically occupying the system ground state while other levels deplete
Mueller  _et al._ (2006). When non-Hermiticity is added to the picture, which
is the case for photon and polariton condensates, this macroscopic occupation
is no longer exclusive to just the ground state but instead other available
higher energy levels might become macroscopically occupied, depending on their
internal losses and coupling to other levels leading to competition over the
condensate gain.

While many-particle two-level systems such as lasers and condensates can give
great insight to fundamental physics of coupled systems, they also form the
solid-state implementation of the qubit, like superconducting circuits, for
quantum information processing. Indeed, the possibility of using Bose-Einstein
condensates for quantum computation was recently proposed Byrnes  _et al._
(2012), followed by some exciting work focused on exciton-polariton
condensates due to their easy optical manipulation Demirchyan  _et al._
(2014); Cuevas  _et al._ (2018); Ghosh and Liew (2020); Xue  _et al._ (2019).
Moreover, non-Hermitian two-level systems of polariton condensates could
potentially model the open Bose-Hubbard dimer Lledó  _et al._ (2019),
$\mathcal{PT}$-symmetric physics Bender  _et al._ (2007); Chestnov  _et al._
(2016), and nonlinear Josephson effects such as macroscopic self-trapping
Ostrovskaya  _et al._ (2000). Recently, they have been theoretically proposed
as clock generators under resonant laser excitation Leblanc  _et al._ (2020),
and to realize propagating parity domain walls in extended condensates
Sigurdsson  _et al._ (2015). It is therefore of quite some interest to be able
to generate and control macroscopically occupied dissipative two-level systems
to investigate the aforementioned possibilities. Such control also opens new
perspectives to study synchronization in extended two-level bosonic systems
where adjacent traps possess different level occupation.

Our study is then motivated by the following question: Can a two-level
dissipative condensate system be designed with tunable energy spacing and
competitive population dynamics which determine the ratio of particles between
the excited state and the ground state? In other words, can the condensate
Bloch vector be tuned to have arbitrary polar angle on the Bloch sphere and
precession frequency? Here, we answer this question in the affirmative. We
show experimentally and theoretically that a non-resonantly excited exciton-
polariton condensate, possessing two energy levels from its laser induced trap
Askitopoulos  _et al._ (2013), displays a gradual change in the relative level
occupation number as a function of excitation power while the laser geometry
determines the level splitting. The condensate population dynamics are found
to follow exactly the competitive Lotka-Volterra equations due to cross-
saturation effects coupling the levels together. Our modeling provides a
transparent description of the physics at play and strengthen the connection
between dissipative polariton condensates and laser systems. We additionally
demonstrate the scalability of our system by measuring the coherent coupling
between two two-level condensates by using two distinct laser beams. There, we
observe that the dissipative nature of their tunneling current synchronizes
the systems to maximizes their mutual gain.

##  II Results

The dissipative nature of polariton condensates requires continuous pumping to
compensate for (mainly optical) losses and maintain the condensate beyond its
intrinsic $\sim\mathrm{ps}$ lifetime. Optical excitation of polariton
condensates by means of a non-resonant laser beam is interlinked with the
generation of an incoherent reservoir of excitons, which acts as both, optical
gain and a repulsive potential for polaritons. The spatial distribution of
this particle interaction-induced potential is inherently linked to the
spatial geometry of the excitation pump profile. Design of potential
landscapes by means of optical pump laser beam-shaping is at the heart of
research for polaritonic devices and applications. While the use of widely-
spaced tightly-focused pump laser beams (spot size $\approx 2\;\mathrm{\upmu
m}$) leads to the excitation of high-energy and radially expanding condensates
Ohadi  _et al._ (2016); Töpfer  _et al._ (2020), on the contrary, the use of
closely-spaced pump beam geometries can show condensation in low-energy and
spatially confined polariton modes Tosi  _et al._ (2012); Cristofolini  _et
al._ (2013); Ohadi  _et al._ (2015). In particular, it has been shown that a
non-resonant circular or elliptical annular pump profile creates a near-
parabolic potential landscape for polaritons and facilitates trapped polariton
condensation in states approximately given by the eigenstates $\Psi_{nm}$ of
the two-dimensional quantum harmonic oscillator (HO) Askitopoulos  _et al._
(2015, 2018). Change in experimental parameters, such as the diameter or the
pump power of the annular excitation profile, modifies the system’s resonance
conditions and allows for condensation in different HO states characterized by
their indices $n$ and $m$.

In the following, we demonstrate polariton condensation in an elliptically-
shaped potential favoring condensation in the lowest two HO modes. Nonlinear
coupling between the two energetically split modes (cross-saturation effects)
is shown to yield a limit cycle state. This non-stationary state involves
linear harmonic oscillations of the condensates center of mass $\braket{x}$ in
analogy to the harmonic oscillations of a mechanical pendulum for small angles
of displacement (see Fig. 1).

Figure 1: Schematic illustration of harmonic oscillations of a polariton
condensate in an optically generated elliptical trapping potential. Dual-mode
condensation into the two energetically lowest modes $\Psi_{00}$ and
$\Psi_{10}$ with energy splitting $\hbar\Omega$ leads to an oscillatory motion
of the condensate’s center of mass $\braket{x}\propto\cos{(\Omega t)}$ in
analogy to the harmonic motion of a mechanical pendulum for small displacement
angles.

###  II.1 Competition of two condensate modes

By shaping the non-resonant continuous wave excitation laser profile into an
elliptical annulus (see Methods) with diameters of $14.2\;\mathrm{\upmu m}$
and $10.6\;\mathrm{\upmu m}$ in $x$\- and $y$-direction, respectively, we lift
the degeneracy of the two excited HO states $\Psi_{10}$ and $\Psi_{01}$
Askitopoulos  _et al._ (2018). Here, the trap’s elongation in x-direction
leads to a larger net gain for the $\Psi_{10}$ state as for the $\Psi_{01}$
state, thus, favoring condensation into the former state. In Fig. 2(a) we show
the measured real-space polariton photoluminescence (PL) for the system
excited with pump power $P=1.1P_{\mathrm{thr}}$, where $P_{\mathrm{thr}}$ is
the system’s condensation threshold pump power. We confirm single-mode
condensation in the $\Psi_{10}$ state by spectrally resolving the emission
along the major axis of the elliptical trap ($y=0$) as illustrated in Fig.
2(d). The observed blueshift of $\approx 1.1\;\mathrm{meV}$ above the lower
polariton ground state energy is a result of the optically induced background
(confinement) potential. By increasing the excitation pump power to
$P=2.1P_{\mathrm{thr}}$ we observe a transition of the system from single-mode
excited state $\Psi_{10}$ occupation to single-mode ground state $\Psi_{00}$
occupation as shown in Figs. 2(c) and 2(f). Interestingly, in the transition
region, i.e. $1.1P_{\mathrm{thr}}<P<2.1P_{\mathrm{thr}}$, we observe the
coexistence of both modes. Such dual-mode operation with approximately equal
occupation of both modes is realized at $P=1.3P_{\mathrm{thr}}$ and shown in
Figs. 2(b) and 2(e). The continuous transition of excited state to ground
state condensation with increasing pump power $P$ is illustrated in Figs. 2(g)
and 2(h) showing the (spatially integrated) spectra and normalised mode
emission intensities (blue and red curves), respectively. We point out that
the mode linewidths for the spectra illustrated in Figs. 2(d-g) are limited by
the spectral resolution $\Delta E\approx 0.02\;\mathrm{meV}$ of our
experimental setup. We separately estimate a linewidth of $\approx
0.004\;\mathrm{meV}$ for the ground state condensate at
$P=2.1P_{\mathrm{thr}}$ from analysis of the condensate’s auto-correlation
function (see Appendix B).

While the energy levels of both modes $E_{00}$ and $E_{10}$ experience a
continuous blueshift with increasing excitation pump power
$P>1.1P_{\mathrm{thr}}$ we find that - within our experimental resolution -
there is no change in energy splitting $\hbar\Omega=E_{10}-E_{00}\approx
0.22\;\mathrm{meV}$ (see black line in Fig. 2(h)). Because the dominant
contribution to the energy splitting $\hbar\Omega$ between ground and first
excited state is given by the optically induced confinement potential, we
conclude that there is no significant change in the potential landscape for
increasing pump power $P$ above condensation threshold.

Figure 2: Two-mode competition in a trapped polariton condensate. (a-c) Real
space and (d-f) energy-resolved real space photoluminescence along the major
axis ($y=0$) of the elliptical annular pump beam (red-dashed line) for varying
excitation pump power $P$ above condensation threshold. (g) Integrated spectra
for increasing excitation pump power $P$. (h) Spectral weights of the ground-
state mode (red) and the first excited state (blue) and their energy splitting
$\hbar\Omega$ (black). The experimental spectral resolution
($0.02\;\mathrm{meV}$) is illustrated in gray-shaded. Energy scale in (d-g) is
scaled as the blueshift with respect to the lower polariton ground state mode
below threshold. (i) Relative mode occupations $N_{0,1}/N_{\mathrm{tot}}$ with
$N_{\mathrm{tot}}=(N_{0}+N_{1})$ for two weakly coupled modes described in
Eqs. (3a) and (3b) as a function of their relative net gain $p_{0}/p_{1}$.

Motivated by our experimental results we describe the condensate wave function
$\Psi(\mathbf{r},t)$ in a two-mode model

$\Psi(\mathbf{r},t)=\psi_{0}(t)\Psi_{00}(\mathbf{r})+\psi_{1}(t)\Psi_{10}(\mathbf{r}),$
(1)

where $\Psi_{00,10}$ represent the orthogonal and normalized spatial field
distributions of the system’s ground- and first-excited state. The complex-
valued envelopes $\psi_{0,1}$ resemble occupation numbers
$N_{0,1}=|\psi_{0,1}|^{2}$ and phases $\theta_{0,1}=\arg(\psi_{0,1})$ of each
mode, respectively. One can derive the nonlinear-coupled two-level system (see
Appendix A)

$\displaystyle i\dot{\psi}_{0}$
$\displaystyle=\left[\omega_{0}+ip_{0}+(\alpha_{0}-i\sigma_{0})|\psi_{0}|^{2}+(\beta-i\kappa)|\psi_{1}|^{2}\right]\psi_{0},$
(2) $\displaystyle i\dot{\psi}_{1}$
$\displaystyle=\left[\omega_{1}+ip_{1}+(\alpha_{1}-i\sigma_{1})|\psi_{1}|^{2}+(\beta-i\kappa)|\psi_{0}|^{2}\right]\psi_{1}.$

Here, for each mode $j=0,1$ the real-valued parameters denote eigenfrequency
$\omega_{j}$, net gain $p_{j}$, self- and cross- nonlinearities $\alpha_{j}$
and $\beta$, as well as the self- and cross-saturation terms $\sigma_{j}$ and
$\kappa$. Equations (2) can analogously be rewritten into four coupled
differential equations for the occupation numbers $N_{0,1}$ and phases
$\theta_{0,1}$ of each mode, i.e.

$\displaystyle\dot{N}_{0}$
$\displaystyle=2\cdot\left(p_{0}-\sigma_{0}N_{0}-\kappa N_{1}\right)N_{0},$
(3a) $\displaystyle\dot{N}_{1}$
$\displaystyle=2\cdot\left(p_{1}-\sigma_{1}N_{1}-\kappa N_{0}\right)N_{1},$
(3b) $\displaystyle\dot{\theta}_{0}$
$\displaystyle=\omega_{0}+\alpha_{0}N_{0}+\beta N_{1},$ (3c)
$\displaystyle\dot{\theta}_{1}$
$\displaystyle=\omega_{1}+\alpha_{1}N_{1}+\beta N_{0}.$ (3d)

Equations (3a) and (3b), which are decoupled from the phase dynamics, have
previously been theoretically studied in the context of mode competition in
polariton condensates Eastham (2008); Pukrop and Schumacher (2020). They
represent the competitive Lotka-Volterra equations, and as such describe
competition between two interacting species populations Murray (2007), as well
as the competition between two laser modes Sargent III  _et al._ (1974). The
relative strength of competition between the two modes sharing the same gain
medium is given by the parameter $C=\kappa^{2}/\sigma_{1}\sigma_{2}$. While in
the so-called strong coupling regime ($C>1$) cross-saturation effects dominate
the system and one species will quench the population of the second entity, in
the following, we restrict our analysis to the experimentally relevant weak
coupling regime ($C<1$) which facilitates the coexistence of both species.
Considering configurations where at least one mode is pumped above threshold,
i.e. $p_{0,1}>0$ , there exist three non-trivial equilibrium points
($\dot{N}_{0,1}=0$) given by

$\displaystyle(\mathrm{I})$ $\displaystyle N_{0}=0,$ $\displaystyle
N_{1}=p_{1}/\sigma_{1},$ (4) $\displaystyle(\mathrm{II})$ $\displaystyle
N_{0}=\frac{\sigma_{1}p_{0}-\kappa p_{1}}{\sigma_{0}\sigma_{1}-\kappa^{2}},$
$\displaystyle N_{1}=\frac{\sigma_{0}p_{1}-\kappa
p_{0}}{\sigma_{0}\sigma_{1}-\kappa^{2}},$ $\displaystyle(\mathrm{III})$
$\displaystyle N_{0}=p_{0}/\sigma_{0},$ $\displaystyle N_{1}=0.$

However, only one of these points can be stable for the same parameters and
stability of each solution $(\mathrm{I}-\mathrm{III})$ is determined by the
system’s Jacobian which yields the conditions

$\displaystyle(\mathrm{I})$ $\displaystyle p_{0}<p_{1}\kappa/\sigma_{1},$
$\displaystyle p_{1}>0,$ (5) $\displaystyle(\mathrm{II})$ $\displaystyle
p_{0}>p_{1}\kappa/\sigma_{1},$ $\displaystyle p_{1}>p_{0}\kappa/\sigma_{0},$
$\displaystyle(\mathrm{III})$ $\displaystyle p_{0}>0,$ $\displaystyle
p_{1}<p_{0}\kappa/\sigma_{0}.$

While two of these points $(\mathrm{I},\mathrm{III})$ represent single-mode
polariton emission, one point $(\mathrm{II})$ corresponds to dual-mode
emission with energies $E_{0,1}=\hbar\dot{\theta}_{0,1}$ given by Eqs. (3c)
and (3d). In Fig. 2(i) we show the equilibrium occupation numbers $N_{0}$ and
$N_{1}$ for increasing relative net gain $p_{0}/p_{1}$ between the two modes.
The continuous transition from single-mode emission in the second mode (I) to
single-mode emission in the first mode (III) is interleaved with a region of
dual-mode emission (II). The physical reasons for the increase of
$p_{0}/p_{1}$ in experiment as a function of laser power can be understood
from two different mechanisms: First, as pump power is increased further above
threshold the number of particles in the condensate increases which blueshifts
the excited state out of the trap making it more and more lossy. Second, the
growing number of particles in the reservoir enables polaritons to relax more
efficiently in energy Wouters and Carusotto (2010); Askitopoulos  _et al._
(2015) which makes $p_{0}$ grow faster than $p_{1}$.

We calculate the transient behavior of the competitive Lotka-Volterra system
towards its fixed point attractors for all three scenarios (I-III) as shown in
Figs. 3(a-e). Here, for each configuration, we also plot the system’s Lyapunov
function $\mathcal{L}(N_{0},N_{1})$ in false-colorscale which can be written
in quadratic form MacArthur (1970),

$\mathcal{L}=-2p_{0}N_{0}-2p_{1}N_{1}+\sigma_{0}N_{0}^{2}+\sigma_{1}N_{1}^{2}+2\kappa
N_{0}N_{1}.$ (6)

We point out that the Lyapunov function $\mathcal{L}$ satisfies the condition
$\dot{\mathcal{L}}\leq 0$ and that the system’s fixed point attractors
correspond to the minima of $\mathcal{L}(N_{0},N_{1})$ within its domain
$\left\\{N_{0}\geq 0,N_{1}\geq 0\right\\}$. We stress that the observed change
in mode populations here is very different from the mode switching previously
reported in Sun  _et al._ (2018) which was attributed to linear
(deterministic) polariton physics where the laser parameters changed which
quantum mode had the lowest condensation threshold.

The presented configurations with stable mode occupation numbers $N_{0,1}$ are
solutions of the nonlinear coupled two-level system
$\boldsymbol{\psi}(t)=(\psi_{1}(t),\psi_{0}(t))^{T}$ described in Eqs. (2).
Hereby, the constant emission energy $E_{0,1}=\hbar\dot{\theta}_{0,1}$ of each
mode is blue-shifted and determined in Eqs. (3c) and (3d). These two-level
states $\boldsymbol{\psi}(t)$ describe orbits on the Bloch sphere spanned by
the three-component vector
$\mathbf{S}=(\boldsymbol{\psi}^{\dagger}\boldsymbol{\sigma}\boldsymbol{\psi})/(\boldsymbol{\psi}^{\dagger}\boldsymbol{\psi})$
with Pauli matrix vector written $\boldsymbol{\sigma}$. One can show that

$\displaystyle S_{x}(t)$
$\displaystyle=\frac{2\sqrt{N_{1}N_{0}}}{N_{1}+N_{0}}\cos(\Omega t+\delta),$
(7) $\displaystyle S_{y}(t)$
$\displaystyle=\frac{2\sqrt{N_{1}N_{0}}}{N_{1}+N_{0}}\sin(\Omega t+\delta),$
$\displaystyle S_{z}(t)$ $\displaystyle=\frac{N_{1}-N_{0}}{N_{1}+N_{0}}.$

where $\delta$ is a stochastic phase-offset chosen during the spontaneous
condensation event. In Figs. 3(f-j) we illustrate the orbits corresponding to
the intensity equilibrium points depicted in Figs. 3(a-e). Condensate
configurations with single-mode occupation in ground or excited state yield
fixed points located at the poles of the Bloch-sphere as shown in (f) and (j).
The periodic orbits in (g-i) for dual-mode operation are centered around the
$S_{z}$-axis and their period is determined by the modes’ energy splitting,
i.e. $T=2\pi\Omega^{-1}$. Modification of the relative occupation numbers (as
demonstrated in Fig. 2(h)) allows control over the position of the orbit’s
$S_{z}$-position. Furthermore, as we demonstrate in Appendix C, the period $T$
of these orbits is optically tuneable in the range of $10-20\;\mathrm{ps}$
through spatial modulation of the pump laser profile.

Figure 3: Numerical analysis of two weakly coupled condensate modes with
occupation numbers $N_{0}$ and $N_{1}$ described in Eqs. (3a) and 3b. (a-e)
Phase portraits and Lyapunov potentials $\mathcal{L}$ (false-colorscale, Eq.
(6)) for increasing relative net gain $p_{0}/p_{1}$ from left to right. Black
circles depict fixed point attractors for vanishing occupation $N_{0}=0$ of
the ground state (a), dual-mode operation $N_{0,1}>0$ (b-d) and vanishing
occupation $N_{1}=0$ of the excited state (e). (f-j) Representation of the
intensity equilibrium points on the Bloch sphere spanned by the three-
component vector $\mathbf{S}=(S_{x},S_{y},S_{z})$ described in Eq. (7). (k-o)
Projection of the orbits on the Bloch sphere (red) onto the two-dimensional
phase-space (blue) spanned by the condensate’s center of mass $\langle
x\rangle$ and momentum expectation value $\langle p\rangle$ illustrates the
regimes in which the trapped condensate is stationary (k,o) and harmonically
oscillating (l-n). Numerical parameters:
$\sigma_{0}=\sigma_{1}=1\;\mathrm{ps}^{-1}$,
$p_{1}=\kappa=0.5\;\mathrm{ps}^{-1}$ and
$p_{0}=(0.25\;\mathrm{ps}^{-1},0.31\;\mathrm{ps}^{-1},0.5\;\mathrm{ps}^{-1},0.81\;\mathrm{ps}^{-1},1\;\mathrm{ps}^{-1})$
for (a-e) and (f-j), respectively.

###  II.2 Harmonic oscillations of a trapped condensate

There is a direct relation between the periodic orbits appearing on the two-
level system’s Bloch-sphere described in Eq. (7) and the harmonic oscillations
of the condensate density $|\Psi(\mathbf{r},t)|^{2}$ illustrated in Fig. 1. In
particular, the expectation values of the condensate’s center of mass along
the horizontal (major) axis in real space $\braket{x}$ and momentum space
$\braket{p}$ are given by

$\displaystyle\braket{x}(t)$ $\displaystyle=x_{0}\cdot S_{x}(t),$ (8a)
$\displaystyle\braket{p}(t)$ $\displaystyle=p_{0}\cdot S_{y}(t),$ (8b)

where the scaling factors $x_{0}=\bra{\Psi_{00}}x\ket{\Psi_{10}}$ and
$p_{0}=\bra{\Psi_{00}}\hbar\partial_{x}\ket{\Psi_{10}}$ are dependent on the
spatial distribution of the two competing modes. The trajectories of the two-
level system on the Bloch sphere in Figs. 3(f-j) can be projected onto the
phase-plane spanned by space $\braket{x}$ and momentum $\braket{p}$
expectation values as illustrated in Figs. 3(k-o). Here, single-mode
condensation in either excited state (k) or ground state (o) both correspond
to a stationary condensate center of mass located at
$\braket{x}=\braket{p}=0$. The periodic orbits appearing in Figs. 3(l-n)
depict harmonic oscillations of the condensate’s space $\braket{x}$ and
momentum $\braket{p}$ expectation values in analogy to the phase-plane orbits
of a mechanical undamped oscillator.

Although we cannot directly assess the condensate’s ultrafast (ps-timescale)
oscillations, we indirectly resolve them through time-correlation
measurements. Hereby, we interfere the real-space polariton photoluminescence
$\Psi(\mathbf{r},t)$ with a retro-reflected and time-shifted version
$\Psi(-\mathbf{r},t+\tau)$ of itself and extract the system’s first-order
correlation function (see Methods)

$g^{(1)}(-\mathbf{r},\mathbf{r};\tau)=\frac{\langle\Psi^{*}(\mathbf{r},t)\Psi(-\mathbf{r},t+\tau)\rangle}{\sqrt{\langle|\Psi(\mathbf{r},t)|^{2}\rangle\langle|\Psi(-\mathbf{r},t)|^{2}\rangle}},$
(9)

where $\langle...\rangle$ denotes time averaging. The coherence function (Eq.
(9)) is a complex-valued and normalized measure for the first-order
correlations between the two signals $\Psi(\mathbf{r},t)$ and
$\Psi(-\mathbf{r},t+\tau)$. When the signals are composed of two frequency
components with frequency-detuning $\Omega$ we expect periodic modulation of
the coherence function $g^{(1)}(-\mathbf{r},\mathbf{r};\tau)$ at the same
frequency $\Omega$ according to the Wiener-Khinchin theorem.

In the following we unravel the time-correlations of the trapped condensate
system presented in Fig. 2(b,e) with approximately equal occupation of ground
and first excited state. Recorded interference patterns and extracted modulus
$|g^{(1)}|$ and argument $\arg(g^{(1)})$ of the correlation function are
illustrated in Figs. 4(a,c,e) for no relative time-delay $\tau=0$ between the
two signals. The equal-time correlations reveal vanishing coherence
($|g^{(1)}(-\mathbf{r},\mathbf{r};\tau=0)$) along the vertical lines
$x\approx\pm 1.9\;\mathrm{\upmu m}$. This is a direct result of the mixing of
two states with opposite parity and equal contribution to the condensate
density, $N_{0}|\Psi_{00}(\mathbf{r})|^{2}=N_{1}|\Psi_{10}(\mathbf{r})|^{2}$.
On the contrary, at small distances $|x|\ll 1.9\;\mathrm{\upmu m}$, where the
condensate density is mainly formed by particles in the ground state, we find
a highly correlated signal $g^{(1)}\approx 1$ resembling the symmetric nature
(even parity) of the ground state $\Psi_{00}$. At larger distances $|x|\gg
1.9\;\mathrm{\upmu m}$ we find an anti-correlated signal $g^{(1)}\approx-1$ as
a result of the anti-symmetric field distribution (odd parity) of the excited
state $\Psi_{10}$. The composition of even and odd parity states is a
necessary condition for a non-vanishing oscillation amplitude
$x_{0}=\bra{\Psi_{00}}x\ket{\Psi_{10}}$ of the condensate’s center of mass.

Introducing a temporal delay of half an oscillation period
$\tau=\pi/\Omega\approx 10\;\mathrm{ps}$ we measure an interference pattern as
illustrated in Fig. 4(b) with corresponding modulus and phase of the
correlation function $g^{(1)}(-\mathbf{r},\mathbf{r};\tau=10\;\mathrm{ps})$
shown in Figs. 4(d) and (f). Here, we have applied a constant offset to the
argument of the correlation function to yield 0 phase at the origin
$\mathbf{r}=0$, i.e. we measure from a co-rotating reference frame. We find
large correlation between spatially mirrored positions
$g^{(1)}(-\mathbf{r},\mathbf{r};\tau=10\;\mathrm{ps})\approx 1$ throughout the
whole system in agreement with the condensate’s spatial density oscillation.
Further, in Figs. 4(g) and (h) we plot the extracted modulus and argument of
the correlation function $g^{(1)}(-x,x;\tau)$ along the trap major axis $y=0$
for varying temporal delay $\tau$. Periodic disappearance and revival of
coherence at $x=\pm 1.9\;\mathrm{\upmu m}$ (see Fig. 4(i)) with frequency
$\Omega/2\pi\approx 52\;\mathrm{GHz}$ reveal the coherent nature of the
present harmonic oscillations persisting over more than 20 oscillations ($\sim
400\;\mathrm{ps}$).

Figure 4: First-order correlations of a dual-mode polariton condensate in an
elliptical trap. (a,b) Recorded interference of the real-space polariton
photoluminescence $\Psi(\mathbf{r},t)$ and its retro-reflected and time-
shifted version $\Psi(-\mathbf{r},t+\tau)$ for (a) $\tau=0$ and (b)
$\tau=10\;\mathrm{ps}$. The annular pump profile is shown as a red-dashed
ellipses. (c,d) Magnitude and (e,f) phase of the digitally reconstructed
first-order coherence function $g^{(1)}(-\mathbf{r},\mathbf{r};\tau)$ for
(c,e) $\tau=0$ and (d,f) $\tau=10\;\mathrm{ps}$, respectively. Temporal
evolution of (g) magnitude and (h) phase of the coherence function
$g^{(1)}(-x,x;\tau)$ along the horizontal slice $y=0$. For better visibility
an offset to the phase pattern has been applied for each time step $\tau$ to
yield 0 phase at the center $r=0$. (i) Harmonic oscillation of the coherence
function $|g^{(1)}(-x,x;\tau)|$ with maximum oscillation-amplitude at
$x=1.9\;\mathrm{\upmu m}$. The condensate is pumped at $1.3$ times the
condensation threshold. Scale bar in (a) corresponds to $4\;\mathrm{\upmu m}$
and applies to (a-f).

###  II.3 Synchronization of two polariton oscillators

In a next step we investigate a system of two spatially separated polariton
condensates, each of them representing a trapped (dual-mode) oscillator. As we
schematically illustrate in Fig. 5(a) for small separation distances between
the two traps a finite coupling rate $\mathcal{V}_{0,1}$ between the system’s
ground states $\Psi_{00}^{(1,2)}$ and excited states $\Psi_{10}^{(1,2)}$,
where the superscript denotes different traps, leads to phase-locking of the
condensates’ periodic orbits and, thus, to synchronized oscillatory motion of
the centers of mass $\braket{x^{(1,2)}}$ of both condensates.

Figure 5: Synchronization of two polariton oscillators. (a) Schematic showing
the coupling between two trapped and spatially separated dual-mode condensates
and the resulting synchronization of their periodic orbits on the Bloch
sphere. (b) Measured real-space condensate photoluminescence excited by two
elliptical annular pump profiles (major axis $\approx 14.7\;\mathrm{\upmu m}$,
minor axis $\approx 10.4\;\mathrm{\upmu m}$) shown as red-dashed lines. (c)
Corresponding energy-resolved emission along symmetry axis $y=0$. (d) Modulus
of the first-order coherence function $|g^{(1)}(-x,x;\tau)|$ showing
correlations between the centers of the two trapped condensates. (e) Periodic
disappearance and revival of coherence $|g^{(1)}(\tau)|$ at
$x=7\;\mathrm{\upmu m}$. (f) Energy-resolved momentum space emission along
$k_{y}=0$. White dashed line depicts the lower polariton dispersion below
condensation threshold. (g) Extracted momentum-space mode profiles for
interfering ground states $\Psi^{(1,2)}_{00}$ (red) and excited states
$\Psi^{(1,2)}_{10}$ (blue) at energies indicated by arrows in (g). Black
dashed lines display the calculated mode profiles in the absence of
interference, i.e. for incoherent modes.

For our experiment we generate two identical elliptical annular pump profiles
with similar dimensions as those presented in Figs. 2 and 4. The two trap
centers are displaced by $d\approx 20.7\;\mathrm{\upmu m}$ and excitation at
pump power $P\approx 1.3P_{\mathrm{thr}}$ leads to approximately equal
occupation of ground state $\Psi^{(1,2)}_{00}$ and excited state
$\Psi^{(1,2)}_{10}$ in each trap $(1)$ and $(2)$, respectively. In Fig. 5(b)
we show the measured real-space condensate PL and in Fig. 5(c) we plot the
spectrally-resolved emission along the symmetry axis $y=0$. We note that -
within the spectral resolution of our experiment - the energies of the ground
and first excited modes are identical in both traps, where each trap has
energy splitting $\hbar\Omega\approx 0.20\;\mathrm{meV}$. Coherence between
the two spatially separated dual-mode condensates is confirmed by
interferometric cross-correlation measurements. In particular, the measured
modulus of the first-order correlation function $|g^{(1)}(-x,x;\tau)|$ along
the symmetry axis $y=0$ is shown in Fig. 5(d). The illustrated spatial range
$5\;\mathrm{\upmu m}<x<15\;\mathrm{\upmu m}$ contains the mutual coherence
between both traps centered at $x\approx\pm 10\;\mathrm{\upmu m}$. In analogy
to our results for a single oscillating condensate displayed in Fig. 4(g) we
observe spatially-distributed time-periodic modulation of the coherence
function demonstrating synchronized spatial oscillations of both trapped
condensates. A profile of the temporal coherence function with maximum
oscillation amplitude at $x=7\;\mathrm{\upmu m}$ is displayed in Fig. 5(e)
showing periodic disappearance and revival of coherence with frequency $\sim
50\;\mathrm{GHz}$. We point out that the measured coherence function displayed
in Fig. 5(e) with vanishing correlation at $\tau=0$ reveals phase-locking of
the two oscillating condensate center of masses $\braket{x^{(1,2)}}$ with
vanishing phase difference (in-phase synchronization).

Coherence between the two condensate traps and phase-locking of their
respective oscillating centers of mass $\braket{x^{(1,2)}}$ is further
confirmed by our measurement of the system’s energy-resolved momentum space
emission shown in Fig. 5(f). Here, we resolve (far-field) interference for
both condensate modes with extracted intensity profiles shown in Fig. 5(g) for
ground state (red) and first excited state (blue). For comparison we plot the
calculated un-modulated profiles (black-dashed) that resemble incoherent
emission between both condensate traps. The experimentally observed far-field
intensity modulation with destructive interference at $k_{x}=0$ reveals anti-
phase synchronization between the two spatially separated condensates at both
energy levels, i.e. $\Delta\theta_{0}=\theta^{(2)}_{0}-\theta^{(1)}_{0}=\pi$
and $\Delta\theta_{1}=\theta^{(2)}_{1}-\theta^{(1)}_{1}=\pi$. One can easily
show that phase-locking with equal phase-differences
$\Delta\theta_{0}=\Delta\theta_{1}$ of the coupled HO modes resembles in-phase
synchronization of the two condensate oscillators as schematically illustrated
for the two Bloch spheres in Fig. 5(a).

In the following, we describe the system’s wave function in the tight-binding
approach as a superposition of two spatially separated dual-mode condensates

$\Psi(\mathbf{r},t)=\sum_{j=1,2}\psi^{(j)}_{0}(t)\Psi^{(j)}_{00}(\mathbf{r})+\psi^{(j)}_{1}(t)\Psi^{(j)}_{10}(\mathbf{r}).$
(10)

Synchronization of the two coupled condensate oscillators can be modeled by
introducing a linear coupling $\mathcal{V}_{0,1}$ term between each of the
trapped HO modes. In general this coupling term consists of both dissipative
and non-dissipative parts, i.e. $\mathcal{V}_{0,1}=i\gamma_{0,1}+J_{0,1}$
Aleiner  _et al._ (2012). Using the short form notation
$i\dot{\psi}^{(1,2)}_{0,1}=H^{(1,2)}_{0,1}\psi^{(1,2)}_{0,1}$ where $H_{0,1}$
is a function describing the two-level dynamics of an isolated condensate
oscillator (i.e., square brackets of Eqs. (2)) we can now write the 4 coupled
equations

$\displaystyle i\dot{\psi}^{(1)}_{0}$
$\displaystyle=H^{(1)}_{0}\psi^{(1)}_{0}+\left(i\gamma_{0}+J_{0}\right)\psi^{(2)}_{0},$
(11a) $\displaystyle i\dot{\psi}^{(2)}_{0}$
$\displaystyle=H^{(2)}_{0}\psi^{(2)}_{0}+\left(i\gamma_{0}+J_{0}\right)\psi^{(1)}_{0},$
(11b) $\displaystyle i\dot{\psi}^{(1)}_{1}$
$\displaystyle=H^{(1)}_{1}\psi^{(1)}_{1}+\left(i\gamma_{1}+J_{1}\right)\psi^{(2)}_{1},$
(11c) $\displaystyle i\dot{\psi}^{(2)}_{1}$
$\displaystyle=H^{(2)}_{1}\psi^{(2)}_{1}+\left(i\gamma_{1}+J_{1}\right)\psi^{(1)}_{1},$
(11d)

to describe the dynamics of each of the mode amplitudes $\psi^{(1,2)}_{0,1}$
constituting the total condensate wave function in Eq. (10). It is known that
dissipative coupling with $\gamma_{0,1}<0$ ($\gamma_{0,1}>0$) facilitates
phase-locking of each trapped condensate mode with phase difference
$\Delta\theta_{0,1}=\theta^{(2)}_{0,1}-\theta^{(1)}_{0,1}=\pi(0)$ Aleiner  _et
al._ (2012). Indeed, direct substitution of the in-phase or the anti-phase
ansatz $\psi_{0,1}^{(1)}=\pm\psi_{0,1}^{(2)}$ respectively results in simply
the same set of coupled equations as given by Eqs. (2) with shifted
$\omega_{0,1}\to\omega_{0,1}\pm J_{0,1}$ and $p_{0,1}\to
p_{0,1}\pm\gamma_{0,1}$. A rigorous stability analysis of these and more
exotic solutions is beyond the scope of the current study. While the far-field
emission pattern of our experimentally realized coupled-trap system displayed
in Fig. 5(f) indicates dissipative coupling with $\gamma_{0,1}<0$, other
synchronization patterns are expected when changing the separation distance
between the two traps Harrison  _et al._ (2020).

##  III Discussion

In this work we have demonstrated optical generation, tunability and coupling
of macroscopic two-level systems realized in trapped polariton condensates.
Non-resonant excitation using annular beam profiles enables optical control
over the mode-composition of each trapped condensate. As we have shown, dual-
mode condensation comprising the two energetically-lowest trapped modes
results in nearly terahertz (ps-timescale) spatial density oscillations in
close analogy to the dynamics of a simple mechanical oscillator. Amplitude and
frequency of such a condensate oscillator are both tuneable through modulation
of the excitation laser beam. Dissipative coupling between two closely spaced
condensate traps, each of which represents a (dual-mode) oscillator,
facilitates phase-locking of their spatial oscillations. This effect
represents a quantum analogue of Huygen’s clock synchronization realized by
spatially oscillating coherent many-particle states. Our work demonstrates an
interconnectable element for future polaritonic circuitry, and lays out
directions for polaritonic applications utilizing coupled macroscopic two-
level systems.

##  IV Acknowledgments

The authors acknowledge the support of the Skoltech NGP Program (Skoltech-MIT
joint project), the UK’s Engineering and Physical Sciences Research Council
(grant EP/M025330/1 on Hybrid Polaritonics), and the RFBR projects No.
20-52-12026 (jointly with DFG) and No. 20-02-00919. H.S. acknowledges
hospitality provided by the University of Iceland.

##  V Methods

Experiments are conducted on a strain-compensated semiconductor microcavity
with InGaAs quantum wells held in a cold finger cryostat at a temperature
$T\approx 6\;\mathrm{K}$. The continuous wave Gaussian excitation laser
(circular polarization, $\lambda\approx 785\;\mathrm{nm}$) is modulated using
a spatial light modulator in the Fourier plane of the optical setup and
focused with (elliptical) annular beam shape onto the microcavity using an
$\mathrm{NA}=0.4$ microscope objective lens. The full-width-at-half-maximum of
the laser profile in radial direction is $\approx 1.5\;\mathrm{\upmu m}$. An
acousto-optic modulator transforms the continuous wave signal into square wave
packets at $5\;\%$ duty cycle and $10\;\mathrm{kHz}$ repetition rate to avoid
heating of the sample. We operate at a negative exciton-photon detuning of
$\Delta\approx 5.5\;\mathrm{meV}$. Resulting polariton photoluminescence
($\lambda\approx 857\;\mathrm{nm}$) is collected in reflection geometry,
filtered from the excitation laser using a long-pass filter and imaged
simultaneously onto two different charge-coupled device (CCD) sensors. One
sensor is placed at the output of a $750\;\mathrm{mm}$ spectrometer (equipped
with a $1800\;\mathrm{grooves}/\mathrm{mm}$ grating) such that through
combination of optical lenses we can choose to image spectrally-resolved
momentum-space (far-field) or real-space (near-field) photoluminescence. The
second sensor is placed at the output of a modified Mach-Zehnder
interferometer, in which one mirror is replaced by a retro-reflector mounted
on a translation stage. Here, we chose to image the interference of the real-
space (near-field) polariton photoluminescence with its retro-reflected and
time-shifted distribution, i.e.
$|\Psi(\mathbf{r,t})+\exp{(i\mathbf{k}\mathbf{r})}\Psi(-\mathbf{r},t+\tau)|^{2}$.
The wavevector $\mathbf{k}$ is controlled by transverse displacements of the
retro-reflected beam and determines the orientation and periodicity of the
interference fringes. Recorded interferograms are analyzed by means of off-
axis digital holography to obtain magnitude and phase of the first-order
correlation function $g^{(1)}(\mathbf{r},-\mathbf{r},\tau)$, where $\tau$ is
the temporal shift controlled by the relative path-length difference of the
two interferometer arms.

##  APPENDIX A COUPLED CONDENSATE MODES IN A TRAP

The dynamics of the macroscopic wave function $\Psi(\mathbf{r},t)$ of the
driven-dissipative polariton condensate can be described in mean field
formalism Wouters and Carusotto (2007) as

$i\dot{\Psi}=\left[\frac{-\hbar\nabla^{2}}{2m}+\alpha|\Psi|^{2}-\frac{i\gamma_{c}}{2}+(g+i)n\right]\Psi.$
(12)

Here, intrinsic non-linearity and dissipation of polaritons are described with
parameters $\alpha$ and $\gamma_{c}$. A pump induced background reservoir
$n(\mathbf{r},t)$ provides both gain and a repulsive potential for the
condensate wave function with interaction parameter $g>0$. In the following we
assume that the reservoir can be removed from the dynamics of the condensate
and only the simplest particle saturation mechanism is considered, i.e.
$n(\mathbf{r},t)\approx P(\mathbf{r})(1-R|\Psi(\mathbf{r},t)|^{2})$ where
$R^{-1}$ is the saturation density. Assuming that the modulation to the
reservoir from the condensate is small in the range of investigated pump
powers and corresponding condensate densities, $1\gg
R|\Psi(\mathbf{r},t)|^{2}$, the condensate trapping potential is approximately
defined by the term $gP(\mathbf{r})$ where $P(\mathbf{r})$ denotes the
constant, non-resonant and annular laser intensity profile. The eigenmodes of
the two-dimensional trap satisfy

$\left[\frac{-\hbar\nabla^{2}}{2m}+gP(\mathbf{r})\right]\Psi_{nm}(\mathbf{r})=\omega_{nm}\Psi_{nm}(\mathbf{r}).$
(13)

and are taken to be normalized and real-valued. Since we operate in a regime,
in which only the trap’s ground state $\Psi_{00}$ and first excited state
$\Psi_{10}$ are populated, we project the total wave function $\Psi$ onto a
truncated Hilbert space

$\Psi(\mathbf{r},t)=\psi_{0}(t)\Psi_{00}(\mathbf{r})+\psi_{1}(t)\Psi_{10}(\mathbf{r}),$
(14)

where $\braket{\Psi_{i0}}{\Psi_{j0}}=\delta_{ij}$. Therefore, the problem can
be described as an effective two-level system following the equations of
motion

$\displaystyle i\dot{\psi}_{0,1}=$
$\displaystyle\left[\omega_{0,1}+ip_{0,1}+(\alpha_{0,1}-i\sigma_{0,1})|\psi_{0,1}|^{2}\right.$
(15) $\displaystyle\left.+(\beta-i\kappa)|\psi_{1,0}|^{2}\right]\psi_{0,1}.$

Here, we have neglected nonlinear off-resonant mixing terms between the modes
which vanish in the time average of our continuous-wave experiment. Each mode
$j={1,2}$ is driven by its net gain
$p_{j}=\bra{\Psi_{j0}}P\ket{\Psi_{j0}}-\gamma_{c}/2$, and stabilized by both
its self-saturation term $\sigma_{j}$ and cross-saturation term $\kappa$.
These non-linear damping terms are determined by the overlap integrals

$\displaystyle\sigma_{j}$ $\displaystyle=\int
P(\mathbf{r})\Psi_{j0}^{4}(\mathbf{r})\mathrm{d}\mathbf{r},$ (16a)
$\displaystyle\kappa$ $\displaystyle=2\int
P(\mathbf{r})\Psi_{00}^{2}(\mathbf{r})\Psi_{10}^{2}(\mathbf{r})\mathrm{d}\mathbf{r}.$
(16b)

Nonlinear blue-shift of each mode due to Kerr and cross-Kerr effects are
included in the integrals

$\displaystyle\alpha_{j}$
$\displaystyle=\int\Psi_{j0}^{4}(\mathbf{r})\mathrm{d}\mathbf{r},$ (17a)
$\displaystyle\beta$
$\displaystyle=2\int\Psi_{00}^{2}(\mathbf{r})\Psi_{10}^{2}(\mathbf{r})\mathrm{d}\mathbf{r}.$
(17b)

##  APPENDIX B CONDENSATE LINEWIDTH

Figure 6: First-order coherence function $|g^{(1)}(\tau)|$ of the ground state
$\Psi_{00}$ mode excited at pump power $P=2.1P_{\mathrm{thr}}$ in an
elliptical trap. An exponential fit (red line) of the temporal decay of
coherence yields a $1/e$ coherence time of $t_{\mathrm{c}}=333\;\mathrm{ps}$.
Inset shows the measured interference of the condensates real-space emission
with its retro-reflected version at $\tau=0$ time delay.

We investigate the linewidth of a single-mode trapped polariton condensate
through analysis of the temporal auto-correlation function of its emission. We
choose an elliptical annular pump with major axis $\approx 14.2\;\mathrm{\upmu
m}$ and minor axis $\approx 10.6\;\mathrm{\upmu m}$ to excite a single-mode
ground state $\Psi_{00}$ condensate above condensation threshold at
$P=2.1P_{\mathrm{thr}}$. Real-space and energy-resolved photoluminescence of
this condensate are presented in the main text in Figs. 2(c) and (f). The
measured decay of temporal coherence
$|g^{(1)}(\mathbf{r}=0,\mathbf{r}=0;\tau)|$ extracted at the Gaussian-shaped
condensate center $\mathbf{r}=0$ is shown in Fig. 6. We find a single
exponential decay of coherence with $1/e$ decay time of
$t_{\mathrm{c}}=333\;\mathrm{ps}$. Thus, for a Lorentzian lineshape we
calculate a (full-width-at-half-maximum) condensate linewidth of $\Delta
E=2\hbar/\tau_{\mathrm{c}}\approx 4\;\mathrm{\upmu eV}$, which is below our
experimental spectral resolution $\Delta\approx 0.02\;\mathrm{meV}$.

##  APPENDIX C TUNING OF OSCILLATION FREQUENCY

Figure 7: Measured energy splitting $\hbar\Omega$ between the lowest two modes
$\Psi_{00}$ and $\Psi_{10}$ in an elliptical trap with diameter $d$ of the
major axis. Inset depicts a typical pump laser profile. A power-law decay
$\Omega\propto 1/d$ (red line) is shown as a guide for the eye.

Control over the oscillation frequency of the trapped condensate system is
achieved by tuning of the energy splitting $\hbar\Omega=E_{10}-E_{00}$ between
the two competing HO modes. A change in spatial dimension of the elliptical
annular pump beam profile directly modifies the near-parabolic polaritonic
potential landscape, and, thus, allows to control the HO mode energy
splitting. In Fig. 7 we show the measured energy splitting $\hbar\Omega$
(left) between $\Psi_{00}$ and $\Psi_{10}$ modes and the corresponding
calculated oscillation frequency $(2\pi)^{-1}\Omega$ (right) for varying
diameter $d$ of the elliptical trap pump profile along its (horizontal) major
axis. In our experiment we simultaneously modify the vertical diameter of the
ellipse to guarantee two-mode condensation of ground state $\Psi_{00}$ and
excited state $\Psi_{10}$, i.e. we quench the population of the $\Psi_{01}$
mode. Hereby, the eccentricity $e$ of the elliptical pump profile stays the
range of $\approx 0.6-0.8$ for all points shown in Fig. 7. A typical
excitation laser profile is illustrated in the inset of Fig. 7. We
experimentally demonstrate tunability of the condensate’s spatial oscillations
in the frequency range of $50\;\mathrm{GHz}-100\;\mathrm{GHz}$, which
corresponds to oscillation periods of $20\;\mathrm{ps}-10\;\mathrm{ps}$.

## References

  * Rabinovich  _et al._ (2008) M. I. Rabinovich, R. Huerta, P. Varona, and V. S. Afraimovich, Transient cognitive dynamics, metastability, and decision making, PLOS Computational Biology 4, 1 (2008).
  * Goodwin (1982) R. M. Goodwin, A growth cycle, in  _Essays in economic dynamics_ (Springer, 1982) pp. 165–170.
  * Hacinliyan  _et al._ (2010) A. S. Hacinliyan, I. Kusbeyzi, and O. O. Aybar, Approximate solutions of maxwell bloch equations and possible lotka volterra type behavior, Nonlinear Dynamics 62, 17 (2010).
  * Mueller  _et al._ (2006) E. J. Mueller, T.-L. Ho, M. Ueda, and G. Baym, Fragmentation of bose-einstein condensates, Phys. Rev. A 74, 033612 (2006).
  * Byrnes  _et al._ (2012) T. Byrnes, K. Wen, and Y. Yamamoto, Macroscopic quantum computation using bose-einstein condensates, Phys. Rev. A 85, 040306 (2012).
  * Demirchyan  _et al._ (2014) S. S. Demirchyan, I. Y. Chestnov, A. P. Alodjants, M. M. Glazov, and A. V. Kavokin, Qubits based on polariton rabi oscillators, Phys. Rev. Lett. 112, 196403 (2014).
  * Cuevas  _et al._ (2018) Á. Cuevas, J. C. López Carreño, B. Silva, M. De Giorgi, D. G. Suárez-Forero, C. Sánchez Muñoz, A. Fieramosca, F. Cardano, L. Marrucci, V. Tasco, G. Biasiol, E. del Valle, L. Dominici, D. Ballarini, G. Gigli, P. Mataloni, F. P. Laussy, F. Sciarrino, and D. Sanvitto, First observation of the quantized exciton-polariton field and effect of interactions on a single polariton, Science Advances 4, 10.1126/sciadv.aao6814 (2018).
  * Ghosh and Liew (2020) S. Ghosh and T. C. H. Liew, Quantum computing with exciton-polariton condensates, npj Quantum Information 6, 16 (2020).
  * Xue  _et al._ (2019) Y. Xue, I. Chestnov, E. Sedov, E. Kiktenko, A. Fedorov, S. Schumacher, X. Ma, and A. Kavokin, Split-ring polariton condensates as macroscopic two-level quantum systems (2019), arXiv:1907.00383 [cond-mat.mes-hall] .
  * Lledó  _et al._ (2019) C. Lledó, T. K. Mavrogordatos, and M. H. Szymańska, Driven bose-hubbard dimer under nonlocal dissipation: A bistable time crystal, Phys. Rev. B 100, 054303 (2019).
  * Bender  _et al._ (2007) C. M. Bender, D. C. Brody, H. F. Jones, and B. K. Meister, Faster than hermitian quantum mechanics, Phys. Rev. Lett. 98, 040403 (2007).
  * Chestnov  _et al._ (2016) I. Y. Chestnov, S. S. Demirchyan, A. P. Alodjants, Y. G. Rubo, and A. V. Kavokin, Permanent rabi oscillations in coupled exciton-photon systems with pt -symmetry, Scientific Reports 6, 19551 (2016).
  * Ostrovskaya  _et al._ (2000) E. A. Ostrovskaya, Y. S. Kivshar, M. Lisak, B. Hall, F. Cattani, and D. Anderson, Coupled-mode theory for bose-einstein condensates, Phys. Rev. A 61, 031601 (2000).
  * Leblanc  _et al._ (2020) C. Leblanc, G. Malpuech, and D. D. Solnyshkov, High-frequency exciton-polariton clock generator, Phys. Rev. B 101, 115418 (2020).
  * Sigurdsson  _et al._ (2015) H. Sigurdsson, I. A. Shelykh, and T. C. H. Liew, Switching waves in multilevel incoherently driven polariton condensates, Phys. Rev. B 92, 195409 (2015).
  * Askitopoulos  _et al._ (2013) A. Askitopoulos, H. Ohadi, A. V. Kavokin, Z. Hatzopoulos, P. G. Savvidis, and P. G. Lagoudakis, Polariton condensation in an optically induced two-dimensional potential, Phys. Rev. B 88, 041308 (2013).
  * Ohadi  _et al._ (2016) H. Ohadi, R. L. Gregory, T. Freegarde, Y. G. Rubo, A. V. Kavokin, N. G. Berloff, and P. G. Lagoudakis, Nontrivial phase coupling in polariton multiplets, Phys. Rev. X 6, 031032 (2016).
  * Töpfer  _et al._ (2020) J. D. Töpfer, H. Sigurdsson, L. Pickup, and P. G. Lagoudakis, Time-delay polaritonics, Communications Physics 3, 1 (2020).
  * Tosi  _et al._ (2012) G. Tosi, G. Christmann, N. Berloff, P. Tsotsis, T. Gao, Z. Hatzopoulos, P. Savvidis, and J. Baumberg, Sculpting oscillators with light within a nonlinear quantum fluid, Nature physics 8, 190 (2012).
  * Cristofolini  _et al._ (2013) P. Cristofolini, A. Dreismann, G. Christmann, G. Franchetti, N. G. Berloff, P. Tsotsis, Z. Hatzopoulos, P. G. Savvidis, and J. J. Baumberg, Optical superfluid phase transitions and trapping of polariton condensates, Phys. Rev. Lett. 110, 186403 (2013).
  * Ohadi  _et al._ (2015) H. Ohadi, A. Dreismann, Y. G. Rubo, F. Pinsker, Y. del Valle-Inclan Redondo, S. I. Tsintzos, Z. Hatzopoulos, P. G. Savvidis, and J. J. Baumberg, Spontaneous spin bifurcations and ferromagnetic phase transitions in a spinor exciton-polariton condensate, Phys. Rev. X 5, 031002 (2015).
  * Askitopoulos  _et al._ (2015) A. Askitopoulos, T. C. H. Liew, H. Ohadi, Z. Hatzopoulos, P. G. Savvidis, and P. G. Lagoudakis, Robust platform for engineering pure-quantum-state transitions in polariton condensates, Phys. Rev. B 92, 035305 (2015).
  * Askitopoulos  _et al._ (2018) A. Askitopoulos, A. V. Nalitov, E. S. Sedov, L. Pickup, E. D. Cherotchenko, Z. Hatzopoulos, P. G. Savvidis, A. V. Kavokin, and P. G. Lagoudakis, All-optical quantum fluid spin beam splitter, Phys. Rev. B 97, 235303 (2018).
  * Eastham (2008) P. R. Eastham, Mode locking and mode competition in a nonequilibrium solid-state condensate, Phys. Rev. B 78, 035319 (2008).
  * Pukrop and Schumacher (2020) M. Pukrop and S. Schumacher, Externally controlled lotka-volterra dynamics in a linearly polarized polariton fluid, Phys. Rev. E 101, 012207 (2020).
  * Murray (2007) J. D. Murray,  _Mathematical biology: I. An introduction_ , Vol. 17 (Springer Science & Business Media, 2007).
  * Sargent III  _et al._ (1974) M. Sargent III, M. Scully, and W. Lamb Jr,  _Laser Physics_ (Addison-Wesley Publishing Co., Reading, Massachusetts, 1974).
  * Wouters and Carusotto (2010) M. Wouters and I. Carusotto, Superfluidity and critical velocities in nonequilibrium bose-einstein condensates, Phys. Rev. Lett. 105, 020602 (2010).
  * MacArthur (1970) R. MacArthur, Species packing and competitive equilibrium for many species, Theoretical population biology 1, 1 (1970).
  * Sun  _et al._ (2018) Y. Sun, Y. Yoon, S. Khan, L. Ge, M. Steger, L. N. Pfeiffer, K. West, H. E. Türeci, D. W. Snoke, and K. A. Nelson, Stable switching among high-order modes in polariton condensates, Phys. Rev. B 97, 045303 (2018).
  * Aleiner  _et al._ (2012) I. L. Aleiner, B. L. Altshuler, and Y. G. Rubo, Radiative coupling and weak lasing of exciton-polariton condensates, Phys. Rev. B 85, 121301 (2012).
  * Harrison  _et al._ (2020) S. L. Harrison, H. Sigurdsson, and P. G. Lagoudakis, Synchronization in optically trapped polariton stuart-landau networks, Phys. Rev. B 101, 155402 (2020).
  * Wouters and Carusotto (2007) M. Wouters and I. Carusotto, Excitations in a nonequilibrium bose-einstein condensate of exciton polaritons, Phys. Rev. Lett. 99, 140402 (2007).


# Testing for Overfitting

James Schmidt111Department of Applied Mathematics and Statistics, email:
<EMAIL_ADDRESS>
Johns Hopkins University
Johns Hopkins University Applied Physics Laboratory

###### Abstract

High complexity models are notorious in machine learning for overfitting, a
phenomenon in which models well represent data but fail to generalize an
underlying data generating process. A typical procedure for circumventing
overfitting computes empirical risk on a holdout set and halts once (or flags
that/when) it begins to increase. Such practice often helps in outputting a
well-generalizing model, but justification for why it works is primarily
heuristic.

We discuss the overfitting problem and explain why standard asymptotic and
concentration results do not hold for evaluation with training data. We then
proceed to introduce and argue for a hypothesis test by means of which both
model performance may be evaluated using training data, and overfitting
quantitatively defined and detected. We rely on said concentration bounds
which guarantee that empirical means should, with high probability,
approximate their true mean to conclude that they should approximate each
other. We stipulate conditions under which this test is valid, describe how
the test may be used for identifying overfitting, articulate a further nuance
according to which distributional shift may be flagged, and highlight an
alternative notion of learning which usefully captures generalization in the
absence of uniform PAC guarantees.

##  1 Introduction

Supervised machine learning is severely underdetermined: a finite labeled data
set is used to search a function space for an appropriate model fitting both
the data and “from where the data comes.” While the full function space is
often at least two infinite orders of magnitude greater than the data,
practitioners usually restrict search to a hypothesis class that is
parametrized as a finite dimensional space. If this hypothesis class is too
restricted, the search may output a model which fails to represent or
approximate the data well enough; if, on the other hand, the class is too
rich, the output model may represent the data too well, in that the model
fails to represent the underlying distribution from which data is drawn.
Generally, this tradeoff between underfitting and overfitting, respectively,
is asymmetric: a model which fits data may (and hopefully does) still
generalize to the underlying distribution, while a model which underfits data
usually does not fit the distribution. Stated differently, underfitting is
detectable in the course of performance evaluation while overfitting cannot be
identified by performance on the training data alone ([3]).

To mitigate the aspect blindness of training data performance to overfitting,
standard practice sets aside a holdout set disjoint from training and computes
performance separately. Thus, training a model ordinarily incorporates two
distinct steps: 1. optimization with training data to fit (model to) data and
2. verification of generalization by evaluating performance on holdout data.
While vague heuristics motivating this two-step procedure abound in the
literature and research community, rigorous statistical rationale less
ubiquitously accompany justification of its use. Moreover, this two-step
process facially treats training data and holdout data as altogether different
kinds of things, with different tasks and different intended uses. As such,
separating the conclusions we draw from training data and holdout data
threatens to undermine the original impetus according to which training data
is used for training in the first place, namely that optimization with respect
to training data should thereby optimize an expectation (generalization). We
explain the reasons for this paradox, and propose a solution that translates
into a statistical test which may be deployed for both defining and
identifying overfitting, using modified Law of Large Numbers (LLN) intuition
that empirical means should approximate their expectation.

In section 2, we review requisite background for the supervised learning
problem, discuss the problem with training data, how it relates to
overfitting, and why we would still like to use model performance on training
data to contribute to assessing generalization. In section 3, we detail the
statistical test for achieving this end, and give commentary on how this test
clarifies the meaning of overfitting. We point out how the test validates
generalization even absent strong but restrictive (e.g. PAC) learnability
guarantees; we also introduce a weaker but still rich notion of learnability.
We end with some plots in section 4 illustrating the use of the test in
simulation.

##  2 Technical Background

###  2.1 Supervised Machine Learning

The setting for a supervised machine learning problem starts with the
following data:

  1. 1.

a joint probability space
$(\mathcal{X}\times\mathcal{Y},\mathbb{P}_{\mathcal{X}\times\mathcal{Y}})$,222We
leave implicit the $\sigma$-algebra of measurable sets and suppose that
anything we try measuring is indeed $\mathbb{P}$-measurable.

  2. 2.

labeled data
$\mathsf{S}=\big{(}(x_{1},y_{1}),\ldots,(x_{m},y_{m})\big{)}\in(\mathcal{X}\times\mathcal{Y})^{\omega}\vcentcolon=\displaystyle\bigcup_{m\in\mathbb{N}}(\mathcal{X}\times\mathcal{Y})^{m}$,

  3. 3.

a hypothesis class $\mathcal{H}\subset\mathcal{Y}^{\mathcal{X}}$ of functions
$\tilde{y}:\mathcal{X}\rightarrow\mathcal{Y}$,333The notation
$\mathcal{Y}^{\mathcal{X}}$ denotes the set
$\big{\\{}\tilde{y}:\mathcal{X}\rightarrow\mathcal{Y}\big{\\}}$ of
unstructured functions with domain $\mathcal{X}$ and codomain $\mathcal{Y}$.
Of course, we require $\mathcal{H}$ to consist only of measurable such
functions. usually finite dimensional, elements $\tilde{y}\in\mathcal{H}$ of
which are called models, and

  4. 4.

a cost function generator
$c:\mathcal{H}\rightarrow\mathbb{R}^{\mathcal{X}\times\mathcal{Y}}$ mapping a
model $\tilde{y}$ to random variable
$c_{\tilde{y}}:\mathcal{X}\times\mathcal{Y}\rightarrow\mathbb{R}$, whose
output $c_{\tilde{y}}(x,y)$ on input $(x,y)$ is a measure of fit between
prediction $\tilde{y}(x)$ and label $y$.

The goal is to concoct an algorithm
$\hat{y}_{(\cdot)}:(\mathcal{X}\times\mathcal{Y})^{\omega}\rightarrow\mathcal{H}$
which outputs a model $\hat{y}_{\mathsf{S}}$ with small expected cost

$\mathbb{E}(c_{\hat{y}_{\mathsf{S}}})\approx\displaystyle\inf_{\tilde{y}\in\mathcal{H}}\mathbb{E}(c_{\tilde{y}}),$

having some guarantees of approximation performance in probability. The
measure $\mathbb{P}_{\mathcal{X}\times\mathcal{Y}}$ generating data
$(x_{i},y_{i})$ is usually unknown, and data $\mathsf{S}$ is used to proxy
approximate expectation and to optimize the expected risk function

$\begin{array}[]{lrcl}\mathbb{E}(c_{(\cdot)}):&\mathcal{H}&\rightarrow&\mathbb{R}\\\
&\tilde{y}&\mapsto&\mathbb{E}(c_{\tilde{y}}).\end{array}$ (1)

The standard algorithm for this optimization is empirical risk minimization,
namely

$\hat{y}_{\mathsf{S}}\in\displaystyle\arg\min_{\tilde{y}\in\mathcal{H}}e_{\mathsf{S}}(\tilde{y}),$
(2)

where empirical risk is defined as

$e_{\mathsf{S}}(\tilde{y})\vcentcolon=\frac{1}{|\mathsf{S}|}\displaystyle\sum_{(x,y)\in\mathsf{S}}c_{\tilde{y}}(x,y).$
(3)

Law of Large Numbers intuition suggests that

$e_{\mathsf{S}}(\tilde{y})\approx\mathbb{E}(c_{\tilde{y}})$ (4)

when $|\mathsf{S}|$ is large, so supposing as much, an output
$\hat{y}_{\mathsf{S}}$ of eq. 2 may be hoped to be a close approximation of
the true goal, in the sense that

$e_{\mathsf{S}}(\hat{y}_{\mathsf{S}})\approx\displaystyle\inf_{\tilde{y}\in\mathcal{H}}\mathbb{E}(c_{\tilde{y}}).$
(5)

To the extent that a model $\tilde{y}\in\mathcal{H}$ (approximately) satisfies
approximation (4), we say that the model generalizes
($\varepsilon$-generalizes if the error in approximation is bounded by
$\varepsilon$), and to the extent that models in $\mathcal{H}$ can be
guaranteed to generalize optimality
$\inf_{\tilde{y}}\mathbb{E}(c_{\tilde{y}})$ (5), we say that $\mathcal{H}$ is
some kind of learnable. The familiar and formal notion of probably
approximately correct (PAC) learnability, for example, extends guarantees of
concentration bounds to an optimization (over $\mathcal{H}$) context, and
defines $\mathcal{H}$ to be PAC learnable if there is a sample complexity
$\mu:(0,1)^{2}\rightarrow\mathbb{N}$ for which $\hat{y}_{\mathsf{S}}$ may be
guaranteed to $\varepsilon$-generalize with at least $1-\delta$ probability as
long as $|\mathsf{S}|>\mu(\varepsilon,\delta)$ ([9], [4]).444Explicitly, if
$m>\mu(\varepsilon,\delta)$ then
$\mathbb{P}_{(\mathcal{X}\times\mathcal{Y})^{m}}\left(\big{|}\mathbb{E}(c_{\hat{y}_{(\cdot)}})-\inf_{\tilde{y}\in\mathcal{H}}\mathbb{E}(c_{\tilde{y}})\big{|}>\varepsilon\right)<\delta$.
Strictly speaking, PAC learnability only requires the existence of an
algorithm
$\hat{y}:(\mathcal{X}\times\mathcal{Y})^{\omega}\rightarrow\mathcal{H}$
satisfying this bound, not necessarily that empirical risk minimization is it.
Properly quantifying the character and richness of $\mathcal{H}$ (as captured,
e.g., by VC dimension) demarcates learnability conditions, and various
theoretical results exist providing such guarantees.

###  2.2 Overfitting and Generalization

Absent formal learnability guarantees, it turns out that LLN reasoning is not
sufficient for ensuring generalization. The reasons are multifarious but
substantively turn around currying ([6, §2.3]) of the cost function generator
$c:\mathcal{H}\rightarrow\mathbb{R}^{(\mathcal{X}\times\mathcal{Y})}$. The
notion of currying reflects pre-fixing arguments of a multivariable function
to generate a function of fewer variables, and casting the learning objective
in this formalism is helpful for understanding the overfitting problem. For a
fixed model $\tilde{y}\in\mathcal{H}$, the map
$c_{\tilde{y}}:\mathcal{X}\times\mathcal{Y}\rightarrow\mathbb{R}$ is a random
variable, and therefore defines a measure $\mathbb{P}_{\mathbb{R}}$ on
$\mathbb{R}$ by
$\mathbb{P}_{\mathbb{R}}([a,b])\vcentcolon=\mathbb{P}_{\mathcal{X}\times\mathcal{Y}}\big{(}c_{\tilde{y}}^{-1}([a,b])\big{)}$.
This means, among other things, given data
$\mathsf{S}=\big{(}(x_{1},y_{1}),\ldots,(x_{m},y_{m})\big{)}\sim_{iid}\mathbb{P}_{\mathcal{X}\times\mathcal{Y}}$,
that
$c_{\tilde{y}}(\mathsf{S})\vcentcolon=\big{(}c_{\tilde{y}}(x_{1},y_{1}),\ldots,c_{\tilde{y}}(x_{m},y_{m})\big{)}\sim_{iid}\mathbb{P}_{\mathbb{R}}$.
And independence invites valid conclusions of various concentration results.

Figure 1: Slicing the Cost Function
$c:\mathcal{H}\times(\mathcal{X}\times\mathcal{Y})^{\omega}\rightarrow\mathbb{R}$

Searching over a function space, in the supervised machine learning setting,
adds complications to otherwise innocuous independence conclusions. For the
learning algorithm
$\hat{y}:(\mathcal{X}\times\mathcal{Y})^{\omega}\rightarrow\mathcal{H}$ first
takes data $\mathsf{S}\in(\mathcal{X}\times\mathcal{Y})^{\omega}$ in search of
a certain minimum with respect to this data. Given different data, the
algorithm outputs a different model. The curried cost generator, by contrast,
$c_{(\cdot)}:\mathcal{H}\rightarrow\mathbb{R}^{\mathcal{X}\times\mathcal{Y}}$
defines an empirical risk generator
$e(\cdot):\mathcal{H}\rightarrow\mathbb{R}^{(\mathcal{X}\times\mathcal{Y})^{\omega}}$
defined by sending $\tilde{y}\mapsto e_{(\cdot)}(\tilde{y})$, the latter of
which is defined by mapping data
$\mathsf{S}\in(\mathcal{X}\times\mathcal{Y})^{\omega}$ to
$e_{\mathsf{S}}(\tilde{y})$ (eq. 3), and with respect to which LLN reasoning
and the like may properly apply. The learning optimization procedure, however,
flips the currying around: fixing data $(x,y)\in\mathcal{X}\times\mathcal{Y}$,
we have a cost on models $c_{(\cdot)}(x,y):\mathcal{H}\rightarrow\mathbb{R}$
defined by $\tilde{y}\mapsto c_{\tilde{y}}(x,y)$, which extends to empirical
risk $e_{\mathsf{S}}(\cdot):\mathcal{H}\rightarrow\mathbb{R}$ mapping model
$\tilde{y}\mapsto e_{\mathsf{S}}(\tilde{y})$, instantiating the curried
function
$e_{(\cdot)}:\big{(}\mathcal{X}\times\mathcal{Y}\big{)}^{\omega}\rightarrow\mathbb{R}^{\mathcal{H}}$.555The
reversal of roles in subscripts between $c$ and $e$ is unfortunate, but
otherwise reflective of the primary purpose of each function, namely that
$c_{\tilde{y}}$ measures performance of model $\tilde{y}$ on a datapoint
$(x,y)\in\mathcal{X}\times\mathcal{Y}$ while $e_{\mathsf{S}}$ measures
empirical risk of fixed data on a model $\tilde{y}\in\mathcal{H}$.

Order of operations matter. Consider uncurried versions
${\mathcal{H}\times(\mathcal{X}\times\mathcal{Y})^{\omega}}$${\mathbb{R}}$$\scriptstyle{c_{(\cdot)}(\cdot)}$$\scriptstyle{e_{(\cdot)}(\cdot)}$
of the cost and empirical risk functions—thereby deprioritizing either data
$\mathsf{S}\in(\mathcal{X}\times\mathcal{Y})^{\omega}$ or model
$\tilde{y}\in\mathcal{H}$—and define

$\Gamma(\hat{y})\vcentcolon=\big{\\{}(\tilde{y},\mathsf{S})\in\mathcal{H}\times(\mathcal{X}\times\mathcal{Y})^{\omega}:\,\tilde{y}=\hat{y}_{\mathsf{S}}\big{\\}},$

the pullback of diagram
${\mathcal{H}\times(\mathcal{X}\times\mathcal{Y})^{\omega}}$${\mathcal{H}}$$\scriptstyle{\hat{y}_{(\cdot)}\circ\pi_{2}}$$\scriptstyle{id_{\mathcal{H}}\circ\pi_{1}}$
(see fig. 1). Evaluation of empirical performance for a model
$\hat{y}_{\mathsf{S}}$ using training data $\mathsf{S}$ lives in
$\Gamma(\hat{y})$, which elucidates why performance of training data is
“aspect blind” to overfitting: Law of Large Numbers reasoning does not apply
in this regime.

Consider that a sequence of datasets

$\begin{array}[]{lcl}\mathsf{S_{1}}&=&(x_{1},y_{1})\in(\mathcal{X}\times\mathcal{Y})^{1},\\\
\mathsf{S_{2}}&=&\big{(}\mathsf{S_{1}},(x_{2},y_{2})\big{)}\in(\mathcal{X}\times\mathcal{Y})^{2},\\\
&\vdots&\\\
\mathsf{S}_{m}&=&\big{(}\mathsf{S}_{m-1},(x_{m},y_{m})\big{)}\in(\mathcal{X}\times\mathcal{Y})^{m},\\\
&\vdots&\end{array}$

with each $\mathsf{S}_{j}\sim_{iid}\mathbb{P}_{\mathcal{X}\times\mathcal{Y}}$,
induces a sequence of models

$\hat{y}_{\mathsf{S}_{1}},\hat{y}_{\mathsf{S}_{2}},\ldots,\hat{y}_{\mathsf{S}_{m}},\ldots\in\mathcal{H}.$

The sequence of models consequently induces a sequence of finite sequences of
costs

$\begin{array}[]{lcl}c_{\hat{y}_{\mathsf{S}_{1}}}(\mathsf{S}_{1})&=&c_{\hat{y}_{\mathsf{S}_{1}}}(x_{1},y_{1})\in\mathbb{R},\\\
c_{\hat{y}_{\mathsf{S}_{2}}}(\mathsf{S}_{2})&=&\big{(}c_{\hat{y}_{\mathsf{S}_{2}}}(x_{1},y_{1}),c_{\hat{y}_{\mathsf{S}_{2}}}(x_{2},y_{2})\big{)}\in\mathbb{R}^{2},\\\
&\vdots&\\\
c_{\hat{y}_{\mathsf{S}_{m}}}(\mathsf{S}_{m})&=&\big{(}c_{\hat{y}_{\mathsf{S}_{m}}}(x_{1},y_{1}),\ldots,c_{\hat{y}_{\mathsf{S}_{m}}}(x_{m},y_{m})\big{)}\in\mathbb{R}^{m},\\\
&\vdots&\end{array}$

which clearly is not guaranteed to be iid, unless miraculously the cost
functions

$c_{\hat{y}_{\mathsf{S}_{1}}},c_{\hat{y}_{\mathsf{S}_{2}}},\ldots,c_{\hat{y}_{\mathsf{S}_{m}}},\ldots$

all induce the same measure
$\mathbb{P}_{c_{\hat{y}}(\mathcal{X}\times\mathcal{Y})}$ on $\mathbb{R}$, for
which there is no apriori reason to suppose. For each such
$\hat{y}_{\mathsf{S}_{j}}$, LLN still holds in the sense that there is, for
$\varepsilon,\delta>0$, a number $m_{j}>0$ for which

$\mathbb{P}_{(\mathcal{X}\times\mathcal{Y})^{m}}\left(\big{|}e_{(\cdot)}(\hat{y}_{\mathsf{S}_{j}})-\mathbb{E}(c_{\hat{y}_{\mathsf{S}_{j}}})\big{|}>\varepsilon\right)<\delta$

whenever $m>m_{j}$. Still, there is no reason to suppose that $m_{j}<j$, and
even if it were, small probability events still exist: the search over
$\mathcal{H}$ in the supervised learning setting incentivizes discovery of
such events, c.f. section 3.2.

One may place an appropriate measure $\mathbb{P}_{\Gamma}$ on
$\Gamma(\hat{y})$—in fact, this pullback naturally inherits from measures on
$(\mathcal{X}\times\mathcal{Y})^{\omega}$—and certainly iid samples
$\mathsf{S}\in(\mathcal{X}\times\mathcal{Y})^{\omega}$ induce iid samples
$\hat{y}_{\mathsf{S}}\sim_{iid}\mathbb{P}_{\Gamma}$, but statements of events
on this set do not extend to iid conditions on sequences of costs. For such,
we must fix and isolate our attention to slices
$\\{\hat{y}_{\mathsf{S}}\\}\times(\mathcal{X}\times\mathcal{Y})^{\omega}\subset\Gamma(\hat{y})$
(see fig. 2), for which sample
$\mathsf{S}^{\prime}=\big{(}(x_{1}^{\prime},y_{1}^{\prime}),\ldots,(x_{k},y_{k}^{\prime})\big{)}\sim_{iid}\mathbb{P}_{\mathcal{X}\times\mathcal{Y}}$
induces truly independent and identically distributed sample
$c_{\hat{y}_{\mathsf{S}}}(\mathsf{S}^{\prime})\vcentcolon=\big{(}c_{\hat{y}_{\mathsf{S}}}(x_{1}^{\prime},y_{1}^{\prime}),\ldots,c_{\hat{y}_{\mathsf{S}}}(x_{k}^{\prime},y_{k}^{\prime})\big{)}\sim_{iid}\mathbb{P}_{c_{\hat{y}_{\mathsf{S}}}(\mathcal{X}\times\mathcal{Y})}$.

Figure 2: Sequences of models
$\hat{y}_{\mathsf{S_{1}}},\hat{y}_{\mathsf{S_{2}}},\ldots,\hat{y}_{\mathsf{S_{m}}},\ldots$

While we therefore cannot rely on empirical risk
$e_{\mathsf{S}}(\hat{y}_{\mathsf{S}})$ by itself to reflect generalization
performance, we may in concert with
$e_{\mathsf{S}^{\prime}}(\hat{y}_{\mathsf{S}})$ for some other data
$\mathsf{S}^{\prime}\in(\mathcal{X}\times\mathcal{Y})^{\omega}$, usually
called a holdout or validation set. Typically performance at each training
stage is evaluated on the holdout set, and early stopping conditions verify
that validation performance continues to improve [5]. An onset of validation
performance degradation can be interpreted as indication of overfitting.
Illustrations of overfitting in the literature (e.g. [1], [3], [10]) display
performance on training data compared with performance on holdout data, often
parameterized by model complexity or training step ([8], [7]).

While discussions of overfitting many times consider generalization against
model complexity, and therefore present performance across models, we
introduce a test, for a fixed model, based on classic concentration
inequalities, with respect to which overfitting may be quantitatively defined,
relying on comparison of validation performance to training set performance.
We reason that because model construction uses and depends on (minimization
with respect to $\tilde{y}$ of) $e_{\mathsf{S}}(\tilde{y})$, we ought to be
able to conclude performance with it. In fact, comparison against empirical
risk $e_{\mathsf{S}}(\hat{y}_{\mathsf{S}})$ provides an anchor against which
we may draw rigorous statistical conclusions. The test we provide amounts to
much of the same as common stopping criteria for training, though the grounds
we give are both grounded in the math and provide threads for distinguishing
causes of error.

##  3 Detecting Overfitting

###  3.1 The Test

We consider only the case where cost
$c_{(\cdot)}:\mathcal{H}\times(\mathcal{X}\times\mathcal{Y})\rightarrow\mathbb{R}$
is bounded as $c_{\hat{y}_{\mathsf{S}}}\subset[0,1]$, such as most
classification problems or restricted classes of regression problems. In this
case, Hoeffding-like bounds abound and we expect that

$\mathbb{P}_{(\mathcal{X}\times\mathcal{Y})^{k}}\left(\left|\mathbb{E}(c_{\hat{y}_{\mathsf{S}}})-e_{(\cdot)}(\hat{y}_{\mathsf{S}})\right|>\varepsilon\right)<2e^{-2\varepsilon^{2}k}.$
(6)

In other words, for independently and identically distributed sampled data
$\mathsf{S}^{\prime}\in(\mathcal{X}\times\mathcal{Y})^{k}$,
$e_{\mathsf{S}^{\prime}}(\hat{y}_{\mathsf{S}})\approx\mathbb{E}(c_{\hat{y}_{\mathsf{S}}})\pm\varepsilon$
with probability at least $1-e^{-2\varepsilon^{2}k}$.666The fact that
$\mathbb{E}(c_{(\cdot)})$ and $e(\cdot)$ both take $\hat{y}_{\mathsf{S}}$ as
argument is irrelevant: the bound holds for any $\tilde{y}\in\mathcal{H}$.
While $\mathsf{S}\in(\mathcal{X}\times\mathcal{Y})^{m}$ is also drawn
independently, by assumption, we cannot quite conclude the same of
$e_{\mathsf{S}}(c_{\hat{y}_{\mathsf{S}}})$ because (as discussed above) with
respect to the $c_{\hat{y}_{\mathsf{S}}}$-induced measure on $\mathbb{R}$, the
sequence
$(c_{\hat{y}_{\mathsf{S}}}(x_{1},y_{1}),\ldots,c_{\hat{y}_{\mathsf{S}}}(x_{m},y_{m}))$
is not. We may, however, suppose that a consequence of independence holds,
namely that

$|\mathbb{E}(c_{\hat{y}_{\mathsf{S}}})-e_{\mathsf{S}}(\hat{y}_{\mathsf{S}})|<\varepsilon/2,$
(7)

and use this (possibly counterfactual) supposition to test its truth. While
possibly counterintuitive, a bound of the form in (7) is exactly what we
desire from a generalizing model $\hat{y}_{\mathsf{S}}$.

We first collect some definitions.

######  Definition 3.1.

Let $\mathsf{S}\sim_{iid}\mathbb{P}_{\mathcal{X}\times\mathcal{Y}}$ and
$\hat{y}_{\mathsf{S}}\in\mathcal{H}$. We say that $\hat{y}_{\mathsf{S}}$
$\varepsilon$-overfits $\mathsf{S}$ if

$e_{\mathsf{S}}(\hat{y}_{\mathsf{S}})<\mathbb{E}(c_{\hat{y}_{\mathsf{S}}})-\varepsilon.$

Similarly, $\hat{y}_{\mathsf{S}}$ $\varepsilon$-underfits $\mathsf{S}$ if
$e_{\mathsf{S}}(\hat{y}_{\mathsf{S}})>\mathbb{E}(c_{\hat{y}_{\mathsf{S}}})+\varepsilon$.
Finally, $\hat{y}_{\mathsf{S}}$ $\varepsilon$-generalizes model
$\hat{y}_{\mathsf{S}}$ neither $\varepsilon$-overfits nor
$\varepsilon$-underfits $\mathsf{S}$.

######  Proposition 3.1 (Test for Overfitting).

Suppose that model $\hat{y}_{\mathsf{S}}$ $\varepsilon/2$-generalizes
(definition 3.1, inequality (7)). Then

$\mathbb{P}_{(\mathcal{X}\times\mathcal{Y})^{k}}\left(\big{|}e_{\mathsf{S}}(\hat{y}_{\mathsf{S}})-e_{\mathsf{S}^{\prime}}(\hat{y}_{\mathsf{S}})\big{|}>\varepsilon\right)\leq
2e^{-\frac{\varepsilon^{2}k}{2}}.$ (8)

Therefore, the null hypothesis that trained model $\hat{y}_{\mathsf{S}}$
$\frac{\varepsilon}{2}$-generalizes may be tested using probability bound eq.
8.

###### Proof.

Since

$\begin{array}[]{ll}\big{|}e_{\mathsf{S}}(\hat{y}_{\mathsf{S}})-e_{\mathsf{S}^{\prime}}(\hat{y}_{\mathsf{S}})\big{|}&=\big{|}e_{\mathsf{S}}(\hat{y}_{\mathsf{S}})-\mathbb{E}(c_{\hat{y}_{\mathsf{S}}})+\mathbb{E}(c_{\hat{y}_{\mathsf{S}}})-e_{\mathsf{S}^{\prime}}(\hat{y}_{\mathsf{S}})\big{|}\\\
&\leq\big{|}e_{\mathsf{S}}(\hat{y}_{\mathsf{S}})-\mathbb{E}(c_{\hat{y}_{\mathsf{S}}})\big{|}+\big{|}\mathbb{E}(c_{\hat{y}_{\mathsf{S}}})-e_{\mathsf{S}^{\prime}}(\hat{y}_{\mathsf{S}})\big{|}\\\
&<\frac{\varepsilon}{2}+\big{|}\mathbb{E}(c_{\hat{y}_{\mathsf{S}}})-e_{\mathsf{S}^{\prime}}(\hat{y}_{\mathsf{S}})\big{|}\end{array}$

we conclude

$\left\\{\big{|}e_{\mathsf{S}}(\hat{y}_{\mathsf{S}})-e_{(\cdot)}(\hat{y}_{\mathsf{S}})\big{|}>\varepsilon\right\\}\subseteq\left\\{\big{|}\mathbb{E}(c_{\hat{y}_{\mathsf{S}}})-e_{(\cdot)}(\hat{y}_{\mathsf{S}})\big{|}>\frac{\varepsilon}{2}\right\\}.$

Inclusion of events implies inequality of measures, and we apply Hoeffding
(inequality (9)) to bound the right hand side probability
$\mathbb{P}\left(\big{|}\mathbb{E}(c_{\hat{y}_{\mathsf{S}}})-e_{(\cdot)}(\hat{y}_{\mathsf{S}})\big{|}>\frac{\varepsilon}{2}\right).$
∎

Notice that use of holdout data for evaluation by itself provides an absolute
approximation of performance, while in tandem with training data, we gain
quantified (un)certainty specifically about generalization. Finally, the
probability in (8) depends on the size of validation data, but not on the size
of training data. This conclusion is correct: while we would like more
training data to correlate with higher likelihood of performance, the problem
in section 2.2 indicates that such intuition may not find a straightforward
grounding in probability. Presumably, one may be less inclined to hypothesize
satisfactory model performance when training with little data. The intuition
finds security in PAC learnability, absent which there is no obvious
guaranteed connection between size of (training) data and performance; we
discuss this issue further in section 3.3.

###  3.2 Interpreting the Output

Overfitting is a heuristic notion which suggests a model has fit the data and
not the distribution which generated it. On closer inspection, however, the
test we propose does not provide indication of only overfitting. In fact, the
supposition of generalization is one with respect to a certain (fixed)
distribution; this test thus additionally assumes that the test data
$\mathsf{S}^{\prime}\sim_{iid}\mathbb{P}_{\mathcal{X}\times\mathcal{Y}}$ as
well. It may not be. For there may be some form of distributional shift
according to which
$\mathsf{S}^{\prime}\sim_{iid}\mathbb{P}^{\prime}_{\mathcal{X}\times\mathcal{Y}}$,
in which case we cannot guarantee the bound in (8), at least not if the
expectation $\mathbb{E}(c_{\hat{y}_{\mathsf{S}}})$ is computed with respect to
the original measure $d\mathbb{P}_{\mathcal{X}\times\mathcal{Y}}$. In other
words, instantiation of event
$\left\\{\big{|}e_{\mathsf{S}}(\hat{y}_{\mathsf{S}})-e_{(\cdot)}(\hat{y}_{\mathsf{S}})\big{|}>\varepsilon\right\\}$
by inequality
$|e_{\mathsf{S}}(\hat{y}_{\mathsf{S}})-e_{\mathsf{S}^{\prime}}(\hat{y}_{\mathsf{S}})|>\varepsilon$
may suggest:

  1. 1.

an unlikely sample $\mathsf{S}^{\prime}$ was received (all the hypotheses
hold),

  2. 2.

$\hat{y}_{\mathsf{S}}$ does not generalize
$\mathbb{P}_{\mathcal{X}\times\mathcal{Y}}$ with respect to $c_{(\cdot)}$
(overfitting), or

  3. 3.

$\mathsf{S}^{\prime}\not\sim_{iid}\mathbb{P}_{\mathcal{X}\times\mathcal{Y}}$
(possible distributional shift).

It is important when running a statistical test to respect the scope of what
it purports to evaluate: namely, if a set of assumptions hold—in this case 1.
that $\hat{y}_{\mathsf{S}}$ $\frac{\varepsilon}{2}$-generalizes (eq. 7) and 2.
$\mathsf{S}^{\prime}\sim_{iid}\mathbb{P}_{\mathcal{X}\times\mathcal{Y}}$—then
the probability that a certain kind of event occurs is bounded by some value
which is explicitly calculable. Realization of the unlikely and unlucky event
by $\mathsf{S}^{\prime}$ can either mean $\mathsf{S}^{\prime}$ really is
unlucky or that one of the assumptions fails.

Finally, while this test is expressed with respect to the cost function
$c_{\tilde{y}}$ or $c_{\hat{y}_{\mathsf{S}}}$, it need not be so limited. In
fact, any map $f:\mathcal{X}\times\mathcal{Y}\rightarrow\mathbb{R}$ may be
used to probe the distribution, substituting the appropriate concentration
inequality depending on the range of $f$. When
$f(\mathcal{X}\times\mathcal{Y})$ is bounded, we may rely on a version of
Hoeffding, which converges exponentially. Subsequent work will investigate the
use of random projections to examine distribution shift and uncertainty
quantification, as a means of testing to eliminate or isolate the above
obfuscating condition #3.

###  3.3 Loosening Uniform Bounds

Figure 3: Model $o$ overfits.

We conclude with commentary on the merits of this test. The bound in eq. 8 is
perhaps unsurprising and at first glance offers little upside beyond the
performance guarantee as provided by eq. 6, which ensures approximation of
empirical mean (of holdout data) to the true mean. Indeed, one may argue that
overfitting, in the sense of definition 3.1, induces little cause for concern:
as long as performance $\mathbb{E}(c_{\hat{y}_{\mathsf{S}}})$ is “good
enough,” (as approximated by $e_{\mathsf{S}^{\prime}}(\hat{y}_{\mathsf{S}})$)
it may not particularly matter whether or that training data performance
matches a model’s generalization performance. On the other hand, guarantees of
the sort which PAC learnability provides ensure that the output of a training
algorithm is near optimal in a hypothesis class. In the presence of
overfitting, one may not know whether better than ‘good enough’ is achievable.
Generalization with training data provides confidence that empirical risk
minimization (2) approximately realizes risk minimization (5) in the absence
of uniform (PAC) guarantees. The test is a workable mechanism for checking
that there is little gap between performance a hypothesis class may achieve on
data and on the data’s distribution.

Consider, for example, fig. 3 which compares level sets for
$\mathbb{E}(c_{(\cdot)})$ and $e_{\mathsf{S}}(\cdot)$. Learnability, as
described by uniform convergence and notions of representability (c.f. [9,
§4.1]), guarantees that these profiles roughly track each other, which is
sufficient for generalization of output model $\hat{y}_{\mathsf{S}}$: if the
value of $\mathbb{E}(c_{(\cdot)})$ and $e_{\mathsf{S}}(\cdot)$ are roughly
approximate everywhere in $\mathcal{H}$, then they certainly are at a
particular point. On the other hand, learnability objectives ultimately seek
generalization of the output, namely that $\mathbb{E}(c_{(\cdot)})$ and
$e_{\mathsf{S}}(\cdot)$ are roughly approximate at $\hat{y}_{\mathsf{S}}$; how
they compare in other regions of $\mathcal{H}$ may be immaterial.

We underscore the point. PAC results guarantee not only that an algorithm will
return an optimal (in the hypothesis class) model, but that the sample
complexity with respect to which the algorithm is expected to reliably work is
independent of distribution. Guarantees of this form are helpful in providing
confidence ahead of time that the learning endeavor is not misguided. On the
other hand, practitioners often engage in the tackling the learning problem
irrespective of knowledge or other assurances that their class is PAC
learnable. Moreover, PAC learnability does not cover the intermediate case
that some distributions may require a larger sample complexity (some tasks are
harder to learn than others), and that there may be no uniform bound over all
measures, even if there are some over subsets. Still, assurance that the
output of training generalizes does not require that the hypothesis class be
PAC learnable, i.e. that uniform bounds hold. Rather: uniform bounds, when
they exist, provide a conceptual framework and analytic setting wherein a
class of results may be generated, in the absence of which, we would
nevertheless like to be able to say something.

We collect this commentary into a definition.

######  Definition 3.2.

Let $\mathcal{P}$ be a collection of probability measures on
$\mathcal{X}\times\mathcal{Y}$. We say that a hypothesis class
$\mathcal{H}\subset\mathcal{Y}^{\mathcal{X}}$ is $\mathcal{P}$-learnable if
there is sample complexity $\mu:(0,1)\rightarrow\mathbb{N}$ and algorithm
$\hat{y}:(\mathcal{X}\times\mathcal{Y})^{\omega}\rightarrow\mathcal{H}$ for
which
$\mathbb{P}_{(\mathcal{X}\times\mathcal{Y})^{m}}\left(\mathbb{E}(c_{\hat{y}_{(\cdot)}})-\inf_{\tilde{y}\in\mathcal{H}}\mathbb{E}(c_{\tilde{y}})>\varepsilon\right)<\delta$
whenever $m>\mu(\varepsilon,\delta)$ and
$\mathbb{P}_{\mathcal{X}\times\mathcal{Y}}\in\mathcal{P}$. We say that
$\mathcal{H}$ is $\mathbb{P}_{\mathcal{X}\times\mathcal{Y}}$-learnable if it
is $\\{\mathbb{P}_{\mathcal{X}\times\mathcal{Y}}\\}$-learnable.

Whereas PAC learnability typically guarantees (approximate) optimality
independent of measure, $\mathcal{P}$-learnability expressly restricts
measures for which $\mathcal{H}$ is (uniformly) suited to learn. While PAC
learnability is powerful in providing guarantees absent prior assumptions
about the measure from which data is drawn, this generality also inhibits the
usefulness of prior knowledge: for a class is PAC learnable or not
irrespective of what is known regarding the data. One might imagine, for
example, that assurance data is subgaussian may be relevant, in that
subgaussianity demarcates a class of measures for which a hypothesis class is
suited to fit.

##  4 In Simulation

Code implementing this test can be found at
https://github.com/schmidttgenstein/qudost.git. In each of the following
illustrations, we plot empirical performance (accuracy for classification
problems) with respect to both training and holdout data on the left, and the
absolute difference on the right. These curves are plotted against training
epoch, and each pair uses a different size for $\mathsf{S}^{\prime}$ with
respect to which either probability or precision depend (appendix A). We fix
the precision in each case, denoted by the dashed red line in the right
figures, and report the probability bound a lá proposition 3.1. Per
Hoeffding’s inequality (9), this precision may be made finer with more data.
Worth noting that the test in proposition 3.1 does not intrinsically relate to
early stopping: a model may overfit and cease to overfit at various epochs in
training (see, e.g., fig. 5).

Results in fig. 4 and fig. 5 use generated data and a multilayer perceptron
for binary classification. Results in fig. 6 and fig. 7 use a simple convnet
on MNIST data.

###  4.1 Simulated Data

Figure 4: $k=15,000$

Figure 5: $k=50,000$

###  4.2 MNIST

Figure 6: MNIST $k=1000$

Figure 7: MNIST $k=6000$

##  Appendix A Hoeffding’s Inequality for Statistical Hypothesis Testing

Hoeffding’s inequality gives a probability bound for independent sample
$\mathsf{S}=(x_{1},\ldots,x_{m})\sim_{iid}\mathbb{P}_{\mathcal{X}}$ when
$\mathcal{X}=[0,1]$, namely:

$\mathbb{P}_{\mathcal{X}}\left(\left|\displaystyle\int_{\mathcal{X}}xd\mathbb{P}_{\mathcal{X}}(x)-\frac{1}{|\mathsf{S}|}\displaystyle\sum_{x\in\mathsf{S}}x\right|>\varepsilon\right)<2e^{-2\varepsilon^{2}|\mathsf{S}|}.$
(9)

Therefore, given any two of confidence specification $\delta\in(0,1)$, data
set sized $|\mathsf{S}|=m$, and precision bound $\varepsilon\in(0,1)$, one may
readily solve for the third.

Proof of its verity and other applications may be found in various probability
texts ([2], [9], [4]).

## References

  * [1] J. Bilmes.  Underfitting and overfitting in machine learning.  UW ECE Course Notes, 2020.
  * [2] S. Boucheron, G. Lugosi, and P. Massart.  Concentration Inequalities: A Nonasymptotic Theory of Independence.  Oxford University Press, 02 2013.
  * [3] B. Ghojogh and M. Crowley.  The theory behind overfitting, cross validation, regularization, bagging, and boosting: Tutorial.  arXiv [stat.ML], 2019.
  * [4] M. Mohri, A. Rostamizadeh, and A. Talwalkar.  Foundations of Machine Learning.  MIT Press, 2018.
  * [5] L. Prechelt.  Early Stopping — But When?, pages 53–67.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2012.
  * [6] E. Riehl.  Category Theory in Context.  Dover, 2016.
  * [7] R. Roelofs.  Measuring generalization and overfitting in machine learning.  2019\.
  * [8] R. Roelofs, V. Shankar, B. Recht, S. Fridovich-Keil, M. Hardt, J. Miller, and L. Schmidt.  A meta-analysis of overfitting in machine learning.  In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
  * [9] S. Shalev-Shwartz and S. Ben-David.  Understanding Machine Learning.  Cambridge University Press, 2014.
  * [10] M. Valdenegro-Toro and M. Sabatelli.  Machine learning students overfit to overfitting.  arXiv [cs.LG].