Prospects for thermalization of microwave-shielded ultracold molecules

We study anisotropic thermalization in dilute gases of microwave shielded polar molecular fermions. For collision energies above the threshold regime, we find that thermalization is suppressed due to a strong preference for forward scattering and a reduction in total cross section with energy, significantly reducing the efficiency of evaporative cooling. We perform close-coupling calculations on the effective potential energy surface derived by Deng et al. [Phys. Rev. Lett. 130, 183001 (2023)], to obtain accurate 2-body elastic differential cross sections across a range of collision energies. We use Gaussian process regression to obtain a global representation of the differential cross section, over a wide range of collision angles and energies. The route to equilibrium is then analyzed with cross-dimensional rethermalization experiments, quantified by a measure of collisional efficiency toward achieving thermalization.

We study anisotropic thermalization in dilute gases of microwave shielded polar molecular fermions.For collision energies above the threshold regime, we find that thermalization is suppressed due to a strong preference for forward scattering and a reduction in total cross section with energy, significantly reducing the efficiency of evaporative cooling.We perform close-coupling calculations on the effective potential energy surface derived by Deng et al. [Phys. Rev. Lett. 130, 183001 (2023)], to obtain accurate 2-body elastic differential cross sections across a range of collision energies.We use Gaussian process regression to obtain a global representation of the differential cross section, over a wide range of collision angles and energies.The route to equilibrium is then analyzed with cross-dimensional rethermalization experiments, quantified by a measure of collisional efficiency toward achieving thermalization.
The ever growing interest in quantum control of polar molecules motivates the cooling of molecular gases to unprecedented cold temperatures [1][2][3][4][5].In bulk gases, reaching such temperatures can be accomplished through evaporative cooling [6], a process which throws away energetic molecules and leverages collisions to rethermalize the remaining, less energetic, distribution.Understanding and controlling 2-body scattering for thermalization is, therefore, of great importance for ultracold experiments.To this end, the exciting advent of collisional shielding with external fields has permitted a large suppression of 2-body losses between molecules [7][8][9][10][11].Thermalization relies instead on the elastic cross section, which is generally dependent on the field-induced dipoledipole interaction and their energy of approach.Of particular interest to this Letter is collisional shielding with microwave fields [12][13][14][15], recently achieved at several labs around the world [16][17][18][19].
In analogous gases of magnetic atoms with comparatively small dipole moments, dipolar scattering remains close-to-threshold [20] at the ultracold but nodegenerate temperatures of T ∼ 100 nK [21][22][23][24].For dipoles, threshold scattering occurs when the collision energy is much lower than the dipole energy E dd , in which case the scattering cross section becomes energy independent [25] with a universal analytic form [26]. Numerical studies of thermalization are made much simpler at universality, since collisions can be sampled regardless of collision energy [27,28].However, this convenience is lost with the polar molecular gases of interest here.Take for instance a gas of fermionic 23 Na 40 K, as we will concern ourselves with in this study.This species has a large intrinsic dipole moment of d = 2.72 D, so that even ultracold temperatures have majority of collisions occurring away from threshold with an energy dependent cross section.
FIG. 1.A log-log plot of T vs N during a forced evaporation protocol.The plot compares the evaporation trajectory for microwave shielded 23 Na 40 K when scattering is realistic and non-threshold (solid black curve), to the artificial case of threshold scattering (dashed red curve).Both 1 and 2-body losses are assumed negligible and ignored here.
By continuously lowering the energetic depth of the confining potential U (t) = U 0 exp(−t/τ ) over a time interval τ , highly energetic molecules are forced to evaporate away, lowering the number of molecules along with the gas temperature as shown in Fig. 1.For the plot, Eq. ( 1) is solved by taking evaporation to occur with an initial trap depth U 0 /k B = 4 µK over τ = 0.5 s, in a harmonic trap with mean frequency ω = 2π × 100 Hz, starting at temperature T 0 = 400 nK and molecule number N 0 = 20, 000.The evaporation efficiency, defined as the slope of T vs N on a log-log scale, is governed by the thermalization rate γ th .The figure shows efficient cooling for the low-energy threshold cross sections (dashed red curve), and significantly less efficient cooling for the realistic cross sections (solid black curve).The remainder of this Letter provides the microscopic mechanisms that lead to this dramatic difference, and efficient theoretical tools we employ to obtain these conclusions.
Shielded collisions-Central to this study, are collisions that occur between molecules shielded by circularly polarized microwaves [15].The resulting potential energy surface between two such molecules is conveniently described by a single effective potential [31]: where r = (r, θ, ϕ) is the relative position between the two colliding molecules, Ê is the axis along which the dipoles are effectively aligned, is the effective molecular dipole moment and Here ∆ and Ω are the detuning and Rabi frequency respectively, of the microwaves.A y = 0 slice of the effective microwave shielding interaction potential is plotted in the inset of Fig. 2. Notably, the long-range 1/r 3 tail of V eff (r) is almost identical to that of point dipole particles, modified only by an overall minus sign.As a result, the close-to-threshold elastic cross sections for microwave shielded molecules are identical to those for point dipoles.
It is natural to introduce units based on the reduced mass µ, dipole length and dipole energy: respectively.Threshold scattering is then expected to occur for collision energies E ≪ E dd .With the microwave parameters ∆ = 2π × 15 MHz and ∆ = 2π × 9.5 MHz, which will be assumed in what follows, the molecules see a dipole length of a d ≈ 3900a 0 , corresponding to a dipole energy of E dd /k B ≈ 360 nK.Therefore, temperatures comparable to E dd /k B are insufficient to keep molecular scattering in the threshold regime [25].Moreover, since the dipole energy scales as E dd ∼ d −4 , larger dipoles require much lower temperatures to achieve universal dipolar threshold scattering as alluded to earlier.
Away from threshold, the integral cross section σ in the presence of microwave shielding (dashed black curve), develops a nontrivial energy dependence that clearly differs from that of plain point dipoles (dotted blue curve) as illustrated in Fig. 2. The plotted cross sections were obtained from close-coupling calculations logarithmically spaced in energy, with a universal loss short-range boundary condition [32] (see Supplementary Material for further details).Away from threshold at E ≈ E dd , the microwave shielded integral cross section does not deviate much from its value at threshold (solid red line in Fig. 2).But the differential cross section could still have its anisotropy changed substantially, which is what ultimately affects thermalization [26].For a study of both non-threshold differential scattering and its implications to thermalization in nondegenerate Fermi gases, we take its nonequilibrium evolution as governed the Boltzmann transport equation [33].Formulated in this way, numerical solutions treat the molecular positions and momenta as classical variables, while collisions can be efficiently computed by means of Monte Carlo sampling [27,34].But on the fly close-coupling calculations would be too expensive for such sampling over a broad range of collision energies and angles.Instead, we propose the following.
Gaussian process fitting-At a given collision energy, the elastic differential cross section D el , is a function of the dipole alignment axis Ê, and the relative ingoing and outgoing momentum vectors ℏk and ℏk ′ , respectively.Collectively, we refer to this set of parameters as β.By first performing close-coupling calculations at several well chosen collision energies E = ℏ 2 k 2 /(2µ) [35], we can use the resultant scattering data to infer an M -dimensional continuous hypersurface that approximates D el , with a Gaussian process (GP) model [36][37][38].
GP regression is a machine learning technique used to interpolate discrete data points, stitching them together to form a continuous global surface.To do so, a GP assumes that D el (β) evaluated any 2 nearby points in its coordinate space, β i and β j , are Gaussian distributed with a covariance given in terms of a function K(β i , β j ), called the kernel.A parameterized functional form for the kernel is chosen prior to the surface fitting process, reducing the task of combing through an infinite space of possible functions that best match the data, to a minimization over the kernel parameters.This minimization step is referred to as training the GP model.Several symmetries in the differential cross section help to reduce the computational load of training slightly.Rotated into the frame where Ê points along the z axis, which we refer to as the dipole-frame, the unique hypersurface regions effectively live in an M = 4 dimensional space, with coordinates β = (E, η, θ s , ϕ s ).As defined, η = cos −1 k • Ê is the angle between the dipole and incident relative momentum directions, where it is convenient to select k to lie in its x, z plane.The angles θ s and ϕ s , denote the inclination and azimuthal scattering angles respectively, in this frame.Doing so, the differential cross section possesses the symmetry Consequently, we only need to specify the differential cross section for angles within the domain η, θ s , ϕ s ∈ [0, π], to fully describe its global structure.More details of the appropriate frame transformations are provided in Supplementary Material.
To perform the interpolation with GP regression, we utilize the Matérn- 5  2 kernel [39], which is better able to capture the sharp jumps in a non-smooth function, over higher-order differentiable kernels such as the radial basis function.This kernel contains a parameter w that sets a length scale over which features of the data vary in coordinate space, that is optimized during the model training process.This kernel is typically not ideal for periodic input data, so we make the periodicity of the angles (η, θ s , ϕ s ) explicitly known to the GP model by training it with the cosine of these angles, instead of the angles themselves.Furthermore, log 10 (E/E dd ) is fed into the GP model in place of E, to reduce the disparity in fitting domains between each coordinate of β.The GP model is trained over the range log 10 (E/E dd ) = −6 to 2, corresponding to collision energies of E/k B ≈ 0.36 pK to 36 µK.After training on ∼ 10, 000 samples of D el (E, η, θ s , ϕ s ), the resulting GP fit obtains a meansquared error of ≈ 0.5% against the close-coupling calculations [40], which we take as an accurate representation of the actual cross section.
In Fig. 3, we plot the total cross section σ(E, η) = D el (E, η, Ω s )dΩ s , at various collision energies.There is a marked variation in the η dependence, indicating a higher tendency for side-to-side collisions (η = 90 • ) over head-to-tail ones (η = 0 • ) at higher energies.To highlight the dominant anisotropic scattering process, Fig. 3 also provides plots of the differential cross section at η = 45 • , the approximate angle at which σ is maximal.As energy increases from subplots (a) to (d), the scattered angle dependence of D el becomes biased toward forward scattering, reducing the effectiveness of collisions for thermalization as discussed below.Alphabetic labels in Fig. 3 consistently correspond to the collision energies: (b) E = 0.2E dd , (c) E = 2E dd and (d) E = 20E dd .The Born approximated cross sections at threshold [26] are labeled with (a).
Collisional thermalization-Fast and easy access to the accurate differential cross section via its GP model now permits accurate theoretical investigations of nondegenerate gas dynamics.More specifically, we are concerned here with a gas' route to thermal equilibrium.A common experiment for such analysis is cross-dimensional rethermalization [41], in which a harmonically trapped gas is excited along one axis, then left alone to re-equilibrate from collisions.
We present results in terms of the temperatures along each axis i, defined in the presence of a harmonic trap as T i = (⟨p 2 i ⟩/m + mω 2 i ⟨q 2 i ⟩)/2, where ⟨. ..⟩ = d 3 qd 3 pf (q, p)(. ..) denotes a phase space average over the phase space distribution f in molecular positions q and momenta p, while ω i are the harmonic trapping frequencies.As is usual in cross-dimensional rethermalization, we consider an excitation of axis i then proceed to measure the thermalization rate along axis j.This is modeled by taking axis i to have an initial out-ofequilibrium temperature T i = T 0 + δ i /k B , with a perturbance in energy δ i , while the the other 2 axes are simply at initial temperature T 0 .
FIG. 4. Plot of εij as a function of the dipole tilt angle Θ, for all 6 unique configurations (subplots a to f) of the excitation axis i, and measured thermalization axis j.The solid red curves are the analytic εij results derived with the Born approximated cross section at threshold, whereas the dashed-dotted curves are those from Monte Carlo integration using the GP interpolated cross sections, at temperatures T = 10 nK (black), T = 100 nK (dark gray), T = 400 nK (gray), and T = 1 µK (light gray).The dashed blue lines are the efficiency for purely p-wave collisions, εp = 1/4.1.
In the case of a dilute gas, the relaxation of T j follows an exponential decay in time, whose rate γ ij is related to the standard collision rate γ coll , by a proportionality factor ε ij = γ ij /γ coll .As defined, the quantity ε ij is the inverse of the so-called number of collisions per rethermalization [41,42], a measure of thermalization common to the literature [10,17,18].We opt to utilize its inverse instead as it is the more natural definition to discuss efficiency of evaporative cooling.Usually defined as γ coll = ⟨n⟩⟨σv r ⟩ with phase space averaged number density ⟨n⟩ and 2-body elastic rate ⟨σv r ⟩, ε ij represents the efficiency of each non-threshold collision toward thermalization of the gas.This collisional efficiency is formally cast in terms of the integral where ∆κ 2 i = κ ′2 i − κ 2 i is the collisional change in adimensional relative momenta κ = p r (mk B T 0 ) −1/2 , α ij = 3/2 if i = j, and α ij = −3 otherwise (see Supplementary Materials).The integral above has been evaluated analytically in the threshold scattering regime [28], both for identical dipolar fermions and bosons.
Evidently from Eq. ( 5), ε ij is symmetric in its indices which leaves only 6 unique configurations of i and j.As-serting the dipoles lie in the x, z-plane and tilted with angle Θ = cos −1 Ê • ẑ, we compute Eq. ( 5) with Monte Carlo integration [43] and plot the results in Fig. 4.Each subplot (a to f) shows a different (i, j) configuration, within which, ε ij is plotted against the dipole tilt angle Θ as dashed curves, for the temperatures T = 10 nK (black), T = 100 nK (dark gray), T = 400 nK (gray) and T = 1 µK (light gray).Interestingly, the ε ij terms involving excitation or rethermalization along y essentially lose their dependence on Θ around 400 nK, beyond which collisions are less efficient than even nondipolar pwave scattering (dashed blue line in Fig. 4) [44] for all Θ.This decrease can be intuited by looking at the differential cross section around η = 45 • , around which the total cross section is maximal.As evidenced from the subplots of D el in Fig. 3, forward scattering is favored at higher collision energies, limiting momentum transfer between axes and therefore, also the efficiency of collisions toward rethermalization.Preferential forward scattering is what ultimately leads to the reduction in evaporation efficiency, earlier described and seen in Fig. 1.There, the rate of thermalization was approximated by the average γ th = γ coll i,j ε ij /9, as is expected for evaporation along all 3-dimensions.The dipoles were assumed aligned along Θ = 90 • , and γ th interpolated over several temperatures to solve Eq. ( 1).
Realistically, forced evaporation by trap depth lowering tends to occur primarily along 1 direction, reducing the evaporation efficiency in the presence of molecular losses [45].The resulting out-of-equilibrium momentum distribution from single axis evaporation will be much like that in cross-dimensional rethermalization experiments, where an anisotropic collisional efficiency could now be used to your advantage.For instance, near unity collisional efficiency is achieved in the threshold regime with ε xz specifically at Θ = 45 • .Optimal evaporation protocols could thus be engineered by varying the molecular dipole orientation relative to the axis of evaporation.We leave such investigations to a future work.
Outlook and conclusions-By constructing a GP model of the elastic differential cross section between microwave shielded polar molecular fermions, we have found that non-threshold collisions can greatly diminish the efficacy of collisions toward thermalization of a nondegenerate gas.It is thus prudent to perform evaporation in the threshold regime, with the caveat that Pauli blocking in fermions would also lower the collisional efficiency below the Fermi temperature [21].If deployed in direct simulation Monte Carlo solvers [27,28,34], this GP model could also permit accurate dynamical studies in the Fermi degenerate or hydrodynamic regimes.The latter is motivated by restrictions of ε ij , only being able to describe thermalization in dilute samples.With larger molecular dipoles at densities required to achieve quantum degeneracy, the collision rate is far exceeded by the mean trapping frequency, demanding equilibration of trapped dipolar gases be treated within a hydrodynamic framework [46][47][48][49].The method of GP interpolation proposed here could similarly be applied to DC field shielded molecules [50]  For 2 polar molecules scattering of the effective potential V eff (r) provided in the main text, scattering solutions can be obtained by first expanding the wavefunction in the basis where Y ℓ,m ℓ (θ, ϕ) are spherical harmonics, and u E,ℓ,m ℓ (r) are solutions to the radial time-independent Schrödinger equation: Above, k 2 = 2µE/ℏ 2 is the collision wavenumber, and the explicit matrix elements ⟨ℓ, Numerical scattering solutions associated to Eq. ( 2), require picking a consistent convention when referencing the associated scattering matrices.We present our adopted convention as follows.First defining the matrices and the fundamental set of radial wavefunction solutions U (r; E), Eq. ( 2) can be recast as the compact system of equations: In principle, these equations can be solved numerically at a given collision energy E, by propagating the log derivative matrix from r = 0 to r → ∞.In practice, however, propagating to ∞ is not possible so we only do so up to r = r match , then match Y (r) to the asymptotic solutions where the distant colliders no longer interact.Moreover, we side step the issue of singularities at the origin by imposing a short-range boundary condition by starting the propagation at a minimum radius r = r min , then initializing the diagonal log-derivative matrix there as [32,52] that assumes universal short-range loss.This boundary condition prevents dipolar scattering resonances [53,54] which simplifies our current study.Propagation is done with an adaptive radial step size version of Johnson's algorithm [55], utilizing where a 0 is the Bohr radius and L is the largest value of ℓ utilized in the calculation.Typically, we utilize L = 121 or as many as is required for numerical convergence.The asymptotic solutions to Eq. ( 2) arise by considering the domain where r is much larger than the range of the potential, so that Eq. ( 2) is well approximated as This asymptotic radial equation is solved by the 2 independent solutions (up to arbitrary normalization): where j ℓ (kr) and n ℓ (kr) are the spherical Bessel and Neumann functions respectively.Then defining the matrices arbitrary solutions to Eq. ( 8), and in fact Eq. ( 2), can be written as where K is the reactance matrix that is responsible for matching the numerical scattering solutions U to the asymptotic solutions in Eq. ( 10) at r = r match .In particular, the off-diagonal elements of K provide information on the channel couplings that arise due to the interaction potential for a given incident collision channel.The matrix N is relevant only for normalization.
The reactance matrix can be written in terms of the logarithmic derivative via from which, we can then compute the other scattering matrices via the relations [25] The scattering matrices above permit us to evaluate the scattering amplitude, noting that m ℓ remains a good quantum number in these collision (App.I A), as which gives the appropriately antisymmetrized elastic differential cross section [26] via total cross section [52] and integral total cross section A.

Matrix elements of effective potential
To perform the scattering calculations on the effective single-channel microwave shielded potential energy surface, we are required to compute the ⟨ℓ, m| V eff (r) |ℓ ′ , m ′ ℓ ⟩ matrix elements.We list these elements explicitly in this section.The matrix elements for V dd (r) are given as while the matrix elements for V 6 (r) are given as

II. FRAME TRANSFORMATIONS FOR GAUSSIAN PROCESS FITTING
For efficient GP fitting of the elastic differential cross section, it is optimal to choose a coordinate frame whereby the symmetries are most conveniently handled.Naively, the differential cross section during a two-body collision involves 3 unit vectors: 1) the dipole alignment axis Ê, 2) the incident relative momentum k and 3) the outgoing relative momentum k ′ , therefore requiring 6 angular coordinates.This description is the case in a lab-frame (LF), where without loss of generality, we define it such that the dipole axis lies in its x, z-plane ÊLF = (sin Θ, 0, cos Θ) T , and the other 2 vectors are given in terms of spherical coordinates as However, we can also define a dipole-frame (DF) which utilizes the dipole alignment direction as its z-axis, Ê = ẑDF , while its x axis is aligned to the plane in which both Ê and k lie, so that The remaining xDF axis is then obtained with the cross product xDF = ŷDF × ẑDF .In the event where Ê and k coincide, we simply choose and also ŷDF = Ê × xDF .The differential cross section only cares about the relative angle between k and Ê, but not the vectors themselves.The dipole frame allows a convenient handling of this fact, we can simply write k = (sin η, 0, cos η) T where η = cos −1 k • Ê.So to obtain the post-collision relative momentum in the dipole-frame, we can construct the required rotation matrix R(DF ← LF), by the method of direction cosines The outbound relative collision vector is then given in dipole-frame as k′ (DF) = R(DF ← LF) k′ (LF).We then denote inclination and azimuthal scattering angles in the dipole frame as θ s and ϕ s respectively.Furthermore, the dipole frame as defined makes the differential cross section symmetric about the x, z-plane, only requiring us to specify ϕ s within the range [0, π].The entire differential cross section can then be obtained by specifying its value in the appropriate energy interval, and for η, θ s , ϕ s ∈ [0, π].

III. THE EIKONAL APPROXIMATION
At collision energies much larger than E dd , the scattering becomes semiclassical with the total cross section well approximated by the Eikonal approximation [25].Within this approximation, the scattering amplitude, considering on the 1/r 3 long-range tail of V eff , is given by is the momentum transfer vector, b is the impact parameter and tildes denote adimensional quantities normalized by the relevant dipole units (see the main text).
From the scattering amplitude, we can compute the total cross section using the optical theorem which requires evaluation of the scattering amplitude at forward scattering q = 0: where ϵ m = 1 if m = 0, and ϵ m = 2 if m > 0, and we utilized the substitution ℓ = kb.Then plugging the forward scattering amplitude into the optical theorem gives which averaged over incident directions, gives identical to the formula obtained for point-dipole scattering in Ref. [25].The result above is applicable to distinguishable dipoles, and so may not be quantitatively accurate in describing our study of scattering between identical fermions.Nevertheless, it serves to provide a useful visual aid to the expectation energy scaling, and seems to actually give rather favorable quantitative agreement.

IV. DERIVING THE COLLISIONAL EFFICIENCY TOWARD THERMALIZATION
Obtaining the form of N ij in the main text requires formulation of the Enskog equations.To do so, we define the phase space averaged quantity ⟨χ i ⟩ = k B (T i − T eq ), which quantifies the system's deviation from its equilibration temperature, T eq .Then multiplying the Boltzmann equation [56] by χ i and integrating it over phase space variables [28], we derived Enskog equations that govern the relaxation of ⟨χ j ⟩: where ⟨n⟩ is the average number density, c r (p r , t) is the distribution of relative momentum p r , and ∆χ ≡ χ ′ +χ ′ 1 −χ−χ 1 denotes the amount by which χ changes during a collision event.
Taken only perturbatively from equilibrium along axis i, Eq. 28a is approximated by the decay law C[χ j ] ≈ −γ ij ⟨χ j ⟩, which results in the short-time relation we move to center of mass and relative momentum coordinates, P = (p + p 1 )/2 and p r = p − p 1 respectively, so that the change in χ i is given as As for the relative momentum distribution, we Taylor expand it at t = 0 with respect to δ i /k B T 0 to get

FIG. 2 .
FIG. 2. Energy dependence of the angular averaged total cross section σ between microwave shielded 23 Na 40 K (black dashed line).The energy dependence clearly differs from the total cross section between fermionic point dipoles (dotted blue curve).For comparison, we plot the low energy Born and high energy Eikonal approximations with solid red lines.The inset shows a y = 0 slice of the effective microwave shielding interaction potential, with Rabi frequency Ω = 2π × 15 MHz and microwave detuning ∆ = 2π × 9.5 MHz.The shielding core is depicted as a white patch surrounding the coordinate origin which saturates the colorbar at V eff > 200E dd .Coordinate axes are plotted in units of 10 3 Bohr radii a0.

FIG. 3 .
FIG. 3. The central plot shows the total cross section as a function of the incident collision angle, obtained from (a) the Born approximation (red dashed curve), and from GP interpolation (solid curves) for 3 different collision energies: (b) E = 0.2E dd (black), (c) E = 2E dd (gray) and (d) E = 20E dd (light gray).In alphabetical correspondence, are angular plots of the differential cross section (in units of a 2 d ) in subplots with the respective collision energies, assuming dipoles pointing along Ê = ẑ and incident collision angle η = 45 • lying in the x, z-plane.Subplot (d) uses a smaller domain for clarity of presentation.

2 d 3 = k B ⟨n⟩ 2 d 3
identifying γ ij as the thermalization rate.Now considering the collision integralC[χ j ] = ⟨n⟩ p r m p r c r (p r , t) d 2 Ω ′ D el ∆χ j , p r m p r c r (p r , t) d 2 Ω ′ D el ∆T j ,