# Multidimensional heterogeneity learning for count value tensor data with
applications to field goal attempt analysis of NBA players

Guanyu Hu, Yishu Xue, and Weining Shen  Department of Statistics, University
of Missouri - Columbia Google Inc Department of Statistics, University of
California, Irvine

###### Abstract

We propose a multidimensional tensor clustering approach for studying how
professional basketball players’ shooting patterns vary over court locations
and game time. Unlike most existing methods that only study continuous-valued
tensors or have to assume the same cluster structure along different tensor
directions, we propose a Bayesian nonparametric model that deals with count-
valued tensors and projects the heterogeneity among players onto tensor
dimensions while allowing cluster structures to be different over directions.
Our method is fully probabilistic; hence allows simultaneous inference on both
the number of clusters and the cluster configurations. We present an efficient
posterior sampling method and establish the large-sample convergence
properties for the posterior distribution. Simulation studies have
demonstrated an excellent empirical performance of the proposed method.
Finally, an application to shot chart data collected from 191 NBA players
during the 2017–2018 regular season is conducted and reveals several
interesting insights for basketball analytics.

##  1 Introduction

There is a rapid growth in sports analytics over the recent years thanks to
the fast development in game tracking technologies (Albert et al., 2017). New
statistics and machine learning methods are largely needed to address a range
of questions from game outcome prediction to player performance evaluation
(Cervone et al., 2014; Franks et al., 2015; Cervone et al., 2016), and to in-
game strategy planning (Fernandez and Bornn, 2018; Sandholtz et al., 2020).

Our focus in this paper is to study the National Basketball Association (NBA)
players’ shooting patterns, in particular, how they change over shooting
locations and game time (e.g., first quarter versus clutch time). In
professional basketball research, shooting pattern (or field goal attempt) is
a fundamental metric since it often provides valuable insights to players,
coaches and managers, e.g., players and coaches will be able to obtain a
better understanding of their current shooting choices and hence develop
further offensive/defensive plans accordingly, while managers can make better
data-informed decisions on player recruitment. In the literature, it is common
to employ spatial and spatial-temporal models (Reich et al., 2006; Miller et
al., 2014; Jiao et al., 2020) to study the spatial and temporal dependence in
the field goal attempt data. The key novelty of our paper is to introduce a
new Bayesian tensor multidimensional clustering method that studies the
heterogeneity of shooting patterns among players.

Our starting point is to recognize that the field goal attempt data in
basketball games enjoys a natural tensor structure. For example, we can divide
the basketball half court into regions under the polar coordinate system and
then summarize the number of field goal attempts over each region during each
of the four quarters of the game as entries in a three-way tensor, where each
tensor direction corresponds to the shooting distance, angle, and game time
(quarter). One advantage of considering a tensor representation is that the
spatial-temporal dependence structure is automatically considered as part of
the tensor structure. Compared to the most existing works that rely on
Cartesian coordinate system (Miller et al., 2014; Jiao et al., 2020; Yin et
al., 2020; Hu et al., 2020; Yin et al., 2020), our approach certainly makes
more sense since shooting angle and distance are two important factors that
affect the shooting selection of professional players (Reich et al., 2006).
Studying the change of shooting patterns over time is also meaningful, e.g.,
Stephen Curry has made a comparable fewer number of attempts in fourth quarter
during regular season simply because Gold State Worriers has often established
a significant lead at the end of the third quarter during the 2017-2018
regular season.

Tensor models have received a great deal of attention in machine learning and
statistics literature (Kolda and Bader, 2009; Sun et al., 2014; Bi et al.,
2021). For tensor clustering problems, most existing work (Sun and Li, 2019;
Chi et al., 2020; Mai et al., 2021) either works with a single tensor or
assumes that the clustering structure is the same across different tensor
directions. It is also common that they only consider tensors with entries
that takes continuous or binary values. In terms of model estimation,
clustering is often based on solving a regularized optimization problem, which
requires pre-specifying the number of clusters or choosing the cluster number
based on certain ad hoc criteria. There are some recent papers on Bayesian
tensor models, e.g., Guhaniyogi et al. (2017); Spencer et al. (2019);
Guhaniyogi (2020). However, all of them are in the regression context, hence
cannot be directly applied to solve our problem. Our proposed approach differs
with aforementioned methods in several ways. First, we consider a flexible
multidimensional tensor clustering problem, which allows different clustering
structures over tensor directions. We believe this is a meaningful relaxation
in many applications, e.g., basketball players’ shooting choice may differ
significantly in terms of shooting distance, angle and game time depending on
players’ position, shooting preference and role in the team. Moreover, we
focus on count-valued tensors (rather than continuous-valued tensors) for the
obvious reason that the number of shot attempts is the main outcome of
interest in our application. Thirdly, our model is fully probabilistic, which
allows an easier interpretation compared to optimization-based method. In
particular, we consider a Bayesian nonparametric model under the mixture of
finite mixtures framework (MFM; Miller and Harrison, 2018), which allows
simultaneous estimation and inference for the number of clusters and the
associated clustering configuration for each direction (e.g. distance, angle,
and quarter). We show that the posterior inference can be conducted
efficiently for our method and demonstrate its excellent empirical performance
via simulation studies. We also establish posterior convergence properties for
our method. Finally, our proposed method reveals several interesting data-
driven insights of field goal attempts data which are useful for professional
basketball players, coaches, and managers.

Main contributions and broad impact: (1) We are among the first (to the best
of our knowledge) to introduce tensor models to sports analytics. It is our
hope that this paper can contribute to promoting more use of tensor methods in
different sport applications. (2) We develop a novel multidimensional tensor
clustering approach that allows different clustering structures over tensor
directions and handles count-valued data. The proposed method is fully
Bayesian, which renders convenient inference on the number of clusters and the
clustering structure. (3) We provide a large-sample theoretical justification
for our method by showing posterior consistency for the cluster number and
contraction rate for the mixing distributions. These results are new for
Bayesian tensor models.

##  2 Method

###  2.1 Probabilistic multi-dimensional tensor clustering

We treat the shooting chart data as three-way tensors and discuss a multi-
dimensional clustering approach in this section. Note that each direction of
the tensor represents the shooting angle, distance to the basket, and one of
the four quarters in the game. Our proposed method can be conveniently
extended to study general multi-way tensor data as well. Let $\bm{Y}$ be a
$p_{1}\times p_{2}\times p_{3}$ tensor with each element $Y_{ijk}$ only taking
count values for $i=1,\ldots,p_{1};j=1,\ldots,p_{2};k=1,\ldots,p_{3}$. It is
natural to consider a Poisson distribution with a mean parameter represented
as a rank-one tensor, that is,

$\displaystyle\bm{Y}\sim\text{Poisson}(\bm{\gamma}_{1}\circ\bm{\gamma}_{2}\circ\bm{\gamma}_{3}),$
(1)

where $\circ$ denotes the outer product between two vectors, and
$\bm{\gamma}_{1}\in\mathbb{R}_{+}^{p_{1}},\bm{\gamma}_{2}\in\mathbb{R}_{+}^{p_{2}},\bm{\gamma}_{3}\in\mathbb{R}_{+}^{p_{3}}$.
Model (1) can also be viewed as a Poisson regression model where the mean
parameter corresponds to an analysis of variance (ANOVA) model with main
effects only, that is,
$\log\text{E}(Y_{ijk})=\log\gamma_{1,i}+\log\gamma_{2,j}+\log\gamma_{3,k}$ for
$1\leq i\leq p_{1},1\leq j\leq p_{2},1\leq k\leq p_{3}$. By ignoring the
interaction effects, the number of parameters is effectively reduced from
$p_{1}p_{2}p_{3}$ to $(p_{1}+p_{2}+p_{3})$. Hence it renders parsimonious
parameter estimation and easy interpretation in our NBA application study,
i.e., the main effects
$\log\bm{\gamma}_{1},\log\bm{\gamma}_{2},\log\bm{\gamma}_{3}$ correspond to
the additive effect of shooting distance, angle, and game time (quarter).

In order to learn the multidimensional heterogeneity pattern, we consider
three independent mixture of finite mixtures (MFM; Miller and Harrison, 2018)
priors on $\bm{\gamma}_{1},\bm{\gamma}_{2},\bm{\gamma}_{3}$ such that the
clustering pattern in those three directions can be learned separately. Here
we present a brief introduction to MFM without getting into more details.
Given $n$ observations, we consider $z_{1},\ldots,z_{n}$ as their clustering
labels, e.g., $z_{1}=z_{2}=z_{4}$ would mean that observations $1,2,4$ belong
to the same cluster. Then the MFM prior can be expressed as

$\displaystyle K\sim p(\cdot),\quad(\pi_{1},\ldots,\pi_{K})\mid
K\sim\mbox{Dir}(\gamma,\ldots,\gamma),\quad z_{i}\mid
K,\pi\sim\sum_{h=1}^{K}\pi_{h}\delta_{h},\quad i=1,\ldots,n,$ (2)

where $K$ is the number of clusters, $(\pi_{1},\ldots,\pi_{K})$ are associated
cluster weights and $\sum_{h=1}^{K}\pi_{h}\delta_{h}$ is the mixture
distribution with $\delta_{h}$ being a point-mass at $h$. Under the Bayesian
framework, all those three quantities are random and hence are assigned with
prior distributions, i.e., we use $p(\cdot)$, which is a proper probability
mass function on $\mathbb{N}_{+}$, as a prior on $K$, and a Dirichlet
distribution on the mixture weights. Compared to the Chinese restaurant
process (CRP), the probability of introducing a new table (cluster) for MFM is
slowed down by a factor of $V_{n}(t+1)/V_{n}(t)$, which allows for a model-
based pruning of the tiny extraneous clusters. Here the coefficient $V_{n}(t)$
is defined as

$\displaystyle\begin{split}V_{n}(t)&=\sum_{n=1}^{+\infty}\dfrac{k_{(t)}}{(\gamma
k)^{(n)}}p(k),\end{split}$

where $k_{(t)}=k(k-1)\ldots(k-t+1)$, and $(\gamma k)^{(n)}={\gamma k}(\gamma
k+1)\ldots(\gamma k+n-1)$, and $\gamma$ is the hyperparameter in the Dirichlet
prior for the weights. The conditional distributions of $z_{i},i=2,\ldots,n$
under (2) can be defined in a Pólya urn scheme similar to CRP:

$\displaystyle P(z_{i}=c\mid
z_{1},\ldots,z_{i-1})\propto\begin{cases}\absolutevalue{c}+\gamma,&\text{at an
existing table labeled}\,c\\\
V_{n}(t+1)/V_{n}(t)\gamma,&\text{if}~{}\,c\,\text{~{}is a new
table}\end{cases},$ (3)

with $t$ being the number of existing clusters.

Now back to the shooting chart data, as we propose to use three independent
MFM priors for clustering shooting distance, angle, and game time, our final
model can be presented in the following hierarchical structure,

$\begin{split}&K_{\ell}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}p_{K},~{}~{}\ell=1,2,3,\\\
&\bm{\pi}_{\ell}=(\pi_{\ell,1},\ldots,\pi_{\ell,K_{\ell}})\mid
K_{\ell}\sim\text{Dir}(\nu,\ldots,\nu),~{}~{},\nu>0,~{}\ell=1,2,3,\\\
&\log{\bm{\gamma}_{\ell,1}},\ldots,\log{\bm{\gamma}_{\ell,K_{\ell}}}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}\text{MVN}_{p_{\ell}}(\bm{0},\Sigma_{\ell}),~{}~{}\ell=1,2,3,\\\
&\Sigma_{\ell}=\sigma^{2}_{\ell}(I_{\ell}-\rho_{\ell}\bm{W}_{\ell}),~{}~{}\ell=1,2,3,\\\
&\sigma^{2}_{\ell}\sim\text{Gamma}(a,b),~{}~{}\ell=1,2,3,\\\
&\rho_{\ell}\sim\text{Unif}(c_{1\ell},c_{2\ell}),~{}~{}\ell=1,2,3,\\\
&P(z_{i\ell}=j\mid\bm{\pi}_{\ell},K_{\ell})=\pi_{j\ell},~{}~{}\ell=1,2,3,~{}~{}j=1,\ldots,K_{\ell},\\\
&\bm{Y}_{i}\sim\text{Poisson}(\bm{\gamma}_{1,z_{i1}}\circ\bm{\gamma}_{2,z_{i2}}\circ\bm{\gamma}_{3,z_{i3}}),~{}i=1,\ldots,n,\end{split}$
(4)

where the main effects for distance, angle, and period are modeled by
multivariate normal distributions and their covariance matrices involve
adjacency matrices, denoted by $\bm{W}_{l}$’s, that are used for incorporating
the potential spatial and temporal correlation information (e.g., two shooting
locations are next to each other). To ensure those covariance matrices
$\Sigma_{\ell}$ are positive definite, we introduce $c_{1\ell}$ and
$c_{2\ell}$ as the reciprocals of minimum and maximum eigenvalues of
$\bm{W}_{l}$, respectively. For the prior $p_{K}$ on the number of clusters,
we consider a truncated Poisson$(1)$ following the recommendations in Miller
and Harrison (2018) and Geng et al. (2019).

Our multidimensional clustering model in (4) sits between two extremes. One is
the usual tensor clustering model that assumes the same cluster structure
across different directions, which certainly is more restrictive compared to
ours. The other is to marginally cluster over each of the tensor directions
and solve multiple clustering problems independently, which does not fully
utilize the tensor structural information. Our method combines the attractive
features from both sides by allowing cluster structures to be different over
directions while borrowing information to improve the estimation efficiency.

Our model in (4) can be viewed as a Bayesian mixture of rank-one tensor
models. Compared to the frequenstist work on tensor clustering (Sun and Li,
2019; Chi et al., 2020), where a Tucker decomposition is usually utilized and
the choice of the rank relies heavily on pre-specification or certain model
selection criteria, our approach is capable of automatically determining the
rank while quantifying the uncertainty in rank selection. Moreover, our method
is fully probabilistic; hence each mixture component is easy to interpret in
practice.

###  2.2 Theoretical Properties

Next we study the theoretical properties for the posterior distribution
obtained from model (4). For convenience, we define three mixing measures
$G_{\ell}=\sum_{i=1}^{K_{\ell}}\pi_{\ell,i}\delta(\gamma_{\ell,i})$ for
$\ell=1,2,3$, where $\delta(\cdot)$ is the point mass measure. In other words,
$G_{1},G_{2},G_{3}$ represent the clustering structures and associated
parameters along each of the three directions in the tensor. In order to
establish the posterior contraction results, we consider a refined parameter
space $\bm{\Theta^{*}}$ defined as
$\cup_{k_{1},k_{2},k_{3}=1}^{\infty}\bm{\Theta_{k}^{*}}$ for
$\bm{k}=(k_{1},k_{2},k_{3})$, where $\bm{\Theta_{k}^{*}}$ is a compact
parameter space for all model parameters including mixture weights and main
effects given a fixed cluster number for each direction, i.e.,
$K_{1}=k_{1},K_{2}=k_{2},$ and $K_{3}=k_{3}$. More precisely, we define
$\bm{\Theta_{k}^{*}}$ as

$\displaystyle\Big{\\{}$
$\displaystyle\pi_{\ell,i}\in(\epsilon,1-\epsilon)~{}~{}\text{for
every}~{}i=1,\ldots,k_{\ell},\ell=1,2,3,~{}\sum_{j}^{k_{\ell}}\pi_{\ell,j}=1~{}~{}\text{for
every}~{}\ell=1,2,3,$ $\displaystyle\gamma_{\ell,i}\in(-M,M)~{}~{}\text{for
every}~{}i=1,\ldots,k_{\ell},\ell=1,2,3,\Big{\\}},$

where $\epsilon$ and $M$ are some pre-specified positive constants. For any
two mixing measures $G_{1}=\sum_{i=1}^{k}p_{i}\delta(\gamma_{i})$ and
$G_{2}=\sum_{j=1}^{k^{\prime}}p_{j}^{\prime}\delta(\gamma_{j})$, we define
their Wasserstein distance as
$W(G_{1},G_{2})=\inf_{q\in\mathcal{Q}}\sum_{i,j}q_{ij}|\gamma_{i}-\gamma_{j}|$,
where $\mathcal{Q}$ denotes the collection of joint discrete distribution on
the space of $\\{1,\ldots,k\\}\times\\{1,\ldots,k^{\prime}\\}$ and $q_{ij}$ is
the probability being associated with $(i,j)$-element and it satisfies the
constraint that $\sum_{i=1}^{k}q_{ij}=p_{j}^{\prime}$ and
$\sum_{j=1}^{k^{\prime}}q_{ij}=p_{i}$, for every $i=1,\ldots,k$ and
$j=1,\ldots,k^{\prime}$.

For $\ell=1,2,3$, let $K_{\ell}^{0}$ and $G_{\ell}^{0}$ be the true number of
clusters and true mixing measure along direction $\ell$. Also let $P_{0}$ be
the associated joint probability measure. Then the following theorem
establishes the posterior consistency and contraction rate for the cluster
number and mixing measure. The proof is given in the Supplementary Material;
and it is based on the results for Bayesian mixture models in Guha et al.
(2021).

######  Theorem 1.

Let $\Pi_{n}(\cdot\mid\bm{Y}_{1},\ldots,\bm{Y}_{n})$ be the posterior
distribution obtained from (4) given i.i.d. observations
$\bm{Y}_{1},\ldots,\bm{Y}_{n}$. Suppose that the true parameters belong to
$\bm{\Theta^{*}}$. Then for each of $\ell=1,2,3$, we have

$\displaystyle\Pi_{n}\left\\{K_{\ell}=K_{\ell}^{0}\mid\bm{Y}_{1},\ldots,\bm{Y}_{n})\right\\}\rightarrow
1,~{}\text{and}~{}~{}\Pi_{n}\left\\{(W(G_{\ell},G_{\ell}^{0})\lesssim(\log
n/n)^{-1/2}\mid\bm{Y}_{1},\ldots,\bm{Y}_{n})\right\\}\rightarrow 1,$

almost surely under $P_{0}$ as $n\rightarrow\infty$.

Theorem 1 shows that as sample size $n\rightarrow\infty$, our proposed
Bayesian model is capable of correctly identifying the unknown number of
clusters along each of the tensor directions with posterior probability
tending to one. Moreover, the latent clustering structure (e.g., cluster
membership) can also be consistently recovered. The contraction rate for
$G_{\ell}$ is nearly parametric with an additional logarithmic factor, which
is common in the Bayesian asymptotic literature (Guha et al., 2021). The
assumption of a compact parameter space $\bm{\Theta^{*}}$ is needed to rule
out extreme scenarios, for example, when some mixture probabilities are
extremely close to $0$, for which it becomes very challenging to distinguish
between our model and a sub-model without these small mixture components. In
practice, this assumption is often satisfied since we can always restrict the
modeling parameters to take values within a pre-specified range, e.g.,
assuming cluster probability to be at least $\epsilon$ for some small
$\epsilon$ values such as $.0001\%$. Our results can also be extended to
general multi-way tensors as long as the independent MFM priors are used for
each direction.

###  2.3 Bayesian Inference

We discuss the posterior sampling scheme for our model. For the MFM prior, we
use the stick-breaking (Sethuraman, 1994) approximation to reconstruct

$K_{\ell}\sim
p_{K},\bm{\pi}_{\ell}=(\pi_{1\ell},\ldots,\pi_{K_{\ell}\ell})\mid
K_{\ell}\sim\text{Dir}(\nu,\ldots,\nu),~{}~{}\ell=1,2,3,$

as follows for each of $\ell=1,2,3$,

  * •

Step 1. Generate
$\eta_{1},\eta_{2},\cdots\overset{\text{iid}}{\sim}\text{Exp}(\psi_{\ell})$,

  * •

Step 2. Let $K_{\ell}=\min\\{j:\sum_{k=1}^{j}\eta_{k}\geq 1\\}$,

  * •

Step 3. Set $\pi_{h\ell}=\eta_{h}$, for $h=1,\cdots,K_{\ell}-1$,

  * •

Step 4. Set $\pi_{h\ell}=1-\sum_{h=1}^{K_{\ell}-1}\pi_{h}$,

where we choose $(K_{\ell}-1)\sim\mbox{Poisson}(\psi_{\ell})$ and $\nu=1$.
Based on the stick-breaking reparameterization, we obtain a similar
hierarchical model as the Dirichlet process mixture model in Ishwaran and
James (2001) when we choose a sufficiently large dimension $T$ for
$\bm{\pi}_{\ell}$ and set the last $T-K_{\ell}$ elements to be zero. Due to
the lack of available analytical form for the posterior distribution of
$\gamma$’s, we employ the MCMC sampling algorithm to sample from the posterior
distribution, and then obtain the posterior estimates of the unknown
parameters. Computation is facilitated by the nimble (de Valpine et al., 2017)
package in R (R Core Team, 2013).

To determine the final clustering configuration based on post-burn-in
iterations, we use the Dahl’s method (Dahl, 2006). The main idea is to obtain
a clustering configuration that best represents the posterior samples based on
comparing the “pairwise similarity” between different cluster structures. The
procedure can be described as follows. First, at MCMC iteration $t$, based on
the $n$-dimensional vector $(z_{1}^{(t)},\ldots,z_{n}^{(t)})$ for the latent
clusters, a membership matrix $M^{(t)}$ consisting of 0 and 1’s can be
obtained, where $M^{(t)}(i,j)=M^{(t)}(j,i)=1(z_{i}^{(t)}=z_{j}^{(t)})$. Next,
the membership matrices are averaged over all post-burn-in iterations to get a
matrix of pairwise similarity, $\overline{M}=\sum_{t=1}^{T}M^{(t)}/T$, where
$T$ denotes the total number of iterations. Finally, the iteration that has
the smallest element-wise Euclidean distance from $\overline{M}$ is taken as
the inferred clustering configuration, i.e., with $t^{*}$ being

$t^{*}=\operatorname*{arg\,min}_{t}\sum_{i=1}^{n}\sum_{j=1}^{n}\left(M^{(t)}(i,j)-\overline{M}(i,j)\right)^{2},$

and the final inferred configuration is obtaind as
$(z_{1}^{t^{*}},\ldots,z_{n}^{t^{*}})$.

##  3 Simulation

To evaluate the performance of the proposed model, simulation studies are
performed on generated data sets with a total of 150 players. We consider two
simulation settings. For the first setting, we consider a three-angle pattern
and three-distance partition of the court, i.e., the court is divided into $9$
parts based on combinations of distance and angle. Two clusters of size 75 are
set for angle, distance, and quarter, respectively. The patterns for angle and
group are visualized in Figure 8 of the Supplementary Material. For quarter
group 1, we choose $\bm{\gamma}_{3}=(-1,-1,-1,-1)^{\top}$; and for quarter
group 2, $\bm{\gamma}_{3}=(-0.5,-2,-0.5,-2)^{\top}$. For the second simulation
setting, we consider a finer partition of the court, including 11 angles, 12
distances, and again two quarter patterns in the same way as design 1. The
true number of clusters is 3 (each cluster cluster of size 50) for the angle,
3 (each cluster of size 50) for the distance, and 2 (each of size 75) for the
quarter. The angle and group patterns are visualized in Figure 9 of the
Supplementary Material. Under both settings, for each piece of the partitioned
court, the corresponding number of shots is generated using the associated
$\gamma_{1},\gamma_{2}$ and $\gamma_{3}$ based on the last line of Equation 4.
The proposed multidimensional clustering approach is then applied to fit the
generated data; and this procedure is repeated for 100 times for each setting.
All the computations were performed on a computing server (256GB RAM, with 8
AMD Opteron 6276 processors, operating at 2.3 GHz, with 8 processing cores in
each) and the running time was within twelve hours.

Figure 1: Simulation results: Rand index boxplots for angle, distance, and
quarter over 100 Monte-Carlo replicates.

To evaluate the clustering performance on each of the tensor directions, we
use the Rand index (RI; Rand, 1971), which is a commonly used metric that
measures the concordance between two clustering schemes. Taking values between
$0$ and $1$, a larger value of RI indicates a higher agreement. To evaluate
whether the true number of clusters is correctly inferred, we also examine the
total number of clusters inferred in each replicate over each of the three
directions.

We consider two competing methods, $K$-means (function kmeans() in R) and
density-based spatial clustering (DBSCAN; implemented in fpc, Hennig, 2020).
To make a fair comparison, we use the number of clusters obtained by our
method for $K$-means. For DBSCAN, as the method depends on a pre-specified
“reachability distance”, we use four candidate values, 25, 50, 75, and 100;
and we denote the methods as DBSCAN-25,$\ldots,$ DBSCAN-100 for the rest of
this paper. Both methods are applied to each of the three directions in an
independent manner. In other words, for 150 simulated players, we sum out the
distance and quarter directions, and obtain 150 11-dimensional count vectors
for clustering.

We summarize the rand indexes from 100 replicates as boxplots in Figure 1 and
also report the average RI in Table 1. From the results, we find a clear
advantage of our method over $K$-means. Compared to DBSCAN, our advantage is
not obvious under the simple setting, Design 1; but becomes significantly
better under Design 2. We also note that the performance of DBSCAN is quite
sensitive to the choice of the reachability distance, e.g., DBSCAN-25 has the
worst performance for all three directions under Design 2, but not for Design
1. Our method, on the other hand, manages to achieve a reasonably high average
Rand Index for different tensor directions under both simulation designs.
These results highlight the benefit of incorporating the tensor structure and
borrowing information from other directions by our method.

Figure 2: Simulation results: Histograms of cluster numbers over 100 Monte-Carlo replicates for angle, distance, and quarter under two designs. Correct estimate of $K$ is marked in orange color.  Table 1: Simulation results: Average Rand Index over 100 Monte-Carlo replicates under two simulation designs for the proposed method and two competing methods. Design | Method | Angle | Distance | Quarter
---|---|---|---|---
Design 1 | Proposed | 0.999 | 0.987 | 0.938
|  $K$-means | 0.720 | 0.712 | 0.985
| DBSCAN-25 | 0.996 | 0.998 | 0.982
| DBSCAN-50 | 1.000 | 1.000 | 1.000
| DBSCAN-75 | 1.000 | 1.000 | 1.000
| DBSCAN-100 | 1.000 | 1.000 | 1.000
Design 2 | Proposed | 0.826 | 0.851 | 0.903
|  $K$-means | 0.773 | 0.786 | 0.851
| DBSCAN-25 | 0.499 | 0.667 | 0.554
| DBSCAN-50 | 0.687 | 0.706 | 0.759
| DBSCAN-75 | 0.720 | 0.720 | 0.776
| DBSCAN-100 | 0.720 | 0.720 | 0.776

We also present the histogram of the estimated number of clusters from 100
replicates in Figure 2. Under Design 1, where there are only two clusters on
each direction, the proposed method manages to correctly estimate the cluster
number $98\%,94\%$, and $87\%$ of times for angle, distance, and quarter,
respectively. In Design 2, with a finer partition of the court, it becomes
harder to infer the number of clusters, and the percentage of correct
estimation reduces to $56\%,60\%$, and $61\%$.

##  4 NBA Data Analysis

We consider the shot attempts made by players during the 2017–2018 NBA regular
season excluding the overtime period. Rookie year players who started their
NBA career in 2017 are excluded. We also exclude players who made very few
number of shots in that season, e.g., due to long-term injury. Shots that were
made at negative degrees (under the polar coordinate system) are also
excluded. At the end, the dataset that we study consists of 122,001 shot
attempts made by 191 players, with Aron Baynes bottoming the list with 317,
and Russell Westbrook topping the list with 1356 shot attempts.

We consider the polar coordinate representation of shot attempts in a similar
way with Reich et al. (2006). We treat the basket as origin and partition the
angle (from 0 to $\pi$) into 11 equal sections. In terms of the shooting
distance, we partition it into 12 sections, with the first 11 be designed so
that the areas of all sectors and annular sectors are the same. The remaining
9 areas correspond to the remaining areas on the offensive half court. The
partition scheme is illustrated in Figure 7 of the Supplementary Material.
Compared to the partition scheme in Figure 2 of Reich et al. (2006), where the
annular sectors only covered the regions near the three-point line, we extend
the annular sectors because of the current trend of making three-point shots
among NBA players, e.g., Stephen Curry and Damian Lillard. For each player, we
further divide the number of shot attempts by four game quarters for each
court partition, and end up with a $11\times 12\times 4$-dimensional tensor.
In Figure 3, we choose three players, Bradley Beal, LeBron James, and Russel
Westbrook, and present their shot charts for demonstration. Some interesting
patterns can be observed from the plots, e.g., LeBron James makes more shots
facing the basket, and Russel Westbrook makes fewer shot attempts in the
fourth quarter on average.

Figure 3: Visualization of shot count tensors for Bradley Beal, LeBron James,
and Russel Westbrook.

We apply the proposed multidimensional heterogeneity learning approach on the
collected tensor data from 191 players. The same neighborhood matrices
$W_{1}$, $W_{2}$, and $W_{3}$ from the simulation studies are used. We
consider a MCMC chain of length 10,000 and and a thinning interval of 2,
resulting in a total of 5,000 posterior samples. We then discard the first
2,000 as burn-in and use the rest of 3,000 samples to obtain the final
clustering configuration using Dahl’s method as described in Section 2.3.

We obtain two clusters of sizes 71 and 120 over the angle direction as shown
in Figure 4. While it can be seen that players in both clusters make more
shots when facing the basket, those in cluster 1 also make a fair amount of
shots at the two wings, as well as the corners. Players in cluster 2, however,
mostly shoot in the region facing the basket and its immediate neighbors.
Compared to those in cluster 1, they make less corner shots, as the estimated
$\bm{\gamma}_{1}$ is almost 0 in the two regions on each side. Representative
players for the two clusters are, respectively, James Harden and John Wall.
Their shot charts are given in Figure 10 of the Supplementary Material.

Figure 4: Visualization of $\bm{\gamma}_{1}$ estimates for two shooting angle
clusters.

Shooting patterns in terms of distance to basket have been clustered into
three groups as visualized in Figure 5. Players in the cluster 1 have two hot
regions: near the basket, and beyond the three point line. Point guards and
shooting guards (small forwards) make the majority of this cluster (75
players), with representative players such as Kyrie Irving and Stephen Curry.
Compared with cluster 1, the 90 players in cluster 2 tend to shoot less beyond
the three point line, but make more perimeter shots. A representative player
for this cluster is Russell Westbrook. Finally, in cluster 3, most of the 26
players only shoot in regions that are closest to the basket, such as DeAndre
Jordan and Clint Capela. Most of their shots are slam dunks and alley-oops.
Some other players (e.g., Fred VanVleet and Jonas Valanciunas) in cluster 3,
although also making perimeter shots and three-pointers, rely heavily on lay-
ups. We pick one representative player from each cluster and present their
shooting charts in Figure 11 of the Supplementary Material.

Figure 5: Visualization of $\bm{\gamma}_{2}$ estimates for the three shooting
distance clusters.

Finally, the two clusters for quarters are visualized in Figure 6. In cluster
1, players make more shots in quarters 1 and 3, and less shots in quarters 2
and 4. Most players in this cluster are leading players in their teams, and
they often take breaks during the second quarter. In the fourth quarter,
leading players may also take breaks if their teams lead or fall behind by
wide margins. Stephen Curry, Kevin Durant and Paul George are in this cluster.
In cluster 2, the distribution of shots across four quarters is more even than
that in cluster 1, and on average the estimated $\bm{\gamma}_{3}$ is
relatively smaller. The cluster sizes are 91 and 100, which indicate these two
patterns are similarly prevalent among the players that we studied. We pick
Anthony Davis and Chris Paul as two representative players from two clusters
and present their shooting charts over four quarters in Figure 12 of the
Supplementary Material.

Figure 6: Visualization of $\bm{\gamma}_{3}$ estimates for two game quarter
clusters.

##  5 Discussion

We propose a new multidimensional tensor clustering method in this paper and
demonstrate its utility by studying how shooting patterns distribute over
court locations and game time among different players. Our method is
applicable to many other sports such as football and baseball, where it is
natural to formulate and model the multi-way array data. The proposed method
also applies to other applications such as imaging analysis and recommender
systems.

Several future work directions remain open. First, our method is based on
Poisson distribution for outcomes in the tensor; and it is of interest to
generalize this assumption by considering other types of distributions such as
zero-inflated Poisson and continuous distributions. Incorporating sparsity in
tensor models is another interesting direction that will allows us to deal
with high-dimensional tensors. From applications point of view, it is of
interest to analyze and compare the shooting patterns between different
periods of games/seasons, e.g., regular season versus playoff games, and
before-pandemic versus 2020 NBA bubble seasons.

## References

  * Albert et al. (2017) Albert, J., M. E. Glickman, T. B. Swartz, and R. H. Koning (2017).  Handbook of statistical methods and analyses in sports.  CRC Press.
  * Bi et al. (2021) Bi, X., X. Tang, Y. Yuan, Y. Zhang, and A. Qu (2021).  Tensors in statistics.  Annual review of statistics and its application 8, 345–368.
  * Cervone et al. (2014) Cervone, D., A. D’Amour, L. Bornn, and K. Goldsberry (2014).  Pointwise: Predicting points and valuing decisions in real time with nba optical tracking data.  In Proceedings of the 8th MIT Sloan Sports Analytics Conference, Boston, MA, USA, Volume 28, pp. 3.
  * Cervone et al. (2016) Cervone, D., A. D’Amour, L. Bornn, and K. Goldsberry (2016).  A multiresolution stochastic process model for predicting basketball possession outcomes.  Journal of the American Statistical Association 111(514), 585–599.
  * Chi et al. (2020) Chi, E. C., B. R. Gaines, W. W. Sun, H. Zhou, and J. Yang (2020).  Provable convex co-clustering of tensors.  Journal of Machine Learning Research 21(214), 1–58.
  * Dahl (2006) Dahl, D. B. (2006).  Model-based clustering for expression data via a Dirichlet process mixture model.  In M. V. Kim-Anh Do, Peter Müller (Ed.), Bayesian Inference for Gene Expression and Proteomics, Volume 4, pp. 201–218. Cambridge University Press.
  * de Valpine et al. (2017) de Valpine, P., D. Turek, C. J. Paciorek, C. Anderson-Bergman, D. T. Lang, and R. Bodik (2017).  Programming with models: writing statistical algorithms for general model structures with NIMBLE.  Journal of Computational and Graphical Statistics 26(2), 403–413.
  * Fernandez and Bornn (2018) Fernandez, J. and L. Bornn (2018).  Wide open spaces: A statistical technique for measuring space creation in professional soccer.  In Sloan Sports Analytics Conference, Volume 2018.
  * Franks et al. (2015) Franks, A., A. Miller, L. Bornn, K. Goldsberry, et al. (2015).  Characterizing the spatial structure of defensive skill in professional basketball.  The Annals of Applied Statistics 9(1), 94–121.
  * Geng et al. (2019) Geng, J., A. Bhattacharya, and D. Pati (2019).  Probabilistic community detection with unknown number of communities.  Journal of the American Statistical Association 114(526), 893–905.
  * Guha et al. (2021) Guha, A., N. Ho, and X. Nguyen (2021).  On posterior contraction of parameters and interpretability in bayesian mixture modeling.  Bernoulli 27(4), 2159–2188.
  * Guhaniyogi (2020) Guhaniyogi, R. (2020).  Bayesian methods for tensor regression.  Wiley StatsRef: Statistics Reference Online, 1–18.
  * Guhaniyogi et al. (2017) Guhaniyogi, R., S. Qamar, and D. B. Dunson (2017).  Bayesian tensor regression.  The Journal of Machine Learning Research 18(1), 2733–2763.
  * Hennig (2020) Hennig, C. (2020).  fpc: Flexible Procedures for Clustering.  R package version 2.2-9.
  * Hu et al. (2020) Hu, G., H.-C. Yang, and Y. Xue (2020).  Bayesian group learning for shot selection of professional basketball players.  Stat, e324.
  * Ishwaran and James (2001) Ishwaran, H. and L. F. James (2001).  Gibbs sampling methods for stick-breaking priors.  Journal of the American Statistical Association 96(453), 161–173.
  * Jiao et al. (2020) Jiao, J., G. Hu, and J. Yan (2020).  A Bayesian joint model for spatial point processes with application to basketball shot chart.  Journal of Quantitative Analysis in Sports.  Forthcoming.
  * Kolda and Bader (2009) Kolda, T. G. and B. W. Bader (2009).  Tensor decompositions and applications.  SIAM review 51(3), 455–500.
  * Mai et al. (2021) Mai, Q., X. Zhang, Y. Pan, and K. Deng (2021).  A doubly-enhanced EM algorithm for model-based tensor clustering.  Journal of the American Statistical Association (just-accepted), 1–44.
  * Miller et al. (2014) Miller, A., L. Bornn, R. Adams, and K. Goldsberry (2014).  Factorized point process intensities: A spatial analysis of professional basketball.  In E. P. Xing and T. Jebara (Eds.), Proceedings of the 31st International Conference on Machine Learning, Volume 32 of Proceedings of Machine Learning Research, Bejing, China, pp. 235–243. PMLR.
  * Miller and Harrison (2018) Miller, J. W. and M. T. Harrison (2018).  Mixture models with a prior on the number of components.  Journal of the American Statistical Association 113(521), 340–356.
  * R Core Team (2013) R Core Team (2013).  R: A Language and Environment for Statistical Computing.  Vienna, Austria: R Foundation for Statistical Computing.
  * Rand (1971) Rand, W. M. (1971).  Objective criteria for the evaluation of clustering methods.  Journal of the American Statistical Association 66(336), 846–850.
  * Reich et al. (2006) Reich, B. J., J. S. Hodges, B. P. Carlin, and A. M. Reich (2006).  A spatial analysis of basketball shot chart data.  The American Statistician 60(1), 3–12.
  * Sandholtz et al. (2020) Sandholtz, N., J. Mortensen, and L. Bornn (2020).  Measuring spatial allocative efficiency in basketball.  Journal of Quantitative Analysis in Sports 16(4), 271–289.
  * Sethuraman (1994) Sethuraman, J. (1994).  A constructive definition of Dirichlet priors.  Statistica Sinica, 639–650.
  * Spencer et al. (2019) Spencer, D., R. Guhaniyogi, and R. Prado (2019).  Bayesian mixed effect sparse tensor response regression model with joint estimation of activation and connectivity.  arXiv preprint arXiv:1904.00148.
  * Sun et al. (2014) Sun, W. W., B. Hao, and L. Li (2014).  Tensors in modern statistical learning.  Wiley StatsRef: Statistics Reference Online, 1–25.
  * Sun and Li (2019) Sun, W. W. and L. Li (2019).  Dynamic tensor clustering.  Journal of the American Statistical Association 114(528), 1894–1907.
  * Yin et al. (2020) Yin, F., G. Hu, and W. Shen (2020).  Analysis of professional basketball field goal attempts via a Bayesian matrix clustering approach.  arXiv preprint arXiv:2010.08495.
  * Yin et al. (2020) Yin, F., J. Jiao, G. Hu, and J. Yan (2020).  Bayesian nonparametric estimation for point processes with spatial homogeneity: A spatial analysis of NBA shot locations.  arXiv preprint arXiv:2011.11178.

## Supplemental Material

We present additional figures referenced in the main paper.

Figure 7: Illustration of the partition scheme imposed on the court. Figure 8:
Visualization for $\bm{\gamma}_{1}$ and $\bm{\gamma}_{2}$ in the first
simulation setting. Figure 9: Visualization for $\bm{\gamma}_{1}$ and
$\bm{\gamma}_{2}$ in the second simulation setting. Figure 10: Shot charts for
representative players from the two angle clusters. Figure 11: Shot charts for
representative players from the three distance clusters. Figure 12: Shot
charts for representative players from the two quarter clusters. Figure 13:
Visualization of
$\widehat{\gamma}_{1}+\widehat{\gamma}_{2}+\widehat{\gamma}_{3}$ for players
in different distance and angle groups and quarter group 1. Figure 14:
Visualization of
$\widehat{\gamma}_{1}+\widehat{\gamma}_{2}+\widehat{\gamma}_{3}$ for players
in different distance and angle groups and quarter group 2.


# Lattice Polaron in a Bose-Einstein Condensate of Hard-Core Bosons

Moroni Santiago-García  Instituto Nacional de Astrofísica, Óptica y
Electrónica, Calle Luis Enrique Erro No.1 Santa María Tonantzintla, Puebla CP
72840, Mexico  Shunashi G. Castillo-López  Arturo Camacho-Guardian
<EMAIL_ADDRESS>Instituto de Física, Universidad Nacional Autónoma de
México, Apartado Postal 20-364, Ciudad de México C.P. 01000, Mexico

###### Abstract

Lattice polarons, quasiparticles arising from the interaction between an
impurity and its surrounding bosonic environment confined to a lattice system,
have emerged as a platform for generating complex few-body states, probing
many-body phenomena, and addressing long-standing problems in physics. In this
study, we employ a variational ansatz to investigate the quasiparticle and
spectral properties of an impurity coupled to a condensate gas of hard-core
bosons in a two-dimensional optical lattice. Our findings demonstrate that the
polaron features can be tuned by adjusting the filling factor of the bath,
revealing intriguing polaron characteristics in the strongly interacting
regime. These results offer valuable insights for lattice polaron experiments
with ultracold gases and can serve as a guide for new experiments in emergent
quantum devices, such as moiré materials, where optical excitations can be
described in terms of hard-core bosons.

##  I Introduction

Ultracold gases have served as a robust platform for quantum simulation of
exotic many-body physics Bloch  _et al._ (2012); Georgescu  _et al._ (2014);
Gross and Bloch (2017); Schäfer  _et al._ (2020). This permitted the
realization of quantum analogs of well-understood phenomena Greiner  _et al._
(2002); Tarruell and Sanchez-Palencia (2018), as well as the exploration of
physics beyond the accessible regimes in condensed matter physics. Nowadays,
the versatility of these systems allows for addressing concise proposals for
long-standing open problems such as high-Tc superconductivity. Bose polaron
physics in ultracold gases has attracted much attention given this context,
and it has dramatically stimulated new theoretical and numerical approaches
Rath and Schmidt (2013); Ardila and Giorgini (2015); Shchadilova  _et al._
(2016); Grusdt  _et al._ (2018a); Christensen  _et al._ (2015); Christianen
_et al._ (2024); Levinsen  _et al._ (2021, 2015); Sun  _et al._ (2017); Sun
and Cui (2017); Christianen  _et al._ (2022a, b); Yoshida  _et al._ (2018);
Massignan  _et al._ (2021); Yegovtsev  _et al._ (2022); Skou  _et al._
(2021a); Cayla  _et al._ (2023); Nielsen and Ardila (2019); Guenther  _et al._
(2021); Bighin  _et al._ (2022); Camacho-Guardian (2023); Liu (2024); Isaule
(2024) to understanding the experimental realization of this phenomenon far
beyond its original formulation Jørgensen  _et al._ (2016); Hu  _et al._
(2016); Ardila  _et al._ (2019); Yan  _et al._ (2020); Skou  _et al._ (2021b,
2022); Morgen  _et al._ (2023); Etrych  _et al._ (2024).

In optical lattices, impurity physics has renewed interest in probing the
Mott-insulator to superfluid transition Colussi  _et al._ (2023), topological
phases Grusdt  _et al._ (2016); Camacho-Guardian  _et al._ (2019); Pimenov
_et al._ (2021), band geometry Pimenov (2024), magnetic polarons Koepsell  _et
al._ (2021); Grusdt  _et al._ (2018b); Wang  _et al._ (2021); Nielsen  _et
al._ (2021); Nielsen (2022), few-body physics Ding  _et al._ (2023), non-
equilibrium dynamics Isaule  _et al._ (2024), and polaron physics in strongly
correlated Fermi-Hubbard models Amelio and Goldman (2024). The study of
lattice polarons is further motivated by the advances in quantum gas
microscopy, which enables the imaging of individual atoms Bakr  _et al._
(2009, 2010), providing intricate spatial details of quantum states that
complement traditional spectroscopic information Koepsell  _et al._ (2021);
Hilker  _et al._ (2017).

On the other hand, recent developments with two-dimensional van der Waals
heterostructures place quantum gases into new territories where polaron
physics and Bose-Fermi mixtures are being realized with control and tunability
Ruiz-Tijerina (2023) with excitons and charge carriers (electrons and holes).
Polaron physics has demonstrated to be a powerful tool to sense correlated
phases of matter Shimazaki  _et al._ (2021); Smoleński  _et al._ (2021); Mazza
and Amaricci (2022); Amelio  _et al._ (2023); Julku  _et al._ (2023). In these
experiments, Bose-Fermi Hubbard systems may arise as a consequence of an
emergent moiré potential Kennes  _et al._ (2021); Huang  _et al._ (2022); Tang
_et al._ (2020). In multilayers, spatially indirect excitons (electron and
hole sitting in different layers) can arise, imprinting strong dipole-dipole
interactions, which can effectively lead to a Hubbard model sensitive to the
intrinsic quasi-bosonic character of the excitations Huang  _et al._ (2023).
Indeed, for relevant experiments, the excitons in a moiré superlattice may
behave as hard-core bosons, which have been predicted to exhibit superfluid
Remez and Cooper (2022) and supersolid phases Julku (2022).

Figure 1: (a) Schematic representation of the state describing a Bose-Einstein
condensate (BEC) of hard-core bosons in terms of the Bloch sphere, the angle
$\theta$ determines the filling factor of the BEC. (b) Cartoon of the impurity
(blue ball) coupled to the BEC of hard-core bosons (pink balls).

Motivated by the progress with ultracold gases and the relevance in the new
quantum materials, in this article, we study the strongly interacting impurity
in a Bose-Einstein condensate (BEC) of hard-core bosons in the atomic context.
For this purpose, we employ a variational ansatz to describe the spectral and
quasiparticle properties of an impurity embedded in a gas of hard-core bosons.
The character of the variational ansatz allows us to understand the coupling
between the collective excitations of the hard-core gas and the impurity
beyond the so-called Fröhlich-like Hamiltonian, which is only valid for weak
impurity-boson interactions. Our results show that the quasiparticle
properties and the emergence of the polaron branches can be tuned with the
underlying properties of the bath and find that for low densities, the
quasiparticle properties are governed by the beyond Fröhlich Hamiltonian.

Our results offer a benchmark for experiments with ultracold gases, where
interactions between the bath and impurities can be tuned on demand.
Furthermore, our study is motivated by recent experiments in van der Waals
heterostructures, where moiré excitons and electrons can be tightly confined
to a moiré superlattice and highly population Bose-Fermi mixture can be
created to realize polaron physics. Our theoretical framework may provide
valuable guidance for understanding these experimental systems.

##  II Model

Consider a single impurity coupled to a two-dimensional gas of bosons confined
in a square optical lattice of $N_{s}$ sites, as illustrated in Fig. 1. The
Hamiltonian of the system is given by

$\displaystyle\hat{H}$ $\displaystyle=$
$\displaystyle-\sum_{i=A,B}\sum_{\langle\mathbf{r},\mathbf{r}^{\prime}\rangle}\left[(t_{i}\hat{c}^{\dagger}_{\mathbf{r},i}\hat{c}_{\mathbf{r}^{\prime},i}+\text{h.c.})-\delta_{\mathbf{r},\mathbf{r}^{\prime}}\mu_{i}\hat{n}_{\mathbf{r},i}\right]$
(1) $\displaystyle+$
$\displaystyle\frac{U_{BB}}{2}\sum_{\mathbf{r}}\hat{c}^{\dagger}_{\mathbf{r},B}\hat{c}^{\dagger}_{\mathbf{r},B}\hat{c}_{\mathbf{r},B}\hat{c}_{\mathbf{r},B}$
$\displaystyle+$ $\displaystyle
U_{AB}\sum_{\mathbf{r}}\hat{n}_{\mathbf{r},B}\hat{n}_{\mathbf{r},A},$

where the creation (annihilation) operation of the A/B atoms is denoted by
$\hat{c}^{\dagger}_{\mathbf{r},A/B}$ ($\hat{c}_{\mathbf{r},A/B}$). We use
index $B$ for the two-dimensional Bose gas, while index $A$ corresponds to the
impurity. The tunneling coefficients are given by $t_{i}$, the chemical
potentials by $\mu_{i}$, and we consider an on-site boson-boson interaction,
$U_{BB}$, much larger than the tunneling, $t_{B}$, so we assume that the $B$
atoms can effectively be regarded as hard-core bosons. The impurity boson
interaction, $U_{AB}$, is also assumed to be short-range. Finally,
$\hat{n}_{\mathbf{r},i}=\hat{c}^{\dagger}_{\mathbf{r},i}\hat{c}_{\mathbf{r},i}$
corresponds to the number operators.

Considering the majority $B$ atoms as hard-core bosons implies explicitly
forbidding double occupancy of sites. To obtain the ground state of the $B$
gas and its collective excitations, we follow Refs. Bernardet  _et al._
(2002); Santiago-García and Camacho-Guardian (2023), and we only briefly
detail their proceeding. The idea is to account for the forbidden double
occupancy of the $B$ atoms by mapping the bosonic field operators into spin
1/2 operators; that is, we start by replacing
$\hat{c}^{\dagger}_{\mathbf{r},B}\rightarrow\hat{S}^{\dagger}_{\mathbf{r}},\
\text{and}\,\hat{c}_{\mathbf{r},B}\rightarrow\hat{S}_{\mathbf{r}}$, where
$\hat{S}^{\dagger}_{\mathbf{r}}$ ($\hat{S}_{\mathbf{r}}$) creates
(annihilates) a spin $+\frac{1}{2}$ ($-\frac{1}{2}$) at location
${\mathbf{r}}$. For $U_{AB}=0$, the ground state of the system can be written
as

$|\Phi_{0}\rangle=\hat{c}^{\dagger}_{\mathbf{k}=0,A}|\psi(\theta)\rangle_{\text{BEC}},$
(2)

this considers an impurity at the bottom of the band on top of the gas of
hard-core bosons, whose many-body wave-function is given by

$\displaystyle|\psi(\theta)\rangle_{\text{BEC}}=\prod_{\mathbf{r}}\left[\sin\left(\frac{\theta}{2}\right)+\cos\left(\frac{\theta}{2}\right)\hat{S}^{\dagger}_{\mathbf{r}}\right]|0\rangle.$
(3)

Here, the angle $\theta$ is determined by a variational approach Bernardet
_et al._ (2002): $\cos\theta=\mu_{B}/4t_{B}$, which defines the filling factor
$n_{B}=(\cos\theta+1)/2$.

To understand the coupling between the impurity and the majority bosons, we
write the Hamiltonian in Eq. (1) in terms of the collective excitations of the
hard-core BEC. In turn, in momentum space, the Hamiltonian of the system is
given by the sum of the following terms

$\hat{H}=\sum_{\mathbf{k}}[\omega(\mathbf{k})\hat{\gamma}_{\mathbf{k}}^{\dagger}\hat{\gamma}_{\mathbf{k}}+\epsilon_{A}(\mathbf{k})\hat{n}_{\mathbf{k},A}]+\hat{H}_{AB},$
(4)

where $\hat{\gamma}_{\mathbf{k}}^{\dagger}$ and $\hat{\gamma}_{\mathbf{k}}$
are the creation and annihilation bosonic operators describing the collective
excitations of the BEC, respectively, with
$\omega(\mathbf{k})=2\sqrt{\alpha^{2}_{\mathbf{k}}-\beta^{2}_{\mathbf{k}}}$
and

$\alpha_{\mathbf{k}}=\frac{1}{2}\left[\frac{\epsilon_{\mathbf{k}}}{2}(\cos^{2}\theta+1)+4t_{B}\right],$
(5a) $\beta_{\mathbf{k}}=-\frac{1}{4}\sin^{2}\theta\epsilon_{\mathbf{k}}.$
(5b)

Here, $\epsilon_{\mathbf{k}}=-2t_{B}[\cos{(k_{x}a)}+\cos{(k_{y}a})]$, where
$a$ is the lattice constant. The coupling between the impurity and the
excitations of the BEC is described by $\hat{H}_{AB}$ as follows

$\displaystyle\hat{H}_{AB}=\frac{U_{AB}}{2}(1+\cos{\theta})\sum_{\mathbf{r}}\hat{n}_{\mathbf{r},A}$
(6)
$\displaystyle-\frac{U_{AB}}{\sqrt{N_{s}}}\frac{\sin\theta}{2}\sum_{\mathbf{k,q}}(u_{\mathbf{q}}-v_{\mathbf{q}})(\hat{\gamma}^{\dagger}_{-\mathbf{q}}+\hat{\gamma}_{\mathbf{q}})\hat{c}^{\dagger}_{\mathbf{k}+\mathbf{q},A}\hat{c}_{\mathbf{k},A}$
$\displaystyle-\frac{1}{N_{s}}U_{AB}\cos{\theta}\sum_{\mathbf{k,q,k^{{}^{\prime}}}}\big{[}u_{\mathbf{k}^{{}^{\prime}}+\mathbf{q}}u_{\mathbf{k}^{{}^{\prime}}}\hat{\gamma}^{\dagger}_{\mathbf{k}^{{}^{\prime}}+\mathbf{q}}\hat{\gamma}_{\mathbf{k}^{{}^{\prime}}}$
$\displaystyle-
u_{\mathbf{k}^{{}^{\prime}}+\mathbf{q}}v_{\mathbf{k}^{{}^{\prime}}}\hat{\gamma}^{\dagger}_{\mathbf{k}^{{}^{\prime}}+\mathbf{q}}\hat{\gamma}^{\dagger}_{-\mathbf{k}^{{}^{\prime}}}-v_{\mathbf{k}^{{}^{\prime}}+\mathbf{q}}u_{\mathbf{k}^{{}^{\prime}}}\hat{\gamma}_{\mathbf{-}\mathbf{(}\mathbf{k}^{{}^{\prime}}+\mathbf{q})}\hat{\gamma}_{\mathbf{k}^{{}^{\prime}}}$
$\displaystyle+v_{\mathbf{k}^{{}^{\prime}}+\mathbf{q}}v_{\mathbf{k}^{{}^{\prime}}}\hat{\gamma}_{\mathbf{-}(\mathbf{k}^{{}^{\prime}}+\mathbf{q})}\hat{\gamma}^{\dagger}_{-\mathbf{k}^{{}^{\prime}}}\big{]}\hat{c}^{\dagger}_{\mathbf{k}-\mathbf{q},A}\hat{c}_{\mathbf{k},A},$

where the coherence factors of the BEC, $u_{\mathbf{k}}$ and $v_{\mathbf{k}}$,
are defined as

$u_{\mathbf{k}}=\sqrt{\frac{1}{2}+\frac{\alpha_{k}}{2\sqrt{\alpha_{k}^{2}-\beta_{k}^{2}}}},$
(7a)
$v_{\mathbf{k}}=\sqrt{-\frac{1}{2}+\frac{\alpha_{k}}{2\sqrt{\alpha_{k}^{2}-\beta_{k}^{2}}}}.$
(7b)

The impurity-boson Hamiltonian of Eq. (6) consists of three contributions: i)
The mean-field term in the first line, ii) the Fröhlich-like interaction where
a single collective excitation can be either absorbed or emitted (second
line), and iii) beyond Fröhlich contributions (third to fifth lines), which
involve processes with two collective modes. Notice that for
$\theta\rightarrow\pi$, the interaction Hamiltonian is completely governed by
the beyond Fröhlich term.

Next, we propose a Chevy-like ansatz Chevy (2006); Li and Das Sarma (2014);
Ding  _et al._ (2023) for the ground state as

$|\Psi\rangle=\left[\phi_{\mathbf{0}}\hat{c}^{\dagger}_{\mathbf{k}=0,A}+\sum_{\mathbf{q}}\phi_{\mathbf{q}}\hat{c}^{\dagger}_{\mathbf{q},A}\hat{\gamma}^{\dagger}_{\mathbf{-q}}\right]|\psi(\theta)\rangle_{\text{BEC}},$
(8)

with the variational parameters $\phi_{\mathbf{0}}$ and $\phi_{\mathbf{q}}$,
obtained from the variational principle
$\delta\langle\Psi|\hat{H}-E|\Psi\rangle/\delta\phi_{\mathbf{0}}^{*}=0$ and
$\delta\langle\Psi|\hat{H}-E|\Psi\rangle/\delta\phi_{{\mathbf{p}}}^{*}=0$,
respectively. This ansatz leads to the following set of equations:

$\epsilon_{MF}\phi_{\mathbf{0}}-\frac{U_{AB}}{2\sqrt{N_{s}}}\sin{\theta}\sum_{\mathbf{p}}\phi_{\mathbf{p}}(u_{\mathbf{p}}-v_{\mathbf{p}})=E\phi_{\mathbf{0}},$
(9a)
$\displaystyle\left[\omega(\mathbf{p})+\epsilon^{(A)}_{\mathbf{p}}+\epsilon_{MF}\right]\phi_{\mathbf{p}}-\frac{U_{AB}}{2\sqrt{N_{s}}}\sin{\theta}(u_{\mathbf{p}}-v_{\mathbf{p}})\phi_{\mathbf{0}}$
$\displaystyle-\frac{1}{N_{s}}U_{AB}\cos{\theta}\sum_{\mathbf{p^{{}^{\prime}}}}(u_{\mathbf{p}}u_{\mathbf{p^{{}^{\prime}}}}+v_{\mathbf{p}}v_{\mathbf{p^{{}^{\prime}}}})\phi_{\mathbf{p^{{}^{\prime}}}}=E\phi_{\mathbf{p}},$
(9b)

where $\epsilon_{MF}=n_{B}U_{AB}$ is the mean-field energy and
$\epsilon^{(A)}_{\mathbf{p}}=-2t_{A}[\cos{(k_{x}a)}+\cos{(k_{y}a})-2]$ is the
impurity dispersion.

For the two-dimensional Bose polaron in homogeneous confinement, Chevy’s
ansatz in Ref. Nakano  _et al._ (2024) has provided a very accurate
description compared to quantum Monte Carlo calculations Ardila  _et al._
(2020), and it is equivalent to the non-self-consistent T-matrix approximation
Cárdenas-Castillo and Camacho-Guardian (2022). Furthermore, the limitations of
Chevy’s ansatz for weakly interacting BEC are related to the formation of
clusters around the impurity Levinsen  _et al._ (2021), clusters that, due to
the bosonic nature of the atoms forming the BEC, can be arbitrarily large. The
formation of large polaron clouds cannot be correctly captured with this
ansatz. Here, however, we work in the opposite regime where these clusters are
intrinsically prevented due to the hard-core nature of the bosons. Therefore,
we expect this ansatz to be even more reliable with the hard-core constraint
for the bath bosons.

##  III Lattice polarons

We can now investigate the spectral features of the impurity coupled to the
hard-core gas and the properties of the resulting lattice polaron. In
particular, we perform calculations of the quasiparticle residue and the
spectral function. The polaron energy is simply given by $E$ since we set the
ground state energy of the BEC and the impurity to zero The quasiparticle
residue, $Z$, is determined by the squared overlap between the interacting
$|\Psi\rangle$ and the non-interacting $|\Phi_{0}\rangle$ states, the latter
corresponding to the case where the impurity sits on top of the BEC without
creating correlations, see Eq. (2):

$\displaystyle Z=|\langle\Phi_{0}|\Psi\rangle|^{2}.$ (10)

On the other hand, the spectral function of the impurity is defined by the
expression

$A(\mathbf{k}=0,\omega)=-\frac{1}{\pi}\mathbf{Im}\bigg{[}\sum_{n=0}\frac{|\langle\Phi_{0}|\Psi_{n}\rangle|^{2}}{\omega-
E_{n}+i\eta}\bigg{]},$ (11)

where $\\{|\Psi_{n}\rangle\\}$ is the set of eigenstates obtained from the
diagonalization of Eqs. (9a) and (9b), and $E_{n}$ denotes the corresponding
eigenenergies. For visibility purposes, we add a small imaginary number,
$i\eta$, to the spectral function. As we focus on the zero-momentum impurity,
we simply denote $A(\mathbf{k}=0,\omega)=A(\omega).$

Figure 2: Zero-momentum polaron energies as a function of the angle $\theta$
for $U_{AB}=-10.0$ (red line), $U_{AB}=-5.0$ (blue line), $U_{AB}=-1.0$ (black
line), $U_{AB}=1.0$ (orange line) and $U_{AB}=5.0$ (purple line).

First, we study the zero-momentum polaron energy, $E/t$, as a function of the
angle $\theta$ (which varies the filling factor $n_{B}$) for several values of
the interaction strength $U_{AB}.$ This is shown in Fig. 2, where we show for
attractive impurity-boson interactions: $U_{AB}=-10.0$ (red line),
$U_{AB}=-5.0$ (blue line), and $U_{AB}=-1.0$ (black line) as well as for
repulsive interactions $U_{AB}=1.0$ (dashed orange line) and $U_{AB}=5.0$
(dot-dashed purple line). Unless explicitly stated, all of our simulations are
performed for a squared two-dimensional lattice of $151\times 151$ sites. We
take $t_{A}=t_{B}=t.$

For attractive interactions in Fig. 2, we observe for small angles, i.e.,
close to unit filling factors, that the polaron energy simply saturates at
$U_{AB}$. Physically, this is very intuitive: when the bosons fill the
lattice, there is no longer a continuum of scattering processes, and the
energy of adding an impurity is simply given by the energy cost of putting a
boson and an impurity together, which is $U_{AB}.$ With increasing angle, that
is, decreasing the filling factor, scattering processes start to be allowed,
then, the energy of the polaron decreases smoothly with the angle $\theta.$
Interestingly, we find that when $\theta\rightarrow\pi$, that is,
$n_{B}\rightarrow 0$, the energy of the polaron does not tend to zero. This is
somewhat surprising because in the limit $n_{B}\rightarrow 0$, the problem can
be regarded as a $B$ atom placed in an empty lattice, and one would expect
that the energy of the polaron to approach zero as $n_{B}\rightarrow 0$.
However, this is a partial description and it is not entirely valid in the
strong coupling regime. To further understand this unexpected behavior, note
that in Eq. (6), the mean-field and Fröhlich terms of the Hamiltonian vanish
when $\theta\rightarrow\pi$, while the term beyond Fröhlich remains finite
(last three lines in Eq. (6)). This term increases with $U_{AB}$, leading to a
non-zero polaron energy for $n_{B}\rightarrow 0$. In this limit, both the
mean-field and Fröhlich contributions become negligible, and the polaron state
is completely dominated by the beyond Fröhlich term of the Hamiltonian, this
is a genuine result of strong interactions.

For repulsive interactions, observe that close to unit filling $(\theta=0)$,
the energy of the polaron saturates to $E\sim U_{AB}$, as shown in Fig. 2.
With increasing angle (decreasing filling factor), the polaron energy smoothly
decreases and tends to zero as the filling factor approaches zero. Our
findings suggest that for strong coupling, the repulsive and attractive
polaron behave differently in the limit $n_{B}\rightarrow 0$.

To understand the quasiparticle character of the polaron states discussed
above, we now turn to the polaron residue as defined in Eq. (10). In Fig. 3,
we show the residue $Z$ as a function of the angle $\theta$ and the
interaction strength $U_{AB}/t$.

Figure 3: Quasiparticle residue, $Z$, for the zero-momentum polaron as a
function of the interaction strength, $U_{AB}$, and the angle $\theta$.

For weak attractive interactions, we find that the polaron is well-defined
with $Z\sim 1$. As the interaction strength increases, the polaron remains a
well-defined quasiparticle for small angles (close to unit filling factors).
In contrast, for small filling factors ($\theta\rightarrow\pi$), the residue
of the polaron starts dropping to zero. That is, the states we discussed
previously, where the energy of the ground-state remains non-zero in the limit
$n_{B}\rightarrow 0$, have indeed very small residue and are consequently ill-
defined quasiparticles. Note that in this same limit ($n_{B}\rightarrow 0$),
the residue of the repulsive polaron is close to one, meaning that in this
case, this branch remains as a well-defined quasiparticle, contrary to the
attractive polaron. On the other hand, we obtain that close to unit filling
and for strong coupling, the attractive polaron remains well-defined whereas
the repulsive polaron cedes its spectral weight and is no longer a well-
defined quasiparticle. One should remark that the breakdown of the
quasiparticle picture occurrs only for strong interactions. For weak impurity-
boson interactions, we obtain that the quasiparticle picture holds valid for
all filling factors, including the limit cases of $n_{B}\rightarrow 0$ and
$n_{B}\rightarrow 1.$

Since our approach preserves the sum-rule $\int A(\omega)d\omega=1$, if the
polaron loses its residue, the spectral weight has to be distributed in high-
energy excitations. To understand how the spectral weight is transferred into
high-energy excitations, we plot in Fig. 4 the spectral function, $A(\omega)$,
as a function of the frequency $\omega/t$ for fixed $U_{AB}/t=-5$ and several
filling factor values.

Figure 4: Spectral function, $A(\omega)$, as a function of the frequency
$\omega/t$ with $U_{AB}=-5$ for $\theta=\pi/4$ (red line), $\theta=3\pi/4$
(blue line), $\theta=7\pi/8$ (black line) and $\theta=\pi-0.05$ (purple line).
With increasing angle (decreasing the filling factor), we observe that the
polaron becomes ill-defined and cedes its spectral weight to the continuum of
collective excitations.

Fig. 4 shows the evolution of the polaron peak with decreasing density
(increasing $\theta$) for strong impurity-boson interactions. For large
filling factors (small angles), the polaron peak retains most of the spectral
weight, as shown by the red line in Fig. 4; in this case, only the
quasiparticle peak at $\omega/t\approx U_{AB}=-5$ is visible in the spectral
function. By decreasing the filling factor (increasing $\theta$), we observe a
reduction of polaron peak, as indicated by the blue and black curves for
$\theta=3\pi/4$ and $\theta=7\pi/8$, respectively. As the polaron peak
shrinks, its spectral weight is transferred into high-energy excitations,
these excitations are incoherent and consist of the continuum of excitations
of our ansatz in Eq. (8); i.e., this continuum is formed by one collective
excitation of the BEC and an impurity state with opposite momentum. Close to
$\theta=\pi$, the polaron peak has disappeared and the quasiparticle residue
vanish, distributing its spectral weight in the continuum of excitations which
become more visible (purple line).

Due to the confinement of the BEC and the impurity within an optical lattice,
the continuum of excitations is bounded from above and below, and its width is
given by $W=4t[2+\sqrt{2(1+\cos^{2}{\theta})}]$, which depends on $\theta$.
Note that the fading of the polaron in the strong coupling regime is a result
that only the beyond Fröhlich term in the Hamiltonian of Eq. (6) can capture.

##  IV Conclusions

In this work, we have studied the problem of an impurity coupled to a Bose-
Einstein condensate of hard-core bosons in a two-dimensional optical lattice.
Employing a variational ansatz to describe the formation of the Bose polaron
for strong coupling, we unveiled the quasiparticle properties of the zero-
momentum impurity and have shown the interplay between the strong impurity-
boson interactions and the inherent collective excitations of the hard-core
bosons which tune the quasiparticle features.

The ability to probe lattice systems site-by-site using quantum microscopy has
increased interest in the study of lattice polarons, which may serve as probes
to measure quantum phase transitions Colussi  _et al._ (2023), topological and
geometric features Grusdt  _et al._ (2016); Camacho-Guardian  _et al._ (2019);
Pimenov  _et al._ (2021); Pimenov (2024), as well as following spatially non-
equilibrium dynamics Isaule  _et al._ (2024). Entering into the strongly
interacting regime may lead to the emergence of new few-body and many-body
states Ding  _et al._ (2023).

Also, the recent experiments with van der Waals materials where excitons in
moiré superlattices behave as hard-core bosons make our study relevant for
moiré polarons, moiré polaritons, and probing strongly correlated phases with
lattice polarons.

##  V Acknowledgments

We thank David Ruiz-Tijerina for valuable discussions and comments to our
manuscript. M. S. G acknowledges to Consejo Nacional de Humanidades, Ciencias
y Tecnología (CONAHCYT) for the scholarship provided. A.C.G. acknowledges
financial support from UNAM DGAPA PAPIIT Grant No. IA101923, and UNAM DGAPA
PAPIME Grants No. PE101223 and No. PIIF23. S.G.C-L acknowledges financial
support from UNAM DGAPA PAPIIT Grant No. TA100724.

Data availability statement.- The data that support the findings of this study
are openly available at the following URL/DOI:
https://doi.org/10.5281/zenodo.10841578

## References

  * Bloch  _et al._ (2012) Immanuel Bloch, Jean Dalibard, and Sylvain Nascimbene, “Quantum simulations with ultracold quantum gases,” Nature Physics 8, 267–276 (2012).
  * Georgescu  _et al._ (2014) I. M. Georgescu, S. Ashhab, and Franco Nori, “Quantum simulation,” Rev. Mod. Phys. 86, 153–185 (2014).
  * Gross and Bloch (2017) Christian Gross and Immanuel Bloch, “Quantum simulations with ultracold atoms in optical lattices,” Science 357, 995–1001 (2017).
  * Schäfer  _et al._ (2020) Florian Schäfer, Takeshi Fukuhara, Seiji Sugawa, Yosuke Takasu, and Yoshiro Takahashi, “Tools for quantum simulation with ultracold atoms in optical lattices,” Nature Reviews Physics 2, 411–425 (2020).
  * Greiner  _et al._ (2002) Markus Greiner, Olaf Mandel, Tilman Esslinger, Theodor W Hänsch, and Immanuel Bloch, “Quantum phase transition from a superfluid to a mott insulator in a gas of ultracold atoms,” nature 415, 39–44 (2002).
  * Tarruell and Sanchez-Palencia (2018) Leticia Tarruell and Laurent Sanchez-Palencia, “Quantum simulation of the hubbard model with ultracold fermions in optical lattices,” Comptes Rendus Physique 19, 365–393 (2018).
  * Rath and Schmidt (2013) Steffen Patrick Rath and Richard Schmidt, “Field-theoretical study of the bose polaron,” Physical Review A 88, 053632 (2013).
  * Ardila and Giorgini (2015) LA Peña Ardila and S Giorgini, “Impurity in a bose-einstein condensate: Study of the attractive and repulsive branch using quantum monte carlo methods,” Physical Review A 92, 033612 (2015).
  * Shchadilova  _et al._ (2016) Yulia E. Shchadilova, Richard Schmidt, Fabian Grusdt, and Eugene Demler, “Quantum dynamics of ultracold bose polarons,” Phys. Rev. Lett. 117, 113002 (2016).
  * Grusdt  _et al._ (2018a) Fabian Grusdt, Kushal Seetharam, Yulia Shchadilova, and Eugene Demler, “Strong-coupling bose polarons out of equilibrium: Dynamical renormalization-group approach,” Physical Review A 97, 033612 (2018a).
  * Christensen  _et al._ (2015) Rasmus Søgaard Christensen, Jesper Levinsen, and Georg M. Bruun, “Quasiparticle properties of a mobile impurity in a bose-einstein condensate,” Phys. Rev. Lett. 115, 160401 (2015).
  * Christianen  _et al._ (2024) Arthur Christianen, J. Ignacio Cirac, and Richard Schmidt, “Phase diagram for strong-coupling Bose polarons,” SciPost Phys. 16, 067 (2024).
  * Levinsen  _et al._ (2021) Jesper Levinsen, Luis A Peña Ardila, Shuhei M Yoshida, and Meera M Parish, “Quantum behavior of a heavy impurity strongly coupled to a bose gas,” Physical Review Letters 127, 033401 (2021).
  * Levinsen  _et al._ (2015) Jesper Levinsen, Meera M Parish, and Georg M Bruun, “Impurity in a bose-einstein condensate and the efimov effect,” Physical Review Letters 115, 125302 (2015).
  * Sun  _et al._ (2017) Mingyuan Sun, Hui Zhai, and Xiaoling Cui, “Visualizing the efimov correlation in bose polarons,” Phys. Rev. Lett. 119, 013401 (2017).
  * Sun and Cui (2017) Mingyuan Sun and Xiaoling Cui, “Enhancing the efimov correlation in bose polarons with large mass imbalance,” Phys. Rev. A 96, 022707 (2017).
  * Christianen  _et al._ (2022a) Arthur Christianen, J Ignacio Cirac, and Richard Schmidt, “Chemistry of a light impurity in a bose-einstein condensate,” Physical Review Letters 128, 183401 (2022a).
  * Christianen  _et al._ (2022b) Arthur Christianen, J Ignacio Cirac, and Richard Schmidt, “Bose polaron and the efimov effect: A gaussian-state approach,” Physical Review A 105, 053302 (2022b).
  * Yoshida  _et al._ (2018) Shuhei M. Yoshida, Shimpei Endo, Jesper Levinsen, and Meera M. Parish, “Universality of an impurity in a bose-einstein condensate,” Phys. Rev. X 8, 011024 (2018).
  * Massignan  _et al._ (2021) Pietro Massignan, Nikolay Yegovtsev, and Victor Gurarie, “Universal aspects of a strongly interacting impurity in a dilute bose condensate,” Phys. Rev. Lett. 126, 123403 (2021).
  * Yegovtsev  _et al._ (2022) Nikolay Yegovtsev, Pietro Massignan, and Victor Gurarie, “Strongly interacting impurities in a dilute bose condensate,” Phys. Rev. A 106, 033305 (2022).
  * Skou  _et al._ (2021a) Magnus G. Skou, Thomas G. Skov, Nils B. Jørgensen, and Jan J. Arlt, “Initial dynamics of quantum impurities in a bose–einstein condensate,” Atoms 9 (2021a), 10.3390/atoms9020022.
  * Cayla  _et al._ (2023) Hugo Cayla, Pietro Massignan, Thierry Giamarchi, Alain Aspect, Christoph I Westbrook, and David Clément, “Observation of 1/k 4-tails after expansion of bose-einstein condensates with impurities,” Physical Review Letters 130, 153401 (2023).
  * Nielsen and Ardila (2019) K Knakkergaard Nielsen and LA Peña Ardila, “Critical slowdown of non-equilibrium polaron dynamics,” New Journal of Physics 21, 043014 (2019).
  * Guenther  _et al._ (2021) Nils-Eric Guenther, Richard Schmidt, Georg M. Bruun, Victor Gurarie, and Pietro Massignan, “Mobile impurity in a bose-einstein condensate and the orthogonality catastrophe,” Phys. Rev. A 103, 013317 (2021).
  * Bighin  _et al._ (2022) G. Bighin, A. Burchianti, F. Minardi, and T. Macrì, “Impurity in a heteronuclear two-component bose mixture,” Phys. Rev. A 106, 023301 (2022).
  * Camacho-Guardian (2023) A. Camacho-Guardian, “Polaritons for testing the universality of an impurity in a bose-einstein condensate,” Phys. Rev. A 108, L021303 (2023).
  * Liu (2024) Ning Liu, “Weak and Strong Coupling Polarons in Binary Bose-Einstein Condensates,” arXiv e-prints , arXiv:2401.11808 (2024), arXiv:2401.11808 [cond-mat.quant-gas] .
  * Isaule (2024) Felipe Isaule, “Functional renormalisation group approach to the finite-temperature bose polaron,” arXiv preprint arXiv:2402.04197 (2024).
  * Jørgensen  _et al._ (2016) Nils B Jørgensen, Lars Wacker, Kristoffer T Skalmstang, Meera M Parish, Jesper Levinsen, Rasmus S Christensen, Georg M Bruun, and Jan J Arlt, “Observation of attractive and repulsive polarons in a bose-einstein condensate,” Physical review letters 117, 055302 (2016).
  * Hu  _et al._ (2016) Ming-Guang Hu, Michael J Van de Graaff, Dhruv Kedar, John P Corson, Eric A Cornell, and Deborah S Jin, “Bose polarons in the strongly interacting regime,” Physical review letters 117, 055301 (2016).
  * Ardila  _et al._ (2019) LA Peña Ardila, NB Jørgensen, T Pohl, S Giorgini, GM Bruun, and JJ Arlt, “Analyzing a bose polaron across resonant interactions,” Physical Review A 99, 063607 (2019).
  * Yan  _et al._ (2020) Zoe Z Yan, Yiqi Ni, Carsten Robens, and Martin W Zwierlein, “Bose polarons near quantum criticality,” Science 368, 190–194 (2020).
  * Skou  _et al._ (2021b) Magnus G Skou, Thomas G Skov, Nils B Jørgensen, Kristian K Nielsen, Arturo Camacho-Guardian, Thomas Pohl, Georg M Bruun, and Jan J Arlt, “Non-equilibrium quantum dynamics and formation of the bose polaron,” Nature Physics 17, 731–735 (2021b).
  * Skou  _et al._ (2022) Magnus G. Skou, Kristian K. Nielsen, Thomas G. Skov, Andreas M. Morgen, Nils B. Jørgensen, Arturo Camacho-Guardian, Thomas Pohl, Georg M. Bruun, and Jan J. Arlt, “Life and death of the bose polaron,” Phys. Rev. Res. 4, 043093 (2022).
  * Morgen  _et al._ (2023) AM Morgen, SS Balling, K Knakkergaard Nielsen, T Pohl, GM Bruun, and JJ Arlt, “Quantum beat spectroscopy of repulsive bose polarons,” arXiv (2023).
  * Etrych  _et al._ (2024) Jiří Etrych, Gevorg Martirosyan, Alec Cao, Christopher J Ho, Zoran Hadzibabic, and Christoph Eigen, “Universal quantum dynamics of bose polarons,” arXiv preprint arXiv:2402.14816 (2024).
  * Colussi  _et al._ (2023) VE Colussi, F Caleffi, C Menotti, and A Recati, “Lattice polarons across the superfluid to mott insulator transition,” Physical Review Letters 130, 173002 (2023).
  * Grusdt  _et al._ (2016) Fabian Grusdt, Norman Y Yao, D Abanin, Michael Fleischhauer, and E Demler, “Interferometric measurements of many-body topological invariants using mobile impurities,” Nature communications 7, 11994 (2016).
  * Camacho-Guardian  _et al._ (2019) A. Camacho-Guardian, N. Goldman, P. Massignan, and G. M. Bruun, “Dropping an impurity into a chern insulator: A polaron view on topological matter,” Phys. Rev. B 99, 081105 (2019).
  * Pimenov  _et al._ (2021) D. Pimenov, A. Camacho-Guardian, N. Goldman, P. Massignan, G. M. Bruun, and M. Goldstein, “Topological transport of mobile impurities,” Phys. Rev. B 103, 245106 (2021).
  * Pimenov (2024) Dimitri Pimenov, “Polaron spectra and edge singularities for correlated flat bands,” arXiv preprint arXiv:2401.15155 (2024).
  * Koepsell  _et al._ (2021) Joannis Koepsell, Dominik Bourgund, Pimonpan Sompet, Sarah Hirthe, Annabelle Bohrdt, Yao Wang, Fabian Grusdt, Eugene Demler, Guillaume Salomon, Christian Gross, _et al._ , “Microscopic evolution of doped mott insulators from polaronic metal to fermi liquid,” Science 374, 82–86 (2021).
  * Grusdt  _et al._ (2018b) F. Grusdt, M. Kánasz-Nagy, A. Bohrdt, C. S. Chiu, G. Ji, M. Greiner, D. Greif, and E. Demler, “Parton theory of magnetic polarons: Mesonic resonances and signatures in dynamics,” Phys. Rev. X 8, 011046 (2018b).
  * Wang  _et al._ (2021) Yao Wang, Annabelle Bohrdt, Shuhan Ding, Joannis Koepsell, Eugene Demler, and Fabian Grusdt, “Higher-order spin-hole correlations around a localized charge impurity,” Phys. Rev. Res. 3, 033204 (2021).
  * Nielsen  _et al._ (2021) K. K. Nielsen, M. A. Bastarrachea-Magnani, T. Pohl, and G. M. Bruun, “Spatial structure of magnetic polarons in strongly interacting antiferromagnets,” Phys. Rev. B 104, 155136 (2021).
  * Nielsen (2022) K. Knakkergaard Nielsen, “Exact nonequilibrium hole dynamics, magnetic polarons, and string excitations in antiferromagnetic bethe lattices,” Phys. Rev. B 106, 115144 (2022).
  * Ding  _et al._ (2023) Shanshan Ding, G. A. Domínguez-Castro, Aleksi Julku, Arturo Camacho-Guardian, and Georg M. Bruun, “Polarons and bipolarons in a two-dimensional square lattice,” SciPost Phys. 14, 143 (2023).
  * Isaule  _et al._ (2024) Felipe Isaule, Abel Rojo-Francàs, and Bruno Juliá-Díaz, “Bound impurities in a one-dimensional bose lattice gas: low-energy properties and quench-induced dynamics,” arXiv preprint arXiv:2402.03070 (2024).
  * Amelio and Goldman (2024) Ivan Amelio and Nathan Goldman, “Polaron spectroscopy of interacting fermi systems: Insights from exact diagonalization,” SciPost Physics 16, 056 (2024).
  * Bakr  _et al._ (2009) Waseem S Bakr, Jonathon I Gillen, Amy Peng, Simon Fölling, and Markus Greiner, “A quantum gas microscope for detecting single atoms in a hubbard-regime optical lattice,” Nature 462, 74–77 (2009).
  * Bakr  _et al._ (2010) Waseem S Bakr, Amy Peng, M Eric Tai, Ruichao Ma, Jonathan Simon, Jonathon I Gillen, Simon Foelling, Lode Pollet, and Markus Greiner, “Probing the superfluid–to–mott insulator transition at the single-atom level,” Science 329, 547–550 (2010).
  * Hilker  _et al._ (2017) Timon A Hilker, Guillaume Salomon, Fabian Grusdt, Ahmed Omran, Martin Boll, Eugene Demler, Immanuel Bloch, and Christian Gross, “Revealing hidden antiferromagnetic correlations in doped hubbard chains via string correlators,” Science 357, 484–487 (2017).
  * Ruiz-Tijerina (2023) David A Ruiz-Tijerina, “Bose–fermi mixtures in 2d solid-state superstructures,” Nature Materials 22, 153–154 (2023).
  * Shimazaki  _et al._ (2021) Yuya Shimazaki, Clemens Kuhlenkamp, Ido Schwartz, Tomasz Smoleński, Kenji Watanabe, Takashi Taniguchi, Martin Kroner, Richard Schmidt, Michael Knap, and Ata ç Imamoğlu, “Optical signatures of periodic charge distribution in a mott-like correlated insulator state,” Phys. Rev. X 11, 021027 (2021).
  * Smoleński  _et al._ (2021) Tomasz Smoleński, Pavel E Dolgirev, Clemens Kuhlenkamp, Alexander Popert, Yuya Shimazaki, Patrick Back, Xiaobo Lu, Martin Kroner, Kenji Watanabe, Takashi Taniguchi, _et al._ , “Signatures of wigner crystal of electrons in a monolayer semiconductor,” Nature 595, 53–57 (2021).
  * Mazza and Amaricci (2022) Giacomo Mazza and Adriano Amaricci, “Strongly correlated exciton-polarons in twisted homobilayer heterostructures,” Phys. Rev. B 106, L241104 (2022).
  * Amelio  _et al._ (2023) Ivan Amelio, N. D. Drummond, Eugene Demler, Richard Schmidt, and Atac Imamoglu, “Polaron spectroscopy of a bilayer excitonic insulator,” Phys. Rev. B 107, 155303 (2023).
  * Julku  _et al._ (2023) Aleksi Julku, Shanshan Ding, and Georg M Bruun, “An exciton interacting with a moire lattice–polarons, strings, and optical probing of spin correlations,” arXiv preprint arXiv:2311.03482 (2023).
  * Kennes  _et al._ (2021) Dante M Kennes, Martin Claassen, Lede Xian, Antoine Georges, Andrew J Millis, James Hone, Cory R Dean, DN Basov, Abhay N Pasupathy, and Angel Rubio, “Moiré heterostructures as a condensed-matter quantum simulator,” Nature Physics 17, 155–163 (2021).
  * Huang  _et al._ (2022) Di Huang, Junho Choi, Chih-Kang Shih, and Xiaoqin Li, “Excitons in semiconductor moiré superlattices,” Nature Nanotechnology 17, 227–238 (2022).
  * Tang  _et al._ (2020) Yanhao Tang, Lizhong Li, Tingxin Li, Yang Xu, Song Liu, Katayun Barmak, Kenji Watanabe, Takashi Taniguchi, Allan H MacDonald, Jie Shan, _et al._ , “Simulation of hubbard model physics in wse2/ws2 moiré superlattices,” Nature 579, 353–358 (2020).
  * Huang  _et al._ (2023) Tsung-Sheng Huang, Peter Lunts, and Mohammad Hafezi, “Non-bosonic moiré excitons,” arXiv preprint arXiv:2310.19931 (2023).
  * Remez and Cooper (2022) Benjamin Remez and Nigel R. Cooper, “Leaky exciton condensates in transition metal dichalcogenide moiré bilayers,” Phys. Rev. Res. 4, L022042 (2022).
  * Julku (2022) Aleksi Julku, “Nonlocal interactions and supersolidity of moiré excitons,” Physical Review B 106, 035406 (2022).
  * Bernardet  _et al._ (2002) K Bernardet, GG Batrouni, J-L Meunier, Guido Schmid, Matthias Troyer, and A Dorneich, “Analytical and numerical study of hardcore bosons in two dimensions,” Physical Review B 65, 104519 (2002).
  * Santiago-García and Camacho-Guardian (2023) Moroni Santiago-García and Arturo Camacho-Guardian, “Collective excitations of a bose–einstein condensate of hard-core bosons and their mediated interactions: from two-body bound states to mediated superfluidity,” New Journal of Physics 25, 093032 (2023).
  * Chevy (2006) F. Chevy, “Universal phase diagram of a strongly interacting fermi gas with unbalanced spin populations,” Phys. Rev. A 74, 063628 (2006).
  * Li and Das Sarma (2014) Weiran Li and S. Das Sarma, “Variational study of polarons in bose-einstein condensates,” Phys. Rev. A 90, 013618 (2014).
  * Nakano  _et al._ (2024) Yasufumi Nakano, Meera M. Parish, and Jesper Levinsen, “Variational approach to the two-dimensional bose polaron,” Phys. Rev. A 109, 013325 (2024).
  * Ardila  _et al._ (2020) L. A. Peña Ardila, G. E. Astrakharchik, and S. Giorgini, “Strong coupling bose polarons in a two-dimensional gas,” Phys. Rev. Res. 2, 023405 (2020).
  * Cárdenas-Castillo and Camacho-Guardian (2022) Luis Fernando Cárdenas-Castillo and Arturo Camacho-Guardian, “Strongly interacting bose polarons in two-dimensional atomic gases and quantum fluids of polaritons,” Atoms 11, 3 (2022).