Dynamical crystal creation with polar molecules or Rydberg atoms in optical lattices

We investigate the dynamical formation of crystalline states with systems of polar molecules or Rydberg atoms loaded into a deep optical lattice. External fields in these systems can be used to couple the atoms or molecules between two internal states: one that is weakly interacting and one that exhibits a strong dipole-dipole interaction. By appropriate time variation of the external fields, we show that it is possible to produce crystalline states of the strongly interacting states with high filling fractions chosen via the parameters of the coupling. We study the coherent dynamics of this process in one dimension (1D) using a modified form of the time-evolving block decimation (TEBD) algorithm, and obtain crystalline states for system sizes and parameters corresponding to realistic experimental configurations. For polar molecules these crystalline states will be long-lived, assisting in a characterization of the state via the measurement of correlation functions. We also show that as the coupling strength increases in the model, the crystalline order is broken. This is characterized in 1D by a change in density-density correlation functions, which decay to a constant in the crystalline regime, but show different regions of exponential and algebraic decay for larger coupling strengths.


I. INTRODUCTION AND OVERVIEW
Recent experimental progress in controlled excitation of Rydberg atoms from cold gases [1][2][3][4][5][6] and in the production of ultracold polar molecules [7] has opened new opportunities for the study of systems with dipole-dipole interactions. Key features of these systems include the ability to control these interactions by varying the external fields, as well as the strong dependence of the interaction strength on the internal state of the atoms or molecule. An important example of this dependence is the Rydberg blockade mechanism [8,9] in which an atom in an excited Rydberg state can prevent nearby atoms from being excited to Rydberg states, because the strong interaction between the two excited atoms makes subsequent excitations non-resonant. This mechanism has been recently demonstrated in experiments [10,11], and has important potential applications to the production of fast quantum gates with two or several neutral atoms [8,9,12,13]. Moreover, the blockade produces rich and intriguing collective dynamics when atoms in a dense gas are excited to Rydberg states [1][2][3][4][5][6][14][15][16]. This cooperative behaviour is manifest in a suppression of the proportion of excited atoms with increasing density (due to the presence of more atoms within the blockade radius), and also in an enhanced laser coupling to collective Rydberg excitations.
In the context of Rydberg atoms, there has been a lot of interest in the possible production of crystalline states, where a small fraction of atoms in excited Rydberg states arrange themselves in a regular array [17][18][19][20]. In particular, when ground and excited states are coupled by a laser field, the resulting effective models correspond to crystalline states in the limit of weak coupling, with the fraction of excited atoms determined by a competition between the interaction energy and the detuning of the laser field from the atomic transition.
Motivated by this experimental and theoretical progress, we consider here a collection of either polar molecules or neutral atoms loaded into a deep optical lattice, e.g., with one particle per site. We assume that tunnelling in the lattice is negligible on the timescale of the experiment. Using dressed rotational states for the polar molecules or Rydberg states for the atoms we show how a microscopic model can be realized where two internal states of the particles are accessed, one of which exhibits a large dipole-dipole interaction, whilst the other is weakly interacting. We show that this system allows the realisation of crystalline formations in the internal states, with a periodicity related to the lattice period by a rational fraction. The regularity of the system makes it possible to realize crystals with high filling factors of strongly interacting states. In addition, the use of polar molecules makes it possible for the crystalline states to be long-lived, aiding in the measurement of characteristic correlation functions. We show that in a 1D tube, increasing the coupling strength between the internal states at a fixed detuning will break the crystalline order. This is characterized in 1D by a change in the density-density correlation (DDC) functions, which decay to a constant in the crystalline regime, but show different regions of exponential and algebraic decay for larger coupling strengths.
We study in detail the dynamical preparation of states with crystalline order of different periodicities in the internal states by an appropriate time-variation of microwave or laser couplings. As was shown in previous studies for Rydberg atoms in an optical lattice modelled by nearest-neighbour interactions [14][15][16], it is not possible to produce crystalline states simply by switching on the coupling field. Instead, this is made possible using an adiabatic variation of the parameters in the effective model. Using an extended version of the Time-Evolving Block Decimation (TEBD) algorithm [21,22], we calcu-late the coherent time-dependent many-body dynamics of the excitations. We show that for system sizes and parameters corresponding to realistic experimental configurations with either polar molecules or Rydberg atoms we obtain fidelities > 99% compared with the ideal target states, i.e., the ground states of the effective model. Whilst the fidelity could be reduced by incoherent processes, e.g. decay of excited states, we estimate that these rates are small, and that the basic structure of the crystalline state will be robust. Note that the time-dependent preparation of crystalline phases here follows the same philosophy as that for production of Rydberg crystals in reference [20]. In focussing on the case of atoms or molecules initially trapped in a lattice, we concentrate here on the case with a high proportion of particles in the strongly interacting case. This is in contrast to reference [20], which treats the case of a dilute gas of Rydberg excitations that is more relevant for a frozen Rydberg gas not trapped in a lattice. Our study of the correlation functions in the crystalline and non-crystalline phases is, however, also relevant for the case of [20].
The remainder of the paper is organized as follows: in section II, we introduce the effective long-range spinmodel that describes both the Rydberg atoms and polar molecules in appropriate parameter regimes. Then we briefly discuss in section III extensions to the TEBD algorithm that we use to study the model. In section IV, we analyse the ground states of this effective model and identify a target state of our time-dependent excitation process. The ground states are characterized as belonging to different quantum phases via the characteristic behaviour of DDC functions. In section V, we discuss a scheme to prepare crystalline states with high fidelity by an appropriate time-dependence of the external coupling field parameters. In section VI, details of two different physical realizations of the effective model, with polar molecules (section VI A) and Rydberg atoms (section VI B) are provided. Finally, section VII, provides the conclusion and outlook.

II. THE SYSTEM
We consider particles loaded into an optical lattice in a regime where the lattice is sufficiently deep that tunnelling can be neglected on the timescale of the experiment. We assume that the particles have two internal states, one of which exhibits a strong dipole-dipole interaction, whilst the other is weakly interacting. This can be realized for realistic experimental parameters either with 1. Polar molecules in an external electric field, where a microwave coupling of rotational states can give rize to dressed states with different dipole-dipole interactions, or 2. Neutral atoms in an external electric field coupled via a laser to excited Rydberg states, which in con-trast to the ground states exhibit a strong dipoledipole interaction.
Below we summarize the key points of each of these implementations. Further details for polar molecules and neutral atoms coupled to excited Rydberg states are given in section VI A and section VI B respectively.

A. Polar molecules
Polar molecules in their electronic-vibrational groundstate manifold possess a number of features that make them appealing to design and implement specific spinmodels. In addition to their electric structure, allowing e.g. for optical trapping, they have a rich internal (rotational) symmetry, with level spacings in the GHz regime, and negligible decoherence. They possess a permanent dipole-moment d 0 (typically on the order of a few Debye [23,24]) along the internuclear axis, which makes it possible to strongly couple their rotational excitations via static and microwave fields and gives rise to dipole-dipole interactions between different molecules.
FIG. 1. Schematic illustration of a system of polar molecules loaded into subsequent sites of an 1D optical lattice with site separation a. The molecules are aligned by a static electric field in ez-direction, E dc , and exposed to two time-dependent microwave fields E z 1 (t), E z 2 (t), both linearly polarized along the z-axis.
We consider a 1D optical lattice (lattice spacing a) along the e x axis with one polar molecule trapped at each lattice site in the lowest vibrational state φ k (x) (see figure 1), which can be realized either by loading an insulator state of polar molecules, or beginning deep in a Mott Insulator state with two atoms per site and associating them to form polar molecules [25]. Without external fields, the rotational states of the molecules do not possess a net dipole moment. However, finite dipole moments along the z-axis in the rotational states can be induced by applying a static external electric field E dc in that direction. The key idea of our implementation is that a suitable choice of E z 2 (t) makes it possible to cancel the induced dipole moment in one of the dressed states, and we denote the complete state of an atom in this dressed state at site k as |α k . The single molecule ground state, which obtains a strong dipole moment from the static electric field, is denoted as |β k . In the case that E dc and E z 2 (t) are properly tuned, two molecules at distinct sites k and l that are in the states |β k and |β l lead to an energy shift due to the dipole-dipole interaction potential of V dd /|k − l| 3 . In section VI A we give example experimental parameters showing that an energy shift V dd ≈ 0.7d 2 0 /4π 0 a 3 ≈ 10 kHz can be achieved. Note that we can therefore assume that the molecules remain in their motional ground state, since the band separation (trapping frequency) ω T is typically of the order of many tens of kHz up to 100 kHz in a realistic scenario. Finally, the states |β k and |α k on all sites k are coupled by a weak microwave field E z 1 (t) with Rabi frequency Ω and detuning δ, which gives rise to the effective microscopic many-body model that is described in section II C below.

B. Neutral atoms coupled to Rydberg states
Neutral atoms in an external electric field exhibit strong dipole-dipole interactions when in certain Rydberg states with relatively large principal quantum number. For alkali atoms in states with principle quantum number n = 17, it is possible to obtain a dipole moment of more than one kiloDebye (see section VI B for details). By finding isolated Rydberg states, the internal dynamics of each atom can be reduced to a two level system, where weakly interacting atoms in the ground state are coupled via a laser to these states with large dipole-dipole interactions.
We assume that either atoms are loaded into a deep optical lattice with tight radial confinement, so that we have one atom per site along a 1D tube, or that we have weak radial confinement and many atoms per site (as occurs, e.g., when an optical lattice is applied along only one direction). In each case, we can assume that the ground-state atoms are confined to a single modefunction at each given lattice site, φ k (x). In the case of a single atom per site, this is the Wannier function for the lowest Bloch band, which is analogous to the case of the polar molecules described above. For many atoms per site, interactions can populate higher bands, and φ k (x) will depend on the interactions. However, we assume that the timescale for excitation to the Rydberg state is much shorter than timescales corresponding to the atomic motion, and that therefore we can neglect population of any modefunction other than the initial modefunction φ k (x) for either the ground or Rydberg states. We then assume that we can have at most one excitation per lattice site, which in the case of large single-site populations is made possible by the Rydberg blockade mechanism, when we ensure that the Blockade radius is larger than the occupied region at any single lattice site (see section VI B). This makes it possible to describe the system with an effective spin-1/2 model, assigning a state on each lattice site to the existence or non-existence of an atom in an excited Rydberg state on that site, labelled for site k by the states |β k and |α k , respectively.
In comparison with polar molecules, the dipole-dipole interaction of the strongly interacting state for Rydberg atoms is much larger. If we choose a principal quantum number n = 14, then V dd ≈ 7 GHz (assuming a lattice spacing of a = 400 nm, see section VI B for more details). This corresponds to much faster experimental timescales. However, the polar molecules are much longer lived than atoms in excited Rydberg states, and open up different opportunities for experimental observation of the crystalline states we discuss here.

C. The effective model
In either the polar molecule or neutral atom case we can derive the same effective microscopic model for the system in a frame rotating with the frequency of the field coupling the two internal states. We assume that at most one particle per site can be in the strongly interacting state, and represent the two possible states (weakly and strongly interacting) at each lattice site as two spin states. We label these states for site k as {|β k , |α k }, where |β k is the state exhibiting strong dipole-dipole interactions, and |α k the non-interacting state. The Hamiltonian can then be formulated as ( ≡ 1) withn k ≡ |β β| k being the local "number operator" for the interacting state at site k. The Pauli matrices are defined viaσ x k ≡ |β α| k + |α β| k andσ z k ≡ |β β| k − |α α| k , where the latter is related to the number operators vian k = (σ z k + 1)/2. Ω denotes the effective Rabi frequency and δ the detuning of the laser, which excites the particles. V dd is the dipole-dipole interaction energy shift between neighbouring sites. Note that we assume the particles to remain in the motional states φ k (x), and also neglect corrections to the |k−m| −3 decay due to the finite width of φ k (x) (see section VI for more details).
Below we will first study the ground state properties of this model, which will give us target states for the dynamical preparation process. We will then study the time dependence of this system starting from a manybody state where all of the particles are weakly interacting ( k |α k ), and investigate the adiabatic preparation of crystalline states as the parameters are varied.
It is interesting to note that for an infinite system, the Hamiltonian can be written aŝ where we have removed constant terms, and ζ(3) ≡ k k −3 ≈ 1.202 denotes Apéry's constant [27]. Note that in the case of δ = ζ(3)V dd , the system is equivalent to an Ising model in a transverse field with dipolar interactions. Ground states of this model have been also studied in [28].

III. NUMERICAL METHOD
To calculate the many-body dynamics and ground states of Hamiltonian (1), we make use of the TEBD algorithm as introduced in [21,22]. This algorithm makes possible the near-exact integration of the Schrödinger equation for 1D lattice and spin Hamiltonians based on a matrix product state ansatz, and has been applied to a range of such Hamiltonians with next-neighbour interactions. There are also applications of matrix product state algorithms to dissipative systems [29][30][31], and these methods have been incorporated within density matrix renormalization group algorithms [32,33]. There is also a strong effort to extend these ideas to higher dimensions [34].
The modification in our work is the extension of these methods to the treatments of finite-range interactions. For the purposes of this study, we include long-range interactions up to an arbitrary distance of l sites, by implementing a routine for swapping site indices in the matrix product state representation. This adds a computational cost of O(l) basic operations compared to the standard TEBD calculation. Here we work with finite systems of size N , and always properly represent interactions over the full length of the system, i.e. l = N −1. This method of swapping also allows us to consider periodic boundary conditions for the dipole-dipole interactions. Note that we always perform convergence tests in the matrix product state bond-dimension χ and other numerical parameters [35]. For larger system sizes, swapping sites can become inefficient, in part due to the higher required values of χ. To implement interaction terms over substantially larger distances than in this study, more efficient algorithms should be used, e.g., algorithms involving use of matrix product operators [36][37][38].

IV. GROUND STATES OF THE EFFECTIVE MODEL
We begin by studying the key ground state properties of the effective model of Hamiltonian (1). This will allow to us to determine in which parameter regimes crystalline order will appear, and to characterize the states that will be the target states of the dynamical preparation process discussed next in section V.
The Hamiltonian we study here has strong similarities to another model treated in several recent studies, in which the dipole-dipole interactions were replaced with van der Waals interactions arising for certain Rydberg systems [16][17][18][19][20]. In [17] it was shown that for Ω = 0 that system reduces to a classical spin model exhibiting a second order phase transition from a paramagnetic to a crystalline phase at δ = 0, and we observe similar behaviour here. Due to the dipolar long range character we expect crystalline states with different periodicities of the excitations to the strongly interacting state to appear for δ > 0 and small values of Ω. The detuning plays the A sketch of the key features of the phase diagram for the system with Hamiltonian (1). For δ < 0 and Ω/V dd 1, the ground state involves all of the atoms in the weakly interacting state in order to minimize the interaction energy (marked by the blue lobe). As the detuning δ > 0 is increased, it becomes favourable to add excitations to the system, but this competes with an increasing interaction energy. This competition gives rise to crystalline phases with periodicities determined by this competition (marked by the brown lobes). For large Ω/V dd , the crystalline order is broken, but some particles exist in excited states (white area).
role of a chemical potential for strongly interacting states (or external magnetic field in the spin model), and it becomes favourable to excite states with a dipole-dipole interaction on the lattice. However, the long-range interactions will compete energetically with the detuning, leading to crystal periodicities that decrease with increasing δ/V dd . Below we show that as the coupling Ω/V dd is increased, the crystalline order is broken, in favour of a paramagnetic phase in the spin language, with no long range density-density ordering. These key features of the phase diagram are sketched in figure 2.
We compute the ground state of equation (1) using the TEBD algorithm in several different parameter regimes. As a first step, we look at the total number of particles in a strongly interacting state |β , In figure 3 we present a plot of N β for several detunings δ and effective Rabi frequencies Ω. We plot N β on a grid with spacings ∆δ = 0.125V dd and ∆Ω = 0.025V dd . As expected, for small Ω/V dd there is essentially no occupation of the strongly interacting states in regions with negative detuning, because a state with no |β k excitations is energetically favoured. In contrast, a positive detuning leads to a reduction of the total energy when occupations of strongly interacting states are added to the system, and in regions with large positive detuning we find states with occupations of |β k on all lattice sites. In the region between we observe varying excitation numbers depending on the choice of parameters. In particular, there is a large plateau in the region 0.5V dd δ 2V dd and Ω 0.1V dd , in which N β remains close to N/2 = 25. In the plot the contours of fillings slightly smaller (N β = 23) and larger (N β = 27) than N/2 are drawn as lines on the δ-Ω plane. Note that the N/2 plateau lobe is centred around a value slightly larger than δ = 1, which is approximately the point where the second term in equation (2) vanishes and our system becomes an effective Ising spin model.
The half filling plateau disappears with increasing Ω and we find a linear interpolation between regions with zero and full filling in this regime. This is consistent with a breakdown of the crystalline structure, and we note that in the limit Ω → ∞ for finite δ we expect the filling to be constant at N/2 since the state N k (|β k − |α k ) / √ 2 becomes the ground state of the system.
Note that the grid is too coarse-grained to resolve plateaus of smaller fractions of the filling. However, we observe several kinks in the particle number at small values of Ω/V dd . We show below, e.g., that the feature at δ ≈ 0.25V dd corresponds to a crystalline state with filling N/3. This will verify the lobe-like structure as depicted in figure 2, however we find that the plateau for filling N/3 is significantly smaller, at values of Ω ≈ 10 −2 V dd . For values of δ ∼ 2.25V dd we also observe some kinks, which correspond to regular patterns in the lattice with two subsequent sites being in the strongly interacting state, followed by one site being in the state |α k .
In order to illustrate the claim that the plateaux seem to belong to crystalline states, we can investigate the fluctuation in the total occupation number of the strongly interacting state, In figure 4 we show ∆N β for δ = 1.1V dd and 0.05V dd ≤ Ω ≤ 5V dd . We find that ∆N β decreases significantly with decreasing Ω and a crossover behaviour is centred around For Ω 0.2 a crystalline state with nearly no fluctuations ∆N β 3 N = 50 is present. For large Ω, ∆N β converges to a value close to N/4 = 12.5, which is the correct result for finite δ and Ω → ∞, i.e. where N k (|β k − |α k ) / √ 2 becomes the exact ground state of the system. Note that the quantity ∆N β is related to the Mandel Q-factor Q ≡ ∆N β /N β − 1, which has been measured in recent experiments for the number of Rydberg excitations in small ensembles of cold atoms [39] and quantifies the deviation from a Poissonian number statistic [40].
To study ground states in more detail for different parameter regimes, and demonstrate clearly the crystalline behaviour we evaluate DDC functions, defined as In a finite system with open boundary conditions, the exact values of G k will differ from the site i where the DDC is evaluated. However, the decay behaviour of the correlation functions for large k should be independent of the site index i in a sufficiently large system. To distinguish ground state phases by their decay behaviour, we will first calculate ground states of a system with periodic boundary conditions, in which the DDC becomes site independent. We then drop the site index and write G k .
In figure 5(a) we show a characteristic DDC in the parameter regime where the filling of the lattice is approximately N/2, choosing a detuning of δ = 1.1V dd , The insets show the modulus |G k | on a logarithmic scale. Here we show results for a system with 50 sites and periodic boundary conditions. We find states with a marked crystalline structure of periodicity 2 in panel (a) and a phase transition at 0.2V dd < Ω < 0.3V dd in terms of the decay behaviour of |G k |. In contrast to panel (a), in (b) |G k | decays exponentially. In the case of Ω = 0.4V dd the long-range decay behaviour (k 10) follows a power-law decay ∝ Ck −3 (with C ≈ 1.7 × 10 −2 ).
with both Ω = 0.1V dd and Ω = 0.2V dd . The inset in figure 5(a) shows the modulus of the correlation function, |G k | on a logarithmic scale. In the regime of the half filling plateau we indeed find crystalline structure with a pronounced zigzag pattern. The long-range crystalline order is indicated by the decay of |G k | to a non-zero value as a function of the site index k. We find that for increasing N the value of the constant to which |G k | decays becomes independent of the system size. We note, already from the comparison in figure 5(a) that as Ω/V dd increases the strength of the long range crystalline order decreases. If we continue to increase Ω/V dd , we find that this order is broken entirely, as shown in figure 5(b), with the same parameters, but for Rabi frequencies Ω = 0.3V dd and Ω = 0.4V dd . In these cases we also observe peaks with periodicity 2, however, these occur only for very small values of k 5. As k increases, |G k | decays to zero, indicating the breakdown of the longrange order. This appears as a phase transition from crystalline order to essentially a "paramagnetic phase", in which the correlations typically decay exponentially. Here we note a special feature of the correlation functions due to the long-range interactions, however, which is most clearly visible in the results for Ω = 0.4V dd in figure 5(b). There we observe that the long-range order has been replaced by exponential decay at short distances (k 10), but that this changes from exponential to a power law-decay. By fitting a function we find that in this region |G k | ≈ C/k 3 , with a constant C. This is the same power law decay as the interaction term, as is consistent with similar results derived for other longrange spin models [28,41]. Physically, the correlations initially decay exponentially due to the primary influence on a given particle being from its nearest neighbours (analogous to exponentially decaying correlations in the paramagnetic phase of an Ising model). However, once correlations arising from nearest neighbour interactions become sufficiently small, long-distance correlations will be dominated by effects of the direct interaction between distant particles. Thus, long range correlations are expected to decay at large distances with the same power as the interaction [41].
We find similar behaviour for crystalline states at other densities, including occupation periodicities larger than two. When the detuning is decreased from δ ≈ 0.5V dd it becomes less favourable to occupy |β k states in the system and therefore in the case of small Ω, crystalline phases with larger spacing between these occupations appear in order to reduce the dipole-dipole interaction energy. For example, in figure 6 we show the DDC indicating an occupation of the states |β k on every third site for Ω ∼ 10 −2 V dd and δ = 0.25V dd . To make this periodicity clearer, we present results for a system with N = 60 sites. As in figure 5, in figure 6 we observe the characteristic long range order in the DDC, combined with a marked peak on every third lattice site. Again we find a phase transition away from crystalline order, with an exponential decay of |G k | exhibited for larger values of Ω. Note that the crystalline phase with 1/3 filling occurs in a much smaller region of parameter space than the crystalline phase with 1/2 filling.
In what follows, we are interested in time-dependent creation of these crystalline states in the case of polar molecules or Rydberg atoms loaded into an optical lattice. These systems will typically occupy N ∼ 50 lattice sites, with open boundary conditions generated by the unoccupied ends of the chain. Thus, rather than attempting to derive the full phase diagram for an infinite system, we focus on the case of open boundary conditions relevant to the experiments. Although the correlation functions are modified by the boundary conditions, we show below that the characteristic signatures of different phases can still be clearly identified.
In figure 7 we show the DDC evaluated from site num- with open boundary conditions. We choose the same parameters as were shown in figures 5(a) and (b), including crystalline phases with periodicity 2. Despite obvious effects of the open boundary conditions we can clearly distinguish crystalline behaviour (Ω = 0.1V dd , 0.2V dd ) from the case were exponentially decaying density-density correlations are present (Ω = 0.3V dd , 0.4V dd ). While in the crystalline phase we find very little decay for increasing k until the effects of the boundary become apparent (k 20), in the case of Ω 0.3 the correlations rapidly decay below values of 10 −5 . As in figure 5 we still observe a long-range power law decay for k 10 proportional to k −3 .
In relation to figure 2, the results presented in figures 5-7 give evidence for the crystalline phases at filling factors 1/2 and 1/3 on the lattice, and show for each that there exist transitions as Ω is increased to non-crystalline phases. This justifies that the lobes are finite along the Ω k . The results shown are for the same parameter regime as in figure 5. Despite strong boundary effects, a crystalline phase is still distinguishable from a phase with such structure via the behaviour over short ranges k 10. In the crystalline phase the long-range decay follows a power law behaviour proportional to Ck −3 (C ≈ 1.7 × 10 −2 is a constant).
axis, and shows that the boundary depicted corresponds to a phase transition.
Note that the DDC functions can be measured experimentally both for Rydberg atoms and for polar molecules. Detection of these correlation functions is discussed in more detail in section VI.

V. DYNAMICAL CREATION OF CRYSTALLINE STATES
A. Time-dependent state preparation We now discuss the dynamical creation of crystalline states in our system. The natural initial state for the experiments corresponds to having all of the particles in the non-interacting states |α k , which would be the ground states for neutral atoms, or a single dressed rotational state for polar molecules. In terms of the effective model, this state is the ground state of the Hamiltonian (1) for large negative detuning and small Ω. The question is now whether one can time-dependently vary the laser parameters Ω and δ, i.e. find a trajectory in the Ω-δ-plane, so that this initial state is adiabatically transferred to a crystalline state of the form discussed in section IV.
In order to guide our discussion of this question, we first consider a small model system consisting of 8 sites, where exact diagonalization of the Hamiltonian (1) is straightforward. We can then estimate possible paths for the adiabatic transfer based on the size of the energy gap ∆E between the ground state and excited states. During the adiabatic ramp, we would like ∆E to remain as large as possible to suppress non-adiabatic transitions to excited states of the effective model. In figure 8 we show a shaded plot of the energy gap between the ground state and the first excited state as a function of Ω and δ, with lines of constant energy gap marked in the plot. Close to Ω = 0 we find a large region where the energy gap assumes small values and where one can recognize lobelike structures. Especially visible are lobes at values of δ ≈ 2.25V dd . As mentioned before these correspond to states with regular patterns of two subsquent sites being in the strongly interacting state. In contrast to the lobe for filling fraction N/3 of the strongly interacting state at δ ≈ 0.25V dd , these are enhanced due to the larger increase in interaction energy when adding single occupations of the strongly interacting state to the system. Therefore, these lobes are more robust when increasing Ω. In general, the system is too small to observe lobes of filling fractions N/3 and smaller and these are only visible as small kinks for δ ≈ 0.25V dd . With increasing Rabi frequency Ω the energy gap increases. The requirement of having at any instance of time a sufficiently large energy gap is most easily fulfilled if one chooses to "steer" a path in parameter space around those regions where the gap is the smallest. A simple path that seems to achieve this is, for example, to first increase Ω at constant large negative detuning δ and then increase δ at a constant large Ω. Both of these two steps can be achieved by changing the parameters δ and Ω rapidly, due to the large gap ∆E present at all times. Afterwards, Ω can be decreased to a small value on a line of constant δ. This path consisting of these three segments is marked on figure 8 by an orange arrow. Now we would like to test these insights by investigating the time-dependent propagation of the system along this path in parameter space. In order to investigate the scaling with the system size, we study systems of size 20, 30, 40 and 50 lattice sites. We note that the large sizes here are already typical of likely experimental setups. As an example, we choose to prepare a crystalline state with excitation periodicity 2, which is the ground state of the our model for the parameters Ω = 0.05V dd and δ = 1.1V dd . We perform real time simulations for a path consisting of three segments, as sketched in figure 8. We start at Ω = 0.05V dd and δ = −3V dd . There, the inner product between the ground state of the Hamiltonian and the state in which all atoms are in the |α k state is sufficiently large for all system sizes. The segments of the path are: 1. We increase Ω from 0.05V dd to 0.5V dd at δ = −3 with a variation rate ∆Ω/∆t = 0.05V 2 dd .
Due to the large gaps in the energy spectrum it is possible for segments 1 and 2 to obtain states with a high fidelity F ≡ | Ψ G (Ω, δ)|Ψ(t) | 2 ∼ 1 in the comparison of our evolved state, |Ψ(t) , with the ground state of the effective model with the corresponding parameters, |Ψ G (Ω, δ) on that path. The most sensitive part of this process is the third segment. In figure 9 we show the fidelity F of the adiabatically evolved state along segment 3 as a function of Ω(t) and for various variation rates ∆Ω/∆t. On N = 30 lattice sites, we find that for the three rates ∆Ω/∆t = 5 × 10 −2 V dd , 1 × 10 −2 V dd , and 5×10 −3 V dd , this fidelity drops below values of 80%, which occurs when crossing the phase transition from figure 5, i.e. in the regime 0.2V dd < Ω < 0.3V dd . At this point the ground state gap becomes too small and excitations to higher excited states analogous to Landau-Zener tunnelling processes occur. However when using a rate of ∆Ω/∆t = 1 × 10 −3 V dd we find that the fraction of the state that is lost into higher excitations remains reasonably small and the final fidelity with the crystalline state at Ω = 0.05V dd is larger than F ≈ 99%. Thus, we have shown that a crystalline state can be produced with very high fidelity in a finite system of 30 sites. The timescale which is required for all three segments is of the order of 500/V dd , and we will show that this is experimentally realisable for our two physical implementations below.
In an experiment with either polar molecules or neutral atoms coupled to excited Rydberg states in an optical lattice, the final state could be characterized by measuring DDC functions. These can be most easily detected directly by in situ imaging experiments [42]. In the case of polar molecules, the crystalline structure will be visible in noise-correlation measurements of state-selective momentum distributions of molecules released from the lattice [43][44][45]. In these correlation measurements, the periodicity of the crystalline structure would be present as a peak at the corresponding reciprocal lattice vector. In the case of Rydberg excitations, direct measurement of DDC functions could also be made by imaging the sample on a channel plate detector [20]. An important question is: to what extent the same fidelities we have observed for a system size of 30 lattice sites can also be achieved for larger system sizes. In general we expect the energy gap between ground and excited states to become smaller with increasing system size, and we expect the measure of the fidelity to become more sensitive due to the exponentially increasing size of the Hilbert space. To analyse how much this affects the final fidelity, in figure 10 we plot the fidelity achieved after ramping the laser parameters through all segments of our preparation procedure, as a function of the time required for the final segment 3. We consider system sizes of 20, 30, 40 and 50 sites.
We find that for all system sizes the fidelity is > 94% when choosing a variation rate slower than ∆Ω/∆t = 1 × 10 −3 V 2 dd , which corresponds to an operation time of approximately 450/V dd on segment 3. In the case of an even slower rate of ∆Ω/∆t = 5 × 10 −4 V 2 dd which requires a time of approximately 900/V dd , the fidelity becomes larger than 98%, even in the case of a system with 50 lattice sites. This fidelity is very high, especially given that the Hilbert space dimension is 2 50 in the latter case.

B. Comparison with experimental parameters
Let us compare the time required to reach a final crystalline phase with experimentally accessible timescales, which are in general limited by the finite lifetime of the strongly interacting state. In the scheme presented above, at the slowest variation rate of ∆Ω/∆t = 5 × 10 −4 V 2 dd on segment 3, a total operation time of the or- der 1000/V dd is required. Note that this time can be further dramatically decreased by choosing shorter paths or by not utilising a constant ramping rate but rather reducing it with time while approaching the critical point of the phase transition.
In the case of the polar molecule implementation (see section VI A for more details), typical parameters give V dd ≈ 10 kHz, assuming a permanent dipole moment d 0 = 1 Debye and a lattice spacing of a = 400 nm. This translates to a time for the complete preparation of approximately 0.1 s, which is significantly shorter than recently measured lifetimes of molecules in an optical lattice of 8 s [46]. By choosing polar molecules with larger dipole moments like LiCs (d 0 ≈ 5.5 Debye [24]), the preparation time can be further reduced down to the order of a few ms.
In the case of the Rydberg atom implementation (see section VI B for more details), if we choose a principal quantum number n = 14, then V dd ≈ 7 GHz (assuming a lattice spacing of a = 400 nm). Hence, the total operation time corresponds to approximately 140 ns in the experiment. This value is well below the typical lifetime of a n = 14 Rydberg level of Lithium, which can be found to be approximately 2.2 µs [47]. Note that these estimates neglect the effects of blackbody radiation, which can redistribute population amongst excited levels. For Rydberg atoms, in general the lifetime increases with increasing n and furthermore the required operation time decreases as n −4 . However, the disadvantage of large principal quantum numbers is the small level distances ∆ inter and ∆ intra (see figure 13), which both decrease with n. This limits the experimentally accessible values of δ. In the case n = 14 (a = 400 nm) at an electric field strength of F el = 100 kV/m, we can estimate by using the results known from the hydrogen atom ∆ inter ≈ 4.6 × 10 3 V dd and ∆ intra ≈ 90V dd . The field strength is thereby well below the Inglis-Teller limit, which is in this case 319 kV/m [48]. This justifies the validity of our two-level model, derived in section VI B, in this regime.
This comparison to the lifetime of a single Rydberg atom is oversimplified, as in the presence of N β excited atoms, the rate of single decay events will increase roughly proportional to N β [49]. However, as the system size (and hence N β ) becomes larger, the effect of a single decay on the crystalline structure should be reduced. Whilst the many-body state fidelity will be very sensitive to a single decay, DDC functions should be robust, at least on length scales smaller than those separating places where spontaneous emissions have occurred. The effect on the many-body state of single decay events also depends on when during the excitation process they occur. A full treatment of this excitation process including decay events could be achieved by studying the dynamics of a master equation including spontaneous emission events.
Therefore and from the results in figure 10 we conclude that the adiabatic passage should be experimentally realisable for both our implementations of the effective spinmodel with a size of ∼ 50 sites. For larger system sizes, longer timescales will be required to reach as high fidelities as we have obtained here. However, the properties of the final states will typically be determined in an experiment by measurements of the correlation functions. In contrast to the state fidelity, correlation functions (together with the physical character of the state) become more robust to small non-adiabaticities (especially those resulting in localized defects) in larger systems.

A. Polar molecules in an optical lattice
In the following we describe how to implement Hamiltonian (1) with a system of polar molecules. We consider N molecules confined to the lowest motional level on subsequent sites of a deep 1D optical lattice, with modefunction φ k (x) on site k. We assume that the lattice is sufficiently deep that tunnelling between wells of the lattice is strongly suppressed on typical experimental timescales. In particular we consider molecules in their electronic vibrational ground-state manifold with a closed shell structure 1 Σ, which possess permanent dipole moments (e.g. 40 K 87 Rb like in recent experiments [7]). The setup is depicted in figure 1. We use a strong and a weak microwave field linearly polarized along the z-axis and an additional static field aligned along the same direction. Given the fact that we consider samples of the order of a few to hundred µm one can neglect the spatial variation of the respective static and microwave field along the optical lattice axis, which is aligned in the e x -direction.
The low energy spectrum of a single molecule at site k in our setup is well described by rotational excitations, governed by a rigid rotor Hamiltonian coupled to the fields via an electric dipole interaction [50][51][52][53] Here, B is the rotational constant typically of the order of tens of GHz ( ≡ 1) [26],Ĵ k is the angular momentum operator,d z k the z-component of the dipole operator and E dc , E z 1 (t) = E z 1 exp(−iω 1 t) + c.c., and E z 2 (t) = E z 2 exp(−iω 2 t) + c.c. the static, weak and strong external microwave fields, respectively.
By applying electric fields one lifts the degeneracy of the state manifolds with J > 0 and induces a finite dipole moment in the dressed rotational states. Since the fields are aligned along the quantisation axis, m J is conserved, and we focus on states with m J = 0 in our setup, which are connected to the ground state. The induced dipole moments of two molecules at distinct sites k and l then lead to an energy shift due to the dipole-dipole interaction potential.
In order to describe this, we need to consider many molecules on the lattice. Molecules in an internal state denoted σ can then be represented by field operatorŝ Ψ σ (x), obeying the standard bosonic or fermionic commutation relations. We consider a situation where no motional modes other than φ k (x) for each site k are populated during the experiment, which is valid provided that the band separation ω T is substantially larger than the laser coupling parameters and interaction energies. Thus, we can expand the field operators aŝ and identify individual molecules occupying each mode described by the operatorsĝ σ k . Since all rotational states we consider here have m J = 0, they are only coupled by an interaction proportional tod z kd z l , which can be integrated over the modefunctions in each lattice site. The resulting many-body Hamiltonian of a system for N A strong linearly-polarized microwave field E z 2 (t) couples the states |1 dc and |2 dc the Rabi-frequency Ω2 and detuning ∆2. Fine-tuning the microwave-field (with respect to the static field) makes it possible to cancel the permanent dipolemoment in one of the two dressed states, |α+ and |α− , cf. dα → 0. A second (weaker) microwave field E z 1 (t) finally couples the ground state |β ≡ |0 dc to the dressed state with vanishing dipole |α with an effective Rabi-frequency Ω1 and detuning ∆1.

molecules can be written aŝ
with lattice spacing a, and where V dd ≈ 1/(8π 0 a 3 ). Note that we neglected corrections to the shape of the interactions, which are of the order O(l H /(a|k − l|)), where l H denotes the characteristic harmonic oscillator length for the optical potential well on each site.
A suitable choice of external electric fields then makes it possible to strongly couple/dress locally two (effective) rotational levels, such that only one level has a strong induced dipole moment, while at the same time the transition dipole moment remains negligibly small. This results in dipole-dipole interactions in equation (8) reducing to the diagonal form of equation. (1).
We now detail our appropriate choice of the electric fields to obtain the desired dressed states. The key idea is to apply the static field E dc to induce a finite dipole moment in the rotational states, leading to the dressed states |J dc , defined by [BĴ 2 k −d z k E dc ]|J dc = E J |J dc (see figure 11(b), and to then couple the two states |1 dc and |2 dc by a strong microwave field E z 2 (t), such that the dipole moment vanishes in one of the resulting dressed states. The field E z 2 (t) is tuned to the transition between the two states with Rabi frequency Ω 2 ≡ 2|d z k |1 dc E z 2 and detuning ∆ 2 ≡ ω 2 − (E 2 − E 1 ) (see figure 11). Diagonalizing the corresponding two-level Hamiltonian and making the rotating wave approximation leads to the two dressed states |α ± k = cos(θ/2)|1 dc − sin(θ/2)|2 dc with the dressing angle θ ≡ arctan(−2Ω 2 /∆ 2 ) (0 < θ < 2π), where the state |α − k (|α + k ) corresponds to θ < π (θ > π). The corresponding energies of the dressed states are E ± α = E 1 − ∆ 2 /2 ± ∆ 2 2 + 4Ω 2 2 /2. Finally, we couple the state |α k ≡ |α − k and the state |β k ≡ |0 dc by the weak electric field with Rabi frequency Ω 1 ≡ α|d z k |β k E z 1 and detuning ∆ 1 ≡ ω 1 − (E α − E 0 ). Thereby we assume that the field E z 1 (t) is much weaker than E z 2 (t), so that the states |α and |β are unaffected from E z 1 (t). Going to the rotating frame and making the rotating wave approximation the manybody Hamiltonian becomes time-independent and readŝ withσ x ≡ |β α| k + |α β| k andn k ≡ |β β| k . In the last term of equation (9),V d is the dipole-dipole interaction matrix in the rotating frame, which in the two-molecule basis {|α k |α l , |α k |β l , |β k |α l , |β k |β l } iŝ Here, the non-vanishing elements ofV d can be evaluated in the rotating wave approximation from the bare singleparticle dipole-operator elements d ij ≡ i|d z k |j dc : v αα αα = cos 2 (θ/2)d 11 + sin 2 (θ/2)d 22 2 + 2 sin 2 (θ/2) cos 2 (θ/2)d 2

12
(11) v αβ αβ = d 00 cos 2 (θ/2)d 11 + sin 2 (θ/2)d 22 (12) v βα αβ = cos 2 (θ/2)d 2 10 + sin 2 (θ/2)d 2 Our goal is to find optimal values of the microwave dressing angle θ and static field strength E dc so that the desired Hamiltonian [equation (1)] is realized with as small an error as possible. Though from equation (13) we see that there is no parameter choice where v βα αβ is exactly zero, we can find regimes where the model is implemented up to very small errors. At this point we note that subtracting a multiple of the identity from V d amounts to an overall shift of the zero of energy, and thus we can equivalently writê where δv ≡ v ββ ββ − v αβ αβ . In Figure 12 we plot the matrix elements δv, |v αα αα − v αβ αβ |, and v βα αβ as a function of the microwave field angle θ and the electric field strengths up to E dc = 20d 0 /B. In figure 12(a) we show the magnitude of the interaction energy for the strongly interacting state, δv. Figure 12(b), and figure 12(c) show the ratio of the other interaction matrix elements to δv, i.e. |v αα αα − v αβ αβ |/δv, and |v βα αβ |/δv. In panel (d) we plot the (normalized) total deviation from an exact implementation of the Hamiltonian [equation (1)]. This is quantified as δV d /δv, given by We find that for field strengths E dc ∼ 14B/d 0 , fine tuning the microwave field at θ ∼ 0.9π makes it possible to obtain deviations δV d ∼ 10 −3 δv, while the interaction strength for the strongly interacting state is δv ≈ 0.65d 2 0 (see contour line in panel (a) ). Thus within small errors we obtain the model interaction with a typical timescale of V dd ≡ 2δvV (1) dd ≈ δv/4π 0 a 3 ≈ 10 kHz for typical parameters (a = 400 nm, d 0 = 0.8 Debye). An alternative approach using circularly polarized coupling fields is also possible, which can give rise to an implementation of the model for appropriately tuned parameters [54].

B. Neutral atoms coupled to Rydberg states
In this section we show how to implement spin Hamiltonian (1) by making use of neutral atoms coupled to excited Rydberg states. We consider alkali atoms confined in an optical lattice in a single tube, with inter-site distance a oriented along the x-direction. Perpendicular to the lattice a homogeneous electric field, F = F el e z is applied, together with a laser which couples the ground state |g of each atom to a specific chosen Rydberg-Stark state.

Dipole-Dipole Interaction
While the atomic ground state is only slightly shifted by the external electric field, the level structure of highlying Rydberg states is strongly altered compared to the field-free case. Nevertheless, the azimuthal symmetry of the atomic Hamiltonian is preserved such that m (the quantum number of the z-component of the electronic orbital angular momentum L z ) remains a good quantum number [48]. States with large quantum defect, i.e. usually s-and p-states, sustain a second-order Stark shift while the degenerate states with higher angular momentum exhibit a linear Stark effect, reminiscent of the hydrogen atom [55]. The derivative of the energy with respect to the electric field strength directly yields the electric dipole moment of the given state. While the precise level structure and therefore the dipole moment of the Stark states depends on the element under consideration the exact results known from the hydrogen atom usually provide good estimates for the quantities of interest. Low angular momentum states with quantum defect exhibit a second order Stark effect. In our scheme we couple the ground state |g of an atom to the Rydberg state |dn , which exhibits the strongest Stark shift within the given n-manifold. The magnetic field strength and the laser parameters are set such that ∆intra ∆inter, ∆intra |δ| and ∆intra Ω0.
A sketch of a Rydberg-Stark level structure is shown in figure 13. For sufficiently small field strengths, i.e. below the Inglis-Teller limit where the atom becomes ionized, the levels are grouped in manifolds which can be labelled by the principal quantum number n. In our study we focus on the energetically highest state of such a manifold whose azimuthal quantum number is m = 0. Employing the formula for hydrogen, the electric dipole moment of this state, which we denote as |d n , evaluates to with e and a 0 being the charge of the electron and the Bohr radius, respectively. Two atoms excited to the state |d n k and |d n l at distinct sites of the lattice k and l, are expected to interact strongly via the dipole-dipole interaction potential, as already for n = 17 equation (17) gives a dipole moment of more than one kiloDebye. We intend to work in a regime in which the interaction can be treated within first order perturbation carried out in single particle product states. To this end it is important to note that the dipole-dipole interaction commutes with L z1 + L z2 , and therefore only couples atomic pair states for which the azimuthal quantum numbers obey m 1 + m 2 = m 1 + m 2 [8]. We assume that all of these couplings are far off-resonant, which can be assured by properly tuning ∆ intra and ∆ inter using the electric field. Moreover, accidental resonances can be avoided by restricting oneself to low principal quantum numbers n since here the two-particle energy spectrum possesses large gaps which make near-resonant excitations unlikely. Employing d n ⊥ e x the interaction Hamiltonian between two atoms being in the product state |d n k |d n l can then be written aŝ with V dd ≈ 1/(8π 0 a 3 ), when the confinement length in each lattice site is much smaller than the lattice spacing.
We now turn to the description of the interaction between the laser and the atom. We employ a two-level approximation assuming that the ground state |g k is coupled only to the state |d n k . This is justified as long as the modulus of the laser detuning δ and the Rabi frequency Ω 0 are much smaller than the energy gaps ∆ intra and ∆ inter . Making the rotating wave approximation, the interaction of the laser with an atom located on the k-th site can be written as ( = 1) with Ω 0 and δ being the Rabi frequency and the detuning, respectively.

Many-body model
We now consider the many-body physics when many atoms in an optical lattice are simultaneously coupled to the state |d n k . We assume that either atoms are loaded into a optical lattice with tight radial confinement, so that we have one atom per site along a 1D tube, or that we have weak radial confinement and many atoms per site. In each case, we can assume that the ground-state atoms are confined to a single modefunction at each given lattice site, φ k (x). In the case of a single atom per site, this is the Wannier function for the lowest Bloch band, and in the case of many (weakly-interacting) atoms it is the solution to the Gross-Pitaevskii equation in the mean-field limit.
In order to formulate the Hamiltonian of the manybody system, we introduce field operatorsΨ g (x) and Ψ r (x), corresponding respectively to atoms in the ground and Rydberg states, and obeying the standard bosonic (or fermionic) commutation relations. We assume that the timescale for excitation to the Rydberg state is much shorter than timescales corresponding to the atomic motion, and that therefore we can neglect population of any modefunction other than the φ k (x) for either the ground or Rydberg states. In this case, we can approximately expand the field operators in terms of mode operators corresponding to each of these modes, where both the initial field operators and the mode operators,ĝ † k andr † k , obey the appropriate bosonic commutation or fermionic anti-commutation relations, depending on which atoms are used in the experiment. We can then write the Hamiltonian for a lattice with N sites aŝ In the last term we have again neglected corrections to the shape of the interaction potential, which are of the order O(l H /(a|k − l|)), with l H the characteristic harmonic oscillator length for the deep optical potential well on each site.
The quantity V on dd is relevant only in the case that we initially have many atoms per site, and gives the energy offset due to the dipole-dipole interaction if two Rydberg atoms are excited on the same site. We will assume that V on dd V dd /|k − l| 3 for k = l and Ω 0 , (N δ) V on dd , so that the occupation of a single site by two Rydberg atoms is far off-resonant and therefore the excitation of such a state is highly improbable. This phenomenon, which is known as the Rydberg blockade, makes it possible to restrict the maximum number of Rydberg atoms per site to one, such thatn k can only assume the two eigenvalues 0 and 1. Secondly, to arrive at a Hamiltonian governing solely the dynamics of the Rydberg excitations, we eliminate the creation and annihilation operatorsĝ † k andĝ k in equation (22). In the case where we have an initial ground state with a single atom per site, this is simply a relabelling of the state with a single ground state atom per site as the vacuum state of the operators r k . In the case where the mode function at each site is occupied by the same macroscopic number N g 1, on the other hand, we make the replacementĝ k ,ĝ † k → N g , in which case we define the effective Rabi frequency as Ω = Ω 0 N g . The model we derive is exact under the stated assumptions when we have a fixed number of atoms per site (e.g., beginning deep in the Mott Insulator regime for bosons or in a band insulator for Fermions), which is the ideal experimental realization. An approximation to the model can also be obtained in a weakly interacting gas, but with some inhomogeneity in Ω for different lattice sites, both due to position-dependence of the mean density (due to a trapping potential), and also due to number fluctuations between sites. In production of crystalline states via the adiabatic ramping process discussed below, we expect that the ramp is somewhat robust against these fluctuations, especially for high filling factor crystalline states. However, this depends on the details of the realization and the inhomogeneity in Ω.
Thus, we arrive at the spin Hamiltonian describing excitation of Rydberg states, where the Pauli matrixσ x k ≡r † k +r k has been defined. This Hamiltonian is equivalent to our model system from equation (1), with the states corresponding to |α k ≡ |g k and |β k ≡ |d n k .
The scheme presented in section V achieves a crystalline phase of Rydberg excitations on very short timescales. However, whilst the preparation can be very fast compared to the lifetime of the Rydberg states, this lifetime will pose an ultimate limit to the stability of the many-body state that is produced. One way around this is to transfer the excited Rydberg states to alternative stable ground states. However, an important issue arises in this case because of the momentum kicks that arise due to the interactions between excited atoms, which can create very significant excitations in the motional states of the atoms. This is especially an issue towards the end of the chain, where the atoms will experience a force primarily in the direction away from the centre of the chain. These kicks can be estimated classically from the force acting due to the dipole-dipole interaction potential between excited Rydberg state. For crystalline states one can estimate that the net momentum kick during the crystal preparation time decreases as k R ∝ 1/(f i) 4 , where i denotes the i-th excited Rydberg atom from the boundary and f the periodicity of the crystal, i.e. with N/f filling. Quantitatively, transitions to higher bound states are given by the size of matrix elements w n (x)w 0 (x)e ik R x dx, where w n (n > 0) denotes the higher band Wannier modefunctions. Thus the transition effects become small in the case that k R a 0 , where a 0 denotes the oscillator length, is small. The oscillator length is typically of the order of tens of nm and we can estimate the momentum kick of the boundary excitation for a N/2 filling crystal as roughly 1/nm. Thus, this will be a problem for atoms at the end of the chain in the N/2 case, but decreases rapidly for atoms away from the edge, or a larger periodicity in the crystal. One way to circumvent this problem for the case where all lattice sites are singly occupied would be to use induced dipoles rather than direct coupling to Rydberg levels from the outset, as discussed in reference [56]. The idea is to make use of two stable ground states, one of which is coupled off-resonantly to a Rydberg level, producing an induced dipole with a strength depending on the parameters of the coupling and the static dipole of the Rydberg level. In this way, the coupling between ground and Rydberg states in our model is replaced by the coupling between two stable ground states. This can be achieved, e.g., via a two-photon transition. This results in much smaller values of Ω and δ (and hence also the interaction strength between excited atoms), which could even be made comparable to or less than the band gap in the optical lattice. The atoms would then always remain trapped in the lattice. The disadvantage of this method would be substantially increased times for the total adiabatic process, however, the ratio of the preparation time to spontaneous emission lifetime of the state could always be made large.
This method would have the added benefit of a much smaller linewidth of the coupling field for the transition between interacting and non-interacting states. Especially for crystal periodicities other than two lattice sites, the narrow regions of detuning δ/V dd for which the phases are found implies a strong condition on the stability of the laser detuning in an experiment (and hence on the linewidth of the exciting laser). For example, a crystalline state with an excitation on every fourth site can be found in the region of Ω/V dd = 4 × 10 −3 for 0.06 < δ/V dd < 0.09. Using V dd ≈ 7 GHz from the previous example, this means that the range of detuning values in which states of the character is found is approximately 200 MHz, and we require that the laser linewidth is smaller than this value. This is realistic, however for crystalline phases with even larger spacings between the excitations, the requirements become even stronger. If we use dressed dipoles, then we would be using a twophoton transition to couple between the states, for which the detuning can be controlled very accurately.

VII. CONCLUSION AND OUTLOOK
In conclusion, using adiabatic sweeps of microwave or laser parameters, it is possible to form crystalline phases of strongly interacting states in either polar molecules or neutral atoms coupled to Rydberg states in an optical lattice. The effective model Hamiltonian exhibits a rich ground state phase diagram, with the crystalline order being destroyed as the coupling between states is increased. For typical experimental size scales and parameters we showed that preparation of crystalline structures can be achieved on reasonable timescales, and the structures can be detected by measurement of characteristic density-density correlation functions.
The results we have obtained here motivate further study of the phase diagram with its lobe-like structure of different crystalline phases (as indicated in figure 8). This also opens the opportunity to extend 1D time-dependent simulations to the study of crystal formation in regimes of significantly larger spacing than is studied here (e.g., using matrix product operator methods [29,30,[36][37][38]) and also to probe time-dependent transitions between the different phases. This will also allow connections to current work in which Rydberg excitations of randomly distributed atoms are studied [20].