### Sample preparation

Table of Contents

To create our molecular samples, we first prepared a density-matched double-degenerate mixture of ^{23}Na and ^{40}K atoms. The atoms were subsequently associated to weakly bound molecules by means of a magnetic Feshbach resonance. Finally, the molecules were transferred to their absolute ground state via STIRAP. Details about the preparation process are described in refs. ^{18,43}. At the beginning of the measurements described in the main text, the molecules were trapped by the 1,064-nm and the 1,550-nm beam shown in Fig. 1a at a d.c. magnetic field of 72.35 G.

For the measurements of the collision rates, the microwave transition strength and to characterize the one-body loss, we worked with thermal molecules and sometimes reduced the molecule number to suppress interactions. For the collision rate measurements, the trap frequencies were (*ω*_{x}, *ω*_{y}, *ω*_{z}) = 2π × (67, 99, 244) Hz. For the evaporation, however, we started with near-degenerate molecules at (*ω*_{x}, *ω*_{y}, *ω*_{z}) = 2π × (45, 67, 157) Hz and ended up, for example, at (*ω*_{x}, *ω*_{y}, *ω*_{z}) = 2π × (52, 72, 157) Hz in case I or at (*ω*_{x}, *ω*_{y}, *ω*_{z}) = 2π × (42, 56, 99) Hz in case II and case III (Fig. 3).

To measure the cross-dimensional thermalization, we heated the weakly bound molecules along the vertical direction after we separated them from unbound atoms and before STIRAP was applied. For this purpose, we used parametric heating by modulating the intensity of the 1,064-nm beam at twice the vertical trap frequency.

### Microwave-field generation

It is essential that the phase noise of the microwave source does not induce transitions between the dressed states. We generated the microwave with a vector signal generator (Keysight E8267D). The microwave passes through a voltage-controlled attenuator (General Microwave D1954) before it is amplified with a 10-W power amplifier (KUHNE electronic KU PA 510590 – 10 A). At 10-MHz carrier offset, we measured −150 dBc Hz^{−1} phase-noise density from the signal generator and no significant enhancement from the amplifier. The microwave is emitted by a five-turn helical antenna (customized by Causemann Flugmodellbau) whose top end is about 2.2 cm away from the molecular sample.

With the voltage-controlled attenuator, we can adiabatically prepare the molecules in the dressed state by ramping the power attenuation linearly within 100 μs over a range of 65 dB.

### Imaging and thermometry

To image the molecules, we transferred them back into the non-dressed absolute ground state by ramping down the microwave power. Subsequently, the dipole traps were turned off and return STIRAP pulses were applied to bring the molecules back into the weakly bound state. After time of flight, typically 10 ms, the atoms were dissociated by ramping the magnetic field back over the Feshbach resonance. The magnetic field has to cross the Feshbach resonance slowly to minimize the release energy. In the end, the dissociated molecules were imaged by absorption imaging. We estimated that the derived temperature of the molecular sample could be overestimated by about 7 nK owing to the residual release energy. It is noted that the values of *T* and *T*/*T*_{F} reported in the main text do not account for the release energy.

To obtain the temperature of the molecular sample, we fit the absorption images with the Fermi–Dirac distribution

$${n}_{{\rm{FD}}}(x,z)={n}_{{\rm{FD}},0}{{\rm{Li}}}_{2}\left[-\zeta \exp \left(-\frac{{x}^{2}}{2{\sigma }_{x}^{2}}-\frac{{z}^{2}}{2{\sigma }_{z}^{2}}\right)\right],$$

(3)

where *n*_{FD,0} is the peak density, Li_{2}(*x*) is the dilogarithmic function, *ζ* is the fugacity and *σ*_{i=x,z} are the cloud widths in the *x* and *z* directions. Given a cloud width σ_{i}, we can calculate the temperature *T*_{i} by

$${\sigma }_{i}=\frac{\sqrt{1+{\omega }_{i}^{2}{t}_{{\rm{TOF}}}^{2}}}{{\omega }_{i}}\sqrt{\frac{{k}_{{\rm{B}}}{T}_{i}}{m}},$$

(4)

where *ω*_{i} is the trapping frequency in the *i*th direction, *t*_{TOF} is the time of flight and *m* is the mass of the molecules. The fugacity can be associated with the ratio of the temperature *T* and the Fermi temperature *T*_{F} with the relation

$${\left(\frac{T}{{T}_{{\rm{F}}}}\right)}^{3}=-\frac{1}{6{{\rm{Li}}}_{3}(\,-\,\zeta )},$$

(5)

where Li_{3}(*x*) is the trilogarithmic function. \({T}_{{\rm{F}}}={(6N)}^{1/3}\hbar \bar{\omega }/{k}_{{\rm{B}}}\) is given by the molecule number *N* and the geometric mean trap frequency \(\bar{\omega }={({\omega }_{x}{\omega }_{y}{\omega }_{z})}^{1/3}\). By rewriting *ζ* and fixing *T*_{F}, we are left with only the fitting parameters *n*_{FD,0}, *T*_{x} and *T*_{z}. We note that the temperature in the direction of the imaging beam *T*_{y} is assumed to be equal to *T*_{x} = *T*_{h}.

In addition, we independently determine the temperatures of the molecular samples from the time-of-flight images by fitting the thermal wings of the cloud to a Gaussian distribution

$${n}_{{\rm{th}}}(x,z)={n}_{{\rm{th}},0}\exp \,(\,-\,\frac{{x}^{2}}{2{\sigma }_{x}^{2}}-\frac{{z}^{2}}{2{\sigma }_{z}^{2}}),$$

(6)

where *n*_{th,0} is the peak density. Similar to ref. ^{33}, we first fit a Gaussian distribution to the whole cloud. We then constrain the Gaussian distribution to the thermal wings of the cloud by excluding a region of 1.5*σ* around the centre of the image. We find that by excluding 1.5*σ*, the ratio of signal to noise allows for the fit to converge for all datasets in Fig. 3a.

The temperatures extracted from fitting the Fermi–Dirac distribution and fitting the Gaussian distribution to the thermal wings are compared in Extended Data Fig. 1.

### Model for elastic and inelastic collisions

The elastic and inelastic collision rate coefficients *β*_{el} and *β*_{in} are experimentally determined from the time evolution of the measured molecule number *N*, the average temperature (2*T*_{h} + *T*_{v})/3 and the differential temperature *T*_{v} − *T*_{h} by numerically solving the differential equations^{19,34}

$$\frac{{\rm{d}}N}{{\rm{d}}t}=\left(-K\frac{2{T}_{{\rm{h}}}+{T}_{{\rm{v}}}}{3}n-{\varGamma }_{1}\right)N,$$

(7)

$$\frac{{\rm{d}}{T}_{{\rm{h}}}}{{\rm{d}}t}=\frac{1}{12}K{T}_{{\rm{v}}}{T}_{{\rm{h}}}n+\frac{{\varGamma }_{{\rm{th}}}}{3}({T}_{{\rm{v}}}-{T}_{{\rm{h}}}),$$

(8)

$$\frac{{\rm{d}}{T}_{{\rm{v}}}}{{\rm{d}}t}=\frac{1}{12}K(2{T}_{{\rm{h}}}-{T}_{{\rm{v}}}){T}_{{\rm{v}}}n-2\frac{{\varGamma }_{{\rm{th}}}}{3}({T}_{{\rm{v}}}-{T}_{{\rm{h}}}),$$

(9)

with the mean density

$$n=\frac{N}{8\sqrt{{{\rm{\pi }}}^{3}{k}_{{\rm{B}}}^{3}{T}_{{\rm{h}}}^{2}{T}_{{\rm{v}}}/{m}^{3}{\bar{\omega }}^{6}}}.$$

(10)

Here, *K* is the temperature-independent two-body loss coefficient, averaged for simplicity over all collision angles, and

$${\varGamma }_{{\rm{th}}}=\frac{n{\sigma }_{{\rm{el}}}v}{{N}_{{\rm{col}}}}$$

(11)

is the rethermalization rate with the elastic scattering cross-section *σ*_{el} and the thermally averaged collision velocity

$$v=\sqrt{16{k}_{{\rm{B}}}(2{T}_{{\rm{h}}}+{T}_{{\rm{v}}})/(3{\rm{\pi }}m)}.$$

(12)

The average number of elastic collisions per rethermalization is taken from ref. ^{48} as

$${N}_{{\rm{col}}}={\bar{{\mathscr{N}}}}_{z}(\varphi )=\frac{112}{45+4\,\cos (2\varphi )-17\,\cos (4\varphi )}$$

(13)

where *ϕ* is the tilt of the dipoles in the trap, which, in our case, corresponds to the tilt of the microwave wave vector with respect to the d.c. magnetic field. Following our characterization of the microwave polarization, we assume \({\bar{{\mathscr{N}}}}_{z}({29}^{\circ })=2.05\).

The anti-evaporation terms, that is, the first terms in equations (8) and (9), assume a linear scaling of the two-body loss rate with temperature. Our calculations predict that this assumption does not hold for small detunings (*Δ* < 2π × 20 MHz), as illustrated in Fig. 2. Our results, however, do not significantly change when we instead assume no temperature dependence in this regime.

Finally, after determining *σ*_{el} and *K*, the elastic and inelastic collision rate coefficients

$${\beta }_{{\rm{el}}}={\sigma }_{{\rm{el}}}v$$

(14)

and

$${\beta }_{{\rm{in}}}=K(2{T}_{{\rm{h}}}+{T}_{{\rm{v}}})/3$$

(15)

are plotted in Fig. 2 assuming a fixed temperature *T* = *T*_{h} = *T*_{v}.

Example data of the loss measurements, performed to determine *β*_{in}, are shown in Extended Data Fig. 2a. At high densities, two-body loss is the dominant contribution, whereas at low densities, the exponential shape of the loss curve shows that one-body effects outweigh inelastic collisions. To limit the number of free-fit parameters, we determine *Γ*_{1} = 1.7(4) Hz in independent measurements at low densities, as shown in Extended Data Fig. 2b. To suppress confounding effects from inelastic collisions, we reduce the initial molecule number to about 2,000 for these measurements. Under these conditions, the 1/e lifetime is 570(100) ms without shielding, which is still mostly limited by residual two-body collisions. Turning on the shielding results in a similar 1/e lifetime of about 590(100) ms. The lifetime reduces to 300(50) ms when a microwave source with a 3 dB higher phase-noise density (Rohde & Schwarz SMF100A) is used. If we isolatethe molecules by loading them into a three-dimensional optical lattice, we measure a lifetime of 8.0(1.2) s in absence of a microwave field, as shown in Extended Data Fig. 2c. Turning on the microwave field (using the Rohde & Schwarz SMF100A) results in a fast exponential decay to about half of the initially detected molecules, followed by a slow exponential decay. Assuming that the particles are isolated on individual lattice sites and that the faster decay is a result of mixing of two dressed states by phase noise of the microwave, we fit the data with the function

$$N(t)={N}_{0}{{\rm{e}}}^{-t/{\tau }_{0}}\left(\frac{1}{2}+\frac{1}{2}{{\rm{e}}}^{-t/{\tau }_{{\rm{MW}}}}\right),$$

(16)

where *N*_{0} is the initial number of molecules. We find a one-body loss time *τ*_{0} = 4.4(1.4) s and a state-mixing time *τ*_{MW} = 210(90) ms, which is in reasonable agreement with the lifetime measurements in the bulk. The slightly faster decay might be caused by molecules in higher bands, leading to residual collisions in the lattice, which are not accounted for in equation (16). The scaling of the lifetime with the microwave phase noise in the bulk and the measurements in the lattice indicate that *Γ*_{1} is currently limited by the noise power spectral density in the dressed-state transition (that is, at around 2π × 10 MHz offset from the carrier), even at a level of −150 dBc Hz^{−1} (ref. ^{38}). In addition, we find that the 1,550-nm light, which is used to trap the molecules in the bulk (Fig. 1a), contributes with about 0.5 Hz to the one-body decay. The underlying loss mechanism is under investigation.