Superradiance in finite quantum systems randomly coupled to continuum

TL;DR

This paper investigates superradiance in finite open quantum systems with random couplings, revealing how decay properties and exceptional points influence the phenomenon, especially near quantum criticality.

Contribution

It introduces a comprehensive analysis of superradiance using random matrix theory and links the effect to the distribution of exceptional points and quantum criticality.

Findings

Superradiance is enhanced at quantum critical points.

Eigenstates separate into short- and long-living due to increasing decay widths.

Robust features are related to the distribution of exceptional points.

Abstract

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.

Equations

$${\hat{H}}^{\left(\{{\lambda }_l\}\right)}=\underset{{\hat{H}}^{(0)}}{\underbrace{\sum _{k=1}^d{E}_k^{(0)}|k⟩⟨k|}}+\sum _{l=1}^n\underset{{\lambda }_l}{\underbrace{({ϵ}_l-i{\gamma }_l)}}|{\varphi }_l⟩⟨{\varphi }_l|,$$

$${\hat{H}}^{(\lambda )}={\hat{H}}^{(0)}+\underset{\lambda }{\underbrace{(ϵ-i\gamma )}}\underset{{\hat{P}}_{\mathrm{D}}}{\underbrace{\sum _{l=1}^n|{\varphi }_l⟩⟨{\varphi }_l|}}.$$

$${\hat{H}}^{(0)}$$

$$=$$

$$E_{\kappa }^{(ϵ-i\gamma )}$$

$${\Gamma }_{\kappa }^{(ϵ-i\gamma )}$$

$$E_{\kappa }^{(ϵ-i\gamma )}$$

$${\Gamma }_{\kappa }^{(ϵ-i\gamma )}$$

$$\frac{1}{2}{N}_{\mathrm{EP}}=n(d-n).$$

$$|\lambda|\!\lesssim\!\frac{S\,d^{2}}{\sqrt{n(d\!-\!n)}}\!\propto\!\left\{\begin{array}[]{ll}d^{3/2}&{\rm for\ }n\!\ll\!d{\ \rm or\ }n_{\bot}\!\ll\!d,\\ d&{\rm for\ }n\!\approx\!n_{\bot}\!\approx\!d/2,\end{array}\right.$$

$${\overline{\mathcal{E}}}^{(\lambda )}\equiv \frac{1}{d}\sum _{\kappa =1}^d{\mathcal{E}}_{\kappa }^{(\lambda )}=\frac{1}{d}\mathrm{Tr}{\hat{H}}^{(\lambda )},$$

$${\overline{{\Delta }^2\mathcal{E}}}^{(\lambda )}$$

$${\overline{E}}^{(ϵ-i\gamma )}-{\overline{E}}^{(0)}$$

$${\overline{\Gamma }}^{(ϵ-i\gamma )}$$

$${\overline{{\Delta }^{2}E}}^{(ϵ-i0)}-{\overline{{\Delta }^{2}E}}^{(0)}=2ϵ\frac{n}{d}A+{ϵ}^2\frac{n}{d}\left(1-\frac{n}{d}\right),$$

$${\overline{{\Delta }^{2}E}}^{(ϵ-i\gamma )}-{\overline{{\Delta }^{2}E}}^{(ϵ-i0)}={\overline{{\Delta }^2\Gamma }}^{(ϵ-i\gamma )}-{\gamma }^2\frac{n}{d}\left(1-\frac{n}{d}\right).$$

$$E_1^{(ϵ-i0)}\le \mathrm{Re}{\mathcal{E}}_{\kappa }^{(ϵ-i\gamma )}$$

$$0\le -\mathrm{Im}{\mathcal{E}}_{\kappa }^{(ϵ-i\gamma )}$$

$${\hat{H}}^{(\lambda )}=\underset{{\hat{H}}^{(0)}}{\underbrace{\left(\begin{array}{cc}{E}_1^{(0)}& 0\\ 0& E_2^{(0)}\end{array}\right)}}+\underset{\lambda }{\underbrace{(ϵ-i\gamma )}}\underset{\left(\begin{array}{c}\cos \vartheta \\ \sin \vartheta \end{array}\right)(\cos \vartheta \sin \vartheta )}{\underbrace{\left(\begin{array}{cc}{\cos }^2\vartheta & \cos \vartheta \sin \vartheta \\ \sin \vartheta \cos \vartheta & {\sin }^2\vartheta \end{array}\right)}},$$

$${\mathcal{E}}_{1,2}^{(\lambda )}\equiv E_{1,2}^{(\lambda )}-i{\Gamma }_{1,2}^{(\lambda )}=\frac{E_1^{(0)}+E_2^{(0)}+\lambda }{2}$$

$$\displaystyle\pm\sqrt{\biggl{(}\frac{E^{(0)}_{1}\!-\!E^{(0)}_{2}}{2}\biggr{)}^{2}+\biggl{(}\frac{\lambda}{2}\biggr{)}^{2}+\lambda\,\frac{E^{(0)}_{1}\!-\!E^{(0)}_{2}}{2}\cos 2\vartheta}$$

$$\lambda_{{\rm EP}\pm}=-\bigl{(}E^{(0)}_{1}-E^{(0)}_{2}\bigr{)}\,e^{\pm 2i\vartheta}\,.$$

$${\hat{H}}^{({\lambda }_{\mathrm{EP}\pm })}\to \left(\begin{array}{cc}{\mathcal{E}}_{\mathrm{EP}\pm }& 1\\ 0& {\mathcal{E}}_{\mathrm{EP}\pm }\end{array}\right),$$

Full text

I Introduction

Imagine a classical system that randomly decays with a rate $\gamma$ whenever it passes a certain decay zone in the phase space. For time $\Delta t$ spent in the zone, $e^{-\gamma\Delta t}$ is the survival probability and ${(1\!-\!e^{-\gamma\Delta 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_{\rm min}$ and $E_{\rm 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 $\Gamma(E)$, which depend on details of dynamics but are proportional to $\gamma$ 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 ${\cal H}$ of state vectors. The quantum decay zone is defined as a certain $n$-dimensional subspace ${{\cal H}_{\rm D}\subset{\cal H}}$, which couples uniformly to an infinite environment with continuous spectrum. Consequently, all states within ${\cal H}_{\rm D}$ have the same decay width $\gamma\!=\!\hbar/(2\tau)$, where $\tau$ stands for the mean lifetime (here we define “energy width” consistently with relation ${\gamma\tau=\frac{1}{2}\hbar}$ instead of a more common definition without factor $\frac{1}{2}$). Since any eigenstate of an independently chosen Hamiltonian has generally a non-vanishing overlap with ${\cal H}_{\rm D}$, the whole system under the Hamiltonian-induced evolution becomes unstable. One can apply a simple description based on the Feshbach theory Fes ; Moi11 , in which the open system is attributed by an effective non-Hermitian Hamiltonian $\hat{H}^{(\gamma)}\!=\!{\hat{H}^{(0)}\!-\!i\gamma\hat{P}_{\rm D}}$, with $\hat{H}^{(0)}$ denoting Hamiltonian of the unperturbed closed system and $\hat{P}_{\rm D}$ a projector to ${\cal H}_{\rm D}$. The eigenstates $|\kappa^{(\gamma)}\rangle$ of $\hat{H}^{(\gamma)}$ for any finite value of $\gamma$ have generally non-zero decay widths $\Gamma^{(\gamma)}_{\kappa}\!=\!\hbar/(2\tau^{(\gamma)}_{\kappa})$. If $\gamma$ grows infinitesimally from zero, $\Gamma^{(\gamma)}_{\kappa}$ of all eigenstates increase proportionally. However, as the width $\gamma$ grows to higher values, the set of eigenstates eventually splits into two groups: The first one contains $n$ states whose widths keep growing as $\Gamma^{(\gamma)}_{\kappa}\!\approx\!\gamma$. In contrast, the widths of remaining ${(d-n)}$ states in the second group overcome the initially increasing trend and return back to zero, ${\Gamma^{(\gamma)}_{\kappa}\to 0}$, with ${\gamma\to\infty}$. For a sufficiently high value of $\gamma$, 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. 1 below.

The above-described phenomenon was probably for the first time pointed out in Ref. Kle85 . Its mechanism was soon related Sok88 ; Sok92 to so-called Dicke superradiance Dic54 , 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. Rot91 ; Vol04 ; Aue11 ; Cel11 ; Liu14 ; Ele14 ; Rot15 ; Gre15 ; Ele17a ; Ele17b ; Jos18 . 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 Liu14 ; Rot15 ; Gre15 ; Jos18 . For detailed reviews with additional references see Refs. Aue11 ; Ele14 ; Rot15 .

Essential insight into the mechanism of the NHSR follows from the mathematics of so-called exceptional points (EPs) Kat66 ; Zir83 . These represent degeneracies of a Hamiltonian $\hat{H}^{(\lambda)}$ with a real discrete spectrum depending on parameter $\lambda$ in complex-extended domain ${\lambda\in\mathbb{C}}$, where $\hat{H}^{(\lambda)}$ is non-Hermitian and its eigenvalues complex. Convergence of EPs to a point $\lambda_{\rm c}$ on the real parameter axis with asymptotically increasing size of the system was shown to trigger quantum phase transitions (QPTs) Hei89 ; Cej07 ; Sin17 ; Str18 . 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 $\gamma$ Ele14 ; Hei98 ; Jun99 .

The role of EPs in the NHSR process was studied in Refs. Ele14 ; Hei98 ; Jun99 . 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. II), demonstrate on it a general effect of EPs (Sec. III), and derive some overall properties of the complex spectrum (Sec. IV). 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. V). At last, we summarize the results (Sec. VI).

II 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 ${\cal H}$, in which $n\in[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. 1. Everywhere we set ${\hbar=1}$.

At first we consider a general Hamiltonian of the form

[TABLE]

where the first term $\hat{H}^{(0)}$ (the unperturbed initial Hamiltonian) describes a closed system with energies ${E_{k}^{(0)}\in\mathbb{R}}$ and orthonormal eigenvectors $|k\rangle$, while the second term (decay-inducing Hamiltonian) defines $n$ decaying states given by orthonormal vectors $|\phi_{l}\rangle$ with complex energies ${\lambda_{l}\in{\mathbb{C}}}$, each composed of a real part ${\epsilon_{l}\in\mathbb{R}}$ and a non-negative decay width ${\gamma_{l}\in\mathbb{R}}$. The orthonormality conditions read as ${\langle k_{1}|k_{2}\rangle=\delta_{k_{1}k_{2}}}$ and ${\langle\phi_{l_{1}}|\phi_{l_{2}}\rangle=\delta_{l_{1}l_{2}}}$, where both vector sets $\{|k\rangle\}_{k=1}^{d}$ and $\{|\phi_{l}\rangle\}_{l=1}^{n}$ are supposed to be incompatible.

To illustrate physical meaning of the schematic Hamiltonian (1), we assume a system of $N$ fermions in a finite Hilbert space generated by $d_{\rm sp}$ single-particle states in a mean-field potential well. Out of these states, $(d_{\rm sp}\!-\!n_{\rm sp})$ are stable, bounded inside the well below the continuum threshold energy, while the remaining $n_{\rm 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\!=\!\binom{d_{\rm sp}}{n_{\rm sp}}$ is spanned by $\binom{d_{\rm sp}-n_{\rm sp}}{N}$ stable and $n\!=\!\binom{d_{\rm sp}}{N}\!-\!\binom{d_{\rm sp}-n_{\rm sp}}{N}$ unstable mean-field configurations. The unstable subspace is generated by $N$-body basis vectors $\{|\phi_{l}\rangle\}_{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 $\{|\phi_{l}\rangle\}_{l=n+1}^{d}$. Now we assume that all fermions are subject to some residual two-body interactions acting between all $d_{\rm 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 $\hat{H}^{(0)}$ which approximately describes the system without tunneling. The eigenbasis $\{|k\rangle\}_{k=1}^{d}$ of $\hat{H}^{(0)}$ is rotated relative to the mean-field basis $\{|\phi_{l}\rangle\}_{l=1}^{d}$. The decay describing part of the Hamiltonian consists of projectors to original unstable states $\{|\phi_{l}\rangle\}_{l=1}^{n}$ with complex coefficients expressing decay widths $\gamma_{l}$ and real-energy corrections $\epsilon_{l}$ depending on the shape of the confining potential barrier. This results in a non-Hermitian Hamiltonian of the form (1).

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

[TABLE]

Here, the decay-inducing part of the Hamiltonian $\lambda\hat{P}_{\rm D}=\hat{H}_{{\rm D}}^{(\lambda)}$ is determined by a single projection operator $\hat{P}_{\rm D}$ and all the decaying states, which form a subspace ${{\cal H}_{\rm D}=\hat{P}_{\rm D}{\cal H}\subset{\cal H}}$, have the same complex energy ${\epsilon-i\gamma=\lambda}$. The imaginary component $\gamma$ expresses the decay width of states in ${\cal H}_{\rm D}$ and the real component $\epsilon$ represents their energy shift. We stress that though the simplified Hamiltonian (2) restricts direct applicability of our model to specific systems, it captures essential features of the more complex Hamiltonian (1) while considerably reducing the number of free parameters.

The decay-inducing term of Hamiltonian (2) can be cast in the form ${\hat{H}_{{\rm D}}^{(\lambda)}=\lambda\hat{P}_{\rm D}+0\hat{P}_{\bot}}$ (with $\hat{P}_{\bot}=\hat{\mathbb{I}}-\hat{P}_{\rm D}$) with two eigenvalues $\lambda$ and 0. The eigenspace associated with eigenvalue $\lambda$ coincides with the decaying subspace ${\cal H}_{\rm D}$ of dimension $n$, while the eigenspace with eigenvalue [math] is the subspace ${{\cal H}_{\bot}=\hat{P}_{\bot}{\cal H}}$ with dimension ${n_{\bot}=d-n}$. Note that the cases ${n=0}$ with ${\hat{H}_{{\rm D}}^{(\lambda)}=0}$ and ${n=d}$ with ${\hat{H}_{{\rm D}}^{(\lambda)}=\lambda\hat{\mathbb{I}}}$ are both trivial and we exclude them.

In our analysis, the unperturbed component $\hat{H}^{(0)}$ of the total Hamiltonian (2) is associated with the Lipkin-Meshkov-Glick model Lip65 . It is build from quasispin operators ${(\hat{J}_{1},\hat{J}_{2},\hat{J}_{3})}$ satisfying the standard angular-momentum commutation relations. The conserved quantity $\hat{J}^{2}$ is fixed at the value $j(j+1)$, so that the Hilbert space ${\cal H}$ spanned by eigenvectors of $\hat{J}_{3}$ has dimension ${d=2j+1}$. To generate the spectrum of energies $E_{k}^{(0)}$, we consider the following three alternative forms:

[TABLE]

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 (3)–(5). We introduce cumulants expressing global properties of the unperturbed spectrum, in particular the average $\overline{E}^{(0)}=\frac{1}{d}\sum_{k=1}^{d}E_{k}^{(0)}$, which represents a central energy of the spectrum, the variance \overline{\Delta^{2}E}^{(0)}=\frac{1}{d}\sum_{k=1}^{d}\bigl{[}E_{k}^{(0)}\!-\!\overline{E}^{(0)}\bigr{]}^{2}, which characterizes a squared spread of the spectrum, and the skewness \overline{\Delta^{3}E}^{(0)}=\frac{1}{d}\sum_{k=1}^{d}\bigl{[}E_{k}^{(0)}\!-\!\overline{E}^{(0)}\bigr{]}^{3}, which quantifies an asymmetry of the spectrum with respect to $\overline{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 $\lambda$.

All Hamiltonians (3)–(5) can be written, using the Holstein-Primakoff transformation Hol40 , in terms of a single pair ${(\hat{q},\hat{p})}$ of conjugate coordinate and momentum operators satisfying in the classical limit a constraint ${q^{2}+p^{2}\leq 2}$ (for a general discussion of this mapping see, e.g., Str19 ). The Hamiltonian $\hat{H}^{(0)}_{\rm 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)\propto q^{2}+{\rm const}}$. Hamiltonians $\hat{H}^{(0)}_{\rm PT1}$ and $\hat{H}^{(0)}_{\rm PT2}$ (“phase transitional” of the first and second kind) describe representative $d$-dimensional spectra of two types of quantum critical systems: $\hat{H}^{(0)}_{\rm 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)\propto 3q^{4}-5q^{2}+{\rm const}}$. The spacing between the low-lying states forming the doublets decreases with dimension as ${E^{(0)}_{k+1}-E^{(0)}_{k}\propto e^{-c_{k}d}}$ (where $c_{k}$ is a constant), while spacing between states outside the doubles remains of order ${E^{(0)}_{k+2}-E^{(0)}_{k+1}\sim O(1)}$. $\hat{H}^{(0)}_{\rm PT2}$ represents a second-order QPT Hamiltonian associated with a quartic potential ${V(q)\propto q^{4}+{\rm const}}$. Its spectrum exhibits a typical power-law cumulation of low-energy levels according to ${E^{(0)}_{k+1}-E^{(0)}_{k}\propto\left(k/d\right)^{1/3}}$. Energy levels of all the three initial Hamiltonians are for a moderate dimension shown in Fig. 2. More details concerning the critical Hamiltonians can be found in Ref. Str18 .

The decay-inducing component of Hamiltonian (2) is represented not by a single fixed operator $\hat{H}_{{\rm D}}^{(\lambda)}$, 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_{\rm R}$ (of order from $10^{1}$ to $10^{4}$) of samples from this ensemble. The orthonormal vectors $\{|\phi_{l}\rangle\}_{l=1}^{n}$ and $\{|\phi_{\bot l^{\prime}}\rangle\}_{l^{\prime}=1}^{n_{\bot}}$ forming in each realization the bases of subspaces ${\cal H}_{\rm D}$ and ${\cal H}_{\bot}$, respectively, result from a random orthogonal transformation of the original eigenbasis $\{|k\rangle\}_{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 Meh04 . In particular, for each realization we perform the following steps: (i) we generate a random $d$-dimensional real symmetric matrix ${\hat{H}_{\rm GOE}\in{\rm GOE}}$, i.e., a matrix with independent elements taken from zero-mean Gaussian distributions with ${\sigma^{2}=2}$ or 1 for diagonal or off-diagonal elements, respectively, (ii) we find an orthogonal matrix $\hat{O}$ transforming $\hat{H}_{\rm GOE}$ to the diagonal form ${\hat{D}=\hat{O}^{\rm T}\hat{H}_{\rm GOE}\hat{O}}$, where ${\hat{O}^{\rm T}=\hat{O}^{{\dagger}}=\hat{O}^{-1}}$ stands for the transpose of $\hat{O}$, (iii) we randomly chose $n$ columns of $\hat{O}$ and associate them with the decaying states $\{|\phi_{l}\rangle\}_{l=1}^{n}$. This procedure is repeated until a large number $N_{\rm R}$ of $\hat{H}_{\rm GOE}$ realizations yields a robust estimate of the analyzed quantities. The averaging over the random-matrix ensemble will be further denoted by angular brackets $\langle\bullet\rangle$, which is in contrast to the bar symbol $\overline{\bullet}$ introduced above for the “non-statistical” averages over the energy levels. Note that orthogonality of the transformation $\hat{O}$, which is inherent in the foundations of GOE, guarantees that matrix elements ${\langle k_{1}|\hat{P}_{\rm D}|k_{2}\rangle=\sum_{l=1}^{n}\langle k_{1}|\phi_{l}\rangle\langle\phi_{l}|k_{2}\rangle=\sum_{l=1}^{n}O_{k_{1}l}O_{k_{2}l}}$ are real, so that we can write: ${{\rm Im}{\langle k_{1}|\hat{H}(\lambda)|k_{2}\rangle}\propto{\rm Im}\lambda=-\gamma}$.

III Bifurcation of decay widths

In the following, we focus on the dependencies of the eigensolutions of the Hamiltonian (2) on parameter ${\lambda=\epsilon-i\gamma}$. The Hamiltonian is a complex matrix symmetric under transposition (but not under the full Hermitian conjugation) with $d$ generally complex eigenvalues ${{\cal E}_{\kappa}^{(\lambda)}=E_{\kappa}^{(\lambda)}-i\Gamma_{\kappa}^{(\lambda)}}$. The real parts $E_{\kappa}^{(\lambda)}$ stand for energies and the imaginary parts ${\Gamma_{\kappa}^{(\lambda)}\geq 0}$ represent decay widths of individual eigenstates enumerated by integer ${\kappa=1,2,\dots,d}$. Because of non-Hermiticity, the eigenstates of $\hat{H}^{(\lambda)}$ must be distinguished to right eigenstates satisfying ${\hat{H}^{(\lambda)}|\kappa_{\rm R}^{(\lambda)}\rangle={\cal E}_{\kappa}^{(\lambda)}|\kappa_{\rm R}^{(\lambda)}\rangle}$, and left eigenstates satisfying $\langle\kappa_{\rm L}^{(\lambda)}|{\hat{H}^{(\lambda)}=\langle\kappa_{\rm L}^{(\lambda)}|{\cal E}_{\kappa}^{(\lambda)}}$. From the transposition symmetry of $\hat{H}^{(\lambda)}$ 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 ${\langle{\kappa_{1}}_{\rm L}^{(\lambda)}|{\kappa_{2}}_{\rm R}^{(\lambda)}\rangle=\delta_{\kappa_{1}\kappa_{2}}}$.

We will assume that $\epsilon\!=\!{\rm Re}\lambda$, is set constant, while $\gamma\!=\!-{\rm Im}\lambda$ is varied within the domain ${\gamma\geq 0}$. We can obviously write $\hat{H}^{(\epsilon-i\gamma)}=\hat{H}^{(\epsilon-i0)}-i\gamma\hat{P}_{\rm D}$, where we explicitly introduce a shifted Hermitian Hamiltonian ${\hat{H}^{(\epsilon-i0)}=\hat{H}^{(0)}+\epsilon\hat{P}_{\rm D}}$ and an anti-Hermitian decay-inducing Hamiltonian ${\hat{H}_{\rm D}^{(0-i\gamma)}=-i\gamma\hat{P}_{\rm D}}$. Evolutions of decay widths for several model settings, averaged over a sample of GOE realizations $N_{\rm R}$, are depicted in Fig. 3. We see that the instability expressed by positive decay widths $\Gamma_{\kappa}^{(\lambda)}$ characterizes all eigenstates of $\hat{H}^{(\lambda)}$ as soon as $\gamma$ deviates from zero to infinitesimally small positive values. Indeed, all eigenstates are expected to have non-zero overlaps with a randomly generated decaying subspace ${\cal H}_{\rm D}$. Using the Hellmann-Feynman formula, which in the non-Hermitian context reads $\frac{d}{d\lambda}{\cal E}_{\kappa}^{(\lambda)}\!=\!{\langle\kappa_{\rm L}^{(\lambda)}|\frac{d}{d\lambda}\hat{H}^{(\lambda)}|\kappa_{\rm R}^{(\lambda)}\rangle}$, we derive the following relations: $\frac{\partial}{\partial\gamma}E_{\kappa}^{(\lambda)}\!=\!{-{\rm Im}\langle\kappa_{\rm L}^{(\lambda)}|\hat{P}_{\rm D}|\kappa_{\rm R}^{(\lambda)}\rangle}$ and $\frac{\partial}{\partial\gamma}\Gamma_{\kappa}^{(\lambda)}\!=\!{{\rm Re}\langle\kappa_{\rm L}^{(\lambda)}|\hat{P}_{\rm D}|\kappa_{\rm R}^{(\lambda)}\rangle}$. These can be used for an estimation of $E_{\kappa}^{(\lambda)}$ and $\Gamma_{\kappa}^{(\lambda)}$ for small values of $\gamma$, where eigenvectors $|\kappa_{\rm R}^{(\lambda)}\rangle$ and $|\kappa_{\rm L}^{(\lambda)}\rangle$ roughly coincide with eigenvectors $|\kappa^{(\epsilon-i0)}\rangle$ of the Hermitian Hamiltonian $\hat{H}^{(\epsilon-i0)}$. We get

[TABLE]

where ${P_{\kappa}=\sum_{l=1}^{n}|\langle\phi_{l}|\kappa^{(\epsilon-i0)}\rangle|^{2}}$ is the probability of identifying the eigenstate $|\kappa^{(\epsilon-i0)}\rangle$ with any of the decaying states $|\phi_{l}\rangle$.

On the other hand, it is clear that for very large values of ${\gamma>0}$ the decay-inducing part of the Hamiltonian dominates, so that the eigenvectors of the full Hamiltonian (2) become approximately those of the second term. Hence for asymptotic $\gamma$ only ${n<d}$ eigenstates (those roughly coinciding with the set ${\{|\phi_{l}\rangle\}_{l=1}^{n}\!\in\!{\cal H}_{\rm D}}$) will yield non-negligible widths $\Gamma_{\kappa}^{(\lambda)}\!\approx\!\gamma$, while the rest of ${n_{\bot}=d-n}$ eigenstates (those approximated by $\{|\phi_{\bot l^{\prime}}\rangle\}_{l^{\prime}=1}^{n_{\bot}}\!\in\!{\cal H}_{\bot}$) will have ${\Gamma_{\kappa}^{(\lambda)}\!\approx\!0}$. Indeed, the Hellmann-Feynman formula gives $\frac{\partial}{\partial\gamma}E_{\kappa}^{(\lambda)}\!\approx\!0$ for all eigenstates, and $\frac{\partial}{\partial\gamma}\Gamma_{\kappa}^{(\lambda)}\!\approx\!1$ or 0, depending on whether the approximate eigenvectors belong to ${\cal H}_{\rm D}$ or ${\cal H}_{\bot}$, respectively. A more systematic treatment makes use of a transformed Hamiltonian $i\gamma^{-1}\hat{H}^{(\epsilon-i\gamma)}=\hat{P}_{\rm D}+i\gamma^{-1}\hat{H}^{(\epsilon-i0)}$ and a perturbative expansion of its complex eigenvalues in small parameter $\gamma^{-1}$. This yields:

[TABLE]

where $c_{\kappa}$ are some coefficients. The log-log dependencies in Fig. 3 support these conclusions.

Having understood the evolution of the spectrum for very small and very large values of $\gamma$, we ask what drives the transition between these limiting regimes. The region in parameter $\gamma$ where the bifurcation of decay widths takes place is scrutinized under a specific magnifying glass in Fig. 4. The vertical axis of this figure shows the log-log slope $\chi_{\kappa}^{(\lambda)}\!=\!\frac{d}{d\log\gamma}\log\Gamma_{\kappa}^{(\lambda)}\!=\!(\gamma/\Gamma_{\kappa}^{(\lambda)})\frac{d}{d\gamma}\Gamma_{\kappa}^{(\lambda)}$ of individual decay widths from the $\epsilon\!=\!0$ panel of Fig. 3. The value $\chi_{\kappa}^{(\lambda)}\!=\!+1$ indicates the $\Gamma_{\kappa}^{(\lambda)}\!\propto\!\gamma$ behavior, while $\chi_{\kappa}^{(\lambda)}\!=\!-1$ implies $\Gamma_{\kappa}^{(\lambda)}\!\propto\!\gamma^{-1}$. We observe that the transitional region between these values is a relatively narrow interval of $\gamma$, for the selected values of $d$ and $n$ roughly given by ${1\lesssim\gamma\lesssim 20}$. In this region, the slopes in Fig. 4 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 $\hat{H}^{(\lambda)}$ is a point $\lambda_{\rm EP}$ of the complex parameter space where two (or more) complex eigenvalues ${\cal E}_{\kappa}^{(\lambda)}$ and ${\cal E}_{\kappa^{\prime}}^{(\lambda)}$ become degenerate. In a typical case, the complex degeneracy has a different character than an ordinary degeneracy of a Hermitian Hamiltonian. First, $\hat{H}^{(\lambda)}$ at ${\lambda=\lambda_{\rm EP}}$ fails to provide a complete system of left and right eigenvectors, but the single eigenvector associated with both degenerate levels becomes self-orthogonal: ${\langle\kappa_{\rm L}^{(\lambda_{\rm EP})}|\kappa_{\rm R}^{(\lambda_{\rm EP})}\rangle=0}$. Second, the EP is not just a conical intersection of the two energy surfaces but a complex square-root singularity. The path in ${\lambda\in\mathbb{C}}$ satisfying the condition ${{\rm Im}\,{\cal E}_{\kappa}^{(\lambda)}={\rm Im}\,{\cal E}_{\kappa^{\prime}}^{(\lambda)}}$ is terminated at $\lambda_{\rm EP}$, where it smoothly connects to herefrom issuing path defined by ${{\rm Re}\,{\cal E}_{\kappa}^{(\lambda)}={\rm Re}\,{\cal E}_{\kappa^{\prime}}^{(\lambda)}}$. On the respective sides of the EP, the real and imaginary parts of ${{\cal E}_{\kappa}^{(\lambda)}\!-\!{\cal E}_{\kappa^{\prime}}^{(\lambda)}}$ 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 $\lambda$, 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\times 2$ matrices. We sketch such a trivial ${d=2}$ version of our model in Appendix A (cf. Ref. Hei98 ). It is shown that setting ${\epsilon={\rm Re}\lambda_{\rm EP}}$, so that the path with increasing $\gamma$ crosses the EP of the Hamiltonian at ${\gamma=|{\rm Im}\lambda_{\rm EP}|}$, we observe a sudden bifurcation of the widths $\Gamma_{1}^{(\lambda)}$ and $\Gamma_{2}^{(\lambda)}$ from a common value equal to $\frac{1}{2}\gamma$ 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 $\epsilon$ from ${\rm Re}\lambda_{\rm EP}$, the decay widths of both levels differ already at small $\gamma$ so that the ratio between the widths of short- and long-living state increases with ${|\epsilon\!-\!{\rm Re}\lambda_{\rm EP}|}$. At the beginning, for $\gamma$ not exceeding $|{\rm Im}\lambda_{\rm EP}|$, both widths grow proportionally to $\gamma$. However, for larger $\gamma$, somewhere at $\gamma\gtrsim|{\rm Im}\lambda_{\rm EP}|$, the smaller of both widths turns back and starts decreasing as $\gamma^{-1}$. This is a smooth but still clear realization of the NHSR effect. Examples of these scenarios, depicted in Fig. 11 below, can be compared to less trivial but rather analogous behavior for ${d>2}$ shown in Fig. 3.

Focusing on these ${d>2}$ cases, we first note that a symmetric matrix of dimension $d$ depending linearly on parameter ${\lambda\in{\mathbb{C}}}$ has in general $\frac{1}{2}d(d\!-\!1)$ pairs of complex-conjugate EPs. However, in case of Hamiltonian (2) a large part of EP pairs migrates to infinity, while at finite values of $\lambda$ we observe a reduced number of pairs:

[TABLE]

This is $\lesssim 50\,\%$ of the full number of pairs (the maximal fraction is reached for ${n=\frac{1}{2}d}$). Figure 5 shows distributions of these EPs in ${\lambda\in{\mathbb{C}}}$ for the three initial Hamiltonians (3)–(5) and for various sizes of the decaying subspace ${\cal H}_{\rm D}$. The distribution is averaged over a large number of GOE realizations of ${\cal H}_{\rm D}$.

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

[TABLE]

where S\!=\!\bigl{(}\overline{\Delta^{2}E}^{(0)}/d^{2}\bigr{)}^{1/2} is the linear spread of the unperturbed spectrum divided by the dimension (e.g., for $H^{(0)}_{\rm HO}$ and $H^{(0)}_{\rm PT2}$ with unit average spacings we get $S\!\approx\!0.3$). The scaling formula (13) follows from the analysis of Ref. Str18 (cf. Zir83 ) which showed that a maximum of the EP distribution for a general perturbed Hamiltonian ${{\hat{H}}^{(0)}+\lambda\hat{V}}$ is achieved when spectral variances (spreads of the spectra) of both terms $H^{(0)}$ and $\lambda\hat{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. 3 (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 $\gamma$ value where the longer-living group turns back to the $\Gamma\!\to\!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 $\gamma$ to the increasing-$\Gamma$ direction) from a dynamic bi-modality (existence of two peaks moving to opposite $\Gamma\!\to\!0$ and $\Gamma\!\to\!\infty$ 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 $\epsilon$ either exceeds the upper peripheral region of the EP distribution projected to the ${\rm Re}\lambda$ axis, or undershoots its lower peripheral region. Analogously, the dynamic bimodality occurs when the value of $\gamma$ passes the upper periphery of the EP distribution projected of to the ${\rm Im}\lambda$ axis. However, it needs to be stressed that details of the formation of both kinds of bi-modalities for $\lambda$ 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. 6 several transient distributions of decay widths for the same model setup as used in the middle row of Fig. 5 (${d=2n=16}$). Each panel depicts a distribution of widths averaged over a large number of GOE realizations for a fixed parameter point ${\lambda=\epsilon-i\gamma}$ together with arrows, which indicate the direction in which a given part of the distribution evolves with increasing $\gamma$. Note that some of the distributions represent $\gamma$=const. cuts of Figs. 3 and 4. The fixed parameter points in Fig. 6 are organized so that panel (a) is associated with $\lambda$ 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 ${\lambda\in{\mathbb{C}}}$ 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 $\gamma$ is not large enough. The dynamic bi-modality is triggered only when $\gamma$ passes the upper periphery of the projected EP distribution.

The distribution of EPs in Fig. 5 for ${n=\frac{1}{2}d}$, which is the case corresponding to Fig. 6, is approximately symmetric under rotations around ${\lambda=0}$. However, the distributions for $n\!=\!1$ and $n=d\!-\!1$ with ${H^{(0)}=H^{(0)}_{\rm PT1}}$ and $H^{(0)}_{\rm PT2}$ exhibit strong asymmetries under rotations and even under the mirror transformation $\epsilon\leftrightarrow-\epsilon$. In these cases, the NHSR dynamics sensitively depends on the sign of $\epsilon$. This is illustrated in Fig. 7, where we show the decay width distributions at two points with ${\rm Re}\lambda\!=\!\pm\epsilon$ and ${\rm Im}\lambda$ small for $H^{(0)}_{\rm PT2}$. Both points belong to the regions within the principal parts of the EP distributions in the rightmost column of Fig. 5. We see that in the $n\!=\!1$ and $n=d\!-\!1$ cases the decay width distributions have significantly different forms for both $\pm\epsilon$ values, in agreement with the observed asymmetries of the respective EP distributions. In both these cases, the width distribution corresponding to the $\epsilon$ value closer to the main maximum of the EP distribution shows a more developed static bi-modality. However, the width distributions for $\pm\epsilon$ differ also in the $n\!=\!\frac{1}{2}d$ 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. Str18 ).

IV Global properties of the complex spectrum

Employing the trace-based method described e.g. in Ref. Str18 (see also Wol80 ), one can derive expressions for cumulants of the whole eigenvalue spectrum of a general (parameter-dependent) Hamiltonian $\hat{H}^{(\lambda)}$. In particular, the average of all eigenvalues is determined from

[TABLE]

and the eigenvalue variance is given by

[TABLE]

For the non-Hermitian Hamiltonian (2), the averages and variances related to complex energies ${\cal E}_{\kappa}^{(\lambda)}$ shall be evaluated separately for real and imaginary parts. From Eq. (14) we get easy formulas for both types of averages,

[TABLE]

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

[TABLE]

where {A=\frac{1}{n}\sum_{k=1}^{d}\bigl{(}\!E_{k}^{(0)}\!\!-\!\overline{E}^{(0)}\bigr{)}\langle k|\hat{P}_{\rm D}|k\rangle}, and

[TABLE]

Formulas (16) and (17) imply linear dependences of the average real energy and the average decay width on $\epsilon$ and $\gamma$, respectively. Formula (18) expresses a quadratic dependence of the real energy variance on $\epsilon$ for ${\gamma=0}$, while formula (19) captures a specific relation between the variances of real energies and decay widths for variable $\gamma$ and constant $\epsilon$.

The linear term of the quadratic dependence of the real energy variance in Eq. (18) 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)}\!\!-\!\overline{E}^{(0)})}$ calculated with normalized weight factors ${w_{k}=\frac{1}{n}\sum_{l=1}^{n}|\langle\phi_{l}|k\rangle|^{2}}$ proportional to the overlap probability of the respective unperturbed eigenstate with ${\cal H}_{\rm D}$. This asymmetry determines the point ${\epsilon_{\rm min}=-Ad/(d\!-\!n)}$ of a minimal quadratic spread $\overline{\Delta^{2}E}^{(\epsilon-i0)}$ of the real spectrum along the real $\lambda$ 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 $\lambda$ 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 $\langle A\rangle\!=\!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 $\langle\Delta^{2}A\rangle$ and $\langle\Delta^{3}A\rangle$ are correlated with the respective cumulants of the unperturbed energy spectrum $\overline{\Delta^{2}E}^{(0)}$ and $\overline{\Delta^{3}E}^{(0)}$. This qualitatively explains the above-mentioned $\epsilon\leftrightarrow-\epsilon$ asymmetries of EP distributions in Fig. 5 for ${H^{(0)}\!=\!H^{(0)}_{\rm PT1}}$ and $H^{(0)}_{\rm PT2}$, whose spectra are apparently skewed (cf. Fig. 2).

We now turn to the consequences of Eq. (19) for the NHSR dynamics. The left-hand side of this formula represents the change of variance of the real spectrum between $\gamma\!=\!0$ and $\gamma\!>\!0$ points with the same $\epsilon$. The right-hand side is the difference between the actual variance of the decay widths at ${\lambda=\epsilon-i\gamma}$ and the variance of a two-component ensemble consisting of $n$ widths equal to $\gamma$ and ${n_{\bot}=d-n}$ widths equal to 0. According to Ref. Wol80 , the spectrum of Hamiltonian (2) satisfies the following inequalities:

[TABLE]

The bounds $E_{1}^{(\epsilon-i0)}$ and $E_{d}^{(\epsilon-i0)}$ in (20) are the lowest and the highest eigenvalues of $\hat{H}^{(\epsilon-i0)}$, and similarly 0 and $\gamma$ in (21) are the lowest and highest eigenvalues of ${\gamma\hat{P}_{\rm D}}$. Inequality (21) implies that the right-hand side of Eq. (19) is semi-negative. The limiting zero value corresponds to ${\gamma=0}$, while for increasing $\gamma$ we expect negative values that remain of order ${\sim O(1)}$ even in the asymptotic case. Hence we see that the increase of the decay rate $\gamma$ 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 ${\gamma\to\infty}$ limit. As shown in Fig. 8, the maximal contraction of the real energy spectrum, reaching the asymptotic value of 50 % of the squared spread of the original spectrum (${\approx 71\,\%}$ of the linear spread), is obtained for ${n=\frac{1}{2}d}$.

V Effects of criticality

As seen in Fig. 6, the distributions of decay widths for various choices (3), (4) and (5) of the initial Hamiltonian $\hat{H}^{(0)}$ are rather close to each other for moderate dimensions $d$. However, the differences get much more significant as $d$ grows. In Fig. 9 we show an analog of Fig. 6 for ${d=1024}$. In comparison to the previous figure, the parameter coordinates ${[\epsilon,\gamma]}$ corresponding to panels (a)–(i) were scaled roughly according to Eq. (13). We immediately see that the choice of $\hat{H}^{(0)}$ plays a substantial role. In general, the first- and second-order QPT Hamiltonians $\hat{H}^{(0)}_{\rm PT1}$ and $\hat{H}^{(0)}_{\rm PT2}$ result at various fixed values of $\lambda$ in much more developed bi-modal structures of the width distribution than the harmonic oscillator Hamiltonian $\hat{H}^{(0)}_{\rm HO}$. Interestingly, comparing pairs of panels (c)-(c’), (f)-(f’) and (i)-(i’) in Fig. 9, we observe clear differences between decay width distributions for both critical Hamiltonians at mirror-conjugate values of ${\rm Re}\lambda$. These again reflect the $-\epsilon\!\leftrightarrow\!+\epsilon$ asymmetries discussed in connection with Fig. 7(b) and Eq. (18).

The differences of decay width distributions in Fig. 9 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 $\gamma$, but more substantially relies on the reversed evolution of both modes as $\gamma$ grows to large values and on the convergence of the long-living mode to $\Gamma\!=\!0$. Whether the criticality of $\hat{H}^{(0)}$ implies some enhancement of the latter process is studied in Fig. 10.

This figure compares the smallest decay widths $\Gamma_{\rm min}^{(\lambda)}(\rm HO)$ and $\Gamma_{\rm min}^{(\lambda)}(\rm PT2)$ obtained by minimization over the entire set of all $d$ levels for the harmonic oscillator Hamiltonian $\hat{H}^{(0)}_{\rm HO}$ and the second-order QPT Hamiltonian $\hat{H}^{(0)}_{\rm PT2}$. The imaginary part of $\lambda$ is taken very large, $\gamma\!=\!16d$, so that we are far above the superradiant transition to decreasing width regime, the real part of $\lambda$ is set to $\epsilon\!=\!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\!=\!\frac{1}{2}d$. The minimal decay widths and their ratios were averaged over several GOE realizations, yielding the dots in Fig. 10. The error bars indicate standard deviations of the distributions obtained. Note that the numbers of realizations $N_{\rm 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 $\hat{H}^{(0)}_{\rm PT1}$, for which the parity doublets of levels become too close with increasing $d$.

We see in the inset of Fig. 10 that for $d\!\gtrsim\!2^{10}$ both GOE averages of the minimal decay width scale roughly as $\Gamma_{\rm min}^{(\lambda)}\!\propto\!d$. This is consistent with the choice $\gamma\!=\!16d$ since, with respect to the scaling formula (13), the pattern in Fig. 3 remains roughly invariant under the change of dimension if variables $\gamma$ and $\Gamma$ are replaced by $\gamma/d$ and $\Gamma/d$, respectively.

The most important information follows from the ratio between both widths for large $d$. A plot of the GOE average $\langle\Gamma_{\rm min}^{(\lambda)}({\rm PT2})/\Gamma_{\rm min}^{(\lambda)}({\rm HO})\rangle$ at $\lambda\!=\!0\!-\!i(16d)$ is shown in the main panel of Fig. 10. It suggests that the asymptotic value of the average ratio is roughly equal to 0.55. Since in the large-$\gamma$ domain the evolution of widths is described by Eq. (11), hence ${\Gamma_{\rm min}^{(\lambda)}\approx c_{\rm min}\gamma^{-1}}$, the above numerical value is assigned also to the average ratio $\langle c_{\rm min}({\rm PT2})/c_{\rm min}({\rm HO})\rangle$ of the respective coefficients. It is clear that the particular value 0.55 would not apply to other choices of $\epsilon$ (cf. Fig. 9), 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 $\hat{H}^{(0)}_{\rm 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 $\lambda\!=\!0$ and some of the levels exhibit the crossover to the ${\Gamma_{\kappa}^{(\lambda)}\propto\gamma^{-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 $\hat{H}^{(0)}_{\rm PT1}$, if it is made available for large-$d$ numerical simulations.

VI 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 $\Gamma\!\to\!\infty$ and $\Gamma\!\to\!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. 11), 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. 3) we observe laminar flows of real energies and decay widths described by the respective pair of equations (6)–(7) or (8)–(11).

(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. 4).

(iv) Essential information follows from projections of the EP distribution (cf. Fig. 5) to real and imaginary axes of the coupling parameter $\lambda\!=\!\epsilon\!-\!i\gamma$. Considering trajectories in the $\lambda$ plane with fixed $\epsilon$, we distinguish two basic ESQPT scenarios (see Fig. 3): First, if $\epsilon$ is within the support of the projected EP distribution on the real axis, the decay widths of individual levels for small $\gamma$ form a single-mode distribution. As $\gamma$ increases, the distribution moves to increasing values of $\Gamma$ and eventually splits into the increasing and decreasing branches when $\gamma$ gets outside the support of the projected EP distribution on the imaginary axis. Second, if $\epsilon$ is outside the support of the real projection of the EP distribution, the decay-width distribution is bi-modal already at small values of $\gamma$. However, the splitting into the increasing and decreasing branches again happens only when $\gamma$ gets out of the imaginary projection of the EP distribution. These conclusions are supported by examples of $\Gamma$ distributions in Figs. 6, 7 and 9, 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. (19)]. The opening of the system always leads to narrowing of the real energy spectrum (see Fig. 8).

(vi) Details of the NHSR dynamics depend on the structure of real energy spectrum of the $\gamma\!=\!0$ Hamiltonian. First, the formation of bi-modal structures in the decay width distribution depends on the value of $\epsilon$ in a way that goes beyond the overall EP distribution (see examples in Fig. 7). Second, even the asymptotic-$\gamma$ behavior carries traces of the $\gamma\!=\!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 $\lambda$ axis, which shifts the crossover of some states to the $\Gamma\!\to\!0$ stage of evolution to smaller values of parameter $\gamma$ and makes the stabilization process of these states more advanced (see Fig. 10). 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.

Appendix A Two-dimensional model

Hamiltonian (2) in the ${d=2}$ case reads

[TABLE]

where $\vartheta$ (assumingly ${\neq 0,\pi}$ since these special values would yield trivial solutions) is an angle determining the single decaying state $|\phi_{1}\rangle$. Note that for ${d=2}$ only the case of ${n=1}$ is non-trivial. Two complex eigenvalues

[TABLE]

become degenerate for $\lambda$ equal to

[TABLE]

These are complex conjugate EPs of the matrix (22), which both represent square-root singularities of the function ${{\cal E}_{1}^{(\lambda)}-{\cal E}_{2}^{(\lambda)}=\sqrt{(\lambda-\lambda_{{\rm EP}+})(\lambda-\lambda_{{\rm EP}-})}}$. At $\lambda=\lambda_{{\rm EP}\pm}$, the Hamiltonian can undergo a similarity transformation to the Jordan form

[TABLE]

where ${{\cal E}_{{\rm EP}\pm}={\cal E}_{1}^{(\lambda_{{\rm EP}\pm})}={\cal E}_{2}^{(\lambda_{{\rm EP}\pm})}}$. This means that two right eigenvectors $|1_{\rm R}^{(\lambda)}\rangle$ and $|2_{\rm R}^{(\lambda)}\rangle$ associated with the eigenvalues ${\cal E}_{1}^{(\lambda)}$ and ${\cal E}_{2}^{(\lambda)}$ for ${\lambda\neq\lambda_{{\rm EP}\pm}}$ contract to a single one $|1_{\rm R}^{(\lambda_{{\rm EP}\pm})}\rangle$ at the EP, and the same happens to left eigenvectors. The unique eigenvector at the EP is self-orthogonal in the sense ${\langle 1_{\rm L}^{(\lambda_{{\rm EP}\pm})}|1_{\rm R}^{(\lambda_{{\rm EP}\pm})}\rangle=0}$.

Setting the real part $\epsilon$ of parameter $\lambda$ to {\rm Re}\lambda_{{\rm EP}\pm}=\bigl{(}E^{(0)}_{2}\!-\!E^{(0)}_{1}\bigr{)}\cos 2\vartheta and varying the imaginary part along a half-line ${\gamma\in[0,\infty)}$, we proceed in opening the system along a trajectory passing the EP located in the lower complex half-plane of $\lambda$. It follows from Eq. (23) that for $\gamma$ increasing from zero towards the EP absolute value, the widths $\Gamma_{1}^{(\lambda)}$ and $\Gamma_{2}^{(\lambda)}$ are both equal to $\frac{1}{2}\gamma$, whereas the real energies $E_{1}^{(\lambda)}$ and $E_{2}^{(\lambda)}$ differ and collapse, according to the square-root formula, to the EP degeneracy as \gamma\to|{\rm Im}\lambda_{\rm EP\pm}|={|\bigl{(}E^{(0)}_{2}\!-\!E^{(0)}_{1}\bigr{)}\sin 2\vartheta|}. For $\gamma$ increasing further from ${|{\rm Im}\lambda_{\rm EP\pm}|}$, the real energies stay degenerate whereas the decay widths exhibit the square-root bifurcation to a short- and long-living asymptotic forms $\Gamma_{\rm short}^{(\lambda)}\!=\!\gamma$ and $\Gamma_{\rm long}^{(\lambda)}\!\propto\!\gamma^{-1}$. This critical (non-analytic) form of the superradiant scenario is illustrated in Fig. 11 for model settings in which $\lambda_{{\rm EP}\pm}=\pm i$.

For an increasing distance of $\epsilon$ from ${\rm Re}\lambda_{{\rm EP}\pm}$, the sharpness of the NHSR transition decreases, but its principal features remain preserved. These cases are also illustrated in Fig. 11. We see that shorter- and longer-living states exist already at small values of $\gamma$ where their widths increase linearly. The ratio between the larger and smaller width grows with the distance from ${\rm Re}\lambda_{{\rm EP}\pm}$. As $\gamma$ reaches values near or above ${\gamma\approx|{\rm Im}\lambda_{{\rm EP}\pm}|}$, one of the widths makes a crossover to the $\propto\gamma^{-1}$ behavior. At the same time, evolution of real energies with $\gamma$ increasing across the EP region for $\epsilon\!\neq\!{\rm Re}\lambda_{{\rm EP}\pm}$ shows a smooth reduction of the spacing $|E_{1}^{(\lambda)}\!-\!E_{2}^{(\lambda)}|$.