Quantum simulation of a 2D quasicrystal with cold atoms

We describe a way to obtain a two-dimensional quasiperiodic tiling with eight-fold symmetry using cold atoms. A series of such optical tilings, related by scale transformations, is obtained for a series of specific values of the chemical potential of the atoms. A theoretical model for the optical system is described and compared with that of the well-known cut-and-project method for the Ammann-Beenker tiling. This type of cold atom structure should allow the simulation of several important lattice models for interacting quantum particles and spins in quasicrystals.


Introduction
With the discovery of the first quasicrystals [1], the quest began for, on the one hand, new quasiperiodic systems with better characterization of structural properties, and on the other hand, for theoretical methods to handle these systems. One of the goals of experiment has been, in particular, obtaining a single component quasicrystal, in the hope of finding direct relationships between its physical and geometrical properties. This may, we hope, become possible in cold atom systems. Cold atoms in optical lattices have been used to simulate quantum behavior of periodic crystals but not, thus far, of quasiperiodic tilings. Cold gases, of cesium, rubidium or potassium atoms for example, are used as quantum simulators for a great variety of systems. As compared to real solid state systems, cold atom systems represent ideal systems, in which model parameters can be tuned at will. Improvements in experimental techniques has resulted in an explosion of experimental simulations of condensed matter physics, quantum information and quantum optics models. Ultracold quantum gases in optical potentials therefore provide an exciting possibility towards the goal of synthesizing a perfect one-component quasicrystal.
The behavior of electrons in quasicrystals remains insufficiently understood, especially in dimensions higher than one. Numerical investigations have thus far been limited by the amount of computational time needed. Much effort has been devoted to trying to understand tight-binding models on simple tilings, as a first step towards the description of more realistic systems. It is thus natural to ask what possibilities exist for realizing a quasiperiodic tiling by trapping cold atoms in an optical potential. The advantages of such a system, if realized, are manifold. First, it would become possible to fabricate samples using a single atomic species -a significant simplification compared to real quasicrystals. It would be possible to directly observe the quantum states of the atoms, described by a tight-binding Hamiltonian. Finally and importantly, many of the parameters of the hopping and interaction Hamiltonians could be tuned. Many body properties in quasiperiodic systems could be properly studied under controlled conditions. The amount of disorder could be tuned as desired, Submitted to Crystals, pages 1 -12 www.mdpi.com/journal/crystals to mimic experimental situations. An experimental set-up to realize a two dimensional tiling was proposed in [2,3]. It was shown that one can obtain an eight-fold tiling bearing a close relation to the Ammann-Beenker tiling [4]. We will present the experimental system, the theoretical model and explain how a perfect quasiperiodic tiling can be obtained by introducing additional small interactions between atoms. This paper begins with a description of the experimental set-up. We then introduce a 4D description of the optical tiling which is naturally suggested by the experimental geometry. We next show the relation between the optical tilings and the perfect Ammann Beenker tiling, and describe how to transform the former into the other. In conclusion, some directions for simulating important theoretical models are suggested.

Experimental set-up
Optical trapping of atoms by laser generated potentials has led to the artificial realization in experiments of many different kinds of models on lattices, in which particles can move and interact via experimentally controlled interactions [5]. Atoms can be trapped by laser light thanks to the dipole force acting on an atom due to the Stark effect in an off-resonance electric field. Under suitable assumptions concerning the decay rate Γ of the excited state (for more details see the review ), it can be shown that the potential energy landscape seen by the atom has the form where v 0 is a constant and I( r) is the average value of the intensity at the position r. Denoting the electric field by E, I( r) = | E( r)| 2 . The sign of the prefactor, v 0 , depends on the polarizability of the atom, and this can be positive or negative depending on the frequency of the laser, ω L , relative to the resonance frequency ω 0 for the atom [6]. The constant v 0 depends, among other parameters, on the laser detuning and is positive for blue-detuning (ω L > ω 0 ) and negative for red-detuning (ω L < ω 0 ). In other words, atoms experience a net force directed towards nodes (antinodes) of the laser intensity pattern if the detuning parameter ∆ = (ω L − ω 0 ) is positive (resp. negative). In this paper we study this latter case, corresponding to atoms being attracted to maxima of the laser intensity pattern I. We note that, for large detuning, there will appear corrections to this interaction potential due to non-conservative processes and we will neglect these. Using this type of optical potential it has been possible to simulate models defined on a large variety of periodic structures including the square, triangular, kagome and other lattices. We will now consider a quasiperiodic structure with eight-fold symmetry obtained by this method. We consider a region where standing waves have been set up using four laser beams oriented at 45 • angles in the xy plane, as shown in Fig.1. The wavelengths, λ, are the same for all the beams, as are their amplitudes. We will consider the situation where all the polarizations are perpendicular to this plane, allowing the amplitudes to sum up to an absolute maximum. The alternative choice, of using in-plane polarizations yields smaller maxima of amplitudes, smaller contrast, and would therefore be less efficient in trapping particles. For the case of four standing waves, the intensity is given by where r = (x, y) is the position vector of the points lying in the plane of the lasers. The four wave vectors are given by with n = 1, .., 4, where k = 2π/λ. Notice that the four beams can have arbitrary different phase-shifts φ n . As long as the relative phase shifts are maintained at some fixed arbitrary values, these phase shifts do not change the nature of the structures obtained, as discussed later. The function I( r) of Eq.1 is quasiperiodic since there is no integer relationship between the four wave vectors k n . The intensity landscape obtained for a random choice of the phases φ n has a complex structure of maxima, minima and saddlepoints. The intensities of the peaks, or local maxima, have a range of values with the upper bound I max = 16I 0 . The potential energy thus has a minimum value of V min = −16V 0 with V 0 = |v 0 |I 0 .
The cold atoms in this region are attracted to local maxima of I(r) -i.e. the local minima of the potential energy. If one now introduces a finite density of atoms, depending on the chemical potential, only maxima corresponding to intensities bigger than a cut-off I c will be occupied by an atom. We will neglect fluctuations due to finite temperature and trapping of atoms in metastable configurations, and instead focus on the ideal structures one expects to find, as a function of I c .
Figs.2 show the optical tilings obtained for three representative values of the intensity cut-off: I c /I 0 = 10.8, 15 and 15.82. The edge-length of the tiles can be seen to increase by a discrete scale factor, the irrational number α = 1 + √ 2, also called the silver mean. The tilings are shown superposed on the intensity profile, represented by a shaded plot (dark shades for small intensity). These are examples of the type of structure that we will refer to as the optical quasicrystal (OT). As can be seen from the shape of the tiles, it is closely related to the standard Ammann-Beenker or octagonal tiling (ABT) [4,7] composed of squares and 45 • rhombuses. For larger and larger values of the cutoff, approaching 16I 0 , only the largest maxima will be occupied, and the atomic density in the xy plane correspondingly decreases. As the cutoff takes on successive special values in the series given in the next section, the lattice vectors increase by powers of the irrational number α.

A four dimensional model for the optical quasicrystal
This section summarizes the model of the OT which was introduced in [2,3], along with some additional details, and clarifications. Our aim here will be to give a mathematical description of the optical tilings shown in Fig.2, and calculate quantities such as the lattice vectors in terms of the cutoff intensity I c . In this model, the only dimensionful parameters are λ = 2π/k the wavelength of the laser light, which sets the length-scale, and V 0 , which sets the energy scale in the problem. All other relevant lengths and energies can then be specified in terms of these two quantities.
Let us return to the intensity of the laser beams given by the function Eq.1. Since the four vectors k n have no rational relationships I( r) is a quasiperiodic function in the sense of H. Bohr [8] and A.S. Besicovitch [9]. The Fourier transform is readily obtained by expanding cosines so that the spectrum is the finite set of all combinations ± k n ± k m . The 8-fold symmetry of I( r) follows from the remark that a rotation γ (see [3] Fig. 3). I( r) is also invariant by the 2-fold symmetry σ w.r.t. the x-axis, which maps { k 1 , crystallographic in 2D, it is crystallographic in the 4D space R 4 : Lifting γ and σ to R 4 gives the integer matrices in the standard basis { ε 1 , ε 2 , ε 3 , ε 4 }. They satisfy Γ 8 = Σ 2 = I and ΣΓΣ −1 = Γ −1 . This 4D representation of C 8v is reducible since R 4 is the Cartesian product of two planes P and P which are orthogonal and invariant. While the restriction of Γ to P is the previous γ rotation, the restriction to P is a rotation γ of 3π/4. One can choose orthonormal bases { e x , e y } in P (the "physical space" ), and { e x , e y } in P (the "perpendicular space") so that the orthogonal projections e n = π( ε n ) and e n = π ( ε n ), all of norm 1/ √ 2, are as shown in Fig. 3.
e y , e x , e y } bases of R 4 . The transformation is given by the following rotation R: The wave vectors k n are the projections in P of orthogonal 4D vectors K n = K ε n = ( k n , k n ), so that their magnitude is K = √ 2 k. The intensity Eq. 1 can be obtained from the 4D periodic function so that I( r) = I( r, 0) is the restriction of I( R) to the plane P.
The maxima of I( R) occur on points where all cosines equal 1, i.e. on the cubic lattice 2π K Z 4 , and on points where all cosines equal −1, which correspond to the body centers. Let B denote the BCC lattice of Z 4 : coordinates are either all integers or all half integers. Accordingly, the maxima 16 I 0 of I( R) are located on the BCC lattice 2π The vertices of B can be written as R = T x = ∑ x n β n with x = (x 1 , x 2 , x 3 , x 4 ) (where the x n are integers) and where T is the matrix A basis of B is given by the unit vectors β n = T ε n (n = 1, .., 4). One can check that det(T) = 1 2 and that T.C 8v = C 8v .T so that B is invariant with respect to C 8v and belongs to the same Bravais class as the hypercubic lattice Z 4 . Furthermore, T commutes with the 8-fold generator Γ. The four primitive lattice vectors β n of B project onto the set b n in P, and b n in P , as shown in Fig. 3. The norms of these vectors are b n = 1 If the cutoff I c is high enough, the last condition is fulfilled in a set of disjoint domains centered on the BCC lattice √ 2π k B. If the cutoff I c is close to the absolute maximum I m = 16 I 0 , one can substitute the quadratic approximation I( R) ≈ (16 − 8k 2 ( r 2 + r 2 ))I 0 . The domains are close to spheres of radius ρ given by 8k 2 ρ 2 = (16I 0 − I c )/I 0 . Their projections on P and P are close to disks of the same radius.
If I( r) > I c is a local maximum, the point ( r, 0) belongs to a domain centered on a vertex R of the BCC lattice, at a distance bounded by ρ in the quadratic approximation. These points r are close to the BCC lattice, allowing for a comparison with the cut-and-project method as discussed below.
The vertices of the octagonal tiling [7] are the projections r of 4D lattice points R = ( r, r ) such that r belongs to an octagonal window generated by the four vectors b n (see [10] for the cut and project algorithm). The area of this window (see Fig.4 Figure 3. Projections of the bases {ε n } of Z 4 and {β n } of of the BCC lattice in P (left) and P (right). the octagonal tiling enlarges the edge length of the tiles in P by a factor α. while distances in P are reduced by the same factor. Inflated octagonal tilings correspond to selection windows of area W/α 2p with p = 0, 1, 2, ... If V c is low enough and if the areas of D and W are equal (up to inflation), a relationship between both structures can be expected. This equality ensures that the two structures have the same density of points. Using the quadratic approximation of I( R) this condition writes The OT of Fig.2 correspond to p = 1, 2 and 3. The circular windows D are shown in Fig.4 for p = 0 and p = 1, inside the window of the octagonal tiling. The criterion in Eq.5 can be further simplified by noting that, for p > 1, the cosine terms in V can be expanded to second order in the distance ρ. In this limit the selection windows for the OT are approximately circular and the radius for the cutoff value I c is given by k 2 ρ 2 c = √ 2πα −2p . It can be shown that the edges of the tiles of the pth OT have length = √ 2π k |b n |α p . The smallest edge length ≈ 3.81λ is obtained for p = 1.

Fourier transform of the OT
We now turn to the diffraction pattern of a structure obtained when atoms occupy all of the allowed vertices corresponding to a particular value of I c . The general method for finding the Fourier transform and indexation of peaks is discussed for example by A. Katz and D. Gratias in [11]. The structure factor consists of Bragg peaks whose positions are given by the reciprocal lattice of the 4D cubic lattice, and whose peak intensities are given by the Fourier Transform of the selection window. Due to the close similarity of their selection windows, the structure factor of the OT and that of the ABT are expected to be very similar, with peaks in the same positions but with slight differences in their intensities. Here we will confine our attention to the indexing of Bragg peaks. It is interesting to observe that spherical selection windows were introduced by Grimm and Baake in the context of the diffraction patterns of quasiperiodic tilings [12], as a useful approximation in the calculation of the structure factor. For the optical quasicrystal, in contrast, we see that the spherical window turns out to be the correct, physically imposed choice, provided only that the cut-off I c is large enough.
The 4D basis vectors are K n = (2 √ 2π/λ) ε n . After projection in P, one obtains the four lattice vectors of Eq.2 The peaks of the structure factor S( k) are found at positions k = ∑ h n k n with the condition that equivalent to the requirement that ∑ n h n be even. This condition is equivalent to saying that the positions of its Bragg peaks are those of k √ 2F, the reciprocal lattice of π √ 2 k B, after projection into P. The values of the intensities are determined by the Fourier transform of the selection window, which is of octagonal shape for the ABT, and circular shape for the OT. The resulting differences between the two structure factors would show up only under a detailed comparison of peak intensities. Fig.5 shows the eight vectors ± k n , and the numerically calculated intensity function S( k) defined by where N is the number of sites in the sample. The set of intense peaks nearest the origin correspond to the combinations {±1, ±1, 0, 0} and permutations thereof, in accordance with the condition in Eq.7. The structure factors of successive OT for p ≥ 1 will be the same, only rescaled by powers of α −p , to take into account the inflation in real space.

From the optical tiling to the Ammann Beenker tiling
In the preceding discussion, we assumed that when atoms are loaded into this optical potential, they occupy sites of lowest potential energy i.e., highest intensity, at very low temperature, giving rise to the optical tiling. While overall, the intensities of the occupied peaks vary within the range I max > I > I c , one can define different categories of sites. It is easily seen that the typical value of the potential energy depends on the local environment, as follows. For convenience, we will describe the environment of each site by the letter A,B,... as determined by its position in the octagonal acceptance window W (see Fig.4), and for purposes of illustration, we consider the OT for p = 1. Fig.6 is a contour plot of the intensity I( r) in perpendicular space. The value of the potential on a given site depends on its perpendicular coordinate, decreasing with the distance from the origin. For the values of the potential shown, these contours are close to circular. On the other hand, as we saw (in Fig.4), the distance from the origin in perpendicular space coordinate tends to get larger as z decreases. As shown in Fig.6, the high coordination number sites (8 > z > 5) sites correspond to the region inside  We turn now to the problem of structure differences between the Ammann Beenker tiling and the optical tiling. These arise from the difference in shape of their acceptance windows in perpendicular space. Comparing the two acceptance domains one sees that they overlap over most of the region they occupy. This ensures that a large fraction of the selected points are identical in the two tilings. The differences between the two windows arise in the outlying regions. The results in, firstly, some sites which are present in the AB tiling but are "missing" in the OT. This leads to the empty hexagons and larger n-gons that one sees in Figs.2. Secondly, pairs of "twin-sites" separated by a very short distance δ appear in the OT, whereas in the AB tiling only one of the members of these pairs is present, the other being related by a "phason-flip". These two types of differences are shown in detail Fig.7.
Differences of intensity of the occupied and unoccupied sites are too small to be detected visually. For example, in Fig.7(a) the intensities of the twin sites are 10.8V 0 (left site) and 11.1V 0 (right site). Similarly, in Fig.7(b) In sum, the difference in the shape of the selection windows leads to the appearance of inhomogeneities: as compared to the AB tiling, the OT has larger density fluctuations. The local intensity for the members of a pair of twin sites are slightly different, because their distance from the origin in perpendicular space are slightly different, as shown in Fig.8. Similarly, it is easy to show that the "vacant" site inside, say, an empty hexagon has a local intensity that lies just below the cut-off value. Small local fluctuations could lead to depopulating the twin sites on the one hand, and populating the empty sites, on the other. This suggests that a way to remove the observed defects in the OT would be to turn on small short range repulsive interactions between atoms. One can, indeed, tune the local interaction strengths between atoms in cold gases using Feshbach resonances [5]. In the standard Hubbard model for fermions, the interaction energy is nonzero when the two fermions are on a single site. In the extended Hubbard model, the interaction terms concern particles 1 The two defects shown are situated on the same vertical worm of the ABT. on nearest neighbor sites. Both kinds of terms can be simulated in cold atom gases. By varying the relative strength of the hopping and onsite interaction terms, it is possible, for example, to induce a superfluid to Mott insulator transition in a bosonic system [13]. Extensions of the Hubbard model to further near neighbors are reviewed in [14].
Adding a small short range repulsive interaction energy will make it unfavorable to occupy simultaneously both sites of the twin pairs. One of the atoms corresponds to a larger distance from the origin in perpendicular space, and therefore has a slightly higher onsite energy than its twin. It will therefore move to an unoccupied site with similar onsite energy, corresponding to one of the empty polygons. Let us make the reasonable assumption that all pairwise nearest neighbor interactions are negligible, except for those between the atoms on twin-sites. Let the strength of the repulsive energy for a pair of twin sites be U . Increasing U will make it increasingly unfavorable to occupy simultaneously both sites of a twin pair. The U value needed to shift an atom to an empty site can be estimated by the following argument. Consider a pair of twin sites, whose positions in the perpendicular space are known, for example as in Fig.8. We now consider the change in the energy of the system ∆E if one of the twin sites is vacated by the atom, which migrates to a defect "vacant" site. In this move, an atom is removed from the site whose coordinate r 1 lies within the cutoff radius but outside the octagonal window, and it is replaced elsewhere, on a site whose whose coordinate r 2 lies outside the cutoff radius but within the octagonal window. Such a move costs potential energy, and the change can be calculated by using the projected form of the potential energy in perpendicular space: where ρ 1 and ρ 2 are the radial distances of the initial and final positions of the atom in perpendicular space and where we used the quadratic approximation of the potential. One can calculate the maximum value of this difference, by taking ρ 1 to be on the midpoint of an edge of the octagonal window, and ρ 2 on a vertex. The difference of distances can then be readily calculated, as a function of p. One thus finds that ∆V V max < π 2 16 (2 − √ 2)α −2p . In order for such a hop to occur, the gain in onsite potential energy has to be compensated by the reduction of energy due to the destruction of a twin pair, namely, U . This argument provides an upper bound to the value of the repulsive interaction needed to eliminate twin pairs. For p = 1, for example, one gets U /V max ∼ ∆V/V max ≤ 0.06, and the values for larger p are even smaller. This estimate shows that the repulsive interaction needed would be quite small relative to the onsite potential energy term.
To transform the OT into the AB tiling, in sum, one must introduce a small additional nearest neighbor interactions, over and above the optical potential. As we noted, in cold atom systems, such repulsive interactions can be introduced, and experimentally controlled. The experimental set-up with the laser potential and repulsive interactions appears to be a feasible proposition for obtaining a perfect defect-free Ammann Beenker tiling, at least in principle. In practice, there will be problems associated with slow dynamics, trapping in metastable configurations, and thus possibly a certain amount of disorder at finite temperatures.
The dynamics of this system can be described, to good approximation, by a tight-binding model, using the basis set of the Wannier states defined on the local minima of the optical lattice. As we have already pointed out, the values of the diagonal terms (the onsite energies), depend on the local environment (coordination number). The degree of localization of the Wannier functions depends on the depth of the potential minimum, and also varies to a limited extent. As for the hopping amplitude between two neighboring sites, they depend on several factors: on the distance between the sites, on the profile of the potential, and on the degree of localization of the local Wannier orbitals. These amplitudes will therefore have a variation of values for different pairs of sites. To evaluate the size of the spread of values, a detailed calculation is necessary. However, it can be argued as in [3] that the hopping amplitudes are expected to fall off quickly with the distances and that the two important processes correspond to hopping along edges and across the small diagonal of the rhombus.
With the experimental and theoretical caveats mentioned above, we see that it is theoretically possible to simulate tight-binding Hamiltonians for fermions and bosons on a perfect octagonal tiling. This quasiperiodic tiling has been a subject of theoretical investigation since the mid-eighties. To cite only a few representative works, spectral properties and quantum dynamics in tight binding Hamiltonians have been investigated by many authors using a variety of techniques ( [15][16][17][18][19] . Many body effects in quasicrystals are also a question of current interest. The effect of Hubbard interactions [20] and more recently, the fate of a local magnetic impurity [21] in this tiling were examined, motivated by recent experimental work in heavy fermion quasicrystal compounds [22,23]. It should be possible also to simulate and study the 2D quasiperiodic antiferromagnetic Heisenberg spin model, for which theoretical works have predicted novel ground state properties [24]. Effects of disorder, magnetic fields etc could also be systematically studied by means of a cold atom simulation.

Conclusions
We have outlined a method of obtaining a two dimensional quasicrystal namely, the Ammann Beenker tiling, by trapping cold atoms in a laser generated potential. This would allow, for the first time, experimental studies of important theoretical paradigms for quasicrystals. The study of dynamics of fermions or bosons in a perfect 2D quasicrystal, and described by a Hubbard model, is one example. Another example is the experimental realization of the Heisenberg spin model. Disorder and effects due to magnetic perturbations could be investigated under controlled conditions. Progress in these problems would significantly advance our understanding of the thermodynamic and transport phenomena in real quasicrystals.