Single-Particle Entanglement Dynamics in Complex Systems

TL;DR

This paper explores how single-particle entanglement entropy changes in complex systems based on system parameters, revealing a universal functional that governs these dynamics.

Contribution

The paper introduces a universal functional that captures the dynamics of single-particle entanglement entropy in complex systems.

Findings

A single functional of system parameters governs the entropy dynamics.

The entropy is sensitive to collective information content rather than individual system details.

The functional identifies universality classes and connections between quantum states.

Abstract

We analyze the effect of varying system conditions on the single-particle entanglement entropy for an arbitrary eigenstate of a complex system that can be described by a multiparametric Gaussian ensemble. Our theoretical analysis leads to the identification of a single functional of the system parameters that governs the entropy dynamics. This reveals a sensitivity of the entropy to collective information content, characterized by the functional, instead of the individual system details. The functional can further be used to identify the universality classes as well as a deep web of connection underlying different quantum states.

Figures

Distribution of SPEE for Anderson Hamiltonian. The figure illustrates the SPEE distribution for a typical eigenstate at energy E=0 of an Anderson Hamiltonian for four different combinations of systems parameters w,w1,k for a fixed L=12. The combinations are chosen such that Equation (25) gives almost same Λ value for them. The results for three different Λ are displayed in three panels, with both (Λ,〈SA〉) mentioned in the legend and the labels (a), (b), (c), and (d), referring to different combinations of w1 and k: (a) k=1,w1=0.0, (b), k=2,w1=0.0, (c) k=1,w1=1.0, and (d) k=2,w1=1.0. The large deviation tails for the SPEE distribution for small and intermediate Λ values are clearly visible here.

Evolution of SPEE measures for the RP ensemble. The figure illustrates the Λe-governed evolution of the (a) average and (b) variance of SPEE for the RP ensemble (with μ=cNα, where c and α are the free parameters) for a fixed matrix size N=123. For comparison, the inset displays the α-governed evolution for different c values. As indicated by the collapse of curves already in terms of Λe, χ seems to be irrelevant here. This suggests an almost independence of χ from c, which thereby plays an insignificant role in the comparison of RPE with different c values. A consideration of χ is, however, relevant when compared with a different system, i.e., AE (see Figure 10). Here again the statistical uncertainties are small compared to symbol size, thus making error bars irrelevant.

Full text

1. Introduction

The ubiquitous complexity in real systems highlights the concept of entropy as an important tool for understanding behavior in diverse fields, e.g., classical thermodynamics, statistical physics, and information theory. Indeed, with its basic definition based on the randomness and uncertainty in the system, an underlying connection among its various forms, e.g., thermal, statistical, and information, is intuitively expected. This is also corroborated by many recent studies; e.g., the study [1] indicates that a consideration of the relationship between thermal and entanglement entropy can place significant constraints on the ground state entanglement structure for “natural” systems. Another study [2] indicates that the entanglement entropy of a topologically ordered system measures the maximal quasiparticle fluctuations and, as a consequence, corresponds to the thermal entropy of the quasiparticles at infinite temperature on the entanglement boundary. The potential connection with thermal entropy motivates us to analyze an important form of the entanglement entropy, namely ‘single-particle entanglement entropy’ (SPEE), developed to treat the ubiquitous indistinguishable particles.

The standard tools of entanglement rely on the tensor product structure of the state space of a composite quantum system. While such an underlying structure is always present in the case of distinguishable particles, the constraints on the particle statistics in the case of indistinguishable ones confine the available Hilbert space, thus rendering the standard approach inapplicable. To overcome this issue, an alternative tool was proposed in the past: it is based on an isomorphism (nonunique) of the Fock space to the state space of a composite (many modes) quantum system [3]. The idea in turn motivated the concept of single-particle entanglement (SPE): consider a particle confined on a bipartitioned lattice of arbitrary dimensions. The single-particle state of the lattice Hamiltonian can then be expressed as a tensor product of the site occupation number basis states in the second-quantized Fock space. The site entanglement in the two parts of the lattice can then be utilized as a criterion to analyze the aspects of single-particle wave dynamics where traditional criteria such as the inverse participation ratio, etc., are inadequate. Some examples in this context are wave dynamics under anisotropic system conditions, localization-to-delocalization phase transition, etc. Indeed, application of the entanglement statistics for characterization of the phase transitions has been studied in many previous studies [4,5,6,7,8,9], but most of them have been confined to bipartite or multipartite entanglement measures. Comparatively, less attention has been paid to single-particle entanglement analysis [10,11], which, albeit counterintuitive, has been proven to be a genuine resource in quantum information processing [3,12]. In this work, we pursue the analysis further by developing a common mathematical formulation of single-particle entanglement entropy dynamics for a wide range of Hamiltonians; the results are numerically verified for two prototypical Hamiltonians, viz., the three-dimensional (3D) Anderson model, representing a particle moving in a random potential with its hopping confined to some nearest neighbors, and the Rosenzweig–Porter ensemble, representing long-range interactions in the basis space.

The first and foremost requirement for the entanglement analysis of a physical system is the knowledge of its Hamiltonian in the physically motivated basis and thereafter its eigenstates. A realistic system is, however, almost always many-body, consisting of complicated interactions among its subunits. Even if the latter are well-known, e.g., Coulomb interactions, an exact determination of the matrix representation of the Hamiltonian is often not possible; this could occur, for example, due to technical issues involved in determination of the matrix elements (e.g., calculation of the integrals either by a theoretical or numerical route). Incomplete knowledge or error in their determination manifests itself by randomization of the matrix elements, which can, however, be of different types. For example, if a matrix element can be determined up to its average value and variance, the maximum entropy hypothesis predicts its distribution to be Gaussian. In case information is available about higher-order moments, it can lead to non-Gaussian distributions. Similarly, the existence of conservation laws can result in correlated distributions or other constraints on the matrix elements, e.g., column/row sum rule. In addition, not all elements need be random; some of them can be exactly determined. For example, for the Hamiltonian of a tight-binding lattice with on-site disorder and nearest neighbour non-random hopping, the diagonals in the site basis are randomly distributed but the off-diagonals are all non-random (given by the hopping strength). Indeed, the nature and type of the distribution of the matrix elements is sensitive to various system conditions, e.g., symmetry and conservation laws, dimensionality and boundary conditions, disorder, etc., and can vary from one element to another. As a consequence, the Hamiltonian matrix is best represented by a system-dependent random matrix, with some or all elements randomly distributed. The information about the system appears through the distribution parameters of the matrix elements. As expected, the randomness underlying the matrix elements of the Hamiltonian also manifests in its eigenstates, and the distribution of the latter can be derived, in principle, from a multivariable integration over the Hamiltonian ensemble density (i.e., the JPDF of the matrix elements); this is achieved by a transformation of variables from the matrix space to the eigenvalue–eigenstate space and thereafter an integration over all eigenvalues and all other eigenfunctions except one of them. The information in turn leads to statistics of the entanglement measures of the remaining eigenstate.

A multivariable integration of the ensemble density of the Hamiltonian is in general technically complicated. In addition, the interactions among various subunits of a many-body system are in general sensitive to a host of system conditions, both of static as well as dynamic type. A variation in these conditions can lead to changes in the mutual interactions among subunits and thereby the ensemble density. This in turn is expected to manifest in the behavior of a typical many-body state and thereby entanglement measures and can lead to serious technical issues: even if the integration is achieved by some approximations, these are often specific to system parameters and may not be applicable once the parameters are changed. The integration route is therefore not often a viable option and motivates a search for alternative tools that take into account arbitrary variation in the system parameters.

The search for an alternative tool to determine quantum correlations for a general class of Hamiltonians represented by multiparametric Gaussian ensembles [13,14,15] motivates us to seek a common mathematical formulation for the single-particle entanglement statistics. As discussed in [15], the joint probability distribution function (JPDF) of the components of an arbitrary eigenstate (referred to as the state JPDF hereafter for brevity) for Hermitian operators represented by multivariate Gaussian ensembles undergoes a multiparametric evolution as the ensemble parameters vary. For a specific variation in the latter, the evolution describes a diffusion in state space and can be exactly derived from the ensemble density. The diffusion equation has an additional benefit: it provides a common mathematical formulation, referred as the complexity parameter formulation, of the state JPDF for a wide range of Hamiltonians where the system information is contained in a single functional of all system parameters. We apply the above formulation in the present work and thereby derive the evolution equation for the single-particle entanglement measures and their solutions.

We proceed as follows. We begin in Section 2 by introducing the standard definitions of the single-particle entanglement entropy (SPEE). This is followed by Section 3, describing two prototypical examples of a large class of single-particle Hamiltonians for which our formulation is intended. Section 4 revisits the complexity parameter formulation, briefly explaining how the changing system conditions lead to an evolution of the JPDF of the matrix elements and in turn of the eigenfunction components. The JPDF dynamics of the eigenfunction components is later used to obtain the evolution equation of the single-particle entanglement in Section 5. Finally, in Section 6, we present the exact diagonalization results of the examples introduced in Section 3 as a validation of our theoretical claims and conclude in Section 7 with a summary of our results and future directions. The theoretical derivations of our results involve detailed technical background material, already published several times; to keep this paper self-contained but to avoid its defocusing, we have moved these details to the Supplementary Material. For the reader’s convenience, we also include a brief description of the symbols used in this paper in Appendix A.

2. Single-Particle Entanglement Measure

Consider the dynamics of a single particle in a lattice of arbitrary dimension D, bipartitioned into parts A and B, created, e.g., by slicing the 3D lattice horizontally, vertically, or in an arbitrary direction as per requirement (as illustrated in Figure 1 by a schematic diagram). With each subpart associated with an occupation probability of the particle, it can be described by state $[eqn]$ (subpart occupied) or $[eqn]$ (subpart empty). An arbitrary state of the particle can then be written as a product over the occupation probabilities in the two parts. Based on whether the particle occupies the sites in subsystem A or B, two such product states are possible. An eigenstate of the Hamiltonian H can then be written as a superposition of the product states,

$$[eqn]$$

But as the probability $[eqn]$ or $[eqn]$ of the particle lying in subpart A or B depends on the number of sites therein, this implies $[eqn]$ and $[eqn]$ and similarly for B. The state $[eqn]$ can then be written as

$$[eqn]$$

where $[eqn]$ and $[eqn]$ correspond to the occupation or vacancy of site r, respectively.

The formulation for the standard Von Neumann entropy of a pure quantum state is analogous to the one for the Shannon entropy. Based on similar ideas, the single-particle entanglement entropy (SPEE) is defined in terms of the probability of occupation of the subparts. Let $[eqn]$ be the probability of the particle occurring in part A,

$$[eqn]$$

The SPEE can then be defined as

$$[eqn]$$

where $[eqn]$ . In the ergodic limit, for a balanced bipartition implying equal occupancy of both subparts, we have $[eqn]$ and thereby $[eqn]$ . The other extreme limit corresponds to the particle being localized in one subpart, thereby implying $[eqn]$ or $[eqn]$ ; both these possibilities correspond to $[eqn]$ . Consequently, $[eqn]$ varies between 0 and $[eqn]$ as the particle wavefunction starts extending from one subpart to another.

3. Hamiltonian and Ensemble Density

Equation (4) gives the SPEE for an eigenfunction of a single-particle Hamiltonian. However, for cases where the Hamiltonian is best described by an ensemble, it is imperative to consider the distribution of the eigenfunction over the ensemble and therefore SPEE statistics. The latter can be derived if the theoretical formulation for the ensemble density of Hamiltonians representing the particle is available. The natural query arising then is how to find the appropriate as well as mathematically tractable ensemble representing a given Hamiltonian. The next section elucidates our approach to construct the ensembles for two prototypical single-particle Hamiltonians; these are later used in the numerical analysis too.

(i) Anderson Hamiltonian: Consider a single electron moving in a random potential, e.g., a D-dimensional disordered lattice. Within tight-binding approximation, the dynamics can be described by the Hamiltonian [10,11,16,17]

$$[eqn]$$

where $[eqn]$ describes the tunneling amplitude of an electron between sites k and l and the summation $[eqn]$ is over z neighboring sites. Based on the nature of the disorder, the ensemble representing H can be of various types. For example, for the independent Gaussian-distributed site energies $[eqn]$ with variance $[eqn]$ , the probability density of $[eqn]$ can be written as $[eqn]$ . The hopping can be chosen to be isotropic or anisotropic and non-random or random (Gaussian). In the case of non-random, anisotropic hopping $[eqn]$ between nearest neighbours, we have $[eqn]$ if $[eqn]$ form nearest neighbour pairs and $[eqn]$ if $[eqn]$ are not nearest neighbours. The probability density $[eqn]$ of the ensemble can therefore be given as

$$[eqn]$$

with C as a normalization constant; $[eqn]$ denotes the product over connected sites only (depends upon the number of neighbors). Similarly, in the case of random nearest neighbor hopping, we have

$$[eqn]$$

Using the Gaussian limit of the delta function, i.e., $[eqn]$ , Equation (6) and Equation (7) can both be written in the form of a multiparametric Gaussian ensemble,

$$[eqn]$$

where $[eqn]$ for disconnected sites and the notation $[eqn]$ corresponds to the limit $[eqn]$ for all $[eqn]$ -pairs corresponding to the non-random hopping or for disconnected sites.

(ii) Rosenzweig–Porter Ensemble: The Rosenzweig–Porter ensemble (RPE), a prototypical model often used for spectral analysis of the localization → delocalization transition [16,18,19,20], is described by the probability density [16]

$$[eqn]$$

i.e., a Hamiltonian with normally distributed diagonal entries and off-diagonal entries with variance controlled by an external parameter $[eqn]$ :

$$[eqn]$$

We note that a typical Hamiltonian taken from the above ensemble describes, e.g., a long-range isotropic hopping in a lattice with on-site disorder. As discussed in [16], a choice of $[eqn]$ , with c and $[eqn]$ as arbitrary parameters, reveals the rich dynamics of a typical wavefunction of the Hamiltonian; it undergoes a localization → extended non-ergodic → ergodic transition of wavefunction dynamics with changing $[eqn]$ for a fixed N [16]. For $[eqn]$ and 2, the spectral statistics is known to be critical and different from both Poisson (localized limit) and Wigner–Dyson (ergodic limit). The eigenfunction statistics in the bulk of the spectrum for these $[eqn]$ -values is known to be multifractal.

4. Complexity Parameter Formulation of the Hamiltonian Ensemble

Our primary focus in the present work is to study the entanglement statistics of a typical eigenstate of a complex system represented by the Hamiltonian H in the single-particle approximation. As the examples discussed in the previous section indicate, a wide range of Hamiltonians can be represented by an ensemble of N-dimensional real-symmetric matrices H, described by an ensemble density,

$$[eqn]$$

with $[eqn]$ and $[eqn]$ as the arbitrary variance and mean value of the matrix element $[eqn]$ and C as the normalization constant. As the limit $[eqn]$ corresponds to a non-random $[eqn]$ taking a value $[eqn]$ , the above density can model a wide range of ensembles, including sparse as well as banded ones. As mentioned in Section 1, for $[eqn]$ to be an appropriate representation of the Hamiltonian, the ensemble parameters are expected to be sensitive to the system conditions. Some examples of such system parameters are disorder, dimensionality, boundary and topological conditions, system size, etc. Indeed, the ensemble densities in Equations (8) and (9) correspond to special cases of Equation (11).

The ensemble density in Equation (11) is a function of $[eqn]$ matrix elements $[eqn]$ and $[eqn]$ ensemble parameters $[eqn]$ . A variation in system conditions affects various interactions among subunits of the system and changes not only $[eqn]$ but also their uncertainties and thereby distribution parameters. This in turn leads to changes in $[eqn]$ both in the matrix space as well as in the ensemble parameter space; the location of each matrix of the ensemble in the matrix space is changed and thereby the moments $[eqn]$ are averaged over the shifted ensemble. The dynamics give rise to many natural queries, a few of which are mentioned below.

(i) Whether the dynamics of $[eqn]$ in the ensemble parameter space has an exact equivalent in the matrix space. This information is relevant for following reason: We assume that the system under consideration is well-described by the ensemble in Equation (11) at a given instant of time. A change in the system conditions can, however, affect the matrix elements. The question is whether the system can still be described by the same ensemble with time-evolved ensemble parameters? Indeed, from Equation (11), a variation in either the matrix elements or the ensemble parameters can change the ensemble density. For an evolving ensemble to continue representing the system, it is desirable that the evolution of $[eqn]$ in H-matrix space is exactly mimicked by that in $[eqn]$ space, where $[eqn]$ refer to the variance and mean matrices, consisting of $[eqn]$ and $[eqn]$ , respectively.

(ii) Does the dynamics indeed depend on individual details of all ensemble parameters or only respond to their collective influence? This query arises from ample evidence in many areas of complex systems that complexity, irrespective of its origin, is sensitive to very few details or relevant system parameters, often resulting in universality of the physical properties.

(iii) Whether the new ensemble still retains its Gaussian form (requires specific conditions on $[eqn]$ ) and whether it be accessed from the initial ensemble, in the ensemble parameter space, just by varying the parameters $[eqn]$ and $[eqn]$ . On intuitive grounds, this is expected to be achieved only for specific dynamics in the matrix space.

4.1. Diffusion Equation for ρ(H)

Indeed, as discussed in detail in [13] (also in Section SI of the Supplementary Material), a specific drift dynamics of $[eqn]$ in $[eqn]$ -space, generated by an operator $[eqn]$ , can exactly mimic a diffusive dynamics with a finite drift in H-matrix space: $[eqn]$ , where $[eqn]$ . Here $[eqn]$ , $[eqn]$ , and $[eqn]$ and $[eqn]$ is a non-zero but otherwise arbitrary constant giving the variance of the matrix elements at the end of the evolution [13,16].

The above equation describes a multiparametric flow of the matrix elements in a real-symmetric matrix space from an arbitrary initial condition, say $[eqn]$ . As discussed in detail in [13,16], the dynamics of $[eqn]$ in $[eqn]$ -space can be mapped to another M-dimensional parametric space $[eqn]$ , in which $[eqn]$ undergoes a single parametric evolution if following set of conditions is satisfied:

$$[eqn]$$

As is clear from the above, $[eqn]$ is the evolution parameter and $[eqn]$ are the constants of the evolution. To distinguish from the $[eqn]$ -space, hereafter the t-space will be referred as the complexity parameter space; its dimensionality M depends on the number of original ensemble parameters $[eqn]$ undergoing variation, and in general $[eqn]$ [13].

The parameters $[eqn]$ can be obtained by solving the characteristic set of equations

$$[eqn]$$

with $[eqn]$ for $[eqn]$ and $[eqn]$ for $[eqn]$ , with $[eqn]$ referring to the number of non-zero $[eqn]$ .

A particular solution of the above equation can be given as

$$[eqn]$$ $$[eqn]$$

with $[eqn]$ as arbitrary constants subjected to the condition $[eqn]$ . While many solutions satisfying the above condition are possible, the appropriate solution is the one that is also applicable to the initial ensemble (as the constants of evolution are same for the initial ensemble). Indeed, the choice of constants $[eqn]$ in Equation (14) depends on the evolving ensemble parameters. For example, if one of the $[eqn]$ remains zero throughout the evolution and does not participate in the evolution, the corresponding coefficient $[eqn]$ can be set to zero. For the case in which all $[eqn]$ remain non-zero throughout the evolution, we have $[eqn]$ , and it is appropriate to choose $[eqn]$ (for all $[eqn]$ ). The choice fulfills the condition $[eqn]$ , and $[eqn]$ can be given as

$$[eqn]$$

To distinguish it from $[eqn]$ , hereafter $[eqn]$ is symbolically referred to as Y. As Y acts as the evolution parameter while $[eqn]$ are the constants of evolution, the former can appropriately be referred to as the ensemble complexity parameter and the latter as complexity constants.

Being a sum over logarithmic functions of $[eqn]$ and $[eqn]$ , Y can be regarded as an average distribution parameter, i.e., an average uncertainty that governs the statistics as system conditions change. Another definition of Y follows from an alternative derivation discussed in [21], implying it to be the average time-scale associated with accuracy fluctuations due to changing system conditions.

The transformation from the space of ensemble parameters to complexity parameters implies the transformation of $[eqn]$ in Equation (11) to $[eqn]$ (the latter given by Equation (10) in Section SI of the Supplementary Material). The ensembles with different $[eqn]$ sets can in general have different values of $[eqn]$ and therefore are represented by different points in the complexity parameter space. The common mathematical framework given by Equation (12) for the wide range of ensembles represented by $[eqn]$ in Equation (11) suggests an analogous formulation for the eigenvalues and eigenfunctions and their classification into an infinite range of universality classes defined by a single parameter Y for a given set of the constants of evolution. The latter in turn depend on the matrix constraints, e.g., whether the H-matrix is subjected to the same symmetry or conservation laws. The ensembles with the same set of $[eqn]$ (i.e., same type of matrix constraints) lie along the same evolutionary path in the complexity parameter space and can be accessed by a variation in their ensemble parameters.

To achieve our objective to derive a Y-based formulation of the single-particle entanglement statistics of a typical Hamiltonian represented by Equation (11), we now require a similar formulation for the eigenfunction components. This is briefly reviewed below (the detailed derivation is discussed in [15] and also in the Supplementary Material).

4.2. Diffusion Equation for the Joint Probability Density of the Eigenfunctions and Eigenvalues

As the single-particle entanglement measures for a Hamiltonian are defined in terms of the statistics of its eigenfunction components, it is necessary to first derive the Y-governed evolution of the latter. This was derived in [15] from an exact route, i.e., an exact diagonalization of the evolution equation for $[eqn]$ . Here we present an alternative approach (although based on second-order perturbation theory, it can be generalized to non-Gaussian cases of $[eqn]$ too).

Consider an $[eqn]$ real-symmetric matrix H with O as the $[eqn]$ eigenvector matrix of H (an orthogonal matrix satisfying $[eqn]$ ) and E as the $[eqn]$ diagonal matrix of its eigenvalues, $[eqn]$ . A small perturbation $[eqn]$ of H in the matrix space changes its eigenvalues and eigenfunctions. But as discussed above, with H represented by Equation (11), this also results in a complexity parameter-governed evolution of $[eqn]$ and thereby the joint probability density $[eqn]$ of the eigenvalues and eigenfunctions; here $[eqn]$ and $[eqn]$ refer to the sets of all eigenfunctions $[eqn]$ (the columns of the U-matrix) and eigenvalues $[eqn]$ (the diagonal of the E-matrix). The derivation of the Y-governed evolution of $[eqn]$ basically depends on three main steps (a detailed discussion is given in Section SII of the Supplementary Material).

(i) We first determine the first- and second-order moments of $[eqn]$ from Equation (12). A comparison of the latter with the standard Fokker–Planck equation gives $[eqn]$ , $[eqn]$ and $[eqn]$ , with $[eqn]$ and $[eqn]$ for $[eqn]$ . All other averages are of a higher order in $[eqn]$ and can be ignored for a small change in Y.

(ii) The information from step (i) is then used, along with second-order perturbation theory for Hermitian operators, to derive the moments of the eigenvalues and eigenfunctions of H. We have $[eqn]$ and $[eqn]$ , with angular brackets implying conditional ensemble averages with fixed $[eqn]$ . Similarly, the moments for the eigenvalues can be given as, up to the first order of $[eqn]$ , $[eqn]$ and $[eqn]$ . Further, to the first order in $[eqn]$ , the ensemble-averaged correlation between $[eqn]$ and $[eqn]$ is zero (for both $[eqn]$ and 2): $[eqn]$ .

(iii) Based on a similar dynamics for the matrix elements as indicated by Equation (12), we assume a Markovian dynamics for the eigenvalues and eigenfunctions of H. The Y-governed evolution of the joint probability density $[eqn]$ can then be described by a standard Fokker–Planck equation,

$$[eqn]$$

where $[eqn]$ and $[eqn]$ refer to two parts of the Fokker–Planck operator, corresponding to eigenvalues and eigenfunction components, respectively. Here $[eqn]$ is given as

$$[eqn]$$

and $[eqn]$ is

$$[eqn]$$

Note here that $[eqn]$ is subjected to the following boundary condition for $[eqn]$ : $[eqn]$ for $[eqn]$ for $[eqn]$ ; this follows because the higher-order moments of the ensemble density are assumed to be negligible.

A substitution of the moments in Equation (17) followed by the latter’s integration over all undesired variables will then lead to an evolution equation for the joint probability density of the desired combination of eigenfunction components. This is briefly described below (discussed in detail in [15] and also in Section SIII of the Supplementary Material).

4.3. Diffusion Equation for the Joint Probability Density of the Components of a Typical Eigenfunction

The JPDF of the components $[eqn]$ of an eigenstate, say $[eqn]$ of H lying between $[eqn]$ and $[eqn]$ with $[eqn]$ , can be defined as

$$[eqn]$$

where $[eqn]$ is the volume element in the eigenvalue–eigenvector space.

A differentiation of Equation (20) with respect to Y gives $[eqn]$ . Using Equation (17) in the integral and repeated partial integrations over the eigenvalues and eigenfunctions other than $[eqn]$ then leads to (see Section SIII of the Supplementary Material)

$$[eqn]$$

with $[eqn]$ and $[eqn]$ . As is clear from the above, the evolution is now governed by

$$[eqn]$$

with $[eqn]$ as an important system-specific spectral range defined as follows: the eigenvalues at distances more than $[eqn]$ around e are uncorrelated. In general, it is of the order of a few local mean level spacings $[eqn]$ at the energy e of interest: $[eqn]$ . Here $[eqn]$ is the number of eigenvalues in the range, intuitively expected to be related to the inverse participation ratio $[eqn]$ at energy e, i.e., $[eqn]$ . This suggests that

$$[eqn]$$

with $[eqn]$ as the proportionality constant with units of $[eqn]$ and $[eqn]$ as the local mean level spacing at energy $[eqn]$ . In the context of disordered systems, it is believed that $[eqn]$ ∼ $[eqn]$ , with $[eqn]$ as the Thouless energy, and is of the order of $[eqn]$ at energy e. We also recall that for many cases, e.g., a d-dimensional Anderson Hamiltonian with nearest neighbor hopping, $[eqn]$ is a ratio of system size N and localization volume $[eqn]$ , with $[eqn]$ as the average localization length at the energy e.

The solution of Equation (21) gives the JPDF of the eigenfunction components for a fixed $[eqn]$ and for an arbitrary initial condition. For finite system sizes, changing system parameters can change $[eqn]$ continuously between 0 and ∞; this implies, for a finite N, the existence of an infinite range of universality classes of the JPDF characterized by continuous values of $[eqn]$ between 0 and ∞. But as $[eqn]$ itself is N-dependent, a variation in system conditions in an infinite size limit leads to, in general, $[eqn]$ or ∞, with the solution of Equation (21) approaching the two end points, i.e., the initial state and Wigner–Dyson ensembles, respectively [15]. For specific system conditions, however, a subtle competition between $[eqn]$ and $[eqn]$ may render $[eqn]$ independent of N; referring to the corresponding $[eqn]$ value as $[eqn]$ , it would then remain finite even in the $[eqn]$ limit: $[eqn]$ [15]. The latter in turn implies, in an infinite size limit, the existence of a new intermediate statistics at $[eqn]$ , different from the two end points (i.e., $[eqn]$ and ∞) and can therefore be referred to as critical. Indeed, based on the complexity of the system, more than one set of system parameters may exist, resulting in different $[eqn]$ and thereby different critical statistics intermediate between the initial state and Wigner–Dyson ensembles. In contrast to finite N, Equation (21) therefore indicates, in the $[eqn]$ limit, the existence of discrete universality classes for the eigenfunction statistics, each characterized by a distinct $[eqn]$ . The explicit appearance of size N in Equation (21) (in addition to its implicit appearance through $[eqn]$ ), however, suggests that the statistics at $[eqn]$ for finite N is different from that of infinite N.

A detailed study in [13,16] of the eigenvalue dynamics of the ensemble (11) with changing ensemble parameters revealed the existence of a single-parameter formulation of the spectral statistics. Referred as the spectral complexity parameter, it is defined as $[eqn]$ [13,15]. Due to the N-dependence of both $[eqn]$ and $[eqn]$ , the spectral statistics shows a finite size scaling as well as a critical behavior at energy e if the system conditions lead to an N-independent, non-zero, finite $[eqn]$ , referred to as $[eqn]$ , that is, if the size dependence of $[eqn]$ is canceled by $[eqn]$ [16]. Equivalently,

$$[eqn]$$

As the correlations among eigenvalues and eigenfunctions are known to be significant in the critical regime, it is useful to express $[eqn]$ in terms of $[eqn]$ ,

$$[eqn]$$

Here $[eqn]$ is again a function of $[eqn]$ as well as energy e (see Section SIII of the Supplementary Material),

$$[eqn]$$

To distinguish it from $[eqn]$ and $[eqn]$ , $[eqn]$ hereafter is referred to as the strength complexity parameter.

The relation in Equation (25) gives an important insight: it suggests the existence of a finite size scaling of the wavefunction statistics too. The latter in turn manifests through a multifractal behavior of the eigenfunctions at the critical point, with scaling exponents (referred to as critical exponents or multifractal dimensions) dependent on the system parameters [15,22]. As $[eqn]$ is sensitive to system specifics, the critical (multifractal) exponents can vary from system to system. The occurrence of a critical point and corresponding multifractal behavior of the eigenstates is, however, not a necessary feature of all infinite-size complex systems; based on the $[eqn]$ -formulation, it requires a specific set of system conditions conspiring with each other and leading to a size-independent $[eqn]$ (a more detailed discussion of $[eqn]$ is included in the Supplementary Material).

While the derivation of Equation (21) presented above is based on the second-order perturbation approach for Hermitian matrices and the Markovian assumption, it can be derived exactly by a direct integration of Equation (12) over all eigenvalues and eigenfunctions other than $[eqn]$ [15]. The alternative approach presented here, however, lends further credence to our complexity parameter formulation.

5. Evolution of the Single-Particle Entanglement with System Parameters

Writing components $[eqn]$ of $[eqn]$ as $[eqn]$ (for notational simplification as well as to retain same notations as used in [15]), we now have $[eqn]$ and $[eqn]$ . Equation (21) can now be used to derive the first two moments of the SPEE as follows.

5.1. Dynamics of SPEE Average

With $[eqn]$ defined in Equation (4), its ensemble average can now be obtained by averaging over the ensemble density $[eqn]$ of the eigenfunction $[eqn]$ ,

$$[eqn]$$

where $[eqn]$ . The above definition, along with Equation (21) implying the $[eqn]$ -dependence of $[eqn]$ , also implies a similar dependence for $[eqn]$ , and its evolution can be derived from Equation (21). The steps are as follows. From Equation (21), we have

$$[eqn]$$

A substitution of Equation (21) on the right side of the above equation leads to $[eqn]$ , where

$$[eqn]$$ $$[eqn]$$

The above integrals can be solved by repeated partial integrations, but this requires terms dependent on how $[eqn]$ responds to a small change in the eigenfunction components. Here we give the required rates of change for the analysis,

$$[eqn]$$

and

$$[eqn]$$

With repeated partial integration of $[eqn]$ and $[eqn]$ and subsequent use of the relations in Equations (31) and (32), Equation (28) can be further simplified (discussed in detail in Section SV of the Supplementary Material). As our focus here is on the large N limit, it can be rewritten as

$$[eqn]$$

with $[eqn]$ and $[eqn]$ . A general solution of the above equation can now be given as

$$[eqn]$$

with $[eqn]$ as the constant of integration to be determined from the initial condition. Noting that either $[eqn]$ or $[eqn]$ in the separability limit, the terms $[eqn]$ and $[eqn]$ in the integral in Equation (34) diverge at the lower integration limit $[eqn]$ if the initial state is separable. Choosing a weak separability limit, with $[eqn]$ or $[eqn]$ with $[eqn]$ , as the initial condition, however, gives $[eqn]$ for $[eqn]$ . Equation (34) then implies $[eqn]$ . Using $[eqn]$ , this gives $[eqn]$ in the limit $[eqn]$ .

As is clear from Equation (33), the determination of the $[eqn]$ -dependence of $[eqn]$ requires prior knowledge of the similar dependence of $[eqn]$ and $[eqn]$ . But this in turn requires knowledge of the negative moments of $[eqn]$ and leads to a set of hierarchical equations with no available solution (this can be seen by following the similar steps used for the derivation of $[eqn]$ ). To overcome this technical difficulty, the only option available at this stage is to consider a crude approximation, i.e., replace $[eqn]$ and similarly for $[eqn]$ in Equation (33) (the approximation is often used for partition function studies of systems with annealed randomness). The $[eqn]$ -dependence of $[eqn]$ can now be determined, proceeding as in the $[eqn]$ -case. By definition, we have

$$[eqn]$$ $$[eqn]$$

where $[eqn]$ . Substitution of Equation (21) in the above equation gives

$$[eqn]$$

Further, by repeated partial integration and realizing that the non-zero terms after differentiation in $[eqn]$ are such that $[eqn]$ , we get

$$[eqn]$$

As a check, we note that, for the balanced bipartition case $[eqn]$ , the above equation gives expected limiting behavior $[eqn]$ as $[eqn]$ . Assuming the separability limit as the initial condition, the solution of Equation (38) can be given as

$$[eqn]$$

With $[eqn]$ , the above equation gives $[eqn]$ and $[eqn]$ ∼1 and $[eqn]$ for $[eqn]$ and $[eqn]$ , respectively, thus indicating a rapid transition as a function of $[eqn]$ with $[eqn]$ as the point of inflection. A smooth crossover around this point can then be seen as a function of $[eqn]$ ; here we recall that $[eqn]$ , defined in Equation (26), also changes with $[eqn]$ .

Substitution of $[eqn]$ and $[eqn]$ in Equation (33) now gives

$$[eqn]$$

with $[eqn]$ defined in Equation (26). As $[eqn]$ depends on the system parameters, taking different values in localization and delocalization regimes, this in turn leads to two possible solutions. But as is clear from the above, a rescaling of $[eqn]$ by $[eqn]$ is expected to map one solution on the other.

As discussed in the next section, while the numerical result for $[eqn]$ has the same form as the one given in Equation (40), it does not exactly match. We believe this discrepancy arises from the approximation $[eqn]$ (based on the numerics for $[eqn]$ discussed in the next section, the approximation seems to be good for large $[eqn]$ ). To rule out the possibility that the discrepancy arises from the approximation $[eqn]$ , we consider the following intuitive approach to improve our theoretical prediction. For small $[eqn]$ ranges, assuming that $[eqn]$ are perturbed slightly from their $[eqn]$ values, i.e., $[eqn]$ with $[eqn]$ , we have $[eqn]$ . For large $[eqn]$ ranges where typically $[eqn]$ , we have $[eqn]$ . This suggests a rapid decay of $[eqn]$ from a divergence to a constant value with increasing $[eqn]$ . Based on the above insights, we conjecture, for the balanced case ( $[eqn]$ ) and for $[eqn]$ -governed evolution of $[eqn]$ from a separable initial state, that

$$[eqn]$$

with $[eqn]$ . The above conjecture is in agreement with the balanced case numerics for $[eqn]$ and $[eqn]$ (Figure 2, numerical details discussed in the next section). As a measure of goodness of fit, we find that the $[eqn]$ value (defined in Appendix B) for the above conjecture in Figure 2 is about $[eqn]$ and the standard error in determining parameters $[eqn]$ and $[eqn]$ is about $[eqn]$ and $[eqn]$ , respectively. Substitution of Equation (41) in Equation (34) gives (for the $[eqn]$ case)

$$[eqn]$$ $$[eqn]$$

with $[eqn]$ as the Beta function and $[eqn]$ as the confluent Hypergeometric function. As discussed in the next section, notwithstanding a good agreement between $[eqn]$ numerics and the conjecture in Equation (41), a comparison of Equation (43) with $[eqn]$ -numerics does not show any significant improvement over Equation (40).

For $[eqn]$ , the evolution approaches the steady state with $[eqn]$ with $[eqn]$ . Equation (40) then gives $[eqn]$ . This is consistent with the prediction based on Equation (4). This also lends credence, indirectly, to the approximation $[eqn]$ and similarly for $[eqn]$ used in Equation (33).

5.2. Dynamics of SPEE Variance

The evolution equation for the variance can similarly be obtained. Here we give the equation for the simple case $[eqn]$ (also used in our numerics),

$$[eqn]$$

where $[eqn]$ . The variance of the SPEE is thus determined by the covariance of $[eqn]$ and $[eqn]$ . Proceeding again as in the case of Equation (33), we obtain

$$[eqn]$$

with $[eqn]$ as the constant of integration to be determined from the initial condition. A choice of the separability condition as the initial condition implies $[eqn]$ for the initial state; Equation (45) then implies $[eqn]$ .

As in the case of $[eqn]$ , here again determination of $[eqn]$ involves negative moments and requires solving a hierarchical set of equations, which is technically complicated. This motivates us again to consider an intuitive approach. As mentioned above, for $[eqn]$ with $[eqn]$ (which is the case for small $[eqn]$ ranges), we have $[eqn]$ and $[eqn]$ . This gives $[eqn]$ ∼ $[eqn]$ . This is also valid for the case with $[eqn]$ . For large $[eqn]$ ranges where typically $[eqn]$ , we have $[eqn]$ and $[eqn]$ ; this again suggests that $[eqn]$ ∼ $[eqn]$ . Based on the above insights, we conjecture that $[eqn]$ as a function of $[eqn]$ has the same mathematical form as $[eqn]$ as a function of $[eqn]$ (with $[eqn]$ ):

$$[eqn]$$

with $[eqn]$ as a rescaling factor. As displayed in Figure 2, the conjecture is corroborated by our numerics for the Anderson case (system details given in Section 6), which gives

$$[eqn]$$ $$[eqn]$$

with $[eqn]$ . As can be checked by a direct substitution, the above functional forms are consistent with Equation (45) too. Here, $[eqn]$ is a rescaling of $[eqn]$ to match the mean with the numerics, and $[eqn]$ and $[eqn]$ for the variance and covariance, respectively. The absolute error in the determination of the normalization factor a from its exact value $[eqn]$ is about $[eqn]$ and $[eqn]$ for the variance and the covariance curves, respectively.

5.3. Relevance of Λ and Complexity Constants

Equations (34) and (45) describe the common mathematical formulation for the SPEE statistics of a typical eigenfunction of the Hamiltonians modeled by Equation (11), with information about the ensemble parameters and thereby system parameters appearing through $[eqn]$ only. This implies that the SPEE statistics respond to changing system parameters in a collective way instead of individually to each of them. The information about crucial system parameters may appear through $[eqn]$ as well as $[eqn]$ . For example, the dimensionality and boundary conditions both influence the basis connectivity, i.e., degree of sparsity of the matrix, which is also reflected in the distribution parameters $[eqn]$ . The dimensionality also affects $[eqn]$ through the average localization length $[eqn]$ (the sensitivity of $[eqn]$ to the system parameters of Anderson Hamiltonians is discussed in detail in [15,16]).

The system size N is another important parameter that affects the evolution of the eigenfunction statistics significantly and thereby the related measures. The dependence on N can, however, vary from one measure to another [15]. Indeed, in contrast to Equation (21), the size N appears in SPEE statistics through $[eqn]$ . For large $[eqn]$ , the statistics is therefore expected to approach the separable state behavior in the limit $[eqn]$ and the maximum SPEE in the limit $[eqn]$ . For $[eqn]$ , the SPEE statistics is predicted to become size-independent and remains different from both the separable state as well as the maximum entangled state, even in an infinite size limit.

Equations (34) and (45) describe the evolution of the SPEE statistics from a separable initial state at $[eqn]$ . While the latter can arise from different sets of initial conditions for different systems, their path of the evolution is fixed by the complexity constants. Thus, in general, SPEE statistics for different Hamiltonians may evolve along different paths in the complexity parameter $[eqn]$ -space. It is, however, possible that two different Hamiltonians, represented in two different bases, still have the same sets of complexity constants $[eqn]$ for both of them. Typical eigenstates of such Hamiltonians will evolve, with $[eqn]$ as the evolution parameter, along the same path in t-space, but again that does not ensure same SPEE statistics for them. As discussed in the previous section, the latter requires an analogy of $[eqn]$ and can further be explained as follows: Consider two different systems labeled as “1” and “2”, both represented by the ensemble (11), with their different system conditions giving rise to different sets of ensemble parameters labeled as “ $[eqn]$ ” and “ $[eqn]$ ”. If the physical conditions in “1” and “2” are such that the relation

$$[eqn]$$

is fulfilled, the statistical behavior of system “1” at energy $[eqn]$ is predicted to be analogous to that of “2” at energy $[eqn]$ . The above analogy, however, need not imply an analogy for previous histories of the two systems (indeed they may show completely different behavior for $[eqn]$ ) or for their future or even at other energies. We emphasize that the above universality is applicable for the Hamiltonian matrices belonging to same matrix constraint class, e.g., real-symmetric or complex Hermitian, etc., which in turn depends on the underlying symmetries and conservation laws (referred to as global constraints) [21]. (We recall that the matrices belonging to same matrix constraint class may belong to ensembles with different ensemble parameters (referred to as ensemble constraints or “local” constraints).

A relevant query in the above context is whether the complexity constants $[eqn]$ can always be chosen the same for different Hamiltonian ensembles (of type represented by Equation (11)) if they are subjected to the same set of global constraints. Intuitively the answer is in the affirmative and can be justified as follows: a physically motivated basis to represent a Hamiltonian is usually based on its global constraints, e.g., symmetry conditions and conservation laws, which the chosen basis is expected to preserve under perturbation of the system. For example, for a time-reversal Hamiltonian, an appropriate basis is the one in which the Hamiltonian matrix is real-symmetric. Different Hamiltonians with a common set of global constraints can then be represented in same/ equivalent basis if their dynamics of interest preserves them. For example, both Anderson and RP Hamiltonians discussed in Section 3 are time-reversal invariant but subjected to no other symmetry constraint. As our interest is in their evolution, which at no stage violates the time-reversal symmetry, both can then be well represented by real-symmetric matrices (this is, however, not the case, e.g., if the Hamiltonian dynamics violates the time-reversal symmetry or is subjected to column constraints). A typical matrix, represented in a fixed basis, in general has many basis constants [23,24], which can then be chosen as the complexity constants $[eqn]$ . This is further explained through two examples (namely Anderson and RP) in Section SV of the Supplementary Material. (It is worth reemphasizing that two different Hamiltonians, both represented by real-symmetric matrices, need not belong to ensemble $[eqn]$ (Equation (11)) with analogous $[eqn]$ sets; indeed the relative matrix elements’ strengths and the ensemble parameters depend on the local constraints, e.g., disorder, interaction, etc. This in turn results in different complexity parameters, e.g., $[eqn]$ , for them, even though their complexity constants are same. While a variation in their respective system conditions may force such Hamiltonians to follow the same evolutionary path in complexity space, their statistics in general need not be same.

6. Numerical Analysis

The theoretical analysis described in the previous section predicts the existence of an infinite range of universality classes of SPEE statistics, characterized by the complexity parameter $[eqn]$ , among the eigenstates if their Hamiltonians are represented by Equation (11) and belong to same symmetry class. The robustness of our theoretical prediction makes it imperative to seek their numerical verification. Here we numerically analyze the two ensembles described by Equations (8) and (9).

Anderson Ensemble (AE): We exactly diagonalize the $[eqn]$ Hamiltonian matrix, represented in the site basis, of a three-dimensional cubic Anderson lattice (Equation (5)) subjected to a periodic boundary condition for two linear lattice sizes $[eqn]$ and 18 (with $[eqn]$ ). The matrix for each N is considered for several on-site disorder strengths w, off-diagonal isotropic random hopping with variance $[eqn]$ (with its mean t fixed to $[eqn]$ ), and the hopping range parameter k. Here $[eqn]$ implies hopping to six nearest neighbours ( $[eqn]$ ) and $[eqn]$ implies hopping to twenty-four next nearest neighbours ( $[eqn]$ ). Throughout this analysis, the hopping rate is kept fixed to $[eqn]$ . Equations (8) and (16) give the ensemble complexity parameter for the case $[eqn]$ , with $[eqn]$ , and $[eqn]$ , with $[eqn]$ . The initial condition of the Hamiltonians in the ensemble is chosen to be the localized phase such that $[eqn]$ . Thus, we have

$$[eqn]$$

with $[eqn]$ as the $[eqn]$ value corresponding to the initial state. The spectral complexity parameter $[eqn]$ for this case can now be calculated from Equation (25), with $[eqn]$ [15] (also discussed in Section SVII of the Supplementary Material). To determine the average localization length $[eqn]$ , here, we use the relation $[eqn]$ (valid in the spectrum bulk near $[eqn]$ ), with $[eqn]$ as the inverse participation ratio (IPR) of the eigenfunction of interest, and $[eqn]$ as the mean level spacing; both $[eqn]$ and $[eqn]$ are determined numerically.

For each set of parameters $[eqn]$ , we consider an ensemble of $[eqn]$ disorder realizations and use the standard shift-invert diagonalization technique to calculate $[eqn]$ of the total eigenfunctions in a small range $[eqn]$ at the energy $[eqn]$ [25,26]. This in turn gives the total $[eqn]$ eigenfunctions for an average analysis of SPEE for a given set of system conditions; the averaging here can then consist of both ensemble as well as spectral averaging. An important point worth emphasizing here regards the choice of spectral range $[eqn]$ : based on our theoretical predictions, $[eqn]$ (Equation (25)) governs the SPEE statistics and its quantitative analogy is necessary to observe the analogy of the statistics for different parametric combinations. But $[eqn]$ is energy-dependent too, which makes it necessary in principle to consider only ensemble averaging. To improve the statistics, however, it may be imperative to consider a spectral averaging over a small optimized range $[eqn]$ such that $[eqn]$ does not vary significantly over the range.

For calculating the SPEE measures, we consider a horizontal bipartition of the lattice through the center (Figure 1) such that $[eqn]$ . Throughout the numerical analysis, we have used log base 2 to calculate the SPEE and its statistics. Furthermore, the figures, except Figure 1, are shown on a semi-log plot. To see the effect of the type of averaging, we analyze the results both ways, i.e., (i) ensemble averaging only and (ii) ensemble–spectral averaging (over various disorder realizations and over about $[eqn]$ of the total eigenstates per realization at $[eqn]$ ) for a fixed set of system parameters. The above analysis is repeated by varying the parameters, which gives us the average and variance of SPEE for many $[eqn]$ values. Our analysis indicates no significant difference in the results obtained by the two types of averaging. This is, however, due to the choice of the eigenfunction for the analysis; we have analyzed the eigenfunction at $[eqn]$ . (The local spectral density and, as a result, $[eqn]$ for AE does not vary significantly within a small energy range around $[eqn]$ [16].) To avoid repetition of almost the same figures, here we include the figures corresponding to ensemble averaging only. The results for $[eqn]$ and its variance are displayed in the top and bottom panels in Figure 3.

To illustrate the relevance of $[eqn]$ as a single parameter governing the SPEE dynamics (instead of $[eqn]$ , which plays that role for spectral statistics [16]), it is important to first understand what happens if the dynamics is studied in terms of one of the system parameters. The insets in the top and bottom panels of Figure 3 display the variation with respect to diagonal disorder w for a fixed lattice size $[eqn]$ (i.e., $[eqn]$ ) but for four different combinations of $[eqn]$ and k. As can be seen from both panels of the figure, the curves corresponding to lattices with nearest neighbour interaction ( $[eqn]$ , denoted by $[eqn]$ ) but with different hopping types (random hopping with off-diagonal variance $[eqn]$ as well as non-random hopping $[eqn]$ both for fixed hopping strength $[eqn]$ ) collapse onto a single curve in terms of $[eqn]$ . The curves, however, collapse to a different curve if the next nearest neighbour interaction ( $[eqn]$ , denoted by $[eqn]$ ) is switched on while keeping the hopping variances the same ( $[eqn]$ and $[eqn]$ ). This is consistent with the $[eqn]$ -based formulation: the SPEE statistics for two different system parameter sets is expected to be analogous only in terms of $[eqn]$ and not $[eqn]$ . As the ensemble-averaged localization length $[eqn]$ and thereby $[eqn]$ for the nearest neighbor differs significantly from that of next neighbour hopping, $[eqn]$ for the two cases is expected to be a different constant, which in turn leads to two different curves. To determine $[eqn]$ , however, prior knowledge of $[eqn]$ for the $[eqn]$ and 2 cases is needed. The numerical analysis for $[eqn]$ gives an analogous $[eqn]$ -variation with w for the two cases (see Figure S1 in the Supplementary Material). Equation (26) then implies $[eqn]$ , with subscripts 1 and 2 implying the variable for cases $[eqn]$ and 2, respectively. Due to a lack of exact information about $[eqn]$ , $[eqn]$ , we conjecture that $[eqn]$ (based on numerical estimation); the latter suggests that $[eqn]$ . Indeed, as can be seen from Figure 4, a rescaling of $[eqn]$ for the $[eqn]$ case by a factor of 12 (keeping all other details the same as in Figure 3) leads to convergence of the corresponding statistics to the $[eqn]$ case, thereby confirming a single-parametric evolution in terms of $[eqn]$ . The display in Figure 3 also provides another insight, i.e., an insensitivity of $[eqn]$ to $[eqn]$ (as revealed by the analogy of $[eqn]$ -governed curves for a fixed k but two different $[eqn]$ values); this insight may be potentially relevant for a theoretical $[eqn]$ -formulation in the future.

In the previous section, we predicted Equation (40) as a theoretical formulation for $[eqn]$ , assuming $[eqn]$ . Figure 4 displays a comparison of the numerics based on the exact SPEE formulation (defined in Equation (4)) with Equation (40). While the latter seems to predict an almost similar functional form for $[eqn]$ , a quantitative deviation between the numerics and approximate theory is clearly visible. The latter is expected because $[eqn]$ seems to be a good approximation only in a large $[eqn]$ range (as displayed in Figure 5 for one set of parameters, we find that $[eqn]$ ). To improve the prediction, we conjectured Equation (41) for $[eqn]$ for the case $[eqn]$ , which in turn led to Equation (43) for $[eqn]$ . A comparison with $[eqn]$ numerics displayed in Figure 2 shows a good agreement with Equation (41) for $[eqn]$ and $[eqn]$ . This motivates us to compare the $[eqn]$ numerics with Equation (43) using the same $[eqn]$ and $[eqn]$ . As displayed in Figure 4, however, Equation (43) does not lead to any significant improvement over Equation (40). We believe that the deviation could be an artifact of computations on the logarithmic scale as well as finite system sizes, while the theoretical predictions are derived in the large N-limit. We also considered two fitted functions for large and small $[eqn]$ ; these results are also displayed in the inset of Figure 4.

In order to understand the effect of system size on SPEE statistics, we consider its dynamics for a larger lattice size $[eqn]$ (i.e., $[eqn]$ ) while keeping the other parameter details the same as in Figure 3. However this time, since we generated only 100 samples, we use both ensemble and spectral averaging. The corresponding $[eqn]$ -governed SPEE behavior (for $[eqn]$ ) is displayed in Figure 6. Similarly to Figure 3, the curves again converge to a single curve in the main panel as a function of $[eqn]$ (while a lack of convergence in the inset for w-governed evolution is still visible). A comparison of the display in the main panels of Figure 6 for $[eqn]$ and Figure 4 for $[eqn]$ also indicates an analogy of their statistics as a function of $[eqn]$ (also displayed in Figure 7 for a single set of system conditions).

While Figure 4 and Figure 6 validate the $[eqn]$ formulation of the SPEE statistics for four different $[eqn]$ combinations, the computational limitations permit us to analyze each combination only for two fixed system sizes. It is therefore natural to query the finite size effects on the $[eqn]$ formulation and whether it remains valid even in an infinite size limit. We recall that theoretically a finite $[eqn]$ value in the limit $[eqn]$ implies a critical wavefunction statistics different from the initial state and the ergodic limit (discussed in Section 4.3). We also recall, with $[eqn]$ given below Equation (50), that $[eqn]$ for a d-dimensional Anderson lattice of linear size L depends on the ratio $[eqn]$ (with $[eqn]$ ), thereby implying a finite size scaling of the spectral statistics (verified in [16]). This in turn suggests, from Equation (25), a finite size scaling of the wavefunction statistics too, although the relevant length scale, say $[eqn]$ , competing with system size L need not be same as the average localization length $[eqn]$ of the eigenfunctions at energy e. Intuitively this can be explained as follows. With usual eigenfunction correlation measures based on the intensities at a few sites, their scaling as a function of $[eqn]$ can be explained as a competition between the tendency to spread over whole basis space and the tendency to localize to a few basis states. In the SPEE case, however, the collective response of the intensities in each subsystem affects the correlations between them. Noting that the intensities in a subsytem need not belong to neighboring basis states, the emergence of new length scales dominating such quantum correlations is therefore not unexpected. While $[eqn]$ is the relevant length scale for the localization-to-delocalization transition in the basis space, $[eqn]$ is expected to play a similar role for the separability-to-maximum entanglement transition between the two subsystems. Indeed, as displayed in Figure 7, a comparison of SPEE statistics as a function of w (for a fixed $[eqn]$ combination), for four different system sizes $[eqn]$ with $[eqn]$ and under the balanced condition $[eqn]$ , indicates a finite size scaling, with a critical disorder $[eqn]$ . Although the curves for different system sizes are quite close here (seemingly an artifact of the balanced condition $[eqn]$ chosen in our analysis), their crossing at $[eqn]$ is clearly visible; the latter therefore marks the critical point of the separability-to-maximum entanglement transition. The determination of $[eqn]$ and its relation with $[eqn]$ requires, however, a detailed theoretical as well as numerical analysis and is currently under study.

Figure 7 reveals another interesting insight: $[eqn]$ behavior for all N values collapse onto a single curve as a function of $[eqn]$ , thereby indicating that the statistics are governed only by $[eqn]$ (similarly for the $[eqn]$ -variance) irrespective of system size. The above analogy for different system sizes (for a fixed set of disorder and hopping parameters) along with those displayed in Figure 3 and Figure 6 (for different disorder and hopping parameters but at a fixed system size) lend further credence to our $[eqn]$ parameter-based theoretical formulation of the SPEE measures, i.e., Equation (34) and Equation (44). Recalling that these equations are applicable for arbitrary $[eqn]$ with $[eqn]$ , it is desirable to numerically establish that the $[eqn]$ -based analogy extends to the case $[eqn]$ too. The collective response of the basis states of one subsystem to another also hints at the importance of relative subsystem sizes in the potential emergence of $[eqn]$ and thereby for the finite size scaling of SPEE. Indeed the existence of such a length scale was reported in the study [4], which was based on a finite size scaling analysis of the von Neumann entropy for engineered bipartite Gaussian states (EPGSs). A finite size scaling of the average SPEE for the $[eqn]$ case was also reported in [11], with one of the subsystems consisting of a single site. A detailed analysis of the case $[eqn]$ and its impact on $[eqn]$ is currently under progress and will be discussed in a subsequent work.

The numerically observed success of the $[eqn]$ -formulation for the first two moments of the SPEE statistics encourages one to seek up to what order of moments it would persist. An appropriate route to answer this query is by directly analyzing the SPEE distribution over the ensemble $[eqn]$ of M Hamiltonians. Although our present theoretical analysis is confined to the average and variance, we numerically analyze the SPEE distribution for an Anderson state at $[eqn]$ and for many combinations of systems parameters $[eqn]$ . To improve the statistics, we consider a spectral averaging of $[eqn]$ over $[eqn]$ eigenfunctions in the neighbourhood of $[eqn]$ . We recall that many different combinations of systems parameters can lead to the same $[eqn]$ value (following from Equation (16) and Equation (25)) and thereby same SPEE statistics. To verify this in case of the SPEE distribution, we again consider four different sets of $[eqn]$ , which lead to approximately the same $[eqn]$ ; the results for three different $[eqn]$ values are displayed in Figure 8. As can be seen from the figure, the distributions for different system parametric combinations almost collapse to one curve if they share approximately the same $[eqn]$ value. (An exact matching, although desirable, requires a detailed analysis of the average localization length $[eqn]$ and therefore transfer matrix analysis, which is time-consuming and will be pursued elsewhere.) The figure also reveals large deviation tails for the SPEE distribution for small and intermediate $[eqn]$ values and appears similar to the one reported for the bipartite entanglement entropy of a typical state of the random-field Heisenberg model [6]. Indeed, the average entanglement entropy of an arbitrary bipartite quantum state with Gaussian-distributed components (with arbitrary variances and mean values) reported in [4,5] also displays a behavior visibly similar to that of the average SPEE discussed in the current work.

RP Ensemble: We next consider the Rosenzweig–Porter Hamiltonian, with its ensemble density described by Equation (9) with $[eqn]$ , where N is the size of the matrix. The initial state of the ensemble is chosen as $[eqn]$ , which corresponds to the localized dynamics of the eigenfunctions and thereby separable state of the two sub-lattices. We again consider the eigenstates at $[eqn]$ and calculate the SPEE statistics for different free parameters c and $[eqn]$ . Equation (16) gives in this case

$$[eqn]$$

Here again we set $[eqn]$ . The spectral complexity parameter $[eqn]$ can now be determined from Equation (25). We note that, in contrast to the AE case, the local mean level spacing $[eqn]$ at energy e is now given as $[eqn]$ (explained in Section SVII of the Supplementary Material), with $[eqn]$ as the ensemble-averaged spectral density at the energy e and $[eqn]$ as the mean level spacing.

The $[eqn]$ dependence of the average of $[eqn]$ and its variance over the ensemble $[eqn]$ is displayed in Figure 9, with the inset showing the corresponding dependence on $[eqn]$ for different c. As evident from the figure, the statistics evolves along different paths for different c values, with $[eqn]$ as the evolution parameter; it follows the same path, however, in terms of $[eqn]$ , evolving rapidly from separable to the maximum limit. But again recalling that the evolution paths for different system parameters are expected to converge in terms of $[eqn]$ , this suggests an almost analogous $[eqn]$ -variation with $[eqn]$ for all c values considered in our analysis.

Distribution of SPEE for Anderson Hamiltonian. The figure illustrates the SPEE distribution for a typical eigenstate at energy E=0 of an Anderson Hamiltonian for four different combinations of systems parameters w,w1,k for a fixed L=12. The combinations are chosen such that Equation (25) gives almost same Λ value for them. The results for three different Λ are displayed in three panels, with both (Λ,〈SA〉) mentioned in the legend and the labels (a), (b), (c), and (d), referring to different combinations of w1 and k: (a) k=1,w1=0.0, (b), k=2,w1=0.0, (c) k=1,w1=1.0, and (d) k=2,w1=1.0. The large deviation tails for the SPEE distribution for small and intermediate Λ values are clearly visible here.

Evolution of SPEE measures for the RP ensemble. The figure illustrates the Λe-governed evolution of the (a) average and (b) variance of SPEE for the RP ensemble (with μ=cNα, where c and α are the free parameters) for a fixed matrix size N=123. For comparison, the inset displays the α-governed evolution for different c values. As indicated by the collapse of curves already in terms of Λe, χ seems to be irrelevant here. This suggests an almost independence of χ from c, which thereby plays an insignificant role in the comparison of RPE with different c values. A consideration of χ is, however, relevant when compared with a different system, i.e., AE (see Figure 10). Here again the statistical uncertainties are small compared to symbol size, thus making error bars irrelevant.

To gain insights into higher-order statistics for the RP case, here again we numerically analyze the SPEE distribution for the eigenstates chosen from the neighborhood of $[eqn]$ and for many $[eqn]$ combinations. The results displayed in Figure 11 indicate, in contrast to the Anderson case, an absence of the large deviation tails for the RP case in the small $[eqn]$ regime. The behavior in the other regimes for the two cases is, however, analogous. Indeed, such behavior is expected because the diffusion equation for the distribution is governed by a single parameter only in the large $[eqn]$ limit. (With our primary focus on the average $[eqn]$ in the current work, the mathematical details for the diffusion equation are omitted here, but it can be derived following similar steps to those given in [5] for engineered bipartite Gaussian states (EPGSs).)

As mentioned near Equation (49), the analogy of $[eqn]$ -governed evolution of SPEE measures is not only predicted for a specific Hamiltonian under different system conditions, but also expected to occur for two different Hamiltonians if both belong to same class of global constraints, e.g., symmetries and conservation laws ([15]). As the matrices in both AE and RP ensembles are real-symmetric, subjected to no additional constraints, the complexity constants for both systems can be chosen same (see the discussion in Section SV of the Supplementary Material). This suggests that, under system conditions leading to $[eqn]$ , the SPEE statistics for AE and RP should be analogous. In addition, the $[eqn]$ -governed paths of their SPEE evolution are theoretically expected to be analogous too. Considering same N for both systems, Equation (26) gives $[eqn]$ and thereby $[eqn]$ ∼ $[eqn]$ (with subscripts “ae” and “rp” implying the measures for AE and RP, respectively). Here again due to lack of knowledge about the constants $[eqn]$ and $[eqn]$ , we determine them by a numerical estimation. (The latter is based on first noting, from Figure 3 and Figure 9, the values of $[eqn]$ and $[eqn]$ that correspond to $[eqn]$ . With the latter equality predicted to occur at $[eqn]$ , this implies the relation $[eqn]$ , which in turn gives $[eqn]$ . Substitution of numerically determined values of $[eqn]$ and $[eqn]$ for both cases now gives the desired ratio.)

Figure 10 displays a comparison of $[eqn]$ statistics for RP ensemble and Anderson ensemble for one set of their system parameters (sufficient as collapse for different sets is already displayed); the good agreement of the two curves in top panel, with $[eqn]$ as the independent variable, and its absence in lower panel, with $[eqn]$ as the independent variable, once again emphasizes the role of $[eqn]$ as the relevant single parameter for SPEE-statistics (instead of $[eqn]$ ). Based on the above analogy, a finite size scaling behavior for the SPEE statistics is expected in the RP case too; this analysis is currently under progress.

The $[eqn]$ -based analogy depicted in Figure 10 corresponds to the eigenstates chosen near energy $[eqn]$ for both AE and RP ensembles. But, as mentioned near Equation (49), the analogy need not be confined to $[eqn]$ only and is theoretically predicted for other energies too (i.e., the two systems at the same or different energies) if the condition in Equation (49) is satisfied. A numerical analysis of the above prediction is currently being pursued in our follow-up work on many-body entanglement and we hope to report it in the near future. (Meanwhile we note that a detailed theoretical as well as numerical verification of $[eqn]$ -based universality for the spectral statistics at different energies of chiral Hamiltonians is discussed in [27].)

7. Conclusions

In this work, we have analyzed the SPEE dynamics in disordered Hamiltonians as the system conditions vary. The theoretical approach presented here is based on the complexity parameter formulation of the eigenfunction statistics of the Hamiltonians represented by multiparametric Gaussian ensembles; the formulation was developed and investigated in detail in [15] based on an exact diagonalization of the diffusion equation for the ensemble density. Hoping to extend the potential applicability of our formulation to non-Gaussian cases too, we have also included here an alternative derivation based on second-order perturbation theory of quantum Hamiltonians and the Markovian dynamics assumption.

Our theory predicts the existence of a common mathematical formulation of the SPEE for the typical eigenstates of a wide range of single-particle Hamiltonians, represented by different multivariate Gaussian ensembles but subjected to the same symmetry constraints and conservation laws. As the ensemble parameters contain information about the system details and are therefore different for ensembles representing different Hamiltonians, the system dependence appears in the formulation collectively through a single parameter $[eqn]$ , a functional of the system parameters. This in turn implies an analogy of the SPEE statistics for different quantum states (same or different Hamiltonians) if they share same $[eqn]$ values and thereby reveals a great deal of universality hidden underneath the world of quantum states. While our SPEE formulation is based on an exact theoretical formulation for the JPDF of the eigenfunction components, the robustness of our claim renders numerical verification desirable. In particular, we have tested out theoretical predictions for many system conditions of two prototypical Hamiltonians, viz., the 3D Anderson model and the Rosenzweig–Porter ensemble. A minor drawback of our analysis at this stage is the approximate information for $[eqn]$ ; while crucial to determine $[eqn]$ , this, however, is a challenging technical problem due to eigenvalue–eigenfunction correlations, with insights seemingly possible only through numerical analysis. We recall that the spectral complexity parameter $[eqn]$ is related to $[eqn]$ through $[eqn]$ . As discussed in detail in [15,22], $[eqn]$ plays an important role in determining the critical spectral statistics and the multifractal exponents, relevant criteria for the localization-to-delocalization transitions in complex systems. The $[eqn]$ dependence of the SPEE statistics (through $[eqn]$ ) can then be used to study the latter’s connection to multifractality exponents; we hope to pursue this study in future.

Our SPEE analysis presented above reveals a potentially useful insight about the entanglement entropy in quantum systems. The prevailing notion in the research community about the SPEE is that it is at best a measure of the quantum correlations with classical counterparts and does not have as mysterious an origin as the bipartite entanglement entropy, which characterizes the correlations with no classical counterpart. A comparison of our SPEE results with those for bipartite entanglement entropy reported in [4,5,28], however, indicates a similar dependence on the complexity parameter (the studies [4,5,28] were based on the bipartite quantum states described by multivariate Gaussian ensembles). A similar dependence is also revealed by our current study of the bipartite entanglement entropy of many-body states; we expect to report it in the near future. This is not surprising: it has been known for last few decades that the complexity, irrespective of its diverse origin, often leads to independence from the system specifics. In addition, our preliminary analysis of the finite size scaling of SPEE measures suggests the existence of a new length scale relevant for separability to the single-particle entanglement transition. With the latter expected to be different from the average localization length that governs the typical localization-to-delocalization transition of the eigenstates, the information may potentially lead to new insights not only into the quantum dynamics of single-particle Hamiltonians but for many-body cases too.

A common mathematical formulation for the entanglement dynamics in single-particle systems has by itself many implications. The first and foremost among them is the existence of an infinite range of universality classes of the SPEE statistics for finite systems, characterized by the strength complexity parameter  $[eqn]$ , with the latter taking continuous values between 0 and ∞. For infinite system sizes, the size dependence of $[eqn]$ , however, results in its rapid approach either to 0 or ∞. If the subsystems are initially chosen to be separable, the SPEE either approaches the separability limit or the maximum entanglement limit. If, however, a subtle conspiracy of the system conditions can lead to size-independent $[eqn]$ for an arbitrary system size, or a non-zero, finite $[eqn]$ even in an infinite size limit, the statistics is then predicted to be different from both limits. In general, in an infinite size limit, only discrete values of non-zero, finite $[eqn]$ are possible; this implies only a partial entanglement of the two subsystems at the energy under consideration. Statistically, however, this entanglement behavior need not be unique and may be achieved by a change in system conditions as well as energy for a fixed pair of subsystems or a change of subsystems with or without other changes. Alternatively stated, the universality implied by the $[eqn]$ formulation reveals the lack of sensitivity of the SPEE to specific system details; it depends only on their collective information. A potential extension of this formulation for multi-particle systems is also important and thus desirable. The entanglement can be then studied in the traditional sense, where a partition consists of a fraction of the total number of particles. An interesting path to pursue will be to develop a single-parametric formulation of the entanglement dynamics for such cases and thereby seek the universality classes among bipartite and multipartite entanglement. This work is currently in progress.