Stationary quantum vortex street in a driven-dissipative quantum fluid of light

TL;DR

This paper explores the formation of a novel class of density-phase defects, including vortex streets, in a driven-dissipative 2D quantum fluid of light, revealing stable configurations and potential for complex potential implementations.

Contribution

It introduces the concept of stationary quantum vortex streets in a driven-dissipative quantum fluid of light, demonstrating their formation, stability, and potential for complex potential engineering.

Findings

Vortex streets form in 2D quantum fluid of light under bistability.

Dark soliton molecules evolve into vortex street configurations.

Flexible potential design enables maze-like vortex connections.

Abstract

We investigate the formation of a new class of density-phase defects in a resonantly driven 2D quantum fluid of light. The system bistability allows the formation of low density regions containing density-phase singularities confined between high density regions. We show that in 1D channels, an odd (1-3) or even (2-4) number of dark solitons form parallel to the channel axis in order to accommodate the phase constraint induced by the pumps in the barriers. These soliton molecules are typically unstable and evolve toward stationary symmetric or anti-symmetric arrays of vortex streets straightforwardly observable in \emph{cw} experiments. The flexibility of this photonic platform allows implementing more complicated potentials such as maze-like channels, with the vortex streets connecting the entrances and thus solving the maze.

Equations

$$i\hslash \frac{\partial \psi }{\partial t}=\left[-\frac{{\hslash }^2{\nabla }^2}{2m}-i\Gamma +g{\left|\psi \right|}^2\right]\psi +(S+P)e^{-i{\omega }_{0}t},$$

$$S_c\approx \frac{2{\left(\hslash {\omega }_0\right)}^{3/2}}{3\sqrt{3}{g}^{1/2}}$$

$$v\sim 2\frac{S-S_c}{S_c}\frac{\xi }{\tau },$$

$$L(\mathbf{r})u(\mathbf{r})+g\psi (\mathbf{r}{)}^{2}v(\mathbf{r})$$

$$L(\mathbf{r})v(\mathbf{r})+g({\psi }^{\ast }(\mathbf{r}){)}^{2}u(\mathbf{r})$$

$$\left(-\hslash {\omega }_0-\frac{{\hslash }^2}{2m}\Delta +g{\left|{\psi }_0\right|}^2-i{\Gamma }_0\right){\psi }_0+S=0$$

$$m_0\frac{d^{2}x}{d{t}^2}=F(x)$$

$$F(x)=g{x}^3-\delta x+S$$

$$U(x)=-\int F(x)dx=-\frac{g{x}^4}{4}+\frac{\delta x^2}{2}-Sx$$

$$U(x_1)=U(x_3),$$

$$\underset{x_1}{\overset{x_3}{\int }}F(x)dx=0$$

$$S_c=\frac{2{\delta }^{3/2}}{3\sqrt{3}{g}^{1/2}}$$

$$\frac{S_c}{S_{min}}=\frac{2}{3\sqrt{3}}\frac{\delta }{\Gamma }$$

$$\begin{gathered} \frac{S_c({\omega }_0+\Delta {\omega }_0,\Gamma +\Delta \Gamma ,g+\Delta g)}{S_c({\omega }_0,\Gamma ,g)}=1+ \\ {c}_{\omega }\frac{\Delta {\omega }_0}{{\omega }_0}+c_{\Gamma }\frac{\Delta \Gamma }{\Gamma }+c_g\frac{\Delta g}{g}. \end{gathered}$$

$$S_c=\sqrt{\frac{\hslash {\omega }_0}{g}}\sqrt{{\left(\frac{3\hslash {\omega }_0}{32}\right)}^2+{\Gamma }^2}$$

$$v=\frac{\partial n}{\partial t}\frac{\Delta x}{\Delta n}$$

$$\frac{\partial \psi }{\partial t}=\frac{S-S_c}{i\hslash }$$

$$\psi \left(t\right)=\frac{S-S_c}{i\hslash }t+\sqrt{n_0}{e}^{i{\varphi }_0}$$

$$\frac{\partial n}{\partial t}=\frac{\partial }{\partial t}{\left|\psi \left(t\right)\right|}^2\approx 2\sqrt{n_0}\frac{S-S_c}{\hslash }\frac{1}{\sqrt{2}}\left(1-\chi \frac{S-S_c}{S_c}\right)$$

$$v=\sqrt{2}\sqrt{n_0}\frac{S-S_c}{\hslash }\frac{\Delta x}{\Delta n}$$

$$v\sim 2\frac{S-S_c}{S_c}\frac{\xi }{\tau }$$

$$i\hslash \frac{\partial \psi }{\partial t}=-\frac{{\hslash }^2}{2m}\Delta \psi +g{\left|\psi \right|}^2\psi -i\Gamma \psi +P(x,y)e^{-i{\omega }_{0}t}$$

$$\hslash {\omega }_0{\psi }_0=-\frac{{\hslash }^2}{2m}\Delta {\psi }_0+g|{\psi }_0{|}^2{\psi }_0-i\Gamma {\psi }_0+P$$

$$\left({\begin{array}[]{*{20}{c}}{-\frac{{{\hbar^{2}}}}{{2m}}\left({\frac{{{\partial^{2}}}}{{\partial{x^{2}}}}-k_{y}^{2}}\right)+2g{{\left|{{\psi_{0}}}\right|}^{2}}-i\Gamma-\hbar{\omega_{0}}-\hbar\omega}&{g\psi_{0}^{2}}\\ {g\psi_{0}^{*2}}&{-\frac{{{\hbar^{2}}}}{{2m}}\left({\frac{{{\partial^{2}}}}{{\partial{x^{2}}}}-k_{y}^{2}}\right)+2g{{\left|{{\psi_{0}}}\right|}^{2}}-i\Gamma-\hbar{\omega_{0}}+\hbar\omega}\end{array}}\right)\left({\begin{array}[]{*{20}{c}}A\\ B\end{array}}\right)=\left({\begin{array}[]{*{20}{c}}0\\ 0\end{array}}\right)$$

$$-\frac{{\hslash }^2}{2m}\frac{{\partial }^{2}A}{\partial x^2}+g|{\psi }_0{|}^{2}A-\hslash {\omega }_{0}A\approx -\frac{\hslash {\omega }_0}{2}A$$

$${\left(\frac{{\hslash }^2{k}_y^2}{2m}+g{\left|{\psi }_0\right|}^2-i\Gamma -\frac{\hslash {\omega }_0}{2}\right)}^2-g^2{\left|{\psi }_0\right|}^4={\hslash }^2{\omega }^2$$

$${\Gamma }_{max}=\frac{\hslash {\omega }_0}{2}-\Gamma$$

$$k_y^{\ast }=\sqrt{\frac{m{\omega }_0}{\hslash }}$$

$$k_y^{\ast }=\frac{1}{\xi \sqrt{2}}$$

$${\hslash }^2{\omega }^2={\left(g|{\psi }_0(x){|}^2-E_n\right)}^2-g^2|{\psi }_0(x){|}^4$$

$$\Im (\hslash \omega )=E_n$$

Full text

A fluid is called quantum when it exhibits quantum-mechanical effects at a macroscopic scale. For bosons, it occurs when many particles can be described by a single-particle wave function. This collective behaviour can arise spontaneously when particles undergo a phase transition towards a quantum-coherent state such as a superconducting state, a superfluid, or a Bose-Einstein condensate (BEC). Cavity exciton-polaritons kavokin2 (polaritons) are 2D photonic modes interacting via their excitonic parts. Their quantum coherence can spontaneously occur through a BEC process kasprzak2006bose , but a unique feature of this photonic system is that such coherence can also be imprinted by a resonant laser and be preserved substantially longer than their lifetime wertz2010spontaneous . The high control of the injected flow combined with a direct optical access to the full wave function (amplitude, phase, space, time) make this platform very attractive to study quantum fluid physics RevModPhys.85.299 . A typical example which revealed the potential of this system is the observation of oblique dark solitons pigeon2011Hydrodynamic ; amo2011polariton ; Hivet , which form when a supersonic quantum fluid hits a defect (initially proposed in 2006 PhysRevLett.97.180405 for atomic BECs). The 2D solitons forming behind the defect remain stable because the transverse ”snake instability”KIVSHAR2000117 ; Dutton2001 , making 2D solitons normally unstable, is carried away by the supersonic flow making the soliton effectively 1D PhysRevLett.100.160402 ; Walker2017 . Such supersonic flow is energetically unstable, but polariton flows can be efficiently decoupled from thermal relaxation wertz2010spontaneous ; Wertz2012 , which made possible the observation of oblique dark solitons amo2011polariton ; Hivet .

An interesting regime occurs if the fluid velocity is decreased just below the speed of sound. In such a case, the subsonic flow still interacts with the defect exhibiting a local acceleration. This leads to the formation of quantum vortex streets composed of vortex-antivortex pairs. This quantum version of von Karman vortex streets can be understood as the decay of the oblique solitons via the snake instability when the stabilization by the supersonic flow is lots. The creation of vortex-antivortex pairs has been reported in time-resolved pulsed experiments, both in polaritons sanvitto2011all ; nardin2011hydrodynamic and atomic quantum fluids Kwon2016 . Similar proposals were made for non-resonantly pumped polaritons Smirnov2014 ; Liew2015 . However, the study of the snake instability dynamics leading to quantum vortex streets requires both cw excitation and time resolution, and it remained elusive so far. In a recent theoretical work, it was proposed to improve this scheme by sustaining the propagating flow against radiative decay by using a support laser covering the whole sample pigeon2017sustained . Interestingly, this configuration demonstrates original density-phase defects. This pump-support scheme is not limited to the study of the flow scattering on defects, but can be used in a more general frame to create and study a large variety of topological defects Sich2012 ; Tercas2014 ; Sich2018 . Topological defects in the driven-dissipative case can be stationary only when the support laser intensity falls in the bistability loop of the non-linear system Baas2004Opticalbistability ; Yulin2008 ; Zhang2012 ; pigeon2017sustained , where the density can be either low or high, depending on the laser absorption. Stationary phase defects exist in low-density regions, where the phase is not fixed, because most particles are not directly injected by the laser but diffuse from higher-density regions. The control of the spatial distribution of intensity and phase allows to realize various confining potentials, such as 1D channels, 0D traps PhysRevB.89.134501 , or circuits made by the combination of both.

In this work, we show that the creation of narrow 1D low-density channels surrounded by high-density regions leads to the formation of stationary dark soliton molecules. These molecules exhibit snake instabilities, leading to the formation of coupled chains of vortex-antivortex pairs – stationary vortex streets, stabilized by the confining potential. We finally show that the pump distribution can be used to create a maze. Vortex streets retreat from the dead ends, while the entrance and the exit of the maze remain connected, implementing an efficient analog all-optical maze solving algorithm.

The resonantly pumped microcavity is modelled by the standard driven-dissipative Gross-Pitaevskii equation, formally equivalent to the Lugiato-Lefever equation Lugiato1987 . We neglect the polarization degree of freedom, the non-parabolicity of the polariton dispersion, and any thermal effects Vishnevsky2012 ; Stepanov2018 . The equation reads:

[TABLE]

where $\Gamma=\hbar/(2\tau)$ is the polariton decay rate ($\tau=15$ ps), $m=8\times 10^{-5}m_{0}$ is the polariton mass ($m_{0}$ is the free electron mass), $g=5~{}\mu$eV$\mu$m2 is the polariton-polariton interaction constant. The detuning between the ground state and the pump laser is $\omega_{0}=0.14~{}\mathrm{meV}/\hbar$. The support $S$ and the pump $P(\mathbf{r})$ are at normal incidence.

We first consider spatially homogeneous pumping (support only). The bistability loop obtained for our parameters is shown in Fig. 1(a). Next, we add a half-space pump ($x<30$ $\mu$m) switching the system to the higher branch of the bistability loop, whereas the other half-space remains on the lower branch. The use of a spatially inhomogeneous pump $P(\mathbf{r})$ allows to control the pumping intensity in the high- and low-density regions independently.

Figure 1(b) shows the two regions separated by a domain wall (DW). The region under the pump shows a large intensity and a fixed phase. The intensity at the DW decays within one healing length $\xi=\hbar/\sqrt{2gnm}$ ($n$ is the density) and then exhibits small periodic oscillations of intensity. The DW is stable against the development of instabilities along $Y$, but can propagate along $X$. Such DW propagation has been previously considered for polaritons kavokin2008neurons ; Amo2010Light and in general for switching waves in optics Rozanov1997 ; Ganne2001 ; odent2016optical and beyond. The velocity $v$ of the DW computed numerically as a function of support intensity is shown in Fig. 1(a) (black). As expected for this class of differential equations, $v$ is linear in $S-S_{c}$, where the critical value of the support $S_{c}$ is given, in the limit $\Gamma\to 0$, by the Maxwell construction Maxwell1875 ; suppl

[TABLE]

When the support is larger than $S_{c}$, the DW propagates to the right with a velocity

[TABLE]

where $\xi\approx 1.8\mu$m, and the high intensity region expands to the whole space. For support values below $S_{c}$ the high-density region shrinks and the DW stops at the boundary of the pump ($x=30\mu$m). Around $S_{c}$, spatially localized solutions of the Gross-Pitaevskii equation bifurcate under the form of dark solitons multiplets Parra2016 .

Next, we consider a second high-density region with its boundary parallel to the first one, defining an all-optically controlled confining potential. Similar configuration, but without the support beam, has been studied in Aioi2013 . We start by considering high density regions with the same phase and a fixed channel width $L=23$ $\mu$m ($\sim 13\xi$ of the high density region). Fig. 2 is computed at $S=0.25S_{c}$ for 2 values of $P$. The first column presents the stationary intensity distribution with 2 dark solitons in the channel. The system is effectively 1D, since it is translationally invariant along $Y$ (periodic boundary conditions). Dark solitons are anti-symmetric states with a $\pi$ phase shift. The phase constraints imposed by the high-density regions therefore only allow an even number of solitons. These soliton multiplets are however unstable with respect to the development of instability along $Y$ for a large range of parameters. The second column of Fig. 2 shows the imaginary part of the energy of the weak excitations versus their longitudinal wave vector $k_{y}$, obtained from the Bogoliubov-de Gennes equations:

[TABLE]

where $L(\mathbf{r})=-\hbar^{2}\nabla^{2}/2m+2g|\psi(\mathbf{r})|^{2}-\hbar\omega_{0}-i\Gamma$. Comparing with Ref. morgan2013stability , the chemical potential $\mu$ for conservative system is replaced by the laser frequency in the driven-dissipative case carusotto2004probing ; solnyshkov2008dispersion . Eqs. (4) represent an eigenvalue problem for $\omega$. The translational invariance along $Y$ allows to replace $-\nabla^{2}\rightarrow k_{y}^{2}-\partial^{2}/\partial x^{2},$ where $k_{y}$ is the wave vector of perturbations. Positive imaginary part of the energy leads to the development of modulational instability (the snake instability, well known in conservative condensates Kuznetsov1988 ). The maximal instability wave vector can be estimated suppl as $k_{y}^{*}=1/(\sqrt{2}\xi)$. The lowest energy mode in the double-well potential formed by the soliton pair is a symmetric bound state, whereas the highest mode is an anti-symmetric anti-bound state. However, the symmetry of the patterns in the 3rd column of Fig. 2 is inversed, because each soliton is anti-symmetric by itself ($\pi$ phase jump). An antisymmetric superposition of solitons is therefore a symmetric function, which gives rise to the symmetric pattern observed in the 2nd line of Fig. 2. A more detailed study of the mode competition is shown in suppl , section II. The mode with a higher imaginary part develops faster (marked by arrows in Fig. 2). The 2D modulational instabilities in the numerical simulations can be triggered by any noise or fluctuations breaking the translational symmetry along the $Y$-axis. Here, we consider a weak Gaussian disorder with a correlation length of $2~{}\mu$m and an amplitude $\gamma=0.01$ meV. The third and fourth columns of Fig. 2 show the intensity and the phase of the stationary wave function after the development of the instability. The precise realization of the disorder determines the positioning of the pattern along $Y$, but does not affect the shape, at least if the disorder amplitude is sufficiently small: $\gamma\ll\hbar\omega_{0}$. Additional simulations shown in suppl confirm that these stationary patterns are accessible considering realistic disorder parameters up to 0.1 meV. In all cases, the solitons break into two vortex anti-vortex chains, which can be seen as stationary vortex streets. 2D analysis confirms the stability of the final patterns. It means that the snake instability develops, but is then frozen by the presence of the confining potential.

Figure 3 shows a phase diagram obtained from the stability analysis versus pump and support intensities (constant $L$ and $\hbar\omega_{0}$). The phase in the upper right corner (purple) corresponds to high density in the channel (no solitons or vortices). It occurs for support values typically smaller than $S_{c}$, because of the particle flow coming from the two DWs instead of one. This regime is qualitatively similar with the one of the polariton neuron picture kavokin2008neurons . The dark grey region corresponds to a non-stationary steady state (limit cycle), at least in the conditions of our simulations, namely without energy relaxation and for sufficiently low disorder. This phase shows a pair of breathing solitons oscillating in time (see movie suppl ). The small light-grey domain corresponds to a lattice of four solitons. This occurs for small $P$ and large $S$, so weak transverse flows which favours the soliton lattice stability. The next phase located at the bottom left corner corresponds to the collapse of 4 solitons into a symmetric pair of vortex chains (see Fig. S1 suppl ). The two next phases located above the blue line correspond to the collapse of 2 solitons into symmetric and anti-symmetric vortex chains respectively (as in Fig. 2(a,b)). A tiny domain (lime-green) exhibits the collapse of 4 solitons into an anti-symmetric pair of vortex chains. The false color scale of the figure shows the maximal instability wave vector $k_{y}^{*}$ (except for the non-stationary and stable phases). The anti-symmetric solutions have a twice larger period than the symmetric ones. The $k_{y}^{*}$ gradient within a given phase is relatively small, which means that the patterns visibility should not be strongly affected by pump/support intensity fluctuations in a real experiment. Disorder broadens the transitions between the phases, but the core regions remain well defined (see suppl for details).

Another interesting possibility offered by this driven-dissipative system is to tune the relative phase of the pump between the walls, working at zero support. In such a case, an odd number of solitons forms suppl , decaying into the same number of chains of vortex-antivortex pairs. This tuning of the soliton number by varying the relative phase between the pumps is a generalization of Goblot2016 to 2D. Here, the 2D character of the system allows instabilities along $Y$. In fact, modern optical techniques allow creating any shape of confining potential, such as various 0D traps PhysRevB.89.134501 , or graphs, connecting 0D sites with 1D channels.

The geometry we address now is a maze of 1D channels [Fig. 4(a,b)]. For a proper value of $(S,P)$, immediately after the jump of the walls on the upper branch, the maze is filled with solitons. However, the dead ends represent a configuration different from that of Fig. 2: the heads of the vortex streets start to withdraw. Fig  4(a) shows the intensity distribution 20 ps after driving pump and support are switched on (blue circle in Fig. 3 in the non-stationary phase). The heads are moving as symbolized by the arrows. Fig  4(b) shows the final intensity distribution ($t=1$ ns), where the street only connects the two exits of the maze. A zoom on the vortex street head is shown in Fig. 4(c). This head is a DW, but the conditions for its motion are different from the lateral DW in an open 1D channel we considered previously. Indeed, the motion of the head-DW along $X$ is facilitated by the confinement and the flow along $Y$ coming from the two lateral walls. When the head-DW arrives at a cross-roads, an open 1D channel geometry is restored and the DW stops. Fig. 4(d) shows the critical support intensity for the motion of the two types of DW (dead/open ends – blue/red). Both decrease when the corridor becomes narrower, tending to zero for $L\approx 14$ $\mu$m. This is the minimal width of a finite quantum well (of $0.14$ meV depth, given by the laser detuning), where the two-node state can exist for such particle mass. For a wide open channel (red), $S$ can even exceed $S_{c}$, because the quantum pressure of the multinode structure prevents the DWs from meeting each other (see suppl for details). As expected, there is a substantial support range where the head-DW propagates, whereas an open 1D corridor remains stable. Within this range, this configuration represents an optical maze solver (see also a supplementary movie for the dynamics suppl ). The head-DW motion can be affected by disorder, but being an extended object, it is less sensitive than vortices which easily pin on defects Lagoudakis2009 . Numerical simulations suppl show that disorder trapping can be avoided by working at higher detunings $\hbar\omega_{0}>3\gamma$ and, if necessary, shorter lifetimes, making the effects observable in realistic systems.

All-optical maze solving is important for the large interdisciplinary field of analog graph solving algorithms Shannon1951 ; Steinbock1995 ; Caruso2016 ; Berloff2017 . The solving time is determined by the velocity $v$ of the head-DW (see suppl for details). In the worst case, the length of the dead end is $NL$, where $N$ is the number of cells in the maze (or vertices in the graph) and $L$ is the width of a channel. If $Z$ is the overall system size, the maximal number of cells is $N=Z^{2}/L^{2}$. The solving time is therefore $t=NL/v$. The best among the other maze solving (pathfinding) algorithms Even2011 such as the Depth-First Search also exhibit the worst-case complexity of $O(N)$, but the practical advantage of the present analog implementation is the small value of the prefactor $L/v\sim 0.5$ ns: the high velocity $v$ reduces the solving time, allowing such an analog maze solver to outperform a modern PC which needs hundreds of clock ticks to check a single cell.

To conclude, we have shown that the modulational instability can be controlled and stabilized in a driven-dissipative polariton system allowing the on-demand formation of soliton molecules and vortex streets exhibiting a particularly rich phase diagram. Non-stationary regimes can be used for fast analog maze solving.

I Supplemental Materials

In this Supplemental Material, we provide additional details on the results shown in the main text. Section I discusses the behavior of isolated domain walls. Section II is devoted to the stability of soliton structures in channels. Section III provides additional results for channels with $\pi$ phase difference. Section IV discusses the behavior of solitons in a maze. Section V describes the two supplementary movies.

II Stationary and propagating domain walls

The behavior of bistable systems and, in particular, spatially inhomogeneous solutions present in such systems were a subject of studies for a long time. An overview on optical systems can be found in Ref. Rozanov1997 .

We start by looking for a stationary spatially inhomogeneous solution where a part of the system is at the upper bistability branch, whereas another part is at the lower bistability branch. These two parts are separated by a domain wall (or switching wall). This problem can be solved analytically for negligibly small dissipation $\Gamma$ (as compared with the laser detuning $\delta=\hbar\omega_{0}$), and then the solution can be generalized to non-zero $\Gamma$.

The stationary driven-dissipative Gross-Pitaevskii equation that needs to be satisfied by the wave function reads:

[TABLE]

where we are going to neglect the term $\Gamma$ at first. Once it is neglected, the wave function $\psi$ and the pumping (support) $S$ can be assumed to take only real values without loss of generality (since all coefficients are real). The equation (5) can therefore be rewritten as

[TABLE]

which is a Newton’s equation of motion for a material point ($x=\psi$) with a mass $m_{0}=\hbar^{2}/2m$ under the effect of a position-dependent ”force”

[TABLE]

to which one may attribute a ”potential”

[TABLE]

The two maxima of this potential located at the coordinates $x_{1}$ and $x_{3}$ correspond to the two stable domains (high density and low density), while the minimum located at $x_{2}$ corresponds to the inaccessible part of the bistability curve. The system is stationary only if the values of the effective potentials at the two maxima are exactly the same:

[TABLE]

otherwise the domain wall starts to propagate. Indeed, a material point should start its motion at one maximum and finish at the other maximum, and for this the two maxima have to be at the same potential height. The points $x_{1}$ and $x_{3}$, corresponding to the extrema of $U(x)$, can be found analytically from the cubic equation $F(x)=0$, and the condition (9) can be rewritten as

[TABLE]

corresponding to the Maxwell construction in thermodynamics Maxwell1875 . Solving this equation analytically for the unknown $S$ gives finally:

[TABLE]

and the ratio with respect to the pumping required for bistability $S_{min}=\Gamma_{0}\sqrt{\hbar\omega/g}$ is given by

[TABLE]

These analytical expressions presented in the main text can be used as estimates for the required support level. A more precise result can be obtained only numerically, as discussed below.

For non-negligible $\Gamma$, all terms in the equation (5) become comparable. The real terms (the kinetic and the interaction energy, the detuning) are of the order of $\delta=\hbar\omega_{0}$, while the only imaginary term is $\Gamma$ (imaginary part of kinetic energy is small with respect to $\Gamma$). Thus, the expression for the critical pumping $S_{c}$ at $\Gamma$ comparable with $\hbar\omega_{0}$ can be sought by replacing of the first three terms in (5) by $\hbar\omega_{0}$ with some coefficient. This coefficient can be found from numerical simulations by the small variation of parameters $\omega_{0}\rightarrow\omega_{0}+\Delta\omega_{0}$, $\Gamma\rightarrow\Gamma+\Delta\Gamma$, and $g\rightarrow g+\Delta g$, which allows obtaining the coefficients ($c_{\omega}\approx 0.75$, $c_{\Gamma}\approx 0.75$, $c_{g}=-0.5$) in the Taylor expansion:

[TABLE]

As a net result, one obtains

[TABLE]

One sees that this equation has the same Taylor expansion as Eq. (13) for actual values of system parameters.

For $S\neq S_{c}$, one of the domains becomes more favorable than the other, and the domain wall starts to propagate. Assuming that the density changes linearly with position across the domain wall (which is valid in the vicinity of the inflection point of this wall), we find the following expression for the speed of the wall:

[TABLE]

where $\Delta x$ is the width of the domain wall and $\Delta n$ is the difference of the densities in the high and low density regions.

To find the derivative $\partial n/\partial t$, we define $n_{0}$ as the density at the inflection point which corresponds to $S_{c}$. The only contribution into $\partial\psi/\partial t$ at $S\neq S_{c}$ can come from the difference $S-S_{c}$ (other terms in the time-dependent driven-dissipative Gross-Pitaevskii equation cancel):

[TABLE]

which allows to write the solution

[TABLE]

This expression strongly depends on the phase $\phi_{0}$ of the wave function at the inflection point (with respect to the phase of the pump). In a homogeneous system, the phase is given by $\tan\phi=-\Gamma/(gn-\hbar\omega_{0})$. In the low-density region, $\phi\to 0$, whereas in the high-density region $\phi\to\pi/2$ just above the threshold density $gn\approx\hbar\omega_{0}$. We assume that $\phi_{0}$ takes an intermediate value $\phi_{0}=\pi/4$, and we can also assume that it changes linearly with $S-S_{c}$ for small deviations from $S_{c}$, in which case the time derivative of the density can be found as:

[TABLE]

For $S\approx S_{c}$, the second-order correction can be neglected, but it starts to become important for larger deviation of the pumping, creating an important difference between the speed of the wall in the case of shrinking or expanding high-density domain.

In the simplest case, the speed can be found as

[TABLE]

By estimating $\sqrt{\Delta n}\approx\sqrt{n_{0}}\approx S_{c}/\Gamma$ and taking $\Delta x=\sqrt{2}\xi$ one rewrites Eq. (20) as

[TABLE]

The expression (20) for the values used in numerical simulations gives $v=12$ $\mu$m/ns, very close to the numerical value of 30 $\mu$m/ns for $(S-S_{c})/S_{c}=0.05$.

III Stability of driven-dissipative solitons in 2D

We analyze the soliton stability in the 2D driven-dissipative configuration using the Bogoliubov-de Gennes equations for weak excitations of the condensate. Let $\psi$ be the non-trivial solitonic solution of the driven-dissipative Gross-Pitaevskii equation written on the polariton basis

[TABLE]

We can write $\psi=\psi_{0}(x,y)e^{-i\omega_{0}t}$, and the perturbed solution is $\psi_{0}(x,y)+A(x)e^{i(k_{y}y-\omega t)}+B^{*}(x)e^{-i(k_{y}y-\omega t)}$, where $\omega$ is the perturbation frequency relative to the laser frequency $\omega_{0}$ and $\psi_{0}(x,y)$ is the solution of the stationary equation

[TABLE]

Linearizing (21), we obtain the Bogoliubov-de Gennes equations for $A$ and $B$:

[TABLE]

where we are interested in the bogolon states confined along the $X$ direction.

III.1 Maximal instability wave vector for a single soliton

For a single state quantized in a potential trap formed by the profile of a single dark soliton, its energy (which is not yet the bogolon frequency) can be estimated as Pitaevskii :

[TABLE]

which supposes that the system is above the bistability threshold and that the soliton is almost dark. The diagonalization of the matrix (23) thus gives an equation

[TABLE]

which allows to estimate the value of the imaginary part of $\omega$. Indeed, if we suppose $\Gamma\ll\hbar\omega_{0}$, the solution for $\omega$ becomes imaginary if the real part of (25) is negative. Maximal imaginary part is achieved if the first square is zero and the second square is maximized, which is obtained at a point $x_{0}$ where $g|\psi_{0}(x_{0})|^{2}=\hbar\omega_{0}/2$, determining the maximal possible positive imaginary part of $\hbar\omega$ as

[TABLE]

In this case, the soliton is clearly always unstable, because the negative contribution to the imaginary part is much smaller than the positive one, as required for the bistability. We can conclude that 2D solitons in the driven-dissipative configuration remain unstable with respect to small perturbations, at least if they are obtained at the upper bistability branch.

This also provides a criterion for the maximal wavevector $k_{y}$ of the unstable region:

[TABLE]

which is allowed at a different point $x_{1}$, where $g|\psi_{0}(x_{1})|^{2}\approx 0$. At the upper part of the bistability curve, where $\hbar\omega_{0}=gn$, this expression can be linked with the inverse of the healing length of the quantum fluid $\xi$:

[TABLE]

This expression is given in the main text.

III.2 Mode competition

If several quantized states are present in the trap with the energies $E_{n}$, the highest energy state forms the bogolon with the highest imaginary part. Indeed, the imaginary part of $\omega$ from the simplified equation

[TABLE]

is due to the second term, which is negative. The minimal value of the first term is zero. It is achieved when $g|\psi_{0}(x)|^{2}=E_{n}$, which determines the value of $g|\psi_{0}(x)|^{2}$, and the imaginary part is therefore simply

[TABLE]

measured from the energy of the bare states. Our conclusion on the mode competition in the trap is therefore that the highest mode in the trap always wins (developing first the instability), as soon as it is sufficiently well localized in it. A more exact result can only be obtained by numerical solution. In numerical analysis, we start by finding the stationary solution of the GPE $\psi_{0}(x)$, which is then plugged into (23). The imaginary parts of the bogolon energy $\hbar\omega$ found by this procedure are plotted in Fig. 2 of the main text, and the winning mode is shown with an arrow.

Figure 5 shows an additional panel with the same elements as Fig. 2 of the main text, corresponding to the lowest pumping case of the phase diagram in Fig. 3 of the main text (orange point), when 4 solitons are formed in the stationary unperturbed regime. These 4 solitons, however, decay into 2 vortex streets once the modulational instability develops.

III.3 Effects of disorder

In experiments, the patterns arising from the development of the instability are affected in space by the pinning induced by the disorder PhysRevB.88.201303 .

Large disorder should lead to a smoothing of the phase diagram. To test the robustness of the phases discussed in the main text with respect to the disorder, we have plotted an alternative phase diagram (Fig. 6) where the disorder amplitude is 0.1 meV (as compared with 0.01 meV in the main text). When disorder is present, the patterns are not as regular as in a quasi-ideal case. One can nevertheless discriminate a symmetric and an anti-symmetric pattern using the standard deviation of the X coordinate of the solitons mass center $X_{c}(y)$. This deviation reads as $\left<\Delta^{2}\right>=\left<(X_{c}(y)-\left<X_{c}(y)\right>_{y})^{2}\right>_{y}$, where $\left<\right>_{y}$ means averaging over the $Y$ coordinate. For symmetric patterns, the center is nearly aligned vertically, while for the anti-symmetric patterns the center undergoes sinusoidal trajectory and thus the deviation is larger. Computing or measuring experimentally this quantity allows to clearly distinguish different phases, as shown in the insets of Fig. 6. The limit cycle with oscillating solitons can still persist in the presence of disorder. This motion results in their broadening and ”shallowing” as clearly seen in the right inset of Fig. 6. Importantly, the disorder can also suppress the motion of the solitons and thus the limit cycle phase region becomes much smaller with stronger disorder. The diagram was calculated with randomly generated noise for each pair ($S,P$), which enhances the sampling.

IV Solitons under $\pi$ phase shift

In the main text, we show the results obtained when the phase of the pumping laser is homogeneous in space, and only the density profile is varying (allowing to obtain the high-density walls). In this subsection, we present additional results concerning the formation of solitons and their stability for a $\pi$ phase difference between the pump at the walls. In this case, no support is used (otherwise it would exhibit different interference with the two pumping lasers of different phase).

Figure 7 shows the results obtained in this configuration, with top and bottom rows corresponding to two different values of the pump $P$. The first column shows the spatial density profile with 3 or 1 solitons, depending on the distance available for them because of the pump-induced broadening of the walls. The second column shows the imaginary part of the energy obtained from the Bogoliubov-de Gennes analysis described above. Both curves exhibit a maximum with positive imaginary part, confirming the existence of modulational instability. The final stage of the development of this instability is shown in the 3rd column: it exhibits either 3 vortex chains or a single vortex chain. Finally, the phase distribution shown in the 4th column confirms the formation of vortices and anti-vortices evidenced by the density shown in the 3rd column.

V Solitons in a maze

In Figs. 1, 2, 3 of the main text, we were considering a system infinite in the $Y$ direction (implemented by periodic boundary conditions). As soon as we consider a dead end instead of an infinite channel (see Fig. 4(c) of the main text), the ”head” of the pair of solitons comes into play. This head can also be considered as a vertical domain wall between the domains. This wall is between a high-density region (the wall limiting the dead end) and a low-density region (the region with solitons). There are now 3 domain walls in the system: two of them (horizontal) are opposite to each other, whereas one (vertical) can propagate, depending on the conditions. The extra kinetic energy appearing because of the variation of the wavefunction in the transverse direction (vertical), absent in the infinite system (for a very large channel), changes the conditions of the local bistability loop. The flows toward the low density region are also larger and both effects together (kinetic energy and flows) make the low-density regime impossible for the same parameters for which it was possible in the infinite system. The domain wall starts to propagate, leading to the expansion of the high-density region, erasing progressively the solitons in the channel.

Figure 8 plots the velocity of the ”soliton head” (blue line) as a function of the Support in the 24 $\mu$m channel for the same value of the pump which was used in the maze solving regime (see Fig. 3). One sees a behavior very similar to the ”pure” domain wall from Fig. 1 (replotted as a black line), but shifted to much smaller values of $S$. This shift is due to the confinement of the domain wall in the transverse direction (vertical) and to the extra flows from the walls of the channel (shown in Fig. 4(c) of the main text).

To elucidate the contribution of these two effects, we plot an additional phase diagram (Fig. 9), showing for a single channel having a dead end the boundary between the two cases (solitons or high density in the channel) versus the channel width $L$ and the support amplitude $S$. The boundary between the two cases in shown in Fig. 4(d) of the main text. It is well described by an inverse square function with a single-parameter fit:

[TABLE]

giving $L_{0}=13.9\pm 0.1$ $\mu$m. This is the narrowest channel width where a pair of solitons can exist without support. The inverse square dependence is typical for a quantity depending on the confinement energy, and the value of $L\approx 14$ $\mu$m corresponds to the minimal width of a finite quantum well (with its depth of $0.14$ meV given by the laser detuning) where the two-node state can be confined for such particle mass (accounting for the finite broadening due to the lifetime). For a wider channel, the kinetic energy of the two-node state is lower, and the support and the flows shift it upwards via the interactions. The inverse square law holds well up to $L=30$ $\mu$m, which is larger than the channel width used in the main text ($25$ $\mu$m). For even larger channels, the dependence becomes more complicated, because the low-density solution can contain more than 2 solitons, depending on $L$. The detailed study of this dependence is beyond the scope of the present work.

One could imagine that the motion of the soliton heads should lead to the switching of the whole maze structure to the high-density regime. However, this does not occur, because when the domain wall associated with a soliton edge disappearing from a dead end arrives to a T-junction where there is another soliton, the situation becomes fully equivalent to Fig. 2 of the main text: the T-junction with one channel in the high-density regime is actually not a T-junction any more, but just a channel with homogeneous walls (red curve in Fig. 4(d)). We conclude that the solitons can only disappear from the dead ends, whereas solitons connecting the exits of the maze persist, thus allowing to solve the maze.

In the main text, we also plot the threshold support intensity versus the channel width for an open-end system, which exhibits higher values than for the dead end case as previously explained. Interestingly, our simulations show that in a large open-end channel the threshold support can even be larger than $S_{c}$. Indeed, in such case the lateral walls start to propagate towards the center where a 4-soliton structure forms. This structure is quite stable and produces a quantum pressure, which counteracts the pressure of the interactions and blocks the motion of the walls even if $S$ is slightly larger than $S_{c}$. We remind that the fact that the threshold for the open end channels is higher than for a dead end channel, is what allows the head domain wall to withdraw from from the dead ends, but to stop at the crossroads, and therefore the maze solving process.

V.1 Effects of disorder on the maze solving property

As discussed in the main text, in presence of a large disorder the detuning has to be increased. Fig. 10 is computed for a realistic disorder of 0.1 meV. To avoid disorder disturbance, we modified all parameters (detuning, lifetime, pump and support intensity), in order for the kinetic and interaction energies to be at least 3 times larger than the disorder amplitude. One should note that the used parameters are perfectly accessible experimentally. As one can see, the maze remains well defined and is solved, despite the disorder value.

VI Supplementary Movies

The supplementary video file video.avi shows an example of the solution of a large maze with time. The image size is $1024~{}\mu$m$\times 1024~{}\mu$m. The maze solving time is 16 ns. We see that the vortex streets withdraw from the dead ends and only the vortex street connecting the two exits of the maze remains.

An additional movie video2.avi shows the solution of a small maze (similar to that of Fig. 4 of the main text) over time, but with a large disorder amplitude, as shown in Fig. 10.