{'fraction_non_alphanumeric': 0.033329429675606044, 'fraction_numerical': 0.0115938634500527, 'mean_word_length': 4.502770975641191, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Theory of vortex-lattice melting in a one-dimensional optical lattice

We investigate quantum and temperature fluctuations of a vortex lattice in a one-dimensional optical lattice. We discuss in particular the Bloch bands of the Tkachenko modes and calculate the correlation function of the vortex positions along the direction of the optical lattice. Because of the small number of particles in the pancake Bose-Einstein condensates at every site of the optical lattice, finite-size effects become very important. Moreover, the fluctuations in the vortex positions are inhomogeneous due to the inhomogeneous density. As a result, the melting of the lattice occurs from the outside inwards. However, tunneling between neighboring pancakes substantially reduces the inhomogeneity as well as the size of the fluctuations. On the other hand, nonzero temperatures increase the size of the fluctuations dramatically. We calculate the crossover temperature from quantum melting to classical melting. We also investigate melting in the presence of a quartic radial potential, where a liquid can form in the center instead of at the outer edge of the pancake Bose-Einstein condensates.


I. INTRODUCTION
Since the onset of experiments on Bose-Einstein condensates, vortices have attracted a lot of attention. When the Bose-Einstein condensate is rotated faster than some critical rotation frequency Ω c , a vortex appears in the gas [1,2,3]. Upon increasing the rotation frequency further, the number of vortices increases [4,5,6,7,8], and the vortices order themselves in a hexagonal Abrikosov lattice [9,10,11,12,13,14,15]. If the rotation frequency is increased even further, the very rapidly rotating ultracold bosonic gases have been predicted to form highly-correlated quantum states [16,17,18,19,20]. In these states, the Bose-Einstein condensate has been completely depleted by quantum fluctuations, and quantum liquids appear with excitations that can carry fractional statistics. Some of these states have been identified with (bosonic) fractional quantum Hall states [19,21,22], such as the Laughlin state [23], the Moore-Read state [24] and various Read-Rezayi states [25,26].
In this article we study the quantum fluctuations of vortices in a one-dimensional optical lattice. Optical lattices provide a powerful tool in manipulating a Bose-Einstein condensate. By using a three-dimensional optical lattice, the Bose-Hubbard model [27] was experimentally realized and the superfluid-Mott insulator transition was observed [28]. Moreover, a two-dimensional optical lattice was used to create one-dimensional Bosegases and to study the crossover between a superfluid Bose gas and the "fermionized" Tonks gas [29,30]. Theoretically, the combination of a two-dimensional optical lattice and a single vortex was predicted to give interesting effects around the superfluid-Mott insulator transition [31]. Optical lattices can also be used to provide a pinning potential for vortices [32,33,34,35,36], and experiments on this topic are ongoing [37]. The physics of a single vortex line in a one-dimensional optical lattice is recently extensively studied [38,39,40,41,42,43]. By putting fermions in the vortex core, this system can even be used to create a superstring in the laboratory [44].
A one-dimensional optical lattice divides the Bose-Einstein condensate into a stack of two-dimensional pancake condensates that are weakly coupled by tunneling as schematically shown in Fig. 1. An important consequence of this setup is that the modes of the on-site two-dimensional vortex lattice form Bloch bands as a function of the axial momentum. In this paper we pay special attention to the Tkachenko modes [45], which are the lowest-lying modes of the vortex lattice. Recently, these modes have been investigated for a single twodimensional vortex lattice both experimentally [13,14] and theoretically [46,47,48,49,50,51,52,53].
The number of particles in a single pancake is much smaller than in a Bose-Einstein condensate in a harmonic trap, and hence the fluctuations are much larger. Therefore, this system is a promising candidate to reach the quantum Hall regime. The requirement for that is that the ratio ν = N/N v of the number of atoms N and the number of vortices N v is smaller than a critical value ν c . The ratio ν plays the role of the filling factor and estimates for the critical ν c are typically around 8 [19,54]. However, observed filling factors are up till now always greater than 100, where almost perfect hexagonal lattices form and no sign of melting can be seen [14]. These experiments are carried out with Bose-Einstein condensates consisting of typically 10 5 particles, whereas the maximum number of vortices observed is around 300. Decreasing the number of particles results in loss of signal, whereas the number of vortices is limited by the rotation frequency that has to be smaller than the transverse trapping frequency. Adding a quartic potential, which stabilizes the condensate also for rotation frequencies higher than the transverse trap frequency, has until now not improved this situation [15], although it has opened up the possibility of forming a giant vortex in the center of the cloud [55,56,57]. Applying a one-dimensional optical lattice produces pancake condensates each containing typically 1000 particles, such that quantum fluctuations are strongly enhanced. Experimental signal is still conserved because of the combined signal of all the pancake. Moreover, because of the small number of particles in each pancake shaped Bose-Einstein condensate, finite-size effects become very pronounced in this setup. In particular, the critical filling factor for the melting of the lattice ν c changes compared to the homogeneous situation. As a further consequence, melting is not expected to occur homogeneously but starts at the outside and then gradually moves inwards as the rotation speed increases [58]. Therefore, phase coexistence is expected, where a vortex crystal is surrounded by a vortex liquid.
The optical lattice gives also the exciting possibility to study the dimensional crossover between twodimensional melting and three-dimensional melting, by varying the coupling between the pancakes. This system also exhibits an intruiging similarity to the layered structure of the high-T c superconductors [59,60]. Recently, the density profiles for quantum Hall liquids in this geometry have been calculated [61]. Also the static properties of the lattice in a double and multilayer geometry have been investigated [62,63]. Finally, classical melting between shells of vortices in a single two-dimensional Bose-Einstein condensate has also been studied recently [64].
In this article we study this interesting physics by expanding on our previous work [65]. The paper is organized as follows. In Sec. II we present the theoretical foundations of our work and derive the effective lagrangian for the vortex fluctuations. In Sec. III we discuss the excitation spectrum and we pay in particular attention to the Bloch bands of the Tkachenko modes. In Sec. IV we present result on the quantum melting, and in Sec. V we take into account the effects of a nonzero temperature. We calculate the crossover temperature from quantum melting to classical melting. In Sec. VI we present results on the addition of a quartic potential. We end up with conclusions in Sec. VII.

II. EFFECTIVE LAGRANGIAN
In this section we derive the effective lagrangian for the vortex fluctuations. Using the ansatz that the condensates wavefunctions are part of the lowest Landau level, we find the classical groundstates and expand the action up to second order in the fluctuations.

A. Tight-binding approximation
We closely follow the approach in Ref. [40] for a single vortex, extending it to the case of a vortex lattice. Our starting point is a cigar-shaped Bose-Einstein condensate trapped by the potential where m is the atomic mass and ω r and ω z are the radial and axial trapping frequencies, respectively. Since we assume a cigar-shaped trap, we have that ω z ≪ ω ⊥ and we neglect the axial trapping in the rest of this work. Experimentally this can also be realized by using endcap lasers, which results in a box-like trapping potential in axial direction [66]. In addition, the Bose-Einstein condensate experiences a one-dimensional optical lattice where V z is the lattice depth and λ is the wavelength of the laser light. The lattice potential splits the condensate into N s two-dimensional condensates with a pancake shape. Each of the pancake Bose-Einstein condensates contains on average N atoms, so the total number of atoms is N s N . We consider a deep optical lattice, such that its depth is larger than the chemical potential of the two-dimensional condensate, but still we assume that there is full coherence across the condensate array. This means in particular that the lattice potential should not be so deep as to induce a superfluid-Mott insulator transition. Typically the required lattice depth to reach the superfluid-Mott insulator transition in a threedimensional lattice with one particle per lattice site is of the order of 10E z , where E z is the recoil energy of an atom after absorbing a photon from the laser beam and given by In a one-dimensional lattice the number of atoms in each lattice site is typically much larger than in a threedimensional lattice and the transition into the insulating state requires a much deeper lattice [67]. The superfluid-Mott insulator transition in a one-dimensional optical lattice has recently been observed [68], but also in that case the number of particles per lattice site is of the order of one.
The action describing the system in the rotating frame is given by S = dt d 3 x L(x, t), where the Lagrange density is given by (4) Here Ψ(x, t) is the Bose-Einstein condensate wavefunction, and Ω is the rotation frequency. Moreover, is the angular momentum operator, and g = 4π 2 a/m is the interaction strength, with a the three-dimensional scattering length.
Since we assume a deep optical lattice potential, we can perform a tight-binding approximation and write where i labels the lattice sites and z i = iλ/2 is the position of the ith site. The wavefunction in the z direction Φ(z) is chosen to be the lowest Wannier function of the optical lattice. Because of the deep optical lattice this wavefunction is well approximated by the ground-state wavefunction of the harmonic approximation to the lattice potential near the lattice minimum. The frequency associated with this harmonic trap is and the wavefunction Φ(z) is given by where ℓ L = /mω L . Writing the Lagrange density as L = T − E, the timederivative part of the Langrage density becomes Terms coupling wavefunctions on neighboring sites do not appear in this kinetic part of the lagrangian, because of the orthogonality of the Wannier functions on different sites. The energy part of the Lagrange density can be written as with and ij denotes that the sum is taken over nearestneighboring sites. We have defined the effective twodimensional coupling strength and the hopping amplitude Using the Gaussian wavefunction from Eq. (8) to calculate this hopping amplitude underestimates the hopping amplitude, since this approximation is only good in the vicinity of the center of the wells of the optical lattice, and not in the classically forbidden regions where the overlap between the neighboring wavefunctions is maximal. Therefore, we use the result which is exact for a deep optical lattice [69]. The parameters g ′ and t are both fully determined by the microscopic details of the atoms and the optical lattice.

B. Lowest Landau level approximation
Melting is only expected for a Bose-Einstein condensate that is weakly interacting. The transverse wavefunction can then be taken to be part of the lowest Landau level [70,71]. This implies that we describe the Bose-Eintein condensate as a compressible fluid. Thus we consider wavefunctions of the form where w = (x + iy)/ℓ, w ni (t) = (x ni (t) + iy ni (t))/ℓ and x ni (t) = (x ni (t), y ni (t)) is the position of the n th vortex at site i. The vortex positions are the dynamical variables and are therefore time dependent. Here ℓ is the "magnetic length", which is normally identified with the radial harmonic length ℓ ⊥ = /mω ⊥ . The validity of the lowest Landau level approximation can then in meanfield theory be estimated to be limited by the condition which is recently compared with exact diagonalization results for small numbers of particles [73]. To increase the validity of this study, we use ℓ as a variational parameter instead of fixed it to the radial harmonic length, such that our results are also valid for stronger interactions and slower rotation [74]. The associated frequency is ω = /mℓ 2 . In practice, it turns out that ω is always well approximated by the rotation frequency Ω.
Making use of the fact that the wavefunctions φ n (w) = w n e −|w| 2 /2 / √ πℓ 2 n! form a complete and orthonormal basis for the lowest Landau level wavefunctions, it is easy to derive that within these approximations we have that where the density n i (r) is a function of the vortex positions x ni (t).
From now on distances are rescaled by ℓ, frequencies are scaled by the radial trapping frequency ω ⊥ , and we define a dimensionless interaction strength by means of The on-site part of the energy functional can then be written as: The lowest Landau level wavefunctions and, therefore, also the atomic density, are fully determined by the location of the vortices. To consider the quantum mechanics of the vortex lattice we, therefore, replace the functional integral over the condensate wavefunctions by a path integral over the vortex positions. This involves a non-trivial Jacobian, which does not change the results presented here, because we always consider the case that N ≫ 1. In the calculations we take the scattering length of 87 Rb, λ = 700 nm and V z /E z = 16, which gives U = 25N . The qualitative features, however, do not depend on the value of U .

C. Classical Abrikosov lattice
To determine the quadratic fluctuations around the Abrikosov lattice, we first have to find the classical groundstate. We calculated this groundstate for up to 37 vortices. The number of vortices in the condensate increases with the rotation frequency. For small numbers of vortices, the groundstate is distorted from the hexagonal lattice [75,76]. In general, there are also vortices far outside the condensate. When there are more than 18 vortices in the condensate, there is generally one vortex in the center, while the other vortices order themselves in rings of multiples of six vortices. However, when a new ring of vortices enters the condensate, there is an instability towards an elliptic shape deformation. This is related to the elliptic shape deformation that occurs before a single vortex enters the condensate, and that has been investigated before theoretically [77,78,79,80] and has also been observed [81]. This shape instability plays an important role in the classical melting of the vortex lattice, as we show lateron. Pictures of the classical groundstate are given in Fig. 2 for fixed interaction U and different rotation frequencies. The number of vortices within the condensate as a function of rotation frequency is plotted in Fig. 3, where we limited our study to vortex lattices that exhibit hexagonal symmetry around the origin. In the fluctuation calculation we only consider the regime where configurations consisting of one vortex in the center surrounded by rings of six vortices are stable. In that case, the coarse-grained atomic density is well approximated by a Thomas-Fermi profile [82].

D. Fluctuations
Next we study the quadratic fluctuations by expanding the action up to second order in the fluctuations u ni = x ni − x ni . This yields an action of the form where u i ≡ (. . . , u ni , . . .) is the total displacement vector of all the point vortices on site i. The matrices T, E and J depend on Ω, U , and the classical lattice positions x ni , and are formally given by We also expand the density in the fluctuations by means of n i (r) = n 0 (r) + n u ni · n n (r) + nm u ni · n nm (r) · u mi 1 + n u ni · n ni + nm u ni · n nm · u mi , where n 0 (r) is the equilibrium particle density associated with the classical lattice. We have defined the vectors which form the dipole densities and are associated with linear variations around the classical solution, and the tensors which form the quadrupole densities and are associated with the quadratic variations around the classical solution. Moreover, n n = d 2 r n n (r) (25)  replacements Ω Nv

FIG. 3: The number of vortices
Nv that is within the condensate as a function of the rotation frequency Ω. Depending on the rotation frequency, the groundstate is either a vortex lattice with a vortex in the center surrounded by rings of a multiple of six vortices such that the number of vortices can be written as 1 + 6nr, or a configuration where there is no vortex in the center and the vortices order themselves in rings of multiples of three vortices, such that the total number of vortices in the condensate is a multiple of three. For Ω/ω ⊥ > .9, the first situation is the groundstate (except for a small region, where the Bose-Einstein condensate becomes elliptically deformed) and we only calculate the fluctuations for this case. For this plot the condensate radius was defined as the radius where the angular averaged density drops below 0.003. and n nm = d 2 r n nm (r), such that n i (r) is always normalized to 1. The density n 0 (r), and the tensors n n (r), n nm (r) are independent of time and the layer index i, since only the fluctuations in the vortex positions u i (t) are taken as the dynamic variables. Making use of the fact that the following expressions can be derived for the dipole density For the quadrupole density for a single vortex we get n xnyn (r) = 0, n ynxn (r) = 0, whereas for two different vortices we find All matrices in the action for the fluctuations can be completely expressed in terms of these functions. From Eq. (19) we read off that (n n n m (r) + n m n n (r)) + n nm (r) (31) +2πωU (3n n n m − 2n nm n 2 0 (r) − 2 (n n n m (r)n 0 (r) + n m n n (r)n 0 (r)) + n n (r)n m (r) + 2n nm (r)n 0 (r) .
From a straightforward derivation we obtain also T xnym = n xn n xm + n yn n ym 4 − n xnxm + n ynym 2 , The calculation is simplified considerably by making use of the hexagonal symmetry of the vortex lattice. Because of this symmetry the number of matrix elements to be calculated for each matrix is less then v , when all matrix elements are calculated independently.

E. Diagonalization
To diagonalize this action along the z axis, we perform a Fourier transformation to obtain Finally, we completely diagonalize this action by a transformation that is normalized such that the action becomes where ω α (k) are the mode frequencies of the vortex lattice. This means that the v kα , where k labels the momentum in the z direction and α labels the mode, correspond to bosonic operators with commutation relation This allows us to calculate the expectation value for the fluctuations in the vortex positions, but also for the correlations between the various point vortices. Note that these results were obtained in the lowest Landau level approximation, which corresponds to the compressible limit.

III. TKACHENKO MODES
The Tkachenko modes are almost transverse modes of the vortex lattice. In a harmonic trap with cylindrical symmetry they become modes which are almost angular. In the radial direction their spectrum is discretized, because of the finite lattice size. The number of radial Tkachenko modes equals the number of rings in which the vortices have ordered themselves. For 37 vortices 6 Tkachenko modes can be identified. A close comparison with continuum theory for a finite-size system, where also a discrete spectrum was found [52], is possible but beyond the scope of this work. Moreover, the Tkachenko modes also have a dispersion in the z direction. Without the optical lattice some aspects of these modes were recently investigated [83]. For typical parameters this dispersion is plotted in Fig. 4, while the modes are displayed in Fig. 5. As can be seen in the latter figure, some of the vortices at the edge deviate from having a purely angular motion. As is clearly visible, there is one gapless mode, which is linear at long wavelengths, while the other modes are roughly just tight-binding-like. Moreover, various avoided cross- Mode 1 is the Goldstone mode and corresponds to a pure rotation, while the modes 2-6 are gapped and ordered along increasing gap. They were previously called s-band modes in Ref. [13]. As is visible, some of the vortices at the edge deviate from having a purely angular motion. ings between these modes are clearly visible. The gapless mode is the Goldstone mode associated with the spontaneously broken rotational O(2)-symmetry due to the presence of the vortex lattice. When the tunneling rate is very small, the gapped modes have exactly a tightbinding dispersion and the gapless mode gets a dispersion proportional to sin(kλ/4). This can be understood by observing that in this case the modes are decoupled and the hamiltonian for the gapless mode reduces to the Josephson hamiltonian Writing p φ = −i∂ φ , the action in momentum space thus reads after a quadratic expansion of the Josephson energy from which we deduce the dispersion (39) It is interesting to note that for a small rotation frequency, which implies a small vortex lattice, the Tkachenko modes are not the lowest-lying modes. For U = 10, a Tkachenko mode becomes the lowest-lying gapped mode when Ω > 0.978, but there are many modes in between the second and the third Tkachenko mode. This confirms the expectation that increasing the vortex lattice will bring down the Tkachenko spectrum more and more.

IV. QUANTUM MELTING
In this section we study vortex-lattice melting due to quantum fluctuations. We apply the Lindemann criterion to estimate the position of the melting transition. We study the influence of tunneling between the pancake condensates and compare with a local-density approximation. By looking at correlations between the vortices, we can distinguish between various phases, which are partially melted.

A. Single-layer geometry
Quantum fluctuations of the vortices ultimately result in melting of the vortex lattice. To decide whether or not the lattice is melted, we use the Lindemann criterion, which in this inhomogeneous situation has to be applied locally. The Lindemann criterion means that the lattice is melted, when where ∆ ni is average distance to the neighboring vortices. The critical value c L , is known as the Lindemann parameter. This parameter is a priori unknown, but values found from comparison with Monte Carlo simulations are typically in the range c L = 0.1 − 0.3. The value c L = 0.1 was recently found from elaborate calculations for a triangular vortex lattice in high-T c superconductors, which also compared well to experiments in that case [84]. We, therefore, use this value in our calculations. Note that the results we obtain depend quantitatively on the value of the Lindemann parameter. Changing this value shifts the curves, but the qualitative features remain the same. Due to the inhomogeneity of our system, we have to apply this criterion locally. Because the coarse-grained particle density decreases with the distance to the origin, vortices on the outside are already melted, while the inner part of the crystal remains solid [58]. Therefore, a crystal phase in the inside coexists with a liquid phase on the outside. In Fig. 7 we compute the radius of the crystal phase R cr normalized to the condensate radius R, as a function of the rotation frequency for fixed a number of particles and a fixed interaction strength U , but for various hopping strengths t. Also here we define the condensate radius R as the radius for which the angularly averaged density drops below 0.003. The crystal radius R cr is defined as the radius of the innermost ring of vortices that is melted according to the Lindemann criterion. When according to this definition R cr > R, we set the crystal radius equal to the condensate radius, i.e., R cr = R. The ratio R cr /R shows discrete steps because of the ring-like structure in which the vortices order themselves.

B. Local-density approximation
We compare this with a simple local density calculation, where the criterion N/N v = n(r)/n v (r) = 6 [19] or N/N v = n(r)/n v (r) = 8 [54] is applied locally, by making use of a Thomas-Fermi density profile to describe the coarse-grained atomic density. Substituting the Thomas-Fermi profile where R TF is the Thomas-Fermi radius, in the on-site energy in Eq. (19), we obtain We have introduced here also the Abrikosov parameter β = n 2 (r)/( n(r)) 2 ≃ 1.1596 for the hexagonal lattice.
Minimizing for R TF gives which on substitution gives for the on-site energy This is minimized for which sets the variational parameter ℓ. The Thomas-Fermi radius is then given by In lowest order the vortex density is given by Solving now we get for the crystal radius normalized to the Thomas Fermi radius the result For comparison this line is plotted in Fig. 7. As can be seen there are important finite-size corrections to the local-density approximation. The local-density approximation fails to take into account the discrete nature of the vortex positions. Moreover, it predicts the total melting of the crystal at considerable higher rotation frequencies than the exact answer.
In principle there are corrections to the vortex density due to the Thomas-Fermi profile of the particle density. A better approximation near the center of the trap is [70] n v (r) = 1 However, this vortex density becomes zero well within the Thomas-Fermi radius. As a result the ratio n(r)/n v (r) diverges and is always bigger than 8. Higher-order contributes should be taken into account to solve this problem. Instead we compare with the constant vortex density which corresponds to a Gaussian distribution of the atomic density with radius R TF . Straightforward derivation gives then This shifts the curves for the local-density theory in Fig.  7 a little bit upwards, and makes the comparison even less favourable.

C. Multi-layer geometry
When the tunneling between pancakes is turned on, the fluctuations are also coupled in the axial direction. This decreases the fluctuations in the vortex displacements because the stiffness of the vortices increases. Hence melting occurs for higher rotation frequencies, as is visible in Fig. 7. In Fig. 8 we show the crystal radius for fixed rotation frequency and increasing hopping amplitude. To determine the presence of crystalline order in the axial direction, we calculate the correlation function e iq·(uni−unj) , which is related to the static structure factor. For our gaussian theory, this reduces to For the central vortex the correlation function u 0i u 0j decays exponentially, which signals long-range order. For the other vortices [q·(u ni −u nj )] 2 grows as a logarithm, such that e − [q·(uni−unj)] 2 decays algebraically. This is in agreement with the expectation for a one-dimensional system at zero temperature. In Fig. 9 we show the axial correlation function for one of the vortices outside the center.
Since the lattice spacing λ/2 is the only length scale in the z direction, it also determines the axial correlation length. The power of the algebraic decay scales with 1/ √ t. This can be understood when we assume that the gapless Tkachenko mode gives the dominant contribution to the axial correlation. Using the effective hamiltonian for this mode from Eq. (37), we obtain the following expression: where we rescaled k , in units of the lattice spacing λ/2, r n is the distance of the vortex to the origin and we used that 1 − cos(k(i − j)) = 2 sin 2 (k(i − j)/2). Only the first term in the brackets is divergent for k → 0, so we neglect the second term. We make the approximation sin(k/2) ≃ k/2. When i − j is large we can approximate the rapidly oscillating function sin 2 (k(i − j)/2) by its average 1/2, except for a region near k = 0. We obtain then where the constant C 1 = γ + log(8π) + Ci(8π) ≃ 3.80296.
Using now the fact that E J ∝ t, we conclude that indeed the power of the algebraic decay of e − (q·(uni−unj)) 2 /2 scales with 1/ √ t. For nonzero temperatures we obtain a linear rise of the correlation function (u ni − u nj ) 2 , which can be understood from the same argument, since then we have to add the Bose-Einstein distribution, which for very low temperatures can be approximated by When we again rescale the momentum by the lattice constant λ/2 and approximate sin k/2 ≃ k/2, we observe that the 1/k-factor in the intergrand of Eq. (53) is replaced by 1/k 2 . This gives that up to a constant (u in − u jn ) 2 ∝ i − j and e − (q·(uni−unj)) 2 /2 decays exponentially.

D. Correlations
Melting for small arrays of electrons in quantum dots that are composed of two or three shells [85,86,87,88], and of vortex lattice shells [64,89] is predicted to occur in two stages. First the oriental order between different shells is destroyed, while after that the radial order is destroyed. We also find this behavior for the vortex lattice. For that we compute the correlation function between various vortices and decompose it in radial and angular fluctuations. A shell structure is found, where vortices with almost (but not necessarily exactly) the same distance group together. Between these shells angular fluctuations dominate, while within the shell angular fluctuations are suppressed and radial fluctuations dominate. This can be understood if we assume that the Tkachenko modes dominate the fluctuations, because they leave the rings intact, but change the relative angle between the rings. To quantify this, we introduce the following order parameter to measure fluctuations in the relative angle of neighboring shells for neighboring vortices n and m on different shells, where r n is the distance of vortex n to the origin and only the angular part of the fluctuations in the displacement field is considered. Let N max v be the maximum number of vortices on either of the two shells. The shells are decoupled with respect to each other if with c L = 0.1 the same Lindemann parameter. Furthermore we introduce the correlation function where ∆ mn is the distance between the neighboring vortices m and n and we split this up in a radial and an- gular component. Moreover we distinguish between the case when n and m are on the same shell, which we denote by Γ rad r and Γ ang r , or on different shells, denoted by Γ rad rr ′ and Γ ang rr ′ . We use the definition that order is destroyed if the order parameters exceeds the Lindemann parameter, i.e., Γ > c 2 L = 0.01. When the fluctuations are small, the vortex lattice is in the crystal phase. This phase is defined as Γ nm = Γ rad nm + Γ ang nm < c 2 L , for n and m on the same shell and on different shells. In particular, there is positional order of the radii of the different shells. However, we can distinguish between the case that there exists oriental order between the shells, in which case the relative angle is locked and Γ SS nm < 2π N max v c 2 L . This we call the solid phase (S). When the oriental order between the shells is destroyed we call this phase the solid of solid shells (SS). The melted phase is charactarized by Γ nm = Γ rad nm + Γ ang nm > c 2 L , and in particular the positional order in the radii of the different shells is destroyed. Within this phase it is still possible to have a well-defined angle between vortices on the same shell, i.e., Γ ang r < c 2 L . This we call the hexatic liquid (HL), whereas when this order is destroyed we call it the liquid (L). These phases are schematically indicated in Fig. 10.
These criteria should again be applied locally because of the inhomogeneous density, and we define the phase boundaries in analogy to the previous definition as the radius of the innermost ring that is partly melted. For our parameters it turns out that the angle between neighboring vortices on the same ring is always well-defined, but we can identify the transition between the solid, the solid of solid shells, and the hexatic liquid. The result of this calculation is presented in Fig. 11. In agreement with known theory for vortex lattices, the hexatic symmetry is a very robust phenomenon.

V. CLASSICAL MELTING
The Bose-Einstein temperature for a two-dimensional noninteracting Bose gas in a harmonic trap is given by where ζ(x) is the Riemann zeta function and ζ(2) ≃ 1.28. When the temperatures are much lower than this temperature, there is no thermal cloud present in the gas and we can easily extend our analysis to this regime. This is experimentally relevant, because the zero-temperature limit is difficult to reach [90]. However, according to the Mermin-Wagner theorem, a two-dimensional crystal cannot exist for nonzero temperatures. For an infinite hexagonal vortex lattice this can be seen as follows. Let u(x, t) be the displacement field of the vortex lattice. We define as the longitudinal and transverse fluctuations in the displacement field, respectively. For long wavelengths the action for the vortex lattice is then given by wheren is the average atomic density and c 1 and c 66 are the compressional modulus and the shear modulus of the vortex lattice, respectively. We have used that the vortex-vortex interaction decays faster than 1/r 2 , since the condensate is assumed to be weakly interacting. From this action we read off that the dispersion of the (almost transverse) Tkachenko mode is quadratic, i.e., The fluctuations in the vortex positions can then be calculated as which is a converging integral if we realize that the integration is over the first Brillouin zone. However, for nonzero temperature, we have to add the Bose-Einstein distribution 1/(e βω(k) − 1), which for long wavelengths can be approximated by k B T /ω(k). The temperature fluctuations are therefore given by which diverges logarithmically in the infrared. Therefore, the usual Lindemann criterion always predicts melting for infinite isolated two-dimensional vortex lattices. For a finite system there is a natural infrared cut-off of the divergence, but still there are large contributions from low-lying modes. Any finite amount of tunneling, however, will turn the system into a three-dimensional system, where the Mermin-Wagner theorem does not apply any more. If k denotes the momentum in the z direction, we can then define u L (q, k, t) = (q x u x (q, k, t) + q y u y (q, k, t)/q, (65) u T (q, k, t) = (q x u y (q, k, t) − q y u x (q, k, t))/q. (66) Tunneling results in an additional term in the action that for long wavelengths is given by where J is an effective hopping parameter. This still gives a quadratic dispersion, but with an anisotropic mass, i.e., ω(q, k) = 2 (c 1 q 2 + Jk 2 )(c 66 q 2 + Jk 2 )/n.
The fluctuations can therefore be calculated as which clearly converges in the infrared and the fluctuations remain finite. For a system of coupled pancake Bose-Einstein condensates we can therefore still use the usual Lindemann criterion. In Fig. 12 the vortex-crystal radius is plotted as a function of the rotation frequency for a fixed and strong tunneling and various temperatures. In the regime Ω/ω ⊥ ∈ [.938, .941] there is a dynamical instability towards elliptic shape deformation. This is indicated by a shaded region in the plotted figures. We did not calculate the fluctuations in this regime. This is related to the elliptic shape deformation that occurs before a single vortex enters the condensate, and that has been investigated before theoretically [77,78,79,80] and has also been observed [81]. Because the unstable mode crosses zero, this causes huge fluctuations in the neighborhood of this instability. In Fig. 13 we calculate the temperature for which quantum fluctuations of the vortex crystal are equal to the temperature fluctuations, which defines the crossover temperature between quantum and classical melting. The temperature should be chosen much lower to observe the effects of quantum melting. In Fig. 14 the vortex-crystal radius is plotted as a function of temperature for a fixed rotation frequency. When tunneling between sites is suppressed, the correlation between neighboring vortices where ∆ mn is the distance between the vortex n and m at site i, has been proposed as an appropriate order parameter, with unchanged Lindemann parameter [49,   Fig. 15. Note that a true solid phase of the vortex lattice is not observable in these phase diagrams, not even for the lowest temperature of T = T c /100.

VI. ANHARMONIC RADIAL CONFINEMENT
The inclusion of a quartic potential, in addition to the usual harmonic potential, has attracted interest for several reasons. The quartic potential plays an important role in the understanding of vortex nucleation [94,95]. For our purposes, however, it is more important that in this way rotation frequencies larger than the trap frequency can be applied. This causes the density to be lower in the center of the trap, and eventually gives rise to the formation of giant vortices, i.e., multiple quantized vortices, in the center of the trap [55,56,57]. Experiments in this setup were performed and it was observed that the vortex lattice became disordered, but no giant vortex was observed [15]. Giant vortex formation was observed by artificially removing the inner part of a fast rotating condensate by means of a laser beam [11]. Shape oscillations of a vortex lattice in an anharmonic potential were also studied experimentally [96]. Theoretically, the phase diagram of this system was studied intensively to identify the parameter space where giant vortex formation can take place [97,98,99]. Other aspects that are studied are the dynamics of forming the giant vortex [100,101], oscillations of the vortex lattice [102], aspects on observation [103,104], and stability of quantum fluc- tuations [105].
It is straightforward to extend our analysis to this case. We first of all add the quartic potential  where ℓ ⊥ is the harmonic length associated with the radial trapping. We again use lowest Landau level wavefunctions and make the ansatz that there is one vortex in the center and that the rest of the vortices order themselves in rings of multiples of six. The resulting density profiles for 37 vortices are plotted in Fig. 16 as a function of the rotation frequency. We compare these density profiles with the Thomas-Fermi-like density profiles of the form that come from minimizing the on-site energy functional without paying attention to the lowest Landau level constraint. When the minus sign is chosen there appears a hole in the density, and the condensate shape is annular with inner radius R 1 and outer radius R 2 . We compare the values coming from this ansatz with the inner and outer radius coming from the lowest Landau level densities in Fig. 17. They agree quite well. We also calculate the quantum fluctuations of the vortices. This is only possible in a limited regime, where the ansatz of a central vortex surrounded by rings of six vortices is dynamically stable. In this regime, the inner part of the vortex crystal is melted, although the particle density is nonzero there. We define the liquid radius R L as the radius of the innermost vortex that is part of the vortex crystal. In the region Ω/ω ⊥ ∈ [1.155, 1.25] the ansatz is stable against fluctuations. We find that the liquid radius in this region is almost constant and given by R L /R ≃ .58. There is a second radius which separates the vortex crystal from the vortex liquid that is at the outside of the condensate, which corresponds to the crystal radius R cr defined before. In this region this radius is always equal to the condensate radius or R cr /R = 1.0.

VII. CONCLUSIONS AND OUTLOOK
In this article we derived the theory of vortex fluctuations in a one-dimensional optical lattice. We showed that in this configuration the modes of the vortex lattice get a dispersion in the axial direction. In particular we discussed the Bloch bands of the Tkachenko modes. Based on the Lindemann criterion we studied the melting of the vortex lattice. Because of the inhomogeneous density the melting occurs from the outside inwards. Comparison with a local density-approximation yields important finite-size corrections, because the local-density theory does not take into account the discrete nature of the vortices and overestimates the rotation frequency that is needed for total melting of the vortex lattice. Tunneling between the sites of the optical lattices decreases the fluctuations and causes freezing of the vortex lattice. Temperatures on the other hand increases the fluctuations considerably. The crossover temperature from classical to thermal melting is very low, which makes it an experimentally difficult problem to see the effects of quantum fluctuations. By looking into the correlations between the vortices, several phases can be distinguished, where part of the order is destroyed by quantum fluctuations.
In the one-dimensional optical lattice there is a clear experimental signal for both the fluctuations in the vortex position and the liquid. The fluctuations can be measured in analogy with the situation of a single vortex in an optical lattice. By imaging in axial direction one will see the vortex cores distributed themselves in a gaussian distribution around their equilibrium position, from which the size of the fluctuations can be extracted [40]. In the liquid the vortices are no longer individually visible, which is a clear distinction from the vortex crystal. An interesting question is whether the liquid will completely restore the rotation symmetry that is broken by the presence of the vortex lattice, or the liquid is partly pinned because of the interaction between the vortex liquid and the vortex crystal.
It remains a challenging problem to describe the coexisting crystal-liquid. This will allow to decide on the occurrence of melting based on energy considerations and thus shed more light on the accuracy of the application of the Lindemann criterion in this inhomogeneous situation. This also applies to the phases where partial melting takes place. It is not only important to find a good description of these phases themselves, but also the predicted phase coexistence between them raises interesting new questions.