# Singularities in nearly-uniform 1D condensates due to quantum diffusion

C. L. Baldwin  Joint Quantum Institute, NIST/University of Maryland, College
Park, Maryland 20742, USA  P. Bienias  Joint Quantum Institute,
NIST/University of Maryland, College Park, Maryland 20742, USA  Joint Center
for Quantum Information and Computer Science, NIST/University of Maryland,
College Park, Maryland 20742, USA.  A. V. Gorshkov  Joint Quantum Institute,
NIST/University of Maryland, College Park, Maryland 20742, USA  Joint Center
for Quantum Information and Computer Science, NIST/University of Maryland,
College Park, Maryland 20742, USA.  M. J. Gullans  Joint Center for Quantum
Information and Computer Science, NIST/University of Maryland, College Park,
Maryland 20742, USA.  M. Maghrebi  Department of Physics and Astronomy,
Michigan State University, East Lansing, Michigan 48824, USA

###### Abstract

Dissipative systems often exhibit wavelength-dependent loss rates. One
prominent example is Rydberg polaritons formed by electromagnetically-induced
transparency, which have long been a leading candidate for studying the
physics of interacting photons and also hold promise as a platform for quantum
information. In this system, dissipation is in the form of quantum diffusion,
i.e., proportional to $k^{2}$ ($k$ being the wavevector) and vanishing at long
wavelengths as $k\to 0$. Here, we show that one-dimensional condensates
subject to this type of loss are unstable to long-wavelength density
fluctuations in an unusual manner: after a prolonged period in which the
condensate appears to relax to a uniform state, local depleted regions quickly
form and spread ballistically throughout the system. We connect this behavior
to the leading-order equation for the nearly-uniform condensate—a dispersive
analogue to the Kardar-Parisi-Zhang (KPZ) equation—which develops
singularities in finite time. Furthermore, we show that the wavefronts of the
depleted regions are described by purely dissipative solitons within a pair of
hydrodynamic equations, with no counterpart in lossless condensates. We close
by discussing conditions under which such singularities and the resulting
solitons can be physically realized.

Dissipative systems are typically described by a constant dissipation rate,
yet many physical platforms are instead subject to momentum-dependent losses.
A prominent example is Rydberg systems, which have received much interest as a
platform for quantum nonlinear optics [1, 2, 3] and quantum information
processing/simulation [4, 5, 6, 7, 8, 9, 10, 11, 12]. The polaritons that form
under the condition of electromagnetically-induced transparency (EIT) [13, 14,
15] undergo quantum diffusion, i.e., a one-body loss rate $\Gamma_{k}\propto
k^{2}$ [2]. A similar form of dissipation occurs in bosonic atoms driven by
two coherent laser beams [14]. This type of loss can be realized in arrays of
microwave resonators as well by coupling the cavity modes to qubits [16, 17].

In a many-body system, momentum-dependent loss can have drastic consequences,
from dissipatively stabilizing condensates [18] to producing exotic critical
or correlated states [19, 20, 17, 21]. These advances notwithstanding, many
consequences of momentum-dependent loss remain undiscovered.

Figure 1: (Top) Time evolution of a small Gaussian density perturbation.
(Bottom) Snapshots of the density profile (normalized by the initial density
$\rho_{0}$) at the times indicated by the dashed white lines in the top panel.
Simulation parameters are the dissipation strength $\lambda=2.0$, height of
the initial Gaussian $h=0.1$, its width $w=15\xi$, the spatial discretization
$\Delta x=0.2\xi$, and the time step $\Delta t=0.1\tau$. $\xi$ and $\tau$
define the coherence length/time scales.

In this paper, we investigate a driven-dissipative condensate in 1D subject to
one-body loss $\Gamma_{k}\sim\lambda k^{2}$. We show that when perturbed from
uniformity, this system exhibits a striking instability, best demonstrated by
the example in Fig. 1. Shown is the density profile of a condensate as a
function of time, obtained by numerical simulation of the Gross-Pitaevskii
equation (details to be explained below). The condensate initially has a
slight localized excess of particles. The excess density begins to spread
throughout the system, and the condensate appears to relax to a uniform state.
However, after a significant delay, the density quickly drops to zero in
certain regions, forming fronts which move ballistically and eventually
consume the entire condensate.

We show that the onset of instability can be attributed to the long-wavelength
equation for the phase of the nearly-uniform condensate. Whereas driven-
dissipative condensates with $k$-independent loss are typically described by
the Kardar-Parisi-Zhang (KPZ) equation [22, 23, 24, 25, 26], a well-known
nonlinear diffusion equation, here we find an analogous nonlinear wave
equation which we refer to as “dispersive KPZ”. Little is known about the
dispersive KPZ equation, at least in the physics literature, but a surprising
feature of the latter is that generic solutions diverge in finite time [27].
We show that this singularity corresponds to the sudden depletion of the
condensate.

The dynamics following formation of the depleted regions can no longer be
described by dispersive KPZ, for which solutions simply do not exist beyond
the singularity time. We thus derive a more general pair of hydrodynamic
equations, and identify soliton solutions which accurately describe the shape
and motion of the fronts seen in Fig. 1. As will be clear, these solitons are
exclusive to dissipative condensates, and in fact, their core size diverges in
the limit of vanishing dissipation, $\lambda\rightarrow 0$.

Dissipative Gross-Pitaevskii equation.—We consider a one-dimensional gas of
particles with contact interactions and single-body loss
$\Gamma_{k}\sim\lambda k^{2}$. Formally, the system is described by the
quantum master equation ($\hbar=1$)

$\displaystyle\partial_{t}\rho$
$\displaystyle=-i\big{(}\hat{H}_{\textrm{eff}}\rho-\rho\hat{H}_{\textrm{eff}}^{{\dagger}}\big{)}+\int\textrm{d}x\,\frac{\lambda}{m}(\partial_{x}\hat{\psi})\rho(\partial_{x}\hat{\psi}^{{\dagger}}),$
(1)
$\displaystyle\hat{H}_{\textrm{eff}}=\int\textrm{d}x\left[\frac{1-i\lambda}{2m}(\partial_{x}\hat{\psi}^{{\dagger}})(\partial_{x}\hat{\psi})+U\hat{\psi}^{{\dagger}2}\hat{\psi}^{2}\right],$
(2)

where $\rho$ is the density matrix of the system and
$\hat{\psi}^{{\dagger}}(x)$ creates a bosonic particle at position $x$. Here
$m$ is the mass and $U>0$ governs the strength of interactions.

Following the standard procedure, e.g., as in Refs. [28, 29], we first derive
the semiclassical equation of motion for the condensate wavefunction
$\psi(x)$, valid at large densities $\rho_{0}\gg mU$:

$i\partial_{t}\psi+\frac{1-i\lambda}{2m}\partial_{x}^{2}\psi-2U\big{|}\psi\big{|}^{2}\psi=0.$
(3)

Equation (3) is quite similar to the standard Gross-Pitaevskii (GP) equation,
with the only difference being that the coefficient of the kinetic term is
complex. Therefore, any spatial variation of the wavefunction leads to
dissipation. We shall focus on the dynamics of a nearly-uniform condensate.
For concreteness, we use initial conditions of the form

$\psi(x,0)=\sqrt{\rho_{0}}\left(1+he^{-\frac{x^{2}}{w^{2}}}\right)^{\frac{1}{2}}.$
(4)

We have confirmed that the conclusions of this paper hold for other initial
conditions as well (sinusoidal perturbations, random density/phase
fluctuations, etc.).

The natural length scale of Eq. (3) is the healing length
$\xi\equiv\sqrt{1/mU\rho_{0}}$, and the natural time scale is $\tau\equiv
m\xi^{2}$. The remaining dimensionless parameters are the dissipation strength
$\lambda$, the magnitude of the density perturbation $h$, and the width of the
density perturbation $w/\xi$.

Figure 2: Demonstration of the scaling form for the singularity time
$\tau_{\textrm{sing}}$. The solid line is a power-law $z^{-2}$, drawn for
comparison. $\tau_{\textrm{sing}}$ is plotted for Gaussian density
perturbations of various heights $h$ and widths $w$ (see the inset), and
various dissipation strengths $\lambda\in\\{0.05,0.1,0.15,0.2,0.25,0.3\\}$.
System size is $L=10000\xi$.

Fig. 1, showcased earlier, displays a representative simulation of Eq. (3)
using the initial profile in Eq. (4). The behavior is highly non-trivial—a
prolonged period during which the condensate is nearly uniform is followed by
the sudden appearance and subsequent spread of fully depleted regions. We
refer to the sudden depletion as a “singularity”. While the density profile is
strictly analytic as a function of time, the long-wavelength equation derived
below exhibits a genuine singularity which acts as a precursor to the
condensate depletion.

For concreteness, let us define $\tau_{\textrm{sing}}$ as the time when
$\rho(x,t)\equiv|\psi(x,t)|^{2}$ first drops below $\rho_{0}/2$ at some
position $x$, i.e., the first time at which
$\textrm{min}_{x}\rho(x,t)<\rho_{0}/2$. Figure 2 plots $\tau_{\textrm{sing}}$
for multiple choices of $\lambda$ and initial conditions. A clear scaling form
is seen:

$\frac{\tau_{\textrm{sing}}(\lambda,w,h)}{\tau}\sim\frac{w}{\xi}{\mathcal{F}}\left(\frac{\lambda
wh}{\xi}\right),$ (5)

where the scaling function appears to fall off as ${\mathcal{F}}(z)\sim
z^{-2}$ for $z\lesssim 1$. Such algebraic dependence implies that the
underlying instability is fundamentally different from nucleation, where a
metastable state tunnels into a true equilibrium state, for which the decay
rate would be exponentially suppressed at small fluctuations/perturbations.
The instability reported here is governed by a different mechanism that
follows from the long-wavelength description of the condensate.

Dispersive KPZ equation.—To derive the long-wavelength effective equation for
the nearly-uniform condensate, starting from Eq. (3), we: i) write
$\psi(x,t)=\sqrt{\rho_{0}+\Delta\rho(x,t)}e^{i\theta(x,t)}$, assuming
$\Delta\rho\ll\rho_{0}$; and ii) retain only the terms in the GP equation
which are both lowest-order in $\Delta\rho/\rho_{0}$ and most relevant at long
wavelengths. The calculation is given in the Supplemental Material (SM) [30].
The end result is

$\frac{1}{c^{2}}\partial_{t}^{2}\theta=\partial_{x}^{2}\theta+\lambda(\partial_{x}\theta)^{2},$
(6)

with $c\equiv\sqrt{2}\xi/\tau$ being a velocity scale that characterizes the
“speed of sound”. The density variation in turn comes out to be
$\Delta\rho=-({\rho_{0}\tau}/{2})\partial_{t}\theta$.

Figure 3: Comparison of solution to GP equation (red lines) against solutions
to approximate equations for time evolution of a Gaussian density
perturbation. Times are indicated in each panel. (Upper three) Comparison to
dispersive KPZ, shown in blue. (Bottom) Comparison to soliton given by Eqs.
(13) and (14), shown in green, with the constant $z_{0}$ chosen to match the
center of the front. Simulation details: dissipation strength is
$\lambda=0.4$, height of the initial Gaussian is $h=0.05$, width is
$w=200\xi$, spatial step size is $\Delta x=0.2\xi$, temporal step size is
$\Delta t=0.1\tau$.

Equation (6) is quite similar to the (noiseless) KPZ equation which has
emerged in generic dissipative condensates [22, 23, 24, 25], except for the
second time derivative on the left-hand side, which results in a wave-like
equation with $c$ defining a causal “light cone” [27]. Being a nonlinear wave
equation, we refer to Eq. (6) as the “dispersive KPZ” equation. Much less is
known about dispersive KPZ than its diffusive counterpart [31, 32, 33, 34],
but one established result is that under certain conditions, solutions to the
dispersive KPZ equation—as well as a larger class of nonlinear hyperbolic
equations—diverge in finite time [27]. On physical grounds, this is due to the
absence of any damping term such as $\partial_{t}\theta$ which could
counteract the nonlinear growth. We have confirmed this divergence through
numerical simulation of Eq. (6).

Figure 3 compares the solution of dispersive KPZ to the solution of the GP
equation for a representative example. We see that: i) the two agree extremely
well for as long as $\Delta\rho/\rho_{0}$ is everywhere small; and ii)
development of the singularity in dispersive KPZ coincides with the depletion
of the condensate. For this reason, we equate the singularity time with
$\tau_{\textrm{sing}}$ 111Unlike the singularity time in dispersive KPZ, which
is a well-defined instant, our definition of $\tau_{\textrm{sing}}$ in the GP
equation is somewhat arbitrary due to setting a threshold at $\rho_{0}/2$ (see
the discussion above Eq. (5)). But in the limit where
$\tau_{\textrm{sing}}\gg\tau$, using any other fraction instead of $1/2$ would
change $\tau_{\textrm{sing}}$ only by a subleading amount..

The scaling form of $\tau_{\textrm{sing}}$ given in Eq. (5) then follows from
the scaling of solutions to dispersive KPZ. Suppose that, just as in the
simulations above, initially $\theta(x,0)=0$ and
$\Delta\rho(x,0)=-\frac{\rho_{0}\tau}{2}\partial_{t}\theta(x,0)$ is of the
form

$\Delta\rho(x,0)=\frac{\rho_{0}h}{2}{\mathcal{G}}\left(\frac{x}{w}\right),$
(7)

for some dimensionless function ${\mathcal{G}}(y)$. Defining $y\equiv x/w$,
$s\equiv tc/w$, $\phi(y,s)\equiv\lambda\theta(x,t)$, the dispersive KPZ
equation together with the initial conditions can be written as

$\begin{gathered}\partial_{s}^{2}\phi=\partial_{y}^{2}\phi+(\partial_{y}\phi)^{2},\\\
\phi(y,0)=0,\quad\partial_{s}\phi(y,0)=-\frac{\lambda
wh}{\sqrt{2}\xi}{\mathcal{G}}(y).\end{gathered}$ (8)

The only dimensionless parameter here is $\lambda wh/\xi$, hence the scaling
form in Eq. (5). This gives further evidence for the applicability of
dispersive KPZ 222For more complicated initial conditions,
$\tau_{\textrm{sing}}$ will not obey so simple a scaling form..

Unfortunately, the dispersive KPZ equation does not have a general analytic
solution (although a solvable special case is given in the SM [30]). Thus let
us briefly discuss an analogous but simpler equation that exhibits similar
features:

$\partial_{t}\tilde{\theta}=\partial_{x}\tilde{\theta}+\lambda\tilde{\theta}^{2},$
(9)

which, in dimensionless coordinates, describes a left-moving wave with an
additional nonlinear term (with $\tilde{\theta}$ roughly mimicking
$\partial_{x}\theta$ 333To see the analogy with dispersive KPZ, define (again
setting $\tau=\xi=1$ for simplicity)
$\theta_{1}\equiv\partial_{t}\theta,\theta_{2}\equiv\partial_{x}\theta$.
Dispersive KPZ can be written as the pair of first-order equations
$\partial_{t}\theta_{1}=\partial_{x}\theta_{2}+\lambda\theta_{2}^{2}$,
$\partial_{t}\theta_{2}=\partial_{x}\theta_{1}$.). It is trivial to solve this
equation by transforming to the frame moving alongside the wave: along the
path $x(t)=x_{0}-t$, Eq. (9) simply becomes
$\textrm{d}\tilde{\theta}/\textrm{d}t=\lambda\tilde{\theta}^{2}$. Thus the
general solution is

$\tilde{\theta}(x_{0}-t,t)=\frac{\tilde{\theta}_{0}}{1-\lambda\tilde{\theta}_{0}t},$
(10)

where $\tilde{\theta}_{0}\equiv\tilde{\theta}(x_{0},0)$. We see that, unless
$\tilde{\theta}(x_{0},0)$ is everywhere negative, $\tilde{\theta}(x,t)$ will
diverge in finite time, regardless of the precise shape of the initial
condition. The same phenomenon occurs in the setting of the dispersive KPZ
equation. Note that this behavior is much more drastic than a linear
instability, where the amplitude would grow exponentially but nonetheless be
finite at any finite time.

Hydrodynamic equations.—For times greater than $\tau_{\textrm{sing}}$, the
dispersive KPZ equation clearly cannot describe the evolution of the
condensate. Thus we derive a pair of hydrodynamic equations which no longer
assume $\Delta\rho\ll\rho_{0}$, only requiring that the relevant length and
time scales still be larger than $\xi$ and $\tau$, respectively. We follow the
standard procedure for quantum fluids by describing the wavefunction in terms
of the density $\rho(x,t)$ and velocity field
$v(x,t)\equiv\partial_{x}\theta(x,t)/m$ [38, 39]. The resulting hydrodynamic
equations become (see the SM for details [30])

$\displaystyle\partial_{t}\rho+\partial_{x}\big{(}\rho v\big{)}=-\lambda m\rho
v^{2},$ (11)
$\displaystyle\partial_{t}v+v\partial_{x}v=-2mU^{2}\partial_{x}\rho.$ (12)

Equation (11) is the analogue of the continuity equation, with the additional
feature that the density is depleted in regions of nonzero velocity. Equation
(12) is the standard Euler equation for an incompressible fluid, with the
pressure given by $P(\rho)=mU^{2}\rho^{2}$ [hence the right-hand side can be
written as $-\rho^{-1}\partial_{x}P(\rho)$] [39].

One can confirm by direct substitution that the above equations admit soliton
solutions—$\rho(x,t)=\rho(x-ut)$, $v(x,t)=v(x-ut)$—for any velocity $u$ such
that $|u|\geq c$ ($c$ being defined as before). Supersonic solitons are likely
unstable, therefore we focus on the case $u=c$, where the soliton moves
rightward at the speed of sound. In terms of $z\equiv x-ct$, we obtain [30]

$\frac{\rho(z)}{\rho_{0}}=1+\frac{v(z)}{c}-\frac{v(z)^{2}}{2c^{2}},\;\;\;\;\frac{v(z)}{c}=f^{-1}\big{[}\frac{\sqrt{2}\lambda}{3\xi}(z_{0}-z)\big{]},$
(13)

where $z_{0}$ is a constant which fixes the center of the soliton and $f^{-1}$
is the inverse of the function

$\displaystyle f(y)=\log{(-y)}$
$\displaystyle-\frac{2\sqrt{3}+3}{6}\log{(\sqrt{3}-1+y)}$ (14)
$\displaystyle+\frac{2\sqrt{3}-3}{6}\log{(\sqrt{3}+1-y)}.$

Note that the density $\rho$ approaches $\rho_{0}$ as $z\to\infty$, while it
vanishes as $z\to-\infty$. The fronts observed in our simulations of the GP
equation agree well with Eq. (13) (the left-moving fronts are easily related
to the above by symmetry). A representative comparison is shown in the bottom
panel of Fig. 3.

These solitons are quite different from those in the dissipation-free GP
equation [$\rho\sim\rho_{0}\tanh^{2}{(z/\sqrt{2}\xi)}$] [38]. Most
importantly, the dissipative solitons have a core size $\xi/\lambda$ (as
opposed to simply $\xi$), which diverges in the limit of vanishing
dissipation, $\lambda\rightarrow 0$. This is consistent with the fact that
these solitons originate from an instability which occurs only in the presence
of quantum diffusion.

Physical realizations.—Let us briefly comment on potential physical
realizations of this phenomenon. As noted above, one possible platform is
Rydberg polaritons via electromagnetically-induced transparency (EIT), formed
when an incoming photon hybridizes with a long-lived Rydberg state through a
lossy intermediate state [40, 41, 4]. At precisely zero momentum, the
polariton is a superposition of Rydberg state and photon with exactly zero
amplitude on the lossy state, and hence is essentially lossless. The deviation
from resonance at small but finite $k$ leads to the $k^{2}$ loss and results
in a diffusion-like term [2]. Furthermore, at low energies, we can neglect
scattering into other modes, leaving Eq. (36) as the effective many-body
Hamiltonian.

While the interaction between polaritons is generically complex-valued as
well, we have confirmed that it is possible to tune microscopic parameters so
that the effective two-body loss rate vanishes while the one-body ($k^{2}$)
loss remains significant; see the SM [30]. Thus the instability reported here
may be observable in Rydberg polariton systems, although the parameter regime
in which $\tau_{\textrm{sing}}\gg\tau$ (where the “singularity” is sharpest)
would necessitate a long atomic medium. Running-wave cavities may provide a
feasible alternative to the long free-space lengths.

An alternate realization could come from a 1D cloud of bosonic atoms driven by
two coherent lasers under EIT condition. With one beam orthogonal to the
atomic gas and the other parallel, detuning (proportional to the atomic wave
vector $k$) due to the Doppler shift leads to diffusion-like dynamics [14]. In
order to ensure that the contact interaction does not itself cause losses, one
would have to properly choose the states involved and tune interactions, e.g.,
with a magnetic field [42]. Finally, microcavity arrays [18, 43] provide
another platform where $k^{2}$ loss can be realized [16, 17]. However, it may
be challenging to engineer coherent interactions and diffusive terms
simultaneously.

Conclusion.—We have shown that 1D driven-dissipative condensates for which
quantum diffusion is the dominant source of dissipation suffer from a peculiar
instability to local density perturbations. The condensate relaxes towards
uniform density until a time $\tau_{\textrm{sing}}$—much larger than the
natural timescale $\tau$—after which certain regions quickly deplete and form
fronts which then spread throughout the condensate. We have traced this
behavior to the long-wavelength effective equation for the phase of the
condensate, a nonlinear wave equation which we refer to as the “dispersive
KPZ” equation. Solutions to dispersive KPZ can diverge at finite times, and we
have observed that the singularity in the long-wavelength description
coincides with depletion of the condensate. We have further derived a pair of
hydrodynamic equations for the condensate that accurately describe the
dynamics even beyond the onset of instability. Interestingly, the fronts are
described by non-standard soliton solutions that emerge solely due to
dissipation.

From a mathematical perspective, it has long been known that the solutions to
nonlinear wave equations can diverge, or more generally become nonanalytic
[44, 45]. It is interesting to note that whereas the divergence is often seen
as an unphysical mathematical pathology, here it corresponds to a genuine
physical phenomenon. Coincidentally, Ref. [27] even comments: “there is, to
our knowledge, no direct application of [dispersive KPZ] to a physical
problem”. The situation discussed here—condensates undergoing quantum
diffusion—provides such an application, the first such to our knowledge.

Many directions for future work remain. First of all, it is desirable to go
beyond the semiclassical limit and investigate the strongly-interacting
quantum regime. A step in this direction would be to include noise terms in
Eqs. (3) and (6) [28, 29]. Most studies of the traditional KPZ equation do
include a noise term, as it is the competition between noise and nonlinearity
which leads to novel scaling properties [46, 47], and so it is natural to ask
whether dispersive KPZ has its own distinct scaling behavior. Furthermore, the
traditional solitons and hydrodynamic behavior of condensates have been well-
studied [38, 48, 39], whereas we have only scratched the surface of the
present equations. Finally, further scrutiny of different physical
realizations is worthwhile. While we have intentionally kept our analysis
theoretical and abstract, more systematic investigations are needed to assess
the feasibility of any specific implementation.

Acknowledgments.—The authors would like to thank J. V. Porto for informative
discussions. This research was performed while C.L.B. held an NRC Research
Associateship award at the National Institute of Standards and Technology.
P.B. and A.V.G. acknowledge funding by the AFOSR, AFOSR MURI, DoE ASCR Quantum
Testbed Pathfinder program (award No. DE-SC0019040), U.S. Department of Energy
Award No. DE-SC0019449, DoE ASCR Accelerated Research in Quantum Computing
program (award No. DE-SC0020312), NSF PFCQC program, and ARO MURI. M.M.
acknowledges support from NSF under Grant No. DMR-1912799, the Air Force
Office of Scientific Research (AFOSR) under award number FA9550-20-1-0073 as
well as the start-up funding from Michigan State University.

## References

  * Gorshkov  _et al._ [2011] A. V. Gorshkov, J. Otterbach, M. Fleischhauer, T. Pohl, and M. D. Lukin, Photon-photon interactions via Rydberg blockade, Phys. Rev. Lett. 107, 133602 (2011).
  * Peyronel  _et al._ [2012] T. Peyronel, O. Firstenberg, Q.-Y. Liang, S. Hofferberth, A. V. Gorshkov, T. Pohl, M. D. Lukin, and V. Vuletić, Quantum nonlinear optics with single photons enabled by strongly interacting atoms, Nature 488, 57 (2012).
  * Dudin and Kuzmich [2012] Y. O. Dudin and A. Kuzmich, Strongly interacting Rydberg excitations of a cold atomic gas, Science 336, 887 (2012).
  * Friedler  _et al._ [2005] I. Friedler, D. Petrosyan, M. Fleischhauer, and G. Kurizki, Long-range interactions and entanglement of slow single-photon pulses, Phys. Rev. A 72, 043803 (2005).
  * Saffman  _et al._ [2010] M. Saffman, T. G. Walker, and K. Mølmer, Quantum information with Rydberg atoms, Rev. Mod. Phys. 82, 2313 (2010).
  * Weimer  _et al._ [2010] H. Weimer, M. Müller, I. Lesanovsky, P. Zoller, and H. P. Büchler, A Rydberg quantum simulator, Nat. Phys. 6, 382 (2010).
  * Li  _et al._ [2013] L. Li, Y. O. Dudin, and A. Kuzmich, Entanglement between light and an optical atomic excitation, Nature 498, 466 (2013).
  * Maxwell  _et al._ [2013] D. Maxwell, D. J. Szwer, D. Paredes-Barato, H. Busche, J. D. Pritchard, A. Gauguet, K. J. Weatherill, M. P. A. Jones, and C. S. Adams, Storage and control of optical photons using Rydberg polaritons, Phys. Rev. Lett. 110, 103001 (2013).
  * Otterbach  _et al._ [2013] J. Otterbach, M. Moos, D. Muth, and M. Fleischhauer, Wigner crystallization of single photons in cold Rydberg ensembles, Phys. Rev. Lett. 111, 113001 (2013).
  * Gorniaczyk  _et al._ [2016] H. Gorniaczyk, C. Tresp, P. Bienias, A. Paris-Mandoki, W. Li, I. Mirgorodskiy, H. P. Büchler, I. Lesanovsky, and S. Hofferberth, Enhancement of Rydberg-mediated single-photon nonlinearities by electrically tuned Förster resonances, Nat. Commun. 7, 12480 (2016).
  * Tiarks  _et al._ [2016] D. Tiarks, S. Schmidt, G. Rempe, and S. Dürr, Optical $\pi$ phase shift created with a single-photon pulse, Sci. Adv. 2, e1600036 (2016).
  * Bienias and Büchler [2020] P. Bienias and H. P. Büchler, Two photon conditional phase gate based on Rydberg slow light polaritons, J. Phys. B 53, 054003 (2020).
  * Harris  _et al._ [1990] S. E. Harris, J. E. Field, and A. Imamoğlu, Nonlinear optical processes using electromagnetically induced transparency, Phys. Rev. Lett. 64, 1107 (1990).
  * Fleischhauer  _et al._ [2005] M. Fleischhauer, A. Imamoglu, and J. P. Marangos, Electromagnetically induced transparency: Optics in coherent media, Rev. Mod. Phys. 77, 633 (2005).
  * Mohapatra  _et al._ [2007] A. K. Mohapatra, T. R. Jackson, and C. S. Adams, Coherent optical detection of highly excited Rydberg states using electromagnetically induced transparency, Phys. Rev. Lett. 98, 113003 (2007).
  * Marcos  _et al._ [2012] D. Marcos, A. Tomadin, S. Diehl, and P. Rabl, Photon condensation in circuit quantum electrodynamics by engineered dissipation, New J. Phys. 14, 055005 (2012).
  * Marino and Diehl [2016] J. Marino and S. Diehl, Driven Markovian quantum criticality, Phys. Rev. Lett. 116, 070407 (2016).
  * Diehl  _et al._ [2008] S. Diehl, A. Micheli, A. Kantian, B. Kraus, H. P. Büchler, and P. Zoller, Quantum states and phases in driven open quantum systems with cold atoms, Nat. Phys. 4, 878 (2008).
  * Eisert and Prosen [2010] J. Eisert and T. Prosen, Noise-driven quantum criticality, arXiv:1012.5013 (2010).
  * Höning  _et al._ [2012] M. Höning, M. Moos, and M. Fleischhauer, Critical exponents of steady-state phase transitions in fermionic lattice models, Phys. Rev. A 86, 013606 (2012).
  * Zeuthen  _et al._ [2017] E. Zeuthen, M. J. Gullans, M. F. Maghrebi, and A. V. Gorshkov, Correlated photon dynamics in dissipative Rydberg media, Phys. Rev. Lett. 119, 043602 (2017).
  * Kuramoto [1984] Y. Kuramoto,  _Chemical Oscillations, Waves, and Turbulence_ (Springer Series in Synergetics, 1984).
  * Grinstein  _et al._ [1993] G. Grinstein, D. Mukamel, R. Seidin, and C. H. Bennett, Temporally periodic phases and kinetic roughening, Phys. Rev. Lett. 70, 3607 (1993).
  * Grinstein  _et al._ [1996] G. Grinstein, C. Jayaprakash, and R. Pandit, Conjectures about phase turbulence in the complex Ginzburg-Landau equation, Physica D 90, 96 (1996).
  * Altman  _et al._ [2015] E. Altman, L. M. Sieberer, L. Chen, S. Diehl, and J. Toner, Two-dimensional superfluidity of exciton polaritons requires strong anisotropy, Phys. Rev. X 5, 011017 (2015).
  * Maghrebi [2017] M. F. Maghrebi, Fragile fate of driven-dissipative XY phase in two dimensions, Phys. Rev. B 96, 174304 (2017).
  * Escudero [2007] C. Escudero, Blow-up of the hyperbolic Burgers equation, J. Stat. Phys. 127, 327 (2007).
  * Kamenev [2011] A. Kamenev,  _Field Theory of Non-Equilibrium Systems_ (CUP, 2011).
  * Sieberer  _et al._ [2016] L. M. Sieberer, M. Buchhold, and S. Diehl, Keldysh field theory for driven open quantum systems, Rep. Prog. Phys. 79, 096001 (2016).
  * [30] See the Supplementary Material at [url] for the following: S1) Derivation of dispersive KPZ from dissipative GP equation; S2) Example of singularity in a solution to dispersive KPZ; S3) Derivation of hydrodynamic equations and soliton solutions; S4) Calculation of one-body and two-body loss rates for Rydberg polaritons. See also Refs. [27, 39, 49, 38, 50, 51, 52].
  * Souplet [1995] P. Souplet, Nonexistence of global solutions to some differential inequalities of the second order and applications, Port. Math. 52, 289 (1995).
  * Makarenko  _et al._ [1997] A. S. Makarenko, M. N. Moskalkov, and S. P. Levkov, On blow-up solutions in turbulence, Phys. Lett. A 235, 391 (1997).
  * Liu and Natalini [2001] H. Liu and R. Natalini, Longtime diffusive behavior of solutions to a hyperbolic relaxation system, Asymptot. Anal. 25, 21 (2001).
  * Orive and Zuazua [2006] R. Orive and E. Zuazua, Long-time behavior of solutions to a nonlinear hyperbolic relaxation system, J. Differ. Equ. 228, 17 (2006).
  * Note [1] Unlike the singularity time in dispersive KPZ, which is a well-defined instant, our definition of $\tau_{\textrm{sing}}$ in the GP equation is somewhat arbitrary due to setting a threshold at $\rho_{0}/2$ (see the discussion above Eq. (5)). But in the limit where $\tau_{\textrm{sing}}\gg\tau$, using any other fraction instead of $1/2$ would change $\tau_{\textrm{sing}}$ only by a subleading amount.
  * Note [2] For more complicated initial conditions, $\tau_{\textrm{sing}}$ will not obey so simple a scaling form.
  * Note [3] To see the analogy with dispersive KPZ, define (again setting $\tau=\xi=1$ for simplicity) $\theta_{1}\equiv\partial_{t}\theta,\theta_{2}\equiv\partial_{x}\theta$. Dispersive KPZ can be written as the pair of first-order equations $\partial_{t}\theta_{1}=\partial_{x}\theta_{2}+\lambda\theta_{2}^{2}$, $\partial_{t}\theta_{2}=\partial_{x}\theta_{1}$.
  * Pethick and Smith [2002] C. J. Pethick and H. Smith,  _Bose-Einstein Condensation in Dilute Gases_ (CUP, 2002).
  * Cazalilla  _et al._ [2011] M. A. Cazalilla, R. Citro, T. Giamarchi, E. Orignac, and M. Rigol, One dimensional bosons: From condensed matter systems to ultracold gases, Rev. Mod. Phys. 83, 1405 (2011).
  * Fleischhauer and Lukin [2000] M. Fleischhauer and M. D. Lukin, Dark-state polaritons in electromagnetically induced transparency, Phys. Rev. Lett. 84, 5094 (2000).
  * Lukin  _et al._ [2001] M. D. Lukin, M. Fleischhauer, R. Cote, L. M. Duan, D. Jaksch, J. I. Cirac, and P. Zoller, Dipole blockade and quantum information processing in mesoscopic atomic ensembles, Phys. Rev. Lett. 87, 037901 (2001).
  * Bloch  _et al._ [2008] I. Bloch, J. Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Rev. Mod. Phys. 80, 885 (2008).
  * Houck  _et al._ [2012] A. A. Houck, H. E. Türeci, and J. Koch, On-chip quantum simulation with superconducting circuits, Nat. Phys. 8, 292 (2012).
  * John [1990] F. John,  _Nonlinear Wave Equations_ (AMS, 1990).
  * Alinhac [1995] S. Alinhac,  _Blowup for Nonlinear Hyperbolic Equations_ (Birkhauser, 1995).
  * Kardar  _et al._ [1986] M. Kardar, G. Parisi, and Y.-C. Zhang, Dynamic scaling of growing interfaces, Phys. Rev. Lett. 56, 889 (1986).
  * Halpin-Healy and Zhang [1995] T. Halpin-Healy and Y.-C. Zhang, Kinetic roughening phenomena, stochastic growth, directed polymers and all that. Aspects of multidisciplinary statistical mechanics, Phys. Rep. 254, 215 (1995).
  * Leggett [2006] A. J. Leggett,  _Quantum Liquids_ (OUP, 2006).
  * Wazwaz [2002] A.-M. Wazwaz,  _Partial Differential Equations_ (CRC Press, 2002).
  * Gullans  _et al._ [2016] M. J. Gullans, J. D. Thompson, Y. Wang, Q. Y. Liang, V. Vuletić, M. D. Lukin, and A. V. Gorshkov, Effective field theory for Rydberg polaritons, Phys. Rev. Lett. 117, 113601 (2016).
  * Bienias  _et al._ [2014] P. Bienias, S. Choi, O. Firstenberg, M. F. Maghrebi, M. Gullans, M. D. Lukin, A. V. Gorshkov, and H. P. Büchler, Scattering resonances and bound states for strongly interacting Rydberg polaritons, Phys. Rev. A 90, 053804 (2014).
  * Barlette  _et al._ [2000] V. E. Barlette, M. M. Leite, and S. K. Adhikari, Quantum scattering in one dimension, Eur. J. Phys. 21, 435 (2000).

## Supplemental Material to “Singularities in nearly-uniform 1D condensates
due to quantum diffusion”

The contents of this Supplemental Material are as follows. In Sec. I, we show
how the dispersive KPZ equation can be derived from the dissipative GP
equation using two approximations: nearly-uniform density and long
wavelengths. In Sec. II, we give a special case of a solution to dispersive
KPZ which exhibits the finite-time singularity. In Sec. III, we derive the
pair of hydrodynamic equations from dissipative GP by making only the long-
wavelength approximation, and then identify soliton solutions to them.
Finally, in Sec. IV, we calculate the effective one-body and two-body loss
rates in terms of microscopic parameters for Rydberg polariton systems, and
demonstrate that it is possible to tune the two-body loss rate to zero.

##  I Derivation of dispersive KPZ from dissipative Gross-Pitaevskii

Here we derive the dispersive KPZ equation, Eq. (6) in the main text. Starting
from the dissipative Gross-Pitaevskii (GP) equation, reproduced here:

$i\partial_{t}\psi+\frac{1-i\lambda}{2m}\partial_{x}^{2}\psi-2U\big{|}\psi\big{|}^{2}\psi=0,$
(15)

let us first convert to a pair of equations for the density $\rho(x,t)$ and
phase $\theta(x,t)$ via $\psi(x,t)\equiv\sqrt{\rho(x,t)}e^{i\theta(x,t)}$:

$\partial_{t}\rho+\frac{1}{m}\partial_{x}\big{(}\rho\partial_{x}\theta\big{)}=\frac{\lambda}{m}\big{(}\sqrt{\rho}\partial_{x}^{2}\sqrt{\rho}-\rho(\partial_{x}\theta)^{2}\big{)},$
(16)
$\partial_{t}\theta+\frac{1}{2m}\big{(}\partial_{x}\theta\big{)}^{2}-\frac{1}{2m\sqrt{\rho}}\partial_{x}^{2}\sqrt{\rho}+2U\rho=\frac{\lambda}{2m\rho}\partial_{x}\big{(}\rho\partial_{x}\theta\big{)}.$
(17)

We next make two approximations:

  * •

The density can be written $\rho(x,t)=\rho_{0}+\Delta\rho(x,t)$ with
$\Delta\rho(x,t)\ll\rho_{0}$, and terms sub-leading in $\Delta\rho/\rho_{0}$
can be neglected. Note that merely by writing $\rho=\rho_{0}+\Delta\rho$, Eq.
(17) acquires a constant term $2U\rho_{0}$ which we can trivially remove by
taking $\theta\rightarrow\theta-2U\rho_{0}t$.

  * •

All terms which are sub-leading at long wavelengths and low frequencies (in
the sense to be defined momentarily) can be neglected.

The first approximation leads to the following (after defining
$\xi\equiv\sqrt{1/mU\rho_{0}}$, $\tau\equiv m\xi^{2}$ as in the main text):

$\partial_{t}\Delta\rho+\frac{\rho_{0}\xi^{2}}{\tau}\partial_{x}^{2}\theta+\frac{\xi^{2}}{\tau}\big{(}\partial_{x}\Delta\rho\big{)}\big{(}\partial_{x}\theta\big{)}=\frac{\lambda\xi^{2}}{2\tau}\partial_{x}^{2}\Delta\rho-\frac{\lambda\rho_{0}\xi^{2}}{\tau}\big{(}\partial_{x}\theta\big{)}^{2},$
(18)
$\partial_{t}\theta+\frac{\xi^{2}}{2\tau}\big{(}\partial_{x}\theta\big{)}^{2}-\frac{\xi^{2}}{4\rho_{0}\tau}\partial_{x}^{2}\Delta\rho+\frac{2}{\rho_{0}\tau}\Delta\rho=\frac{\lambda\xi^{2}}{2\tau}\partial_{x}^{2}\theta+\frac{\lambda\xi^{2}}{2\rho_{0}\tau}\big{(}\partial_{x}\Delta\rho\big{)}\big{(}\partial_{x}\theta\big{)}.$
(19)

For the second, we heuristically consider $\partial_{x}$ to be of order
$k\ll\xi^{-1}$, then neglect all terms which are necessarily sub-leading under
only the assumptions that $\xi k$ and $\Delta\rho/\rho_{0}$ are small (we are
not making any assumptions for $\theta$). For example,
$\xi^{2}\partial_{x}^{2}\Delta\rho\sim\xi^{2}k^{2}\Delta\rho$ can be neglected
compared to $\Delta\rho$ since $\xi k\ll 1$, and
$\partial_{x}^{2}\Delta\rho\sim k^{2}\Delta\rho$ can be neglected compared to
$\rho_{0}(\partial_{x}\theta)^{2}\sim\rho_{0}k^{2}$ since
$\Delta\rho\ll\rho_{0}$, but both $\Delta\rho/\rho_{0}$ and
$\xi^{2}\partial_{x}^{2}\theta\sim\xi^{2}k^{2}$ must be kept since we have not
assumed anything about how $\Delta\rho/\rho_{0}$ compares to $\xi k$. One then
obtains

$\partial_{t}\Delta\rho=-\frac{\rho_{0}\xi^{2}}{\tau}\partial_{x}^{2}\theta-\frac{\lambda\rho_{0}\xi^{2}}{\tau}\big{(}\partial_{x}\theta\big{)}^{2},$
(20)
$\partial_{t}\theta=-\frac{\xi^{2}}{2\tau}\big{(}\partial_{x}\theta\big{)}^{2}+\frac{\lambda\xi^{2}}{2\tau}\partial_{x}^{2}\theta-\frac{2}{\rho_{0}\tau}\Delta\rho.$
(21)

Now taking a time derivative of Eq. (21) together with Eq. (20) yields a
closed equation for $\theta$. Also note that the right-hand side become small
in the $k\rightarrow 0$, $\Delta\rho\rightarrow 0$ limit. Thus we can treat
$\tau\partial_{t}\sim\tau\omega$ as another small quantity, and neglect terms
which are sub-leading by factors of $\tau\omega$. We finally arrive at the
equation for $\theta$:

$\partial_{t}^{2}\theta=\frac{2\xi^{2}}{\tau^{2}}\partial_{x}^{2}\theta+\frac{2\lambda\xi^{2}}{\tau^{2}}\big{(}\partial_{x}\theta\big{)}^{2}.$
(22)

This is precisely the dispersive KPZ equation given in the main text [Eq.
(6)]. The excellent agreement in comparisons to the GP equation further
confirms its validity.

##  II Example of a diverging solution to dispersive KPZ

Here we give a specific solution to the dispersive KPZ equation which diverges
in finite time. Although the equation does not have a general analytic
solution, it turns out that one family is simply parabolas with time-dependent
coefficients:

$\theta(x,t)=a(t)x^{2}+b(t).$ (23)

The coefficients $a(t)$ and $b(t)$ are prescribed initial values and
derivatives at $t=0$, and the dispersive KPZ equation determines their
subsequent values.

The fact that $\theta(x,t)$ diverges as $x\rightarrow\pm\infty$ makes Eq. (23)
somewhat unphysical, but we will explain how these results imply singularities
for more physical field profiles. For now, merely note that $\theta(x,0)$ is
finite for any fixed value of $x$. We shall show that
$\theta(x,t)\rightarrow\infty$, still at fixed $x$, as $t$ approaches some
finite $\tau_{\textrm{sing}}$.

Inserting Eq. (23) into Eq. (22), we find a closed pair of equations for
$a(t)$ and $b(t)$:

$\frac{\textrm{d}^{2}}{\textrm{d}t^{2}}a(t)=\frac{8\lambda\xi^{2}}{\tau^{2}}a(t)^{2},\quad\frac{\textrm{d}^{2}}{\textrm{d}t^{2}}b(t)=\frac{4\xi^{2}}{\tau^{2}}a(t).$
(24)

The equation for $a(t)$, which gives the curvature of the $\theta$ field, is
equivalent to the equation of motion for a particle in a cubic potential
$U(a)\equiv-8\lambda\xi^{2}a^{3}/3\tau^{2}$. Thus unless the “energy”
$E\equiv\frac{1}{2}(\textrm{d}a/\textrm{d}t)^{2}+U(a)$ is exactly 0, the
particle will “roll downhill” to $a=+\infty$. The time required for $a$ to
diverge is given by

$\tau_{\textrm{sing}}=\int_{a(0)}^{\infty}\textrm{d}a\frac{\textrm{sgn}\big{[}\textrm{d}a/\textrm{d}t\big{]}}{\sqrt{2\big{(}E-U(a)\big{)}}}.$
(25)

Since the integrand goes as $a^{-3/2}$ at large $a$, the integral converges
and thus $\tau_{\textrm{sing}}$ is finite.

The singularity of the parabolic solution actually implies that a family of
initially bounded $\theta$ profiles is singular as well. The initial profile
need only have an interval of sufficiently long length in which it is given by
Eq. (23). This is because dispersive KPZ, like the linear wave equation, has a
finite “speed of light”: the value of $\theta(x_{0},t_{0})$ depends only on
$\theta(x,t)$ for spacetime points in the (past) light cone $|x-x_{0}|\leq
c(t_{0}-t)$, where $c\equiv\sqrt{2}\xi/\tau$; a proof of such causal behavior
can be found in Ref. [27]. As a result, any $\theta(x)$ which is a parabola in
a region of length $2L\equiv 2c\tau_{\textrm{sing}}$, where
$\tau_{\textrm{sing}}$ is given by Eq. (25), will become singular in no more
than time $\tau_{\textrm{sing}}$, regardless of the shape of $\theta(x,0)$
outside the region.

##  III Derivation of hydrodynamic equations from dissipative Gross-Pitaevskii

Here we first derive the hydrodynamic equations in the main text, Eqs. (11)
and (12), and then construct soliton solutions to them. For the derivation, we
no longer assume $\Delta\rho\ll\rho_{0}$ but continue to make the long-
wavelength approximation. In particular, assume that spatial variations in
$\rho$ and the velocity $v\equiv\partial_{x}\theta/m$ are on a significantly
longer length scale than $\xi$. Thus in Eq. (16), written as

$\partial_{t}\rho+\partial_{x}\big{(}\rho
v\big{)}=\frac{\lambda}{m}\sqrt{\rho}\partial_{x}^{2}\sqrt{\rho}-m\lambda\rho
v^{2},$ (26)

the first term on the right-hand side, called the “quantum pressure” [39], can
be dropped since it contains two spatial derivatives. Equation (17) in turn
takes the form

$\partial_{t}v+v\partial_{x}v-\frac{1}{2m^{2}}\partial_{x}\left(\frac{1}{\sqrt{\rho}}\partial_{x}^{2}\sqrt{\rho}\right)+\frac{2U}{m}\partial_{x}\rho=\frac{\lambda}{2m}\partial_{x}\left(\frac{1}{\rho}\partial_{x}\big{(}\rho
v\big{)}\right).$ (27)

Dropping terms containing two and three derivatives, we find Eqs. (11) and
(12) in the main text. Anticipating a solution where the density approaches a
constant $\rho_{0}$ asymptotically and recalling
$\xi\equiv\sqrt{1/mU\rho_{0}}$, $\tau\equiv m\xi^{2}$,
$c\equiv\sqrt{2}\xi/\tau$, these equations can be cast as

$\displaystyle\partial_{t}\rho+\partial_{x}\big{(}\rho v\big{)}$
$\displaystyle=-\frac{\sqrt{2}\lambda}{\xi c}\rho v^{2},$ (28)
$\displaystyle\partial_{t}v+v\partial_{x}v$
$\displaystyle=-\frac{c^{2}}{\rho_{0}}\partial_{x}\rho.$ (29)

Figure 4: Velocity profile of solitons, for three choices of $u/c$. Plotted is
$z$ as a function of $v(z)$, in accordance with Eq. (35). The solution is
valid only when $v(z)$ is single-valued for all $z$ (indicated by a green
background). For certain $u$ and/or ranges of $v$, the solution is not single-
valued and thus invalid (indicated by red).

Soliton solutions to the above can be found by making the ansatz
$\rho(x,t)=\rho(x-ut)$, $v(x,t)=v(x-ut)$, with $u$ to be determined. Denote
$z\equiv x-ut$. Eq. (29) becomes

$-u\frac{\textrm{d}v}{\textrm{d}z}+v\frac{\textrm{d}v}{\textrm{d}z}=-\frac{c^{2}}{\rho_{0}}\frac{\textrm{d}\rho}{\textrm{d}z},$
(30)

which can be integrated to give

$\rho(z)=\rho_{0}\left(1+\frac{u}{c^{2}}v(z)-\frac{1}{2c^{2}}v(z)^{2}\right).$
(31)

To determine the integration constant, we have used the fact that in the
region where the condensate is uniform and stationary, i.e., where
$v(z)\rightarrow 0$, we should have $\rho(z)\rightarrow\rho_{0}$. Since we are
only interested in solutions with $\rho(z)\geq 0$, the allowed range of $v(z)$
is given by

$u-\sqrt{u^{2}+2c^{2}}\leq v(z)\leq u+\sqrt{u^{2}+2c^{2}}.$ (32)

Substituting Eq. (31) into Eq. (28) gives

$\left(1-\frac{u^{2}}{c^{2}}+\frac{3u}{c^{2}}v(z)-\frac{3}{2c^{2}}v(z)^{2}\right)\frac{\textrm{d}v}{\textrm{d}z}=-\frac{\sqrt{2}\lambda}{\xi
c}\left(v(z)^{2}+\frac{u}{c^{2}}v(z)^{3}-\frac{1}{2c^{2}}v(z)^{4}\right).$
(33)

We have the partial fraction decomposition

$\displaystyle\frac{2c^{2}-2u^{2}+6uv-3v^{2}}{2c^{2}v^{2}+2uv^{3}-v^{4}}=$
$\displaystyle\left(2+\frac{u^{2}}{c^{2}}\right)\frac{u}{c^{2}}\frac{1}{v}+\left(1-\frac{u^{2}}{c^{2}}\right)\frac{1}{v^{2}}$
(34) $\displaystyle+$
$\displaystyle\,\frac{\sqrt{u^{2}+2c^{2}}}{\big{(}u+\sqrt{u^{2}+2c^{2}}\big{)}^{2}}\frac{1}{v-u-\sqrt{u^{2}+2c^{2}}}-\frac{\sqrt{u^{2}+2c^{2}}}{\big{(}u-\sqrt{u^{2}+2c^{2}}\big{)}^{2}}\frac{1}{v-u+\sqrt{u^{2}+2c^{2}}}.$

Thus Eq. (33) is integrated to give

$\displaystyle\frac{\sqrt{2}\lambda}{\xi c}\big{(}z-z_{0}\big{)}=$
$\displaystyle\left(1-\frac{u^{2}}{c^{2}}\right)\frac{1}{v(z)}-\left(2+\frac{u^{2}}{c^{2}}\right)\frac{u}{c^{2}}\log{\big{|}v(z)\big{|}}$
(35) $\displaystyle-$
$\displaystyle\,\frac{\sqrt{u^{2}+2c^{2}}}{\big{(}u+\sqrt{u^{2}+2c^{2}}\big{)}^{2}}\log{\Big{|}v(z)-u-\sqrt{u^{2}+2c^{2}}\Big{|}}+\frac{\sqrt{u^{2}+2c^{2}}}{\big{(}u-\sqrt{u^{2}+2c^{2}}\big{)}^{2}}\log{\Big{|}v(z)-u+\sqrt{u^{2}+2c^{2}}\Big{|}},$

where $z_{0}$ is the integration constant. Inverting Eq. (35) gives the
velocity profile $v(z)$.

Not every choice of $u$ yields a valid solution, however. In order for Eq.
(35) to be invertible and $v(z)$ to be defined for all $z$, the curve $z(v)$
must be monotonic and range from $-\infty$ to $\infty$ over an interval of
$v$. Figure 4 gives some examples using $u>0$ (the case $u<0$ is related by
symmetry—take $v\rightarrow-v$, $z\rightarrow-z$ to obtain the solution for
$-u$). We see that $u=c/2$ does not have any valid solution—we have confirmed
that this is the case for any $|u|<c$. Furthermore, the cases $u=c$ and
$u=6c/5$ only have valid solutions for $v<0$. One can confirm that the same is
true for any $u>c$ (e.g., by setting $\textrm{d}z/\textrm{d}v=0$ in Eq.
(33)—there will be two positive solutions for any $u>c$, hence the $z(v)$
curve is qualitatively similar to the rightmost panel of Fig. 4).

To summarize, soliton solutions exist only for $|u|\geq c$ and with the
condensate moving in the opposite direction ($v(z)<0$ for $u>0$ and $v(z)>0$
for $u<0$). In numerical simulations, however, solitons with $|u|>c$ do not
appear, which can be argued on several grounds. First, the long-wavelength
dispersive KPZ equation is causal [27], hence the wavefront cannot cut into
the condensate – in a region still described by the latter equation – at a
speed faster than $c$. Second, such supersonic solitons, while emerging in
many (even relativistic) field theories [49], are not expected to emerge for
physical initial conditions (for example, when the perturbation away from a
uniform condensate is confined to a finite region). Furthermore, even if they
do emerge, they are likely to be unstable due to mechanisms such as the Landau
criterion for the onset of dissipation in a superfluid [38].

##  IV  Effective contact interactions between dark-state Rydberg polaritons

Here we calculate the effective one-body and two-body loss rates for a Rydberg
polariton system given by the Hamiltonian

$H=\int\textrm{d}z\begin{pmatrix}\psi_{e}(z)\\\ \psi_{p}(z)\\\
\psi_{s}(z)\end{pmatrix}^{{\dagger}}\begin{pmatrix}-ic\partial_{z}&g&0\\\
g&-\Delta&\Omega\\\ 0&\Omega&0\end{pmatrix}\begin{pmatrix}\psi_{e}(z)\\\
\psi_{p}(z)\\\
\psi_{s}(z)\end{pmatrix}+\frac{1}{2}\int\textrm{d}z\textrm{d}z^{\prime}\frac{C_{6}}{\big{|}z-z^{\prime}\big{|}^{6}}\psi_{s}(z)^{{\dagger}}\psi_{s}(z^{\prime})^{{\dagger}}\psi_{s}(z^{\prime})\psi_{s}(z),$
(36)

where $\psi_{e}(z)$, $\psi_{p}(z)$, and $\psi_{s}(z)$ are the bosonic field
operators for respectively a photon, atomic $p$ state, and atomic Rydberg
state at position $z$. $c$ is the speed of light, $g$ is the atom-photon
coupling, $\Delta\equiv\delta+i\gamma$ is the complex detuning which takes
into account the finite linewidth $2\gamma$ of the $p$ state, and $2\Omega$ is
the control field Rabi frequency. We shall show that it is possible to tune
parameters so that the two-body loss rate vanishes while the one-body loss
rate (going as $k^{2}$) remains finite.

We shall refer to Refs. [50, 51], and only present the relevant equations
here. First, for $g\gg|\Delta|,\Omega,\gamma$ and small $k$, we can focus on
the middle branch (the dark state) resulting from diagonalization of the non-
interacting part of the Hamiltonian. The effective (complex) mass $m$ of the
dark-state polariton is given by

$m=-\frac{g^{4}}{2\Delta c^{2}\Omega^{2}},$ (37)

hence in the notation of the main text, the one-body ($k^{2}$) dissipation
strength is simply $\lambda=\gamma/|\delta|$ for $\delta<0$. The effective
interaction potential between polaritons is

$U(r)=\frac{1}{\chi}\frac{1}{\frac{|\chi C_{6}|}{\chi
C_{6}}r^{6}/r_{b}^{6}-1},$ (38)

where the blockade radius is $r_{b}\equiv|C_{6}\chi|^{1/6}$, and

$\chi\equiv\frac{1}{2\Delta}-\frac{\Delta}{2\Omega^{2}}.$ (39)

For coherent interactions, the two-body problem can be solved analytically in
the limit $\kappa r_{b}|\Delta/2\Omega^{2}\chi|\ll 1$ [51, 50]. The analysis
carries over to the more general dissipative case considered here. Assuming
weak interactions, we can approximate the effective interaction potential by a
square well of width $r_{b}$ and depth $u_{0}^{2}$ chosen to match

$\vspace{-0em}u_{0}^{2}r_{b}=-\int_{0}^{\infty}\text{d}r\,mU(r)=-\frac{\pi}{3}r_{b}\kappa^{2}\frac{|\Delta|^{2}\left(\frac{\Delta}{\Delta^{2}-\Omega^{2}}\right)^{5/6}}{\Delta\Omega^{1/3}\sqrt[6]{\left|\frac{1}{\Delta}-\frac{\Delta}{\Omega^{2}}\right|}},\vspace{-0.em}$
(40)

where $\kappa=g^{2}/c|\Delta|$. In this case, we can find a complete analytic
solution to the scattering phase for all $k$ [52]. The contact interaction
between polaritons in a dilute regime can be characterized by the strength
$-1/ma$ (denoted $U$ in the main text) with

$\vspace{-0em}a=r_{b}+\frac{1}{u_{0}\tan(u_{0}r_{b})}.\vspace{-0.em}$ (41)

More generally, however, one must solve the two-body problem numerically. In
either event, by choosing appropriate values for $\kappa r_{b}$ and $\delta$,
one can hope to access a regime in which $ma\in\mathbb{R}$, i.e., interactions
are lossless.

Figure 5: Strength of contact interactions $-1/ma$ (denoted $U$ in the main
text) between two dark-state polaritons for $\delta/\gamma=-10$,
$\Omega/|\Delta|=6$, and $C_{6}>0$. Solid blue and dashed orange curves are
respectively the real and imaginary parts within the square-well
approximation, while dotted red and dot-dashed green curves are those of the
full numerical solution. Using the numerical solution, we see a zero-crossing
of Im$(1/ma)$ for $\kappa r_{b}\approx 5.5$.

We have found that such a regime exists, and that furthermore: i) this occurs
when $\Omega>|\Delta|$, and thus the effective interactions are repulsive
[51]; and ii) $|\delta|\gg\gamma$, and thus the effective mass has a small
(but nonzero) imaginary part. These are the ideal conditions for the phenomena
observed in the main text.

Figure 5 shows a representative example, where $\delta/\gamma=-10$,
$\Omega/|\Delta|=6$, and $C_{6}>0$. We see the desired zero-crossing of
Im$(1/ma)$ for $\kappa r_{b}\approx 5.5$. We have also compared the analytic
square-well results ($-1/ma_{\textrm{sq}}$) against the numerical solution of
the two-polariton scattering problem ($-1/ma_{\textrm{num}}$). As expected,
the two agree at small $\kappa r_{b}$. The full numerical solution is required
to identify the zero-crossing of the imaginary part.