text:
\section{Introduction}

Despite some interest in clustering of normally distributed
data \cite{Nielsen2009}, most studies focus on models for
clustering   such as mixture models or classification-tree based
methods \cite{Wasserman2010,Bishop2006}, and few consider
the (usually unknown) distribution of the underlying data.
More recently the distribution of data has been considered through
clustering non-normal distributions or assuming distributions with
outliers (for example see
Refs \cite{Punzo2020,Dotto2019,Bagnato2017,Lin2014,Lee2013,Contreras2012}).
However, multivariate normal distributions commonly arise,
in particular describing the distribution of maximum
likelihood estimates (MLEs) of parameterised models.
As discussed below, there are many advantages to
fitting a parameterised model to data and clustering the
fitted parameters - the original intended application of this work.

Section \ref{sectionFit} outlines the advantages of fitting
a parameterised model prior to clustering, an approach that
presently is rarely used \cite{Kim2008,Kim2017,Park2019}.
Section \ref{sectionCluster} will describe a simple model
for clustering multivariate normal distributions.
It assumes that clusters have items with the same
mean, and that the likelihood is determined by the independent
likelihood of membership of each individual cluster.
No assumption is made for the distribution of clusters.
The exactly calculated log-likelihood is equivalent to an
exact calculation of the Bayesian Information
Criterion (BIC) \cite{Wasserman2010,Schwarz1978,Djuric1998,Cavanaugh1999}.
It includes a weighted sum of squares term that penalises poor fits,
and a term to capture model complexity that has similarities to the equivalent
model-complexity terms in
AIC \cite{Wasserman2010,Akaike1973}
and BIC \cite{Wasserman2010,Schwarz1978,Djuric1998,Cavanaugh1999}.
For recent progress on approximate calculations of BIC and a review of
related literature, see \cite{Teklehaymanot2018}.
Section \ref{prior} generalises the calculation in
Sections \ref{sectionCluster} and \ref{secLL} with a flat prior, to
a normally distributed prior with zero-mean but arbitrary covariance.
Section \ref{limits} considers some important limits.
Section \ref{example} provides numerical examples.
Section \ref{sectionDiscuss} discusses
the results and highlights topics for future work.

\section{Clustering of parametric models}\label{sectionFit}

Most parametric models use maximum likelihood estimates (MLEs) to
obtain normally distributed estimates of their parameters,
specified by a mean $\mu$ and a
covariance matrix $\Sigma$ \cite{Wasserman2010}.
This suggests an alternative to directly clustering labelled data, that
involves firstly fitting a model to capture relevant information, and then
clustering the normally distributed estimates.
For example, we could fit a survival model for the presence of disease as
a function of known risk-factors and confounders associated with disease
risk, and then cluster the normally distributed estimates for risk-factor
associations.
There are several advantages to this approach:
\begin{enumerate}

	\item{The distribution of underlying data can be unknown, but
		under reasonable regularity conditions, MLE estimates
		are normally distributed \cite{Wasserman2010,Hardle2015}.}

    \item{Estimates can be stratified and multiply-adjusted
     for known confounders, helping to extract the information of interest
     from potentially noisy
     data \cite{Wasserman2010,Collett2014}.}

    \item{The fitted models can be more interpretable and familiar
    to the scientific community. For example,
    proportional hazards models are commonly used by medical researchers.}

    \item{Marginalisation \cite{Wasserman2010,Bishop2006} can be used
    to cluster by
    parameter subsets of greater interest, e.g. risk factors as
    opposed to confounders, or minimum versus maximum quantiles.}

    \item{By fitting a model, we can build-in
     prior knowledge through the model.}

\end{enumerate}
These benefits are becoming recognised \cite{Kim2008,Kim2017,Park2019},
with a similar approach being used to detect changes in gene expression by
clustering Fourier series coefficients \cite{Kim2008,Kim2017}.
Clustering parameters of linear-models such as a Fourier series, are examples of
clustering the normally-distributed MLE estimates of parameterised models.
Here we consider the general problem of clustering multivariate normals,
to determine both the number and membership of clusters.

\section{Clustering multivariate normal distributions}\label{sectionCluster}

The derivation below is equivalent to an exact calculation of
the Bayesian information criterion for the model, where the
normally distributed MLEs for the mean and covariance
$\{ (\hat{\mu}_i,\hat{\Sigma}_i) \}$ are the data, and
the model will be the proposed clusters with an equal mean
within each cluster.
We are interested in the likelihood of a
specific partition with clusters labeled by $f(i)$.
The model assumes,
\begin{equation}\label{muEq}
\hat{\mu}_i \sim N \left( \mu_{f(i)}, \hat{\Sigma}_i \right)
\end{equation}
with $\mu_{f(i)}=\mu_{f(j)}$ iff $f(i)=f(j)$.
Writing $x_i = \hat{\mu}_i$,
$\Gamma_i = \hat{\Sigma}_i^{-1}$,
$x=\{x_i\}$, $\Gamma=\{\Gamma_i\}$, and
using Bayes theorem \cite{Wasserman2010,Bishop2006},
the probability of a specific partition
$G=\{G_g\}$ with means $\mu=\{\mu_g\}$ and $g=1..m$, is,
\begin{equation}\label{bayes1}
\begin{array}{ll}
P \left( G | x, \Gamma \right)
&\propto
P \left( x | G, \Gamma \right)
P\left( G  | \Gamma \right)
\\
&=
\int_{-\infty}^{\infty} d\mu_1 ... \int_{-\infty}^{\infty} d\mu_m
P \left( x, \mu |
G, \Gamma \right)
P\left( G | \Gamma \right)
\\
&=
\int_{-\infty}^{\infty} d\mu_1 ... \int_{-\infty}^{\infty} d\mu_m
P \left( x | \mu,
G, \Gamma \right)
P \left( \mu
| G, \Gamma \right)
P\left( G | \Gamma \right)
\end{array}
\end{equation}
For independent normally distributed $x_i$ with covariance
$\Gamma_i^{-1}$, from cluster $g$ with mean $\mu_g$,
\begin{equation}
P\left( x | \mu, G, \Gamma \right)
= \Pi_{g=1}^m \Pi_{i\in G_g}
P \left( x_i | \mu_g, \Gamma_i \right)
\end{equation}
with,
\begin{equation}\label{toFact}
P\left( x_i | \mu_g, \Gamma_i \right)
=
\frac{1}{\sqrt{(2\pi)^p |\Gamma_i^{-1}|}}
\exp \left( - \frac{1}{2} (x_i-\mu_g)^T
\Gamma_i (x_i-\mu_g) \right)
\end{equation}
Assuming a model with
$P( \mu | G, \Gamma )=\Pi_{g=1}^m P(\mu_g)$
and $P(G|\Gamma)=P(G)$,
then using Eqs. \ref{bayes1}-\ref{toFact} and factorising terms,
\begin{equation}\label{intmuEq}
\begin{array}{c}
P\left( G| x, \Gamma \right) \propto
P(G)
\left( \Pi_{g=1}^m \Pi_{i\in G_g} \frac{1}{\sqrt{(2\pi)^p|\Gamma_i^{-1}|}} \right)\times
\\    \Pi_{g=1}^m \int_{-\infty}^{\infty} d\mu_g
P(\mu_g)
\exp \left(
-\frac{1}{2}  \sum_{i\in G_g}
(x_i- \mu_g)^T \Gamma_i (x_i-\mu_g) \right)
\end{array}
\end{equation}

In principle, the choice of $P(G)$ depends on the application.
There could be applications where
$P(G)$ should be proportional to the number of partitions with
$m$ clusters and $n_g$ members in each cluster (the Stirling number of
the 2nd kind).
This would be analogous to the binomial distribution, where only the total
numbers in each of two possible clusters are known.
However, the members are all identifiable, similarly to a Bernoulli distribution.
Here $P(G)$ is regarded as constant, but the number of partitions may
need to be considered by MCMC schemes to generate correctly
distributed samples, or in physical applications where noise influences
cluster formation.

\section{Evaluating the likelihood - a uniform prior}\label{secLL}

The (prior) distribution $P(\mu_g)$ will usually be unknown.
The following Section \ref{prior} shows how
a normal distribution for $P(\mu_g)$ with
$\mu_g\sim N(0,\Gamma_0^{-1})$, can be incorporated into the analysis.
Section \ref{limits} shows that a
``uniform prior'' cannot in general be regarded as a limiting case
of the normal prior $N(0,\Gamma_0^{-1})$,
for the example with
$\Gamma_0=I/\sigma^2$ where $I$ is the identity matrix, in the limit where
$\sigma^2 \rightarrow \infty$.
Here for simplicity of presentation, we firstly consider a uniform prior.

To integrate over $\mu$ write $\Gamma_i=\Sigma_i^{-1}$ and note
that $\Sigma_i$ and their inverses $\Gamma_i$ are symmetric, and use this to write,
\begin{equation}\label{sseq}
\begin{array}{l}
    \sum_i (x_i-\mu)^T  \Gamma_i (x_i-\mu)
    \\ =
    \left( \mu - \left( \sum_i \Gamma_i \right)^{-1} \sum_i \Gamma_i x_i \right)^T
    \left( \sum_i \Gamma_i \right)
    \left( \mu - \left( \sum_i \Gamma_i \right)^{-1} \sum_i \Gamma_i x_i \right)
    \\+ \sum_i x_i^T \Gamma_i x_i
    - \left( \sum_i x_i^T \Gamma_i \right)
    \left( \sum_i \Gamma_i \right)^{-1}
    \left( \sum_i \Gamma_i x_i \right)
\end{array}
\end{equation}
The terms involving $\mu$ in \ref{sseq} will factorise in
Eq. \ref{intmuEq}, and lead to Gaussian integrals that
integrate to give functions of $\{\Gamma_i\}$
that are independent of $\{ x_i \}$.
The remaining terms will be,
\begin{equation}\label{remEq}
\begin{array}{l}
    \sum_i x_i^T \Gamma_i x_i
    - \left( \sum_i x_i^T \Gamma_i \right)
    \left( \sum_k \Gamma_k \right)^{-1}
    \left( \sum_j \Gamma_j x_j \right)
    \\=
    \sum_{i,j} x_i^T
    \Gamma_i \left( \sum_k \Gamma_k \right)^{-1} \Gamma_j
    x_i
    - \sum_{i,j} x_i^T
    \Gamma_i \left( \sum_k \Gamma_k \right)^{-1} \Gamma_j
    x_j
\end{array}
\end{equation}
Because $\Sigma_i$ and their inverses $\Gamma_i$ are symmetric, then $C_{ij}=\Gamma_i \left( \sum_k \Gamma_k \right)^{-1}\Gamma_j$
has $C_{ij}=C_{ji}^T$, as can be seen by taking the transpose of $C_{ij}$.
Using $C_{ij}=C_{ji}^T$, $a^T b=b^T a$ for vectors $a$ and $b$,  and relabeling the indices $i$ and $j$,
\begin{equation}\label{Cij}
    \begin{array}{l}
    \frac{1}{2} \sum_{i,j} (x_i-x_j)^T C_{ij} (x_i-x_j)
    \\=
    \frac{1}{2} \sum_{i,j} x_i^T C_{ij} x_i
    +\frac{1}{2} \sum_{i,j} x_j^T C_{ij} x_j
    -\frac{1}{2} \sum_{i,j} x_i^T C_{ij} x_j
    -\frac{1}{2} \sum_{i,j} x_j^T C_{ij} x_i
    \\=
    \sum_{i,j} x_i^T C_{ij} x_i
    -\frac{1}{2} \sum_{i,j} x_i^T C_{ij} x_j
    -\frac{1}{2} \sum_{i,j} x_j^T C_{ji}^T x_i
    \\=
    \sum_{i,j} x_i^T C_{ij} x_i - x_i^T C_{ij} x_j
    \end{array}
\end{equation}
Hence using Eqs. \ref{remEq} and \ref{Cij} we have,
\begin{equation}\label{ssEqf}
    \begin{array}{l}
    \sum_i x_i^T \Gamma_i x_i
    - \left( \sum_i x_i^T \Gamma_i \right)
    \left( \sum_k \Gamma_k \right)^{-1}
    \left( \sum_j \Gamma_j x_j \right)
    \\=
    \frac{1}{2} \sum_{i,j} (x_i-x_j)^T
    \Gamma_i \left( \sum_k \Gamma_k \right)^{-1}  \Gamma_j
    (x_i-x_j)
   \end{array}
\end{equation}
where the sums over $i$, $j$, and $k$ will range over elements in cluster $g$.

Eq. \ref{ssEqf} is an intuitively reasonable result, that after
marginalisation with respect to the means $\mu_g$, the likelihood
will be determined from differences in the estimated means within
each cluster with an inverse covariance matrix
$\Gamma_i (\sum_{k \in G_g} \Gamma_k )^{-1} \Gamma_j$, that is weighted
by the inverse covariances of $x_i$ and $x_j$, and the inverse of the
average-inverse-covariance of their cluster.
Note that the extra factor of $1/2$ prevents double counting in the sum
over $i$ and $j$, with $(1/2)\sum_{i,j} F(i,j)=\sum_i\sum_{j\geq i} F(i,j)$
for any  $F(i,j)$ with $F(i,i)=0$ and is symmetric with respect to exchange of
$i$ and $j$, so that the sum involves exactly one term for
every unique combination of $i$ and $j$ in cluster $g$.
Using Eqs. \ref{sseq} and \ref{ssEqf}, and integrating over $\mu_g$ gives,
\begin{equation}\label{gmuEq}
\begin{array}{c}
    \int_{-\infty}^{\infty} d\mu_g
    \exp \left(
    -\frac{1}{2}  \sum_{i\in G_g}
    (x_i- \mu_g)^T \Gamma_i (x_i-\mu_g) \right)
    \\=
    \left| 2\pi \left( \sum_{k \in G_g} \Gamma_k \right)^{-1} \right|^{1/2} \times
    \\
    \exp \left( - \frac{1}{4}
    \sum_{i,j \in G_g}
    (x_i- x_j)^T \Gamma_i
    \left( \sum_{k\in G_g} \Gamma_k \right)^{-1}
    \Gamma_j (x_i-x_j) \right)
\end{array}
\end{equation}
Using Eq. \ref{intmuEq} and \ref{gmuEq},  the log-likelihood has,
\begin{equation}\label{logLik}
    \begin{array}{ll}
    \log\left( P\left( G | x, \Gamma  \right)  \right) + C =
    &
    - \frac{1}{2} \sum_i \log  \left| 2\pi \Gamma_i^{-1} \right|
    \\&
    + \frac{1}{2} \sum_{g=1}^m
    \log \left|2\pi
    \left(\sum_{k\in G_g} \Gamma_k \right)^{-1}
    \right|
    \\&- \frac{1}{4}
    \sum_{g=1}^m
    \sum_{i,j \in G_g}
    (x_i- x_j)^T \Gamma_i \left( \sum_{k\in G_g} \Gamma_k \right)^{-1}  \Gamma_j (x_i-x_j)
\end{array}
\end{equation}
where $C$ is a constant that ensures correct normalisation of the probability
distribution.

The first term on the right side of Eq. \ref{logLik} is
independent of the clustering model, the final term involves a sum of
squares that measures the goodness
of fit of the data to the model, and the middle term captures
the model's complexity.
For $\mu$ with dimension $p$, the determinants in the middle
term can be considered as a product of $p$ eigenvalues, leading
to a sum of $p$ logarithms of the eigenvalues for each cluster.
Those $p\times m$ terms are analogous to the number of free
parameters in the AIC, $m$ clusters with $p$ dimensions give
$m\times p$ parameters.
An estimated covariance is roughly proportional to the number
of data used to estimate it \cite{Wasserman2010,Hardle2015}, which
brings a slightly more
complicated dependence on the data into Eq. \ref{logLik},
that has some similarities to the $\log(n)$ term in BIC for
the number of data points.

Note that the sum of squares term involves
$\left( \sum_{k\in G_g} \Gamma_k \right)$, that
uses information from all members of the cluster.
As a result, clusters that minimise Eq. \ref{logLik} will
in general differ from clusters that use pair-wise
distance based measures.
Therefore clusters that minimise Eq. \ref{logLik} will
typically differ from pair-wise distance-based e.g.
hierarchical clustering models.

We can use Eq. \ref{logLik} to evaluate and compare the log-likelihoods
of clusterings proposed for example by heirarchical clustering methods.
Alternately we can directly maximise Eq. \ref{logLik}
to obtain an MLE.
Because the likelihood factorises in terms of the clusters,
if Metropolis MCMC were used to generate a sample of clusters,
then the change in log-likelihood at each iteration is determined
solely by the change in log-likelihood of the clusters whose membership changes.

\section{A normal prior}\label{prior}

As mentioned earlier, if we had taken a normally distributed
prior for $P(\mu_g)$ with $\mu_g\sim N(0,\Gamma_0^{-1})$ then
there would be an extra term $\mu^T \Gamma_0 \mu$ on the
left side of Eq. \ref{sseq}, causing the $\sum_{i} \Gamma_i$
terms on the right side to be replaced by
$\Gamma_0+\sum_{i} \Gamma_i$, with,
\begin{equation}
\begin{array}{l}
\sum_i (x_i-\mu)^T  \Gamma_i (x_i-\mu) + \mu^T \Gamma_0 \mu
\\ =
\left( \mu - \left( \Gamma_0 + \sum_i \Gamma_i \right)^{-1} \sum_i \Gamma_i x_i \right)^T
\left( \Gamma_0 + \sum_i \Gamma_i \right)
\left( \mu - \left( \Gamma_0 + \sum_i \Gamma_i \right)^{-1} \sum_i \Gamma_i x_i \right)
\\+ \sum_i x_i^T \Gamma_i x_i
- \left( \sum_i x_i^T \Gamma_i \right)
\left( \Gamma_0 + \sum_i \Gamma_i \right)^{-1}
\left( \sum_i \Gamma_i x_i \right)
\end{array}
\end{equation}
Noting that,
\begin{equation}
\begin{array}{rl}
\sum_i x_i^T \Gamma_i x_i =
&
\sum_{i} x_i^T \Gamma_i
\left( \Gamma_0 + \sum_k \Gamma_k \right)^{-1}
\left( \Gamma_0 + \sum_j \Gamma_j \right)
x_i
\\
=&
\sum_{i} x_i^T \Gamma_i
\left( \Gamma_0 + \sum_k \Gamma_k \right)^{-1}
\Gamma_0 x_i
\\+&
\sum_{i,j} x_i^T \Gamma_i
\left( \Gamma_0 + \sum_k \Gamma_k \right)^{-1}
\Gamma_j x_i
\end{array}
\end{equation}
and using this with Eq. \ref{Cij}, with
$C_{ij}=\Gamma_i \left( \Gamma_0 + \sum_k \Gamma_k \right)^{-1} \Gamma_j$, we get,
\begin{equation}\label{Eqrhs}
\begin{array}{l}
\sum_i x_i^T \Gamma_i x_i
- \left( \sum_i x_i^T \Gamma_i \right)
\left( \Gamma_0 + \sum_k \Gamma_k \right)^{-1}
\left( \sum_j \Gamma_j x_j \right)
\\=
\sum_{i} x_i^T \Gamma_i
\left( \Gamma_0 +  \sum_k \Gamma_k \right)^{-1}
\Gamma_0 x_i
\\
+
\frac{1}{2} \sum_{i,j} (x_i-x_j)^T
\Gamma_i \left( \Gamma_0 +  \sum_k \Gamma_k \right)^{-1}  \Gamma_j
(x_i-x_j)
\end{array}
\end{equation}
The first term on the right of Eq. \ref{Eqrhs} can alternately be written in a symmetrical form with,
\begin{equation}
\begin{array}{l}
\sum_i x_i^T \Gamma_i x_i
- \left( \sum_i x_i^T \Gamma_i \right)
\left( \Gamma_0 + \sum_k \Gamma_k \right)^{-1}
\left( \sum_j \Gamma_j x_j \right)
\\=
\frac{1}{2}
\sum_{i} x_i^T \left( \Gamma_i + \Gamma_0 \right)^T
\left( \Gamma_0 +  \sum_k \Gamma_k \right)^{-1}
\left( \Gamma_0 + \Gamma_i \right) x_i
\\
-
\frac{1}{2}
\sum_{i} x_i^T \Gamma_i^T
\left( \Gamma_0 +  \sum_k \Gamma_k \right)^{-1}
\Gamma_i x_i
-
\frac{1}{2}
\sum_{i} x_i^T \Gamma_0^T
\left( \Gamma_0 +  \sum_k \Gamma_k \right)^{-1}
\Gamma_0 x_i
\\
+
\frac{1}{2} \sum_{i,j} (x_i-x_j)^T
\Gamma_i \left( \Gamma_0 +  \sum_k \Gamma_k \right)^{-1}  \Gamma_j
(x_i-x_j)
\end{array}
\end{equation}
Recalling that,
\begin{equation}
P(\mu_g)=\frac{1}{ \sqrt{\left| 2\pi \Gamma_0^{-1} \right| } }
\exp \left\{ -\frac{1}{2} \mu_g^T \Gamma_0 \mu_g \right\}
\end{equation}
and that there is one for each of the $m$ clusters,
then after marginalisation over each $\mu_g$,
the log-likelihood will have,
\begin{equation}\label{logLik2}
\begin{array}{ll}
\log\left( P\left( G | x, \Gamma  \right)  \right) + C =
&
- \frac{m}{2}        \log  \left| 2\pi \Gamma_0^{-1} \right|
- \frac{1}{2} \sum_i \log  \left| 2\pi \Gamma_i^{-1} \right|
\\
&+ \frac{1}{2} \sum_{g=1}^m
\log \left|2\pi
\left( \Gamma_0 + \sum_{k\in G_g} \Gamma_k \right)^{-1}
\right|
\\&- \frac{1}{4}
\sum_{g=1}^m
\sum_{i,j \in G_g}
(x_i- x_j)^T \Gamma_i \left( \Gamma_0 + \sum_{k\in G_g} \Gamma_k \right)^{-1}  \Gamma_j (x_i-x_j)
\\&-
\frac{1}{2}
\sum_{g=1}^m
\sum_{i\in G_g}
x_i^T \Gamma_i \left( \Gamma_0 + \sum_{k\in G_g} \Gamma_k \right)^{-1}  \Gamma_0 x_i
\end{array}
\end{equation}
where $C$ is a constant that ensures correct normalisation of the probability
distribution.

\section{Important limits}\label{limits}

A flat prior cannot be considered as a limiting case of a normal prior with,
for example, $\Gamma_0^{-1} = \sigma_0^2I$ and $\sigma_0^2 \rightarrow \infty$.
This can be seen from Eqs. \ref{intmuEq} and Eq. \ref{logLik2}.
Considering Eq. \ref{logLik2},
the first term is $-(m/2)p\log(2\pi \sigma_0^2)$ and tends to $-\infty$, the
2nd term is independent of $\Gamma_0$, the 3rd term tends to
$\log ( 2\pi (\sum_{k\in G_g} \Gamma_k){^{-1}})$,
the 4th term has
$( \Gamma_0 + \sum_{k\in G_g} \Gamma_k )^{-1}
\rightarrow ( \sum_{k\in G_g} \Gamma_k)^{-1}$,
and the final term is proportional to $\Gamma_0=I/\sigma_0^2$ and tends to zero.
The behaviour of the 1st term that normalises the factors of $P(\mu_g)$ can be
understood from Eq. \ref{intmuEq}.
As $\sigma_0^2$ becomes larger, $P(\mu_g=0)$ must become increasingly small to ensure
that $P(\mu_g)$ is correctly normalised, and there is one factor
of $P(\mu_g)$ per cluster.
For this example with $\Gamma_0^{-1}=\sigma_0^2 I$, the first term provides a penalty that
is proportional to the number of free parameters $m\times p$.
For the alternative limit with $\sigma_0^2 \rightarrow 0$, then the first and third
terms cancel, but the final term diverges at a rate proportional to $1/\sigma_0^2$.

Often data without covariances are clustered.
To explore this limit,
take $\Gamma_i^{-1} = \sigma^2 I$ and let $\sigma \rightarrow 0$.
The final term diverges most rapidly ($\sim 1/\sigma^2$), and has,
\begin{equation}
(x_i- x_j)^T
\Gamma_i
\left( \Gamma_0 + \sum_{k\in G_g} \Gamma_k \right)^{-1}
\Gamma_j (x_i-x_j)
\\
\rightarrow \frac{1}{\sigma^2} \frac{(x_i-x_j)^T(x_i-x_j)}{n_g}
\end{equation}
where $n_g$ is the number of items in group $g$.
This can be written as,
\begin{equation}
- \frac{1}{4}
\sum_{g=1}^m
\sum_{i,j \in G_g}
(x_i- x_j)^T \Gamma_i \left( \Gamma_0 + \sum_{k\in G_g} \Gamma_k \right)^{-1}  \Gamma_j (x_i-x_j)
\\
\rightarrow
-\frac{1}{2\sigma^2}
\sum_{g=1}^m
\frac{1}{2n_g}
\sum_{i,j \in G_g}
(x_i-x_j)^T(x_i-x_j)
\end{equation}
showing that as the variance in the data becomes increasingly small,
the log-likelihood is minimised by the minimum sum of within-group
sum of squares.

\section{Example}\label{example}

To illustrate the differences between minimising Eqs. \ref{logLik}, \ref{logLik2}, and
traditional statistical tests, hierarchical clustering was used with the
ward D2 algorithm and the Battacharyya distance to generate clusters, and the
optimum number of clusters was compared for the different approaches.
For a more traditional statistical test,
consider the null hypothesis of equal means in each group.
Data from the same group $g$ have $x_i \sim N(\mu_g,\Gamma_i^{-1})$, so
summing over all pairs in all groups we have,
\begin{equation}
\frac{1}{2} \sum_{g=1}^m \sum_{i,j\in G_g}
\left( x_i - x_j \right)^t \left( \Gamma_i^{-1} + \Gamma_j^{-1} \right)^{-1}
\left( x_i - x_j \right)
\sim \chi^2_q
\end{equation}
where $q=\sum_{g=1}^m n_g(n_g-1)/2$ and $n_g$ is the number of data points in
cluster $g$, so that $n_g(n_g-1)/2$
is the number of distinct pairs of different diseases
in each cluster and $q$ is the total sum of disease pairs.
This provides a statistical test for the null hypothesis of equal means
in each group, that we can compare with the clustering results from minimising
Eqs. \ref{logLik} or \ref{logLik2}.

Data $\{x_i\}$ were simulated by firstly sampling
cluster means $\mu_g \sim N(0,S_0)$ and covariance
matrices $\{S_g\}$ (see Appendix \ref{theApp} for details),
then sampling $x_i \sim N(\mu_g, S_g)$.
The $\{x_i, S_g\}$ could represent fitted MLEs for example.
Sampled means $\mu_g$ were omitted unless they had a
statistically significant difference from 0 after an FDR multiple testing
adjustment of a multivariate $\chi^2$ test, using the covariance $S_g$
of a proposed group.
Noise was added to $S_g$, to give $S_i$ prior to clustering,
allowing tests of the sensitivity of results to noise that may exist
in fitted parameters.
Alternately, we could have added noise to the $S_g$ prior to sampling $\{x_i\}$, to
explore the sensitivity of results to underlying differences in the covariances.
Examples are in figure \ref{fig1}, with details in Appendix \ref{theApp}.

The simulations considered a flat prior and a normal prior with  $\sigma_0^2 = 1$.
Data with 10 parameters were sampled, giving a prior estimate of $E[\mu_g^2] = 10$.
Examples were considered for a high signal with $E[\mu_g^2]\simeq 20$, and a
lower signal with $E[\mu_g^2]\simeq 5$.
Clustering tended to fall into two types depending on the
signal to noise ratio, as characterised by
$E[\mu_g^2]/E[x_i^2] \sim Tr(S_0)/Tr(E[S_g])$.
For a high signal to noise ratio with $E[\mu_g^2]/E[x_i^2]\simeq 8.0$,
minimising Eqs. \ref{logLik} and \ref{logLik2} gave similar results
as choosing the minimum number of groups for which $P(\chi^2_q) \geq 0.05$,
and all select a similar number of clusters to that being simulated.
In all cases the flat prior had a minimum close to
the fewest clusters with $p>0.05$, but the minimum is very shallow, which would
lead to very wide confidence intervals.
For a lower signal to noise ratio with
$E[\mu_g^2]/E[x_i^2]\simeq 2.0$, all methods
selected fewer clusters than were sampled,
but this is much more pronounced for Eq. \ref{logLik2}.
Although the minimum in Eq. \ref{logLik2} tended to be for fewer clusters
than the minimum of Eq. \ref{logLik}, both
occur near the ``elbow'' in $\chi^2_q$,
offering a less subjective form of ``elbow criterion''.

\begin{figure}[ht]
	\centering
	\vspace*{-50mm}

	\hspace*{-8mm}
	\includegraphics[width=1.1\linewidth,angle=0]{clusterPlot}

	\vspace*{-32mm}
	\caption[]{

		We clustered 300 items with 10 parameters each, from 50 clusters
		with between 2 and 14 items.

		The top row has the highest signal to noise ratio, with $E[\mu_g^2]/E[x_i^2]=8.0$ (blue), compared to $2.0$ (red and orange).

		Its minimum log-likelihood coincided with the fewest clusters
		with p$>$0.05, at 50 clusters.

		For lower signal to noise ratios, both methods underestimated the number of
		clusters.

		The bottom row (orange) added

		noise to the covariances,

		with $S_i = S_g+ (1/2) \mbox{mean}(S_g) \mbox{outer}(u,u)$, and $u\sim U(0,1)$,
		further increasing the p-value based underestimation of
		the number of clusters.

		For all cases, the minimum log-likelihoods remain near the elbow in $\chi^2$.

	}
	\label{fig1}
\end{figure}

\section{Discussion}\label{sectionDiscuss}

The results presented in Sections \ref{secLL} and \ref{prior} are immediately
applicable to estimated data with covariances, where e.g. hierarchical
clustering needs an objective method to select the number of clusters.
They provide intuitively reasonable results when the data are noisy and statistical
tests are of limited use, providing a similar, but less subjective alternative to
the elbow criterion.

The existence of a log-likelihood for clustering opens a range
of fascinating problems.
Firstly, the log-likelihood can be maximised to determine the optimum
cluster, or clusters, and coded algorithms are needed to do this.
Because the MLE will change by discrete amounts when elements are
placed in different clusters, gradient-based methods to find MLEs
may need to be modified or replaced.
More interestingly, the existence of a likelihood for both the number and membership
of clusters allows the possibility of determining a confidence set for the maximum  likelihood estimate, although it may take further thought about how to use and interpret it.
A first step will be to relabel clusters to identify equivalent
clusterings (obtained by permuting the order of the clusters),
e.g. we can label items $1$ to $n$, and then
order within each cluster from smallest to largest, before
ordering clusters by the smallest item number in each cluster to produce a unique
description of equivalent clusterings.
Confidence sets will need to be characterised to describe the
similarity between clusters, such as the number of clustered pairs
common to all clusterings.
Although it is unclear how best to characterise a confidence set,
it should be relatively easy to test whether a particular clustering's
log-likelihood lies within an MLE's 95\% confidence interval.

Looking ahead, the model might be extended to clusters whose members are
assumed to have both the same mean and covariance.
In the present model the covariance matrices are solely used to
assess the likelihood of having the same mean.
It might be possible to evaluate the BIC exactly for this, or other situations,
and to explore the accuracy of AIC and BIC approximations.

\section{Acknowledgements}
Anthony Webster is supported by a fellowship from the
Nuffield Department of Population Health, University of Oxford, UK.

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}
The Efimov effect, as predicted in 1970\;\cite{Efimov1970,Efimov1973}, has attracted broad interest in atomic and nuclear physics\;\cite{Kraemer2006,Thomas2010,Berninger2011,Huang2014,Zaccanti2009,Bloom2013,Hu2014,Maier2015,Kato2017,Pires2014,Tung2014Dec,HUlmanis2016,hafner2017role,
Barontini2009,blume2012few,Greene2017universal,Jensen2004structure,Kolganova2011,gong2019compensating,yao2018combined}. When the two-body s-wave scattering length $a$ is tuned to a large value compared to the characteristic range $r_{0}$ of the two-body interaction potential, an infinite series of trimer states are formed. These Efimov trimer states follow a discrete scale invariance: $E_{n} = \lambda^{2}E_{n+1}$, where $n$ denotes the $n$th Efimov state with three pairs of interactions in the resonant limit. $\lambda = e^{\pi/s_{0}}$ is the scaling constant, where $s_0$ depends on the number of resonant interactions as well as the quantum statistics and mass ratio of a trimer's constituents\;\cite{Eric2006,Naidon2017,IncaoMass2006,Incao2005}.

In ultracold gases, where two-body scattering length $a$ can be tuned to any value using Feshbach resonances, the Efimov effects can be observed by measuring the three-body event rate constant, which is resonantly enhanced when an Efimov trimer state approaches zero binding energy at a negative $a$ denoted $a_{\scriptscriptstyle-}$. On the other side of the Feshbach resonance, interference effects
lead to minima in the three-body rate coefficient at a positive $a$ denoted $a_{\scriptscriptstyle+}$ and the scattering length where an Efimov trimer state degenerates with the atom-molecule threshold denoted $a_{\scriptscriptstyle*}$\;\cite{Naidon2017,D_Incao_2018}. These Efimov features have been observed in experiments with ultracold homonuclear Bose, three-component Fermi, and heteronuclear Bose-Fermi gases\;\cite{Kraemer2006,Nakajima2010,Lompe2010,Pollack2009,Knoop2009,Floerchinger2009,Zaccanti2009,Bloom2013,Hu2014,Maier2015,Kato2017,Pires2014,Tung2014Dec,HUlmanis2016,hafner2017role}. The universal relationship
between these features for the case of three identical bosons was first studied by Braaten and Hammer in Ref.\;\cite{Eric2006} and later expanded upon in Ref.\;\cite{Gogolin2008Apr}. Therefore, with these universal relationships, the ground Efimov state energy $E_{0}$ and the first Efimov resonance position can fully determine the Efimov spectrum; these factors are usually considered as the three-body parameters (3BP). Recent measurements of Efimov resonances for a number of ultracold atom species have revealed an unexpected universality, in which the first Efimov resonance $a^{\scriptscriptstyle(1)}_{\scriptscriptstyle-}$ is determined by the van der Waals (vdW) length $r_{\scriptscriptstyle\textsl{vdW}}$: $a^{\scriptscriptstyle(1)}_{\scriptscriptstyle-}=-(8.9 \pm 1.8)r_{\scriptscriptstyle\textsl{vdW}}$\;\cite{D_Incao_2018,Thomas2010,Berninger2011,Huang2014,Ottenstein2008,Scott2009,Chin2010}, where $r_{\scriptscriptstyle\textsl{vdW}}=(2\mu_{2b}C_{6})^{1/4}/2$, which characterizes the length scale between two neutral atoms interacting via the van der Waals interaction $-C_{6}/r^{6}$ with two-body reduced mass $\mu_{2b}$. The calculations from theories based on a single-channel van der Waals model confirmed this vdW universality for broad Feshbach resonances\;\cite{WangJia2012,WangYujun2012,NaidonMar2014,NaidonAug2014,Huang2014Oct}. Similarly, universality of the 3BP has been proposed for heteronuclear atomic systems such as BBX that consist of two heavy atoms B resonantly interacting with a light one X, with an additional dependence on the intraspecies interaction between the atoms B\;\cite{WangYujun2012}.

There are several experimental groups working on Efimov unfavored systems of K-Rb admixtures\;\cite{Barontini2009,Bloom2013,Hu2014,Wacker2016KRb}. The initial experiment in $^{41}$K-$^{87}$Rb Bose-Bose mixtures reported an unexpected nonuniversal behavior of Efimov resonances $a_{\scriptscriptstyle-}=-246\,a_{0}$\;\cite{Barontini2009}, where $a_{0}$ denotes Bohr radius. It thus prompted several further experiments in $^{40}$K-$^{87}$Rb Bose-Fermi mixtures, which all surprisingly showed an Efimov-type resonance around $a_{\scriptscriptstyle*}=230(30)\,a_{0}$ for collisions between atoms and dimers but not for atomic three-body recombination\;\cite{Zirbel2008,Bloom2013,Kato2017}. Moreover, the predicted universal relationships for the
location of Efimov features due to the same trimer state are $|a_{\scriptscriptstyle-}|/a_{\scriptscriptstyle*} = 1.96 $\;\cite{Helfrich2010}. Given the observed
resonance at $a_{\scriptscriptstyle*}=230(30)\,a_{0}$, these relations yield $a_{\scriptscriptstyle-}=-450\,a_{0}$, which is
excluded by the measurements mentioned above. Therefore, a breakdown of universal relationships exists between Efimov resonance and interference features in K-Rb mixtures under the zero-range theory, which assumes $a_{\scriptscriptstyle\textsl{BB}}\rightarrow 0$.

In experiments with ultracold mixtures of Cesium and Lithium atoms, the absence of ground state resonance for positive Cs-Cs scattering length is observed, which suggests that two Efimov branches may exist for $a_{\scriptscriptstyle\textsl{BB}}>0$ ($a_{\scriptscriptstyle\textsl{BB}}$ represents the B-B scattering length) in the BBX system, with different microscopic mechanisms governing the value of the 3BP\;\cite{Tung2014Dec,HUlmanis2016,hafner2017role}. The existence of two Efimov branches may
help explain the breakdown of universal relationships between resonance and interference features in K-Rb mixtures. The atom-dimer features in the RbRb + K channel would shed new light on universal few-body scenarios in mixed systems. Ref.\;\cite{Acharya2016} has derived a set of universal relations that connect the different Efimov features, including the effects of the finite interspecies effective range and the finite intraspecies scattering length, which requires two three-body observables as input parameters. For the $^{87}$Rb-$^{87}$Rb-$^{40}$K system, the $a_{\scriptscriptstyle*}=230(30)\,a_{0}$ resonance is a reliable parameter that has been observed by several experiments. Thus, data of $a_{\scriptscriptstyle+}$ or $a_{\scriptscriptstyle-}$ is important for elucidating the universal relationships between different Efimov features as well as the role of finite-range corrections. The Efimov resonance $a_{\scriptscriptstyle-}< -30000\,a_{0}$ was suggested by Wang \textit{et al.}\;\cite{WangYujun2012}, yet is far outside the range of previous measurements. It is also very difficult for theory to obtain the three-body recombination rate at such a large s-wave scattering length for the numerical stability. Thus, we focus on the positive side of the Rb-K Feshbach resonance, where interference effects lead to minima in the three-body rate coefficient at $a_{\scriptscriptstyle+}$. The prediction that the possible recombination minima at $a_{\scriptscriptstyle+}=2800\,a_{0}$\;\cite{WangYujun2012} for $^{40}$K-$^{87}$Rb-$^{87}$Rb was excluded by the late experiment\;\cite{Bloom2013}, and the analysis suggested that $a_{\scriptscriptstyle+}$ should be shifted to a higher scattering length value.

In this paper, we calculate the Efimov features in $^{87}$Rb-$^{40}$K Bose-Fermi mixtures by investigating the three-body recombination and $^{40}$K-$^{87}$Rb$^{87}$Rb elastic collision processes with $^{87}$Rb-$^{87}$Rb interaction fixed at $a_{\scriptscriptstyle\textsl{RbRb}} = 100\,a_{0}$, which is close to the experimental condition. We use the R-matrix propagation method in the hyperspherical coordinate frame based on the Lennard-Jones model potential for the interactions between atoms. The hyperradius is divided into two regions: short distances, where the slow-variable-discretization (SVD)\;\cite{Tolstikhin1996SVD} method is applied to overcome the numerical difficulties at sharp nonadiabatic avoided crossings; and large distances, where the traditional adiabatic hyperspherical method is employed to avoid the large memory and central processing unit time needed in SVD. Then, the $\underline{\mathcal{R}}$ matrix was propagated from short to large distances. The scattering properties are obtained through the $\underline{\mathcal{S}}$ matrix by matching the $\underline{\mathcal{R}}$ matrix with asymptotic functions and boundary conditions\;\cite{WangJia2011}.

The paper is organized as follows.
In Sec. II, our calculation method and all necessary formulas for calculations are presented.
In Sec. III, we discuss the results and the universal relationship between the different Efimov features.
Finally, we provide a brief summary. Atomic units are
applied throughout the paper unless stated otherwise. \\

\section{Theoretical formalism}
In this work, we calculate the collision properties of the $^{87}$Rb-$^{87}$Rb-$^{40}$K system in $J^{\Pi}=0^{+}$ symmetry. $m_{i}$(i=1,2,3) represent the masses of the three atoms $^{87}$Rb, $^{87}$Rb and $^{40}$K, respectively. We employ Delves's hyperspherical coordinates and introduce the mass-scaled Jacobi coordinates. The first Jacobi vector $\vec{\rho}_{1}$ is chosen to be the vector from atom 1 to atom 2, with reduced mass $\mu_{1}$, and the second Jacobi vector $\vec{\rho}_{2}$ goes from the diatom center of mass to the third atom, with reduced mass $\mu_{2}$. $\theta$ is the angle
between $\vec{\rho}_{1}$ and $\vec{\rho}_{2}$. The hyperradius $R$ and hyperangle $\phi$ are defined as\\
\begin{equation}
\label{1}
\mu R^{2}=\mu_{1}\rho_{1}^{2}+\mu_{2}\rho_{2}^{2} ,
\end{equation}
and\\
\begin{equation}
\label{2}
\tan\phi=\sqrt{\frac{\mu_{2}}{\mu_{1}}}\frac{\rho_{2}}{\rho_{1}},\;\; 0 \leq\phi\leq\frac{\pi}{2},
\end{equation}
respectively, where $R$ is the only coordinate with the dimension of length and represents the overall size of the three-body system. $\theta$, $\phi$, and three Euler angles $(\alpha, \beta, \gamma)$ are represented collectively by $\Omega$ $[\Omega \equiv (\theta, \phi, \alpha, \beta, \gamma)]$, which describe the rotation of the plane that contains the three particles. $\mu$ is an arbitrary scaling factor that is chosen as $\mu=\sqrt{\mu_{1}\mu_{2}}$ in our calculation.\\

In hyperspherical coordinates, the Schr$\mathrm{\ddot{o}}$dinger equation can be written after rescaling the three-body wave function $\Psi_{\upsilon'}$ as $\psi_{\upsilon'}(R;\theta,\phi)=\Psi_{\upsilon'}(R;\theta,\phi) R^{5/2} \sin\phi \cos\phi$:\\
\begin{equation}
\label{3}
\left[ -\frac{1}{2\mu}\frac{d^{2}}{dR^{2}}+\left( \frac{\Lambda^{2}-\frac{1}{4}}{2\mu R^{2}}+V(R;\theta,\phi)\right) \right] \psi_{\upsilon'}(R;\Omega)=E\psi_{\upsilon'}(R;\Omega),
\end{equation}
where $\Lambda^{2}$ is the squared ``grand angular momentum operator'', and its expression is as given in Ref.\;\cite{lin1995}. The three-body interaction potential $V(R;\theta,\phi)$ is expressed as a sum of three pairwise two-body interaction potentials $v(r_{ij})$:\\
\begin{equation}
\label{4}
V(R;\theta,\phi)=v(r_{12})+v(r_{13})+v(r_{23})\,,
\end{equation}
where $r_{ij}$ is the interparticle distance. Here, we use the Lennard-Jones potential with a vdW tail to model the interactions between two atoms \;\cite{WangJia2012}, which was proven to be an excellent model potential to explore van der Waals universality in Efimov physics\;\cite{WangYujun2012,WangJia2012,HUlmanis2016,NaidonAug2014}. The potential is expressed in the form of \\
\begin{equation}
 \label{5}
 v(r_{ij})=-\frac{C_{6,ij}}{r_{ij}^{6}}\left[1-\frac{1}{2}\left(\frac{\lambda_{ij}}{r_{ij}}\right)^{6}\right],
 \end{equation}
where the $i$ and $j$ indices in $r_{ij}$ label particles $i$ and $j$. The two-body scattering length $a_{ij}$ is changed by adjusting the $\lambda_{ij}$, and $C_{6,ij}$ is the dispersion coefficient. The values of $C_{6,RbRb}$ and $C_{6,RbK}$ used here are 4698 from Ref.\;\cite{C6RbRb2014} and 4106.5 from Ref.\;\cite{C6RbK1999}, respectively.\\

The three-body wave function $\psi_{\upsilon^{'}}$ can be expanded with the complete, orthonormal set of angular wave function $\Phi_{\nu}$ and radial wave functions $F_{\nu\upsilon'}$ as\\
\begin{equation}
\label{6}
\psi_{\upsilon'}(R;\varOmega)=\sum\limits_{\nu=0}^{\infty}F_{\nu\upsilon'}(R)\Phi_{\nu}(R;\varOmega).
\end{equation}
The adiabatic potentials $U_{\nu}(R)$ and channel functions $\Phi_{\nu}(R;\varOmega)$ at fixed $R$ can be obtained by solving the following adiabatic eigenvalue equation:\\
\begin{equation}
\label{7}
\left( \frac{\Lambda^{2}-\frac{1}{4}}{2\mu R^{2}}+V(R;\theta,\phi)\right) \Phi_{\nu}(R;\varOmega)=U_{\nu}(R) \Phi_{\nu}(R;\varOmega).
\end{equation}

Then, we solve Eq.\;(\ref{3}) using the R-matrix propagation method. First, we divide the hyperradius into $(N-1)$ intervals with a set of grid points $R_{1}<R_{2}< ...<R_{N}$ according to the method of Ref.\;\cite{yao2019}. At a short distance, we use the SVD method to solve Eq.\;(\ref{3}) in the interval $[R_{i},R_{i+1}]$. The total wave function $\psi_{\upsilon'}(R;\Omega)$ is expanded in terms of the discrete variable representation (DVR) basis $\pi_{i}$ and the channel functions $\Phi_{\nu}(R;\Omega)$ as\\
\begin{equation}
\label{8}
\psi_{\upsilon'}(R;\Omega)=\sum\limits_{i}^{N_{DVR}}\sum\limits_{\nu}^{N_{chan}}C_{i\mu}^{\upsilon'}\pi_{i}(R)\Phi_{\nu}(R;\Omega),\\
\end{equation}
where $N_{DVR}$ is the number of DVR basis functions and $N_{chan}$ is the number of included channel functions. Inserting Eq.\;(\ref{8}) into Eq.\;(\ref{3}), the three-body Schr$\mathrm{\ddot{o}}$dinger equation can be rewritten as:
\begin{equation}
\label{9}
\sum\limits_{j}^{N_{DVR}}\sum\limits_{\mu}^{N_{chan}}T_{ij}O_{i\nu,j\mu}C_{i\nu}^{\upsilon'}+U_{\nu}(R_{i})C_{i\nu}^{\upsilon'}=E^{\upsilon'}C_{i\nu}^{\upsilon'},\\
\end{equation}
where the kinetic energy matrix elements can be expressed as\\
\begin{equation}
\label{10}
T_{ij}=\frac{1}{2\mu}\int_{R_{i}}^{R_{i+1}}\frac{d}{dR}\pi_{i}(R)\frac{d}{dR}\pi_{j}(R)dR,\\
\end{equation}
$R_{i}$ and $R_{i+1}$ are the boundaries of the calculation box, and the overlap matrix elements between the adiabatic channels defined at different quadrature points are:\\
\begin{equation}
\label{11}
O_{i\nu,j\mu}=\langle\Phi_{\nu}(R_{i};\Omega)|\Phi_{\mu}(R_{j};\Omega)\rangle.\\
\end{equation}

At large distances, we use the traditional adiabatic hyperspherical method to solve Eq.\;(\ref{3}). Substituting Eq.\;(\ref{6}) into Eq.\;(\ref{3}), a set of coupled ordinary differential equations are obtained:\\
\begin{equation}
\label{12}
\left[-\frac{1}{2\mu}\frac{d^{2}}{dR^{2}}+U_{\nu}(R)-E\right]F_{\nu,\upsilon'}(R)-\frac{1}{2\mu}\sum\limits_{\mu}\left[2P_{\mu\nu}(R)\frac{d}{dR}+Q_{\mu\nu}(R)\right]F_{\mu\upsilon'}(R)=0,\\
\end{equation}
where the nonadiabatic couplings $P_{\mu\nu}$ and $Q_{\mu\nu}$ that control the inelastic transition and the width of resonance supported by adiabatic potential $U_{\nu}(R)$ are \\
\begin{equation}
\label{13}
P_{\mu\nu}(R)=\int\Phi_{\mu}(R;\Omega)^{*}\frac{\partial}{\partial R}\Phi_{\nu}(R;\Omega)d\Omega,\\
\end{equation}
and\\
\begin{equation}
\label{14}
Q_{\mu\nu}(R)=\int\Phi_{\mu}(R;\Omega)^{*}\frac{\partial^{2}}{\partial R^{2}}\Phi_{\nu}(R;\Omega)d\Omega.\\
\end{equation}
respectively.

The effective potentials are usually more physical than adiabatic potentials and can thus give physical asymptotic behaviors of the system at large R, which are defined as:\\
\begin{equation}
\label{15}
W_{\nu\nu}(R)=U_{\nu}(R)-\frac{\hbar^{2}}{2\mu}P_{\nu\nu}^{2}(R).\\
\end{equation}

During scattering calculation, the $\underline{\mathcal{R}}$ matrix is a fundamental quantity that can be subsequently used to determine the scattering matrix $\underline{\mathcal{S}}$, which is the main goal of the study. We use the $R$-matrix propagation method over interval $[a_{1},a_{2}]$ to calculate the relative $\underline{\mathcal{R}}$ matrix from $a_1$ to $a_2$ as follows:
\begin{equation}
\label{16}
\underline{\mathcal{R}}(a_2)=\underline{\mathcal{R}}_{22}-\underline{\mathcal{R}}_{21}\left[\underline{\mathcal{R}}_{11}+\underline{\mathcal{R}}(a_1)\right]^{-1}\underline{\mathcal{R}}_{12}.
\end{equation}

After we have the $\underline{\mathcal{R}}$ matrix at large distances, the physical scattering matrix $\underline{\mathcal{S}}$ and the relative reaction matrix $\underline{\mathcal{K}}$ can be determined by applying asymptotic boundary conditions:\\
\begin{equation}
\label{17}
\underline{\mathcal{K}}=(\underline{f}-\underline{f}'\underline{\mathcal{R}})(\underline{g}-\underline{g}'\underline{\mathcal{R}})^{-1},
\end{equation}
\begin{equation}
\label{18}
\underline{\mathcal{S}}=(\underline{1}+i \underline{\mathcal{K}})(\underline{1}-i\underline{\mathcal{K}})^{-1}.\\
\end{equation}
where $\underline{f}$, $\underline{f}'$, $\underline{g}$, and $\underline{g}'$ are diagonal matrices with matrix elements given by energy-normalized asymptotic solutions $f_{\nu\nu'}$ and $g_{\nu\nu'}$ and derivatives $f'_{\nu\nu'}$ and $g'_{\nu\nu'}$, respectively. $f_{\nu\nu'}$ and $g_{\nu\nu'}$ are described in terms of spherical Bessel and Neumann functions: \\
\begin{equation}
\label{19}
f_{\nu\nu'}=\sqrt{\frac{2\mu k_{\nu}}{\pi}}R j_{l_{\nu}}(k_{\nu}R)\delta_{\nu\nu'},\\
\end{equation}
and\\
\begin{equation}
\label{20}
g_{\nu\nu'}=\sqrt{\frac{2\mu k_{\nu}}{\pi}}R n_{l_{\nu}}(k_{\nu}R)\delta_{\nu\nu'}.\\
\end{equation}
respectively. For the $\nu$th recombination channel, $l_\nu$ is the angular momentum of the third atom relative to the dimer and is given by Eq.\;(\ref{21}), and $k_{\nu}=\sqrt{2 \mu\left(E-E_{2 b}\right)}$.
For the $\nu$th entrance channel, $l_{\nu}=\lambda_{\nu}+3/2$, and $k_{\nu}=\sqrt{2 \mu E}$.\\
\section{Results and Discussion}

\subsection{Atom-dimer elastic scattering cross sections}
In this section, we focus on the case where $^{87}$Rb-$^{87}$Rb has a finite positive scattering length $100\,a_{0}$ and the Rb-K scattering length is varied at the negative side of the Feshbach resonance. Figure\;\ref{fig1a} shows the hyperspherical potential curves for the $^{87}$Rb-$^{87}$Rb-$^{40}$K system with $a_{\scriptscriptstyle\textsl{RbK}}=-450\,a_{0}$ and $a_{\scriptscriptstyle\textsl{RbRb}}=100\,a_{0}$. The black line in Fig.\;\ref{fig1} denotes the recombination channel (labeled by f), which corresponds to a K atom relative to the RbRb dimer. The effective potentials for the recombination channels have asymptotic behavior at large R as,\\
\begin{equation}
 \label{21}
 W=\frac{l_{f}(l_{f}+1)}{2\mu R^{2}}+E_{2b},\\
 \end{equation}
where $E_{2b}$ is the two-body bound state energy and $l_{f}$ is the corresponding angular momentum of the third atom relative to the dimer. The dashed lines represent the three-body continuum channels (labeled i), which corresponds to the case in which the three atoms are far away from each other as $R \to \infty$. The corresponding potential curves behave as\\
\begin{equation}
 \label{22}
 W=\frac{\lambda(\lambda+4)+15/4}{2\mu R^{2}}.\\
 \end{equation}
Here, $\lambda(\lambda+4)$ is the eigenvalue of the grand angular momentum operator $\Lambda^{2}$, and $\lambda$ are nonnegative integer values, which are restricted by the requirements of permutation symmetry.\\
\begin{figure}[htbp]
  \centering
\subfigure{
  \includegraphics[width=7.8cm]{1a.eps}
  \label{fig1a}
}
\subfigure{
  \includegraphics[width=7.8cm]{1b.eps}
  \label{fig1b}
}
\caption{(Color online) The $J^{\Pi}=0^{+}$ adiabatic hyperspherical potential curves U$_{\nu}$(R) for $^{87}$Rb-$^{87}$Rb-$^{40}$K with $a_{\scriptscriptstyle\textsl{RbRb}}=100\,a_{0}$. The values of $\lambda$ indicate the asymptotic behavior of the potential curves representing the three-body continuum channels, as given in Eq.\;(\ref{22}). (a) The potential curves obtained with $a_{\scriptscriptstyle\textsl{RbK}}=-450\,a_{0}$. The lowest adiabatic potential curve denotes a recombination channel that corresponds to an RbRb dimer and a K atom. (b) The potential curves with $a_{\scriptscriptstyle\textsl{RbK}}=178\,a_{0}$. The lowest adiabatic potential curve denotes the RbRb + K channel, and the red solid line is the RbK + Rb channel.}
  \label{fig1}
\end{figure}
In these calculations, the basis sets $N_{\theta}=168$ and $N_{\phi}=434$ are chosen, and the potential curves have at least six significant digits.
We denote the atom-dimer scattering length by $a_{ad}$. The cross-section for atom-dimer scattering is expressed in terms of the $S$ matrix as\\
\begin{equation}
\label{23}
\sigma_{2}=\sum\limits_{J,\Pi}\sigma_{2}^{J,\Pi}=\sum\limits_{J,\Pi}\frac{(2J+1)\pi}{k_{ad}^{2}}|S_{f\leftarrow i}^{J,\Pi}-1|^{2},\\
\end{equation}
where $\sigma_{2}^{J,\Pi}$ is the partial atom-dimer elastic scattering cross section corresponding to the $J^{\Pi}$ symmetry, $k_{ad}=\sqrt{2\mu_{ad}(E-E_{d})}$ is the atom-dimer wavenumber, $\mu_{ad}$ is the atom-dimer reduced mass, and $E_{d}$ is the binding energy of the weakly bound molecular state. The threshold behavior of $\sigma_{2}^{J,\Pi}$ is determined solely by $J$ and follows the standard Wigner threshold law\;\cite{Esry2001Dec}\\
\begin{equation}
\label{24}
\sigma_{2}^{J,\Pi} \sim (E-E_{d})^{2J}.\\
\end{equation}
The atom-dimer scattering length is defined as\\
\begin{equation}
\label{25}
a_{ad}=-\lim\limits_{k_{ad} \to 0}\frac{\tan\delta_{0}}{k_{ad}},\\
\end{equation}
where $\delta_{0}$ is the phase shift for atom-dimer elastic scattering and is related to the diagonal $S$-matrix element by the formula\\
\begin{equation}
\label{26}
S_{f\leftarrow i}^{0+}=\exp{(2i\delta_{0})}.\\
\end{equation}
We typically use $15$ channels and $349$ sectors distributed from $R = 18\; a_{0}$ to $R = 3100\; a_{0}$. In the matching procedure of this section, the asymptotic wave functions are expressed in rotated Jacobi coordinates, which improve the convergence of observables as a function of matching distance\;\cite{cyzhao2022}. The convergence of the scattering length with respect to the final matching distance $R_{m}$ is shown in Table\;\ref{tableI} when the Rb-K scattering length is $-450\,a_{0}$. The atom-dimer scattering length values reach converged values at $R_{m}$ = 2974 $a_{0}$.

\begin{table*}
\caption{\label{t4} Convergence of the scattering length between K and RbRb $a_{ad}$ with respect to the matching distances $R_{m}$. The Rb-K and Rb-Rb scattering length is $a_{\scriptscriptstyle\textsl{RbK}}=-450\,a_{0}$ and $a_{\scriptscriptstyle\textsl{RbRb}}=100 \,a_{0}$, respectively.}
\begin{ruledtabular}
\begin{tabular}{cllllllll}
&\multicolumn{1}{c}{$R_{m}$}&\multicolumn{1}{c}{2543}&\multicolumn{1}{c}{2752}&\multicolumn{1}{c}{2974} \\
\hline
& $a_{ad}$ &467.02  &467.34 &467.36 \\
\end{tabular}
\end{ruledtabular}
\label{tableI}
\end{table*}

In Fig.\;\ref{fig2} we present cross sections for $^{40}$K + $^{87}$Rb$^{87}$Rb elastic scattering as a function of the collision energy $E$ with $a_{\scriptscriptstyle\textsl{RbRb}}=100\,a_{0}$ and $a_{\scriptscriptstyle\textsl{RbK}}=-450\,a_{0}$. The data show that $\sigma_{2}^{0+}$ obeys the threshold behavior as $\sigma_{2}^{0+} \propto E^{0}$ and approaches the constant value $ 4\pi a_{ad}^{2}$ in the ultracold limit.\\

\begin{figure}[htbp]
 \centering
 \includegraphics[width=0.6\linewidth]{2.eps}
\caption{Cross sections $\sigma_{2}$ for $^{40}$K + $^{87}$Rb$^{87}$Rb elastic scattering as a function of the collision energy $E$ in $J^{\Pi}=0^{+}$ symmetry with $a_{\scriptscriptstyle\textsl{RbRb}}=100\,a_{0}$ and $a_{\scriptscriptstyle\textsl{RbK}}=-450\,a_{0}$.}
\label{fig2}
\end{figure}

Fig.\;\ref{fig3a} shows the cross section for $^{40}$K + $^{87}$Rb$^{87}$Rb elastic scattering as a function of $^{87}$Rb-$^{40}$K scattering length at the negative side of the Feshbach resonance. A peak at approximately $a_{\scriptscriptstyle\textsl{RbK}}= - 155\,a_{0}$ is shown. We interpret the peak as being caused by the appearance of a three-body bound state at the atom-dimer threshold ($^{40}$K + $^{87}$Rb$^{87}$Rb). Owing
to the presence of the near-threshold trimer state, the atom-dimer scattering length is expected to diverge. As illustrated in Fig.\;\ref{fig3b}, the $^{40}$K + $^{87}$Rb$^{87}$Rb elastic scattering lengths $a_{ad}$ are divergent at the same Rb-K scattering length value. A schematic diagram of the relevant three-body energy spectrum as a function of Rb-K scattering length is shown in Fig.\;\ref{fig3c}, illustrating the energies of trimer states (the blue solid curve) and atom-dimer thresholds (the red dashed curve). The energy spectrum gives clear signatures at the threshold locations, which are consistent with the peaks of the cross section as well as the divergence position of the atom-dimer scattering length. Thus, intriguing questions are whether the obtained three-body bound state supported by the lowest hyperspherical potential curve is an Efimov state, and whether it is related to the observed atom-dimer resonance on the positive Rb-K scattering length (or emerges from it). It is known that Efimov states have universal properties that generally require scattering lengths much larger than $r_{\scriptscriptstyle\textsl{vdW}}$ and energies much smaller than $E_{\scriptscriptstyle\textsl{vdW}}$. The Rb-Rb-K system can be characterized by a length $r_{\scriptscriptstyle\textsl{vdW,RbRb}}=83\,a_{0}$ and a corresponding energy $E_{\scriptscriptstyle\textsl{vdW,RbRb}}= 9.25\times 10^{-10}$ a.u. The three-body state is located at $|a_{\scriptscriptstyle\textsl{vdW,RbK}}|= 155\, a_{0}< 2\, r_{\scriptscriptstyle\textsl{vdW,RbRb}} $, with a binding energy larger than $E_{\scriptscriptstyle\textsl{vdW,RbRb}}$ by one order of magnitude. Thus, the three-body state is not an Efimov state, nor is it related to the Efimov resonance in the positive Rb-K scattering length. Analogous to a typical two-body Feshbach resonance, such a three-body bound state could offer the possibility to tune the atom-dimer interaction from attractive to repulsive with the atomic two-body scattering length tuned at small range adjustment.

\begin{figure}[htbp]
\centering
\subfigure{
\includegraphics[width=5.0cm]{3a.eps}
\label{fig3a}
}
\subfigure{
\includegraphics[width=5.0cm]{3b.eps}
\label{fig3b}
}
\subfigure{
\includegraphics[width=5.4cm]{3c.eps}
\label{fig3c}
}
\caption{(Color online)
The $^{40}$K + $^{87}$Rb$^{87}$Rb elastic scattering observables with the $J^{\Pi}=0^{+}$ symmetry as a function of Rb-K scattering lengths when $a_{\scriptscriptstyle\textsl{RbRb}} = 100\,a_{0}$.
All the data are obtained at low collision energy $E=100$ nk. (a) Cross sections: the red dashed line shows the position of the peak; (b) the atom-dimer scattering lengths; and the red dashed line shows the divergence position; (c) three-body energy spectrum. The red dashed line is the RbRb energy, and the solid line is the three-body state energy. The inset shows the Rb-K scattering length at which the three-body state disappears.}
\label{fig3}
\end{figure}

\subsection{Three-body recombination rates}
We now consider three-body recombination, which is widely used to reveal Efimov physics in cold atoms. The event rate constant for three-body recombination is \\
\begin{equation}
\label{27}
K_{3}^{(f\leftarrow i)}=\sum\limits_{J,\Pi}K_{3}^{J\Pi}=2!\sum\limits_{J,\Pi}\sum\limits_{f,i}\frac{32(2J+1)\pi^{2}}{\mu k^{4}}|S_{f \leftarrow i}^{J,\Pi}|^{2}.\\
\end{equation}
Here, $K_{3}^{J\Pi}$ is the partial recombination rate corresponding to the $J^{\Pi}$ symmetry, and $k=(2\mu E)^{1/2}$ is the hyperradial wavenumber in the incident channel. $S_{f \leftarrow i}^{J,\Pi}$ represents the scattering matrix element from the initial three-body continuum channels $i$ to the final atom-dimer channel $f$ for the $J^{\Pi}$ symmetry. The factor (2!) derives from the number of indistinguishable bosonic particles. The asymptotic form of the effective hyperradial potentials in $R\to \infty$ allows for generalization of Wigner's threshold law\;\cite{Esry2001Dec,Incao2005Jun} for recombination (i.e., the low-energy behavior of the recombination rates), which leads to the partial recombination rates $K_{3}^{J\Pi}$ satisfying\\
\begin{equation}
\label{28}
K_{3}^{J\Pi}\propto E^{\lambda_{min}}
\end{equation}
near threshold, where $\lambda_{min}$ is the minimum value of $\lambda$ in Eq.\;(\ref{22}).

We first study the three-body recombination rates $K_{3}^{(f\leftarrow i)}$ for the $^{87}$Rb$^{87}$Rb$^{40}$K system with $J^{\Pi}=0^{+}$ symmetry on the positive Rb-K scattering length side. Figure\;\ref{fig1b} shows the hyperspherical potential curves for the $^{87}$Rb-$^{87}$Rb-$^{40}$K system, in which the Rb-K scattering length is $a_{\scriptscriptstyle\textsl{RbK}}=178\,a_{0}$ and Rb-Rb interacts via the s-wave scattering length $a_{\scriptscriptstyle\textsl{RbRb}}=100\,a_{0}$. It is shown that there are two recombination channels: one being the weakly bound $^{87}$Rb$^{40}$K + $^{87}$Rb channel and the other being the deeply bound channel $^{40}$K + $^{87}$Rb$^{87}$Rb. Figure\;\ref{fig4} shows the three-body recombination rate $K_{3}^{(f\leftarrow i)}$ as a function of the collision energy $E$ for the $^{87}$Rb$^{87}$Rb$^{40}$K system with $J^{\Pi}=0^{+}$ symmetry. At lower collision energies, the recombination rates $K_{3}^{(f\leftarrow i)}$ are constant and follow the Wigner threshold law prediction, $K_{3}^{(f \leftarrow i)}\propto E^{0}$.\\

\begin{figure}
 \centering
   \includegraphics[width=0.6\linewidth]{4K.eps}
\caption{(Color online) The three-body recombination rate $K_{3}^{(f \leftarrow i)}$ as a function of the collision energy $E$ for the $^{87}$Rb$^{87}$Rb$^{40}$K system in the $J^{\Pi}=0^{+}$ symmetry with $a_{\scriptscriptstyle\textsl{RbRb}}=100\,a_{0}$ and $a_{\scriptscriptstyle\textsl{RbK}}=178\,a_{0}$. The black and red solid lines indicate three-body recombination rates for the sum of different incident channels $i$ to particular recombination channels $f=1,2$, respectively. The blue solid line indicates the sum of the above two three-body recombination rates. The recombination rates $K_{3}^{(f \leftarrow i)}$ follow the Wigner threshold law prediction, $K_{3}^{(f \leftarrow i)}\propto E^{0}$ in the low-energy region.}
   \label{fig4}
\end{figure}

Figure\;\ref{fig5a} shows the $a_{\scriptscriptstyle\textsl{RbK}}$ dependence of the three-body recombination rates at positive Rb-K scattering lengths with the Rb-Rb interaction fixed at $a_{\scriptscriptstyle\textsl{RbRb}}=100\,a_{0}$. The data are obtained at a collision energy of $E=100$ nk. The filled black squares represent recombination into the deeply bound channel RbRb + K, red solid circles represent recombination into the weakly bound channel RbK + Rb, and blue solid triangles depict their sums. The recombination rates of the deep channel (RbRb + K) are close to those of the shallow channel (RbK + Rb) at small Rb-K scattering lengths. Then, the recombination rates of the deep channel begin to decrease when the Rb-K scattering lengths are tuned larger. A minima appears at approximately $a_{\scriptscriptstyle\textsl{RbK}}=1641\,a_{0}$ in the deep recombination channel. However, this recombination minima has no effect on the total recombination rates since it is smaller than that into the weakly bound channel by at least 3 orders of magnitude. An Efimov recombination minima for total rates arises at approximately $a_{\scriptscriptstyle\textsl{RbK}} = 3638\,a_{0}$. For $a > 0$, the three-body recombination minimum is a well-known feature of Efimov physics\;\cite{Zaccanti2009,Pollack2009}, which is explained as the destructive interference between two different decay pathways.
Figure\;\ref{fig5b} shows the dependence of Efimov recombination minima on the collision energy. The recombination minima is at $a_{\scriptscriptstyle+} = 3264\;a_{0}$ with the collision energy of $E = 50$ nk and is $ a_{\scriptscriptstyle+} = 3638\,a_{0} $ for $E = 100$ nk. However, this minimum is less evident at the collision energy of $E = 300$ nk.

\begin{figure}[htbp]
\centering
\subfigure{
\includegraphics[width=7.8cm]{5a.eps}
\label{fig5a}
}
\subfigure{
\includegraphics[width=7.8cm]{5b.eps}
\label{fig5b}
}
\caption{(Color online) (a) The $J^{\Pi} = 0^{+}$ three-body recombination rate $K_{3}^{(f \leftarrow i)}$ is shown as a function of Rb-K scattering length $a_{\scriptscriptstyle\textsl{RbK}}$ for $E=100$ nk. The Rb-Rb interaction is fixed at $a_{\scriptscriptstyle\textsl{RbRb}}=100\,a_{0}$. The filled squares (black) represent recombination of the deep channel Rb + Rb+ K $\to$ RbRb + K, the filled circles (red) represent recombination of the shallow channel Rb + Rb + K $\to$ RbK + Rb, and the filled triangles (blue) depict the sum of the channels. The local minima is at $a_{\scriptscriptstyle\textsl{RbK}}=1641\,a_{0}$ in the deep recombination channel. An Efimov recombination minima on the total rates is shown at $a_{\scriptscriptstyle\textsl{RbK}} = 3638 \,a_{0}$. (b) Magnifying the position of Efimov recombination minima for different collision energies.}
\label{fig5}
\end{figure}

Aside from calculating the positions of the three-body recombination minima, it is desirable to compare the actual shapes of the $a$-dependence of the loss rate with the experimental measurements.
To date, the three-body recombination rate in $^{87}$Rb-$^{40}$K mixtures has been measured by Ref.\;\cite{Bloom2013,Hu2014}.

In Figure.\;\ref{fig6} we show our numerically calculated three-body recombination rates $K_{3}^{(f \leftarrow i)}$ for the $^{87}$Rb$^{87}$Rb$^{40}$K system compared with experimental data. The solid circles show experimental data obtained from Ref.\;\cite{Bloom2013}, and the solid squares show numerical results at a collision energy of $E=100$ nk. The solid line shows an $a^{4}$ dependence. At positive scattering lengths, our numerical results are multiplied by 8, and at negative scattering lengths, our numerical results are multiplied by 105. These scaled numerical results show reasonable agreement with the experimental results. The expected $a^{4}$ dependence of recombination rates is shown for both positive and negative scattering lengths. The experimental and numerical data both show that there are no Efimov features at negative $a$ between $-200\,a_{0}$ and $-3000\,a_{0}$ or at positive $a$ between $200\,a_{0}$ and $3000\,a_{0}$. A possible Efimov recombination minima at $a_{\scriptscriptstyle\textsl{RbK}}=3638\, a_{0}$ is displayed in our numerical calculations.

\begin{figure}
 \centering
   \includegraphics[width=0.9\linewidth]{6K.eps}
\caption{(Color online) The three-body recombination rates $K_{3}^{(f \leftarrow i)}$ for the $^{87}$Rb$^{87}$Rb$^{40}$K system as functions of $^{87}$Rb-$^{40}$K scattering lengths with $J^{\Pi}=0^{+}$ symmetry. The solid circles are experimental data obtained from Ref.\;\cite{Bloom2013}, and the solid squares are numerical results at a collision energy of $E=100$ nk. The solid line shows an $a^{4}$ dependence.}
   \label{fig6}
\end{figure}

\subsection{Discussion of the universal relationship between the different Efimov features }
The universal relations between various Efimov features have been predicted in the zero-range limit by Helfrich
\textit{et al.}\;\cite{Helfrich2010} using an effective-field-theory framework. The predicted universal relationship for the
location of Efimov features due to the same trimer state is $|a_{\scriptscriptstyle-}|/a_{\scriptscriptstyle*}=1.96$. In $^{40}$K-$^{87}$Rb mixtures, an atom-dimer resonance was observed at $a_{\scriptscriptstyle*}= 230\,a_0$ \;\cite{Zirbel2008,Bloom2013,Kato2017}. If the atom-dimer resonances are associated with the lowest Efimov state, universal scaling predicts $a_{\scriptscriptstyle-}=-450\,a_{0}$, which is excluded by several measurements. The breakdown of the universal relation is expected to be associated with the fact that the universal relationship does not consider the effect of finite effective ranges and nonresonant RbRb interactions. Thus, a set of universal relations that connects the different Efimov features, which include the effects of the finite interspecies effective range and the finite intraspecies scattering length, was derived by Ref.\;\cite{Acharya2016} (Eq. (25) of the paper). This equation requires two Efimov features as inputs.

For the $^{87}$Rb-$^{87}$Rb-$^{40}$K system, $\lambda=122.7$, and the other parameters are provided in Ref.\;\cite{Acharya2016}.
Taking $a_{\scriptscriptstyle*}=230\,a_{0}$, $a_{\scriptscriptstyle+}=3638\,a_{0}$($E=100$ nk), $a_{22}=a_{\scriptscriptstyle\textsl{RbRb}}=100\,a_{0}$, and $r_{0}=r_{\scriptscriptstyle \textsl{vdW,RbK}}=72\,a_{0}$ as inputs, the first Efimov resonance is at $a_{\scriptscriptstyle-}=-52679\,a_{0}$ predicted through Eq. (25) of the paper. For the collision energy of $E=50$ nk, the recombination minima predicted by us resides at $a_{\scriptscriptstyle+}=3264\, a_{0}$, and the first Efimov resonance would be at $a_{\scriptscriptstyle-}=-48315\,a_{0}$ predicted by the universal relation. These results are consistent with the predicted value of $a_{\scriptscriptstyle-}<-30000\,a_{0}$ from Wang \textit{et al.}\;\cite{WangYujun2012} and outside current experiment reach. Therefore, the universal relation including the finite range effects and the nonresonant intraspecies scattering length can give reasonable predictions, and no breakdown universal relations occur.
\section{Summary}
In this paper, we studied the atom-dimer elastic scattering cross sections, scattering lengths and three-body recombination rates in $^{87}$Rb-$^{40}$K mixtures with total angular momentum $J=0$. We use the R-matrix propagation method in the hyperspherical coordinate frame based on the Lennard-Jones model potential with a vdW tail for the interactions between atoms. The finite intraspecies scattering length is considered with the $^{87}$Rb-$^{87}$Rb s-wave scattering length fixed at $a_{\scriptscriptstyle\textsl{RbRb}} = 100\,a_{0}$.

To make clear whether there is an Efimov state disappearing at the $^{40}$K + $^{87}$Rb$^{87}$Rb threshold, we investigate the $^{40}$K + $^{87}$Rb$^{87}$Rb elastic scattering at a negative $^{87}$Rb-$^{40}$K scattering length. We calculate the scattering lengths and cross sections in the range from $-100\,a_{0} $ to $-1000\,a_{0}$. An atom-dimer elastic scattering resonance at $a_{\scriptscriptstyle\textsl{RbK}} = - 155\,a_0$ is found, which is due to the presence of the near-threshold trimer state at the atom-dimer threshold ($^{40}$K + $^{87}$Rb$^{87}$Rb). The three-body state is located at $|a_{\scriptscriptstyle\textsl{RbK}}|= 155\,a_{0}< 2 r_{\scriptscriptstyle \textsl{vdW,RbRb}} $, and its binding energy is larger than $E_{\scriptscriptstyle \textsl{vdW,RbRb}}$ by one order of magnitude. Thus, the three-body state is not an Efimov state, and it is not related to the previously observed atom-dimer Efimov resonance at $a_{\scriptscriptstyle*} = 230(30)\,a_{0}$.

We further investigate the three-body recombination, which is widely used to reveal Efimov physics in cold atoms. Our calculations find an Efimov recombination minima at a positive interspecies scattering length. Taking the position of the recombination minima and the experimentally observed atom-dimer resonance $a_{\scriptscriptstyle*} = 230(30)\,a_{0}$ as inputs, the universal relations accounting for the finite range effects and finite intraspecify scattering length are tested. An Efimov resonance arising at the three-body threshold $a_{\scriptscriptstyle-}=-52679\,a_{0}$($E=100$ nk) is predicted, in agreement with previous theoretical predictions.
Testing shows that the main reason for the breakdown of universal relations in $^{40}$K-$^{87}$Rb predicted by the zero-range to be finite range effects and the intraspecify scattering length. This study also demonstrates the importance of finite range effects in $^{40}$K-$^{87}$Rb mixtures. Finally, the present results help elucidate the universality of three-body parameters for heteronuclear systems with finite intraspecies scattering lengths.

\section{Acknowledgments}
We thank C. H. Greene, P. G. Panos and Chen Ji for helpful discussions.
Hui-Li Han was supported by the National Natural Science Foundation of China under Grant No. 11874391.
Ting-Yun Shi was supported by the Strategic Priority Research Program of the Chinese Academy of Sciences under Grant No. XDB21030300.