# Many-body Chemical Reactions in a Quantum Degenerate Gas

Zhendong Zhang1 Shu Nagata1 Kaixuan Yao1 Cheng Chin1 1James Franck Institute,
Enrico Fermi Institute and Department of Physics, University of Chicago,
Chicago, Illinois 60637, USA

###### Abstract

Chemical reactions in the quantum degenerate regime are described by mixing of
matterwave fields. Quantum coherence and bosonic enhancement are two unique
features of many-body reactions involving bosonic reactants and products. Such
collective reactions of chemicals, dubbed “super-chemistry”, is an elusive
goal in quantum chemistry research. Here we report the observation of coherent
and collective reactive coupling between Bose condensed atoms and molecules
near a Feshbach resonance. Starting from an atomic condensate, the reaction
begins with a rapid formation of molecules, followed by oscillations of their
populations in the equilibration process. Faster oscillations are observed in
samples with higher densities, indicating bosonic enhancement. We present a
quantum field model which describes the dynamics well and identifies three-
body recombination as the dominant reaction process. Our findings exemplify
the highly sought-after quantum many-body chemistry and offer a new paradigm
for the control of quantum chemical reactions.

Ultracold atoms and molecules form an ideal platform toward controlling
chemical reactions at the level of single internal and external quantum
states. Ultracold molecules can be prepared in an individual internal state
by, e.g., magneto- Julienne2006_2  and photoassociation Julienne2006_1  of
ultracold atoms and direct laser cooling Shuman2010 . External motion of
molecules can be constrained by loading them into optical lattices Yan2013  or
tweezers Cairncross2021 . These experiments lead to the realization of state-
to-state ultracold chemistry Joschka2017 ; Rui2017 ; Liu2021 .

A number of experiments on cold molecules have reached the regime of quantum
degeneracy, which promise new forms of molecular quantum matter and reaction
dynamics. For instance, molecular Bose-Einstein condensates (BECs) formed in
atomic Fermi gases have stimulated tremendous interest in the BEC-BCS
(Bardeen-Cooper-Schrieffer) crossover Qijin2005 ; Giorgini2008 . Degenerate
fermionic molecules are created by magneto-association of bosonic and
fermionic atoms and optical transitions to the lowest rovibrational state
Marco2019 . Here quantum degeneracy suppresses chemical reactions due to the
fermion anti-bunching effect Marco2019 .

Figure 1: Reactive coupling between atomic and molecular quantum fields. a,
Bose condensed atoms described by a single wavefunction $\psi_{a}$ are coupled
to molecules condensed in the state $\psi_{m}$. The coupling synthesizes and
decomposes molecules. Wavy lines represent dissipation. b, We introduce a
reaction potential $V$ to describe the many-body dynamics of the atomic and
molecular fields. A pure sample of atoms or molecules first relaxes towards
lower potential and then equilibrates near the potential minimum. Due to
bosonic stimulation, the potential scales as $V\propto N^{\alpha}$, where $N$
is the total particle number and $\alpha$ is the reaction order, see text.

Recently, molecular BECs are realized based on atomic BECs near a Fesbach
resonance Zhang2021 . The reactive coupling between condensed atoms and
molecules promises a new regime of quantum chemistry, dubbed ?quantum super-
chemistry?, which highlights the coherent coupling of macroscopic matterwaves
and Bose stimulation of the reaction process Heinzen2000 ; Malla2021 . A key
feature of the coherence is the collective oscillations between the reactant
and product populations. Because of Bose statistics, enhancement of the
reaction dynamics is anticipated to significantly modify the branching ratio
Moore2002 .

At quantum degeneracy, reaction dynamics come fundamentally from mixing of the
matterwave fields of reactants and products. For instance, consider Feshbach
coupling which converts two atoms into one molecule and vice versa, described
by the chemical equation:
$\mathrm{A}+\mathrm{A}\longleftrightarrow\mathrm{A}_{2}$. In a quantum gas,
the reaction is described by the many-body Hamiltonian with reaction order
$\alpha=3$:

$\displaystyle\hat{H}=\varepsilon_{m}\hat{\psi}_{m}^{\dagger}\hat{\psi}_{m}+g_{2}\hat{\psi}_{m}^{\dagger}\hat{\psi}_{a}^{2}+g_{2}\hat{\psi}_{a}^{\dagger
2}\hat{\psi}_{m},$ (1)

where $\hat{\psi}_{a}$ ($\hat{\psi}_{m}$) is the atomic (molecular) field
operator, $g_{2}$ is the Feshbach coupling strength and $\varepsilon_{m}$ is
the energy of one bare molecule relative to two bare atoms. Here we define the
reaction order $\alpha$ as the maximum number of field operators in the
reaction terms.

Another prominent example that couples ultracold atoms and molecules is three-
body recombination, where three colliding atoms are converted into a diatomic
molecule and another atom, and vice versa. This process is described by the
chemical equation
$\mathrm{A}+\mathrm{A}+\mathrm{A}\longleftrightarrow\mathrm{A}_{2}+\mathrm{A}$.
At quantum degeneracy, the recombination process can resonantly couple atomic
and molecular fields as

$\displaystyle\hat{h}_{\mathrm{3}}=g_{3}\hat{\psi}_{m}^{\dagger}\hat{\psi}_{a}^{\dagger}\hat{\psi}_{a}^{3}+g_{3}\hat{\psi}_{a}^{\dagger
3}\hat{\psi}_{a}\hat{\psi}_{m},$ (2)

where $g_{3}$ is the recombination coupling strength. Here the reaction order
is $\alpha=5$.

To understand the dynamics of the coupled quantum fields, we present the
following picture. We show that the molecular population
$\hat{N}_{m}=\hat{\psi}_{m}^{\dagger}\hat{\psi}_{m}$ follows the form of an
“energy conservation” law as

$\displaystyle\frac{\hbar^{2}}{2}\dot{\hat{N}}_{m}^{2}+\hat{V}(\hat{N}_{m})=const.,$
(3)

where $\hbar^{2}\dot{\hat{N}}_{m}^{2}/2$ resembles the kinetic energy and we
introduce the many-body reaction potential
$\hat{V}=[\hat{N}_{m},\hat{H}]^{2}/2+const.$ (see Supplement). In this
picture, the system tends towards lower potential. Quantum fluctuations of the
nonlinear field coupling, however, can effectively damp the dynamics of the
populations Vardi2001 ; Richter_2015 . In experiments, damping can also come
from inelastic scattering and coupling to a thermal field. Thus one expects
that the system first relaxes towards the potential minimum, and then
equilibrates near the minimum with small amplitude coherent oscillations (see
Fig. 1). In the thermodynamic limit with total particle number $N\gg 1$, the
reaction potential and the oscillation frequency near the minimum scale with
the particle number as $V\propto N^{\alpha}$ and $\omega_{0}\propto
N^{\alpha/2-1}$. The dependence on the particle number signals bosonic
enhancement of the reaction dynamics Heinzen2000 ; Vardi2001 .

In this paper, we report the observation of coherent and Bose stimulated
reactions between Bose condensed Cs atoms and Cs2 molecules. The reaction is
initiated by tuning the magnetic field near a narrow $g-$wave Feshbach
resonance, which couples scattering atoms and diatomic molecules in a single
high-lying rovibrational state (See Supplement). Near the resonance, atomic
and molecular populations quickly relax toward a dynamical equilibrium,
followed by coherent oscillations between atoms and molecules in the
equilibration process. We show that the oscillation frequency strongly depends
on the particle number. From the dependence, we conclude that three-body
recombination is the dominant reaction process that couples the atomic and
molecular fields near the Feshbach resonance.

Figure 2: Comparison of molecule formation rate in classical and quantum
degenerate regimes. a, Dynamics of molecule formation in an atomic gas after
quenching the magnetic field $3(2)$ mG above the Feshbach resonance at
$B_{0}=19.849(2)$ G. The solid lines are fits to the data in the initial
growth stage for extraction of the molecule formation rate
$\dot{N}_{\mathrm{m}}$. b, The extracted molecule formation rate coefficient
$\beta$ above and below the critical temperature $T_{c}$. The red line is a
fit to the data based on the classical kinetic theory prediction
$\beta=b_{\mathrm{Cl}}\Gamma_{\mathrm{Cl}}$, from which we obtain the
classical branching ratio $b_{\mathrm{Cl}}=7(1)\%$ (see text). The blue line
fits the data in the quantum regime with
$\beta=b_{\mathrm{Q}}\Gamma_{\mathrm{Q}}$, which gives the quantum branching
ratio $b_{\mathrm{Q}}=3.9(3)\%$ (see text). The inset shows the rate
coefficient normalized to the classical gas expectation
$\Gamma_{\mathrm{Cl}}$. In panel a error bars represent one standard deviation
of the mean, estimated from 4-8 measurements. In panel b error bars represent
95% confidence intervals of the mean. Figure 3: Coherent reaction dynamics in
quantum gases of atoms and molecules across a Feshbach resonance. a, Evolution
of atomic and molecular populations in an atomic BEC after quenching the
magnetic field $2(1)$ mG below the resonance. Solid lines are fits to capture
the dynamics in relaxation and equilibration processes. b, Effective
temperatures determined from time-of-flight measurements of the atoms and
molecules at the same field. Solid lines are guides to the eye. Insets in
panels a and b are sample images of atoms and molecules after the time-of-
flight (see Supplement). c, Loss rate of atoms immediately after the quench.
Solid (empty) circles represent samples prepared below (above) the resonance.
Green line is a Lorentzian fit with center at $B_{0}$ and width $\Delta B$.
Magenta solid (dashed) line is a fit near and below (above) the resonance
based on $\gamma_{a}=\gamma_{0}/[1+|(B-B_{0})/\delta B|^{\epsilon_{\pm}}]$
(see text), from which we obtain the exponent $\epsilon_{-}=6(2)$
($\epsilon_{+}=2.9(4)$) below (above) the resonance. Inset is a zoomed-in view
near the resonance. d, Oscillation frequency of molecular populations from
atomic samples at mean BEC density of $2.9\times 10^{13}\mathrm{cm}^{-3}$ and
BEC fraction of $80\%$ (red) and $60\%$ (purple). Solid lines are empirical
fits based on $\omega=\sqrt{\delta\mu^{2}(B-B_{m})^{2}+\omega_{0}^{2}}$, where
$B_{m}$ and $\omega_{0}$ are fitting parameters. The values of $B_{m}$ from
the fits are consistent with the resonance position $B_{0}$ within our
measurement uncertainty. Dashed lines are the asymptotes
$\omega=|\delta\mu(B-B_{m})|$. Data in panels a and b are averages of 3-4
measurements, and error bars represent one standard deviation of the mean.
Data in panels c and d are obtained from the fits, see Supplement, and error
bars represent 95% confidence intervals.

Our experiment starts with an ultracold Bose gas of $6\times 10^{3}$ to
$5\times 10^{5}$ cesium atoms in an optical trap. The atoms can form a pure
BEC either in a three-dimensional (3D) harmonic potential or a 2D square well
potential Clark2017 . We induce the reaction by switching the magnetic field
near the g-wave Feshbach resonance, which can convert an atomic BEC into a
molecular BEC Zhang2021 . We determine the resonance position
$B_{0}=19.849(2)$ G, resonance width $\Delta B=8.3(5)$ mG, and the relative
magnetic moment $\delta\mu=h\times 0.76(3)$ MHz/G, where $h$ is the Planck
constant, from measurements of the molecular binding energy
$\varepsilon_{m}\approx\delta\mu(B-B_{0})$ and the scattering length (see
Supplement). After the reaction, we decouple the atoms and molecules by
quickly tuning the magnetic field far off the resonance and image each
independently Zhang2021 .

To show that chemical reactions follow different rules in a degenerate quantum
gas versus in a normal gas, we compare the molecule production rate for
samples prepared above and below the BEC critical temperature $T_{c}$. We
extract the molecule production rate coefficient
$\beta=\dot{N}_{m}/N_{0}n_{0}$ right after the magnetic field switch, where
$N_{0}$ and $n_{0}$ are the initial total atom number and mean atomic density,
respectively, see Fig. 2.

The measured molecule formation rate shows distinct behaviour in the two
regimes. In a thermal gas with temperature $T>T_{c}$, the molecule formation
rate is $\beta=b_{\mathrm{Cl}}\Gamma_{\mathrm{Cl}}$, where
$\Gamma_{\mathrm{Cl}}$ and $b_{\mathrm{Cl}}$ are the classical atomic
collision rate coefficient and the branching ratio into the molecular state,
respectively. Near the resonance, the collision rate coefficient is unitarity
limited as $\Gamma_{\mathrm{Cl}}=8h^{2}(\pi^{3}m^{3}k_{\mathrm{B}}T)^{-1/2}$
Chin2010 , where $m$ is the atomic mass and $k_{\mathrm{B}}$ is the Boltzmann
constant. Our measurement in the thermal regime is consistent with the
$T^{-1/2}$ scaling. From the fit we extract the branching ratio
$b_{\mathrm{Cl}}=7(1)\%$ (see Fig. 2b).

Figure 4: Bose-enhanced atom-molecule reaction dynamics on Feshbach resonance.
a, Molecules formed in atomic BECs with different initial atom numbers $N_{0}$
following the magnetic field quench to $1(1)$ mG below the Feshbach resonance.
Solid lines are fits to the data. b, The extracted oscillation frequencies for
different initial atom number $N_{0}$. The red (magenta) solid line is a power
law fit with exponent given by the two-(three-)body model. The blue dashed
line is a power law fit with a free exponent, which yields the scaling
$\omega_{0}\propto N_{0}^{0.7(2)}$. In the inset, we show the frequency
dependence on the mean atomic density $n_{0}$ and the associated fits that
yields the scaling $\omega_{0}\propto n_{0}^{1.7(4)}$. c, The reaction
potential $V\equiv V_{3}\approx-4g_{3}^{2}N^{5}f_{m}(1-f_{m})^{4}$ of the
three-body process described in Eq. (2) for different total particle number
$N$ (green solid lines). In the thermodynamic limit $N\gg 1$, the minimum
occurs at molecule fraction $f_{m}\equiv 2N_{m}/N=1/5$. d, Evolution of
molecule fraction $f_{m}$ for different initial atom number $N_{0}$. The mean
molecule fractions in the equilibration phase $f_{m}=19(2)\%$ and
$f_{m}=21(1)\%$ are consistent with the predicted minimum position of $V_{3}$
at $f_{m}=1/5$. Here the uncertainties represent $95\%$ confidence intervals.
Data in panels a and d are averages of 3-5 measurements and error bars
represent one standard deviation of the mean. Error bars in panel b represent
95% confidence intervals.

Entering the quantum degenerate regime $T<T_{c}$, we observe a steep drop in
the rate coefficient, see Fig. 2b. At low temperatures, we model the rate
coefficient as $\beta=b_{\mathrm{Q}}\Gamma_{\mathrm{Q}}$, where
$b_{\mathrm{Q}}$ and $\Gamma_{\mathrm{Q}}=4\sqrt{3}(h/m)(6\pi^{5}n)^{-1/3}$
are the branching ratio and the rate coefficient predicted by the universal
theory in the quantum regime Jin2014 ; Eismann2016 . The model fits the
measurement well, and we extract the branching ratio to be
$b_{\mathrm{Q}}=3.9(3)\%$. The sharp transition of molecule formation rate
around the critical temperature $T_{c}$ indicates different laws in the
classical and quantum degenerate regimes.

Close examination of molecule formation dynamics in atomic BECs reveals
additional interesting features of quantum many-body reactions. To understand
the underlying reaction processes, we study the atom loss rate
$\gamma_{a}=-\dot{N}_{a}/N_{0}$ right after switching the magnetic field,
where $N_{a}$ is the atom number (see Fig. 3c). Far from the resonance
$|B-B_{0}|\gg\Delta B$, atoms decay slowly and the loss rate follows a
symmetric Lorentzian profile centered at the resonance
$\gamma_{a}\propto(B-B_{0})^{-2}$. We attribute the lineshape to the Feshbach
coupling Chin2010 .

Near the resonance, the loss rate greatly exceeds the expectation from the
Lorentzian profile. This rapid atom loss only lasts for a few 100 $\mu$s and
is accompanied with fast molecule production and heating of both atoms and
molecules. We identify this fast process as the relaxation dynamics described
in Fig. 1b. In order to characterize the enhanced reaction rate, we fit the
loss rate near the resonance as $\gamma_{a}\propto[1+|(B-B_{0})/\delta
B|^{\epsilon_{\pm}}]^{-1}$ from which we extract the exponents
$\epsilon_{+}=2.9(4)$ above the resonance and $\epsilon_{-}=6(2)$ below the
resonance. The exponents $\epsilon_{\pm}$ larger than 2 are consistent with
the enhanced atom loss near the resonance beyond the Lorentzian profile.

The relaxation dynamics stem from three-body recombination, evidenced by the
fast heating of both species in the relaxation phase, see Fig. 3b Greene2017 .
In addition, the measured exponent $\epsilon_{+}=2.9(4)$ from the enhanced
atom loss is consistent with the predicted value of $3.5$ for three-body
recombination near a narrow Feshbach resonance Petrov2004 . We attribute the
even larger exponent $\epsilon_{-}=6(2)$ below the resonance to bosonic
enhancement of the three-body process.

Following the relaxation, both atomic and molecular populations oscillate for
several ms before they slowly decay over a much longer time scale (see Fig.
3a,b and Supplement). The oscillation is consistent with the equilibration
dynamics near the reaction potential minimum described in Fig. 1b. The
frequency $\omega$ of the oscillation depends on the magnetic field and is
well fit to $\omega=\sqrt{\varepsilon_{m}^{2}/\hbar^{2}+\omega_{0}^{2}}$ (see
Fig. 3d). Far from the resonance, the frequency approaches the molecular
binding energy $|\varepsilon_{m}|$. On resonance with $\varepsilon_{m}=0$, the
frequency $\omega=\omega_{0}$ is given by the collective reactive coupling
between the atomic and molecular fields. Large-amplitude oscillations are also
observed in samples with magnetic field modulation, see Supplement.

To demonstrate the many-body nature of the reactive coupling, we probe the
atom-molecule oscillations right on the Feshbach resonance with different
initial atom number $N_{0}$ and mean density $n_{0}$. After quenching the
magnetic field, we observe that samples with higher populations and densities
display faster oscillations, see Fig. 4a. Fitting the data, we obtain the
scaling $\omega_{0}\propto N_{0}^{0.7(2)}$ and $\omega_{0}\propto
n_{0}^{1.7(4)}$, see Fig. 4b. Note that the two scalings are linked by
$n_{0}\propto N_{0}^{2/5}$ for BECs in a harmonic trap Pethick2008 . The
particle number dependence of the reactive coupling supports bosonic
enhancement of the reaction process.

The scaling with respect to the particle number also reveals the underlying
reaction mechanism. For the two-body process described in Eq. (1), we derive
the effective potential $V_{2}=-g_{2}^{2}N^{3}f_{m}(1-f_{m})^{2}+O(N^{2})$,
where $f_{m}=2N_{m}/N$ is the molecule fraction, from which the resonant
oscillation frequency is calculated to be $\omega_{0}\propto N^{1/5}$ in a
harmonic trap (see Supplement). For the three-body recombination process
described in Eq. (2), the effective potential is
$V_{3}=-g_{3}^{2}N^{5}f_{m}(1-f_{m})^{4}+O(N^{4})$, which yields the scaling
$\omega_{0}\propto N^{3/5}$. Our measurement agrees well with the three-body
model, see Fig. 4b.

Moreover, we find the molecule fraction oscillates around 20(1)% in the
equilibration phase, which is consistent with the minimum position of the
reaction potential $V_{3}$ at $f_{m}=1/5$ (see Fig. 4c,d). Two-body Feshbach
process, on the other hand, predicts a different minimum of $V_{2}$ at
$f_{m}=1/3$.

To conclude, we observe collective many-body chemical reactions in an atomic
BEC near a Feshbach resonance. The dynamics are well described by a quantum
field model derived from three-body recombination. In particular, the coherent
oscillations of atomic and molecular fields in the equilibration phase support
quantum coherence and Bose enhancement of the reaction process. The
observation of coherent and collective chemical reactions in the quantum
degenerate regime paves the way to explore the interplay between many-body
physics and ultracold chemistry.

Acknowledgement We thank P. Julienne, K. Levin, D. Mazziotti, D. DeMille and
K.-K. Ni for helpful discussions. We thank K. Patel and L. Weiss for carefully
reading the manuscript. We thank J. Jachinowski for experimental assistance
and carefully reading the manuscript. This work was supported by the National
Science Foundation under Grant No. PHY1511696 and PHY-2103542, by the Air
Force Office of Scientific Research under award number FA9550-21-1-0447. Z.Z.
is supported by the Grainger Graduate Fellowship. S.N. acknowledges support
from the Takenaka Scholarship Foundation.

## References

  * (1) Köhler, T., Góral, K. & Julienne, P. S.  Production of cold molecules via magnetically tunable Feshbach resonances.  _Rev. Mod. Phys._ 78, 1311–1361 (2006).
  * (2) Jones, K. M., Tiesinga, E., Lett, P. D. & Julienne, P. S.  Ultracold photoassociation spectroscopy: Long-range molecules and atomic scattering.  _Rev. Mod. Phys._ 78, 483–535 (2006).
  * (3) Shuman, E. S., Barry, J. F. & DeMille, D.  Laser cooling of a diatomic molecule.  _Nature_ 820–823 (2010).
  * (4) Yan, B. _et al._ Observation of dipolar spin-exchange interactions with lattice-confined polar molecules.  _Nature_ 521–525 (2013).
  * (5) Cairncross, W. B. _et al._ Assembly of a rovibrational ground state molecule in an optical tweezer.  _Phys. Rev. Lett._ 126, 123402 (2021).
  * (6) Wolf, J. _et al._ State-to-state chemistry for three-body recombination in an ultracold rubidium gas.  _Science_ 921–924 (2017).
  * (7) Rui, J. _et al._ Controlled state-to-state atom-exchange reaction in an ultracold atom–dimer mixture.  _Nature Physics_ 699–703 (2017).
  * (8) Liu, Y. _et al._ Precision test of statistical dynamics with state-to-state ultracold chemistry.  _Nature_ 379–384 (2021).
  * (9) Chen, Q., Stajic, J., Tan, S. & Levin, K.  BCS–BEC crossover: From high temperature superconductors to ultracold superfluids.  _Physics Reports_ 1–88 (2005).
  * (10) Giorgini, S., Pitaevskii, L. P. & Stringari, S.  Theory of ultracold atomic Fermi gases.  _Rev. Mod. Phys._ 80, 1215–1274 (2008).
  * (11) Marco, L. D. _et al._ A degenerate Fermi gas of polar molecules.  _Science_ 853–856 (2019).
  * (12) Zhang, Z., Chen, L., Yao, K.-X. & Chin, C.  Transition from an atomic to a molecular Bose–Einstein condensate.  _Nature_ 708–711 (2021).
  * (13) Heinzen, D. J., Wynar, R., Drummond, P. D. & Kheruntsyan, K. V.  Superchemistry: Dynamics of coupled atomic and molecular Bose-Einstein condensates.  _Phys. Rev. Lett._ 84, 5029–5033 (2000).
  * (14) Malla, R. K., Chernyak, V. Y., Sun, C. & Sinitsyn, N. A.  Coherent reaction between molecular and atomic Bose-Einstein condensates: Integrable model.  _Phys. Rev. Lett._ 129, 033201 (2022).
  * (15) Moore, M. G. & Vardi, A.  Bose-enhanced chemistry: Amplification of selectivity in the dissociation of molecular Bose-Einstein condensates.  _Phys. Rev. Lett._ 88, 160402 (2002).
  * (16) Vardi, A., Yurovsky, V. A. & Anglin, J. R.  Quantum effects on the dynamics of a two-mode atom-molecule Bose-Einstein condensate.  _Phys. Rev. A_ 64, 063611 (2001).
  * (17) Richter, F. _et al._ Ultracold chemistry and its reaction kinetics.  _New Journal of Physics_ 055005 (2015).
  * (18) Clark, L. W., Gaj, A., Feng, L. & Chin, C.  Collective emission of matter-wave jets from driven Bose–Einstein condensates.  _Nature_ 356–359 (2017).
  * (19) Chin, C., Grimm, R., Julienne, P. & Tiesinga, E.  Feshbach resonances in ultracold gases.  _Rev. Mod. Phys._ 82, 1225–1286 (2010).
  * (20) Greene, C. H., Giannakeas, P. & Pérez-Ríos, J.  Universal few-body physics and cluster formation.  _Rev. Mod. Phys._ 89, 035006 (2017).
  * (21) Makotyn, P., Klauss, C. E., Goldberger, D. L., Cornell, E. A. & Jin, D. S.  Universal dynamics of a degenerate unitary Bose gas.  _Nature Physics_ 116–119 (2014).
  * (22) Eismann, U. _et al._ Universal loss dynamics in a unitary Bose gas.  _Phys. Rev. X_ 6, 021025 (2016).
  * (23) Petrov, D. S.  Three-boson problem near a narrow Feshbach resonance.  _Phys. Rev. Lett._ 93, 143201 (2004).
  * (24) Chin, C. & Grimm, R.  Thermal equilibrium and efficient evaporation of an ultracold atom-molecule mixture.  _Phys. Rev. A_ 69, 033612 (2004).
  * (25) Pethick, C. J. & Smith, H.  _Bose–Einstein condensation in dilute gases_ (Cambridge university press, 2008).
  * (26) Hung, C.-L., Zhang, X., Gemelke, N. & Chin, C.  Accelerating evaporative cooling of atoms into Bose-Einstein condensation in optical traps.  _Phys. Rev. A_ 78, 011604 (2008).
  * (27) Chin, C. _et al._ Observation of Feshbach-like resonances in collisions between ultracold molecules.  _Phys. Rev. Lett._ 94, 123201 (2005).
  * (28) Zhang, Z.  _Coherent dynamics and reactions in atomic and molecular Bose-Einstein condensates_ (The University of Chicago, 2022).
  * (29) Thompson, S. T., Hodby, E. & Wieman, C. E.  Ultracold molecule production via a resonant oscillating magnetic field.  _Phys. Rev. Lett._ 95, 190404 (2005).
  * (30) Lange, A. D. _et al._ Determination of atomic scattering lengths from measurements of molecular binding energies near Feshbach resonances.  _Phys. Rev. A_ 79, 013622 (2009).
  * (31) Hung, C.-L.  _In situ probing of two-dimensional quantum gases_ (The University of Chicago, 2011).
  * (32) Hung, C.-L., Zhang, X., Gemelke, N. & Chin, C.  Observation of scale invariance and universality in two-dimensional Bose gases.  _Nature_ 236–239 2011).
  * (33) Castin, Y. & Dum, R.  Bose-Einstein condensates in time dependent traps.  _Phys. Rev. Lett._ 77, 5315–5319 (1996).
  * (34) Berninger, M. _et al._ Feshbach resonances, weakly bound molecular states, and coupled-channel potentials for cesium at high magnetic fields.  _Phys. Rev. A_ 87, 032517 (2013).
  * (35) Mark, M. J., Meinert, F., Lauber, K. & Nagerl, H.-C.  Mott-insulator-aided detection of ultra-narrow Feshbach resonances.  _SciPost Phys._ 5, 55 (2018).
  * (36) Barnett, S. M. & Knight, P. L.  Thermofield analysis of squeezing and statistical mixtures in quantum optics.  _J. Opt. Soc. Am. B_ 467–479 (1985).
  * (37) Kumar, S. & Mehta, C. L.  Theory of the interaction of a single-mode resonant radiation field with $n$ two-level atoms.  _Phys. Rev. A_ 21, 1573–1588 (1980).
  * (38) Ketterle, W., Durfee, D. S. & Stamper-Kurn, D.  Making, probing and understanding Bose-Einstein condensates.  _arXiv preprint cond-mat/9904034_ (1999).
  * (39) Claussen, N. R. _et al._ Very-high-precision bound-state spectroscopy near a ${}^{85}\mathrm{Rb}$ Feshbach resonance.  _Phys. Rev. A_ 67, 060701 (2003).
  * (40) Liu, B., Fu, L.-B. & Liu, J.  Shapiro-like resonance in ultracold molecule production via an oscillating magnetic field.  _Phys. Rev. A_ 81, 013602 (2010).
  * (41) Barone, A. & Paterno, G.  _Physics and applications of the Josephson effect_ , vol. 1 (Wiley Online Library, 1982).

Supplementary Material

##  I Experimental Procedure

Our experiment starts with a ultracold Bose gas of 6,000 to 470,000 133Cs
atoms at a temperature of 2 to 232 nK in a 3D harmonic trap. We tune the
temperature and atom number by changing the trap depth at the end of
evaporation process Hung2008 . The harmonic trap frequencies are
($\omega_{x}$,$\omega_{y}$,$\omega_{z}$) = $2\pi\times$(24,13,74) to
$2\pi\times$(36,15,91) Hz. The atoms are polarized into the hyperfine ground
state $|F=3,m_{\mathrm{F}}=3\rangle$, where $F$ and $m_{\mathrm{F}}$ are
quantum numbers for the total spin and its projection along the magnetic field
direction, respectively. The narrow g-wave Feshbach resonance couples
$\mathrm{Cs}$ atoms into $\mathrm{Cs}_{2}$ molecules at
$|f=4,m_{f}=4;l=4,m_{l}=2\rangle$, where $f$ and $l$ represent quantum numbers
for the sum of the spins of two individual atoms and the orbital angular
momentum of a molecule, $m_{f}$ and $m_{l}$ are projections of $f$ and $l$
along the magnetic field direction Chin2005 .

To induce the molecule formation dynamics, we quench the magnetic field close
to the resonance position $B_{0}$ from 19.5 G where the samples are prepared.
After holding for variable times, we switch the field back to either 19.5 G or
17.17 G to decouple atoms and molecules. We can image the remaining atoms at
this field by absorption imaging. We can also wait for the remaining atoms to
fly away after a resonant light pulse and image the molecules by jumping the
field up to 20.4 G to dissociate them into atoms and then image the atoms from
the dissociation Zhang2021 . For the atom loss measurements shown in Fig. 3c,
BECs with $\sim$40,000 atoms are transferred from the harmonic trap to a 2D
flat-bottomed optical potential before we quench the field to different values
near the resonance Clark2017 . For the rest of the data shown in Figs. 2-4, we
start from atomic samples in the 3D harmonic dipole trap.

To measure the temperature of atoms or molecules (e.g. as is shown in Fig.
3b), we release them into a horizontally isotropic harmonic trap for a quarter
of the trap period, which converts the particle distribution from the real
space to the momentum space zhang2022coherent . We extract the temperature $T$
by fitting the momentum distribution with the condensate around zero momentum
excluded using the Gaussian function
$n(k_{r})=n(0)\exp{[-\hbar^{2}k_{r}^{2}/(2mk_{\mathrm{B}}T)]}$, where $k_{r}$
is the radial wave number and $k_{\mathrm{B}}$ is the Boltzmann constant.

##  II Determination of the Feshbach resonance position and width

To determine the position of the narrow g-wave Feshbach resonance in our
system, we perform measurements of molecular binding energy at different
offset magnetic fields using magnetic field modulation spectroscopy
Thompson2005 ; Lange2009  and find the field value where the binding energy
reaches zero.

Figure S1: Bound state energy diagram for cesium atoms in the hyperfine ground
state $|F=3,m_{\mathrm{F}}=3\rangle$ and molecular energy measurement near the
g-wave Feshbach resonance around 20 G using modulation spectroscopy. a, Energy
diagram for $\mathrm{Cs}_{2}$ molecular states close to the atomic scattering
continuum adapted from Fig. 22 in Ref. Chin2010 . b, Molecular energy
$\varepsilon_{\mathrm{m}}$ obtained from modulation spectroscopy at different
offset magnetic fields. The solid line is a linear fit which reaches 0 at
$B_{0}=19.849(1)$ G.

We start with atomic BECs at 23 nK prepared at 19.5 G. Then we quench the
field to an offset value $B_{\mathrm{dc}}$ near the resonance and
simultaneously modulate the field sinusoidally with an amplitude
$B_{\mathrm{ac}}=5$ mG for 5 ms. We scan the modulation frequency and measure
the spectrum of the remaining atom number. From the atom loss peak of the
spectrum due to the conversion from atoms into molecules, we extract the
resonant frequency that corresponds to molecular energy at an offset magnetic
field $B_{\mathrm{dc}}$ near the g-wave Feshbach resonance zhang2022coherent ;
Thompson2005 , see Fig. S1. We have confirmed that the resulting atom loss
peak position is not sensitive to the modulation amplitude and modulation
time. A linear fit to the data in Fig. S1 gives the resonance position
$B_{0}=19.849(1)$ G where the molecular energy goes to zero. The slope of the
linear fit gives the magnetic moment difference between two bare atoms and one
bare molecule as $\delta\mu=h\times 0.76(3)$ MHz/G, which is consistent with
Ref. Chin2005 . We emphasize that for the narrow resonance we are using, the
molecular energy approaches zero quadratically only within a small fraction of
the resonance width. Our linear fit to the molecular energy data
underestimates the resonance position by $\sim 0.3$ mG based on our
calculation using the resonance width from the following scattering length
measurements Chin2010 . The systematic error of our calibration of the
absolute magnetic field is less than 20 mG. Throughout this work, we perform
the magnetic field calibration based on the same procedure to ensure a
constant systematic error.

Next we measure the s-wave scattering length near the resonance to obtain the
resonance width. Here the scattering length is inferred from the expansion of
a quasi-2D BEC prepared with trap frequencies
$(\omega_{x},\omega_{y},\omega_{z})=2\pi\times(11,13,895)$ Hz. During the
expansion, the mean field interaction energy is converted into kinetic energy.
We first prepare the BEC at an initial magnetic field $B_{i}=20.481$ G or
$19.498$ G where the scattering is $a_{i}$. The column density distribution of
atoms in the Thomas-Fermi regime is Hung2011 :

$\displaystyle
n(x,y)=[\mu-\frac{1}{2}m\omega_{x}^{2}x^{2}-\frac{1}{2}m\omega_{y}^{2}y^{2}]/g_{\mathrm{2D}},$
(S1)

where $g_{\mathrm{2D}}=(\hbar^{2}/m)\sqrt{8\pi}a_{i}/l_{z}$ is the coupling
strength, $l_{z}=\sqrt{\hbar/m\omega_{z}}$ is the harmonic oscillator length
in the tightly confined z direction and
$\mu=\sqrt{g_{\mathrm{2D}}Nm\omega_{x}\omega_{y}/\pi}$ is the chemical
potential determined by $g_{\mathrm{2D}}$, the total atom number $N$ and the
initial trap frequencies $\omega_{x}$ and $\omega_{y}$ in the horizontal
plane.

Figure S2: Scattering length measurement near the narrow g-wave Feshbach
resonance by time-of-flight. a, Atomic density distributions after 20 ms time-
of-flight at different magnetic fields near the Feshbach resonance. The images
with $B<19.865$ G ($B>19.865$ G) come from initial BECs prepared below (above)
the Feshbach resonance. b, Scattering length extracted from the Thomas-Fermi
radii in the time-of-flight images, see text. The circular (diamond) data
points come from initial BECs prepared below (above) the resonance. The solid
line is a fit to the data excluding the points at
$19.858~{}\mathrm{G}<B<19.909~{}\mathrm{G}$ based on Eq. (S4), from which we
obtain the resonance width $\Delta B=8.3(5)$ mG. The points at
$19.855~{}\mathrm{G}<B<19.909~{}\mathrm{G}$ are excluded because of the
heating effect near the resonance. c, Total atom number extracted from the
time-of-flight images.

Then we quench the magnetic field to a different value $B_{f}$ where the
scattering length is $a_{f}$ and simultaneously switch off the harmonic trap
in the horizontal plane. According to Ref. Castin1996 , the dynamics of a BEC
after the release follow a simple dilation with scaling parameters
$\lambda_{x}(t)$ and $\lambda_{y}(t)$, which determine the density
distribution at time $t$ as:

$\displaystyle
n(x,y,t)=\frac{\mu-\frac{1}{2}m\omega_{x}^{2}x^{2}/\lambda_{x}^{2}(t)-\frac{1}{2}m\omega_{y}^{2}y^{2}/\lambda_{y}^{2}(t)}{g_{\mathrm{2D}}\lambda_{x}(t)\lambda_{y}(t)},$
(S2)

where the scaling parameters evolve according to:

$\displaystyle\ddot{\lambda}_{x}(t)$
$\displaystyle=\frac{a_{f}}{a_{i}}\frac{\omega_{x}^{2}}{\lambda_{x}^{2}(t)\lambda_{y}(t)}$
$\displaystyle\ddot{\lambda}_{y}(t)$
$\displaystyle=\frac{a_{f}}{a_{i}}\frac{\omega_{y}^{2}}{\lambda_{y}^{2}(t)\lambda_{x}(t)}.$
(S3)

We scan the magnetic field and measure the Thomas-Fermi radii
$R_{j}=\sqrt{2\mu\lambda_{j}^{2}(t)/m\omega^{2}_{j}}$ where $j=x,y$ after 20
ms expansion. Eventually we extract $a_{f}$ based on its one-to-one
correspondence to the Thomas-Fermi radii according to Eq. (S3). The results
are summarized in Fig. S2 and we fit the scattering length data using the
formula Chin2010

$\displaystyle a(B)=a_{\mathrm{bg}}[1+\eta(B-B_{0})](1-\frac{\Delta
B}{B-B_{0}}),$ (S4)

where we obtain the resonance width $\Delta B=8.3(5)$ mG, the resonance
position $B_{0}=19.861(1)$ G, the background scattering length on resonance
$a_{\mathrm{bg}}=163(1)a_{0}$ and the slope of the background scattering
length $\eta=0.31(2)$/G. The background scattering length $a_{\mathrm{bg}}$
and the slope $\eta$ are consistent with Ref. Berninger2013  and the resonance
width $\Delta B$ is consistent with Ref. Mark2018 , where a different method
is used. The fitted resonance position deviates from that in the binding
energy measurement by $\sim$10 mG, which we attribute to the heating of atoms
near the resonance. The binding energy measurement, however, suffers less from
the heating issue Thompson2005 . Throughout the whole paper, we adopt the
resonance position $B_{0}=19.849(1)$ G from the binding energy measurement.

##  III Many-body chemical reactions in quantum degenerate regime

Here we derive the effective potential $V$ for three-body and two-body
processes and calculate the minimum position of $V$ and the oscillation
frequency of the system near the minimum. We start from the many-body
Hamiltonian for the three-body recombination process
$\mathrm{A}+\mathrm{A}+\mathrm{A}\longleftrightarrow\mathrm{A}_{2}+\mathrm{A}$
near a Feshbach resonance,

$\displaystyle\hat{H}=\varepsilon_{m}\hat{\psi}_{m}^{\dagger}\hat{\psi}_{m}+g_{3}(\hat{\psi}_{m}^{\dagger}\hat{\psi}_{a}^{\dagger}\hat{\psi}_{a}^{3}+\hat{\psi}_{a}^{\dagger
3}\hat{\psi}_{a}\hat{\psi}_{m}),$ (S5)

where we only consider the ground states for atoms and molecules since their
associated coupling terms dominate over others that involve excited states,
due to the macroscopic population in the ground state of atoms that we begin
with.

To characterize the evolution of the system, we derive the equation of motion
for the atomic population $\hat{N}_{a}=\hat{\psi}_{a}^{\dagger}\hat{\psi}_{a}$
and molecular population $\hat{N}_{m}=\hat{\psi}_{m}^{\dagger}\hat{\psi}_{m}$,
which only depends on the number operators and conserved quantities of the
system as Barnett1985 ; Kumar1980

$\displaystyle\frac{\hbar^{2}}{2}\dot{\hat{N}}_{m}^{2}=\frac{\hbar^{2}}{8}\dot{\hat{N}}_{a}^{2}=-\frac{1}{2}[\hat{N}_{m},\hat{H}]^{2}=-\frac{1}{2}(\hat{H}-\varepsilon_{m}\hat{N}_{m})^{2}+g_{3}^{2}\hat{N}_{a}[(\hat{N}_{m}+1)(\hat{N}_{a}-1)(\hat{N}_{a}-2)^{2}+\hat{N}_{m}\hat{N}_{a}(\hat{N}_{a}+1)(\hat{N}_{a}+2)],$
(S6)

where $\hat{H}$ is the total energy. We introduce the many-body reaction
potential $\hat{V}_{3}$ through
$\hbar^{2}\dot{\hat{N}}_{m}^{2}/2+\hat{V}(\hat{N}_{m})=const.$ as shown in Eq.
(3) and obtain it under the boundary condition
$\langle\hat{V}_{3}\rangle|_{N_{m}=0,N_{a}=N}=0$ as:

$\displaystyle\hat{V}_{3}$
$\displaystyle=\frac{1}{2}[\hat{N}_{m},\hat{H}]^{2}-\frac{1}{2}[\hat{N}_{m},\hat{H}]^{2}|_{\hat{N}_{m}=0,\hat{N}_{a}=\hat{N}}$
$\displaystyle=-\varepsilon_{m}\hat{N}_{m}(\hat{H}-\frac{\varepsilon_{m}}{2}\hat{N}_{m})-g_{3}^{2}P_{3}(\hat{N}-2\hat{N}_{m},\hat{N}_{m})+g^{2}\hat{N}(\hat{N}-1)(\hat{N}-2)^{2},$
(S7)

where the three-body polynomial is
$P_{3}(x,y)=x[(y+1)(x-1)(x-2)^{2}+yx(x+1)(x+2)]$ and
$\hat{N}=\hat{N}_{a}+2\hat{N}_{m}$ is the conserved total population that
commutes with the Hamiltonian.

Figure S3: The reaction potential $V_{3}$ for the three-body process and the
oscillation frequency $\omega$ near the minimum of $V_{3}$ at different
detuning $\epsilon$. a, The reaction potential calculated based on Eq. (S9) at
$\epsilon=0$ (blue), $\epsilon=0.5$ (orange) and $\epsilon=1$ (yellow). The
region below the dashed line at $V_{3}=0$ is where the system is allowed to
reach based on the conservation law in Eq. (3). b, The oscillation frequency
$\omega$ at different detunings $\varepsilon_{m}$ (solid line). The dashed
line represents the asymptote $|\varepsilon_{m}|$ in the large detuning limit.
On resonance $\epsilon=0$, the frequency is $\omega=16g_{3}N^{3/2}/5$.

In our experiment, the atoms and molecules have macroscopic population in the
reaction process $N_{a},N_{m}\gg 1$. We may replace all the operators by their
expectation values and express the reaction potential in terms of the total
particle number $N$ and the molecular fraction $f_{m}=2N_{m}/N$ as

$\displaystyle V_{3}$
$\displaystyle=-g_{3}^{2}Nf_{m}\bigg{\\{}[(1-f_{m})^{4}-\frac{\epsilon^{2}}{8}f_{m}]N^{4}+(\frac{\epsilon\tilde{H}}{2}+2f_{m}^{3}-7f_{m}^{2}+9f_{m}-5)N^{3}$
$\displaystyle+5(2f_{m}^{2}-5f_{m}+4)N^{2}+2(5f_{m}-9)N+4\bigg{\\}},$ (S8)

where the dimensionless detuning and total energy are
$\epsilon=\varepsilon_{m}/g_{3}N^{3/2}$ and $\tilde{H}=H/g_{3}N^{3/2}$,
respectively. In the thermodynamic limit where the particle number $N\gg 1$,
the potential is reduced to:

$\displaystyle V_{3}=-g_{3}^{2}N^{5}f_{m}[(1-f_{m})^{4}-\epsilon^{2}f_{m}/8].$
(S9)

The potential curves at different detunings are shown in Fig. S3a. The
potential minimum $f_{m0}$ satisfies $\partial_{f_{m}}V_{3}|_{f_{m0}}=0$. On
resonance $\epsilon=0$, the molecule fraction at the minimum is

$\displaystyle f_{m0}=\frac{1}{5},$ (S10)

and the oscillation frequency $\omega_{0}$ of the system around the minimum
determined by the curvature of the potential $V_{3}$ at $f_{m0}$ is

$\displaystyle\omega_{0}=\frac{16}{5}g_{3}N^{3/2}.$ (S11)

At finite detuning, the minimum position $f_{m0}$ and the oscillation
frequency $\omega$ can be solved numerically from
$(1-f_{m0})^{3}(1-5f_{m0})-\epsilon^{2}f_{m0}/4=0$ and
$\omega=g_{3}N^{3/2}\sqrt{16(1-f_{m0})^{2}(2-5f_{m0})+\epsilon^{2}}$, see Fig
S3b. In the large detuning limit $\epsilon\gg 1$, the frequency approaches the
absolute value of the molecular energy $|\varepsilon_{m}|$.

For the nonuniform Thomas-Fermi density distribution of a BEC in a harmonic
trap, the coupling strength $g_{3}$ depends on the particle number and is
given by $g_{3}=g_{3}^{\prime}/\Omega^{3/2}$, where $g_{3}^{\prime}$ is a
coupling constant, the effective trap volume is
$\Omega=(14\pi/15^{2/5})(Na\bar{a}^{4})^{3/5}$ and the oscillator length
$\bar{a}$ is determined by the trap frequencies as
$\bar{a}=\sqrt{\hbar/m(\omega_{x}\omega_{y}\omega_{z})^{1/3}}$ Pethick2008 .
Thus on resonance the frequency $\omega_{0}\propto N^{3/5}$. Note that in our
experiment, we measure the dependence of the oscillation frequency
$\omega_{0}$ on the initial total atom number $N_{0}$ and we assume that the
atom number in the equilibration phase is proportional to $N_{0}$. By fitting
the data in Fig. 4b, we extract the three-body coupling constant
$g_{3}^{\prime}/h=5.5(5)\times 10^{-18}\mathrm{cm}^{9/2}/\mathrm{s}$.

Similarly, for the Hamiltonian in Eq. (1) for the two-body process
$\mathrm{A}+\mathrm{A}\longleftrightarrow\mathrm{A}_{2}$, the exact effective
potential $\hat{V}_{2}$ under the boundary condition
$\langle\hat{V}_{2}\rangle|_{N_{m}=0,N_{a}=N}=0$ is derived to be

$\displaystyle\hat{V}_{2}=\varepsilon_{m}\hat{N}_{m}(\frac{1}{2}\varepsilon_{m}\hat{N}_{m}-\hat{H})-2g_{2}^{2}\hat{N}_{m}[(\hat{N}-2\hat{N}_{m})^{2}-\hat{N}+2].$
(S12)

After replacing the operators by their expectation values we obtain

$\displaystyle V_{2}=-g_{2}^{2}Nf_{m}\bigg{\\{}[(1-f_{m})^{2}-\epsilon^{\prime
2}f_{m}/8]N^{2}+(\epsilon^{\prime}\tilde{H}^{\prime}/2-1)N+2\bigg{\\}},$ (S13)

where the dimensionless detuning and total energy are
$\epsilon^{\prime}=\varepsilon_{m}/g_{2}\sqrt{N}$ and
$\tilde{H}^{\prime}=H/g_{2}\sqrt{N}$ for the two-body process. In the
thermodynamic limit $N\gg 1$, the potential is reduced to:

$\displaystyle V_{2}=-g_{2}^{2}N^{3}f_{m}[(1-f_{m})^{2}-\epsilon^{\prime
2}f_{m}/8],$ (S14)

based on which we derive the analytical expressions for the minimum position
and the oscillation frequency near the minimum as:

$\displaystyle f_{m0}$ $\displaystyle=\frac{16+\epsilon^{\prime
2}-\sqrt{\epsilon^{\prime 4}+32\epsilon^{\prime 2}+64}}{24},$
$\displaystyle\omega$ $\displaystyle=g_{2}\sqrt{N}(\epsilon^{\prime
4}+32\epsilon^{\prime 2}+64)^{1/4}.$ (S15)

On resonance $\epsilon^{\prime}=0$, we have

$\displaystyle f_{m0}$ $\displaystyle=\frac{1}{3},$ $\displaystyle\omega_{0}$
$\displaystyle=2\sqrt{2}g_{2}\sqrt{N}.$ (S16)

For the nonuniform density distribution, the frequency is
$\omega_{0}\propto\sqrt{N/\Omega}=N^{1/5}$.

In general, for reactions of order $\alpha$ where $\alpha$ is the maximum
number of field operators in the terms of the reaction Hamiltonian, the
reaction potential is $V\propto N^{\alpha}$ in the thermodynamic limit for a
uniform system where each mode is macroscopically occupied. Then the
oscillation frequency of the system around the potential minimum $f_{m0}$ is
determined by the curvature of the potential as
$\omega_{0}\propto\sqrt{\partial^{2}V/\partial N^{2}|_{f_{m0}}}\propto
N^{\alpha/2-1}$.

##  IV Evaluation and theory of the molecule formation rate coefficient

The molecule formation rate coefficient $\beta$ presented in Fig. 2b is
evaluated based on $\beta=\dot{N}_{m}/N_{0}n_{0}$. Apart from the growth rate
of the molecule fraction $f_{m}=2N_{m}/N_{0}$ shown in Fig. 2a, we also need
to evaluate the initial mean atomic density $n_{0}$ to calculate $\beta$. We
model the atomic density distribution for the general case where both normal
and BEC components exist as ketterle1999

$\displaystyle
n(\mathbb{r},t)=n_{n}(\mathbb{r},t)+n_{\mathrm{B}}(\mathbb{r},t)=n_{np}g_{3/2}(\prod_{i=1}^{3}e^{-r_{i}^{2}/\sigma_{i}^{2}})+n_{\mathrm{B}p}(1-\sum_{i=1}^{3}\frac{r_{i}^{2}}{R_{i}^{2}})\theta(1-\sum_{i=1}^{3}\frac{r_{i}^{2}}{R_{i}^{2}}),$
(S17)

where $(r_{1},r_{2},r_{3})=(x,y,z)$ are spatial coordinates, $\sigma_{i}$ and
$R_{i}$ are the widths of the normal and BEC components,
$g_{3/2}(x)=\sum_{k=1}^{\infty}x^{k}/k^{3/2}$ is the polylogarithm function
and $\theta(x)$ is the Heaviside step function. Our vertical absorption
imaging integrates over the z axis and we obtain the 2D column density as

$\displaystyle
n(x,y,t)=\int_{-\infty}^{\infty}n(\mathbb{r},t)dz=\sqrt{\pi}n_{np}\sigma_{z}g_{2}(\prod_{i=1}^{2}e^{-r_{i}^{2}/\sigma_{i}^{2}})+\frac{4}{3}n_{\mathrm{B}p}R_{z}(1-\frac{x^{2}}{R_{x}^{2}}-\frac{y^{2}}{R_{y}^{2}})^{3/2}\theta(1-\frac{x^{2}}{R_{x}^{2}}-\frac{y^{2}}{R_{y}^{2}}).$
(S18)

We use Eq. (S18) to fit 2D column density distributions from our experiments
and extract the widths $\sigma_{x}$, $\sigma_{y}$, $R_{x}$ and $R_{y}$ and
total particle number $N_{n}$ and $N_{\mathrm{B}}$ in the normal and BEC
components. We use the scaling to calculate the cloud widths in the z
direction as $\sigma_{z}=(\omega_{y}/\omega_{z})\sigma_{y}$ and
$R_{z}=(\omega_{y}/\omega_{z})R_{y}$ based on our independent measurements of
trap frequencies $\omega_{y}$ and $\omega_{z}$. Then the 3D peak densities are
evaluated as

$\displaystyle
n_{np}=\frac{N_{n}}{g_{3}(1)\pi^{3/2}\bar{\sigma}^{3}},~{}n_{\mathrm{B}p}=\frac{15N_{\mathrm{B}}}{8\pi\bar{R}^{3}},$
(S19)

where $\bar{\sigma}=(\sigma_{x}\sigma_{y}\sigma_{z})^{1/3}$ and
$\bar{R}=(R_{x}R_{y}R_{z})^{1/3}$. Eventually, with all the parameters in Eq.
(S17) determined, we evaluate the initial mean 3D density as

$\displaystyle n_{0}=\frac{\int n^{2}(\mathbb{r},0)d^{3}\mathbb{r}}{\int
n(\mathbb{r},0)d^{3}\mathbb{r}},$ (S20)

as shown in Fig. S4a.

To capture the sharp transition of the rate coefficient around the critical
temperature $T_{c}$ as shown in Fig. 2b, we develop the following two-
component model incorporating both the thermal and condensed atoms. There are
three possible collision processes leading to molecule formation, which are
represented by the following rate equations involving molecular density
$m(\mathbb{r},t)$ and atomic density $n(\mathbb{r},t)$ with the subscripts $n$
and $\mathrm{B}$ indicating normal and BEC components, respectively.

1\. Molecule formation from collisions between condensed atoms:

$\dot{m}_{\mathrm{B}}(\mathbb{r},t)=\beta_{\mathrm{B}}n^{2}_{\mathrm{B}}(\mathbb{r},t)$
(S21)

with the rate coefficient given by
$\beta_{\mathrm{B}}=b_{\mathrm{Q}}\Gamma_{\mathrm{Q}}$, where the collision
rate is
$\Gamma_{\mathrm{Q}}=\frac{4\sqrt{3}}{(6\pi^{5})^{1/3}}\frac{h}{m}\frac{\langle
n^{2/3}\rangle}{\langle n\rangle}$, $\langle n\rangle=\int
n^{2}(\textbf{r},t)d^{3}\textbf{r}/\int n(\textbf{r},t)d^{3}\textbf{r}$ and
$\langle n^{2/3}\rangle=\int n^{5/3}(\textbf{r},t)d^{3}\textbf{r}/\int
n(\textbf{r},t)d^{3}\textbf{r}$.

2\. Molecule formation from collisions between thermal atoms:

$\dot{m}_{n}(\mathbb{r},t)=\beta_{n}n^{2}_{n}(\mathbb{r},t),$ (S22)

with the rate coefficient given by
$\beta_{n}=b_{\mathrm{Cl}}\Gamma_{\mathrm{Cl}}$, where the collision rate is
$\Gamma_{\mathrm{Cl}}=8h^{2}(\pi^{3}m^{3}k_{\mathrm{B}}T)^{-1/2}$.

3\. Molecule formation from collisions between one BEC atom and one thermal
atom:

$\dot{m}_{n\mathrm{B}}(\mathbb{r},t)=2\beta_{n\mathrm{B}}n_{\mathrm{B}}(\mathbb{r},t)n_{n}(\mathbb{r},t)=2\beta_{n}n_{\mathrm{B}}(\mathbb{r},t)n_{n}(\mathbb{r},t),$
(S23)

where we assume
$\beta_{n\mathrm{B}}=\beta_{n}=b_{\mathrm{Cl}}\Gamma_{\mathrm{Cl}}$.

Figure S4: Measurement of the mean atomic density $n_{0}$ and fitting the
molecule formation rate coefficient by the two-component model. a, Initial 3D
mean atomic density as a function of temperature, evaluated based on Eq.
(S20). The solid line is a guide to the eye. b, Molecule formation rate
coefficient fit by Eq. (S25) with the branching ratio $b_{\mathrm{Q}}$ as a
fitting parameter (green solid line).

The total molecular formation rate is given by the summation of Eqs.
(S21)-(S23):

$\displaystyle\dot{m}(\textbf{r},t)=\beta_{\mathrm{B}}n_{\mathrm{B}}^{2}(\textbf{r},t)+\beta_{n}[n_{n}(\textbf{r},t)+2n_{\mathrm{B}}(\textbf{r},t)]n_{n}(\textbf{r},t)=\beta
n(\textbf{r},t).$ (S24)

After the spatial integration, we obtain the molecular formation rate
coefficient in the mixture

$\displaystyle\beta=w_{\mathrm{B}}\beta_{\mathrm{B}}+(1-w_{\mathrm{B}})\beta_{n},$
(S25)

where the weight of the BEC component is

$\displaystyle w_{\mathrm{B}}=\frac{\int
n^{2}_{\mathrm{B}}(\mathbb{r},t)d^{3}\mathbb{r}}{\int
n^{2}(\mathbb{r},t)d^{3}\mathbb{r}}.$ (S26)

For a uniform system, the weight is reduced to
$w_{\mathrm{B}}=f^{2}_{\mathrm{BEC}}$ where
$f_{\mathrm{BEC}}=N_{\mathrm{B}}/N$ is the BEC fraction. Since we perform the
rate coefficient measurements in a harmonic trap, the density distributions
are nonuniform and the weight is enhanced to be
$w_{\mathrm{B}}>f^{2}_{\mathrm{BEC}}$. Our two-component model captures the
transition of the measured rate coefficients around $T_{c}$ very well, see
Fig. S4b.

##  V Extraction of molecule oscillation frequency and atom loss rate

We use the following function to fit the data in the equilibration phase at
$1~{}\mathrm{ms}<t<3~{}\mathrm{ms}$ shown in Fig. 3 for the extraction of
molecule oscillation frequencies Claussen2003 :

$\displaystyle N_{m}(t)=N_{m}(0)-\gamma_{1}t+\Delta
N_{m}e^{-\gamma_{2}t}\sin(\omega t+\phi),$ (S27)

where the fitting parameters are the molecule number $N_{m}(0)$ extrapolated
to time $t=0$, decay rates $\gamma_{1}$ and $\gamma_{2}$, oscillation
amplitude $\Delta N_{m}$, oscillation frequency $\omega$ and phase $\phi$.
Here the two decay rates $\gamma_{1}$ and $\gamma_{2}$ characterize the decay
of molecule number and the damping of molecule oscillation amplitude, which
are generally different.

For the data shown in Fig. 4, we fit the data at
$0.3~{}\mathrm{ms}<t<3~{}\mathrm{ms}$ using the function

$\displaystyle N_{m}(t)=e^{-\gamma t}[N_{m}(0)+\Delta N_{m}\sin(\omega
t+\phi)],$ (S28)

where we find the single decay rate $\gamma$ is enough to describe the data
very well. For each fit, we subtract from the time $t$ a delay time of 0.15 ms
due to the finite speed of our magnetic field switch.

Figure S5: Examples of atomic density evolution in a 2D flat-bottomed optical
potential for the data presented in Fig. 3c. For data below the resonance,
BECs are initially prepared at 19.5 G and magnetic field is quenched to values
between 0.05 and 1 G (panel a) and between 5 and 50 mG (panel c) below the
resonance. Relaxation and equilibration phases are marked with different
background colors in panel c. For data above the resonance, BECs are initially
prepared at 20.4 G and magnetic field is quenched to values between 0.1 and 1
G (panel b) and between 10 and 50 mG (panel d) above the resonance. Solid
lines are fits for extracting the atom loss rates, see text.

To prevent the fits from getting stuck in a local optimum, we vary the initial
guess of the frequency $\omega$ for the fits and use the result that has the
minimum root mean square error.

For the atom loss rate measurement shown in Fig. 3c, we present example time
traces of the averaged atomic density in the 2D flat-bottomed trap in Fig. S5.
Far from the resonance, see Fig. S5a and b, the atomic density decays slowly
and we fit the data using

$\displaystyle n_{a}(t)=n_{a}(0)e^{-\gamma_{a}t},$ (S29)

where $n_{a}(0)$ is the initial atomic density and $\gamma_{a}$ is the atom
loss rate. The fit is applied to the data above half of the initial density.

Below and near the resonance, see Fig. S5c, we find the density first decays
rapidly and then settles around some equilibrium value before a slow decay
kicks in at a longer time scale than 3 ms. In this case, we use the following
fit function:

$\displaystyle
n_{a}(t)=n_{a}(0)\bigg{\\{}\theta(t_{0}-t)+[(1-s)e^{-\gamma_{a}(t-t_{0})}+s]\theta(t-t_{0})\bigg{\\}},$
(S30)

where $t_{0}$ is the time when the decay begins and $s$ represents the
fractional density the system settles to after the initial fast decay. On the
other hand, above and near the resonance, see Fig. S5d, the data is fit well
by a single exponential decay as in Eq. (S29).

Figure S6: Amplification of coherent molecule oscillations through magnetic
field modulation. a, Schematic diagram of molecule formation near the Feshbach
resonance with additional sinusoidal magnetic field modulation, where the
black wavy arrow represents the RF photon from the modulation. b, Evolution of
the molecular population after quench to a static magnetic field 2(1) mG below
the resonance, where the contrast $C$ of the molecule oscillation is $9(3)\%$.
c, Time traces of molecule number with magnetic field modulation at 0.85
(red), 1.71 (green) and 3.42 kHz (magenta) with modulation amplitude
$B_{\mathrm{ac}}\approx 4$ mG around the same field as in panel b. The solid
lines are fits based on Eq. (S27). The contrast is defined as
$C=N_{m}(0)/\Delta N_{m}$, see text. d, Frequency of the molecule oscillations
extracted from the fits in panel c. Blue solid line is a linear fit without an
offset, which gives a slope 0.97(3). The purple dashed line represents the
molecule oscillation frequency in panel b.

##  VI Amplification of the molecule oscillation through magnetic field
modulation

We observe large amplitude molecule number oscillations by applying additional
external driving to atomic BECs near the Feshbach resonance. The coupled
atomic and molecular BECs in our system effectively form a bosonic Josephson
junction Liu2010 , in close analogy to Cooper-pair Josephson junctions in
superconducting devices barone1982 . Inspired by the Shapiro resonance effect
where adding a small resonant ac component to an applied voltage enhances dc
tunneling current in the superconducting Josephson junction, we modulate the
magnetic field at a frequency close to the free molecule oscillation frequency
at a static field, with the hope to facilitate the reaction. At 2 mG below the
resonance, the molecule number oscillates around 4,000 at 1.4 kHz, with a
contrast of $\sim$9% (see Fig. S6b). After adding a magnetic field modulation
with an amplitude of 4 mG around the same field, we find the oscillation
contrast increases by a factor of 3 to 4, see Fig. S6c. We also see clear
damping in the oscillations with the driving field, likely due to additional
heating introduced by the driving.