thesis_berceanu.txt

Superfluidity, the ability of a fluid to flow without apparent
viscosity, is one of the most striking consequences of collective
quantum coherence, with manifestations ranging from metastability of
supercurrents in multiply connected geometries to the appearance of
quantized vortices, or the existence of a critical velocity for
frictionless flow when scattering against a defect. While
traditionally investigated in equilibrium systems, like liquid
He and ultracold atomic gases, experimental advances in
nonlinear optics, in particular regarding microcavity
exciton-polaritons, paved the way for studying superfluid-related
phenomena in a driven-dissipative framework.

Equally exciting is the possibility of realising topological phases of
matter, such as the integer or fractional quantum Hall states, outside
of traditional electronic systems. Photonics experiments in
driven-dissipative resonator arrays, in particular, offer a high
degree of controllability and tunability, as well as unprecedented
experimental access to the eigenstates and energy spectrum.

This thesis reports on hydrodynamic effects, as well as topological
properties, of driven-dissipative systems. In particular, we analyze
the superfluid-like behaviour of microcavity exciton-polaritons, as
well as the momentum-space topology of coupled resonator arrays.

Microcavity exciton-polaritons are quasiparticles resulting from the
mixing of excitons (bound electron-hole pairs) and photons confined
inside semiconductor microcavities. While polariton fluids have been
shown to display collective coherence, the connection between the
various manifestations of superfluid behaviour is more involved
compared to equilibrium systems. In this manuscipt, we consider both
the case of a single-fluid pump-only configuration, as well as the
three-fluid optical parametric oscillator regime that results from
parametric scattering of the pump to the signal and idler states. In
both cases, we look at the response of the moving polaritons
scattering against a weak static defect present in the microcavity.

For the single fluid, we evaluate analytically the drag exerted by the
fluid on the defect. For low fluid velocities, the pump frequency
classifies the collective excitation spectra in three different
categories: linear, diffusive-like and gapped. We show that both the
linear and diffusive-like cases share a qualitatively similar
crossover of the drag from the subsonic to the supersonic regime as a
function of the fluid velocity, with a critical velocity given by the
speed of sound found for the linear regime. In contrast, for gapped
spectra, we find that the critical velocity exceeds the speed of
sound. In all cases, we show that the residual drag in the subcritical
regime is caused by the nonequilibrium nature of the system. Also,
well below the critical velocity, the drag varies linearly with the
polariton lifetime, in agreement with previous numerical studies.

The optical parametric oscillator regime presents an additional
challenge, as one is dealing with three coupled fluids. The
spontaneous macroscopic coherence following the phase locking of the
signal and idler fluids has been already shown to be responsible for
their simultaneous quantized flow metastability. We find that the
modulations generated by the defect in each fluid are not only
determined by its associated scattering ring in momentum space, but
each component displays additional rings because of the cross-talk
with the other components imposed by nonlinear and parametric
processes. We single out three factors determining which one of these
rings has the biggest influence on each fluid response: the coupling
strength between the three fluids, the resonance of the ring with the
polariton dispersion, and the values of each fluid group velocity and
lifetime together establishing how far each modulation can propagate
from the defect.  For the typical conditions of parametric scattering,
the pump is in the supercritical regime, so the signal and idler will
show the modulations of the pump, meaning none of the three states
manifests superfluid behaviour. However, the signal appears to flow
without friction in the experimental study, because the three factors
mentioned above conspire to reduce the amplitude of its modulations
below currently detectable levels.

Driven-dissipative systems can show interesting phenomena also without
interactions, stemming from the nontrivial topology of their energy
bands. In the final part of this thesis, we present a realistic
proposal for an optical experiment using state-of-the-art coupled
resonator arrays. We study theoretically the driven-dissipative
Harper-Hofstadter model on a square lattice in the presence of a weak
harmonic trap. Without pumping and losses, the eigenstates of this
system can be understood, under certain approximations, as
momentum-space toroidal Landau levels, where the Berry curvature, a
geometrical property of an energy band, acts like a momentum-space
magnetic field. We show how key features of these eigenstates can be
observed in the steady-state of the driven-dissipative system under a
monochromatic coherent drive. We also show that momentum-space Landau
levels would have clear signatures in spectroscopic measurements in
such experiments, and we discuss the insights gained in this way into
geometrical energy bands and particles in magnetic fields.








List of Publications
List of Publications

The following articles have been published in the context of this thesis:









  Onset and dynamics of vortex-antivortex pairs in polariton optical parametric oscillator superfluids,
  G. Tosi, F. M. Marchetti, D. Sanvitto, C. Anton, M. H. Szymanska, A. C. Berceanu, C. Tejedor, L. Marrucci, A. Lemaitre, J. Bloch, and L. Vina,
  Phys. Rev. Lett. 107, 036401 (2011).

  Drag in a resonantly driven polariton fluid,

  A. C. Berceanu, E. Cancellieri, and F. M. Marchetti,

  J. Phys.: Condens. Matter 24, 235802 (2012) [Chapter ].

  Multicomponent polariton superfluidity in the optical parametric oscillator regime,
  A. C. Berceanu, L. Dominici, I. Carusotto, D. Ballarini,

  E. Cancellieri, G. Gigli, M. H. Szymanska, D. Sanvitto,
 and F. M. Marchetti,
  Phys. Rev. B 92, 035307 (2015) [Chapter ].

  Momentum-space Landau levels in driven-dissipative cavity arrays,

  A. C. Berceanu, H. M. Price, T. Ozawa, and I. Carusotto,

  Phys. Rev. A 93, 013827 (2016) [Chapter ].


Conference proceedings:


  Momentum-space Landau levels in arrays of coupled ring resonators,

  H. M. Price, A. C. Berceanu, T. Ozawa, and I. Carusotto,

  Proc. SPIE 9762, 97620W (2016).









Preface
Preface

How birds fly together
Systems far from thermal equilibrium frequently show novel features
when compared to their equilibrium counterparts. As an everyday
example, consider a group of birds displaying long-range ordered
behaviour, manifested by forming a flock under certain
conditions. This behaviour can be modeled by introducing a time step
rule, such that each individual bird in a group determines its next
direction on each time step by averaging the directions of its
neighbours and adding some random noise on top of
that . It can be shown that, in the limit of the
velocity going to zero, the model reduces to the XY model in two
dimensions, where the spin is represented by the bird velocity. Since
the 2D XY model does not spontaneously break the symmetry at any
finite temperature (as justified by the Mermin-Wagner theorem), one
can show that the appearance of the long-range ordered phase is a
direct consequence of nonequilibrium aspects of the model. In a
nutshell, the neighbours of one particular bird will be different at
different times, depending on the velocity field. This gives rise to a
time-dependent variable-ranged interaction, which can stabilize the
ordered phase.

Driven-dissipative systems
Non-equilibrium driven-dissipative photonic systems, such as
polaritons in semiconductor microcavities or arrays of coupled optical
resonators, have recently attracted a lot of interest due to the
possibility of observing quantum phenomena which normally require very
low temperatures and/or intense magnetic fields and are traditionally
restricted to the domain of solid-state systems or ultracold atomic
gases. Besides being highly tunable, these optical systems facilitate
direct experimental access to observables such as the wavefunction or
energy spectrum, all at room temperature. In particular, microcavity
polaritons have allowed the observation of collective hydrodynamic
phenomena, ranging from frictionless flow around a small defect to the
formation of quantized vortices and dark solitons at the surface of
large impenetrable obstacles , while
ring-resonator arrays coupled to artificial magnetic fields have
recently allowed engineering topological edge states robust to
disorder .


Polaritons
Microcavity polaritons are quasiparticles resulting from the strong
coupling of cavity photons and quantum well excitons 
, and have the prerogative of being both easy to
manipulate, via an external laser, and detect, via the light escaping
from the cavity . The finite polariton lifetime
establishes the system as intrinsically out of equilibrium: an
external pump is needed to continuously replenish the cavity of
polaritons, that quickly, on a scale of tens of picoseconds,
escape. The pumping can be done resonantly, close to the polariton
energy dipersion, or non-resonantly.

Incoherent pumping: polariton BEC
The landmark observation  of polariton
condensation in 2006, as shown in Fig. , was achieved
using a non-resonant experimental setup. In the experiments, the
system was incoherently excited by a laser beam tuned at a very high
energy. Relaxation of the excess energy 
 lead to a population of the cavity polariton states
and, above a certain laser power threshold, to Bose-Einstein
condensation into the lowest polariton state.

    
    
    Experimental observation of polariton Bose-Einstein
    condensation. A sharp and intense peak corresponding to the lowest
    momentum state is formed in the center of the far-field emission [top
    panels], with increasing pump power (left to right). The
    corresponding energy-resolved emission [bottom panels] show that
    above the condensation threshold, the emission comes almost entirely
    from the lowest energy state situated at the bottom of the polariton
    dispersion. From Ref. .
    
  

Coherent pumping: OPO
Due to their energy dispersion and strong nonlinearity inherited from
the excitonic component, polaritons resonantly injected by the
external laser into the pump state with a suitable wavevector and
energy can undergo coherent stimulated scattering into two conjugate
states , called the signal and
the idler, in a process known as optical parametric oscillator
(OPO). Since their first realisation 
, the
interest in microcavity optical parametric phenomena has involved
several fields of fundamental and applicative
research 
. 


      
    
    Top panels: Simulated (right panel) and observed (left panel)
    edge state propagation in a ring-resonator array.  The light enters
    from one corner and exits from the other, as signaled by the
    arrows. Adapted from Ref. .  Bottom
      panels: Simulated edge state propagation in a polariton system.
    Intensities of the exciton (left panel) and photon (right panel)
    fields obtained when pumping a lattice at the center of its lower
    edge. From Ref. .
    
  
Superfluidity
The superfluid properties of a resonantly pumped polariton quantum
fluid, both in the pump-only configuration (without parametric
scattering) as well as the OPO regime, have been actively investigated
experimentally, as well as theoretically . In
particular, a supression of scattering in the pump-only case was
observed  below a critical velocity, similar to what
has been predicted by the Landau criterion for equilibrium superfluid
condensates.

While in equilibrium condensates different aspects of superfluidity
are typically closely related , this is no longer
true in a non-equilibrium context.  Independent of the pumping scheme,
the driving and the polariton finite lifetime force one to reconsider
the meaning of superfluid behaviour, when the spectrum of collective
excitations is complex rather than real, raising question about the
applicability of a Landau criterion .

An additional complexity characterises the OPO regime, namely, the
simultaneous presence of three oscillation frequencies and momenta for
pump, signal and idler correspondingly increases up to 12 (or 6,
depending on the approximation used to describe the system) the number
of collective excitation branches . Note that from
the experimental point of view, pioneering
experiments  have observed a ballistic nonspreading
propagation of signal/idler polariton wavepackets in a triggered-OPO
configuration, as well as demonstrating the existence and
metastability of vortex configurations  in the
signal and idler.
 

Topology in polaritons
Very recently, the polariton community started to explore topological
effects in polariton lattices 

. The idea is to break time reversal symmetry by
the application of a strong external magnetic field, giving rise to
energy bands with nontrivial topology. In particular,
Ref.  proposes an exciton-photon coupling with
a winding phase in momentum space, giving rise to polaritonic bands
with chiral edge modes that allow unidirectional propagation,
protected against backscattering. The bottom panels of
Fig.  show the exciton and photon fields
traveling together as a unique polaritonic counter-clockwise chiral
edge mode. In this context, dissipation in the form of photon losses
provides a coupling to the outside continuum of modes and, as such,
can be used to detect or inject edge modes.

It would be interesting to study the feasability of coupled
micropillars for this type of physics, in view of the recent
experimental findings of a linear graphene-like
dispersion  as well as edge
states  in honeycomb micropillar lattices.  

On the other hand, topologically protected edge states have already
been experimentaly observed  in the case of
silicon-based optical resonator arrays (see top panels of
Fig. ), due to recent advances in creating
synthetic gauge fields, which have opened new horizons for simulating
topological phases of matter also with neutral particles, such as
photons  (or ultracold
atoms 
). Rather than simply replicating previous
measurements, experiments with synthetic gauge fields allow for
unprecedented access to properties such as the eigenstates or
eigenspectrum, while the tunability and controllability of these
experiments offer the prospect of simulating novel physics.


Contents of this thesis
Thesis layout
This manuscript is split in two parts: part I is a detailed review of
the basic concepts, including both the theoretical formalism, as well
as the relevant experiments, necessary for understanding part II,
which presents the three main works published as part of my PhD
(Chapters ,  and ). The two
systems chosen for ilustrating the basic physical concepts are
ultracold atomic gases (Chapter ) and microcavity
exciton-polaritons (Chapter ). We now give a brief
description of the content of each of the chapters.

In Chapter , we describe the phenomenon of
Bose-Einstein condensation, and introduce the Gross-Pitaevskii
equation which is extended in Chapter  to the case
of polariton condensates. We then review the linear response of a
moving atomic condensate to a weak stationary defect, leading to
discussion on superfluidity. This scattering problem is the
equilibrium conterpart of the problems studied in
Chapters  and  in the context of polaritons
in the resonantly pumped and optical parametric oscillator regimes,
respectively. Finally, we introduce the Harper-Hofstadter model, which
is the backbone of Chapter .

In Chapter , we introduce microcavity
exciton-polaritons, and their theoretical description in terms of a
driven-dissipative Gross-Pitaevskii equation. We discuss both the
resonantly pumped system which is used in Chapter , as
well as the optical parametric oscillator of Chapter ,
and we make the connection to the superfluid-related phenomena
previously introduced in Chapter .

In Chapter , we study the scattering of a resonantly
pumped polariton condensate against a static defect of the
microcavity, for the simplified case of small fluid velocities, when
the polariton dispersion can be considered quadratic. After discussing
the effects of dissipation on the Bogoliubov excitation spectra, we
find that the finite polariton lifetime also affects the drag exerted
by the condensate on the defect. In particular, we show that there is
a nonzero drag force, entirely due to the out-of-equilibrium nature of
the system, even in the "superfluid regime". Finally, we
characterise the behaviour of the drag as a function of the condensate
velocity and polariton lifetime.

In Chapter , we present a theoretical and experimental
study of the same scattering problem, now in the context of an optical
parametric oscillator consisting of three coupled condensates: the
pump, the signal and the idler. Apart from using the linear response
analysis first introduced in Chapter  and then
extended to the case of a pump-only condensate in
Chapter , we also numerically solve the full
driven-dissipative Gross-Pitaevskii equation. We show that, while the
modulation patterns are present in all three condensates, their
relative amplitudes depend on various factors. In particular, for the
typical experimental conditions that favour parametric scattering, the
signal has undetectable modulations, while the pump and idler show the
same response.

In Chapter , we add a harmonic trap to the
Harper-Hofstadter lattice model of Chapter . We
discuss how the eigenstates of the new Hamiltonian can be seen as
momentum-space Landau levels, by making use of the nontrivial topology
of the Harper-Hofstadter energy bands. We then extend the model to
include driving and dissipation, and finally present a proposal for
the physical implementation of the model in state-of-the-art
driven-dissipative photonic systems.

The source code belonging to Chapters , and  is available online .








Condensed-matter systems

Ultracold atomic gases

Laser cooling  allowed
achieving temperatures in the micro-Kelvin regime, and eventually led
to the realization of optical lattices . It
also paved the way for more powerful cooling techniques, such as
evaporative cooling, which made possible the Bose-Einstein
condensation of dilute atomic
gases .

Ultracold gases have been at the forefront of simulating quantum
phenomena with analogs throughout physics, from nonlinear optics to
condensed matter systems . In particular,
their link to quantum simulation of condensed matter phenomena becomes
obvious when adding optical lattice
potentials , combined with synthetic magnetic
fields .

In this Chapter, we review the physics of ultracold atomic gases, with
an emphasis on their link to hydrodynamic effects such as
superfluidity, as well as their connection to traditional solid state
lattice models such as the celebrated Harper-Hofstadter model which
originally describes the single-particle physics of band electrons in
intense magnetic fields.

Chapter organization
This Chapter is organized as follows: in Section  we
present a short history of Bose-Einstein condensation, describe its
main features and state its formal definition in terms of the
Penrose-Onsager criterion, leading to the weakly-interacting Bose gas
paradigm and the Gross-Pitaevskii equation (Sec. ), an
essential theoretical tool for the mean-field description of atomic
condensates. In Section  we introduce the
linear response formalism, which proves useful for interpreting
experiments where a weak perturbation is applied to the condensate. We
employ this formalism in order to study a scattering problem
concerning the flow of a BEC in the presence of a weak static defect
(Sec. ). In this context, we review
Bogoliubov's excitation spectrum and its associated Landau criterion
for superfluidity, followed by a detailed discussion on superfluidity
and related phenomena, such as quantized vortices, in
Sec. . We briefly touch on the subject of
synthetic gauge fields for neutral atoms, before investigating the
properties of BECs in periodic potentials created by optical lattices
(Sec. ). Finally, we combine the concepts of
synthetic gauge fields and optical lattices in
Section , where we show the main features of the
Harper-Hofstadter model, which was recently realized using
atomic gases.




Bose-Einstein condensation


BEC in noninteracting gas
In 1925, Albert Einstein (prompted by the earlier work of the indian
polyglot Satyendra Nath Bose) considered what would happen to a
non-interacting bosonic gas of non-relativistic particles in the
thermodynamic limit, as one lowers the temperature. He predicted the
phenomenon we now call Bose-Einstein condensation (BEC),
namely a phase transition to a new state of matter, in which a finite
fraction of all the particles would occupy the same single-particle
state. The transition occurs at fixed density below a critical
temperature  or, alternatively, at fixed temperature, above a
critical density. In particular, if we take  neutral particles in a
cubic box of volume , then they would predominantly occupy the
zero-wavevector state , and the critical temperature would
be 


with the density  and ,  being the particle
mass and Boltzmann's constant, respectively.

At its core, BEC is a paradigm of quantum statistical mechanics,
stemming from the indistinguishability of elementary particles and the
Bose-Einstein statistics that they obey. One can hand-wavingly deduce
the critical temperature (or critical density) where quantum
degeneracy would start playing a role in a many-body system, by
arguing that the thermal de Broglie wavelength should be comparable to
or greater than the inter-particle distance (which in our case is
 on average) . Apart from the
numerical prefactor, we get the same answer as
Eq. eq:Tc3D. While one may argue that elementary massive
bosons do not exist, it is worth emphasizing that indistinguishability
only plays a role when there is a finite probability for exchange
processes to occur between the particles. In that sense, all
odd-isotope alkali atoms under relevant experimental conditions (see
below) effectively behave as bosons: their many-body wavefunction is
symmetric under the exchange of any two such atoms.

BEC in interacting system
Interestingly enough, BEC was considered by many at the time to be a
pathological behaviour of the non-interacting gas, which would resolve
once interactions were properly accounted for. In fact, it is well
known that the ideal Bose gas has infinite compresibility.  This
pathology is cured by introducing a weak repulsive interaction between
bosons, a regime where BEC survives, as we will see next.

One-body density matrix
Following Leggett, we characterise each of the  particles (assumed
spinless, for simplicity) by a position vector , with the label
 running from 1 to . Any pure state  of the (now interacting)
system - which can be also subjected to an external potential - can
be described at time  by the many-body wavefunction
. Therefore, the most general state
of the system (also called mixed state) can be written as a
superposition of pure orthonormal states  with different weights
. The single-particle density matrix
 represents the probability
amplitude, at time , of finding a specific particle at position
, multiplied by the amplitude of finding it at ,
and averaged over the positions of all the other particles:

  _1(,^,t) & N _s p_s d_2d_3d_N
_s^(,_2,,_N,t)_s(^,_2,,_N,t)

& = _i n_i(t) _i^(,t) _i(^,t)

where in the second line we have re-written the density matrix in
diagonal form, introducing its eigenvalues  and eigenvectors
, which form a complete orthonormal set at any time  (here
 labels a good quantum number of the problem, i.e. momentum in a
translationally-symmetric situation).

Penrose-Onsager criterion for BEC
We are now ready to state the Penrose-Onsager criterion for
condensation, first formulated in 1956: if at any given time  it is
possible to find a complete orthonormal basis of states of
 such that one and only one of these states has an
eigenvalue of order  (the rest being of order 1), then we say the
system exhibits BEC. One should note that this definition only applies
to "simple" BEC, as opposed to the "fragmented" case (of no
concern to us here), where two or more of the eigenvalues of the
one-body density matrix are of order .

Condensate wavefunction and the order parameter
We denote the single macroscopic eigenvalue of the density matrix by
, and its corresponding eigenfunction by
.  is called the condensate
  wavefunction and the  particles occupying it the
condensate, while the ratio  is the condensate
  fraction. It is not necessarily true that , even at zero
temperature. Also note that, while  behaves as a
single-particle Schrodinger wavefunction, it is generally not an
eigenfunction of the single-particle part of the Hamiltonian, or of
any other simple operator for that matter, other than
. Another useful quantity frequently found in the
literature is the so-called order parameter,
. We see that, while
 is normalized to 1,  will be normalized to .

No-go theorem for lower dimensionality
It is worth mentioning the existence of a theorem due to
Hohenberg , stating that, in the thermodynamic
limit, BEC cannot occur at a finite temperature in any system moving
freely in space in less than three dimensions, irespective of the
existence and/or sign of the interparticle interactions, as thermal
fluctuations would destroy the condensate. Note that this theorem,
however, only applies under equilibrium conditions, the nonequilibrium
case still being an open question.  Furthermore, there is no general
proof that a realistic system of interacting bosonic particles must
show BEC, even at zero temperature - the solid phase of He
constitutes an obvious counter-example.

Experimental proof
Most gases, with the notable exception of He, are solids at the
densities and temperatures predicted by Eq. eq:Tc3D. That is
why it took no less than 70 years between Einstein's original paper
and the first experimental observation of BEC in an atomic gas. In
1995, the group of Eric Cornell and Carl Wieman succesfully condensed
a cloud of Rb atoms  (closely followed by
the group of Wolfgang Ketterle at MIT with Na
atoms ), by first bringing the system to a
very low density, and then cooling it fast enough to prevent any
recombination processes that would have lead to the formation of the
solid phase. While other odd-isotope alkali elements, especially
Na or Li, are also routinely used in experiments, the
first non-alkali atom to be cooled into the BEC phase was hydrogen
H. Due to the extreme diluteness of these systems
( atoms/cm), the typical range of  is from
 nK to a few K. Achieving such ultra-low temperatures
stimulated the development of novel experimental techniques, such as
magnetic/laser trapping and evaporative cooling of atoms.

Diagnostic techniques
The very low densities of alkali gases also limit the range of
available diagnostic techniques. The most commonly employed method in
BEC experiments is optical absorption imaging, where one shines a
laser on the gas and detects the percentage of transmitted power. The
image is usually taken after removing the trap and allowing the gas to
expand. This gives information about the gas density as a function of
coordinates and time, with a spatial resolution of a few m. In
stark contrast to liquid He, density-related information seems
to be sufficient for most practical purposes.

Diluteness/weak interaction condition


Gross-Pitaevskii equation

Short history and utility of GP equation
As it turns out, many of the experimental results in ultracold gases
can be interpreted on the basis of a single equation for the
condensate wavefunction . This equation, first derived
in 1961 independently by Eugene Gross and Lev Pitaevskii, was
originally intended as a phenomenological description of quantum
vortices in the superfluid phase of liquid He, below the lambda
point. Since liquid helium is a strongly interacting system however,
the GP equation turned out to be much better suited to alkali gases.
Before giving the concrete formulation of the GP equation, we must
first explore the nature of the inter-atomic interactions.

s-wave scattering length
In dilute systems, the inter-atomic distance  is on
the order of 1000 , while the range  of the inter-atomic
potential, namely the extent of the last bound state of the van der
Waals interaction, is about 50-100 . As , the
probability of three-atom colissions is substantially diminished. This
justifies limiting ourselves to a binary (instead of three-body or
more) scattering problem: consider two atoms, separated by a relative
distance  and interacting in three dimensions through a potential
. We can therefore decouple their center-of-mass motion from
their relative one and write a Schrodinger equation for the
scattering states .

For temperatures below , the thermal de Broglie wavelength
 (as mentioned in Sec. ), meaning all
significantly occupied states will have small wavevectors,
. This directly translates to a low relative kinetic
energy, and hence small relative wavevectors, for the scattering
problem outlined above. However, we know from scattering theory that
the probability for two atoms, with relative angular momentum ,
of being separated by a distance  is proportional to
, therefore essentially negligible in the limit
. That is of course, unless , meaning their
relative state is -wave, which is what we will assume from now
on. Since , one can use the asymptotic expression for
, which only depends on the scattering amplitude. At small
wavevectors, this amplitude can be safely replaced by the
-wave scattering length , which will encapsulate all
the interaction effects on macroscopic properties of the atomic gas.

One can now replace the two-body potential  with an effective
interaction, , provided it gives the same
scattering length. The limit of small wavevectors prompts us to only
consider the lowest Fourier component of , equivalent
in real space to a contact interaction(Technically, one
  should also include a regularizing part in order to remove any 
  divergencies of the wavefunction.)
, where we have introduced the
interaction coupling constant , whose value can be calculated using
the first-order Born approximation  



The scattering length  therefore becomes the small parameter of
the theory of weakly-interacting ultracold gases, and the validity of
the Born approximation rests on the following two conditions


  k a_s & 1

    a_s & ^-1/3
Eq. eq:diluteness-condition is called the "diluteness
condition", and it paves the way to various mean-field approaches,
such as the GP equation.


Derivation of Gross Pitaevskii equation
Formally, the GP equation corresponds to the lowest-order expansion in
 of the more exact Bogoliubov theory. However, we will try to
give a hand-waving justification of it for the zero-temperature
case. At , all  particles are in the condensate, therefore one
could neglect all inter-particle correlations and introduce the
simplest (Hartree-Fock) ansatz, expressing the ground state many-body
wavefunction  in the symmetrized form



As mentioned in Sec. , the single-particle state 
(now occupied by all the bosons) obeys a Schrodinger-like
equation, to which we must add the energy of the effective binary
interactions. In mean-field, these interactions contribute the
equivalent of a one-particle potential term proportional to
.  Together with the kinetic part,
this results in the nonlinear equation(We have tacitly
  assumed that  is large enough, such that .)



where we have also included an external potential , normally
used to model harmonic trapping of the gas. Note that
Eq. eq:TDGP is valid for physics occuring over distances much
larger than the scattering length , which in turn must be smaller
than the typical range  of the potential .

Caveats of TDGP
Eq. eq:TDGP is the time-dependent Gross-Pitaevskii (TDGP)
equation, a mean-field result where the condensate wavefunction
 must be calculated self-consistently. It is important to
emphasize that the TDGP equation is also valid at nonzero temperatures
, provided that the density of non-condensed particles is
much smaller than the condensate density. In that case, the condensate
number  is smaller (but still on the order of) the total particle
number . Finally, one must note that the nonlinearity of
Eq. eq:TDGP builds a bridge connecting BEC to nonlinear
optics, where a similar relation is used, under the name of
nonlinear Schrodinger equation.

Time-independent GP
In case the external potential  does not depend explicitly on time,
the stationary solutions of Eq. eq:TDGP evolve with a trivial
phase factor . This yields a time-independent GP
equation for  (we set  from here on)


where  is chemical potential of the gas, the energy required to
add one more particle to the system.(Technically, it is the
  Lagrange multiplier associated to the conservation of particle
  number , and can be shown to be very close to the actual
  chemical potential in the thermodynamic limit. )


Linear response theory

Following loosely the formalism presented in Ref. , we now let  be the steady state solution to the GP equation in
the time-independent trapping potential 



with the GP Hamiltonian defined as



This Hamiltonian describes a bosonic condensate of  particles
with contact interactions quantified by , and chemical potential
.  Now consider adding a small time-dependent perturbation on top
of the trap, giving . We are
interested in the response of the condensate to this perturbation.
For weak perturbations, we can perform a linearization of the GP
equation Eq. eq:GP-atoms around the stationary solution
 - an approach known in the literature as the "linear
response" formalism. The condensate wavefunction 
evolves according to



We assume a small deviation of the wavefunction from its initial
steady state



such that we can expand Eq. eq:GP-atoms-evolution and keep only linear terms in  and . We get



Note that Eq. eq:GP-atoms-lin is not strictly linear due to
the coupling of  to . To restore
linearity, we consider the functions  and
 as being independent and write the linear system



where we have introduced the linear operator



and the source term .  Note that
 is a non-Hermitian operator!

We now consider the eigenvalue equation for the operator 



with  being the right eigenvector and  its
corresponding eigenvalue



Similarly, we also introduce the left eigenvector, obeying
,
and the orthonormality condition
. 

Notice that  and  are
connected by the unitary transformation(Note that this holds
  as long as the Hamiltonian  only contains real
  terms.)



where  is the third Pauli
matrix. We say that  is -Hermitian, meaning
that one can define a new scalar product
 with a
different signature, such that

. The operator  is
usually called the metric operator, and, not suprisingly in our case,
it is the same as the one of the scalar Klein-Gordon equation. A
pseudo-Hermitian operator usually also posesses antilinear symmetries,
and as we will see below this is also the case for
. Interestingly, for operators with a real spectrum,
it can be shown that one can define another metric , which
guarantees a positive-definite inner product, or, in other words,
 (provided  of
course). This can be used to formulate a probabilistic quantum theory
for the new wave-functions  and . For the general
theory and properties of pseudo-Hermitian operators, we point the
interested reader to Ref. .


Using Eq. eq:symmetry-1, we get the general form of the left
eigen-vectors as



with  a normalization factor.  We can chose
 and group the eigenvalues of 
into 3 families, according to the quantity



We therefore have: the "" family, corresponding to , the
"" family, such that  and the "" family, with
.








We are now ready to write the completeness relation



Using Eq. eq:completeness, we can decompose any column vector
as(The modes in the "" family do not appear in this
  expansion as their components live in the space orthogonal to the one
  of our solution.)

  l_1l_2 = _k "+" family [u_kl_1 - v_kl_2]u_kv_k

  + _k "-" family [v_kl_2 - u_kl_1]u_kv_k

There is now a further symmetry of  that we can
exploit in our problem, a sort of time-reversal "spin"-flip
symmetry, namely



where , with
 the first Pauli matrix and 
the complex conjugation antilinear operator. This results in a
duality between the "" family with eigenvectors  and
energy  and the "" family with eigenvectors
 and energy
.

We can now finally project Eq. eq:GP-atoms-system onto the
eigenvectors of . Using the above-mentioned duality and
Eq. eq:decomposition, we get



with the complex amplitudes  satisfying 



where we introduced




Cherenkov emission of Bogoliubov excitations


We now turn to applying the formalism developed in
Sec.  to a concrete physical example, namely
a flowing condensate scattering against a static
defect . The BEC(We integrate all
  density profiles along the  direction, resulting in an effective
  2-dimensional description.) is therefore in a state with
well-defined momentum, described by the plane wave



and a chemical potential . Since we have
no trap, , and Eq. eq:GP-atoms produces the
equation of state



where we have introduced the condensate density
. 

We now introduce a weak perturbation in the form of a static localized
defect potential , which can represent
for instance a laser spot depleting a small area of the condensate, as
shown in Fig. .

    
    
    Density profiles of an expanding BEC hitting a stationary defect
created by the repulsive potential of a blue-detuned laser beam. The
condensate has different speeds in the two panels, moving roughly
twice as fast in the right-panel. Notice the Mach cone formed behind
the defect, which gets narrower as the condesate moves faster. From
Ref. .
    

    
    
    Top panels: Bogoliubov dispersion
Eq. eq:bogoliubov. The dotted lines indicate the 
 plane.
    Middle panels: Locus  of intersection of the 2D
dispersion with the  plane. Green arrows are normal
to , while the dashed lines indicate the Cherenkov
cone.
    Bottom panels: Real-space density modulation, with a
-defect at . Dashed lines show the Mach cone. Left
column panels are for , and right column for 
. From Ref. .
    

Using Eq. eq:atom-MF, the GP Hamiltonian becomes
 and the
source term
. We now
get the linear operator for our problem in the form



Notice that, due to the presence of the off-diagonal exponential
terms,  does not commute with the momentum operator, which is
the generator of the spatial translation group. Luckily, however, we
can restore translational invariance by a simple unitary
transformation, as shown below.

Using the standard commutation relations, one can show that, for a
constant wavevector , the unitary operator(The hat
  symbol denotes operators in the relevant Hilbert space.)



performs a translation in momentum space,
, with the
ket  representing a single particle state with
wavevector  such that
. Using the
definitions above, one can easily obtain the commutator



This allows us to rewrite the following expressions
  
    T^(k_0)kT(k_0)& = k - k_0I

    T(k_0)kT^(k_0)& = k + k_0I  
  
We now recognize the two exponentials in Eq. eq:ourL as being
the real-space representation of  and its
hermitian conjugate. This motivates us to define the following unitary
operator



such that a unitary transformation of our operator  now restores
translational symmetry. Indeed, one can see that



where we have made use of Eqs. eq:products and we have written
 in a base-independent representation.  In the subspace of
momentum eigenstates , we can write the (right-)eigenvalue
equation corresponding to Eq. eq:translated-L as



where we have recovered the matrix representation of
Eq. eq:LGP, and introduced the notation



Here  labels the 2 different eigenmodes, and we defined
the condensate speed .

Notice that the  mode has only one eigenvector. However, one can
safely exclude it as this mode does not imply energy or momentum
transport. Excluding the  point, one can then solve
Eq. eq:right-eigen, obtaining the celebrated Bogoliubov
excitation spectrum



with , as before. Here  represents the momentum of
the quasiparticle excitation with respect to the momentum  of the
condensate. Note that the complex amplitudes  and
 only depend on the absolute value of , while the
(real) spectrum of Eq. eq:translated-L, ,
is the Bogoliubov spectrum with an additional Galilean boost
.

We can now further simplify the problem. As can be seen from
Eq. eq:symmetry-2, the 2 eigen-families  and 
are linked by a duality, stemming from the 
symmetry(This can be actually formally proven after defining
  the parity and time-reversal operators corresponding to our
  problem. For details, see Ref. .) of the
Bogoliubov operator .  We therefore drop the subscript 
and make the convention that
.

Bogoliubov spectrum discussion
The Bogoliubov spectrum Eq. eq:bogoliubov is shown in the top
row of Fig. , for two different values of
, which sets the slope of the dotted lines (indicating the
 plane). In the non-interacting case  and the
spectrum reduces to a simple parabola characterising a free
particle. For repulsive interactions,  and we can distinguish
two qualitatively different domains, after first introducing the
typical length scale of the problem, called the healing
  length. The healing length  is a measure of the distance over
which the condensate density recovers its equilibrium value 
when forced to vary away from this value. For example, in a box the
boundary conditions fix the density to zero at the positions of the
walls. In mathematical terms,




We first explore the domain of small momenta, , which is
characterised by a linear behaviour of the dispersion,
, that implies the propagation of low-energy
excitations in the form of sound waves, with a
velocity(The sound velocity here is measured in the
  condensate rest-frame ().)  given by



The Landau criterion for superfluidity determines the maximum velocity at which a weak impurity can travel
through the condensate without dissipating energy. In order for it to
dissipate energy, such an impurity must be able to create
quasiparticle excitations in the condensate. Conservation of energy
and momentum then results in a critical velocity



below which no dissipation can occur. In our case, this velocity is
precisely equal to the speed of sound, . Furthermore, the
two situations, the one of a particle moving through the condensate,
or of the condensate moving against a fixed defect, are physically
equivalent, being connected by a Galilean transformation. We can
therefore conclude that we must have  in order to
observe any propagating perturbation, otherwise for  the
superfluid will remain unperturbed. Before moving on, it is worth
noting that the Landau criterion has some asociated caveats.  One is
assuming that the only excitations are density excitations,
phonons. Thus one is, for example, neglecting the nucleation of
vortices by a macroscopic defect with a size comparable to the healing
length, which would lower the effective critical velocity. Vortices
would furthermore also break the translational invariance along the
transverse directions, an invariance that we already made use of. The
second caveat is that quantum fluctuations are also neglected. As seen
in Ref. , they could lead to nonzero
dissipation even at sub-sonic speeds.

The second domain of interest is the one of large momenta,
. Looking at Eq. eq:right-eigen, one notices that
the only -dependent parts of  are the
diagonal terms. These terms have opposite sign, so the off-diagonal
coupling between  and  becomes highly off-resonant at
large . Completely neglecting it gives the free-particle-like
shifted parabola , with
 and .

The discontinuity in the spectrum at zero momentum can be explained by
the fact that the diagonal and off-diagonal parts of 
 are equal in absolute value.

In order to obtain the defect-induced density perturbation, we must
now also determine the eigenvectors of the problem. Making use of
Eq. eq:symmetry-1, we can act with  on the
eigenstates of  to obtain the ones of
. This finally leads us to a biorthonormal
basis

, containing 4 basis vectors



which fulfill the orthonormality condition



and the completeness relation



provided of course that we normalize in such a way that
. In this basis, the spectral
decomposition of Eq. eq:translated-L is the diagonal form



The concrete form of , coupled with the
normalization condition Eq. eq:norms, determines the
eigenvectors of the "+" family up to a phase factor. Indeed, one can
choose , with 



Furthermore, in case the Hamiltonian doesn't contain any time-reversal
symmetry-breaking terms, one can chose  and  to be real
quantities, without loss of generality.

It is now straightforward to solve the linearized evolution equation
Eq. eq:GP-atoms-system. In particular, for a localized static
defect potential , quasiparticle modes
at all wavevectors  are excited.(If one wishes to
  selectively excite a pair of modes, one can use a periodic
  potential, following, for example, Ref. .) The
source term is time-independent, and hence the quasi-particle
amplitudes of Eq. eq:amplitudes-bk have the simple form
. Equivalently, we can obtain the
defect-induced perturbation of the wavefunction from its initial
steady state by directly inverting Eq. eq:bidecomposition and
then reversing the unitary transformation that was applied to obtain
Eq. eq:translated-L. Whichever route we take, it is clear that
the final answer will have a resonant structure, containing
 in the denominator, hence the dominant modes will be the
ones that satisfy .

Landau causality rule
We must also mention the subject of adiabatic switching at
this point. Adiabatic switching is neccesary to causally distinguish
the past from the future, making sure that the time  is
prior to that of the occurence of any cause (in our case the defect)
giving rise to the effect (pertubation of the condensate density). The
way to achieve this in practice is by shifting the real poles of the
Bogoliubov dispersion Eq. eq:bogoliubov into the lower half of
the complex plane by an infinitesimal amount,
. This corresponds to a
weak damping of the plane wave solution and ensures that no Bogoliubov
excitations were present at .

Geometrical construction
We show the defect-induced perturbation of the condensate density in
the bottom row of Fig. , for two distinct
values of the condensate speed . Rather than giving its full
analytical expression, it is more instructive to present a geometrical
construction, detailed in Ref. , that can shed
light on the main features of the condensate response.  We have
already identified the importance of the poles of ; the
solutions of  can be visualized if one
plots the intersection of the Bogoliubov dispersion surface
 with the  plane. The locus of this
intersection is a closed curve that we will denote by  and
that is plotted in the middle row of
Fig. . Making the connection with the Landau
criterion presented earlier, we can identify two regimes.  For small
velocities  we have the sub-sonic regime, where
 is a single point at the origin. Consequently, one can
observe superfluid-like behaviour with no propagating density
modulation. The modulation will stay localized in the vecinity of the
defect, and, in the Galilean-equivalent problem of a particle moving
through the condensate, it would renormalize the particle
mass. 
As we gradually increases the condensate speed, and hence the slope of
the  plane, this plane will touch the surface of the
dispersion relation when , marking the entry into the
super-sonic (dissipative) regime. Further increasing the
speed will increase the size of , as can be seen in
Fig. . The green arrows, orthogonal to the
curve  at each point, represent the group velocity of the
Bogoliubov mode at that particular -value. They show the direction
of propagation of the density perturbation away from the defect, up to
infinity. It is the interference of these propagating modes that we
see as the real-space density pattern. Note that the locus 
has two distinct regions, inherited from the linear and quadratic
domains of the dispersion Eq. eq:bogoliubov.

The linear region of , close to the origin, is characterised
by essentially the same physics as the Cherenkov effect in
non-dispersive media . The Charenkov effect
consists in the emission of electromagnetic radiation by a charged
particle moving relativistically through a dielectric medium at a
velocity higher than the (phase) velocity of light in that medium. The
emission is concentrated into a Cherenkov cone in momentum
space, of aperture , where
.  The higher the particle
speed  with respect to the speed of light , the wider the angle
 between the emission and the direction of motion. The
Cherenkov cone is depicted by the dashed lines in the middle panels of
Fig. , and the angle  in our case is
of course given by . Due to the momentum-space
singularity at the origin, the group velocity has a jump, defining a
whole region of space where no phonons are emitted. This corresponds
in real-space to the so-called Mach cone, named in analogy to
the cone that is created by a super-sonic aircraft (since we are
dealing with terminology, it is worth noting that the ratio 
goes by the name of Mach number). The Mach cone is depicted
by dashed lines in the bottom panels of Fig. :
as its aperture  is quantified by , that
means that, the faster the fluid, the narrower the cone will be.

The high-momentum, rounded region of , which corresponds to
the quadratic single-particle-like dispersion, has no equivalent in
Cherenkov physics. It is responsible for the hyperbolic-like
wavefronts(The wavefronts are actually parabolic in the
  non-interacting limit of .), emitted in the positive
 direction (upstream). The physical origin of these
rounded waves lies in the interference between the coherent matter
wave of the BEC and the wave scattered off the defect. It it worth
noting that the curve  also helps one determine the spacing
between the emitted wavefronts, which is inversely proportional to the
value of the momentum at the particular point on  that
corresponds to the wave propagation direction. In practice, we see for
example that the spacing along , in the positive 
direction, is wider in the left-bottom panel than in the right one.

Aside from being just a useful geometrical construction, the locus
 can actually be observed in scattering experiments, both in
the context of atomic condensates, as well as for microcavity
polaritons (Chapter ), where it takes
the name of Rayleigh scattering ring.

Before closing this Section, it is worth making the connection between
the Cherenkov waves shown in the bottom panels of
Fig.  and the more mundane example of surface
waves created at the interface between a fluid layer and a gas. We
first need to define the capillary length
, with  the surface
tension of the fluid-gas interface,  the gravitational constant,
and  the fluid density. Now, if the height  of the fluid
layer in question is smaller than , the fluid
dispersion is dominated by capillary effects, and its form is similar
to Eq. eq:bogoliubov (see Ref. ). But we
have just seen that the form of the dispersion (more precisely, its
intersection with the  plane) determines the shape of the
ougoing waves. This means that, for very shallow fluids, we will
observe similar wavefronts to the ones of
Fig. . More concretely, for the water-air
interface, a thickness of  millimeter and a speed of 
cm/s should do the trick!

Superfluidity

Short history
The history of superfluidity is intrinsically linked to that of helium,
the liquefaction of which was first achieved in 1908 by
Kammerlingh-Onnes. This ushered in a new era of low-temperature
physics, and it was soon realized that the thermal expansion
coefficient of liquid He had a discontinuity at the so-called
lambda-(-)point. This lead to the classification of liquid
helium in two distinct phases, above (He-I) and below (He-II) the
critical temperature . It was in this context
that, 30 years later, Kapitza in Moscow and Allen and Misener in
Cambridge investigated the flow of He-II through a narrow opening
between two large containers. Observing that the fluid essentially
behaves as having no viscosity(In fact, formally one cannot
  even define viscosity in this case, as the ratio of the mass flow to
  the pressure differential diverges.), Kapitza introduced the term
"superfluidity" in analogy to the superconductivity observed in
charged systems. It turned out to be an inspired analogy, as the
prevalent view of the scientific community nowadays is that the two
are essentially the same phenomenon, stemming from BEC (in the neutral
case), respectively Cooper-pairing in the charged case (arguably a
kind of pseudo-BEC). The connection between superfluidity and BEC was
first made by Fritz London, who noted that  was not too
far off from the critical temperature one would predict by applying
Eq. eq:Tc3D to a noninteracting gas of helium atoms. Of
course, liquid helium is anything but a noninteracting gas, so in that
sense it is rather amusing that this connection was made. The modern
conception of superfluidity is that it represents a collection of
phenomena, most of them related to flow properties, and the initial
observations of Kapitza can in fact be broken down into such
"elementary" ingredients, as we shall see.

Two-fluid model
A great leap in the understanding of superfluidity was made with the
help of Lev Landau's phenomenological two-fluid
  model . In a nutshell, one can think of the
problem in terms of two liquids, one formed by the condesate occupying
the single-particle state  and flowing without friction, and
the other behaving as a normal fluid. To be fair, it must be said that
Landau actually opposed the whole notion of BEC, and formulated his
model without it. We must also note that the superfluid and normal
components in He-II are not the condensed and noncondensed atoms, as
the whole liquid is superfluid at zero temperature, while only around
10 of the atoms are in the BEC state. That being said, the model
describes He-II in thermodynamic equilibrium by introducing two
independent velocities ( and ), and associated
densities ( and ), which are
temperature-dependent. One can then express the total mass current
density  and kinetic energy of flow
 as 
  J() & = _s()v_s() + _n()v_n()

  Q() & = _s()v_s^2()2 + _n()v_n^2()2
with . We can futher define the
superfluid () and normal () fractions as
 and , with
. Furthermore, we have the limits

)
  f_s(T) &[T 0]) 1

  f_s(T) &[T T_)]) 0

It is in this sense that one is dealing with an intuitive picture of
He-II as a mixture of two independent, interpenetrating components,
each associated with its own mass current and flow energy.


Superfluid velocity)
We now look at the  atoms condensed in the state , in
order to investigate the origin of the superfluid velocity
. One can express the condesate wavefunction in the Madelung
form(We explicitly re-introduced  in the equations of
this Section to show their quantum nature.)



with  being the particle
density(We again emphasize that, in general, the superfluid
  density  is not equal to the density of condensed
  atoms.)  and  the phase (in dimensions of an action), both
real functions.

The quantum mechanical definition of the probability current density
of the condensate yields (see Appendix )

  j_s(,t) & = N_0(t)2im[_0^(,t)_0(,t) - _0(,t)_0^(,t)]

& = N_0(t) _0(,t)^2 1m (,t)

with , as before, representing the particle mass. The ratio
 will then be the condensate velocity
(called superfluid velocity in the literature for the
historical reasons detailed above)




Onsager-Feynman quantization condition
As the curl of a gradient is zero, the flow given by
Eq. eq:supervelocity is irrotational, in other words the
vorticity



provided, of course, that  is defined in all space (i.e.
the magnitude of  is finite). If that is not the case, one can
consider a contour  around the region where  vanishes, and
the condition that  be single-valued then gives the famous
Onsager-Feynman quantization



with  being the windinding
  number, the number of turns made by the phase in the complex plane,
as we go around the integration contour . Thus the circulation of
the irrotational flow is quantized, with the quantum of circulation
being  cms.

Two effects intro
We are now in the position to discuss the two elementary ingredients
which together account for the original obsevation of superfluidity by
Kapitza. Following Ref. , we consider a multiply
connected geometry in the form of a hollow cylindrical pipe with solid
walls shown in Fig. .

    
    
    The multiply connected toroidal geometry, of radius , used to
    showcase superfluid behaviour. The anullus, which contains He
    atoms, rotates at an angular velocity .
    
  
Assuming the pipe of radius  is filled with liquid helium, and
rotates with angular velocity , one can then obtain the
(temperature-dependent) orbital angular momentum of the system 
by minimizing the effective free energy 


Apart from the already defined quantities, we have introduced the
moment of inertia , and defined the quantum unit of
angular velocity, . Finally,  is the
nearest integer to , expressed as



One needs to compare Eq. eq:HF-angularmom with the angular
momentum of a normal liquid under the same circumstances, which is
given by . Keeping this in
mind, we can now introduce the two qualitatively different
experiments.


Persistent currents



The first experiment starts by rotating the torus at a large angular
velocity , on the high-temperature side of the
-point. Once the fluid equilibrates with the motion of the
container walls, we then cool the system through ,
maintaining the same angular velocity . Finally, we stop the
container from rotating, and measure the angular momentum.



After a short while, proportional to the viscosity of the normal
component, its velocity  becomes zero. However, the
superfluid part of the liquid experiences no friction with the pipe
walls, and will therefore continue circulating indefinately! Since we
are not dealing with the ground state of the system, but with an
excited state with astronomical lifetime, we call this phenomenon
metastability of supercurrents ("persistent currents" is
another commonly found term).



However, metastability of superflow is not a direct consequence of the
BEC state - the underlying mechanism is one of topological origin:
consider the contour  in Eq. eq:feynman to run around the
torus of Fig. . The current in the ring will be
proportional to the circulation Eq. eq:feynman. Since the
latter is quantized, the windinding number  is conserved, hence
so is the (super)current!



Similar arguments explain the existence of so-called topological
defects in superfluid systems. Such defects are singularities of the
condensate wavefunction, for instance, vortex lines (or rings, if the
line closes in on itself) associated with the superfluid
component. The conservation laws presented above make the vortex a
highly stable configuration, as opposed to the case of a normal fluid,
where it is short-lived. In this sense, the circulating currents which
constitute the vortex are "persistent", corresponding to metastable
supercurrents. Furthermoe, a superfluid system cannot sustain bulk
vorticity due to Eq. eq:vorticity, and the vortices are
quantized according to Eq. eq:feynman.

In the final state with the container at rest, , so the
first term of Eq. eq:HF-angularmom dissapears. The second term
survives, as the system cannot change its winding number. In the limit
of large  the angular momentum will then be approximately given
by . We see that one can reversibly
increase (or decrease) the angular momentum  by varying the
temperature (provided, of course, that it is always kept below
).



Since the superfluid circulates despite any roughness of the container
walls, we can reasonably expect that an object moving though a
stationary fluid also experiences reduced friction. This led to a
series of experiments with charged ions being accelerated at various
speeds by an applied external electric field 
. It was seen that, in accordance with the Landau
criterion presented in Sec. , there exists
a critical speed below which the ions experience reduced viscosity
(and actually no viscosity at all in the  limit).



HF effect
We now turn to the second experiment, which is arguably one of the
fundamental defining properties of superfluidity. The effect was first
predicted by Fritz London and then observed by G. Hess and W. Fairbank
in 1967 in the context of liquid helium, hence the name Hess-Fairbank
(HF) effect.



As before, we cool the sample while rotating with a lower angular
velocity . We wait for the system to equilibrate,
and then measure the angular momentum, this time without stopping the
rotation of the container.



The value of  will be given by Eq. eq:HF-angularmom,
instead of the classical  (that is why the HF effect is also
called nonclassical rotational inertia, NCRI).

We can now distinguish two different regimes. The first one, for
, implies , therefore
, so the superfluid component contributes
nothing to the current. At ,  and we get
the spectacular result that the liquid completely ceases to rotate
along with the container!

In the general case, the second term of Eq. eq:HF-angularmom
survives. If we again consider the zero-temperature limit, and plot
the angular momentum  versus angular velocity
, we get a series of plateaus of
width  at positions
, and so on. The centers
of these plateaus all fall on the "classical" line
.



A hand-waving explanation, for the simplest case of a noninteracting
atomic gas, goes as follows. In a noncondensed system rotating with
angular velocity , the effective single-particle energies are
shifted by an amount , with  the
angular momentum quantum number of state . The Bose-Einstein
function then redistributes the particles to accomodate this shift,
producing a "classical" angular momentum. In a BEC however, all
 atoms must be in the same state , and posess the same value
of . Their only allowed contribution to  is then
. For small , the condensate is in the
 state, and so it does not rotate with the container! This
argument can be formalized, so one can show the HF effect to be
impossible inside the framework of classical statistical mechanics,
much like the Bohr-Van Leeuwen theorem shows magnetism to be a purely
quantum phenomenon. Note that, while the HF effect can be (at least
hand-wavingly) explained without including interactions, these are
crucial to any theory of metastable currents.


Meissner effect
As already explained in the introductory paragraph of this Section,
there is a deep underlying connection between superfluidity in a
neutral system and superconductivity in a charged one. In fact, the
analog of the HF effect in a charged system is the well-known Meissner
effect.

We now consider the same topology of Fig. , and replace
the container by a crystal lattice and the liquid helium (or atomic
gas for that matter) by the electrons inside this lattice. If we now
keep the crystal at rest and produce a tangential vector potential
, this will generate a flux
 inside the conducting ring
(usually called an "Aharonov-Bohm"
flux). The definition of the "superfluid"
velocity Eq. eq:supervelocity now needs to be revised
according to the gauge-coupling recipe, which consists of replacing
the canonical momentum  by the kinematic one:
. As a result, we get



where the factors of 2 appear because the macroscopic wavefunction is
now that of the center of mass of a Cooper-pair of electrons, so 
and  are both multiplied by 2. The resulting Onsager-Feynman
quantization Eq. eq:feynman then becomes



where  is the superconducting flux
  quantum. If the (mass) supercurrent is given, as before, by
, then this new quantization condition
results in a non-zero current even for . Indeed, for
 (compare to the regime
 discussed for the HF effect), the
 state is realized and the electrical current  is
given by



This proves that, while in a normal conductor the flux  would
not generate any circulating current, a superconducting ring would
respond with a diamagnetic (note the minus sign in
Eq. eq:diamagnet) electrical current proportional(The
  proportionality constant is called the susceptibility.) to the
vector potential . Coupling this with Maxwell's equations
immediately results in the well-known Meissner effect, namely the
expulsion of an applied magnetic field inside a superconducting
sample.


Synthetic gauge field
We now proceed to a closer analysis of the connection between the
charged and neutral case. The Hamiltonian  for a particle of charge
 and mass , moving in a vector potential , and in the
presence of an external potential , is given by
Eq. eq:charged-ham. In our problem,  is imposed by the
toroidal geometry of the container, which is for now assumed to be
stationary in the inertial frame of the stars. On the other hand, the
Hamiltonian , of a neutral particle with the same mass,
but this time in the frame that rotates with the container at an
angular velocity  is shown in
Eq. eq:neutral-ham,


    H  = 12m(p - qA())^2 & + V()

    H^  = 12m(p^-m    ^)^2 & + V(^) - 12m(    ^)^2
with  being the coordinate (and  the
canonical momentum) in the rotating frame. Making use of the
cylindrical symmetry, one can decompose  as ,
and therefore .

We now see that a neutral system in the rotating frame, apart from a
centrifugal term ,
is formally identical to a charged system in the rest frame, with the
correspondence 
. This connection can be used to make a neutral particle
experience an "effective" electromagnetic vector potential. Its
gauge-coupling to this potential is what eventually led to the
celebrated concept of synthetic gauge field.

In fact, such artificial, rotation-induced, gauge fields have already
been proposed in the context of ultracold atomic gases, in order to
asess their normal and superfluid
fractions .


BEC Experiments

    
    
    Scanning a blue-detuned laser beam with velocity , though the
    cigar-shaped atomic condensate (left panel), creates a density hole
    (right panel) larger than the healing length. The friction acting on
    this hole can then be used to asess the Landau critical velocity. From
    Ref. .
    
  
This naturally leads us to the topic of experimental evidence for
superfluid behaviour. The answer, however, will not be a clear-cut
one, as we have seen that it actually depends on which aspects of
superfluidity one is interested in, such as NCRI, persistent currents,
quantized vortices or reduced friction on impurities. Historically,
liquid He-II was the default system for exploring such physics,
however, new systems such as ultracold gases and, more recently,
microcavity polaritons (the subject of Chapter )
have also entered the picture. In the case of liquid helium,
superfluidity is relatively easy to detect, while BEC can only be
asessed through circumstantial evidence (such as neutron scattering),
as strong interactions oppose any large density changes. Ultracold
gases, on the other hand, present a spectacular BEC onset that can be
observed in the density profile, while superfluidity proves somewhat
more mysterious.

Regarding metastable supercurrents in a simply connected geometry,
Chevy and coworkers  succeded in stirring an atomic
condensate, and then measuring its angular momentum using the
precession of certain collective modes.  Other studies were also shown
to produce a metastable single-vortex state (and even vortex lattices)
that persist, despite not being the ground state of the system (in
fact, it can be shown that the atomic BEC ground state is
characterised by the absence of any macroscopic current flow).
Interferometric measurements are commonly employed to see the
circulation around these vortex cores.  However, the lifetime of the
vortices is limited by how long it takes before they escape the
magnetic trap. To circumvent this limitation, persistent currents in
multiply connected geometries (such as toroidal traps) were also
investigated, for instance in Refs. and .

As far as the Landau criterion is concerned, an experiment similar to
the one accelerating charged ions through liquid helium was reported
in Ref. . The authors used a Raman laser to
transfer a fraction of the condensed atoms into a different
(hyperfine) state. As a result, those atoms no longer experienced any
magnetic trapping, but only the effect of gravity. The mobility of the
free-falling atoms was then computed, revealing information about the
friction  acting on them.  A fairly sharp
transition with a sudden increase in friction was observed, at
velocities close to the local speed of sound.

In the opposite limit of macroscopic "impurities" (large compared to
the healing length), Ref.  details the use of a
blue-detuned laser beam, in order to create a density hole which is
then moved around the atomic condensate (see Fig. ).  By
measuring the density difference immediately ahead and behind the
hole, one can deduce the frictional force acting on it (see
Appendix ). This experiment resulted in a
critical velocity about ten times smaller than the speed of sound,
presumably due to vortex creation by the macroscopic defect.


Optical lattices

Laser trapping
An electric dipole moment  placed in an external
electric field  aquires an energy
. If
an atom is placed in a light field, the oscillating electric field of
the latter induces an electric dipole moment
 on the atom,
with  the free-space permittivity and  the atomic
polarizability . Integrating, we obtain the
interaction energy of the atom with the field as

. A classical model of the atom as a conducting
sphere of radius  yields a polarizability of . In
practice however, atomic transitions complicate this simple picture,
as the polarizability  depends on the driving frequency of the
light field. In the framework of the Drude
model , one has



with the resonant frequency , damping time  and
 the static polarizability. The frequency-dependent
interaction energy  thus obtained is sometimes
called the ac Stark effect . If
the light frequency is smaller than that of the atomic resonance
, the induced dipole is in phase with the electric
field, as , and the force on the atom will
point in the direction of the square-field maxima (the situation is of
course reversed if instead ). This effect was first
used in Ref.  to trap a BEC of sodium atoms
using a focused infrared laser beam (the dominant transition of Na is
the familiar  yellow emission). Interestingly, the
same physics can be used to create an optical lattice, as explained
below.

Optical lattice
Consider a (plane wave) laser beam propagating along the  axis,
with an associated electric field
,
which reflects from a mirror situated at . The interference
between the reflected and original beam produces a standing wave
,
with an intensity

. This results in an effective periodic lattice for
the atoms, with lattice spacing , where
. Such a one-dimensional lattice was first
demonstrated in Ref. . By tuning the strength of
the laser field, one can adjust the lattice depth, and by changing the
laser wavelength, one can modify the lattice constant. Adding more
laser beams can produce two- or even three-dimensional lattices, in a
variety of different geometries. For example, the simplest way to
create a 2D square lattice would be to superimpose a second lattice in
the  direction, producing an optical potential of depth
:



The beam along  should be orthogonally polarized with respect to
the one along , in order to avoid interference effects.

Band structure
The eigen-problem of a particle moving in a periodic potential given
by Eq. eq:square-potential is well-known in the context of
solid-state band theory , as the electrons
inside a solid also move in such a potential, produced in that case by
the crystal ions. The energy eigenstates are the Bloch waves
, where
 is the quasimomentum, restricted to the first
Brillouin zone (BZ)

 and 
represents the band index. The functions
 are periodic in both the  and
 direction, with period . To see how these bands
arise, one can expand both the wavefunction  as
well as the potential  in a Fourier series, making use of
their periodicity. Inserting this expansion into the time-independent
Schrodinger equation, we get a linear system of equations for the
Fourier coefficients. For fixed , we obtain
different(Truncating the Fourier expansion to a finite number
  of terms  results in  bands.) eigenenergies,
. These eigenenergies form bands (called Bloch
  bands) when  is varied inside the BZ. Each band as a certain
width in energy (unsurprisingly called the bandwidth) and consecutive
bands are separated by gaps, meaning certain energies are not
allowed. Furthermore, each eigenenergy has an associated
eigenfunction, completely determined by the set of Fourier
coefficients with index . The solution of course depends on
the value of , and we can distinguish two limiting cases,
depending on the ratio of the bandwidths to the energy gaps. In the
weak-potential limit, the bandgaps are much smaller than the
bandwidths, and the eigenenergies have a pronounced dependence on
quasimomentum. The opposite limit, which is the one we want to focus
on, is the regime of deep lattices, also called the
tight-binding (TB) limit. The band dispersion is now flatter
compared to the previous case, and the bandwidth/gap ratio is much
smaller.

Tight binding
For low temperatures in the TB limit, only the lowest energy band
 is occupied, so we subsequently drop the band
index. Since the potential wells are so deep, the lowest band is
adequately described by considering localized wavefunctions at each
well, i.e. at the minima of Eq. eq:square-potential. Suppose
now that a function  exists such that the Bloch waves are
given by



where  is a normalization factor and
 are the well
positions. Eq. eq:wanier represents a wave,
, whose phase and amplitude at each well is
given by the localized function , called the Wannier
  function . It can be shown that Wannier
functions centered at different wells are orthogonal to each other and
form a complete basis. The Wannier functions are thus a unitary
transformation of the Bloch functions, suitable for the case when
particle dynamics is described by inter-well tunneling (hopping)
between nearest-neighbour wells. 

In the case of a BEC, when the number of ultracold atoms per lattice
well is small, the "discrete structure" of the condensate becomes
relevant and we must replace the GP description Eq. eq:TDGP we
have used so far with a lattice model, such as the Bose-Hubbard
model , as suggested in
Ref. . For a condensate confined in the optical
potential Eq. eq:square-potential, with an additional harmonic
trap , the spinless Bose-Hubbard Hamiltonian in 2D can be
written as 


where  is the on-site interaction strength and
 is the energy offset due to the
harmonic trap. Futhermore,
 denotes the
number operator for site  and 
() is the bosonic creation (destruction) operator for
that site. Finally, the kinetic part



can be diagonalized to yield the Bloch dispersion



where the tunneling amplitude  is an exponentially decreasing
function of  . The characteristic tunneling
energy scale is given by the bandwidth, which in this case is .

Hamiltonian Eq. eq:BH-ham supports a zero-temperature quantum
phase transition between a superfluid and a Mott-insulating phase, at
a critical value of the ratio . This transition was
experimentally observed for the case of a bosonic BEC in a 3D optical
lattice, as reported in Ref. . Increasing the
lattice depth , the authors went from a gapless, superfluid
regime dominated by the kinetic term Eq. eq:kinetic-part to a
gapped Mott-insulator state, where the atoms are localized and the
on-site interaction term in Eq. eq:BH-ham is dominant.

Harper-Hofstadter model
Landau levels
In Sec.  we introduced the quantized Bloch
bands, characterising particles in a periodic potential. A similar
quantization also occurs in the case of a charged particle moving
perpendicular to a constant uniform magnetic field, and the allowed
energy bands are called Landau
  levels . In a magnetic field
, the relevant energy and length scales for a
particle with charge  and mass  are set by the cyclotron
  frequency  and magnetic length
, respectively. In this Section we want to
explore the fate of the Bloch bands in presence of a magnetic
field. Alternatively, one could also start from a Landau quantized
system and perturb it via a weak periodic
potential , the two approaches
having distinct validity ranges.

Hofstadter butterfly
One of the simplest TB models for Bloch electrons in uniform magnetic
fields is the Harper-Hofstadter
  model . Its
relevant dimensionless parameter is the ratio of flux  through a
lattice cell to one Dirac flux quantum(As opposed to
  the superconducting flux quantum of Sec. .)




Physically,  reflects(Not to be confused with the
  atomic polarizability, also denoted by .) the
commensurability of the lattice period  and the magnetic length
, as it can be recast in the form
.  The resulting
single-particle energy spectrum of the model depends on the
rationality of . If , for  and  co-prime
integers, the Bloch band splits into  distinct magnetic
subbands. For  irrational, the spectrum consists instead of an
infinite number of energy levels, forming a Cantor set. The union of
all allowed energies as a function of  forms a self-similar
fractal, shown in Fig.  and known as
Hofstadter’s butterfly . In the
continuum limit (, with fixed ), 
and the bands group into clusters that can be identified with the
Landau fan predicted in the absence of a lattice potential.

    
    
    Hofstadter's butterfly, representing the internal structure of the
    lowest Bloch band in a uniform magnetic field. For rational values of 
, the spectrum has  magnetic subbands, which become
    increasingly narrow as  gets larger. Adapted from
    Ref. .
    
  

HH experiments
Solid-state experiments in the Hofstadter regime are notoriously
difficult, since the inter-atomic spacing of a crystal lattice is a
few AAngstroms, and achieving a comparable magnetic length would
require a field of about  Tesla. Signatures of the Hofstadter
bands have instead been seen in artificial superlattices with
periodicities of tens of
nanometers . In the context
of neutral ultracold atoms, synthetic gauge potentials paved the way
for the exploration of similar
physics . For instance,
rotating a harmonically trapped condensate with a high angular
velocity (comparable to the radial trap frequency) allowed BEC
experiments  to enter the lowest Landau
level (LLL) regime. Rotating 2D optical lattices have also been
realized , but the first implementation
of the Harper-Hofstadter (HH) model with ultracold
atoms  employed a different
method  of generating the artificial
magnetic field, namely laser-assisted
  tunneling  based on a linear superlattice
potential . In a nutshell, tunneling
between lattice sites required atoms to absorb and emit photons
produced by external lasers. The interaction with said lasers
generated a net phase change equivalent to the flux of a magnetic
field. The total phase accumulated by an atom when performing a loop
around a lattice cell (plaquette) was  in these
experiments.

Aharonov-Bohm phase
In quantum mechanics, the phase of the wavefunction of a charged
particle moving between positions  and  is changed by
the presence of a magnetic field 
by 


where the integral is taken along an oriented path from point "1" to
point "2". Therefore,  and, if we consider a
smooth closed integration contour  instead of two
separate points, the total phase change (or Aharonov-Bohm
  phase) along  is



where  is the magnetic flux through a smooth surface
 bounded by the curve  with positive
orientation, according to Stoke's theorem.

Peierls phases
On the 2D square lattice in Fig. , one can define
the magnetic phases
 along the
 and  directions in a similar manner, as

    
    
    2D square lattice pierced by a homogeneous static magnetic field
    . The lattice constant is , and the curved
    arrows indicate the field-induced hopping phases.
    
  

  ^x_m,n &= e__m,n^_m+1,n dA() &
  ^y_m,n &= e__m,n^_m,n+1 dA()
where the integrals are now taken along the links joining neighbouring
sites. Using Eq. eq:AB-closed, the flux  through a
plaquette can be related to the total phase 
accumulated on a contour enclosing the plaquette in a counterclockwise
sense



Note that  is gauge-independent, unlike the
Peierls phases .



HH Hamiltonian
To obtain the HH Hamiltonian, we introduce a static uniform magnetic
field into the single-band Bose-Hubbard model Eq. eq:BH-ham of
Sec. . We follow Hofstadter's original
treatment, and neglect on-site interactions(For a treatment
  that includes interactions in the Bogoliubov approximation, see
  Ref. .), as well as the harmonic
trapping. The kinetic term Eq. eq:kinetic-part is then
modified according to Peierls' substitution  to
include the hopping phases 



These phases act to "frustrate" the hopping, such that for
noninteger  it is impossible to minimize the kinetic energy
simultaneously on all lattice links . The
Hamiltonian is no longer invariant under translations by a single
lattice site in any direction, because fixing a gauge for the vector
potential  reduces the physical symmetry of the
lattice. In fact, it can be shown  that a
translation of  is equivalent to a gauge transformation,
suggesting that new translation operators that would commute with
Eq. eq:HH-hamiltonian can be constructed as a combination of
translation and gauge transformation.


Magnetic translation operators
We define the magnetic translation operators (MTOs)
  as



  &= _m,n a^_m+1,na_m,ne^i^x_m,n &
  &= _m,n a^_m,n+1a_m,ne^i^y_m,n

where the phases  are to be determined by imposing
, leading
to 

  ^x_m,n &= ^x_m,n + 2n &
  ^y_m,n &= ^y_m,n - 2m 

While the specific form of the MTOs depends on the choice of gauge,
their commutator is gauge independent



and vanishes only if  is an integer, a situation equivalent to
the zero-field case. In order to proceed, we point out that
Eq. eq:TxTy has a transparent physical interpretation. If we
act on the single-particle state




we see that performing a counterclockwise loop around a unit cell
results in the accumulation of a phase proportional to the flux inside
the cell. This can be generalized to a macrocell containing
 original lattice cells, to deduce that



For  the commutator  vanishes
if  is an integer. The smallest possible macrocell is given
by  and is called the magnetic unit cell. For this
choice of unit cell,  form a complete
set of commuting operators, so we can find simultaneous eigenstates,
labeled by the quasimomentum  which is now restricted to
the magnetic Brillouin zone (MBZ)

.


HH spectrum

    
    
    Top panels: Single-particle energy spectrum (left)
     of the HH Hamiltonian Eq. eq:HH-hamiltonian
    for  and dispersion of the two lowest bands, 
    (middle) and  (right) in (part of) the MBZ.  Bottom
      panel: Energy gap (continuous line) 

    (here  denotes an average over the MBZ)
    and bandwidth of  (dashed line, ) for 
. The dotted vertical line marks the location of the  case
    shown above.
    
  
To obtain the spectrum of , we need to solve its
associated Schrodinger equation,
. We expand the eigenstates
 in the single-particle basis
 and obtain an equation
for the coefficients 

  -J (e^-i^x_m,n_m+1,n+e^i^x_m-1,n_m-1,n

   + e^-i^y_m,n_m,n+1 + e^i^y_m,n-1_m,n-1) = E_m,n

We fix the Landau gauge  for simplicity,
so  and we choose the magnetic unit
cell(In the symmetric gauge, the magnetic unit cell is
  larger, containing  sites.)
, for which the commuting MTOs read


  T^q_x &= _m,n a^_m+q,n a_m,n &
  T_y &= _m,n a^_m,n+1 a_m,n

justifying the plane wave ansatz(Here  is defined in the
  full BZ.)   subject to
the periodic boundary conditions . Substituting
this ansatz into Eq. eq:HH-schrodinger yields the
Harper equation 


Calculating the single-particle spectrum of
Eq. eq:HH-hamiltonian is therefore equivalent to solving the
eigenvalue equation

,  for the  matrix  defined
as(If ,  and
  .)



For integer  we recover the lowest Bloch band
Eq. eq:tb-dispersion, while for  this band
splits into  bands with dispersion , 
. 

Spectrum discussion
In the top-left panel of Fig.  we show an
example of the spectrum computed by diagonalizing  for
. The  level is at the centre of the middle band,
and there are  gaps separating the bands. For even values of ,
the spectrum is symmetric with respect to  and the two central
bands touch at zero energy. Around these degeneracy points, the
dispersion is linear , while in all other
cases it is quadratic near its minimum. The  case is
particularly interesting by analogy to graphene-related
physics , as there are two inequivalent
"Dirac points" and the Hamiltonian preserves time-reversal
symmetry. Inspecting the middle and right panels we note that, despite
the reduced symmetry of , the dispersion has the full
 rotational symmetry of the underlying square lattice. One can
formally define rotation and reflection operators which, together with
the MTOs, constitute the full projective symmetry
  group .  The bottom panel of
Fig.  shows the evolution of the lowest
energy gap , as well as the bandwidth  of the lowest
band as function of  (for the case ). The lowest band
gets flatter, as the continuum () LLL
eigenstates are the solutions of the lowest band of the HH
model , and  decreases like
, as the number of gaps is  (excluding the band
bottom and band top).


HH topology
The energy bands of the HH model are topologically nontrivial, as
characterized by their non-zero Chern numbers - topological
invariants that are linked to the quantization of the Hall
conductivity in the integer quantum Hall (QH)
effect . In fact, solutions of the Harper equation
Eq. eq:Harper were also studied in connection with the QH
effect . We make use of these topological
properties in Chapter , where we consider the HH model
in the presence of an external harmonic trap in a dissipative system.



Conservation laws for the GP field

In this Appendix we use (classical) field theory to derive the
conservation laws associated with the time-dependent Gross-Pitaevskii
Eq. eq:TDGP. We start by recasting the TDGP using Einstein's
summation convention, with the notation 
 and  (
)



This equation can be deduced by extremalization of the action
, where the Lagrangian density
 is given by



For compactness, let ,( is the speed of
  light, but it is just a constant for our purposes, and will drop out
  from all final results.)  and we can define the canonical
energy-momentum tensor , following the usual
field-theory prescription 


where the indices  run from 1 to 4. As the complex scalar
field  (and ) satisfy the Euler-Lagrange
equations  resulting in Eq. eq:GP-field
(and its complex conjugate), we have



from which energy and momentum conservation immediately
follow.(Conservation of angular momentum imposes the
  additional requirement that .)

To see this explicitly, we split the temporal and spatial parts of
Eq. eq:conservation-laws to obtain


  _t T_4 4 + i c _k T_k 4 & = _t L

  _t T_4 l + i c _k T_k l & = i c _l L
Direct subsitution into Eq. eq:en-mom-tensor shows that
, where the energy density  is given by



while , with the energy flux density



Regarding Eq. mom-cons, we get , with the mass
current density (compare to Eq. eq:bec-current)

and the momentum flux density tensor 



Combining all the above, one can rewrite
Eq. eq:conservation-laws as


  _t H + _k S_k & = N_0^2 _t U

  _t J_l  + _k _k l & = -N_0^2 _l U
We now focus on Eq. eq:mom-cons-pi, which shows momentum
conservation. For a more transparent physical interpretation, one can
make use of the Madelung transformation Eq. eq:madelung and
integrate over a volume , obtaining



where we made use of the divergence theorem for tensor
fields  in order to express the volume integral as an
integral over the surface  which encloses . Here 
is the outward-pointing unit normal to  at each point, and we also
introduced the total momentum of the field,
 (with  the velocity
Eq. eq:supervelocity). The hydrodynamic form of the momentum
flux density then reads



with the pressure
. We
can now write Eq. eq:gauss in vector form



and see that the vector between square brackets is the amount of
momentum per unit time and area "flowing" through a surface
orthogonal to .

If we now consider a stationary regime (so that the last term of
Eq. eq:vectorial-flow drops out), and assume the potential 
describes some fixed obstacle (such as a cylinder) inside , then by
definition  the "drag" force acting on a surface
element of the obstacle is the momentum flux though this element. The
total force, of course, will be given by the surface integral of the
momentum flux, or equivalently, by the volume integral of the
potential gradient in Eq. eq:vectorial-flow. While in
conservative systems this drag is caused by pressure gradients across
the obstacle , nonconservative ones such as
the resonantly pumped polariton fluid of Chapter  also
have a viscous contribution to the momentum flux tensor, giving rise
to Navier-Stokes-type physics.







Microcavity exciton-polaritons

In Chapter  we have seen, among other things, that
interactions are essential for the appearance of superfluidity, and
the existence of a BEC is not sufficient. We now want to explore
similar physics in the context of nonlinear optical systems, as their
nonlinear polarization is able to mediate appreciable interactions
between photons . One way to achieve the strong
interactions necessary for collective fluid-like behavior of the
many-photon system is by coupling the photons with matter degrees of
freedom, resulting in new quasiparticles called
polaritons .

The connection between BEC and spontaneous coherence effects in
nonlinear optical systems became apparent starting with the
experiments on exciton-polaritons in semiconductor
microcavities 
. However, the steady state in such devices typically
results from a dynamical balance between driving and dissipation, in
contrast to the thermal equilibrium achieved in ultracold atomic
gases .

In this Chapter we review the physics of microcavity exciton
polaritons, with an emphasis on their hydrodynamics properties, and in
particular the new effects that their non-equilibrium nature
introduces.

Chapter organization
This Chapter is organized as follows: we start from a microscopic
description of a semiconductor microcavity (Sec. )
in terms of electrons, holes and cavity photons
(Sec. ), before switching to a higher level
weakly-interacting Bose model for polaritons in
Sections  and . We
discuss various ways of pumping the system, as well as the important
roles played by decay (Sec. ) and disorder
(Sec. ), leading to a driven-dissipative
Gross-Pitaevskii equation at the mean-field level
(Sec. ). We investigate the optical limiter and
bistable regimes of resonantly pumped polariton condensates, by means
of a single-state (pump only) ansatz (Sec. ) for the
polaritonic GP equation. Finally, we consider a three-state solution,
the optical parametric oscillator regime (Sec. ) obtained
by parametric scattering of pump polaritons into the signal and idler
states. In Section , we also review some
superfluid-related phenomena in an OPO context, such as quantized
vortices, and the  symmetry breaking resulting in a Goldstone
mode of the OPO Bogoliubov spectrum.


Model system



Consider the following (simplified) physical system, at zero
temperature: we have an optical planar microcavity containing photons
confined along the growth direction (say ) by its (identical)
metalic mirrors, separated by a distance . At position
, the cavity contains a 2D quantum well (QW), which is
a thin semiconductor layer sandwitched in between another alloy acting
as a barrier. The setup is shown in Fig. (for an in-depth introduction to the physics of semiconductor
microcavities, see Ref. ). The photons create
electron-hole hydrogenlike bound pairs (excitons) confined by the
barrier inside the QW. Exactly describing the system of electrons,
holes and photons, taking into account their Coulomb interaction and
the disorder created by defects inside the QW, as well as
imperfections present in the cavity mirrors is a daunting task that
would provide little physical insight into this system. We therefore
use a series of approximations and construct an effective theory that
is more physical.

    
    
    Pictorial representation of microcavity polaritons, as electron-hole pairs (excitons) inside a semiconductor quantum well (QW), 
    coupled with the photons trapped by the dielectric Bragg mirrors of the cavity. 
    Momentum conservation dictates that a cavity field component of momentum  must decay into external radiation of frequency , emitted at an angle  satisfying .
    From Ref. .
    
  
Microscopic description


We start by writing down the kinetic energy of the cavity photons



where we've made use of the translational symmetry of the problem in
the  plane and labeled the bosonic creation (annihilation)
operators  () by the in-plane momentum
. These operators satisfy the standard commutation relations

  [a_, a_^^] & =_,^

  [a_, a_^] & =[a_^, a_^^] = 0

The confinement along the  direction leads to the quantization of
the photon momentum , with  indexing
the longitudinal modes. Fig.  depicts the
photonic field corresponding to the  mode as a standing wave,
representing the measurable electric field inside the cavity. This
field can be expressed in terms of photonic operators by
 H.c., where the photon mode wavefunction
 has the sinusoidal shape 




The cavity photon dispersion for each of these modes is therefore
, with  the speed of light
in the semiconductor. Expanding the square root for , one
gets . We see that cavity
photons, as opposed to free-space photons, have a quadratic dispersion
with an effective mass . For typical cavities used in
experiments,  is of the order of  (free) electron
masses  (see Tab. ).

We now turn our attention to the semiconductor QW. For simplicity,
consider a spin-polarized direct bandgap semiconductor such as GaAs,
so that we can neglect additional spin degrees of freedom. We further
assume that the dispersions of single particle states in the
conduction and valence bands have the quadratic form
 and
, with
 the semiconductor bandgap and  and  the
effective masses for electrons and (heavy) holes. Excitons have a
typical "extension"  (see
Tab.  for typical values for a GaAs-based
microcavity) called the (2D) exciton Bohr radius, with  the
static dielectric constant,  the electron charge and
 their reduced mass. Their binding
energy is called the exciton Rydberg, defined as
 (see
Tab. ).

We further define  and  as the operators which
create an electron in the empty conduction band, respectively a hole
in the filled valence band, and which obey the Fermi anticommutation
rules

  c_, c_^^ & =_,^

  c_, c_^ & = c_^, c_^^ = 0
and similar for . We can therefore write the electronic
Hamiltonian in terms of these new operators
as 


Here we have introduced the electron and hole densities,
 and
. The
matrix element of the Coulomb potential
 depends on cavity quantization area
, but this dependence drops out when passing from discrete sums
over states to integrals over the momenta .

Finally, the last ingredient of our microscopic description concerns
the interaction between electrons and photons. Making use of the
dipole approximation, we have



with the strength

written in the dipole gauge, where  is the inter-band
dipole matrix element.


Effective Hamiltonian


As mentioned in Sec. , we are interested in an
approximate description of the problem where we can consider the
excitons as fundamental quasiparticle excitations from the ground
state of the semiconductor and treat the Coulomb term in
Eq. eq:coulomb-hamiltonian as an effective exciton-exciton
interaction. We then couple the resulting excitons to light, and
obtain the so-called weakly interacting boson model for
polaritons. The name stems from the fact that, at moderate
electron-hole densities, such that the inter-exciton distance is much
larger than , one can assume excitons to behave esentially
as (composite) bosons , and therefore
describe them using the Bose creation and annihilation operators
 and . Making an Usui
transformation  and truncating the interaction
terms at fourth order results in the following three-part effective
Hamiltonian 


where we have defined

  H_0 & =_(a_^, b_^)_C(k)_R_R_X(k)a_b_

  H_XX & =12_,^,qU_- ^,qb_+ q^b_^ - q^b_^b_

  H_XC^sat & =-_R_satA_,^,q(b_^-q^b_+ q^b_a_^ + H.c.) contains the kinetic energy of the excitons and cavity photons,
as well as their harmonic coupling, ie. the conversion of an exciton
to a cavity photon at the Rabi frequency , which depends on
the overlap between the wavefunctions of the exciton and photon
through the  oscillator strength surface density  of
the excitonic transition



where  is the maximum amplitude of the cavity photon
electric field, at one of its antinodes from
Eq. eq:electric-field. It is worth noting that both
Eq. eq:dipolar-inter and Eq. eq:polariton-ham make use
of the rotating-wave approximation, neglecting antiresonant terms
which do not conserve the number of excitations, such as
 and . This is
justified so long as , which
is normally the case, as can be seen in Tab. ,
which reports typical parameters of a GaAs-based microcavity with
embedded QWs.

Eq. eq:polariton-ham can be diagonalized by means of a unitary
transformation of the form 


The normal modes of  are called upper and lower
  polaritons (UP and LP) and correspond to the Bose operators
 and , with eigenenergies

  
    
    
    Dispersion relation for the two polaritonic branches (UP, LP) as a
    function of emission angle  of photons outside the cavity,
    which is directly related to the momentum  of polaritons inside
    the cavity. The exciton and cavity photon dispersions are depicted
    with dashed (red) lines. Due to their light mass, the polaritons have
    a very sharp dispersion compared to excitons. Also note the typical
    anticrossing form characteristic of Hermitian coupling. Parameters:
     eV,  meV, 
    and  meV.
    
  


where the  signs refer to the UP and LP branches and the Hopfield
coefficients appearing in Eq. eq:hopfield-trans are

  X_k & =[1 + (_R_LP(k) - _C(k))^2]^-1/2

  C_k & =-[1 + (_LP(k) -
        _C(k)_R)^2]^-1/2The two branches of the polaritonic dispersion relation are shown in
Fig. , together with the dispersions of
the cavity photons and excitons. One can therefore think of polaritons
as being massive photons dressed by matter excitations.

Note that the 1s-exciton energy can be approximated as
, with the exciton mass
 of the order of the electron mass and
. Since, as we have already
noted,  is several orders of magnitude larger than the cavity
photon mass , one can in practice neglect the momentum
dependence of the excitonic dispersion, considering it as flat
compared to the cavity photon one. Furthermore, we can measure
energies starting from  and denote the detuning between
exciton and photon bands as 
. For zero detuning (), the splitting between
the two polariton branches is . In order to tune 
in experiments, one tipically builds the cavity mirrors with a
wedge .



 quantifies the effective repulsive exciton-exciton
interaction, with a momentum dependent strength
 that can be calculated from the
Coulomb exchange term in the Born approximation .  A
further approximation that one can make for wave vectors much smaller
than  is to reduce  to
a contact potential

, as shown in
Ref. . As a result,  physically represents the
interaction between two excitons in the same single-particle momentum
eigenstate .

Finally, the last term of Eq. eq:total-ham,
, can be interpreted as a saturation
effect due to the Pauli exclusion principle underlying the fermionic
character of the excitons . The exciton saturation
density  generally
depends on the specific shape of the internal wavefunction of the
exciton .  For the problems we will be treating in
the remainder of this manuscript it can be considered large enough to
justify neglecting this anharmonic contribution to the total
Hamiltonian Eq. eq:total-ham (see Tab. ).

Before concluding this Section, one should mention the validity
limitations of the weakly interacting boson model for polaritons
reviewed here. Its derivation implies an expansion of the original
fermionic operators of Eq. eq:coulomb-hamiltonian in powers of
bosonic operators, having as small parameter the number of excitons
per Bohr radius . Truncating this series at fourth order
results in an upper bound on the exciton density of around
. At higher densities, an electron-hole plasma-like
state is formed, with a a qualitatively different
behaviour .  We see that this model is
particularly useful for the case when the dominant interaction is the
Coulomb repulsion between the excitons.

Lower polaritons

Assuming the Rabi energy  to dominate over the kinetic and
interaction energies we can neglect the upper polaritons and project
Eq. eq:total-ham on the basis of lower-polariton states
only. Employing the approximations discussed above, we obtain the
lower-polariton Hamiltonian



where the strength of the repulsive interaction between lower
polaritons now depends on in-plane momentum through the Hopfield
coefficients(The interaction strength  can usually be
  rescaled to 1 by a simple transformation.)



This effective interaction, originating from the binary Coulomb
scattering of the excitons, is the one responsible for the collective
behaviour of the polaritons, and in particular superfluid
hydrodynamics, an aspect which we will explore in the following two Chapters.

Physically, the quantities  and  are the
excitonic and photonic fractions of the lower polariton mode and the
effective mass of the lower polaritons, , can be
computed from the curvature of their dispersion
Eq. eq:polariton-dispersion by



which at normal incidence  reduces to 



if we neglect the flat exciton dispersion in the limit 
. For zero detuning, we see that  and the
light and matter content of lower polaritons are each equal to .

As long as the interesting physics is limited to small (compared to
) wavevectors , close to the bottom
of the LP branch, one can safely approximate the LP dispersion as
being parabolic, by expanding Eq. eq:polariton-dispersion in a
Taylor series around  and truncating at second order. We
obtain the expression



with an effective mass given by Eq. eq:LP-mass. In practice,
many interesting experimental conditions can be adequately described
by limiting ourselves to the LP branch with a quadratic dispersion.


Pumping and dissipation

We have so far assumed that the mirrors at the two ends of the optical
cavity in Fig.  are perfect. In reality of
course, this is not the case, and photons continuously escape from the
cavity. In other words, decay is a built-in feature of this
system and, rather than being a hindrance, it is actually useful as it
allows one to extract information about the system. By looking at
the near-field (far-field) emission from the cavity, one can measure
the real (momentum) space photon density. The polaritonic spectrum can
also be obtained in this way, as a function of the angle  of
emission of a photon out of the cavity, as shown in
Fig. .

Since the system is lossy, one needs to drive it by means of an
incoming laser beam (called the pump), thus replenishing the photons
and achieving a stationary state. It is important to note, however,
that this state should not be confused with thermal equilibrium, as we
are dealing with intrinsic out-of-equilibrium physics
here. Experiments have so far explored two main ways of pumping these
cavities: resonantly (coherent) and non-resonantly (incoherent). We
will not concern ourselves with the latter in this manuscript, but we
point the interested reader to the excelent review in
Ref.  (and references therein) for details.

Resonant pumping is done by means of a continuous-wave coherent laser
beam which is transmitted through one of the mirrors, at energies and
momenta close to the LP or UP. The coherence of the beam is then
inherited by the polariton fluid.  As already mentioned, conservation
of in-plane momentum and photon frequency dictates that a laser wave
coming in at an angle  and frequency  resonantly
excites a microcavity mode with wavevector
 (see
Fig. ). As a result, we see that, besides
balancing the losses, the pump can also be used to tune the system to
a particular desired state, injecting polaritons of a given momentum
and energy. In the following we will consider that the pumping is done
quasi-resonantly, close to the bottom of the LP branch, by using a
homogeneous planewave pump of amplitude , wavevector ,
frequency 


Depending on the pump parameters (specifically the pumping angle,
which in turn determines ), one can distinguish two qualitatively
different regimes: 1) the pump-only state where the only stable
configuration of the system is the mode at  and 2) the regime
where polaritons at  undergo stimulated scattering and populate
two other modes. Indeed, if the pumping angle is such that the
wavevector it excites lies above the inflection point of the LP
dispersion, the system can enter a new regime called optical
  parametric oscillator (see Sec.  for details about
OPO).

To quantitatively account for the driving and decay, one can use
the input-output formalism that was extended to planar microcavities
in Ref. . As far as photonic losses are concerned, we
must first note that, apart from the radiative decay channels already
mentioned, there is also a nonradiative contribution due to photon
absorbtion inside the cavity. The total cavity photon decay rate
 can therefore be expressed as the sum
, where the radiative
decay rate is 


for normal incidence, where  is the speed of light in vacuum and 
the mirror transmittivity.

The excitons from the QW are of course also subject to decay
processes. While their radiative decay can only take place via one of
the cavity modes (due to strong light-matter coupling(In
  general, one is in the strong coupling regime when  is
  greater than the combined losses in the system.)), we must include
nonradiative recombination processes as well as the additional
dephasing  caused by interactions with carriers inside
the well and variations in well thickness. All this results in an
effective decay rate  for excitons, which is generally
weaker than the photonic one quantified by . Typical values
of both rates for GaAs-based microcavities are shown in
Tab. .


Disorder effects


It is noteworthy to remark that Hamiltonian Eq. eq:total-ham
does not include any disorder. We have already mentioned that there
are two distinct disorder classes in this problem. One is the
excitonic disorder, acting on length scales of around 10 nanometers,
which tipically affects the exciton oscillator strengths
, but not the spatial polaritonic density. This type of
disorder stems from variations in the QW thickness as well as alloy
imperfections and is typically not strong enough to dissociate
excitons. Neglecting excitonic disorder also implies that each
localized exciton state couples to precisely one extended photon
state. The resulting polaritons are then nothing but the coherent
superposition of these excitonic and photonic states.

Photonic disorder, on the other hand, acts on lengthscales
 on the order of microns (see
Tab. ) and therefore can disrupt the polariton
density profile. Photonic disorder is normally caused by the presence
of imperfections inherent to the growth process of the cavity
mirrors. These defects (as well as any layer mismatch) act as
scattering centers for the polaritons, and can be included in the
Hamiltonian as an additional external scalar potential. The simplest
example is that of a spatially localized defect of strength 
located at positon , which can be approximated by the potential




To get a more quantitative idea of the physics involved, it is worth
describing realistic samples typically used in state-of-the-art
experiments. Instead of simple metallic mirrors, such microcavities
are usually bound by a series of up to around 20 distributed Bragg
reflectors (DBRs), which are alternating dielectric layers of
thickness  and different refractive indices, achieving a
total reflectivity of over  . 

Multiple CdTe (or GaAs, GaN, ZnO) QWs a few nanometers thick are
placed inside the cavity of length m, at
the antinodes (maxima) of the photonic electric field, in order to
maximize the coupling between excitons and photons (see
Eq. eq:omega-R, which shows maximum coupling for
). When  distinct QWs are placed
in the cavity at the electric field maxima, the light-matter coupling
 is enhanced by a factor . Furthermore, the cavity
cutoff frequency  is normally chosen close to the exciton
frequency .

Typical values of the relevant parameters for GaAs-based
microcavities are shown in Tab. . 


    
  Typical parameters of a GaAs-based microcavity with
    embedded QWs. Left column shows properties of a QW, while the
    right one refers to the microcavity. For these values, the ratio
    between saturation and Coulomb interaction is small,
    . From
    Ref. .
  


Mean-field description

As it turns out, a useful approach for tackling the
hydrodynamics-related problems that we will encounter in the following
Chapters is the so-called mean-field approximation. In the context of
conservative, weakly interacting Bose gases at thermal equilibrium,
this approach leads to the celebrated Gross-Pitaevskii
equation , and in our case it will lead to its
driven-dissipative counterpart.

If the Rabi energy is larger than the pump detuning from the bottom of
the LP branch  and also larger than
the decay rates  and , we can start from the lower
polariton Hamiltonian of Sec. .  The idea
is that one replaces the lower polariton operators  in
Eq. eq:LP-ham with their expectation values
, where  is now
interpreted as a classical field (order parameter representing lower
polariton polarization) describing the state with momentum label
. The mean-field equation for  is then obtained from
the Heisenberg equation of motion
 by
replacing all the operators with their corresponding expectation
values. This mean-field treatment is of course exact as far as the
pumping and the losses are concerned, as well as for any two-operator
product such as the polariton kinetic energy. However, the effective
lower polariton interactions are approximated following an approach
similar to the one pioneered by Bogoliubov in the context of atomic
BECs, namely



We obtain a generalized Gross-Pitaevskii equation (GPE) in momentum
space under coherent driving, which includes losses and describes the
time evolution of 

  i_t_ = [_LP(k) - i2 _k]_ + __1,_2 g_,_1,_2 ^__1 + _2 -  __2 __1

  + C_ _q C_q V_d(- q) _q
  + C_ f_p (-i _pt) _,_p
Here  is the
(lower) polariton loss rate,  is the Fourier transform of
the external potential Eq. eq:defect-pot and the effective
polariton-polariton interaction strength is given by



Going beyond the mean-field approximation to include quantum and
thermal fluctuations of the (quantum) field  is not
an easy task in general, see Ref.  for some
approaches such as using the semiclassical Wigner representation.


In order to obtain the corresponding form of Eq. eq:mom-GPE in
real-space, we introduce the Fourier transform of the classical
wavefunction  by means of(Replacing the discrete
  sums over  states by continuous integrals is justified in the
  limit of a large system.)



Knowing that Fourier transforms take multiplication to convolution, we
immediately get the real-space equation for 
  i _t () = [ _LP(-i ) - i2 (-i ) ] ()
 + d_1,2,3 g(_1 - , - _2, - _3) ^(_1) (_2) (_3)

+ d_1,2 C(- _1) C(_1 - _2) V_d(_1) (_2) + C(_p) F(r,t)
Note that the kinetic energy operator has a particular dependence on
the gradient instead of the usual Laplacian, reflecting the
non-parabolic LP dispersion.  Furthermore, the nonlocal interaction
term in the second line of Eq. eq:convoluted-GPE has the form
of the most general third order nonlinear term that one can come up
with, if  were an arbitrary function. This freedom is
somewhat restrained in practice by the fact that  is related to
Eq. eq:LP-inter by



and we have also introduced the functions , which are
the Fourier transforms of the Hopfield coefficients . These
coefficients do not have a very pronounced momentum dependence, as in
fact for ,  is between  and . If we
neglect altogether the Hopfield coefficients in
Eq. eq:convoluted-GPE and consider the case of parabolic lower
polariton dispersion, we finally obtain the real-space
driven-dissipative GPE in its simplified form



Here we have also neglected the momentum dependence of the polariton
loss rate, justified by the fact that, at large momenta, the
diminishing radiative width is normally compensated by non-radiative
broadening.  To make a stronger connection with nonlinear optics, we
note that a similar equation holds for cavities containing a nonlinear
medium, with the role of  being taken by the 
susceptibility of the medium . This equation
allows for an ab-initio description of the coherently pumped system,
unlike the more complicated case of incoherent pumping where one
generally also needs a phenomenological description of the various
relaxation processes .

Before moving on, we note that one can also describe the dynamics of
polaritons scattering against a defect using the two-component
GPE  for the coupled exciton and
cavity fields :



where



Eq. eq:2-component-GPE is usually solved numerically, as we
shall see in Chapter .

Equation of state

Because we drive the system with the coherent plane-wave continuous
pump Eq. eq:pwpump, we can look for stationary state solutions
with similar oscillatory behaviour



In absence of an external potential, inserting this ansatz into
Eq. eq:mom-GPE, we obtain the equation of state for lower
polaritons



where we employ the simplified notation
, ,
 and
. One
immediately notices that, unlike the case of the atomic mean-field
Eq. eq:atom-MF, the phase of the polariton field amplitude
 is set by the pump, therefore explicitly breaking the 
gauge symmetry associated with global rotations of the condensate
phase. This implies that, unlike the atomic BECs, no Goldstone branch
will be present in the polaritonic excitation spectrum. Note also that
the oscillation frequency of the condensate wavefunction is fixed by
the pump frequency . In Chapter , we will
compare the main features of the quasiparticle excitation spectrum for
polaritons to equilibrium systems like cold atomic gases in the BEC
regime. Furthermore, a locked phase clearly prevents the formation of
phase dislocations, such as vortices and solitons. For this reason, it
has been suggested  and experimentally
realised  that the defect can be located just outside
the pump spot, leaving the condensate phase unbound and allowing
hydrodynamic nucleation of vortices, vortex-antivortex pairs, arrays
of vortices, and solitons to be observed when the fluid collides with
the extended defect. Similarly, nucleation of vortices in the wake of
the obstacle has been observed in pulsed experiments
 .



    
  
    
    Pump energy blueshift of Eq. eq:pump-mf-sq as a function
    of pump laser intensity in the optical limiter (left panel, 
 meV) and bistable (right panel, 
 meV) regimes. The arrows in the right panel show
    the hysteretic jumps between the lower and upper branches at different
    intensities, while the dotted middle branch is single-mode
    unstable. The black dot marks a saddle node bifurcation.
    
  

The repulsive polaritonic interaction is responsible for a blueshift
of the LP dispersion .  Introducing the pump energy blueshift
, we can rewrite
Eq. eq:pump-mf in absolute value as



where we have also defined the pump laser intensity
. One can now distinguish two
qualitatively different regimes, based on the value of the bare pump
detuning . The first regime, called
optical limiter is for
 and it is
plotted in the left panel of Fig. . We see that the
pump blueshift increases monotonically with increasing pump power. The
increase is sublinear, as the blueshift moves the lower polariton
dispersion out of resonance with the pump.

The second regime
, shown in the
right panel of Fig. , is called optical
  bistability and is no longer monotonic. The blueshift-intensity
curve has a characteric S-shape, and shows hysteresis: increasing the
pump power, we get to the first turning point (filled circle in the
figure) where the blueshift suddenly jumps from the low to the high
branch, while the reverse process of decreasing the pump power shows a
jump at a second turning point on the S curve. Defining an
interaction-renormalized pump detuning
,
we see that this quantity is positive up to the second turning point
on the S curve, and then changes sign, becoming negative on the higher
branch.

The middle branch of the S curve, with negative slope, is dynamically
unstable, while the two turning points where the stable and unstable
solutions meet are called saddle node bifurcations.  To see
this, one must consider small perturbations of wave-vector  of
the form



where . Substituting this into
Eq. eq:mom-GPE and linearising with respect to the amplitudes
 and , we obtain an eigenvalue problem from which the two
(complex) eigenvalues  can be determined in a straightforward
manner. 

Dynamical stability is then ensured when
 for any wave-vector .  We can
therefore distinguish two main types of instabilities, depending on
the value of : for  (the degenerate case),
the system shows a Kerr single-mode (or pump-only)
instability which only involves the pump wave-vector, while the
general case of  implies a parametric
instability, signaling the onset of the optical parametric oscillator
(OPO) regime. This regime involves two separate modes (apart from the
pump), namely the so-called signal and idler states
at wavevectors  and  respectively,
as we will discuss at large in Sec. .



OPO regime

    
    
    The driving (pumping) of a swing by a person, using the standing
    position, shown as an (albeit contrived) example of a degenerate
    parametric amplifier. From Ref. .
    
  
To understand why the OPO regime arises, one has to first look at the
form of the Heisenberg equation of motion Eq. eq:mom-GPE. This
equation describes the coupled dynamics of an arbitrary state ,
given an applied optical field which excites polaritons with in-plane
wavevector . As a result of polariton-polariton interactions,
the coherent occupation at generic wavevectors  and  can
scatter into  and , a process one can
represent as



As these states further scatter into others by the same matter
wave-mixing mechanism, we get a cascade of so-called satellite states,
with diminishing population.  

Out of the multitude of possible scattering channels, we limit our
description to the one describing the process
, which is
precisely the parametric scattering of two pump polaritons into the
signal and idler states that we just introduced in
Sec. . As it turns out, due to the peculiar shape
of the lower polariton dispersion relation
Eq. eq:polariton-dispersion, this scattering process can occur
in a completely resonant way, with conservation of both in-plane
momentum as well as energy, expressed by the phase-matching conditions
 and  (we
have introduced here the signal and idler frequencies). The parametric
scattering condition that has to be satisfied then reads



We must further emphasize the importance of the unique form of the LP
dispersion in Eq. eq:conservations, as it allows efficient
parametric coupling on a relatively broad range of pump momenta.
Parametric scattering would not be possible for particles with a
quadratic dispersion: in the atomic case for instance, one must employ
the use of an optical lattice in order to deform the dispersion and
promote parametric amplification . 

While this may all seem far-removed from everyday experience, we
remind the reader that parametric amplifiers and oscillators are
commonly found in mechanical (as well as electronic) systems. For
instance, one can consider a person on a swing in the standing
position, as shown in Fig.  (and described in
detail in Ref. ). In this particular example, we
are dealing with a degenerate parametric amplifier, as the signal and
idler have the same frequency and are represented by the swing. The
person driving the swing then acts as the pump, and one can see that
his (or hers) up-down motion has a periodicity which is twice that of
the swing, therefore  in this case.

A schematic depiction of the polariton parametric scattering process
with the signal state located at zero momentum and the pump close to
the inflection point of the LP dispersion is shown in
Fig. . It is clear that for a given
signal wavevector, Eq. eq:conservations then uniquely
determines the pump and idler states.

  
    
    
    Upper (UP) and lower (LP) polariton dispersions are shown as solid
    (black) lines, while the original exciton and cavity photon are
    represented by (red) dashed lines, as a function of the laser
    wavevector . The red circle represents the pump polaritons being
    parametrically scattered into the signal (upper blue triangle) and
    idler (lower green triangle) states. The scattering process, depicted
    by the curved arrows, is resonant, conserving both energy and
    momentum. Parameters: , , 
,  meV and 
.
    
  
The fact that  is close to  is not incidental, but was
observed in multiple experiments performed in the OPO regime, as well as numerical simulations of the
full GP equation . An intuitive explanation was
put forward , based on the
interaction-induced blueshift of the LP dispersion.

We now turn our attention to the mechanism explaining the onset of the
parametric scattering instability in the simple case of the optical
limiter regime shown in the left panel of Fig. . As
already mentioned, the polariton-polariton repulsive interaction
causes a blueshift of the LP dispersion proportional to the pump mode
population . This shift also affects the signal and
idler frequencies , bringing them in
resonance with the pump frequency . The pump population (and
hence the shift) of the optical limiter is a monotonically increasing
function of laser intensity, quantified by
Eq. eq:pump-mf-sq. As a result, there is a certain minimum
threshold intensity where the resonance starts and the
pump-only solution becomes unstable towards parametric
scattering. Further increasing the laser power, we hit an upper
threshold, where the blueshift becomes too large and the resonance is
lost, forbiding parametric oscillation. This mechanism defines a
certain window of pump intensities where the parametric instability
can arise, the region depicted by the dashed line in
Fig. . Inside this window, the pump-only solution still
exists, but it is unstable towards OPO.

  
    
    
    Signal  ([blue] upper triangles), pump  ([red] circles), and
    idler  ([green] lower triangles) mean-field energy blueshifts as a
    function of pump laser intensity in the optical limiter regime. The
    dashed line is the pump blueshift in the pump-only solution
    Eq. eq:pump-mf-sq, while the vertical dotted line marks the
    pump intensity used in Fig. . Parameters are
     meV.
    
  

The simplest ansatz one can make for describing the OPO regime above
threshold is to consider the three separate states (signal, pump and
idler), each with its own population, wavevector and frequency, in the
form 


where energy and momentum conservation dictate that
 and . As
before, we insert this ansatz into Eq. eq:mom-GPE and we
obtain three coupled complex nonlinear equations for 
which determine the dynamics of the signal, pump and idler. It is
important to note that these equations become closed only if we
neglect the multiple scattering processes mentioned at the start of
this Section and consider only the conversion of two pump polaritons
into a signal and an idler polariton. We will not give the full form
of these OPO equations of state, but instead refer the reader to
Ref. , where a complete discussion of the
oscillation threshold and ensuing dynamics is
given. Fig.  shows a particular solution of the equations
for the optical limiter case. One can see the smooth growth of the
signal population (and subsequent depletion of the pump) inside the
OPO (dashed) region. It is of course never implied that the whole OPO
window is dynamically stable, and in fact a stability analysis needs
to be performed by adding small fluctuations on top of the (now) three
coupled states, in a manner similar to
Eq. eq:pump-fluctuations, following
Ref. . The complete treatment will be presented in
Chapter .

Before moving on, it is worth noting that a limitation of the ansatz
Eq. eq:3-state-ansatz is that it does not allow solving the
so-called "selection" problem of determining the value of 
above threshold. As a result, this value has to be considered as an
input parameter to our problem. Futhermore, similarly to the pump-only
case and inherent in any mean-field approach, we are of course
neglecting any quantum fluctuations of the three macroscopically
occupied modes. Last but not least, we are also neglecting the
so-called TE-TM (transverse electric - transverse magnetic)
splitting , which can act as an effective
spin-orbit interaction term . This approximation
allows one to work with a single spin state, considering the circular
polarization of the pump to be fully transmitted to the signal and
ider as well.

We now turn our attention to a separate discussion, concerning the
phases of the three OPO states.  As before, the phase of the pump is
fixed by the outside laser, and now the sum of the signal and idler
phases is given by the matching condition
. Their phase difference however is a
free parameter, and as such it is spontaneously chosen by the system
every time the threshold is crossed. In summary, the OPO equations of
motion for the three states posess a  symmetry, corresponding to
a simultaneous and opposite phase rotation of the signal and idler,
that is


  
    _s& _s e^i  

    _i & _i e^-i 
  
This intrinsic symmetry is spontaneously broken above threshold, a
random value of the phase of the signal (and hence idler) being
selected at each realization of the experiment. 

We can thus see the onset of OPO as a phase transition in
out-of-equilibrium conditions, signaling the spontaneous macroscopic
occupation of the signal and idler states. Below threshold, the
incoherent signal and idler emission starts being stimulated and then
becomes coherent once the threshold is crossed. These coherence
properties were studied both experimentally  as well as
theoretically, by means of quantum Monte Carlo
methods . The order parameter of the transition is
the matter polarization, the analog of magnetization in the
ferromagnetic phase transition. At the critical point, noise
fluctuations are amplified and a macroscopic polarization is
spontaneously created. 

The phase transition can be first or second order, depending on the
pump parameters - both behaviors have been observed in
experiments . In the optical
limiter case detailed above, both pump and signal are continuous
functions of the incident laser power, and the transition is second
order as the signal increases smoothly. In the pattern formation
language, one clasifies the associated bifurcation as being of the
supercritical Hopf type. If one considers the bistable pump
regime of the right panel of Fig.  however, things
get more complicated, as the bifurcation can now be of the
subcritical Hopf type and hence the associated OPO initiates
discontinuously . 

It is relatively well known that in an infinite uniform
two-dimensional equilibrium system, thermal fluctuations at any
strictly nonzero temperature are strong enough to destroy the fully
ordered BEC state. The associated phase transition in that case is
instead described by the so-called Berezinskii-Kosterlitz-Thouless
(BKT) theory, and is driven by interactions as opposed to being a
purely statistical phenomenon. As the polariton system is out of
equilibrium (and in practice no experimental system is truly infinite
anyway), the exact nature of the OPO transition in this case has not
yet been fully elucidated, although it seems that BKT-like physics has
been both predicted  and
observed  under incoherent pumping.

The breaking of the  symmetry in an atomic gas below the BEC
critical point is related to the appearance of a soft Goldstone mode
in the form of zero sound. Similar physics arises below the Curie
temperature in a ferromagnet, where the magnon excitations (spin
waves) are the corresponding Goldstone bosons due to the spontaneous
breaking of the initial rotational symmetry. Since the resonantly
pumped polariton system above threshold spontaneously breaks a 
symmetry, Goldtone's theorem predicts a massless mode corresponding in
this case to a spatial twist of the signal-idler phases. Indeed, such
a mode has been identified in Ref.  by calculating
the quasiparticle excitation spectrum on top of the three-state OPO
anzatz Eq. eq:3-state-ansatz. The Bogoliubov spectrum in
question is shown in Fig. , and the Goldstone mode
represented by the solid heavy line labeled "G".

  
    
    
    Real (left) and imaginary (right) parts of the Bogoliubov spectrum
    of collective excitations on top of the OPO stationary state as a
    function of . The chosen pump intensity is
    marked with a dotted line in Fig. . The breaking of the
     symmetry corresponding to the simultaneous phase rotation of
    the signal and idler is manifest in the presence of the Goldstone mode
    G.
    
  

 
                
    
    Top panels: Integrated photon emission intensities for the
    signal  (red dots), pump  (blue squares) and idler  (green
    triangles) as a function of the renormalized pump intensity (left) and
    real-space signal emission corresponding to positions 1-4, including
    random photonic disorder, shown as contour lines (right).
    Middle panels: Photonic component of the OPO spectrum (log
    scale) as a function of energy and momentum (left) and integrated in
    momentum (right). Note that the emission is -like in energy.
    Bottom panels: Full photonic emission (left) and filtered
    real-space emission of signal, pump and idler (last 3 panels). From
    Ref. .
    
  
Note that both the real and imaginary parts of the frequency of this
gapless mode tend to zero in the long wavelength limit
. Close to this point, the real part has a
non-zero slope due to the finite wavevector of the injected pump
polaritons, while the imaginary part has a parabolic shape. This means
that, as opposed to the BEC case, a localized perturbation in the OPO
regime does not freely propagate as a sound wave. Instead, one should
observe a diffusive-like behaviour, where the localized signal-idler
phase perturbation decays and is then dragged along by the flow of the
pump. Another important difference compared to the BEC spectrum is the
absence of a singularity of the Goldstone branch around the
 point.

Further insight into OPO physics for realistic experimental conditions
can be gained by numerically solving the full 2-component GP
equation eq:2-component-GPE for the coupled exciton and cavity
photon fields, as first done in Ref. .  A
similar approach, reviewed in detail in Ref. ,
will be employed in Chapter , where we study
multicomponent superfluidity in the OPO regime.
Using a continuous wave pump with a top-hat profile, one injects
polaritons at a specific wavevector, close to the inflection point of
the LP dispersion. Above a certain threshold pump power, random
numerical noise starts getting amplified and the OPO instability sets
in, populating the three states. After the system reaches a steady
state, one can obtain the time evolution of the photon and exciton
fields  and  both in real and momentum space. Since
experiments measure the photonic component of the polariton field,
that is what is plotted in Fig. . Filtering in a
narrow cone in momentum (or equivalently, in energy), one can obtain
separately the emission of the signal, pump and idler states, while
the energy-momentum spectrum can be calculated as the discrete Fourier
transform of the total photonic wavefunction recorded at equally
spaced time-intervals.

Inspecting Fig. , we see that the threshold is a
continuous one in this case, with the OPO first switching on in a
narrow region in the center of the system and then extending to the
size of the whole pump spot before finally surviving only on its edge
when pumping is increased even further. Interestingly, the OPO
transition is signaled by the appearance of a moving striped pattern
in the polariton density (see bottom left panel), spontaneously
breaking the translational symmetry of the system. This pattern is
caused by the interference of the coherent emission of the three
states, and is similar to the Rayleigh-Benard instability in the
field of fluid dynamics. The threshold point in the OPO case is where
the stimulated scattering rate exceeds the losses, while the
Benard cells form when the temperature gradient (playing the role
of pumping) exceeds the viscosity (dissipation). In these simulations,
the pump (and hence also signal and idler) wavevectors are oriented
along the  direction, , so the fringes
appear along the  axis. They are however usually not visible
in experiments, due to the time integration of the photonic emission
in most experimental setups. Seen in this new light, the OPO phase
transition has an associated Goldstone mode that would correspond to a
rigid translation of the spatial stripe pattern.

The photonic component of the OPO spectrum is shown in the middle
pannel of Fig. . Integrating the spectrum over all
momenta, one sees that, as a result of being in the steady state, the
emission is -like in energy, while it is broadened in momentum
due to the finite size of the pump. The satellite states, equally
spaced in energy and less populated, are also visible in the middle
left panel. They are due to secondary scattering processes, from the
signal/idler into the pump and secondary-signal/secondary-idler, and
so on.

Finally, the filtered emission from the signal, pump and idler states
is shown in the lower panels of Fig. . Since we are
in the steady state, the density profiles are time-independent and
the idler emission is less intense due to the smaller radiative decay
of polaritons at large momenta. One important difference from cold
atomic gases in the BEC regime is the presence of nonvanishing
currents, which can flow across the microcavity even in the ground
state in a polariton system. These currents can be seen in the second
panel for the signal state, and are due the complex interaction
between nonlinearities, dissipation and parametric gain. Note that we
are not refering here to the uniform flow caused by the signal being
at finite momentum , as in fact this dominant current is
substracted from the plot.

OPO superfluidity

    
    
    Snapshots of the time evolution of the phase (upper row) and the
    density (bottom row) of an OPO polariton condensate triggered by a
    pulsed probe carrying a  vortex. From Ref. .
    
  

While still on the topic of persistent currents, we can make the
connection to Sec. , and discuss
superfluid-related phenomena in the OPO regime. A prime example of
this is the stability of quantized vortices, reported by Sanvitto and
coworkers in Ref. . They used a pulsed, narrow
probe beam in a Gauss-Laguerre state in order to imprint a finite
angular momentum on the OPO condensate previously created by means of
a wide pump spot. In a similar fashion to the rotating BEC experiments
performed with ultracold gases, the system responded by forming a
quantized vortex, which eventually drifted out of the condensate, due
to the simply-connected geometry (see
Fig. ). The vortex was visible in both
the signal and idler state, and its lifetime was much longer than the
average polariton lifetime.
In a similar direction,
Ref.  presented a theoretical study predicting
the formation of spontaneous quantized vortices in the OPO regime,
provided the pump spot is spatially narrow. These vortices were shown
to drift to their stable equilibrium positions, caried by the
nontrivial OPO currents discussed in the previous paragraph. Further
work by Tosi et al.  investigated the appearance and
dynamics of vortex-antivortex pairs created by a pulsed probe, much
narrower than the pump spot. It was shown that, while both vortices
and antivortices form along the edge of the probe laser, they form at
different relative positions, as determined by the steady-state
currents. We will deepen our investigation of OPO superfluidity in
Chapter , where we look at the flow of the signal, pump
and idler states against a small defect.

Finally, departing somewhat from the OPO regime, it was
shown in Ref.  that, despite not strictly
obeying the Landau criterion, the superfluid density of incoherently
pumped polaritons need not vanish. However, due to their
nonequilibrium nature, the normal fraction would remain finite even at
zero temperature. Last but not least, the same publication proposes
using artificial gauge fields (generated by applying a real magnetic
field in an anti-Helmholtz configuration) for measuring the superfluid
and normal components of the system.








Applications












Drag in a coherently-driven polariton fluid

In a conservative quantum liquid flowing past a small defect, the
Landau criterion for superfluidity (presented in
Sec. ) links the onset of dissipation at a
critical fluid velocity with the shape of the fluid collective
excitation spectrum . In particular, for weakly
interacting Bose gases, the dispersion of the low-energy excitation
modes being linear implies that the critical velocity for superflow
coincides with the speed of sound . Clearly, this is strictly
correct only for vanishingly small
perturbations , while for a defect with
finite size and strength, the critical velocity can be smaller than
 , due to vortex creation by
the macroscopic defect.

    
    
    Experimental images of the real- (top row) and momentum- (bottom
    row) space polariton density extracted from the near/far-field light
    emitted from the cavity. The different columns correspond to
    increasing values of the polariton density, from left to right. For
    the highest density, polariton superfluidity is apparent as a
    suppression of the real-space modulation (panel III) accompanied by
    the collapse of the Rayleigh scattering ring (panel VI).
    
    From Ref. .
    
  

However, even for perturbatively weak defects, in out-of-equilibrium
systems, where the spectrum of excitations is complex, the validity of
the Landau criterion has to be
questioned . In the
particular case of coherently driven polaritons in the pump-only
configuration, it has been predicted ,
and later observed , that scattering is suppressed at
either strong enough pump powers or small enough flow velocities (see
Fig. ). Yet, on closer scrutiny, it has been shown
that, despite the apparent validity of the Landau criterion, the
system always experiences a residual drag force, even in the limit of
asymptotically large densities  or small
velocities. This result has been proven by numerically solving the
Gross-Pitaevskii equation describing the resonantly-driven polariton
system in presence of a non-perturbative extended defect. Here, the
drag force exerted by the defect on the fluid has been shown to
display a smooth crossover from the subsonic to the supersonic regime,
similar to what it has been found in the case of non-resonantly pumped
polaritons . In this Chapter, we find an even
richer phenomenology for the dependence of the drag force on the fluid
velocity and two different kinds of crossovers from the sub- to the
supercritical regime. Furthermore, we show that the origin of the
residual drag force, which, in agreement with
Ref. , lies in the polariton lifetime only, can
be demonstrated even within a linear response approximation.

More specifically, we apply the linear response theory described in
Sec.  for the case of an equilibrium
condensate and extended here to a driven-dissipative fluid in order to
analytically evaluate the drag force exerted by the coherently driven
polariton fluid in the pump-only configuration on a point-like
defect. To simplify the formalism, we restrict our analysis to the
case of resonant pumping close to the bottom of the lower polariton
dispersion, where the dispersion is quadratic. Here, the properties of
the collective excitation spectrum have been shown to be uniquely
determined by three parameters only : the fluid
velocity , the interaction-renormalised pump detuning ,
and the polariton lifetime . In particular, the sign of the
detuning  determines three qualitatively different types of
spectra: linear for , diffusive-like for ,
and gapped for .

For both cases of linear and diffusive spectra, we find a
qualitatively similar behaviour of the drag force as a function of the
fluid velocity : In particular, the drag displays a crossover
from a subsonic or superfluid regime - characterised by the absence
of quasiparticle excitations - to a supersonic regime - where
Cherenkov-like waves, similar to those of
Sec. , are generated by the defect and
propagate into the fluid. The crossover becomes sharper for increasing
polariton lifetimes  and displays the typical threshold
behaviour for  with a critical velocity given by the
speed of sound of the linear regime, , exactly as for weakly
interacting equilibrium superfluids (in the case of perturbatively
weak defects). This behaviour is similar to the one predicted for
polariton superfluids excited non-resonantly ,
where the spectrum in that case is diffusive-like.

However, for gapped spectra at , we find that the
critical velocity governing the drag crossover exceeds the speed of
sound, , and we determine an analytical expression of 
as a function of the detuning . Furthermore, for ,
the drag has a threshold-like behaviour qualitatively different from
the one of weakly interacting equilibrium superfluids, with the drag
jumping discontinuously from zero to a finite value at .

We evaluate the drag as a function of the polariton lifetime 
and find for all three cases that: In the supercritical regime,
, the lifetime tends to suppress the propagation of the
Cherenkov waves away from the defect and therefore to suppress the
drag. Instead, well in the subcritical regime, , we find
that the residual drag goes linearly to zero with the polariton
lifetime , in agreement to what  was found in
Ref. , by making use of a non-perturbative
numerical analysis for a finite size defect. Similar to
Ref. , here, we do also find that the residual
drag in the subcritical regime can be explained in terms of an
asymmetric perturbation induced in the fluid by the defect in the
direction of the fluid velocity.

This Chapter is structured as follows: In Sec.  we
briefly review the linear approximation scheme and extend it to the
case of a resonantly pumped polariton fluid. We classify the three
types of collective excitation spectra in the simplified case of
excitation close to the bottom of the LP dispersion in
Sec. . In Sec.  we derive the drag force
and characterise the crossover from the subsonic to the supersonic
regime in the three cases of zero, positive and negative detuning. In
this section, we also evaluate the drag as a function of the polariton
lifetime, interpreting therefore the results of
Ref. .  Brief conclusions are drawn is
Sec. .


Linear response

As explained in Chapter , the description of
cavity polaritons resonantly excited by an external laser is usually
formulated in terms of a classical non-linear Schrodinger equation
(or Gross-Pitaevskii equation) for the LP field 
(see Eq. eq:rspace-GPE):



The LP dispersion is expressed in terms of the photon
 and exciton
 energies, the photon mass , and the Rabi splitting
 (see Eq. eq:polariton-dispersion):



Because polaritons continuously decay at a rate  (see
Sec. ), the cavity is replenished by a continuous
wave resonant pump  at a wavevector  (we will
later assume  directed along the -direction,
) and frequency  (see
Eq. eq:pwpump):


Note that, as discussed in
Secs. -,
Eq. eq:basic is a simplified description of the polariton
system. It implies that the interaction nonlinearities are small
enough not to mix the lower and upper polariton branches and that we
pump far from the UP dispersion. Moreover, starting from a formulation
in terms of coupled exciton and photon fields, the polariton lifetime
would be momentum dependent and, similarly, the polariton-polariton
interaction strength  is not contact-like as instead assumed in
Eq. eq:basic. However, as shown in Appendix ,
these simplifications do not affect our results qualitatively, rather,
they allow us to write them in terms of simpler
expressions. Furthermore, we have checked that, whenever the system is
excited near the bottom of the lower polariton dispersion, the results
for the drag force reported in Sec.  coincide with the
ones obtained using the photon-exciton coupled field description
Eq. eq:2-component-GPE introduced in
Chapter .

The potential  in Eq. eq:dispe describes a
defect, which can be either naturally present in the cavity
mirror  or it can be created by an additional
laser . Later on, we will assume the defect to be
point-like  and weak, so that we can
apply the linear response approximation .

As detailed in Sec. , one divides the
response of the LP field in a mean-field component 
corresponding to the case when the perturbing potential is absent, and
a fluctuation part  reflecting the linear
response of the system to the perturbing potential (see
Eq. eq:ansatz-atoms):



By substituting eq:mfield into eq:basic, we obtain a
mean-field equation and by retaining only the linear terms in the
fluctuation field and the defect potential, the following first order
equation in  (the equivalent of
Eq. eq:GP-atoms-system):



where the operator  is given by (compare to
Eq. eq:ourL for the atomic case):



with 
. We solved the complex cubic mean-field
equation eq:pump-mf for  in
Sec. . Here, we want to study the response of the
system to the presence of the defect and how different behaviours of
the onset of dissipation can be described in terms of the different
excitation spectra one can get for polaritons resonantly pumped close
to the bottom of the LP dispersion.


 


Collective excitation spectra for the subsonic (thick
solid [black] line at , with )
and supersonic (dashed [red] line at ) regimes and for an
interaction-renormalised pump detuning  (a,
b),  (c, d),  (e, f) and
 (g, h). Real parts of the spectra are
plotted in the left panels and the corresponding imaginary parts in
the right panels for  - note that in our
description the spectrum imaginary parts do not depend on the fluid
velocity .




Spectrum of collective excitations

    
  
    
    Momentum-space response (top row) 
 (arb. units) and normalised real-space
    wavefunction (bottom row) 
    corresponding to the nondiffusive spectra (a)-(d) in
    Fig. . System parameters: ,
    , and  is 
 (left column),  (middle
    column) and 0 (right column).
    
  
    
  
    
    Momentum-space response (top row) 
 (arb. units) and normalised real-space
    wavefunction (bottom row) 
    corresponding to the diffusive spectra (e)-(h) in
    Fig. . Polariton lifetime 
, and , 
    (left column); ,  (middle
    column) and ,  (right
    column).
    
  
The spectrum of the collective excitations can be obtained by
diagonalising the operator  in the momentum space
representation, following the approach of
Sec.  (see Eq. eq:translated-L)



where, . The description of the
spectrum simplifies in the case when the pumping is close to the
bottom of the LP dispersion, that can be approximated as parabolic
(see Eq. eq:parabolic-disp)



where  is the fluid velocity, and
 is the LP mass of Eq. eq:LP-mass. This
simplification allows one to describe the complex spectrum in terms of
three parameters only, namely the fluid velocity , the
interaction-renormalised pump detuning (defined in
Sec. )



and the LP lifetime . One observation we can immediately make
is that, unlike the atomic case, one will no longer have a sharp
corner(Unless, of course, .)  in the spectrum
at , as that feature depended on the equality (up
to a sign) between the diagonal and off-diagonal elements of
.

The elementary excitation spectrum of coherently-driven polaritons
reads (compare this to the spectrum given by
Eq. eq:boosted-bogoliubov):



where . If energies
are measured in units of the mean-field energy blue-shift 
 (we will use the notation 
 and ), then the
fluid velocity  is measured in units of the speed of sound 
. In order to make connection with the current
experiments, note that, for blue-shifts in the range 
 meV, typical values of the speed of sound  are
 m/s. Similarly, for common values of the LP mass,
the range in momenta in Fig.  comes of the order of
 m.


    
    
    Hermitian (left column) and anti-Hermitian (right column) coupling
    between two damped oscillators of energies , , decay rates
    ,  and interaction energy . Top row shows the
    coupling matrices, while middle (bottom) row shows the evolution of
    the real (complex) eigenvalues, with increasing . From
    Ref. .
    
  

The spectrum eq:spect can be classified according to the sign
of the interaction-renormalised pump detuning
  - see
Fig. . For  (panels [a,b]), the
real part of the spectrum lacks the sonic behaviour at small
, and shows a gap that increases with
, while the imaginary part is determined by the
polariton lifetime  only.

If one applies the Landau criterion for the real part
of the spectrum only, then one finds a critical velocity



always larger than the speed of sound for . 

If the fluid velocity is subcritical,  (see the [black] solid
lines in Fig. (a)), then no quasiparticles can
be excited and thus, for infinitely living polaritons ,
the fluid would experience no drag when scattering against the
defect. 

For the case of , the  branch
touches the  plane in one point (see left column of
Fig. ), resulting in a localized perturbation
around the defect, with an additional stripe pattern of wavevector
, caused by the interference of the pump with the momentum of
the scattered state.

For supercritical velocities instead,  see the [red] dashed
lines in Fig. (a), one expects dissipation in
the form of radiation of Cherenkov-like waves from the defect into the
fluid. In the supercritical regime, the set of wavevectors 
for which  form a closed curve in
-space with no singularity of the derivative (see middle
column of Fig. ). As explained in the
geometrical model presented in Sec. , this
means the radiation can be emitted in all possible directions around
the defect. This, as we will see in the next section, will imply that
the drag force for  goes abruptly, rather than
continuously, from zero at  to a finite value at
.

The spectrum gap closes to zero in the resonant situation at
, when the two branches  touch at
 (panels [c,d] of Fig. ).

Here, the real part of the spectrum displays the standard sonic
  dispersion at small wavevectors (as for the weakly interacting
bosonic gases of Chapter ) with a slope given by
. The imaginary part, as in the previous case, is
constant and equal to . It is clear therefore that in this
case, when , one recovers the equilibrium results valid
for weakly interacting gases ,
where the critical velocity for superfluidity equals the speed of
sound, , and the drag displays a threshold-like
behaviour. Here, in the supersonic regime , the closed
curve  has instead a singularity,
resulting in the standard Mach cone of aperture , 
, inside which radiation from the defect cannot be
emitted  (see right column of
Fig. , as well as Fig. of Sec. ).

For , the real parts of the two Bogoliubov branches
cross and stick together, giving rise to flat regions near the
crossing points. This is due to level attraction caused by the
anti-Hermitian coupling of Eq. eq:opell, and is accompanied
by a splitting of the imaginary parts, as illustrated in
Fig. . We further distinguish two separate
cases,  and .

In the first case (panels [e,f]), there is only one flat region in
momentum-space, and the intensity of the Rayleigh scattering is
amplified on a segment situated at  and oriented
parallel to the  axis, as shown in the left column of
Fig. . A consequence of this (parametric)
amplification is the long shadow seen in the real-space
image . As the imaginary part only has one peak
(corresponding to the pump), this is the precursor of the Kerr
single-mode instability described in Sec. .

The second case (panels [g,h]) has a different momentum-space
topology, with two flat regions, producing two distinct peaks in the
imaginary part. As soon as these become positive, one enters the
regime of parametric oscillation (see Sec. ), with the
signal and idler momenta situated at the two maxima. The corresponding
far-field image plotted in the middle column of
Fig.  shows two pronounced peaks on the line
passing through , the interference of which results in
a near-field "zebra-like" pattern . For very small
fluid velocities (see right column of Fig. ), the
two Bogoliubov branches cross on a ring of wavevectors, giving rise to
cylindrical wavefronts .

We note that spectra [e]-[h] have no correspondence in equilibrium
systems, because a finite polariton lifetime  is needed in
order to insure stability, . Following
the literature , we refer to these spectra as
diffusive-like.

We also note that, for these spectra, even if considering only the
real part of the collective excitation spectrum, as soon as the fluid
is in motion, dissipation in the form of waves is possible.

However, we will see that, similar to the case of non-resonantly
pumped polaritons , when decreasing  (and
accordingly  in order to have stable solutions), this
situation continuously connects to the case where a threshold-like
behaviour with  was found.


We will see in the next section how these different spectra imply only
two qualitatively different types of crossover of the drag force as a
function of the fluid velocity, for either  or
 pump detunings.


 


Drag force  as a function of the fluid velocity
 for different values of the pump detuning :
 (a),  (b), and 
(c), and for different values of the polariton lifetime - here, we
use the notation  and 
.



Drag force

The steady state response of the system to a static and weak defect
can be evaluated starting from Eq. eq:linre:

*

For a point-like defect, this can be written in momentum space as:

*

while the other component  can be
obtained by complex conjugation and by substituting
. The drag force exerted by the defect on the
fluid is given by the expectation value of the operator
 over the condensate
wavefunction :



This definition is justified for conservative systems in
Appendix , using the momentum-flux tensor. In
the steady state linear response regime, we obtain:


  F_d = g_V dk(2)^2 ik
  [_p^* _s (k + k_p) + _p     _s^* (k_p - k)]


  = 2g_V^2_p^2 dk(2)^2 ik
    (k)_+ (k)_- (k)
    
The drag is clearly oriented along the fluid velocity ,
i.e., . If , then the
integral in Eq. eq:dragf is finite only if poles exist when
, i.e., when quasiparticles can be
excited, in agreement with the Landau criterion. For finite polariton
lifetimes, however, it is clear that the integral will always be
different from zero for .


We now analyse the behaviour of the drag force as a function of the
fluid velocity for the three (, , and
) different spectra illustrated in the previous section.





Drag force  as a function of the inverse polariton
lifetime  in the (a) subcritical regime
() and (b) supercritical regime (). In both
cases we have fixed  ()
but these results are qualitatively similar for any other value of the
pump detuning. We plot in the right panels the normalised real-space
wavefunction  for two specific
values of  and .


For the linear spectrum, at , in the equilibrium
limit, , we recover for the drag the known result of
weakly interacting Bose gases in two
dimensions :



with a threshold-like behaviour at a critical fluid velocity equal to
the speed of sound . This limiting result is plotted as a bold
gray line in the panels (b,c) of Fig. . For
 and finite lifetimes , we find a smooth crossover
from the subsonic to the supersonic regime, with the drag being closer
to the equilibrium threshold behaviour for decreasing  (see
Fig. (b)). A finite lifetime tends to increase the
value of the drag in the subsonic region , giving place
to a residual drag force, similar to what was found in the numerical
simulations of Ref.  (see
Fig. ). Instead, in the supersonic region
, the finite lifetime tends to decrease the value of the
drag. 

In the case of diffusive-like spectra at  the
situation is qualitatively very similar to the resonant case (see
Fig. (c)), with the difference that now, in order to
have stable solutions, we can decrease the value of the lifetime only
by decreasing accordingly also the value of the pump detuning
. The crossover for both  and 
is also qualitatively very similar to the case of non-resonantly
pumped polaritons , where the spectrum of
excitation is diffusive-like. In Ref. , a similar
approach was used to analytically derive the drag for nonresonantly
pumped polaritons in 1D, finding a continuous crossover from a regime
dominated by viscous drag (of Stokes type) to one dominated by wave
resistance. Futhermore, the authors proved that it was not possible to
separate the viscous and wave-resistance components of the drag.
Shortly after the publication of Ref.  (on which
this Chapter is based), Van Regemortel et
al.  found that the parametric amplification
mechanism described in Sec.  can increase the
scattering of particles in the flow direction (for low condensate
speeds), leading to a negative drag force, i.e. a force
directed opposite to the flow direction. To see how this comes about,
one can look at the right column of Fig. : the top
panel shows that the average momentum of the scattered modes is in the
positive  direction. This, in turn, leads to a pileup of
the fluid density behind the defect (), as can be seen in the
bottom panel.

    
    
    Top row: Photon density along , perturbed by a
    circular potential of radius  m and height 110 meV (centered
    at the origin), in the subcritical regime. The pump momentum is 
 m, and the panels (from left to right) correspond to
    decreasing exciton () and photon () decay rates
     1.3, 0.44 and 0.011 meV. 
    
    
    Bottom row: Residual drag force as a function of decay
    rate. The left panel compares two different values of , namely
    0.7 (black solid line) and 1 m (red dashed line). The
    right panel shows two different potential heights, of -22 (blue
    dashed line) and -2.2 meV (red dashed line), with 
    m.
    
    In all plots, the pump is 0.44 meV blue-detuned above the bare LP
    branch . Adapted from
    Ref. .
    
  
In the case of gapped spectra, the situation is, however,
qualitatively different (see Fig. (a)). For infinitely
living polaritons, , the drag force can also be
evaluated analytically and its expression is similar to
Eq. eq:drag0, but with a critical velocity larger than the
speed of sound, the expression of which is given in
Eq. eq:criti:



Therefore now the drag experiences a jump for , rather than a
continuous threshold as for the resonant case . As already
mentioned in the previous section, this discontinuous behaviour of the
drag for the gapped spectra is connected to the fact that, as soon as
quasiparticles can be excited by the defect at ,
Cherenkov-like waves can be immediately emitted in all directions,
rather than being restricted in a region outside the Mach cone as
before. For , the cone was gradually closing with
increasing  fluid velocity.

Both the increase of the value of the drag in the subcritical region
as a function of the polariton lifetime and the decrease in the
supercritical region, are behaviours common to all the types of
spectra. We plot the drag force as a function of  in
Fig. , for two values of the fluid velocity  and a
specific value of the pump detuning , though we have checked
that the following results are generic. For , we find that
the residual drag is a finite-lifetime effect only and, well below the
critical velocity, the drag force goes linearly to zero for
. This is in agreement with the results of
Ref. , where the full GP equation for the
coupled exciton and photon fields was solved numerically for a
finite-size defect: see the bottom row of Fig.  for
the numerical drag in the limit of asymptotically large densities.

In the resonant case , the slope of the
drag for  can be evaluated analytically starting from the
expression eq:dragf:

*

The residual drag in the subsonic regime is an effect of the
broadening of the quasi-particles energies: Even when the spectrum
real part does not allow any scattering against the defect (e.g., for
), the broadening produces some scattering close to
the defect. This results in a perturbation of the fluid around the
defect, asymmetric in the direction of the fluid velocity (see panel
(a) of Fig. ), similar to what was obtained in
Ref.  (see top row of
Fig. ). Instead, in the supersonic regime, the drag
force is weaker in the non-equilibrium case with respect to the
equilibrium one. This is caused by the finite lifetime tending to
suppress the propagation of the Cherenkov waves away from the defect,
as shown in panel (b) of Fig. .


Conclusions and discussion

To conclude, we have analysed the linear response to a weak defect of
resonantly pumped polaritons in the pump-only state, and we have been
able to determine two different kinds of threshold-like behaviours for
the drag force as a function of the fluid velocity. In the case of
either zero or positive pump detuning, one can continuously connect to
the case of equilibrium weakly interacting gases of
Chapter , where the drag displays a continuous
threshold with a critical velocity equal to the speed of
sound. However, for negative pump detuning, where the spectrum of
excitations is gapped, the drag shows a discontinuity with a critical
velocity larger than the speed of sound. In this sense, the case of
coherently driven microcavity polaritons in the pump-only
configuration displays a richer phenomenology than the case of
nonresonantly pumped polariton superfluids. We have also seen that the
absence of a long-range wake does not imply the absence of
dissipation, as a residual drag force due to the finite polariton
lifetime is always present in the system. In this sense, one can say
that we are not dealing with superfluid behaviour in a strict sense.


GP equation for the LP branch

If one starts from a description of polaritons in terms of separate
exciton and cavity photon fields, a rotation into the LP and UP basis,
followed by neglecting the occupancy of the upper polariton branch, as
explained in detail in Chapter , results in the
following Gross-Pitaevskii equation (compare to
Eq. eq:mom-GPE) for the LP field in momentum space

 :


  i_t _LP,k = f_p e^-i_p t
  _k,k_p + [_LP () - i    ()/2]_LP,k +


  _k_1, k_2 g_k, k_1, k_2
  ^*_LP,k_1 + k_2-k _LP,k_1
  _LP,k_2 + C_ _k_1 V_d,k - k_1 _LP,k_1 C__1

where  is the effective LP
decay rate,



is the interaction strength, and where
.

In these expressions, the coefficients (see Eqs. eq:hopfield-X
and eq:hopfield-C) 



are the Hopfield coefficients used to diagonalise the free polariton
Hamiltonian. We want here to justify the simplified description done
in Eq. eq:basic. If we follow the linear response expansion as
in eq:mfield, the operator  in momentum
space analogous to eq:opell reads as:



where now 

. It is easy to show that
the eigenvalues of this operator coincide with our approximated
expressions eq:spect in the limit of ,
when , 
 and when we can simply
rename  and 
. It is interesting to note that, even if we would
retain the linear terms in  in the
expansion of , this would result in
a renormalisation of the fluid velocity  in the
expression eq:spect which takes into account the blue-shift of
the LP dispersion due to the interaction.


















Polariton superfluidity in the OPO regime

















Our findings
This Chapter presents a joint theoretical and experimental study of an
OPO configuration (see Fig. ) where a wide and
steady-state condensate hits a stationary localized defect in the
microcavity.  Contrary to the
criterion for quantized flow metastability for which the signal and
idler display simultaneous locked responses, we find that their
scattering properties when the OPO hits a static defect are different.
In particular we investigate the scattering properties of all three
fluids, the pump, the signal and the idler, in both real and momentum
space. We find that the modulations generated by the defect in each
fluid are not only determined by its associated Rayleigh scattering
ring, but each component displays additional rings because of the
cross-talk with the other components imposed by nonlinear and
parametric processes.  We single out three factors determining which
one of these rings has the biggest influence on each fluid response:
the coupling strength between the three OPO states, the resonance of
the ring with the blue-shifted LP dispersion, and the values of each
fluid group velocity and lifetime together establishing how far each
modulation can propagate from the defect.  The concurrence of these
effects implies that the idler strongly scatters, inheriting the same
modulations as the pump, while the modulations due to its own ring can
propagate only very close to the defect and cannot be
appreciated. However, the modulations in the signal are strongly
suppressed, and not at all visible in experiments, because the slope
of the polariton dispersion in its low momentum component brings all
Rayleigh rings coming from pump and idler out of resonance.

    
    
    Pictorial representation of the optical parametric oscillator
    (OPO) state inside a planar microcavity. Cavity modes of the incoming
    external laser [black arrow] couple to the excitons formed by electron
    [e, solid spheres] - hole [h, hollow spheres] pairs
    inside the semiconductor quantum well (QW), creating the pump [red]
    polaritons, which coherently scatter into the signal [blue] and idler
    [green] states. All three types of polaritons decay into external
    radiation, emitted through the top distributed Bragg reflector (DBR)
    of the cavity.
    
  
There is always scattering due to the pump
Note that the kinematic conditions for OPO are incompatible with the
pump and idler being in the subsonic regime. Thus, the coupling
between the three components always implies some degree of scattering
in the signal. In practice, the small value of the signal momentum
strongly suppresses its visible modulations, as confirmed by the
experimental observations.

Model

Exciton-photon GPE
As explained in Chapter , the dynamics of
polaritons in the OPO regime and their hydrodynamic properties when
scattering against a defect can be described via a classical
driven-dissipative non-linear Gross-Pitaevskii equation (GPE) for the
coupled exciton and cavity fields
 :



The dispersive - and -fields decay at a rate  and
are coupled by the Rabi splitting , while the nonlinearity
is regulated by the exciton coupling strength :



We describe the defect via a potential  acting on the
photonic component; this can either be a defect in the cavity mirror
or a localized laser field .

In the conservative, homogeneous, and linear regime [
], the eigenvalues of  are given by the
LP and UP energies Eq. eq:polariton-dispersion.

The cavity is driven by a continuous-wave laser field 

 into the OPO regime: Here, polaritons are continuously injected
into the pump state with frequency  and momentum
, and, above a pump strength threshold, they undergo coherent
stimulated scattering into the signal  and
idler  states. 


 

OPO mean-field blueshifts and fluctuation
  Rayleigh rings in the linear-response scheme for homogeneous
  pumping. Left panel: signal  ([blue] upper triangles), pump 
  ([red] circles), and idler  ([green] lower triangles) mean-field
  energy blueshifts  (in units of 
) vs the rescaled pump intensity  (in
  units of ) in the optical limiter regime. Parameters are
   meV, zero cavity-exciton detuning, 
 meV,  meV, 
, , and .  The shaded
  area is stable OPO region, while the vertical dashed line
  corresponds to the pump power value chosen for plotting the right
  panel. Right panel: blueshifted LP dispersion eq:blues with
  superimposed Rayleigh curves 
 evaluated within the linear-response
  approximation (same symbols as left panel). The two rings
  corresponding to the signal state, 
, are shrunk to zero because .


LP GPE in momentum space
As a first step, it is useful to get insight into the system behaviour
in the simple case of a homogeneous pump of strength
. A numerical study of the coupled
equations eq:gpequ for the more realistic case of a
finite-size top-hat pump profile  will be
presented in Sec. .

To further simplify our analysis, we assume here that the UP
dispersion does not get populated by parametric scattering processes
and thus, by means of the Hopfield coefficients
Eqs. eq:hopfield-X and eq:hopfield-C, we project the
GPE eq:gpequ onto the LP
component 
, where

. As explained in detail in
Sec. , we obtain Eq. eq:mom-GPE:


  i_t _k^ = [LP(k) -
    i_k2]_k^ +
  C_k _q C_q V_d(k - q)
  _q^
 + _k_1, k_2 g_k,
    k_1, k_2 ^*_k_1 + k_2-k
  _k_1^ _k_2^ + f_p(t)
  _k,k_p

Here, 
is the effective LP decay rate, the interaction strength is given by

, and the pumping term is
given by . Note
that the dependence on the exciton-exciton interaction strength 
can be removed by rescaling both the LP field
 and the pump
strength , something we will do later
on to all effects, working in terms of energy blueshifts and the
rescaled pump intensity.


Linear-response theory


Spectrum of collective excitations. Cut at  of the
  real part of the quasiparticle energy dispersion
   plotted versus 
. The spectrum is evaluated within the linear approximation
  scheme and the parameters are the same ones used for
  Fig. .

Linear response ansatz
As already explained in Sec. , in the limit where the
homogeneously pumped system is only weakly perturbed by the external
potential , we apply a linear-response
analysis : The LP field is expanded around
the mean-field terms for the three  OPO
states  (see Eqs. eq:pump-fluctuations
and eq:3-state-ansatz)



where . Eq. eq:efflp is
expanded linearly in both the fluctuation terms,
 and , as well as
the defect potential.  At zeroth order, the three complex uniform
mean-field equations can be solved to obtain the dependence of the
signal, pump and idler energy blueshifts,
 on the system
parameters . A typical behaviour of 
as a function of the rescaled pump intensity
 in the optical limiter
regime is plotted in the left panel of Fig. .

At first order, one obtains six coupled equations diagonal in momentum
space 


for the six-component vector

 and for the potential part,


. To understand the origin of the factor  in
Eq. eq:fluct, it is sufficient to consider one component only,
and start from Eq. eq:mfield used in Sec. . The
connection becomes clear once we write the fluctuation part as

 and compare
it to Eq. eq:pump-fluctuations of Sec. .


Bogoliubov matrix
In  eq:fluct we have only kept the terms oscillating at the
frequencies  and neglected the other terms in the
expansion (i.e., ), which are
oscillating at frequencies far from the LP band, and thus with
negligible amplitudes.

In the particlelike and the holelike channels, the Bogoliubov matrix
determining the spectrum of excitations can be written
as 


where the three OPO states components are


  (M_k^)_mn &= [_LP(k_m + k)
    - _m - i_k_m + k2]
  _m,n
   + 2_q,t=1^3
  g_k_m+k,k_n+k,k_t _q^*
  _t^ _m+q,n+t


  (Q_k^)_mn &= _q,t=1^3
  g_k_m+k,k_q,k_t _q^ _t^
  _m+n,q+t

OPO spectrum
In absence of a defect potential (),
Eq. eq:fluct is the eigenvalue equation for the spectrum of
excitations of a homogeneous OPO, i.e.,
. We plot in
Fig.  a typical collective dispersion (here we consider
the same system parameters as the ones used in Fig. ),
by plotting the real part of the Bogoliubov matrix eigenvalues
 as a function of 
(cut at ). The spectrum has six branches,
, labeled by  and
. Even though these degrees of freedom are mixed together, at
large momenta, one recovers the LP dispersions shifted by the three
states' energies and momenta, i.e.,



where  () corresponds to the  () particlelike (holelike)
branch.

The OPO solution is stable (shaded area in Fig. ) as
far as . Note also that,
owing to the  symmetry of the OPO equations of motion (see
Sec. ), no restoring force would oppose a rotation of the
signal-idler phases, meaning that the generator of such global
rotations,
, is an
eigenvector of  with eigenvalue
 .

Rayleigh rings
The shape of the patterns, or Cherenkov-like waves, resulting from the
elastic scattering of the OPO 3-fluids against the static ()
defect can be determined starting from the spectrum, and in particular
evaluating the closed curves  in
-space, or "Rayleigh rings"  defined by
the condition . Even if
they do not appear to be relevant here, note that the presence of a
non-vanishing imaginary part of the excitation spectrum
 introduces some
complications: Even in the absence of any Rayleigh ring, the drag
force can be non-vanishing (as we saw in Chapter  for
the pump-only case) and the standard Landau criterion may fail to
identify a critical velocity .

The modulations propagate with a direction
 orthogonal to each curve
, a pattern wavelength given by
the corresponding , and a group velocity

, where


 determines the distance, at any given direction
, over which the perturbation
extends away from the defect. For a single fluid under a coherent
pump, the qualitative shape of the modulation pattern generated in the
fluid by the defect is mostly determined by the excitation
spectrum .



Blueshifted dispersion and rings for ks=0
For OPO, the spectrum of excitation on top of each of the three,
, states (see Fig. ) generates six identical
Rayleigh rings  for the three
states. Note that the Rayleigh rings can be found by finding the
intersections .

The Rayleigh rings for the OPO conditions specified in
Fig.  are clearly visible in the right panels of
Fig.  (b), where we plot the -space
photoluminescence filtered at the energy of each state, i.e.,


.

We have chosen here a -like defect potential,
, but we have checked that our results do not
depend on the specific shape of the defect potential: In particular,
we have also considered the response to defects with smooth
Gaussian-like profiles, whose effect is only to partially weaken the
upstream modulations in real space.


As detailed in Sec. , the mean-field OPO ansatz assumes
the signal momentum  to be an input parameter, rather than
solving the full "selection problem" . This is
in contrast to both the full numerical solution of the GP
Eq. eq:gpequ (Sec. ) as well as
experiments (Sec. ), where the signal
momentum is automatically selected (typically close to the bottom of
the LP band) by parametric scattering processes. The reason why
parametric scattering processes select a signal with a momentum close
to zero, already very close to the lower pump thresold for OPO, is
still awaiting an explaination.

linresp


  
  
  Linear response of the three OPO states to a static
  -defect for  [column (a)],  [column (b)] and
   m [column (c)]. The parameters are given in
  Fig. , and m.  Column
    (a): Top row is the analogue of Fig. , for negative
  signal momentum. The three bottom rows show the rescaled emission in
  real  (left panels in linear
  scale) and momentum space 
  (right panels in logarithmic scale) filtered at the energies of the
  signal (top), pump (middle) and idler (bottom).  Column (c): Same
  as column (a), for a positive signal momentum.  Column (b):
  Vanishing signal momentum; the corresponding mean-field blueshifts and
  Rayleigh rings are plotted in Fig. . The inset shows
  the Gaussian-filtered (see main text) signal emission.

In particular, this cannot be addressed within a spatially homogeneous
approximation where the three mean-field solutions for pump, signal,
and idler states are described by plane waves. However, one can
show  that, within the same mean-field
approximation scheme, when increasing the pump power towards the upper
threshold for OPO, the blue-shift of the LP polariton dispersion due
to the increasing mean-field polariton density, causes the signal
momentum to converge towards zero .

In the following, we take advantage of the broad range of possible
choices for  (and thus ) at mean-field
level , and analyze three distinct cases:
negative ( m), vanishing
( m) and positive ( m)
signal momenta.





Vanishing signal momentum

For the OPO conditions considered here, the signal momentum is at
, and thus only four of the six Rayleigh rings are
present. The same rings are also plotted in the right panel of
Fig. , shifted at each of the three OPO states'
momentum ,  and
energies .

It is important to note that, even though the three OPO states have
locked responses because they display the same spectrum of
excitations, only one of the rings 
 is the most resonant at 
with the interaction blueshifted LP dispersion,



where  are the
mean-field energy blueshifts (measured in Fig.  in
units of ).  This implies that the
most visible modulation for each fluid should be the most resonant
one, with superimposed weaker modulations coming from the other two
state rings.

Signal for ks=0
In the specific case of Fig. , the signal is at 
 and thus produces no rings in momentum space. The other four
rings are very far from being resonant with the blueshifted LP
dispersion eq:blues at , and thus the signal
displays only an extremely weak modulation coming from the next
closer ring, which is the one associated with the pump state,
.

We estimate that the signal modulation amplitudes are roughly  of
the average signal intensity and about a factor of 10 times weaker
than the modulation amplitudes in the pump fluid.

To show that the signal has weak modulations coming from
the pump, we apply a Gaussian filter to the real space images (see the
inset of Fig.  (b)).

This consists of the following procedure. The original data for the
real space profile  are convoluted with a Gaussian
kernel , obtaining a new profile,

, where short wavelengths features are smoothened out. The
convoluted image  is then subtracted from the
original data, giving , and
effectively filtering out all long wavelength details.

This procedure reveals that indeed the pump imprints its modulations
also into the signal, even though these are extremely weak, thus
leaving the signal basically insensitive to the presence of the
defect.

Pump and idler for ks=0
Pump and idler states are each mostly resonant with their own
rings, i.e.,  at 
and  at ,
respectively. Thus one should then observe two superimposed
modulations in both pump and idler filtered emissions, the stronger
one for each being the most resonant one.

However, the modulations associated with the idler only propagate very
close to the defect, at an average distance
m before getting
damped, and thus are not clearly visible. For the OPO conditions
considered, this is due to the small idler group velocity
, as the dispersion is almost
excitonic at the idler energy.

Conclusion for ks=0
We can conclude that, for the typical OPO condition with a signal at
, considered in Figs. and  (b), the signal fluid does not show modulations
and the extremely weak scattering inherited from the pump state can be
appreciated only after a Gaussian filtering procedure of the image. In
contrast, the idler has a locked response to the one of the pump
state.

Finite signal momentum


Real space signal, pump and idler OPO one-dimensional
  filtered profiles in the linear-response approximation. OPO filtered
  emissions along the  direction, 
.  A signal at 
 m [panel (a)] corresponds to Fig.   (a), a signal at  m [panel (b)] corresponds to
  Fig.  (b) and a signal at  m
  [panel (c)] corresponds to Fig.  (c). In each panel we
  plot the filtered profiles of signal, pump, and idler, while the
  horizontal gray dashed lines represent the values of the mean-field
  emission without a defect.

Given the freedom of choice for the signal momentum  close
to the lower OPO threshold, we consider here two additional cases,
that could not be studied neither experimentally, nor within a full
numerical approach, but that instead we can easily analize within the
linear-response theory.

In particular, we have left fixed the pumping conditions
 m and  meV and
considered two opposite situations.

In the first case, the signal has a finite and positive momentum
 m and thus the idler is at low momentum,
 m. The results are shown in
Fig.  (c).

Here we see that all six Rayleigh rings are clearly visible and, in
addition, as the idler is at lower momentum compared to the case
considered in Sec. , and thus its
dispersion steeper, the idler group velocity is large enough to
appreciate the modulation of the Rayleigh ring associated to this
state. As a result, each of the three filtered OPO emissions exhibits
as the strongest modulation the one coming from its own Rayleigh ring,
including the signal which is now at finite and large momentum. In
this case, the OPO response of each filtered state profile looks
completely independent from the other, as if we were pumping each
state independently.

In the second case, shown in Fig.  (a), the signal is
finite and negative,  m, and the idler is now
at very large momentum,  m, where its
dispersion is very exciton-like and flat, and thus the idler has a
very small group velocity and its own modulations are visible only
very close to the defect. For this case, we can appreciate in the
idler filtered profile overlapped modulations from all the three state
Rayleigh rings (note that because the signal is at negative momentum,
its modulations have an opposite direction compared to the ones of
pump and idler), while in the signal we can mostly see the signal long
wavelength modulations and only very weakly the pump one.



Real space profiles of three uncoupled signal, pump and idler
  fluids. The rescaled profiles 
  are obtained by setting all off-diagonal couplings in
  Eq. eq:bogoliubov-opo to zero, resulting in three uncoupled
  signal (top row), pump (middle row), and idler (bottom row)
  fluids. The three columns correspond to the three different cases
  analysed within the linear-response approximation: the case
  of a signal at  m (left column) corresponds to
  the same conditions as Fig.  (a), a signal at 
 m (middle column) corresponds to Fig.   (b), and a signal at  m (right column)
  corresponds to Fig.  (c).


Discussion
We can compare the different modulation strengths of the three OPO
profiles by looking at the color bars plotted next to the profiles. In
order to better compare them on the same plot, we show in
Fig.  the one-dimensional OPO filtered emissions along
the  direction, rescaled by the mean-field solution  in
absence of the defect, for the three cases analysed above. While for
the OPO conditions with a signal at  m (middle
panel, corresponding to Fig.  (b)), the imprinted
modulations from the pump are hardly visible, and can only be
appreciated after a Gaussian filter manipulation, both OPO cases with
a finite momentum signal result in modulations in the signal with an
amplitude of the same order of magnitude of both pump and idler
fluids.

In Fig.  we instead show the real space profiles
obtained by setting all off-diagonal couplings in the Bogoliubov
matrix  of Eq. eq:bogoliubov-opo to
zero, resulting in three uncoupled signal (top row), pump (middle
row), and idler (bottom row) fluids. This underlines the importance of
the coupling between the three fluids in the three different OPO
regimes we analysed within the linear-response approximation. In
particular, for the OPO conditions such that the signal has a finite
and positive momentum  (right column of
Fig. , which corresponds to the conditions shown in
Fig.  (c)), the coupling has little effect and the
three fluids respond to the defect in practice in a independent way
(each Rayleigh ring influences its own fluid). The OPO condition with
a finite and negative momentum , which corresponds to the
conditions shown in Fig.  (a), are shown in the left
column of Fig. : Here, we can see that the coupling
between the three fluids plays a role. The modulations of the pump can
be appreciated both in the signal profile (though weakly), as well as
in the idler profile. In the idler fluid its own modulations can only
propagate very close to the defect and thus the only really
appreciable modulations are the ones inherited from the pump. Finally,
for the experimentally relevant case of an OPO with a signal at zero
momentum (middle column of Fig. , which corresponds to
the conditions shown in Fig.  (b)), there is also a
role played by the coupling, though, as already thoroughly analysed,
the modulations in the signal inherited from the pump can only be
appreciated after a Gaussian filtering manipulation.

Finally note that, for the OPO conditions of
Sec. , as well as for the finite momenta
cases considered in Sec. , the subsonic to
supersonic crossover of the pump-only state  happens at
pump intensities well above the region of stability of OPO - the
shaded gray regions of Fig.  and Fig. (a) and (c). Thus, it is not possible to study a case where the pump
is already subsonic and at the same time promotes stimulated
scattering.

Numerical analysis


OPO spectrum obtained by full numerics. Left
  panel: Photonic component of the OPO spectrum in presence of a
  point-like defect,  (logarithmic scale),
  as a function of the rescaled energy  versus
  the -component of momentum  (cut at ) for a top-hat
  pump (see text for the space profile and parameter values), with
  intensity  above the OPO threshold, pump
  wave-vector  m in the -direction and
   meV. The symbols indicate the signal
  ([blue] upper triangle), pump ([red] circle), and idler ([green]
  lower triangle) energies, as well as the two momenta  on each
  state Rayleigh ring at . Note that the logarithmic scale
  results in a fictitious broadening in energy of the spectrum, which
  is in reality -like (see right panel). The bare LP
  dispersion, including its broadening due to finite lifetime, is
  plotted as a shaded grey region, while the bare UP dispersion as a
  (black) dot-dashed line. Right panel: Momentum integrated spectrum,
   (linear scale) as a
  function of the rescaled energy , where it can
  be clearly appreciated that the emission is -like in
  energy.

The results obtained within the linear-response approximation are
additionally confirmed by an exact full numerical analysis of the
classical driven-dissipative non-linear GPE eq:gpequ for the
coupled exciton and cavity fields  in the case
of a finite size pump. Eq. eq:gpequ is solved numerically on a
2D grid of  points and a separation of
 m (i.e., in a box 
 m m) by using a 5-order
adaptive-step Runge-Kutta algorithm. Convergence has been checked both
with respect the resolution in space  as well as in momentum
, without  as well as in
presence of the defect.

The same approach has been already used in the literature, for a
review, see Refs. ).

As for the system parameters, we have considered a LP dispersion at
zero photon-exciton detuning, , a
dispersionless excitonic spectrum, 
and a quadratic dispersion for photons
, with the photon mass
, where  is the bare electron
mass. The LP dispersion Eq. eq:polariton-dispersion is
characterised by a Rabi splitting  meV. Furthermore,
the exciton and cavity decay rates are fixed to
 meV.

For the defect we choose a -like potential



where its location  is fixed at one of the 
points of the grid.

Note that in a finite-size OPO, local currents lead to inhomogeneous
OPO profiles inside the pump spot, despite the external pump having
a top-hat profile with a completely flat inner
region  - as shown later, this
can be observed in the filtered OPO profiles evaluated in absence of a
defect, shown as dashed lines in Fig. .

We have thus chosen the defect location so that it lies in the
smoothest and most homogeneous part of the OPO profiles.

Also we have checked that our results do not qualitatively depend on
the strength  (nor on the sign) of the defect potential, as far
as this does not exceed a critical value above which it destabilises
the OPO steady-state regime.

The pump, 
, has a smoothen and rotationally
symmetric top-hat profile, 

 with strength 
 meV/m and parameters m,
m.

We pump at  m in the -direction, 
, and at  meV, i.e., roughly
 meV above the bare LP dispersion. By increasing the pump
strength , we find the threshold  above which
OPO switches on, leading to two conjugate signal and idler states.

We then fix the pump strength just above this threshold (
), where we find a steady state OPO solution which is
stable (see Ref.  for further details). In
absence of the defect, this condition corresponds to a signal state at
 m and  meV and an
idler at  m and 
 meV. 

It is interesting to note that already very close to the lower pump
power threshold for OPO, the selected signal momentum is very close to
zero. This contrasts with what one obtains in the linear approximation
scheme, where instead just above the lower OPO threshold there exists
a broad interval of permitted values for  (and thus ) - we
have already mentioned this "selection problem" for parametric
scattering in Sec. .

We evaluate the time dependent full numerical solution
of eq:gpequ , until a steady state
regime is reached. Here, both its Fourier transform to momentum
 and energies  can be evaluated numerically.

We plot on the left panel of Fig.  a cut at 
of the photonic component of the OPO spectrum in presence of a
point-like defect, , as a function of the
rescaled energy  versus the -component of
momentum  (cut at ). In the right panel we plot instead
the corresponding momentum integrated spectrum, 
.

Here, we can clearly see that the presence of the defect does not
modify the fact that the OPO emission for the OPO signal ([blue] upper
triangle), pump ([red] circle), and idler ([green] lower triangle)
states has a completely flat dispersion in energy, thus indicating
that a stable steady state OPO solution has been reached. Note that in
the spectrum map of the left panel of Fig. , the
logarithmic scale results in a fictitious broadening in
energy. However, from the integrated spectrum plotted in linear scale
in the right panel of Fig.  one can clearly appreciate
that this emission is -like, exactly as it happens for the
homogeneous OPO case .




Full numerical responses to a static defect of
  the three OPO states in real and momentum space. Filtered OPO
  emissions (signal [top panels], pump [middle], and idler [bottom])
  in real space  (left panels in
  linear scale) and momentum space 
 (right panel in logarithmic scale)
  obtained by a full numerical evaluation of eq:gpequ. For the
  top left panel of the signal space emission, Gaussian filtering is
  applied to enhance the short wavelength modulations of this state,
  revealing that the modulations corresponding to the pump state are
  also imprinted (though weakly) into the signal. The symbols indicate
  the pump ring diameter extracted from fitting the upstream
  modulations and resulting in a density-wave wavevector coinciding
  with that of the pump,  m.


Real space signal, pump and idler OPO one-dimensional
  filtered profiles derived from finite size numerics with and without
  a defect. Rescaled OPO filtered emissions along the  direction,
  , obtained by
  numerically solving the GPE Eq. eq:gpequ of
  Sec. . While the dashed lines represent the
  filtered emissions of signal (top panel), pump (middle) and idler
  (bottom) for a top-hat pump without a defect, the solid lines are
  the same OPO conditions but now for a defect positioned at
   m corresponding to the vertical
  dotted lines. The system parameters are the same ones as those of
  Fig. .

Thus the effect of the defect is to induce only elastic (i.e., at the
same energy) scattering; now the three OPO states emit each on its own
Rayleigh ring (given each by the symbols on the left panel of
Fig.  which represent the rings at a cut for
). This makes it rather difficult to extract the separated
signal, pump and idler profiles by filtering in momentum, as done
previously for the homogeneous case, but it still allows to filter
those profiles very efficiently in energy. In fact, because the
emission is -like, it is enough to fix a single value of the
energy  to the one of the three states ,
thus extracting the filtered profiles either in real space
 or in momentum space
 - we have however checked that
integrating in a narrow energy window around  does not
quantitatively change the results.

The results of the above described filtering are shown in
Fig. , where real-space emissions
 are plotted in the left panels, while
the ones in momentum space  are
plotted in the right panels. We observe a very similar phenomenology
to that one obtained in the linear approximation shown in
Fig.  (b). The signal now is at slightly negative
values of momenta  m, thus implying a very small
Rayleigh ring associated with this state. Thus we observe that only
the modulations associated with the pump are the ones that are weakly
imprinted in the signal state and that can be observed by means of a
Gaussian filtering (inset of the top-left panel). We have fitted the
upstream wave crests and obtained the same modulation wavevector as
the pump one ([blue] upper triangles). Similar to the linear-response
case, we also find here that the most visible perturbation in the
emission filtered at the idler energy is the one due to the pump
Rayleigh ring. As before, the modulations due to the idler Rayleigh
ring cannot propagate far from the defect because of the small group
velocity associated with this state.


In Fig.  we instead plot the corresponding
one-dimensional profiles in the  direction both in presence
(solid line) and without (dashed line) a defect. Here, we can observe
that, even if the pump has a top-hat flat profile, as also commented
previously, in absence of a defect, the finite-size OPO is
characterised by inhomogeneous profiles of signal, pump and idler
because of localised currents. Furthermore, we observe that the
presence of a defect induces strong modulations in pump and idler.

Notice that the modulation imprinted by the pump into the signal is
hardly visible in the top panel of Fig.  without the
Gaussian filtering.


Experiments



    
  Experimental OPO spectrum (lower panel) and filtered
    emissions of signal (top), pump (middle), and idler (bottom) in
    presence of a structural defect. The six panels show the filtered
    emission profiles in real space . A Gaussian
    filtering to enhance the short wavelength modulations is applied
    in the right column. The extracted wave crests from the idler
    (yellow contours in the bottom panel) are also superimposed on the
    pump profile (middle) by applying a  phase-shift.

We now turn to the experimental analysis, where a continuous-wave
laser is used to drive a high quality () GaAs microcavity
sample into the OPO regime - details on the sample can be found in
Refs. .

The polariton dispersion is characterised by a Rabi splitting
 meV and the exciton energy is
 meV, while the cavity-exciton detuning is
slightly negative,  meV. The pump is at m and
 meV, and, at pump powers -times
above threshold, an OPO appears, with signal at small wavevector
m and  meV, and
idler at m and
 meV.

The defect used in the sample is a localized inhomogeneity naturally
present in the cavity mirror. Note that the exact location of the
defect can be extracted from the emission spectrum and is indicated
with a dot (orange) symbol in the profiles of Fig. .

To filter the emission at the three states' energies,
, and to obtain 2D spatial maps for the three
OPO states, one uses a spectrometer and, at a fixed position ,
one obtains the intensity emission as a function of energy and
position, . By changing  one can build the
full emission spectrum as a function of energy and 2D position,
. The filtered emission for each OPO state is
obtained from the integrals

, with  meV. The results
are shown in Fig.  for, respectively, the signal (top
panel), pump (middle) and the idler (bottom) profiles. Energy and
momentum of the three OPO states are labeled with a [blue] upper
triangle (signal), a [red] circle (pump), and a [green] lower triangle
(idler), while the localised state, clearly visible just below the
bottom of the LP dispersion, is indicated with the symbol . The
bare LP dispersion is extracted from an off-resonant low pump power
measurement, as well as the emission of the exciton reservoir (X) and
that of the UP dispersion (each in a different scale).

The signal profile shows no appreciable modulations around the defect
locations, nor could any  be observed after applying a Gaussian
filtering procedure to the image.

In contrast, in agreement with the theoretical results, both filtered
profiles of the pump and idler show the same Cherenkov-like pattern. We
extracted the wave crests from the idler profile ([yellow] contours in
the bottom panel) and superimposed them on the pump profile (middle
panel) with an added  phase-shift, revealing that the only
modulations visible in the idler state are the ones coming from the
pump state.


Conclusions

To conclude, we have presented a joint theoretical and experimental
study of the superfluid properties of a nonequilibrium condensate of
polaritons in the so-called optical parametric oscillator
configuration by studying the scattering against a static defect.

We have found that, while the signal is basically free from
modulations, the pump and idler lock to the same response. We have
highlighted the role of the coupling between the OPO components due to
nonlinear and parametric processes. These are responsible for the
transfer of the spatial modulations from one component to the
other. This process is most visible in the clear spatial modulation
pattern that is induced by the nonsuperfluid pump onto the idler,
while the same modulations are only extremely weakly transferred into
the signal, because of its low characteristic wavevector, so much that
experimentally cannot be resolved.

The main features of the real- and momentum-space emission patterns
are understood in terms of Rayleigh scattering rings for each
component and a characteristic propagation length from the defect; the
rings are then transferred to the other components by nonlinear and
parametric processes.

Much interest has been recently devoted to aspects related to
algebraic order  and superfluid
response  in drive-dissipative polariton
condensates. Our theoretical and experimental results further stress
the complexities and richness involved when looking for superfluid
behavior in nonequilibrium multicomponent condensates such as the ones
obtained in the OPO regime.





















Landau levels in driven-dissipative cavity arrays

In 1984, Michael Berry showed that the adiabatic evolution of energy
eigenfunctions with respect to a time-dependent Hamiltonian contains a
phase of geometrical origin, commonly known today by the name of
"Berry's phase".  Although Berry's seminal paper gave
the example of a spinor's evolution under slowly changing magnetic
field, geometric phases with similar
origin have been encountered in a plethora of physical contexts,
ranging from hydrodynamics to quantum field theory, the quantum Hall
effect and topological insulators.

In a condensed matter context, the eigenstates of electrons in a
periodic lattice can be labeled by a band index  and the crystal
momentum . The Berry-phase physics arising when adiabatically
transitioning between neighbouring wavevectors was appreciated only
recently. This phenomenon can explain properties such as electric
polarization, anomalous Hall conductivity or the quantization of
conductance in the integer quantum Hall effect.

In the presence of a (synthetic) gauge field, the eigenstates making
up an energy band can have nontrivial geometrical properties, as
encoded in the Berry connection and Berry curvature 
. Understanding the geometry of eigenstates in a
band is of great importance, not least because the integral of the
Berry curvature over the 2D Brillouin zone (BZ) gives the first Chern
number: the topological invariant responsible for the integer quantum
Hall effect . Consequently, there has been much work in
recent years to develop new techniques with which to probe the
properties of energy bands in photonics and ultracold gases. For
example, the Berry curvature can be measured in the semiclassical
dynamics of a wavepacket in an optical lattice 

 or in photon transport in a cavity
array . In all these cases, the physics can be most
naturally understood by recognising that the Berry curvature acts like
a magnetic field in momentum space 
.

The analogy between Berry curvature and magnetism is most powerful
when a geometrical energy band is subjected to an additional weak
harmonic potential . Then, in the effective
momentum-space Hamiltonian, a harmonic potential acts like the kinetic
energy of a particle in real space. Just as the physical momentum
 is the sum of the canonical momentum  and
the magnetic vector potential  in the magnetic
Hamiltonian, so the physical position ,
is given by the canonical position  and the Berry connection
 of band  in the effective
Hamiltonian 

. The Berry curvature,
, is then like a
momentum-space magnetic field. For certain models, this analogy leads
to a clear analytical understanding of single-particle
dynamics 
. In particular, we will focus on the small-flux
limit of the Harper-Hofstadter
Hamiltonian , which
is a model that has recently been realized in a multitude of
experimental configurations, ranging from ultracold
gases ,
solid state superlattices and silicon photonics  to classical systems
such as coupled pendula  and oscillating
circuits . As first shown in
Ref. , the eigenstates of this model in the
presence of a harmonic trap are toroidal Landau levels in momentum
space. Not only would an observation of these states constitute the
first exploration of analogue magnetic states in momentum space, but
also the first experimental study of magnetism on a torus.
 
Most theoretical works on momentum-space Landau levels have focused on
conservative dynamics ,
while photonic systems naturally include driving and
dissipation . In this Chapter, we present a
realistic experimental proposal for the observation of these states in
a driven-dissipative 2D lattice of cavities, such as the array of
coupled silicon ring resonators of Ref. ,
where link resonators were used to simulate a synthetic gauge field
for photons. We combine this set-up with a harmonic potential,
introduced, for example, by a spatial modulation of the resonator
size. We demonstrate numerically that the main features of
momentum-space Landau levels will be observable spectroscopically in
this system for realistic parameters.

In this Chapter, we also emphasize how the inherent driving and
dissipation in photonics can be a key advantage in probing properties
that are otherwise inaccessible.  Firstly, the spectroscopic
measurements discussed here are sensitive to the absolute energy of a
state. From this, we show how to extract the energy shift due to the
off-diagonal matrix elements of the Berry connection
 relating eigenstates in different bands 
and . Only very recently has the first measurement of such effects
been reported in ultracold atomic
gases , and the approach
used in this experiment would be difficult to apply to a photonics
set-up. The scheme we present may therefore be useful for the
characterisation of energy bands in topologically-nontrivial photonic
systems.

Secondly, since the photon steady-state depends on the overlap between
the (observable) spatial amplitude profile of the drive and of the
eigenstates , the observables will depend on
the phase of the eigenfunctions and thus on the specific synthetic
magnetic gauge that is implemented in a given experimental realization
of the Harper-Hofstadter Hamiltonian using a synthetic gauge field. We
note that a related gauge-sensitivity has also recently been of much
interest in ultracold gases in suitably designed time-of-flight
experiments 
. Experiments on synthetic magnetic fields therefore
present the opportunity of straightforwardly probing gauge-dependent
physics.

This Chapter is organized as follows: in Section  we
add a harmonic trap to the HH model introduced in
Sec. . We explain how the eigenstates can be
understood as momentum-space Landau levels in Section
, before we discuss the breakdown of
approximations in Section , focusing on the
energy shift from the off-diagonal matrix elements of the Berry
connection. Then in Section , we add
driving and dissipation to the model. In Section ,
we show numerical results highlighting gauge-dependent effects, before
presenting a viable proposal for a photonics-based experiment in
Section . Finally, we draw conclusions in Section
.


The trapped Harper-Hofstadter Model

In this Chapter, we study the Harper-Hofstadter Hamiltonian
 in the presence of an external harmonic trap. The full
tight-binding Hamiltonian  of this system is

 H &= H_0+12_m,n[(m-m_0)^2+(n-n_0)^2]a_m,n^a_m,n 

H_0 &= -J_m,n(e^i
_m,n^xa_m+1,n^a_m,n +e^i
_m,n^ya_m,n+1^a_m,n) +
H.c. 
where  is the real hopping amplitude and 
() are the creation (annihilation) operators for a
particle on a square lattice at site . The harmonic trap is of
strength  and is centered at a position  which, in
general, need not coincide with a lattice site. Throughout, the
lattice spacing is set equal to one.

As explained in Sec. , the hopping phases
 are the Peierls phases gained by
a charged particle hopping in the presence of a perpendicular magnetic
field . The sum of
the phases around a square plaquette of the lattice is therefore equal
to , where  is the number of magnetic flux quanta
through the plaquette (with ). For neutral particles, such
as photons, these phases can be imposed artificially to simulate the
effects of magnetism, for example, by inserting link resonators into
an array of silicon ring resonators as mentioned
above .

Although the sum of phases around a plaquette is set by the external
(synthetic) flux, the exact form of the hopping phases themselves
depends on the choice of magnetic gauge. In the Landau gauge, for
example,  such that only the hopping
amplitude along one direction is modified. Conversely, in the
symmetric gauge,  and so hopping
terms along both  and  are affected, preserving the 
rotational invariance of the lattice. This gauge-dependence of the
hopping phases is reflected in the spatial profile of the phase of an
eigenstate of . In a photonics experiment, this phase is
an observable quantity as the intensity response of a system to a
given external driving is determined by the overlap of the spatial
amplitude distribution of the pump with the eigenstates. Such
experiments will therefore be sensitive to the synthetic magnetic
gauge as we discuss in Section .

Toroidal Landau Levels in Momentum
Space
Having introduced the full Hamiltonian in Eq. eq:model, we now
review how the eigenstates of this model in an appropriate limit can
be understood as toroidal Landau levels in momentum
space . Throughout the following discussion,
we assume that the trap is centered at the origin .

We begin from the eigenstates of the Harper-Hofstadter model

, where  is the energy dispersion of
band  at crystal momentum . As the spatially-dependent
hopping phases in  break translational invariance, new
magnetic translation operators must be introduced to define a larger
magnetic unit cell, containing an integer number of magnetic flux
quanta . Then
translational symmetry is restored and Bloch's theorem can be applied
to write the eigenstates as 
, where
 is the periodic Bloch function and  is the
number of lattice sites. Thanks to the new larger unit cell, the
crystal momentum and the periodic Bloch functions here are defined in
the smaller magnetic Brillouin zone (MBZ). For example, hereafter, we
take the number of flux quanta per plaquette to be of the form
, where  is an integer. Then the magnetic unit cell can
be chosen to be  times larger than the original unit cell, while
the MBZ is  times smaller than the original BZ.

Adding the harmonic trap breaks all translational symmetry of the
lattice, but we can use the eigenstates of  as a basis
in which to expand the new wave function 

. Substituting this expansion into the full
Schrodinger equation 
, it can be shown that the expansion coefficients
 satisfy :

 i _t _n(k) = E_n(k) +
2_n^',n^"(_n,n^'i
_k + A_n,n^'(k))
 (_n^',n^"i_k +
A_n^',n^"(k)) _n^"(k)  where 
 is the
matrix-valued Berry connection.

To proceed, we consider the harmonic trap to be sufficiently weak
compared to the bandgap that we can make a single-band
approximation . This assumes that only one
coefficient  is non-negligible and that the external trap does
not significantly mix different energy bands. Then Eq. eq:first
reduces to

 where we have introduced the effective momentum-space
Hamiltonian



For the moment we focus on the first two terms; we discuss the role of
the last term, which comes from the off-diagonal matrix elements of
the Berry connection, in detail in the next subsection.  As can be
seen, there is a close analogy between the first two terms in the
momentum-space Hamiltonian and that of a charged particle in an
electromagnetic field in real space:



In this analogy, the role of the particle mass  is played by
, while the scalar potential  is replaced by
the energy band dispersion  and the magnetic vector
potential  by the intra-band Berry connection
. We note that both the magnetic
vector potential and the Berry connection are gauge-dependent
quantities. We hereafter refer to the gauge choice for the Berry
connection as the Berry gauge, and the gauge choice for a real-space
magnetic vector potential as the magnetic gauge. From the Berry
connection, we can also define the geometrical Berry curvature

, which acts like a momentum-space magnetic field
.

The topology of momentum space also plays a crucial role here, as the
MBZ is topologically equivalent to a torus. One important consequence
of this is of course that the integral of Berry curvature over the
whole MBZ is quantised in units of the first Chern number
. In the above analogy with magnetism, this means that
the particle is confined to move on the surface of a torus, while the
Chern number counts the number of magnetic monopoles contained
inside .

The above analogy with magnetism is particularly powerful because
there are natural limits for our model Eq. eq:model in which the
eigenstates of Eq. eq:mag and hence of Eq. eq:dual are
known analytically 
. We will focus on the flat-band limit in which the
bandwidth is much smaller than the trapping energy; for the energy
bands of  with , this assumption improves
as  decreases or as  increases. In this limit, we can
firstly approximate , so that
the first term is analogous to the kinetic energy of a particle in a
uniform magnetic field. Secondly, we can approximate 
 so that the second term of  is
just a constant energy shift. Hence, the corresponding eigenstates can
be understood as toroidal Landau levels in momentum space. We note
that the opposite limit, in which the trapping energy is small
compared to the bandwidth, also yields very interesting physics
including the realisation of a Harper-Hofstatder model in momentum
space .

As shown in Ref. , the momentum-space toroidal
Landau levels form semi-infinite ladders of states:

 where we have introduced the Landau level quantum
number  , and where 
 can be
recognised as the analogue of the cyclotron frequency 
. Again, we note that here we have neglected the contribution
from the last term in Eq. eq:dual, as this will be discussed in
the next subsection. As can be shown from the Diophantine equation for
the Hall conductivity, for odd values of , the Chern number of all
bands except the middle band is 
 . Then as the Chern number is
related to the uniform Berry curvature as 
, where  is the
MBZ area, the Berry curvature is given by 
 .

While the above spectrum does not directly depend on the toroidal
topology of momentum space, the topology does enter into the
eigenstate degeneracy, which is equal to , as well as
into the analytical form of the eigenstates in the MBZ. For example,
for the bands with , the eigenstates can be
written as 
 _(k) &= N_^l__n
_j = - ^ e^- i k_y j  e^ - ( k_x + j
l__n^2 ) ^2 / 2 l__n^2
 &H_( k_x
/ l__n + j l__n)

 N_^l__n &= (
2/q 2^! 2 l__n^2
)^1/2

where  are the Hermite polynomials and 
 is the analogue of the magnetic length. Here we
have taken the Berry gauge to have a Landau form
 parallel to
the  unit vector in momentum space. We have also
assumed that the MBZ is of length  in one momentum direction,
and  in the other direction, corresponding to a magnetic
unit cell of  plaquettes containing one flux quantum. While this
choice of MBZ is valid in any magnetic gauge of the underlying
Harper-Hofstadter model, it is particularly natural when the hopping
phases in  are in the Landau magnetic gauge 
, as we shall discuss further below.


The Berry connection
We now study the effects of the last term in the momentum-space
Hamiltonian Eq. eq:dual which comes from the off-diagonal matrix
elements of the Berry connection :
 E_n(k) &2_n^'n
A_n,n^'(k)^2 
 &=
2_n^'n
u_n,k_kH_0(k)u_n^',k^2[E_n^'(k)
- E_n(k)]^2
   This can be recognised as a momentum-space counterpart of
the real-space geometrical scalar potential previously studied in
atomic systems . In these
systems, the scalar potential arises from real-space Berry
connections, which can be created, for example, by using
spatially-dependent optical fields to dress the atoms.
 
In the flat-band limit, we have checked numerically that we can
approximate the momentum-space geometrical scalar potential as 
 and so this contributes only a uniform
constant energy shift. This term was not considered in previous
works , as these works focused
on systems such as ultracold atomic gases, where the absolute energy
is not easily experimentally observable. As we will discuss in the
next section, spectroscopic measurements in photonics are sensitive to
the energy of a state, therefore such corrections may be extracted
experimentally. We note that an analogous effect has also been derived
for the effective momentum-space magnetic Hamiltonian for a trapped
particle in an ideal flat band , and predicted
for the frequency spectrum associated with excitonic states in
transition metal dichalcogenides .

To see under what conditions the off-diagonal elements of the Berry
connection are relevant, we compare the energy 
obtained from an exact numerical diagonalization of Eq. eq:model
with the analytical eigenenergy  predicted by
Eq. eq:ladders. We focus on the lowest ladder of states,
associated with band  in Eq. eq:ladders, at energies which
are below the onset of the second ladder around energy . This
allow us to easily identify which numerical eigenvalue should
correspond to which Landau level quantum
number .

We introduce two dimensionless parameters  and
 to quantify deviations between the numerics and
analytics. The former represents the error in the "zero-point
energy", and is defined as the energy of the lowest numerical state
relative to the analytical  Landau level. The latter is the
level spacing error, which we define as the difference between the
numerical energy spacing between two neighbouring states and the
analytical spacing between states with  and  quantum
numbers.

Considering first the "zero-point energy" error ,
the analytical energy of the  Landau level is given from
Eq. eq:ladders by:



where we have used that  and where
we calculate the uniform energy shift
 as the average band-energy
over the MBZ. This definition generalises our flat-band approximation
to account for the non-zero bandwidth of the lowest band.

We then define the dimensionless parameter  as
 This dimensionless error is plotted with a dashed line
as a function of  for  in the top panel of
Fig. . At small , there is a large bandgap between the
lowest two Harper-Hofstadter energy bands, but the lowest band also
has a large bandwidth. In this regime, the single-band approximation
is reasonable, while the flat-band approximation breaks down leading
to large errors. This limit requires a different analytical approach
as previously presented in Ref. .  To account for
the error at large , we include the shift from the off-diagonal
matrix elements of the Berry connection Eq. eq:shift. We
calculate this shift numerically from the eigenstates of the
Harper-Hofstadter model, and we incorporate it into a second parameter
 This is plotted as a solid line in the top panel of
Fig. . As can be seen, this shift dramatically reduces
the error in the zero point energy at large . We also calculate
this shift considering only the effects of band mixing with the second
lowest band ; this is indistinguishable on this scale from the
full shift. This can be understood from the dependence on the bandgaps
in Eq. eq:shift, which shows that the contributions of
high-energy bands are suppressed. As discussed further in
Section , it would be possible experimentally to extract
the energy of the lowest state; this could constitute the first direct
measurement of the effects of the off-diagonal matrix elements of the
Berry connection in a photonics system.

    
  Top panel: "Zero-point energy" error, with (solid
curve, ) and without (dashed line,
) the shift from the off-diagonal matrix elements
of the Berry connection. Including just the first term in the sum
Eq. eq:shift gives an identical curve to the one of
. The chosen trap strength is 
.  Bottom panel: Level-spacing error, for the same
trap strength, considering  (solid curve), 
(dashed curve),  (dashed-dotted curve) and 
(dotted curve).
  
We turn now to the level-spacing error . This can
be expressed as
 where we have used that the analytical level spacing
from Eq. eq:ladders is simply .  We plot
the level-spacing error in the bottom panel of Fig.  for
 and 4. As can be seen here, there is a large
variation in the errors at small  due to the large
bandwidth . On the other hand, we see that
, for , where the flat-band
approximation improves. In this regime, the level-spacing error is
much smaller than the zero-point error. This is because when the shift
from the off-diagonal matrix elements of the Berry connection
Eq. eq:shift is approximately uniform over the MBZ at large
, it just acts as a uniform energy shift on all the states in a
ladder with band index . Consequently, this shift drops out of the
level spacing error between states, leaving only higher-order
band-mixing terms. From perturbation theory, it is expected that
mixing with other bands leads to a negative energy shift on states in
the lowest band, and indeed this can be seen in both
 and  in the small
negative errors found at large .

Driving and dissipation
We now include in our model the driving and dissipation that are an
integral part of the proposed photonics experiment. We assume there
are uniform and local losses characterized by a loss rate ,
and that the pump is monochromatic, with frequency  and a
spatial profile . Following the treatment of
Ref. , we replace the bosonic creation and
annihilation operators with their expectation values, as can be
justified for a noninteracting system. The steady state evolution of
the photon-field amplitude in a cavity then follows that of the pump
as . Combining Hamiltonian
evolution with pumping and losses, one arrives at a set of linear
coupled equations that can be solved numerically for the
steady-state:

 f_m,n
=J[e^-i_m,n^xa_m+1,n+e^i_m-1,n^xa_m-1,n
. 

. +e^-i_m,n^ya_m,n+1+e^i_m,n-1^ya_m,n-1]

 +[_0+i12((m-m_0)^2+(n-n_0)^2)]a_m,n
 where we have reintroduced the position of the harmonic
trap centre , although unless otherwise specified we set
 in our simulations.

The expectation values  correspond to the number of
photons at site , whereas the intensity spectrum is given by
their total sum  as a function of pump
frequency . These observables can be directly related to the
eigenstates of the Hamiltonian in Eq. eq:model. Firstly, the
different eigenmodes of a driven-dissipative system will appear as
peaks in the transmission and/or absorption spectra under a coherent
pump . The resonance peaks will be broadened
by the decay rate , while the area of the peaks will depend on
the overlap between the spatial amplitude profile of the pump and the
underlying eigenstate of  at that energy.

Secondly, when the pump frequency is set on resonance with a given
mode, the intensity profiles in both real- and momentum-space
reproduce the wave function of that
mode . This corresponds respectively to
measuring the near-field and far-field spatial emission of photons
from the cavity array. We note that the far-field emission is simply
the Fourier-transform of the real-space wave function and so will be a
function of crystal momentum defined in the full BZ. To reach the MBZ,
a further processing step is required; for example, if the
Harper-Hofstadter Hamiltonian is in the Landau gauge and if we choose
a magnetic unit cell of  plaquettes along , the
appropriate transformation takes a particularly simple
form :



where  is the wave function coefficient
in the MBZ, while  is that in the original
BZ. In this expression,  is an integer, while
 is the magnetic reciprocal lattice
vector, where the factor of  is due to the enlarged magnetic unit
cell. We note that for other magnetic gauges or for other choices of
the magnetic unit cell, this transformation will in general be more
complicated. In this sense, we call this choice of magnetic unit cell,
a "natural" choice when the Harper-Hofstadter Hamiltonian is in the
Landau gauge.  In the rest of the Chapter, we denote the momentum in
the original BZ as , and that in the MBZ as
.

Results and discussion

Pumping and gauge-dependent effects

    
  Intensity spectra for different pumping conditions: (a)
pumping the single site (5,5), (b) pumping with a Gaussian profile
centered at site (5,5) with width , (c) homogeneous pumping
across all lattice sites and (d) pumping with a random phase across
all lattice sites. These results were obtained by numerical solving
Eq. eq:linear_problem for the steady-state in a lattice of 
 sites, with , 
 and .  Black (solid) curves correspond to
using the Landau gauge while green (dashed) ones to the symmetric
gauge. The dotted vertical lines (with labels indicating the value of
) mark the states which were selected for later analysis. The
spectra in panels (a) and (d) are identical for both gauges.
  
As introduced above, spectroscopic measurements can be used in a
driven-dissipative photonics experiment to study the trapped
Harper-Hofstadter model and hence toroidal Landau levels in momentum
space. In this section, we focus on the effects of the pumping,
exploring how different pumping schemes excite the eigenstates with
different weights. We find that such spectroscopic measurements are
sensitive to the underlying synthetic magnetic gauge chosen in a given
implementation of the Harper-Hofstadter Hamiltonian.

To best illustrate these gauge-dependent effects, we present the
results of numerically solving Eq. eq:linear_problem for the
steady-state in a large lattice of  sites,
with ,  and .  These
parameters are chosen to highlight the key features of different
pumping schemes; we will present numerical results for a more
realistic experimental system in Section .  The
numerical code was written in Julia and is available online, as supplemental material to
Ref. .

The intensity spectrum of the steady-state as a function of pump
frequency is shown in Fig. , where we compare
results for both the Landau and symmetric gauge for four pumping
schemes, discussed in turn below. For simplicity we limit ourselves to
pump frequencies located between the two lowest-lying
Harper-Hofstadter bands of the untrapped system. This allows us to
focus only on states within the first ladder of the trapped system
(Eq. eq:ladders with ). At higher energies, the clear
identification of states is more difficult as more than one ladder of
toroidal Landau levels can overlap, as shown, for example, in
Section .

Single-site pumping-The first and simplest case that we
consider is that of pumping a single site 
 at an off-center lattice site. These results are shown
in Fig.  panel (a), where the uniform energy
spacing of the toroidal Landau levels can be clearly observed. For
this pumping scheme, we find no significant differences between the
spectra for the Landau or symmetric gauge. This is to be expected as
changing the gauge is equivalent to changing the relative phase
between different sites, but as we are only pumping one site, this
phase difference is unimportant.

Instead, for both magnetic gauges, we see that the peak height is very
small for low energy states, rising to a maximum as energy increases,
before decreasing once more. This behaviour can be understand by
considering the form of the real space wavefunctions of .
In real space, the eigenstates are rings of finite width which
increase in radius as the energy increases (as can be seen in
Fig. ). Analytically, we can predict how the
ring radius scales with energy by remembering that the term
 in Eq. eq:ladders is the
momentum-space kinetic energy  where
 is the
physical position operator in the lowest band. From this, we deduce
that , as can be confirmed
numerically. Therefore, if one pumps an off-center site, there will
only be a limited range of rings that will have radii that will
overlap with the pump spot and so be excited. Here, we have set the
pump spot to be at position , and from the above scaling, the
toroidal Landau level that best overlaps with this pump will have a
quantum number , which is in good agreement with the
numerical results shown in Fig.  (a).

     
  Real space reconstruction of the states , 6, 15
and 26 of Fig. (a).
  
Gaussian pumping- We now consider a Gaussian pump as the next
logical step up in complexity from a single-site pump. This has the
form 
, and we choose  and for the pump centre to be at
, as for single-site pumping. The results are shown
in Fig.  (b). The main effect is, as
expected, that more states become visible in both the low- (smaller
) and high-energy (larger ) sections of the
spectrum. This is because the pump has a greater spatial width and so
overlaps with a larger range of real-space eigenstates.

However, we can also see that the intensity spectrum now depends on
the underlying magnetic gauge, as the phase of the eigenstates is
important. In particular, more high-energy peaks can be seen for the
Landau gauge than for the symmetric gauge. This can be most easily
understood by noting that in momentum space, the symmetric-gauge
states also have a ring-like structure (see bottom panel of
Fig. ), where the ring radius increases with
. To see this, we note that, in the symmetric gauge, the
real-space wavefunctions have a phase which winds around the ring as
 where  is the polar angle around the
ring. This phase-winding sets the radius of the rings in momentum
space as ; a scaling that can be
confirmed numerically and seen in Fig. , bottom
panel. (The white spot close to the edges of the rings in these
figures is due to destructive interference with the pump.)  As the
Fourier transform of the Gaussian pump is again a Gaussian, it follows
that only a limited range of low-energy symmetric-gauge momentum-space
states will have a good overlap with the pump. The high energy portion
of the spectrum is therefore washed out compared to its Landau gauge
counterpart, where states have higher amplitude close to the centre of
the BZ and so better overlap with the pump.

    
  Momentum space reconstruction of the eigenstates. Top
row: states corresponding to , 2, 4 and 6 in
Fig. (c), using the Landau gauge.
Bottom row: states corresponding to , 1, 9 and 20 in
Fig. (b), using the symmetric gauge.
  
Before continuing, we also note that for sufficiently large values of
 the symmetric-gauge rings in momentum space will increase to
the point where they touch the BZ boundaries. When this occurs,
self-interference patterns appear in the wave function as shown for
example in Fig. . The extra ring-like structures
appearing for  in Fig.  are due to
the close proximity of states pertaining to other ladders with 
. Note that in order to excite such high-energy states, we have
used a pump with a homogeneous amplitude and a random onsite phase, as
will be presented in the fourth pumping case below.

Homogeneous pumping with uniform phase- If we now take the
limit of a very wide Gaussian, we reach a homogeneous pump profile
extended over all lattice sites. The results for this pumping scheme
are shown in panel (c) of Fig. . Now 
, and we see that the intensity spectrum is strongly
magnetic-gauge dependent. In the Landau gauge, firstly, there are
visible peaks for only half of the states. This can be understood by
noting that a homogeneous pump in real space is a  function in
momentum space centered in the middle of the BZ. If we consider the
Landau-gauge eigenstates in the full BZ, as shown in the top row of
Fig. , we see that the states with an even value
of  have an even number of nodes, with a lobe at the BZ
center. Conversely, the states with odd values of  have an odd
number of nodes, including one at the BZ center. (These feature can be
related back to the properties of the Hermite polynomials in the
analytical eigenstates in the MBZ Eq. eq:chi.) Consequently,
only states with even values of  have a good overlap with the
pump, and the intensity spectrum contains half the expected peaks, now
separated by twice the toroidal Landau level energy spacing.

In the symmetric gauge, secondly, we find only one out of every four
states for homogeneous pumping, as can be seen in the inset of
Fig.  (c). This is due to the fact that, on a
square lattice, the angular momentum is conserved modulo 4, respecting
the 4-fold rotational symmetry. The peak intensity gets smaller for
larger  because of the diminishing overlap of the localized
central pump with the increasing momentum-space ring discussed above.

     
  Momentum space reconstruction of the eigenstates in the
full BZ, using the symmetric gauge and homogeneous pumping with a
random on-site phase. Top row: states corresponding to 
.  Bottom row: states corresponding to 
. Parameters are the same as in
Fig. .
  

Homogeneous pumping with a random on-site phase - As the
fourth scheme, we consider a pump with a uniform amplitude over the
lattice but a random site-dependent phase :
. The phases are chosen from a random
uniform distribution, and have values in the interval . The
bottom panel of Fig.  was obtained by
averaging over 100 distinct realizations of these random phases. This
results in a relatively even intensity distribution, for both gauges,
where we can associate a peak to each toroidal Landau level in this
energy window. While such a pumping scheme would therefore be the best
way to excite all the eigenstates and to fully probe the
momentum-space physics, we note that this would also be difficult to
achieve in an experiment.


     
  Momentum space reconstruction of the eigenstates in the
full BZ, using the Landau (top row) and symmetric gauge (middle and
bottom row) and a spatially homogeneous pump with a random on-site
phase. We have considered the state  for different trap
positions . For the top and middle rows, we have (from
left to right): (0,0) (trap in the center), (2,0), (5.5,0) and (11,0),
whereas for the bottom row we chose the positions (0,2), (0,5.5),
(0,11) and (11,11). Parameters are the same as in
Fig. .
  

     
  Top row: Intensity spectrum for a lattice of
     sites, with ,  and
    .  Second row: Profile of the  state
    in real space (left), momentum space (center) and the population
    over bands in the MBZ (right) for the conservative system.
    Third row: Reconstruction of the  mode
    wavefunction in real space (left), momentum space (center) and the
    population over bands in the MBZ (right) obtained in the
    driven-dissipative system by pumping at , on
    resonance with the desired mode (black dotted line in the top
    panel).  The -like pump at (0,5) is visible as a dark
    square in the left panel.  Bottom row: Slice along the
     line in the MBZ (solid black line), 
    analytical prediction of Eq. eq:chi 
    (blue dotted line) and population over bands for the
    nondissipative system (dashed orange line).
  
Before continuing, we give a final example of an interesting
gauge-dependent effect that could be studied experimentally in this
system. Unlike the physics discussed above, this is not directly
related to the pumping but instead to the behaviour of the wave
function under a change in the centre of the harmonic trap 
. As derived in Ref.  , moving the harmonic
trap in space changes the boundary conditions on the wave function in
the MBZ. We note that although this derivation was made explicitly for
the magnetic Landau gauge in the MBZ, numerically we observe here that
this physics is also seen in the full BZ in both gauges.  As shown
in Fig. , a shift in the harmonic trap centre in
one direction shifts the observed momentum-space pattern in the
perpendicular direction. This behaviour can be understood as a
realisation of Laughlin's Gedankenexperiment for the quantum
Hall effect but now in momentum space . As we
observe, the momentum-space wave function returns to itself after the
harmonic trap has been moved  lattice sites for the magnetic Landau
gauge but  lattice sites for the magnetic symmetric gauge,
reflecting the underlying translational symmetry of  in
the two different gauges.


Results for realistic experimental parameters

We now present numerical results for system parameters within current
experimental reach, to demonstrate that the essential characteristics
of toroidal Landau levels could be probed experimentally for the first
time in photonics. We choose a small lattice of only 
sites, with losses of ; this loss rate is in the same
range as those present in the experiment of
Ref. . Such a large loss rate broadens the
peaks in the intensity spectrum, making closely-spaced eigenenergies
harder to resolve. From Eq. eq:ladders, we see that the level
spacing is given by , and so we can
increase the energy spacing by applying a stronger harmonic potential,
chosen here as . Increasing the strength of the
harmonic trap improves our flat-band approximation, but weakens the
single-band approximation. To compensate for this, we consider a
larger value of , for which the larger band-gap 
 reduces band-mixing effects.

As in the experiment of Ref. , we work in the
Landau gauge for the Harper-Hofstadter Hamiltonian, with hopping
phases given by . In order to model the
experimental pumping scheme where light was injected into a single
resonator at the edge of the system via an external integrated
waveguide , we consider a localized pump on a
single site situated on the upper border of the system at
.  The corresponding intensity spectrum is shown in
the 1st row of Fig. . Apart from the expected broadening
due to larger losses, the peaks observed correspond well to the
expected eigenenergies.  As discussed above, single-site pumping
limits the number of visible peaks, as the heights of the peaks at
low-energies are suppressed due to the poor overlap of the real-space
eigenstate with the pump position. However, as this pumping scheme is
closest to that used in experiments, we emphasise that even in this
case, enough peaks can be observed to extract quantitative
measurements of the toroidal Landau level spacing. The orange vertical
(dashed) lines show the first ladder of eigenstates of 
from Eq. eq:model. We also note that here for frequencies
larger than , we also start to see states from the second
ladder  (see Eq. eq:ladders), which are
depicted as green vertical dash-dotted lines. Their proximity to the
first ladder states means they cannot be easily resolved as separate
peaks in the dissipative spectrum.

Setting the pump frequency at the energy indicated by the black dotted
line, we plot the numerical near- and far-field emission in the left
and center panels of the 3rd row of Fig. . This
corresponds to the wave function in real space and in the full BZ,
respectively. By applying the transformation in Eq. eq:trans, we
can also map the wave function in the full BZ to the population over
bands in the MBZ, as shown in the right panel of the 3rd row of
Fig. . For comparison, we plot these quantities for the
corresponding numerical eigenstate of 
Eq. eq:model in the second row of Fig. , for
which there is no pumping and dissipation.

We find very good qualitative agreement between the numerical results
in the MBZ and the analytical toroidal Landau level Eq. eq:chi
with , as expected. We can make a quantitative comparison with
this analytical eigenstate by taking a slice along the dash-dotted
vertical lines () in the right column of rows 2 and 3;
these cuts are shown in the bottom panel of Fig.  along
with a dotted blue curve indicating the analytical eigenstate. As can
be seen, there is excellent agreement between the numerics without
driving and dissipation and the analytical result. We have checked
that reducing  makes this fit even better, pointing towards
band-mixing effects. Introducing pumping and dissipation distorts the
eigenstate, but many characteristic features are still clearly
observable.

It is particularly interesting to note that in the driven-dissipative
steady-state in real space, shown in the left panel of the 3rd row of
Fig. , the photon distribution breaks the rotational
symmetry of the ring eigenstate. While this can be physically
understood as a decaying cyclotron orbit with an inverse lifetime set
by , in terms of eigenmodes the exponential decay (and more
generally the breaking of the rotational symmetry) results from the
interference of several modes which overlap in frequency due to the
relatively large value of .  In the same way that real-space
Landau levels give rise to real-space cyclotron orbits under the
effect of the magnetic field, the observation of momentum-space Landau
levels can provide clear evidence of a cyclotron orbit in momentum
space under the effect of the Berry curvature, whose effect is indeed
that of a momentum-space magnetic field.

Finally, we briefly summarize how one can practically measure the
contribution  from the off-diagonal matrix elements of the
Berry connection (see Eq. eq:shift) from the intensity
spectrum. Starting from an experimental spectrum, one first needs to
select a particular peak and determine its  label by comparing
the MBZ reconstruction with the analytical result. The distance
between two neighbouring peaks gives the level spacing 
. Finally, to separate the shift  from the
Harper-Hofstadter ground state energy  in Eq. eq:ladders,
one can make use of the fact that the former depends on the trap
strength , while the latter does not. Preparing two otherwise
identical samples with different trap strengths and subtracting the
ground state energy will then allow for a direct measurement of the
contribution from the off-diagonal matrix elements of the Berry
connection.

Conclusion

In conclusion, we have shown that the observation of toroidal Landau
levels in momentum space is within experimental reach for
state-of-the-art driven-dissipative photonic systems. Our proposal
combines the recent realisation of the Harper-Hofstatder model in an
array of silicon-based coupled ring resonators in
Ref. , with a harmonic potential, which could
be introduced through a spatial modulation of the resonator size. We
have presented numerical results to show that even for very small
lattices, in the presence of driving and strong dissipation, key
characteristics of the toroidal Landau levels can still be
extracted. This would be a first direct investigation of analogue
magnetic eigenstates in momentum space.

We have also emphasised that the proposed photonics experiment would
be able to highlight a momentum-space analog of the cyclotron motion
as well as to measure the energy shift due to the off-diagonal matrix
elements of the Berry connection, which, as these are inter-band
geometrical properties, are hard to access by other means. We have
also discussed how the spectroscopic measurements presented here are
sensitive to the specific synthetic magnetic gauge implemented in an
experiment.

Finally, an interesting outlook would be to include the effect of
photon-photon interactions in the model, as the degenerate ground
states predicted in  for a weakly-interacting
trapped Harper-Hofstadter model may lead to interesting nonlinear
dynamical features. In the longer run, when the synthetic gauge field
is combined with strong interactions, one can hope to observe the
hallmarks of fractional quantum Hall
physics .








Conclusions
Conclusions


Summary of the manuscript
The main theme of this thesis was the study of well-established
phenomena, namely superfluidity and magnetism, in a novel context of
driven-dissipative photonic systems, which are subject to losses that
need to be compensated by continuous pumping. After introducing the
basic concepts in the setting of ultracold atomic gases in thermal
equilibrium, we looked at how the effects of driving and dissipation
changed this picture. In particular, we investigated the response of
microcavity exciton-polaritons scattering against a stationary defect,
with an eye on superfluid-like effects, and then turned to the
momentum-space magnetism of resonator arrays.

Polaritons
In the case of polaritons, we looked at the single-fluid pump-only
regime, as well as the three coupled fluids of the optical parametric
oscillator regime created by the parametric scattering of the pump to
the signal and idler states. We found that none of the two regimes
displayed frictionless flow in the strict sense of the word. The
pump-only fluid showed a density modulation localized close to the
defect even in the "superfluid regime", resulting in a residual drag
force that is entirely due to the finite polariton lifetime. On the
other hand, we showed that the optical parametric oscillator regime,
in the optical limiter configuration, always displays propagating
density modulations in the pump, signal and idler, and therefore
violates the definition of frictionless flow in a stronger, conceptual
way. Futhermore, we determined two distinct types of threshold-like
behaviour of the drag force as a function of fluid velocity in the
pump-only case, one of which has no direct analog in equilibrium
weakly-interacting atomic gases. In the optical parametric oscillator
case, we singled out three factors which together determine the
amplitude of the density modulations, and explained why, for typical
experimental conditions, the signal state may appear superfluid.

Harper Hofstadter
In the last part of the manuscript, we have discussed how
momentum-space Landau levels could be observed experimentally in
driven-dissipative photonic systems, by breaking time-reversal
invariance by means of a synthetic gauge field and creating
topologically nontrivial energy bands. We have presented a realistic
proposal for engineering these states by combining a harmonic trap
with an artificial magnetic field for photons in a two-dimensional
ring resonator array. An observation of momentum-space Landau levels
would be the first realisation of magnetism on a toroidal surface. We
have also demonstrated that the main properties of momentum-space
Landau levels may be observed spectroscopically for systems with
realistic experimental parameters. Furthermore, spectroscopic
measurements will be able to access the absolute energy of
eigenstates, measuring novel geometrical features of energy bands,
such as a contribution from the inter-band Berry connection. The
pumping and losses present in these experiments also open the way
towards studying analogue cyclotron orbits in momentum space.


Outlook
In the future, more in-depth studies of polariton superfluidity,
especially regarding the optical parametric oscillator regime, are
needed. In experiments, one could look at multiply-connected
geometries, and on the theoretical side, use nonequilibrium Keldysh
field theory, which was already successfully employed in the case of
incoherently pumped systems, in order to rigorously compute the
superfluid and normal fractions. As far as the resonator arrays are
concerned, it would be interesting to include photon-photon
interaction in future studies, paving the way to fractional quantum
Hall physics in a driven-dissipative setting.