# On the density of complex eigenvalues of Wigner reaction matrix in a
disordered or chaotic system with absorption

Yan V. Fyodorov 1
1 King’s College London, Department of Mathematics, London WC2R 2LS, United
Kingdom

###### Abstract

In an absorptive system the Wigner reaction $K-$matrix (directly related to
the impedance matrix in acoustic or electromagnetic wave scattering) is non-
selfadjoint, hence its eigenvalues are complex. The most interesting regime
arises when the absorption, taken into account as an imaginary part of the
spectral parameter, is of the order of the mean level spacing. I show how to
derive the mean density of the complex eigenvalues for reflection problems in
disordered or chaotic systems with broken time-reversal invariance. The
computations are done in the framework of nonlinear $\sigma-$ model approach,
assuming fixed $M$ and $N\to\infty$. Some explicit formulas are provided for
zero-dimensional quantum chaotic system as well as for a semi-infinite
quasi-1D system with fully operative Anderson localization.

##  1 Introduction

Consider the problem of wave scattering from a piece of random medium confined
to a spatial domain ${\cal D}$ and described by a self-adjoint Hamiltonian
$H$, e.g.

$H=\sum_{{\bf r}\in\mathbf{\Lambda}}\,V({\bf r})\left|{\bf
r}\right\rangle\left\langle{\bf r}\right|+\sum_{{\bf r}\sim{\bf
r}^{\prime}}\left(t_{{\bf r}{\bf r}^{\prime}}\left.|{\bf
r}\right\rangle\left\langle{\bf r}^{\prime}\right|+t_{{\bf r}^{\prime}{\bf
r}}\left|{\bf r}^{\prime}\right\rangle\left\langle{\bf r}\right|\right),$ (1)

where the second sum runs over nearest neighbours on a lattice
$\mathbf{\Lambda}$ (assumed to be confined to the domain ${\cal D}$). The
parameters $t_{{\bf r}^{\prime}{\bf r}}$ are in general complex satisfying
$t^{*}_{{\bf r}{\bf r}^{\prime}}=t_{{\bf r}^{\prime}{\bf r}}$ to ensure the
Hermiticity of the Hamiltonian: $H=H^{\dagger}$, where we use $t^{*}$ to
denote complex conjugation of $t$ and $H^{\dagger}$ for Hermitian conjugation
of $H$. The disordered nature of the medium is taken into account by choosing
the on-site potentials $V({\bf r})$ and/or hopping parameters $t_{{\bf
r}^{\prime}{\bf r}}$ to be random variables. Such construction is known in the
literature as the Anderson model, and provides the paradigmatic framework to
study single-particle localization phenomena. Note that the form (1) can be
used also for modelling a quantum particle motion on any graph ${\bf
r}\in\mathfrak{G}$ , with $t_{{\bf r}{\bf r}^{\prime}}$ being the elements of
the adjacency matrix of the graph $\mathfrak{G}$

The tight-binding representation is convenient as it allows one to think of
such a Hamiltonian as described by a large $N\times N$ random matrix $H$, with
$N$ being the number of sites in the lattice or graph. Alternatively, one may
think of its continuum analogue, $H=-\frac{\hbar^{2}}{2m}\nabla^{2}+V({\bf
r}),\quad{\bf r}\in{\cal D}$, with the appropriate (e.g. Dirichlet) conditions
at the boundary of ${\cal D}$. In fact, under appropriate conditions the
essentially random nature of wave scattering can be generated by an
irregularly shaped boundary of the domain ${\cal D}$, without any intrinsic
potential disorder. This is the standard case in the so-called wave billiards,
the paradigmatic toy systems to study effects of quantum or wave chaos, see
e.g. [1, 2, 3, 4]. In such a case the famous Bohigas-Giannoni-Schmidt
conjecture [5] allows to describe universal features of such systems
efficiently by replacing the Hamiltonian $H$ with a random $N\times N$ matrix
from Gaussian Ensembles: Gaussian Orthogonal (GOE), Gaussian Symplectic (GSE),
or Gaussian Unitary (GUE), depending on the presence or absence of time-
reversal symmetry (and/or other relevant symmetries) in the system.

A very convenient framework for describing scattering of classical or quantum
waves from the disordered or chaotic medium has been formulated in Ref. [6],
see e.g. [7] for more detail. Within such a framework, which is frequently
called in the literature the “Heidelberg model”, one constructs the unitary
$M\times M$ energy-dependent scattering matrix $S(E)$ describing scattering of
waves incident on a random medium at some energy $E$ and then exiting it via
$M$ open scattering channels, numbered by $c=1,\ldots,M$, see a sketch below.
Unitarity reflects the flux conservation: the vectors ${\bf
a}=(a_{1},\ldots,a_{M})$ of incoming and ${\bf b}=(b_{1},\ldots,b_{M})$ of
outgoing amplitudes are linearly related via ${\bf b}=S(E)\,{\bf a}$ and have
the same norm.

.

Figure 1: A sketch of a chaotic wave scattering from a region schematically
represented by a cavity and assumed to contain a random medium inside. An
operator governing wave dynamics in such a cavity decoupled from the channels
is assumed to be effectively described by a large random matrix $H$. An
infinite lead is assumed to support $M$ propagating channels in the considered
energy range, and is coupled to the cavity region by a matrix/operator $W$.
The ensuing $M\times M$ unitary scattering matrix $S$ can be related to $H$
and $W$ in the framework of the Heidelberg approach, and is given by Eq. (2).

The relation between $S(E)$ and the medium Hamiltonian $H$ is then provided by
the following expression

$S(E)=\frac{{\bf 1}-i\mathbf{K}}{{\bf
1}+i\mathbf{K}},\qquad\mathbf{K}(E)=W^{\dagger}\frac{1}{E-H}W$ (2)

where columns $W_{c}$ ($c=1,\ldots,M$) of an $N\times M$ matrix $W$ of
coupling amplitudes to $M$ open scattering channels can be taken as fixed
vectors satisfying the orthogonality condition:

$\sum_{i=1}^{N}W^{*}_{ci}W_{bi}=\gamma_{c}\delta_{cb},$ (3)

with $\gamma_{c}>0$ $\forall c=1,\ldots,M$ determining the “bare” strength of
coupling of a given channel to the scattering system. The resulting $M\times
M$ Hermitian matrix $\mathbf{K}(E)$ is known in the literature as the Wigner
reaction $K-$ matrix. It is Hermiticity of $\mathbf{K}$ which implies
$S-$matrix unitarity, hence implies the flux conservation. Note that the
Wigner $K-$matrix is experimentally measurable in microwave scattering
systems, as it is directly related to the systems impedance matrix, see e.g.
[8, 9, 10].

One of serious challenges related to theoretical description of scattering
characteristics is however related to the fact that experimentally measured
quantities suffer from inevitable energy losses (absorption), e.g. due to
damping in resonator walls and other imperfections. Such losses violate
unitarity of the scattering matrix and are important for interpretation of
experiments, hence considerable efforts were directed towards incorporating
them into the Heidelberg approach [11]. At the level of the model (2) the
losses can be taken into account by allowing the spectral parameter $E$ to
have finite imaginary part by replacing $E\rightarrow E+i\alpha\in\mathbb{C}$
with some $\alpha>0$. This replacement violates Hermiticity of the Wigner
matrix $\mathbf{K}(E+i\alpha)$; in particular entries of $\mathbf{K}$ become
now complex even for real symmetric choice of $H$ and real $W$. The most
interesting, difficult and experimentally relevant regime occurs when
absorption parameter $\alpha$ is comparable with the mean separation
$\Delta(E)$ between neighbouring eigenvalues of the wave-chaotic Hamiltonian
$H$. For example, if one uses the Gaussian random matrix model for $H$,
normalized to have the mean eigenvalue density given by Wigner semicircle
$\nu(E)=1/(2\pi)\sqrt{4-E^{2}}$ in a finite interval $|E|<2$, one has
$\Delta(E)=(\nu(E)N)^{-1}$ as $N\to\infty$. The statistics of the real and
imaginary parts of $K-$matrix entries in such a regime was the subject of a
considerable number of theoretical papers [12, 13, 14, 15, 16] and by now
well-understood and measured experimentally with good precision for systems
with preserved time-reversal invariance in microwave cavities [17, 8, 9, 10]
and microwave simulations of quantum graphs [18, 19, 20, 21]. More recently
experimental results for K-matrices in systems with broken time-reversal
invariance [22, 23] and eventually symplectic symmetry [24] have been also
reported.

In the present paper we will be interested in yet another characteristics of
the non-Hermitian Wigner matrix $\mathbf{K}(E+i\alpha)$, the mean density of
its complex eigenvalues $K_{c}=\mathop{\rm Re}K_{c}-i\mathop{\rm Im}K_{c}$,
$\forall c=1,\ldots,M$ defined as

$\rho_{M}(u,v;y)=\left\langle\sum_{c=1}^{M}\delta\left(u-\mathop{\rm
Re}K_{c}\right)\delta\left(v-\mathop{\rm Im}K_{a}\right)\right\rangle,$ (4)

where we suppressed the energy dependence for simplicity, indicating instead
explicit dependence on the appropriately scaled absorption parameter
$y=\frac{2\pi\alpha}{\Delta}$. Here and henceforth the angular brackets
$\left\langle\ldots\right\rangle$ indicate the averaging over ensemble of
random Hamiltonians $H$. Note that by selecting the coupling vectors $W_{c}$
coinciding with the first $M$ basis vectors in $N-$dimensional space, i.e.
$W_{1}={\bf e}_{1}=(1,0,\ldots,0),\,W_{2}={\bf e}_{2}=(0,1,0,\ldots,0)$, etc.
converts the $K$-matrix to $M\times M$ top left corner of the $N\times N$
resolvent matrix $(E+i\alpha-H)^{-1}$. Physically, this corresponds to $M$
perfectly coupled channels, attached to first $M$ sites. From that angle we
aim to characterize the eigenvalue density for the corner resolvent minor at
complex values of the spectral parameter, an interesting and potentially rich
mathematical problem. We are not aware of any systematic studies in that
direction.

Note that in fully chaotic, zero-dimensional system the positions of channel
attachement do not play any role due to inherent ergodicity. In a more general
non-ergodic situation, which may arise due to presence of Anderson
localization phenomena, one may think of such an arrangement as corresponding
to a wave reflection problem. In such setting the density Eq.(4) has appeared
recently in the paper [25] as an important quantity facilitating computation
of the mean density of S-matrix poles, also known as resonances, in the
complex energy plane. The latter density is experimentally measurable in wave-
chaotic system [26, 27] and is a subject of long-standing theoretical
interest, see e.g. [28, 29, 30, 31, 32, 33, 34]. Clearly, the density Eq.(4)
is also experimentally measurable in principle, provided accurate experimental
data can be sampled for the whole $K-$matrix. The paper [25] included, without
a proper derivation, an explicit expression for such density, valid for a
general class disordered systems with broken time-reversal invariance, namely
for those which can be mapped on the so-called supersymmetric nonlinear
$\sigma-$model, see [35] and discussions in [25] for more information. The
present paper aims to fill in that gap by providing a detailed derivation,
which involves several steps only relatively sketchy described in the
available literature.

To begin with, for the special simplest case $M=1$ the $K-$matrix consists of
a single element, and finding the density Eq.(4) is equivalent to computing
the joint probability density of the real and complex parts of such an
element. Such density has been originally addressed in [36] via quite a
tediuos calculations in the $\sigma-$model approximation. A much more
efficient approach has been proposed later in [14], see also an account in
[11]. Our goal in this paper is to show how to generalize that approach to any
number of open channels $M$, on an example of systems with broken time-
reversal invariance. Along these lines we also try to elucidate some features
of the method which were omitted in the exposition of [14, 11].

##  2 Derivation of the main results

###  2.1 General exposition of the method

Given two real parameters $p\in\mathbb{R}$ and $q>0$ we start with defining
the following object:

$C_{\alpha}(p,q):=\left\langle
Tr\left(z-\mathbf{K}(E+i\alpha)\right)^{-1}Tr\left(\overline{z}-\mathbf{K}(E-i\alpha)\right)^{-1}\right\rangle,$
(5)

where we denoted $z=p+iq,\,\overline{z}=p-iq$ and assumed the real energy $E$
and the absorption parameter $\alpha>0$ to be fixed. As eigenvalues of the
matrices $\mathbf{K}(E+i\alpha)$ and $\mathbf{K}(E-i\alpha)$ are complex
conjugates of each other, one can write each trace in terms of
$K_{c}=\mathop{\rm Re}K_{c}-i\mathop{\rm Im}K_{c}:=u_{c}-iv_{c}$, with
$v_{c}>0$, representing Eq.(5) as a sum of diagonal and off-diagonal
contributions:

$C_{\alpha}(p,q)=C^{(diag)}_{\alpha}(p,q)+C^{(off)}_{\alpha}(p,q)$

where

$C^{(diag)}_{\alpha}(p,q):=\left\langle\sum_{c=1}^{M}\frac{1}{|z-K_{c}|^{2}}\right\rangle=\left\langle\sum_{c=1}^{M}\frac{1}{(p-u_{c})^{2}+(q+v_{c})^{2}}\right\rangle,$
(6)

and

$C^{(off)}_{\alpha}(p,q):=\left\langle\sum_{c\neq
c^{\prime}}^{M}\frac{1}{(z-K_{c})(\overline{z}-\overline{K_{c^{\prime}}}}\right\rangle$
(7) $=\left\langle\sum_{c\neq
c^{\prime}}^{M}\frac{1}{u_{c}-u_{c^{\prime}}-i(2q+v_{c}+v_{c^{\prime}})}\left(\frac{1}{p-u_{c}+i(q+v_{c})}-\frac{1}{p-u_{c^{\prime}}-i(q+v_{c^{\prime}})}\right)\right\rangle$

At the next step let us introduce the Fourier-transform in variable $p$:

$\tilde{C}_{\alpha}(k,q):=\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{ipk}C_{\alpha}(p,q)\,dp.$

Taking into account $q>0,v_{c}\geq 0,\forall c$ we get

$\tilde{C}^{(diag)}_{\alpha}(k,q)=\left\langle\frac{1}{2}\sum_{c=1}^{M}\frac{e^{iku_{c}-|k|(q+v_{c})}}{q+v_{c}}\right\rangle$
(8)

whereas the Fourier transformed off-diagonal part reads as

$\tilde{C}^{(off)}_{\alpha}(k,q)=\left\langle\sum_{c\neq
c^{\prime}}^{M}\frac{-i}{u_{c}-u_{c^{\prime}}-i(2q+v_{c}+v_{c^{\prime}})}\right.$
(9)
$\left.\times\left(\theta(-k)e^{iku_{c}+k(q+v_{c})}+\theta(k)e^{iku_{c^{\prime}}-k(q+v_{c^{\prime}})}\right)\right\rangle$

where $\theta(k)=1$ for $k\geq 0$ and zero otherwise.

The next step is to continue analytically in the parameter $q$ from positive
real values to the whole complex plane slit along the negative real line:
$q=-v,v>0$, and evaluate the jump across the slit, defined as

$\delta\tilde{C}_{\alpha}(k,v>0):=\lim_{\epsilon\to
0}\left(\tilde{C}_{\alpha}(k,-v-i\epsilon)-\tilde{C}_{\alpha}(k,-v+i\epsilon)\right)$
(10)

For the diagonal part one finds after straightforward algebra

$\delta\tilde{C}^{(diag)}_{\alpha}(k,v>0)=i\lim_{\epsilon\to
0}\left\langle\sum_{c=1}^{M}e^{iku_{c}+|k|(v-v_{c})}\left(\frac{\epsilon\cos(\epsilon|k|)}{\epsilon^{2}+(v-v_{c})^{2}}-\frac{\sin(\epsilon|k|)(v-v_{c})}{\epsilon^{2}+(v-v_{c})^{2}}\right)\right\rangle$
(11)

which upon using

$\lim_{\epsilon\to
0}\left(\frac{\epsilon\cos(\epsilon|k|)}{\epsilon^{2}+(v-v_{c})^{2}}-\frac{\sin(\epsilon|k|)(v-v_{c})}{\epsilon^{2}+(v-v_{c})^{2}}\right)=\pi\delta(v-v_{c})$

reduces the diagonal contribution to

$\delta\tilde{C}^{(diag)}_{\alpha}(k,v>0)=i\pi\left\langle\sum_{c=1}^{M}e^{iku_{c}}\,\delta(v-v_{c})\right\rangle$
(12)

At the same time straightforward computations show that assuming that the
eigenvalues of the $K-$matrix are all distinct, i.e. $u_{c}-iv_{c}\neq
u_{c}^{\prime}-iv_{c^{\prime}}$ for $c\neq c^{\prime}$, the off-diagonal part
does not generate any nonvanishing jump across the slit at $q=-v,v>0$, that is
$\delta\tilde{C}^{(off)}_{\alpha}(k,v>0)=0$. Finally, applying in Eq.(13) the
inverse Fourier transform in the variable $k$ and comparing with the
definition Eq.(4) provides the expression for the density of complex
eigenvalues of the $K-$matrix in the form

$\rho_{M}(u,v;y)=\frac{1}{2i\pi^{2}}\int_{-\infty}^{\infty}e^{-iku}\delta\tilde{C}_{\alpha=y\Delta/2\pi}(k,v>0)\,dk.$
(13)

In this way the problem of computing the density $\rho_{M}(u,v;y)$ is reduced
to ability to evaluate explicitly the correlation function $C_{\alpha}(p,q>0)$
in Eq.(5) and perform the required Fourier transforms and jump evaluation.
Below we show how this program is executed for those disordered or chaotic
systems with broken time-reversal invariance which can be mapped onto the
corresponding nonlinear $\sigma-$model.

###  2.2 Computations for systems with broken time-reversal invariance

Referring the interested reader to [25] and references therein for a detailed
discussion of physical assumptions behind such mapping, we just mention here
that it provides the most powerful and systematic approaches to addressing
universal single particle features of wave propagation in a disordered medium,
including Anderson localization phenomena. Developed in the seminal works by
Efetov [35] building on earlier ideas of Wegner [37] the model is defined by
specifying a weight function $e^{-{\cal S}[Q]}$, with the action ${\cal S}[Q]$
describing interaction between supermatrices $Q({\bf r})$ (i.e. matrices with
Grassmann/anticommuting/ fermionic and ordinary/commuting/bosonic entries)
associated to every site ${\bf r}\in\tilde{\mathbf{\Lambda}}$ located on an
auxiliary lattice $\tilde{\mathbf{\Lambda}}$. The size of supermatrices
involved depends on the underlying symmetries of the Hamiltonian $H$, and in
the simplest case of the Hamiltonians with fully broken time-reversal
symmetry, denoted in the standard nomenclature as class A with Dyson parameter
$\beta=2$, the supermatrices are of the size $4\times 4$. Physically such
model provides, in a certain sense, a coarse-grained description of the
original microscopic Anderson model or its continuous equivalent, with non-
universal features on scales smaller than the mean-free path $l$ being
effectively integrated out. In such a picture every (super)matrix $Q({\bf r})$
associated to a single lattice site in $\tilde{\mathbf{\Lambda}}$ ”lumps
together” behaviour of the microscopic model on scales of the order of the
mean-free path $l$. From this point of view the billiards in the quantum
chaotic regime, where essentially $l$ is of the same order as the system
length $L$, are effectively characterized by nonlinear $\sigma-$models with a
single matrix $Q$ without any spatial dependence. Such limit is traditionally
called ”zero-dimensional”. At the same time all effects of the Anderson
localization require considering extended lattices of interacting
$Q-$matrices.

One of the central objects of such theory turns out to be the so called “order
parameter function” (OPF) $F_{\bf r}(Q)$ which is formally defined [38] by
integrating the weight $e^{-{\cal S}[Q]}$ over all but one supermatrix $Q({\bf
r})$. Due to global symmetries of the action, the OPF can be shown to actually
depend only on a few real Cartan variables parametrizing $Q$ matrices. In
particular, for systems with broken time-reversal symmetry one has $F_{\bf
r}(Q):={\cal F}(\lambda,\lambda_{1})$, with $\lambda\in[-1,1]$ and
$\lambda_{1}\in[1,\infty]$ being the compact and non-compact coordinates,
respectively (we omitted spatial dependence on ${\bf r}$ for brevity). Note
that the OPF characterizes the closed system which (in the absence of
absorption) conserves the number of particles, whereas allowing
particles/waves at a given energy to be sent via the lead to the random medium
and then collecting the reflected waves renders the medium open. However, if
one makes an assumption of “locality” of the lead, whose transverse extent is
assumed to be much smaller than the mean-free path $l$ in the disordered
medium, makes the coupling to it effectively point-wise at the level of
$\sigma$-model description. Still, even such point-wise lead may support
arbitrary many propagation channels $M$, though we will be always assuming $M$
remaining negligible to the number of sites in the underlying microscopic
lattice $\mathbf{\Lambda}$.

The power of nonlinear $\sigma$-model description in our case lies in our
ability to provide an explicit representation for the correlation function
$C_{\alpha}(p,q)$ defined in Eq.(5) in terms of the OPF ${\cal
F}(\lambda,\lambda_{1})$ at the point of lead attachment. For systems with
broken time-reversal invariance such computation has been already performed in
[7], albeit formally only in the ”zero-dimensional” limit, with OPF taking an
especially simple form ${\cal
F}(\lambda,\lambda_{1})=e^{-y(\lambda_{1}-\lambda)}$, where as before
$y=2\pi\alpha/\Delta$ is the effective absorption parameter. It is however
straightforward to adapt the calculation for arbitrary nonlinear sigma-model,
see Appendix B of [11], the result being given by the sum of two
contributions, the disconnected one

$C_{\alpha}^{(disc)}(p,q)={\small\sum_{c=1}^{M}\frac{1}{p-\gamma_{c}\frac{E}{2}-i(q+\pi\nu(E)\gamma_{c})}\sum_{b=1}^{M}\frac{1}{p-\gamma_{c}\frac{E}{2}+i(q+\pi\nu(E)\gamma_{b})}}$
(14)

and the connected one

$C_{\alpha}^{(con)}(p,q)=\int_{-1}^{1}d\lambda\int_{1}^{\infty}d\lambda_{1}\frac{{\cal
F}(\lambda,\lambda_{1})}{(\lambda_{1}-\lambda)^{2}}\,{\cal
R}_{M}(p,q|\lambda,\lambda_{1})$ (15)

where the last factor in Eq.(15) is given by

${\cal R}_{M}(p,q|\lambda,\lambda_{1}):=\mathbf{\cal
L}_{p,q}\prod_{c=1}^{M}\frac{(p-\gamma_{c}\frac{E}{2})^{2}+q^{2}+2\pi\nu(E)\gamma_{c}\lambda+(\pi\nu(E)\gamma_{c})^{2}}{(p-\gamma_{c}\frac{E}{2})^{2}+q^{2}+2\pi\nu(E)\gamma_{c}\lambda_{1}+(\pi\nu(E)\gamma_{c})^{2}}$
(16)

with the coupling coefficients $\gamma_{c}$ defined in Eq.(3) and the
differential operator $\mathbf{\cal
L}_{p,q}:=\frac{1}{4}\left(\frac{\partial^{2}}{\partial
p^{2}}+\frac{\partial^{2}}{\partial q^{2}}\right)$.

These expressions provide the basis for implementing the analytic continuation
procedure described above. For simplicity we consider below explicitly only
the case $E=0$, so that $\pi\nu(E)=1$, and largely concentrate on the
simplest, yet important case of equivalent channels:
$\gamma_{c}=\gamma,\,\forall c=1,\ldots,M$ (see however Eq.(32) for two non-
equivalent channels). The analytic continuation procedure for the disconnected
part amounts to a straightforward repetition of our derivation of Eq.(13) and
yields $\rho^{(disc)}(u,v)=\sum_{c=1}^{M}\delta(u)\delta(v-\gamma_{c})$. The
connected contribution to the density is much less trivial and we consider it
below.

One starts with rewriting Eq.(16) in the form

${\cal R}_{M}(p,q|\lambda,\lambda_{1})=\mathbf{\cal
L}_{p,q}\left(1-\frac{2q\gamma(\lambda_{1}-\lambda)}{p^{2}+q^{2}+2\gamma_{c}\lambda_{1}+\gamma_{c}^{2}}\right)^{M},$
(17)

which after expanding the binomial reduces to

${\cal
R}_{M}(p,q|\lambda,\lambda_{1})=-\sum_{l=1}^{M}\left(\begin{array}[]{c}M\\\
l\end{array}\right)\frac{(\lambda_{1}-\lambda)^{l}}{(l-1)!}\frac{\partial^{l-1}}{\partial\lambda_{1}^{l-1}}\mathbf{\cal
L}_{p,q}\frac{2q\gamma}{p^{2}+q^{2}+2\gamma_{c}\lambda_{1}+\gamma_{c}^{2}}.$
(18)

The latter form makes it an easy task to perform the Fourier transform in the
variable $p$ assuming $q>0$, which essentially amounts to making in Eq.(18)
the replacement

$\mathbf{\cal
L}_{p,q}\frac{2q\gamma}{p^{2}+q^{2}+2\gamma_{c}\lambda_{1}+\gamma_{c}^{2}}\,\,\longrightarrow\,\,\phi(k,q)=\frac{\pi\gamma}{2}\left(\frac{\partial^{2}}{\partial
q^{2}}-k^{2}\right)\,q\frac{e^{-|k|\sqrt{q^{2}+2\gamma\lambda_{1}q+\gamma^{2}}}}{\sqrt{q^{2}+2\gamma\lambda_{1}q+\gamma^{2}}}$

Following the procedures described in Eq(10) we now continue analytically in
the parameter $q$ from positive real values to the whole complex plane slit
along the negative real line: $q=-v,v>0$, and evaluate the associated jump
across the slit

$\delta\phi(k,v>0):=\lim_{\epsilon\to
0}\left(\phi(k,-v-i\epsilon)-\phi(k,-v+i\epsilon)\right)$ (19)

which is easily found to be equal to

${\small\delta\phi(k,v>0)=\pi\gamma\left(\frac{\partial^{2}}{\partial
v^{2}}-k^{2}\right)v\frac{\cos{k\sqrt{2\gamma\lambda_{1}v-v^{2}-\gamma^{2}}}}{\sqrt{2\gamma\lambda_{1}v-v^{2}-\gamma^{2}}}\,\theta(2\gamma\lambda_{1}v-v^{2}-\gamma^{2}).}$
(20)

Straightforward inversion of the Fourier-transform in the variable $k$
converts the above into

$\delta\phi(u,v>0)=\frac{\gamma}{2}\left(\frac{\partial^{2}}{\partial
u^{2}}+\frac{\partial^{2}}{\partial v^{2}}\right)v$
$\times\frac{\delta\left(u-\sqrt{2\gamma\lambda_{1}v-v^{2}-\gamma^{2}}\right)+\delta\left(u+\sqrt{2\gamma\lambda_{1}v-v^{2}-\gamma^{2}}\right)}{\sqrt{2\gamma\lambda_{1}v-v^{2}-\gamma^{2}}}\,\theta(2\gamma\lambda_{1}v-v^{2}-\gamma^{2}).$
$=\frac{1}{2}\left(\frac{\partial^{2}}{\partial
u^{2}}+\frac{\partial^{2}}{\partial
v^{2}}\right)\delta\left(\lambda_{1}-x_{\gamma}\right),\quad
x_{\gamma}:=\frac{u^{2}+v^{2}+\gamma^{2}}{2\gamma\,v}.$ (21)

Next we trade the derivatives over $\lambda_{1}$ for those over $x_{\gamma}$
by the identity

$\frac{\partial^{l-1}}{\partial\lambda_{1}^{l-1}}\delta\left(\lambda_{1}-x_{\gamma}\right)=(-1)^{l-1}\frac{\partial^{l-1}}{\partial
x_{\gamma}^{l-1}}\delta\left(\lambda_{1}-x_{\gamma}\right)$

and in this way arrive at replacing Eq.(18) with

$\tilde{{\cal
R}}_{M}(u,v|\lambda,\lambda_{1})=-\frac{1}{2}\left(\frac{\partial^{2}}{\partial
u^{2}}+\frac{\partial^{2}}{\partial
v^{2}}\right)\sum_{l=1}^{M}\left(\begin{array}[]{c}M\\\
l\end{array}\right)\frac{(-1)^{l-1}(\lambda_{1}-\lambda)^{l}}{(l-1)!}\frac{\partial^{l-1}}{\partial
x_{\gamma}^{l-1}}\delta\left(\lambda_{1}-x_{\gamma}\right).$ (22)

With this Eqs. (13) and Eq.(15) imply the density of $K-$matrix eigenvalues
via

$\rho^{(con)}(u,v)=\frac{1}{2\pi}\int_{-1}^{1}d\lambda\int_{1}^{\infty}d\lambda_{1}\frac{{\cal
F}(\lambda,\lambda_{1})}{(\lambda_{1}-\lambda)^{2}}\,\tilde{{\cal
R}}_{M}(u,v|\lambda,\lambda_{1})$ (23)

which upon substituting (22) into it and changing the order of integrations
yields

$\rho^{(con)}(u,v)=\frac{1}{4\pi}\left(\frac{\partial^{2}}{\partial
u^{2}}+\frac{\partial^{2}}{\partial
v^{2}}\right)\int_{-1}^{1}d\lambda\,G_{M}(\lambda|x_{\gamma}).$ (24)

Here we denoted

$G_{M}(\lambda|x_{\gamma}):=\sum_{l=1}^{M}\left(\begin{array}[]{c}M\\\
l\end{array}\right)\frac{(-1)^{l-1}}{(l-1)!}\frac{\partial^{l-1}}{\partial
x_{\gamma}^{l-1}}\left[(x_{\gamma}-\lambda)^{l}T(\lambda|x_{\gamma})\right],$
(25)

with

$T(\lambda|x_{\gamma})=\frac{{\cal
F}(\lambda,x_{\gamma})}{(x_{\gamma}-\lambda)^{2}}.$ (26)

Applying the Leibnitz formula

$\frac{\partial^{l-1}}{\partial
x_{\gamma}^{l-1}}\left[(x_{\gamma}-\lambda)^{l}T(\lambda|x_{\gamma})\right]=\sum_{k=0}^{l-1}\left(\begin{array}[]{c}l-1\\\
k\end{array}\right)\frac{l!}{(k+1)!}(x_{\gamma}-\lambda)^{k+1}\frac{\partial^{k}}{\partial
x_{\gamma}^{k}}T(\lambda|x_{\gamma})$

and substituting it back to Eq.(25) one may change the order of summation as

$G_{M}(\lambda|x_{\gamma})=\sum_{l=1}^{M}A_{l}\sum_{k=0}^{l-1}B_{k,l}V_{k}=\sum_{k=0}^{M-1}V_{k}\sum_{l=k+1}^{M}A_{l}B_{k,l}$

with $V_{k}:=(x_{\gamma}-\lambda)^{k+1}\frac{\partial^{k}}{\partial
x_{\gamma}^{k}}T(\lambda|x_{\gamma})$ and

$A_{l}:=\left(\begin{array}[]{c}M\\\
l\end{array}\right)\frac{(-1)^{l-1}}{(l-1)!},\quad
B_{k,l}:=\left(\begin{array}[]{c}l-1\\\ k\end{array}\right)\frac{l!}{(k+1)!}.$

This gives

$\sum_{l=k+1}^{M}A_{l}B_{k,l}=\frac{M!}{k!(k+1)!}\sum_{l=k+1}^{M}\frac{(-1)^{l-1}}{(M-l)!}\frac{1}{(l-k-1)!}$
$=\frac{(-1)^{k}\,M!}{(M-k-1)!k!(k+1)!}\sum_{n=0}^{M-k-1}(-1)^{n}\left(\begin{array}[]{c}M-k-1\\\
n\end{array}\right)=\frac{(-1)^{M-1}}{(M-1)!}\delta_{k,M-1}$

using the Kronecker symbol $\delta_{k,k^{\prime}}$, since the sum over $n$ is
vanishing for all $0\leq k<M-1$, and is equal to unity at $k=M-1$.

As the result we get the final expression for the connected part of the mean
density of $K-$matrix eigenvalues in the form

$\rho^{(con)}_{M}(u,v)=\frac{1}{4\pi}\frac{(-1)^{M-1}}{(M-1)!}\left(\frac{\partial^{2}}{\partial
u^{2}}+\frac{\partial^{2}}{\partial
v^{2}}\right)\int_{-1}^{1}d\lambda\,(x_{\gamma}-\lambda)^{M}\,\frac{\partial^{M-1}}{\partial
x_{\gamma}^{M-1}}\frac{{\cal
F}(\lambda,x_{\gamma})}{(x_{\gamma}-\lambda)^{2}}$ (27)

A few remarks are here in order which help to properly interpret and
appreciate the content of Eq.(27).
Remark 1. Recalling from Eq.(21) that

$x_{\gamma}=\frac{u^{2}+v^{2}+\gamma^{2}}{2\gamma\,v}\equiv\frac{u^{2}}{2\gamma\,v}+\frac{1}{2}\left(\frac{v}{\gamma}+\frac{\gamma}{v}\right)\geq
1$ (28)

one may straightforwardly check that for any smooth enough function $\Phi(x)$
holds

$\left(\frac{\partial^{2}}{\partial u^{2}}+\frac{\partial^{2}}{\partial
v^{2}}\right)\Phi(x_{\gamma})=\frac{1}{v^{2}}\frac{\partial}{\partial
x_{\gamma}}(x_{\gamma}^{2}-1)\frac{\partial}{\partial
x_{\gamma}}\Phi(x_{\gamma}),\quad x_{\gamma}>1.$ (29)

This was exactly the form used to represent the density in [25].
There is however a subtlety in Eq.(27) related with its content at
$x_{\gamma}\to 1$. In our derivation we tacitly assumed $x_{\gamma}>1$.
However, a more careful analysis shows that the integral in the right-hand
side of Eq.(27) should be pre-multiplied with the step-function factor
$\theta(x_{\gamma}-1)$ arising as the result of performing integration over
$\lambda_{1}\in[1,\infty)$ with the factor $\delta(\lambda_{1}-x_{\gamma})$.
Presence of such a seemingly innocent $\theta-$factor has however important
consequences: when acted upon with the differential operator in the right-hand
side of Eq.(29) it generates the $\delta$-function factors exactly cancelling
the contribution from the disconnected part,
$\rho^{(disc)}(u,v)=\sum_{c=1}^{M}\delta(u)\delta(v-\gamma_{c})$. As a result,
the formula Eq.(27) as it is written (i.e. without $\theta-$factor) in fact
gives the full, properly normalized, eigenvalue density for the $K-$matrix in
absorptive systems. A similar mechanism of cancellation of $\delta-$terms has
been first noticed in [39], and we explain in the Appendix A how it works in
our case using the simplest case of $M=1$ as an example.

Remark 2. With the hindsight, one may notice that one could have arrived to
the same expression Eq.(27) by a much simpler procedure. Namely, by defining

$\tilde{x}:=\frac{p^{2}+q^{2}+\gamma^{2}}{2\gamma\,q}$ (30)

rewrite Eq.(17) in the form

${\cal
R}_{M}(p,q|\lambda,\lambda_{1})=\frac{1}{4q^{2}}\frac{\partial}{\partial\tilde{x}}(\tilde{x}^{2}-1)\frac{\partial}{\partial\tilde{x}}\left(\frac{\tilde{x}+\lambda}{\tilde{x}+\lambda_{1}}\right)^{M},$
(31)

Then simply replace $u\to p,q\to-v-i0$ implying $\tilde{x}\to-x_{\gamma}+i0$
and calculate the associated jump across the cut using

$\mathop{\rm
Im}\left(\frac{-x_{\gamma}+i0+\lambda}{-x_{\gamma}+i0+\lambda_{1}}\right)^{M}=(x_{\gamma}-\lambda)^{M}\frac{(-1)^{M-1}}{(M-1)!}\frac{\partial^{M-1}}{\partial
x_{\gamma}^{M-1}}\mathop{\rm Im}\frac{1}{x_{\gamma}-\lambda_{1}-i0}$
$=\pi(x_{\gamma}-\lambda)^{M}\frac{(-1)^{M-1}}{(M-1)!}\frac{\partial^{M-1}}{\partial
x_{\gamma}^{M-1}}\delta(x_{\gamma}-\lambda_{1})$

Such recipe was exactly one employed for $M=1$ in [14], though without a
proper explanation provided there or in the review [11].

Armed with such a recipe, one can easily apply it to the case of non-
equivalent channels. General formulas look in that case quite complicated, but
in the simplest case of two non-equivalent channels with coupling constants
$\gamma_{1}\neq\gamma_{2}$ one gets a relatively compact expression:

$\rho_{\gamma_{1},\gamma_{2}}(u,v)=\frac{1}{4\pi}\left(\frac{\partial^{2}}{\partial
u^{2}}+\frac{\partial^{2}}{\partial
v^{2}}\right)\left\\{-\int_{-1}^{1}\,d\lambda\,\frac{{\cal
F}(\lambda,x_{1})-{\cal F}(\lambda,x_{2})}{x_{1}-x_{2}}\right.$ (32)
$\left.+\int_{-1}^{1}\,d\lambda\left[\frac{{\cal
F}(\lambda,x_{1})}{x_{2}-\lambda_{2}}+\frac{{\cal
F}(\lambda,x_{2})}{x_{1}-\lambda_{1}}\right]\right\\},$

where we defined

$x_{1}=\frac{u^{2}+v^{2}+\gamma_{1}^{2}}{2\gamma_{1}\,v},\quad
x_{2}=\frac{u^{2}+v^{2}+\gamma_{2}^{2}}{2\gamma_{2}\,v}.$ (33)

Remark 3. It is clear that performing further analysis of Eq.(27) hinges on
our ability to have a good understanding of the OPF ${\cal F}(\lambda,x)$ for
the closed counterpart of the scattering system, which in general also depends
on the (appropriately normalized) absorption parameter $\alpha$. Such
knowledge is currently available mainly in two cases (i) the ”zero-
dimensional” limit, with OPF taking an especially simple form ${\cal
F}^{(0d)}(\lambda,x)=e^{-y(x-\lambda)}$, where as before $y=2\pi\alpha/\Delta$
and (ii) in a (semi) infinite quasi-one dimensional wire, see the sketch
below, of length $L\to\infty$, with one edge closed for the waves and second
edge attached to an infinite waveguide with $M$ propagating channels.

Figure 2: A sketch of the “quasi-1D” model. The left part in grey represents
an infinite-length ideal lead supporting $M$ propagating modes. The disordered
part is of a finite length $L$ and contains finite concentration of random
impurities inside.

Such wire is characterized by a classical microscopic diffusion constant $D$
related to the localization length $\xi$ of quantum wave problem as
$\xi=2\pi\nu D$, with $\nu$ being as before the mean eigenvalue density at a
given energy. Note that mathematically such wires can be modelled by a large
banded random matrix [40, 41]. In such a system the OPF at points close to its
edges has been originally found in [42] and takes the following form in terms
of the modified Bessel functions $I_{p}(z),K_{p}(z)$:

${\cal F}^{(1d)}(\lambda,x)=K_{0}(a)bI_{1}(b)+aK_{1}(a)I_{0}(b),$ (34)

with

$a=\kappa\sqrt{(x+1)/2},\qquad b=\kappa\sqrt{(\lambda+1)/2},$ (35)

where the parameter $\kappa$ is related to the absorption $\alpha$ as

$\kappa=\sqrt{8\alpha/\Delta_{\xi}},$ (36)

with an important energy scale
$\Delta_{\xi}=(4\pi^{2}D\nu^{2})^{-1}=D/\xi^{2}$ giving the mean level spacing
in the quasi-one dimensional wires whose length $L$ equal to the localization
length $\xi$.

In the ”zero-dimensional” limit, due to a simple form of the OPF one can
relatively straightforwardly perform the required integrations and
differentiations in Eq.(27) and get the explicit formulas, which we present
below for the simplest cases $M=1$ and $M=2$ of equivalent channels:

$\rho_{0d,M=1}(u,v)=\frac{1}{2\pi
v^{2}}e^{-x_{\gamma}}\left(y\cosh{y}-\sinh{y}(1-yx_{\gamma})\right)$ (37)

and

$\rho_{0d,M=2}(u,v)=\frac{1}{2\pi
v^{2}}e^{-x_{\gamma}}\left(\sinh{y}[y(x^{2}_{\gamma}-1)-2x_{\gamma}]+2[y\cosh{y}-\sinh(1-yx_{\gamma})]\right).$
(38)

with the same definition of $x_{\gamma}$, Eq.(28). The formula equivalent to
Eq.(37) appeared already in the literature, see Eq.(5) in [13], the two
channel case seems to be new. As to the quasi-$1D$ system of infinite length,
it turns out that again the results can be found explicitly in the general
case. Below we present it only for the simplest case of a single attached
channel, when the density aquires quite an elegant form after manipulations
outlined in the Appendix B to this paper:

$\rho_{1d,M=1}(u,v)=\frac{1}{2\pi v^{2}}P_{0}(x_{\gamma}),\quad
P_{0}(x)=\frac{\kappa^{2}}{4}\left[I_{2}(\kappa)K_{0}\left(\kappa\sqrt{\frac{x+1}{2}}\right)\right.$
(39)
$\left.+I_{1}(\kappa)\sqrt{\frac{x+1}{2}}K_{1}\left(\kappa\sqrt{\frac{x+1}{2}}\right)\right]$

As is shown in [13], for $M=1$ and $\gamma=1$ the variable $r=(x-1)/(x+1)$ is
nothing else but the modulus of the reflection coefficient, which in the
absorptive system is smaller than one. Correspondingly, the function
$P_{0}(x)$ in Eq.(39) provides the distribution for $x$, hence for $r$, in a
single-channel quasi-$1D$ system with absorption. This complements a result
for the same geometry in the case of no absorption inside the sample, but for
the second edge of the sample being in contact with perfectly absorbing lead,
see eqs. (12)-(13) in [12]. Note also that it is not difficult to integrate
further out the variable $u$, getting an explicit formula for the distribution
of variable $v$, known as the local density of states, corresponding to
locations close to the sample edge. The latter is an important characteristic
of disordered single-particle systems, see [36, 43, 44].

In conclusion, we derived the mean density of complex eigenvalues for random
Wigner reaction $K-$matrices for absorptive disordered or chaotic systems with
broken time-reversal invariance, in the sigma-model approximation. Extension
of these results to systems with preserved time-reversal invariance (and then
eventually symplectic symmetry) is certainly possible along similar lines,
generalizing $M=1$ results presented in [11]. These subjects are left for
future publications.

Acknowledgments: This research has been supported by the EPSRC Grant
EP/V002473/1 “Random Hessians and Jacobians: theory and applications”.

Appendix A: cancellation of the disconnected part.
Our starting point is the formula Eq.(27) with included $\theta-$factor,
specified for simplicity and transparency to the case of a single channel
$M=1$ and $\gamma=1$, so that $x_{\gamma=1}=x$. We write it in the form

$\rho^{(con)}_{M=1}(u,v)=\frac{1}{4\pi v^{2}}{\cal
L}_{x}\left[\theta(x-1)\Phi(x)\right],\quad\Phi(x)=\int_{-1}^{1}\,\frac{{\cal
F}(\lambda,x)}{(x-\lambda)}\,d\lambda,$ (40)

where we introduced the differential operator

${\cal L}_{x}:=\frac{\partial}{\partial x}(x^{2}-1)\frac{\partial}{\partial
x}.$ (41)

Straightforward differentiation then gives

${\cal
L}_{x}\left[\theta(x-1)\Phi(x)\right]=\delta(x-1)\left[2x\Phi(x)+2(x^{2}-1)\Phi^{\prime}(x)\right]$
(42) $+\delta^{\prime}(x-1)\left[(x^{2}-1)\Phi(x)\right]+\theta(x-1){\cal
L}_{x}\Phi(x).$

Further using the integration by parts identity

$\delta^{\prime}(x-1)\left[(x^{2}-1)\Phi(x)\right]=-\delta(x-1)\frac{d}{dx}\left[(x^{2}-1)\Phi(x)\right]=-\delta(x-1)\left[2x\Phi(x)+(x^{2}-1)\Phi^{\prime}(x)\right]$

we conclude that

${\cal
L}_{x}\left[\theta(x-1)\Phi(x)\right]=\delta(x-1)\left[(x^{2}-1)\Phi^{\prime}(x)\right]+\theta(x-1){\cal
L}_{x}\Phi(x)$ (43)

so it remains to evaluate $\lim_{x\to
1}\left[(x^{2}-1)\Phi^{\prime}(x)\right]$. To this end we notice that it can
be generally shown that $\lim_{x\to 1}{\cal F}(\lambda,x)=1$, hence from
Eq.(40) we have $\Phi(x\to
1)\approx\int_{-1}^{1}\,\frac{1}{(x-\lambda)}\,d\lambda=\ln{\frac{x+1}{x-1}}$,
which immediately implies $\lim_{x\to
1}\left[(x^{2}-1)\Phi^{\prime}(x)\right]=-2$. This gives the singular
contribution to the density Eq.(40) in terms of the variables $u,v$ given by

$-\frac{2}{4\pi
v^{2}}\delta\left(\frac{u^{2}+v^{2}+1}{2v}-1\right)=-\frac{1}{\pi
v}\delta\left(u^{2}+(v-1)^{2}\right)=-\delta(u)\delta(v-1)$

which exactly cancels the contribution from the disconnected part.

Appendix B
In this appendix we show how Eq.(34) when substituted to Eq.(27) implies
Eq.(39). Throughout this appendix we again use $x_{\gamma}=x$ and ${\cal
L}_{x}:=\frac{\partial}{\partial x}(x^{2}-1)\frac{\partial}{\partial x}$.
First of all, we use the identity (43) from the paper [25], which claims that

$\frac{\partial}{\partial\kappa}{\cal
F}^{(1d)}(\lambda,x)=-\frac{\kappa}{2}(x-\lambda)K_{0}\left(\kappa\sqrt{\frac{x+1}{2}}\right)I_{0}\left(\kappa\sqrt{\frac{\lambda+1}{2}}\right)$
(44)

By differentiating both sides of Eq.(27) over $\kappa$ and using Eq.(44) in
the right-hand side yields

$\frac{\partial}{\partial\kappa}\rho_{1d,M=1}(u,v)=-\frac{1}{8\pi v^{2}}{\cal
L}_{x}\left[\kappa\,K_{0}\left(\kappa\sqrt{\frac{x+1}{2}}\right)\int_{-1}^{1}I_{0}\left(\kappa\sqrt{\frac{\lambda+1}{2}}\right)\,d\lambda\right]$
(45)

and after performing the integral by substitution $\lambda=2z^{2}-1,z\in[0,1]$
find that

$\frac{\partial}{\partial\kappa}\rho^{(con)}_{1d,M=1}(u,v)=-\frac{1}{2\pi
v^{2}}{\cal
L}_{x}\left[K_{0}\left(\kappa\sqrt{\frac{x+1}{2}}\right)I_{1}\left(\kappa\right)\right]$
(46) $=-\frac{1}{2\pi v^{2}}\frac{\partial}{\partial
x}\sqrt{\frac{x+1}{2}}\left[\frac{1-x}{2}\,\kappa
I_{1}\left(\kappa\right)\,K_{1}\left(\kappa\sqrt{\frac{x+1}{2}}\right)\right]$
(47)

At the next step we employ the following identity (c.f. 5.54 in p.624 of
[45]):

$\frac{1-x}{2}\,\kappa
I_{1}\left(\kappa\right)\,K_{1}\left(\kappa\sqrt{\frac{x+1}{2}}\right)=\frac{\partial}{\partial\kappa}\left[\kappa
I_{2}\left(\kappa\right)K_{1}\left(\kappa\sqrt{\frac{x+1}{2}}\right)\right.$
(48) $\left.+\sqrt{\frac{x+1}{2}}\kappa
I_{1}\left(\kappa\right)K_{2}\left(\kappa\sqrt{\frac{x+1}{2}}\right)\right]$

Using the fact that $\rho^{(con)}_{1d,M=1}(u,v)\to 0$ as $\kappa\to\infty$ we
then may conclude that Eq.(47) and Eq.(48) together imply

$\rho_{1d,M=1}(u,v)=-\frac{1}{2\pi v^{2}}\frac{\partial}{\partial
x}\sqrt{\frac{x+1}{2}}\left\\{\kappa
I_{2}\left(\kappa\right)K_{1}\left(\kappa\sqrt{\frac{x+1}{2}}\right)+\sqrt{\frac{x+1}{2}}\kappa
I_{1}\left(\kappa\right)K_{2}\left(\kappa\sqrt{\frac{x+1}{2}}\right)\right\\}$
$=-\frac{1}{2\pi
v^{2}}\left\\{\frac{I_{1}\left(\kappa\right)}{\kappa}\frac{\partial}{\partial
x}\left[\kappa^{2}\frac{x+1}{2}K_{2}\left(\kappa\sqrt{\frac{x+1}{2}}\right)\right]+I_{2}\left(\kappa\right)\frac{\partial}{\partial
x}\left[\kappa\sqrt{\frac{x+1}{2}}K_{1}\left(\kappa\sqrt{\frac{x+1}{2}}\right)\right]\right\\}$

Finally introducing in the above the variable $z=\kappa\sqrt{\frac{x+1}{2}}$,
using the chain rule and the identity (see 8.846.14 in [45])

$\frac{d}{dz}\left(z^{p}K_{p}(z)\right)=-z^{p}K_{p-1}(z)$

allows to bring the density $\rho^{(con)}_{1d,M=1}(u,v)$ to the final form
Eq.(39).

## References

  * [1] U. Kuhl, O. Legrand, and F. Mortessagne. Microwave experiments using open chaotic cavities in the realm of the effective Hamiltonian formalism. Fortschr. Phys., 61, 414–419 (2013)
  * [2] G. Gradoni, J.-H. Yeh, B. Xiao, T.-M. Antonsen , S.-M. Anlage, and E. Ott. Predicting the statistics of wave transport through chaotic cavities by the Random Coupling Model: a review and recent progress. Wave Motion, 51, 606–621 (2014)
  * [3] B. Dietz and A. Richter. Quantum and wave dynamical chaos in superconducting microwave billiards. Chaos, 25, 097601 (2015)
  * [4] H. Cao and J. Wiersig. Dielectric microcavities: Model systems for wave chaos and non-Hermitian physics.” Rev. Mod. Phys., 87, 61–111 (2015)
  * [5] O. Bohigas, M.-J. Giannoni and C. Schmit. Characterization of chaotic quantum spectra and universality of level fluctuation laws. Phys. Rev. Lett., 52, Issue 1, 1–5 (1984)
  * [6] J. J. M. Verbaarschot, H.A. Weidenmüller and M.R. Zirnbauer. Grassmann integration in stochastic quantum physics: the case of compound-nucleus scattering. Phys. Rep., 129, 367–438 (1985)
  * [7] Y. V. Fyodorov and H.-J. Sommers. Statistics of resonance poles, phase shifts and time delays in quantum chaotic scattering: Random matrix approach for systems with broken time-reversal invariance. J. Math. Phys. 38, Issue 4, 1918–1981 (1997)
  * [8] S. Hemmady, X. Zheng, E. Ott, T.M. Antonsen, S.M. Anlage. Universal impedance fluctuations in wave chaotic systems. Phys. Rev. Lett. 94, 014102 (2005)
  * [9] S. Hemmady, X. Zheng, J. Hart, T. M. Antonsen, E. Ott, and S. M. Anlage, Universal properties of two-port scattering, impedance, and admittance matrices of wave-chaotic systems. Phys. Rev. E 74, 036213 (2006)
  * [10] X. Zheng, S. Hemmady, T. M. Antonsen, Jr., S. M. Anlage, and E. Ott. Characterization of fluctuations of impedance and scattering matrices in wave chaotic scattering. Phys. Rev. E 73, 046208 (2006).
  * [11] Y. V. Fyodorov, D. V. Savin and H.-J. Sommers. Scattering, reflection and impedance of waves in chaotic and disordered systems with absorption. J. Phys. A: Math. Gen. 38, Issue 49, 10731–10760 (2005)
  * [12] Y.V. Fyodorov. Induced vs. spontaneous breakdown of S-matrix unitarity: probability of no return in quantum chaotic and disordered systems. JETP Letters 78, 250–254 (2003).
  * [13] Y.V. Fyodorov, D.V. Savin. Statistics of impedance, local density of states, and refection in quantum chaotic systems with absorption. JETP Letters 80, 725 –729 (2004).
  * [14] D.V. Savin, H.-J. Sommers, and Y.V. Fyodorov. Universal statistics of the local Green’s function in wave chaotic systems with absorption. JETP Letters 82, 544 – 548 (2005).
  * [15] D.V. Savin. Y. V. Fyodorov, and H.-J. Sommers. Correlation Functions of Impedance and Scattering Matrix Elements in Chaotic Absorbing Cavities. Acata Physica Polonica A 109, 53–64 (2006)
  * [16] S.B. Fedeli, Y.V. Fyodorov. Statistics of off-diagonal entries of Wigner K-matrix for chaotic wave systems with absorption. J Physics A: Math. Theor. 53 (16), 165701 (2020)
  * [17] R.A. Mendez-Sanchez, U. Kuhl, M. Barth, C.V. Lewenkopf, and H.-J. Stöckmann. Distribution of reflection coefficients in absorbing chaotic microwave cavities. Phys. Rev. Lett. 91, 174102 (2003).
  * [18] O. Hul, S. Bauch, P. Pakonski, N. Savytskyy, K. Zyczkowski, and L. Sirko. Experimental simulation of quantum graphs by microwave networks. Phys. Rev. E 69, 056205 (2004).
  * [19] O. Hul, O. Tymoshchuk, S. Bauch, P.M. Koch, and L. Sirko. Experimental investigation of Wigner’s reaction matrix for irregular graphs with absorption. J. Phys. A 38, 10489 (2005).
  * [20] M. Lawniczak, O. Hul, S. Bauch, P. Seba, and L. Sirko. Experimental and numerical investigation of the reflection coefficient and the distributions of Wigner’s reaction matrix for irregular graphs with absorption. Phys. Rev. E 77, 056210 (2008).
  * [21] M. Lawniczak, S. Bauch, O. Hul, and L. Sirko. Experimental investigation of the enhancement factor for microwave irregular networks with preserved and broken time reversal symmetry in the presence of absorption. Phys. Rev. E 81, 046204 (2010)
  * [22] M. Lawniczak and L. Sirko. Investigation of the diagonal elements of the Wigner’s reaction matrix for networks with violated time reversal invariance. Scientific Reports 9, Article number: 5630 (2019)
  * [23] M. Lawniczak, B. van Tiggelen, and L. Sirko. Experimental investigation of distributions of the off-diagonal elements of the scattering matrix and Wigner’s K-matrix for networks with broken time reversal invariance. Phys. Rev. E 102, 052214 (2020)
  * [24] M. Lawniczak, A. Akhshani, O. Farooq, M. Bialous, S. Bauch, B. Dietz, and L. Sirko. Distributions of the Wigner reaction matrix for microwave networks with symplectic symmetry in the presence of absorption. Phys. Rev. E 107, 024203 (2023)
  * [25] Y.V. Fyodorov, M. Skovortsov, and K.Tikhonov. Resonances in a single-lead reflection from a disordered medium: $\sigma$-model approach. arXiv:2211.03376
  * [26] U. Kuhl, R. Höhmann, J Main, and H.-J. Stöckmann. Resonance Widths in Open Microwave Cavities Studied by Harmonic Inversion. Phys. Rev. Lett. 100, 254101 (2008)
  * [27] L. Chen, S.M. Anlage and Y.V. Fyodorov. Statistics of Complex Wigner Time Delays as a counter of S-matrix poles: Theory and Experiment. Phys. Rev. Lett. 127 (20), 204101 (2021)
  * [28] V.V. Sokolov and V.G. Zelevinsky. Dynamics and statistics of unstable quantum states. Nucl. Phys. A 504, Issue 3, 562–588(1989)
  * [29] F. Haake, F. Izrailev, N. Lehmann, D Saher, and H.-J. Sommers. Statistics of complex levels of random matrices for decaying systems. Z. Phys. B: Cond. Matt. 88, 359-370 (1992)
  * [30] Y. V. Fyodorov and H.-J. Sommers. Statistics of S-matrix poles in few-channel chaotic scattering: Crossover from isolated to overlapping resonances. JETP Lett., 63 1026–1030 (1996).
  * [31] Y. V. Fyodorov, B. A. Khoruzhenko. Systematic Analytical Approach to Correlation Functions of Resonances in Quantum Chaotic Scattering. Phys. Rev. Lett., 83, 65 –68 (1999)
  * [32] H.-J. Sommers, Y. V. Fyodorov, and M. Titov. S-matrix poles for chaotic quantum systems as eigenvalues of complex symmetric random matrices: from isolated to overlapping resonances. J. Phys. A: Math. Gen., 32, L77 (1999)
  * [33] H. Schomerus, K.M. Frahm, M. Patra, C.W.J. Beenakker. Quantum limit of the laser line width in chaotic cavities and statistics of residues of scattering matrix poles. Physica A, 278, 469–496 (2000)
  * [34] Y.V. Fyodorov and B. Mehlig. Statistics of resonances and nonorthogonal eigenfunctions in a model for single-channel chaotic scattering. Phys Rev E 66 (4), 045202(R) (2002)
  * [35] K.B. Efetov. Supersymmetry in disorder and chaos. (Cambridge University Press, Cambridge, 1997)
  * [36] A.D. Mirlin, Y.V. Fyodorov. Statistical properties of one-point Green functions in disordered systems and critical behavior near the Anderson transition. J. de Phys. I 4 (5), 655–673 (1994)
  * [37] F. Wegner. The mobility edge problem: continuous symmetry and a conjecture. Z. Phys. B 35, 207–- 210 (1979)
  * [38] M.R. Zimbauer. Anderson localization and non-linear sigma model with graded symmetry. Nucl. Phys. B 265, Issue 2, 375–408 (1986)
  * [39] Y.V. Fyodorov, I. Rozhkov and J.R. Weaver. Statistics of transmitted power in multichannel dissipative ergodic structures. Phys. Rev. E 68, 016204 (2003).
  * [40] Y.V. Fyodorov and A.D. Mirlin. Scaling properties of localization in random band matrices: a $\sigma$-model approach. Phys. Rev. Lett. 67 (1991), no. 18, 2405–-2409 (1991)
  * [41] Y.V. Fyodorov and A.D. Mirlin. Statistical properties of eigenfunctions of random quasi 1D one-particle Hamiltonians. Int. J. Mod. Phys. B 8 (27), 3795–3842 (1994)
  * [42] M. A. Skvortsov and P. M. Ostrovsky. Local correlations of different eigenfunctions in a disordered wire. Pis’ma Zh. Eksp. Teor. Fiz. 85, 79 (2007) [JETP Lett. 85, 72 (2007)].
  * [43] A.D. Mirlin and Y.V. Fyodorov. Distribution of the local densities of states, order parameter function, and critical behavior near the Anderson transition. Phys. Rev. Lett. 72, 526, (1994)
  * [44] G Schubert, Jens Schleede, K. Byczuk, H. Fehske and D. Vollhardt. Distribution of the local density of states as a criterion for Anderson localization: numerically exact results for various lattices in two and three dimensions. Phys. Rev. B 81, 155106 (2010)
  * [45] L.S. Gradshteyn, I.M. Ryzhik. Table of Integrals, Series, and Products. Ed. A. Jeffrey $\&$ D. Zwillinger. 6th Ed. (Academic Press, 2000).


# Long-range magnetic order induced surface state in GdBi and DyBi

Yevhen Kushnirenko<EMAIL_ADDRESS>Division of Materials Science and
Engineering, Ames Laboratory, Ames, Iowa 50011, USA  Department of Physics and
Astronomy, Iowa State University, Ames, Iowa 50011, USA  Lin-Lin Wang
Division of Materials Science and Engineering, Ames Laboratory, Ames, Iowa
50011, USA  Zhuoqi Li  Division of Materials Science and Engineering, Ames
Laboratory, Ames, Iowa 50011, USA  Department of Physics and Astronomy, Iowa
State University, Ames, Iowa 50011, USA  Brinda Kuthanazhi  Division of
Materials Science and Engineering, Ames Laboratory, Ames, Iowa 50011, USA
Department of Physics and Astronomy, Iowa State University, Ames, Iowa 50011,
USA  Benjamin Schrunk  Division of Materials Science and Engineering, Ames
Laboratory, Ames, Iowa 50011, USA  Evan O’Leary  Division of Materials Science
and Engineering, Ames Laboratory, Ames, Iowa 50011, USA  Department of Physics
and Astronomy, Iowa State University, Ames, Iowa 50011, USA  Andrew Eaton
Division of Materials Science and Engineering, Ames Laboratory, Ames, Iowa
50011, USA  Department of Physics and Astronomy, Iowa State University, Ames,
Iowa 50011, USA  P. C. Canfield  Division of Materials Science and
Engineering, Ames Laboratory, Ames, Iowa 50011, USA  Department of Physics and
Astronomy, Iowa State University, Ames, Iowa 50011, USA  Adam Kaminski
<EMAIL_ADDRESS>Division of Materials Science and Engineering, Ames
Laboratory, Ames, Iowa 50011, USA  Department of Physics and Astronomy, Iowa
State University, Ames, Iowa 50011, USA

###### Abstract

The recent discovery of unconventional surface-state pairs, which give rise to
Fermi arcs and spin textures, in antiferromagnetically ordered rare-earth
monopnictides attracted the interest in these materials. We use angle-resolved
photoemission spectroscopy (ARPES) measurements in conjunction with density
functional theory (DFT) calculations to investigate the evolution of the
electronic structure of GdBi and DyBi. We find that new surface states,
including a Dirac cone, emerge in the AFM state. However, they are located
along a direction in momentum space that is different than what was found in
NdSb, NdBi, and CeBi [1, 2, 3, 4, 5]. The observed changes in the electronic
structure are consistent with the presence of AFM-II-A type order.

††preprint: APS/123-QED

##  I Introduction

The rare-earth monopnictides family of materials crystallizes in a simple NaCl
structure. Many members of this family order antiferromagnetically (AFM) [6,
7, 8, 9], and some of them have complex phase diagrams with multiple AFM
phases [10, 11, 12, 13, 14]. Recent predictions [15, 16, 17, 18, 19, 20] of
the existence of topological states that result from the inversion between 6p
band of pnictogen element and 5d band of rare-earth element in some of these
materials inspired more experimental studies, including angle-resolved
photoemission spectroscopy (ARPES). These studies demonstrated not only
topological Dirac states [21, 22, 23, 24, 25, 26, 27, 28] but also the
development of interesting surface states (SS) in NdBi [1, 2, 3], NdSb [1, 4,
5], and CeBi [1, 14] upon AFM transition. These states undergo Kaminski-
Canfield (KC) band splitting that leads to the formation of spin textured
Fermi arcs.

GdBi and DiBi were also shown to have an AFM phase. The Neel temperature is
$T_{N}$ = 27.5 K for GdBi [7, 29, 30] and $T_{N}$ = 11 K for DyBi [6, 31, 32].
The AFM structure of these materials is different from that of NdBi, NdSb, and
CeBi. Neutron diffraction studies [7, 33] have shown, that in both GdBi and
DyBi the magnetic moments are aligned in type-II AFM structure. However,
neutron diffraction measurements cannot distinguish AFM-II-A (Fig. 2(b)) from
AFM-II-B (Fig. 2(c)) order [34]. Similarly to several other rare-earth
monopnictides, the DFT calculations for GdBi in the paramagnetic (PM) phase
[17, 25, 30] predict a band inversion between Bi 6p and Gd 5d bands along
$\Gamma-X$ direction. This band inversion should lead to the formation of one
SS Dirac cone at the $\overline{\Gamma}$-point and two cones at
$\overline{M}$-point of the 2D Brillouin zone (BZ) [24, 3]. An earlier
synchrotron-based ARPES measurements [25] demonstrated the presence of these
Dirac states.

In this work, we investigate the evolution of GdBi and DyBi electronic
structures upon cooling through their AFM transitions using high-resolution
laser-based ARPES. We demonstrate the development of Dirac states in the
magnetically ordered phase near the $\overline{X}$-point.

##  II Experimental details

Single crystals of GdBi and DyBi were grown out of In flux [35]. The elements
with an initial stoichiometry of R4Bi4In96 (R = Gd, Dy) were put into a
fritted alumina crucible [36] (sold as Canfield crucible set [37]) and sealed
in a fused silica tube under partial pressure of argon. The prepared ampoules
were heated to 1150∘ C over 5 hours and held there for 5 hours. This was
followed by a slow cooling to 700∘ C (for GdBi) and 400∘ C (for DyBi) over 100
hours and decanting of the excess flux through the frit using a
centrifuge.[38] The cubic crystals obtained were stored and handled in a
glovebox under nitrogen atmosphere.

Most of the ARPES data were collected using vacuum ultraviolet (VUV) laser
ARPES spectrometer that consists of a Scienta DA30 electron analyzer, tunable
picosecond Ti:Sapphire oscillator paired with fourth-harmonic generator [39]
and 3.5 eV continuous laser paired with second-harmonic generator. Data were
collected with 6.2-7 eV photon energy. Angular resolution was set at $\sim$
0.1∘ and 1∘, along and perpendicular to the direction of the analyzer slit,
respectively, and the energy resolution was set at 2 meV. The VUV laser beam
was set to vertical polarization. The diameter of the photon beam on the
sample was $\sim 15\,\mu$m. The measurements at photon energy of 21.2 eV were
carried out using R8000 analyzer and GammaData helium discharge lamp with
custom focusing optics. The diameter of the photon beam on the sample was
$\sim 1$mm. Samples were cleaved in-situ along (001) plane, usually producing
very flat, mirror-like surfaces. The measurements were performed at a base
pressure lower than 2$\times$10-11 Torr.

Figure 1: (a) Results of magnetization measurements ($d(\chi T)/dT$) of GdBi
and DyBi. (b) Fermi surface maps of NdBi measured in the AFM state. (c) and
(d) Fermi surface maps of GdBi measured in the PM and AFM phases,
respectively. (e) and (f) Fermi surface maps of DyBi measured in the PM and
AFM phases, respectively. (g) and (h) Band dispersion measured in GdBi in the
PM and AFM phases, respectively, along $\overline{\Gamma}-\overline{X}$
direction (marked as #1 in (c)). (i) and (j) Band dispersion measured in GdBi
in the PM and AFM phases, respectively, along $\overline{\Gamma}-\overline{M}$
direction (marked as #2 in (c)). The cuts along #3 and #4 are shown in Fig. 3.
The green arrows in (b, d, f) mark SS states that appear in the AFM phase. The
red arrows in (d) mark the star-shaped bulk pocket, which appears in the AFM
phase.

##  III Results and Discussion

Temperature-dependent magnetic susceptibility measurements made in H = 1 kOe
were made on GdBi and NdBi. The $d(\chi T)/dT$ data (fig. 1a) show that there
is clear AFM ordering below 27.3 K and 11 K for GdBi and DyBi, respectively.
These results are in agreement with earlier reported $T_{N}$ values. [7, 29,
30, 6, 31, 32]

Figure 2: Plot of calculated FS using DFT. (a) 2D projections of bulk Fermi
surface for PM phases. (b) and (c) schematic representation of AFM-II-A and
AFM-II-B order, respectively. (d) and (e) 2D projections of bulk+surface Fermi
surface for AFM-II-A and AFM-II-B order, respectively. In type AFM-II-A, the
four face-centered cubic sublattices of antiparallel spins are correlated such
that on each (111) plane, the spins are parallel, and the spins on two
neighboring (111) planes are opposite. In type AFM-II-B, the four face-
centered cubic sublattices divide the full face-centered cubic lattice in such
a way that no two nearest neighboring ions are in the same sublattice. [34]

In Fig. 1(c), we show a Fermi surface map of GdBi measured in the paramagnetic
state. This Fermi surface is similar to the ones of other rare-earth
monopnictides. It consists of several pockets around the center of the BZ,
formed mostly by bulk states. This result agrees with DFT calculations [Fig.
2(a)]. Besides these pockets, the calculations predict the existence of two
pockets formed by bulk states around the corner of the BZ. One of these
pockets is partially seen in our experimental data. In the AFM state, we
observe the development of additional states along
$\overline{\Gamma}-\overline{X}$ direction. These states form a small sharp
pocket (marked with green arrow) and a bigger broad star-shaped pocket (marked
with red arrows) on the FS (Fig. 1(d)). Both of these pockets are centered
around $\overline{X}$-point. The size and shape of the bigger pocket resemble
the outer bulk pocket located around the $\overline{\Gamma}$-point. The
comparison of the cuts taken along $\overline{\Gamma}-\overline{X}$ from PM
and AFM data sets (Fig. 1(j, h)) show that the additional small pocket is
formed by linear dispersion. The bands located near $\overline{\Gamma}$-point
do not show any changes upon the magnetic transition. The cuts taken along
$\overline{\Gamma}-\overline{M}$ direction (Fig. 1(g, j)) also show no
qualitative change upon the magnetic transition. In Fig. 1(e, f), we show
Fermi surface maps measured from DyBi in the PM and AFM states, respectively.
These maps also demonstrate the development of additional states near the
$\overline{X}$-point in the AFM state.

Figure 3: Band dispersion in AFM state of GdBi. (a)-(c) Band dispersion
measured along $\overline{\Gamma}-\overline{Y}$ direction (marked as #3 in
Fig. 1(d)) using 6.7, 7, and 6.2 eV photons (left) and corresponding second
derivative plots (right). These three spectra were measured at T = 6.5, 4.5,
5.5 K, respectively. (d) MDCs taken from (a) at 382 and 408 meV binding energy
(marked with arrows of corresponding colors). (e) and (f) Band dispersion
measured along $\overline{\Gamma}-\overline{Y}$ direction in the first and the
second BZ using 21.2 eV light from He-lamp at T = 6.5 K. (g) Band dispersions
extracted from plots (a),(d)-(f) by fitting MDCs. Points at zero-momentum were
extracted by fitting EDCs. (h) Band dispersion measured at T = 6.5 K along the
$\overline{X}-\overline{M}$ direction (marked as #4 in Fig 1.(d)).

These results are different from other rare-earth monopnictides [1, 2, 3, 4,
5, 14]. For comparison, we demonstrate a dataset measured from NdBi in AFM
state in Fig. 1(b). These data do not indicate the presence of any states
around the $\overline{X}$-point but instead show the development of Fermi arcs
and elliptical pockets along $\overline{\Gamma}-\overline{M}$ direction.

The spectrum in Fig. 3(a) is measured in the AFM phase along
$\overline{\Gamma}-\overline{Y}$ direction (cut #3 marked in Fig. 1(d)). It
demonstrates bands that, at first glance, resembles a Dirac cone. However, a
closer inspection reveals that its upper part, instead of merging at one point
with the lower part, crosses it, and the top of the lower band is rounded. The
crossing can be even better seen in the second derivative plot and MDCs in
Fig. 3(d). The MDC extracted through the top of the hole-like band(E = -382
meV) has two small peaks on the background of a large, broad peak. The
additional peaks correspond to the upper band. The MDC extracted at E = -408
meV (below the top of the holelike band) has three peaks present. The outer
two peaks correspond to the hole-like dispersion, and the central peak
corresponds to the bottom of the upper band. The upper band is linear and less
intense than hole-like dispersion. Thus, we can assume that they are parts of
the Dirac cone, the lower part of which cannot be properly distinguished
because of the strong background of the hole-like pocket. Only some hints of
the lower part of the Dirac cone can be seen on the second derivative plot. To
demonstrate the SS nature of the linear dispersions, we measured spectra
similar to the one in Fig. 3(a) using light with different photon energy (see
Fig. 3(b,c)). The linear dispersion shows no changes with photon energy. We
fitted the shape of linear dispersions in the 6.7 eV spectrum and plotted the
results over the second derivative plots of 6.2 and 7.0 eV spectra as dashed
lines. The curves perfectly match the shape of the linear dispersions in these
spectra. At the same time, the hole-like dispersion shows a substantial
$k_{z}$-dependence: its top shifts by 65 meV between 6.2 and 7.0 eV spectra,
and thus, it is definitely a bulk band, while the linear dispersion has SS
origin.

For a better understanding of the electronic structure and, in particular, to
prove the existence of the lower part of the Dirac cone, we performed
measurements using a He-lamp ($h\nu$ = 21.2 eV). Fig. 3(e, f) demonstrate the
spectra measured in $\overline{\Gamma}-\overline{Y}$ direction in the first
and the second BZ, respectively. The lower part of the Dirac cone can be seen
in both spectra. It is better seen in Fig. 3(e): here, we see two linear
dispersions that cross at E = -406 meV and a blob of intensity above them (E =
-356 meV), which corresponds to the top of the holelike bulk band. For further
analysis, we extracted band dispersion from spectra in Fig. 3(a, e, f) by
fitting the experimental data (see Fig. 6 in the appendix for details). The
results of these fits are shown in Fig. 3(g). The linear dispersions extracted
from 6.7 and 21.2 eV spectra match each other very well, which again proves
that they are Dirac surface states.

Figure 4: Temperature dependence of the electronic structure of GdBi. (a)
Spectra measured along the $\overline{X}-\overline{M}$ direction at several
temperatures. (b) EDCs taken from data in panel (a) at zero-momentum (marked
with dashed line). (c) Temperature dependence of the area under the main peak
from the EDCs in (b) after subtracting the area at 30 K.

In Fig. 3(h), we show a spectrum measured through $\overline{X}$-point (cut #4
in Fig 1.(d)) in the AFM state. This spectrum has a sharp Dirac cone on the
background of broad features. Strictly speaking, this spectrum is measured not
along a straight cut for a constant value of $k_{x}$ component of momentum
because of the geometry of the ARPES experiment. Because of this, in order to
extract the real shape of the bands near the $\overline{X}$-point, several
spectra should be analyzed (see SM for details). The shape of the upper part
of the cone (cut through its center) extracted this way is shown in Fig. 3(g)
with a black line. This shape matches the shape of the Dirac cone at the
$\overline{\Gamma}$-point very well. This indicates that the bands that appear
near the $\overline{X}$-point in the AFM state are the result of the folding
of the bands from the center of the BZ. Another consequence of measuring at
high angles of emission is that the photoemission experiment is expected to
become less sensitive to the bulk states. Thus, the SS should be more
pronounced. And indeed, the Diraclike feature folded from the
$\overline{\Gamma}$-point becomes dominant near the $\overline{X}$-point. This
is additional evidence that this Dirac-like feature at the
$\overline{X}$-point and, as a consequence, the Dirac-like feature at the
$\overline{\Gamma}$-point are originating from SS. Despite low bulk
sensitivity, some bulk features can still be identified in Fig. 3(h). We have
extracted the shape of the band that crosses the Fermi level at $\sim$0.13
$\AA^{-1}$ and compared it with the shape of the outermost band in Fig. 3(a).
The results are shown in Fig. 3(g). These two band dispersions match each
other reasonably well, which again demonstrates that the bands the
$\overline{X}$-point are the result of band folding. The small difference
between bulk bands in Fig. 3(a, h) can be explained by the $k_{z}$ dispersion
of this bulk band ($k_{z}$ component also depends on the angle of emission).

Figure 5: Band dispersion of DyBi. (a) and (b) Spectrum measured along
$\overline{\Gamma}-\overline{Y}$ direction at T = 16 K and the corresponding
second derivative plot. (c) and (d) Band dispersion measured T = 6.6 K through
the first and second $\overline{\Gamma}$-point using 21.2 eV light from He-
lamp. (e) and (f) Spectrum measured near the $\overline{X}$-point in AFM and
PM state, respectively.

To analyze the evolution of the electronic structure near the
$\overline{X}$-point with temperature, we measured a series of spectra through
the $\overline{X}$-point at several temperatures below and above $T_{N}$. Five
of these spectra are shown in Fig. 4(a). All spectra measured in the AFM state
have the same set of bands. These bands do not change their shape or location
with temperature. This can be better seen in Fig. 4(b), which demonstrates
EDCs taken through the band center. The peak which corresponds to the Dirac
point does not sift with temperature. However, the intensity of all bands at
these spectra tends to decrease with increasing temperature, and all of them
disappear above $T_{N}$ in the PM state. In Fig. 4(c), we plot a temperature
dependence of the area under the EDC peak after subtracting the area above
$T_{N}$. This plot demonstrates that the intensity of the bands folded to the
area near the $\overline{X}$-point follows the AFM order parameter.

The folding of the bands from the center of the BZ to the region near the
$\overline{X}$-point is predicted by the DFT calculations for the case of AFM-
II-A order [Fig. 2(d)] agrees with our experimental results. Based on our
data, we can exclude the scenario of AFM-II-B ordering. In this case, the
folding of the bands from the $\overline{M}$-point to the center of the BZ is
expected (Fig. 2(e)). This contradicts our experimental results, which show no
changes at the $\overline{\Gamma}$-point upon AFM transition (see Fig. 1(g-j)
and Fig. 7 in the appendix).

As was already mentioned, the FS of DyBi in the AFM phase (Fig. 1(f)) is
similar to the one of GdBi (Fig. 1(d)). For further analysis, we show the
Spectrum measured along $\overline{\Gamma}-\overline{Y}$ direction (Fig. 5(a))
and the corresponding second derivative plot (Fig. 5(b)). This spectrum has
the same set of bands as Fig. 4(a-c), including a cone-like feature that
crosses a more intense hole-like band. The spectra measured using 21.2 eV
photons along $\overline{\Gamma}-\overline{Y}$ direction in the first and the
second BZ are shown in Fig. 5(c, d), respectively. These spectra are similar
to the spectra measured from GdBi. The separation between the lower part of
the cone and the blob of intensity, which we associate with the top of the
bulk hole-like dispersion, is more clearly seen in DyBi than in GdBi.

Fig. 5(e, f) demonstrate spectra measured near the $\overline{X}$-point in the
AFM and PM state, respectively. In this experiment, we could not reach the
$\overline{X}$-point where the upper and the lower parts of the cone are
expected to touch each other. Nevertheless, in Fig. 5(e), one can see the
upper parts of the cone and broad bulk features that are absent in the PM
phase. This indicates that, in DyBi the bands from the center of the BZ are
folded to the X-point in the AFM state in the same way as it happens in GdBi.

##  IV CONCLUSIONS

Our ARPES measurements demonstrate the emergence of new surface states near
the $\overline{X}$-point in GdBi and DyBi upon a transition to AFM phase.
These states are different from those observed in rare-earth monopnictides
with AFM-I order (NdSb, NdBi, and CeBi). The observed additional states are
likely the result of the band folding from the center of the BZ. These states
include a Dirac cone, which was observed in both $\overline{\Gamma}$\- and
$\overline{X}$-points of the BZ. The electronic structure near the center of
the BZ does not change upon AFM transition. The observed results are
consistent with AFM-II-A ordering in these materials.

Figure 6: Supplementary figure. (a-c) tree projections of a set of points
extracted by fitting 13 spectra measured near the $\overline{X}$-point. The
red curve shows the share of the band dispersion along the
$\overline{X}-\overline{M}$ cut extracted by fitting this set of points. (d)
correspondence between the experimental data and the results of the fit. (e)
set of points extracted by EDC and MDC fitting of the spectrum in Fig. 3(a).
(f) set of points extracted by EDC and MDC fitting of the spectra in Fig.
3(e,f).

##  Appendix A Band shape fitting ($\overline{X}$-point)

The measurement of the band structure near the $\overline{X}$-point with the
6.7 eV laser has to be performed at large take-off angles of the
photoelectron: $\Theta\sim$ 71 degrees to reach the Dirac point. At such
angles, the fact that the $k_{y}$ is not equal in all parts of the spectrum
cannot be neglected. Thus, in order to extract the shape of the bands along a
straight cut at a constant value of $k_{y}$, we need to compile data from
several closely located cuts.

The shape of the Dirac cone near the $\overline{X}$-point was extracted by
fitting 13 spectra measured at angles from 51 to 75 degrees. The set of points
expected by fitting MDC peaks is shown in Fig. 6(a-c). In the next step, we
fit this set of points with the following function:
$E=E_{0}+A_{1}R+A_{2}R+A_{3}R^{3}$, where
$R=\sqrt{(k_{x}-k_{x0})^{2}+(k_{y}-k_{y0})^{2}}(1+Bcos(4\Theta))$, and
$\Theta=arctan(k_{x}-k_{x0})/(k_{y}-k_{y0}))$. Here $Bcos(4\Theta)$ component
represents a deviation of the pocket shape from a perfect circle. The fit
gives the following values of the coefficients: $E_{0}=$ 245$\pm$3; $k_{x0}=$
0.6809$\pm$0.0001; $k_{y0}=$ 0.008$\pm$0.0001; $A_{1}=$ 3.87$\pm$0.18;
$A_{2}=$ 5.4$\pm$3.4; $A_{3}=$ 19$\pm$19; $B=$ -0.046$\pm$0.011.

Fig. 6(d), which shows a correspondence between the experimental data and the
results of the fit, and the colors in Fig. 6(c), which show a difference
between the experimental data and the results of the fit for different parts
of the band, indicate good quality of the fit. If we extend the function,
which describes the band shape, to the lower cone by changing the sign of
effective radius, the result would also match the experimental data in that
energy range (see red markers in Fig. 6(d)). The fit coefficients for non-
linear components of the fit ($A_{2}R^{2}$ and $A_{3}R^{3}$) are comparable to
the fit error. Also, these components are insignificant, which can be visually
seen from the shape of the red curve in Fig. 6(b), which represents the
results of the fit at $k_{x}$ = $k_{x0}$. All this indicates that the observed
band is indeed a Dirac cone.

##  Appendix B Band shape fitting ($\overline{\Gamma}$-point)

In Fig. 6(e, f), we demonstrate the intermediate steps of extracting the shape
of the bands in Fig. 2(a,d, e). The markers at $k_{x}\neq 0$ represent the
results of fitting MDC peaks, and the markers $k_{x}=0$ represent the results
of fitting EDC peaks. The lines represent the results of fitting these points.
The resulting shape of the upper band from the spectrum measured with the
laser was fitted with a line, and the shape of the lower band was fitted with
a polynomial of degree four. The shape of the band from the spectra measured
with the He-lamp was fitted with a line.

Figure 7: GdBi temperature dependence: additional data. (a) Spectra measured
along $\overline{\Gamma}-\overline{Y}$ direction using 21.2 eV light (top) and
corresponding second derivative plots (bottom). (b) Spectra measured along
$\overline{\Gamma}-\overline{Y}$ direction using 7.0 eV light.

## References

  * Kushnirenko  _et al._ [2022] Y. Kushnirenko, B. Schrunk, B. Kuthanazhi, L.-L. Wang, J. Ahn, E. O’Leary, A. Eaton, S. L. Bud’ko, R.-J. Slager, P. C. Canfield, and A. Kaminski, Phys. Rev. B 106, 115112 (2022).
  * Schrunk  _et al._ [2022] B. Schrunk, Y. Kushnirenko, B. Kuthanazhi, J. Ahn, L.-L. Wang, E. O’Leary, K. Lee, A. Eaton, A. Fedorov, R. Lou, V. Voroshnin, O. J. Clark, J. Sánchez-Barriga, S. L. Bud’ko, R.-J. Slager, P. C. Canfield, and A. Kaminski, Nature 603, 610–615 (2022).
  * Honma  _et al._ [2023a] A. Honma, D. Takane, S. Souma, K. Yamauchi, Y. Wang, K. Nakayama, K. Sugawara, M. Kitamura, K. Horiba, H. Kumigashira, _et al._ , Nature Communications 14, 7396 (2023a).
  * Kushnirenko  _et al._ [2023] Y. Kushnirenko, B. Kuthanazhi, L.-L. Wang, B. Schrunk, E. O’Leary, A. Eaton, P. C. Canfield, and A. Kaminski, Physical Review B 108, 115102 (2023).
  * Honma  _et al._ [2023b] A. Honma, D. Takane, S. Souma, Y. Wang, K. Nakayama, M. Kitamura, K. Horiba, H. Kumigashira, T. Takahashi, Y. Ando, _et al._ , Physical Review B 108, 115118 (2023b).
  * Tsuchida and Wallace [1965] T. Tsuchida and W. E. Wallace, The Journal of Chemical Physics 43, 2087 (1965).
  * Nereson and Arnold [1971] N. Nereson and G. Arnold, Journal of Applied Physics 42, 1625 (1971).
  * Bartholin  _et al._ [1979] H. Bartholin, P. Burlet, S. Quezel, J. Rossat-Mignod, and O. Vogt, Le Journal de Physique Colloques 40, C5 (1979).
  * Kumigashira  _et al._ [1996] H. Kumigashira, S.-H. Yang, T. Yokoya, A. Chainani, T. Takahashi, A. Uesawa, T. Suzuki, O. Sakai, and Y. Kaneta, Phys. Rev. B 54, 9341 (1996).
  * Rossat-Mignod  _et al._ [1983] J. Rossat-Mignod, P. Burlet, S. Quezel, J. Effantin, D. Delacôte, H. Bartholin, O. Vogt, and D. Ravot, Journal of Magnetism and Magnetic Materials 31, 398 (1983).
  * Wiener and Canfield [2000] T. Wiener and P. Canfield, Journal of Alloys and Compounds 303, 505 (2000).
  * Kuroda  _et al._ [2020] K. Kuroda, Y. Arai, N. Rezaei, S. Kunisada, S. Sakuragi, M. Alaei, Y. Kinoshita, C. Bareille, R. Noguchi, M. Nakayama, S. Akebi, M. Sakano, K. Kawaguchi, M. Arita, S. Ideta, K. Tanaka, H. Kitazawa, K. Okazaki, M. Tokunaga, Y. Haga, S. Shin, H. S. Suzuki, R. Arita, and T. Kondo, Nature Communications 11, 1 (2020).
  * Kuthanazhi  _et al._ [2022] B. Kuthanazhi, N. H. Jo, L. Xiang, S. L. Bud’ko, and P. C. Canfield, Philosophical Magazine 102, 542 (2022).
  * Kushnirenko [2024] Y. Kushnirenko, Unpublished (2024).
  * Guo  _et al._ [2017] C. Guo, C. Cao, M. Smidman, F. Wu, Y. Zhang, F. Steglich, F.-C. Zhang, and H. Yuan, Npj Quantum Materials 2, (2017).
  * Li  _et al._ [2017] Z. Li, D.-D. Xu, S.-Y. Ning, H. Su, T. Iitaka, T. Tohyama, and J.-X. Zhang, International Journal of Modern Physics B 31, 1750217 (2017).
  * Duan, Xu  _et al._ [2018] Duan, Xu, Wu, Fan, Chen, Jia, Zhang, Peiran, Liu, Yang, Yuan, Huiqiu, and Cao, Chao, Communications Physics 1, 71 (2018).
  * Huang  _et al._ [2020] Z. Huang, C. Lane, C. Cao, G.-X. Zhi, Y. Liu, C. E. Matt, B. Kuthanazhi, P. C. Canfield, D. Yarotski, A. J. Taylor, and J.-X. Zhu, Phys. Rev. B 102, 235167 (2020).
  * Zeng  _et al._ [2015] M. Zeng, C. Fang, G. Chang, Y.-A. Chen, T. Hsieh, A. Bansil, H. Lin, and L. Fu, arXiv preprint arXiv:1504.03492 (2015).
  * Wang  _et al._ [2023] L.-L. Wang, J. Ahn, R.-J. Slager, Y. Kushnirenko, B. Ueland, A. Sapkota, B. Schrunk, B. Kuthanazhi, R. McQueeney, P. Canfield, and A. Kaminski, Communications Physics (2023).
  * Zeng  _et al._ [2016] L.-K. Zeng, R. Lou, D.-S. Wu, Q. Xu, P.-J. Guo, L.-Y. Kong, Y.-G. Zhong, J.-Z. Ma, B.-B. Fu, P. Richard, _et al._ , Physical review letters 117, 127204 (2016).
  * Niu  _et al._ [2016] X. Niu, D. Xu, Y. Bai, Q. Song, X. Shen, B. Xie, Z. Sun, Y. Huang, D. Peets, and D. Feng, Physical Review B 94, 165163 (2016).
  * Lou  _et al._ [2017] R. Lou, B.-B. Fu, Q. Xu, P.-J. Guo, L.-Y. Kong, L.-K. Zeng, J.-Z. Ma, P. Richard, C. Fang, Y.-B. Huang, _et al._ , Physical Review B 95, 115140 (2017).
  * Kuroda  _et al._ [2018] K. Kuroda, M. Ochi, H. S. Suzuki, M. Hirayama, M. Nakayama, R. Noguchi, C. Bareille, S. Akebi, S. Kunisada, T. Muro, M. D. Watson, H. Kitazawa, Y. Haga, T. K. Kim, M. Hoesch, S. Shin, R. Arita, and T. Kondo, Phys. Rev. Lett. 120, 086402 (2018).
  * Li  _et al._ [2018] P. Li, Z. Wu, F. Wu, C. Cao, C. Guo, Y. Wu, Y. Liu, Z. Sun, C.-M. Cheng, D.-S. Lin, _et al._ , Physical Review B 98, 085103 (2018).
  * Oinuma  _et al._ [2019] H. Oinuma, S. Souma, K. Nakayama, K. Horiba, H. Kumigashira, M. Yoshida, A. Ochiai, T. Takahashi, and T. Sato, Phys. Rev. B 100, 125122 (2019).
  * Sakhya  _et al._ [2022a] A. P. Sakhya, S. Kumar, A. Pramanik, R. P. Pandeya, R. Verma, B. Singh, S. Datta, S. Sasmal, R. Mondal, E. F. Schwier, _et al._ , Physical Review B 106, 085132 (2022a).
  * Sakhya  _et al._ [2022b] A. P. Sakhya, B. Wang, F. Kabir, C.-Y. Huang, M. M. Hosen, B. Singh, S. Regmi, G. Dhakal, K. Dimitri, M. Sprague, R. Smith, E. D. Bauer, F. Ronning, A. Bansil, and M. Neupane, Phys. Rev. B 106, 235119 (2022b).
  * Li  _et al._ [1996] D. Li, Y. Haga, H. Shida, T. Suzuki, and Y. Kwon, Physical Review B 54, 10483 (1996).
  * Dwari  _et al._ [2023] G. Dwari, S. Sasmal, S. Dan, B. Maity, V. Saini, R. Kulkarni, S. Banik, R. Verma, B. Singh, and A. Thamizhavel, Physical Review B 107, 235117 (2023).
  * Hulliger [1980] F. Hulliger, Journal of Magnetism and Magnetic Materials 15, 1243 (1980).
  * Zhao  _et al._ [2022] J. Zhao, A. Rahman, W. Liu, F. Meng, L. Ling, Y. Han, C. Xi, W. Tong, L. Xu, Z. Tian, _et al._ , Physical Review B 106, 224412 (2022).
  * Nereson and Struebing [1972] N. Nereson and V. Struebing, in  _AIP Conference Proceedings_ , Vol. 5 (American Institute of Physics, 1972) pp. 1385–1389.
  * Li [1955] Y.-Y. Li, Physical Review 100, 627 (1955).
  * Canfield and Fisk [1992] P. C. Canfield and Z. Fisk, Philosophical Magazine Part B 65, 1117 (1992).
  * Canfield  _et al._ [2016] P. C. Canfield, T. Kong, U. S. Kaluarachchi, and N. H. Jo, Philosophical Magazine 96, 84 (2016).
  * [37] URL:, https://lspceramics.com/canfield-crucible-sets-2/ .
  * Canfield [2019] P. C. Canfield, Reports on Progress in Physics 83, 016501 (2019).
  * Jiang  _et al._ [2014] R. Jiang, D. Mou, Y. Wu, L. Huang, C. D. McMillen, J. Kolis, H. G. Giesber, J. J. Egan, and A. Kaminski, Review of Scientific Instruments 85, 033902 (2014).