\section{Introduction}
\label{Introduction}
 The study of congruences of Fourier coefficients of modular forms
 is an  interesting area of research in number theory. There were several important works in the literature on the congruences
 of Fourier coefficients (cf.~\cite{Rib84},~\cite{Hid85}), but we are particularly focusing on~\cite{CKO05}, where the authors studied
 the $p$-divisibility properties of Fourier coefficients of modular forms for $\mathrm{SL}_2(\mathbb{Z})$.  As an application, they have retrieved the results of Hatada, Hida on non-ordinary primes.
 Later, these results were generalized by El-Guindy
 for cusp forms of level $2,3,5,7 $ and $13$ (cf.~\cite{Gui07}).

The study of sign changes of Fourier coefficients of modular forms is also an active area of  research in modern number theory (cf.~\cite{Mur83},~\cite{KS06}). In~\cite{Sie69}, Siegel obtained an upper bound for the first sign change of Fourier coefficients of modular forms for $\mathrm{SL}_2(\mathbb{Z})$ by studying the space of linear relations among the initial Fourier coefficients of modular forms for $\mathrm{SL}_2(\mathbb{Z})$. A general result on the linear relations were
obtained by Choie, Kohnen, and Ono for $\mathrm{SL}_2(\mathbb{Z})$ (cf.~\cite{CKO05}), by El-Guindy for cusp forms of level $2,3,5,7 $ and $13$ (cf.~\cite{Gui07}).

 It is of real interest  to answer  similar questions  for the coefficients in the $u$-series expansion
 of Drinfeld modular forms. Such a study has been initiated
 by Choi. In~\cite{Cho08}, she studied ($T^q-T$)-divisibility properties of the coefficients  of Drinfeld modular forms of any weight, any type for $\mathrm{GL}_2(A)$, and determined all the linear relations between the initial coefficients  of Drinfeld modular forms of trivial type.
 In~\cite{Kaz08}, Kazalicki also obtained similar results for $\mathrm{GL}_2(A)$.

 For the level $\Gamma_0(T)$, these results are only known for Drinfeld modular forms of any weight, trivial type.
 In~\cite{Cho09}, Choi studied $(T^q-T)$-divisibility properties of the coefficient of  Drinfeld modular forms
 of any weight, trivial type for $\Gamma_0(T)$. In \textit{loc. cit.}, she determined all the linear relations between the initial coefficients in the $u$-series expansion of Drinfeld modular forms of trivial type, level $\Gamma_0(T)$.

In this article, we continue this study for $\Gamma_0(T)$ by generalizing the results of~\cite{Cho09} for
Drinfeld modular forms of any weight, \textbf{arbitrary type},  level $\Gamma_0(T)$.
Throughout the article, we fix to use the following notations:
\begin{itemize}
\item $p$ is an odd prime number and $q=p^r$ for some $r \in \mathbb{N}$.
\item $k\in \mathbb{N}$ and $l\in \mathbb{Z}/(q-1)\mathbb{Z}$ such that $k\equiv 2l \pmod {q-1}$.
      Let $0 \leq l \leq q-2$ be a lift of $l \in \mathbb{Z}/(q-1)\mathbb{Z}$. By abuse of notation,
      we continue to write $l$ for the integer as well as its class. Then, we define
      $r_{k, l} := \frac{k-2l}{q-1}$  and $r_{k, l,N} := r_{k,l} + N + 1$, where $N$ is a non-negative integer. Recall that $\dim M_{k,l}(\Gamma_0(T))=1+r_{k,l}$ (cf. \cite[Proposition 4.1]{DK}).
\end{itemize}
\subsection{An overview of the article }
The article is organized as follows. In \S \ref{Basic Theory} we recall the basic theory of Drinfeld modular forms. In \S \ref{Main_Section_1} we study congruences of coefficients and generalize \cite[Theorem 3.4]{Cho09} for Drinfeld modular forms of any weight, \textbf{arbitrary type},  level $\Gamma_0(T)$.
Finally, in \S \ref{Sec_Linear Relations} we determine all the  linear relations between the initial coefficients in $u$-series expansion of Drinfeld modular forms of any weight,  \textbf{arbitrary type} for  $\Gamma_0(T)$ which generalizes \cite[Theorem 4.1]{Cho09}.

\section{Basic theory of Drinfeld modular forms}
\label{Basic Theory}

In this section, we shall recall some basic theory of Drinfeld modular forms
(see \cite{Gos80},\cite{Gos80a}, \\ \cite{Gek88},\cite{GR96} for more details).

Let $\mathbb{F}_q$ denote the finite field of order $q$. Then, we set $A=\mathbb{F}_q[T]$
and
$K=\mathbb{F}_q(T)$.
Let $K_\infty=\mathbb{F}_q((\frac{1}{T}))$ be the completion of $K$
with respect to the infinite place $\infty$ (corresponding to $\frac{1}{T}$-adic valuation) and denote by $C$ the completion of an algebraic closure of $K_{\infty}$.
Let $L=\tilde{\pi}A \subseteq C$ be the $A$-lattice of rank $1$,
corresponding to the rank $1$ Drinfeld module given by $\rho_T=TX+X^q,$ where $\tilde{\pi}\in K_\infty(\sqrt[q-1]{-T})$ is defined up to a $(q-1)$-th root of unity.
The Drinfeld upper half-plane $\Omega= C-K_\infty$
has a rigid analytic structure.
The group $\mathrm{GL}_2(K_\infty)$ acts on $\Omega$ via fractional linear transformations.
  Any $x\in K_\infty^\times$ has the unique expression
$x= \zeta_x\big(\frac{1}{T} \big)^{v_\infty(x)}u_x,$
where $\zeta_x\in \mathbb{F}_q^\times$, and $v_\infty(u_x-1)> 0$ ($v_\infty$ is the valuation at $\infty$).
For any $k\in \mathbb{N},l\in \mathbb{Z}/(q-1)\mathbb{Z}, \gamma=\psmat{a}{b}{c}{d}\in \mathrm{GL}_2(K_{\infty})$, and $f:\Omega \longrightarrow C$,  we define
$f|_{k,l} \gamma := \zeta_{\det\gamma}^{l}\big(\frac{\det \gamma}{\zeta_{\det(\gamma)}} \big)^{k/2}(cz+d)^{-k}f(\gamma z).$
For an ideal $\mathfrak{n} \subseteq A$, let $\Gamma_0(\mathfrak{n})$  denote
the congruence subgroup $\Gamma_0(\mathfrak{n}) :=\{\psmat{a}{b}{c}{d}\in \mathrm{GL}_2(A) : c\in \mathfrak{n} \}.$
\begin{dfn}
A rigid holomorphic (resp., meromorphic) function $f:\Omega \longrightarrow C$ is said to be a holomorphic (resp., meromorphic) Drinfeld modular form for $\Gamma_0(\mathfrak{n})$
of weight $k$ and type $l$ if
\begin{enumerate}
\item $f|_{k,l}\gamma= f$ , $\forall \gamma\in \Gamma_0(\mathfrak{n})$,
\item $f$ is holomorphic (resp., meromorphic) at the cusps of $\Gamma_0(\mathfrak{n})$.
\end{enumerate}
\end{dfn}
Every holomorphic (resp., meromorphic) Drinfeld modular form has an unique expansion at the cusp $\infty$ w.r.t the parameter $u(z) :=  \frac{1}{e_L(\tilde{\pi}z)},$ where $e_L(z):= z{\prod_{\substack{0 \ne \lambda \in L }}}(1-\frac{z}{\lambda})$. Let $M_{k,l}(\Gamma_0(\mathfrak{n}))$ denote the $C$-vector space of holomorphic
Drinfeld modular forms of weight $k$ and type $l$ for $\Gamma_0(\mathfrak{n})$. If $f\in M_{k,l}(\Gamma_0(\mathfrak{n}))$ vanishes at the cusps of $\Gamma_0(\mathfrak{n})$, then we say that $f$ is a Drinfeld cusp form.
Any $f\in M_{k,l}(\Gamma_0(\mathfrak{n}))$ has $u$-series expansion at $\infty$ of the form
$\sum_{0 \leq \ i \equiv l \mod (q-1)}a_f(i)u^{i}$
(similarly, for meromorphic function $f$, the index $i$ starts from a negative integer).

By definition, if $k\not \equiv 2l \pmod {q-1}$ then $M_{k,l}(\Gamma_0(\mathfrak{n}))=\{0\}$. Hence, we always with the vector space $M_{k,l}(\Gamma_0(\mathfrak{n}))$ where $k\in \mathbb{N}$ and $l\in \mathbb{Z}/(q-1)\mathbb{Z}$ such that $k\equiv 2l \pmod {q-1}$.

We now give some useful examples of Drinfeld modular forms.
\begin{exple}[\cite{Gos80}, \cite{Gek88}]
\label{Eisenstein Series}
Let $d\in \mathbb{N}$. For $z\in \Omega$, the function
\begin{equation*}
g_d(z) := (-1)^{d+1}\tilde{\pi}^{1-q^d}L_d \sum_{\substack{a,b\in \mathbb{F}_q[T] \\ (a,b)\ne (0,0)}} \frac{1}{(az+b)^{q^d-1}}
\end{equation*}
is a Drinfeld modular form of weight $q^d-1$ and type $0$ for $\mathrm{GL}_2(A)$,
where $\tilde{\pi}$ is the Carlitz period and $L_d:=(T^q-T)\ldots(T^{q^d}-T)$ is the least common multiple of all monics of degree $d$.
We refer $g_d$ as an Eisenstein series and it does not vanish at $\infty$.
\end{exple}

\begin{exple}[Poincar\'e series]
\label{Poincare Series}
In~\cite{Gek88}, Gekeler defined the Poincar\'e series as follows:
$
 h(z) = \sum_{\gamma\in H\char`\\ \mathrm{GL}_2(A)} \frac{\det \gamma . u(\gamma z)}{(cz+d)^{q+1}},
$
where $H=\big\{\psmat{*}{*}{0}{1}\in \mathrm{GL}_2(A)\big\}$ and $\gamma = \psmat{a}{b}{c}{d}\in \mathrm{GL}_2(A)$.
Then $h$ is a cusp form of weight $q+1$, type $1$ for $\mathrm{GL}_2(A)$.
The $u$-series expansion of $h$ at $\infty$ is given by $ -u -u^{(q-1)^2+1}+\ldots$.
By the properties of $\Delta$-function in~\cite[Page 228]{Gek86a}, we deduce  that $h$ vanishes exactly once (resp., $q$-times) at $\infty$ (resp., at $0$) as a Drinfeld modular form of level $\Gamma_0(T)$.

\end{exple}

We end this section by introducing an important function $E$. In~\cite{Gek88}, Gekeler defined the function $E(z):= \frac{1}{\tilde{\pi}} \sum_{\substack{a\in \mathbb{F}_q[T] \\ a \ \mathrm{monic}}} ( \sum_{b\in \mathbb{F}_q[T]} \frac{a}{az+b} )$
 which is analogous to the Eisenstein series of weight $2$ over $\mathbb{Q}$. Though $E$ is not modular but we can use it to construct a Drinfeld modular form
\begin{equation}
\label{E_T}
E_T(z) := E(z)- TE(Tz) \in M_{2,1}(\Gamma_0(T)).
\end{equation}
The $u$-series expansion of $E_T$ at $\infty$ is given by $u-Tu^q+\ldots$
(cf.~\cite[Proposition 4.3]{DK21} for a detailed discussion about this function).

\section{Congruences for coefficients of Drinfeld Modular forms}
\label{Main_Section_1}

In this section, we generalize~\cite[Theorem 3.4]{Cho09} to Drinfeld modular forms $f$ of level $\Gamma_0(T)$, arbitrary type.
 We start by introducing the modular forms $\Delta_T$ and $\Delta_W$.
Recall that $0$ and $\infty$ are the only cusps of the Drinfeld modular curve $X_0(T):= \overline{\Gamma_0(T)\char`\\ \Omega}$ and the operator $W_T:=\psmat{0}{-1}{T}{0}$ permutes the cusps.
Consider the functions
$ \Delta_T(z) := \frac{g_1(Tz)-g_1(z)}{T^q-T}\ \mrm{and} \
\Delta_W(z) := \frac{T^qg_1(Tz)-Tg_1(z)}{T^q-T}=-T^{\frac{q+1}{2}}\Delta_T|_{q-1,0}W_T$.
Note that,  $\Delta_T$, $\Delta_W\in M_{q-1,0}(\Gamma_0(T))$ and their $u$-expansions are given by
$\Delta_T =u^{q-1}-u^{q(q-1)}+\ldots\in A[[u]]$ and
$\Delta_W = 1+ Tu^{q-1} - T^qu^{q(q-1)} + \ldots \in A[[u]]$.

\begin{prop}\cite[Proposition 4.3]{DK}
\label{IMP_Prop}
Let $\Delta_T, \Delta_W$ and $E_T$ be as defined before.
\begin{enumerate}
\item \label{P1}  \label{P2}
       The modular form $\Delta_T$(resp., $\Delta_W$) vanishes $q-1$ times at  $\infty$ (resp., at  $0$) and non-zero on $\Omega \cup \{0\}$ (resp.,  on $\Omega \cup \{\infty\}$). Hence,
       $\Delta_T$ and $\Delta_W$ are algebraically independent.
 \item \label{P3}
        The set $S:= \{\Delta_W^{r_{k,l}}E_T^l, \Delta_W^{r_{k,l}-1}\Delta_TE_T^l, \ldots, \Delta_W\Delta_T^{r_{k,l}-1} E_T^l, \Delta_T^{r_{k,l}}E_T^l \}$ forms a basis for the $C$-vector space $M_{k,l}(\Gamma_0(T))$.
 \item  \label{P4}
        For any $k \in \mathbb{N}$, the mapping $\eta: M_{k-2l,0}(\Gamma_0(T)) \longrightarrow M_{k,l}(\Gamma_0(T))$
         defined by $f \longrightarrow fE_T^l $
         is an isomorphism.
 \item \label{P5} \label{P6}
        We have $E_T^{q-1}= \Delta_W\Delta_T$.  In particular, the function $E_T\in M_{2,1}(\Gamma_0(T))$ vanishes exactly once at the cusps $0,\infty$ and non-vanishing elsewhere.
\item \label{P7} $h(z)= -\Delta_W(z)E_T(z)$.
\end{enumerate}
\end{prop}

In~\cite[Theorem 3.4]{Cho09}, Choi multiplied $f$ with $E$ to get the congruence for the coefficients of $fE$ modulo $T^{q}-T$, which imply  the congruences for the coefficients of $f$ (cf. Corollary~\cite[Corollary 3.5]{Cho09}). Although the methodology of the proof of the main result is similar, the novelty in our work is to multiply $f$ with $E_T^{q-l}$ to get the congruence for the coefficients of $fE_T^{q-l}$ modulo $T^{q^d}-T$ and deduce the congruences for the coefficients of $f$.

Before we state our main result of this section, let us recall some important results which are useful in the proof. First, we shall recall the following Theorem
(cf.~\cite[Theorem 7.14.2]{Har77}).

\begin{thm}[Residue Theorem]
\label{Residue Theorem}
Let $\mathfrak{n}$ be an ideal of $A$.
For any $1$-form $\omega$ on $X_0(\mathfrak{n}):= \overline{\Gamma_0(\mathfrak{n})\char`\\ \Omega}$, we have $\sum_{p\in X_0(\mathfrak{n})}\mathrm{Res}_p\omega=0.$
\end{thm}

\begin{prop}
\label{Residue at infinity}
Let $G$ be a meromorphic Drinfeld modular form of weight $2$, type $1$ for $\Gamma_0(T)$.
\begin{enumerate}
\item  \label{residue P1}
 If $G$ is holomorphic at $P\in \Gamma_0(T)\char`\\ \Omega$, then $\mathrm{Res}_P\ G(z)dz=0$.
\item \label{residue P2}
If the $u$-series expansion of $G$ at $\infty$ is given by
      $G(z)=\sum_{i\geq -n} a_{G,\infty} (i(q-1)+1)u^{i(q-1)+1},$ then
      $\mathrm{Res}_\infty G(z)dz=-\frac{a_{G,\infty} (1)}{\tilde{\pi}}$,
      where $\tilde{\pi}$ is the Carlitz period.
\item \label{residue P3}
If the $u_0$-series expansion of $G|_{2,1}\psmat{0}{-1}{1}{0}$ at $\infty$ is given by $\sum_{i\geq -n} a_{G,0}(i(q-1)+1)u_0^{i(q-1)+1}$ where $u_0:=u(z/T)$,
      then $\mathrm{Res}_0\ G(z)dz=-\frac{T}{\tilde{\pi}}a_{G,0}(1)$.
      Moreover, if the order of vanishing of $G$ at $0$ is at least $2$, then $\mathrm{Res}_0\ G(z)dz=0$.
\end{enumerate}
\end{prop}

\begin{proof}
\begin{enumerate}
\item By~\cite[\S 2.10]{GR96}, if $G$ is holomorphic at $P\in \Gamma_0(T)\char`\\ \Omega$,
      then $G(z)dz$ is also holomorphic at $P$. Hence $\mathrm{Res}_P\ G(z)dz=0$.
\item The parameter at $\infty$ is given by  $u=\frac{1}{e_L(\tilde{\pi}z)}$ which implies
      $dz=-\frac{1}{\tilde{\pi}} u^{-2}du$. Thus $G(z)dz = - \frac{1}{\tilde{\pi}} \cdot \sum_{i\geq -n} a_{G,\infty} (i(q-1)+1)u^{i(q-1)-1}du$. The coefficient of $u^{-1}$ gives the required result.
\item The parameter at $0$  is given by $u_0(z): =\frac{1}{e_L(\tilde{\pi}z/T)}$ which implies that
      $dz=-\frac{T}{\tilde{\pi}} u_0^{-2}du_0$. Thus
      $\big(G(z)|_{2,1}\psmat{0}{-1}{1}{0} \big)dz =- \frac{T}{\tilde{\pi}} \cdot \sum_{i\geq -n} a_{G,0}(i(q-1)+1)u_0^{i(q-1)-1}du_0$.
      Since the matrix $\psmat{0}{-1}{1}{0}$ permutes the cusps $0$ and $\infty$,
      by comparing the coefficient of $u_0^{-1}$ on both sides,
      we get  $\mathrm{Res}_0\ G(z)dz=-\frac{T}{\tilde{\pi}}a_{G,0}(1)$.
      If the order of vanishing of $G$ at $0$ is $\geq 2$, then $G(z)|_{2,1}\psmat{0}{-1}{1}{0} =\sum_{i\geq 1} a_{G,0}(i(q-1)+1)u_0^{i(q-1)+1}$, which implies $\mathrm{Res}_0\ G(z)dz=0$.
      \end{enumerate}
\end{proof}
Now, we are in a position to state and prove the main result of this section.
\begin{thm}
\label{MT_1}
Let $f\in M_{k,l}(\Gamma_0(T))$ be a non-zero Drinfeld modular form
such that the $u$-series expansion at $\infty$ belongs to $A[[u]]$.
Suppose the $u$-series expansion of $fE_T^{q-l}$  at $\infty$ is given by
$fE_T^{q-l}= \sum_{j\geq 0} a_{fE_T^{q-l}}(j(q-1)+1)u^{j(q-1)+1}$ where $E_T$ is as in~\eqref{E_T}.
Fix $d \in \mathbb{N}$. Let  $b \in \mathbb{N}$ such that $r_{k, l}+2+a\  \frac{q^d-1}{q-1}=p^b$ holds
for some integer $a \geq 0$.
Then the following congruence
$$a_{fE_T^{q-l}}(p^b(q-1)+1)\equiv 0 \pmod {T^{q^d}-T}$$
holds.
Moreover, if $a$ can chosen to be $0$, then $a_{fE_T^{q-l}}(p^b(q-1)+1)=0$.
\end{thm}

\begin{proof}

For any integer $a \geq 0$, consider the function $G(a):=\frac{g_d^a}{\Delta_T^{a.\frac{q^d-1}{q-1}}}.\frac{h}{\Delta_T^{r_{k, l}+1}}.\frac{f}{E_T^l}$.
\begin{itemize}
 \item  Since $g_d$ is non-vanishing at $\infty$ and $\Delta_T$ vanishes only at $\infty$,
        the function $\frac{g_d^a}{\Delta_T^{a.\frac{q^d-1}{q-1}}}$ is holomorphic on $\Omega \cup \{0\}$ and the possible pole is  only at $\infty$.
 \item  The function $h$ vanishes exactly once (resp., $q$-times) at $\infty$ (resp., at $0$) as
        a Drinfeld modular form for $\Gamma_0(T)$. The function $\Delta_T$
        vanishes only at $\infty$ and the order of vanishing at $\infty$ is $q-1$ (cf. Proposition~\ref{IMP_Prop}(\ref{P1})). This implies that the function $\frac{h}{\Delta_T^{r_{k, l}+1}}$ is holomorphic on $\Omega \cup \{0\}$ with order of vanishing at $0$ is $q$ and has a pole only at  $\infty$.
 \item  By Proposition \ref{IMP_Prop}(\ref{P4}) we get $\frac{f}{E_T^l}\in M_{k-2l,0}(\Gamma_0(T))$ is a holomorphic on $\Omega\cup \{0,\infty\}$.
\end{itemize}
Using the above properties, we get  the function $G(a)$ is a meromorphic Drinfeld modular form of weight $2$ and type $1$ for $\Gamma_0(T)$. Moreover, the function $G(a)$ is holomorphic on $\Omega\cup \{0\}$ with order of vanishing at least $q$ at $0$ and has a pole only at $\infty$. In particular,
by Proposition \ref{Residue at infinity}, we obtain
\begin{itemize}
 \item $\mathrm{Res}_P\ G(a)(z)dz=0$ \ for $P\in \Gamma_0(T)\char`\\ \Omega$, $\mathrm{Res}_0\ G(a)(z)dz=0$,
 \item $\mathrm{Res}_\infty\ G(a)(z)dz=-\frac{a_{G(a)}(1)}{\tilde{\pi}}$,
  where $a_{G(a)}(1)$ is the coefficient of $u$ in the $u$-series expansion of $G$ at $\infty$.
\end{itemize}
By Theorem~\ref{Residue Theorem}, we get  $\mathrm{Res}_\infty G(a)(z)dz= 0$
and hence $a_{G(a)}(1)=0$.
By Proposition~\ref{IMP_Prop}((\ref{P5})\&(\ref{P7})) we have $h= -\frac{E_T^q}{\Delta_T}$. Consequently,
$G(a) = -\frac{g_d^aE_T^{q-l}f}{\Delta_T^{r_{k, l}+a.\frac{q^d-1}{q-1}+2}}.$
Since $g_d\equiv 1 \pmod {T^{q^d}-T}$, we get
 \begin{equation}
 \label{congruence mod product of degree d}
-G(a) = \frac{g_d^aE_T^{q-l}f}{\Delta_T^{r_{k, l}+a.\frac{q^d-1}{q-1}+2}} \equiv
      \frac{E_T^{q-l}f}{\Delta_T^{r_{k, l}+a.\frac{q^d-1}{q-1}+2}}\pmod {T^{q^d}-T}.
\end{equation}
Now, we  begin the proof of the theorem.   By hypothesis, let $b\in \mathbb{N}$  such that
$r_{k, l}+a.\frac{q^d-1}{q-1}+2 =p^b $
holds for some $a \geq 0$. With this choice of $a$, ~\eqref{congruence mod product of degree d} becomes
\begin{equation}
 \label{derived congruence mod product of degree d}
-G(a) =\frac{g_d^aE_T^{q-l}f}{\Delta_T^{p^b}} \equiv
      \frac{E_T^{q-l}f}{\Delta_T^{p^b}}\pmod {T^{q^d}-T}.
\end{equation}
Since $\frac{1}{\Delta_T^{p^b}}= u^{-p^b(q-1)}(1-u^{p^b(q-1)^2}+\ldots)$, an easy computation shows that the coefficient of $u$ in the $u$-series expansion of $\frac{E_T^{q-l}f}{\Delta_T^{p^b}}$
at  $\infty$ is $a_{fE_T^{q-l}}(p^b(q-1)+1)$.  By~\eqref{derived congruence mod product of degree d}, we obtain
$a_{fE_T^{q-l}}(p^b(q-1)+1) \equiv - \pmod {T^{q^d}-T}.$
Since $a_{G(a)}(1)=0$, we are done.
\end{proof}
\begin{remark}
When $d=1$, the existence of $a\geq 0$ with $r_{k,l}+2+a=p^b$ is automatic if $b \gg 0$.
Theorem~\ref{MT_1} is a generalization of~\cite[Theorem 3.4]{Cho09}
from trivial type to arbitrary type.
\end{remark}

Now, we have a corollary of Theorem~\ref{MT_1}.
\begin{cor}
\label{Choi_Corollary}
Let $f\in M_{k,l}(\Gamma_0(T))$ be a Drinfeld modular form such that the $u$-series expansion at $\infty$ is given by
$f=\sum_{j\geq 0} a_f(j(q-1)+l)u^{j(q-1)+l}\in A[[u]]$. If $p^\alpha \mid l$ for some $\alpha \in \mathbb{N}$ and
 $m \leq \alpha$ be a natural number such that $p^m>r_{k, l}+1$, then the congruence
$$a_f((p^m-1)(q-1)+l)\equiv 0 \pmod {T^q-T} \ holds.$$
\end{cor}

\begin{remark}
Corollary~\cite[Corollary 3.5]{Cho09} is a special case of Corollary~\ref{Choi_Corollary} for $l=0$ and $m=1$. In fact,  the assumption $q > p+1$ in Corollary $\it{loc. cit.}$ is redundant.
\end{remark}
\begin{proof}[Proof of Corollary~\ref{Choi_Corollary}]
Let $q-l=p^\alpha x$ for some $x\in \mathbb{N}$. Now
\begin{align*}
&\sum_{n\geq 0} a_{fE_T^{p^\alpha x}}( n(q-1)+1)u^{n(q-1)+1}=fE_T^{p^\alpha x}\\
&= \sum_{j\geq 0} a_f(j(q-1)+l)u^{j(q-1)+l}. \Big(\sum_{i\geq 0} (a_{E_T}(i(q-1)+1))^{p^\alpha}u^{p^\alpha i(q-1)+p^\alpha}\Big)^x.
\end{align*}
By comparing the coefficients of $u^{p^m(q-1)+1}$ on both sides, we get
$$a_{fE_T^{q-l}}(p^m(q-1)+1)= a_f((p^m-1)(q-1)+l)a_{E_T}(1)^{p^\alpha x}=a_f((p^m-1)(q-1)+l),$$
since $a_{E_T}(1)=1$. The first equality follows from the fact that $p^\alpha (q-1)+p^\alpha>p^m(q-1)+1$.
Now, the result follows from Theorem~\ref{MT_1} by taking $d=1,b=m$.
\end{proof}

Now, we shall give an example satisfying Corollary~\ref{Choi_Corollary}.
\begin{exple}
Let $p=3, q=3^2, l=6$ and $k=12$, the assumptions of Corollary \ref{Choi_Corollary} are satisfied and the space $M_{12,6}(\Gamma_0(T))$ is generated by $E_T^6$.
Hence, by Corollary~\ref{Choi_Corollary}, we get that
$a_{E_T^6}((p-1)(q-1)+6) \equiv 0\ (T^q-T).$ In fact, a straight forward calculation shows that  $a_{E_T^{pm}}((p-1)(q-1)+pm)=0$ for any $m\in \mathbb{N}$ with $pm<q-1$.
\end{exple}
We can also produce infinitely many examples satisfying Corollary~\ref{Choi_Corollary} as follows:
\begin{exple}
Suppose $p|l$. Consider the function $g_1^{p-2}E_T^l\in M_{(q-1)(p-2)+2l,l}(\Gamma_0(T))$.
Since $p> r_{(q-1)(p-2)+2l, l}+1$,  we get $a_{g_1^{p-2}E_T^l}((p-1)(q-1)+l)\equiv 0 \pmod {T^q-T}$
by Corollary~\ref{Choi_Corollary}.
\end{exple}

\section{linear relations between the initial Fourier coefficients}
\label{Sec_Linear Relations}
In this section, we determine all the  linear relations between the initial coefficients in $u$-series expansion of Drinfeld modular forms of any weight,  arbitrary type for  $\Gamma_0(T)$.
For any integer $i \geq 0$, consider the map
$$a_i^*: M_{k,l}(\Gamma_0(T)) \rightarrow C \  \ \mathrm{defined \ by} \
f \rightarrow a_f(i(q-1)+l),$$
where $f=\sum_{0\leq j\equiv l \pmod {q-1}}^\infty a_f(j)u^{j}.$
Then $a_i^* \in (M_{k,l}(\Gamma_0(T)))^*$, the dual of $M_{k,l}(\Gamma_0(T))$. The elements of $S$ in Proposition~\ref{IMP_Prop}(\ref{P3}) satisfy $a_i^*(\Delta_W^{r_{k,l}-j}\Delta_T^j E_T^l)=0$ for $i<j$, which implies
\begin{prop}
\label{Basis for the dual}
The set $\{a_i^*\}_{i=0}^{r_{k,l}}$ forms a basis for the $C$-vector space
$(M_{k,l}(\Gamma_0(T)))^*$.
\end{prop}

For any integer $N\geq 0$, consider the surjective map
$$\psi_{k,l;N} : C^{r_{k, l,N}+1} \longrightarrow (M_{k,l}(\Gamma_0(T)))^* \ \ \mathrm{defined \ by} \
(c_0,c_1,\ldots,c_{r_{k, l,N}}) \longrightarrow \sum_{i=0}^{r_{k, l,N}}c_i a_i^*.$$
\begin{dfn}
\label{definition of L_k,l;N(Gamma_0T)}
For $k \in \mathbb{N}$, $l\in \mathbb{Z}/(q-1)\mathbb{Z}$ such that $k\equiv 2l \pmod {q-1}$, we define
$L_{k,l;N}(\Gamma_0(T)):= \mrm{Ker}(\psi_{k,l;N})$. By Proposition~\ref{Basis for the dual}, $\dim_C L_{k,l;N}(\Gamma_0(T))= N+1$.
\end{dfn}

Now, we state the main result of this section, which is a generalization of~\cite[Theorem 4.1]{Cho09} from trivial type to arbitrary type.  For each Drinfeld modular form $g\in M_{N(q-1)+2l,l}(\Gamma_0(T))$, we define the elements $b(k,l,N,g;i)$, $c(k,l,N,g;i)\in C$ by the $u$-series expansion
\begin{equation}
\label{required function for making weight 2 type 1 form}
\frac{hg}{\Delta_T^{r_{k, l,N}}E_T^{2l}}=\sum_{i=0}^{r_{k, l,N}}b(k,l,N,g;i)u^{-i(q-1)+1-l}+ \sum_{i=1}^\infty c(k,l,N,g;i)u^{i(q-1)+1-l}.
\end{equation}

\begin{thm}
\label{MT_2}
The map $\phi_{k,l;N}: M_{N(q-1)+2l,l}(\Gamma_0(T)) \longrightarrow L_{k,l;N}(\Gamma_0(T))$ defined by
$\phi_{k,l;N}(g)  = (b(k,l,N,g;0),b(k,l,N,g;1),\ldots, b(k,l,N,g;r_{k, l,N}) )$
is an isomorphism of $C$-vector spaces.
\end{thm}

\begin{proof}
First we  show that the image of $\phi_{k,l;N}$ belongs to $L_{k,l;N}(\Gamma_0(T))$.
Let $g\in M_{N(q-1)+2l,l}(\Gamma_0(T))$.  Then, for any $f=\sum_{j=0}^\infty a_f(j(q-1)+l)u^{j(q-1)+l} \in M_{k,l}(\Gamma_0(T))$, consider the function
$G:= \frac{hg}{\Delta_T^{r_{k, l,N}}E_T^{2l}} \cdot f=
\frac{h}{\Delta_T^{r_{k, l,N}}}\cdot \frac{g}{E_T^l} \cdot \frac{f}{E_T^l}$.

\begin{itemize}
 \item The function $h$ vanishes exactly once (resp., $q$-times) at $\infty$ (resp., at $0$) and $\Delta_T$ vanishes only at $\infty$ and the order of vanishing is $q-1$
       (Proposition~\ref{IMP_Prop}(\ref{P1})). This implies that $\frac{h}{\Delta_T^{r_{k, l,N}}}$ is holomorphic on $\Omega \cup \{0\}$ with order of vanishing $q$ at $0$ and a pole only at $\infty$.
 \item By Proposition~\ref{IMP_Prop}(\ref{P4}), $\frac{g}{E_T^l}\in M_{N(q-1),0}(\Gamma_0(T))$ and $\frac{f}{E_T^l}\in M_{k-2l,0}(\Gamma_0(T))$.

\end{itemize}
These properties imply that the function $G$ is a meromorphic Drinfeld modular form of weight $2$ and type $1$ for $\Gamma_0(T)$. Moreover, the function $G$ is holomorphic on $\Omega\cup\{0\}$ with order of vanishing at least $q$ at $0$ and has a pole only at $\infty$. Note that, the coefficient of $u$ in the $u$-series expansion of $G$ is $\sum_{i=0}^{r_{k, l,N}}b(k,l,N,g;i)a_f(i(q-1)+l)$, which is equal to $0$
by Theorem~\ref{Residue Theorem} and Proposition~\ref{Residue at infinity}.
Thus the image of $\phi_{k,l;N}$ belongs to $L_{k,l;N}(\Gamma_0(T))$.
Clearly, $\phi_{k,l;N}$  is linear.

Now, we show that $\phi_{k,l;N}$ is injective. Suppose that $\phi_{k,l;N}(g)=0$ for some $g\in M_{N(q-1)+2l,l}(\Gamma_0(T))$. For any $f\in M_{k,l}(\Gamma_0(T))$, by~\eqref{required function for making weight 2 type 1 form}, we get
\begin{equation}
\frac{hg}{\Delta_T^{r_{k, l,N}}E_T^{2l}}.f= \sum_{i=1}^\infty c(k,l,N,g;i)u^{i(q-1)+1-l}\cdot\sum_{j=0}^\infty a_f(j(q-1)+l)u^{j(q-1)+l}.
\end{equation}
Therefore $\frac{hgf}{\Delta_T^{r_{k, l,N}}E_T^{2l}}$ is a doubly cuspidal holomorphic Drinfeld modular form of weight $2$, type $1$ for $\Gamma_0(T)$. Since the modular curve $X_0(T)$ has genus $0$, the function
$\frac{hgf}{\Delta_T^{r_{k, l,N}}E_T^{2l}}$ is identically zero.  This forces $g = 0$.
Hence the map $\phi_{k,l;N}$ is injective. Since the dimensions of
$L_{k,l;N}(\Gamma_0(T))$, $M_{N(q-1)+2l,l}(\Gamma_0(T))$ are equal,
the map $\phi_{k,l;N}$ is an isomorphism.

\end{proof}
The Theorem above and  Proposition~\ref{IMP_Prop}((\ref{P3}) and (\ref{P4})), we get
\begin{cor}
The $C$-vector spaces $L_{k,l;N}(\Gamma_0(T))$ and $L_{k-2l,0;N}(\Gamma_0(T))$ are isomorphic.
\end{cor}
We conclude the article with a remark that one can  extend \cite[Theorem 4.1]{Cho08} for $\mathrm{GL}_2(A)$ from trivial type to arbitrary type by replacing the function $\frac{g_1^{q-\alpha}}{h^{(q-1)(r+N-1)}}u$ in \cite[Page 99]{Cho08} by $\frac{g_1^{q-\beta}}{h^{(q-1)(\dim_C  M_{k,l}(\mathrm{GL}_2(A))+N)+2l-1}}f$,  where $f\in M_{N(q^2-1)+l(q+1),l}(\mathrm{GL}_2(A))$ and $\beta = \frac{k-(q+1)l}{q-1}+(1-\dim_C M_{k,l}(\mathrm{GL}_2(A)))(q+1)$.

\bibliographystyle{plain, abbrv}


\section{Introduction}

A hot cup of tea, left on the counter, will eventually cool to room temperature, losing all information encoded in its initial conditions. A similar process of thermalisation is also expected in generic isolated quantum systems. At present, this is best understood through the eigenstate thermalisation hypothesis (ETH) \cite{DAl-16}, which conjectures that an isolated quantum system themalises because its individual energy eigenstates $\ket{E}$ appear thermal with respect to expectation values of realistic observables $\hat{\mathcal{O}}$ (that is, $\bra{E}\hat{\mathcal{O}}\ket{E} \approx \mathcal{O}_\text{th}(E)$ where $\mathcal{O}_\text{th}(E)$ is the thermal expectation value at the energy $E$).

It is well established that some systems can fail to thermalise. Integrable models, for example, have a large number of local integrals-of-motion that constrain the dynamics so that the state space cannot be explored ergodically. From the perspective of the ETH, this failure to thermalise in quantum integrable models is attributed to rare non-thermal, i.e., $\bra{E}\hat{\mathcal{O}}\ket{E} \not\approx \mathcal{O}_\text{th}(E)$, energy eigenstates $\ket{E}$ of its Hamiltonian \cite{Bir-10}. Similar behaviour is also expected in many-body localised (MBL) models, where disorder induced local integrals-of-motion lead to an effective integrability \cite{Gor-05,Bas-06,Pal-10}.

For non-integrable systems, anomalous thermalisation can arise due, for example, to pre-thermalisation \cite{Aba-17, Els-17, Nul-20}, strong zero-modes \cite{Kel-15, Fen-16, Mor-17}, or metastability \cite{Lan-18, Pan-20}. Recently, however, it was also discovered that some non-integrable systems can fail to thermalise due to rare non-thermal eigenstates called quantum many-body scars (QMBS) \cite{Shi-17, Mou-18a, Mou-18b, Khe-19, Ho-19, Shi-20, Tur-18a, Tur-18b, Shi-19, Lin-19, Cho-19, Khe-20, Doo-20b, Doo-21a}. One of the first and most prominent examples of QMBS is in the so-called PXP Hamiltonian \cite{Tur-18a, Tur-18b, Shi-19, Lin-19, Cho-19}, which was used to model an experiment on a chain of interacting Rydberg atoms \cite{Ber-17, Les-11}. A defining feature of the PXP Hamiltonian is a constraint that makes certain states inaccessible to the dynamics. Similar dynamical constraints also play a central role in a systematic route to constructing QMBS via the ``embedding method'' of Ref. \cite{Shi-17}, and in the anomalously slow thermalisation observed in Refs. \cite{Lan-18,Pan-20}.

In this paper, our aim is to further explore the role of such local dynamical constraints in ETH-violation. To this end, in section \ref{sec:model} we construct a family of spin chain models that generalises the PXP model. The key feature here is a discrete variable $\ell$ that controls the strength of a dynamical constraint. For the strongest non-trivial constraint ($\ell = 2$) we recover the PXP Hamiltonian. Increasing $\ell$ weakens the constraint, in the sense that fewer states are excluded from the dynamics.

In section \ref{sec:results}, we show that the scarring (i.e., the number of non-thermal states) becomes more extreme when the constraint is weakened. This is directly observable in the behaviour of the eigenstate expectation values $\langle E | \hat{\mathcal{O}} | E \rangle$ and in their fluctuations around their micro-canonical average. We argue that this behaviour occurs because, as $\ell$ increases, the number of energy eigenstates $\ket{E}$ that can evade the constraint grows dramatically. On the other hand, for finite $\ell$, it remains the case that most eigenstates $\ket{E}$ do not evade the constraint. Thus the generic behaviour of the bulk spectrum remains distinctly thermal.

\section{Model}\label{sec:model}

We consider a one-dimensional chain of $L$ spin-$1/2$ particles with the Hamiltonian: \begin{equation} \hat{H}_{\ell} = \hat{P}_{\ell} \hat{H}_0 \hat{P}_{\ell} , \qquad \hat{H}_0 = \frac{\Omega}{2}\sum_{i=1}^L \hat{\sigma}_i^x , \label{eq:H_l} \end{equation} where the projector: \begin{equation} \hat{P}_{\ell} = \prod_{i} \left( \hat{\mathbb{I}} - \ket{\uparrow_i\uparrow_{i+1}\hdots\uparrow_{i+\ell-1}}\bra{\uparrow_i\uparrow_{i+1}\hdots\uparrow_{i+\ell-1}} \right) , \label{eq:P} \end{equation} enforces a dynamical constraint of radius $\ell \in \{ 1,2,...,L \}$. Here, $\{ \ket{\uparrow_i}, \ket{\downarrow_i} \}$ is a basis for the spin-$1/2$ subsystem at site $i$ of the chain and $\hat{\sigma}_i^x = \ket{\uparrow_i}\bra{\downarrow_i} + \ket{\downarrow_i}\bra{\uparrow_i}$.

The constraint can be seen as excluding any states with $\ell$ consecutive $\uparrow$-states from the dynamics, giving a constrained state space that is smaller than the full $2^L$-dimensional unconstrained state space. For example, if $\ell = 3$, states such as $\ket{\hdots \uparrow\uparrow\uparrow \hdots}$ are annihilated by the projector $\hat{P}_{\ell = 3}$, and are therefore trivial zero-energy eigenstates of the Hamiltonian $\hat{H}_{\ell = 3}$ and may be neglected. The PXP Hamiltonian is a special case for $\ell = 2$.

In Appendix \ref{app:underlying_H} we show that our Hamiltonian Eq. \ref{eq:H_l} may be derived as an effective model, starting from an underlying Hamiltonian with strong $\ell$-body interactions between the spin-1/2 particles on the chain. Since $\ell$-body interactions may be considered unphysical for large $\ell$ we also show in Appendix \ref{app:underlying_H} that our Hamiltonian can be mapped to a more physically plausible underlying model with $2$-body interactions between spin-XX particles. We also note that the constraints in our Hamiltonian Eq. \ref{eq:H_l} can, in principle, be generated by repeated fast projective measurements, as was done experimentally in Ref. \cite{Bar-15} for $\ell = L$.

Our constrained Hamiltonian $\hat{H}_{\ell}$ has several symmetries. A reflection $i \to L - i + 1$ of site indices about the midpoint of the chain leaves the Hamiltonian invariant, implying that it has a spatial parity symmetry. Also, if there are periodic boundary conditions, the Hamiltonian has a translation symmetry under the transformation $i \to i+1$ of the site index. In the following, we can restrict to a symmetry sector of $\hat{H}_\ell$ to reduce the numerical cost of calculating its eigensystem.

\section{Results}\label{sec:results}

\subsection{Increasing $\ell$ weakens the dynamical constraint}

Let $\mathcal{D}_{L,\ell}$ denote the dimension of the constrained state space of the Hamiltonian $\hat{H}_\ell$ in a chain of length $L$. The strongest possible constraint is for $\ell = 1$, since the constrained Hilbert space dimension is $\mathcal{D}_{L,1} = 1$ and only the state $|\downarrow\rangle^{\otimes L}$ is allowed. The constraint becomes weaker as $\ell$ increases, in the sense that fewer and fewer states are excluded from the constrained state space. The weakest possible constraint is for $\ell = L$, since only a single state $|\uparrow\rangle^{\otimes L}$ is excluded and the constrained state space dimension is $\mathcal{D}_{L,L} = 2^L - 1$.

More generally, for a spin chain with open boundary conditions (OBC) the dimension $\mathcal{D}_{L,\ell}^\text{OBC}$ of the constrained subspace is given by the recurrence relation (see Appendix \ref{app:dimension} for the proof): \begin{equation} \mathcal{D}_{L,\ell}^\text{OBC} = \mathcal{D}_{L-1,\ell}^\text{OBC} + \mathcal{D}_{L-2,\ell}^\text{OBC} + \hdots + \mathcal{D}_{L-\ell,\ell}^\text{OBC}  \label{eq:D_OBC} \end{equation} with the initial conditions $\mathcal{D}_{L,L}^\text{OBC} = 2^L - 1$ and $\mathcal{D}_{L,\ell}^\text{OBC} = 2^L$ if $L<\ell$. This recurrence relation defines a generalised Fibonacci sequence, in which the next number in the sequence is obtained as the sum of the previous $\ell$ numbers. For instance, $\mathcal{D}_{L,2}^\text{OBC}$ is the ($L+2$)'th number in the usual Fibonacci sequence \cite{Les-12a, Tur-18a}, while $\mathcal{D}_{L,3}^\text{OBC}$ is the ($L+3$)'th number in the so-called tribonacci sequence.

\begin{figure}[t]

  \includegraphics[width=\columnwidth]{constrained_dim.png}
  \caption{(a) for $L \gg 1$ the constrained state space dimension increases exponentially in the system size $\mathcal{D}_{L,\ell} \sim (d_\ell)^L$ (solid lines for OBC and dotted lines for PBC are almost indistinguishable in the figure). For OBC one can show that $d_\ell$ is the solution to the equation $\ell = \log\left( 2 - d_\ell \right)^{-1} / \log d_\ell$, which is plotted in (b) (see Ref. \cite[p.~101]{Gar-61}). For $\ell = 2$ this gives $d = (1+\sqrt{5})/2$, the golden ratio, while for $\ell \to \infty$ we have $d \to 2$.}
\label{fig:constrained_dim}
\end{figure}

For periodic boundary conditions (PBC), the constrained Hilbert space dimension has the slightly more complicated form (again, see Appendix \ref{app:dimension} for the proof): \begin{eqnarray} \mathcal{D}_{L,\ell}^\text{PBC} &=& \mathcal{D}_{L-1,\ell}^\text{OBC} + \mathcal{D}_{L-3,\ell}^\text{OBC} + 2 \mathcal{D}_{L-4,\ell}^\text{OBC} \nonumber \\ && + 3 \mathcal{D}_{L-5,\ell}^\text{OBC} + \hdots + (\ell - 1)\mathcal{D}_{L-\ell-1,\ell}^\text{OBC} . \label{eq:D_PBC} \end{eqnarray}

In Fig. \ref{fig:constrained_dim}(a) we plot $\mathcal{D}_{L,\ell}$ as a function of the system size $L$, for both OBC and PBC. We see that for $L \gg 1$ the constrained state space dimension grows exponentially in the system size $\mathcal{D}_{L,\ell} \sim d_\ell^L$, for some number $d_\ell$. Fig. \ref{fig:constrained_dim}(b) shows $d_\ell$ as a function of the constraint radius $\ell$. As $\ell$ increases, $d_\ell$ also increases, with $d_\ell \to 2$ in the $\ell \to \infty$ limit. We interpret this as meaning that increasing $\ell$ weakens the dynamical constraint, and that the constraint becomes negligible as $\ell\to\infty$.

\subsection{The finite-$\ell$ constraint breaks integrability}

\begin{figure}[h]

  \includegraphics[width=0.9\columnwidth]{spectrum_b.png}
  \caption{Left column: the distribution of energy level spacings $p(s_\alpha/\bar{s})$ in the zero-momentum and even reflection-parity symmetry sector of $\hat{H}_\ell$. Level statistics are calculated in an energy window $\Lambda_E = [E- \Delta E, E+\Delta E]$, shown in the gray shaded region in the corresponding figure in the right column. The energy window $\Lambda_E$ is chosen with $\Delta E = 0.45\Omega$ and is centred at $E = 0.5 \Omega$ if $L$ is even or $E = \Omega$ if $L$ is odd. This ensures that $\Lambda_E$ excludes the degenerate peaks when $\ell \geq 4$, but captures the states between two peaks. Right column: The density of states of $\hat{H}_\ell$ in its zero-momentum and even reflection-parity symmetry sector. The sharp peaks correspond to degeneracies in the spectrum.}
\label{fig:spectrum}
\end{figure}

The unconstrained Hamiltonian $\hat{H}_0 = \frac{\Omega}{2}\sum_i \hat{\sigma}_i^x$ is non-interacting and is clearly integrable. However, the constraint $\hat{P}_\ell$ breaks the integrability for finite $\ell$. This is shown numerically in the left column of Fig. \ref{fig:spectrum}, where we plot the normalised level spacing distribution $p(s_\alpha/\bar{s})$, where $s_\alpha = E_{\alpha + 1} - E_\alpha$ is the spacing between two consecutive energy eigenvalues in a narrow energy window, and $\bar{s}$ is the mean level spacing in the energy window. We note that in calculating $p(s_\alpha/\bar{s})$ we restrict to the even-reflection and zero-momentum (we have assumed PBC) symmetry sector of $\hat{H}_\ell$. For $\ell \in \{2,3,\hdots,6 \}$ the distribution shows clear level-repulsion, and is close to the distribution of the Gaussian orthogonal ensemble (GOE) of random matrices, as expected for a non-integrable Hamiltonian.

Another standard numerical test of integrability is to compute the $r$-value \cite{Oga-07}, defined as $\langle r_\alpha \rangle = \langle \text{min}(s_\alpha,s_{\alpha+1})/\text{max}(s_\alpha,s_{\alpha+1}) \rangle$ where the average $\langle\bullet\rangle$ is taken over all eigenvalues in the symmetry sector.  In Fig. \ref{fig:r_statistics} we see that the (slightly perturbed   \footnote{We use a slightly perturbed Hamiltonian  $\hat{H}_0 = \frac{\Omega}{2}\sum_{i} \hat{\sigma}_i^x + \lambda \sum_{i} (\hat{\sigma}_i^x \hat{\sigma}_{i+1}^z + \hat{\sigma}_i^z \hat{\sigma}_{i+1}^x), \quad |\lambda| \ll |\Omega| $  to calculate the $r$-values in Fig. \ref{fig:r_statistics}. This is because the $r$-value is not well defined when there are a large number of degeneracies in the spectrum. the perturbation above has the effect of breaking these degeneracies by a small amount so that the $r$-value gives meaningful results.})  $\hat{H}_0$ is integrable, but that the constraint $\hat{H}_\ell = \hat{P}_\ell \hat{H}_0 \hat{P}_\ell$ breaks the integrability for $\ell \in \{ 2,3,\hdots,7 \}$. In fact, Fig. \ref{fig:r_statistics} shows that the integrability is broken more strongly for $\ell \geq 3$ than for the PXP model ($\ell = 2$), in the sense that the $r$-value approaches the GOE value more quickly as the symmetry sector dimension increases.

\begin{figure}

  \includegraphics[width=\columnwidth]{r_values.pdf}
  \caption{The $r$-value is a standard numerical test of integrability. Here, we plot the plot the $r$-value as a function of the dimension $\tilde{\mathcal{D}}_{L,\ell}$ of the zero-momentum, even-reflection symmetry sector of $\hat{H}_\ell = \hat{P}_\ell \hat{H}_0 \hat{P}_\ell$, where $\hat{H}_0 = \frac{\Omega}{2}\sum_{i=1}^L \hat{\sigma}_i^x + \lambda \sum_{i} (\hat{\sigma}_i^x \hat{\sigma}_{i+1}^z + \hat{\sigma}_i^z \hat{\sigma}_{i+1}^x)$ is a (slightly perturbed, $\lambda = 0.01 \times \Omega$) unconstrained Hamiltonian. For $\ell = 0$ we have $\langle r_\alpha \rangle \to 0.386$ as the system size increases, as expected for an integrable model. For $\ell \geq 2$, however, it tends towards the value $\langle r_\alpha \rangle \approx 0.53$ that is expected for nonintegrable systems (GOE).}
  \label{fig:r_statistics}
\end{figure}

\subsection{Degeneracies in the $\ell \geq 4$ constrained Hamiltonian}

The eigenstates of the unconstrained Hamiltonian $\hat{H}_0 = \frac{\Omega}{2}\sum_i \hat{\sigma}_i^x$ are the product states $\{ \ket{x_1,x_2,...,x_L} \}$, where each $|x_i\rangle \in \{ \ket{\pm} \}$ is a $\hat{\sigma}^x$-eigenstate. The eigenvalues $E_\alpha = \frac{\Omega}{2} \sum_{i=1}^L x_i $ are highly degenerate, and are integer multiples of $\Omega$ (for even $L$), or half-integer multiples of $\Omega$ (for odd $L$). This is illustrated in Fig. \ref{fig:spectrum}(a), where a histogram shows the density of states in the zero-momentum, even-reflection symmetry sector. The degeneracies are indicated in the histogram by sharp peaks at integer/half-integer multiples of $\Omega$. Of course, any linear combination of degenerate eigenstates is also a valid eigenstate. For example, the Dicke states $\ket{j,m_x}$, defined as simultaneous eigenstates of $\hat{J}^x \equiv \frac{1}{2} \sum_{j=1}^L \hat{\sigma}^x_j $ [with eigenvalue $m_x \in \{-j,-j+1,\hdots,j\}$] and $\hat{J}^2 \equiv (\hat{J}^x)^2 + (\hat{J}^y)^2 + (\hat{J}^z)^2$ [with eigenvalue $j(j+1)$], are also valid eigenstates of the unconstrained Hamiltonian.

\begin{figure}[b]

  \includegraphics[width=\columnwidth]{degeneracies.png}
  \caption{For $\ell \geq 4$ and $L$ even (odd) the degeneracy at integer (half-integer) multiples of $\Omega$ increases exponentially in $L$. This resembles the degeneracy of the unconstrained $\ell = 0$ model (the purple trianglular markers). [Plotted data is restricted to the even-reflection and zero-momentum symmetry sector (PBC are assumed).]}
\label{fig:degeneracies}
\end{figure}

The right column of Fig. \ref{fig:spectrum} also shows the density of states for the constrained Hamiltonian $\hat{H}_\ell$ with $\ell \in \{ 2,3,...,6 \}$ in the zero-momentum (assuming PBC), even-reflection symmetry sector. For $\ell \in \{ 2, 3 \}$ [Fig. \ref{fig:spectrum}(i,k)] the density of states is approximately a smooth Gaussian \cite{Tur-18a}. As the constraint is weakened to $\ell \geq 4$, however, we begin to see deviation from the smooth density of states. This is due to the appearance of sharp peaks in the middle of the spectrum at integer or half-integer multiples of $\Omega$ [see Fig. \ref{fig:spectrum}(g,e,c)]. We find numerically that each of these peaks is degenerate, with the number of degenerate states increasing exponentially in the chain length $L$ (see Fig. \ref{fig:degeneracies}). This is reminiscent of the degeneracies in the $\ell=0$ unconstrained model, which are also exponentially growing in system size (see the purple markers in Fig. \ref{fig:degeneracies}). For the unconstrained model the degeneracies are associated with the conserved magnetization $\hat{J}^x$ and the associated quantum number $m_x$. However, for any non-zero $\ell$ this symmetry is broken and the degeneracies at integer/half-integer multiples of $\Omega$ are not associated with any obvious local conserved quantities. In the next section we will show that the highly degenerate states in the spectrum of $\hat{H}_{\ell \geq 4}$ are in fact many-body scar states.

\subsection{Extreme quantum scarring}

Let $\ket{E_\alpha}$ and $E_\alpha$ be the eigenstates and eigenvalues of $\hat{H}_\ell$. For an observable $\hat{\mathcal{O}}$, and a microcanonical energy window $\Lambda_E \equiv [E - \Delta E, E + \Delta E]$, we define the microcanonical average as $\overline{\mathcal{O}}_{\Lambda_E} \equiv \mathcal{N}_{\Lambda_E}^{-1} \sum_{E_\alpha \in \Lambda_E} \mathcal{O}_{\alpha\alpha}$, where $\mathcal{O}_{\alpha\alpha} = \bra{E_\alpha}\hat{\mathcal{O}}\ket{E_\alpha}$ are the eigenstate expectation values (EEVs) of the observable and $\mathcal{N}_{\Lambda_E}$ is the number of states in $\Lambda_E$. With this definition one can formulate a strong and a weak version of the ETH.

For strong-ETH the condition $\mathcal{O}_{\alpha\alpha} \approx \overline{\mathcal{O}}_{\Lambda_E}$ must be satisfied by \emph{all} eigenstate expectation values in the energy window. More precisely, the strong-ETH is satisfied if the quantity
\begin{equation}
 I_s = \max_{E_\alpha \in \Lambda_E} | \mathcal{O}_{\alpha\alpha} -  \overline{\mathcal{O}}_{\Lambda_E}|
\end{equation}
 vanishes in the thermodynamic limit.

 For weak-ETH it is enough that \emph{most} eigenstates in the energy window satisfy $\mathcal{O}_{\alpha\alpha} \approx \overline{\mathcal{O}}_{\Lambda_E}$ \cite{Mor-18}. In other words, the mean fluctuation around the microcanonical average
 \begin{equation}
 I_w = [ \mathcal{N}_{\Lambda_E}^{-1} \sum_{E_\alpha \in \Lambda_E} ( \mathcal{O}_{\alpha\alpha} -  \overline{\mathcal{O}}_{\Lambda_E} )^2 ]^{1/2}
 \end{equation}
 should vanish in the thermodynamic limit.

There is now considerable numerical evidence suggesting that the weak-ETH holds for generic integrable and non-integrable quantum systems \cite{Mor-18}: For non-integrable systems $I_w$ decays exponentially with system size $L$ \cite{Ste-14b,Beu-14,Kim-14,Ike-15}, while for generic integrable systems, $I_w$ has a power-law decay with system size $L$ \cite{Bir-10,Ike-13,Alb-15}. Since quantum many-body scars are strong-ETH violating states, but in non-integrable systems, they can be identified through exponential decay of $I_w$ and a non-decaying $I_s$.

\begin{figure

  \includegraphics[width=\columnwidth]{EEVs.png}
  \caption{The left column shows eigenstate expectation values (EEVs) of the observable $\hat{\mathcal{O}} = \hat{J}^z/L$, for constraint radius $\ell = 0,6,5,4,3,2$ (from top to bottom). The color scale indicates the density of datapoints. The right column shows the ETH indicators $I_w$ (red) and $I_s$ (blue) plotted against the dimension $\mathcal{D}_{L,\ell}^{k=0}$ of the zero-momentum ($k=0$) symmetry sector. The indicators are calculated in an energy window $\Lambda_E = [E-\Delta E, E + \Delta E]$, with $\Delta E = 0.5\Omega$ and centred at $E = 1.0\Omega$ ($E=1.5\Omega$ ) if $L$ is even (odd) (e.g., the shaded gray regions in the left column). We see that the weak-ETH is satisfied for the constrained models ($I_w$ decays as a power law in $\mathcal{D}$, i.e., exponentially in system size $L$), but the strong-ETH is violated ($I_s$ does not decay with system size). [Note that we have used $I_s = \max_{\Lambda_E} (\mathcal{O}_{\alpha\alpha} - \mathcal{O}_{\Lambda_E})$ instead of $I_s = \max_{\Lambda_E} |\mathcal{O}_{\alpha\alpha} - \mathcal{O}_{\Lambda_E}|$, since this gives cleaner results and is still a valid indicator of strong-ETH violation.]}
\label{fig:EEVs}
\end{figure}

In the left column of Fig. \ref{fig:EEVs} we plot the eigenstate expectation values (EEVs) for the observable $\hat{J}^z = \frac{1}{2}\sum_i \hat{\sigma}_i^z$ in the zero-momentum symmetry sector of $\hat{H}_\ell$ (again, we have assumed PBC). In a narrow energy window $\Lambda_E$ away from the edges of the spectrum (the grey shaded areas in the left column of Fig. \ref{fig:EEVs}) we see that the EEVs are, for the most part, concentrated around the microcanonical average (marked by a black $\mathsf{x}$). However, there are some eigenstates in $\Lambda_E$ for which the EEV deviates significantly from the microcanonical average, giving a relatively large value of the strong-ETH indicator $I_s$. In the right column of Fig. \ref{fig:EEVs} the blue lines show that $I_s$ does not decay as the system size increases, verifying that the strong-ETH is indeed violated for these eigenstates. The special eigenstates that prevent the decay of $I_s$ are regarded as quantum many-body scars.

From Fig. \ref{fig:EEVs}(h), we notice that for $\ell = 3$ there are two scar eigenstates with a particularly pronounced violation of the strong-ETH. We find that, for any odd $L \geq 7$, these two scar states appear in the zero-momentum sector with \emph{exactly} $E_\alpha = \pm\sqrt{7}/2$ and $J^z_{\alpha\alpha} = -3/7$ (up to numerical precision) \footnote{We note that, for $L$ odd, as we increase $L \to L+2$ this special state alternates between the $(k,p) = (0,1)$ and the $(k,p)=(0,-1)$ sectors, where $k$ denotes momentum and $p$ denotes reflection-parity.}.  We have also found numerically that these two special states are eigenstates of the total angular momentum operator $\hat{J}^2 = (\hat{J}^x)^2 + (\hat{J}^y)^2 + (\hat{J}^z)^2$, with the total angular momentum $j = 3/2$.

For $\ell \geq 4$ the scar states that are clearly visible in Figs. \ref{fig:EEVs}(b,d,f) at $J^z_{\alpha\alpha} = 0$ are exactly the degenerate eigenstates that were discussed in the previous section, with energies that are integer/half-integer multiples of $\Omega$ for $L$ even/odd. A closer examination of these degenerate integer/half-integer energy eigenstates shows that they are Dicke states $\ket{j,m_x}$, with $m_x = E_\alpha / \Omega$ and $j \ll N/2$. In other words, they are eigenstates of the unconstrained Hamiltonian $\hat{H}_{\ell = 0}$ as well as the constrained Hamiltonian $\hat{H}_{\ell \geq 4}$. This suggests that at $\ell = 4$ the constraint has been weakened sufficiently that some of the eigenstates of $\hat{H}_0$ can avoid the constraint completely. We note that $\bra{j,m_x} \hat{J}^z \ket{j,m_x} = 0$ for Dicke states $\ket{j,m_x}$, which explains why $J^z_{\alpha\alpha} = 0$ for these scar states in Figs. \ref{fig:EEVs}(b,d,f).

A clearly visible trend in the left column of Fig. \ref{fig:EEVs} is that the scarring becomes more severe as $\ell \geq 3$ increases. It is natural to ask whether the scarring can become so extreme that -- not only is the strong-ETH violated -- but the weak-ETH is also violated. This might occur if, for example, there are always a similar number of Dicke scar states and thermal states in the energy window $\Lambda_E$, preventing the decay of the weak-ETH indicator $I_w$ as the system size increases. In the right column of Fig. \ref{fig:EEVs}, however, the red lines show that the weak indicator decays as $I_w \sim \mathcal{D}^{-\gamma}$ for $\ell \in \{2,3,4,5,6 \}$, where $\mathcal{D}$ is the zero-momentum symmetry sector dimension. This suggests that, even though the number of scars is increasing exponentially, in the large $L$ limit they still only make up a negligible fraction of the total number of eigenstates in the energy window.

\begin{figure}[t]

  \includegraphics[width=\columnwidth]{I_w_decay_exponent.png}
  \caption{As shown in the right column of Fig. \ref{fig:EEVs} above, the weak ETH-indicator decays as $I_w \sim \mathcal{D}^{-\gamma}$. Here we plot the inverse of the decay exponent. We see here that for $4 \leq \ell \leq 8$ it closely fits the form $\gamma = 1/(a\ell + b)$, for constants $a$ and $b$. Error bars are calculated from the error in fitting to the red markers in the right column of Fig. \ref{fig:EEVs}.}
\label{fig:I_w_decay_exponent}
\end{figure}

However, it is also clear that for $\ell \geq 3$ the decay exponent $\gamma$ becomes smaller as $\ell$ increases. In Fig. \ref{fig:I_w_decay_exponent} we plot $\gamma$ as a function of $\ell$ for $\ell \geq 4$. It appears to be well approximated by the line $\gamma = 1/(a\ell + b)$ for real constants $a$ and $b$. This trend suggests that the $I_w \sim \mathcal{D}^{-\gamma}$ scaling persists for finite $\ell$, although the decay exponent can be very small for large $\ell$. This slow decay of $I_w$ with system size reflects the high degree of scarring for large finite values of $\ell$.

\section{Conclusion}\label{sec:conclusion}

In this paper we investigate the phenomenon of many-body scarring in quantum systems with dynamical constraints. We construct a class of spin-chain models for which the radius of a local dynamical constraint is given by a discrete variable $\ell$. Increasing $\ell$ weakens the constraint by decreasing the dimension of the subspace that is excluded from the dynamics. We have shown that increasing $\ell$ also corresponds to more extreme quantum many-body scarring, with the number of scar states increasing exponentially in system size for $\ell \geq 4$. Another example of a model exhibiting an exponentially increasing number of scars was given in Ref. \cite{Shi-17}. However, our model differs in our ability to adjust the severity of scarring with our parameter $\ell$.

To the best of our knowledge, our model $\hat{H}_\ell$ is the first that allows the degree of quantum scarring to be tuned. Since the scarring becomes more extreme with increasing $\ell$ it would be interesting to test our results for larger values, and to explore the ultimate limit of extreme scarring. With our current numerical approach (exact diagonalisation of dense matrices) this is challenging due to finite-size effects becoming more significant for larger values of $\ell$. In the future larger system sizes may be accessed, however, by using the ``shift-invert'' algorithm with sparse matrices, and focussing on a narrow energy window at finite density \cite{Beu-14, Lui-15}.

\emph{Note added:} During the completion of this work we became aware that the same model is also studied in Ref. \cite{Des-21}. In that work, the focus is on the adjacency graph corresponding to $\hat{H}_\ell$ (and other constrained models), which is then used to provide interesting insights into the origin of many-body wave function revivals in the dynamics.

\begin{acknowledgments}
We would like to thank Masudul Haque, Ian Jubb, Kevin Kavanagh, and Luuk Coopmans for helpful discussions. This work was funded by Science Foundation Ireland through Career Development Award No. 15/CDA/3240. SD also acknowledges financial support from the SFI-EPRC joint project QuamNESS. We thank the DJEI/DES/SFI/HEA Irish Centre for High-End Computing (ICHEC) for the provision of computational facilities.
\end{acknowledgments}


\section{\label{sec:level1}Introduction}

Describing the propagation of a foreign particle through a medium---the polaron problem---is a fundamental
physics problem, with a wide variety of applications \cite{Devreese2009, Kuper1963}.
The problem is especially non-trivial for quantum systems at low temperatures and
in reduced dimensions, which enhance quantum fluctuations. Recent advances in
experimental techniques in the field of ultracold gases allow efficient
realization, control and measurements of strongly interacting quantum systems,
including two-component gases with large population imbalance, of which a single
a polaron problem is a limiting case \cite{Palzer2009, Catani2012, Meinert2017, Fukuhara2013}.

The motion of the impurity injected into a one-dimensional (1D) degenerate zero-temperature
quantum gas strongly depends on the ratio of the impurity mass, and the mass of the host particles.
This is already seen from the semiclassical analysis of the kinetic Boltzmann equation,
which predicts that the terminal velocity of the impurity depends on the mass
ratio \cite{GamayunLychkovskiy2014}. For equal masses, the Boltzmann equation
approach fails unless multiple coherent scattering processes are taken into account
\cite{Gamayun2014}.

For equal masses, however, an alternative numerical approach is based on a Bethe Ansatz
solvable model where the host gas consists of non-interacting fermions
and impurity-host interaction is a short-range delta-function repulsion (see Eq.\ \eqref{h0} below).
This model, a special case of a Gaudin-Yang model \cite{Gaudin1967, Yang1967}, is integrable
and the exact spectrum can be computed by solving a system of Bethe Ansatz equations \cite{McGuire1965}.
Explicit expressions for matrix elements for several observables are known \cite{Zotos1993, Mathy2012},
and some physically relevant observables can be computed numerically via summations
over subsets of Bethe Ansatz states \cite{Mathy2012, Burovski2014, Robinson2016, Gamayun2018}.

In this paper, we consider the case where the impurity mass differs from that of the host particles.
While the exact Bethe Ansatz solution is not available, we construct an exact diagonalization
procedure using a suitably truncated basis of the Bethe Ansatz states of an integrable model with equal masses
(so that the off-diagonal matrix elements are related to the mass difference).
Applying our procedure to a mesoscopic system of and impurity and $N=5$ host particles, we study
the evolution of the spectrum with the mass ratio and observe clear signatures of level repulsion
and a redistribution of the spectral weight for $m/M \neq 1$.
Two notes are in order. First, we note that our approach is non-perturbative and does
not rely on either a coupling constant or a mass difference being small. Second,
our numerical procedure is controllable: the truncation error is directly related to the
saturation of the sum rule related to the completeness of the full basis of the Hilbert space.

Our approach is inspired by the so-called Truncated Spectrum Approach (TSA), which is
used to describe perturbations of conformal field theories by relevant operator
(see e.g. Ref.\ \cite{James2018} for a review). In this approach, the subset of eigenstates
that forms a basis for the perturbed Hamiltonian is truncated by their energy values,
which makes the method tailored to capture low-energy physics. The dependence of the
energy cut-off is regulated by the corresponding renormalization group.
In contrast to the TSA approach, we rely on the saturation of a sum rule; this
is similar to the application of
the TSA for the quench description in the Lieb-Liniger model \cite{Caux2012}.
(See {Fig.\ \ref{fig:rho_convergence}} below for a comparison of the ordering strategies.)

It is also instructive to compare our technique to alternative numerical approaches, which can
be extended to the unequal mass problem. One popular approach is to directly simulate
real-time dynamics using time-dependent density-matrix renormalization
group \cite{Knap2014, Massel2013, Peotta2013,Ozaki2015}. Our approach is complementary, as
it is based on a different strategy of approximating the low-lying part of the spectrum.
Directly solving a many-body Schroedinger equation was used in Ref. \cite{Petrov2011} for
studying the response to a parametric modulation of the coupling constant for
both equal and unequal mass systems. This approach is obviously superior for studying
few-particle systems, but quickly becomes impractical for larger number of particles.

The rest of the paper is organized as follows. In Sec.\ \ref{sec:model} we detail the microscopic
model, and discuss our numerical approaches for studying both integrable, equal-mass
(Sec.\ \ref{subsec:eig_int}) and non-integrable, unequal-mass (Sec.\ \ref{subsec:eig_non_int})
cases. Sec.\ \ref{sec:results} presents results of our numerical experiments, and
Sec.\ \ref{sec:concl} concludes the paper.

\section{Model and method}
\label{sec:model}

We consider an impurity of mass $M$ injected into the 1D gas of free spinless fermions
(equivalently, Tonks-Girardeau bosons) at $T=0$.
The host gas contains $N$ identical fermions of mass $m$. The impurity interacts
with the host particles via the short-range delta-function repulsion. Then the
Hamiltonian in the first quantization reads

\begin{equation}
\widehat{H} = \frac{\widehat{P}^2}{2M} + \sum_{i=1}^N \frac{\widehat{p_i}^2}{2m} + \frac{g}{m} \sum_{j=1}^N \delta(X - x_j) \,.
\label{ham_orig}
\end{equation}

Here $X$ is the coordinate of the impurity and $x_j$ is the coordinate of the $j$-th host
particle, $\widehat{P}$ and $\widehat{p_j}$ are the momenta operators of the impurity and  $j$-th
host particle, respectively; $g > 0$ is the strength of the delta-function coupling between
host particles and the impurity. In Eq.\ \eqref{ham_orig} and below we use the
units where the Planck's constant $\hbar=1$, and we also set $2m=1$.

For equal masses, $M=m$, the model \eqref{ham_orig} admits the exact solution
via Bethe Ansatz \cite{McGuire1965}. For unequal masses, $m\neq M$, the model is no longer integrable, and the Bethe
Ansatz machinery of \cite{McGuire1965} breaks down.

We identically rewrite \eqref{ham_orig} to separate the integrable part, $\widehat{H_0}$,
which is nothing but Eq.\ \eqref{ham_orig} with $m=M$:

\begin{gather}
\widehat{H} = \widehat{H_0} + \left( \frac{1}{2M} - \frac{1}{2m} \right) \widehat{P^2} \,,  \label{h1}\\[1em]
\widehat{H_0} = \frac{\widehat{P}^2}{2m} + \sum_{i=1}^N \frac{\widehat{p_i}^2}{2m} + \frac{g}{m} \sum_{j=1}^N \delta(X - x_j) \label{h0}\,,
\end{gather}

Note that the second term in Eq.\ \eqref{h1} is off-diagonal in the basis of eigenstates of Eq.\ \eqref{h0}.

\subsection{Eigenstates of the integrable model}
\label{subsec:eig_int}

In this subsection, we discuss the eigenstates of the integrable model
\eqref{h0}. Following Refs. \cite{McGuire1965, Mathy2012}, we consider $N$ host
particles with periodic boundary conditions on a ring of circumference $L$.
Consider the spectrum of the Hamiltonian \eqref{h0} with a given total momentum,
$Q$. Let $| \psi \rangle$ is an eigenstate of $\widehat{H_0}$ with the energy eigenvalue
$\varepsilon_\psi$: $\widehat{H_0}|\psi\rangle = \varepsilon_\psi |\psi\rangle$.
The eigenstates are fully characterized by a set of $N+1$ rapidities,
$z_j$, $j=1, \dots, N+1$. Rapidities are roots of the Bethe equation \cite{McGuire1965},
\begin{equation}
1 / \tan{z_j} = a z_j - c \;,
\label{eq:ctg}
\end{equation}
where $a = 8 / g L$ and the free term $c$ ensures the conservation of the total
momentum of the system,
\begin{equation}
Q = \frac{2}{L}\sum_{j=1}^{N+1} z_j \;.
\label{eq:Q}
\end{equation}

Due to periodicity of the cotangent function, solutions of the Bethe equation
\eqref{eq:ctg} can be written as
\begin{equation}
z_j = \pi n_j - \delta_j\;,
\label{eq:phases}
\end{equation}
where $n_j$ are integers and phase shifts $\delta_j \in [0, \pi)$. Therefore,
an eigenstate of Eq.\ \eqref{h0} with a given value of the total momentum $Q$
is completely defined by specifying an integer partition $\vec{\lambda}$, i.e. a set of
$N+1$ distinct integers, $n_j$,
\begin{equation}
\vec{\lambda} = \{n_j \in \mathbb{Z}\;,\quad j=1,\cdots, N+1\} \;.
\label{eq:partition}
\end{equation}
Without loss of generality, we assume the partitions to be ordered, $n_1 < n_2 < \cdots n_{N+1}$.

The constructive procedure of finding an eigenstate of \eqref{h0} with a given
total momentum $Q$ is then as follows \cite{McGuire1965, Mathy2012}:
(i) fix a partition $\vec{\lambda}$, Eq.\ \eqref{eq:partition}, and (ii) solve
a system of $N+2$ nonlinear equations: $N+1$ Bethe equations \eqref{eq:ctg}
define the phase shifts, $\delta_j$, as functions of $c$; the latter is fixed by
solving Eq.\ \eqref{eq:Q} given the total momentum $Q$.

Given a partition $\vec{\lambda}$, and its corresponding rapidities, $z_j$, the
energy eigenvalue is
\begin{equation}
\varepsilon_\psi = \frac{1}{2m}\frac{4}{L^2}\sum_{j=1}^{N+1} z_j^2 \;,
\label{eq:e0}
\end{equation}
and the wave functions in the coordinate representation can be constructed explicitly
as certain Slater-type determinants \cite{Mathy2012} of size $N+1$.

Given the initial state at $t=0$, $|\Psi_0\rangle$, being the product state
of an impurity moving with the momentum $Q$ and the Fermi sea of the host fermions,
the completeness relation for the eigenstates of the Hamiltonian \eqref{h0} reads
\begin{equation}
\sum_{\psi} |\langle \Psi_0 | \psi \rangle |^2 = 1 \,.
\label{eq:sum_rule}
\end{equation}
Note that the size of the Hilbert space is infinite even for a finite number of
host particles $N$. The challenge is to find a relevant subset of eigenstates,

\begin{equation}
\mathcal{S} = \{ \psi : \sum_{\psi} |\langle \Psi_0 | \psi \rangle |^2 = 1 -\epsilon \}\;,
\label{eq:subspace}
\end{equation}
which saturates the sum rule to a predefined accuracy $\epsilon$.
The relevance of the various subsets of eigenstates for correlation functions
in the Lieb-Liniger model was recently under active exploration \cite{Panfil2021, Panfil2021_1}.

\subsubsection{Numerical approaches}
\label{sec:numerical}

Several approaches for generating representative subsets $\mathcal{S}$,
Eq.\ \eqref{eq:subspace}, have been devised. Numerical approaches are most
transparently described using the (formal) analogy between Bethe Ansatz
rapidities and spinless (pseudo)fermions. Specifically, rapidities, $z_j$,
can be formally interpreted as quasi-momenta of $N+1$ spinless pseudo-fermions
(see e.g. \cite{Caux2009}). Then, an arbitrary partition
$\vec{\lambda}$ is uniquely specified by listing
pseudo particle-hole pairs relative to the pseudo Fermi sea, which is nothing but
a set of $N+1$ consecutive integers
$$
\vec{\lambda}_0 = \{-(N+1)/2, \cdots, (N-1)/2 \}\;.
$$
For integrable models, it is typical that contributions from states with small
number of pseudo particle-hole pairs dominate \cite{Caux2009} and thus
numerics can be expected to converge with respect to taking into account
states with an increasing number of pseudo particle-hole pairs. In what follows
we call these pseudo particle-hole pairs ``excitations''.

Ref.\ \cite{Mathy2012} performed a brute-force enumeration of up to $10^5$
states with up to three excitations. When states are ordered by the absolute value of
the overlap with the in-state, $\langle \Psi_0 | \psi \rangle$, the saturation of the sum rule,
Eq.\ \eqref{eq:subspace} with the number of states, $N_s$, is reported to be a power law,
$\epsilon \sim N_s^\alpha$ with some value of the exponent $\alpha$.

Ref.\ \cite{Burovski2014} developed an alternative, \emph{stochastic} enumeration
approach of constructing a subset of states, Eq.\ \eqref{eq:subspace}, where the states
are generated via a Markov process in the space of partitions $\vec{\lambda}$,
based on the Metropolis algorithm \cite{Metropolis} with the transition probabilities
proportional to $|\langle \Psi_0 | \psi \rangle|^2$. This way, the process automatically finds
states with largest contributions to the sum rule \eqref{eq:sum_rule}. In practice,
the set of updates where we only change two values in a partition $\vec{\lambda}$
are sufficient for a quick convergence of the sum rule, \eqref{eq:subspace}.

The resulting eigenstates of Eq.\ \eqref{h0} can be classified into several families.
In the thermodynamic limit, $N\to\infty$, the sum rule \eqref{eq:sum_rule} is dominated
\cite{Burovski2014}
by the single-parameter family of states $\mathcal{S}_k$, of the states
\begin{equation}
\vec{\lambda}_k = \{ -(N+1)/2, \cdots, (N-1)/2, k \}\;,
\label{eq:dominant_family}
\end{equation}
where the pseudo-hole is fixed at $(N+1)/2$ and the pseudo-particle is located
at $n_{N+1} = k$. Note that the momentum conservation, \eqref{eq:Q}, constraints
$\mathcal{S}_k$---in fact, any single-excitation family with a fixed pseudo-hole
position ---to contain $O(N)$ states.

At any finite $N$, the contribution of other states is non-negligible. This includes
other single-excitation states, the pseudo Fermi sea, $\vec{\lambda}_0$, and multiply
excited states (see the top panel of Fig.\ \ref{fig:overlaps}).

\subsection{Eigenstates of the non-integrable model}
\label{subsec:eig_non_int}

We now turn our attention to the full Hamiltonian, Eq.\ \eqref{ham_orig}.
In the rest of the paper, we use Latin symbols ($|f\rangle$ etc) to label its eigenstates and
Greek symbols ($|\psi\rangle$ etc) to label eigenstates of the integrable model \eqref{h0}.

Consider a complete basis of the integrable Hamiltonian, $\widehat{H_0}$, Eq.\ \eqref{h0}.
In this basis, the full Hamiltonian, \eqref{h1}-\eqref{h0} has off-diagonal elements due to
the kinetic energy of the impurity, $\propto\widehat{P}^2$. Formally the matrix elements of
$\widehat{H}$ are given by
\begin{equation}
\langle \phi | \widehat{H} | \psi\rangle = \varepsilon_\psi \delta_{\phi\psi} +%
 \sum_{\alpha \in \mathcal{S}} \langle \phi | \widehat{P} | \alpha\rangle \langle \alpha | \widehat{P} | \psi\rangle \left(\frac{1}{2M} - \frac{1}{2m} \right)  \;,
\label{h_matrix}
\end{equation}
where $|\psi\rangle$ is an eigenstate of $\widehat{H}_0$ with the energy eigenvalue $\varepsilon_\psi$, and we
used the completeness relation $\sum_\alpha |\alpha\rangle \langle \alpha | = 1$.

Eq.\ \eqref{h_matrix} is exact if all eigenstates are included. We, instead, consider
a subspace, Eq.\ \eqref{eq:subspace}, of the full Hilbert space, and diagonalize
Eq.\ \eqref{h_matrix} \emph{in that subspace}.

The resulting eigenvalues of the truncated Hamiltonian \eqref{h_matrix}, $\varepsilon_f$, approximate
the exact eigenvalues of $\eqref{ham_orig}$---and the accuracy of
the approximation is controlled by the value of $\epsilon$, i.e. by the total missing
weight of the discarded eigenstates of the integrable model. We also note that the
total discarded weight is invariant under the diagonalization procedure, which is nothing but
a unitary transformation between two orthonormal bases, $\{ | \psi \rangle \}$ and $\{ | f\rangle \}$, thus
$\sum_{\psi\in\mathcal{S}} |\langle \psi | \Psi_0 \rangle |^2 = \sum_{f} \langle f | \Psi_0 \rangle |^2$.

To summarize our numeric procedure, we (i) enumerate a set of states $\mathcal{S}$
 \eqref{eq:subspace} with a predefined level $\epsilon$ using the Bethe Ansatz machinery
of the previous section; (ii) construct the truncated matrix of the Hamiltonian \eqref{h_matrix},
using the explicit expressions for the matrix elements $\langle \phi | \widehat{P} | \psi\rangle$ with
$\phi, \psi \in \mathcal{S}$, which were derived in Ref. \cite{Mathy2012}; and (iii) diagonalize
the resulting matrix numerically.

\section{Results and discussion}
\label{sec:results}

We first study the convergence of the sum rule, \eqref{eq:sum_rule} , with the
size of the subspace \eqref{eq:subspace}, i.e. the number of states, $N_s$, included into the $\mathcal{S}$.
The dependence of $\epsilon$ on $N_s$ is not unique and depends on which states are
included, not only on their number. Fig. \ref{fig:rho_convergence} shows an illustrative example
where we take $N=5$, $Q = 1.2 k_F$ (here $k_F = \pi N/L$ is the Fermi momentum of $N$ fermions)
and the dimensionless coupling constant $\gamma = g/(N/L) = 4.2$. We generate $3\times 10^3$
eigenstates of Eq.\ \eqref{h0} with up to three excitations using the stochastic enumeration
algorithm of Sec.\ \ref{sec:numerical}, which, collectively, saturate the sum rule \eqref{eq:subspace}
with $\epsilon = 2\times 10^{-5}$.

Fig.\ \ref{fig:rho_convergence} shows partial sums
as a function of $N_s$ for two strategies: where the states are included according
to their energy in the ascending order (red squares), and where the states are included
according their overlap with the $t=0$ state, $|\langle \psi | \Psi_0\rangle|^2$, in the descending order (open circles).
Clearly, the latter strategy results in much faster convergence of Eq.\ \eqref{eq:subspace}.
In this case, the dependence of $\epsilon$ on $N_s$ has three qualitatively different regimes.
For $10^2 \lesssim N_s \lesssim 10^3$, the convergence can be approximated by a power law
(cf Ref.\ \cite{Mathy2012}), $\epsilon \propto 1/N_s^a$ with $a=1.9$. For smaller
values of $N_s$, the data can be also approximated by a power law, with a smaller exponent, $a\sim 1.2$.
For very large values of $N_s$ (i.e., $\epsilon \lesssim 5\cdot 10^{-5}$),
Fig.\ \ref{fig:rho_convergence} clearly shows convergence which is slower than a power law.
We checked that qualitatively this three-regime behavior with the power-law dependence
with an exponent $a\sim 2$ crossing over to a slower-then-exponential convergence
at large $N_s$ (equivalently, small $\epsilon$), is typical for a wide range of
couplings and larger number of particles.

\begin{figure}[htb]
\includegraphics[width=0.99\columnwidth, keepaspectratio=True]{rho_s_N5L21m3.png} \\%
\caption{\label{fig:rho_convergence} Convergence of the sum rule \eqref{eq:sum_rule} with $N_s$, the number of states
included into the subspace \eqref{eq:subspace}. Open circles show partial sums
of $N_s$ states ordered by the overlap $|\langle \psi | \Psi_0\rangle|^2$, and filled squares
show partial sums where states are ordered by the energy, Eq.\ \eqref{eq:e0}. Dashed
line illustrates the power-law $\propto N_s^{a}$ with $a=-1.2$ and the dash-dotted line is the
power-law dependence $\propto N_s^{a}$ with $a = -1.9$. See text for discussion.
}
\end{figure}

Fig.\ \ref{fig:overlaps}(top) shows individual states of Fig.\ \ref{fig:rho_convergence}.
Specifically, for each state we show its overlap, $|\langle \psi | \Psi_0\rangle|^2$, vs energy $\varepsilon_\psi$,
relative to the energy of the in-state, $E_\mathrm{in} = \langle \Psi_0 | \widehat{H}_0 |\Psi_0\rangle$. We also
indicate special single-excitation families, similar to \eqref{eq:dominant_family},
which are parameterized by the position of the pseudo-particle, $k$, at a fixed location
of the pseudo-hole, $l$. For a finite number of host particles ($N=5$ here), these
families are not separated from the rest of the spectrum---the separation increases
with the number of particles $N$ (cf Fig. 2 of Ref.\ \cite{Burovski2014}).
At finite $N$, doubly excited states have a finite contribution to Eq.\ \eqref{eq:sum_rule},
including a set of states with the overlap $\sim 10^{-4}$, which are visually separated from the
lower-overlap part in Fig.\ \ref{fig:overlaps}(top). The contribution of
three-excitation states to the sum rule, \eqref{eq:sum_rule}, is minor: $\sim 600$ states
present in Figs.\ \ref{fig:rho_convergence}-\ref{fig:overlaps}(top) contribute $\sim 3\times 10^{-5}$
in total.

We now turn our attention to constructing the eigenstates of the full model,
Eq.\ \eqref{h1}--\eqref{h0}. We construct the full Hamiltonian matrix \eqref{h_matrix}
in the subspace $\mathcal{S}$, which is spanned by the states shown in Figs. \ref{fig:rho_convergence}
and \ref{fig:overlaps}, and numerically diagonalize it using standard LAPACK
routines. This procedure generates the set, $\mathcal{S}_f$, of eigenstates
with energy eigenvalues $\varepsilon_f$, and orthonormal eigenvectors, $|f\rangle$,
of the Hamiltonian \eqref{ham_orig} in the basis of eigenstates of Eq. \eqref{h0}:
$| f \rangle = \sum_{\psi\in\mathcal{S}} | \psi \rangle \langle \psi| f \rangle$.

\begin{figure}[hbt]
\includegraphics[width=0.99\columnwidth, keepaspectratio=True]{overlaps_N5L21_integrable.png} \\%
\includegraphics[width=0.99\columnwidth, keepaspectratio=True]{overlaps_N5L21_diag.png} \\%
\caption{\label{fig:overlaps} (top panel) Overlaps, $|\langle \psi | \Psi_0\rangle|^2$ as a function
of energy, $\varepsilon_\psi$, for eigenstates of the integrable model \eqref{h0}.
Open squares, triangles and diamonds show single-excitation families, similar to
\eqref{eq:dominant_family} with a pseudo-hole fixed at $l=2$, $1$, and $0$, respectively.\\
(bottom panel) Overlaps, $|\langle f | \Psi_0\rangle|^2$ as a function
of energy, $\varepsilon_f$, for eigenstates of the non-integrable model \eqref{h1}--\eqref{h0},
with $m/M = 1.2$, obtained by diagonalization of Eq.\ \eqref{h_matrix} in the
subspace $\mathcal{S}$ spanned by the states shown in Fig.\ \ref{fig:rho_convergence}.
In both panels we only show states with overlaps larger then $10^{-7}$ and energies
smaller then $11 E_\mathrm{in}$, for clarity.
See text for discussion.
}
\end{figure}

The bottom panel of Fig.\ \ref{fig:overlaps} shows the result of the diagonalization
procedure in the same subspace as Fig. \ref{fig:rho_convergence}. The overlaps with
the in-state are computed via $\langle \Psi_0 | f \rangle = \sum_{\psi\in\mathcal{S}} \langle \Psi_0 | \psi \rangle \langle \psi| f \rangle$.
Here we take the ratio of masses of host particles and the impurity, $m/M = 1.2$.
The overall structure is similar for other mass ratios.

While the total weight is invariant under the unitary transformation,
$\sum_{\psi\in\mathcal{S}} |\langle \psi|\Psi_0 \rangle|^2 = \sum_{f\in\mathcal{S}_f} |\langle f|\Psi_0 \rangle|^2$, it is clear
from Fig.\ \ref{fig:overlaps} that having $M \neq m$ drastically changes the
structure of relative weights of individual states. For a non-integrable model,
\eqref{h1}--\eqref{h0}, the overlaps  with in-state as a function of energy
for the majority of states cluster in a vicinity of a certain master curve,
in a stark contrast to the states of the integrable model, \eqref{h0}.
This situation resembles the typical behavior of thermalization in a generic
isolated quantum system \cite{Rigol2008}. This approach was recently used to describe
observables for the Holstein polaron, which is very close to our case \cite{Jansen2019}.
The eigenstates thermalization hypothesis was clearly demonstrated in the sense
that diagonal matrix elements are functions of the energy only.
We see that in our case for $m\neq M$  similar trend is present even for the
overlaps with the in-state.  We believe that this result might be useful for
the computations of correlation functions.

\begin{figure}[hbt]
\includegraphics[width=0.99\columnwidth, keepaspectratio=True]{NNLS_hist.png} \\%
\caption{\label{fig:NNSD} Nearest neighbor level spacing distributions for
the integrable model, $m=M$ (solid blue line), and non-integrable models with a
heavy impurity, $m/M = 0.8$ and a light impurity, $m/M = 1.2$ (dashed lines).
To set the energy scale, level spacing is given relative to the in-state energy $E_\mathrm{in}$
of the integrable model.
See text for discussion.
}
\end{figure}

To further probe the effect of $M\neq m$, we consider statistics of energy levels. For an
integrable case, $M = m$, the levels are expected to be uncorrelated (the system
is then said to be regular) \cite{Berry1977}, while for $M \neq m$ spectral
statistics is expected to be chaotic. Fig.\ \ref{fig:NNSD} shows the distribution of energy level spacing,
$\Delta\varepsilon = \varepsilon_{j+1} - \varepsilon_{j}$ for both integrable
and non-integrable cases. Here we take the subspace $\mathcal{S}$ of
Fig.\ \ref{fig:rho_convergence}, which contains $3\times 10^3$ states, and diagonalize the
Hamiltonian matrix \eqref{h_matrix} for $m/M = 1.2$ and $m/M = 0.8$. A clear feature of
Fig.\ \ref{fig:NNSD} is a reduction of distribution for
$\Delta\epsilon / E_\mathrm{in} \lesssim 10^{-3}$ for $M\neq m$---which signals
level repulsion, typical for non-integrable systems \cite{Bohigas1984}. We also note that
the shape of the distribution distinguishes between integrable and non-integrable
systems, but is similar for both light or heavy impurity ($M < m$ or $M > m$).

\begin{figure}[htb]
\includegraphics[width=0.99\columnwidth, keepaspectratio=True]{ene_mr.png} \\%
\caption{\label{fig:ene_mr} Evolution of six lowest-lying energy states with the
mass ratio, $m/M$. To set the energy scale, we scale the energies relative to
the in-state of the integrable model, $E_\mathrm{in}$.
See text for discussion.
}
\end{figure}

Fig.\ \ref{fig:ene_mr} illustrates the dependence of the ground state and five lowest
excited states with the mass ratio, $0.75 \leqslant m/M \leqslant 1.25$. The ground state is non-degenerate and
monotonic with $m/M$; excited states show the Landau-Zener behavior in the vicinity
of $m/M = 1$.

\section{Conclusions and outlook}
\label{sec:concl}

We study a model which captures the physics of an impurity injected into a degenerate 1D
Fermi gas. For equal masses of the impurity and the host particles, the model admits
an exact Bethe Ansatz solution, which we use for studying the non-integrable model
with unequal masses. We construct an exact diagonalization procedure using
truncated bases of the Bethe Ansatz states of the $m=M$ model. Our method is
non-perturbative---we do not rely on either a coupling constant or the mass difference
being small---and controllable: the accuracy of the truncation of the Hilbert space
is controlled by the ``missing weight'' of the discarded configurations, which
is, in turn, controlled by choosing the set of Bethe Ansatz states we include into
the diagonalization.
We illustrate our numerical procedure on a system of $5+1$ particles and compare
spectral properties and redistribution of spectral weights for $m=M$ and $m\neq M$.

Comparing our approach to alternative numerical methods, we note that typically
exact diagonalization and matrix-product state simulations work with
lattice models. (For parabolic traps, a truncated basis of harmonic oscillator
states can also be used \cite{Grass2015}.)
Our approach is, however, fully off-lattice---it can of course
be extended to lattice models which allow constructing a Bethe Ansatz basis to
build the diagonalization procedure on.

For few-body systems off-lattice, the Schroedinger equation can be directly solved
numerically \cite{Petrov2011}. However, this approach is limited to $N < 4$, while
our procedure can be directly extended to mesoscopic systems of a few tens of
particles, as typical in current experiments with mixtures of ultracold quantum
gases \cite{Palzer2009, Catani2012, Meinert2017}. Furthermore, our numerical results
can be useful for benchmarking and extending various approximate analytical and semi-analytical
approaches, see e.g. \cite{Doggen2013,Mehta2014,Dehkharghani2015,Petkovic2015,Das2018,Heitmann2020,Kwasniok2020,Sorout2020,Brauneis2021, Levinsen2021, Vidana2021,Doi2021,Dolgirev2021}.

\textit{Acknowledgments.---} E.B. acknowledges support of the Academic Fund Program
at the NRU Higher School of Economics (HSE) in 2020-2021 (grant No 20-01-025).
O.G. acknowledges support from the Polish National Agency for Academic Exchange (NAWA) through the Grant No. PPN/ULM/2020/1/00247