Liquid crystal phases of two-dimensional dipolar gases and Berezinskii-Kosterlitz-Thouless melting

Liquid crystals are phases of matter intermediate between crystals and liquids. Whereas classical liquid crystals have been known for a long time and are used in electro-optical displays, much less is known about their quantum counterparts. There is growing evidence that quantum liquid crystals play a central role in many electron systems including high temperature superconductors, but a quantitative understanding is lacking due to disorder and other complications. Here, we analyse the quantum phase diagram of a two-dimensional dipolar gas, which exhibits stripe, nematic and supersolid phases. We calculate the stiffness constants determining the stability of the nematic and stripe phases, and the melting of the stripes set by the proliferation of topological defects is analysed microscopically. Our results for the critical temperatures of these phases demonstrate that a controlled study of the interplay between quantum liquid and superfluid phases is within experimental reach for the first time, using dipolar gases.

is the Fermi energy of a 2D non-interacting gas with areal density n 0 , the system is effectively 2D with the dipoles frozen in the harmonic oscillator ground state in the z direction. An external field aligns the dipoles so that their dipole moment d is perpendicular to the y-axis and forms an angle Θ with the z-axis. The dipole-dipole interaction is θ ( , ) = ( − ) /( + ) / V r z D r z 1 3 cos rd for magnetic ones. The strength of the interaction is determined by the dimensionless parameter π = / g mD k 4 3 F 2 , and the degree of anisotropy is controlled by the tilting angle Θ. The system is rotationally symmetric for Θ = 0 and becomes more anisotropic with increasing Θ. Above a critical interaction strength (Θ) g c , it is predicted to form density stripes at = T 0, where the density exhibits periodic modulations of the form Here,  q e c y is the wave vector of the stripes, and / n 2 1 and u are, respectively, the amplitude and the phase of the complex stripe order parameter ψ = ( )/ n iu exp 2 1 1 . The density modulation is formed along the y-direction so as to minimise the interaction energy. The system thus exhibits liquid-like correlations along the x-direction and crystalline correlations along the y-direction. This phase has been predicted by Hartree-Fock theory [14][15][16][17]42 , density-functional theory 19 , and by a variant of the so-called STLS method 18 . Remarkably, Hartree-Fock and density-functional theory predict essentially the same critical coupling strength ( ) .  g 0 06 c for stripe formation at Θ = 0, whereas the STLS method obtains a somewhat higher value. In general, since the formation of stripes occurs at relatively strong coupling it is difficult to make quantitatively accurate predictions, and one will probably have to wait for experiments to sort out the precise physics. For  π Θ . 0 23 , the system is predicted to become a p-wave superfluid 20 , which for strong enough coupling can coexist with the stripe order forming a supersolid 21 . Quantum Monte Carlo simulations also predict the dipoles to form a triangular Wigner crystal at very large values of  g 10 for Θ = 0 24,25 in analogy with the case of bosonic dipoles [43][44][45] . This very strong coupling regime is outside the scope of the present paper.
Stripe phase at finite T and effective XY model. Since the stripe phase breaks translational invariance along the y-direction, it is a quantum analog of a classical smectic liquid crystal 46,47 . Indeed, the system has a manifold of equivalent ground states distinguished only by a constant factor u, which specifies the position of the stripes along the y-direction. Consequently, there are low energy collective excitations associated with a spatially dependent phase ( , ) u x y . Moreover, since a change from u to π + u 2 returns the system to the same ground state, it follows that the low energy degrees of freedom of the stripe phase are described by a 2D anisotropic XY model. Specifically, the simplest form of the elastic free energy congruent with the symmetry of the system is given by Berezinskii-Kosterlitz-Thouless melting. As the stripe phase is described by the XY model, it exhibits algebraic long-range order at sufficiently low temperatures and it melts via the Berezinskii-Kosterlitz-Thouless mechanism due to the proliferation of topological defects [48][49][50][51] . In the case of the stripe phase, the topological defects are dislocations. The phase field for a single dislocation of charge = ± , ± … where the path of the integration encloses the core of the dislocation. The presence of such a dislocation corresponds to inserting Q extra stripes to the left (right) of the dislocation for > Q 0 ( < ) Q 0 . The energy of a single dislocation consists of a core part E c , and a part that scales logarithmically with the size of the system. Pairs of bound dislocations with opposite charges = ± Q 1 ( > Q 1 are energetically suppressed), however, have a finite energy even for an infinite system size and can be thermally excited in the stripe phase. This is due to the fact that the phase fields of the oppositely-charged dislocations cancel at large distances, which results in merely a local disturbance of the density. In Fig. 1, we illustrate dislocation pairs with opposite charges = ± Q 1 centered at ( , ) x y 1 1 and ( , ) x y 2 2 respectively. The stripe amplitude is suppressed in the core regions of the defects due to the large energy cost associated with ∇ ∝ / u r 1 , where r is the distance to the core. From the rescaling → / ⊥ x B B x it follows that the energy of a vertically displaced dislocation pair distance δ apart is the same as that of a pair displaced horizontally by the distance δ / ⊥ B B . Since > ⊥ B B as we will demonstrate below, this shows that the dislocation pairs along the x-direction are more tightly bound than those along the y-direction.
The spontaneous thermal excitation of bound dislocation pairs decreases the elastic coefficients at a macroscopic scale. The softening of the effective stiffness constant B can be calculated from the well-known renormalisation group equations as described in the methods section. At a critical temperature T c st , the renormalised elastic coefficient B R drops to zero by a sudden jump of magnitude π / T 2 c st . This disappearance of elastic rigidity signals the melting of the density stripes.
Calculation of bare stiffness constants. We now turn to a microscopic calculation of the "bare" stiffness constants B and ⊥ B unrenormalised by dislocation pairs. The relevant thermodynamic quantity is the free energy of the system ( ) F q , which depends on the stripe wave vector q. Any non-uniform phase fluctuation increases the free energy by an amount given by (2) for long wave lengths. To extract the elastic coefficients ⊥ B and B , we consider two specific distortions: an infinitesimal rotation and an infinitesimal compression/expansion of the stripes away from the equilibrium configuration, as illustrated in Fig. 2.
These distortions are described by the phase field δ ( ) = − ⋅ u r q r, where δ δ = ⋅ q q e y for the compression and δ δ = ⋅ ⊥ q q e x for the rotation. They are thus equivalent to a variation of the stripe vector δ = + q q q c . Inserting the phase fluctuations into Eq. (2), we obtain the increment of the free energy ( ) where A is the area of the system, and we have used the equilibrium condition ∇ ( ) | = F q 0 q c . We thus find The interaction energy per particle due to stripe formation scales as ( / ) − / D n n n 2 0 3 2 1 0 2 . Assuming that the interaction energy is dominant, we find that the elastic coefficient B scales as ε ( / ) n n g F 1 0 2 for a fixed Θ. The magnitude of B can be further reduced by a geometrical factor depending on Θ, since the system becomes rotationally symmetric for Θ = 0, as we shall discuss below. x y 1 1 and ( , ) The dashed lines indicate the position of the density maxima. In order to microscopically calculate the bare stiffness constants, we employ Hartree-Fock mean-field theory for the free energy, writing µ ( ) where Ω MF is the mean-field thermodynamic potential given by μ is the chemical potential and N is the total number of particles. The quasiparticle energies are ε jk , where = , , j 1 2 is the band index and k is restricted to the first Brillouin zone of the 1D periodic potential set up by the stripes. We subtract the interaction energy E int to avoid double counting. The details of this calculation are given in the methods section.
In Fig. 3, we plot the bare elastic coefficients obtained from this approach as a function of temperature for = g 1 and π Θ = . 0 28 , for which the system has a large stripe amplitude / .  n n 0 6 1 0 at low temperatures. In order to minimize finite size effects, we determine the elastic coefficients by fitting a parabolic curve to the free energy in the vicinity of = q q c in accordance with (3), instead of performing a numerical differentiation following (4). This is illustrated in the insets of Fig. 3. This procedure allows us to obtain numerically accurate values for the elastic coefficients. From Fig. 3, we see that both elastic coefficients decrease with increased temperature. This is expected since thermal excitations of quasi-particles reduce the stripe amplitude and thus their rigidity. We also find that ⊥  B B , which suggests that compressing/expanding the stripes costs more energy than a rotation. This difference in magnitude becomes even more profound for small Θ when ⊥ B is strongly suppressed by the weak anisotropy of the system. Finally we note that for = g 1 and Θ π = .
0 28 , the system is in fact predicted to have additional superfluid pairing at = T 0 21 . However, as demonstrated in ref. 21, the superfluid order has negligible effects on the stripe formation, and it can thus be safely neglected when analysing the elastic properties of the stripes.
In Fig. 4, we plot the bare elastic constants as a function of the tilting angle Θ for = g 1 and ε = .
T 0 01 F . The elastic constant B depends non-monotonically on Θ, first decreasing and then increasing exhibiting a minimum at π Θ .  0 24 . This is consistent with the mean-field phase diagram, which shows that the stripe formation is  somewhat suppressed for intermediate values of Θ 21,42 . To illustrate this, we plot as an inset the stripe amplitude n 1 as a function of Θ; we see that it exhibits the same non-monotonic behaviour as B . It is an intriguing question whether this non-monotonic behaviour is an artifact of the Hartree-Fock approximation or whether it is a real physical effect. In comparison to this behaviour, Fig. 4 shows that ⊥ B increases monotonically in Θ. In particular, we have → ⊥ B 0 for Θ → 0 as shown in detail in the inset. This reflects that the system is rotationally symmetric for Θ = 0 such that a rotation of the stripes costs no energy.
Renormalised stiffness constants and stripe melting. The bare elastic constants obtained from the mean-field theory can now be used as initial values in the RG equations to determined the renormalised elastic constants. We also need the dislocation core energy, which must scale as ~B. Therefore, we write where κ is a constant of order unity. In Fig. 5, we plot the renormalised elastic coefficient B R as a function of temperature, obtained by solving (16) Fig. 6 as a function of Θ for = g 1 and κ = 3. It increases rapidly with Θ, indicating that the degree of anisotropy of the system increases such that the stripes become more rigid. An extrapolation of our calculation for = g 1 and κ = 3 shows that ε .  0 3 . The critical temperature also increases with the coupling strength, scaling as ε ( / )T B n n g c F st 1 0 2 0 . We note that in addition to the explicit linear dependence on g, the T c st can further increase with the coupling strength through the dependence on n 1 . Our results show that in order to observe the stripe phase and the associated BKT physics with dipoles, it is preferable to choose a large tilting angle in addition to having a large dipole moment. However, the tilting angle cannot exceed π Θ .  0 3 above which the system exhibits a density collapse for large coupling strengths 20,21 . We note that the critical temperature increases with increasing core energy κ, since it becomes more energetically costly to create dislocation pairs. A microscopic determination of κ for the stripe phase is unfortunately very complicated and beyond the scope of the present paper. We have therefore adopted a pragmatic approach simply choosing the value κ = 3 in Fig. 6, which is intermediate between the values of two microscopic models: for the XY model on a lattice one has κ π = /2 2 52 , whereas BCS theory yields κ π = / 3 for the 2D superfluid transition, if one equates the core energy with the loss of condensation energy inside a radius given by the BCS coherence length 53,54 . Melting of supersolid phase. The system exhibits p-wave pairing for Θ > ( / ) arcsin 2 3 20 , which can coexist with stripe order for > (Θ) g g c at = T 0 21 . We now determine the critical temperature T c sf for the superfluid transition. The 2D superfluid transition is in principle also determined by the BKT mechanism, where the topological defects are now vortices. For weak pairing, however, the mean-field BCS theory in fact gives a good estimate of the transition temperature. We thus determine the critical temperature by solving the linearised gap equation Here ∆ k is the gap parameter and ( , − ′) V V k k is the effective interaction between the quasiparticles in the stripe phase with energy dispersion ξ k measured from the Fermi surface. The details of this calculation are given in the methods section. The critical temperature obtained from this calculation is shown in Fig. 6 for = g 1 and for several tilting angles. This mean-field result gives an upper bound to the critical temperature, but since ε  T c F sf we expect that a more detailed BKT calculation yield only slightly smaller values. This should be contrasted with the melting of the stripes, where an estimate of the critical temperature from a vanishing stripe order would give a much higher value compared to the BKT calculation. This can be seen from Fig. 3, which shows that the mean-field elastic coefficients remain large up to ε = .
T 0 12 F . Thus, it is crucial to use the BKT theory to analyse the stripe melting.
Using a simple p-wave ansatz for the gap parameter φ ∆ ∆  cos k , where φ is the polar angle of the wave vector k, one can obtain an approximate solution for the critical temperature as ) ( where C is a constant related to an effective momentum cutoff in the integral in (7). We find that the data obtained from solving (7) numerically are in fact very well described by (8) with .  C 0 4. Quantum nematic phase for Θ = 0. Figure 6 shows that the critical temperature for the stripe phase vanishes as Θ → 0. This is a direct consequence of the rotational symmetry rendering = ⊥ B 0 for Θ = 0. In this case, the system is no longer described by the XY model. Instead, an appropriate expression for the elastic energy of stripe fluctuations is 55 where λ is a length comparable to the stripe spacing. Dislocations again play an important role in determining the finite temperature properties of the system described by (9). In contrast to the Θ > 0 case, however, single dislocations now have a finite energy and can be thermally excited. When the presence of the free dislocations is taken into account, a system described by (9) is predicted to be in a nematic phase for < < T T 0 c n , and in an isotropic liquid phase for > T T c n55 . In the nematic phase, the translational order exists only within a length scale ξ D , which is determined by the density of the free dislocations. The stripe orientations, averaged over the length scale ξ D , are however algebraically correlated. As a crude physical picture, one can think of the nematic phase as blobs of stripe order of area ξ D 2 , which are all oriented more or less in the same direction, but which are not positionally correlated with each other. The nematic phase is in this sense analogous to the 2D hexatic phase of a crystal, which exhibits bond orientational order but no long-range translational order 47,56,57 . A quantum hexatic phase was recently predicted to exist in 2D dipolar gases for very strong coupling  g 27 22,23 . The results presented here point out the intriguing possibility to realise a quantum version of the nematic phase with dipoles for smaller coupling strengths. We expect the critical temperature T c n for the melting of the quantum nematic phase to scale as B . However, a quantitative calculation of the critical temperature for the dipolar system requires knowledge of The system is in the stripe phase below the Kosterlitz-Thouless melting temperature which is calculated taking κ = 3 for the core energy. For Θ = 0, the striped phase melts at = + T 0 into a nematic phase with long range orientational order but no translational order. The nematic phase melts into an isotropic liquid at a temperature  T B c n indicated by the blue cross. For  Θ ( / ) arcsin 2 3 the system is in a supersolid phase at = T 0 with both stripe and superfluid order. The transition temperatures of this phase calculated for various tilting angles are indicated by the diamonds. The dashed curve is a fit to the data by (8).
The system collapses for  π Θ . 0 3 . the parameter λ, whose determination is beyond our current theoretical framework. In Fig. 6, we have indicated the critical temperature T c n using a somewhat smaller value than the bare B due to renormalisation effects. We note that Fermi surface deformation 58 leading to a nematic phase 59 has been predicted for fermionic dipoles in 3D.

Discussions
An important question concerns whether the critical temperature for the predicted quantum liquid crystal phases is within experimental reach. As an example, let us consider a recent experiment reporting the trapping of chemically stable 23 Na 40 K molecules in their ground state close to quantum degeneracy. The group obtained an induced dipole moment of = .
d 0 8 Debye and a maximum 3D density of = . × n 2 5 10 3D 11 cm −3 10 . Estimating a corresponding 2D areal density as = / n n 0 3D 2 3 , these values correspond to .  g 0 57. This coupling strength can be increased by reaching a larger fraction of the permanent electric dipole moment of 23 Na 40 K, which is 2.72 Debye 60 , or by increasing the density of the gas. Since the critical temperature for the nematic and the stripe phases both scale as ε ( / ) n n g F 1 0 2 , this indicates that the quantum liquid crystal physics discussed in this paper is within experimental reach, once dipolar gases can be cooled down significantly below their Fermi temperature.
The formation of stripe and superfluid order can be observed as correlation peaks in time-of-flight (TOF) experiments 21,42 . One can also detect the stripes directly as density modulations, either after TOF or in situ, provided that the experimental resolution is sufficiently high. Observing the proliferation of dislocations would directly confirm the microscopic mechanism behind the BKT transition.
An interesting question is how the presence of a harmonic trapping potential in the xy plane will influence the results presented here. In the case of a 2D superfluid Fermi gas, recent experiments combined with Monte-Carlo simulations show that the BKT transition survives the presence of a harmonic potential 41 . The exponent η describing the power law decay of the correlation function ( ) G r was however found to be significantly larger than its value for a homogenous system. We speculate that the same will be the case for the striped system considered here, since the elastic constants scale as εg F and therefore decrease with decreasing density, so that a trap average will lead to a larger η. We note that the present system allows for the measurement of the local stiffness constants near the centre of the trap where the system is nearly homogenous, simply by observing the local stripe fluctuations if the experimental resolution is sufficiently high. Alternatively, one can avoid the complications due to a harmonic trapping potential altogether by implementing the box shaped potentials, which have recently been realised experimentally 61,62 .
Finally, we would like to mention a recent fixed note Monte-Carlo calculation which suggests that the stripe phase is not the ground state for Θ = 0 for any coupling strength 24 . We speculate that this result, which contradicts that of refs [14][15][16][17][18]19,42, is due to the approximate nature of the calculation combined with the fragility of the striped phase, which melts at any non-zero temperature for Θ = 0, as shown by our results.
In summary, we analysed the phase diagram of a 2D dipolar gases, which exhibits stripe, nematic and supersolid phases corresponding to the breaking of translational, rotational and gauge symmetry. For a non-zero tilting angle Θ, the low energy degrees of freedom of the striped phase are described by an anisotropic 2D XY model. We calculated the stiffness constants corresponding to a rotation and a compression/expansion of the stripes microscopically. This should be contrasted with electron systems, where such stiffness constants are often simply unknown parameters of the theory. The stripes were shown to melt via the Berezinskii-Kosterlitz-Thouless mechanism due to the proliferation of dislocations, and we obtained the melting temperature by solving the relevant renormalisation group equations. We also calculated the critical temperature of the supersolid phase. For Θ = 0, the striped phase is stable only at = T 0, which melts into a nematic phase for arbitrarily small temperatures. Our analysis of the melting temperatures demonstrated that these phases should be within experimental reach. An observation of them would constitute a major breakthrough in our understanding of the interplay between liquid crystal and superfluid order in low-dimensional many-body systems.

Methods
Mean-field theory of stripe formation. The mean-field Hamiltonian that takes into account the possibility of stripe formation with a wave vector q is given by 42 where ˆ † c k creates a dipole with momentum k, ε k is the single particle Hartree-Fock energy and h k is a real off-diagonal element defined by The quasi-2D interaction in Fourier space is obtained by averaging the interaction over the harmonic oscillator ground state in the z direction. This gives (up to an irrelevant constant term) 63 Scientific RepoRts | 6:19038 | DOI: 10.1038/srep19038 where Θ is the orientation angle of the dipole (see Fig. 7) and φ is the polar angle of k. We diagonalise the mean-field Hamiltonian by generalising the method described in refs 21,42 to an arbitrary stripe vector q. This yields the Hamiltonian Here γ = ∑ , + +Û c jk j G k G k G annihilates a quasiparticle with energy ε jk , where = , , j 1 2 is the band index, = , = , ± , l l G q 0 1 is the reciprocal lattice vector and k is restricted to the first Brillouin zone of the 1D periodic potential set up by the stripes. We can then calculate the mean-field free energy as µ where Ω MF is the mean-field thermodynamic potential given by (5)  are the scale-dependent stiffness constant and dislocation fugacity respectively. They both decrease with increasing l as the renormalisation due to dislocation pairs at larger length scales are included via the solution of (16). The initial values of ( ) K 0 and ( ) y 0 are the bare (local) values unrenormalised by dislocation pairs, which we calculate microscopically as described in the text. At a critical temperature T c st , the long range renormalised elastic coefficient