Spontaneous crystallization of light and ultracold atoms

Coherent scattering of light from ultracold atoms involves an exchange of energy and momentum introducing a wealth of non-linear dynamical phenomena. As a prominent example particles can spontaneously form stationary periodic configurations which simultaneously maximize the light scattering and minimize the atomic potential energy in the emerging optical lattice. Such self-ordering effects resulting in periodic lattices via bimodal symmetry breaking have been experimentally observed with cold gases and Bose-Einstein condensates (BECs) inside an optical resonator. Here we study a new regime of periodic pattern formation for an atomic BEC in free space, driven by far off-resonant counterpropagating and non-interfering lasers of orthogonal polarization. In contrast to previous works, no spatial light modes are preselected by any boundary conditions and the transition from homogeneous to periodic order amounts to a crystallization of both light and ultracold atoms breaking a continuous translational symmetry. In the crystallized state the BEC acquires a phase similar to a supersolid with an emergent intrinsic length scale whereas the light-field forms an optical lattice allowing phononic excitations via collective back scattering, which are gapped due to the infinte-range interactions. The studied system constitutes a novel configuration allowing the simulation of synthetic solid state systems with ultracold atoms including long-range phonon dynamics.


I. INTRODUCTION
For a gas of point-like particles off resonantly illuminated by coherent light, the individual dipoles oscillate in phase, each emitting radiation in a characteristic pattern.When several particles contribute to the scattering, the corresponding amplitudes interfere, which leads to a strongly angle-dependent scattering distribution [1][2][3].In addition, if the motional degree of freedom is relevant on the considered time scales, any high field seeking particle will be drawn towards the corresponding local light field maxima, where in turn light scattering is enhanced.This directional energy and momentum transfer between the gas and the field leads to an instability resulting in density fluctuations and potentially also in the formation of an ordered pattern.While for a room temperature gas this typically occurs only at very high pump powers [4][5][6], it can become important for very strong scatterers as larger nano-or microparticles [7][8][9][10][11][12][13]. The stringent threshold conditions can be relaxed by laser cooling the gas to temperatures well below the mK-range as well as by recycling the scattered light in optical resonators.In this case much weaker forces and thus lower light power is needed to a create substantial back-action effect of the scattered light onto the particles.This back-action was predicted to lead to roton-like instabilities and spatial bunching even at moderate pump powers, as observed in several configurations [14][15][16][17][18][19][20][21][22][23][24][25].
A relevant question is thus whether these instabilities can in some cases lead to the formation of a stable crystalline phase in the steady state of such driven, dissi- * stefan.ostermann@uibk.ac.at pative systems.The first and simplest instance of such crystals is the self-ordered phase of transversally driven atoms in optical resonators [26][27][28][29], with the corresponding transition observable also as a quantum phase transition at zero temperature [30,31].It has been shown recently that a similar phase is also realizable in longitudinally pumped ring-cavities [32].
While this self-ordered phase shows some aspects shared by standard crystals like a roton-like mode [33], other characteristic features like the breaking of a continuous translational symmetry and a crystal spacing which is not externally fixed are both missing, since the resonator mirrors select a single electromagnetic mode.In order to include such features one necessarily needs to couple the particles to several electromagnetic modes, ideally a continuum.This is the case in onedimensional tapered optical nanofibers [34,35] or confocal cavities [36], where transversally driven atoms are predicted to spontaneously break the continuous symmetry into a crystal phase.The existence of a continuum of electromagnetic modes opens up the possibility for photons to crystallize, as it was studied with light propagating under Electromagnetically-Induced-Transparency conditions through a nonlinear medium [37,38].
In this work, we propose and characterise a novel crystalline phase of light and ultracold atoms.We consider a mirror symmetric and translation invariant setup as it is depicted in Fig. 1.It involves an elongated Bose-Einstein condensate (BEC) longitudinally illuminated by two counter propagating Gaussian beams far detuned from any atomic resonance.The beams have either orthogonal polarization or a sufficiently large frequency difference to suppress any interference effects.Above a finite driving intensity both atoms and light break a continuous translational symmetry leading to pattern for- An elongated BEC interacting with two counterpropagating, non-interfering laser beams of orthogonal polarization.The two beams are far detuned from any atomic resonance in order to avoid mixing between the two polarizations.Both polarizations are assumed to be equivalent with respect to the considered atomic transition, the latter thus involving a spherically (or at least cylindrically) symmetric ground-state.Alternatively to the use of two different polarizations, sufficiently different frequencies of the two counterpropagating lasers can be chosen.
mation with an intrinsically defined lattice spacing determined by the polarizability and density of the gas.The resulting state corresponds to a supersolid BEC trapped in an emerging optical lattice, the latter showing collective phononic excitations.The appearance of an emergent length scale in combination with lattice phononsi.e. the appearance of a crystal of light -is a crucial difference to configurations where the drive is transverse to the direction in which the system organises [34][35][36].
A useful property of the chosen geometry is that ample information about the coupled system dynamics can be retrieved from the reflected light fields in a completely non-invasive manner.The present study opens a new direction in (ultra)cold atom-lattice physics, naturally including long range phonon-type interactions and realtime non-destructive monitoring.

II. MODEL
We consider a trapped atomic BEC interacting with the electromagnetic (EM) field driven by two far offresonant, counterpropagating, orthogonally polarized laser beams, as depicted in Fig. 1.In the dispersive regime considered below, the EM field provides an optical potential for the BEC (see Eq. ( 1)), while the BEC significantly modifies the refractive index (see Eq. ( 3)), thus both field and matter are dynamical quantities.
The BEC is treated within the Gross-Pitaevskii (GP) mean-field approximation [39], whereby the condensate wave function satisfies the equation where m denotes the particle's mass, g c is the effective swave atom-atom interaction strength, and N is the atom number .For computational simplicity we assume the BEC to be confined by an extra transverse trapping potential V trap (x, y, z) such that the dynamics along the y and z axis is negligible.Therefore, the BEC wave function ψ is assumed to be in the ground state of the transverse trap with characteristic size d y = d z = √ A, where A denotes the BEC cross section.Such a quasi one-dimensional treatment is eligible if the BEC's chemical potential µ is much smaller than the characteristic transverse trap frequency: µ ω y,z .The wavefunction satisfies the normalization condition: dx|ψ(x, t)| 2 = 1.
The total optical potential for the BEC has two contributions: representing the static trapping potential V trap and the longitudinal (along x) optical potential V opt determined by the dynamical part of the injected and scattered EM field (see Eq. ( 5)).The latter consists of two far offresonant fields with orthogonal polarizations driven from the left (L) and right (R) side of the BEC as depicted in Fig. 1.The two polarization components of the field satisfy the Helmholtz Eq. ( 3) The atoms inside the BEC are described as linearly polarizable particles with a scalar polarizability α where the imaginary part is negligibly small, i. e. spontaneous emission of the atoms is neglected.This corresponds to the assumption that the driving laser frequency ω l is sufficiently far detuned form any atomic resonance to prevent substantial internal excitation.This avoids spontaneous emission and thus mixing of the two counterpropagating EM components via Raman scattering as it is used for near resonant polarization gradient cooling may be neglected.
While for spin-polarized atoms the polarizability is field direction dependent in general, we assume the same polarizability for both polarizations orthogonal to the laser axis being the quantization axis.This corresponds to transitions from a spherically (or at least cylindrically) symmetric atomic ground-state.The impinging laser fields from left and right are approximated by plane waves so that we can write the EM field components as E L,R (x, t) = E L,R (x)e iω l t + c.c. e L,R with the orthogonality condition e L • e R = 0.As the light transit time through the sample is negligible compared to all other time scales, the propagation delay of the EM field is adiabatically eliminated and the two field envelops (L for the field from left and R for the field from right) satisfy the Helmholtz equations with the wavenumber k 0 of the incoming beams and the susceptibility χ(x) of the BEC.This susceptibility depends on the condensate's density and is given by where ψ(x, t) is the solution of Eq. ( 1).The directionality of the field propagation in the Helmholtz equations ( 3) is defined by the boundary conditions, according to which the L-component has a finite imposed amplitude on the left end of the system and the R-component has such on the right end (see also Appendix B).
As soon as one knows the spatial distribution of the electric fields, one can calculate the optical potential for the atoms via Inserting the optical potential (5) into Eq.( 1) leaves us with the set of three coupled differential equations: the GP Eq. ( 1) and the two Helmholtz Eq. ( 3), describing the nonlinear dynamics of our system.The degree of nonlinearity resulting from the atom-light coupling is quantified by the dimensionless constant ζ defined as where n = N/AL is the three-dimensional density of the homogeneous BEC with L its characteristic extension along x.Due to the adiabatic approximation involved in the Helmholtz equation, the EM fields depend only parametrically on time through the dynamical refractive index set by the BEC density.Due to the orthogonality of the two chosen polarizations there is no interference between the two counterpropagating components of the EM fields.Therefore, the optical potential (5) only depends on the absolute value squared of the fields.This important feature guarantees the translation invariance of the setup along the x direction nevertheless maintaining a mirror symmetric setup.Indeed, since we are driving with plane-wave lasers, as long as the BEC density is homogeneous, the EM fields E L,R (x) in Eq. (3) are also plane waves, leading to a translation invariant optical potential (5).This invariance with respect to continuous translations is spontaneously broken above a finite driving intensity, as discussed in section III.In the resulting crystalline phase, the lattice constant is intrinsically established as it is discussed in section IV.This is due to the fact that no specific modes are selected and the fields can counterpropagate independently.

III. DYNAMICAL INSTABILITY TOWARDS CRYSTALLISATION
As already mentioned above, due to the orthogonality of the polarizations of the two injected counterpropagating laser fields the particles do not feel any longitudinal optical forces.Naively, one could thus expect the BEC to remain unperturbed independently of the pump intensity.In this section we show that this is actually not the case, as above a particular threshold driving strength small density fluctuations lead to backscattering of light which in turn amplifies these fluctuations.This leads to an instability towards crystallization in the longitudinal direction.The latter can be described by considering the collective excitation spectrum of the system for a spatially homogeneous density distribution of the BEC ψ 0 (x, t) = 1/ √ L with the corresponding propagating field solution of Eq. ( 3).These are plane waves of the form E (0) L,R = C exp(±ik eff x), with the modified wavenumber where C is a real number fixed by the driving strength.
The spectrum is obtained by linearizing the coupled equations ( 1), (3) Here δψ and δE are small deviations from the stationary solutions ψ 0 or E (0) L,R and µ is the BEC chemical potential (we refer to Appendix A for any details).This yields Here I L,R denotes the intensity (in W/m 2 ) of the incoming light which we have chosen to be equal from left and right.
The above analytical expression Eq. ( 8) is very useful to understand some essential features of the atom-light interaction in the present setup and in particular the nature of the crystallisation transition.Apart from the last term, we recognize the known form of the Bogoliubov spectrum of interacting BECs [39], with the linear-in-q behavior corresponding to phononic excitations at low q.The last term on the other hand is the only one resulting from the atom-light interactions.The first thing to note is that its denominator vanishes at q = ±2k eff , which tells immediately that the modified wavenumber (7) sets the favoured momentum for the appearance of the instability.However, the vanishing of the denominator is compensated by the diverging BEC length L, at every finite atom number N (note that ζ ∼ N ).The limit L → ∞ of Eq. ( 8) actually has to be taken, since the stationary plane-wave solution which we linearized only makes sense for a homogeneous and infinite atomic medium, so that the edges may be neglected.This indeed allows us to neglect the reflection of the incident wave by the change in refractive index at the BEC edges.Such finite-size effects, included in the numerical solutions described in section IV, become irrelevant for large systems, as we demonstrate below.
One way to obtain the proper result for Eq. ( 8) in the limit L → ∞ is to consider that for any finite L the allowed momenta q take only quantized values as multiples of 2π/L.Before taking the limit L → ∞ it is instructive to compute the spectrum (8) for fixed finite L, as it is shown in Fig. 2 for different values of I L,R .One recognizes a gap opening at q = 2k eff for any finite I L,R , i.e. any finite driving strength.The spectrum develops a minimum at the finite momentum q = 2k eff + (2π)/L, which corresponds to a roton minimum in the language commonly adopted for standard crystal formation [40].It constitutes a generalization to continuous-symmetry breaking of the roton-like instability observed with a BEC in an optical cavity [41].In a similar manner as in standard crystals, the crystallisation threshold can be calculated by finding the drive intensity at which the roton energy approaches zero.This leads to the threshold-condition ω 2k eff +(2π)/L = 0. We are now in the position to take the limit L → ∞.In doing this we note that we have to keep the atom number N constant in order to get a finite critical drive strength.Otherwise, if we perform the standard thermodynamic limit N/L = const., the energy of the system diverges and the crystallisation threshold vanishes.This divergence is an artefact of our model in which the light-mediated atomatom interaction is of infinite range since the EM field is adiabatically adapting to the BEC configuration.The inclusion of the dynamics of the EM field (retardation effects) would introduce a finite range and thus eliminate the divergence in the energy.Still, the resulting range is expected to be larger than the typical BEC size L so than our calculation should be valid for any realistic sys- tem size.Taking the L → ∞ limit we thus get the critical driving intensity where we introduced the recoil energy Note that in the L → ∞ limit with constant N the BEC becomes more and more dilute, which renders the direct atom-atom coupling ∼ g c eventually irrelevant.In Fig. 3 the analytical expression ( 9) is compared with numerically estimated thresholds for large system sizes (see section IV).We find full agreement between the linear instability threshold and the numerical threshold found by studying the imaginary time evolution of eqns.( 1) and (3).This numerical approach to finite-sized systems is described next.

IV. CRYSTAL OF LIGHT AND ATOMS
After showing that the homogeneous system is unstable above a certain driving intensity, we are going to show that a stable crystalline phase is reached and study its properties by numerically solving the coupled GP (1) and Helmholtz (3) equations.We perform an imaginary time evolution of the system (1)-(3), i. e. replace t → iτ , which yields the ground state of the system for long enough evolution times.For a detailed description of the numerical methods we refer to Appendix B.
To determine the crystal transition point as a function of driving intensity we compute the total reflectivity of the BEC with respect to the intensity of either one of the incident beams, which we again take to be equal.For large enough system sizes, a clear threshold behavior is visible at a critical driving intensity, whereby the reflectivity grows from essentially zero with almost infinite slope, cf.Fig. 4. The hereby found critical intensity is in  2 perfect agreement with the analytical result(see Fig. 3).As mentioned already in the previous section, finite-size effects manifest due to the presence of the edges of the BEC.In the calculations described in this section and in section VI, there is no further trapping potential along x and the BEC is confined within a box of size L, so that the BEC has sharp edges for the light impinging at x = 0 and x = L (see Appendix B for more details).In section VII we add an harmonic trap along x and show that the qualitative behavior is the same as described here.The BEC edges create a quick increase of the refractive index which induces a small amount of reflection of the incoming beam.As apparent from Fig. 4 this reflection is irrelevant for large system sizes compared to the reflection present in the crystalline phase.
The large light reflection above threshold is due to the appearance of a large spatial modulation of the BEC, forming the density grating shown in Fig. 5(a).This corresponds to a continuous symmetry breaking at the threshold leading to a crystalline phase which for the phase-coherent BEC implies supersolid order.Each peak in the density grating reflects the incoming light, resulting in a damped modulation of the intensity of each polarization component across the condensate, as shown in Fig. 5(b).While the modulation of each component's intensity I L,R is damped across the system, the modulation of the total intensity I tot = I L + I R is not damped, resulting in a periodic optical-lattice potential for the BEC which matches its density grating.
An important feature of the optical lattice emerging in the crystalline phase is the intrinsic character of the lattice spacing, which is not fixed externally but rather set by the BEC density and atom polarizability.This is a clear difference with respect to the self-ordering in optical resonators where the spacing is externally fixed by the cavity mirrors [26]; and also to the case of self-ordering of transversally driven atoms coupled to the continuum of modes of optical fibers, where the spacing is fixed by the driving frequency and fiber dispersion [34,35].As anticipated in section III, the appearance of the roton-like instability at the characteristic momentum 2k eff leads to the following prediction for the emergent lattice spacing: The emergent spacing is always smaller than the one in vacuum λ 0 /2.This feature can be qualitatively reproduced also within a toy-model where the medium is approximated by a set of beam splitters [42].This typically small but nonetheless crucial effect is also present when using counterpropagating beams with equal polarization and is essential for atom trapping in optical lattices [43].Would the atoms indeed be trapped with the vacuum spacing λ 0 /2, the EM field would be perfectly reflected and no standing wave could actually be formed and thus no trapping be possible.It is only through the slight renormalization d < λ 0 /2 that perfect reflection is avoided.What our scheme with orthogonally polarized counterpropagating beams allows is to make the small renormalization of d coincide with the appearance of a large density modulation out of a homogeneous phase, i. e. a crystallisation.The existence of an intrinsic lattice spacing in the crystalline phase implies as well the presence of phononic excitations of the the lattice, as discussed in the next section.

V. EXCITATIONS OF THE CRYSTAL: PHONONS
Further insight in the properties of the atom-light crystal is provided by analyzing its excitation spectrum.As done in Section III, we linearize the coupled system of Eqs.(1),(3).However, now the perturbation is performed around the symmetry-broken stationary solution.The result is presented in Fig. 6 for a driving intensity slightly above threshold.Details of the calculation are given in the Appendix A. Since translation-invariance is broken, the matrices describing the linear system are not diagonal in momentum space requiring a discretization of the position(momentum) continuum.Moreover, while the total light intensity and atom density are periodic, the intensity of each polarization component is not, due to accumulated reflection along the density grating, introducing the decaying evelope shown in Fig. 5(b).This prevents the use of the quasi-momentum to label the excitation modes.
In Fig. 6 we labelled the eigenvalues based on their dominant momentum component q max , extracted from the corresponding eigenvector.This allows to split the spectrum into three regions separated by gaps at q max = k eff and q max = 2k eff .
The gap at q max = k eff opens up for I L,R > I L,R c due to the appearance of an optical lattice potential for the atoms with a π/k eff periodicity.It separates the two bands which, slightly above threshold, are characterized by eigenvectors with a clearly dominant momentum component (see left and middle inset in Fig. 6).
On the other hand, the gap at 2k eff is the same one appearing in the homogeneous phase (see Fig. 2).As discussed in Section III, at the critical drive intensity I L,R c the gap is such that the energy of the mode with momentum q = 2k eff + 2π/L (momentum is still a good quantum number for I L,R ≤ I L,R c ) vanishes.Out of this zeroenergy mode at 2k eff (not resolved with the discretization of Fig. 6), and beyond the critical point: I L,R > I L,R c , the lattice phonon branch develops for q max > 2k eff .The momentum distribution of the lattice-phonon eigenvectors is characterized by the splitting of the single peak at 2k eff into two neighbouring peaks (see rightmost inset of Fig. 6).The phonon wavelength is set by the distance between the two nearby maxima appearing in the momentum distribution.This generates the slow beating in coordinate space.With a finite system size L, the longest wavelength is of the order of L.
Moreover, the lattice-phonon branch is gapped, in the sense that its lowest energy mode at q max slightly above 2k eff has a finite energy, as visible in Fig. 6.More impor-tantly, this gap remains finite in the thermodynamic limit L → ∞.We can estimate the size of the lattice-phonon gap close to threshold by using Eq. ( 8) and computing the energy of the mode next to the zero-energy mode.This yelds which takes the value ∆ 2 ph 8E 2 rec in the thermodynamic limit L → ∞ with N = const..As discussed in section III, in this limit I L,R c remains finite while n → 0 and k eff → k 0 .Another choice of thermodynamic limit is possible: L, N → ∞ with n = const.,where I L,R c → 0 and the gap is still given by (11).The existence of an energy gap for lattice phonons is due to the long-range nature of the interactions, as it can be already predicted within a classical model of interacting point-like particles [44].From a more general field-theoretical perspective, some of the gapless Goldstone modes expected from the continuous-symmetry breaking can indeed disappear (i.e.become gapped) due to the long range of the interactions, as it for instance happens to the longitudinal phonons of a three-dimensional Wigner crystal [45].As long as retardation effects can be neglected our interactions will be infinite-ranged, the lattice-phonons gapped, and thus quantum/thermal fluctuations will not destroy crystalline order even in truly one dimension [46].
The existence of lattice phonons among the collective excitations is confirmed by numerical simulations of the real-time dynamics of the system, as it is described in the next section.

VI. CRYSTALLISATION DYNAMICS AFTER A QUENCH
In this section we investigate the real time dynamics of the system by directly solving eqns.( 1) and ( 3).This allows us to analyse the crystallization dynamics after a sudden turn on (quench) of the pump laser strength from zero to a value above threshold at t = 0.The corresponding time evolution of the BEC reflectivity, kinetic energy, as well as the evolution of the BEC density and total light intensity are shown in Figs.7 and 9.
As apparent from the behavior of the reflectivity and kinetic energy E kin (t) = dx 2 |∂ x ψ| 2 /2m, the crystalline order is reached after a few inverse recoil frequencies, after which both quantities perform oscillations about a finite value.These residual oscillations are triggered by the energy gained by the system upon forming the density grating together with the optical lattice.The reason why this effect takes in a prominent role in the studied case is found by looking at Fig. 8, which shows the zoom into two peaks of the intensity distribution of the crystal.One recognizes that the maxima of the intensity distributions of the two fields coming from left and right (blue dots in Fig. 8) do not coincide with the maximum of the total intensity distribution (black dot in Fig. 8) at which the atoms are trapped.Therefore, the trapped atoms feel a strong field gradient for each single component because they do not sit at the maxima of the two counterpropagating fields, as it would be for example the case in optical lattices.This leads to a large coupling between the two counterpropagating fields and the atoms, leading to strong long-range interactions inducing collective excitations.The corresponding dynamics of the BEC density and the total light intensity is shown in Fig. 9.As one can see from the solid lines marking the evolution of the intensity maxima, they start at a lattice spacing of λ 0 /2 and move closer together in time reaching the emergent spacing d.In addition, we see the presence of residual oscillations about the crystalline order.In particular, the light intensity shows both compression modes, modulating the amplitude of the optical lattice in time, and phonons, modulating the spacing.The latter are clearly visible from the dynamics of the intensity maxima shown in Fig. 10.Since we are neglecting retardation of the fields, the energy can be redistributed among the collective degrees of freedom but not dissipated.Initially, for ω rec t ∼ 1, mostly compression modes are excited.Subsequently part of the energy stored in compression modes is transferred to lattice phonons for ω rec t > ∼ 5.In Fig. 10 we see a single-frequency oscillation of the intensity maxima, the latter moving almost in phase.This indeed corresponds to a low-wavelength lattice phonon, which becomes occupied for long enough times.As discussed in section V, the longest wavelength is of the order of the system size L, consistently with the almost in-phase oscillations of Fig. 10.
As discussed in the previous section, lattice-phonons have a finite gap.They can efficently be excited in a quench experiment provided the energy available for collective excitations is large enough compared to ∆ ph (see Eq. ( 11)).

VII. EXPERIMENTAL IMPLEMENTATION WITH ULTRACOLD BOSONS
BECs with high densities and a controlled shape trapped in optical dipole traps are currently available in many laboratories.In principle, the setups normally employed are already very close to the one needed to study the crystallization effects presented in this work.In the following, we will discuss the conditions needed to study our model in realistic experimental conditions, as well as the required parameter regime for observing the crystallization.Let us remark that the basic physics underlying the crystallization transition discussed here does not rely on the atoms being Bose-condensed.This phenomenon could in principle also be observed with thermal clouds or fermionic gases.Apart from the fundamentally very interesting feature of supersolidity, the practical advantage of a BEC with respect to a thermal cloud resides in its high density and low temperature, both decreasing the required laser power.On the other hand, for degenerate fermi gases, one could expect a strong dependence on the ratio between Fermi momentum and lattice constant [47][48][49][50].
We start by noting that using single beam optical traps can also lead to heating instabilities but never generate a stationary lattice [17].Similarly, operating very close to an atomic resonance has been shown to generate instabilities and a short time formation of an optical lattice structures via so called end-fire modes [16].As this requires significant atomic excitation, it involves fast transverse acceleration with heating and destruction of the BEC.This is prevented in our model by an improved geometry and much larger atom-field detuning.
Our model ( 1)-( 3) is essentially 1D, which relies on the assumption that both, the atoms and the light move and propagate essentially unidirectionally along x.In practice this can be implemented by using a transverse trapping of the atoms tight enough to freeze out the dynamics along y, z.With harmonic trapping potentials this amounts to the requirement that ω ho y,z is sufficiently larger than the BEC chemical potential µ.Here we still describe the one-dimensional BEC using the GP equation, which requires the atom density to be large enough to be in the mean-field regime [39].The enforcement of unidirectional propagation of light is more demanding since an appreciable amount of diffraction out of the BEC axis would be present inducing propagation also perpendicular to x. Apart from the use of hollow-core optical fibers around the BEC [51], one option available in many laboratories today is using a two-dimensional array of tubes with spacing comparable with the wavelength of the light.This arrangement would generically produce destructive interference between the transverse field components diffracted from different tubes, so that if the latter are long enough only the forward propagation along the tube axis would remain.In this configuration, each tube will act equally while the field propagates inside a medium with a refractive index given by the sum of the contributions from each tube.Indeed, since all tubes share the same backreflected field there is a natural synchronization of the different tube lattices.
In any experimental realization a trap to confine the BEC along x will also be present.In addition, the two laser intensities might differ to some extend due to experimental inaccuracies.As an exemplary case we study the crystallization as in section IV but add a harmonic trapping potential V ext (x) = Etrap 2 x 2 /λ 2 0 and chose different pump intensities I l = I r .It can be seen from Fig. 11 that the qualitative features of the crystalline phase remain the same as in the homogeneous case.The only difference is the parabolic envelope for the density as well as for the light intensity distribution and the shift of the distribution towards the direction of the higher intensity.The threshold behaviour remains similar to the one presented in Fig. 4 with the only difference being an increase of the threshold intensity.A useful feature of the considered configuration is that the crystallization process can be observed in real time by looking at the amount of reflected light, since the transmitted part of the counterpropagating beam can be separated from the reflected part having orthogonal polarization.
In order to choose the most suitable atomic transition, pump detuning and power, as well as BEC parameters like density and extension, one must consider the following constraints: we need to have i) a low enough critical driving strength (9), which depends on the detuning ∆ a and spontaneous emission γ through the real part of po-  larizability Reα ∼ γ/∆ a , reading and at the same time ii) a low enough BEC heating rate, which at the critical power reads with n A = N/A being the surface density of the medium with respect to the light propagation.From Eqs. ( 12) and ( 13) one sees that the crystallisation is easier achieved before the BEC is heated up if we increase the BEC surface density n A .There is no favorable scaling neither with detuning ∆ a nor with the linewidth γ, since both heating rate (13) and critical power (12) scale with γ 2 /∆ 2 a .For commonly employed transitions like the Rb or Cs D lines, the required laser power is easily achieved, but the heating rate can become a problem at too low densities due to the required laser powers and detunings.For instance, taking N = 10 6 atoms confined over a transverse cross section A ∼ 5 × 5µm 2 and λ 0 ∼ µm, we estimate a required power I c ∼ W/cm 2 with a heating rate Γ heat ∼ 10 Hz for the Rubidium 780nm line with a detuning ∆ a = 100GHz as well as for the Cesium D2 line with a detuning ∆ a = 20GHz.Such a heating rate would still allow to observe the crystal fromation since, as we see in Fig. 9, this process takes place on the inverse recoil time scale, which is of the order of ms.

VIII. CONCLUSIONS AND OUTLOOK
We predict that in suitable geometries roton instabilities originating from non-linear free space atom light interactions can be tailored to generate stationary crystalline states.They involve an optical lattice showing an emergent spacing and phononic excitations, trapping the atoms at the intensity maxima.
The required translation invariant, mirror symmetric geometry can be realized using two orthogonal polarization degrees of freedom or frequency shifted counterpropagating beams.We estimate that the dynamics studied in this work should be accessible in already existing experimental setups on large quasi-1D Bose-Einstein condensates.Actually, in comparison with standard crossed beam dipole traps, one simply has to adapt and control the polarizations of the trapping lasers and choose suitable detunings.The ordering process should be easily observable not only by measuring the atomic distributions but directly by looking at the reflected light from the condensate.This non-destructive measurement allows for a real-time monitoring of the dynamics.
Our results open up an intriguing new direction in quantum simulations with ultracold atoms in optical lattices, where the latter are enriched by the presence of collective phononic excitations resulting from the spontaneous crystallisation of light.In this spirit, the application of our approach to two-dimensions and the inclusions of retardation effects as well as quantum fluctuations constitute the natural extension of this study.
The last line of equation (A13) corresponds to the stationary GP equation and therefore it vanishe, as it defines how the chemical potential is related to the field amplitude and the particle particle interaction g c , namely via where I tot (k) and n 0 (k) are the Fourier transforms of the total intensity distribution and the BEC density.Besides we defined the additional matrices V L,R (k, k ) := Let us now define the spinor Ψ(q) := (ψ(q), ψ * (q)) T where ψ(q) defines a single momentum component of ψ ψ ψ arbitrary initial condition (B1) and (B2), leading to well defined fields at the boundaries allowing for an estimation of the right hand amplitudes C R and D R .Hence, R and T can be calculated via (B3) and (B4).As soon as we know the initial conditions we can find the solution of the Helmholtz equation via spatial integration performed by a fourth order Runge-Kutta solver.
The solution of the HH equation is then used to calculate the optical potential (5).The time evolution of the GP equation with the newly found potential is then calculated by using a split step method.Note that the HH equation has to be solved within each time step resulting in a modified potential for the next time step in the GP equation.The time evolution is finished as soon as the system is found in a stationary state.

Figure 1 .
Figure 1.Schematic representation of the considered setup.An elongated BEC interacting with two counterpropagating, non-interfering laser beams of orthogonal polarization.The two beams are far detuned from any atomic resonance in order to avoid mixing between the two polarizations.Both polarizations are assumed to be equivalent with respect to the considered atomic transition, the latter thus involving a spherically (or at least cylindrically) symmetric ground-state.Alternatively to the use of two different polarizations, sufficiently different frequencies of the two counterpropagating lasers can be chosen.

2 Figure 4 .
Figure 4. Dependence of the reflection coefficient of the BEC on the incoming field amplitudes for different atom-field couplings ζ = 0.1 (dashed red) and ζ = 0.2 (solid blue).The remaining parameters are the same as in Fig. 2

Figure 5 .
Figure 5. (a) Crystal ground state for ζ = 0.1, I l = Ir = 200 and (b) corresponding intensity distribution for the field from left (green) and right (red).The solid blue line depicts the sum of both intensities.A zoom into the yellow shaded region can be found in Fig. 8.The remaining parameters are gcN/Aλ0 = Erec and L = 10λ0.

Figure 6 .
Figure 6.Excitation spectrum of the atom-light crystal.The blue points are the eigenvalues of the GP and Helmholtz equations linearized about the crystalline stationary state (see Appendix A).The numerical diagonalization is performed with a momentum-space discretization dq = 2π/L.The parameters are the same as in Fig. 2 except for L = 50 and a fixed drive intensity I L,R = 50 (slightly above threshold).qmax is the momentum corresponding to the largest component of the eigenvector of each eigenvalue.The insets show examples of eigenvectors (unormalized probability in momentum space) for three different eigenvalues representative of each region of the spectrum, from left to right: λ0qmax = 1.38 < λ0k eff , λ0k eff < λ0qmax = 11.9 < 2λ0k eff , and λ0qmax = 12.7 > 2λ0k eff .The latter region corresponds to lattice phonons, characterized by a two symmetric pairs of peaks about a finite momentum.This phononic branch: qmax > 2k eff has a gap ∆ ph .Its analytical estimate in Eq. (11) yelds ∆ ph 2 √ 2Erec, in reasonable agreement with the numerical data.

2 Figure 7 .Figure 8 .
Figure 7. (a) Real-time evolution of the kinetic energy for ζ = 0.1, I l = Ir = 100, gcN = 1.(b) Real-time evolution of the reflection coefficient for the same parameters as in figure (a).The solid black line shows the mean value of the corresponding functions.

Figure 9 .
Figure 9. Real time dynamics of the (a) BEC density distribution and b) the total light intensity for the same parameters as in Fig. 5.The solid black lines in figure (b) show the time evolution of the intensity maxima.

Figure 10 .
Figure 10.Real time evolution of the maxima of the intensity distribution as it is shown in Fig. 9.To simplify the comparison between the single curves the maxima positions were shifted so that they all start at x = 0. Fig.(a) shows the total time evolution where one can clearly recognize collective phonon-like excitations of the lattice after ωrect > ∼ 5. Fig.(b) shows the zoom into the yellow marked area in fig.(a) in order to demonstrate the slight dephasing between the oscillations of the maxima.All parameters are chosen as in Fig. 5.

Figure 11 .
Figure 11.a) Crystal ground state for and b) corresponding intensity distribution for the field from left (green) and right (red) for the same parameters as in Fig. 5 with an additional external potential Vext(x) = E trap 2 x 2 /λ 2 0 with Etrap = 1.0Erec and for different pump intensities from left and right I l = 200 and Ir = 150.The solid blue line depicts the sum of both intensities.