# Achieving Multi-beam Gain in Intelligent Reflecting Surface Assisted
Wireless Energy Transfer

Chi Qiu, Qingqing Wu, Meng Hua, Xinrong Guan, Yuan Wu C. Qiu, Q. Wu, M. Hua
and Y. Wu are with the State Key Laboratory of Internet of Things for Smart
City, University of Macau, Macau, 999078, China (email:
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>yuanwu@um.edu.mo). X. Guan is with the College of Communications Engineering,
Army Engineering University of PLA, Nanjing 210007, China (e-mail:
guanxr@aliyun.com).

###### Abstract

Intelligent reflecting surface (IRS) is a promising technology to boost the
efficiency of wireless energy transfer (WET) systems. However, for a multiuser
WET system, simultaneous multi-beam energy transmission is generally required
to achieve the maximum performance, which may not be implemented by using the
IRS having only a single set of coefficients. As a result, it remains unknowns
how to exploit the IRS to approach such a performance upper bound. To answer
this question, we aim to maximize the total harvested energy of a multiuser
WET system subject to the user fairness constraints and the non-linear energy
harvesting model. We first consider the static IRS beamforming scheme, which
shows that the optimal IRS reflection matrix obtained by applying semidefinite
relaxation is indeed of high rank in general as the number of energy receivers
(ERs) increases, due to which the resulting rank-one solution by applying
Gaussian Randomization may lead to significant loss. To achieve the multi-beam
gain, we then propose a general time-division based novel framework by
exploiting the IRS’s dynamic passive beamforming. Moreover, it is able to
achieve a good balance between the system performance and complexity by
controlling the number of IRS shift patterns. Finally, we also propose a time-
division multiple access (TDMA) based passive beamforming design for
performance comparison. Simulation results demonstrate the necessity of multi-
beam transmission and the superiority of the proposed dynamic IRS beamforming
scheme over existing schemes.

###### Index Terms:

IRS, WET, dynamic beamforming, resource allocation.

##  I Introduction

Intelligent reflecting surface (IRS) has been proposed as a promising solution
to reshape the wireless propagation environment by tuning the reflecting
elements [1]. This thus motivates recent research on IRS-aided wireless energy
transfer (WET) [2] and other related applications [3, 4, 5, 6]. In particular,
without the need of any complicated signal processing units, IRS enables
intelligent signal reflection over a large aperture to compensate the severe
wireless channel attenuation over long distance for improving WET efficiency
[7].

Fig. 1: An IRS-aided WET system.

Fig. 2: Transmission protocol of the proposed dynamic IRS beamforming design.

On the other hand, in multiuser WET system, how to fairly enhance energy
transmission efficiency is also an important issue, which has been studied in
prior works [8, 9]. It was shown in [8] that the multi-beam energy
transmission is the optimal design for the multi-antenna energy transmitter
(ET) to maximize the minimum energy harvested by all the energy receivers
(ERs). It also reveals that transmitting only one energy beam with adjustable
weights over time achieves the same optimal performance as conventional multi-
beam design. A similar phenomenon has also been revealed in [9] under non-
linear energy harvesting (EH) model. Recall that IRS is only implemented with
a single set of coefficients, which, however, is not capable of providing the
optimal so-called multi-beam transmission in WET. The above results motivate
us to design a series of IRS single beam patterns within the channel coherence
time to approach the maximum performance of simultaneous multi-beam
transmission, and thus to fully unleash the potential of IRS.

To this end, it still remains an open question that whether there exists
performance gain from such time-sharing of IRS reflection patterns and thus
needed to be exploited. In this letter, we consider a multiuser IRS-aided WET
system. To avoid severe resource allocation mismatches and significant
performance degradation for practical systems, the non-linear EH model is
adopted. We first study a special case, where only one reflection pattern is
allowed during the whole WET. It is found that the rank of the optimal relaxed
IRS reflection coefficient matrix by applying semidefinite relaxation (SDR)
indeed grows with the number of ERs in the system, which suggests that a
substantial performance will be incurred if a rank-one solution is obtained.
To mimic the simultaneous multi-beam transmission, we propose a general
framework for dynamic IRS beamforming, where the IRS is allowed to adjust its
phase-shift vectors for an arbitrary number of times based on the rank of the
reflection matrix from the above to maximize the total harvested energy.
Particularly, the advantage of the conventional multiple beams are harnessed
via a series of time-varying IRS’s reflection patterns. For computational
complexity reduction and performance comparison, we then propose a time-
division multiple access (TDMA) based beamforming scheme. Specifically, each
phase-shift vector is exclusively designed to maximize each ER’s harvested
energy in each time slot. To solve the formulated direct-current (DC) energy
maximization problems for the above cases, we propose efficient algorithms by
exploiting successive convex approximation (SCA) technique. Numerical results
demonstrate the superiority of the proposed dynamic IRS beamforming scheme,
and reveal the upside of time-sharing principle.

##  II System Model

As shown in Fig. 1, we consider an IRS assisted WET system, where an IRS with
$N$ passive reflecting elements is deployed to assist the WET between a
single-antenna ET and a set of single-antenna ERs, denoted by the set
$\mathcal{K}=\\{1,\dots,K\\}$. A quasi-static flat fading channel is
considered, and the total available transmission time is given by $T$. We
assume that the CSI is perfectly known at the ET by using the state-of-art
channel estimation approaches to achieve the system performance upper bound
[3]. Denote $h_{d,k}\in\mathbb{C}^{1\times 1}$,
$\mathbf{g}\in\mathbb{C}^{N\times 1}$, and
$\mathbf{h}_{r,k}^{H}\in\mathbb{C}^{1\times N}$, $\forall k\in\mathcal{K}$, as
the baseband equivalent channels from the ET to ER $k$, from the ET to the
IRS, and from the IRS to ER $k$, respectively. As shown in Fig. 2, a novel
transmission protocol of the dynamic IRS beamforming design is proposed. To be
specific, each channel coherence time $T$ is divided into $J$ time slots, each
with duration of $\tau_{j}$, $j\in{\cal J}=\\{1,\dots,J\\}$. In each time
slot, we assume that the IRS adjusts its reflection pattern individually and
ET adjusts its transmit power. To maximize the reflected signal power by the
IRS, the reflection amplitude is set to be one. Hence, the reflection
coefficient matrix of the IRS in the $j$th time slot is given by
$\mathbf{\Phi}_{j}=\operatorname{diag}\\{e^{j\theta_{j,1}},\dots,e^{j\theta_{j,N}}\\}\in\mathbb{C}^{N\times
N}$, where $\theta_{j,n}\in[0,2\pi),n\in\mathcal{N}=\\{1,\dots,N\\}$ denotes
the phase shift of the $n$th IRS reflecting element.

The RF-power received at ER $k$ in time slot $j$ is given by

$P^{r}_{k,j}=P^{t}_{j}|\mathbf{h}_{r,k}^{H}\mathbf{\Phi}_{j}\mathbf{g}+h_{d,k}^{H}|^{2}=P^{t}_{j}|\mathbf{q}_{k}^{H}\bm{\theta}_{j}+h_{d,k}^{H}|^{2},$
(1)

where $P^{t}_{j}$ represents the allocated transmit power at ET in time slot
$j$, $\mathbf{q}_{k}^{H}=\mathbf{h}_{r,k}^{H}\operatorname{diag}(\mathbf{g})$
and $\bm{\theta}_{j}=[e^{j\theta_{j,1}},\dots,e^{j\theta_{j,N}}]^{T}$ denotes
the IRS passive beamforming vector. Here, we adopt a practical non-linear EH
model to convert the RF power to DC power [10]. As such, the harvested DC
power by ER $k$ in the $j$th time slot is given by

$\Phi_{{k,j}}=\frac{X_{k}}{1+e^{-a_{k}(P^{r}_{k,j}-b_{k})}}-Y_{k},\vspace{-0.2cm}$
(2)

where $X_{k}=M_{k}\left(1+e^{a_{k}b_{k}}\right)/e^{a_{k}b_{k}}$,
$Y_{k}=\frac{M_{k}}{e^{a_{k}b_{k}}}$, $a_{k}$, $b_{k}$ and $M_{k}$, represent
the ER $k$’s circuit configurations.

The objective of this letter is to maximize the transferred energy to all ERs
subject to energy fairness constraints by jointly optimizing the IRS passive
beamforming and resource allocation. We first study _static IRS beamforming_
scheme, then propose the _dynamic IRS beamforming_ scheme, and finally study
the _TDMA-based IRS beamforming_ scheme, which are discussed in details as
follows.

###  II-A WET with Static IRS Beamforming

((a)) Average rank of $\bm{\Theta}$ versus $K$ with $N=100$.

((b)) The 10 largest average eigenvalues of $\bm{\Theta}$ with $K=60$.

Fig. 3: Illustration of properties of $\bm{\Theta}$.

In this case, the phase-shift pattern of the IRS remains unchanged throughout
WET. Let $\bm{\theta}$ denote the phase-shift vector for static IRS
beamforming scheme. Thus, the transmit power will also remain constant,
denoted by $P$. Thus, the harvested DC power by ER $k$ is
$\Phi_{{k}}=\frac{X_{k}}{1+e^{-a_{k}(P|\mathbf{q}_{k}^{H}\bm{\theta}+h_{d,k}^{H}|^{2}-b_{k})}}-Y_{k}$.
The corresponding optimization problem is formulated as

$\displaystyle\max_{P,\bm{\theta},e}$ $\displaystyle\ e$ (3a) s.t.
$\displaystyle T\Phi_{k}\geq\alpha_{k}e,\forall k\in\mathcal{K},$ (3b)
$\displaystyle PT\leq E_{\operatorname{tot}},0\leq P\leq
P_{\operatorname{max}},$ (3c) $\displaystyle|[\bm{\theta}]_{n}|=1,\forall
n\in\mathcal{N}.$ (3d)

Constraint (3b) ensures the fairness among all ERs, with $\alpha_{k}\geq 0$,
denoting the target ratio of the $k$th ER’s harvested energy to total
harvested energy with $\sum_{k=1}^{K}\alpha_{k}=1$. (3c) denotes the total
energy and maximum power constraint, with $E_{\operatorname{tot}}$ denoting
the energy budget and $P_{\operatorname{max}}$ representing the maximum
allowed transmit power. (3d) stands for the unit-modulus constraint of the IRS
phase-shift.

Note that the above problem is non-convex, we propose an SDR-based algorithm
to solve it. Specifically, to fully utilize the available energy at the ET, it
can be easily verified that the optimal transmit power is
$P=E_{\operatorname{tot}}/{T}$, and thus constraint (3c) can be dropped. Let
$|\mathbf{q}_{k}^{H}\bm{\theta}_{j}+h_{d,k}^{H}|=|\bar{\mathbf{q}}_{k}^{H}\bar{\bm{\theta}}_{j}|$,
where $\bar{\bm{\theta}}_{j}=[\bm{\theta}_{j}^{H}\ 1]^{H}$ and
$\bar{\mathbf{q}}_{k}^{H}=[\mathbf{q}_{k}^{H}\ h_{d,k}^{H}]$. Then, we can
rewrite $P^{r}_{k,j}$ in (1) as
$P^{t}_{j}\bar{\bm{\theta}}_{j}^{H}\mathbf{Q}_{k}\bar{\bm{\theta}}_{j}$.
Define
$\mathbf{Q}_{k}=\bar{\mathbf{q}}_{k}\bar{\mathbf{q}}^{H}_{k}=\left[\begin{array}[]{ll}\mathbf{q}_{k}\mathbf{q}^{H}_{k}&\mathbf{q}_{k}h_{d,k}^{H}\\\
h_{d,k}\mathbf{q}_{k}^{H}&h_{d,k}h_{d,k}^{H}\end{array}\right]$ and
$\bm{\Theta}=\bar{\bm{\theta}}\bar{\bm{\theta}}^{H}$, which needs to satisfy
$\bm{\Theta}\succeq 0$ and $\operatorname{rank}(\bm{\Theta})=1$. Since
constraint $\operatorname{rank}(\bm{\Theta})=1$ is non-convex, we can apply
SDR technique to relax it. As a result, problem (3) is reduced as

$\displaystyle\max_{\bm{\Theta},e}$ $\displaystyle\ e$ (4a) s.t.
$\displaystyle P\operatorname{tr}(\bm{Q}_{k}\mathbf{\Theta})\geq
b_{k}-\frac{1}{a_{k}}\operatorname{ln}\left(\frac{X_{k}}{\frac{\alpha_{k}e}{T}+Y_{k}}\right),\forall
k\in\mathcal{K},$ (4b) $\displaystyle\bm{\Theta}\succeq 0,$ (4c)
$\displaystyle[\bm{\Theta}]_{n,n}=1,\forall n\in\mathcal{N}.$ (4d)

Although problem (4) is still non-convex, we observe that with fixed $e$, it
is reduced to a linear feasibility check problem. As such, the global optimal
solution can be obtained via a bisection/one-dimensional search over $e$.
Thus, the solution to problem (4) serves as an upper bound to problem (3).
Since the obtained $\bm{\Theta}$ is not guaranteed to be a rank-one matrix,
the Gaussian Randomization (GR) approach can be applied to obtain the sub-
optimal solution. As aforementioned, the rank of the IRS reflection
coefficient matrix is of interest. For illustration, Fig. 3(a) shows the
average rank of $\bm{\Theta}$ versus the number of ERs and Fig. 3(b) shows the
distribution of the 10 largest average eigenvalues with $K=60$ over $100$
channel realizations, both under the same system setup as will be specified
later in Section V. The rank of $\bm{\Theta}$ is calculated based on the
number of eigenvalues exceeds $2\%$ of the maximum one, from which we can see
that as $K$ increases, the rank of $\bm{\Theta}$ also increases, which will
result in performance degradation when a rank-one solution is reconstructed by
applying GR approach. Motivated by this, we propose _dynamic IRS beamforming_
and _TDMA-based IRS beamforming_ schemes, elaborated in the next section.

##  III Dynamic Beamforming Based Multi-beam Design

In this section, we propose a general time-division based novel framework by
exploiting the IRS’s dynamic passive beamforming. Then, for performance
comparison, we also propose a TDMA-based passive beamforming design.

###  III-A Dynamic IRS Beamforming

In this case, the IRS reflection pattern is allowed to adjust for $J$ times,
where $J=\operatorname{rank}(\bm{\Theta})$ from Section II-A. As such, the IRS
generates multi-beams in a time division manner. As introduced early, in this
design, we jointly optimize the IRS passive beamforming, time allocation and
power allocation. Accordingly, the problem can be formulated as follows

$\displaystyle\max_{\\{P^{t}_{j}\\},\\{\tau_{j}\\},\\{\bm{\theta}_{j}\\},e}e$
(5a) $\displaystyle\text{ s.t.
}~{}\sum_{j=1}^{J}\tau_{j}\Phi_{{k,j}}\geq\alpha_{k}e,\forall
k\in\mathcal{K},$ (5b) $\displaystyle\qquad\sum_{j=1}^{J}\tau_{j}\leq
T,\tau_{j}\geq 0,\forall j\in\mathcal{J},$ (5c)
$\displaystyle\qquad\sum_{j=1}^{J}P^{t}_{j}\tau_{j}\leq
E_{\operatorname{tot}},0\leq P^{t}_{j}\leq P_{\operatorname{max}},\forall
j\in\mathcal{J},$ (5d) $\displaystyle\qquad|[\bm{\theta}_{j}]_{n}|=1,\forall
n\in\mathcal{N},\forall j\in\mathcal{J}.$ (5e)

Constraint (5c) represents the total and non-zero time constraint.

###  III-B TDMA-Based IRS Beamforming

To reduce computational complexity, the TDMA-based IRS beamforming scheme is
proposed. In this design, $T$ is divided into $K$ time slots. In any time slot
$k$, the IRS designs its phase-shift pattern to maximize the harvested energy
at ER $k$, with
$\bm{\theta}_{k}^{*}=e^{j(arg(h_{d,k}^{H})-arg(\bm{q}_{k}^{H}))}$. In this
case, only time and power allocations are left to be optimized. Then the
harvested DC power by $k$-th ER is give by
$\Phi_{{k,j}}^{*}=\frac{X_{k}}{1+e^{-a_{k}(P^{t}_{j}|\mathbf{q}_{k}^{H}\bm{\theta}_{j}^{*}\\!+\\!h_{d,k}^{H}|^{2}-b_{k})}}-Y_{k}$.
Thus, the problem is reduced to

$\displaystyle\max_{\\{P^{t}_{j}\\},\\{\tau_{j}\\},e}e$ (6a)
$\displaystyle\text{ s.t.
}~{}\sum_{j=1}^{K}\tau_{j}\Phi_{k,j}^{*}\geq\alpha_{k}e,\forall
k\in\mathcal{K},$ (6b) $\displaystyle\qquad\eqref{4c},\eqref{4d}.$ (6c)

The above problem formulations in (5) and (6) are both difficult to solve due
to the non-convex constraints and highly coupled variables. In the following
section, we propose an efficient algorithm, based on SCA techniques to solve
them suboptimally.

##  IV Proposed solutions

In this section, efficient algorithms to solve problems (5) and (6) are
proposed.

###  IV-A Dynamic IRS Beamforming

Recall that
$|\mathbf{q}_{k}^{H}\bm{\theta}+h_{d,k}^{H}|^{2}=\bar{\bm{\theta}}^{H}\mathbf{Q}_{k}\bar{\bm{\theta}}$,
and by introducing slack variables
$z=\\{z_{k,j},j\in\mathcal{J},k\in\mathcal{K}\\}$, constraint (5b) can be
equivalently transformed to

$\sum_{j=1}^{J}\tau_{j}\left(\frac{X_{k}}{1+z_{k,j}}-Y_{k}\right)\geq\alpha_{k}e,\forall
k\in\mathcal{K},\vspace{-0.2cm}$ (7)

with additional constraint as

$e^{-a_{k}\left(P^{t}_{j}\bar{\bm{\theta}}_{j}^{H}\mathbf{Q}_{k}\bar{\bm{\theta}}_{j}-b_{k}\right)}\leq
z_{k,j},\forall j\in\mathcal{J},\forall k\in\mathcal{K}.\vspace{-0.2cm}$ (8)

Note that the inequality constraint (8) is always satisfied with equality for
the optimal solution. This can be proved by contradiction. Suppose that the
optimal solution is obtained with inequality in (8), the value of $z_{k,j}$
can always be decreased to obtain a larger objective value, which is
contradictory to the assumption. It can be observed that constraints (7) and
(8) are still non-convex. However, the left-hand side (LHS) of (7) is jointly
convex with respect to (w.r.t). $\sqrt{\tau_{j}}$ and $z_{k,j}$. Recall that
any convex function is lower-bounded by its first-order Taylor expansion.
Thus, the SCA technique is applied. Therefore, for any local points
$\tau_{j}^{r}$ and $z_{k,j}^{r}$ in the $r$th iteration, we have

$\displaystyle\frac{\left(\sqrt{\tau_{j}}\right)^{2}}{1+z_{k,j}}$
$\displaystyle\geq\frac{\left(\sqrt{\tau_{j}^{r}}\right)^{2}}{1+z_{k,j}^{r}}-\frac{\left(\sqrt{\tau_{j}^{r}}\right)^{2}}{\left(1+z_{k,j}^{r}\right)^{2}}\left(z_{k,j}-z_{k,j}^{r}\right)$
$\displaystyle+\frac{2\sqrt{\tau_{j}^{r}}}{1+z_{k,j}^{r}}\left(\sqrt{\tau_{j}}-\sqrt{\tau_{j}^{r}}\right)\triangleq
f_{k,j}^{lb}\left(z_{k,j},\tau_{j}\right),$ (9)

which can be readily checked that $f_{k,j}^{lb}$ is jointly convex w.r.t.
$z_{k,j}$ and $\tau_{j}$. Then, to tackle the non-convexity of (8), we first
rewrite it into a more simplified form as

$P^{t}_{j}\bar{\bm{\theta}}_{j}^{H}\mathbf{Q}_{k}\bar{\bm{\theta}}_{j}\geq
b_{k}-\frac{\ln z_{k,j}}{a_{k}},\forall j\in\mathcal{J},\forall
k\in\mathcal{K}.\vspace{-0.2cm}$ (10)

Note that the optimization variables $P^{t}_{j}$ and $\bar{\bm{\theta}}_{j}$
are coupled in the LHS of (10). Let
$P^{t}_{j}\bar{\bm{\theta}}_{j}^{H}\mathbf{Q}_{k}\bar{\bm{\theta}}_{j}=\bar{\bm{\theta}}_{j}^{H}\mathbf{Q}_{k}\bar{\bm{\theta}}_{j}/\frac{1}{P^{t}_{j}}$.
It is observed that although it is not jointly convex w.r.t.
$\bar{\bm{\theta}}_{j}$ and $P^{t}_{j}$, it is jointly convex w.r.t.
$\bar{\bm{\theta}}_{j}$ and $\frac{1}{P^{t}_{j}}$, which motivates us to apply
SCA to tackle it. Specifically, with any feasible points
$\bar{\bm{\theta}}^{r}_{j}$ and $P^{t,r}_{j}$, the lower bound of
$\bar{\bm{\theta}}_{j}^{H}\mathbf{Q}_{k}\bar{\bm{\theta}}_{j}/\frac{1}{P^{t}_{j}}$
can be obtained as

$\displaystyle\frac{\bar{\bm{\theta}}_{j}^{H}\mathbf{Q}_{k}\bar{\bm{\theta}}_{j}}{\frac{1}{P^{t}_{j}}}\geq
2P^{t,r}_{j}\operatorname{Re}\\{\bar{\bm{\theta}}^{H}_{j}\mathbf{Q}_{k}\bar{\bm{\theta}}^{r}_{j}\\}-P^{t,r}_{j}\bar{\bm{\theta}}^{rH}_{j}\mathbf{Q}_{k}\bar{\bm{\theta}}^{r}_{j}$
$\displaystyle-\bar{\bm{\theta}}^{rH}_{j}\mathbf{Q}_{k}\bar{\bm{\theta}}^{r}_{j}/\frac{1}{P^{t,r2}_{j}}\left(\frac{1}{P^{t}_{j}}-\frac{1}{P^{t,r}_{j}}\right)\triangleq
g_{k,j}^{lb}\left(\bar{\bm{\theta}}_{j},P^{t}_{j}\right),$ (11)

which is jointly convex w.r.t. $\bar{\bm{\theta}}_{j}$ and $P^{t}_{j}$.

Note that constraint (5d) is non-convex due to its coupled variables
$P^{t}_{j}$ and $\tau_{j}$ in its LHS. We first rewrite
$\sum_{j=1}^{J}P^{t}_{j}\tau_{j}$ as
$\frac{1}{4}\sum_{j=1}^{J}((P^{t}_{j}+\tau_{j})^{2}-(P^{t}_{j}-\tau_{j})^{2})$.
Since the term $-(P^{t}_{j}-\tau_{j})^{2}$ is jointly concave w.r.t.
$P^{t}_{j}$ and $\tau_{j}$, its upper bound can be obtained by applying the
SCA technique with given initial points $P^{t,r}_{j}$ and $\tau^{r}_{j}$,
yielding

$\displaystyle-(P^{t}_{j}-\tau_{j})^{2}$
$\displaystyle\leq-(P^{t,r}_{j}-\tau^{r}_{j})^{2}-2(P^{t,r}_{j}-\tau^{r}_{j})(P^{t}_{j}-P^{t,r}_{j})$
$\displaystyle+2(P^{t,r}_{j}-\tau^{r}_{j})(\tau_{j}-\tau^{r}_{j})\triangleq\eta^{ub}_{j}\left(P^{t}_{j},\tau_{j}\right).$
(12)

Finally, for the non-convex unit-modulus constraint (5e), we relax it into a
convex form, given by

$|[\bm{\theta}_{j}]_{n}|\leq 1,\forall n\in\mathcal{N},\forall
j\in\mathcal{J}.$ (13)

As a result, given points $z_{k,j}^{r}$, $\tau_{j}^{r}$,
$\bar{\bm{\theta}}^{r}_{j}$ and $P^{t,r}_{j}$, problem (5) can be approximated
as

$\displaystyle\max_{\\{P^{t}_{j}\\},\\{\tau_{j}\\},\\{\bm{\theta}_{j}\\},e}e$
(14a) $\displaystyle\text{ s.t.
}~{}\sum_{j=1}^{J}\left(X_{k}f_{k,j}^{lb}\left(z_{k,j},\tau_{j}\right)-Y_{k}\tau_{j}\right)\geq\alpha_{k}e,\forall
k\in\mathcal{K},$ (14b) $\displaystyle\qquad
g_{k,j}^{lb}\left(\bar{\bm{\theta}}_{j},P^{t}_{j}\right)\geq b_{k}-\frac{\ln
z_{k,j}}{a_{k}},\forall j\in\mathcal{J},\forall k\in\mathcal{K},$ (14c)
$\displaystyle\qquad\frac{1}{4}\sum_{j}^{J}((P^{t}_{j}+\tau_{j})^{2}+\eta^{ub}_{j}\left(P^{t}_{j},\tau_{j}\right))\leq
E_{\operatorname{tot}},$ (14d) $\displaystyle\qquad 0\leq P^{t}_{j}\leq
P_{\operatorname{max}},\forall j\in\mathcal{J},$ (14e)
$\displaystyle\qquad\eqref{4c},\eqref{12},$ (14f)

which is convex and can be optimally solved by the interior-point method [11].
In addition, since the obtained phase shifts may not satisfy the unit-modulus
constraint, the feasible IRS phase-shift vectors are reconstructed by

$[\bm{\theta}_{j}]_{n}^{*}=\frac{[\bm{\theta}_{j}]_{n}}{|[\bm{\theta}_{j}]_{n}|},\forall
n\in\mathcal{N},\forall j\in\mathcal{J}.$ (15)

Finally, time and power allocations are updated by solving (14) based on the
reconstructed phase shifts. It follows that the obtained objective value
serves as a lower bound of that of problem (5).

###  IV-B TDMA-Based IRS Beamforming

Similar to problem (14), with the given points $z_{k,j}^{r}$, $\tau_{j}^{r}$
and $P^{t,r}_{j}$ in the $r$th iteration, problem (6) can be approximated as

$\displaystyle\max_{\\{P^{t}_{j}\\},\\{\tau_{j}\\},e}e$ (16a)
$\displaystyle\text{ s.t.
}P^{t}_{j}\bar{\bm{\theta}}_{j}^{*H}\mathbf{Q}_{k}\bar{\bm{\theta}}_{j}^{*}\geq
b_{k}-\frac{\ln z_{k,j}}{a_{k}},\forall j\in\mathcal{J},\forall
k\in\mathcal{K}.$ (16b) $\displaystyle\qquad 0\leq P^{t}_{j}\leq
P_{\operatorname{max}},\forall j\in\mathcal{J},$ (16c)
$\displaystyle\qquad\eqref{4c},\eqref{14b},\eqref{14d},$ (16d)

which can be solved by using the interior-point method [11]. Similarly, the
optimal objective value of problem (16) serves as a lower bound of that of
problem (6).

###  IV-C Complexity Analysis

The computational complexity of the two schemes is analyzed as follows. For
fairness comparison, the same $J$ is assumed. Specifically, the dynamic IRS
beamforming scheme needs to optimize more optimization variables and its
complexity is given by $\mathcal{O}\left((2J+NJ)^{3.5}\right)$ based on the
analytical results in [12], whereas the low complexity TDMA-based beamforming
scheme mainly needs to solve for problem (5) with only resource allocation,
which results in complexity of $\mathcal{O}\left((2J)^{3.5}\right)$.

##  V Numerical Results

In this section, we provide numerical results to validate the effectiveness of
our proposed scheme. In the simulation, a three-dimensional coordinate setup
is considered, where the ET is located at $(0,0,0)$ meter (m) and the IRS is
located at $(30,0,5)$ m. The ERs are randomly distributed within a circle
centered at $(30,0,0)$ m with radius of $5$ m. Rician fading channel is
considered for the ET-IRS and IRS-ER links with a Rician factor of $3$ dB,
while Rayleigh fading channel is considered for the ET-ER link. The large-
scale path loss model is given by $PL(d)=c_{0}\left(d/d_{0}\right)^{-\alpha}$,
where $c_{0}=-30$ dB is the path loss at a reference distance of $d_{0}=1$ m,
and path loss exponent $\alpha$ being $2.2$ for the ET-IRS and IRS-ER
channels, while $3.6$ for ET-ER channel. We assume all ERs have the same
circuit configurations with non-linear EH model parameters given by
$a_{k}=150$, $b_{k}=0.014$ and $M_{k}=0.024$ [13]. In addition, the ET antenna
gain and ER antenna gain are given by $10$ dBi and $3$ dBi, respectively.
Other parameters are set as $N=100$, $\alpha_{k}=1/K,\forall k,$
$E_{\operatorname{tot}}=10$ Joule (J), $P_{\operatorname{max}}=46$ dBm, $T=1$
s, $d_{0}=1$ m and $\epsilon=10^{-3}$. We consider the following schemes for
performance evaluation: 1) Upper bound: problem (3) solved by SDR-based
algorithm in Section II-A without applying GR; 2) Static IRS, GR: GR is
applied to recover the rank-one $\bm{\theta}$ based on the solution of SDR
[3]; 3) Static IRS, SCA: problem (3) solved by the SCA-based algorithm [6]; 4)
Dynamic IRS: problem (5) solved by the SCA-based algorithm proposed in Section
IV-A; 5) TDMA-based IRS beamforming: problem (6) solved in Section IV-B; 6)
Without IRS.

Fig. 4 shows the total harvested energy of all ERs versus the number of ERs.
It is observed that the IRS-aided WET system outperforms that without IRS,
which demonstrates the benefits of integrating IRS into WET system.
Particularly, the proposed dynamic IRS beamforming scheme achieves the best
performance among all schemes. This is excepted since the IRS is able to
achieve the multi-beam gain by proactively generating a series of time-varying
beam patterns. It is worth pointing out that the total harvested energy grows
as $K$ increases for all schemes, except for the GR approach. This is because
that the performance of GR approach significantly degrades from the high rank
of matrix $\bm{\Theta}$ as $K$ increases, which is consistent with analysis in
Section II-A. As such, even the proposed TDMA-based scheme with much lower
complexity shows superiority over it. Furthermore, though the performance of
the SCA-based scheme improves to some extent as compared to the GR-based
scheme, it suffers from significant performance loss, which becomes more
prominent when $K$ is larger. This is fundamentally because it is essentially
still a single-beam based transmission, which thus demonstrates the
effectiveness of our proposed dynamic IRS beamforming scheme. It approaches
the upper bound with a performance loss of less than 9%. As compared with
other benchmark schemes, for example, with $K=60$, the total harvested energy
obtained by dynamic IRS beamforming is about $2.4\times 10^{-5}$ J, which
achieves performance gains of 27.5% higher than the SCA-based static
beamforming design, respectively. This thus reveals the advantage of properly
designed joint phase-shift patterns optimization and resource allocation,
highlighting the necessity of multi-beam transmission.

Fig. 4: Total harvested energy versus $K$.

Fig. 5: Total harvested energy versus $J$ for dynamic IRS beamforming scheme.

Fig. 5 studies the total harvested energy of all ERs versus the number of IRS
reflection patterns for the proposed dynamic IRS beamforming scheme under
different system parameters. It is observed that the total harvested energy
increases as $J$ increases, and tends to be saturated. This is expected since
by dynamically adjusting IRS’s passive beamforming, the system is able to
mimic the simultaneous multi-beam transmission through IRS and thus increase
the system performance. Furthermore, when the number of IRS reflection
patterns, i.e., $J$ exceeds the rank of $\bm{\Theta}$ (as shown in Fig. 3(a)),
the performance gain from dynamic IRS beamforming is maximized. Thus the
performance improvement becomes negligible when $J$ becomes larger. In
addition, as shown in Fig. 3(b), the matrix $\bm{\Theta}$ is dominated by the
first few eigenvalues. As such, even if the average rank of $\bm{\Theta}=4.5$
with $K=60$, the system readily achieves a good performance with $2$ or $3$
reflection patterns. This implies that in practice, $J$ can be properly chosen
to balance between the system performance and the signaling overhead incurred
by sending IRS reflection patterns from the ET to IRS.

##  VI Conclusion

In this letter, we investigated the joint passive beamforming and resource
allocation optimization for an IRS-aided multiuser WET system under the
practical non-linear EH model. We first studied static IRS beamforming scheme,
from which we observed that a high rank IRS reflection matrix is generally
needed to optimize the multiuser WET with user fairness consideration. To
mimic the simultaneous multi-beam transmission with IRS single beam patterns,
we then propose the dynamic IRS beamforming scheme and the low-complexity
TDMA-based IRS beamforming scheme. The total harvested energy maximization
problems were formulated, and efficient algorithms based on SDR and SCA
techniques were proposed. Simulation results reveal that the time-division
based dynamic IRS beamforming design can improve the WET efficiency and thus
provide new insight on the passive beamforming design in multiuser WET
systems.

## References

  * [1] Q. Wu and R. Zhang, “Beamforming optimization for wireless network aided by intelligent reflecting surface with discrete phase shifts,” _IEEE Trans. on Commun._ , vol. 68, no. 3, pp. 1838–1851, Mar. 2020.
  * [2] S. Bi, C. K. Ho, and R. Zhang, “Wireless powered communication: opportunities and challenges,” _IEEE Commun. Mag._ , vol. 53, no. 4, pp. 117–125, Apr. 2015.
  * [3] Q. Wu and R. Zhang, “Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 11, pp. 5394–5409, Nov. 2019.
  * [4] S. Hu _et al._ , “Robust and secure sum-rate maximization for multiuser MISO downlink systems with self-sustainable IRS,” _IEEE Trans. on Commun._ , vol. 69, no. 10, pp. 7032–7049, Oct. 2021.
  * [5] S. Zhang _et al._ , “Airis: Artificial intelligence enhanced signal processing in reconfigurable intelligent surface communications,” _China Commun._ , vol. 18, no. 7, pp. 158–171, Jul. 2021.
  * [6] M. Hua and Q. Wu, “Joint dynamic passive beamforming and resource allocation for IRS-aided full-duplex WPCN,” _IEEE Trans. on Wireless Commun._ , pp. 1–1, Aug. 2021.
  * [7] Q. Wu _et al._ , “Intelligent reflecting surface-aided wireless energy and information transmission: An overview,” _Proceedings of the IEEE_ , vol. 110, no. 1, pp. 150–170, Jan. 2022.
  * [8] J. Xu, S. Bi, and R. Zhang, “Multiuser MIMO wireless energy transfer with coexisting opportunistic communication,” _IEEE Wireless Commun. Lett._ , vol. 4, no. 3, pp. 273–276, Jun. 2015.
  * [9] G. Ma _et al._ , “Time-division energy beamforming for multiuser wireless power transfer with non-linear energy harvesting,” _IEEE Wireless Commun. Lett._ , vol. 10, no. 1, pp. 53–57, Jan. 2021.
  * [10] E. Boshkovska _et al._ , “Practical non-linear energy harvesting model and resource allocation for SWIPT systems,” _IEEE Commun. Lett._ , vol. 19, no. 12, pp. 2082–2085, Dec. 2015.
  * [11] S. Boyd and L. Vandenberghe, _Convex optimization_. Cambridge university press, 2004.
  * [12] K.-Y. Wang _et al._ , “Outage constrained robust transmit optimization for multiuser MISO downlinks: Tractable approximations by conic optimization,” _IEEE Trans. on Signal Process._ , vol. 62, no. 21, pp. 5690–5705, Nov. 2014\.
  * [13] S. Zargari _et al._ , “Max-min fair energy-efficient beamforming design for intelligent reflecting surface-aided SWIPT systems with non-linear energy harvesting model,” _IEEE Trans. on Veh. Technol._ , vol. 70, no. 6, pp. 5848–5864, Jun. 2021.


# Spectrum of Feshbach resonances in NaLi $+$ Na collisions

Juliana J. Park<EMAIL_ADDRESS>Research Laboratory of Electronics, MIT-
Harvard Center for Ultracold Atoms, Department of Physics, Massachusetts
Institute of Technology, Cambridge, Massachusetts 02139, USA  Hyungmok Son
Research Laboratory of Electronics, MIT-Harvard Center for Ultracold Atoms,
Department of Physics, Massachusetts Institute of Technology, Cambridge,
Massachusetts 02139, USA  Department of Physics, Harvard University,
Cambridge, Massachusetts 02138, USA  Yu-Kun Lu  Research Laboratory of
Electronics, MIT-Harvard Center for Ultracold Atoms, Department of Physics,
Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA
Tijs Karman  Institute for Molecules and Materials, Radboud University,
Heijendaalseweg 135, 6525 AJ Nijmegen, the Netherlands  Marcin Gronowski
Faculty of Physics, University of Warsaw, Pasteura 5, 02-093 Warsaw, Poland
Michał Tomza  Faculty of Physics, University of Warsaw, Pasteura 5, 02-093
Warsaw, Poland  Alan O. Jamison  Institute for Quantum Computing and
Department of Physics & Astronomy, University of Waterloo, Waterloo, Ontario
N2L 3G1, Canada  Wolfgang Ketterle  Research Laboratory of Electronics, MIT-
Harvard Center for Ultracold Atoms, Department of Physics, Massachusetts
Institute of Technology, Cambridge, Massachusetts 02139, USA

###### Abstract

Collisional resonances of molecules can offer a deeper understanding of
interaction potentials and collision complexes, and allow control of chemical
reactions. Here, we experimentally map out the spectrum of Feshbach resonances
in collisions between ultracold triplet ro-vibrational ground-state NaLi
molecules and Na atoms over a range of 1400 G. Preparation of the spin-
stretched state puts the system initially into the non-reactive quartet
potential. A total of 25 resonances are observed, in agreement with quantum-
chemistry calculations using a coupled-channels approach. Although the theory
cannot predict the positions of resonances, it can account for several
experimental findings and provide unprecedented insight into the nature and
couplings of ultracold, strongly interacting complexes. Previous work has
addressed only weakly bound complexes. We show that the main coupling
mechanism results from spin-rotation and spin-spin couplings in combination
with the anisotropic atom-molecule interaction, and that the collisional
complexes which support the resonances have a size of 30-40 $a_{0}$. This
study illustrates the potential of a combined experimental and theoretical
approach.

##  I Introduction

Collisional resonances that are electromagnetically tunable have become an
established tool for modifying interactions between ultracold atoms and are
the key for many applications, from magnetic association of loosely bound
molecules to quantum simulations [1, 2]. For ultracold molecular systems,
tunable collisional resonances can control chemical reactions [3] and also
provide microscopic information about interaction potentials and collision
complexes.

In the case of cold collisions between alkali-metal atoms, the number of
resonant states remains typically small, and the resonances are usually
tractable. However, in the case of cold collisions involving molecules, due to
strong and anisotropic interactions, ro-vibrational excited states can also
contribute to resonant states, and therefore resonances themselves may not be
well separated and are difficult to identify. Due to the large density of
states of molecular systems [4, 5] it has been a challenge to perform rigorous
scattering calculations and different methods of approximations have become an
active field [6, 7]. Despite these efforts, due to the extreme sensitivity of
the low-temperature observables to details of potential energy surfaces [8,
9], there has been very little success in predicting Feshbach resonances in
molecular collisions.

Rather than performing exact quantum calculations, one feasible alternative
approach is the use of simple statistical short-range models while treating
the physics of long-range scattering within multichannel quantum defect
theory. This has been pursued for the description of molecular resonances [4,
10, 5, 11]. However, the validity of statistical short-range models is
controversial and, therefore, detailed experimental and further theoretical
studies are needed.

Previously, only resonances supported by long-range states of collision
complexes could be assigned to specific quantum states since the quantum
numbers of the separated atoms and molecules are approximately preserved.
Resonant states in collisions involving Feshbach molecules [12, 13, 14, 15]
and in collisions between ${}^{40}\rm{K}$ and ground state
${}^{23}\rm{Na}^{40}\rm{K}$($X^{1}\Sigma^{+}$) [16, 17] were successfully
analyzed in this way. Here we are addressing the challenge of more strongly
interacting complexes by using, for the first time, a combined experimental
and theoretical quantum chemistry approach towards this goal. The agreement
between experiment and theory validates approximations in the calculations,
and the theoretical results can then be used to assign quantum numbers to
resonances and to identify the microscopic mechanism of the resonant
couplings. To keep the problem tractable, we focus on collisions between
triplet ro-vibrational ground-state NaLi and Na, both prepared in the
maximally spin-stretched state. Collisions of
${}^{23}\rm{Na}^{6}\rm{Li}$($a^{3}\Sigma^{+}$) with $\rm{Na}$ are generally
chemically reactive (${}^{23}\rm{Na}^{6}\rm{Li}$($a^{3}\Sigma^{+}$) $+$
$\rm{Na}\rightarrow\rm{Na}_{2}$$(X^{1}\Sigma_{g}^{+})+\rm{Li}+\rm{heat}$).
However, the chemical reaction in the fully spin-polarized atom–molecule
system is strongly suppressed due to the approximate conservation of the total
spin [18, 19], and we can regard NaLi $+$ Na in the quartet potential as
approximately chemically stable. The system is suitable for modeling molecular
scattering resonances because of the relatively small density of states and
number of electrons, and therefore more accurate quantum calculations are
feasible compared to other heavier molecular systems. Although the observed
resonances involve strongly bound complexes, we can make complete assignments
of spin states and total mechanical angular momentum and identify the relevant
coupling terms.

In previous work, we have analyzed the lineshape of one very strong resonance
in collisions between ${}^{23}\rm{Na}$ and
${}^{23}\rm{Na}^{6}\rm{Li}$($a^{3}\Sigma^{+}$) [18]. Here we report a study of
collisions in both spin-stretched states over a magnetic field range of over
1400 G and report more than 20 new Feshbach resonances. This now allows us to
draw conclusions about typical features of collisional complexes. The direct
comparison with quantum chemistry calculations provides major new insight into
the interaction potential and coupling mechanism of the strongly interacting
collision complex.

##  II Experimental protocol and results

Fig. 1: Collisional loss spectrum of NaLi molecules with Na atoms as a
function of magnetic field. Spectra were recorded for both the upper (a) and
lower (b) stretched hyperfine states. Shown is the normalized number of NaLi
molecules left after sweeping down the bias field by 13 G for 200 ms. The
number of Na atoms in each pancake was around 115, and the initial number of
NaLi molecules was about 35. A total of 25 resonances, indicated by red lines,
were observed: 8 in the upper spin-polarized mixture and 17 in the lower spin-
polarized mixture. Each data point represents three to eight measurements.
Error bars correspond to one standard deviation of the mean. Grey lines are a
guide to the eye, obtained by interpolating data with a piecewise cubic
polynomial. Collision channel :
$\mathbf{\ket{7/2,7/2}_{\rm{NaLi}}+\ket{2,2}_{\rm{Na}}}$
---
| $B_{0}\;(\rm{G})$ | $\Delta B\;(\rm{G})$ |  $\beta\;(\rm{cm}^{3}\rm{s}^{-1}$)
1 | 203.7(2) | 5.3(7) | $1.5(6)\times 10^{-11}$
2 | 237.9(1) | 6.3(4) | $7.1(1)\times 10^{-11}$
3 | 372.0(7) | 20(3) | $2.3(5)\times 10^{-11}$
4 | 657.8(3) | 5.2(8) | $1.0(1)\times 10^{-10}$
5 | 884.3(3) | 10(2) | $6(1)\times 10^{-11}$
6 | 978.2(2) | 4.9(3) | $7.4(5)\times 10^{-10}$
7 | 1029.7(3) | 10(2) | $1.5(3)\times 10^{-11}$
8 | 1135.2(8) | 15(3) | $3.0(5)\times 10^{-11}$
Collision channel : $\mathbf{\ket{7/2,-7/2}_{\rm{NaLi}}+\ket{2,-2}_{\rm{Na}}}$
| $B_{0}\;(\rm{G})$ | $\Delta B\;(\rm{G})$ |  $\beta\;(\rm{cm}^{3}\rm{s}^{-1}$)
1 | 82.5(1) | 1.3(2) | $1.0(3)\times 10^{-10}$
2 | 145(1) | 18(4) | $5.3(9)\times 10^{-11}$
3 | 309(1) | 13(4) | $1.9(4)\times 10^{-11}$
4 | 361.82(5) | 1.0(2) | $4.5(8)\times 10^{-11}$
5 | 449.8(2) | 3.1(5) | $6(2)\times 10^{-11}$
6 | 478(3) | 35(13) | $3.4(6)\times 10^{-11}$
7 | 561.8(2) | 2.4(6) | $8(2)\times 10^{-11}$
8 | 590.7(2) | 9(1) | $8(2)\times 10^{-11}$
9 | 642.3(5) | 3(2) | $2.9(4)\times 10^{-11}$
10 | 661.5(2) | 3.9(9) | $2.4(5)\times 10^{-11}$
11 | 722.5(1) | 2.7(3) | $8(4)\times 10^{-11}$
12 | 860.8(2) | 2.1(6) | $2.4(6)\times 10^{-10}$
13 | 963.3(2) | 6(1) | $2(1)\times 10^{-10}$
14 | 1030.1(2) | 3.3(8) | $9(3)\times 10^{-11}$
15 | 1083.3(3) | 4(1) | $5(1)\times 10^{-11}$
16 | 1176.3(3) | 3.2(9) | $8(2)\times 10^{-11}$
17 | 1269.2(1) | 2.6(3) | $3.3(5)\times 10^{-11}$
Table 1: Observed Feshbach resonances. Resonance positions and widths are
obtained by Lorentzian fits to loss features, except for resonances 4 to 6 for
the incoming collision channel
$\ket{7/2,7/2}_{\rm{NaLi}}+\ket{2,2}_{\rm{Na}}$. For these resonances, the
position and width are determined by Lorentzian fits to the field-dependent
loss rates [18].

We prepare a mixture of triplet ground-state
${}^{23}\text{Na}^{6}\text{Li}$$(a^{3}\Sigma^{+})$ molecules and
${}^{23}\text{Na}$ atoms in the spin-polarized quartet potential. There are
two possible states: the “upper stretched state”, where all nuclear and
electron spins are aligned to the bias magnetic field direction
($\ket{F,M_{F}}_{\rm{NaLi}}+\ket{F,M_{F}}_{\rm{Na}}=\ket{7/2,7/2}_{\rm{NaLi}}+\ket{2,2}_{\rm{Na}}$)
and the “lower stretched state”, where all nuclear and electron spins are
anti-aligned to the field direction
($\ket{7/2,-7/2}_{\rm{NaLi}}+\ket{2,-2}_{\rm{NaLi}}$). Here, $F$ is the
quantum number for the total spin (electron and nuclear) and $M_{F}$ is the
$B$-field projection of $F$. The molecule and atom mixture in the upper
stretched state with typical numbers of ${{\sim}}3{\times}10^{4}$ and
${{\sim}}3{\times}10^{5}$, respectively, is produced in a $1596$-nm one-
dimensional optical lattice following the method described in Ref. [20, 18].
The lower stretched state is produced by coherent transfer from the sample in
the upper stretched state using a magnetic field sweep in the presence of
radio frequency waves [21]. For this process, the bias field is dropped from
$745$ G, where the upper stretched state is prepared, to a low field of around
$8$ G in $15$ ms. After state preparation, the bias field is ramped to a
target value in $15\;\text{ms}$. Collisional lifetimes of the atom-molecule
mixtures are determined by holding the sample for a variable time at the
target magnetic field.

The loss of NaLi molecules with Na atoms is measured as a function of the bias
field for both spin-polarized states. First, we performed a coarse search by
sweeping the bias field over a range of 12.6 G during 200 ms and recording the
number of remaining molecules, normalized to the number without the field
sweep. This procedure was repeated over a range from near zero to 1420 G in
steps of 8.76 G, with the results shown in Fig. 1. We identified 8 resonances
in the upper stretched state and 17 resonances in the lower stretched state
indicated with red vertical lines in Fig. 1.

Near the resonances found in the coarse scan, finer scans were performed with
sweep range and step size of around 1 G.

Fig. 2: Feshbach resonance near 884 G. This resonance was observed for NaLi +
Na in the upper stretched state. (a) Normalized NaLi number as a function of
bias field after a hold time of 100 ms. (b) Loss rate of NaLi molecules
obtained from decay curves as shown in the inset at 884.3 G. Data values and
error bars in (a) and inset of (b) represent the average and one standard
deviation of the mean, estimated from statistical errors of three to six
measurements respectively. Error bars in the main plot of (b) represent one
standard deviation of a fitted decay parameter. The blue lines are the fits to
a Lorentzian function plus a background with a linear slope.

Using the data from the fine scans, each loss feature was fit to a Lorentzian
with a slope accounting for background loss and nearby resonances [an example
is shown in Fig. 2(a)]. For all resonances, the thus obtained peak positions
and widths are listed in Table 1. Strong losses lead to a broadened lineshape.
Therefore, for three resonances (labelled 4, 5 and 6 in the upper stretched
state), we also determined the atom loss as a function of hold time and
determined loss rates as a function of the bias field. These results were used
for a more accurate determination of the resonance position and width, as
shown in Fig. 2(b). We confirmed that the width obtained from observed loss
features can be broader than the width from the loss rate measurement. Fig. 2
illustrates this for resonance 5 in the upper stretched state which is near
884 G. The width obtained from a Lorentzian fit to the loss feature is 17(1)
G, whereas the loss rate measurements give a width of 10(2) G.

For all resonances, we determined peak loss rate constants $\beta$ by
recording decay curves as a function of the hold time near the center of the
loss features. They are summarized in the last column of Table 1. The
determination of loss rate constants requires knowledge of the densities of
sodium atoms overlapped with molecules. Instead of absolute calibration, we
follow Ref. [18] and compare measured decay rates with the decay rate of the
mixture in a non-stretched spin state which occurs at the universal rate [22].
This approach was validated in Ref. [18].

Away from resonances, the observed background loss rates are more than an
order of magnitude smaller than the universal loss rate constant which for Na
$+$ NaLi $s$-wave collisions is $1.7\times 10^{-10}\rm{cm}^{3}\rm{s}^{-1}$.
More specifically, at low field near 0.5 G, the upper and lower stretched Na
$+$ NaLi loss rate constants are $7.5(2.2)\times
10^{-12}\rm{cm}^{3}\rm{s}^{-1}$ and $6.7(2.0)\times
10^{-12}\rm{cm}^{3}\rm{s}^{-1}$, respectively, and at high field (near 540 G),
they are $4.5(1.4)\times 10^{-12}\rm{cm}^{3}\rm{s}^{-1}$ and $3.5(1.0)\times
10^{-12}\rm{cm}^{3}\rm{s}^{-1}$, respectively. Within the accuracy of
measurement, the two stretched states have the same loss rate constants, with
lower background loss rates at higher magnetic fields.

##  III Analysis

In this section, we summarize the experimental findings and draw some
conclusions. First, we have observed similar rates of background loss and
Feshbach enhanced losses for both stretched states. This implies that Zeeman
energies do not play a major role and that dipolar relaxation is not the
dominant decay mechanism. Instead, the decay is probably dominated by shorter-
range chemical reactions or inelastic loss, which are expected to be similar
for both spin-stretched initial states. Regarding the number of resonances, we
have observed more resonances for the lower stretched state. In the range up
to 1400 G, we observed 8 resonances in the upper stretched state and 17 in the
lower stretched state. Intuitively this observation can be rationalized that
for the lower stretched state, the non-initial Zeeman states are closed
channels that can support additional resonance states. Surprisingly, our
quantum chemistry calculations below show that the difference is instead
caused by fast decay of some upper stretched state resonances that become too
broad to be observed.

Next, we discuss the widths and spacings between resonances. The linewidths of
the resonances range from about 1 G to about 30 G. If the linewidths are
interpreted as due to the finite lifetime of the resonant state, this implies
atom-molecule complexes with lifetimes in the range of 10 to 350 ns. The
average spacing between resonances is around 100 G. One possible
interpretation of the spacing is rotational structure of the intermediate
complex. The spacing between resonances of the order of 100 G corresponds to
single spin-flip energy differences of 280 MHz. A rotational constant of this
value corresponds to the moment of inertia of a complex with a size (i.e.,
atom/molecule separation) of around 30 $a_{0}$. Note that the rotationally
excited states of the NaLi triplet molecules with $N=1$ are at 9.25 GHz and
$N=2$ are at 27.75 GHz [23] corresponding to a double spinflip energy at 1650
G and 5000 G. In the accompanying paper, we show that this may lead to strong
resonant features [24], but these occur far outside the magnetic field range
studied here.

Finally, we can perform a statistical analysis of the distribution of nearest-
neighbor resonance spacings. The observed distribution follows the Wigner-
Dyson distribution [25, 26], which is a feature of a chaotic system (see
appendix). However, statistical conclusions from only 17 resonances are only
tentative.

##  IV Coupled-channels calculations

With the coupled-channels calculations we will try to answer the following
questions: (1) What is the mechanism of coupling between the initial
scattering channel and the loss channels for the background loss? (2) What are
the dominant interactions that are responsible for the resonances? (3) Why are
there more observed resonances for the lower stretched state? (4) What quantum
numbers describe the collisional complex? (5) What is the size of the
collisional complex?

Fig. 3:  Calculated spectrum of Feshbach resonances. Plotted are the
calculated resonances for the upper and lower stretched states for
$N_{\mathrm{max}}=30$ and $\lambda=-0.02$. The comparison with Fig. 1 shows
that most experimental features (except for resonance positions) are
reproduced by the calculations.

The main result is shown in Fig. 3. A calculated spectrum reproduces most of
the observed features (magnitude of loss rate, widths and number of
resonances, more resonances for the lower stretched states), but cannot
predict the observed resonance positions. This illustrates the power (and
limitations) of state of the art quantum chemistry calculations.

These calculations are fully quantum mechanical coupled-channel calculations
including the electronic interactions, the Zeeman interaction with the
magnetic field, the magnetic dipole-dipole interaction, and the spin-rotation
and spin-spin couplings. The electronic interaction is partitioned into two-
body interactions that are taken from experiment [27, 28, 29], and a non-
additive three-body interaction that is calculated with state-of-the-art
coupled-cluster methods including single, double, and triple excitations with
large Gaussian basis sets extrapolated to the complete basis set limit that
are described in the accompanying paper [24]. The interaction-induced
variation of spin-rotation and spin-rotation couplings is neglected. We assume
that the molecular bond length is fixed at the equilibrium position of the
triplet potential, which naively cannot explain the chemical reactions that
form singlet NaLi or Na2 molecules at the low-spin potential. In our coupled-
channel calculations, we model these by imposing an absorbing boundary
condition at $R=4.5~{}a_{0}$, which can be reached on the low-spin potential,
but not on the high-spin potential, which is highly repulsive at these short
distances.

For simplicity, we start by ignoring hyperfine and vibrational degrees of
freedom. The scattering wavefunction is expanded in the basis of fully coupled
channel functions of the form

$\ket{(NL)J(s_{\mathrm{mol}},s_{\mathrm{atom}})S;\mathcal{J}\mathcal{M}}=\sum_{M_{J},M_{S}}{\braket{JM_{J}SM_{S}}{\mathcal{J\mathcal{M}}}}\\\
\times\ket{(NL)JM_{J}}\ket{(s_{\mathrm{mol}}\ s_{\mathrm{atom}})SM_{S}}$ (1)

where $\braket{JM_{J}SM_{S}}{\mathcal{J\mathcal{M}}}$ is a Clebsh-Gordan
coefficient. The quantum number $N$ represents the rotational angular momentum
of the NaLi molecule, and $L$ the angular momentum associated with the
relative motion of the atom and molecule. $N$ and $L$ are Clebsch-Gordan
coupled to a total mechanical angular momentum $J$ with $z$-component $M_{J}$.
Similarly, $s_{\mathrm{mol}}=0$ or $1$ denotes the NaLi molecular electronic
spin and $s_{\mathrm{atom}}={}^{1}/_{2}$ is the atomic electronic spin, and
$S$ the total electronic spin with $B$-field projection $M_{S}$. In the
coupled basis, $J$ and $S$ are subsequently coupled to a total angular
momentum $\mathcal{J}$ and a magnetic field projection
$\mathcal{M}=M_{J}+M_{S}$. $\mathcal{M}$ is strictly conserved, whereas, for a
large enough magnetic field, $M_{S}$ becomes a good quantum number and
therefore $M_{J}=\mathcal{M}-M_{S}$ is also a good quantum number. Due to the
large singlet-triplet splitting in the NaLi molecule, $s_{\mathrm{mol}}=0$ or
1 is also a good quantum number. For a separated atom and molecule,
$m_{s_{\mathrm{mol}}}$ and $m_{s_{\mathrm{atom}}}$ would separately become
good quantum numbers, but at chemically relevant distances the exchange
splitting between the doublet and quartet interaction potentials is dominant,
so that $S=1/2$ and $3/2$ are good quantum numbers. Hence, we can effectively
consider each $\ket{S\ M_{S}}$ state separately, with only perturbatively weak
couplings between them. For each of these spin channels, there are strong and
anisotropic interactions that couple different $N$ and $L$ channels but
conserve $J$ and $M_{J}$. Since we ignore nuclear spin, the initial channel
corresponds to $s$-wave collisions in the upper spin-stretched ground state,
$\ket{(NL)JM_{J}}\ket{SM_{S}}=\ket{(00)00}\ket{3/2\;3/2}$, or the lower spin-
stretched state $\ket{(00)00}\ket{3/2,\;-3/2}$.

First, we consider the sensitivity of the scattering rates to the interaction
potential shown in Figure 4(a). Here, we scale by a factor $1+\lambda$ the
non-additive three-body part of the interaction potential, that is, the part
that is computed _ab initio_ and is uncertain up to an conservatively
estimated 5 %. By modifying the potential by only 0.1 % we find that the
resonances start to shift so that realistically their positions are completely
undetermined, and when the scaling reaches several percent we tune across
magnetic field-independent resonances, indicating that the background
scattering length is undetermined. Next, shown in Figure 4(b), we again scale
the three-body interaction but now only for the low-spin doublet potential,
leaving the high-spin quartet potential unchanged. In this case, we find that
several of the resonances are now completely independent of the scaling of the
low-spin potential up to $\lambda=0.1$. This implies that no stable resonance
states are supported by the chemically reactive low-spin potential and that
all the predicted (and observed) resonances originate in the quartet
potential.

The analysis above indicates that the _ab initio_ prediction of resonance
positions is beyond the capability of state-of-the-art calculations. Although
the coupled-channels calculations cannot predict the positions of the
resonances, we can still use these calculations as a “numerical experiment” to
investigate the nature of the resonance states, the coupling mechanisms, and
the observed differences between the two spin-stretched states.

Fig. 4:  Calculated loss rates as a function of magnetic field. Calculations
are done with scaling of the spin-independent three-body interaction by a
factor $1+\lambda$ for $N_{\rm{max}}=6$. (a) Scaling of the three-body
interaction for both the doublet and quartet state. The figure shows that
scaling the three body interactions within their uncertainty of at least
several percent dramatically changes the loss rates: The positions of several
Feshbach resonances change, and several $B$-independent resonances appear when
the $\lambda$-scaling creates a resonance in the initial spin-stretched
potential, i.e., a bound state near zero energy. The conclusion is that the
prediction of resonance positions and background loss requires knowledge of
the interaction potentials to an accuracy better than can be achieved by _ab
initio_ calculations. (b) Scaling of the three-body interaction only for the
doublet potential. In this case, the resonances follow vertical lines and are
independent of the scaling up to $\lambda=0.1$, implying that the resonances
originate in the quartet potential.

We perform coupled-channels scattering calculations with the interactions
scaled by $1+\lambda$ and analyze the _typical_ behavior observed for
different $\lambda$ between $-0.1$ and $+0.1$. Representative magnetic field
scans can be seen in Fig. 5(a). We observe approximately 10 resonances for the
lower spin-stretched state and only around 5 resonances in the upper spin-
stretched state, respectively. This is in qualitative agreement with the
experiment which observes 17 and 8 resonances, respectively. In the companion
paper we show that based on the density of states we would expect to see
approximately 10 resonances for either spin-stretched state [24]. We
furthermore show that including molecular vibrations would increase the
density of state by approximately 50 %, which can partially explain the lower
number of observed resonances compared to experiment. We show that including
hyperfine interactions increases the density of resonances somewhat. With
these effects in mind, one could claim almost quantitative agreement with
experiment regarding the density of resonances.

One may expect that the lower spin-stretched state supports more magnetically
tunable resonances because the non-initial Zeeman states correspond to closed
channels – and hence can support Feshbach resonances – even in the rotational
ground state, whereas closed channels for the upper spin-stretched state occur
only for excited rotational states. We investigate this in our coupled-
channels calculations by artificially excluding the non-initial Zeeman states
in the rotational ground state, see Fig. 5(b). Somewhat surprisingly we find
that _excluding_ channels from the calculation does not reduce the number of
resonances for the lower spin-stretched state, but rather increases the
observed number of resonances for the upper spin-stretched state, where the
excluded channels correspond to asymptotically open channels. In the presence
of these open channels some of the resonances decay rapidly by spin-spin
coupling to lower Zeeman levels (both in the doublet and quartet potentials)
such that they are not resolved, leading to a lower number of observable
resonances. This explains the observed qualitative difference between the
upper and lower spin-stretched states.

Fig. 5:  Calculated spectrum of Feshbach resonances. Figure (a) shows results
for the upper and lower spin-stretched states in orange and blue solid lines
(same data as in Figure 3). The lower spin-stretched state typically shows
around 10 resonances below 1500 G, whereas the upper spin-stretched state
shows around half as many resonances. This qualitative difference between the
two states is also observed experimentally. We investigate whether the
anisotropy of the electronic interaction and the spin-spin and spin-rotation
coupling act perturbatively, by scaling down these couplings by a factor of
two and scaling the resulting cross section up by a factor four. Agreement
with the full calculation indicates the spin-spin and spin-rotation coupling
act perturbatively, whereas the interaction anisotropy does not. The coupling
mechanism however does involve the anisotropy as turning this off entirely
produces a much smaller cross section dominated by the magnetic dipole-dipole
interaction (dash-dotted line). When both anisotropy and dipole-dipole are
turned off, the calculated cross section is zero. (b) Removal of Zeeman
states. The qualitative differences in the number of resonances for the upper
and lower spin-stretched states disappear when we _exclude_ from the
calculation channels corresponding to Zeeman states in the rotational ground
state, except for the initial one.

Finally, we investigate numerically the coupling mechanism that gives rise to
the observed resonances. As argued above, each resonance state can essentially
be assigned a molecular and a total electron spin quantum number
$s_{\mathrm{mol}}=1$ and $S={}^{3}/_{2}$, as only the non-reactive quartet
spin state supports stable resonance states. The resonance is magnetically
tunable only if the Zeeman state $M_{S}$ changes. To couple states with
$\Delta M_{S}\neq 0$ a spin-dependent coupling must be involved through spin-
rotation and spin-spin coupling.

From the tensor rank of these couplings we can determine they couple states
with the selection rules, $J=0\rightarrow 1$ and $|\Delta M_{S}|\leq 1$ (for
spin-rotation coupling), and $J=0\rightarrow 2$ and $|\Delta M_{S}|\leq 2$
(for spin-spin coupling), respectively. The spin-spin couplings and spin-
rotation couplings arise from terms
$(\hat{s}_{\mathrm{mol}}\cdot\hat{r})(\hat{s}_{\mathrm{mol}}\cdot\hat{r})$ and
$\hat{s}_{\mathrm{mol}}\cdot\hat{N}$ for the NaLi molecule, where $\hat{r}$
points along the molecular axis and $\hat{s}_{\mathrm{mol}}$ is the molecular
electronic spin. Since they depend on the orientation of the molecular axis,
they cause exchange of spin angular momentum with the molecular rotation $N$.

Since $J$ is approximately a good quantum number, both mechanisms give rise to
distinct and independent resonances. Since the density of states of the
collision complex scales with $2J+1$ and the differential magnetic moment is
higher for larger $|\Delta M_{S}|$ transitions, we conclude that most –
approximately 5/6 – of the resonances are due to spin-spin coupling, and the
remaining 1/6 is due to spin-rotation coupling. Spin-spin coupling does not
occur for ${}^{2}\Sigma$ molecules, and hence is a somewhat unique coupling
mechanism for NaLi in the triplet ground state (${}^{3}\Sigma$).

To confirm the role of the spin-dependent interactions in the coupling
mechanism we have performed coupled-channels calculations where we reduced
these couplings by a factor two. The resulting cross section multiplied by
four is shown as the crosses in Fig. 5(a). The agreement with the full
calculation indicates scaling with the square of the coupling strength that is
expected for perturbatively weak spin-dependent couplings. Note that on
resonance, the dependence on the coupling strength is not an overall scaling,
as the coupling strength also determines the resonance widths. For smaller
coupling, the peaks are narrower and higher. Those results confirm that we can
fully assign the resonances approximately good quantum numbers
$s_{\mathrm{mol}}=1$, $S={}^{3}/_{2}$, each resonance also has definite
$M_{S}$ constrained by $\Delta M_{S}\leq 1$ and $\leq 2$ selection rules, and
$J=1$ or $2$, for spin-rotation and spin-spin coupling, respectively, whereas
the $N$ and $L$ quantum numbers are strongly mixed due to the anisotropic
interaction at short range. Figure 5(a) also implies that the loss mechanisms
are the same for resonant losses and background losses. For the value of
$\lambda$ chosen in Fig. 5(a), the background loss is higher for the upper
spin-stretched state, but for most values of $\lambda$, they are similar, as
observed in the experiment [24].

The coupling mechanism involving the spin-rotation and spin-spin coupling also
requires an anisotropic interaction potential. The physical picture is that
the anisotropic interaction with the atom can reorient the molecule, and
because the spin is coupled to the molecular axis by spin-rotation and spin-
spin coupling, this can lead to Zeeman transitions. To confirm this picture we
have performed calculations that exclude interaction anisotropy, shown as the
dash-dotted line in Fig. 5(a). The resulting cross section is much smaller and
results from long-range Zeeman relaxation by the magnetic dipole-dipole
coupling. If both the magnetic dipole-dipole interaction and the interaction
anisotropy are switched off, the cross section in our model vanishes. The role
of the interaction anisotropy is non-perturbative, however, as can be seen
from comparison between the solid and dotted line in Fig. 5(a), which compares
the cross section from the full calculation to four times the cross section
obtained with the interaction anisotropy halved. The crucial and non-
perturbative role of the anisotropic electronic interaction implies that the
spectrum of resonances cannot be described by a simplified model that accounts
only for the isotropic long-range $R^{-6}$ interaction, contrary to previous
observations of ultracold atom-molecule resonances [17]. In summary: The
observed losses arise from the interplay of spin-spin and spin-rotation
coupling with the anisotropic interactions which connect the incoming non-
reactive quartet states to doublet states that can decay to singlet NaLi or
Na2 molecules. This coupling can occur directly, or via first coupling to non-
stretched quartet states.

We make a direct comparison to the scattering calculations with the density of
states computed quantum mechanically using the same channel basis as used in
the scattering calculations in the accompanying paper [24]. The total number
of states from the quasiclassical estimate is close to the number of
resonances from the scattering calculations, but not in perfect agreement with
it, because of the light masses and relatively weak interactions in the spin-
stretched state. We find that most of the 3-body states are between 20 and 40
$a_{0}$ which is in agreement with the complex size estimated from the simple
interpretation of the spacing between resonances as rotational energy
splitting of the collision complex. This size is short compared to the range
of the van der Waals potential, but much longer than the short-range
interactions - the minimum of the potential wells is around 10 $a_{0}$. The
resonances depend both on short-range and long-range physics: They are
supported by the long-range potential, but require the anisotropic short-range
interactions as a coupling mechanism.

The calculated widths of the resonances vary between 1 and 30 G, in
qualitative agreement with the experiment. Those widths reflect the lifetime
since there are no other broadening mechanism in our simulations. Experimental
broadening such as finite resolution of magnetic field strength and magnetic
field inhomogeneity is very small (estimated to be around 100 mG near 1400 G).

##  V Conclusion and Outlook

We have provided a large-scale map of Feshbach resonances between ultracold
atoms and molecules. We have chosen collisions between NaLi in the triplet
ground state and Na for the two fully spin-stretched initial collision
channels in the quartet potential. Although the NaLi $+$ Na mixture is
chemically reactive, there is no reactivity in the quartet potential, i.e.
inelastic collisions require a spinflip to the doublet potential [18, 19]. The
possibility of metastable collision complexes makes this system promising for
the study of Feshbach resonances. Indeed, we have observed 8 and 17 resonances
within a ${\sim}$1400 G range in the upper and lower stretched states,
respectively.

We have compared the experimental results to full coupled-channel calculations
for the NaLi$+$Na mixture. Even state-of-the-art quantum chemistry
calculations cannot predict the position of resonances because of the
uncertainty in the interaction potentials. However, they can provide a deep
understanding of the relevant states and their couplings. Our simulations have
shown that the resonance states are 3-body complexes in the quartet state with
rotational excitation, either of NaLi (quantum number $N$), or rotation of Na
around NaLi (quantum number $L$). These complexes have a typical size of 30-40
$a_{0}$. The input state and the complex state involve different $M_{S}$,
which is necessary for a magnetic Feshbach resonance. Molecular eigenstates in
different Zeeman levels can be coupled by a strong anisotropic electronic
interaction since the different molecular Zeeman states have different
rotational state decompositions mainly due to spin-spin coupling and also
spin-rotation coupling. Note that the spin-rotation and spin-spin coupling
depend on the orientation of the molecular axis and, therefore, provide
tensorial couplings between the spin and the mechanical angular momentum.
Approximately 5/6 of the resonances can be attributed to spin-spin coupling to
states with total mechanical angular momentum $J=2$, and the remaining 1/6 is
attributed to spin-rotation coupling and assigned $J=1$.

The loss of molecules, either the background loss or via Feshbach resonances,
occurs through transition to the chemically reactive doublet potential, either
directly or through lower Zeeman states of the quartet potential (for the
upper stretched state). These transitions involve an interplay of the
anistotropic electronic interaction with spin-spin or spin-rotation coupling
in the NaLi molecules. Dipolar transitions are considerably weaker.

Our coupled-channels calculations explain why the number of observed
resonances is smaller for the upper stretched state. The upper stretched state
supports a similar number of resonances, but several of them are washed out by
the rapid decay of the intermediate complexes to lower-lying Zeeman states via
spin-spin coupling.

The agreement between experiment and calculations validates the assumptions
made in the calculations. At least for a semi-quantitative analysis, it is
sufficient to neglect hyperfine interactions and vibrationally excited state
of the molecules, and use electronic potentials at fixed nuclear separation
for the molecule.

The physical mechanisms we have identified in the Na $+$
NaLi($a^{3}\Sigma^{+}$) mixture are unique to molecules in the triplet or
higher-spin state. It would be interesting to study other bi-alkali triplet
molecules and check if their collision dynamics is dominated by similar
couplings. The only other case where atom-molecule Feshbach resonances have
been studied involved NaK in the singlet ground state [16, 17]. For this
system, the observed Feshbach resonances were caused by long-range van der
Waals interactions [17, 30]. These long-range states then follow regular
patterns dictated by the quantum defect theory for a simple isotropic long-
range $R^{-6}$ interaction [31]. The resonances described here, however,
behave quite differently as the coupling mechanism involves perturbatively the
spin-spin and spin-rotation coupling, and the strong non-perturbative
anisotropic interaction at short range. Although the resulting resonances are
supported by the long-range interaction, their positions are sensitive to this
anisotropic short-range interactions and hence cannot be described by the
regular patterns predicted by quantum defect theory.

Our work here has focused on the non-reactive quartet potential to support
metastable collisional complexes. However, recent studies have shown that
collisional resonances and collisional complexes should also occur in highly
reactive systems. A single $p$-wave Feshbach resonance has been observed in
collisions between NaLi molecules in the triplet ground state [21]. This
system has no barrier for reactions. Losses much smaller than the universal
loss rate have been observed in s-wave collisions between magnetically trapped
NaLi molecules in different hyperfine states of the triplet state [32], a
system which should be highly reactive. Losses far below the universal rate
imply a finite reflection probability at short range and should lead to
Feshbach resonances when the magnetic field tuning creates strong interference
between reflections at short and long range [18, 33].

These studies emphasize that our understanding of reactivity and metastable
collisional complexes is incomplete, and further experimental explorations are
required. Further studies should involve NaLi in various hyperfine states and
with various collision partners, as well as considering other bi-alkalis. As
suggested in Ref. [30], for non-stretched states, the Fermi contact
interaction between electron and nuclear spin provides another mechanism for
Feshbach resonances. Such broader studies are required to find out which
mechanisms are universal, and which occur only in specific systems.

###### Acknowledgements.

The MIT work was supported by the NSF through the Center for Ultracold Atoms
and Grant No. 1506369 and from the Air Force Office of Scientific Research
(MURI, Grant No. FA9550-21-1-0069). J. J. P. and H. S. acknowledge additional
support from the Samsung Scholarship. M. G. and M. T. were supported by the
National Science Centre Poland (Sonata Bis Grant no. 2020/38/E/ST2/00564) and
the PL-Grid Infrastructure (computational grant no. PLG/2020/014342).

##  Appendix A Statistical analysis

A statistical analysis of the separations between resonances can provide
insight into the nature of the resonant states. We quantify the resonance
statistics of 17 resonances from Na+NaLi collisions in the lower stretched
state provided in Table 1 using the Brody parameter $\eta$, which is a
standard measure of chaos. For non-chaotic systems, in which resonances have
no correlations, the distribution of nearest-neighbor spacings is given by the
Poisson distribution, $P_{\text{P}}(s)=e^{-s}$. On the other hand, for chaotic
systems, which emerge when the mean spacing between bound states is comparable
to the coupling strength, repulsion between energy levels occurs. In this
regime, the distribution of nearest-neighbor spacings is given by a Wigner-
Dyson distribution, $P_{\text{WD}}(s)=\frac{\pi}{2}se^{-\frac{\pi}{4}s^{2}}$
[25, 26]. The two distributions are smoothly interpolated through the Brody
parameter $\eta$ as $P_{\text{B}}(s)=As^{\eta}e^{-\frac{A}{\eta+1}\cdot
s^{\eta+1}}$, known as the Brody distribution, where
$A=(\eta+1)\cdot\Gamma\left(\frac{\eta+2}{\eta+1}\right)^{\eta+1}$ [34, 35].
Here, $\eta=1$ and $\eta=0$ lead to the Wigner-Dyson and Poisson
distributions, respectively. The cumulative probability function of the Brody
distribution is given as

$F_{B}(s)=(\eta+1)\cdot[1-e^{-\alpha s^{\eta+1}}].$ (2)

Fig. 6 shows the cumulative probability of resonance spacing of the 17
resonances from Na$+$NaLi collisions in the lower stretched state. The best
fit of the data to Eq.2 gives $\eta=1.1(1)$, which shows the statistical
signature of quantum chaos.

Fig. 6: Cumulative probability of resonance spacings. Shown are the
experimental results for Na+NaLi collisions in the lower stretched state.
Black line is the best fit to the cumulative probability function of the Brody
distribution (Eq.2) and the blue (red) dashed line is the cumulative
distribution with $\eta=1$ ($\eta=0$).  Fig. 7:  Distribution of scaled
reduced widths for scattering resonances. This plot is based on the data for
the lower stretched state Na+NaLi mixture listed in Table 1. The black line
shows the Porter-Thomas distribution of scaled reduced widths.

However, statistical conclusions from only 17 resonances are tentative. On the
basis of our quantum-scattering calculations, the Wigner-Dyson statistics is
not expected because each resonance can be assigned by $M_{S}$, $J$, and
$M_{J}$ quantum numbers (resonances with different values for these quantum
numbers do not affect each other). In addition, it is difficult to conclude
that the broad trend of the distribution of resonance widths follows the
Porter-Thomas distribution of resonance widths,
$P_{\text{PT}}(\bar{\gamma})=\bar{\gamma}^{-1/2}e^{-\bar{\gamma}/2}$, which is
also a statistical character of quantum chaos. Here, the scaled reduced
widths, $\bar{\gamma}=\gamma_{n}^{0}/\langle\gamma_{n}^{0}\rangle$ where
$\gamma_{n}^{0}$ is the resonance width normalized by the resonance energy,
and $\langle\gamma_{n}^{0}\rangle$ is the average of $\gamma_{n}^{0}$ [36].
Fig. 7 shows the histogram of the scaled reduced resonance width distribution
and the Porter-Thomas distribution. Due to the small statistical sample, it is
difficult to compare the broad trend of the distribution with the Porter-
Thomas distribution. Nevertheless, Fig. 6 represents a purely empirical
analysis of experimental data, which should find a theoretical explanation
(which is probably not quantum chaos).

## References

  * Chin  _et al._ [2010] C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Feshbach resonances in ultracold gases, Reviews of Modern Physics 82, 1225 (2010).
  * Bloch  _et al._ [2012] I. Bloch, J. Dalibard, and S. Nascimbene, Quantum simulations with ultracold quantum gases, Nature Physics 8, 267 (2012).
  * Zhao and Pan [2022] B. Zhao and J.-W. Pan, Quantum control of reactions and collisions at ultralow temperatures, Chemical Society Reviews (2022).
  * Mayle  _et al._ [2012] M. Mayle, B. P. Ruzic, and J. L. Bohn, Statistical aspects of ultracold resonant scattering, Physical Review A 85, 062712 (2012).
  * Christianen  _et al._ [2019] A. Christianen, T. Karman, and G. C. Groenenboom, Quasiclassical method for calculating the density of states of ultracold collision complexes, Physical Review A 100, 032708 (2019).
  * Morita and Tscherbul [2019] M. Morita and T. V. Tscherbul, Restricted basis set coupled-channel calculations on atom-molecule collisions in magnetic fields, The Journal of Chemical Physics 150, 074110 (2019).
  * Koyu  _et al._ [2022] S. Koyu, R. Hermsmeier, and T. V. Tscherbul, Total angular momentum representation for state-to-state quantum scattering of cold molecules in a magnetic field, The Journal of Chemical Physics 156, 034112 (2022).
  * Wallis  _et al._ [2011] A. O. Wallis, E. J. Longdon, P. S. Żuchowski, and J. M. Hutson, The prospects of sympathetic cooling of nh molecules with li atoms, The European Physical Journal D 65, 151 (2011).
  * Morita  _et al._ [2019] M. Morita, R. V. Krems, and T. V. Tscherbul, Universal probability distributions of scattering observables in ultracold molecular collisions, Physical review letters 123, 013401 (2019).
  * Mayle  _et al._ [2013] M. Mayle, G. Quéméner, B. P. Ruzic, and J. L. Bohn, Scattering of ultracold molecules in the highly resonant regime, Physical Review A 87, 012709 (2013).
  * Christianen  _et al._ [2021] A. Christianen, G. C. Groenenboom, and T. Karman, Lossy quantum defect theory of ultracold molecular collisions, Physical Review A 104, 043327 (2021).
  * Chin  _et al._ [2005] C. Chin, T. Kraemer, M. Mark, J. Herbig, P. Waldburger, H.-C. Nägerl, and R. Grimm, Observation of feshbach-like resonances in collisions between ultracold molecules, Physical review letters 94, 123201 (2005).
  * Knoop  _et al._ [2009] S. Knoop, F. Ferlaino, M. Mark, M. Berninger, H. Schöbel, H.-C. Nägerl, and R. Grimm, Observation of an efimov-like trimer resonance in ultracold atom–dimer scattering, Nature Physics 5, 227 (2009).
  * Zenesini  _et al._ [2014] A. Zenesini, B. Huang, M. Berninger, H.-C. Nägerl, F. Ferlaino, and R. Grimm, Resonant atom-dimer collisions in cesium: Testing universality at positive scattering lengths, Physical Review A 90, 022704 (2014).
  * Wang  _et al._ [2019] F. Wang, X. Ye, M. Guo, D. Blume, and D. Wang, Observation of resonant scattering between ultracold heteronuclear feshbach molecules, Physical Review A 100, 042706 (2019).
  * Yang  _et al._ [2019] H. Yang, D.-C. Zhang, L. Liu, Y.-X. Liu, J. Nan, B. Zhao, and J.-W. Pan, Observation of magnetically tunable feshbach resonances in ultracold ${}^{23}\rm{Na}^{40}\rm{K}+^{40}\rm{K}$ collisions, Science 363, 261 (2019).
  * Wang  _et al._ [2021] X.-Y. Wang, M. D. Frye, Z. Su, J. Cao, L. Liu, D.-C. Zhang, H. Yang, J. M. Hutson, B. Zhao, C.-L. Bai, _et al._ , Magnetic feshbach resonances in collisions of ${}^{23}\rm{Na}^{40}\rm{K}$ with ${}^{40}\rm{K}$, New Journal of Physics 23, 115010 (2021).
  * Son  _et al._ [2022] H. Son, J. J. Park, Y.-K. Lu, A. O. Jamison, T. Karman, and W. Ketterle, Control of reactive collisions by quantum interference, Science 375, 1006 (2022).
  * Tomza  _et al._ [2013] M. Tomza, K. W. Madison, R. Moszynski, and R. V. Krems, Chemical reactions of ultracold alkali-metal dimers in the lowest-energy ${}^{3}{\Sigma}$ state, Physical Review A 88, 050701 (2013).
  * Son  _et al._ [2020] H. Son, J. J. Park, W. Ketterle, and A. O. Jamison, Collisional cooling of ultracold molecules, Nature 580, 197 (2020).
  * Park  _et al._ [2023] J. J. Park, Y.-K. Lu, A. O. Jamison, T. V. Tscherbul, and W. Ketterle, A feshbach resonance in collisions between triplet ground-state molecules, Nature 614, 54 (2023).
  * Hermsmeier  _et al._ [2021] R. Hermsmeier, J. Kłos, S. Kotochigova, and T. V. Tscherbul, Quantum spin state selectivity and magnetic tuning of ultracold chemical reactions of triplet alkali-metal dimers with alkali-metal atoms, Physical review letters 127, 103402 (2021).
  * Rvachov  _et al._ [2017] T. M. Rvachov, H. Son, A. T. Sommer, S. Ebadi, J. J. Park, M. W. Zwierlein, W. Ketterle, and A. O. Jamison, Long-lived ultracold molecules with electric and magnetic dipole moments, Phys. Rev. Lett. 119, 143001 (2017).
  * [24] T. Karman, M. Gronowski, M. Tomza, J. J. Park, H. Son, Y.-K. Lu, A. O. Jamison, and W. Ketterle,  _Ab initio_ calculation of the spectrum of feshbach resonances in NaLi $+$ Na collisions, Accompanying paper that is under preparation for submission.
  * Wigner [1993] E. P. Wigner, On a class of analytic functions from the quantum theory of collisions, in  _The Collected Works of Eugene Paul Wigner_ (Springer, 1993) pp. 409–440.
  * Dyson [1962] F. J. Dyson, Statistical theory of the energy levels of complex systems. i, Journal of Mathematical Physics 3, 140 (1962).
  * Knoop  _et al._ [2011] S. Knoop, T. Schuster, R. Scelle, A. Trautmann, J. Appmeier, M. Oberthaler, E. Tiesinga, and E. Tiemann, Feshbach spectroscopy and analysis of the interaction potentials of ultracold sodium, Physical Review A 83, 042704 (2011).
  * Steinke  _et al._ [2012] M. Steinke, H. Knöckel, and E. Tiemann, $({X})1^{1}{\Sigma}^{+}$ state of LiNa studied by fourier-transform spectroscopy, Physical Review A 85, 042720 (2012).
  * Rvachov  _et al._ [2018] T. M. Rvachov, H. Son, J. J. Park, S. Ebadi, M. W. Zwierlein, W. Ketterle, and A. O. Jamison, Two-photon spectroscopy of the nali triplet ground state, Physical Chemistry Chemical Physics 20, 4739 (2018).
  * Frye and Hutson [2021] M. D. Frye and J. M. Hutson, Complexes formed in collisions between ultracold alkali-metal diatomic molecules and atoms, New Journal of Physics 23, 125008 (2021).
  * Frye and Hutson [2022] M. D. Frye and J. M. Hutson, Long-range states and feshbach resonances in collisions between ultracold alkali-metal diatomic molecules and atoms, arXiv preprint arXiv:2212.08030 (2022).
  * Park  _et al._ [2022] J. J. Park, Y.-K. Lu, A. O. Jamison, and W. Ketterle, Magnetic trapping of ultracold molecules at high density, arXiv preprint arXiv:2211.11120 (2022).
  * Julienne  _et al._ [2011] P. S. Julienne, T. M. Hanna, and Z. Idziaszek, Universal ultracold collision rates for polar molecules of two alkali-metal atoms, Phys. Chem. Chem. Phys. 13, 19114 (2011).
  * Brody [1973] T. Brody, A statistical measure for the repulsion of energy levels, Lettere al Nuovo Cimento (1971-1985) 7, 482 (1973).
  * Brody  _et al._ [1981] T. A. Brody, J. Flores, J. B. French, P. Mello, A. Pandey, and S. S. Wong, Random-matrix physics: spectrum and strength fluctuations, Reviews of Modern Physics 53, 385 (1981).
  * Porter and Thomas [1956] C. E. Porter and R. G. Thomas, Fluctuations of nuclear reaction widths, Physical Review 104, 483 (1956).