June 28, 2026
Embedded random matrix ensembles operating in nuclear shell model spaces, with nucleons occupying a finite set of single particle orbits and interacting via a two-body interaction, form the basis for statistical shell model. With sufficiently strong interaction, the level densities in shell model spaces take close to a Gaussian form and transition strength distributions close to a bivariate Gaussian form. In practice, partitioning via spherical configurations (\(\widetilde{m}\)) and angular momentum \(J\) (also isospin where appropriate) are essential. The resulting statistical spectroscopy or statistical shell model was applied successfully in the past in some studies of nuclear level densities, orbit occupancies, \(\beta\)-decay matrix elements and so on. Going beyond these, recently it is recognized that embedded ensembles, in a better approximation, generate in-fact \(q\)-normal form (\(q=1\) gives Gaussian and \(q=0\) Wigner’s semi-circle) for density of eigenvalues, bivariate \(q\)-normal form for transition strengths and conditional \(q\)-normal form for strength functions. These then allow us to develop statistical shell model with \(q\)-normal forms. These new developments in embedded ensembles and statistical shell model are briefly reviewed in this paper. Also described, using some examples, is the role of the \(q\) parameter in generating statistical properties of general quantum many-particle systems.
[orcid=0000-0001-5958-1143]
[orcid=0000-0003-2084-1109]
The subject of statistical shell model (SSM) (also called spectral distribution method or statistical nuclear spectroscopy) started with an article in 1967 by J.B. French where he showed that shell model energy (\(E\)) eigenvalue densities \(\rho(E)\) (i.e. frequency function for \(E\)) are close to Gaussian in form [1]. Also, the partial densities defined over shell model spherical configurations (\(\Gamma = \widetilde{m}\) or \(\widetilde{m} J\) or \(\widetilde{m} JT\) ) are Gaussians with small corrections. Note that, summing the partial densities \(\rho^{\Gamma}(E)\) will generate level densities. These results are obtained using the then developed Rochester-Oak Ridge shell model code [2]. Construction of Gaussians is possible without using shell model matrices as the centroid \(E_c\) and variance \(\sigma^2\) defining a Gaussian will follow from the traces of the shell model Hamiltonian \(H\) and that of \(H^2\) respectively. Similarly, lower order corrections to the Gaussian follow from \(H^3\) and \(H^4\) traces. Note that a Gaussian \(\rho_{{\cal G}}(E) = (1/\sqrt{2\pi}\sigma)\exp -(E-E_c)^2/2\sigma^2\). Statistical spectroscopy or SSM with Gaussians (plus corrections) is possible due to the fact that traces of operators propagate from the defining spaces to the \(m\)-particle spaces with symmetries. Various developments and applications of SSM till about 2003 are described in detail in two books [3], [4] and in a conference proceedings [5]. Some major contributions in the development of SSM in the early years are as follows. (i) As Gaussians extend from \(-\infty\) to \(+\infty\), a method to determine the ground state energy \(E_g\) is needed and Ratcliff produced a method [6] that is used by many groups in the applications of SSM. (ii) Trace propagation principles with configuration traces as an example are described in detail in [7]. (iii) Besides level densities, SSM theory for expectation values of an operator \(K\) follow from the parametric differentiation of the eigenvalue density, i.e. the density generated by \(H +\alpha K\) with \(\alpha \rightarrow 0\). This give a polynomial expansion (polynomials defined by \(\rho(E)\)) for expectation values [8]. More importantly, this shows for example, orbit occupancies (expectation values of orbit number operators) will be essentially linear in energy (polynomial expansion has simple convergence properties) and thus it is possible to calculate with good accuracy ground state orbit occupancies. Examples for occupancies for neutrinoless double beta decay nuclei calculated using SSM and comparison with data are given recently [9]. (iv) SSM also gives a double polynomial expansion for transition strength densities. Given a transition operator \({\cal O}\), the transition strength connecting a initial state with energy \(E_i\) to a final state \(E_f\) is \(\left|\left\langle E_f \mid {\cal O}\mid E_i\right\rangle\right|^2\). Then, transition strength density is the transition strength multiplied by the state densities at the two energies involved. Examination of the convergence of the double polynomial expansion showed that the transition strength densities take bivariate Gaussian form and the correlation coefficient defining this bivariate Gaussian is given by the trace of the operator \(({\cal O}^\dagger H {\cal O}H)\) [10]. Finally, basis for the various distributions mentioned above lie in the operation of random matrix ensembles in shell model spaces and in particular the ensembles that include particle number, fermion character and interactions; see [11] and also [10], [12]. We will discuss this in more detail ahead.
Investigating from 1995 some aspects of spectral distributions using the OXBASH code [13], [14], in light of the developments in quantum chaos in nuclei and other finite quantum many-particle systems, Zelevinsky’s group considered several improvements in the applications of SSM with particular reference to nuclear level densities. These studies include the following. (i) Developed exponential convergence method for much better determination of the ground state and also advocated using shell model codes, where possible, for exact determination of the ground state [15]. (ii) Showed that direct construction of \(\widetilde{m}J\) (where needed \(JT\)) partial densities is much better than using \(\widetilde{m}\) partial densities with Bethe’s spin-cutoff factors for \(J\) projection. Following this and using some of the earlier work on trace propagation, generated formulas for \(\widetilde{m}J\) centroids (trace of \(H\)) and variances (trace of \(H^2\)) that are good for machine calculations. More importantly, recognized that the partial variances should be treated properly when multi-\(\hbar \omega\) excitations are involved [16]. (iii) Developed a good method for exact removal of the center-of-mass spurious states from level densities [17]. (iv) Recognized that it is important to remove the tails of the Gaussians by using finite range Gaussians [18]. Incorporating (i)-(iv), a high-performance algorithm to calculate spin- and parity-dependent nuclear level densities is developed and made available the associated computer codes [18], [19]. These codes are used in tabulating level densities of \((2s1d)\) shell nuclei [20] and also in some applications involving \((2p1f)\) shell nuclei [21]. All the work by Zelevinsky’s group is well summarized in [22]; see also [23].
Although SSM approach is well developed for level densities, orbit occupancies and transition strengths (we will leave aside the improvements still needed in the theory for transition strengths), most important question is about the basis for SSM. Very early, it is recognized by French and Wong [24] and Bohigas and Flores [25] that random matrix theory (RMT) provides the basis but with the inclusion of the two-body nature of the nucleon-nucleon interaction. This gave rise to the introduction of two-body random matrix ensemble (TBRE) and shell model calculations with TBRE gave the ensemble averaged eigenvalue density to be Gaussian [classical Wigner-Dyson Gaussian orthogonal (GOE) and unitary (GUE) ensembles for example give semi-circle form [26]) consistent with the shell model results generated by realistic effective interactions. Note that, in a TBRE the various two-particle \(J\) matrices are represented by independent GOE’s and then the two-body interaction defined by each member is used in shell model codes to generate the \(H\) matrix for a given nucleus (for valence nucleons systems). With more general \(k\)-body interactions and for \(m\) spinless fermions occupying \(N\) number of single particle (sp) states, we have the so called fermionic embedded random matrix ensembles of \(k\)-body interactions FEE(\(k\)). With a GOE in \(k\) particle spaces, we have in \(m\) particle spaces FEGOE(\(k\)). Similarly, a GUE in \(k\) particle spaces generate FEGUE(\(k\)). Mathematical definition of FEGOE(\(k\))/FEGUE(\(k\)) is given in Section 2 and their extension to bosonic ensembles BEGOE(\(k\))/BEGUE(\(k\)) is straightforward. The FEGOE(\(k\)) were investigated analytically for the first time by Mon and French [11] and they showed, using the so called binary correlation approximation, that as \(k\) changes from 1 to \(m\) the eigenvalue density changes from Gaussian form to semi-circle form. We will return to this important result in RMT in Section 3. Though the FEGOE(\(k\)) were investigated further by Rochester group [10], [27], a major further development is in investigating FEGOE with two-body interactions in presence of a mean-field, i.e. FEGOE(1+2). Analysis of FEGOE(1+2) ensembles (then the two-particle matrix generated by a two-body interaction is represented by a GOE and the mean-field one-body part is taken to be fixed in one particle space or it is also represented by an independent GOE in one particle spaces) showed that in the many-particle spaces, as the strength of the two-body interaction increases, there is onset of chaos and finally thermalization. To this end, the change in level statistics (from Poisson to GOE Wigner form), strength functions (changing from Breit-Wigner to Gaussian form), chaos measures such as information entropy and occupancy entropy, fidelity decay etc. are studied using both shell model examples and FEGOE(1+2) systems [also BEGOE(1+2) systems]. For details of these studies see [12], [22], [28]–[31] and references therein. In the shell model studies used are Rochester-Oak Ridge [2], OXBASH [13], NATHAN [28], [32] and BIGSTICK [33] codes. The FEGOE(1+2) results clearly established that nuclear interactions are strong enough to validate SSM principles [23].
Though it is well established that embedded random matrix ensembles (EE)provide the basis for the Gaussian forms used in SSM, more recently using a closer analysis of some of the analytical results for FEGOE(\(k\))/FEGUE(\(k\)), it is recognized that EE indeed generate \(q\)-normal forms [34]–[36]. Therefore, the \(q\)-normal forms, instead of Gaussians, need to be used in SSM [37]. Mathematical details of \(q\)-normal forms are given for example in [38]–[41] and they are briefly discussed in Section 2. This new advancement in EE and SSM owes to the RMT results for the so called SYK model as derived by Verbaarschot and collaborators [42]–[45]. The SYK-RMT model involves Majorana fermions and it is a two parameter embedded random matrix ensemble. These are number of sp states \(N\) and the body rank \(k\) of the interactions. The classical ensemble GOE for example is a one parameter ensemble with the parameter being the matrix dimension \(d\). In comparison, the basic EE(\(k\)) are three parameter ensembles (see Section 2.1) with the parameters being number of sp states \(N\), number of particles (fermions or bosons) \(m\) and the body rank of the interactions \(k\). Verbaarschot et al. [42]–[44] proved analytically that the eigenvalue density for SYK-RMT is \(q\)-normal. This has opened up a new line of investigation of EE(\(k\))’s and their applications to SSM. The purpose of the present article is to present a short review of the results obtained so far for the \(q\) normal forms in SSM. Now we will give a preview.
In Section 2, for completeness, we will introduce EGOE(\(k\))/EGUE(\(k\)) for fermion and boson systems. Also presented are some of the important properties of \(q\)-normal, bivariate \(q\)-normal and conditional \(q\)-normal. In Section 3, binary correlation results for the moments of eigenvalue densities, transition strength densities and strength functions are presented and they show that these distributions follow \(q\)-normal, bivariate \(q\)-normal and conditional \(q\)-normal form respectively. Following these results, in Section 4 results from further investigations of EE(\(k\)) and EE(1+\(k\)) are presented. Section 5 deals with \(q\)-normal and \(q\)-Hermite extensions of SSM. Finally, Section 6 gives conclusions.
Given a system of \(m\) particles (fermions or bosons) distributed in \(N\) degenerate sp states and interacting via \(k\)-body \((1 \leq k \leq m)\) interactions, a class of embedded ensembles are generated by representing the \(k\) particle Hamiltonian by a classical GOE/GUE and then the many-particle Hamiltonian (\(m>k\)) is generated by propagating the \(k\) particle matrix elements to \(m\) particle spaces using the Hilbert space geometry. In other words, \(k\)-particle Hamiltonian is embedded in the \(m\)-particle Hamiltonian and the non-zero \(m\)-particle Hamiltonian matrix elements are appropriate linear combinations of the \(k\)-particle matrix elements. Therefore, these random matrix ensembles are generically called embedded random matrix ensembles or simply embedded ensembles. Due to the \(k\)-body selection rules, many matrix elements of the \(m\)-particle Hamiltonian will be zero unlike in a GOE/GUE.
The random \(k\)-body Hamiltonian operator in second quantized form for a FEGOE/BEGOE (\(\beta=1\)) and FEGUE/BEGUE (\(\beta=2\)) is, \[H(k,\beta) = \displaystyle\sum_{\alpha,\;\gamma} \; v^{\alpha,\gamma}_{k,\beta} \; \psi^\dagger(k; \alpha) \; \psi(k;\gamma) \;. \label{eq1}\tag{1}\] Here, \(\alpha\) and \(\gamma\) are indices denoting \(k\)-particle states (configurations) in occupation number basis. Distributing \(k\) particles (fermions in agreement with Pauli’s exclusion principle or bosons) in \(N\) sp states will generate the complete set of these distinct configurations. Total number of these configurations are \(\binom{N}{k}\) for fermions and \(\binom{N+k-1}{k}\) for bosons. Operators \(\psi^\dagger(k; \alpha)\) and \(\psi(k;\gamma)\) respectively are \(k\)-particle creation and annihilation operators for fermions or bosons, i.e. \(\psi^\dagger(k; \alpha) = \prod_{i=1}^{k} a^\dagger_{\mu_i}\), \(\psi(k;\gamma) = \prod_{i=1}^{k} a_{\mu_i}\) for fermions and \(\psi^\dagger(k; \alpha) = {\cal N}_{\alpha} \; \prod_{i=1}^{k} b^\dagger_{\mu_i}\), \(\psi(k;\gamma) = {\cal N}_{\gamma} \; \prod_{i=1}^{k} b_{\mu_i}\) for bosons. Here, \({\cal N}_{\alpha}\) and \({\cal N}_{\gamma}\) are the factors that guarantees unit normalization of \(k\)-particle bosonic states. The creation and annihilation operators \(a^\dagger_p\) and \(a_q\) for fermions (\(b^\dagger_p\) and \(b_q\) for bosons) satisfy the usual anti-commutation (commutation) relations for fermions (bosons).
In Eq. (1 ), \(v^{\alpha,\;\gamma}_{k,\beta}\) matrix is chosen to be a \(\binom{N}{k}\) [\(\binom{N+k-1}{k}\)] dimensional GOE (for \(\beta=1\)) or GUE (for \(\beta=2\)) in \(k\)-particle spaces for fermion [boson] systems. That means \(v^{\alpha,\;\gamma}_{k,\beta}\) are anti-symmetrized (symmetrized) \(k\)-particle matrix elements for fermions (bosons) chosen to be independent random Gaussian variables with zero mean and variance \[{\overline{v^{\alpha,\gamma}_{k,\beta} \; v^{\alpha^\prime,\gamma^\prime}_{k,\beta}}} = v^2 \; \left( {\delta_{\alpha,\gamma^\prime}} {\delta_{\alpha^\prime,\gamma}} + \delta_{\beta,1} {\delta_{\alpha,\alpha^\prime}} {\delta_{\gamma^\prime,\gamma}} \right) \;. \label{eq2}\tag{2}\] Here and else where in this article, the bar denotes ‘ensemble averaging’ and we choose \(v=1\) without loss of generality. Now, distributing the \(m\) fermions (bosons) in all possible ways in \(N\) sp states generates all many-particle basis states defining \(d_F(N,m)=\binom{N}{m}\) [\(d_B(N,m)=\binom{N+m-1}{m}\)] dimensional Hilbert space. Action of the Hamiltonian operator \(H(k,\beta)\) defined by Eq. (1 ) on the many-particle states generates the FEGOE(\(k\))/FEGUE\((k)\)/BEGOE(\(k\))/BEGUE\((k)\) ensembles in \(m\)-particle spaces. Thus, these EE will have three parameters \((N,m,k)\). However, if we include other degrees of freedom such as spin, then there will be more parameters; see [23], [30], [31] for EE with additional degrees of freedom.
Let us begin with \(q\) numbers \([n]_q\) defined by \[\left[n\right]_q = \displaystyle\frac{1-q^n}{1-q} = 1+q + q^2 + \ldots+q^{n-1}\;. \label{eq3}\tag{3}\] Note that \([n]_{q \rightarrow1}=n\). Similarly \([n]_q! = \displaystyle\Pi^{n}_{j=1} \,[j]_q\) with \([0]_q!=1\). Given these, the \(q\)-normal distribution \(f_{qN}(x|q)\), with \(x\) being a standardized variable (then \(x\) is zero centered with variance unity), is defined as [38]–[40] \[f_{qN}(x|q) = \displaystyle\frac{\displaystyle\sqrt{1-q} \displaystyle\prod_{k^\prime=0}^{\infty} \left(1- q^{k^\prime+1}\right)}{2\pi\,\displaystyle\sqrt{4-(1-q)x^2}}\; \displaystyle\prod_{k^\prime=0}^{\infty} \left[(1+q^{k^\prime})^2 - (1-q) q^{k^\prime} x^2\right]\;. \label{eq4}\tag{4}\] The \(f_{qN}(x|q)\) function is defined in the range with \[S(q) = \left(x_{min} = -\displaystyle\frac{2}{\displaystyle\sqrt{1-q}}\;,\;x_{max}=+\displaystyle\frac{2}{\displaystyle\sqrt{1-q}}\right) \label{eq5}\tag{5}\] and \(q\) takes values \(0\) to \(1\) in all the results presented in this review. Note that \(f_{qN}(x|q) = 0\) outside \(S(q)\) and the integral of \(f_{qN}(x|q)\) over \(S(q)\) is unity, \(\int_{S(q)} f_{qN}(x|q)\,dx =1\). For \(q=1\), taking the limit properly will give \(f_{qN}(x|1)= (1/\sqrt{2\pi})\,\exp-x^2/2\), the Gaussian with \(S(q=1)=(-\infty , \infty)\). Also, \(f_{qN}(x|0)=(1/2\pi) \sqrt{4-x^2}\), the semi-circle with \(S(q)=(-2,2)\). If we put back the centroid \(\epsilon\) and the width \(\sigma\) in \(f_{qN}\), then \(S(q)\) changes to \[S(q:\epsilon,\sigma) = \left(\epsilon -\displaystyle\frac{2\sigma}{\sqrt{1-q}}\;,\;\epsilon +\displaystyle\frac{2\sigma}{\sqrt{1-q}}\right)\;.\] The \(q\)-Hermite polynomials \(He_n(x|q)\), that are orthogonal with \(f_{qN}\) as the weight function, are defined by the recursion relation \[x\,He_n(x|q) = He_{n+1}(x|q) + \left[n\right]_q\,He_{n-1}(x|q) \label{eq6}\tag{6}\] with \(He_0(x|q)=1\) and \(He_{-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). The Orthogonal property of \(He_n(x|q)\)’s is \[\displaystyle\int^{2/\sqrt{1-q}}_{-2/\sqrt{1-q}} He_n(x|q)\,He_m(x|q)\,f_{qN}(x|q)\,dx = \left[n\right]_q!\,\delta_{mn}\;. \label{eq7}\tag{7}\] Using Eq. (7 ), it is easy to derive formulas for the lower order reduced central moments \(\mu_r(q)=\int_{Sq)}\,x^r f_{qN}(x)\,dx\) giving for example (\(\mu_2(q)=1\) by definition), \[\begin{array}{l} \mu_4(q) = 2+q\;,\\ \mu_6(q) = 5+6q+3q^2+q^3\;,\\ \mu_8(q) = 14+28q+28q^2+20q^3+10q^4+4q^5+q^6\;. \end{array}\label{eq8}\tag{8}\] Note that all \(\mu_{2r+1}(q)\) with \(r=1,2,\ldots\) are all zero as \(f_{qN}(x)\) is symmetrical in \(x\). See [46] for various properties of moments, central moments, reduced central moments and cumulants [lowest two cumulants are the shape parameters skewness (\(\gamma_1=\mu_3\)) and excess (\(\gamma_2=\mu_4-3\)) parameters] of probability distributions.
Going further, the bivariate \(q\)-normal distribution \(f_{biv-qN}(x,y|\xi , q)\), normalized to unity with with \(x\) and \(y\) being standardized variables and defined over \(S(q)\) in both \(x\) and \(y\) spaces, is given by [39] \[\begin{array}{l} f_{biv-qN}(x,y|\xi , q) = f_{qN}(x|q) f_{qN}(y|q) h(x,y|\xi , q)\;;\\ \\ h(x,y|\xi ,q) \\ = \displaystyle\prod_{k^\prime=0}^\infty \displaystyle\frac{1-\xi^2 q^{k^\prime}}{ (1-\xi^2 q^{2k^\prime})^2 -(1-q)\,\xi\, q^{k^\prime}\,(1+\xi^2 q^{2k^\prime})\,xy + (1-q)\xi^2 q^{2k^\prime} (x^2 +y^2)}\;, \end{array}\label{eq9}\tag{9}\] where \(\xi\) is the bivariate correlation coefficient. Note that \(f_{qN}(x|q)\) and \(f_{qN}(y|q)\) are the marginal densities of \(f_{biv-qN}\). Bivariate reduced central moments \(\mu_{rs}(q) = \int_{S(q)} x^ry^s f_{biv-qN}(x,y|\xi , q)dxdy\) are symmetrical, i.e. \(\mu_{rs}=\mu_{sr}\) and \(\mu_{rs}=0\) for \(r+s\) odd. Also \(\mu_{20}=\mu_{02}=1\) and \(\xi=\mu_{11}\). The \(f_{biv-qN}\) can be expanded in terms of \(q\)-Hermite polynomials giving [39] (called Poisson-Mehler formula), \[f_{biv-qN}(x,y|\xi , q) = f_{qN}(x|q) f_{qN}(y|q) \left[\displaystyle\sum_{n=0}^{\infty} \displaystyle\frac{\xi^n}{\left[ n\right]_q!}\,He_n(x|q) He_n(y|q)\right]\;. \label{eq10}\tag{10}\] Putting \(\xi=1\) on both sides and operating over \(S(q)\) gives, \[\delta(x-y) = f_{qN}(x|q) \left[\displaystyle\sum_{n=0}^{\infty} \displaystyle\frac{1}{\left[n\right]_q!} He_n(x|q) He_n(y|q)\right] \label{eq11}\tag{11}\] and this formula plays an important role in SSM. More importantly, Eq. (10 ) allows us to derive formulas for \(\mu_{rs}(q)\). For example, the lower order bivariate moments \(\mu_{rs}\) with \(r+s=4\) and \(6\) are, \[\begin{array}{l} \mu_{40}(q)=\mu_{04}(q)=2+q\;,\\ \mu_{31}(q)=\mu_{13}(q)= \xi\,(2+q)\;,\\ \mu_{22}(q) = 1 + \xi^2\,(1+q)\;,\\ \mu_{60}(q)=\mu_{06}(q)=5+6q+3q^2+q^3\;,\\ \mu_{51}(q)=\mu_{15}(q)=\xi\,\mu_{60}(q)\;,\\ \mu_{42}(q)=\mu_{24}(q)=(2+q) +\xi^2\,(3+5q+3q^2+q^3)\;,\\ \mu_{33}(q)=\xi\;(2+q)^2 + \xi^3\,(1+q)(1+q+q^2)\;. \end{array}\label{eq12}\tag{12}\]
Given the bivariate \(q\)-normal \(f_{biv-qN}\), the conditional \(q\)-normal densities (\(f_{CqN}\)) follow easily from Eq. (9 ). Then, with \(h(x,y|\xi ,q)\) from Eq. (9 ), we have \[\begin{array}{l} f_{CqN}(x|y; \xi , q)=f_{qN}(x|q) h(x,y|\xi ,q)\;,\\ f_{CqN}(y|x; \xi , q)=f_{qN}(y|q) h(x,y|\xi ,q)\; \end{array}\label{eq13}\tag{13}\] A very important property of \(f_{CqN}\) that follows from Eq. (10 ) is \[\displaystyle\int_{S(q)} He_n(x|q) f_{CqN}(x|y; \xi ,q) dx = \xi^n He_n(y|q)\;. \label{eq14}\tag{14}\] From Eq. (14 ), it is easy to infer that \(f_{CqN}\) is normalized to unity over \(S(q)\). It is important to recognize that the centroid \(\epsilon(y:\xi ,q)\) of \(f_{CqN}(x|y;\xi ,q)\) is not zero and its variance \(\sigma^2(y:\xi ,q)\) is not unity. Using Eq. (14 ) we have, \[\epsilon(y:\xi ,q) = \xi\,y\;,\;\;\;\sigma^2(y:\xi ,q) = (1-\xi^2) \label{eq15}\tag{15}\] Thus, the centroid is linear in \(y\) and variance independent of \(y\). In addition, by writing \(x^3\) and \(x^4\) in terms of \(q\)-Hermite polynomials and using Eq. (14 ) we have for the skewness and excess the formulas \[\begin{array}{rcl} \gamma_1(y:\xi ,q) & = & -(1-q)\displaystyle\frac{\xi y}{\displaystyle\sqrt{1-\xi^2}}\;,\\ \gamma_2(y:\xi ,q) & = & (q-1) + (1-q)^2 \left[\displaystyle\frac{\xi y}{\displaystyle\sqrt{1-\xi^2}}\right]^2 + (1-q^2)\displaystyle\frac{\xi^2}{(1-\xi^2)}\;. \end{array}\label{eq16}\tag{16}\] Thus, the \(\gamma_1\) and \(\gamma_2\) are zero only when \(q=1\) (then the conditional distribution is a Gaussian and this result is well known [46]). Let us mention that the so called Al-Salam-Chihara polynomials are orthogonal with \(f_{CqN}\) as the weight function [39].
Turning to the eigenvalue densities (also transition strength densities and strength functions or local density of states) generated by EE(\(k\)), we will discuss in this Section results for FEGOE(\(k\))/FEGUE(\(k\)) and point out their extensions to BEGOE(\(k\))/BEGUE(\(k\)) where appropriate. It is important to recognize that the various densities considered in this Section all will be of same form for FEGOE(\(k\)) and FEGUE(\(k\)). Therefore we will be using either of them though the final results apply to both.
The ensemble averaged eigenvalue density for FEGOE(\(k\))/FEGUE(\(k\)) is deduced using the moment method supplemented by numerical verifications in [34]. Firstly, the moments are defined by \(\overline{\left\langle H^p \right\rangle^m}\) with \(H\) given in Eq. (1 ). By definition of the ensembles, the odd moments (p odd) will be zero and this includes the centroid. Formulas for the even moments to order 8 (\(p=8\)) are derived in [11] for FEGOE(\(k\)) and in [47] for FEGUE(\(k\)). Used here is the dilute limit defined by \(N \rightarrow \infty\), \(m \rightarrow \infty\), \(m/N \rightarrow 0\) and \(k/m\) fixed. Then, using the so called binary correlation approximation [11], [30], the reduced moments up to order 8 for FEGUE(\(k\)) are (with the variance or the second moment is given by \(\overline{\left\langle H^2 \right\rangle^m} = \binom{m}{k} \,\overline{\left\langle H^2(k) \right\rangle^k}\)) [47], \[\begin{array}{l} \mu_4(m,k) = 2 + f(m,k,1) \;,\\ \mu_6(m,k) = 5 + 6f(m,k,1) + 3\left[f(m,k,1)\right]^2 + f(m,k,2)f(m,k,1) \;,\\ \mu_8(m,k) = 14 + 28f(m,k,1) + 28\left[f(m,k,1)\right]^2 + 12\left[f(m,k,1)\right]^3 + 8f(m,k,2)f(m,k,1)\;, \\ + 4f(m,k,1) \left[f(m,k,2)\right]^2 + 8 \left[f(m,k,1)\right]^2 f(m,k,2)\\ + f(m,k,1)f(m,k,2)f(m,k,3) + 2\left[f(m,k,1)\right]^2 \Delta\;; \\ f(m,k,r)={\displaystyle\binom{m}{k}}^{-1}\;\displaystyle\binom{m-rk}{k}\;. \end{array}\label{eq17}\tag{17}\] For the formula for \(\Delta\) in Eq. (17 ), see [47]. Comparing these with the FEGOE(\(k\)) formulas given in [11], it is seen that the moments to order 6 for FEGOE(\(k\)) are same as those given in Eq. (17 ) and for \(\mu_8\) only the last term is different. Comparing the formulas in Eqs.(17 ) and (8 ), it is seen that the lower order reduced moments of FEGUE(\(k)\) and FEGOE(\(k\)) will be essentially same those of the \(q\)-normal \(f_{qN}(x)\) if we identify \(q=f(m,k,1)\). The differences in the 6th and 8th moments are verified to be a few percent or less. Significantly, as seen from the formula for \(\mu_4(m,k)\), the \(q\)-parameter for FEGOE(\(k\))/FEGUE(\(k\)) is given by \[q = \mu_4 -2 = {\displaystyle\binom{m}{k}}^{-1}\;\displaystyle\binom{m-k}{k} \sim 1-k^2/m \;. \label{eq18}\tag{18}\] Including finite \(N\) corrections to \(\mu_4\) as given in [48]–[50], a better approximation for the \(q\) parameter for FEGOE(\(k\))/FEGUE(\(k\)) is, \[\begin{array}{l} q(N,m,k) = \displaystyle\binom{N}{m}^{-1} \displaystyle\sum_{\nu=0}^{min(k,m-k)}\; \displaystyle\frac{\Lambda^\nu(N,m,m-k)\;\Lambda^\nu(N,m,k)\;d(g_\nu)}{ \left[\Lambda^0(N,m,k)\right]^2} \,; \\ \Lambda^\nu(N,m,r) = \displaystyle\binom{m-\nu}{r}\;\displaystyle\binom{N-m+r-\nu}{r}\;,\;\;\;d(g_\nu) = \displaystyle\binom{N}{\nu}^2-\displaystyle\binom{N}{\nu-1}^2\;. \end{array}\label{eq19}\tag{19}\] Note that \(\overline{\left\langle H^2\right\rangle^m} = \Lambda^0(N,m,k)\). Numerical calculations show that \(q\) rapidly becomes \(q=0\) as \(k\) value increases. Note that \(q=0\) for GOE/GUE and then correctly we have the semi-circle form for the eigenvalue density. Similarly, for \(q=1\) we have Gaussian form as expected from shell model with large \(m\) value and using realistic interactions, TBRE and also EGOE(2) (then \(k/m << 1\) and \(q \sim 1\)). For example, for \(N=20\) and \(m=8\), Eq. (19 ) gives \(q=0.814\), \(0.417\), \(0.119\), 0.015 and \(0\) for \(k=1\), \(2\), \(3\), \(4\) and \(\ge 6\) respectively. In addition, numerical embedded ensemble results for Gaussian to semi-circle transition are shown in Fig. 1 and they are compared with the curve given by \(f_{qN}(x)\) with \(q\) defined by Eq. (19 ). The agreement is excellent.
Turning to boson systems, here dense limit defined by \(N \rightarrow \infty\), \(m \rightarrow \infty\), \(m/N \rightarrow \infty\) and \(k/m\) fixed is important (dense limit will not exist for fermionic systems). Extension of the FEGOE(\(k\))/FEGUE(\(k\)) results for the moments given by Eq. (17 ) to BEGOE(\(k\))/BEGUE(\(k\)) is available only for \(\mu_4\) and it is proved using the so called \(N \rightarrow -N\) law [50]. It is expected that the \(q\)-normal also applies to BEGOE(\(k\))/BEGUE(\(k\)) in the dense limit and this is well verified by numerical results in Fig. 2. It is important to recognize that the finite \(N\) formula for \(\mu_4\) gives the formulas for \(q\) for BEGOE(\(k\))/BEGUE(\(k\)) [34], \[\begin{array}{l} q(N,m,k) = \displaystyle\binom{N+m-1}{m}^{-1} \displaystyle\sum_{\nu=0}^{\nu_{max}}\; \displaystyle\frac{ \Lambda_B^\nu(N,m,m-k)\;\Lambda_B^\nu(N,m,k) \;d_B(g_\nu)}{\left[\Lambda_B^0(N,m,k)\right]^2}\;;\\ \Lambda_B^\nu(N,m,r) = \displaystyle\binom{m-\nu}{r}\;\displaystyle\binom{N+m+\nu-1}{r}\;,\;\;\; d_B(g_\nu) = \displaystyle\binom{N+\nu-1}{\nu}^2-\displaystyle\binom{N+\nu-2}{\nu-1}^2\;. \end{array}\label{eq20}\tag{20}\] In Figs. 1 and 2, the values for \(q\) given by finite \(N\) formulas are used in constructing \(f_{qN}(x)\). Eq. (20 ) gives for example for a \((N=5, m=10)\) system the \(q\) values to be 0.969, 0.861, 0.664, 0.405, 0.172, 0.045, 0.008, 0 for \(k=1\), 2, 3, 4, 5, 6, 7, \(\geq 8\).
Although we presented the results only for identical fermion (also identical boson) systems, for systems with two types of fermions (for example protons and neutrons) with H preserving the two fermion numbers, it is seen that the \(q\)-normal form extend to the eigenvalue density for two-fermion systems; to establish this, formulas for moments up to sixth order are derived in Ref. [51]. The \(q\)-normal is seen to extend also to two species boson systems. The extension of the \(q\)-normal form proton-neutron systems is clearly important for applications to SSM.
Let us consider a system of \(m\) spinless fermions in \(N\) sp states with the Hamiltonian \(H\) a \(k\)-body operator. Say \(H\) generates the eigenstates \(\left|E\right\rangle\) with \(E\) denoting energy. Now, given a \(t\)-body transition operator \({\cal O}(t)\) acting on an eigenstate \(\left|E_i\right\rangle\) will populate the state \(\left|E_f\right\rangle\) with transition strength given by \(\left|\left\langle E_f \mid {\cal O}\mid E_i\right\rangle\right|^2\). The resulting bivariate transition strength density (normalized to unity) is, \[\rho_{biv-{\cal O}}(E_i,E_f) = \left[\left\langle\left\langle{\cal O}^\dagger{\cal O}\right\rangle\right\rangle^m\right]^{-1}\; \left\langle\left\langle{\cal O}^\dagger\delta(H-E_f) {\cal O}\delta(H-E_i)\right\rangle\right\rangle^m\;. \label{eq21}\tag{21}\] Note that \(\left\langle\left\langle X \right\rangle\right\rangle^m =\sum_E \left\langle E \mid X \mid E\right\rangle\). Our interest is in deriving the statistical law for the transition strength densities as appropriate for nuclei. To this end, \(H\) is represented by FEGOE(\(k\))/FEGUE(\(k\)) and the \({\cal O}\) by an independent FEGOE(\(t\))/FEGUE(\(t\)). Then, formulas for the ensemble averaged (average with respect to both \(H\) and \({\cal O}\) ensembles) bivariate moments \(\mu_{rs}(m,k,t)\) of \(\overline{\rho_{biv-{\cal O}}(E_i,E_f)}\) are derived in [35] for \(r+s \leq 6\) using again the binary correlation approximation. These formulas are given in [10], [12], [52] and quite strikingly, they are close to those given by Eq. (12 ) for \(f_{biv-qN}\). We will describe this briefly below. Let us add that dividing the \(\overline{\rho_{biv-{\cal O}}(E_i,E_f)}\) by the eigenvalue densities at \(E_i\) and \(E_f\) with give the ensemble averaged transition strengths. Also, the bivariate moments are defined by \(\overline{\left\langle{\cal O}^\dagger H^s {\cal O}H^r\right\rangle^m}\). Note that all odd moments (\(r+s\) odd) will be zero. Thus, the centroids are zero and the variances \(\mu_{20} = \mu_{02} = \binom{m}{k}\). Similarly \(\left\langle{\cal O}^\dagger{\cal O}\right\rangle^m = \binom{m}{t}\). Now we will consider all other bivariate moments with \(r+s=2\), 4 and 6.
Firstly, the bivariate correlation coefficient \(\xi=\mu_{11}\) for ensemble averaged transition strength density is \[\mu_{11} = \xi= \left[\binom{m}{k}\right]^{-1}\;\displaystyle\binom{m-t}{k} \;. \label{eq22}\tag{22}\] Similarly, the moments to order 4 and 6 and the \(q\) parameter are (by using equations in [10], [12] and rewriting them in terms of \(\xi\) and \(q\) parameters), \[\begin{array}{rcl} \mu_{40} = \mu_{04} & = & 2 + q\;;\;\;\; q = \left[\binom{m}{k}\right]^{-1}\,\binom{m-k}{k}\;,\\ \mu_{31} = \mu_{13} & = & \xi \mu_{40}\;,\\ \mu_{22} & = & 1 + \xi^2 (1+q)+ \xi (\Delta_0)\;,\\ \mu_{60} = \mu_{06} & = & 5 + 6q + 3 q^2 + q^3 + q (\Delta_1)\;,\\ \mu_{51} = \mu_{15} & = & \xi\,\mu_{60}\;,\\ \mu_{42} = \mu_{24} & = & (2+q) + \xi^2 (3 + 5 q + 3 q^2 + q^3) + \xi (X)\;,\\ \mu_{33} & = & \xi \left[4 + 4q + q^2\right] + \xi^3 \left[ 1 + 2q +2q^2 + q^3\right] + \xi (Z)\;. \end{array}\label{eq23}\tag{23}\] Formulas for \(\Delta_0\), \(\Delta_1\), \(X\) and \(Z\) are given in [10], [12]. Numerical calculations are used to verify that for some typical values of \((N,m,k,t)\), these are indeed \(\sim 0\). With this, by comparing Eq. (23 ) with Eq. (12 ), it is clear that the transition strength densities generated by FEGOE/FEGUE are well represented by the bivariate \(q\)-normal distribution. Thus, by changing \(E_i\) and \(E_f\) by the corresponding standardized variables \(x\) and \(y\) respectively, \[\overline{\rho_{biv-{\cal O}}(x,y)} \rightarrow f_{biv-qN}(x,y)\] with \(\xi\) and \(q\) given by Eqs. (19 ) and (22 ). Let us mention that it is well known in statistics [46] and in random matrix theory [4], [27] that lower order moments generate the form of a probability distribution.
Eq. (19 ) gives the formula for the \(q\) parameter with finite \(N\) corrections. Similarly, formula with finite \(N\) corrections to the correlation coefficient \(\xi\) is [35] \[\xi(t,k) = \displaystyle\sum_{\nu=0}^{min(t,m-k)}\; \displaystyle\frac{\Lambda^\nu(N,m,m-t)\;\Lambda^\nu(N,m,k)\;d(g_\nu)}{\binom{N}{m}\; \Lambda^0(N,m,k)\;\Lambda^0(N,m,t)} \,. \label{eq24}\tag{24}\] Note that the function \(\Lambda\) is defined in Eq. (19 ). For example, for \((N,m)=(20,10)\) the values for \(\xi(1,k)\) are 0.682, 0.559, 0.455, 0.364, 0.152 for \(k=2\), 3, 4, 5 and 8 respectively. Similarly, for \(\xi(2,k)\) they are 0.465, 0.314, 0.21, 0.136 and 0.027 respectively.
Although we have restricted to transition operators that preserve \(m\) (then \(E_i\) and \(E_f\) belong to the same system) in the discussion above, it is also possible to analyze \(\mu_{rs}\) with \(r+s=4\) and \((rs)=(11)\) for beta and neutrinoless double beta decay type operators and also for particle removal operators using the results in [52]. More importantly, they will also give formulas, with finite \(N\) corrections, for \(\xi\) and \(q\) for the transition strength densities generated by these operators. Again, these are important for applications of bivariate \(q\)-normal for transition strength densities in SSM.
Wavefunction structure in finite quantum many-body systems such as atomic nuclei follows from the form of the strength functions. Given the eigenstates expanded in terms of a set of physically motivated basis states, strength functions correspond to the spread of a basis state over the eigenstates. They also correspond to partial densities in SSM; for a more general discussion of strength functions see [12]. By deriving analytical formulas, using EE, for the lowest four moments of the strength functions, it is shown in [36] that the conditional \(q\)-normal\(f_{CqN}\) to a good approximation represents strength functions. For good numerical tests of this result, see [53], [54]. Here below we will describe the EE formulas and important structures they display.
Let us consider a system of \(m\) fermions in \(N\) single particle (sp) states and the Hamiltonian \(H\) for the system is say, \[H = H_0(t) + \lambda V(k) \label{eq25}\tag{25}\] where \(H_0\) is a \(t\)-body operator, \(V\) is a \(k\)-body operator and \(\lambda\) is the strength parameter. We will assume that \(t << k\) and for \(m\) fermions, obviously \(k \leq m\). In many physical applications \(t=1\) with \(H_0\) representing a mean-field one-body Hamiltonian [30]. Now, strength functions form is studied by considering the structure of eigenfunctions of \(H\) expanded in terms of the unperturbed \(H_0\) eigenstates (basis states). Denoting \(\left|\kappa, \alpha\right\rangle\) as the eigenstates of \(H_0\) forming a complete set with \(H_0 \left|\kappa, \alpha\right\rangle= E_\kappa \left|\kappa, \alpha\right\rangle\) and \(\left|E, \beta\right\rangle\) as the eigenstates of \(H\) forming a complete set with \(H\left|E, \beta \right\rangle= E \left|E, \beta \right\rangle\) (with \(\alpha\) and \(\beta\) labeling the respective degeneracies in \(H_0\) and \(H\) spectrum), we can expand the eigenstates of \(H_0\) in the eigenbasis of \(H\) giving, \[\left|\kappa, \alpha\right\rangle= \sum_{E, \beta} C_{\kappa, \alpha}^{E, \beta} \left|E, \beta \right\rangle\;. \label{eq26}\tag{26}\] Now, the strength function \(F_\kappa(E)\) is, \[\begin{array}{rcl} F_\kappa(E) & = & \left\langle\delta(H-E)\right\rangle^{\kappa}\;;\\ & = & \displaystyle\frac{1}{d\rho_1(E_\kappa))} \displaystyle\sum_{\alpha \in \kappa , \beta \in E} \left|C_{\kappa , \alpha}^{E , \beta}\right|^2 \;. \end{array}\label{eq27}\tag{27}\] Here, \(d\) is \(m\) particle space dimension, \(d\rho_1(E_\kappa)\) gives number of \(H_0\) states with same basis state energy \(E_\kappa\) and similarly \(d\rho_2(E)\) gives number of eigenstates of \(H\) with same eigen energy \(E\). Note that \(\rho_2(E) = \left\langle\delta(H-E)\right\rangle^m\) is the eigenvalue density generated by \(H\) and similarly, \(\rho_1(E_\kappa) = \left\langle\delta(H_0 - E_\kappa) \right\rangle^m\) is the eigenvalue density generated by \(H_0\).
In order to derive the form of \(F_\kappa(E)\), formulas for the ensemble averaged lower order moments of \(F_\kappa(E)\) are derived by representing \(H_0(t)\) by FEGOE(\(t\))/FEGUE(\(t\)) and \(V(k)\) by FEGOE(\(k\))/FEGUE(\(k\)) and assume that they are independent. From now on, for brevity, often we will drop \(t\) in \(H_0\) and \(k\) in \(V\). By definition, we have \(\overline{\left\langle H_0\right\rangle^m}=0\) and \(\overline{\left\langle H\right\rangle^m}=0\). These are the centroids of \(\overline{\rho_1(E_\kappa)}\) and \(\overline{\rho_2(E)}\) respectively. Similarly, the corresponding variances are \(\sigma_{H_0}^2 = \overline{\left\langle H_0^2\right\rangle^m}\) and \(\sigma_{H}^2 = \overline{\left\langle H^2\right\rangle^m} =\overline{\left\langle H_0^2\right\rangle^m} + \lambda^2 \overline{\left\langle V^2\right\rangle^m}\). Scaling the eigenvalues \(E\) with their width \(\sigma_H\), the moments of \(F_\kappa(E)\) are given by \[M_r(E_\kappa) = \left[\sigma_H^r\right]^{-1}\;\overline{\left\langle H^r \right\rangle^\kappa}\;. \label{eq28}\tag{28}\] It is important to note that \(\left\langle H_0^p \right\rangle^{\kappa} = E^p_\kappa\) as \(\kappa\) are eigenstates of \(H_0\) with eigenvalues \(E_\kappa\). Therefore, \(\left\langle H^r\right\rangle^\kappa = \left\langle H^r\right\rangle^{E_\kappa}\) is a expectation value and it can be written in terms of the polynomials generated by \(\overline{\rho_1(E_\kappa)}\). For a general operator \(K\), the polynomial expansion to second order (usually this is sufficient) is[4], [8], [36], \[\begin{array}{l} \left\langle K \right\rangle^{\kappa} = \displaystyle\sum_\mu \left\langle K P_\mu(H_0) \right\rangle^m P_\mu(E_\kappa)\;;\\ P_0(x) = 1,\;\;\;P_1(x) = \hat{x},\;\;\;P_2(x) = \displaystyle\frac{(\hat{x})^2-1}{\displaystyle\sqrt{\mu_4-1}}\;. \end{array}\label{eq28a}\tag{29}\] Note that \(\hat{x}=E_{\kappa}/\sigma_{H_0}\) and \(\mu_4\) is the fourth reduced moment of \(\overline{\rho_1(E_\kappa)}\). Before proceeding further, let us mention that the moments \(M_r\) in Eq. (28 ) give the central moments \({\cal M}_r\) to be (for \(r \le 4\)), \[\begin{array}{rcl} {\cal M}_2 & = & M_2 - M_1^2\;,\\ {\cal M}_3 & = & M_3 - 3M_2 M_1 +2M_1^3 \;,\\ {\cal M}_4 & = & M_4 - 4M_3 M_1 +6M_2 M_1^2 -3M_1^4\;. \end{array}\label{eq29}\tag{30}\] Without going into details (see [36] for the details) we will now give the formulas for \(M_r\), \(r=1-4\). These will then give the variance \(\sigma^2(E_\kappa)\), skewness \(\gamma_1(E_\kappa)\) and excess \(\gamma_2(E_\kappa)\) via Eq. (30 ). Firstly the centroid \(M_1(E_\kappa)\) is \[M_1(E_\kappa) = \left[\sigma_H\right]^{-1} \;\overline{\left\langle H \right\rangle^{\kappa}} = \left[\sigma_H\right]^{-1}\;\overline{\left\langle H_0\right\rangle^{\kappa}} = \xi\, {\hat{E}_\kappa}\;. \label{eq30}\tag{31}\] Note that the ensemble average of \(\left\langle(H_0)^r [V(k)]^s \right\rangle^\kappa\) is zero for any \(r\) and any odd \(s\) and this is used above (for \(r=0,s=1\)). More importantly, \(\xi\) is a correlation coefficient and clearly (with \({\hat{E}_\kappa}=E_\kappa/\sigma_{H_0}\)), \[\xi=\sigma_{H_0}/\sigma_H = \displaystyle\sqrt{\displaystyle\frac{\binom{m}{t}}{\binom{m}{t} + \lambda^2 \; {\binom{N}{t}}^{-1}\,\binom{N}{k}\;\binom{m}{k}}}\;. \label{eq31}\tag{32}\] Now, the second moment \(M_2(E_\kappa)\) is (absorbing \(\lambda\) in \(V(k)\) and using \(\sigma_H^2=\sigma_{H_0}^2+\sigma_V^2\)), \[\begin{array}{rcl} M_2(E_\kappa) & = & \left[\sigma^2_H\right]^{-1} \overline{\left\langle H^2 \right\rangle^\kappa} = \left[\sigma^2_H\right]^{-1}\;\overline{\left\langle H^2_0 + V^2 + (H_0V + VH_0) \right\rangle^{\kappa}} = \xi^2 ({\hat{E}_\kappa})^2 + (1-\xi^2) \\ \Rightarrow {\cal M}_2(E_\kappa) & = & (1-\xi^2) \;. \end{array}\label{eq31a}\tag{33}\] Thus, the strength function centroid is linear in \(E_\kappa\) (with slope \(\xi\)) and the variance \({\cal M}_2(E_\kappa)\) is independent of \(E_\kappa\). Now let us consider \(M_3(E_\kappa)\), \[\begin{array}{l} M_3(E_\kappa) = \left[\sigma^3_H\right]^{-1} \;\overline{\left\langle H^3 \right\rangle^{\kappa}} \\ = \left[\sigma^3_H\right]^{-1}\;\overline{\left\langle H_0^3 + V^3 + (H_0 V^2 + V^2H_0) + (H^2_0 V + VH^2_0) + VH_0V + H_0 V H_0 \right\rangle^{\kappa}} \\ = \xi^3 ({\hat{E}_\kappa})^3 + 2 \xi (1-\xi^2){\hat{E}_\kappa}+ \left[\sigma^3_H\right]^{-1}\; \overline{\left\langle V H_0 V\right\rangle^{\kappa}}\;. \end{array}\label{eq32}\tag{34}\] Now evaluating the last term using only the first two terms in Eq. (29 ), and simplifying Eqs. (34 ) and (30 ) will give the remarkable formula [36], \[\mu_3({\hat{E}_\kappa}) = \gamma_1({\hat{E}_\kappa}) = - \displaystyle\frac{\xi \left(1-q^{hv}\right) {\hat{E}_\kappa}}{\displaystyle\sqrt{1-\xi^2}} \;. \label{eq33}\tag{35}\] Here, the \(q^{hv}\) is \[q^{hv} = \displaystyle\frac{\overline{\left\langle H_0 V H_0 V\right\rangle^m}}{\sigma_{H_0}^2 \sigma_V^2} = {\binom{m}{k}}^{-1}\;\binom{m-t}{k}\;. \label{eq34}\tag{36}\] Therefore, the first three moments \(M_1(E_\kappa)\), \({\cal M}_2(E_\kappa)\) and \(\mu_3(E_\kappa)\) of the strength function from EE are same as those of \(f_{CqN}(y)\) with \(y={\hat{E}_\kappa}\) and \(\xi\) and \(q\) given by Eqs. (32 ) and (36 ) respectively; see Eqs. (15 ) and (16 ). Finally, the fourth moment is given by \[M_4(E_\kappa) = \left[\sigma^4_H\right]^{-1} \;\overline{\left\langle H^4 \right\rangle^{\kappa}} = \left[\sigma^4_H\right]^{-1} \;\overline{\left\langle\left(H^2_0 + V^2 + (H_0V + VH_0)\right)^2 \right\rangle^{\kappa}}\;. \label{eq35}\tag{37}\] Using Eq. (29 ) and carrying out simplifications as described in detail in [36], we have \[\mu_4({\hat{E}_\kappa}) \simeq \mu_4^0({\hat{E}_\kappa}) = \left(2+q^{hv}\right) + \displaystyle\frac{\xi^2 ({\hat{E}_\kappa})^2 \left(1-q^{hv}\right)^2 + \xi^2 \left[1-\left(q^{hv}\right)^2\right]}{1-\xi^2}\;. \label{eq36}\tag{38}\] with \(\gamma_2({\hat{E}_\kappa}) = \mu_4({\hat{E}_\kappa})-3\), it is seen by comparing Eq. (38 ) with Eq. (16 ) that the strength functions from EE are well approximated by \(f_{CqN}(E|y;\xi, q^{hv})\) with \(y={\hat{E}_\kappa}\) and \(\xi\) and \(q^{hv}\) given by Eqs. (32 ) and (36 ) respectively. It is useful to note that the finite \(N\) formula for \(q^{hv}\) is [36], \[q^{hv} = \displaystyle\frac{\displaystyle\sum_{\nu=0}^{min(t,m-k)}\, \Lambda^{\nu}(N,m,k)\, \Lambda^{\nu}(N,m,m-t)\,d(g_\nu)}{ \displaystyle\binom{N}{m}\;\Lambda^0(N,m,t)\,\Lambda^0(N,m,k)}\;. \label{eq37}\tag{39}\] Also, the \(q^V\) from \(V(k)\) is given by Eq. (19 ) and the same equation gives \(q^h\) from \(H_0(t)\) by replacing \(k\) by \(t\). Now, some comments are in order. (i) Firstly let us comment on the role of \(\lambda\) in Eq. (25 ). It is well known that for very small values of \(\lambda\), the strength functions will be delta functions at \({\hat{E}_\kappa}\) and as \(\lambda\) increases they will take Breit-Wigner (BW) form. With further increase the BW form changes to \(f_{CqN}\) with \(\lambda\) sufficiently large. Thus, for the strength functions to take \(f_{Cqn}\) form \(\lambda\) need to be large and this form will be good for \(\xi^2=1/2\) (this is thermalization regime) and then \(\sigma^2_{H_0} = \lambda^2 \sigma^2_{V}\); see [23], [30], [31], [36], [53], [54]. An open question is to construct a \(q\)-BW form. (ii) It is seen from Eq. (35 ) that \(\gamma_1({\hat{E}_\kappa})\) is not zero and therefore \(F_\kappa(E)\) from EE are not symmetrical. For \({\hat{E}_\kappa}\) negative \(F_\kappa(E)\) will be skewed in the positive direction and for \({\hat{E}_\kappa}\) positive it will be skewed in the negative direction . (iii) From Eq. (38 ), it is easy to see that \(\gamma_2({\hat{E}_\kappa})=q^{hv}(1-q^{hv})\) for \({\hat{E}_\kappa}=0\) and \(\xi^2=1/2\) and hence in the thermodynamic region it is always positive. (iv) Numerical calculations show that for EE the \(f_{CqN}\) form is certainly good for \(|{\hat{E}_\kappa}| \le 2\) and it may need corrections beyond this (see Section 4.3 for further discussion); (v) Using \(f_{CqN}\) form it is possible to derive EE formulas for the quantum chaos measures information entropy and inverse participation ratio [36], [53], [54]. (vi) Though we do not have BEE formulas for the lower order moments \(M_r({\hat{E}_\kappa})\), it is expected that the \(f_{CqN}\) form applies to BEE in the dense limit (see [53] and Section 4.3).
Ensemble averaged (smoothed) eigenvalue density, transition strength density and strength functions that take \(q\)-normal forms as described in Section 3 are directly useful in SSM (see Section 5). However, the errors for example in using the smooth eigenvalue density can be estimated only if we have the knowledge of two-point and higher order correlations in eigenvalues (topic of correlations in eigenvector components or transition strengths is a far more complex subject [27]). More importantly,level and strength fluctuations as quantified by the two and higher point functions are essential for understanding quantum chaos and thermalization in isolated finite many-particle quantum systems. Dyson defined correlation functions in eigenvalues [55] and the lowest of these is the two-point correlation function (smoothed eigenvalue density is the one-point function). For a random matrix, it is given by the ensemble average of the product of the density of eigenvalues at two eigenvalues say \(E\) and \(E^\prime\). Fluctuation measures such as the number variance and Dyson-Mehta \(\Delta_3\) statistic are defined by the two-point function and so also the variance of the level motion in the ensemble. As follows from the results in Section 3.1, clearly the eigenvalue density generated by a member of FEGOE(\(k\))/FEGUE(\(k\)) can be expanded by starting with the (smoothed) \(q\)-normal form and using the associated \(q\)-Hermite polynomials \(He_\zeta(x|q)\). Covariances \(\overline{S_\zeta S_{\zeta^\prime}}\) (overline representing ensemble average) of the expansion coefficients \(S_\zeta\) with \(\zeta \ge 1\) here determine the two-point function as they are a linear combination of the bivariate moments \(\Sigma_{PQ}\) of the two-point function. Recently [56], using binary correlation approximation formulas for the covariances \(\overline{S_\zeta S_{\zeta^\prime}}\) for low values of \(\zeta\) are derived and they are briefly described in this Section. Let us stress that for the GOE/GUE not just the lowest covariances, but the full two-point function was obtained using the binary correlation approximation in [27], [57]. Using this method or any other, till now the two-point function could not be obtained for FEGOE(\(k\))/FEGUE(\(k\)) with \(k \ne m\) (let us mention that \(k=1\) is special [58]–[60]). Here below, we will restrict to FEGUE(\(k\)) and the results can be extended to FEGOE(\(k\)).
Given the ensemble averaged eigenvalue density \(\overline{\rho(E)}\) of FEGUE(\(k\)) where \(\rho(E)\) is the eigenvalue density for each member of EGUE(\(k\)), integral of \(\rho(E)\) defines the distribution function, \(F(x) = d\,\int_{-\infty}^x \rho(E)\,dE\). Note that \(F(x)\) gives number of levels up to the eigenvalue \(x\) and \(d=\binom{N}{m}\) is the matrix dimension. Now, the two-point correlation function \(S^{\rho}(x,y)\) for the eigenvalues and its integral version \(S^F(x,y)\) are (in this Section \(x\) and \(y\) also denote energies or eigenvalues), \[\begin{array}{rcl} S^{\rho}(x,y) & = & \overline{\rho(x)\;\rho(y)} -{\overline{\rho(x)}}\;{\overline{\rho(y)}}\;,\\ S^F(x,y) & = & d^2\; \displaystyle\int_{-\infty}^x \displaystyle\int_{-\infty}^{y} S^\rho(x^\prime,y^\prime) dx^\prime dy^\prime\;\;=\;\;\overline{F(x)\;F(y)} -{\overline{F(x)}}\;{\overline{F(y)}}\;. \end{array}\label{fl1}\tag{40}\] It is clear that \(S^\rho\) (and \(S^F\)) gives measures for level fluctuations and the simplest two-point measure is the number variance \(\Sigma^2(\overline{n})\). Say, there are \(n\) number of levels between energies \(x\) and \(y\). Then \(n=F(x)-F(y)\) and \(\overline{n}= \overline{F(x)} - \overline{F(y)}\). With these, the number variance \(\Sigma^2(\overline{n}) = \overline{(n-\overline{n})^2}\) is simply, \[\Sigma^2(\overline{n}) = S^F(x,x) + S^F(y,y) -2 S^F(x,y)\;. \label{fl2}\tag{41}\] In addition, the Dyson-Mehta \(\Delta_3\) statistic is related to \(\Sigma^2(\overline{n})\) involving an integral with \(\Sigma^2(r)\) [27]. Also, \(S^F(x,x)\) gives the variance of the fluctuation in a eigenvalue \(E\) measured in units of the local level spacing. Importantly, \(S^F(x,y)\) and \(S^\rho(x,y)\) can be probed or constructed using the bivariate moments \({\widetilde{\Sigma}}_{PQ}\) of \(S^{\rho}(x,y)\), \[{\widetilde{\Sigma}}_{PQ} = \displaystyle\int x^P\, y^Q \,S^{\rho}(x,y) \,dx dy = \overline{\left\langle H^P\right\rangle\left\langle H^Q\right\rangle} -\overline{\left\langle H^P\right\rangle}\;\;\overline{\left\langle H^Q\right\rangle}\;. \label{fl3}\tag{42}\] Also, with \(\Sigma_{PQ} = \overline{\left\langle H^P\right\rangle\left\langle H^Q\right\rangle}\), we have \(\Sigma_{P,0} = \overline{\left\langle H^p\right\rangle^m}\), the \(P\)-th moment of \(\overline{\rho(E)}\).
Proceeding further, first the eigenvalue density \(\rho(E)\) for various members of FEGUE(\(k\)) can be expanded in terms of \(q\)-Hermite polynomials starting with \(q\)-normal giving, \[\rho(E)\,dE = f_{qN}({\hat{E}}|q) \left[1+ \displaystyle\sum_{\zeta \ge 1}^{\infty} S_\zeta\;\frac{He_\zeta({\hat{E}}|q)}{\left[\zeta\right]_q!} \right]\,d{\hat{E}}\;; \;\;{\hat{E}}=(E-E_c)/\sigma\;. \label{fl4}\tag{43}\] Here, \(S_\zeta\) are the expansion coefficients and the \(S_\zeta\) should not be confused with \(S^\rho(x,y)\) used for the two-point function. It is important to recall that the ensemble averaged eigenvalue density \(\overline{\rho(E)}\) for FEGUE(\(k\)) is \(f_{qN}\), the \(q\)-normal. The \(S_\zeta\)’s in Eq. (43 ) are for a given member of the FEGUE(\(k\)) ensemble and it is easy to see that \(\overline{S_\zeta}=0\). Note that Eq. (43 ) gives an expansion for \(S^\rho(x,y)\) in terms of \(q\)-Hermite polynomials (in the reminder of this paper, the symbols \(x\) and \(y\) are standardized variables), \[S^{\rho}(x,y) = f_{qN}(x|q)\,f_{qN}(y|q)\,\displaystyle\sum_{\zeta\,,\,\zeta^\prime= 1}^{\infty} \overline{S_\zeta\,S_{\zeta^\prime}}\;\frac{He_\zeta(x |q)}{\left[\zeta\right]_q!}\;\frac{He_{\zeta^\prime}(y |q)}{\left[\zeta^\prime\right]_q!}\;. \label{fl5}\tag{44}\] Also we have easily, \[\left\langle H^p\right\rangle= \overline{\left\langle H^p\right\rangle} + \displaystyle\sum_{\zeta \ge 1} S_\zeta\,\displaystyle\frac{\sigma^p}{\left[\zeta\right]_q!}\;\displaystyle\int_{S(q)} \,x^p\,f_{qN}(x|q)\,He_\zeta(x|q)\,dx \;. \label{fl6}\tag{45}\] Note that \(\sigma^2=\Sigma_{2,0}=\Sigma_{0,2}\). Now, writing \(x^p\) in terms of \(q\)-Hermite polynomials and using the results in [41] (see also Section 2) we have the important result, \[\begin{array}{l} \left\langle H^p\right\rangle= \overline{\left\langle H^p\right\rangle} + \displaystyle\sum_{\zeta \ge 1} S_\zeta\,\sigma^p\,C_{\frac{p-\zeta}{2},p}(q)\;;\\ \Rightarrow {\hat{\Sigma}}_{PQ} = \displaystyle\frac{{\widetilde{\Sigma}}_{PQ}}{\left[\Sigma_{2,0}\right]^{(P+Q)/2}} = \displaystyle\sum_{\zeta , \zeta^\prime= 1}^{\infty} \overline{S_\zeta\,S_{\zeta^\prime}}\;C_{\frac{P-\zeta}{2},P}(q) C_{\frac{Q-\zeta^\prime}{2},Q}(q)\;. \end{array}\label{fl7}\tag{46}\] The \(\widetilde{\Sigma}_{PQ}\) is defined by Eq. (42 ) and a formula for the \(C_{--}\) factors in Eq. (46 ) is given in [41]. Let us add that \({\hat{\Sigma}}_{PQ}=0\) for \(P+Q\) odd and similarly \(\overline{S_\zeta\,S_{\zeta^\prime}}=0\) for \(\zeta +\zeta^\prime\) is odd. Also, \({\hat{\Sigma}}_{P0}=0\), \({\hat{\Sigma}}_{PQ} = {\hat{\Sigma}}_{QP}\), \(\overline{S_\zeta}=0\) and \(\overline{S_\zeta\,S_{\zeta^\prime}}=\overline{S_{\zeta^\prime}\,S_{\zeta}}\). Using Eq. (46 ) successively with \(P+Q\) increasing from 2, the covariances \(\overline{S_\zeta\,S_{\zeta^\prime}}\) can be written in terms of the moments \({\hat{\Sigma}}_{PQ}\). For example, formulas for \(\zeta+\zeta^\prime\le 6\) are, \[\begin{array}{rcl} \overline{S_1\,S_1} & = & {\hat{\Sigma}}_{11}\;,\;\;\; \overline{S_3\,S_1} = {\hat{\Sigma}}_{31} - C_{13}\,{\hat{\Sigma}}_{11}\;,\;\;\; \overline{S_2\,S_2} = {\hat{\Sigma}}_{22}\;,\\ \overline{S_5\,S_1} & = & {\hat{\Sigma}}_{51} - C_{15}\,\overline{S_3\,S_1} - C_{25}\,\overline{S_1\,S_1}\;,\;\;\; \overline{S_4\,S_2} = {\hat{\Sigma}}_{42} - C_{14}\,\overline{S_2\,S_2}\;,\\ \overline{S_3\,S_3} & = & {\hat{\Sigma}}_{33} - C^2_{13}\,\overline{S_1\,S_1} - 2 C_{13}\,\overline{S_1\,S_3}\;. \end{array}\label{fl8}\tag{47}\] Here, \(C_{1,3}=q+2\), \(C_{1,4}=q^2+2q+3\), \(C_{1,5}=q^3+2q^2+3q+4\) and \(C_{2,5}(q) = q^3 + 3q^2 + 6q + 5\).
Formulas for the moments \(\Sigma_{PQ}\) and hence for \(\hat{\Sigma}_{PQ}\), for a system of \(m\) fermions in \(N\) single particle states are derived, for \(P+Q \le 8\) in [56] using binary correlation approximation. For example, \[\begin{array}{rcl} \Sigma_{1,1} & = & \overline{\left\langle H\right\rangle^m \left\langle H\right\rangle^m} = \displaystyle\frac{1}{d^2} \displaystyle\sum_{\alpha_1, \alpha_2} \overline{H_{\alpha_1 \alpha_1}H_{\alpha_2 \alpha_2}} \\ \Sigma_{2,2} & = & \overline{\left\langle H^2\right\rangle^m \left\langle H^2\right\rangle^m} \\ &= & \left[\overline{\left\langle H^2\right\rangle^m}\right]^2 + 2 \displaystyle\frac{1}{d^2} \displaystyle\sum_{\alpha_1, \alpha_2,\alpha_a,\alpha_b} \overline{H_{\alpha_1 \alpha_2} H_{\alpha_a \alpha_b}}\;\overline{H_{\alpha_2 \alpha_1} H_{\alpha_b \alpha_a}} \;. \end{array}\label{fl9}\tag{48}\] Note that here, \(H_{\alpha \beta}\) are \(H\) matrix elements in \(m\) particle spaces. The expressions in Eq. (48 ) are simplified by applying the Wigner-Racah algebra of \(U(N)\) as described in [50], [52]. Then, we have the finite \(N\) formulas, \[\begin{array}{rcl} \hat{\Sigma}_{1,1} & = & \displaystyle\frac{\Lambda^0(N,m,m-k)}{\binom{N}{m}\,\Lambda^0(N,m,k)}\;, \\ \hat{\Sigma}_{2,2} & = & \displaystyle\frac{2\;\displaystyle\sum_{\nu=0}^k\;\left[\Lambda^\nu(N,m,m-k)\right]^2 d(\nu)}{\left[\binom{N}{m}\,\Lambda^0(N,m,k) \right]^2}\;. \end{array}\label{fl10}\tag{49}\] Proceeding further, in [56] formulas are obtained for \(\hat{\Sigma}_{P,Q}\) with \(P+Q \le 8\). Then, in the asymptotic limit defined by \(N \rightarrow \infty\), \(m \rightarrow \infty\), \(m/N \rightarrow 0\) with \(k\) finite, we have for example for \(P+Q \le 6\) \[\begin{array}{rcl} \hat{\Sigma}_{1,1} & = & \displaystyle\frac{\binom{m}{k}}{\binom{N}{k}^2}\;,\;\;\; \hat{\Sigma}_{3,1} = 3 \displaystyle\frac{\binom{m}{k}}{\binom{N}{k}^2} = 3\,\hat{\Sigma}_{1,1} \;,\;\;\;\hat{\Sigma}_{2,2} = \displaystyle\frac{2}{\binom{N}{k}^2}\;,\\ \hat{\Sigma}_{5,1} & = & 5 \displaystyle\frac{\binom{m}{k}}{\binom{N}{k}^2} \left[ 2 + \,\displaystyle\frac{\binom{m-k}{k}}{\binom{m}{k}}\right] = (10+5q)\hat{\Sigma}_{1,1} \;,\\ \hat{\Sigma}_{4,2} & = & \displaystyle\frac{4}{\binom{N}{k}^2} \left[2 + \displaystyle\frac{\binom{m-k}{k}}{\binom{m}{k}}\right] = (4+2q)\hat{\Sigma}_{2,2} \;,\\ \hat{\Sigma}_{3,3} & = & 9 \hat{\Sigma}_{1,1} + \displaystyle\frac{3}{\binom{m}{k}\,\binom{N}{k}^2} + O\left(\frac{1}{\binom{N}{k}^4}\right) \;. \end{array}\label{fl11}\tag{50}\] Using these in Eq. (47 ) show that for \(q=1\), \(\overline{S_\zeta\,S_{\zeta^\prime}} = \delta_{\zeta \zeta^\prime} \overline{S_\zeta^2}\). This result was obtained in [11] as the assumption there is \(q=1\) for FEGOE(\(k\)) with \(k << m\). However, we see from the formula for \(q\) given before that for any reasonable values of \((N,m)\), the \(q\) value will not be close to \(1\). In addition, we see that the formulas in Eqs. (47 ) and (50 ) (also, those for \(P+Q=8\) given in [56]) agree with the formulas for GUE (i.e. for \(k=m\) or \(q=0\)) given in [27], [57]. Also, they agree with the results for FEGUE(\(k\)) in the \(k^2/m \rightarrow 0\) limit as given in [11], [27] and this corresponds to \(q \rightarrow 1\) . Thus, the formulas in Eq. (50 ) cover the two extreme limits and therefore expected to apply to all \(k\) values. Let us add that direct derivation of asymptotic limit formulas for many other \(\hat{\Sigma}_{PQ}\) for \(P+Q > 8\) may prove to be useful as they will provide systematics for \(\hat{\Sigma}_{PQ}\) and hence for \(\overline{S_iS_j}\). With this, it may be possible to carry out the sum in Eq. (44 ) and obtain the two-point function (or the number variance) for FEGUE(\(k\)) just as it was carried out using the moment method for GOE and GUE in the past [27], [30], [57].
Besides studying level fluctuations as described in the previous subsection, it is of interest to investigate how the distribution of the largest or smallest eigenvalue for FEGOE(\(k\)) / BEGOE(\(k\)) (also for the GUE versions) change as \(q\) varies from \(1\) to \(0\) (i.e. as \(k\) changes from 1 to \(m\)). This belongs to the subject of extreme value statistics (EVS) [61]. It is well known that the classical EVS are classified into Frechet, Gumbel and Weibull distributions [62], [63]. However, for the Wigner-Dyson GOE/GUE ensembles, EVS of the (lowest) largest eigenvalues is described by the celebrated Tracy-Widom (TW) distribution [64]–[66]. This corresponds to the situation with \(k=m\) for FEGOE(\(k\))/ BEGOE(\(k\))/ FEGUE(\(k\))/ BEGUE(\(k\)). Following these, there is a first attempt in [31] to study numerically the lowest eigenvalue distribution (LED) for two-body fermionic and bosonic EGOE and they are found to follow the modified-Gumbel distribution that was used earlier in a Sherrington-Kirkpatrick model study in [67]. As analytical study of EVS for FEGOE(\(k\))/BEGOE(\(k\)) and FEGUE(\(k\))/BEGUE(\(k\)) appears to be intractable, numerical analysis was attempted more recently for FEGOE(\(k\))/FEGUE(\(k\)) in [68] and for BEGOE(\(k\))/BEGUE(\(k\)) in [69]. For the fermionic ensembles, it is seen in the numerical calculations that the LED varies from Gaussian to TW form as \(k\) changes from 1 to \(m\). As attention was paid to the important result that the ensemble averaged eigenvalue density takes \(q\)-normal form for EE in the study of LED using bosonic ensembles in [69], we will describe the results of the bosonic ensembles analysis briefly. Here, in the numerical studies Gaussian, TW and modified-Gumbel forms are employed.
Firstly, the \(q\)-normal distribution for the eigenvalue density shows that the distribution has a cut-off at \(-2 [\beta \Lambda_B^{0}(N,m,k)]^{1/2} / [1-q(N,m,k)]^{1/2}\) at the lower edge when measured with respect to the ensemble averaged eigenvalue centroid; this follows easily from Eq. 5 . For finite \((N,m)\), there will be departures as \(q\)-normal is an asymptotic form for the eigenvalue densities. Therefore, the following parametrization was suggested in [69] for the centroid \(\lambda_c\) of the lowest eigenvalues \(\lambda\) for the various members of a FEGOE(\(k\))/FEGUE(\(k\)). For a given \((N,m,k)\), the ansatz is, \[\lambda_c(N,m,k) = \displaystyle\frac{-2}{\displaystyle\sqrt{1- q(N,m,k)}}{\left[ \beta \Lambda_B^{0}(N,m,k) \right]}^{\alpha}. \label{eq-g1}\tag{51}\] Note that for both BEGOE(\(k\)) and BEGUE(\(k\)), formulas for \(q(N,m,k)\) and for \(\Lambda^0_B(N,m,k)\) follow from Eq. 20 . Also, as mentioned in Section 3 the ensemble averaged variance of the eigenvalue density is \(\Lambda^0_B(N,m,k)\). In the limit \(q\)-normal form is exact, the parameter \(\alpha = 1/2\). Secondly, by applying Eq. 20 , it is easy to see that Eq. 51 gives for \(k=m\) the well-known result for TW for GOE/GUE. Therefore, using the calculated \(\lambda\) values, via least-square procedure, values of \(\alpha\) for various \((N,m,k)\) values are obtained and the results are shown in Fig. 3. For \(0 \leq q \leq 0.75\), the agreement with \(q\)-normal form result that \(\alpha = 0.5\) is almost exact. Though not shown in the figure, for \(q(N,m,k) \mathrel{\mathchoice {\vcenter{\offinterlineskip\halign{\hfil \displaystyle##\hfil\cr>\cr\sim\cr}}} {\vcenter{\offinterlineskip\halign{\hfil\textstyle##\hfil\cr>\cr\sim\cr}}} {\vcenter{\offinterlineskip\halign{ \hfil\scriptstyle##\hfil\cr>\cr\sim\cr}}} {\vcenter{\offinterlineskip\halign{\hfil\scriptscriptstyle## \hfil\cr>\cr\sim\cr}}}}0.8\) the deviations from \(\alpha = 1/2\) are significant and they correspond to \(k = 1\). As the member to member fluctuations in the spectral width are largest for \(k=1\), corrections to the \(k=1\) results are applied as described in [69] and the then the results obtained are shown in the figure. Thus, Eq. (51 ) is a good formula for LED centroid.
The variance for TW distribution is \(D^{-1/6}\) for a \(D\)-dimensional GOE/GUE. This corresponds to BEGOE(\(k\))/BEGUE(\(k\)) with \(k = m\) and \(D = \binom{N+k-1}{k}\). This is also proportional to the spectral width. Therefore, the following parametrization is suggested in [69] for BEGOE(\(k\))/BEGUE(\(k\)), \[\sigma_{\lambda} (N,m,k) = \left[\Lambda^{0}_B(N,m,k) \right]^{\mu_1},\;\;\; \sigma_{\lambda}(N,m,k) = {\left[\Lambda^{0}_B(N,m,k) \right]}^{\mu_2} \binom{N+k-1}{k}^{-1/2}. \label{eq-g2}\tag{52}\] For \(k=m\) case, i.e. for classical Gaussian ensembles or for TW, the exponents \(\mu_1=-1/6\) and \(\mu_2=1/3\) and they give correctly \(\sigma_{\lambda}=D^{-1/6}\). Using the calculated \(\lambda\) values for the same set of \((N,m,k)\) used in Fig. 3, via least-square procedure, the values of \(\mu_1\) are obtained as a function of \(q\) parameter and the results are shown in Fig. 4. For small values of \(q \mathrel{\mathchoice {\vcenter{\offinterlineskip\halign{\hfil \displaystyle##\hfil\cr<\cr\sim\cr}}} {\vcenter{\offinterlineskip\halign{\hfil\textstyle##\hfil\cr<\cr\sim\cr}}} {\vcenter{\offinterlineskip\halign{ \hfil\scriptstyle##\hfil\cr<\cr\sim\cr}}} {\vcenter{\offinterlineskip\halign{\hfil\scriptscriptstyle## \hfil\cr<\cr\sim\cr}}}}0.1\), there is a sharp increase in the values of \(\mu_1\) and then, it essentially become a constant except for \(k=1\) where the values of \(\mu_1\) decrease. Note that \(\mu_1\) becomes positive which implies larger variance compared to that for TW. This trend is common for both FEGOE/FEGUE [68]. Though not shown in the figure, the results for \(\mu_2\) vs \(q\) are essentially similar to those for \(\mu_1\) vs \(q\). Going beyond the centroid and variance of LED, there are no ansatz formulas for the shape parameters, skewness \(S\) and kurtosis \(\kappa\). Some numerical results for these are given in [68], [69]. It is seen that in general, the skewness and kurtosis values are larger than those for TW distribution for intermediate \(q\) values. For example, for BEGOE(\(k\)) for \(0.1 \mathrel{\mathchoice {\vcenter{\offinterlineskip\halign{\hfil \displaystyle##\hfil\cr<\cr\sim\cr}}} {\vcenter{\offinterlineskip\halign{\hfil\textstyle##\hfil\cr<\cr\sim\cr}}} {\vcenter{\offinterlineskip\halign{ \hfil\scriptstyle##\hfil\cr<\cr\sim\cr}}} {\vcenter{\offinterlineskip\halign{\hfil\scriptscriptstyle## \hfil\cr<\cr\sim\cr}}}}q \mathrel{\mathchoice {\vcenter{\offinterlineskip\halign{\hfil \displaystyle##\hfil\cr<\cr\sim\cr}}} {\vcenter{\offinterlineskip\halign{\hfil\textstyle##\hfil\cr<\cr\sim\cr}}} {\vcenter{\offinterlineskip\halign{ \hfil\scriptstyle##\hfil\cr<\cr\sim\cr}}} {\vcenter{\offinterlineskip\halign{\hfil\scriptscriptstyle## \hfil\cr<\cr\sim\cr}}}}0.8\) the \(S/S_{TW}\) decreases from \(\sim 4\) to \(2\) while \(\kappa/\kappa_{TW}\) from \(\sim 1.5\) to \(1\) (note that \(S_{TW} = -0.2935\) and \(\kappa_{TW} = 3.1652\)).
Going beyond the lowest four moments, LED’s obtained numerically are compared with Gaussian (\({\cal G}\)), classical TW and modified Gumbel distributions. Of these three, Gaussian is simplest one and for lowest eigenvalue \(E\) with zero center and unit variance, we have \({\cal G}(E) = (2\pi)^{-1/2}\;\exp (-E^2/2)\). The TW form is given by an integral [64]–[66] and its numerical values are used to obtain a smooth form. The modified Gumbel distribution is given by [62], [63], \[G_{\mu}(E)= w \exp \left[\mu\left(\frac{E-u}{v}\right)-\mu\exp \left( \frac{E-u}{v}\right)\right] \label{gumbel}\tag{53}\] where \(u\) and \(v\) are rescaling parameters and \(w\) is a normalization constant. The results for a 5000 member BEGOE(\(k\)) are given in Fig. 5. These are obtained for \(N = 5\) and \(m = 10\) system with \(k\) varying from 1 to 10. The numerical histograms are computed using the scaled lowest eigenvalues \(\tilde{\lambda} = [\sigma_\lambda(N,m,k)]^{-1}[\lambda - \lambda_c(N,m,k)]\). It is clearly seen from the figure that for \(k=1\), the distributions are close to Gaussian form. However, for \(2 \leq k \leq 6\) the distributions are close to modified Gumbel form and there is clear transition to TW for \(k = 10\). Thus, the LED for BEGOE(\(k\)) changes from Gaussian to modified Gumbel to TW as \(q(N,m,k)\) changes from \(\sim 1\) (\(k = 1\)) to 0 (\(k = m\)). As the results in Figs. 3, 4 and 5 are obtained for modest values of \((N,m)\), these can be further substantiated (more importantly these may lead to analytical formulation) by using larger \((N,m)\) values say with \(N\) going to \(10\) and \(m\) upto \(20\). Secondly, the results reported for fermionic systems [68] need to be re-analyzed using the \(q\) parameter so that the results for fermionic and bosonic systems can be properly compared. Also, calculations for fermionic ensembles with much larger values of \((N,m)\) are needed. All these require much larger computing power and these studies are for future.
Nuclear Hamiltonians consist of a mean-field one-body part \(h(1)\), a residual two-body part \(V(2)\) and a small \(3\)-body part (perhaps also a four-body part) [70]. Thus, \(H=h+V\) with \(h=h(1)\) and \(V=V(2)\) or \(V(2)+V(3)\) or \(V(2)+V(3)+V(4)\). Then, FEGOE(1+2) or FEGOE(1+2+3) or FEGOE(1+2+3+4) need to be analyzed. The FEGOE(1+2) has been analyzed in the past in detail, assuming that the eigenvalue densities are close to Gaussian form, in terms of the additional parameter (relative strength of the two-body part) in the model [30]. Following this, it is of interest to study the more general FEGOE(1+\(k\))[also BEGOE(1+\(k\))], i.e. with \(H\) consisting of the mean-field one-body part and a \(k\)-body interaction, with \(k\) changing from \(2\) to \(m\) and employing \(q\)-normal forms. This will be same as considering \(H\) in Eq. (25 ) with \(t=1\). Here below we will present some results for FEGOE(1+\(k\)) and BEGOE(1+\(k\)). We will return to FEGOE(1+2+3) and FEGOE(1+2+3+4) in Section 5.
Wavefunction structure for complex systems follow by examining strength functions and quantum chaos measures such as number of principle components (NPC) and information entropy. Firstly, for interacting many-fermion systems (extension of all the results in this subsection to dense boson systems is straight forward) using Hamiltonian \(H\), which is a sum of one-body \(h(1)\) and an embedded GOE of \(k\)-body interactions \(V(k)\) with strength \(\lambda\), \[H_{\lambda}(1+k) = h(1) + \lambda\,V(k), \label{eqk1}\tag{54}\] it is easy to see that the eigenvalue density will be \(q\)-normal for any \(\lambda\) value (exceptions may happen for \(h(1)\) generating singular spectra [4]). This result is well verified in a number of numerical examples in [37], [53], [54]. Going further, with \(\lambda=0\), strength functions will be delta functions at the \(h(1)\) basis state energies \(E_\kappa\) with \(\kappa\) denoting the \(m\)-particle \(h(1)\) basis states. Without loss of generality we consider \(h(1)\) defined by a set of single particle energies \(\epsilon_i\) with \(i\) denoting sp states. Then, \(h(1)=\sum_i \epsilon_i n_i\) where \(n_i\) is number operator for the \(i\)-th sp state and say there are \(N\) number of sp states. With \(\lambda=0\), the \(m\) fermion system basis states (\(\kappa\))are given by the distribution of \(m\) fermions in \(N\) sp states. The basis state energies \(E_\kappa=\sum_i \epsilon_i\) with the summation over the occupied orbits for the given \(\kappa\) configuration (in practice, i.e. in numerical calculations, the diagonal energies \(\left\langle\kappa \mid V(k) \mid \kappa\right\rangle\) are added to \(E_\kappa\)). Now, increasing \(\lambda\) value, the delta functions will spread and mix and thus, changing the form of the strength functions \(F_{\kappa}(E)\) where \(E\) are \(H\) eigenvalues in \(m\)-particle spaces. For FEGOE(1+2), it is well established that the delta functions change to BW form after some value of \(\lambda\) and with further sufficient increase in \(\lambda\), the the BW form starts changing to near Gaussian form [23], [30]. Then, with further increase in \(\lambda\) to a value \(\lambda_t\) we reach the region of thermalization where different quantities defining the eigenstate properties such as entropy, strength functions, temperature etc, give the same values irrespective of the defining basis. The value of \(\lambda_t\) is determined by the correlation coefficient \(\xi\) in Eq. (32 ). Then \(\xi^2(\lambda_t) = \sigma^2_{h(1)}/\sigma^2_{h(1) + \lambda\,V(2)}=1/2\) [71]. Assuming that this result extends to FEGOE(1+\(k\)), it follows from the results in Section 3.3 that the strength functions will take \(f_{CqN}\) form when \(\sigma^2_{h(1)}/\sigma^2_{h(1) + \lambda_t\,V(k)} \sim 1/2\) giving \(\lambda_t=\sigma_{h(1)}/\sigma_{V(k)}\). Given the \(\epsilon_i\), it is easy write formulas for moments generated by \(h(1)\) [4], [7] and then it is easy to see that \(\sigma^2_{h(1)} = [N(N-1)]^{-1} m(N-m) \sum_i {\tilde{\epsilon}}_i^2\). Note that \(\tilde{\epsilon}_i = \epsilon_i - \overline{\epsilon}\); \(\overline{\epsilon} = N^{-1} \sum_i \epsilon_i\). With this and using \(\sigma^2_{V(k)} = \Lambda^0(N,m,k)\) as given in Section 3.1, we have \[\lambda_t \sim \displaystyle\sqrt{\displaystyle\frac{m(N-m) \displaystyle\sum_i \tilde{\epsilon}_i^2}{N(N-1)\;\Lambda^0(N,m,k)}}\;. \label{eqk2}\tag{55}\] Thus, ensemble averaged strength functions for \(H_{\lambda}(1+k)\) in many fermion spaces follow \(f_{CqN}\) form for \(\lambda \sim \lambda_t\), i.e. for sufficiently large \(\lambda\) values.
Figure 6 shows the ensemble averaged strength function results obtained for FEGOE(\(1+k\)) with \(m=6\) fermions in \(N=12\) sp states. Here, \(\lambda\) is chosen equal to 0.5, so that the system is in thermalization region for all \(k\). In the ensemble calculations [54], the \(E_i\) and \(E_{\kappa}\) spectra are scaled to have zero centroid and unit width for each member of the ensemble using the eigenvalue distribution and \(E_\kappa\)-energies distribution, respectively. Therefore, \(E \rightarrow \hat{E}\) and \(E_{\kappa}/\sigma_H=\xi\,\hat{E_{\kappa}}\). The ensemble averaged \(F_{\kappa}(E)\) results are shown for \(\hat{E_\kappa} = 0.0, \pm 1.0\), and \(\pm 2.0\) using body rank \(k = 2, 3, 4\) and \(6\). All these strength function histograms \(F_\kappa(E)\) are fitted with \(f_{CqN}(\hat{E}|{\hat{E}_\kappa};\xi,q)\). For each \(k\), the smooth curves in Figure 6 are obtained using ensemble averaged values for \(\xi\) and \(q\). It is clearly seen from the results shown in the figure that the ensemble averaged histograms are in very good agreement with the smooth forms obtained using \(f_{CqN}\). From the results shown in the Figure 6, it is clearly seen that for sufficiently large \(\lambda\), the smoothed strength functions \(F_\kappa(E)\) are very well represented by \(f_{CqN}\) and they make a transition from a near Gaussian form to semi-circle form, as the body rank \(k\) in EGOE(1+\(k\)) changes from 2 to \(m\). Also, \(F_\kappa(E)\) results for \(\hat{E_\kappa}=0\) are symmetric and for \(\hat{E_\kappa} \neq 0\), away from the center of the spectrum, \(F_\kappa(E)\) results are asymmetrical about \(\hat{E}\) as expected from Eq. (35 ).
All the results described above apply also to dense boson systems [53]. This is demonstrated in Figure 7. In this figure, histograms represent ensemble averaged \(F_\kappa(E)\) results for a 250 member BEGOE(1+\(k\)) with \(m=10\) bosons in \(N=5\) sp states and \(\lambda=0.5\). The strength function plots are obtained for \({\hat{E}_\kappa}= 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 (for all \(k\)) exists in the region of thermalization [72]. The histograms, representing BEGOE(1+\(k\)) results of strength functions, are compared with the conditional \(q\)-normal density function as given by, \(F_\kappa(E)= f_{CqN}(x=E|y=E_\kappa;\xi,q)\). The smooth black curves in Figure 7 for each \(k\) are obtained via \(f_{CqN}\) using corresponding ensemble averaged \(\xi\) and \(q\) values. The results in Figure 7 clearly show very good agreement between the numerical histograms and continuous black curves for all body rank \(k\). The \(F_\kappa(E)\) results for \({\hat{E}_\kappa}=0\) are given in Figure 7 clearly demonstrate that the strength functions are symmetric and also exhibit a transition from Gaussian form to semi-circle as \(k\) changes from \(m=2\) to \(m=10\). The smooth form given by \(f_{CqN}\) interpolates this transition very well. Going further, \(F_\kappa(E)\) results for \({\hat{E}_\kappa}\neq 0\) are also shown in Figures 7. One can see that \(F_\kappa(E)\) results are asymmetrical about \(E\) as recognized earlier for bosonic systems in [73]. Also, \(F_\kappa (E)\) are skewed more in the positive direction for \({\hat{E}_\kappa}>0\) and skewed more in the negative direction for \({\hat{E}_\kappa}< 0\) exactly as predicted by Eq. (35 ).
Going beyond strength functions, we will briefly consider the quantum chaos measure NPC(\(E\)) that gives number of basis states (\(h(1)\) states) that make up an eigenstate with energy \(E\). In terms of the \(C\) coefficients in Eq. (26 ), NPC is given by \[NPC(E) = \left\{\displaystyle\sum_\kappa \left|C^E_\kappa\right|^4\right\}^{-1}\;. \label{eqk3}\tag{56}\] Using the formulation developed for FEGOE(1+2) for deriving a formula for NPC in [74] and assuming that this extends to FEGOE(1+\(k\)), a formula in terms of an integral involving strength functions and state densities can be written. Then, translating Eq.(4) in [74] to \(q\)-normal forms and assuming further that all the densities involved have the same \(q\) value, we have \[NPC(E) = \displaystyle\frac{d}{3}\;\left\{\displaystyle\int_{S(q)} d{\hat{E}_\kappa}\,f_{qN}({\hat{E}_\kappa}) \left[f_{qN}({\hat{E}})\right]^{-2}\;\left[f_{CqN}({\hat{E}}|{\hat{E}_\kappa}; \xi,q)\right]^2\right\}^{-1}\;. \label{eqk4}\tag{57}\] Now, writing \(f_{CqN}\) in terms of \(q\)-Hermite polynomials using Eqs. (13 ), (9 ) and (10 ) in that order and then carrying out the \({\hat{E}_\kappa}\) integration using Eq. (7 ) will give the formula, \[NPC(E) = \displaystyle\frac{d}{3}\;\left[h({\hat{E}}, {\hat{E}}|\xi^2 , q)\right]^{-1} \label{eqk5}\tag{58}\] with the \(h\) function given by Eq. (9 ). Eq. (58 ) was given first in [53] with a summation over \(q\)-Hermite polynomials as in Eq. (10 ). For a GOE \(H\), we have \(k=m\) with \(\xi=0\). Then, Eq. (58 ) gives correctly NPC to be \(d/3\) independent of \(E\). With \(f_{CqN}({\hat{E}}| {\hat{E}}\; \xi^2 , q)\) for \(q=1\) reducing to conditional normal form [39] and similarly \(f_{qN}\) reducing to normal form, simplifying \[h({\hat{E}}, {\hat{E}}|\xi^2 , q=1) = \left[f_{qN}({\hat{E}})\right]^{-1}\;f_{CqN}({\hat{E}}| {\hat{E}}; \xi^2 , q=1)\] gives the formula \[NPC(E:q=1) = \displaystyle\frac{d}{3}\;\displaystyle\sqrt{1-\xi^4}\;\exp-\frac{\xi^2}{1+\xi^2} \,{\hat{E}}^2\;. \label{eqk6}\tag{59}\] This is same as the formula derived in [74] where \(q=1\) is assumed. Thus, Eq. (58 ) gives correctly the known results for GOE and the formula in \(q \rightarrow 1\) limit. In addition, it is well verified in a number of numerical ensemble calculations [53], [54] that for FEGOE(1+\(k\)) and BEGOE(1+\(k\)), Eq. (58 ) applies for all \(k\) for sufficiently large values of the \(\lambda\) parameter (\(\lambda > \lambda_t\)) appearing in Eq. (54 ). As an example we show in Fig. 8 some results for FEGOE(1+\(k\)).
Nuclear Hamiltonians, as already mentioned in the previous Section, consist of a mean-field one-body part \(h(1)\), a residual two-body part \(V(2)\) and a small \(3\)-body part (perhaps also a four-body part). Just as SSM was developed and applied using Gaussian forms (see [1-10,15-22] in the past (assuming \(H\) to be one plus two-body), with the \(q\)-normal forms now established to approximate better the state, transition strength and strength functions/partial densities, it is possible to develop SSM with \(q\)-normal forms. It is important to recognize that with \(q\)-normal forms, it is possible to consider \(H\) to be not only one plus two-body but also one plus two plus three-body (may be plus four-body also). Following the first attempt in this direction in Ref. [37], in this Section described are some basic approaches one may adopt using \(q\)-normal distributions and the associated \(q\)-Hermite polynomials in SSM. Note that they will include information about the fourth moment of the density of eigenvalues and the boundedness of the eigenvalue density and other distributions in a natural way in SSM.
Let us begin with nuclear level densities that are by definition statistical quantities and they are important as they are measurable and needed for Astrophysical reaction rates calculations. Given a nucleus with valence protons (say \(m_p\) in number) occupying shell model sp orbits \(j^p_1\), \(j^p_2\), \(\ldots\), and similarly \(m_n\) number of valence neutrons occupying sp orbits \(j^n_1\), \(j^n_2\), \(\ldots\), the \((m_p,m_n)\) space can be decomposed into \(p-n\) configurations \((\widetilde{m_p}, \widetilde{m_n})\) by distributing the nucleons in their respective valence sp orbits. Then the state (eigenvalue) density \(I^{(m_p, m_n)}(E)\) can be written as a sum of the partial densities defined over \((\widetilde{m_p}, \widetilde{m_n})\) configurations giving, \[I^{(m_p, m_n)}(E) = d(m_p,m_n)\,\rho^{(m_p, m_n)}(E) = \displaystyle\sum_{(\widetilde{m_p}, \widetilde{m_n})}\, d(\widetilde{m_p},\widetilde{m_n})\,\rho^{(\widetilde{m_p}, \widetilde{m_n})}(E) = \displaystyle\sum_{(\widetilde{m_p}, \widetilde{m_n})}\,I^{(\widetilde{m_p}, \widetilde{m_n})}(E) \;. \label{eq46ann1}\tag{60}\] Eq. (60 ) is exact and here the \(d\)’s are dimensions, \(\rho\)’s are normalized to unity and \(I\)’s are normalized to the dimensions. Without loss of generality, from now on, we will denote \((\widetilde{m_p}, \widetilde{m_n})\) by \({\widetilde{m}}\), \(d(m_p,m_n)\) by \(d\) and \(\rho^{m_p, m_n)}(E)\) by \(\rho(E)\). The moments \(M_p({\widetilde{m}})\) of \(\rho^{({\widetilde{m}})}(E)\) are \(M_p({\widetilde{m}}) = \left\langle H^p\right\rangle^{{\widetilde{m}}}\) with the centroid \(\epsilon({\widetilde{m}})=\left\langle H\right\rangle^{{\widetilde{m}}}\) and the variance \(\sigma^2({\widetilde{m}}) =\left\langle H^2\right\rangle^{{\widetilde{m}}} -[\epsilon({\widetilde{m}})]^2\). Now, using the \(q\)-normal form for \(\rho^{({\widetilde{m}})}(E)\) we have, assuming \(q\) is independent of \({\widetilde{m}}\), \[I^{(m_p,m_n)}(E) = \displaystyle\sum_{{\widetilde{m}}} d({\widetilde{m}}) \rho^{({\widetilde{m}})}(E) \approx \displaystyle\sum_{{\widetilde{m}}} \displaystyle\frac{d({\widetilde{m}})}{\sigma({\widetilde{m}})}\,f_{qN}^{({\widetilde{m}})}({\hat{E}}({\widetilde{m}})|q) \label{eq46ann2}\tag{61}\] with \({\hat{E}}({\widetilde{m}})=(E-\epsilon({\widetilde{m}}))/\sigma({\widetilde{m}})\). The \(f_{qN}\) is defined over the interval \[\left(\epsilon({\widetilde{m}})-\displaystyle\frac{2\,\sigma({\widetilde{m}})}{\displaystyle\sqrt{1-q}}\;,\;\epsilon({\widetilde{m}})+\displaystyle\frac{2\,\sigma({\widetilde{m}})}{\displaystyle\sqrt{1-q}}\right)\;.\] In practice we can use the EGOE formula for \(q\) for \(p-n\) systems deducing via \(\left\langle H^4\right\rangle^{(m_p,m_n)}\) formula given in [52] and averaging its value for one, two and three-body \(H\). An alternative is to use \(q\) as a free parameter. It is also possible to derive for a realistic \(H\) the exact formula for \(\left\langle H^4\right\rangle^{(m_p,m_n)}\) using the methods given in [3], [7] but this is challenging with three-body forces. As \(\rho^{({\widetilde{m}})}(E)\) is in fact a strength function, a better approximation incorporating the results in Section 3.3 is [37] \[\begin{array}{l} \rho^{({\widetilde{m}})}(E) = \left\{\sigma({\widetilde{m}})\right\}^{-1}\;f^{({\widetilde{m}})}_{qN}({\hat{E}}|q) \left[1+\displaystyle\frac{\gamma_1({\widetilde{m}})}{[3]_q!} He_3({\hat{E}}|q) +\displaystyle\frac{(\gamma_2({\widetilde{m}})+1-q)}{[4]_q!} He_4({\hat{E}}|q)\right]\;;\\ \\ \gamma_1({\widetilde{m}}) \approx -(1-q)\left[\displaystyle\frac{\epsilon({\widetilde{m}})-E_c(m)}{\sigma({\widetilde{m}})}\right]\;,\\ \\ \gamma_2({\widetilde{m}}) = (q-1) + (1-q)^2 \left[\displaystyle\frac{\epsilon({\widetilde{m}})-E_c(m)}{\sigma_{{\widetilde{m}}}}\right]^2 + (1-q^2) \displaystyle\frac{\sigma^2_h(m)}{\sigma^2_V(m)}\;. \end{array}\label{eq46ann3}\tag{62}\] It is possible to use Eqs. (60 )-(62 ) to calculate state densities. Alternatively, by replacing \({\widetilde{m}}\) by \({\widetilde{m}}\,J\) everywhere will give level densities. However, calculating \(J\) dependent centroids, variances and \(q\) value, in particular with 3-body forces, is computationally intensive. Before going further, it is important to mention that Fig. 1 shows that \(f_{qN}\) is good for \(k\)-body interactions and Fig. 6 shows that \(f_{qN}\) is good also for \((1+k)\)-body interactions. However, as we have \((1+2+3)\)-body interactions in nuclei, we show in Fig. 9 an example demonstrating that \(f_{qN}\) is good for these also.
Transition strengths and transition strength sums are in some situations measurable carrying new nuclear structure information and more importantly a theory for these is needed for many applications such as in calculating \(\beta\)-decay rates in Astrophysics, neutrinoless double \(\beta\)-decay transition matrix elements and so on. The SSM theory for transition strengths given in [9], [10], [30] can be extended to incorporate the bivariate \(q\)-normal in place of bivariate Gaussian form. Also the needed bivariate correlation coefficient for different types of transition operators follow from [52]. Some discussion of all these is given in [37].
Given a transition operator \({\cal O}\), a formula for transition strength sum \(\left\langle{\cal O}^\dagger {\cal O}\right\rangle^E\) (i.e. for the sum of the transition strengths originating from an eigenstate of \(H=h+V\) with energy \(E\)) is, to a good approximation given by, \[\left\langle{\cal O}^\dagger{\cal O}\right\rangle^{E} = \displaystyle\sum_{(\widetilde{m_p}, \widetilde{m_n})}\;\displaystyle\frac{I^{(\widetilde{m_p},\widetilde{m_n})}(E)}{I^{(m_p ,m_n)}(E)}\;\left\langle{\cal O}^\dagger{\cal O}\right\rangle^{(\widetilde{m_p},\widetilde{m_n})}\;. \label{eq46ann4}\tag{63}\] Formulas for \(\left\langle{\cal O}^\dagger {\cal O}\right\rangle^{(\widetilde{m_p},\widetilde{m_n})}\) can be written down for a variety of one and two-body operators [4]. With \({\cal O}^\dagger\) a creation operator, Eq. (63 ) gives shell model orbit occupancies and the formula is exact for occupancies. Also, note that \(\left\langle{\cal O}^\dagger{\cal O}\right\rangle^{E} I^{(m_p,m_n)}(E)\) gives transition strength density. In applying Eq. (63 ), the densities \(I\)’s will be replaced by the corresponding \(q\)-normal densities as described in Section 5.1. Besides the non-energy strength sum, also important are the lower order energy weighted strength moments, i.e. the centroid, variance and skewness, of the distribution of strengths originating from an eigenstate with energy E, as they are also measurable in many situations in nuclei. For example, the moments \(M_p(E)\) are given by \[M_p(E) = \left\{\left\langle{\cal O}^\dagger{\cal O}\right\rangle^E\right\}^{-1}\;\displaystyle\sum_{E_f} |\left\langle E_f \mid {\cal O}\mid E\right\rangle|^2\,(E_f)^p \label{eq46ann5}\tag{64}\] and then the centroid \(\epsilon(E)=M_1(E)\) and the variance \(\sigma^2(E) = M_2(E) -(M_1(E))^2\). Similarly, the skewness \(\gamma_1(E)\) is defined via \(M_3(E)\). As these are moments of the conditional density of the bivariate transition strength density, their variation with \(E\) follows from Eqs. (31), (33) and (35). Then, \(\epsilon(E)\) will be linear in \(E\), the variance \(\sigma^2(E)\) is a constant (does not depend on \(E\)) and \(\gamma_1(E)\) will be linear in \(E\) with negative slope. All these results are tested in a numerical example and the results are shown in Figs. 10(a)-(d). In the calculations, used is a EGOE\((1+2+3)\) ensemble with (\(N = 12\), \(m = 6\), \(\lambda_2 = 0.3\) and \(\lambda_3 = 0.2\)). For the transition operator \({\cal O}\), chosen is the one-body operator \(a^\dagger_2 a_9\). In Fig. 10(a), numerical results (histogram and stars) are compared with analytical curves for strength sum density and strength sum. They follow from Eq. (64 ) using \(m\)-particle averages. Then, the strength sum density is a marginal of the bivariate transition strength density and thus, it follows \(q\)-normal distribution and this corresponds to the smooth curve in the figure. Similarly, strength sum is the ratio of strength sum density and state density. As the transition operator is not completely random and the Hamiltonian operator has a fixed one-body part along with a mixture of two and three body rank operators, there is a shift of the centroid of the strength sum density relative to the state density centroid. Going to Figs. 10(b)-(d), it is clearly seen that the strength centroid \(\epsilon(E)\), variance \(\sigma^2(E)\) and skewness \(\gamma_1(E)\) follow from the equations for the moments of the conditional \(q\)-normal distribution; see Eqs. (31), (33) and (35). The agreements with theory are very good with some deviations at the spectrum edges.
In conclusion, it is important to add that the Gaussian form used in the past in SSM is reasonably good as long as the systems considered have sufficiently large number of particles and the Hamiltonian is \((1+2)\)-body. However, with growing knowledge on 3-body (perhaps also 4-body) interactions in nuclei, certainly in future one needs the formulation, with \(q\)-normal forms, as briefly described in this Section. In future, it is important to carry out tests of SSM with \(q\)-normal forms using shell model codes with realistic \((1+2+3)\)-body Hamiltonians. In addition, it is also important to carry out applications to nuclear level densities, astrophysical reaction rates calculations and so on.
Embedded random matrix ensembles, introduced 50 years back, continue to be of interest in nuclear physics in particular and in quantum many-particle physics in general. The EE provide the basis for SSM approach for nuclear structure. More recently (from 2017), following the RMT results, for the so called SYK model involving Majorana fermions, as derived by Verbaarschot and collaborators, a new direction in exploring EE and SSM has opened. Using the formulas for the lower order moments of eigenvalue densities, transition strength densities and strength functions on one hand and a variety of numerical calculations on the other, it is now well established that FEE(\(k\)) and BEE(\(k\)) indeed generate \(q\)-normal forms. These results are described in Sections 3.1-3.3. The \(q\)-normal forms and some of their properties are given in Section 2 for easy reference. Further, presented are some results, in Sections 4.1 and 4.2 showing the role of the \(q\) parameter in level fluctuations in the bulk and also in ground state (lowest eigenvalue) fluctuations. Although analytical results showing \(q\)-normal forms are available only for FEE(\(k\)) and BEE(\(k\)), in nuclear physics applications in particular (for SSM), it is important to consider \(k\)-body interactions in presence of a mean-field one-body term. As shown using numerical examples in Section 4.3, with strong enough \(k\)-body interaction strength, the EE(1+k) also generate \(q\)-normal forms for various densities. Further, the \(q\)-normal forms are also good with the more realistic \((1+2+3)\)-body interactions (see Figs. 9 and 10).
With the \(q\)-normal forms well established, clearly it is necessary to modify SSM formulations used before by replacing Gaussians with \(q\)-normal forms. Some aspects of SSM with \(q\) normal forms and the associated \(q\)-Hermite polynomials are presented in Section 5 with more details given in [37]. In applying SSM with \(q\) normal forms, a technical problem that need to be solved is in extending the codes in [19] to calculate also \(\left\langle H^4\right\rangle^{(\widetilde{m_p}, \widetilde{m_n})J}\) or at least \(\left\langle H^4\right\rangle^{(m_p m_n)J}\). These will give the \(q\) parameter values with realistic nuclear interactions. Another important issue that need to be addressed (this will extend the scope of SSM) is to solve the embedded random matrix that includes multi-\(\hbar\omega\) mixing. For this one has to consider for example the partitioned FEGOE described in Section 13.3 and Fig. 13.3 in [30]. It is expected that these ensembles will give multi-modal distributions. Here, it is important to mention that after the NTSE-2026 meeting, there appeared a preprint presenting a more formal derivation of the \(q\)-normal forms for FEE(\(k\)) and BEE(\(k\)) using a method based on ‘Wick product of non-commuting Gaussian random variables’ [75]. It remains to be seen if this new method gives a solution to the partitioned FEGOE and also solve the two-pint correlation function describing level fluctuations. We hope that this review will lead to further investigations of EE and SSM and their applications. Clearly, these future explorations and application need HPC (High-Performance Computing).
Thanks are due to S. Tomsovic for some useful correspondence. NDC acknowledges financial support from University Research Project No.DR/Dir./26-27/17/Sr No-5.