Superradiance in finite quantum systems randomly coupled to continuum
Pavel Stránský Pavel Cejnar
quant-ph math-ph math.MP






stransky@ipnp.troja.mff.cuni.cz
Institute of Particle and Nuclear Physics, Faculty of Mathematics and Physics, Charles University, V Holešovičkách 2, 18000 Prague, Czech Republic 

cejnar@ipnp.troja.mff.cuni.cz
Institute of Particle and Nuclear Physics, Faculty of Mathematics and Physics, Charles University, V Holešovičkách 2, 18000 Prague, Czech Republic 




We study the effect of superradiance in open quantum systems, i.e., the separation of short- and long-living eigenstates when a certain subspace of states in the Hilbert space acquires an increasing decay width. 
We use several Hamiltonian forms of the initial closed system and generate their coupling to continuum by means of the random matrix theory.
We average the results over a large number of statistical realizations of an effective non-Hermitian Hamiltonian and relate robust features of the superradiance process to the distribution of its exceptional points.
We show that the superradiance effect is enhanced if the initial system is at the point of quantum criticality.



Superradiance in finite quantum systems randomly coupled to continuum
    Pavel Cejnar
    December 30, 2023
=====================================================================





§ INTRODUCTION



Imagine a classical system that randomly decays with a rate γ whenever it passes a certain decay zone in the phase space.
For time Δ t spent in the zone, e^-γΔ t is the survival probability and (1-e^-γΔ t) is the probability for the system to disappear. 
As the phase space is assumed to describe the composite system and not the constituents into which it disintegrates, any phase-space probability density gradually looses its normalization to unity.
The accessibility of the decay zone depends on energy E determined by the system's Hamiltonian.
Minimal and maximal Hamiltonian values E_ min and E_ max within the zone (which is assumed to be compact) determine an energy window in which the system gets unstable. 
Assuming ergodicity of classical motions, all states within the above interval will be characterized by  decay rates Γ(E), which depend on details of dynamics but are proportional to γ for any fixed energy.

An analogous setup in the quantum world has significantly different consequences.
Consider a simplified model of a bound quantum system with d-dimensional Hilbert space H of state vectors.
The quantum decay zone is defined as a certain n-dimensional subspace H_ D⊂ H, which couples uniformly to an infinite environment with continuous spectrum.
Consequently, all states within H_ D have the same decay width γ=ħ/(2τ), where τ stands for the mean lifetime (here we define energy width consistently with relation γτ=1/2ħ instead of a more common definition without factor 1/2).
Since any eigenstate of an independently chosen Hamiltonian has generally a non-vanishing overlap with H_ D, the whole system under the Hamiltonian-induced evolution becomes unstable.
One can apply a simple description based on the Feshbach theory <cit.>, in which the open system is attributed by an effective non-Hermitian Hamiltonian Ĥ^(γ)=Ĥ^(0)-iγP̂_ D, with Ĥ^(0) denoting Hamiltonian of the unperturbed closed system and P̂_ D a projector to H_ D.
The eigenstates |κ^(γ)⟩ of Ĥ^(γ) for any finite value of γ have generally non-zero decay widths Γ^(γ)_κ=ħ/(2τ^(γ)_κ). 
If γ grows infinitesimally from zero, Γ^(γ)_κ of all eigenstates increase proportionally.
However, as the width γ grows to higher values, the set of eigenstates eventually splits into two groups:
The first one contains n states whose widths keep growing as Γ^(γ)_κ≈γ.
In contrast, the widths of remaining (d-n) states in the second group overcome the initially increasing trend and return back to zero, Γ^(γ)_κ→ 0, with γ→∞.
For a sufficiently high value of γ, the original real spectrum 

is reorganized so that it consists of n very wide and (n-d) very narrow states with modified but not too distant real energies.

An example of such behavior is shown in Fig. <ref> below.

The above-described phenomenon was probably for the first time pointed out in Ref. <cit.>.
Its mechanism was soon related <cit.> to so-called Dicke superradiance <cit.>, in which mutual coupling mediated by a common electromagnetic field in a dense ensemble of atoms leads to a collective enhancement and time squeeze of the spontaneous radiation emitted from the ensemble.
In the case of an open system, the role of the mediating field is played by mutual coupling of the unperturbed eigenstates by the decay-inducing part of the full Hamiltonian, and an analog of the superradiant burst is the creation of the group of very short-living (super-radiant) states on a background of long-living (sub-radiant) ones.
To distinguish the latter effect from the original notion of superradiance, we call it non-Hermitian superradiance (NHSR), emphasizing the non-Hermiticity of the model Hamiltonian which captures the coupling of its eigenstates to the continuum. 

The NHSR has notable implications in complex many-body systems, such as nuclei, baryon excitations, atoms, atomic clusters and molecules, open quantum system with gain and loss etc. <cit.>. 
Atomic nuclei, in particular, show neat examples of narrow quasi-stationary states (such as neutron or proton resonances) coexisting with much broader structures (various kinds of doorway states or giant resonances), and the above toy model of resonance trapping provides an elementary background for their description.
It is relevant also in various artificial quantum systems realizable with the aid of recent laboratory quantum simulators <cit.>. 
For detailed reviews with additional references see Refs. <cit.>.

Essential insight into the mechanism of the NHSR follows from the mathematics of so-called exceptional points (EPs) <cit.>. 
These represent degeneracies of a Hamiltonian Ĥ^(λ) with a real discrete spectrum depending on parameter λ in complex-extended domain λ∈ℂ, where Ĥ^(λ) is non-Hermitian and its eigenvalues complex.
Convergence of EPs to a point λ_ c on the real parameter axis with asymptotically increasing size of the system was shown to trigger quantum phase transitions (QPTs) <cit.>.
In the setup of an open quantum system, the EPs are responsible for the NHSR-type redistribution of complex eigenvalues of the Hamiltonian with a running decay rate γ <cit.>.

The role of EPs in the NHSR process was studied in Refs. <cit.>.
In the present paper we extend those studies in several directions:
First, we generate the decay-inducing part of the Hamiltonian by means of the random matrix theory.
Statistical averaging over an ensemble of Hamiltonian realizations yields robust results, washing out any particular correlation between the initial eigenbasis and the set of open states.
Second, we investigate the dependence of the superradiance process on the form of the initial Hamiltonian.
In particular, we connect the two above-mentioned roles of EPs showing that if the initial Hamiltonian is at the quantum critical point, the system exhibits sharper NHSR dynamics. 

The outline of the paper is as follows:
We first describe our statistical model to study the NHSR phenomenon (Sec. <ref>), demonstrate on it a general effect of EPs (Sec. <ref>), and derive some overall properties of the complex spectrum (Sec. <ref>). 
We further analyze special properties of the NHSR arising from criticality of the initial Hamiltonian, in particular from its association with the first- and second-order QPT (Sec. <ref>). 
At last, we summarize the results (Sec. <ref>).




§ MODEL HAMILTONIAN



   

In this section we introduce a simple model of NHSR used in our analysis.
The model works in a finite, d-dimensional Hilbert space H, in which n∈[1,d-1] orthogonal states are supposed to be equally coupled to continuum.
An example of the superradiant separation of short- and long-living states for this model is shown in Fig. <ref>.
Everywhere we set ħ=1.

At first we consider a general Hamiltonian of the form

    Ĥ^({λ_l})=∑_k=1^dE_k^(0)|k⟩⟨k|_Ĥ^(0)+∑_l=1^n(ϵ_l-iγ_l)_λ_l|ϕ_l⟩⟨ϕ_l| ,

where the first term Ĥ^(0) (the unperturbed initial Hamiltonian) describes a closed system with energies E_k^(0)∈ℝ and orthonormal eigenvectors |k⟩, while the second term (decay-inducing Hamiltonian) defines n decaying states given by orthonormal vectors |ϕ_l⟩  with complex energies λ_l∈ℂ, each composed of a real part ϵ_l∈ℝ and a non-negative decay width γ_l∈ℝ.
The orthonormality conditions read as k_1k_2=δ_k_1k_2 and ϕ_l_1ϕ_l_2=δ_l_1l_2, where both vector sets {|k⟩}_k=1^d and {|ϕ_l⟩}_l=1^n are supposed to be incompatible.

To illustrate physical meaning of the schematic Hamiltonian (<ref>), we assume a system of N fermions in a finite Hilbert space generated by d_ sp single-particle states in a mean-field potential well.
Out of these states, (d_ sp-n_ sp)  are stable, bounded inside the well below the continuum threshold energy, while the remaining n_ sp states are quasi-stable, located above the continuum threshold but confined with large lifetimes in the potential well region due to a barrier.
For N fermions, the total Hilbert space of dimension d=d_ spn_ sp is spanned by d_ sp-n_ spN stable and n=d_ spN-d_ sp-n_ spN unstable mean-field configurations.
The unstable subspace is generated by N-body basis vectors {|ϕ_l⟩}_l=1^n in which at least one fermion is in the unstable single-particle state, the stable subspace is generated by the remaining vectors {|ϕ_l⟩}_l=n+1^d. 

Now we assume that all fermions are subject to some residual two-body interactions acting between all d_ sp single-particle states.
The real parts of the unperturbed energies and the residual interaction between individual mean-field states (its matrix elements can be calculated from wave functions restricted to the domain inside the well) jointly form a Hermitian Hamiltonian Ĥ^(0)  which approximately describes the system without tunneling.
The eigenbasis {|k⟩}_k=1^d of Ĥ^(0) is rotated relative to the mean-field basis {|ϕ_l⟩}_l=1^d.
The decay describing part of the Hamiltonian consists of projectors to original unstable states {|ϕ_l⟩}_l=1^n with complex coefficients expressing decay widths γ_l and real-energy corrections ϵ_l depending on the shape of the confining potential barrier.
This results in a non-Hermitian Hamiltonian of the form (<ref>).

In the following, we will use a simplified version of this Hamiltonian, namely

    Ĥ^(λ)=Ĥ^(0)+(ϵ-iγ)_λ∑_l=1^n|ϕ_l⟩⟨ϕ_l|_P̂_ D .

Here, the decay-inducing part of the Hamiltonian λP̂_ D=Ĥ_ D^(λ) is determined by a single projection operator P̂_ D and all the decaying states, which form a subspace H_ D=P̂_ D H⊂ H, have the same complex energy ϵ-iγ=λ.
The imaginary component γ expresses the decay width of states in H_ D and the real component ϵ represents their energy shift. 
We stress that though the simplified Hamiltonian (<ref>) restricts direct applicability of our model to specific systems, it captures essential features of the more complex Hamiltonian (<ref>) while considerably reducing the number of free parameters.  

The decay-inducing term of Hamiltonian (<ref>) can be cast in the form Ĥ_ D^(λ)=λP̂_ D+0P̂_ (with P̂_=𝕀̂-P̂_ D) with two eigenvalues λ and 0.
The eigenspace associated with eigenvalue λ coincides with the decaying subspace H_ D of dimension n, while the eigenspace with eigenvalue 0 is the subspace H_=P̂_ H with dimension n_=d-n. 
Note that the cases n=0 with Ĥ_ D^(λ)=0 and n=d with Ĥ_ D^(λ)=λ𝕀̂ are both trivial and we exclude them.

In our analysis, the unperturbed component Ĥ^(0) of the total Hamiltonian (<ref>) is associated with the Lipkin-Meshkov-Glick model <cit.>. 
It is build from quasispin operators (Ĵ_1,Ĵ_2,Ĵ_3) satisfying the standard angular-momentum commutation relations.
The conserved quantity Ĵ^2 is fixed at the value j(j+1), so that the Hilbert space H spanned by eigenvectors of Ĵ_3 has dimension d=2j+1.
To generate the spectrum of energies E_k^(0), we consider the following three alternative forms:

    Ĥ^(0)   =   Ĥ^(0)_ HO =s_0(j) Ĵ_3+a_0(j)
     ,
       =   Ĥ^(0)_ PT1=s_1(j)(Ĵ_3-3/jĴ_1^2)+a_1(j)
     ,
       =   Ĥ^(0)_ PT2=s_2(j)(Ĵ_3-1/2jĴ_1^2)+a_2(j)
     ,

where s_i(j) and a_i(j) with i=0,1,2 are dimension-dependent scaling and shift constants ensuring invariant bounds E_1^(0)=0 and E_d^(0)=d of the unperturbed spectrum and therefore a unit average spacing between neighboring energy levels in all cases (<ref>)–(<ref>).
We introduce cumulants expressing global properties of the unperturbed spectrum, in particular the average E^(0)=1/d∑_k=1^dE_k^(0), which represents a central energy of the spectrum, the variance Δ^2 E^(0)=1/d∑_k=1^d[E_k^(0)-E^(0)]^2, which characterizes a squared spread of the spectrum, and the skewness Δ^3 E^(0)=1/d∑_k=1^d[E_k^(0)-E^(0)]^3, which quantifies an asymmetry of the spectrum with respect to E^(0).
We note that all results below are independent of an overall shift of the spectrum, while a change of the average spacing between levels (spread of the spectrum) can be compensated by an appropriate rescaling of parameter λ.  

   

All Hamiltonians (<ref>)–(<ref>) can be written, using the Holstein-Primakoff transformation <cit.>, in terms of a single pair (q̂,p̂) of conjugate coordinate and momentum operators satisfying in the classical limit a constraint q^2+p^2≤ 2 (for a general discussion of this mapping see, e.g., <cit.>).
The Hamiltonian Ĥ^(0)_ HO (the abbreviation standing for harmonic oscillator) represents an equidistant spectrum composed of d levels with spacing E^(0)_k+1-E^(0)_k=1 associated with a potential V(q)∝ q^2+ const.
Hamiltonians Ĥ^(0)_ PT1 and Ĥ^(0)_ PT2 (phase transitional of the first and second kind) describe representative d-dimensional spectra of two types of quantum critical systems:
Ĥ^(0)_ PT1 stands for the first-order QPT Hamiltonian with energy spectrum consisting of parity doublets of levels in a degenerate symmetric double-well potential V(q)∝ 3q^4-5q^2+ const. 
The spacing between the low-lying states forming the doublets decreases with dimension as E^(0)_k+1-E^(0)_k∝ e^-c_kd (where c_k is a constant), while spacing between states outside the doubles remains of order E^(0)_k+2-E^(0)_k+1∼ O(1).
Ĥ^(0)_ PT2 represents a second-order QPT Hamiltonian associated with a quartic potential V(q)∝ q^4+ const. 
Its spectrum exhibits a typical power-law cumulation of low-energy levels according to E^(0)_k+1-E^(0)_k∝(k/d)^1/3.
Energy levels of all the three initial Hamiltonians are for a moderate dimension shown in Fig. <ref>. 
More details concerning the critical Hamiltonians can be found in Ref. <cit.>. 

The decay-inducing component of Hamiltonian (<ref>) is represented not by a single fixed operator Ĥ_ D^(λ), but by a suitable statistical ensemble of its possible (in some sense equivalent) realizations. 
The results are obtained by averaging over a large number N_ R (of order from 10^1 to 10^4) of samples from this ensemble.
The orthonormal vectors {|ϕ_l⟩}_l=1^n and {|ϕ_ l'⟩}_l'=1^n_ forming in each realization the bases of subspaces H_ D and H_, respectively, result from a random orthogonal transformation of the original eigenbasis {|k⟩}_k=1^d.

To achieve a completely unbiased (isotropic) generation of these bases, we use the eigenstate components of matrices from the Gaussian orthogonal ensemble (GOE), which has the equivalence of bases in its definition <cit.>.
In particular, for each realization we perform the following steps: (i) we generate a random d-dimensional real symmetric matrix Ĥ_ GOE∈ GOE, i.e., a matrix with independent elements taken from zero-mean Gaussian distributions with σ^2=2 or 1 for diagonal or off-diagonal elements, respectively, (ii) we find an orthogonal matrix Ô transforming Ĥ_ GOE to the diagonal form D̂=Ô^ TĤ_ GOEÔ, where Ô^ T=Ô^†=Ô^-1  stands for the transpose of Ô, (iii) we randomly chose n columns of Ô and associate them with the decaying states {|ϕ_l⟩}_l=1^n.
This procedure is repeated until a large number N_ R of Ĥ_ GOE realizations yields a robust estimate of the analyzed quantities.
The averaging over the random-matrix ensemble will be further denoted by angular brackets ∙, which is in contrast to the bar symbol ∙ introduced above for the non-statistical averages over the energy levels.
Note that orthogonality of the transformation Ô, which is inherent in the foundations of GOE, guarantees that matrix elements k_1P̂_ Dk_2=∑_l=1^nk_1ϕ_lϕ_lk_2=∑_l=1^nO_k_1lO_k_2l are real, so that we can write: Imk_1Ĥ(λ)k_2∝ Imλ=-γ.




§ BIFURCATION OF DECAY WIDTHS



   

   

In the following, we focus on the dependencies of the eigensolutions of the Hamiltonian (<ref>) on parameter λ=ϵ-iγ. 
The Hamiltonian is a complex matrix symmetric under transposition (but not under the full Hermitian conjugation) with d generally complex eigenvalues E_κ^(λ)=E_κ^(λ)-iΓ_κ^(λ). 
The real parts E_κ^(λ) stand for energies and the imaginary parts Γ_κ^(λ)≥ 0 represent decay widths of individual eigenstates enumerated by integer κ=1,2,…,d.
Because of non-Hermiticity, the eigenstates of Ĥ^(λ) must be distinguished to right eigenstates satisfying Ĥ^(λ)|κ_ R^(λ)⟩= E_κ^(λ)|κ_ R^(λ)⟩, and left eigenstates satisfying  ⟨κ_ L^(λ)|Ĥ^(λ)=⟨κ_ L^(λ)| E_κ^(λ).
From the transposition symmetry of Ĥ^(λ) it follows that the left eigenvector is just the matrix transpose of the right one.
It can be shown that both types of eigenvectors satisfy the bi-orthonormality condition κ_1_ L^(λ)κ_2_ R^(λ)=δ_κ_1κ_2.

We will assume that ϵ= Reλ, is set constant, while γ=- Imλ is varied within the domain γ≥0.
We can obviously write Ĥ^(ϵ-iγ)=Ĥ^(ϵ-i0)-iγP̂_ D, where we explicitly introduce a shifted Hermitian Hamiltonian Ĥ^(ϵ-i0)=Ĥ^(0)+ϵP̂_ D and an anti-Hermitian decay-inducing Hamiltonian Ĥ_ D^(0-iγ)=-iγP̂_ D.
Evolutions of decay widths for several model settings, averaged over a sample of GOE realizations N_ R, are depicted in Fig. <ref>.
We see that the instability expressed by positive decay widths Γ_κ^(λ) characterizes all eigenstates of Ĥ^(λ) as soon as γ deviates from zero to infinitesimally small positive values.
Indeed, all eigenstates are expected to have non-zero overlaps with a randomly generated decaying subspace H_ D.
Using the Hellmann-Feynman formula, which in the non-Hermitian context reads d/dλ E_κ^(λ)=κ_ L^(λ)d/dλĤ^(λ)κ_ R^(λ), we derive the following relations: ∂/∂γE_κ^(λ)=- Imκ_ L^(λ)P̂_ Dκ_ R^(λ) and ∂/∂γΓ_κ^(λ)= Reκ_ L^(λ)P̂_ Dκ_ R^(λ).
These can be used for an estimation of E_κ^(λ) and Γ_κ^(λ) for small values of γ, where eigenvectors |κ_ R^(λ)⟩ and |κ_ L^(λ)⟩ roughly coincide with eigenvectors |κ^(ϵ-i0)⟩ of the Hermitian Hamiltonian Ĥ^(ϵ-i0).
We get 

    E_κ^(ϵ-iγ)   ≈   E_κ^(ϵ-i0)+O(γ^2)
     ,
    Γ_κ^(ϵ-iγ)   ≈   γ P_κ+O(γ^2)
     ,

where P_κ=∑_l=1^n|ϕ_lκ^(ϵ-i0)|^2 is the probability of identifying the eigenstate  |κ^(ϵ-i0)⟩ with any of the decaying states |ϕ_l⟩.

On the other hand, it is clear that for very large values of γ>0 the decay-inducing part of the Hamiltonian dominates, so that the eigenvectors of the full Hamiltonian (<ref>) become approximately those of the second term. 
Hence for asymptotic γ only n<d eigenstates (those roughly coinciding with the set {|ϕ_l⟩}_l=1^n∈ H_ D) will yield non-negligible widths Γ_κ^(λ)≈γ, while the rest of n_=d-n eigenstates (those approximated by {|ϕ_ l'⟩}_l'=1^n_∈ H_) will have Γ_κ^(λ)≈0.
Indeed, the Hellmann-Feynman formula gives ∂/∂γE_κ^(λ)≈0 for all eigenstates, and ∂/∂γΓ_κ^(λ)≈1 or 0, depending on whether the approximate eigenvectors belong to H_ D or H_, respectively.
A more systematic treatment makes use of a transformed Hamiltonian iγ^-1Ĥ^(ϵ-iγ)=P̂_ D+iγ^-1Ĥ^(ϵ-i0) and a perturbative expansion of its complex eigenvalues in small parameter γ^-1.
This yields:

    E_κ^(ϵ-iγ)   ≈   E_κ^(ϵ-i∞)+O(γ^-2)
     ,
    Γ_κ^(ϵ-iγ)   ≈   {[     γ+c_κγ^-1+O(γ^-3) for |κ^(ϵ-i∞)⟩∈ H_ D,;       c_κγ^-1+O(γ^-3)   for |κ^(ϵ-i∞)⟩∈ H_, ].

where c_κ are some coefficients.
The log-log dependencies in Fig. <ref> support these conclusions.

Having understood the evolution of the spectrum for very small and very large values of γ, we ask what drives the transition between these limiting regimes.
The region in parameter γ where the bifurcation of decay widths takes place is scrutinized under a specific magnifying glass  in Fig. <ref>.
The vertical axis of this figure shows the log-log slope χ_κ^(λ)=d/dlogγlogΓ_κ^(λ)=(γ/Γ_κ^(λ))d/dγΓ_κ^(λ) of individual decay widths from the ϵ=0 panel of Fig. <ref>.
The value χ_κ^(λ)=+1 indicates the Γ_κ^(λ)∝γ behavior, while χ_κ^(λ)=-1 implies Γ_κ^(λ)∝γ^-1.
We observe that the transitional region between these values is a relatively narrow interval of γ, for the selected values of d and n roughly given by 1≲γ≲ 20.
In this region, the slopes in Fig. <ref> show large fluctuations, indicating a kind of turbulent evolution of individual level widths.
It turns out that the exceptional points play a crucial role in this evolution.

The EP of a parameter-dependent Hamiltonian Ĥ^(λ) is a point λ_ EP of the complex parameter space where two (or more) complex eigenvalues E_κ^(λ) and E_κ'^(λ) become degenerate.
In a typical case, the complex degeneracy has a different character than an ordinary degeneracy of a Hermitian Hamiltonian.
First, Ĥ^(λ) at λ=λ_ EP fails to provide a complete system of left and right eigenvectors, but the single eigenvector associated with both degenerate levels becomes self-orthogonal: κ_ L^(λ_ EP)κ_ R^(λ_ EP)=0.
Second, the EP is not just a conical intersection of the two energy surfaces but a complex square-root singularity.
The path in λ∈ℂ satisfying the condition Im  E_κ^(λ)= Im  E_κ'^(λ) is terminated at  λ_ EP, where it smoothly connects to herefrom issuing path defined by Re  E_κ^(λ)= Re  E_κ'^(λ).
On the respective sides of the EP, the real and imaginary parts of E_κ^(λ)- E_κ'^(λ)  bifurcate according to the square-root dependence, which—with regard to the imaginary part—hints at the NHSR behavior.
A passage in a vicinity of an EP induces an avoided crossing of complex eigenvalues (which may include actual crossings of real energies or decay widths) for a single pair of levels along with fast modifications of the corresponding eigenvectors.
The closer is the EP to the selected trajectory of parameter λ, the sharper are the associated changes.
Multiple EPs induce a turbulent flow of eigenvalues with many avoided crossings and ongoing structural redistributions of eigenvectors, while an absence of EPs near the trajectory implies a laminar flow of eigenvalues and virtual freeze of eigenvectors.

Features of EPs are usually first elucidated on 2× 2 matrices.
We sketch such a trivial d=2 version of our model in Appendix <ref> (cf. Ref. <cit.>).
It is shown that setting ϵ= Reλ_ EP, so that the path with increasing γ crosses the EP of the Hamiltonian at γ=| Imλ_ EP|, we observe a sudden bifurcation of the widths Γ_1^(λ) and Γ_2^(λ) from a common value equal to 1/2γ to distinct—decreasing and increasing—components.
This is a non-analytic realization of the NHSR in the d=2 system.
If we increase the distance of ϵ from Reλ_ EP, the decay widths of both levels differ already at small γ so that the ratio between the widths of short- and long-living state increases with |ϵ- Reλ_ EP|.
At the beginning, for γ not exceeding | Imλ_ EP|, both widths grow proportionally to γ. 
However, for larger γ, somewhere at γ≳| Imλ_ EP|, the smaller of both widths turns back and starts decreasing as γ^-1. 
This is a smooth but still clear realization of the NHSR effect.
Examples of these scenarios, depicted in Fig. <ref> below, can be compared to less trivial but rather analogous behavior for d>2 shown in Fig. <ref>.

   

Focusing on these d>2 cases, we first note that a symmetric matrix of dimension d depending linearly on parameter λ∈ℂ has in general 1/2d(d-1) pairs of complex-conjugate EPs.
However, in case of Hamiltonian (<ref>) a large part of EP pairs  migrates to infinity, while at finite values of λ we observe a reduced number of pairs:

    12N_ EP=n(d-n)
     .
 
This is ≲ 50 % of the full number of pairs (the maximal fraction is reached for n=1/2d).
Figure <ref> shows distributions of these EPs in λ∈ℂ for the three initial Hamiltonians (<ref>)–(<ref>) and for various sizes of the decaying subspace H_ D.
The distribution is averaged over a large number of GOE realizations of H_ D.

The Hilbert space dimension in Fig. <ref>  is rather moderate, d=16, but the distributions of EPs remain qualitatively similar for larger dimensions, except an overall scaling of the absolute value |λ|.  
The essential fraction of the EP distribution is expected to be located within a domain defined by

    |λ|≲S d^2/√(n(d-n))∝{[                d^3/2 for  n≪d   or  n_≪d,;                    d       for  n≈n_≈d/2, ].

where S=(Δ^2 E^(0)/d^2)^1/2 is the linear spread of the unperturbed spectrum divided by the dimension (e.g., for H^(0)_ HO and H^(0)_ PT2 with unit average spacings we get S≈0.3).
The scaling formula (<ref>) follows from the analysis of Ref. <cit.> (cf. <cit.>) which showed that a maximum of the EP distribution for a general perturbed Hamiltonian Ĥ^(0)+λV̂ is achieved when spectral variances (spreads of the spectra) of both terms H^(0) and λV̂ are approximately equal to each other.
The formula is also confirmed by our numerical simulations.

   

Let us proceed to the discussion of the impact of the EP distribution on the NHSR dynamics.
We saw in Fig. <ref> (the lower right panel) that two distinct groups of states with shorter and longer lifetimes can exist already long before the superradiant transition takes place, i.e., below the γ value where the longer-living group turns back to the Γ→0 path.
So in the investigation of the NHSR we have to distinguish a static bi-modality (existence of two peaks in the distribution of decay widths which both move with increasing γ to the increasing-Γ direction) from a dynamic bi-modality (existence of two peaks moving to opposite Γ→0 and Γ→∞ directions).
The EPs play an important—though not fully deterministic—role in the description of both these features.
The static bimodality is achieved if the value of ϵ either exceeds the upper peripheral region of the EP distribution projected to the Reλ axis, or undershoots its lower peripheral region.
Analogously, the dynamic bimodality occurs when the value of γ passes the upper periphery of the EP distribution projected of  to the Imλ axis.
However, it needs to be stressed that details of the formation of both kinds of bi-modalities for λ located still within the principal part of the EP distribution are certainly beyond the above simplified description.

   

To support the claims from the previous paragraph, we show in Fig. <ref> several transient distributions of decay widths for the same model setup as used in the middle row of Fig. <ref> (d=2n=16).
Each panel depicts a distribution of widths averaged over a large number of GOE realizations for a fixed parameter point λ=ϵ-iγ together with arrows, which indicate the direction in which a given part of the distribution evolves with increasing γ. 
Note that some of the distributions represent γ=const. cuts of Figs. <ref> and <ref>.
The fixed parameter points in Fig. <ref> are organized so that panel (a) is associated with λ close to zero, while the pair of panels (b),(c) and further the triples (d),(e),(f) and (g),(h),(i) correspond to increasing distance from the origin in various directions of the λ∈ℂ plane.
We see that the panels within the same group show a comparable degree of splitting of the spectrum to short- and long-living states, but the bi-modality has the static character if γ is not large enough.
The dynamic bi-modality is triggered only when γ passes the upper periphery of the projected EP distribution.
 
The distribution of EPs in Fig. <ref> for n=1/2d, which is the case corresponding to Fig. <ref>, is approximately symmetric under rotations around λ=0.
However, the distributions for n=1 and n=d-1 with H^(0)=H^(0)_ PT1 and H^(0)_ PT2 exhibit strong asymmetries under rotations and even under the mirror transformation ϵ↔-ϵ.
In these cases, the NHSR dynamics sensitively depends on the sign of ϵ. 
This is illustrated in Fig. <ref>, where we show the decay width distributions at two points with Reλ=±ϵ and Imλ small for H^(0)_ PT2. 
Both points belong to the regions within the principal parts of the EP distributions in the rightmost column of Fig. <ref>.
We see that in the n=1 and n=d-1 cases the decay width distributions have significantly different forms for both ±ϵ values, in agreement with the observed asymmetries of the respective EP distributions.
In both these cases, the width distribution corresponding to the ϵ value closer to the main maximum of the EP distribution shows a more developed static bi-modality.
However, the  width distributions for ±ϵ differ also in the n=1/2d case, for which the EP distribution is roughly symmetric.
This indicates that initial stages of the static bi-modality formation transcend the description based on the overall EP distribution, which does non reflect links of individual EPs to specific pairs of levels (cf. Ref. <cit.>).




§ GLOBAL PROPERTIES OF THE COMPLEX SPECTRUM



Employing the trace-based method described e.g. in Ref. <cit.> (see also <cit.>), one can derive expressions for cumulants of the whole eigenvalue spectrum of a general (parameter-dependent) Hamiltonian Ĥ^(λ).
In particular, the average of all eigenvalues is determined from

    E^(λ)≡1/d∑_κ=1^d E_κ^(λ)=1/d  TrĤ^(λ) ,

and the eigenvalue variance is given by 

    Δ^2 E^(λ)   ≡   1/d∑_κ=1^d[ E_κ^(λ)- E^(λ)]^2=
    
       =   1/d  TrĤ^(λ)^2-1/d^2 Tr^2Ĥ^(λ) .

For the non-Hermitian Hamiltonian (<ref>), the averages and variances related to complex energies E_κ^(λ) shall be evaluated separately for real and imaginary parts.
From Eq. (<ref>) we get easy formulas for both types of averages,

    E^(ϵ-iγ)-E^(0)   =   n/dϵ ,
    
    Γ^(ϵ-iγ)   =   n/dγ ,

and from Eq. (<ref>) we obtain the following slightly more complicated relations for the variances:

    Δ^2E^(ϵ-i0)-Δ^2E^(0)=2ϵ n/d A+ϵ^2 n/d(1-n/d)
    ,

where A=1/n∑_k=1^d(E_k^(0)-E^(0))kP̂_ Dk, and

    Δ^2E^(ϵ-iγ)-Δ^2E^(ϵ-i0)=Δ^2Γ^(ϵ-iγ)-γ^2 n/d(1-n/d)
    .

Formulas (<ref>) and (<ref>) imply linear dependences of the average real energy and the average decay width on ϵ and γ, respectively. 
Formula (<ref>) expresses a quadratic dependence of the real energy variance on ϵ for γ=0, while formula (<ref>) captures a specific relation between the variances of real energies and decay widths for variable γ and constant ϵ.

The linear term of the quadratic dependence of the real energy variance in Eq. (<ref>) is proportional to the coefficient A, which quantifies an asymmetry of the unperturbed energy spectrum with respect to the decaying subspace. 
It is evaluated as an average of the energy displacement (E_k^(0)-E^(0)) calculated with normalized weight factors w_k=1/n∑_l=1^n|ϕ_lk|^2 proportional to the overlap probability of the respective unperturbed eigenstate with H_ D.
This asymmetry determines the point ϵ_ min=-Ad/(d-n) of a minimal quadratic spread Δ^2E^(ϵ-i0) of the real spectrum along the real λ axis, i.e., the point of maximal compression of the spectrum.
This point shall roughly correspond to the projection of the center of mass of the EP distribution to the real λ axis.
The coefficient A depends on the statistical realization of the decaying subspace, as well as on the energy spectrum of the initial Hamiltonian.
The GOE average is obviously A=0, but an analytic evaluation of the higher cumulants is hindered by non-trivial correlations of weight factors w_k for various k.
Nevertheless, we can infer that Δ^2A and Δ^3A are correlated with the respective cumulants of the unperturbed energy spectrum Δ^2E^(0) and Δ^3E^(0).
This qualitatively explains the above-mentioned ϵ↔-ϵ asymmetries of EP distributions in Fig. <ref> for H^(0)=H^(0)_ PT1 and H^(0)_ PT2, whose spectra are apparently skewed (cf. Fig. <ref>).

   

We now turn to the consequences of Eq. (<ref>) for the NHSR dynamics.
The left-hand side of this formula represents the change of variance of the real spectrum between γ=0 and γ>0 points with the same ϵ.
The right-hand side is the difference between the actual variance of the decay widths at λ=ϵ-iγ and the variance of a two-component ensemble consisting of n widths equal to γ and n_=d-n widths equal to 0. 

According to Ref. <cit.>, the spectrum of Hamiltonian (<ref>) satisfies the following inequalities:

    E_1^(ϵ-i0)≤ Re  E_κ^(ϵ-iγ)   =   E_κ^(ϵ-iγ)≤ E_d^(ϵ-i0) 
    
    0≤- Im  E_κ^(ϵ-iγ)   =   Γ_κ^(ϵ-iγ)≤γ .

The bounds E_1^(ϵ-i0) and E_d^(ϵ-i0) in (<ref>) are the lowest and the highest eigenvalues of Ĥ^(ϵ-i0), and similarly 0 and γ in (<ref>) are the lowest and highest eigenvalues of  γP̂_ D.
Inequality (<ref>) implies that the right-hand side of Eq. (<ref>) is semi-negative.
The limiting zero value corresponds to γ=0, while for increasing γ we expect negative values that remain of order ∼ O(1) even in the asymptotic case. 
Hence we see that the increase of the decay rate γ always reduces the spread of the real energy spectrum.
This reduction is correlated with the increase of the spread of decay widths of individual eigenstates and survives the γ→∞ limit.
As shown in Fig. <ref>, the maximal contraction of the real energy spectrum, reaching the asymptotic value of 50 % of the squared spread of the original spectrum (≈ 71 % of the linear spread), is obtained for n=1/2d.




§ EFFECTS OF CRITICALITY



   




As seen in Fig. <ref>, the distributions of decay widths for various choices (<ref>), (<ref>) and (<ref>) of the initial Hamiltonian Ĥ^(0) are rather close to each other for moderate dimensions d.
However, the differences get much more significant as d grows.
In Fig. <ref> we show an analog of  Fig. <ref> for d=1024.
In comparison to the previous figure, the parameter coordinates [ϵ,γ] corresponding to panels (a)–(i) were scaled roughly according to Eq. (<ref>).
We immediately see that the choice of Ĥ^(0) plays a substantial role.
In general, the first- and second-order QPT Hamiltonians Ĥ^(0)_ PT1 and Ĥ^(0)_ PT2 result at various fixed values of λ in much more developed bi-modal structures of the width distribution than the harmonic oscillator Hamiltonian Ĥ^(0)_ HO.
Interestingly, comparing pairs of panels (c)-(c'), (f)-(f') and (i)-(i') in Fig. <ref>, we observe clear differences between decay width distributions for both critical Hamiltonians at mirror-conjugate values of Reλ.
These again reflect the -ϵ↔+ϵ asymmetries discussed in connection with Fig. <ref>(b) and Eq. (<ref>).




The differences of decay width distributions in Fig. <ref> may indicate that for the critical Hamiltonians the NHSR dynamics is more advanced than for the non-critical one.
However, as shown above, the NHSR is not defined by a static separation of short- and long-living states at a given γ, but more substantially relies on the reversed evolution of both modes as γ grows to large values and on the convergence of the long-living mode to Γ=0.
Whether the criticality of Ĥ^(0) implies some enhancement of the latter process is studied in Fig. <ref>.

This figure compares the smallest decay widths Γ_ min^(λ)( HO) and Γ_ min^(λ)( PT2) obtained by minimization over the entire set of all d levels for the harmonic oscillator Hamiltonian Ĥ^(0)_ HO and the second-order QPT Hamiltonian Ĥ^(0)_ PT2.
The imaginary part of λ is taken very large, γ=16d, so that we are far above the superradiant transition to decreasing width regime, the real part of λ is set to ϵ=0. 
Results are displayed for several Hilbert space dimensions ranging from d=4 to d=8192, the dimension of the decaying subspace being always n=1/2d.
The minimal decay widths and their ratios were averaged over several GOE realizations, yielding the dots in Fig. <ref>.
The error bars indicate standard deviations of the distributions obtained.
Note that the numbers of realizations N_ R decrease with d for computability reasons, but they are sufficient for the purposes of the present study.
The numerical precision, however, prevents us to do the same calculation for Ĥ^(0)_ PT1, for which the parity doublets of levels become too close with increasing d. 

We see in the inset of Fig. <ref> that for d≳2^10 both GOE averages of the minimal decay width scale roughly as Γ_ min^(λ)∝d.
This is consistent with the choice γ=16d since, with respect to the scaling formula (<ref>), the pattern in  Fig. <ref> remains roughly invariant under the change of dimension if variables γ and Γ are replaced by γ/d and Γ/d, respectively.

The most important information follows from the ratio between both widths for large d.
A plot of the GOE average Γ_ min^(λ)( PT2)/Γ_ min^(λ)( HO) at λ=0-i(16d) is shown in the main panel of Fig. <ref>.
It suggests that the asymptotic value of the average ratio is roughly equal to 0.55.
Since in the large-γ domain the evolution of widths is described by Eq. (<ref>), hence Γ_ min^(λ)≈ c_ minγ^-1, the above numerical value is assigned also to the average ratio c_ min( PT2)/c_ min( HO) of the respective coefficients.
It is clear that the particular value 0.55 would not apply to other choices of ϵ (cf. Fig. <ref>), but the NHSR process is generally more advanced for the second-order critical system than for the non-critical harmonic oscillator.

Let us stress that all initial Hamiltonians in this study have the same spread of spectra (average spacing of levels), so the discussed effect reflects some tinier differences in the distributions of levels.
Indeed, Hamiltonian Ĥ^(0)_ PT2 with the quartic potential exhibits a strong (increasing with d) accumulation of levels near the lowest energy, so a higher occurrence of small spacings between levels.
Therefore, some EPs appear closer to λ=0 and some of the levels exhibit the crossover to the Γ_κ^(λ)∝γ^-1 regime sooner than in the case of the harmonic oscillator. 
This explains the observed acceleration of the NHSR process for the second-order QPT Hamiltonian. 
There is no doubt that a similar (even stronger) effect would be observed for the first-order critical Hamiltonian Ĥ^(0)_ PT1, if it is made available for large-d numerical simulations.




§ SUMMARY



We summarize our main results in the following items:

(i) The non-Hermitian superradiance (i.e., splitting of energy eigenstates of an open quantum system with increasing coupling of a certain subset of states to continuum into the groups of short- and long-living states with decay widths Γ→∞ and Γ→0) is a universal effect, qualitatively independent of the form of the initial Hamiltonian and the choice of the decaying subspace.
The effect is in a rudimentary form present already in a two-level system (Fig. <ref>), which contains the essence of more complicated dependencies for larger dimensions. 
 
(ii) The initial (small-coupling) and asymptotic (large-coupling) stages of the NHSR evolution are understood from elementary perturbative expressions. 
In log-log plots (like Fig. <ref>) we observe laminar flows of real energies and decay widths described by the respective pair of equations (<ref>)–(<ref>) or (<ref>)–(<ref>).

(iii) The intermediate stages of the NHSR evolution are strongly influenced by the distribution of EPs of an effective non-Hermitian Hamiltonian.
At these stages, the real energies and decay widths evolve in a turbulent way, showing numerous collisions (actual or avoided crossings) and large fluctuations of log-log slopes (see Fig. <ref>).

(iv) Essential information follows from projections of the EP distribution (cf. Fig. <ref>) to real and imaginary axes of the coupling parameter λ=ϵ-iγ.
Considering trajectories in the λ plane with fixed ϵ, we distinguish two basic ESQPT scenarios (see Fig. <ref>):
First, if ϵ is within the support of the projected EP distribution on the real axis, the decay widths of individual levels for small γ form a single-mode distribution. 
As γ increases, the distribution moves to increasing values of Γ and eventually splits into the increasing and decreasing branches when γ gets outside the support of the projected EP distribution on the imaginary axis.
Second, if ϵ is outside the support of the real projection of the EP distribution, the decay-width distribution is bi-modal already at small values of γ.
However, the splitting into the increasing and decreasing branches again happens only when γ gets out of the imaginary projection of the EP distribution.
These conclusions are supported by examples of Γ distributions in Figs. <ref>, <ref> and <ref>, but must be taken only as raw simplifications.  

(v) Squared spreads of the spectra of real energies and decay widths are mutually correlated [see Eq. (<ref>)].
The opening of the system always leads to narrowing of the real energy spectrum (see Fig. <ref>).

(vi) Details of the NHSR dynamics depend on the structure of real energy spectrum of the γ=0 Hamiltonian.
First, the formation of bi-modal structures in the decay width distribution depends on the value of ϵ in a way that goes beyond the overall EP distribution (see examples in  Fig. <ref>). 
Second, even the asymptotic-γ behavior carries traces of the γ=0 energy spectrum.
In particular, the presence of small spacings between energy levels is equivalent to a closer approach of some of the EPs to the real λ axis, which shifts the crossover of some states to the Γ→0 stage of evolution to smaller values of parameter γ and makes the stabilization process of these states more advanced (see Fig. <ref>).
In this sense one can observe a speedup of the NHSR effect.
This applies especially to systems at quantum critical points, which contain closely located levels by definition.  




§ ACKNOWLEDGMENTS


We acknowledge funding by the Charles University Project UNCE/SCI/013.




§ TWO-DIMENSIONAL MODEL





Hamiltonian (<ref>) in the d=2 case reads

    Ĥ^(λ)=(E^(0)_1   0
    0   E^(0)_2)_Ĥ^(0)+(ϵ-iγ)_λ(cos^2ϑ   cosϑsinϑ
    sinϑcosϑ   sin^2ϑ)_(cosϑ
    sinϑ)(cosϑ sinϑ) ,

where ϑ (assumingly ≠ 0,π since these special values would yield trivial solutions) is an angle determining the single decaying state |ϕ_1⟩. 
Note that for d=2 only the case of n=1 is non-trivial.
Two complex eigenvalues

    E_1,2^(λ)≡ E_1,2^(λ)-iΓ_1,2^(λ)=E^(0)_1+E^(0)_2+λ/2    
    ±√((E^(0)_1-E^(0)_2/2)^2+(λ/2)^2+λ E^(0)_1-E^(0)_2/2cos 2ϑ)
 
become degenerate for λ equal  to 

    λ_ EP±=-(E^(0)_1-E^(0)_2) e^± 2 iϑ .

These are complex conjugate EPs of the matrix (<ref>), which both represent square-root singularities of the function E_1^(λ)- E_2^(λ)=√((λ-λ_ EP+)(λ-λ_ EP-)).
At λ=λ_ EP±, the Hamiltonian can undergo a similarity transformation to the Jordan form

    Ĥ^(λ_ EP±)→( E_ EP±   1
    0    E_ EP±)
     ,

where E_ EP±= E_1^(λ_ EP±)= E_2^(λ_ EP±).
This means that two right eigenvectors |1_ R^(λ)⟩ and |2_ R^(λ)⟩ associated with the eigenvalues E_1^(λ) and E_2^(λ) for λ≠λ_ EP± contract to a single one |1_ R^(λ_ EP±)⟩ at the EP, and the same happens to left eigenvectors.
The unique eigenvector at the EP is self-orthogonal in the sense 1_ L^(λ_ EP±)1_ R^(λ_ EP±)=0.

Setting the real part ϵ of parameter λ to Reλ_ EP±=(E^(0)_2-E^(0)_1)cos 2ϑ and varying the imaginary part along a half-line γ∈[0,∞), we proceed in opening the system along a trajectory passing the EP located in the lower complex half-plane of λ.
It follows from Eq. (<ref>) that for γ increasing from zero towards the EP absolute value, the widths Γ_1^(λ) and Γ_2^(λ) are both equal to 1/2γ, whereas the real energies E_1^(λ) and E_2^(λ) differ and collapse, according to the square-root formula, to the EP degeneracy  as γ→| Imλ_ EP±|=|(E^(0)_2-E^(0)_1)sin 2ϑ|.
For γ increasing further from  | Imλ_ EP±|, the real energies stay degenerate whereas the decay widths  exhibit the square-root bifurcation to a short- and long-living asymptotic forms Γ_ short^(λ)=γ and Γ_ long^(λ)∝γ^-1.
This critical (non-analytic) form of the superradiant scenario is illustrated in Fig. <ref> for model settings in which λ_ EP±=± i.

For an increasing distance  of ϵ from Reλ_ EP±, the sharpness of the NHSR transition decreases, but its principal features remain preserved.
These cases are also illustrated in Fig. <ref>.
We see that shorter- and longer-living states exist already at small values of γ where their widths increase linearly.
The ratio between the larger and smaller width grows with the distance from Reλ_ EP±. 
As γ reaches values near or above γ≈| Imλ_ EP±|, one of the widths makes a crossover to the ∝γ^-1 behavior.
At the same time, evolution of real energies with γ increasing across the EP region for ϵ≠ Reλ_ EP± shows a smooth reduction of the spacing |E_1^(λ)-E_2^(λ)|.

99
Fes H. Feshbach, Ann. Phys. 5, 537 (1958); 19, 287 (1962). 
Moi11 N. Moiseyev, Non-Hermitian Quantum Mechanics (Cambridge University Press, Cambridge UK, 2011).
Kle85 P. Kleinwächter and I. Rotter, Phys. Rev. C 32, 1742 (1985).
Sok88 V.V. Sokolov and V. Zelevinsky, Phys. Lett. B 202, 10 (1988).
Sok92 V.V. Sokolov and V. Zelevinsky, Ann. Phys. 216, 323 (1992).
Dic54 R.H. Dicke, Phys. Rev. 93, 99 (1954).
Rot91 I. Rotter, Rep. Prog. Phys. 54, 635 (1991).
Vol04 A. Volya and V. Zelevinsky, AIP Conf. Proc. 777 (2005), ed. V. Zelevinsky, p. 229.
Aue11 N. Auerbach and V. Zelevinsky, Rep. Prog. Phys. 74, 106301 (2011).
Cel11 G.L. Celardo, N. Auerbach, F.M. Izrailev, and V.G. Zelevinsky,  Phys. Rev. Lett. 106, 042501 (2011).
Liu14 C. Liu, A. Di Falco, and A. Fratalocchi, Phys. Rev. X 4, 021048 (2014).
Ele14 H. Eleuch and I. Rotter, Eur. Phys. J. D 68, 74 (2014).
Rot15 I. Rotter and J.P. Bird, Rep. Prog. Phys. 78, 114001 (2015).
Gre15 Ya.S. Greenberg and A.A. Shtygashev, Phys. Rev. A 92, 063835 (2015).
Ele17a H. Eleuch and I. Rotter, Phys. Rev. A 95, 022117 (2017).
Ele17b H. Eleuch and I. Rotter, Phys. Rev. E 95, 062109 (2017).
Jos18 S. Joshi and I. Galbraith, Phys. Rev. A 98, 042117 (2018).
Kat66 T. Kato, Perturbation Theory of Linear Operators (Springer, New York, 1966).
Zir83 M. R. Zirnbauer, J.J.M. Verbaarschot, and H.A. Weidenmüller, Nucl. Phys. A 411, 161 (1983).
Hei89 W. D. Heiss, Z. Phys. A 329, 133 (1989).
Cej07 P. Cejnar, S. Heinze, and M. Macek, Phys. Rev. Lett. 99, 100601 (2007).
Sin17 M. Šindelka, L. F. Santos, and N. Moiseyev, Phys. Rev. A 95, 010103(R) (2017).
Str18 P. Stránský, M. Dvořák and P. Cejnar, Phys. Rev. E 97, 012112 (2018).
Hei98 W.D. Heiss, M. Müller, and I. Rotter, Phys. Rev. E 58, 2894 (1998).
Jun99 C. Jung, M. Müller, and I. Rotter, Phys. Rev. E 60, 114 (1999).
Lip65 H. J. Lipkin, N. Meshkov, and A. J. Glick, Nucl. Phys. 62, 188 (1965).
Hol40 T. Holstein and H. Primakoff, Phys. Rev. 58, 1098 (1940).
Str19 M. Macek, P. Stránský, A. Leviatan and P. Cejnar, Phys. Rev. C 99, 064323 (2019). 
Meh04 M.L. Mehta, Random Matrices (Academic, London, 2004). 
Wol80 H. Wolkowicz and G.P.H. Styan, Linear Algebra Appl. 29, 471 (1980).


