# Liquid and crystal phase of dipolar fermions in two dimensions

###### Abstract

The liquid and crystal phase of a single-component Fermi gas with dipolar interactions are investigated using quantum Monte Carlo methods in two spatial dimensions and at zero temperature. The dipoles are oriented by an external field perpendicular to the plane of motion, resulting in a purely repulsive interaction. In the liquid phase we calculate the equation of state as a function of the interaction strength and other relevant properties characterizing the Fermi-liquid behavior: effective mass, discontinuity at the Fermi surface and pair correlation function. In the high density regime we calculate the equation of state of the Wigner crystal phase and the critical density of the liquid to solid first order phase transition. Close to the freezing density we also search for the existence of a stripe phase, but such a phase is never found to be energetically favorable.

###### pacs:

05.30.Fk, 03.75.Hh, 03.75.SsThe recent rapid developments in the field of ultracold dipolar atoms and molecules has opened up new fascinating prospects for investigating many-body effects in quantum degenerate gases with long-range interactions (for a review see, e.g., Ref. Review-Baranov ()). In this respect, single-layer and multi-layer configurations of two-dimensional fermions are particularly intriguing because of the competing interplay, depending on the strength of the dipolar interaction and on the distance between layers, between Fermi liquid behavior, superfluid pairing, crystal order and density-wave instabilities Bruun (); Yamaguchi (); DasSarma (); Santos (); Zinner (); Demler (); Baranov (); Parish (); Shlyapnikov ().

Fermionic molecules of KRb, which can have a strong electric dipole moment, have been created using coherent transfer of Feshbach molecules to their rovibrational ground state Jila1 () and have been brought toward the quantum degenerate regime Jila2 (). Other fermionic molecules are now being actively studied experimentally Zwierlein (); Ketterle (). Atomic species with a large magnetic moment, such as dysprosium, offer a different possibility of realizing degenerate Fermi gases of dipoles that was successfully pursued in the recent experiment reported in Ref. Lev ().

A particularly simple geometrical arrangement of a single-component dipolar Fermi gas in 2D is when the dipoles are oriented perpendicular to the plane of motion by means of a sufficiently strong external field. This configuration has been proven to greatly suppress the chemical reaction rate of molecules, thereby enhancing their lifetime Jila3 (). Here particles at distance interact via a purely repulsive, rotationally symmetric and long range potential. Still the phase diagram at zero temperature is expected to be quite rich: interlayer dimers and a novel BCS-BEC superfluid crossover are predicted in bilayer systems Santos (), while in-plane and out-of-plane density ordered phases are predicted in multilayer systems Zinner (); Demler (). In the case of a single layer, a Fermi liquid with peculiar scattering properties is stable at low density Shlyapnikov () and a Wigner crystal emerges at high density, where the classical potential energy of dipoles largely exceeds their kinetic energy. For intermediate values of the interaction strength an instability at finite wave vector is predicted to set in Yamaguchi (); DasSarma (); Parish (), driving the system to a stripe phase that breaks both rotational and translational symmetry (in the direction perpendicular to the stripes). A similar scenario, involving microemulsion phases (e.g. stripes or bubbles) is expected for the meting of the Wigner crystal at in a 2D Coulomb gasSpivak (). These results are derived within a mean-field approach: an important question concerns the quantitative determination of the phase diagram using more accurate theoretical tools, such as quantum Monte Carlo (QMC) techniques Note1 ().

In this Letter we examine a 2D system of identical fermionic particles of mass that interact with the Hamiltonian

(1) |

where is the distance between particle and and is the intensity of the electric (or magnetic) dipole moment. The strength of the dipolar interaction is conveniently expressed in terms of the dimensionless parameter , where is the characteristic length of the dipole-dipole force and is the Fermi wave vector of the 2D gas determined by the density . The energy scale set by is the Fermi energy . As a function of we investigate the ground-state properties of the Fermi liquid phase including the equation of state, the effective mass, the discontinuity of the momentum distribution at , the pair correlation function and the static structure factor. By comparing the energy of the Fermi liquid (FL) and of the Wigner crystal (WC) phase, we determine the value of the freezing density. Furthermore, in the region of interaction strengths close to freezing, we calculate the energy corresponding to a stripe phase finding that it is never favorable compared to the FL and WC state. The main results on the equation of state are shown in Fig. 1 in units of the Hartree-Fock energy

(2) |

corresponding to the lowest order perturbation expansion of the FL in the interaction parameter Yamaguchi (); Shlyapnikov ().

We use the fixed-node diffusion Monte Carlo method (FN-DMC), a projector technique that, starting from an antisymmetric trial wave function , finds the state having the lowest energy compatible with the many-body nodal surface of that is kept fixed during the calculation. The method provides a rigorous upper bound to the energy of the fermionic ground-state FNDMC ().

Simulations are performed in a box of volume , where we always take . The density is and we use periodic boundary conditions (PBC) in both spatial directions. Since the dipole-dipole interaction is long range, the potential energy contribution to the Hamiltonian, given by the second term in Eq. (1), requires a careful treatment

(3) |

where and label particles in the simulation cell and the vectors correspond to the positions of all images of particle in the array of replicas of the simulation cell. The combination of all images has the same average density of the simulation box and provides a good approximation of the homogeneous medium. In the case of the Coulomb potential the summation in (3) is carried out by means of the Ewald method FNDMC (). For the faster convergent potential the mean interaction energy can be evaluated using the simpler formula

(4) |

where denotes the sum (3) with the constraint and is the contribution from distances larger than assuming a uniform distribution of particles. The cut off length is chosen large enough () to yield an average interaction energy that is independent on , within statistical uncertainty.

Fermi liquid phase: The trial wave function describing the FL phase is assumed of the Jastrow-Slater form

(5) |

Here with are the wave vectors complying with PBC in the square box and is a non-negative function satisfying the boundary condition . The short-range behavior of is of the form , where is the modified Bessel function, and fulfills the cusp condition of the dipole-dipole potential Astra ().

A delicate issue related to QMC calculations of the equation of state is the extrapolation to the thermodynamic limit (TL). In the FL phase apart from the size dependence affecting the potential energy contribution, which we treat using the procedure in Eq. 4, significative shell effects are present in the kinetic energy contribution. We consider closed-shell configurations corresponding to for which the relative error in the energy of the non-interacting gas compared to the TL can be as large as . An extrapolation method based on FL theory is provided by the fitting formula Ceperley1 ()

(6) |

which involves the parameter , determining the inverse effective mass of the particles, and the coefficient of the residual size dependence assumed to be linear in . Here is the QMC output energy of the -particle system with the potential contribution evaluated using Eq. 4. The values of are shown as red symbols in Fig. 2. Their scattered dependence on is considerably suppressed if one accounts for the effective mass, as it is shown by the green symbols corresponding to . A more reliable convergence to the TL is obtained using the method of twist-averaged boundary conditions (TABC) Ceperley2 (). Here the PBC wave vectors of the plane waves in the Slater determinant of Eq. (5) are replaced by , where , are continuous variables in the interval . In the grand canonical implementation of TABC described in Refs. Ceperley2 (); Chiesa () the wave vectors are constrained by and different values of the twist can correspond to different numbers of particles. The number of particles and the energy are obtained from averages over all possible twist angles. With our use of TABC we still find a residual size effect Note2 (). The extrapolation to the TL is performed using Eq. (6) and statistical agreement in between PBC and TABC is obtained for all values of the density (see Fig. 2).

The FN-DMC results of the ground-state energy are reported in Fig. 1. At low density we find good agreement with the result , which was derived in Ref. Shlyapnikov () including corrections to the lowest order expansion (2).

The effective mass obtained from Eq. (6) is shown in Fig. 3 for different values of . At weak coupling our results well reproduce the perturbation expansion reported in Ref. Shlyapnikov (), but for larger couplings the reduction of is less pronounced than the perturbative prediction and approaches the value 0.4 for close to freezing. We also calculate the momentum distribution and using the fit , where is the step function and is a continuous function of , we extract the renormalization factor . Results are reported in Fig. 3. In the inset of Fig. 3 we show in the strongly interacting regime for different numbers of particles. It is worth noticing that finite size effects on this quantity are much less severe here than in the 2D Coulomb gas Holzmann () and are similar to the case of hard-disks investigated in Ref. Drummond ().

We finally calculate the pair correlation function giving the probability of finding two particles at the distance . Results for different values of in the FL phase are shown in Fig. 4. It is interesting to notice that by increasing the interaction strength the short-range repulsion increases and a shell structure starts to appear on approaching the freezing density. The Fourier transform of yields the static structure factor . This quantity can also be calculated directly in the FN-DMC algorithm by evaluating the average of the product of density fluctuation operators . Results are reported in Fig. 5 for both estimators Note3 (). For large values of the direct estimator exhibits a more pronounced peak compared to the smoother Fourier transform at the wave vector corresponding to the lowest non-zero reciprocal lattice vector of the triangular lattice.

Wigner crystal phase: To describe the WC phase we make use of the following trial wave function

(7) |

where the Jastrow correlation term is the same as in the FL phase and the single-particle orbitals in the determinant are constructed with Gaussians, whose width is a variational parameter, centered at the lattice points of the triangular Bravais lattice. In order to enforce PBC, both the number of particles and the box sizes and must be multiples of a primitive cell, which can be chosen as a rectangle of side lengths containing two atoms. Extrapolation to the thermodynamic limit is obtained using a linear fit in over the FN-DMC energies. The results for the WC equation of state are reported in Fig. 1. It is worth noticing that the antisymmetric constraint imposed in the wave function (7) for particle exchange is negligible for the value of the energy. In fact, statistically compatible results are obtained with a node less wave function of distinguishable boltzmanons in agreement with the findings of Ref. Astra (). This behavior is expected at large density, where the energy of the WC phase is given by the result Mora ()

(8) |

obtained by including the contribution from the zero-point motion of phonons to the purely classical interaction energy of a system of dipoles arranged in a triangular Bravais lattice. The above expansion, holding for large , is shown in Fig. 1 and is indeed approached by our QMC results. The difference between the ground state energy of the WC and FL phase is shown in the inset of Fig. 1. From a fit to the equation of state of the two phases we can determine the intersection point at . This value is almost a factor two smaller compared to the critical density Astra (); Pupillo (); Mora () of an equivalent system of bosons having the same mass, density and dipolar strength. This can be understood if one considers that the equation of state of the crystal is practically independent of statistics while the energy of the fermionic fluid is significantly larger than the bosonic one. From the equation of state of the FL and WC phase in the vicinity of the freezing density one can also estimate the width of the region where phase separation occurs driving the first-order liquid to solid transition. By imposing equilibrium of pressure and chemical potential in the two phases, the coexistence region turns out to be , a very small value consistent with a similar finding in the bosonic case Mora (). The pair correlation function and the static structure factor deep in the crystal phase are shown respectively in Fig. 4 and Fig. 5. In particular, exhibits a large peak at corresponding to the lowest non-zero wave vector of the reciprocal lattice.

Stripe phase: A pattern of equally spaced stripes is assumed in the -direction corresponding to the trial wave function

(9) |

where denotes the coordinate of the -th stripe and are the PBC wave vectors in the -direction. The number of fermions is the same in each stripe and once multiplied by the number of stripes determines the overall density in the volume . First, using the variational method, we optimize the width of the stripes and their separation . For the latter quantity we find that is an optimal value corresponding to a square box having . We then perform FN-DMC simulations using PBC with different numbers of particles and we extrapolate to the thermodynamic limit relying on a linear fit. The results are shown in the inset of Fig. 1, where we report the energy difference between the stripe and FL phase. For all values of in the vicinity of the freezing density the stripe phase is never energetically favorable compared to neither the FL nor the WC phase.

In conclusion, we investigated the ground state of a purely repulsive system of dipolar fermions and its liquid to solid transition. Important extensions of this work concern the effect of a tilting angle, making the interaction in the 2D plane anisotropic, and the coupling to a second layer inducing interlayer attraction.

Useful discussions with G. Astrakharchik, M. Holzmann and M. Capone are gratefully aknowledged. This work has been supported by ERC through the QGBE grant. Calculations were performed on the AURORA supercomputer at Fondazione Bruno Kessler.

## References

- (1) M.A. Baranov, Phys. Rep. 464, 71 (2008).
- (2) G.M. Bruun and E. Taylor, Phys. Rev. Lett. 101, 245301 (2008).
- (3) Y. Yamaguchi, T. Sogo, T. Ito and T. Miyakawa, Phys. Rev. A 82, 013643 (2010).
- (4) Kai Sun, Congjun Wu and S. Das Sarma, Phys. Rev. B 82, 075105 (2010).
- (5) A. Pikovski, M. Klawunn, G.V. Shlyapnikov and L. Santos, Phys. Rev. Lett. 105, 215302 (2010).
- (6) N.T. Zinner and G.M. Bruun, Eur. Phys. J. D 65, 133 (2011).
- (7) M. Babadi and E. Demler, Phys. Rev. B 84, 235124 (2011).
- (8) L.M. Sieberer and M.A. Baranov, Phys. Rev. A 84, 063633 (2011).
- (9) M.M. Parish and F.M. Marchetti, Phys. Rev. Lett. 108, 145304 (2012).
- (10) Zhen-Kai Lu and G.V. Shlyapnikov, Phys. Rev. A 85, 023614 (2012).
- (11) K.-K. Ni, S. Ospelkaus, M.H.G. de Miranda, A. Pe’er, B. Neyenhuis, J.J. Zirbel, S. Kotochigova, P.S. Julienne, D.S. Jin and J. Ye, Science 322, 231 (2008).
- (12) S. Ospelkaus, K.-K. Ni, D. Wang, M.H.G. de Miranda, B. Neyenhuis, G. Qumner, P.S. Julienne, J.L. Bohn, D.S. Jin and J. Ye, Science 327, 853 (2010).
- (13) J.W. Park, C.-H. Wu, I. Santiago, T.G. Tiecke, P. Ahmadi and M.W. Zwierlein, Phys. Rev. A 85, 051602(R) (2012).
- (14) M.-S. Heo, T.T. Wang, C.A. Christensen, T.M. Rvachov, D.A. Cotta, J.-H. Choi, Y.-R. Lee and W. Ketterle, preprint, arXiv:1205.5304.
- (15) Mingwu Lu, N.Q. Burdick and B.L. Lev, preprint, arXiv:1202.4444 (2012).
- (16) K.-K. Ni, S. Ospelkaus, D. Wang, G. Quemener, B. Neyenhuis, M.H.G. de Miranda, J.L. Bohn, J. Ye and D.S. Jin, Nature 464, 1324 (2010).
- (17) B. Spivak and S.A. Kivelson, Phys. Rev. B 70, 155114 (2004).
- (18) For the 2D Coulomb gas see B.K. Clark, M. Casula and D.M. Ceperley, Phys. Rev. Lett. 103, 055701 (2009).
- (19) For more details see, e.g., J. Kolorenc̆ and L. Mitas, Rep. Prog. Phys. 74, 026502 (2011).
- (20) G.E. Astrakharchik, J. Boronat, I.L. Kurbakov and Yu. E. Lozovik, Phys. Rev. Lett. 98, 060405 (2007).
- (21) D.M. Ceperley, Phys. Rev. B 18, 3126 (1978).
- (22) C. Lin, F.H. Zong and D.M. Ceperley, Phys. Rev. E 64, 016702 (2001).
- (23) S. Chiesa, D.M. Ceperley, R.M. Martin and M. Holzmann, Phys. Rev. Lett. 97, 076404 (2006).
- (24) The residual kinetic energy size error is too small to allow for a meaningful extraction of the effective mass ratio
- (25) M. Holzmann, B. Bernu, V. Olevano, R.M. Martin and D.M. Ceperley, Phys. Rev. B 79, 041308(R) (2009).
- (26) N.D. Drummond, N.R. Cooper, R.J. Needs and G.V. Shlyapnikov, Phys. Rev. B 83, 195429 (2011).
- (27) The functions , and have been calculated using the extrapolation customary in DMC techniques FNDMC ().
- (28) C. Mora, O. Parcollet and X. Waintal, Phys. Rev. B 76, 064511 (2007).
- (29) H.P. Büchler, E. Demler, M. Lukin, A. Micheli, N. Prokof’ev, G. Pupillo and P. Zoller, Phys. Rev. Lett. 98, 060404 (2007).