# Structure of wavefunction for interacting bosons in mean-field with random
$k$-body interactions

Priyanka Rao  N. D<EMAIL_ADDRESS>Department of
Applied Physics, Faculty of Technology and Engineering, The Maharaja Sayajirao
University of Baroda, Vadodara-390001, India

###### Abstract

Wavefunction structure is analyzed for dense interacting many-boson systems
using Hamiltonian $H$, which is a sum of one-body $h(1)$ and an embedded GOE
of $k$-body interaction $V(k)$ with strength $\lambda$. In the first analysis,
a complete analytical description of the variance of the strength function as
a function of $\lambda$ and $k$ is derived and the marker $\lambda_{t}$
defining thermalization region is obtained. In the strong coupling limit
($\lambda>\lambda_{t}$), the conditional $q$-normal density describes Gaussian
to semi-circle transition in strength functions as body rank $k$ of the
interaction increases. In the second analysis, this interpolating form of the
strength function is utilized to describe the fidelity decay after $k$-body
interaction quench and also to obtain the smooth form for the number of
principal components, a measure of chaos in finite interacting many-particle
systems. The smooth form very well describes embedded ensemble results for all
$k$ values.

##  I Introduction

It is now well established that Random Matrix Theory, due to it’s universality
[1], successfully describes the spectral as well as wavefunction properties of
isolated finite many-particle quantum systems [2]. The spectral statistics
deals only with the energy eigenvalues while the statistical properties
related to the structure of the wavefunctions can reveal different layers of
chaos and hence give profound understanding of various problems in the field
of quantum many-body chaos and thermalization, in isolated finite interacting
particle systems such as atomic nuclei, atoms, mesoscopic systems (quantum
dots, small metallic grains), interacting spin systems modeling quantum
computing core, ultra- cold atoms and quantum black holes with SYK model and
so on [2, 3, 4, 5, 6, 7, 8, 9]. To analyze the wavefunction properties, it is
very crucial to examine the so-called strength functions (also known as local
density of states) in detail, as they give information about how a particular
basis state spreads onto the eigenstates. The chaos measures like number of
principal components (NPC), information entropy, fidelity decay etc. can also
be determined by examining the general features of the strength functions [2].

The statistical properties of isolated finite many-particle quantum systems
investigated by employing random matrix ensembles are generally referred as
Gaussian ensembles (and in particular the Gaussian orthogonal ensemble (GOE))
for $m$-particle system. They involve interaction up to $m$-body in character
and are dominated by the $m$-body interactions. However, constituents of
isolated quantum systems interact via few-body interactions. Hence the concept
of embedded ensemble (EE) of $k$-body interaction, in particular EGOE($k$)
(GOE version of EE($k$)) was introduced by French and co-workers [10, 11].
These models for the particles in a mean-field and interacting via two-body
interactions ($k=2$) and their various extended versions form good models for
understanding various aspects of chaos in interacting particle systems [2] and
they are investigated in detail both for fermion systems (called EGOE(1+2))
[12, 13, 14, 15, 16, 17] as well as boson systems (called BEGOE(1+2) with ’B’
for bosons) [18, 19, 20, 21, 22, 23]. Here, with $m$ particles distributed in
$N$ single particle (sp) states, two limiting situations exist, one is the
dilute limit (defined as $m\rightarrow\infty$, $N\rightarrow\infty$ and
$m/N\rightarrow 0$) and another is the dense limit (defined by
$m\rightarrow\infty$, $N\rightarrow\infty$ and $m/N\rightarrow\infty$). In the
dilute limit, one can expect similar behavior for both fermion and boson
systems while the dense limit is feasible only for boson systems and therefore
the focus was on the dense limit in BEGOE investigations [18, 19, 20, 21, 22,
23, 24]. For EGOE(1+2) in dilute limit and for BEGOE(1+2) in dense limit, as a
function of the two-body interaction strength $\lambda$ (measured in units of
the average spacing between the one-body mean-field sp levels), exhibits three
transition or chaos markers $(\lambda_{C},\lambda_{F},\lambda_{t})$: (a) as
the two-body interaction is turned on, level fluctuations exhibit a transition
from Poisson to GOE at $\lambda=\lambda_{C}$; (b) with further increase in
$\lambda$, the strength functions make a transition from Breit-Wigner (BW)
form to Gaussian form at $\lambda=\lambda_{F}>\lambda_{C}$; and (c) beyond
$\lambda=\lambda_{F}$, there is a region of thermalization around
$\lambda=\lambda_{t}$ where the basis dependent thermodynamic quantities like
entropy behave alike. It is important to note that the transitions mentioned
above are inferred from large number of numerical calculations and they are
well verified to be valid in the bulk part of the spectrum. For further
details see [2] and references there in.

Going beyond two-body interaction, it is seen that the higher body
interactions i.e. $k>2$ play an important role in strongly interacting quantum
systems [25, 26], nuclear physics [27], quantum black holes [7, 28] and
wormholes [29] with SYK model and also in quantum transport in disordered
networks connected by many-body interactions [30, 31, 32]. Therefore, it is
necessary to extend the analysis of EE to higher $k$-body interactions in
order to understand these problems. From the previous studies, it is known
that with EGOE($k$) or (BEGOE($k$)), the eigenvalue density for a system of
$m$ fermions/bosons in $N$ sp states changes from Gaussian form to semi-circle
as $k$ changes from 2 to $m$ [2, 6, 13, 33]. Very recently, $q$-Hermite
polynomials have been employed to study spectral densities of the so-called
SYK model [34, 35] and quantum spin glasses [36], along with studying the
strength functions and fidelity decay (also known as survival or return
probability) in EE, both for fermion as well as boson systems [33]. The smooth
form of eigenvalue density can be given by the so-called $q$-normal
distribution $f_{qN}$ and formulas for parameter $q$ in terms of $m$, $N$ and
$k$ are derived for fermionic and bosonic EE($k$) in [33] which explain the
Gaussian to semi-circle transition in spectral densities, strength functions
and fidelity decay in many-body quantum systems as a function of rank $k$ of
interactions. Recently, the lower-order bivariate reduced moments of the
transition strengths are examined for the action of a transition operator on
the eigenstates generated by EGOE($k$) and it is shown that the ensemble
averaged distribution of transition strengths follows a bivariate $q$-normal
distribution $f_{biv-qN}$ and a formula for NPC in the transition strengths
from a state is obtained [37]. Very recently, analytical formulas for the
lowest four moments of the strength functions for fermion systems modeled by
EGOE(1+$k$) are derived and it is shown that the conditional $q$-normal
density $f_{CqN}$ can be used to represent strength functions in the strong
coupling limit [38]. One can expect similar behavior for isolated finite
interacting boson systems with $k$-body interactions in the dense limit. The
purpose of the present letter is firstly to demonstrate that in strong
coupling domain (in the thermalization region), the strength functions indeed
can be represented by the conditional $q$-normal distribution $f_{CqN}$ in the
dense interacting boson systems interacting via $k$-body interaction.
Secondly, using $f_{CqN}$ form and parameters that enter in this form,
fidelity decay is described in BEGOE(1+$k$) and an analytical formula for NPC
is derived.

The Letter is organized as follows. We briefly introduce BEGOE(1+$k$) and
$q$-Hermite polynomials along with their generating function and conditional
$q$-normal distribution in Section II. The numerical results of the variation
of parameter $q$ as a function of $k$-body interaction strength $\lambda$ in
BEGOE(1+$k$) are presented in Section III. Also the formula of $q$ for
BEGOE($k$) is given for the sake of completeness, even though it is clearly
given in [6, 33]. Further, a complete analytical description of the variance
of the strength function, in terms of the correlation coefficient $\zeta$, for
BEGOE(1+$k$) is given and ($m$,$N$,$k$) dependence of marker $\lambda_{t}$ is
derived. In Section IV, the results for the variation of strength function, in
the strong coupling domain ($\lambda>>\lambda_{t}$), are presented as a
function of body rank $k$ and ensemble averaged results are compared with
smooth forms given by $f_{CqN}$. In Section V the interpolating form $f_{CqN}$
for the strength function is utilized to describe the fidelity decay after
random $k$-body interaction quench in BEGOE(1+$k$) in the thermalization
region. Further, two parameter ($\zeta$ and $q$) analytical formula for NPC is
derived as a function of energy for $k$-body interaction and tested with
numerical embedded ensemble results in Section VI. Finally, the concluding
remarks are given in section VII.

##  II Preliminaries

###  II.1 Embedded bosonic ensembles - BEGOE(1+$k$)

Consider $m$ spinless bosons distributed in $N$ degenerate sp states
interacting via $k$-body ($1\leq k\leq m$) interactions. Distributing these
$m$ bosons in all possible ways in $N$ sp states generates many-particle basis
of dimension $d={N+m-1\choose{m}}$. The $k$-body random Hamiltonian $V(k)$ is
defined as,

$V(k)=\displaystyle\sum_{k_{a},k_{b}}V_{k_{a},k_{b}}B^{\dagger}(k_{a})B(k_{b})\;.$
(1)

Here, operators $B^{\dagger}(k_{a})$ and $B(k_{b})$ are $k$-boson creation and
annihilation operators. They obey the boson commutation relations.
$V_{k_{a},k_{b}}$ are the symmetrized matrix elements of $V(k)$ in the
$k$-particle space with the matrix dimension being $d_{k}={N+k-1\choose k}$.
They are chosen to be randomly distributed independent Gaussian variables with
zero mean and unit variance, in other words, $k$-body Hamiltonian is chosen to
be a GOE. BEGOE($k$) is generated by action of $V(k)$ on the many-particle
basis states. Due to $k$-body nature of interactions, there will be zero
matrix elements in the many-particle Hamiltonian matrix, unlike a GOE. By
construction, we have a GOE for the case $k=m$. For further details about
these ensembles, their extensions and applications, see [2, 39, 40] and
references therein.

In realistic systems, bosons also experience mean-field generated by presence
of other bosons in the system and hence, it is more appropriate to model these
systems by BEGOE($1+k$) defined by,

$H=h(1)+\lambda V(k)$ (2)

Here, the one-body operator $h(1)=\sum_{i=1}^{N}\epsilon_{i}n_{i}$ is
described by fixed sp energies $\epsilon_{i}$; $n_{i}$ is the number operator
for the $i$th sp state. The parameter $\lambda$ represents the strength of the
$k$-body interaction and it is measured in units of the average mean spacing
of the sp energies defining $h(1)$. In this analysis, we have employed fixed
sp energies $\epsilon_{i}=i+1/i$ in defining the mean-field Hamiltonian
$h(1)$. As the dense limit is more interesting for bosons, for numerical
study, we have chosen $N=5$, $m=10$ with space dimensionality of $d=1001$ and
varied $k$ from 2 to $m$. It is now known that in nuclear reactions and
strongly interacting quantum systems $k=2,3,4$ are of physical importance[7,
25, 26]. However for the sake of completeness, to study the generic features
of embedded ensembles and the possibility of higher $k$ becoming prominent, we
address $k=2$ to $m$.

###  II.2 $q$-Hermite polynomials and conditional $q$-normal distribution

The $q$-Hermite polynomials were first introduced by L. J. Rogers in
Mathematics. Consider $q$ numbers $[n]_{q}$ defined as
$\left[n\right]_{q}=(1-q)^{-1}(1-q^{n})$. Then, $[n]_{q\rightarrow 1}=n$, and
$[n]_{q}!=\Pi^{n}_{j=1}[j]_{q}$ with $[0]_{q}!=1$. Now, $q$-Hermite
polynomials $H_{n}(x|q)$ are defined by the recursion relation [41],

$x\,H_{n}(x|q)=H_{n+1}(x|q)+\left[n\right]_{q}\,H_{n-1}(x|q)$ (3)

with $H_{0}(x|q)=1$ and $H_{-1}(x|q)=0$. Note that for $q=1$, the $q$-Hermite
polynomials reduce to normal Hermite polynomials (related to Gaussian) and for
$q=0$ they will reduce to Chebyshev polynomials (related to semi-circle).
Importantly, $q$-Hermite polynomials are orthogonal within the limits $\pm
2/\sqrt{1-q}$, with the $q$-normal distribution $f_{qN}(x|q)$ as the weight
function defined by [37],

$f_{qN}(x|q)=\displaystyle\frac{\sqrt{1-q}}{2\pi\sqrt{4-(1-q)x^{2}}}\displaystyle\prod_{i=0}^{\infty}(1-q^{i+1})[(1+q^{i})^{2}-(1-q)q^{i}x^{2}].$
(4)

Here, $-2/\sqrt{1-q}\leq x\leq 2/\sqrt{1-q}$ and $q\in[0,1]$. Note that
$\int_{s(q)}f_{qN}(x|q)\;dx=1$ over the range
$s(q)=(-2/\sqrt{1-q},2/\sqrt{1-q})$. It is seen that in the limit
$q\rightarrow 1$, $f_{qN}(x|q)$ will take Gaussian form and in the limit $q=0$
semi-circle form. Now the bivariate $q$-normal distribution $f_{biv-
qN}(x,y|\zeta,q)$ is defined as follows [37, 42],

$\begin{array}[]{rcl}f_{biv-
qN}(x,y|\zeta,q)&=&f_{qN}(x|q)f_{CqN}(y|x;\zeta,q)\\\ \\\
&=&f_{qN}(y|q)f_{CqN}(x|y;\zeta,q)\end{array}$ (5)

where $\zeta$ is the bivariate correlation coefficient and the conditional
$q$-normal densities, $f_{CqN}$ can be given as,

$\begin{array}[]{rcl}f_{CqN}(x|y;\zeta,q)&=&f_{qN}(x|q)\;\displaystyle\prod_{i=0}^{\infty}\frac{(1-\zeta^{2}q^{i})}{h(x,y|\zeta,q)};\\\
\\\
f_{CqN}(y|x;\zeta,q)&=&f_{qN}(y|q)\;\displaystyle\prod_{i=0}^{\infty}\frac{(1-\zeta^{2}q^{i})}{h(x,y|\zeta,q)};\\\
\\\ h(x,y|\zeta,q)&=&(1-\zeta^{2}q^{2i})^{2}-(1-q)\zeta
q^{i}(1+\zeta^{2}q^{2i})xy+(1-q)\zeta^{2}(x^{2}+y^{2})q^{2i}.\end{array}$ (6)

The $f_{CqN}$ and $f_{biv-qN}$ are normalized to 1 over the range $s(q)$,
which can be inferred from the following property,

$\int_{s(q)}H_{n}(x|q)f_{CqN}(x|y;\zeta,q)\;dx=\zeta^{n}H_{n}(y|q).$ (7)

The first four moments of the $f_{CqN}$ can be given [38] as,

$\begin{array}[]{rcl}\text{Centroid}&=&\zeta y,\\\ \\\
\text{Variance}&=&1-\zeta^{2}\;,\\\ \\\
\text{Skewness,}\;\gamma_{1}&=&-\displaystyle\frac{\zeta(1-q)y}{\sqrt{1-\zeta^{2}}},\\\
\\\
\text{Excess,}\;\gamma_{2}&=&(q-1)+\displaystyle\frac{\zeta^{2}(1-q)^{2}y^{2}+\zeta^{2}(1-q^{2})}{(1-\zeta^{2})}\;.\end{array}$
(8)

Recently, it is shown that generating function for $q$-Hermite polynomials
describes Gaussian to semi-circle transition in the eigenvalue density as $k$
changes from from $1$ to $m$ in spectral densities using $k$-body EGOE and
their Unitary variants EGUE, both for fermion and boson systems [33]. Very
recently, in the strong coupling domain the lowest four moments of the
strength function for $k$-body fermionic embedded ensemble are obtained and it
is shown that they are essentially same as that of $f_{CqN}$ [38]. Therefore,
one can use $f_{CqN}$ distribution to represent the smooth forms of the
strength functions and analyze the wavefunction structure in quantum many-body
systems with $k$-body interactions. With this, the width of $f_{CqN}$ (and
also of the strength fucntion) is related to the correlation coefficient
$\zeta$ by Eq. (8). In the next section, we will present our results for the
variation of parameter $q$ and the correlation coefficient $\zeta$ as a
function of $k$-body interaction strength $\lambda$ in BEGOE(1+$k$). Also, a
complete analytical description of $\zeta$, in terms of $N,m$,$k$ and
$\lambda$, for BEGOE(1+$k$) is given.

##  III Parameter dependence of $q$ and $\zeta$ : results for BEGOE(1+$k$)

###  III.1 Formula of $q$-parameter

It has already been demonstrated that the state density for EE($k$)(and also
EE(1+$k$)) in general exhibits Gaussian to semi-circle transition as $k$
increases from $1$ to $m$ [17]. This is now well verified in many numerical
calculations and analytical proofs obtained via lower order moments [2, 6, 9,
20, 39, 43]. Figure 1(a) represents ensemble averaged state density obtained
for a 100 member BEGOE(1+$k$) ensemble with $m=10$ bosons distributed in $N=5$
sp states and the body rank of interaction changing from $k$ = 2 to 10. In
these calculations, the eigenvalue spectrum for each member of the ensemble is
first zero centered ($\epsilon_{H}$ is centroid) and scaled to unit width
($\sigma_{H}$ is width) and then the histograms are constructed. The results
clearly display transition in the spectral density from Gaussian to semi-
circle form as $k$ changes from 2 to $m=10$. With $E$ as zero centered and
using $x=E/\sigma_{H}$, the numerical results are compared with the normalized
state density $\rho(E)=d\;f_{qN}(x|q)$ with
$\epsilon_{H}-\frac{2\sigma_{H}}{\sqrt{1-q}}\leq
E\leq\epsilon_{H}+\frac{2\sigma_{H}}{\sqrt{1-q}}$. Here the parameter $q$ is
computed using the formula, valid for BEGOE($k$)(i.e. $H=V(k)$), given in
[33],

$\begin{array}[]{l}q_{V(k)}\sim\displaystyle\binom{N+m-1}{m}^{-1}\displaystyle\sum_{\nu=0}^{\nu_{max}=\min[k,m-k]}\;\displaystyle\frac{X(N,m,k,\nu)\;d(g_{\nu})}{\left[\Lambda^{0}(N,m,k)\right]^{2}}\;;\\\
\\\ X(N,m,k,\nu)=\Lambda^{\nu}(N,m,m-k)\;\Lambda^{\nu}(N,m,k)\;;\\\ \\\
\Lambda^{\nu}(N,m,r)=\displaystyle\binom{m-\nu}{r}\;\displaystyle\binom{N+m+\nu-1}{r}\;,\\\
\\\
d(g_{\nu})=\displaystyle\binom{N+\nu-1}{\nu}^{2}-\displaystyle\binom{N+\nu-2}{\nu-1}^{2}\;.\end{array}$
(9)

In the strong coupling domain, one can also apply Eq.(9) to BEGOE(1+$k$), as
the $k$-body part of the interaction is expected to dominate over one-body
part. One can see that the ensemble averaged results in Figure 1(a) are in
excellent agreement with the smooth forms obtained using $f_{qN}$. With
$\lambda=0$ in Eq.(2) i.e. one-body part $h(1)$ only, the analytical formula
of $q$ for bosons, based on trace propagation method [44], can be given as,

$\begin{array}[]{lcl}q_{h(1)}&=&{\langle h(1)^{4}\rangle}^{m}-2\\\ \\\
&=&\displaystyle{\\{\frac{3(m-1)N(1+N)(1+m+N)}{m(2+N)(3+N)(m+N)}-2\\}}\\\
&&+\displaystyle{\frac{m^{2}+(N+m)^{2}+(N+2m)^{2}}{m(N+m)}\frac{\sum_{i=1}^{N}\tilde{\epsilon_{i}}^{4}}{(\sum_{i=1}^{N}\tilde{\epsilon_{i}}^{2})^{2}}}.\\\
\end{array}$ (10)

Here, ${\langle h(1)^{4}\rangle}^{m}$ is the reduced fourth moment of one-body
part and $\tilde{\epsilon_{i}}$ are the traceless sp energies of $i$’th state.
With $H=h(1)$ and uniform sp energies $\epsilon_{i}=i$, Eq.(10) gives $q=0.71$
for ($m=5,N=10$) and $q=0.68$ for ($m=10,N=5$). While with sp energies
$\epsilon_{i}=i+1/i$, used in the present study, one obtains $q=0.68$ for
($m=5,N=10$) and $q=0.63$ for ($m=10,N=5$). Figure 1(b) shows variation of
$q_{h(1)}$ as a function of $N$ for various values of $m/N$. Here, sp energies
$\epsilon_{i}=i+1/i$ are used. It can be clearly seen that in the dense limit
($m\rightarrow\infty$, $N\rightarrow\infty$ and $m/N\rightarrow\infty$),
$q_{h(1)}\rightarrow 1$. In the dilute limit ($m\rightarrow\infty$,
$N\rightarrow\infty$ and $m/N\rightarrow 0$), similar variation in $q_{h(1)}$
can be observed due to $m\leftrightarrow N$ symmetry between the dense limit
and the dilute limit as identified in [18, 44]. Furthermore, the variation of
parameter $q$ is also studied as the interaction strength $\lambda$ varies in
BEGOE(1+$k$) for a fixed body rank $k$. Here, the ensemble averaged value of
$q$ is computed for a system of 100 member BEGOE(1+$k$) ensemble with $m=10$
bosons in $N=5$ sp states and results are shown in Figure 1(c). $q$ estimates
are also shown in the figure by horizontal marks for $H=h(1)$ and $H=V(k)$ on
left and right vertical axes respectively. One can see that for very small
values of $\lambda$, ensemble averaged $q$ values are found very close to
$q_{h(1)}$ for all body rank $k$. While for a sufficiently large $\lambda$,
where $k$-body part dominates over one-body part and ensemble averaged $q$
values reach corresponding $q_{V(k)}$ given by Eq.(9). From the variation of
ensemble averaged $q$ values in Figure 1(c), one can see that the shape of the
state density takes intermediate form between Gaussian to semi-circle as
$\lambda$ changes in BEGOE(1+$k$) for a fixed $k$. Therefore, the $q$-normal
distribution $f_{qN}$ formula can be used to describe the transition in the
state density with any value of $\lambda$ and $k$ in BEGOE(1+$k$).

|
---|---

Figure 1: (a) Histograms represent the state density vs. normalized energy $E$
results of the spectra of a 100 member BEGOE($1+k$) ensemble with $m=10$
bosons in $N=5$ sp states for different $k$ values. The strength of
interaction $\lambda=0.5$ is chosen and in the plots $\int\rho(E)dE=d$.
Ensemble averaged state density histogram is compared with $q$-normal
distribution (continuous black curves) given by $f_{qN}(x|q)$ with the
corresponding $q$ values given by Eq. (9). (b) $q_{h(1)}$ vs. $N$ for various
values of $m/N$. $q_{h(1)}$ is obtained using Eq. (10) with sp energies
$\epsilon_{i}=i+1/i$. Dense limit curve corresponds to the result with
$m/N=1000$. (c) Ensemble averaged $q$ vs. $\lambda$ for a 100 member
BEGOE(1+$k$) ensemble with $m=10$ bosons in $N=5$ sp states for different $k$
values. The horizontal black mark on left $q$-axis indicates $q$ estimate for
$H=h(1)$ given by Eq. (10), while the colored marks on right $q$-axis
represent the $q$ values, given by Eq. (9), for corresponding $k$-body rank
with $H=V(k)$. See text for more details.

###  III.2 Formula of $\zeta$

The parameter $\zeta$, which is the correlation coefficient between full
Hamiltonian $H$ and the diagonal part $H_{\text{dia}}$ of the full
Hamiltonian, is related to the width $\sigma_{F}$ of the strength functions,
given by,

$\zeta=\sqrt{1-\displaystyle\frac{\sigma_{H_{\text{off-
dia}}}^{2}}{\sigma_{H}^{2}}}=\sqrt{1-\sigma_{F}^{2}},\;\;\;\;\sigma_{F}=\displaystyle\frac{\sigma_{H_{\text{off-
dia}}}}{\sigma_{H}}$ (11)

In the above equation, $\sigma_{H}^{2}$ and $\sigma_{H_{\text{off-dia}}}^{2}$
are variances of the eigenvalue distribution using full Hamiltonian and by
taking all diagonal matrix elements as zero, respectively. Since $\zeta$ and
$\sigma_{F}$ are simply related as $\sigma_{F}^{2}=1-\zeta^{2}$, here the
discussion is in terms of $\zeta$. For BEGOE(1+$k$) ensemble, analytical
expression for $\zeta$ based on the method of trace propagation can be derived
as follows. For $H=V(k)$ i.e. with all sp energies as degenerate, it is known
that [20],

$\begin{array}[]{rcl}\sigma_{H=V(k)}^{2}&=&\displaystyle
T(N,m,k)\binom{N+k-1}{k}^{-1}\;\sum_{\alpha,\beta}\overline{w^{2}_{\alpha\beta}}\;,\\\
\\\ T(N,m,k)&=&\displaystyle\Lambda^{0}(N,m,k)/\binom{N+k-1}{k}\;.\end{array}$
(12)

Here, $\alpha$ and $\beta$ denote $k$-particle states. In $k$-particle space,
the $H$ matrix is GOE. Therefore, the $k$-particle matrix elements
$w_{\alpha\beta}$ are Gaussian random variates with zero mean and unit
variance. The variance of diagonal matrix elements is
$\overline{w^{2}_{\alpha\alpha}}=2$ while that of off-diagonal matrix elements
is $\overline{w^{2}_{\alpha\beta}}=1$ for ($\alpha\neq\beta$). With this,

$\sigma_{H=V(k)}^{2}=T(N,m,k)\;\binom{N+k-1}{k}^{-1}\left\\{2\times\text{no-
dia}+2\times\text{no-offdia}\right\\},$ (13)

here the number of independent diagonal $k$-body matrix elements is ’no-
dia’$=\binom{N+k-1}{k}$ and that of off-diagonal is ’no-
offdia’$=\frac{1}{2}\binom{N+k-1}{k}\\{\binom{N+k-1}{k}-1\\}$. Similarly,
$\sigma_{H_{\text{off-dia}}}$ is given by removing the contribution of
diagonal $k$-body matrix elements from the above equation. Then using Eq.(11)
for $H=V(k)$,

$\zeta^{2}=\frac{4}{{N+k-1\choose k}+1}\;.$ (14)

Here, it can be immediately seen that $\zeta^{2}$ is independent of $m$ for
BEGOE($k$). In the dense limit with $N\rightarrow\infty$ and
$m\rightarrow\infty$, $\sigma_{F}\rightarrow 1$ giving $\zeta\rightarrow 0$ as
was suggested in [21]. Also, with $k<<m$, $\zeta^{2}\propto 1/N^{k}$. Using
$m\leftrightarrow N$ symmetry between the dense limit and the dilute limit
formula [18, 44], we have $\zeta^{2}\propto 1/m^{k}$ in the dilute limit and
this result is in agreement with [38]. Going further, with inclusion of one-
body part defined by the external sp energies ($\epsilon_{i}$), and with
$H=h(1)+\lambda V(k)$, we have

$\begin{array}[]{rcl}\sigma_{H}^{2}&=&\sigma_{h(1)}^{2}+\lambda^{2}\;\sigma_{V(k)}^{2},\\\
\\\
&=&\frac{m(N+m)}{N(N+1)}\;\sum\tilde{\epsilon_{i}}^{2}+\lambda^{2}\;\sigma_{V(k)}^{2}.\end{array}$
(15)

The analytical expression for $\zeta^{2}$ can be given by,

$\zeta^{2}=\frac{\frac{m(N+m)}{N(N+1)}\;\sum\tilde{\epsilon_{i}}^{2}+2\;\lambda^{2}\;T(N,m,k)}{\frac{m(N+m)}{N(N+1)}\;\sum\tilde{\epsilon_{i}}^{2}+\lambda^{2}\;T(N,m,k)\;\\{1+\binom{N+k-1}{k}\\}}\;.$
(16)

In the above equation, the contribution from the diagonal part of $V(k)$ is
also included into the numerator term. The analytical expression for
$\zeta^{2}$ given by Eq.(16) is tested with the numerical ensemble averaged
results obtained using a 100 member BEGOE(1+$k$) ensemble with $(m=10,N=5)$.
The results of $\zeta^{2}$ as a function of $k$-body interaction strength
$\lambda$ for different body rank $k$ are presented in Figure 2. The black
smooth curve in each plot is obtained using Eq.(16) with fixed sp energies
employed in the present study. It can be seen from the results that agreement
between the ensemble averaged values (red solid circles) and the smooth forms
obtained by Eq.(16) is very good for all $k$ values. Small difference with
large $\lambda$, for $k<5$, is due to neglect of induced sp energies. The
contribution of induced sp energies reduces as $\lambda$ and $k$ increases.
One can see from the results shown in Figure 2 that the width of the strength
function is strongly dependent on $\lambda$. For $\lambda\rightarrow 0$,
$\zeta^{2}\rightarrow 1$ for all $k$ and the strength functions are known to
be represented by $\delta$ functions. With increase in $\lambda$
i.e.$\lambda\geq\lambda_{C}$, the strength functions are known to be described
by the Briet-Wigner (Lorentz) form. With further increase in
$\lambda>>\lambda_{F}$, $\zeta^{2}$ goes on decreasing smoothly leading to a
fully chaotic domain giving the Gaussian or semi-circle or intermediate to
Gaussian and semi-circle character of the strength functions depending upon
the values of $\lambda$ and $k$. One can also observe the BW to Gaussian to
semi-circle transition in strength functions by changing both $\lambda$ and
$k$. Therefore, it is possible to have a shape intermediate to BW and semi-
circle for some values of $\lambda$ and $k$ [45].

For two-body interaction, the thermodynamic region $\lambda=\lambda_{t}$ can
be determined using the condition $\zeta^{2}=0.5$ [23, 46]; i.e. the spreading
produced by one-body part and two-body part are equal. Similarly, one can
obtain marker $\lambda_{t}$ for $k$-body interactions in presence of mean
field by considering the spreading produced by one-body part and $k$-body part
equal in Eq.(16). Solving it for $\lambda$, ($m$, $N$, $k$) dependence of
marker $\lambda_{t}$ is given by,

$\lambda_{t}=\sqrt{\frac{m(N+m)\;\sum\tilde{\epsilon_{i}}^{2}}{N(N+1)\Lambda^{0}(N,m,k)(1-3\;\binom{N+k-1}{k}^{-1})}}\;\;.$
(17)

Figure 3 shows the variation of marker $\lambda_{t}$ in dense boson systems
with BEGOE(1+$k$) as a function of $N$ for the fixed sp energies used in the
present study. The results are shown for body rank values $k=2,3$ and $4$, and
with $m/N=2$ and $5$. From the results one can clearly see that $\lambda_{t}$
decreases as the rank of the interaction $k$ increases. Hence, the
thermalization sets in faster as the rank of interaction $k$ increases.

Recently, using $k$-body embedded ensembles both for fermions and bosons, it
is demonstrated that in the thermalization region ($\lambda\geq\lambda_{t}$),
shape of the strength functions changes from Gaussian to semi-circle for the
states close to the center of the spectrum as the rank of the interaction $k$
increases and they can be well represented by $f_{qN}$ form for all $k$ values
in $V(k)$ [33]. The strength functions are symmetrical in $E$ near the center
of the spectrum as is the result with $f_{qN}$. However, it is seen in some
calculations with $k=2$ that the strength functions become asymmetrical in $E$
as one moves away from the center [24]. This feature can be incorporated by
representing strength function using $f_{CqN}$ which can not be generated by
$f_{qN}$. This will be verified with a numerical example in the next section
and more importantly, a single interpolating function $f_{CqN}$, in terms of
parameters $q$ and $\zeta$, is considered for describing Gaussian to semi-
circle transition in the strong coupling domain as the body rank $k$ in
BEGOE(1+$k$) is changed.

Figure 2: Ensemble averaged $\zeta^{2}$ (red solid circles) as a function of
interaction strength $\lambda$, calculated for BEGOE(1+$k$) ensemble with
$N=5,m=10$ example, are shown for different $k$ values. The smooth black
curves are due to Eq.(16) using fixed sp energies $\epsilon_{i}=i+1/i$
employed in the present study. Figure 3: Variation of marker $\lambda_{t}$ as
a function of $N$ for dense boson systems with BEGOE(1+$k$). Results are shown
for various values of ($k$,$m/N$) using Eq.(17).

##  IV Strength function

Given $m$-particle basis state $\left|{\kappa}\right\rangle$, the diagonal
matrix elements of $m$-particle Hamiltonian $H$ are denoted as energy
$\xi_{{\kappa}}$, so that $\xi_{{\kappa}}=\langle{\kappa}|H|{\kappa}\rangle$.
The diagonalization of the full matrix $H$ gives the eigenstates
$\left|E_{i}\right\rangle$ with eigenvalues $E_{i}$, where
$\left|{\kappa}\right\rangle=\sum_{i}C_{{\kappa}}^{i}\left|E_{i}\right\rangle$.
The strength function that corresponds to the state
$\left|{\kappa}\right\rangle$ is defined as
$F_{\xi_{\kappa}}(E)=\sum_{i}{|C_{{\kappa}}^{i}|}^{2}\;\delta(E-E_{i})$. In
the present study, we take the $\left|{\kappa}\right\rangle$ states to be the
eigenstates of $h(1)$. In order to get an ensemble averaged form of the
strength functions, the eigenvalues $E_{i}$ are scaled to have zero centroid
and unit variance for the eigenvalue distribution. The ${\kappa}$-energies,
$\xi_{{\kappa}}$, are also scaled similarly. Now, for each member, all
${|C_{{\kappa}}^{i}|}^{2}$ are summed over the basis states ${\kappa}$ with
energy $\xi$ in the energy window $\xi\pm\Delta$. Then, the ensemble averaged
$F_{\xi}(E)$ vs. $E$ are constructed as histograms by applying the
normalization condition $\int_{s(q)}F_{\xi}(E)\;dE=1$. In Figure 4, histograms
represent ensemble averaged $F_{\xi}(E)$ results for all body rank $k$ values
with $\lambda=0.5$ using a 250 member BEGOE(1+$k$) ensemble with $m=10$ and
$N=5$ system. The strength function plots are obtained for $\xi=0.0,\pm 1.0$
and $\pm 2.0$. The value of $k$-body interaction strength is chosen such that
$\lambda>>\lambda_{t}$, i.e. the system exists in the region of thermalization
[9, 23]. The histograms, representing BEGOE(1+$k$) results of strength
functions, are compared with the conditional $q$-normal density function as
given by,

$F_{\xi}(E)=f_{CqN}(x=E|y=\xi;\zeta,q).$ (18)

The smooth black curves in Figure 4 for each $k$ are obtained via Eq.(18)
using corresponding ensemble averaged $\zeta$ and $q$ values. With
$\lambda>>\lambda_{t}$, $\zeta^{2}<<1/2$, the $q$ value in Eq.(18) can fairly
be given by Eq.(9) [38]. The results in Figure 4 clearly show very good
agreement between the numerical histograms and continuous black curves for all
body rank $k$. The $F_{\xi}(E)$ results for $\xi=0$ are given in Figure 4(a)
which clearly demonstrate that the strength functions are symmetric and also
exhibit a transition from Gaussian form to semi-circle as $k$ changes from $2$
to $m=10$. The smooth form given by Eq.(18) using the conditional $q$-normal
density function interpolates this transition very well. Going further,
$F_{\xi}(E)$ results for $\xi\neq 0$ are shown in Figures 4(b) and 4(c). One
can see that $F_{\xi}(E)$ results are asymmetrical about $E$ as demonstrated
earlier [24]. Also, $F_{\xi}(E)$ are skewed more in the positive direction for
$\xi>0$ and skewed more in the negative direction for $\xi<0$ and their
centroids vary linearly with $\xi$. We have also computed the first four
moments (centroid, variance, skewness ($\gamma_{1}$) and excess
($\gamma_{2}$)) of the strength function results shown in Figure 4 for the
body rank $k$ going from $2$ to $m=10$. Figure 5 represents results for
centroid, $\gamma_{1}$ and $\gamma_{2}$ for various values of $\xi$. As
discussed earlier in Section III, the variance of the strength functions is
independent of $\xi$ and simply related to correlation coefficient; for more
details, see results of $\zeta^{2}$ (Figure 2). From the numerical results
obtained for strength functions (Figure 4) along with results of lower order
moments (Figure 5), one can clearly see that in the thermodynamic domain, the
strength functions of dense interacting many-boson systems, with $k$-body
interaction, follow the conditional $q$-normal distribution $f_{CqN}$. The
results are also consistent with the analytical forms derived in [38].

---

|
---|---

Figure 4: Strength function vs. normalized energy $E$ for a system of $m=10$ bosons in $N=5$ sp states with $\lambda=0.5$ for different $k$ values in BEGOE(1+$k$) ensemble. An ensemble of 250 members is used for each $k$. Strength function plots are obtained for (a) $\xi=0$ (purple histogram) , (b) $\xi=-1.0$ (blue histogram) and $1.0$ (red histogram) and (c) $\xi=-2.0$ (blue histogram) and $2.0$ (red histogram). In the plots $\int F_{\xi}(E)dE=1$. The continuous black curves are due to fitting with $f_{CqN}$ given by Eq. (18) using $q$ and $\zeta$ values obtained by Eq. (9) and Eq. (11), respectively. See text for more details. |  |
---|---|---
Figure 5: Ensemble averaged (a) Centroid, (b) $\gamma_{1}$ and (c)
$\gamma_{2}$ as a function of body rank $k$ for the strength function results
presented in Figure 4. Results are shown for various values of $\xi$.

In the study of thermalization and relaxation dynamics of an isolated finite
quantum system after a random interaction quench, strength functions play an
important role. Having tested that in the thermodynamic region with
$\lambda>>\lambda_{t}$, ensemble averaged strength functions of dense boson
systems with $k$-body interaction can be represented by smooth forms given by
$f_{CqN}$, we will now utilize these interpolating forms, in the coming
sections, to study fidelity decay and NPC in dense boson systems with $k$-body
interaction.

##  V Fidelity decay after an interaction quench

Fidelity decay or return probability of a quantum system after a sudden quench
is an important quantity in the study of relaxation of a complex (chaotic)
system to an equilibrium state. Let’s say the system is prepared in one of the
eigenstates ($\psi(0)=\left|{\kappa}\right\rangle$) of the mean-field
Hamiltonian $H=h(1)$. With the quench at $t=0$ by $\lambda V(k)$, the system
evolves unitarily with respect to $H\rightarrow h(1)+\lambda V(k)$ and the
state changes after time $t$ to
$\psi(t)=\left|{\kappa}(t)\right\rangle=\exp(-iHt)\left|{\kappa}\right\rangle$.
Then, the probability to find the system in it’s initial unperturbed state
after time $t$, called fidelity decay, is given by,

$\begin{array}[]{lll}W_{0}(t)&=&|\left\langle\psi(t)|\psi(0)\right\rangle|^{2}=\left|\sum_{E}\left[C_{k}^{E}\right]^{2}\exp-
iEt\right|^{2}\\\ \\\ &=&\int F_{\xi}(E)\exp-iEt\;dE\\\ \\\
&=&\int_{s(q)}f_{CqN}(E|\xi;\zeta,q)\exp-iEt\;dE\;.\end{array}$ (19)

Thus, fidelity is the Fourier transform in energy of the strength function;
this is valid for times not very short or very long. In the thermalization
region, the form of $F_{\xi}(E)$ is Gaussian for $k=2$ while it is semi-circle
for $k=m$. These two extreme situations are recently studied, both
analytically [47] as well as numerically [48, 49, 50]. The formula for
$W_{0}(t)$ can be given in terms of width of $\lambda V(k)$ scaled by
$\sigma_{H}$. Clearly, following the results of the previous section,
$f_{CqN}$ can be used to obtain $W_{0}(t)$ generated by BEGOE(1+$k$). As
analytical formula for the Fourier transform of $f_{CqN}$ is not available,
therefore we evaluated Eq.(19) numerically. Figure 6 shows results for
$W_{0}(t)$ (red solid circles) for a 100 member BEGOE(1+$k$) ensemble with
$m=10$, $N=5$ and $\lambda=0.5$ for various $k$ values and they are compared
with numerical Fourier transform (black smooth curves) of Eq.(18). Here, we
have used normalized eigenenergies in the computation of $W_{0}$ and therefore
the time $t$ is measured in the units of $1/\sigma_{H}$. It is clear from the
results that the Fourier transform of $f_{CqN}$ describes the short-time
behavior nicely and also captures the positions of the oscillations. The
results generated here are consistent with the reported results in [33],
obtained using $f_{qN}$ form for $F_{\xi}(E)$.

It is known that in the strong interaction domain, the decrease in $W_{0}$
(for $k=2$) follows quadratic in time and this Gaussian decrease can last for
a quite large time and after that an exponential one emerges [51]. The
transition time depends on the ratio of the spectral width and the square of
the second moment of strength fucntion ($\sigma_{F}^{2}$). As here
$\lambda>>\lambda_{t}$, $\zeta^{2}\rightarrow 0$ giving $\sigma_{F}^{2}\approx
1$, $t$ is in $1/\sigma_{H}$ units and the spectral width will be in
$\sigma_{H}$ units. Therefore, the results in Figure 6 describe $W_{0}$ nicely
for short time and the standard exponential decrease for long time for $k=2$
seems absent. The long time behavior of fidelity decay is of great interest as
it is expected that $W_{0}$ surely demonstrates a power-law behavior i.e.
$W_{0}(t)\propto t^{-\gamma}$ with $\gamma\geq 2$ implying thermalization
[52], no matter how fast the decay may initially be. As shown in [52], the
power-law behavior appears due to the fact that the energy spectrum is bounded
from both the ends. This condition is essentially satisfied by $f_{CqN}$.
Therefore, it is important to analyze the long-time behavior of fidelity decay
for embedded ensembles first to establish it’s universality and second to test
whether it can be explained with the use of $f_{CqN}$. These are open
questions.

Figure 6: Fidelity decay $W_{0}(t)$ as a function of time for a 100 member
BEGOE(1+$k$) ensemble with $N=5$ and $m=10$ represented by the red solid
circles; the $\psi(0)$ here corresponds to middle states of $h(1)$ spectrum.
red Here $t$ is measured in the units of $\sigma_{H}^{-1}$. The black smooth
curves are obtained by taking numerical Fourier transform of the strength
functions represented by Eq.(18).

In the study of fidelity decay, strength function with $\xi=0$ is involved.
However, the statistical properties, related to wavefunction structure, namely
NPC and $S^{\text{info}}$ can be written as integrals involving strength
functions over all $\xi$ energies. Very recently, an integral formula for NPC
in the transition strengths from a state as a function of energy for fermionic
EGOE($k$) using the bivariate $q$-normal form is presented in [37]. In the
past, the smooth forms, for NPC and $S^{\text{info}}$, were derived in terms
of energy and correlation coefficient $\zeta$ for two-body interaction [53].
In the next section, we present our results for NPC and $S^{\text{info}}$
using $f_{CqN}$ forms for the strength functions and compare with those for
dense interacting boson systems with $k$-body interaction.

##  VI NPC and Information entropy

The NPC in wavefunction characterizes various layers of chaos in interacting
particle systems [16, 54, 55] and for a system like atomic nuclei, NPC for
transition strengths is a measure of fluctuations in transition strength sums
[37]. For an eigenstate $|E_{i}\rangle$ spread over the basis states
$|{\kappa}\rangle$, with energies
$\xi_{\kappa}=\langle{\kappa}|H|{\kappa}\rangle$, NPC (also known as inverse
participation ratio) is defined as,

$\mbox{NPC}(E)=\left\\{{\displaystyle\sum\limits_{\kappa}{\left|{C_{{\kappa}}^{i}}\right|^{4}}}\right\\}^{-1}$
(20)

NPC essentially gives the number of basis states
$\left.\left|{\kappa}\right.\right\rangle$ that constitute an eigenstate with
energy $E$. The GOE value for NPC is $d/3$. NPC can be studied by examining
the general features of the strength functions $F_{\xi}(E)$. The smooth forms
for NPC$(E)$ can be written as [53],

$\mbox{NPC}(E)=\displaystyle\frac{d}{3}\left\\{\displaystyle\int
d\xi\;\displaystyle\frac{\rho^{H_{\kappa}}(\xi)[F_{\xi}(E)]^{2}}{[\rho^{H}(E)]^{2}}\right\\}^{-1}\;,$
(21)

where $\rho^{H_{\kappa}}(\xi)$ and $\rho^{H}(E)$ are normalized eigenvalue
densities generated by diagonal Hamiltonian $H_{\kappa}$ matrix and full
Hamiltonian $H$ matrix, respectively. Taking $E$ and $\xi$ as zero centered
and scaled by corresponding widths, the above equation can be written in terms
of $f_{qN}$ and $f_{CqN}$ [37, 38],

$\mbox{NPC}(E)=\displaystyle\frac{d}{3}\left\\{\displaystyle\int_{S(q)}d\xi\;\displaystyle\frac{f_{qN}(\xi|q)[f_{CqN}(E|\xi;\zeta,q)]^{2}}{f_{qN}(E|q)}\right\\}^{-1}\;,$
(22)

In general, $q$’s in the above equation need not be same [37, 38]. However, in
the thermalization region, with $\zeta^{2}\leq 1/2$, one can approximate
$\gamma_{2}\approx(q-1)$ in Eq.(8). Then, the formula for $q$ given by Eq.(9)
is valid for $f_{qN}$ as well as for $f_{CqN}$. This is well verified
numerically in Section II. Also, the results of $\gamma_{2}$ in Figure 5(c)
corroborate this claim. With this, it is possible to simplify Eq.(22) using
Eqs.(6) and (7) and a simple two parameter formula, valid in chaotic domain,
for NPC can be written as,

$\mbox{NPC}(E)=\displaystyle\frac{d}{3}\displaystyle\left\\{\sum_{n=0}^{\infty}\frac{\zeta^{2n}}{[n]_{q}!}\,H_{n}^{2}(E|q)\right\\}^{-1},$
(23)

It is easy to see from above formula that NPC($E$) approaches GOE value $d/3$
as $\zeta\rightarrow 0$. Also for $q\rightarrow 1$, $f_{qN}$ and $f_{CqN}$ in
Eq.(22) reduce to Gaussian and then Eq.(23) gives similar results obtained for
$k=2$ in [53]. We have tested this formula with numerical ensemble averaged
BEGOE(1+$k$) results. Figure 7, shows results for ensemble averaged NPC vs.
normalized energy, for a 100 member BEGOE(1+$k$) with $m=10$ and $N=5$ example
for different values of $\lambda$ and $k$. The ensemble averaged NPC values
are shown with red solid circles and continuous lines are obtained using the
theoretical expression given by Eq. (23). One can see from the results that
with fixed $k$ (i) for small value of $\lambda$, where the one-body part of
the interaction is dominating, the numerical NPC values are zero and the
theoretical curve is far away from the numerical results indicating that the
wavefunctions are completely localized (the bottom panels in Figure 7); (ii)
with further increase in $\lambda$, the theoretical estimate for NPC in the
chaotic domain is much above the ensemble averaged curve indicating that the
chaos has not yet set in; (iii) However, with sufficiently large $\lambda$, we
see that the ensemble averaged curve is matching with the theoretical estimate
given by Eq. (23), indicating that system is in chaotic domain corresponding
to the thermalization region given by $\zeta^{2}\sim 1/2$ [23] and the
strength functions $F_{\xi}(E)$ are well represented by conditional $q$ normal
distribution. Again with further increase in $\lambda$ (the top panels in
Figure 7), the match between the theoretical chaotic domain estimate and the
ensemble averaged values is very well in the bulk part of the spectrum
($|E|<2$) for all values of $k$ with deviations near the spectrum tails.
Hence, in the chaotic domain, the energy variation of NPC($E$) using Eq. (23)
is essentially given by two parameters, $\zeta$ and $q$. The results clearly
show that the thermalization sets in faster with increase in the body rank
$k$.

Another statistical quantity normally considered is the information entropy
defined by $S^{\text{info}}(E)=-\sum_{{\kappa}}p_{\kappa}^{i}\ln
p_{\kappa}^{i}=-\sum_{\kappa}|C_{{\kappa}}^{i}|^{2}\ln|C_{{\kappa}}^{i}|^{2}$,
here $p_{\kappa}^{i}$ is the probability of basis state ${\kappa}$ in the
eigenstate at energy $E_{i}$. The localization length, $l_{H}$ is related to
$S^{\text{info}}(E)$ by $l_{H}(E)=\exp{S^{info}(E)}/(0.48d)$. Then the
corresponding embedded ensemble expression for $l_{H}$ involving $F_{\xi}(E)$,
can be written as[53],

$l_{H}(E)=-\displaystyle\int
d\xi\;\displaystyle\frac{F_{\xi}(E)\;\rho^{H_{\kappa}}(\xi)}{\rho^{H}(E)}\ln\left\\{\displaystyle\frac{F_{\xi}(E)}{\rho^{H}(E)}\right\\}\;.$
(24)

Replacing $\rho^{H_{\kappa}}(\xi)$ and $\rho^{H}(E)$ by $f_{qN}$ and
$F_{\xi}(E)$ by $f_{CqN}$, formula for $l_{H}$ valid in chaotic domain is
given by,

$l_{H}(E)=-\displaystyle\int
d\xi\;\displaystyle\frac{f_{CqN}(E|\xi;\zeta,q)f_{qN}(\xi|q)}{f_{qN}(E|q)}\;\ln\left\\{\displaystyle\frac{f_{CqN}(E|\xi;\zeta,q)}{f_{qN}(E|q)}\right\\}\;.$
(25)

Simplifying Eq.(25) for $l_{H}$ is an open problem and therefore, it is
evaluated numerically and results are compared with ensemble averaged
numerical results of BEGOE(1+$k$). Figure 8, shows results for ensemble
averaged $l_{H}$ vs. normalized energy $E$ for a 100 member BEGOE($1+k)$ with
$m=10$ bosons in $N=5$ sp states for different values of $k$. Here, we choose
$k$-body interaction strength $\lambda=1$ so that the system will be in
thermalization region. Numerical embedded ensemble results (red solid circles)
are compared with theoretical estimates (black curves) obtained using Eq.
(25). The $\zeta$ values are shown in the figure. A very good agreement
between numerical results and smooth form is obtained for all values of $k$ in
the bulk of the spectrum with small deviations near the spectrum tails. Hence,
in the chaotic domain, the energy variation of $l_{H}(E)$, with Eq. (25), is
essentially given by conditional $q$ forms for the strength functions.

Figure 7: Ensemble averaged NPC as a function of normalized energy $E$ for a
100 member BEGOE(1+$k$) with $m=10$ interacting bosons in $N=5$ sp states for
different values of $k$. Ensemble averaged BEGOE(1+$k$) results are
represented by solid circles while continuous curves correspond to the
theoretical estimates in the chaotic domain obtained using Eq. (23). The
ensemble averaged $\zeta$ and $q$ values are also given in the figure. GOE
estimate is represented by dotted line in each graph. Figure 8: Ensemble
averaged localization lengths $l_{H}$ vs. normalized energy $E$ for a 100
member BEGOE(1+$k$) with $m=10$ interacting bosons in $N=5$ sp states for
different $k$ values. Here, $\lambda=1$ is chosen for all $k$. Ensemble
averaged BEGOE(1+$k$) results (red solid circles) are compared with the smooth
forms obtained via Eq.(25) involving parameters $\zeta$ and $q$. The ensemble
averaged $\zeta$ values are given in the figure and Eq.(9) is used for $q$
values. Dotted lines in each graph represent GOE estimate.

##  VII Conclusions

In the present work, we have analyzed wavefunction structure of dense many-
body bosonic systems with $k$-body interaction by modeling the Hamiltonian of
these complex systems using BEGOE(1+$k$). We have shown that for dense boson
systems with BEGOE(1+$k$), the $q$-polynomials are used to describe the
transition from Gaussian to semi-circle in the state density as the strength
of the $k$-body interaction increases. A complete analytical description of
the correlation coefficient $\zeta$, which is related to variance of strength
functions, is obtained in terms of $N$,$m$,$k$ and $\lambda$ and it is found
to describe the embedded ensemble results very well for all the values of rank
of interaction $k$. Also, in the dense limit $\zeta\rightarrow 0$. We have
also obtained formula for $\lambda_{t}$ in terms of ($m$, $N$, $k$). Further,
it is shown that in the strong interaction domain ($\lambda>>\lambda_{t}$),
the strength functions make transition from Gaussian to semi-circle as the
rank of interaction $k$ increases in BEGOE(1+$k$) and their smooth forms can
be represented by the $q$-normal distribution function $f_{CqN}$ to describe
this crossover. Moreover, the variation of the lowest four moments of strength
functions computed numerically are in good agreement with the analytical
formulas obtained in [38]. With this, we have first utilized the interpolating
form for strength function $f_{CqN}$ to describe the fidelity decay in dense
boson systems after $k$-body random interaction quench. Secondly, using smooth
forms for $f_{qN}$ and $f_{CqN}$, we have also derived two parameter ($q$ and
$\zeta$) formula for NPC valid in thermalization region and shown that these
smooth forms describe BEGOE(1+$k$) ensemble averaged results very well.
Therefore, the results of this work, along with [33, 37, 38], establish that
the $q$-Hermite polynomials play a very significant role in analyzing many-
body quantum systems interacting via $k$-body interaction. The generic
features explored in this work are important for a complete description of
many-body quantum systems interacting via $k$-body interaction as the nuclear
interactions are now known to have some small 3-body and 4-body parts and
higher body interactions may become prominent in strongly interacting quantum
systems [7, 25, 26].

Following the work in [52], it is interesting to analyze power-law behavior of
fidelity decay for very long time using embedded ensembles with $k$-body
forces as smooth forms of strength functions can be represented by $f_{CqN}$.
Further, as smooth forms for the density of states can be represented by
$f_{qN}$, it is possible to study normal mode decomposition of the density of
states for various $k$ values using $f_{qN}$ [13, 17, 56] and thereby one can
study spectral statistics in strongly interacting quantum systems. This is for
future. It is also known that the strength functions and the entanglement
essentially capture the same information about eigenvector structure [55, 57]
and therefore it is important to study entanglement properties using embedded
ensembles with $k$-body forces.

## Acknowledgements

Thanks are due to Manan Vyas for collaboration in the initial stages of this
work and V. K. B. Kota for many useful discussions. Authors acknowledge
support from Department of Science and Technology(DST), Government of India
[Project No.: EMR/2016/001327]. NDC acknowledges support from the
International Centre for Theoretical Sciences (ICTS) during a visit for
participating in the program - Thermalization, Many body localization and
Hydrodynamics (Code: ICTS/hydrodynamics2019/11).

## References

  * [1] O. Bohigas, M.J. Giannoni, C. Schmit, Phys. Rev. Lett. 52 (1984) 1.
  * [2] V. K. B. Kota, Embedded Random Matrix Ensembles in Quantum Physics, Springer-Verlag, Heidelberg, 2014.
  * [3] A. Polkovnikov, K. Sengupta, A. Silva, Rev. Mod. Phys. 83 (2011) 863.
  * [4] L. D’Alessio, Y. Kafri, A. Polkovnikov, M. Rigol, Adv. Phys. 65 (2016) 239.
  * [5] F. Borgonovi, F.M. Izrailev, L.F. Santos, V.G. Zelevinsky, Phys. Rep. 626 (2016) 1.
  * [6] V. K. B. Kota, N. D. Chavda, Int. J. Mod. Phys. E 27 (2018) 1830001.
  * [7] J. S. Cotler et al., J. High Energy Phys. 05 (2017) 118.
  * [8] J. J. M. Verbaarschot, Quantum chromodynamics, in The Oxford Handbook of Random Matrix Theory edited by G. Akemann, J. Baik, and P. Di. Francesco, Oxford University Press, Oxford, 2011.
  * [9] V. K. B. Kota N. D. Chavda, Entropy 20 (2018) 541.
  * [10] J. B. French, S. S. M. Wong, Phys. Lett. B 33 (1970) 449.
  * [11] O. Bohigas, J. Flores, Phys. Lett. B 34 (1971) 261.
  * [12] V. K. B. Kota, Phys. Rep. 347 (2001) 223.
  * [13] T.A. Brody, J. Flores, J.B. French, P.A. Mello, A. Pandey, S.S.M. Wong, Rev. Mod. Phys. 53 (1981) 385.
  * [14] T. Papenbrock, H. A. Weidenmüller, Rev. Mod. Phys. 79 (2007) 997.
  * [15] V.V. Flambaum, G.F. Gribakin, F.M. Izrailev, Phys. Rev. E 53 (1996) 5729.
  * [16] V.V. Flambaum, F.M. Izrailev, Phys. Rev. E 56 (1997) 5144.
  * [17] K. K. Mon, J. B. French, Ann. Phys. (N.Y.) 95 (1975) 90.
  * [18] K. Patel, M.S. Desai, V. Potbhare, V.K.B. Kota, Phys. Lett. A 275 (2000) 329.
  * [19] T. Asaga, L. Benet, T. Rupp, H. A. Weidenmüller, Eur. Phys. Lett. 56 (2001) 340.
  * [20] T. Asaga, L. Benet, T. Rupp, H.A. Weidenmüller, Ann. Phys. (N.Y.) 298 (2002) 229.
  * [21] N. D. Chavda, V. Potbhare, V. K. B. Kota, Phys. Lett. A 311 (2003) 331.
  * [22] N. D. Chavda, V. Potbhare, V. K. B. Kota, Phys. Lett. A 326 (2004) 47.
  * [23] N.D. Chavda, V.K.B. Kota, V. Potbhare, Phys. Lett. A 376 (2012) 2972.
  * [24] N. D. Chavda, V.K.B. Kota, Ann. Phys. (Berlin) 529 (2017) 1600287.
  * [25] D. W. E. Blatt, B. H. J. McKellar, Phys. Rev. C 11 (1975) 2040.
  * [26] H-W. Hammer, A. Nogga, A. Schwenk, Rev. Mod. Phys. 85 (2013) 197.
  * [27] K. D. Launey, T. Dytrych, J. P. Draayer, Phys. Rev. C 85 (2012) 044003.
  * [28] A.M. Garcia-Garcia, Y.Jia, J.J.M. Verbaarschot, Phys. Rev. D 97 (2018) 106003.
  * [29] A.M. Garcia-Garcia, T. Nosaka, D. Rosa, J.J.M. Verbaarschot, Phys. Rev. D 100 (2019) 026002.
  * [30] A. Ortega, M. Vyas, L. Benet, Ann. Phys. (Berlin) 527 (2015) 748.
  * [31] A. Ortega, T. Stegmann, L. Benet, Phys. Rev. E 94 (2016) 042102.
  * [32] A. Ortega, T. Stegmann, L. Benet, Phys. Rev. E 98 (2018) 012141.
  * [33] M. Vyas, V. K. B. Kota, J. Stat. Mech. Theor. Expt. 10 (2019) 103103.
  * [34] A.M. Garcia-Garcia, J.J.M. Verbaarschot, Phys. Rev. D 94 (2016) 126010.; Phys. Rev. D 96 (2017) 066012.
  * [35] Y.Jia, J.J.M. Verbaarschot, J. High Energy Phys. 07 (2020) 193.
  * [36] L. Erdos, D. Schroder, Math. Phys. Anal. Geom. 17 (2014) 441.
  * [37] V. K. B. Kota, M. Vyas, arXiv:2003.09191v1
  * [38] V. K. B. Kota, M. Vyas, arXiv:2011.05799v1
  * [39] L. Benet, T. Rupp, H. A. Weidenmüller, Ann. Phys. (N.Y.) 292 (2001) 67.
  * [40] M. Vyas, T. H. Seligman, AIP Conf. Proc. 1950 (2018) 030009.
  * [41] M. E. H. Ismail, D. Stanton, G. Viennot, Europ. J. Combinatorics 8 (1987) 379.
  * [42] P. Szablowski, Electronic Journal of Probability 15 (2010) 1296.
  * [43] R. A. Small, S. Müller, Ann. Phys. (N.Y.) 356 (2015) 269.
  * [44] V.K.B. Kota, V. Potbhare, Phys. Rev. C 21 (1980) 2637.
  * [45] B. Lauritzen, P. F. Bortignon, R. A. Broglia, V. G. Zelevinsky, Phys. Rev. Lett. 74 (1995) 5190.
  * [46] D. Angom, S. Ghosh, V.K.B. Kota, Phys. Rev. E 70 (2004) 016209.
  * [47] S. K. Haldar, N. D. Chavda, M. Vyas, V. K. B. Kota, J. Stat. Mech. Theor. Expt. 2016 (2016) 043101.
  * [48] E. J. Torres-Herrera, M. Vyas, L. F. Santos, New J. Phys. 16 (2014) 063010.
  * [49] E.J. Torres-Herrera, J. Karp, M. Tavora, L. F. Santos, Entropy 18 (2016) 359.
  * [50] L. F. Santos, E. J. Torres-Herrera, AIP Conf. Proc. 1912 (2017) 020015.
  * [51] M.Schiulaz, E. J. Torres-Herrera, and L. F. Santos, Phys. Rev. B 99 (2019) 174313.
  * [52] M. Távora, E. J. Torres-Herrera, L. F. Santos, Phys. Rev. A 95 (2017) 013604.
  * [53] V.K.B. Kota, R. Sahu, Phys. Rev. E 64 (2001) 016219.
  * [54] P.G. Silvestrov, Phys. Rev. E 58 (1998) 5629.
  * [55] C. Mejia-Monasterio, J. Richert, T. Rupp, H.A. Weidenmüller, Phys. Rev. Lett. 81 (1998) 5189.
  * [56] R. J. Leclair, R. U. Haq, V. K. B. Kota, and N. D. Chavda, Phys. Lett. A 372 (2008) 4373.
  * [57] W. G. Brown, L. F. Santos, D. J. Starling, and L. Viola, Phys. Rev. E 77 (2008) 021106.