---
abstract: 'We have used an extension of our slow light technique to provide a method for inducing small density defects in a Bose-Einstein condensate. These sub-resolution, micron-sized defects evolve into large amplitude sound waves. We present an experimental observation and theoretical investigation of the resulting breakdown of superfluidity. We observe directly the decay of the narrow density defects into solitons, the onset of the ‘snake’ instability, and the subsequent nucleation of vortices.'
author:
- |
    \
    \
    \
    \
    \
    \
    \
    \
    \
    \
    \
    \
    \
title: 'Observation of Quantum Shock Waves Created with Ultra Compressed Slow Light Pulses in a Bose-Einstein Condensate'
---

Superfluidity in Bose condensed systems represents conditions where frictionless flow occurs because it is energetically impossible to create excitations. When these conditions are not satisfied, various excitations develop, and experiments on superfluid helium, for example, have provided evidence that the nucleation of vortex rings occurs when ions move through the fluid faster than a critical speed ([*1,2*]{}). Under similar conditions, shock waves would occur in a normal fluid ([*3*]{}). Such discontinuities are not allowed in a superfluid and instead topological defects, such as quantized vortices and solitons, are nucleated when the spatial scale of density variations becomes on the order of the healing length. This is the length scale at which the kinetic energy, associated with spatial variations of the macroscopic condensate wave function, becomes on the order of the atom-atom interaction energy ([*2,4*]{}). It is therefore also the minimum length over which the density of a condensate can adjust.

Bose-Einstein condensates (BECs) of alkali atoms ([*5*]{}) have provided a system for the study of superfluidity, which is theoretically more tractable than liquid helium and allows greater experimental control. An exciting recent development is the production of solitons and vortices. Experiments have employed techniques that manipulated the phase of the BEC ([*6-9*]{}), or provided the system with a high angular momentum which makes vortex formation energetically favorable ([*10,11*]{}). However, a direct observation of the formation of vortices via the breakdown of superfluidity has been lacking. Rather, rapid heating from ’stirring’ with a focused laser beam above a critical velocity was observed as indirect evidence of this process ([*12,13*]{}).

In this context, it is natural to ask what would happen if one were to impose a sharp density feature in a BEC with a spatial scale on the order of the healing length. Optical resolution limits have prevented direct creation of this kind of excitation. In this paper we present an experimental demonstration of creation of such defects in sodium Bose-Einstein condensates, using a novel extension of the method of ultra slow light pulse propagation ([*14*]{}) via electromagnetically induced transparency (EIT) ([*15,16*]{}).

Our slow-light setup is described in ([*14*]{}). By illuminating a BEC with a ‘coupling’ laser, we create transparency and low group velocity for a ‘probe’ laser pulse subsequently sent into the cloud. In a geometry where the coupling and probe laser beams propagate at right angles, we can control the propagation of the probe pulse from the side. By introducing a spatial variation of the coupling field, along the pulse propagation direction, we vary the group velocity of the probe pulse across the cloud. Here we accomplish this by blocking half of the coupling beam so it illuminates only the left hand side ($z<0$) of the condensate, setting up a light ‘roadblock’. In this way, we compress the probe pulse to a small spatial region at the center of the BEC, while bypassing the usual bandwidth requirements for slow light ([*17*]{}). The technique produces a short wavelength excitation by suddenly removing a narrow disk of the condensate, with the width of the disk determined by the width of the compressed probe pulse.

We find that this excitation results in short wavelength, large amplitude sound waves that shed off ‘gray’ solitons ([*18-20*]{}), and we make the connection to the formation of shock waves in classical fluids. The ‘snake’ (Kadomtsev-Petviashvili) instability is predicted to cause solitons to decay into vortices ([*21-24*]{}). This has been observed with optical solitons ([*25*]{}) and recently the JILA group predicted and observed that solitons in a BEC decay into vortex rings ([*9*]{}). Here we present a direct observation of the dynamics of the snake instability in a BEC and the subsequent nucleation of vortices. The images of the evolution are compared to numerical propagation of the Gross-Pitaevskii equation in two dimensions.

Details of our Bose-Einstein condensation apparatus can be found in ([*26*]{}). We use condensates with about 1.5 million sodium atoms in the state $|1 \rangle \equiv |3S_{1/2},F=1,M_F=-1
\rangle$ and trapped in a 4-Dee magnet. For the experiments presented here, the magnetic trap has an oscillator frequency $\omega_z=(2\pi) 21$ Hz along its symmetry direction and a frequency $\omega_x = \omega_y = 3.8 \omega_z$ in the transverse directions. The peak density of the condensates is then $9.1
\times 10^{13}~\mathrm{cm}^{-3}$. The temperature is $T \sim
0.5\,T_c$, where $T_c=300$ nK is the critical temperature for condensation, and so the vast majority ($\sim 90 \%$) of the atoms occupy the ground state.

To create slow and spatially localized light pulses, the coupling beam propagates along the $x$ axis (Fig. 2B), is resonant with the $D_1$ transition from the unoccupied hyperfine state $|2 \rangle
\equiv |3S_{1/2},F=2,M_F=-2 \rangle$ to the excited level $ |3
\rangle \equiv |3P_{1/2}, F=2, M_F=-2 \rangle$, and has a Rabi frequency $\Omega_c = (2\pi) 15$ MHz ([*27*]{}). We inject probe pulses along the $z$ axis, resonant with the transition and with peak Rabi frequency $\Omega_p = (2\pi) 2.5$ MHz. The pulses are Gaussian shaped with a $1/e$ half-width of $\tau = 1.3\;\mu$s. With the entire BEC illuminated by the coupling beam, we observe probe pulse delays of $4\;\mu$s for propagation through the condensates, corresponding to group velocities of 18 m/s at the center of the clouds. A pulse with a temporal half-width $\tau$ is spatially compressed from a length $2 \, c \tau$ in vacuum to ([*14,17,28,29*]{}) $$L=2 \tau V_g = 2 \tau \frac{|\Omega_c|^2}{\Gamma f_{13} \sigma_0
n_c}$$

inside the cloud, where $\Gamma = (2 \pi) 10$ MHz is the decay rate of state $|3 \rangle$, $n_c$ is the cloud density, $\sigma_0 = 1.65 \times 10^{-9}\;\mathrm{cm}^2$ is the absorption cross-section for light resonant with a two-level atom, and $f_{13}=1/2$ is the oscillator strength of the $|1 \rangle
\rightarrow |3 \rangle$ transition. The atoms are constantly being driven by the light fields into a dark state, a coherent superposition of the two hyperfine states $|1 \rangle$ and $ |2
\rangle$ ([*15*]{}). In the dark state, the ratio of the two population amplitudes is varying in space and time with the electric field amplitude of the probe pulse as

$$\psi_2 = - \frac{\Omega_p}{\Omega_c}\,\psi_1,$$

where $\psi_1,\;\psi_2$ are the macroscopic condensate wave functions associated with the two states $|1 \rangle$ and $
|2 \rangle$.

For the parameters listed above, the probe pulse is spatially compressed from 0.8 km in free space to only 50 $\mu$m at the center of the cloud, at which point it is completely contained within the atomic medium. The corresponding peak density of atoms in $|2 \rangle$, proportional to $|\psi_2|^2$, is $1/34$ of the total atom density. The $|1 \rangle$ atoms have a corresponding density depression.

From Eqs. 1 and 2, it is clear that to minimize the spatial scale of the density defect, we need to use short pulse widths and low coupling intensities. However, for all the frequency components of the probe pulse to be contained within the transmission window for propagation through the BEC ([*17*]{}), we need a pulse with a temporal width $\tau$ of at least $2 \sqrt{D} \Gamma /
{\Omega_c}^2 \approx 0.3~\mu$s (here $D \approx 520$ is the optical density of a condensate for on-resonance two level transitions) to avoid severe attenuation and distortion. Furthermore, we see from Eq. 2 that to maximize the amplitude of the density depression would favor use of a peak Rabi frequency for the probe of $\Omega_p \sim \Omega_c$. This also severely distorts the pulse.

Both of these distortion effects accumulate as the pulse propagates through large optical densities. This motivated us to introduce a roadblock in the condensate for a light pulse approaching from the left hand ($z<0$) side. By imaging a razor blade onto the right half of the condensate, we ramp the coupling beam from full to zero intensity over the course of a $12\;\mu$m region in the middle of the condensate, determined by the optical resolution of the imaging system. In the illuminated region ($z<0$), our bandwidth and weak-probe requirements are well satisfied and we get undistorted, unattenuated propagation through the first half of the cloud to the high-density, central region of the condensate. As the pulse enters the roadblock region of low coupling intensity, it is slowed down and spatially compressed. The exact shape and size of the defects which are created with this method are dependent on when absorption effects become important.

To accurately model the pulse compression and defect formation, we account for the dynamics of the slow light pulses, the coupling field, and the atoms self-consistently. At sufficiently low temperatures, the dynamics of the two-component condensate can be modelled with coupled Gross-Pitaevskii (GP) equations ([*4,5*]{}). Here we include terms to account for the resonant two-photon light coupling between the two components:

$$\begin{aligned}
i \hbar \frac{\partial}{\partial t}\psi_1 & = &
\left(-\frac{\hbar^2 \nabla^2}{2m} + V_1(\mathbf{r}) + U_{11}
|\psi_1|^2 + U_{12}
|\psi_2|^2 \right) \psi_1 \nonumber \\
& & \hspace{2 cm} - i \frac{\mid{\Omega_p}\mid^2}{2 \Gamma} \psi_1
- i \frac{{\Omega_p}^\ast \Omega_c}{2 \Gamma} \psi_2 \nonumber
- i N_c \sigma_e \hbar \frac{k_{2 \gamma}}{2 m} |\psi_2|^2 \psi_1,  \\
i \hbar \frac{\partial}{\partial t}\psi_2 & = &
\left(-\frac{\hbar^2 (\nabla^2 + i \mathbf{k_{2 \gamma}} \cdot
\nabla )}{2m} + V_2(\mathbf{r}) + U_{22} |\psi_2|^2 + U_{12}
|\psi_1|^2
\right) \psi_2 \nonumber \\
& & \hspace{2 cm}- i \frac{\mid{\Omega_c}\mid^2}{2 \Gamma} \psi_2
- i \frac{\Omega_p {\Omega_c}^\ast}{2 \Gamma} \psi_1 - i N_c
\sigma_e \hbar \frac{k_{2 \gamma}}{2 m} |\psi_1|^2 \psi_2.\end{aligned}$$

Here $V_1(\mathbf{r}) = \frac{1}{2} m {\omega_z}^2(
\lambda^2 (x^2+y^2) + z^2)$, where $m$ is the mass of the sodium atoms, and $\lambda=3.8$. Due to the magnetic moment of atoms in state $|2 \rangle$, $V_2(\mathbf{r})= - 2 V_1(\mathbf{r})$, and atoms in this state are repelled from the trap. The EIT process involves absorption of a probe photon and stimulated emission of a coupling photon, leading to a $4.1$ cm/s recoil velocity. This is described by a term in the second equation, containing $\mathbf{k_{2 \gamma}} = \mathbf{k_p - k_c}$, the difference in wave vectors between the two laser beams. (Here we use a gauge where the recoil momentum is transformed away.) Atom-atom interactions are characterized by the scattering lengths, $a_{ij}$, via $U_{ij} = 4 \pi N_c \hbar^2 a_{ij}/m$, where $a_{11}=2.75\;\mathrm{nm}$, $a_{12}=a_{22} = 1.20\,a_{11}$ ([*30*]{}), and $N_c$ is the total number of condensate atoms. To obtain the light coupling terms in Eq. 3, we have adiabatically eliminated the excited state amplitude $\psi_3$ ([*31*]{}), as the relaxation from spontaneous emission occurs much faster than the light coupling and external atomic dynamics driving $\psi_3$. In our model, atoms in $|3 \rangle$ that spontaneously emit are assumed to be lost from the condensate, which is why the light coupling terms are non-Hermitian. Finally, the last term in each equation accounts for losses due to elastic collisions between high momentum $|2\rangle$ atoms and the nearly stationary $|1\rangle$ atoms ($\sigma_e = 8 \pi {a_{12}}^2$) ([*32*]{}).

The spatial dynamics of the light fields are described classically with Maxwell’s equations in a slowly varying envelope approximation, again using adiabatic elimination of $\psi_3$ :

$$\begin{aligned}
\frac{\partial}{\partial z}\Omega_p & = & - \frac{1}{2}f_{13}
\sigma_0 N_c(\Omega_p |\psi_1|^2 +
\Omega_c {\psi_1}^\ast \psi_2), \nonumber \\
\frac{\partial}{\partial x}\Omega_c & = & - \frac{1}{2}f_{23}
\sigma_0 N_c(\Omega_c |\psi_2 \mid^2 + \Omega_p \psi_1
{\psi_2}^\ast).\end{aligned}$$

In the region where the coupling beam is illuminating the BEC ($z<0$), the light coupling terms dominate the atomic dynamics and solving Eqs. 3 and 4 reduces to Eqs. 1 and 2 above.

We have performed numerical simulations in two dimensions ($x$ and $z$) to track the behavior of the light fields and the atoms. The probe and coupling fields were propagated according to Eq. 4 with a second order Runge-Kutta algorithm ([*33*]{}) and the atomic mean fields were propagated according to Eq. 3, with an Alternating-Direction Implicit variation of the Crank-Nicolson algorithm ([*33,34*]{}). In this way, Eqs. 3 and 4 were solved self-consistently ([*35*]{}). Profiles of the probe pulse intensity along $z$, through $x=0$, are shown in Fig. 1A. As the pulse runs into the roadblock, a dramatic compression of the probe pulse’s spatial length occurs. When the probe pulse enters the low coupling region, the Rabi frequency $|\Omega_p|$ becomes on the order of $|\Omega_c|$. So the density of state $|2 \rangle$ atoms, $N_c |\psi_2|^2$, increases in a narrow region, which is accompanied by a decrease in $N_c |\psi_1|^2$ (Fig. 1B). The half-width of the defect is $2\;\mu$m. As the compression develops, absorption/spontaneous emission events eventually start to remove atoms from the condensate and reduce the probe intensity.

Experimental results are shown in Fig. 2. Fig. 2A is an in-trap image of the original condensate of $|1 \rangle$ atoms, Fig. 2B diagrams the beam geometry, and Fig. 2C shows a series of images of state $|2 \rangle$ atoms as the pulse propagates into the roadblock. The corresponding optical density (OD) profiles along $z$ through $x=0$ are also shown. The OD is defined to be $-
\mathrm{ln}(I/I_0)$, where $(I/I_0)$ is the transmission coefficient. All imaging is done with near resonant laser beams propagating along the $y$-axis, and with a duration of 10 $\mu$s. There is clearly a build-up of a dense, narrow sample of $|2
\rangle$ atoms at the center of the BEC as the pulse propagates to the right. Note that the pulse reaches the roadblock at the top and bottom edges of the cloud before the roadblock is reached at the center, which is a consequence of the transverse variation in the density of the BEC, with the largest density along the center line. After the pulse compression is achieved, we shut off the coupling beam to avoid heating and phase shifts of the atom cloud due to extended exposure to the coupling laser, and the subsequent dynamics of the condensate are observed. (We observed that exposure to the coupling laser alone, for the exposure times used to create defects, causes no excitations of the condensates).

In considering the dynamics resulting from this excitation of a condensate, it should be noted that the roadblock ‘instantaneously’ removes a spatially selected part of $\psi_1$. The entire light compression happens in approximately $15\;\mu\mathrm{s}$. After the pulse is stopped and the coupling laser turned off, the $|2 \rangle$ atoms remaining in the condensate $\psi_2$ have a $4.1$ cm/s recoil and atoms which have undergone absorption and spontaneous emission events have a similarly sized but randomly directed recoil. So the $\psi_2$ component and the other recoiling atoms interact with $\psi_1$ for less than $0.5$ ms before leaving the region. Both of these time scales are short compared with the several millisecond timescale over which most of the subsequent dynamics of $\psi_1$ occur, as discussed in the following.

We first considered the one-dimensional dynamics along the $z$ axis. Snapshots of both condensate components, obtained from numerical propagation in 1D according to Eqs. 3 and 4, are shown at various times after the pulse is stopped at the roadblock (and the coupling laser turned off) (Fig. 3A). In the linear hydrodynamic regime, where the density defect has a relative amplitude $A \ll 1$ and a half-width $\delta \gg \xi$ (here $\xi=
1/\sqrt{8 \pi N_c |\psi_{1}|^2
 a_{11}}$ is the local healing length which is $0.4\;\mu\mathrm{m}$ at the center of the ground state condensate in our experiment), one expects to see two density waves propagating in opposite directions at the local sound speed, $c_s = \sqrt{U_{11} |\psi_1
\mid^2/m}$, as seen experimentally in ([*36*]{}). However, for the parameters used in our experiment, the sound waves are seen to shed off sharp features propagating at lower velocities. Examination of the width, speed, and the phase jump across these features shows that they are gray solitons. Describing the slowly varying background wave function of the condensate with $\psi^{(0)}_{1}\!,$ the wave function in the vicinity of a gray soliton centered at $z_0$ is ([*18-20*]{}):

$$\psi_1(z,t) = \psi^{(0)}_{1}(z,t) \left( i \sqrt{1-\beta^2} +
\beta \; \mathrm{tanh}\left(\frac{\beta}{\sqrt{2}~\xi}
(z-z_0)\right) \right).$$

The dimensionless constant $\beta$ characterizes the ‘grayness’, with $\beta = 1$ corresponding to a stationary soliton with a $100\%$ density depletion. With $\beta \ne 1$, the solitons travel at a fraction of the local sound velocity, $c_s
\sqrt{1-\beta^2}$. As seen in the figure, after a shedding event, the remaining part of the sound wave continues to propagate at a reduced amplitude. Our numerical simulations show that the solitons eventually reach a point where their central density is zero and then oscillate back to the other side, in agreement with the discussions in ([*18,20*]{}).

In Fig. 3A, each of the two sound waves shed off two solitons. By considering the available free energy created by a defect, one finds that, when the defect size is somewhat larger than the healing length and the defect amplitude, $A$, is on the order of unity, the number of solitons that can be created is approximately $\sqrt{A} \delta /(2 \xi)$.

One obtains a simple physical estimate of the conditions necessary for soliton shedding by calculating the difference in sound speed associated with the difference in atom density between the center and back edge of the sound wave. As confirmed by our numerical simulations, this difference leads to development of a steeper back edge and an increasingly sharp jump in the phase of the wave function. This is the analog of shock wave formation from large amplitude sound waves in a classical fluid ([*3*]{}). When the spatial width of the back edge has decreased to the width of a soliton with amplitude $\beta = \sqrt{A/2}$ (according to Eq. 5), such a soliton is shed off the back. It’s subsonic speed causes it to separate from the remaining sound wave. Furthermore, by creating defects with sizes on the order of the healing length, we excite collective modes of the condensate, with wave vectors on the order of the inverse healing length. In this regime, the Bogoliubov dispersion relation is not linear ([*4,5*]{}), and accordingly some of the sound wave will disperse into smaller ripples, as seen in Fig. 3A.

Considering the evolution of a defect of relative density amplitude $A$ and half-width $\delta$ in an otherwise homogenous medium, we estimate that solitons of amplitude $\beta =
\sqrt{A/2}$ will be created after the two resulting sound waves have propagated a distance

$$z_{sol} = \frac{2 \delta}{A}\left( \frac{1 - \frac{\xi}{2
\delta}}{1 - \frac{{\pi}^2 \xi^2}{\delta^2 A}}\right).$$

This is in agreement with our numerical calculations. We conclude that the minimum soliton formation length is obtained for large amplitude defects with a width just a few times the healing length. This dictates the defect width picked in the experiments presented here. Narrower defects disperse, whereas larger defects, comparable to the cloud size, couple to collective, nonlocalized excitations of the condensate.

We explored the soliton formation experimentally by creating defects in a BEC with the light roadblock. We controlled the size of the defects by varying the intensity of the probe pulses, which had a width $\tau = 1.3~\mu$s. OD images of state $|1 \rangle$ condensates are shown (Fig. 4) in one particular case. Immediately after the defect is created, the trap is turned off, and the cloud evolves and expands for 1 ms and 10 ms, respectively. As seen from the 1 ms picture, a single deep defect is formed initially, which results in creation of 5 solitons after 10 ms of condensate dynamics and expansion. The initial defect created in the trap could not be resolved with our imaging system, which has a resolution of 5 $\mu$m. By varying the probe intensity, we find that the number of solitons formed scales linearly with the probe pulse energy, as expected.

To study the stability of the solitons, we first performed 2D numerical simulations of Eqs. 3 and 4. Fig. 3B shows density profiles, $N_c |\psi_1|^2$, obtained for the same parameters as used in Fig. 1B. Again, the profiles are shown at various times after the pulse is compressed and stopped. The deepest soliton (the one closest to the center) is observed to quickly curl and eventually collapse into a vortex pair. The wave function develops a $2 \pi$ phase shift in a small circle around the vortex cores, which shows that they are singly quantized vortices. Also, the core radius is approximately the healing length. (Upon collapse, a small sound wave between the two vortices carries away some of the remaining soliton energy.) This decay can be understood as resulting from variation in propagation speed along the transverse soliton front. As discussed in ([*22-25*]{}), a small deviation will be enhanced by the nonlinearity in the Gross-Pitaevskii equation, and thus, the soliton collapses about the deepest (and therefore slowest) point.

Fig. 5 shows experimental images of state $|1 \rangle$ condensates. After the defect is created, the condensate of $|1
\rangle$ atoms is left in the trap for a varying amount of time (as indicated on the figures). The trap was then abruptly turned off and, $15$ ms after release, we imaged a selected slice of the expanded condensate, with a thickness of 30 $\mu$m along $y$ ([*37*]{}). The release time of $15$ ms is picked large enough that the condensate structures are resolvable with our imaging system ([*38*]{}). The slice was optically pumped from state $|1 \rangle$ to the $|3S_{1/2},\,F=2 \rangle$ manifold for $10\;\mu\mathrm{s}$ before it was imaged with absorption imaging by a laser beam nearly resonant with the transition from $|3S_{1/2},\,F~=~2
\rangle$ to $|3P_{3/2},\,F~=~3 \rangle$. The total pump and imaging time was small enough that no significant motion due to photon recoils occurred during the exposure. The slice was selected at the center of the condensate by placing a slit in the path of the pump beam.

For the data in Fig. 5A, it is seen that the deepest soliton curls as it leaves certain sections behind, and at 1.2 ms it has nucleated vortices. This is a direct observation of the snake instability. In Fig. 5B, at 0.5 ms, the snake instability has caused a complicated curving structure in one of the solitons, and vortices are observed after 2.5 ms. The vortices are seen to persist for many milliseconds and slowly drift towards the condensate edge. We observed them even after $30$ ms of trap dynamics, long enough to study the interaction of vortices with sound waves reflected off the condensate boundaries. Preliminary results, obtained by varying the $y$ position of the imaged condensate slices, indicate a complicated 3D structure of the vortices. In addition, the defect has induced a collective motion of the condensate whereby atoms, originally in the sides of the condensate, attempt to fill in the defect. This leads to a narrow and dense central region, which then slowly relaxes (Fig. 5B).

We performed the experiment with a variety of Rabi frequencies for the probe pulses, and saw nucleation of vortices only for the peak $\Omega_p > (2 \pi) 1.4$ MHz. The free energy of a vortex is substantially smaller near the border of the condensate where the density is smaller, so smaller (and thus lower energy) defects will form vortices very near the condensate edges, seen as ‘notches’ in Figs. 3B and 5.

In conclusion, we have studied and observed how small wavelength excitations cause a breakdown of superfluidity in a BEC. Our results show how localized defects in a superfluid will quite generally either disperse into high frequency ripples or end up in the form of topological defects such as solitons and vortices, and we have obtained an analytic expression for the transition between the two regimes. By varying our experimental parameters, we can create differently sized and shaped defects, and also control the number of defects created, allowing studies of a myriad of effects. Among them are soliton-soliton collisions, more extensive studies of soliton stability, soliton-sound wave collisions, vortex-soliton interactions, vortex dynamics, interaction between vortices, and the interaction between the BEC collective motion and vortices.

References {#references .unnumbered}
==========

1.  R.J. Donnelly, [*Quantized vortices in Helium II*]{} (Cambridge Univ. Press, Cambridge, 1991).

2.  R.M. Bowley, J. Low Temp. Phys. [**87**]{}, 137 (1992).

3.  L.D. Landau and E.M. Lifshitz, [*Fluid Mechanics*]{} (Pergamon Press, New York, 1959).

4.  A.L. Fetter, in [*Bose-Einstein Condensation in Atomic Gases, Proceedings of the International School of Physics Enrico Fermi, Course CXL*]{}, M. Inguscio, S. Stringari, C. Wieman, Eds. (International Organisations Services B.V., Amsterdam, 1999).

5.  For a review, see F. Dalfovo, S. Giorgini, L.P. Pitaevskii, S. Stringari, Rev. Mod. Phys. [**71**]{}, 463 (1999).

6.  M.R. Matthews [*et al.*]{}, Phys. Rev. Lett. [**83**]{}, 2498 (1999).

7.  S. Burger, K. Bongs, S. Dettmer, W. Ertmer, K. Sengstock, Phys. Rev. Lett. [**83**]{}, 5198 (1999).

8.  J. Denschlag [*et al.*]{}, Science [**287**]{}, 97 (2000).

9.  B.P. Anderson [*et al.*]{}, Phys. Rev. Lett. [**86**]{}, 2926 (2001).

10. K.W. Madison, F. Chevy, W. Wohlleben, J.Dalibard, Phys. Rev. Lett. [**84**]{}, 806 (2000).

11. J.R. Abo-Shaeer, C. Raman, J.M. Vogels, W. Ketterle, Science [**292**]{}, 476 (2001).

12. C. Raman [*et al.*]{}, Phys. Rev. Lett. [**83**]{}, 2502 (1999).

13. B. Jackson, J.F. McCann, C.S. Adams, Phys. Rev. A [**61**]{}, 051603(R) (2000).

14. L.V. Hau, S.E. Harris, Z. Dutton, C.H. Behroozi, Nature [**397**]{}, 594 (1999).

15. S.E. Harris, Physics Today [**50**]{}, 36 (1997).

16. M.O. Scully and M.S. Zubairy, [*Quantum Optics*]{} (Cambridge Univ. Press, 1997).

17. S.E. Harris and L.V. Hau, Phys. Rev. Lett. [**82**]{}, 4611 (1999).

18. S.A. Morgan, R.J. Ballagh, K. Burnett, Phys. Rev. A [**55**]{}, 4338 (1997).

19. W.P. Reinhardt and C.W. Clark, J. Phys. B [**30**]{}, L785 (1997).

20. Th. Busch and J.R. Anglin, Phys. Rev. Lett. [**84**]{}, 2298 (2000).

21. B.B. Kadomtsev and V.I. Petviashvili, Sov. Phys. Dokl. [**15**]{}, 539 (1970).

22. C.A. Jones, S.J. Putterman, P.H. Roberts, J. Phys. A [**19**]{}, 2991 (1986).

23. C. Josserand and Y. Pomeau, Europhys. Lett. [**30**]{}, 43 (1995).

24. D.L. Feder, M.S. Pindzola, L.A. Collins, B.I. Schneider, C.W. Clark, Phys. Rev. A [**62**]{}, 053606 (2000).

25. A.V. Mamaev, M. Saffman, A.A. Zozulya, Phys. Rev. Lett. [**76**]{}, 2262 (1996).

26. L.V. Hau [*et al.*]{}, Phys. Rev. A [**58**]{}, R54 (1998).

27. We use the definition $\Omega_p
    = \mathbf{E}_{p0} \cdot \mathbf{d}_{13}, \Omega_c =
    \mathbf{E}_{c0} \cdot \mathbf{d}_{23}$, where $\mathbf{E}_{p0},\mathbf{E}_{c0}$ are the slowly varying electric field amplitudes, and $\mathbf{d}_{13},\mathbf{d}_{23}$ are the electric dipole matrix elements of the transitions.

28. S.E. Harris, J.E. Field, A. Kasapi, Phys. Rev. A [**46**]{}, R29 (1992).

29. C. Liu, Z. Dutton, C.H. Behroozi, L.V. Hau, Nature [**409**]{}, 490 (2001).

30. J.P. Burke, C.H. Greene, J.L Bohn, Phys. Rev. Lett. [**81**]{}, 3355 (1998). A Mathematica notebook is available at http://fermion.colorado.edu/$\sim$chg/Collisions/.

31. J. Javanainen and J. Ruostekoski, Phys. Rev. A [**52**]{} 3033 (1995).

32. Y.B. Band, M.Tribbenbach, J.P. Burke, P.S. Julienne, Phys. Rev. Lett. [**84**]{}, 5462 (2000).

33. W. H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, [*Numerical Recipes in C, Second Edition*]{} (Cambridge University Press, Cambridge, 1992).

34. S.E. Koonin and D.C. Merideth, [*Computational Physics*]{}, (Addison-Wesley, Reading, MA, 1990).

35. For propagation of the GP equation in 1D, we typically used a spatial grid with 4000 points and $dz=0.040\;\mu$m. In 2D simulations, we typically used a $750 \times 750$ grid with $dz=0.21\;\mu$m and $dx=0.057\;\mu$m. To solve the equations self-consistently, we kept track of the wave functions at previous time points and projected forward to second order. Smaller time steps and grid spacing were also used to assure convergence of the results. To mimic the nonlinear interaction strength at the center of a [*three*]{} dimensional cloud we put in an effective condensate radius (calculated with the Thomas-Fermi approximation ([*39*]{})) in the dimensions which were not treated dynamically. In all calculations, the initial condition was the ground state condensate wave function with all atoms in $|1 \rangle$, obtained by propagating a Thomas-Fermi wave function in imaginary time.

36. M.R. Andrews, [*et al.*]{}, Phys. Rev. Lett. [**79**]{}, 553 (1997).

37. M.R. Andrews [*et al.*]{}, Science [**275**]{}, 637 (1997).

38. F. Dalfavo and M. Modugno, Phys. Rev. A [**61**]{}, 023605 (2000).

39. G. Baym and C. Pethick, Phys. Rev. Lett. [**76**]{}, 6 (1996).

40. Supplementary material is available at [*Science*]{} Online at www.sciencemag.org.

41. This work was supported by the Rowland Institute for Science, the Defense Advanced Research Projects Agency, the U.S. Airforce Office of Scientific Research, the U.S. Army Research Office OSD Multidisciplinary University Research Initiative Program, the Harvard Materials Research Science and Engineering Center sponsored by the National Science Foundation, and by the Carlsberg Foundation, Denmark. C.S. is supported by a National Defense Science and Engineering Grant sponsored by the Department of Defense.

Figure Captions {#figure-captions .unnumbered}
===============

**Fig. 1.** **(A)** Compression of a probe pulse at the light ’roadblock’, according to 2D numerical simulations of Eqs. 3 and 4. The solid curves indicate probe intensity profiles along $z$ (at $x=0$), normalized to the peak input intensity. The snapshots are taken at a sequence of times, indicated in the figure, where $t=0$ is defined as the time when the center of the probe pulse enters the BEC from the left. For reference, the atomic density profile of the original condensate is plotted (in arbitrary units) as a dashed curve. The gray shading indicates the relative coupling input intensity as a function of $z$, with white corresponding to full intensity, and the darkest shade of gray corresponding to zero. The spatial turn off of the coupling field is centered at $z=0$ and occurs over $12\;\mu$m, as in the experiment. The number of condensate atoms is $1.2 \times 10^6$ atoms, the peak density is $6.9 \times
10^{13}\; \mathrm{cm}^{-3}$, and the coupling Rabi frequency is $\Omega_c=(2 \pi) 8.0$ MHz. The probe pulse has a peak Rabi frequency of $\Omega_p=(2 \pi)2.5$ MHz and a $1/e$ half-width of $\tau = 1.3 \;\mu$s. **(B)** Creation of a narrow density defect in a BEC. Density profiles of the two condensate components, $N_c|\psi_1|^2$ (dashed) and $N_c|\psi_2|^2$ (solid), along $z$ at $x=0$ for a sequence of times. Note that the $z$ range of the plot is restricted to a narrow region around the roadblock at the cloud center. The densities are normalized relative to the peak density of the original condensate indicated by the red dashed curve. The other curves correspond to times $1
\;\mu$s (green), $4 \;\mu$s (blue), and $14\; \mu$s (black). The width of the probe pulse is $\tau = 4 \;\mu$s and the other parameters are the same as in (A). An animated version is provided in the supplemental material ([*40*]{}) as Animation 1.\
\
**Fig. 2.** **(A)** Experimental in-trap OD image of a typical BEC before illumination by the probe pulse and coupling field. The condensate contains $1.5 \times 10^6$ atoms. The imaging beam was -30 MHz detuned from the $|1 \rangle
\rightarrow |3P_{1/2}, F=2, M_F = -1 \rangle$ transition. **(B)** Top view of the beam configuration used to create and study localized defects in a BEC as discussed in the text. **(C)** Build up of state $|2 \rangle$ atoms at the road block. In-trap OD images (left) show the transfer of atoms from $|1\rangle$ to $|2 \rangle$ as the probe pulse propagates through the condensate and runs into the roadblock. The atoms in $|2\rangle$ were imaged with a laser beam -13 MHz detuned from the $| 2\rangle \rightarrow |3P_{3/2}, F=3, M_F = -2\rangle$ transition. To allow imaging, the probe pulse propagation was stopped at various times, indicated in the figure, by switching the coupling beam off ([*29*]{}). The figures on the right show the corresponding line cuts along the probe propagation direction through the center of the BEC. The probe pulse had a Rabi frequency $\Omega_p=(2 \pi) 2.4$ MHz and the coupling Rabi frequency was $\Omega_c= (2 \pi) 14.6$ MHz.\
\
**Fig. 3.** **(A)** Formation of solitons from a density defect. The plots show results of a 1D numerical simulation of Eqs. 3 and 4. The light roadblock forms a defect and the subsequent formation of solitons is seen. The defect is set up with the same parameters for the light fields as in Fig. 1B. The number of condensate atoms is $N_c=1.2 \times 10^6$ and the peak density is $7.5 \times 10^{13}\;\mathrm{cm}^{-3}$. The times in the plots indicate the evolution time relative to the time when the probe pulse stops at the roadblock and the coupling beam is switched off (at $t=8\;\mu\mathrm{s}$ with $t=0$ as defined in Fig. 1A). The solid and dashed curves show the densities of $|1\rangle$ atoms ($N_c | \psi_1|^2$) and $|2\rangle$ atoms ($N_c |\psi_2|^2$). The phase of $\psi_1$ is shown in each case with a dotted curve (with an arbitrary constant added for graphical clarity). In the first two frames, the $|2\rangle$ atoms quickly leave due to the momentum recoil, leaving a large-amplitude, narrow defect in $\psi_1$ (Because this is a 1D simulation, the momentum kick in the $x$ direction is ignored). **(B)** The snake instability and the nucleation of vortices. The plots show the density $N_c |\psi_1|^2$ from a numerical simulation in 2D, with white corresponding to zero density and black to the peak density ($6.9 \times
10^{13}\;\mathrm{cm}^{-3}$). The parameters are the same as in Fig. 1B and the times indicated are relative to the coupling beam turn-off at $t=21\;\mu\mathrm{s}$. The solitons curl about their deepest point, eventually breaking and forming vortex pairs of opposite circulation (seen first at 3.5 ms). Several vortices are formed and the last frame shows the vortices slowly moving towards the edge of the condensate. At later times, they interact with sound waves which have reflected off the condensate boundaries. Animated versions (Animations 2 and 3) are provided in the supplemental material ([*40*]{}).\
\
**Fig. 4.** Experimental OD images and line cuts (at $x=0$) of a localized defect (top) and the resulting formation of solitons (bottom) in a condensate of $|1 \rangle$ atoms. The imaging beam was detuned $-30$ MHz and $-20$ MHz, respectively, from the $|\,3S_{1/2},\;F=2\rangle \rightarrow |\,3P_{3/2},
F=3\rangle$ transition. Prior to imaging, the atoms were optically pumped to $|\,3S_{1/2}, F=2\rangle$ for $10\;\mu\mathrm{s}$. The probe pulse had a peak Rabi frequency $\Omega_p = (2 \pi) 2.4$ MHz. The coupling laser had a Rabi frequency of $\Omega_c = (2
\pi) 14.9$ MHz, was turned on $6\;\mu\mathrm{s}$ before the probe pulse maximum, and had a duration of $18\;\mu\mathrm{s}$.\
\
**Fig. 5.** Experimental OD images of a $|1\rangle$ condensate, showing development of the snake instability and the nucleation of vortices. In each case, the BEC was allowed to evolve in the trap for a variable amount of time after defect creation. **(A)** The deepest soliton (nearest the condensate center) is observed to curl due to the snake instability and eventually break, nucleating vortices at 1.2 ms. Defects were produced in BECs with $1.9\times 10^6$ atoms by probe pulses with a peak $\Omega_p = (2 \pi) 2.4$ MHz, and a coupling laser with $\Omega_c = (2 \pi)14.6$ MHz. The imaging beam was $-5$ MHz detuned from the $|\,3S_{1/2},\,F=2\rangle \rightarrow
|\,3P_{3/2},\,F=3\rangle$ transition. **(B)** The snake instability and behavior of vortices at later times. The parameters in this series are the same as in (A), except that the peak $\Omega_p = (2\pi) 2.0$ MHz, the number of atoms in the BECs was $1.4 \times 10^6\!,$ and the pictures were taken with the imaging beam on resonance.

Animation Captions {#animation-captions .unnumbered}
==================

**Animation 1.** An animation, based on 2D numerical calculations, showing creation of a narrow density defect in a BEC by the light roadblock. The parameters and conventions are the same as in Fig. 1B. Successive frames are spaced by $1~\mu$s. The solid curve shows the build-up of $|2 \rangle$ atoms as the probe pulse runs into the roadblock, and the dashed curve shows the corresponding depletion of the density of the condensate of $|1
\rangle$ atoms.

**Animation 2.** Animation of 40 ms of BEC dynamics based on 1D numerical simulations. Parameters and conventions are the same as in Fig. 3A, but the phase is not plotted here. Successive frames are spaced by 0.25 ms. The animation shows a narrow density defect in the $|1 \rangle$ condensate decaying into four solitons due to the steepening of the back edge of the sound waves. The high frequency ripples are due to the nonlinear part of the Bogoliubov dispersion curve. When the solitons reach a point where their amplitude, $\beta$, equals unity they turn around. Upon reaching the center of the condensate, they pass through each other unaffected.

**Animation 3.** Animation showing 30 ms of dynamics of the state $|1 \rangle$ condensate, based on 2D numerical simulations. Parameters and conventions are the same as in Fig. 3B. (The plot range of each frame is $96.8~ \mu$m$~\times~
31.2~\mu$m). Successive frames are spaced by 0.4 ms. The narrow density defect in the condensate decays into several solitons and the deepest solitons decay, via the snake instability, into vortices and release their remaining energy as sound waves. The vortices drift slowly, while some of the sound waves reflect off the condensate boundaries and subsequently interact with the vortices.
