A Dielectric Superfluid of Polar Molecules

We show that, under achievable experimental conditions, a Bose-Einstein condensate (BEC) of polar molecules can exhibit dielectric character. In particular, we derive a set of self-consistent mean-field equations that couple the condensate density to its electric dipole field, leading to the emergence of polarization modes that are coupled to the rich quasiparticle spectrum of the condensate. While the usual roton instability is suppressed in this system, the coupling can give rise to a phonon-like instability that is characteristic of a dielectric material with a negative static dielectric function.

The electrical properties of any material are encoded in its dielectric function κ(q, ω), which relates the induced dipole moments within the material to an applied electric field with wave vector q and frequency ω. For its predicted relationship to a variety of novel physical phenomena, the possibility of κ taking on negative values has generated much interest. For the ac case, it is well known that the dielectric function can become negative, for example, at frequencies near an atomic resonance, resulting in the high reflectivity of metals [1]. Similarly, under specially engineered circumstances the magnetic permeability of certain materials can become negative, resulting in a negative index of refraction [2,3].
For a dc field, however, the existence of a negative dielectric function becomes more subtle. Formally, in the infinite-wavelength limit q → 0, a material is stabilized only if its dielectric function, or constant, is non-negative [4]. Such is not necessarily the case for q = 0, however, where the dielectric constant may take on negative values without violating causality [5]. The presence of a negative dielectric constant (κ(q, ω = 0) < 0) is predicted to significantly alter the electrical properties of materials by, for example, giving rise to the attraction of like charges within a material. In this case, the negative dielectric constant is thought to have profound implications for the maximum critical temperature for superconductivity [4]. Experimentally testing such assertions has proven difficult, although the dc limit has been approached in certain materials [6].
In this article, we propose that the realm of negative dielectric constant can be studied in soon-to-be constructed ultracold molecule experiments, in which electrically polar molecules are produced in a Bose-Einstein condensate (BEC). Specifically, we consider molecules whose permanent dipole moments are easily aligned in sufficiently weak applied fields. In this case, the molecules exert fields on one another that are comparable to the applied field, and the BEC acquires dielectric as well as superfluid character. We show that altering either the density of the gas or the strength of the applied field can drive the dielectric constant of the system to be negative, triggering a long-wavelength instability of the condensate. This instability is qualitatively completely different than those that have been studied previously in dipolar BECs with fixed dipole moments, being the phonon [7] and roton instabilities [8]. Additionally, the decay times of this dielectric instability are found to be relatively long, allowing for this state to be observed and studied in detail. Critical in our analysis is the emergence of polarization modes that are coupled to quasiparticle excitations of the BEC. Thus, we show that nontrivial physics emerges in the system of ultracold polar molecules in small fields, and motivate the exploration of this regime, or the regime of the dielectric superfluid, for future experiments.
As a qualitative guide, we first consider a homogeneous dielectric placed in an applied field. This system is characterized by the well-known Clausius-Mossotti relation, which expresses the dielectric constant κ in terms of the polarizability α of the microscopic constituents, where n is the number density of the dielectric [1]. Eq. (1) suggests that large enough values of nα could drive the dielectric constant of the homogeneous material to be negative. As discussed in [4], this is indeed possible for materials with static dielectric functions κ(q, ω = 0) with small but non-zero wave-number (q > 0). Most dielectric materials such as ceramics, glasses and plastics have κ ∼ 1 − 10, mostly due to their large densities [9], while most gases have κ ∼ 1 because their densities are very small. Gases of heteronuclear polar molecules, however, can realize large values of nα due to their possibly large polarizabilities, resulting in very large dielectric constants for these systems even in the dilute gaseous state where n is small. The heteronuclear molecules that we consider here possess permanent dipolar moments d max in their body-fixed frames. In the absence of an applied electric field, the ground state of each molecule is a parity eigenstate and the orientations of the molecular dipole moments average to zero in the laboratory frame. However, the presence of an applied electric field mixes the molecular parity states, and the average dipole orientation becomes nonzero. Because the energy splitting between these states, ∆, can be very small, the resulting molecular polarizabilities can be very large (see table 1 for examples). In particular, this polarizability is distinct from the more familiar polarizability, which is on the order of 1-100 a.u., and arises solely from the electric degrees of freedom. In a weak electric field E, the molecules then acquire a laboratory-frame dipole moment d = αE characterized by the molecular polarizability α ∝ d 2 max /∆, where d max is the limiting value of the molecular dipole moment in large applied fields and ∆ is the splitting of the parity doublet [10]. This polarizability is an intrinsic property of a heteronuclear molecule. The molecular dipole moment in the laboratory frame is therefore uniquely determined by the local electric field via the polarizability α.
While the gas of heteronuclear polar molecules can possess a large polarizability due to the small parity splittings, this feature is quickly lost at finite temperature. Indeed, a thermal gas of polar molecules, at temperatures T such that k B T > ∆, has both parity states equally populated, causing the dipolar character of the system to vanish. However, in an ultracold, Bose-condensed gas the lower molecular state can be dominantly populated, thus realizing high polarizability α. Thus, polar molecules, at sufficiently low temperature, can exhibit polarizability and hence dielectric character, albeit from a different microscopic origin than atoms or non-polar molecules.
To deal with the properties of this inhomogeneous dielectric system, we cannot simply employ the results derived for a homogeneous material. We therefore begin by considering a generalized dielectric formalism. The energy density of a dielectric material with a polarization P(x) induced by an electric field E(x) is given by To account fully for the internal interaction energy of the system, we write the electric field as E(x) = E app +E P (x), where E app is the applied, or external field in the absence of the BEC and E P (x) is the polarization field. The polarization, by definition, is the dipole density P(x) = n(x)d(x), where n(x) is the particle density of the dielectric and d(x) is the local dipole moment of a constituent at point x, or the dipole field. In a dielectric material, the dipoles are induced by the local field and can be expressed, in the linearly polarizable regime, as d(x) = αE(x). This leads to an expression for the polarization of the dielectric, We can now calculate the field produced at a point x due to the presence of this polarization simply by the convolution over the electric field produced by a dipole at a point x ′ , where n is the unit vector in the direction x − x ′ . A self consistent solution to Eqs. (3) and (4) can be used to calculate the energy density of a dielectric material, Eq. (2), given the density n(x) and polarizability α of the constituents. For a Bose-condensed gas of polar molecules in the dilute regime, we interpret the dielectric density as the magnitude squared of the mean-field condensate wave function, n(x) = |Ψ(x)| 2 , where dx|Ψ(x)| 2 = N , and N is the number of molecules in the system. Now, the polarization can be expressed as P(x) = d(x)|Ψ(x)| 2 and the field E P (x) in Eq. (4) can be expressed in terms of the dipole field d(x) and the wave function Ψ(x), The dipole field d(x) gives the dipole moment of a molecule at point x. From Eqs. (3) and (5), it follows that the dipole field obeys the self-consistent expression depending on Ψ(x), where we see that a molecule with polarizability α at point x obtains an induced dipole moment d(x) due to contributions from both the applied field E app and the field produced by the other dipoles in the system. From Eq.
(2), we can now express the energy of this system, including kinetic and trapping energies, as a functional of Ψ(x) and the dipole field d(x). Neglecting any short-range, non-dipolar interactions between the molecules, this energy functional is given by where M is the mass of a polar molecule and U (x) is the trapping potential. The ground state of this system is found by minimizing the total energy with respect to any variation of Ψ ⋆ . Such a task is seemingly straightforward, however, care must be taken when considering the Ψ(x) dependence of the dipole field d(x). Due to the self-consistent nature of Eq. (6), a functional derivative of d(x) with respect to Ψ ⋆ (x) is not straightforward. To handle this relationship, we treat Eq. (6) as a constraint on the system and employ a Lagrange multiplier technique for minimizing (7). We define a modified energy functional Ω[Ψ, d], given by where F[Ψ, d] is defined so that Eq. (6) is satisfied when ∂Ω/∂λ i = 0 and, by construction, λ is a vector Lagrange multiplier with units of electric field. Now, a wave function Ψ(x) that minimizes the total energy of this system is found by setting ∂Ω/∂Ψ ⋆ = 0. Performing this operation, we derive the non-linear Schrödinger Equation, or Gross-Pitaevskii equation (GPE) for Ψ(x), given by where µ = µ E + λ · µ F is the chemical potential of the condensate, whose terms µ E and µ F,i can be calculated by projecting Eq. (9) onto Ψ ⋆ (x) and matching the terms proportional to 1 and λ i , respectively. The components λ i of the vector Lagrange multiplier are found by enforcing ∂Ω/∂d i = 0, and obey the coupled equations where θ i is the angle between n andî. Equations for the components λ j and λ k are found by cyclic permutations of the indices in Eq. (10). Thus, the ground state of the dilute BEC of polar molecules in weak applied fields is described by Eqs. (6), (9) and (10). Interestingly, we see that these equations couple the condensate wave function to the dipole field, so a change in density can result in a change in the dipole moment of the polar molecules, and vice versa. In the strong field limit of fully polarized dipoles, however, d(x) → d max and α → 0. In this case, Eq. (9) becomes identical to the non-local GPE that is typically used to study dipolar BECs [11]. An instructive and relevant case to consider is a gas of polar molecules in a one-dimensional (1D) harmonic trap, where U (x) = 1 2 M ω 2 z z 2 and the applied field is homogeneous in the z-direction, E app =ẑE app . Indeed, such quasi-two dimensional (q2D) geometries are sought after to stabilize dipolar gases against collisional loss [12,13,14] and energetic instabilities due to the attractive part of the dipole-dipole interaction [15]. In this translationally invariant geometry, all transverse (in-plane) components of the electric field due to the dipole polarization will cancel, so we can write E P (x) =ẑE P (x), and subsequently d(x) =ẑd(x). With these simplifications, the GPE (Eq. (9)) can be rewritten as and the equation for the dipole field (Eq. (6)) can be rewritten as where θ z is the angle between n andẑ. In this geometry, the vector Lagrange multiplier has just one non-zero component corresponding to the z-direction, λ =ẑ · λ, where Eqs. (11), (12) and (13) fully describe the BEC of polar molecules in the q2D geometry described above. Thus, we seek solutions for the condensate wave function Ψ(x) and the corresponding dipole field d(x) that self-consistently satisfy these equations. We note that the results presented in the following are nearly quantitatively equivalent to those that follow from neglecting the Ψ(x) dependence of d(x) in the derivation of the GPE, Eq. (9), which amounts to setting λ = 0 in this work. The governing equations for this system are, however, quite general for describing a BEC of polar molecules in weak applied fields, and may be crucial for future studies.
To describe the pure condensate in this q2D geometry, we note that the lowest energy state is that with zero in-plane momentum, Ψ(x) → Ψ(z) = √ n 2D χ(z), where n 2D is an integrated 2D density and χ(z) is the axial wave function of the condensate, normalized to unity. For a condensate wave function with this form, the polarization field, and therefore the dipole field of the condensate, depend only on z. For this case, all of the convolution integrals in Eqs. (11), (12) and (13) can be performed analytically, using the result from Ref. [8] Thus, we see that Eq. (12) reduces to which yields the solution for the condensate dipole field, The form of this result is in stark contrast with the well-known Clausius-Mossotti result (1), where the denominator of Eq. (16) is replaced by 1 − 4π 3 nα, where n is the density of the homogeneous 3D dielectric. Indeed, the Clausius-Mossotti result does not necessarily hold for inhomogeneous systems, like the one we consider here, which has a finite extent in the polarization direction. For this case, the net electric fields of the other dipolar molecules in the system sum to point locally in the direction opposing the applied field. For linearly polarizable molecules, this results in a decrease of the local dipole moment, which is reflected in Eq. (16). The case is the opposite for the homogeneous 3D system, or a quasi-2D system where the applied field lies parallel to the plane of symmetry, or indeed any system that is homogeneous in the direction of the applied or polarization field. In this case, the result [8] can be used to reproduce the Clausius-Mossotti result.
For the case at hand, however, Eq. (16) makes clear the interesting relation between the condensate density and the corresponding dipole field, revealing that the dipole moments of the condensed molecules are smaller where their density is larger. This dipole field is shown in figure 1 for various values of n 2D α|χ(0)| 2 , where we take χ(z) to be a Gaussian of width l z = /M ω z . The widely studied dipolar BECs with fixed dipole moments, trapped in oblate or q2D geometries, are predicted to be unstable to roton-like quasiparticles beyond some critical density or interaction strength set by the dipole moment of the bosons [8,16]. Eq. (16) suggests that this mechanism may be suppressed in the dielectric BEC for larger values of nα as the dipole moments are diminished in the region of higher density.
We explore the stability of this system more quantitatively by investigating the behavior of small fluctuations on top of the ground state, that is, when Ψ(z) → The dipole field corresponding to the condensate mode of the q2D system, given by Eq. (16). For this figure, we take the axial wave function χ(z) to be a Gaussian normalized to unity with width lz. The lines correspond to i) n 2D α|χ(0)| 2 = .046, ii) n 2D α|χ(0)| 2 = 0.70 and iii) n 2D α|χ(0)| 2 = 3.71, and coincide with the points shown in figure 2. For a given polarizability, higher densities reduce the local dipole moment inside the gas. Ψ(z) + δψ(x, t), and δ ≪ 1. Such perturbations to the ground condensed state do not necessarily possess translational symmetry, and a full treatment of the meanfield Eqs. (11) and (12) must be considered. To handle the time-dependence of the fluctuations ψ(x, t), we generalize Eq. (11) to its time-dependent form, where the stationary condensate has the time dependence e −iµt . In this case, µ is replaced by i ∂ t in Eq. (11).
Whereas the usual linearization of the GPE results in a single equation for ψ(x, t), the linearization of Eqs. (11) and (12) results in an equation for ψ(x, t) and an equation for the corresponding linear response of the dipole field. To this end, we take d(z) → d(z) + δσ(x, t), where σ(x, t) describes the linear deviations from the condensate dipole field d(z) due to the fluctuations ψ(x, t). The equation for ψ(x, t) is given by (18) and the equation for σ(x, t), the linear response of the dipole field, is given by Now, a complete description of the small fluctuations of the Bose condensate of polar molecules in the linearly polarizable regime must self-consistently satisfy the coupled Eqs. (18) and (19).
In the q2D geometry, the fluctuations ψ(x, t) take the form of plane-wave Bogoliubov quasiparticles [17], where the u and v amplitudes are normalized to |u| 2 − |v| 2 = 1 and we assume that the lowest-lying quasiparticle excitations occupy the axial condensate wave function χ(z), which amounts to using the single-mode approximation. Typically, inserting this ansatz into Eq. (18) results in a coupled set of Bogoliubov de Gennes (BdG) equations for the frequencies ω [17]. Now, these BdG equations have terms ∝ σ(x, t) in addition to the more familiar terms ∝ ψ(x, t). Using the BdG ansatz (20), the function σ(x, t) spatially decouples into radial plane waves and a function that depends only on z, much like the form of ψ(x, t) in Eq. (20). Thus, a solution to Eq. (19) describes a polarization mode with plane-wave character corresponding to a quasiparticle with in-plane momentum q.
We solve the modified BdG equations by discretizing them on a numeric grid and using the iterative Arnoldi diagonalization method, solving for σ(x, t) in Eq. (19) at each iteration via Gaussian elimination. We take χ(z) to be a Gaussian normalized to unity with width l z , and find no qualitative and little quantitative difference between this Gaussian ansatz and solving for χ(z) exactly in the NLSE (Eq. (9)). The solutions are characterized by the quasiparticle dispersion relation, relating the energy ω of a quasiparticle and its corresponding polarization mode to its momentum q. As is clear from the ansatz (20), any ω with a nonzero imaginary part signifies a dynamic instability of the in-plane homogeneous ground state.
To execute the Gaussian elimination algorithm, we represent the integral in Eq. (19) as a matrix operator in a basis of grid points. This procedure is stable for the momenta that are relevant for characterizing the stability/phase diagram of the system. However, for a given n 2D α|χ(0)| 2 , we find larger wave number(s) for which this matrix is singular and the Gaussian elimination fails. This singularity corresponds to an unphysical divergence in σ(x, t) and is beyond the scope of our model, in which the dipole moments are limited by a physical cutoff d max . The renormalization of this singularity is under active investigation [18]. It does not, however, affect our present conclusions.  (color online) Stability/phase diagram for a) the RbCs system and b) the ThO system in the q2D geometry as a function of applied field Eapp and integrated 2D density n 2D . In the RbCs diagram, the blue region (i) signifies roton instability and the purple and red (ii & iii) shaded regions signify long-wavelength phonon-like instability. In the red region (iii) there is an instability at all small but nonzero momenta, while in the purple region there is an instability at small momenta, but a gap of stability near zero momentum (a long-wavelength instability with roton character). The red region (iii) persists for all log 10 n 2D α|χ(0)| 2 > 0. The black dashed lines are interpolations of the stability/phase boundaries, and the black circles indicate the applied field above which the q2D dipolar BEC has a roton instability when the effect of P(x) is ignored, i.e., when d(x) = αEapp. In the ThO system, the roton instability is completely suppressed and only the phonon-like instabilities are present in the range of applied fields that are permitted in the linearly polarizable regime.
To proceed in characterizing this stability/phase diagram, we consider two bosonic molecular species, RbCs and ThO. The relevant microscopic parameters for these molecules are given in table 1. While RbCs has a relatively small rotational splitting in its ground state (∆ = 7.6 × 10 −8 a.u.), ThO has a very small Λdoublet splitting in its metastable H 3 ∆ 1 state (∆ = 4.6 × 10 −12 a.u.). This results in a very large polarizability for the ThO molecule, and thus enhanced dielectric effects. These molecules remain in the linearly polarizable regime for applied fields E app d max /(2α), corresponding a maximum applied field of about 200 V/cm for RbCs and 12 mV/cm for ThO. We note that both of these species have been produced in experiments [23,24,25], although not yet as quantum degenerate gases.
Figures 2(a) and 2(b) characterize the stability/phase of the RbCs system and the ThO system, respectively, as a function of applied field and integrated density. In 2(a), the blue region (i) shows where the RbCs system has a roton instability due to excitations of finite wavelength ∼ 2πl z [8]. For dipoles of fixed dipole moment, the instability would occur at smaller dipole moment for increasing density, i.e., the boundary of the blue region would be a strictly decreasing function of n 2D α (shown by the black circles). However, in the dielectric BEC this trend reverses beyond a critical point and the rotons become more stable, owing to the diminishing dipole moment in the center of the gas (as anticipated in Eq. (16)). Whereas the roton instability is suppressed in the stability/phase diagram for the RbCs system, it is completely absent from the linearly polarizable ThO system because the large polarizabilities of this molecule serve to suppress the dipole moment in the center of the condensate more rapidly as a function of applied field.
By contrast, the purple and red (ii & iii) regions, present in both diagrams 2(a) and 2(b), signify long-wavelength instabilities. In the purple region (ii), the long- wavelength instability occurs at small but finite momenta, leaving an interval of dynamically stable modes near q = 0. As n 2D α|χ(0)| 2 becomes greater than unity (red region), the instability occurs at all small, non-zero momenta, and looks like a phonon instability. The imaginary parts of these dispersions are given by purple dot-dashed (ii) and the red dashed (iii) lines in figure 3, respectively. This behavior, being independent of E app and occurring at a critical value of nα, is reminiscent of the long-wavelength instability of a homogeneous dielectric with negative static dielectric constant. Indeed, this instability signifies the transition to a negative static dielectric function in the q2D geometry, wherein novel physical behavior, such as the attraction of like electric charges, is expected to occur [4].
We find that the emergence of the dielectric instability is universal, and occurs at n 2D α|χ(0)| 2 = 1 for any linearly polarizable molecule at zero temperature. For RbCs in a trap with axial frequency ω z = 2π × 20kHz, this corresponds to a large critical integrated density of n 2D = 4.02 × 10 12 cm −2 , while for ThO, the critical density is just 7.85 × 10 6 cm −2 because these molecules have much larger polarizabilities. Indeed, the stability/phase diagram of the ThO system in figure 2(b) is identical to that of RbCs in figure 2(a), but with a much smaller cutoff in E app (necessary to remain in the linearly polarizable regime). Thus, the onset of the dielectric instability occurs at n 2D α|χ(0)| 2 = 1 for both species, but the region of roton instability is far suppressed from the ThO diagram due to this molecule's large polarizability.
For the mean-field theory presented in this work to hold, we require that the gas be sufficiently dilute. In the process of mapping these stability/phase diagrams shown in figure 2, we check that the diluteness criteria is satisfied by calculating the ratio of the characteristic interaction length to the average interparticle spacing across the density profile of the condensate, which is characterized by the 2D gas parameter n 2D a 2 dd (z), where a dd (z) = M d 2 (z)/3 2 is the z-dependent characteristic dipole length of a molecule in the system [8], varying with the condensate density profile. We find that this gas parameter is sufficiently small ( 10 −2 ) near the stability/phase boundaries in figure 2 to justify our use of a condensate mean-field, while it can become larger (> 1) for larger densities and applied fields. In this case, two-body correlations may become important in characterizing the ground state of the gas, however, the dielectric properties of the system should share the qualitative features of the mean field ground state. Similarly, we expect the qualitative dielectric properties to persist if the molecules possess large Van der Waals lengths [26], which can approach the mean intermolecular spacing at sufficiently high density.
In a realistic experimental setup, a radial harmonic trap is present in addition to the tight axial trap. Here, long-wavelength quasiparticles manifest as breathing or quadrupole modes. For the ThO condensate, we consider occupations of N = 500, 1000, 2000 molecules in a trap with radial frequency ω ρ = 2π × 200 Hz and the same axial frequency as the q2D system. We employ the local-density approximation (LDA) to solve the NLSE and calculate the breathing and quadrupole mode frequencies. This approximation works well to describe these long-wavelength modes in the dipolar system. We find that, at critical applied fields similar to those in the q2D regime, the breathing and quadrupole modes develop imaginary frequencies, showing that the phonon-like dielectric instability persists in the trapped system, that is, the decay proceeds via the longest wavelength modes available, with a decay time on the order of milliseconds.
To summarize, we have considered a BEC of polar molecules in an applied electric field that is sufficiently weak to keep the dipoles in the linearly polarizable regime. In this regime, the system develops dielectric character, resulting in a suppressed roton instability and the emergence of a dielectric instability due to the development of a negative static dielectric function. While this physics is essential to consider for any weak-field studies of molecular BECs, it also introduces new physics that is accessible due to the possibly very large polarizabilities of heteronuclear molecules.