Ultracold field-linked tetratomic molecules – Nature

Microwave setup

The microwave setup is described in detail in ref. ¹¹. We use a dual-feed waveguide antenna capable of synthesizing arbitrary polarization using two independent controllable feeds. The finite ellipticity and interference between these feeds result in an observed change in Rabi frequency of approximately ±4% when adjusting the relative phase, thereby contributing to the systematic uncertainty.

In terms of the control electronics, we have upgraded the amplifiers to 100 W (Qualwave QPA-5600-5800-18-47, the gate voltage of which is provided by a custom-made linear power supply) to achieve higher Rabi frequencies. Moreover, we have implemented filter cavities to suppress phase noise. Furthermore, we have incorporated a voltage-controlled phase shifter, enabling dynamic control of the relative phase between the two feeds for fine-tuning the microwave ellipticity. To maintain a constant output power while adjusting the ellipticity, we monitor the power in each feed using a power detector and use a feedback control using a voltage-controlled attenuator.

Dimer loss near the field-linked resonance

We experimentally map out the field-linked resonance by measuring the dimer loss. Extended Data Fig. 1 shows the remained dimer number after a 100-ms hold time at Ω = 2π × 29(1) MHz, Δ = 2π × 9.5 MHz, as a function of ellipticity ξ. The loss dip position matches the theoretical resonance position ξ = 4.8°.

Conditions for efficient electroassociation

We experimentally identify the optimum condition for electroassociation. We obtain the tetramer number from the difference between images with and without the tetramer removal process outlined previously. First, we probe the timescale of the tetramer formation. We ramp the ellipticity from ξ = 0(1)° to 8(1)° and vary the ramp speed. As shown in Extended Data Fig. 2a, we observe the formation of tetramers within 0.3(1) ms and subsequently decay because of the finite lifetime. We estimate that the tetramers scatter on average more than once during the association, bringing them close to thermal equilibrium with the remaining dimers.

Next, we investigate the role of quantum degeneracy in efficient electroassociation. For magnetoassociation of Feshbach molecules, it has been shown that a low entropy sample is crucial to achieve high conversion efficiency, because of the improved phase–space overlap between the atoms⁴⁹. Here we vary the degeneracy of our initial dimer samples by changing the final trap depth of the evaporation⁴². We observe an increase in the conversion efficiency η, that is, the fraction of dimers converted into tetramers, with quantum degeneracy of the dimer gas. We achieve a maximum η = 25(2)% conversion efficiency at T = 0.44(1)T_F. Similar to that for magnetoassociation⁴⁹, a maximum unity conversion efficiency is expected at zero temperature.

Association spectra analysis

We determine the binding energy of the tetramers for different target ellipticities (Fig. 2b) and find excellent agreement between the experimental data and coupled-channel calculations without free parameters.

We assume the dimer loss in the modulation spectra to be proportional to the number of formed tetramers. The line shape can be modelled using Fermi’s golden rule⁵⁰

where ν is the modulation frequency and E_b is the binding energy of the tetramer. The function \(g({{\epsilon }}_{{\rm{r}}})\propto {{\rm{e}}}^{-{{\epsilon }}_{{\rm{r}}}/{k}_{{\rm{B}}}T}\) denotes the number of colliding pairs per relative kinetic energy interval dϵ_r. Here, the temperatures T are obtained from the data located away from the association transitions. The function \(F({{\epsilon }}_{{\rm{r}}})\propto \sqrt{{{\epsilon }}_{{\rm{r}}}}{(1+{{\epsilon }}_{{\rm{r}}}/{E}_{{\rm{b}}})}^{-2}\) denotes the Franck–Condon factor F(ϵ_r) between the unbound dimer state and the tetramer state, which we assume to take the same form as for Feshbach molecules⁵⁰. The product F(ϵ_r)g(ϵ_r) is convoluted with a Gaussian distribution with the width σ to account for the linewidth of the tetramer state and the finite energy resolution. The extracted linewidth shows a similar trend with ellipticity as the theoretical linewidth but slightly larger.

Estimation of the elastic scattering rates

We estimate the elastic dipolar scattering rates of dimer–tetramer and tetramer–tetramer collisions. The scattering rate coefficient is given by β = σv, where \(v=\sqrt{8{k}_{{\rm{B}}}T/{{\uppi }}\mu }\) denotes the average relative velocity and σ denotes the cross-section. In the regime of large dipole moment \(E > {\hbar }^{6}/{\mu }^{3}{d}_{1}^{2}{d}_{2}^{2}\), the cross-section σ can be estimated using the semiclassical formula given by⁵¹

Here d₁ and d₂ are the dipole moments of the two colliding particles, μ is the reduced mass and E is the kinetic energy. We neglect the effect of a small ellipticity ξ and estimate the effective dipole moment of the dimers to be \({d}_{0}/\sqrt{12(1+{(\varDelta /\varOmega )}^{2})}\). The dipole moment of tetramers is roughly twice as large as that of dimers. With that, the above formula provides an estimation for the elastic scattering rates to be 9.7 × 10⁻⁹ cm³ s⁻¹ for dimer–tetramer and 1.9 × 10⁻⁸ cm³ s⁻¹ for tetramer–tetramer. This implies that tens of elastic collisions can occur within the lifetime of tetramers.

Lifetime analysis

For the measurements in time of flight, we verify in the absence of tetramers that the two-body loss between dimers is negligible during the hold time. Thus we fit an exponential decay with a constant offset given by the unpaired dimer number N(t) = 2N_Te^−Γt + N_D. The offset N_D is extracted from the data with ellipticity over 8°, in which the number undergoes a fast initial decay and stays constant afterwards.

To investigate the collisional stability of tetramers, we also assess their lifetimes while the dipole trap remains active. Our observations indicate a combined one-body and two-body loss of the detected dimer number, and we confirm that the two-body loss arises from dimer–dimer collisions. Apart from the data near the collisional threshold ξ = 5(1)°, in which in-trap measurements are influenced by thermal dissociation, we do not detect notable additional loss of tetramers in in-trap measurements compared with those in time-of-flight experiments. The deduced inelastic collision rates are consistent with zero within the error bar. We estimate that more than ten elastic collisions can occur throughout the lifetime of tetramers, which suggests that collisions with tetramers are predominantly elastic.

For measurements in a trap, we ramp up the trap depth by 50% simultaneously with the association, to compensate for the force from the inhomogeneous microwave field. The spatially varying microwave changes the dressed state energy, and thus exerts a force on the molecules that lowers the trap depth and leads to additional loss in the trapped lifetime measurements.

We first measure the total number of tetramers and dimers, and then do a comparison measurement in which we remove the tetramers as described in the main text. As shown in Extended Data Fig. 3a, we observe a two-body decay in the dimer number, in contrast to the time-of-flight measurements. To account for this background loss, we first determine the two-body loss rate Γ₂ and the initial dimer number N_D,0 from the comparison measurement and then perform a fit of one-body plus two-body decay in which we fix Γ₂ and N_D,0. The fit function is given by N_D(t) = 2N_T,0e^−Γt + N_D,0/(1 + Γ₂t). Extended Data Fig. 3b,c shows that the tetramer decay in trap and in free space are similar. The extracted decay rates differ by 9(9) × 10¹ Hz, which we use to obtain an upper bound for the inelastic scattering rate coefficients. By assuming that the additional loss is either purely dimer–tetramer or tetramer–tetramer, we estimate the upper bounds of their inelastic collision rate coefficients to be 2(2) × 10⁻¹⁰ cm³ s⁻¹ and 9(9) × 10⁻¹⁰ cm³ s⁻¹, respectively. Both values are consistent with zero within the error bar. Even for the worst-case estimation, the inelastic collision rate coefficients remain orders of magnitude lower than the estimated elastic dipolar scattering rate coefficients.

The lifetime of the long-range field-linked tetramers is much longer than that observed in polyatomic Feshbach molecules, which are either short lived (<1 μs) (ref. ²²) or unstable in the presence of an optical trap³⁷. These features make them a promising candidate for realizing a BEC of polyatomic molecules. Using the resonance at circular polarization, the improved shielding increases the tetramer lifetime to hundreds of milliseconds. As our experiments suggest that they are stable against dimer–tetramer collisions, it seems promising to evaporatively cool tetramers to lower temperatures⁵².

Association timescale analysis

We apply the following double-exponential fit to the tetramer number as a function of ramp time t in Extended Data Fig. 2a

where τ gives the timescale for association and τ_T gives the timescale for tetramer decay. The time t_T ≈ 0.4(t + t_disso) is the time at which the ramp is above the field-linked resonance, which is about a factor of 0.4 of the association time t and the dissociation time t_disso = 0.5 ms. We extract τ = 0.3(1) ms and τ_T = 2(1) ms.

Hyperfine transitions in the modulation spectra

We observe the effects of the hyperfine structure of NaK molecules in the modulation spectra. When we modulate the ellipticity of the microwave by phase modulation, we generate two sidebands that are offset from the carrier by the modulation frequency ν. When ν matches the ground- or excited-state hyperfine splitting of the dimer, a two-photon hyperfine transition occurs. In Extended Data Fig. 4b, we map out the transition spectrum by Landau–Zener sweeps, in which the modulation frequency is ramped from one data point to the next. If a sweep is performed over a hyperfine transition, molecules are transferred to another hyperfine state causing a depletion of the detected number of dimers. We observe three main hyperfine transitions from 2 kHz to 200 kHz and a few weaker ones. We verify that these transitions are not affected by changes in the ellipticity, which confirms that they are not related to the tetramer states. To obtain a clear spectrum, when measuring the dissociation spectrum, we use a small modulation amplitude to minimize power broadening and ensure that we avoid measuring near these transitions.

Tetramer dissociation spectrum analysis

For modulation dissociation, we first create tetramers at ξ = 8(1)° using electroassociation, then modulate the ellipticity for 2 ms to dissociate them. Meanwhile, we turn off the trap to suppress further association of dimers. Afterwards, we remove the remaining tetramers and let the dissociated dimers expand for another 6 ms before absorption imaging.

In addition to the hyperfine transitions mentioned above, the association of background dimers into tetramers also affects the measurement of the dissociation spectrum. However, it is worth noting that the association spectra are considerably narrower than the dissociation spectrum, and their influence can be mitigated by using a small modulation amplitude. To provide evidence for this, we present a comparative measurement in Extended Data Fig. 4a, conducted under identical experimental conditions, except that the ellipticity ramp is as fast as 0.5 μs so that no tetramers are formed. Note that the modulation time is much shorter than for the association spectra in Fig. 2a. The observed constant background in this measurement demonstrates that the frequencies at which we measure the dissociation spectrum remain unaffected by hyperfine transitions or association.

We fit the dissociation spectrum with a dissociation line shape that is similar to the one used to describe the dissociation of Feshbach molecules³⁹

where Θ(ν − E_b/h) is the step function and γ = 20(7) kHz accounts for the broadening of the signal.

Imaging method for the dissociated tetramers

Here we describe the measurement in Fig. 4b–d. We turn off the trap after the electroassociation and image the cloud after 4.5 ms of expansion time. To image the molecules, we ramp the ellipticity back to circular to rapidly dissociate the tetramers in 0.3 ms, then turn off the microwave and reverse the stimulated Raman adiabatic passage to transfer the dimers to the Feshbach molecule state. Finally, we separate the bound atoms using magnetodissociation, directly followed by absorption imaging of the atoms to minimize additional cloud expansion from residual release energy of the tetramer and Feshbach molecule dissociation.

Angular distribution of the dissociation patterns

We average along the radial direction of the dissociation patterns to obtain their angular distribution, as shown in Extended Data Fig. 5. The distribution of the average optical density shows a sinusoidal oscillation, which matches the p-wave symmetry. We extract the orientation angle ϕ₀ by fitting a function proportional to \(1+c\cos (2(\widetilde{\phi }-{\phi }_{0}))\), where \(\widetilde{\phi }\) is the angle relative to the horizontal axis of the image and c accounts for the finite contrast.

Tetramer lifetime at circular polarization

The lifetime of the tetramers can be improved by shifting the field-linked resonance towards circular polarization, in which the microwave shielding is more efficient. With circular polarization, two nearly degenerate tetramer states emerge above the field-linked resonance at Rabi frequency Ω = 2π × 83 MHz and Ω = 2π × 85 MHz, which corresponds to the two p-wave channels with angular momentum projection m = 1 and m = −1, respectively, as shown in Extended Data Fig. 6. For the m = 1 state, the lifetime at binding energy E_b < h × 4 kHz exceeds 100 ms. In comparison, we show the decay rate for ξ = 5° for which the resonance occurs at Ω = 2π × 28 MHz. For the same binding energy, the lifetime is 10 times shorter than that for the m = 1 state because of the smaller Rabi frequency.

Rovibrational excitations of field-linked tetramers

We investigate only the first field-linked bound state in the current experiment. At higher ellipticities and Rabi frequencies, the potential is deep enough to hold more than one bound state, which corresponds to the rovibrational excitation of the tetramers. For vibrational (rotational) excitations, the radial (axial) wavefunction of the constituent dimers has one or more nodes⁵³. These excited field-linked states have more complex structures, which can be probed similarly with microwave-field modulation.

Field-linked states of polyatomic molecules

Here we discuss the applicability of field-linked resonances to complex polyatomic molecules. For molecules in which the dipole moment is orthogonal to one of the axes of inertia, the same calculation can be performed within the corresponding rotational subspace, as shown in ref. ¹⁰ for CaOH and SrOH. For more complex molecules in which the body-frame dipole moment is not orthogonal to any of the three axes of inertia, the microwave can induce the π transition between the ground state and the m_J = 0 rotational excited state. However, this detrimental π coupling can be suppressed by applying a d.c. electric field to shift the m_J = 0 state away from the m_J = ±1 states, so that the microwave can be off-resonant to the π transition, as shown in ref. ⁵⁴. With that, a similar analysis of field-linked resonances can be applied.

Theory

We apply coupled-channel calculations to study the scattering of molecules governed by the Hamiltonian \(\widehat{H}=-{{\nabla }}^{2}/M+{\sum }_{j=1,2}{\widehat{h}}_{{\rm{in}}}(\,j)+V({\bf{r}})\), where the reduced Planck constant ħ = 1.

The dynamics of a single molecule is described by the Hamiltonian \({\widehat{h}}_{{\rm{in}}}={B}_{{\rm{rot}}}{{\bf{J}}}^{2}+\varOmega {{\rm{e}}}^{-{\rm{i}}{\omega }_{0}t}\left|{\xi }_{+}\right\rangle \left\langle 0,0\right|/2+{\rm{h.c.}}\) with the rotational constant B_rot = 2π × 2.822 GHz. Here, we focus only on the lowest rotational manifolds (J = 0 and 1) with the four states |J, M_J⟩ = |0, 0⟩, |1, 0⟩ and |1, ±1⟩, where M_J denotes the projection of angular momentum with respect to the microwave wavevector. The microwave field of frequency ω₀ and the ellipticity angle ξ couples |0, 0⟩ and \(| {\xi }_{+}\rangle \equiv \cos \xi \,| 1,1\rangle +\sin \xi \,| 1,-1\rangle \) with the Rabi frequency Ω. In the interaction picture, the eigenstates of \({\widehat{h}}_{{\rm{in}}}\) are \(| 0\rangle \equiv | 1,0\rangle ,| {\xi }_{-}\rangle \equiv \cos \xi \,| 1,-1\rangle -\sin \xi \,| 1,1\rangle ,| \,+\,\rangle \equiv u| 0,0\rangle +v| {\xi }_{+}\rangle \) and \(| \,-\,\rangle \equiv u| {\xi }_{+}\rangle -v| 0,0\rangle \), and the corresponding eigenenergies are \({E}_{0}={E}_{{\xi }_{-}}=-\varDelta \) and E_± = (−Δ ± Ω_eff)/2, where \(u=\sqrt{(1+\varDelta /{\varOmega }_{{\rm{eff}}})/2}\) and \(v=\sqrt{(1-\varDelta /{\varOmega }_{{\rm{eff}}})/2}\) with Δ > 0 being the blue detuning and \({\varOmega }_{{\rm{eff}}}=\)\(\sqrt{{\varDelta }^{2}+{\varOmega }^{2}}\) the effective Rabi frequency.

The interaction of two molecules V(r) = V_dd(r) + V_vdW(r) contains the dipolar interaction

and the van der Waals interaction −C_vdW/r⁶ (C_vdW = 5 × 10⁵ arbitrary units; ref. ⁵⁵). We can project the Schrödinger equation in the two-molecule symmetric subspace \({{\mathcal{S}}}_{7}\equiv {\{| \alpha \rangle \}}_{\alpha =1}^{7}=\{| +,+\rangle ,{| +,0\rangle }_{s},{| +,{\xi }_{-}\rangle }_{s},\)\({| +,-\rangle }_{s},{| -,0\rangle }_{s},{| -,{\xi }_{-}\rangle }_{s},| -,-\rangle \}\) as \({\sum }_{{\alpha }^{{\prime} }}{\widehat{H}}_{\alpha {\alpha }^{{\prime} }}{\psi }_{{\alpha }^{{\prime} }}({\bf{r}})=E{\psi }_{\alpha }({\bf{r}})\), where \({| i,j\rangle }_{s}=(| i,j\rangle +\)\(| \,j,i\rangle )/\sqrt{2}\) is the symmetrization of \(\left|i,j\right\rangle \). Under the rotating wave approximation, the Hamiltonian reads

where \({{\mathcal{E}}}_{\alpha }=\{0,-\frac{1}{2}(\varDelta +{\varOmega }_{{\rm{eff}}}),-\frac{1}{2}(\varDelta +{\varOmega }_{{\rm{eff}}}),-{\varOmega }_{{\rm{eff}}},-\frac{1}{2}(\varDelta +3{\varOmega }_{{\rm{eff}}}),-\frac{1}{2}(\varDelta +3{\varOmega }_{{\rm{eff}}}),\)\(-2{\varOmega }_{{\rm{eff}}}\}\) are asymptotic energies of seven dressed states with respect to the highest dressed state channel \(\left|1\right\rangle \) and \({V}_{\alpha {\alpha }^{{\prime} }}({\bf{r}})=\left\langle \alpha \right|V({\bf{r}})\left|{\alpha }^{{\prime} }\right\rangle \).

To obtain the binding energy and the decay rate of the tetramer in the dressed state \(\left|1\right\rangle \), we consider a pair of molecules with incident energy \({{\mathcal{E}}}_{2} < E < {{\mathcal{E}}}_{1}\), the angular momentum l and its projection m along the z-direction. We use the log-derivative method⁵⁶ to numerically solve the Schrödinger equation in the angular momentum basis, that is, \({\psi }_{\alpha }({\bf{r}})={\sum }_{lm}{\psi }_{\alpha lm}(r){Y}_{lm}(\widehat{r})/r,\) where the loss induced by the formation of a four-body complex is characterized using the absorption boundary condition at r_a = 48.5a₀. By matching the numerical solution ψ_αlm(r) with the exact wavefunction in the asymptotic region r > R_c, we obtain the scattering amplitudes \({f}_{\alpha lm}^{{\alpha }^{{\prime} }{l}^{{\prime} }{m}^{{\prime} }}\) and the scattering cross sections \({\sigma }_{\alpha lm}^{{\alpha }^{{\prime} }{l}^{{\prime} }{m}^{{\prime} }}\) from the channel (αlm) to the channel \(({\alpha }^{{\prime} }{l}^{{\prime} }{m}^{{\prime} })\). All results are convergent for (l, ∣m∣) > 7 and R_c > 5 × 10⁴a₀. We note that a different position of the absorption boundary (for example, r_a = 32a₀ and r_a = 64a₀) does not affect the result because the wavefunction has a negligible component inside the shielding core.

Without loss of generality, we concentrate on the cross-section \({\sigma }_{210}^{210}\) of the incident and outgoing molecules in the channel (210). When the incident energy is resonant with the tetramer state, a peak appears in the cross-section \({\sigma }_{210}^{210}\), where the width of the peak is the decay rate of the tetramer. The cross-section \({\sigma }_{210}^{210}\) quantitatively agrees with the lineshape

where G(E) = 1/(E − E_b + iΓ/2) is the tetramer propagator, \({k}_{2}=\sqrt{M(E-{{\mathcal{E}}}_{2})}\) and S_bg are the incident momentum and the background scattering amplitude of molecules in the dressed state channel \(\left|2\right\rangle \), respectively. By fitting \({\sigma }_{210}^{210}\) and σ(E), we obtain the binding energy E_b and the decay rate Γ of the tetramer. We remark that for the incident and outgoing molecules in other channels α ≈ 3–7, the propagator G(E) in equation (7) does not change. Therefore, fitting \({\sigma }_{\alpha lm}^{{\alpha }^{{\prime} }{l}^{{\prime} }{m}^{{\prime} }}\) in a different scattering channel leads to the same binding energy E_b and decay rate Γ.

For a tetramer with a small decay rate, its wavefunction ψ_b(r) can be obtained by solving the Schrödinger equation \({H}_{{\rm{eff}}}{\psi }_{{\rm{b}}}({\bf{r}})={\bar{E}}_{{\rm{b}}}{\psi }_{{\rm{b}}}({\bf{r}})\). The single-channel model H_eff = −Δ²/M + V_eff(r) is determined by the effective potential¹⁵

for two molecules in the dressed state channel \(\left|1\right\rangle \), where \({{\mathcal{F}}}_{\xi }(\phi )\,=\)\(-\sin 2\xi \cos 2\phi ,\theta \) and ϕ are the polar and azimuthal angles of r. The strength \({C}_{3}={d}^{2}/\left[48{{\uppi }}{{\epsilon }}_{0}(1+{\delta }_{r}^{2})\right]\) of the dipole–dipole interaction depends only on the relative detuning δ_r = ∣Δ∣/Ω, whereas the C₆ term describes an anisotropic shielding potential that prevents destructive short-range collisions. Using the B-spline algorithm, we obtain the binding energy \({\bar{E}}_{{\rm{b}}}\) and the wavefunction ψ_b(r) ≈ Y₁₋(r)φ₁(r)/r of the first tetramer bound state, where \({Y}_{1-}(r)=({Y}_{11}(r)-{Y}_{1-1}(r))/\sqrt{2}\). The binding energies \({\bar{E}}_{{\rm{b}}}\) and E_b obtained from the single-channel model and the seven-channel scattering calculation agree with each other quantitatively for small ξ and Ω. For the largest ξ and Ω in Fig. 1, the relative error of \({\bar{E}}_{{\rm{b}}}\) is less than 30%. The tetramer wavefunction in the momentum space is the Fourier transform ψ_b(k) = ∫dre^−ik·rψ_b(r)/(2π)^3/2 of ψ_b(r).

For the modulation dissociation, the transition probability p_k to the momentum state k is determined by the coupling strength \({g}_{{\bf{k}}}=\int \,{\rm{d}}{\bf{r}}{\psi }_{{\bf{k}}}^{* }({\bf{r}}){\partial }_{\xi }{V}_{{\rm{eff}}}({\bf{r}}){\psi }_{{\rm{b}}}({\bf{r}})\). Here, ψ_k(r) represents the wavefunction of the scattering state. The coupling strength g_k is primarily influenced by \({Y}_{1-}(\widehat{k})\), which characterizes the angular distribution of ψ_b(k). This dominance arises because ∂_ξV_eff maintains mirror symmetry with respect to the x–y plane. Therefore, by measuring \({p}_{{\bf{k}}}\approx | {Y}_{1-}(\widehat{k}){| }^{2}\), we can effectively probe the angular dependence of the tetramer state in the momentum space.