Dynamical preparation of laser-excited anisotropic Rydberg crystals in 2D optical lattices

We describe the dynamical preparation of anisotropic crystalline phases obtained by laser-exciting ultracold Alkali atoms to Rydberg p-states where they interact via anisotropic van der Waals interactions. We develop a time- dependent variational mean field ansatz to model large, but finite two-dimensional systems in experimentally accessible parameter regimes, and we present numerical simulations to illustrate the dynamical formation of anisotropic Rydberg crystals.


Introduction
Highly excited Rydberg states of atoms have unique properties. This includes the size of the Rydberg orbitals scaling as n 2 , the polarizabilities as n 7 and a long lifetime as n 3 with n the principal quantum number. These properties are also manifest in interactions between Rydberg states, e.g. in van der Waals (vdW) interactions ∝ n r 11 6 , which can be controlled and tuned via external fields. Exciting ground state atoms with a laser to Rydberg states thus provides a means to study many body systems with strong, long-range interactions [1,2]. With the atomic ground state and the Rydberg state defining an effective spin-1 2, we can describe the many body dynamics in terms of a model of interacting spins [3][4][5], reflecting the competition between the laser excitation and vdW interactions, at least in the short time limit where the motion of the atoms can be neglected (the so-called frozen gas regime).
The availability of ultraviolet (UV) laser sources allows now to excite Rydberg p-states in a single photon transition [29] resulting in much larger Rabi frequencies compared to a three-photon process. However, in contrast to the well known s-states, the vdW interactions between two atoms in a Rydberg p-state are anisotropic [30,31]. The goal of this paper is to investigate the quantum phases and their dynamical preparation with a laser pulse protocol for these anisotropic interactions. We are interested in two dimensional (2D) systems with a relatively high density of excitations involving a larger number of atoms, and in particular in the dynamical formation of anisotropic Rydberg crystals. Our studies are performed within a time-dependent variational mean field ansatz, beyond what can be accessed by exact diagonalization (ED) techniques.

Model and method
that the atoms are trapped in a 2D square lattice with exactly one atom per lattice site, as obtained in a Mott insulator phase. The atoms are initially prepared in the ground state, denoted by|↓〉, and coherently excited by a laser to a Rydberg state| ↑ 〉 with Rabi frequency Ω and laser detuning Δ (see figure 1(a)). Two atoms i and k both excited to the Rydberg state|↑〉 and located at positions r i and r k , respectively, interact via vdW interactions 6 , where θ i k , is the angle between their relative vector − r r i k and the z direction of the lattice. These vdW interactions exceed typical ground state interactions of cold atoms by several orders of magnitude. We are interested in a situation where the vdW interaction has a non-trivial angular dependence θ C ( ) i k 6 , . Such an angular dependence arises, for example, in laser excitation to higher angular momentum states, e.g. to Rydbergp-states, as opposed to excitation of s-states where the interactions are isotropic [30]. In the remainder of this paper we will illustrate the anisotropic interactions by explicitly considering the stretched state = 〉 n P m | , 32 j 2 3 2 of Rubidium for which the θ C ( ) i k 6 , is dominated by a term proportional to θ sin i k 4 , [31]. Interactions are therefore much stronger along the x direction compared to the z direction (see figure 1(b)). The atomic physics underlying this interaction will be discussed in detail in section 3 below.
In its simplest form the dynamics of the driven Rydberg gas can be described by an interacting system of pseudospin-1/2 particles

∑ ∑∑
Ωσ Δσ correspond to the local Pauli matrices and = | ↑ 〉 〈 ↑ | P i i is the projection operator on the Rydberg level. We note that in this model atoms are assumed to be pinned to the lattice sites, which is referred to as the frozen gas approximation [2]. For isotropic interactions spin models of this type have been discussed in previous theoretical work [6][7][8][9][10][11], and have been the basis for interpreting experiments [13,14].
The modeling of the laser excited Rydberg gas as a coherent spin dynamics governed by the Hamiltonian (1) is valid for sufficiently short times, typically tens of μs. First, as we noted above, the model (1) ignores the motion of the atoms: laser excited Rydberg atoms are typically not trapped by the optical lattice for the ground state atoms, and there will be (large) mechanical forces associated with the vdW interactions. In addition, Rydberg states have a finite life time, scaling as τ ∼ n 3 τ ∼ n ( ) 5 for low (high) angular momentum states with n the principal quantum number, and black body radiation can drive transitions between different Rydberg states, further decreasing the lifetime by approximately a factor of two [32]. The regime of validity has been analyzed in detail in [33]: there the long time dynamics of laser excited Rydberg gas including motion and dissipation was treated, including the validity of the frozen gas approximation and the cross over regime. (a) Setup: the ground state atoms| ↓ 〉 are placed in a square optical lattice and are excited to a Rydberg state| ↑ 〉 via a homogeneous laser beam with Rabi frequency Ω and detuning Δ. The vdW interaction V between two Rydberg atoms i and k is a function of their relative distance − r r | | i k but also of the angle θ i k , between their relative vector and the magnetic field B which is set along the z direction of the lattice. The details of these interactions in the fine structure manifold n P We emphasize that the various quantum phases predicted by the spin model (1) as a function of the laser parameters and interactions, and their preparation in an experiment, can only be understood in a dynamical way. In an experiment all atoms are initially prepared in their atomic ground state, , which is the ground state of the many-body Hamiltonian (1) for Ω = 0 and Δ < 0. Preparation of the ground state of the spin Hamiltonian (1) for a given Ω and Δ can thus be understood in the sense of adiabatic state preparation, where starting from an initial time t 0 with laser parameters Ω = t ( ) 0 0 and Δ < t ( ) 0 0 we follow the evolution of the many body state for a parameter trajectory to the final time t f with figure 1(c). This dynamical preparation of many-body states and quantum phases of the spin-model (1) in a time-dependent mean field ansatz, in particular in the anisotropic case, will be a central question to be addressed below.
While our focus below will be on the anisotropic model, we find it worthwhile to summarize the basic properties and signatures of the quantum phases (ground states) of the spin model (1) for isotropic interactions. Even for this case, the ground-state phase diagram of the Hamiltonian (1) is rather complicated. In the so-called classical limit, Ω → 0, where all terms in the Hamiltonian (1) commute, the ground-state corresponds to the minimum energy configuration of classical charges on a square lattice interacting via a r 1 6 potential, and Δ serves as a chemical potential [34]. As noted above, for Δ < 0 this corresponds to the state 〉 G | with all atoms in the ground state. For Δ > 0 a finite density of excited Rydberg atoms is energetically favorable and the competition between the laser detuning and the vdW interactions results in a complex crystalline arrangement with a typical distance between excited atoms set by the length scale ℓ Δ ≡  C [ ( )] 6 1 6 . For one dimensional (1D) systems it is known that all possible commensurate crystalline phases form a complete devilʼs staircase [35]. By contrast, the problem of finding the classical ground states in 2D systems is much more complex, as the Rydberg atoms ideally want to form a triangular lattice to maximize their distance, which will compete with the square optical lattice for the setup of figure 1. Although an exact solution of this problem is not known in 2D, we expect a plethora of different crystalline phases in analogy to the 1D case. These different crystalline phases break the lattice symmetries in different ways and are stable over some parts of the phase diagram [34,36,37].
Away from the classical limit, Ω ≠ 0, the crystalline states of Rydberg atoms are expected to be stable for sufficiently small Ω. By increasing Ω quantum fluctuations will eventually melt the crystalline phases and we reach a quantum disordered phase. The nature of the corresponding quantum phase-transition has been studied in 1D [22,23] and remains an open issue in higher dimensions.
Concerning anisotropic interactions, it is natural to expect that the angular dependence of the vdW coefficient, θ C ( ) i k 6 , , is responsible for the presence of an anisotropic crystalline phase at small Ω. In summary, our goal below is to describe the dynamical formation of such crystals in large but finite systems similar to realistic experimental situations, where finite size effects still play an important role, and to compare the final state to the ground state of the system in order to assess the fidelity of the dynamical preparation. To this end, we developed an approach based on a time-dependent variational principle (TDVP) which proved very useful to describe the crystalline states for isotropic as well as anisotropic interactions with a large number of excitations, i.e. in a parameter regime where an exact solution cannot be applied.

Time dependent variational ansatz for many-body systems
In the following we present our variational ansatz and the corresponding equations of motion which we use to describe the dynamical preparation of Rydberg crystals governed by (1). The simplest variational ansatz which is able to describe crystalline states of Rydberg atoms takes the most general product state form where N denotes the number of atoms and the coefficients α i and β i obey the normalization condition . Crystalline phases correspond to states where the probability β | | i 2 to find an atom in the Rydberg state at lattice site i is spatially modulated and its Fourier components serve as an infinite set of order parameters. In contrast, in the quantum disordered phase the Rydberg density is homogeneous and β | | i 2 is equal on all lattice sites i.

Equilibrium properties of the variational ansatz
Before deriving equations of motion for the time dependent variational parameters α t ( ) i and β t ( ) i we discuss equilibrium properties of our variational ansatz. Note that the ansatz (2) captures the exact ground-and excited states of our model Hamiltonian (equation (1)) in the classical limit Ω → 0, where all eigenstates are product states. In the general case Ω > 0, it is an approximation and its validity will be discussed at the end of this subsection. In principle we find the variational ground-state by minimizing the expectation value of the Hamiltonian with respect to the variational parameters α i and β i . In the ground-state these parameters can be chosen to be real and the variational ground-state energy, Note that the number of equations is equal to the number of lattice sites and solving these equations is thus numerically feasible only for finite systems. For this reason we do not attempt to make predictions about the phase diagram of (1) in the thermodynamic limit and rather focus on experimentally relevant systems with a finite but large number of atoms instead. There is one notable exception, however: we can make a statement about the melting transition between the quantum disordered phase at large Ω and the adjacent crystalline phase within our variational (mean-field) approach. In the thermodynamic limit, the quantum disordered phase has a homogeneous Rydberg density = ≡ f p p : R i and we can determine at which point the homogeneous solution becomes unstable to density modulations. Linearizing the mean-field equation (4) in small perturbations around the homogeneous solution are Fourier components of the interaction potential (note that the density f R of Rydberg excitations depends on Ω and Δ). For monotonously decaying potentials some Fourier coefficientsV k are indeed negative and the minimum ofV k in equation (5) thus determines the phase boundary beyond which density modulations form.
We note that an expansion in small density modulations around the homogeneous solution implicitly assumes that the melting transition is continuous. It is possible that this transition could be first order, however. In order to rule out a discontinuous melting transition we minimized the variational energy (3) numerically on a lattice with N = 441 sites and found that the melting transition is indeed continuous.
The momentum k 0 at which the interaction potentialV k is minimal determines the wave-vector at which density modulations form in the crystalline phase. In the isotropic case this minimum is at , where a is the lattice constant of the optical lattice. If one approaches the crystalline phase from the quantum disordered phase, the leading instability is thus always towards a crystalline state with Neel-type order, which breaks a Z 2 lattice symmetry. Only at smaller Ω more complicated crystalline states appear, which are most likely separated by first order phase transitions 4 . Consequently, within our variational approach the quantum phase transition between the disordered and the crystalline phase is always in the Ising universality class for isotropic interactions, independent of Ω and Δ.
For an anisotropic interaction potential with angular dependence θ C ( ) i k 6 , , present for example between states of Rubidium as discussed in section 3, the minimum is at a wave-vector π = a k ( , 0) 0 . Again, if we approach the crystalline phase from the quantum-disordered regime, crystalline order will form only in x-direction with a period of two lattice spacings, whereas no crystalline order is present in z-direction. This transition is again continuous. The system thus decouples into an array of quasi 1D Rydberg gases. Upon further decreasing Ω, we expect a transition to a state with incommensurate, floating crystalline order in z-direction, in analogy to 1D systems [22,23]. The system remains long-range ordered in x-direction, however, and finally settles into a commensurate, genuinely 2D crystalline state at sufficiently small Ω. We leave a detailed investigation of this two-step directional melting transition open for future study.
The phase boundary obtained from (5) is shown in figure 2 for both isotropic as well as anisotropic interactions with the angular dependence represented in figure 1(b). As in the anisotropic case the interactions are much stronger in the x direction compared to the z-direction, the corresponding phase-boundary significantly differs from the isotropic curve and is very close to the one obtained for a 1D system. Note that the mean-field phase boundary has an unphysical re-entrance behavior at negative detunings. This is because our variational ansatz vastly overestimates the ground-state energy in the quantum disordered phase at finite Ω, where pair-correlations are important.
Indeed, our product ansatz of (2) does not describe correlations between local density fluctuations of Rydberg atoms for our variational wave-function. Deep in the crystalline phase these correlations are weak and decay exponentially with distance, however. We can give an upper bound on the strength of such density density correlations and consequently make a statement about the validity of our ansatz by estimating the strength of local density fluctuations Accordingly, density-density correlations are negligible deep in the crystal, where〈 〉 P i is either close to zero or one. As a consequence we expect our ansatz to be valid also at finite Ω as long as we are deep in the crystalline phase. At sufficiently large Ω, where the system enters the quantum disordered phase and quantum correlations become predominant, our ansatz is not a good approximation for the exact ground-state wave function.

Time-dependent variational ansatz and Euler-Lagrange equations
One of the main goals of our paper is to describe the dynamical formation of crystalline states of Rydberg excitations in a large but finite system during a slow change of the laser parameters. For this reason we incorporate our ansatz into a time-dependent variational approach, where the formation of Rydberg excitations is described by the time evolution of the variational coefficients α t , we use the TDVP [39] to derive the equations of motion for the variational coefficients during a slow change of Ω and Δ as in typical dynamical state preparation schemes [8,9,14]. The TDVP states that the time-evolution of the variational coefficients satisfy the least action principle which means that they can be derived using the Euler-Lagrange equations: where L is the Lagrangian atoms (see figure 1(c)). The red dashed line shows the phase boundary for isotropic interactions. For comparison, the black dotted line represents the mean-field phase boundary of a one-dimensional system. leading to a set of N 2 nonlinear coupled equations We note that these equations conserve the norm of the wavefunction for all times, i.e. α β . If Ω and Δ do not evolve in time, the expectation value of the energy Φ Φ = 〈 〉 E H | | is also a conserved quantity. However, this is not the case for a dynamical state preparation and the final energy depends crucially on the parameter trajectory.
In a perfectly adiabatic situation we would obtain the variational ground-state for Ω t ( ) f and Δ t ( ) f at the end of the time evolution. Since our sweep protocols are limited to timescales smaller than the lifetime of the Rydberg state, the preparation will not be perfectly adiabatic and we discuss deviations from adiabaticity in section 4.1. We note that given the phase diagram shown in figure 2 and the typical parameter sweep we consider ( figure 1(b)), the system has to undergo a quantum phase transition from the quantum disordered phase to the crystalline phase at some point in the preparation which may reduce the adiabaticity of the preparation significantly. In the finite systems that we consider in the rest of this work, a finite-size gap is always present which reduces this problem, however.
Finally, we note that our ansatz (2) is particularly suited to study the experimentally relevant situation where the Rydberg laser is switched off at the end of the parameter sweep Ω = t ( ) 0 f , because it captures the exact ground-and excited states of the Hamiltonian (1) in the classical limit Ω = 0, as discussed above. In section 4.1 we also estimate the typical Rabi frequency Ω at which our ansatz fails by comparing our approach to the exact solution of the Schrödinger equation.
In the next section, we explain in detail the implementation of the model Hamiltonian (1)  Rydberg states in order to provide realistic parameters for our numerical section 4.

Anisotropic interactions for Rydberg atoms in p-states
In the following we discuss the derivation of our model Hamiltonian (equation (1)) from a microscopic Hamiltonian Rydberg state. We first focus on the Rydberg manifold and their interactions and then discuss the laser excitations. The first term of (9) MHz/G the Bohr magneton and g j the Lande factor for = j 3 2. Note, that the quantization axis of the corresponding eigenstates is given by the direction of the magnetic field, B, and is aligned in plane along the z-axis, see figure 1. In order to neglect higher order shifts and to prevent mixing between different fine-structure manifolds the energy shifts Δ μ = E gBm m B j z j j have to be much smaller than the fine-structure splitting − E E nP nP 3 2 1 2 of the Rydberg manifolds. Typically, the fine structure splitting is of the order of tens of GHz, e.g. 7.9 GHz for n = 25.
Away from Foerster resonances two laser excited Rydberg atoms dominantly interact via vdW interactions [1,2]. In general, these vdW interactions,V vdW , will mix different Zeeman sublevels 〉 m | j in the nP 3 vdW dd , dd whereV vdW is understood as an operator acting in the manifold of Zeeman sublevels We note, that in the absence of an external magnetic field, → B 0, the new eigenenergies obtained from diagonalizingV vdW are isotropic. Anisotropic vdW interactions can be obtained by lifting the degeneracy between the Zeeman sublevels e.g. with a magnetic field. For distances large enough, such that the off-diagonal vdW coupling matrix elements of (11) are much smaller than the energy splitting between the Zeeman sublevels, it is possible to simply consider interactions between = 〉 nP m | , 32 j 3 2 states and neglect transitions to different m j levels. Typically, for Rydberg p-states the off-diagonal vdW matrix elements are of the same order as the diagonal interaction matrix elements. Pairwise interactions between two atoms excited to the| ↑ 〉 = = 〉 nP m | , 3 2 j 3 3 state are then described by the Hamiltonian 6 the vdW interaction potential, giving rise to the second term of (1). Here, θ = ∢ − B r r ( , ) is the angle between the relative vector and the quantization axis given by the magnetic field direction B (see figure 1).
The second term of (9), H L i ( ) , accounts for the laser excitation of Rubidium 87 Rb atoms prepared in their electronic ground state, which we choose as| Rydberg states. This can be done using a single-photon transition with Rabi frequency Ω, scaling as Ω ∼ − n 3 2 . Using a UV laser source it is possible to obtain Rabi frequencies of several MHz in order to excite Rydberg states around ∼ n 30 with a large density of excitations. The single particle Hamiltonian governing the laser excitation of atom i is is the laser frequency detuned by Δ from the atomic transition (including the magnetic field splitting) and k L is the wave vector of the laser. Unwanted couplings to Zeeman levels with Thus, π μ = C h ( 2) 6.35 MHz m 6 6 and μ = C h (0) 0.269 MHz m 6 6 . The dominant θ ∼sin 4 term arises from dipole-dipole transitions to nS 1 2 states, while residual interactions at θ = 0 and deviations from the θ ∼sin 4 shape originate from couplings to D-channels as discussed in [30]. The lifetime of the Rydberg state is τ μ ≈ 29 s [32]. We emphasize that although we only consider the example of a specific, strongly anisotropic, P 3 2 Rydberg state in this work, our formalism can also be applied to describe the crystals obtained for other Rydberg states. Moreover, considering in particular other values of n or P , 1 2 D 3 2 states, the function θ C ( ) 6 takes the same form as in equation (14) [30], with different coefficients.
We finally note that it is also possible to switch from the anisotropic configuration defined above to an isotropic configuration where the angle θ is fixed to a constant value π 2 . In this configuration, the magnetic field is rotated from the z to the y-direction (see figure 1)

Numerical results
In the following we present our numerical results obtained by propagating the equations of motion (8) along different parameter trajectories Ω Δ t t ( ( ), ( ))in order to describe the dynamical state preparation of Rydberg crystalline phases. To this end we first estimate the domain of validity of our approach based on the TDVP comparing in the case of small systems our results to the ED solution which is obtained from the Schrödinger equation. We then present the results of our approach for large systems ( > N 500) with large densities of Rydberg excitations where ED techniques cannot be applied.

Validity of the variational approach for small systems
It is instructive to start by considering small systems where an exact numerical solution is available which allows us to estimate the validity of our approach. The first situation we have in mind is the classical limit (Ω = 0) of the model (1) for a 1D system where the number of excitations takes the form of the stair case as a function of the number of atoms N (obtained by varying the chain length = − L N a ( 1) , for a fixed lattice spacing a) [35]. As its existence was recently demonstrated experimentally [14] we consider as a first illustration of our approach the same parameters as in [14] with the notable exception that we choose a larger sweep time = t 32 μs in order to describe the dynamical preparation of states which are as close as possible to the ground state of the system. We now test our approach by comparing the exact number of excitations n e for theses parameters, obtained after a truncation of the Hilbert space [8,12,27,42], to the one corresponding to our variational ansatz. The result is shown in figure 3 as a function of the number of atoms N where the laser sweep is represented in the inset. Our approach describes very well the excitation stair case and apart from some defects (for example for N = 8) compares very well with the exact solution (red crosses). We finally emphasize that as the sweep time is increased, our solution converges towards the exact classical ground state, as expected.
Our approach describes the key feature of the 1D system but as it relies on a mean-field approximation of the many-body Hamiltonian (1), its domain of validity may strongly depend on the dimension of the system. As a second illustration of our variational approach we consider, therefore, a small 2D system of N = 16 atoms, in an isotropic configuration where an ED solution based on the truncation of the Hilbert space is still available.
Our goal here is to show the influence of the final Rabi frequency Ω t ( ) f on the validity of our approach. To this end, we consider three different sweeps of the Rabi-frequency Ω and detuning Δ along the paths shown in figure 4 corresponding to three final Rabi frequencies Ω π = t ( ) (2 ) 0, 1, 4 f MHz at a positive detuning (2 ) and Δ π (2 ) as a function of time t used for the dynamical state preparation.  Figure 5 shows final distributions of Rydberg atoms computed using our TDVP approach (upper panel) as well as ED (lower panel) for the three sweeps (a) (b) (c). Note that in contrast to the 1D case, the variational ansatz describes the sweep to the classical limit Ω = t ( ) 0 f perfectly well, even though our approach propagates the wave-function through the non-classical region Ω > 0, where our ansatz is not strictly valid. In a perfectly adiabatic situation corresponding to → ∞ t f , the final state obtained with our time-dependent variational approach would coincide with the variational ground state, which is the exact ground state in this case. The fact that our results for a finite sweep time compare very well with the exact solution suggests that deviations from adiabaticity are negligible. For such a small system size, the competition between the laser excitation and the vdW interactions results in a regular pattern of four Rydberg atoms placed at the corners of the system. We also obtain a good agreement for Ω π = t ( ) (2 ) 1 f MHz and the corresponding pattern is not modified compared to the classical limit. However, for Ω π = t ( ) (2 ) 4 f MHz, our ansatz overestimates the ground state energy considerably, even though the distribution of excitations looks similar to the exact result. It is also instructive to study how the system behaves during the sweep. Figures 6(a) and (b) show a comparison of Rydberg density f R and energy as function of time during the parameter sweep for the three sweep protocols shown in figure 4. Again a substantial difference between TDVP and ED is only visible for sweeps to large final Rabi frequencies.
The results shown in figures 5 and 6(b), (c) allow to assess the validity of our approach for a realistic dynamical state preparation. We are also interested in estimating the typical value of the parameters Ω, Δ where our ansatz can describe the ground state of the model (1), regardless of the details of the dynamical state preparation. To this end we consider a very large sweep time = t 150 f μs to ensure that the equations of motion (8) lead to the formation of the variational ground state whereas the solution obtained with ED results in the exact ground state. We then estimate the regime of validity of the variational ground state as follows: we compute the Rydberg density at the end of the parameter sweep using both TDVP and ED and plot the relative difference of f R between the two approaches as a function of the final Rabi frequency Ω t ( ) f at the end of the sweep. Results are shown in figure 6(c) for four different final detunings Δ. We see that for a final detuning Δ π = t ( ) (2 ) 2 f MHz the difference Δf f R R starts to deviate from zero if the sweep protocol samples Rabi frequencies which are larger than Ω π > t ( ) (2 ) 1 f MHz. Accordingly, for Δ π = (2 ) 2 MHz our variational ansatz is correct as long as Ω π ≲ (2 ) 1 MHz. We note, however, that this criterion was obtained for a small system and potentially depends on system size.

Isotropic Rydberg crystals
Now that we have assessed the regime of validity of our ansatz for a small system and checked in particular that it can quantitatively describe the dynamical preparation of Rydberg crystals in small 2D systems, let us now present our results for large system sizes where an exact numerical treatment is not possible.
We first describe the formation of Rydberg crystals in an isotropic configuration. In analogy to the experimental setup [14], we start from a circular (cookie shaped) distribution of N = 777 ground-state atoms considering the three sweeps path shown in figure 4 keeping the other parameters of the last subsection unchanged.
Let us first comment on our choice of sweep paths (figure 4). In order to prepare a state which is as close as possible to the variational ground-state, it is particularly important to circumvent the region around the critical point Ω Δ = = 0, 0 [7] during the sweep into the crystalline phase. Indeed we found that the energy of the final state increases substantially if the initial negative detuning is chosen too small. On the other hand, if the initial detuning is too large, the length of the sweep path in parameter space is long and the rate of  commensurability issues arise, however, in particular if the average distance between excitations is on the order of a few lattice spacings. We observe that the crystalline structure is strongly pinned by boundary effects in our case and the crystalline structures in the classical limit thus do not resemble those which supposedly exist in the thermodynamic limit [34].

Anisotropic Rydberg crystals
We now describe the preparation of anisotropic Rydberg crystals. In this case, the magnetic field is set along the z-direction of the optical lattice, as shown in figure 1(a), keeping the other parameters such as sweep paths, atom distribution and the Rydberg level unchanged. Accordingly, the interaction between Rydberg atoms is anisotropic and stronger in x-than in z-direction.
Results for the three sweep paths are presented in the upper panel of figure 9. As in the isotropic case, the crystal progressively melts as Ω is increased. Note, however, that the anisotropy is visible in all cases and the crystalline structure melts first in the weakly interacting z-direction while translational symmetry is still broken in the strongly interacting x-direction. Again, we observe that the form of the Rydberg crystal is strongly pinned by boundary effects, similar to the isotropic case.
In the classical limit Ω = 0, we find an anisotropic crystal with an average distance between excitations on the order of Δ ≈ π ⎡ ⎣ ⎤ ⎦  ( ) in the z-direction. Again, the results for the sweep to the classical limit Ω = t ( ) 0 f indicate that the preparation was not perfectly adiabatic. Indeed, the excitation probabilities differ from 0 or 1 at the end of the parameter sweep, as in the case of isotropic interactions. The deviations from adiabaticity are even more pronounced for anisotropic interactions, as we find non-classical excitation probabilities on the order of ∼0.5 in this case. This can be attributed to the fact that the excitation gaps are smaller compared to the case of isotropic interactions, due to the substantially weaker interaction in z-direction.
In order to estimate the fidelity of the dynamic state preparation scheme we plot the distribution of Rydberg atoms obtained after a direct minimization of the ground-state energy within our variational ansatz in the lower panel of figure 9. Comparing this to the distributions obtained after the parameter sweep it is apparent that a number of defects are created due to the not fully adiabatic sweep protocol. Again, the crystalline arrangement is not strongly affected by the rather short sweep time, however.

Conclusion and outlook
In the present work we have developed a time-dependent mean field theory to model the dynamical preparation of anisotropic Rydberg crystals with atoms in 2D optical lattices. In addition we have presented results of numerical simulations relevant for experimentally realistic system sizes, in the limit of patterns with a large number of Rydberg excitations.
We note that the anisotropic character of the vdW interactions has been seen experimentally in a recent Rydberg-blockade experiment involving Rydberg s and d-states [44]. In contrast to the present work, where we considered the anisotropic vdW interactions between single Zeeman levels of the Rydberg states, i.e. in the limit of large Zeeman splitting, in these experiments vdW couplings involving transitions between Zeeman levels can be important. This interplay leads to several new physical phenomena, which will be presented in a future work.