11institutetext: V.Yu. Irkhin and S.V. Streltsov22institutetext: M. N. Mikheev
Institute of Metal Physics, 620108 Ekaterinburg, Russia
22email<EMAIL_ADDRESS>
# S. Vonsovsky and the dawn of the theory of strongly correlated systems

V. Yu. Irkhin  S.V. Streltsov

(Received: / Accepted: date)

###### Abstract

A survey of scientific and social activity of S. V. Vonsovsky is presented.
His fundamental works on the theory of electronic properties and magnetism of
transition and rare-earth metals and their compounds are reviewed. The
historical development of the many-electron polar (preceding the famous
Hubbard model) and $s-d(f)$ exchange (“Vonsovsky-Zener”, “Kondo lattice”)
models is considered.

###### Keywords:

magnetism electron correlations polar model $s-d(f)$ exchange model Kondo
effect

###### pacs:

71.27.+a 75.10.Lp 71.30.+h

Sergey Vasiljevich Vonsovsky (1910-1998) was one of founders of the Institute
of Metal Physics in Sverdlovsk (now Yekaterinburg), as well as the entire Ural
scientific school on the quantum theory of solids and the physics of magnetic
phenomena. Most of the Ural theoretical physicists of the last century
considered him as their teacher.

Vonsovsky was the founder of the famous “Kourovka” - the winter school of
theoretical physics, which has been held in the Urals since 1960 and is known
for its warm and democratic atmosphere. This volume of the Journal of
Superconductivity and Novel Magnetism represents contributions of past and
present Kourovka attendees.

In the present paper, we discuss most important scientific works by Vonsovsky,
especially in the many-electron theory of solids (together with his teacher
Semen Petrovich Shubin). These contributions are not very well known to a
younger generation grown on articles by J. Hubbard and P. W. Anderson (older
colleagues are of course familiar with them), while ideas introduced by S.
Vonsovsky into the field of strongly correlated electron systems were
nonetheless seminal. This is one of the aims of the present paper - to shortly
remind about two of his models: the polar (Section 1), and $s-d(f)$ exchange
(Section 2) models. Some historical background is presented in more detail in
Section 3.

##  1 Polar model

In their works on the polar model 1 ; 2 , Shubin and Vonsovsky stated a goal
to simultaneously describe a wide range of phenomena in solids, including
magnetism and electrical conductivity. Already in 1934-1935 (i.e. nearly 30
years before Hubbard) they formulated a minimal model for interacting
electrons in a solid. This so-called polar model was proposed as a synthesis
of the Heisenberg homeopolar model describing the system of localized moments,
and the Slater method for describing the itinerant-electron system of a metal.
In a crystal with one electron per atom, it means taking into account polar
states –- doubles and holes, i.e. doubly occupied and empty sites, and leads
to a physically transparent picture. In the original formulation of the model,
hoppings of electrons from site to site were taken into account, as well as
all types of interelectron interaction. These papers, which provided the
foundations of many-electron solid-state physics, were carried out by means of
difficult and cumbersome calculations (there was no convenient method of
second quantization at that time).

The polar model was further generalized by Bogolyubov bog , who derived its
Hamiltonian via the successive expansion in terms of the hopping parameters
(loosely speaking the overlap between atomic wave functions) in the second-
quantization representation. In the simplest case of a nondegenerate band,
this Hamiltonian can be written in the form

$\displaystyle{\mathcal{H}}$ $\displaystyle=$
$\displaystyle\sum_{\nu_{1}\neq\nu_{2},\sigma}t_{\nu_{1}\nu_{2}}c_{\nu_{1}\sigma}^{\dagger}c_{\nu_{2}\sigma}$
(1) $\displaystyle+$
$\displaystyle\frac{1}{2}\sum_{\nu_{i}\sigma_{1}\sigma_{2}}I_{\nu_{1}\nu_{2}\nu_{3}\nu_{4}}c_{\nu_{1}\sigma_{1}}^{\dagger}c_{\nu_{2}\sigma_{2}}^{\dagger}c_{\nu_{4}\sigma_{2}}c_{\nu_{3}\sigma_{1}}.$

Here, $c_{\nu\sigma}^{\dagger}$ are the creation operators for an electron at
the site $\nu$ with the spin projection $\sigma$, $t_{\nu_{1}\nu_{2}}$ are the
hopping integrals (matrix elements of one-electron transfer). Thus, first term
is the kinetic energy; being transformed to the reciprocal space it leads to
formation of the band structure. The second term in (1) describes interaction
between electrons set up by the matrix elements
$I_{\nu_{1}\nu_{2}\nu_{3}\nu_{4}}$. One can introduce very different
interactions in this way, e.g.
$V_{\nu_{1}\nu_{2}}=I_{\nu_{1}\nu_{2}\nu_{1}\nu_{2}}$ is the Coulomb
interaction at different sites (responsible, for example, for the charge
ordering), $J_{\nu_{1}\nu_{2}}=-I_{\nu_{1}\nu_{2}\nu_{2}\nu_{1}}$ is the
“direct” exchange interaction (this term is obtained by the transposition
(exchange) of spin subscripts).

The next step was made by Gutzwiller Gutzwiller  and Hubbard H1 ; H3 ; H4 ,
who singled out in the model the most important part of the Coulomb
interaction – the strong repulsion of electrons at the same site
$U=I_{\nu\nu\nu\nu}$. In the case of the non-degenerate band, the Hubbard
Hamiltonian is given by

$\mathcal{H}=\sum_{\mathbf{k}\sigma}\varepsilon_{\mathbf{k}}c_{\mathbf{k}\sigma}^{\dagger}c_{\mathbf{k}\sigma}+U\sum_{i}c_{i\uparrow}^{\dagger}c_{i\uparrow}c_{i\downarrow}^{\dagger}c_{i\downarrow},$
(2)

where $\varepsilon_{\mathbf{k}}$ is the band spectrum. The Hubbard model was
widely used for studying ferromagnetism of itinerant electrons, metal-
insulator transition, and other physical phenomena connected with electron
correlations. In spite of the obvious simplicity, this model retains very rich
physics. Its strict treatment is an extremely complicated and still not solved
(in a general case) problem. Here we only mention that the first term results
in formation of the bands with the width $W$. If the repulsion between
electrons described by the second term is large, $U>W$, there appear the lower
and upper Hubbard bands and the system becomes insulating, as illustrated in
Fig. 1, for a detailed consideration we refer to a more specialized
literature, e.g., Fazekas ; Montrosi ; II1 ; Kivelson .

Figure 1: The band structure due to the first term in (2) is shown in (a). (b)
illustrates formation of lower and upper Hubbard bands in regime of large on-
site repulsion between electrons. The one-dimensional lattice with period $a$
is used for this sketch.

In Ref.H4  Hubbard introduced the many-electron representation of generalized
projection operators:

$X(\Gamma,\Gamma^{\prime})=|\Gamma\rangle\langle\Gamma^{\prime}|,\quad\sum_{\Gamma}X(\Gamma,\Gamma)=1,$
(3)

where $|\Gamma\rangle$ are the exact (many-body) eigenstates of the local part
of the Hamiltonian (on-site interaction and the crystal-field splitting in the
case of a few orbitals), see also Refs. II ; Ovchinnikov . The Hubbard
operators describe transitions between different many-electron terms.

Let us consider what kind of local states should be included into the basis
set to solve the non-degenerate model (2) at half-filling, i.e. with one
electron per site in average. Even in this situation there are not only states
with one electron having spin up $|+\rangle$ or down $|-\rangle$ – so-called
non-polar states (no additional electrons or holes), but also polar states
without electrons $|0\rangle$ (hole) and with two electrons $|2\rangle$
(double) at a given lattice site. Shubin and Vonsovsky argued that the correct
description should include not only one-electron Bloch states obtained by
Fourier transformation of $|+\rangle$ and $|-\rangle$, but also polar states.
This is exactly how correlation effects appear.

Shubin and Vonsovsky used in their works a quasi-classical approximation. It
in fact reduces to the replacement of $X$-operators by c-number functions
$\varphi_{i},\psi_{i},\Psi_{i}$, and $\Phi_{i}$, which determine the
probability amplitude of the $i$th site occupation by zero, one or two
electrons

$\displaystyle
X_{i}(+,0)\rightarrow\varphi_{i}^{\ast}\Psi_{i},X_{i}(2,-)\rightarrow\Phi_{i}^{\ast}\psi_{i},X_{i}(2,0)\rightarrow\Phi_{i}^{\ast}\Psi_{i}$

with a subsidiary normalization

$|\varphi_{i}|^{2}+|\psi_{i}|^{2}+|\Phi_{i}|^{2}+|\Psi_{i}|^{2}=1.$ (4)

In this approximation the total energy can be calculated via the variational
principle with a wave function

$\phi=\prod_{i}\left(\varphi_{i}^{\ast}X_{i}(+,0)+\psi_{i}^{\ast}X_{i}(-,0)+\Phi_{i}^{\ast}X_{i}(2,0)+\Psi_{i}^{\ast}\right)|0\rangle.$
(5)

This function mixes the Fermi- and Bose-type excitations and thereby does not
satisfy the Pauli principle. Nevertheless, the quasi-classical approximation
provides a rough description of some physical phenomena, e.g., of metal-
insulator transition in spirit of the Gutzwiller approximation Gutzwiller
(see Ref. 7 ).

Further on, the many-electron theory was developed by introducing various
representations of auxiliary (“slave”) Bose and Fermi operators for the
$X$-operators (see, e.g., i ). In particular, in connection with the theory of
high-temperature superconductivity, P.W. Anderson 633b  put forward the idea
to separate spin and charge degrees of freedom of an electron by using the
representation

$c_{i\sigma}^{\dagger}=X_{i}(\sigma,0)+\sigma
X_{i}(2,-\sigma)=s_{i\sigma}^{\dagger}e_{i}+\sigma
d_{i}^{\dagger}s_{i-\sigma}.$ (6)

Here $s_{i\sigma}^{\dagger}$ are the creation operators for neutral fermions
(spinons), and $e_{i}^{\dagger}$, $d_{i}^{\dagger}$ for charged spinless
bosons. The physical meaning of such exotic excitations with unusual
statistics can be explained as follows. Let us consider a lattice with one
electron per site with a strong Hubbard repulsion $U$, so that each site is
neutral. In the ground resonance valence bond (RVB) state, each site
participates in one bond. When the bond is broken, there appear two uncoupled
sites, each of them possessing spins 1/2. The corresponding excitations
(spinons) are not charged. At the same time, the vacancy (hole) in the system
carries a charge, but not a spin.

##  2 $s-d(f)$ exchange model

In contrast to the polar model, in the $s-d(f)$ exchange model sd ; sd1  the
existence of localized magnetic moments is not derived, but postulated; in
this case, the subsystems responsible for magnetism and electrical
conductivity are separated: itinerant “$s$-electrons” determine transport
properties, and localized “$d$-electrons” are the source of the magnetic
moments. Initially, the $s-d$ exchange model was constructed for $d$-metals,
where it gives only a semi-phenomenological picture, but its postulates turned
out to be even better justified for systems with $f$-electrons, which are
usually well localized. Later, the $s-d$ model was used to describe a number
of $d$\- and $f$-compounds, for example, so-called “heavy fermion systems”
with giant electronic specific heat. It also turned out to be surprisingly
fruitful for magnetic semiconductors.

In the western papers, the $s-d$ model is sometimes called the Vonsovsky-Zener
model in the honor of the American scientist C. Zener, who in 1951 proposed
Zener  the double-exchange mechanism of ferromagnetism. Basically it was based
on the same conception of two very different types of electrons. Localized
$d$-electrons provide their spins, which can be ferro- or
antiferromagnetically ordered, but contact (in fact, Hund’s intra-atomic)
exchange prefers ferromagnetic order, since in this situation $c$ electrons
propagate easily through the lattice as shown in Fig. 2, see also, e.g.,
Ref.Khomskii . Such concepts were used to theoretically describe giant
magnetoresistance in manganites based on LaMnO3.

Figure 2: Sketch illustrating the Zener’s double-exchange mechanism of
ferromagnetism (FM) based on the idea of itinerant ($s$) and localized ($d$)
electrons.

The Hamiltonian of the $s-d(f)$ model can be presented as

$\displaystyle\mathcal{H}=\sum_{\mathbf{k}\sigma}t_{\mathbf{k}}c_{\mathbf{k}\sigma}^{\dagger}c_{\mathbf{k}\sigma}+\sum_{\mathbf{q}}J_{\mathbf{q}}\mathbf{S}_{-\mathbf{q}}\mathbf{S}_{\mathbf{q}}$
$\displaystyle-I\sum_{i\sigma\sigma^{\prime}}(\mathbf{S}_{i}\mbox{\boldmath$\sigma$}_{\sigma\sigma^{\prime}})c_{i\sigma}^{\dagger}c_{i\sigma^{\prime}},$
(7)

where $\boldmath\sigma$ stands for the Pauli matrices, $\mathbf{S}_{i}$ are
operators of localized spins, summation in the last term runs over sites and
spins, $I$ is the parameter of the $s-d(f)$ exchange interaction, which is
assumed to be contact. Often (for example, in rare-earth metals), the
interaction between the localized spins $J_{\bf{q}}$ is in fact the indirect
Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction, see e.g., Ref. Yosida , via
conduction electrons, which appears in the second order of the perturbation
theory in $I$. However, when constructing perturbation theory, it is
convenient to include $J$ in the zero Hamiltonian.

The fate of the $s-d$ model turned out to be very interesting and difficult,
in a sense even dramatic. In the modern physical literature the $s-d$ model is
frequently called the Kondo (or Kondo lattice) model. This historical
injustice is connected with an exceptional importance of the Kondo effect
kondo , which is, maybe, the most beautiful phenomenon in condensed matter
physics. However, Kondo in its pioneering study of the anomalous magnetic
$s-d(f)$ scattering proceeded just from the second-quantization representation
Hamiltonian, which was proposed by Vonsovsky and Turov sd1 . In fact, the
field of application of the $s-d$ model goes far beyond the systems with a
considerable “Kondo” suppression of localized moments: this describes
successfully many magnetic materials.

Unlike the intra-atomic Coulomb (Hubbard) interaction, the $s-d(f)$
interaction $I$ is, as a rule, not strong. Nevertheless, it leads to essential
effects in the electronic spectrum. From a microscopic viewpoint, the $s-d(f)$
interaction can have different nature. In rare-earth metals and in a number of
their insulating and semiconducting compounds this is the intra-atomic Hund
type exchange, which is ferromagnetic. In metallic systems, frequently the
$s-d(f)$ exchange interaction is an effective one, caused by the hybridization
between the conduction $s$-bands and atomic levels of the $d(f)$-electrons. In
this case, the sign of $I$ is negative (antiferromagnetic).

Such an important physical phenomenon as the Kondo effect appears in higher
orders of perturbation theory: there occurs a singular logarithmic correction
effective (renormalized) coupling parameter in the second order, and the
corresponding contribution to resistivity in the third order. Thus, in the
case of a negative (antiferromagnetic) exchange $I$, the effective exchange
interaction becomes infinite (the infrared divergence leading to anomalies of
electron properties), so that the magnetic scattering leads to a complete
screening of localized moments 558 .

Being formulated already in the first half of XX century, the polar and
$s-d(f)$ exchange models are still used in different forms in the modern
condensed matter physics. They provide a basis for new theoretical concepts
describing physical phenomena (including exotic ones) discovered by
experimentalists. The model approaches which treat effects of strong electron
correlations in $d$\- and $f$-compounds turn out to be very useful from the
point of view of the qualitative microscopic description and provide beautiful
physics.

##  3 Historical background

S.V. Vonsovsky was born and spent his childhood in Tashkent. His mother Sofya
Ivanovna came from an old, but impoverished noble family; maternal grandfather
was a country doctor. Father, Vasily Semenovich, originated from a large
peasant family. After graduating from the gymnasium and the Faculty of Physics
and Mathematics of Moscow University, he went to teach children physics and
mathematics in Turkestan. Here he became the director of the Tashkent Women’s
Gymnasium. Before the revolution, he was an active member of the
Constitutional Democratic Party of People’s Freedom (Cadets), being personally
acquainted with the leader of this party, Petr Milyukov.

Sofya Ivanovna, being an excellent pianist, thought of making her son a
musician. Having retained his love for classical music for life, young Sergey
nevertheless followed the path of his father, choosing physics, and never
regretted it. Being already a famous scientist, he liked to repeat “Physics is
the mother of all sciences”.

S.V. Vonsovsky graduated from Leningrad University in 1932. He gratefully
recalled V.I. Smirnova, O.D. Khvolson, Yu.A. Krutikova, V.A. Fock, P.I.
Lukirsky and his other teachers at the Leningrad State University. After
graduation, he was sent to the Ural Physical-Technical Institute (later
Institute of Metal Physics) in Sverdlovsk, where he began work under the
guidance of a young professor Semen Petrovich Shubin (1908-1938), who became
his main teacher.

Figure 3: Sergey Vasiljevich Vonsovsky, Semen Petrovich Shubin, Lubov
Abramovna Shubina

Soon there appeared the articles by Shubin and Vonsovsky on the polar model of
metals and semiconductors 1 ; 2 . They were printed in the prestigious English
journal Proceedings of the Royal Society and (in more detail) in the Kharkov
journal Phys. Zs. UdSSR published in German. Unfortunately, these paper
remained largely underestimated, being in a way ahead of their time. The ideas
contained in them were then rediscovered and developed abroad (here we should
mention again the famous Hubbard model, which is a special case of the polar
model).

In 1937, Shubin was arrested for his previous participation in the Trotskyist
opposition to Stalinism and died in the Kolyma Gulag camps. Vonsovsky, who
lived in the same house as Shubin, was a witness at his last arrest. Shubin
had two children, and the third child was to be born soon. Vonsovsky took upon
himself all the worries about the family of an untimely deceased friend and
teacher. Lubov Abramovna Shubina (Shatskina) was sister of Lazar Abramovich
Shatskin (1902-1937) – one of the founders of Komsomol, which was also
arrested and killed for his struggle against Stalinism.

It would hardly have been possible to physically survive without the support
of M.N. Mikheev, a long-term director of the Institute of Metal Physics, who
himself did not escape temporary dismissal for a his principled stand.
Vonsovsky was fired from the Institute too, but his colleague Rudolf Janus
helped by employing him to an experimental laboratory for several months.

Since 1939 S.V. Vonsovsky is head of the Department of Theoretical Physics,
later deputy director of the institute. In wartime, Vonsovsky together with
Rudolf Janus carried out defense orders related to the flaw detection of
shells, and worked in Nizhny Tagil. After the war, scientific activity
gradually improved. In addition, Vonsovsky had to pay much attention to
administration. With his active participation there began intensive studies on
neutronography, low-temperature physics, radiation physics, and the theory of
dislocations in the institute.

Vonsovsky’s works in the field of the quantum theory of solids, especially the
many-electron theory of metals and semiconductors, the theory of ferro- and
antiferromagnetism, and superconductivity, received worldwide recognition. In
1946, he proposed the $s-d$ exchange model sd ; sd1  partly based on the work
with Shubin. Vonsovsky wrote:

“In parallel with the polar model, there was development of another – now so-
called $s-d$\- or $s-f$-exchange model of transition metals, which was
invented by Semen Petrovich during his oral discussion with L.D. Landau (S.P.
greatly respected and appreciated him). I finished developing this model with
my employees when Semen Petrovich was no longer with us”.

Together with his students and colleagues Vonsovsky elaborated the theory of
ferromagnetism of alloys, magnetic anisotropy and magnetostriction. His papers
were also related to the theory of superconductivity in transition metals and
alloys, the problem of the coexistence of ferromagnetism and
superconductivity. The fundamental monograph “Magnetism” magn  has become a
reference book for physicists of several generations and has been translated
into many languages. Other monographs by Vonsovsky shur ; res ; el ; iz ; vk
are widely known in Russia and abroad.

Since 1947, Vonsovsky became a professor at the Ural State University, where
he lectured for many years quantum mechanics, solid state physics, and the
theory of magnetism; for a number of years he headed the Department of
Theoretical Physics. Being a leading magnetologist in the country, Vonsovsky
was chairman of the Scientific Council on Magnetism at the USSR Academy of
Sciences for thirty years. Many major Soviet and international conferences on
magnetism (including the well-known Moscow State University conferences) were
organized under his leadership. He closely communicated with many prominent
foreign scientists. From the moment of the formation of the Ural Scientific
Center of the USSR Academy of Sciences (1971) to 1985, Vonsovsky served its
chairman.

Physics of Metals and Metallography journal was established by Vonsovsky in
1955 (he was its permanent editor-in-chief) and became popular in the physical
society of USSR. The scientific and social activities of Vonsovsky were
estimated by many government awards, including the Orders of Lenin (three
times), of the Red Star, of the Red Banner of Labor (twice) and the title of
Hero of Socialist Labor (1969). He was a laureate of the State Prize (1975,
1982), received Vavilov medal (1982), Demidov Prize (1993), was a member of
the German and Polish Academies of Sciences. Vonsovsky was a deputy of the
Supreme Soviet of the Russia Federation (1963-1971).

The most important personal quality Vonsovsky was kindness: he was simply
physically incapable of hurting a person by refusing when he was asked for
something. In Soviet times, he helped his colleagues and students a lot, not
only contributing to the solution of material problems, but also morally,
protecting them with his colossal authority. Administrative and deputy duties
left Vonsovsky little time for science, but the door to his office was always
open. Talking to him was easy and simple: he was attentive to any person,
regardless of his positions and regalia.

In the 1990s Vonsovsky gave more attention to humanitarian activities. He was
one of the founders, rector and honorary president of the Humanitarian
University of Yekaterinburg, where he himself lectured. Sergey Vasilyevich
devoted the last years of his life to work on his memoirs and the textbook
book . Initially Vonsovsky planned to write a textbook on natural sciences for
the Humanitarian University. But the book quickly went beyond this: he felt
the need to convey not only his knowledge of physics, but also the cultural
traditions of his generation. Vonsovsky managed to present many complex
concepts of modern physics in a fascinating form, but without lowering the
level: fundamental symmetries, the problem of the “grand unification” of all
fundamental interactions, the theory of supersymmetry and superstrings, not
forgetting about solid state physics, materials science and magnetism.

In the last chapter of his textbook, Sergey Vasilyevich cited the words of
Andrey Sakharov, whom he treated with great human respect. They probably also
express the scientific-philosophical and ethical position of Vonsovsky
himself: “I support the cosmological hypothesis which states that the
development of the universe is repeated in its basic features an infinite
number of times. In accordance with this, other civilizations, including more
“successful” ones, should exist an infinite number of times on the “preceding”
and the “following” pages of the Book of the Universe. Yet this should not
minimize our sacred endeavors in this world of ours, where, like faint
glimmers of light in the dark, we have emerged for a moment from the
nothingness of dark unconsciousness of material existence. We must make good
the demands of reason and create a life worthy of ourselves and of the goals
we only dimly perceive.”

##  4 Acknowlegments

We are grateful to I. Mazin and S. Ovchinnikov for useful discussions. Authors
appreciate support of the Ministry of Science and Higher Education of the
Russian Federation (theme “Quantum” No. 122021000038-7).

## References

  * (1) S. Schubin and S. Wonsowsky, Proc. Roy. Soc. A 145, 159 (1934).
  * (2) S. Schubin and S. Wonsowsky Zs. Sow. Phys. 7, 292 (1935); 10, 348 (1936).
  * (3) N. N. Bogolyubov, Lectures in Quantum Statistics, in Selected Works, Vol. 2 (Naukova Dumka: Kiev, 1970) [in Russian].
  * (4) M. C. Gutzwiller, Phys. Rev. Lett. 10, 159 (1963).
  * (5) J. Hubbard, Proc. Roy. Soc. A 276, 238 (1963).
  * (6) J. Hubbard, Proc. Roy. Soc. A 281, 401 (1964).
  * (7) J. Hubbard, Proc. Roy. Soc. A 285, 542 (1965).
  * (8) S. V. Vonsovsky, M. I. Katsnelson and A. V. Trefilov, Phys. Met. Metallogr. 76, 343 (1993).
  * (9) P. Fazekas, Lecture Notes on Electron Correlation and Magnetism (World Scientific, 1999).
  * (10) A. Montorsi (ed): The Hubbard Model - A Reprint Volume (World Scientific, Singapore, 1992).
  * (11) D. P. Arovas, E. Berg, S.A. Kivelson, and S. Raghu, Annu. Rev. Condens. Matter Phys. 13, 239 (2021).
  * (12) V. Yu. Irkhin and Yu. P. Irkhin, Electronic Structure, Correlation Effects and Properties of d- and f-metals and their Compounds (Cambridge International Science Publishing, 2007).
  * (13) V. Yu. Irkhin and Yu. P. Irkhin, phys. stat. sol. (b) 183, 9 (1994).
  * (14) S.G. Ovchinnikov, V.V. Valkov, Hubbard Operators in the Theory of Strongly Correlated Electrons (World Scientific, 2004).
  * (15) V. Yu. Irkhin, Phys. Met. Metallogr. 110, 602 (2010).
  * (16) P. W. Anderson, Int. J. Mod. Phys. B 4, 181 (1990).
  * (17) S. V. Vonsovsky, Zh. Eksp. Theor. Fiz. 16, 981 (1946).
  * (18) S. V. Vonsovsky and E.A. Turov, Zh. Eksp. Theor. Fiz. 24, 419 (1953).
  * (19) C. Zener, Phys. Rev. 81, 440 (1951).
  * (20) D.I. Khomskii and S.V. Streltsov, Chemical Reviews 121, 2992 (2021).
  * (21) K. Yosida, Theory of Magnetism (Springer, 2010)
  * (22) J. Kondo, Theory of Dilute Magnetic Alloys, in: Solid State Physics, Vol. 23, ed. F. Seitz, D. Turnbull, and H. Ehrenreich (Academic, New York, 1969), p. 183.
  * (23) A. M. Tsvelick and P. B. Wiegmann, Adv. Phys. 32, 745 (1983).
  * (24) S. V. Vonsovsky, Magnetism (Wiley: New York, 1974).
  * (25) S. V. Vonsovsky and Ya. S. Shur, Ferromagnetism (Moscow, 1948) [in Russian].
  * (26) S. V. Vonsovskii, Ferromagnetic Resonance (Pergamon: Oxford, 1966)
  * (27) S. V. Vonsovsky, Magnetism of Elementary Particles (Mir: Moscow, 1975).
  * (28) S. V. Vonsovsky, Yu. A. Izyumov, E. Z. Kurmaev, Superconductivity of Transition Metals, Their Alloys and Compounds (Springer, 1982).
  * (29) S. V. Vonsovsky, M. I. Katsnelson, Quantum Solid-State Physics (Springer, 1989).
  * (30) S.V. Vonsovsky, Modern Natural-Scientific Picture of the World (UrO RAN: Ekaterinburg, 2005) [in Russian].


# Graph Co-Ranking Approach to Keyphrase Annotation
for Digital Libraries

Adrien Bougouin
LINA, Université de Nantes,
2 rue de la houssinière
BP 92208, 44322 Nantes Cedex 03,
France
<EMAIL_ADDRESS>Florian Boudin
LINA, Université de Nantes,
2 rue de la houssinière
BP 92208, 44322 Nantes Cedex 03,
France
<EMAIL_ADDRESS>
&Béatrice Daille
LINA, Université de Nantes,
2 rue de la houssinière
BP 92208, 44322 Nantes Cedex 03,
France
<EMAIL_ADDRESS>

###### Abstract

Keyphrase annotation is the task of identifying textual units that represent
the main content of a document. Automatic keyphrase annotation methods fall
into two disjoint categories: extraction, which outputs keyphrases directly
from the textual units within the content of the document; and assignment,
which identifies keyphrases controlled by domain-specific vocabularies. Many
keyphrase extraction methods have been proposed, whereas keyphrase assignment
methods have rarely been proposed. However, keyphrase-based document indexing
for Digital Libraries requires both categories of methods. In this paper, we
apply the concept of co-ranking for simultaneously extracting and assigning
keyphrases from a graph combining document topics and domain-specific
keyphrases. This approach leverages the global context — the domain — of the
document to better extract keyphrases from the topics, and the local context —
the document — of domain-specific keyphrases to better assign keyphrases that
cannot be extracted from the document. To the best of our knowledge, this is
the first attempt to simultaneously perform both keyphrase extraction and
assignment. Experiments on three evaluation datasets show statistically
significant improvements in f1-score compared to state-of-the-art results.

##  1 Introduction

Keyphrases are words or phrases that give a synoptic picture of what is
important within a document. They are useful in many tasks such as document
indexing [gutwin1999keyphrasesfordigitallibraries], text categorization
[hulth-megyesi:2006:COLACL] or summarization [litvak2008graphbased]. However,
most documents do not provide keyphrases, and the daily flow of new documents
makes the manual keyphrase annotation impractical. As a consequence, automatic
keyphrase annotation has received special attention in the Natural Language
Processing (NLP) community and many different methods have been proposed
[hasan2014state_of_the_art].

The task of automatic keyphrase annotation consists in identifying the main
concepts, or topics, addressed in a document. Keyphrase annotation methods
fall into two broad categories: keyphrase extraction and keyphrase assignment
methods. Keyphrase extraction methods extract the most important words or
phrases occurring in a document, while assignment methods provide controlled
keyphrases from a domain-specific terminology (controlled vocabulary).

The automatic keyphrase annotation task is often reduced to the sole keyphrase
extraction task. Unlike assignment methods, extraction methods do not require
controlled vocabularies that are costly to create and to maintain.
Furthermore, they are able to identify new concepts. However, extraction
methods often output ill-formed or inappropriate keyphrases
[medelyan2008smalltrainingset], and they produce only keyphrases that actually
occur in the document.

Observations made on manually annotated keyphrases show that human annotators
both extract keyphrases from the content of the document and assign keyphrases
based on their knowledge of the domain [liu2011vocabularygap]. Here, we
propose an approach that mimics this behavior and jointly extracts and assigns
keyphrases. We use two graph representations, one for the document and one for
the domain, and apply a co-ranking algorithm to perform both keyphrase
extraction and assignment in a mutual reinforcing manner.

Along with this approach come two contributions. First, we propose a
reproducible extension to graph-based keyphrase extraction methods for
assigning keyphrases. Second, we circumvent the need of controlled vocabulary
by leveraging the reference keyphrases from training data.

We perform experiments on four datasets of different languages and domains.
Results show that our co-ranking approach significantly improves performance
over state-of-the-art approaches.

The rest of this paper is organized as follows. Section 2 presents the
previous work on automatic keyphrase annotation, section LABEL:sec:topicrankpp
describes our method, section LABEL:sec:experimental_settings presents our
experimental settings and section LABEL:sec:results discusses the results.
Finally, section 3 concludes and shows some directions for future work.

##  2 Related Work

###  2.1 Keyphrase extraction

Keyphrase extraction is the common approach to tackle the automatic keyphrase
annotation task. Previous work proposed many approaches, including binary
classification, statistical ranking and graph-based ranking of keyphrase
candidates. As our approach uses graph-based ranking, we focus on the latter.
For a detailed overview of keyphrase extraction methods, refer to
[hassan2010conundrums, hasan2014state_of_the_art].

Since the seminal work of ?), graph-based ranking approaches to keyphrase
extraction are becoming increasingly popular. The original idea behind these
approaches is to build a graph from the input document and rank its nodes
according to their importance using centrality measures.

In TextRank [mihalcea2004textrank], the input document is represented as a co-
occurrence graph in which nodes are words. Two words are connected by an edge
if they co-occur in a fixed-sized window of words. A random walk algorithm is
used to iteratively rank the words, then generate the keyphrases by
concatenating the most important words.

The random walk algorithm simulates the “voting concept”, or recommendation: a
node is important if it is connected to many other nodes, and if many of those
are important. Thus, let $G=(V,E)$ be an undirected graph with a set of
vertices $V$ and a set of edges $E$, and let $E(v_{i})$ be the set of nodes
connected to node $v_{i}$. The score $S(v_{i})$ of a vertex $v_{i}$ is
computed iteratively until convergence using the following equation:

$\displaystyle S(v_{i})=(1-\lambda)+\lambda\sum_{v_{j}\in
E(v_{i})}{\frac{S(v_{j})}{|E(v_{j})|}}$ (1)

where $\lambda$ is a damping factor defined to $0.85$ by ?).

?) first added edge weights to the random walk and further improved the graph
with co-occurrence information borrowed from similar documents. To extract
keyphrases from a document, they first look for five similar documents, then
use them to add new edges between words within the graph and reinforce the
weight of existing edges.

?) biased multiple graphs with topic probabilities. For each preliminary
retrieved topic, a graph is biased with the conditional probabilities of the
words given the topic. This method performs as many rankings as the number of
topics and gives higher importance scores to high-ranking words for as many
topics as possible. By doing so, ?) increase the topic coverage provided by
the extracted keyphrases.

?) explored further the idea of ?). Rather than using multiple graphs biased
by single topics, they added topics to graph nodes and performed co-ranking
between words and topics. Topics are important if they are connected to
important words and words are important if they are highly connected to
important topics, to important words and to other words. At the end, only
words are used for the keyphrase extraction.

?) introduced TopicRank, which is a method that clusters keyphrase candidates
into topics, ranks these topics, and selects for each important topic one
representative keyphrase. As our work extends that of ?), we present a
detailed description of their method in section LABEL:subsec:topicrank.

Graph-based methods have been well studied in the past few years. However,
none of the proposed approaches have considered the keyphrase assignment
problem nor used a similar representation to analyze and leverage reference
keyphrase relations within training data.

###  2.2 Keyphrase assignment

Keyphrase assignment provides coherent keyphrases for every document of a
specific domain using controlled vocabulary. Dissimilar to keyphrase
extraction, keyphrase assignment also aims to provide keyphrases that do not
occur within the document. This task is more difficult than keyphrase
extraction and has, therefore, seldom been employed for automatic keyphrase
annotation.

?) proposed KEA++, a shallow keyphrase assignment method. KEA++ uses a domain-
specific thesaurus to assign keyphrases to a document. First, keyphrase
candidate are selected among the n-grams of the document. N-grams that do not
match a thesaurus entry are either removed or substituted by a synonym that
matches a thesaurus entry. This candidate selection approach induces a
limitation of keyphrase assignment, refered to as keyphrase indexing by ?),
because it only assigns keyphrases if they occur within the document. Second,
KEA++ exploits the semantic relationships between keyphrase candidates within
the thesaurus as the main feature of a Naive Bayes classifier. By leveraging a
thesaurus, KEA++ achieves better performance than similar methods without a
thesaurus. However, thesauri are not readily available for most domains and
could be too general to perform accurate indexing. The application scenarii of
KEA++ are thus restricted.

####

Keyphrase extraction and keyphrase assignment methods have been studied
separately. Our work is, to our knowledge, the first attempt to perform both
tasks as a distinct one.

##  3 Conclusion and Future Work

In this paper, we have proposed a co-ranking approach to performing keyphrase
extraction and keyphrase assignment jointly. Our method, TopicCoRank, builds
two graphs: one with the document topics and one with controlled keyphrases
(training keyphrases). We design a strategy to unify the two graphs and rank
by importance topics and controlled keyphrases using a co-ranking vote.

We performed experiments on four datasets of different languages and domains.
Results showed that our approach benefits from both controlled keyphrases and
document topics, improving both keyphrase extraction and keyphrase assignment
baselines. TopicCoRank can be used to annotate keyphrases in the way of
professional indexers.

In future work, we will explore unsupervised constructions of the domain graph
(using Wikipedia data, bootstrapping techniques, etc.) and explore more
complex unification strategies in order to more accurately connect the
document and the domain graph.


# Taming the Bloch-Redfield equation: Recovering an accurate Lindblad equation
for general open quantum systems

Diego Fernández de la Pradilla<EMAIL_ADDRESS>Departamento de Física
Teórica de la Materia Condensada and Condensed Matter Physics Center (IFIMAC),
Universidad Autónoma de Madrid, E-28049 Madrid, Spain  Esteban Moreno
<EMAIL_ADDRESS>Departamento de Física Teórica de la Materia Condensada
and Condensed Matter Physics Center (IFIMAC), Universidad Autónoma de Madrid,
E-28049 Madrid, Spain  Johannes Feist<EMAIL_ADDRESS>Departamento de
Física Teórica de la Materia Condensada and Condensed Matter Physics Center
(IFIMAC), Universidad Autónoma de Madrid, E-28049 Madrid, Spain

###### Abstract

Master equations play a pivotal role in investigating open quantum systems. In
particular, the Bloch-Redfield equation stands out due to its relation to a
concrete physical environment. However, without further approximations it does
not lead to a Lindblad master equation that guarantees that the density matrix
stays completely positive, which has raised some concerns regarding its use.
This study builds on previous efforts to transform the Bloch-Redfield
framework into a mathematically robust Lindblad equation, while fully
preserving the effects that are lost within the secular approximation that is
commonly used to guarantee positivity. These previous approaches introduce two
potential deficiencies: the environment-induced energy shift can be non-
Hermitian and some decay rates can be negative, violating the assumptions of
Lindblad’s theorem. Here, we propose and evaluate straightforward solutions to
both problems. Our approach offers an effective and general procedure for
obtaining a Lindblad equation, derived from a concrete physical environment,
while mitigating the unphysical dynamics present in the Bloch-Redfield
equation.

††preprint: APS/123-QED

##  I Introduction

Inevitably, any quantum mechanical system interacts with its surroundings [1].
Therefore, it is not possible to fully understand the dynamics of a system
without properly accounting for the environment. In practice, direct
approaches to calculate the evolution with the Schrödinger equation are far
from feasible, due to the large (even infinite) number of degrees of freedom
in the environment. Consequently, it becomes imperative to develop approximate
methods that identify and preserve the relevant information. Typically, one’s
primary interest is the system, and the system-environment coupling is weak
enough that only negligible correlations appear between the two. In this
context, the framework of open quantum systems [2] provides a description of
the system’s time evolution through so-called quantum master equations (QMEs).
By incorporating the effect of the environment perturbatively, QMEs have
garnered use in most fields of physics, from condensed matter physics and
quantum optics [3, 4, 5, 6] to high-energy physics and cosmology [7, 8].
However, despite their widespread use, QMEs are perturbative in nature and
rely on approximations whose physical meaning and mathematical accuracy have
to be considered. More specifically, different QMEs can depend on varying
degrees of approximation, which can in turn invalidate one in favor of
another.

Lindblad’s theorem shows that any Markovian QME corresponding to a trace-
preserving and completely positive map can be written as a Gorini-Kossakowski-
Sudarshan-Lindblad master equation [9, 10, 2], from here on out called
Lindblad master equation (LME) for conciseness. It is given by

$\dot{\rho}=-\frac{i}{\hbar}[H+\hbar\Delta,\rho]+\sum_{ij}\frac{\gamma_{ij}}{2}\Big{(}-\\{A_{i}^{\dagger}A_{j},\rho\\}+2A_{j}\rho
A_{i}^{\dagger}\Big{)}$ (1)

and guarantees that the density matrix stays physical at all times, i.e., that
it is Hermitian, has trace equal to one, and only contains non-negative
diagonal entries. Here, $\rho$ is the system’s density matrix, $H$ is the bare
system Hamiltonian, $\hbar\Delta$ is the (Hermitian) Hamiltonian induced by
the interaction with the environment, $A_{i}$ are the collapse operators and
$\gamma_{ij}$ is a positive semidefinite matrix known as the Kossakowski
matrix. The first term on the right-hand side describes the unitary evolution
of the system, while the second term describes the dissipation induced by the
environment. The LME is commonly written in a slightly simpler form obtained
by transforming to the eigenbasis of the Kossakowski matrix $\gamma_{ij}$,
which makes the last term diagonal in $i$ and $j$. While many works use ad-hoc
LMEs with dissipation rates and operators chosen manually, this can lead to
unphysical effects, especially if the system consists of coupled subsystems
[11, 12, 13, 14, 15, 16]. A more systematic approach consists in starting from
a microscopic Hamiltonian including both system and bath, and then tracing out
the environmental degrees of freedom. Unfortunately, this procedure does not
yield an LME in general, but another QME known as the Bloch-Redfield equation
(BRE) [2], which does not guarantee positivity of the populations and can thus
lead to unphysical density matrices. To finally obtain an LME, the traditional
approach is to perform the so-called secular approximation on the BRE, in
which terms that oscillate (within the interaction picture) are assumed to
average out and are removed. However, it has been shown that this step is not
always well-justified or accurate [17, 18, 19].

Indeed, the secularization procedure completely eliminates the couplings
between populations and coherences and can therefore miss significant effects.
This is the case if, for instance, the energy spectrum has levels that are
close to each other compared to the energy scale associated to the coupling to
the environment. With the secular approximation, the resulting LME disregards
important interference effects between these states. The BRE, on the other
hand, includes the relevant couplings. Its mathematical flaws, however, have
motivated the development of improved procedures that yield an LME and
overcome the limitations of the secular approximation [20, 21, 22, 23, 24]. In
particular, a recent paper took a simple and effective step towards this goal
[22]. There, the authors propose a prescription to turn the BRE into an LME
that works well for simple system-environment interaction Hamiltonians within
the rotating wave approximation. We further extended those results in [25] by
incorporating counter-rotating terms into more complex system-environment
interactions. Still, the master equation used in [25] can lead to two
undesired possibilities: (i) non-Hermitian, off-diagonal couplings, depending
on the interplay between the energy structure of the system and the resonances
of the bath, and (ii) a non-positive Kossakowski matrix if there are cross-
correlations between bath operators. Either option invalidates the
requirements for an LME, and both are a direct consequence of the application
of the methodology of [22] to complex systems and environments.
Notwithstanding, under conditions of geometrical symmetry and with a careful
treatment of the system level structure, these difficulties can be overcome.
Indeed, we took this approach in [25] to reveal strong Casimir-Polder-induced
state mixing and counterintuitive state protection in a system involving a
hydrogen atom and a dielectric nanoparticle.

The goal of the present article is to present a way to deal with arbitrary
setups by modifying the prescription of [22]. We introduce and compare several
alternatives that resolve the issues of non-Hermitian couplings and non-
positive semidefiniteness (NPS) of the Kossakowski matrix, while maintaining
an accurate description of the dynamics, turning the BRE into a general and
precise LME. On one hand, we find that the non-Hermiticity of the energy shift
can be resolved by using the arithmetic instead of the geometric mean in the
symmetrization procedure for the environment-induced energy shift. On the
other hand, the NPS issue can be resolved effectively by discarding the
negative eigenvalues of the Kossakowski matrix. We explain the details of the
procedures and justify them with model examples that capture the essence of
each problem. The Python library QuTiP [26] has been used throughout this work
for the simulations that illustrate and support our claims. The article is
organized as follows: In section II, we discuss the BRE and its deficiencies
explicitly, along with the explanation of the methodology developed in [22] to
turn it into an LME. The non-Hermiticity issue is addressed in section III,
where two possibilities are compared and, eventually, one is identified as the
best. The analysis of the NPS of the Kossakowski matrix is presented in
section IV. Finally, we summarize the results of our exploration in section V.

##  II Theory

In this section, we write the BRE and elucidate how the method proposed in
[22] transforms it into an LME, subject to certain conditions. Deviation from
these conditions implies that the resulting equation is not an LME, which
motivates the development of the extensions developed in this article.

###  II.1 Bloch-Redfield master equation

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1: Breakdown of population positivity with the BRE. Top panels:
Population of the states of a three-level system, as a function of the
coupling strength, $g/\kappa$, and time, $t\kappa$. Blue and red indicate
negative and positive population, respectively. Bottom panels: Same
populations at $g=\kappa$ (dashed lines in the top panels). The blue region
represents negative, unphysical population, and the dashed line offers a
visual aid to see how the population of the states $|1\rangle$ and $|2\rangle$
add up to more than one.  Figure 2: Side-by-side comparison of the system’s
level structure, $J(\omega)$ and $\lambda(\omega)$. The parameters for the
figure are $\hbar\omega_{M}=1$ eV, $\hbar\kappa=0.1$ eV, $\hbar g=0.1$ eV.

The BRE serves as starting point in numerous derivations of LMEs, due to its
physically motivated derivation. We write it for a general system-bath
interaction $H_{\mathrm{sb}}=\sum_{\alpha}A_{\alpha}B_{\alpha}$, where
$A_{\alpha}$ and $B_{\alpha}$ represent system and bath operators,
respectively:

$\displaystyle\dot{\rho}=-\frac{i}{\hbar}[H$ $\displaystyle,\rho]$
$\displaystyle+\sum_{ij}\bigg{\\{}$
$\displaystyle-i\Big{(}\Lambda_{ij}(\omega_{j})\sigma_{i}^{\dagger}\sigma_{j}\rho-\rho\Lambda_{ij}(\omega_{i})\sigma_{i}^{\dagger}\sigma_{j}\Big{)}$
$\displaystyle+i\Big{(}\Lambda_{ij}(\omega_{j})-\Lambda_{ij}(\omega_{i})\Big{)}\sigma_{j}\rho\sigma_{i}^{\dagger}$
$\displaystyle-\frac{1}{2}\Big{(}\Gamma_{ij}(\omega_{j})\sigma_{i}^{\dagger}\sigma_{j}\rho+\rho\Gamma_{ij}(\omega_{i})\sigma_{i}^{\dagger}\sigma_{j}\Big{)}$
$\displaystyle+\frac{1}{2}\Big{(}\Gamma_{ij}(\omega_{j})+\Gamma_{ij}(\omega_{i})\Big{)}\sigma_{j}\rho\sigma_{i}^{\dagger}\bigg{\\}}.$
(2)

In the above equation, $\rho$ and $H$ represent the system’s density matrix
and Hamiltonian, the indices $i$ and $j$ are combined indices that denote
transitions, with $\sigma_{j}=|n_{j}\rangle\langle m_{j}|$ being the $j$th
transition operator between eigenstates $|n_{j}\rangle$ and $|m_{j}\rangle$ of
$H$, $\omega_{j}=(E_{m_{j}}-E_{n_{j}})/\hbar$ is the $j$th transition
frequency, and the functions $\Gamma_{ij}$ and $\Lambda_{ij}$ are defined as

$\displaystyle\Gamma_{ij}(\omega)$
$\displaystyle=\sum_{\alpha\beta}(A_{\alpha})_{i}^{*}\gamma_{\alpha\beta}(\omega)(A_{\beta})_{j},$
(3a) $\displaystyle\Lambda_{ij}(\omega)$
$\displaystyle=\sum_{\alpha\beta}(A_{\alpha})_{i}^{*}\lambda_{\alpha\beta}(\omega)(A_{\beta})_{j}.$
(3b)

Here, $(A_{\beta})_{j}=\mathrm{Tr}\\{\sigma_{j}^{\dagger}A_{\beta}\\}$ is the
transition matrix element associated with the $j$th transition, and the real
functions $\gamma_{\alpha\beta}$ and $\lambda_{\alpha\beta}$ are given by

$\frac{\gamma_{\alpha\beta}(\omega)}{2}+i\lambda_{\alpha\beta}=\int_{0}^{\infty}\mathrm{d}\tau\
\langle B_{\alpha}(\tau)B_{\beta}(0)\rangle e^{i\omega\tau},$ (4a) where the
imaginary and real parts are related through a Kramers-Kronig-type relation:
$\lambda_{\alpha\beta}(\omega)=\frac{1}{2\pi}\mathrm{P}\int\mathrm{d}\omega^{\prime}\
\frac{\gamma_{\alpha\beta}(\omega^{\prime})}{\omega-\omega^{\prime}}.$ (4b)

The P denotes the principal value and $\bm{\gamma}$ is a non-diagonal and
positive definite tensor-valued function. It is related to the so-called
spectral density $\mathbf{J}$ through $\bm{\gamma}=2\pi\mathbf{J}$, which
characterizes the coupling between the system and the bath.

Roughly speaking, $\Lambda_{ij}$ and $\Gamma_{ij}$ can loosely be associated
with the environment-induced energy shifts and decay rates, respectively.
However, this distinction is not completely clear within the BRE itself, as
the most natural interpretation depends on the specific arrangement of the
equation. For instance, in Eq. (missing)2, we would recover an equation of
Lindblad-form if the following conditions were met: (i) the second line should
become a commutator, (ii) the third line should disappear, (iii) the fourth
line should become an anticommutator, and (iv) the last line should become the
“refilling” or “quantum jump” term, with the same prefactor as the previous
line. Then, the second line would straightforwardly become the energy
correction, and the last two lines would embody Lindblad-like dissipators,
with clearly distinct roles for $\Lambda$ and $\Gamma$. These requirements are
not directly satisfied in Eq. (missing)2 due to $\Lambda_{ij}$ and
$\Gamma_{ij}$ being evaluated at two different frequencies in each line,
$\omega_{j}$ and $\omega_{i}$, but the approach followed in [22] achieves
conditions (i) through (iv) by means of an additional, well-justified
approximation. In this paper we explore the same procedure. Nevertheless, it
should be noted that the conditions above are not the only way to transform
the BRE into a Lindblad-like equation. In fact, with the following commutator
identity,

$aO_{1}O_{2}+bO_{2}O_{1}=\frac{a+b}{2}\\{O_{1},O_{2}\\}+\frac{a-b}{2}[O_{1},O_{2}],$
(5)

Eq. (missing)2 can be reformulated as

$\displaystyle\dot{\rho}$ $\displaystyle=-\frac{i}{\hbar}[H,\rho]$
$\displaystyle+\sum_{ij}\bigg{\\{}-\frac{K^{+}_{ij}(\omega_{j})-K^{-}_{ij}(\omega_{i})}{2}[\sigma_{i}^{\dagger}\sigma_{j},\rho]$
$\displaystyle+\frac{K^{+}_{ij}(\omega_{j})+K^{-}_{ij}(\omega_{i})}{2}\Big{(}-\\{\sigma_{i}^{\dagger}\sigma_{j},\rho\\}+2\sigma_{j}\rho\sigma_{i}^{\dagger}\Big{)}\bigg{\\}},$
(6)

where $K^{\pm}_{ij}(\omega_{j})=\Gamma_{ij}(\omega_{j})/2\pm
i\Lambda_{ij}(\omega_{j})$. In this form, the BRE looks like an LME, with both
$\Lambda$ and $\Gamma$ appearing in the energy correction term and in the
dissipators, blurring their individual role. Note that even written as Eq.
(missing)6, the BRE is not an LME, because the corresponding Kossakowski
matrix is not positive semidefinite.

We conclude this subsection on the BRE by showing in Fig. 1 that Eq.
(missing)2 can actually predict non-positive density matrices. To produce the
figure, we have simulated a simple system composed of a three-level system
coupled to a reservoir with a Lorentzian spectral density, schematically
depicted in Fig. 2. The system eigenstates are $|0\rangle$, $|1\rangle$ and
$|2\rangle$, with two allowed lowering transitions from the excited states
$|1\rangle$ and $|2\rangle$ to $|0\rangle$, characterized by transition
frequencies $\omega_{1}$ and $\omega_{2}$, respectively. Let us denote the
transition operators with non-zero matrix element as

$\displaystyle\sigma_{1}$ $\displaystyle=|0\rangle\langle 1|,$
$\displaystyle\sigma_{2}$ $\displaystyle=|0\rangle\langle 2|,$ (7a)
$\displaystyle\sigma_{1^{\prime}}$
$\displaystyle=\sigma_{1}^{\dagger}=|1\rangle\langle 0|,$
$\displaystyle\sigma_{2^{\prime}}$
$\displaystyle=\sigma_{2}^{\dagger}=|2\rangle\langle 0|,$ (7b)

For concreteness, we set all transition matrix elements equal to each other
and absorb them in the coupling to the environment. Its spectral density is
given by,

$J(\omega)=\frac{g^{2}}{\pi}\frac{\kappa/2}{(\omega-\omega_{M})^{2}+(\kappa/2)^{2}},$
(8a) where $g$ measures the coupling strength, and $\omega_{M}$ and $\kappa$
are the position and width of the peak of the mode described by the
Lorentzian. The corresponding integral $\lambda$ is analytical and expressed
as
$\lambda(\omega)=\mathrm{P}\int_{\mathbb{R}}\mathrm{d}\omega^{\prime}\frac{J(\omega^{\prime})}{\omega-\omega^{\prime}}=g^{2}\frac{\omega-\omega_{M}}{(\omega-\omega_{M})^{2}+(\kappa/2)^{2}}.$
(8b)

The parameters are set to $\hbar\omega_{M}=1$ eV, $\hbar\kappa=0.1$ eV,
$\hbar\omega_{1}=0.75$ eV, $\hbar\omega_{2}=1.35$ eV, and the initial state to
$|\psi(0)\rangle=(|1\rangle-|2\rangle)/\sqrt{2}$. The population dynamics are
depicted in Fig. 1 as a function of the coupling strength and time. In the top
panels, negative and positive values are represented in blue and red,
respectively. The bottom panels are a cut along the dashed lines, at
$g=\kappa$. In the first column, Fig. 1a and Fig. 1d reveal a very prominent
negative population in the groundstate $|0\rangle$, an unmistakable sign of
non-physicality. The preservation of the trace under BRE evolution implies
that the populations of $|1\rangle$ and $|2\rangle$ together must exceed 1.
This compensation is evident in the slight initial bumps observed in Fig. 1e
and Fig. 1f, both above the 0.5 line.

###  II.2 Prescription to obtain an LME

In this subsection, we explain the procedure devised in [22] to turn the BRE
into a Lindblad-like equation. We begin by noting that there are two terms in
every line in the sum of Eq. (missing)2, evaluated at $\omega_{j}$ and
$\omega_{i}$, respectively. The asymmetry in the indices $i$ and $j$ is the
reason why the conditions (i) to (iv) stated in the previous subsection are
not fulfilled. To overcome this difficulty, it was proposed to replace
$\Gamma_{ij}(\omega_{j})$ and $\Gamma_{ij}(\omega_{i})$ by their geometric
mean,
$\tilde{\Gamma}_{ij}=\sqrt{\Gamma_{ij}(\omega_{j})}\sqrt{\Gamma_{ij}(\omega_{i})}$,
and likewise for $\Lambda_{ij}$. This adjustment, whose accuracy is justified
below, removes the difference in the roles of $\omega_{j}$ and $\omega_{i}$.
We designate this prescription as $g\Lambda-g\Gamma$, as both magnitudes are
replaced by geometric averages. To be more exact, [22] deals with systems with
one system-bath interaction, while this explicit formulation, first presented
in [25], accounts for multiple interactions, including cross-correlations
between bath operators.

After applying the replacement of the geometric mean in Eq. (missing)2, the
master equation is

$\displaystyle\dot{\rho}=$
$\displaystyle-\frac{i}{\hbar}[H,\rho]-i\sum_{ij}\tilde{\Lambda}_{ij}[\sigma_{i}^{\dagger}\sigma_{j},\rho]$
$\displaystyle+\sum_{ij}\frac{\tilde{\Gamma}_{ij}}{2}\Big{(}-\\{\sigma_{i}^{\dagger}\sigma_{j},\rho\\}+2\sigma_{j}\rho\sigma_{i}^{\dagger}\Big{)}.$
(9)

The justification for the replacement lies in the consideration of a slowly
varying spectral density, where
$\mathbf{J}(\omega+\Lambda)\simeq\mathbf{J}(\omega)$ and
$\mathbf{J}(\omega+\Gamma)\simeq\mathbf{J}(\omega)$. Under these conditions,
every term in Eq. (missing)2 with
$|\omega_{i}-\omega_{j}|<\max\\{\Gamma_{ij},\Lambda_{ij}\\}$ satisfies
$\Gamma_{ij}(\omega_{i})\simeq\Gamma_{ij}(\omega_{j})\simeq\tilde{\Gamma}_{ij}$
and
$\Lambda_{ij}(\omega_{i})\simeq\Lambda_{ij}(\omega_{j})\simeq\tilde{\Lambda}_{ij}$,
and the replacement is trivially precise. Conversely, when
$|\omega_{i}-\omega_{j}|>\max\\{\Gamma_{ij},\Lambda_{ij}\\}$ the accuracy of
the replacement decreases, but the effect of these terms becomes negligible.
Consequently, the preciseness of the replacement is irrelevant, and these
terms can even be discarded through secularization.

After the replacement, Eq. (missing)9 is an LME if (i) the energy correction,
$\Delta=\sum_{ij}\tilde{\Lambda}_{ij}\sigma_{i}^{\dagger}\sigma_{j}$, is
Hermitian and (ii) the Kossakowski matrix, $\tilde{\Gamma}_{ij}$, is positive
semidefinite. In the examples from [22] both conditions were automatically
satisfied. As for the system explored in [25], the non-Hermitian terms were
removed through an accurate secularization, and the geometric symmetry of the
setup implied a positive semidefinite Kossakowski matrix. In general, however,
these requirements are not guaranteed and an LME is not obtained. In the
following, we elaborate on strategies to address both problems.

##  III Non-Hermiticity

Here, we discuss how non-Hermiticity can appear in the energy correction
Hamiltonian and propose two solutions that ensure Hermiticity.

###  III.1 Description of the problem

As mentioned above, for the master equation to be an LME,
$\Delta=\sum_{ij}\tilde{\Lambda}_{ij}\sigma_{i}^{\dagger}\sigma_{j}$ must be
Hermitian, which is automatically fulfilled when $\tilde{\Lambda}$ is
Hermitian. In principle, it is possible for $\Delta$ to be Hermitian even when
$\tilde{\Lambda}$ is not, but this would require fortuitous cancellations
between system and environment properties that do not occur in general. Thus,
we focus on approaches to convert $\tilde{\Lambda}$ into a Hermitian matrix.

To illustrate the origin of non-Hermitian terms in $\Delta$, we consider the
same three-level system coupled to a Lorentzian bath described before (cf.
Fig. 2). Despite its simplicity, this model is enough to reveal the source of
non-Hermitian terms in $\Delta$. Simultaneously, the corresponding Kossakowski
matrix is positive semidefinite, due to the absence of environment cross-
correlations, allowing us to isolate the non-Hermiticity. According to Eq.
(missing)9, the energy shift for this system is

${\Delta}=\begin{pmatrix}\tilde{\Lambda}_{1^{\prime}1^{\prime}}+\tilde{\Lambda}_{2^{\prime}2^{\prime}}&0&0\\\
0&\tilde{\Lambda}_{11}&\tilde{\Lambda}_{12}\\\
0&\tilde{\Lambda}_{21}&\tilde{\Lambda}_{22}\end{pmatrix},$ (10)

where
$\tilde{\Lambda}_{12}=\tilde{\Lambda}_{21}=\sqrt{\lambda(\omega_{1})}\sqrt{\lambda(\omega_{2})}$.
Notice that $\lambda(\omega)$ changes sign at $\omega_{M}$. Thus, if
$\omega_{1}$ and $\omega_{2}$ lie on opposite sides of $\omega_{M}$, then
$\tilde{\Lambda}_{12}=\tilde{\Lambda}_{21}=i|\tilde{\Lambda}_{12}|$, and
$\Delta$ is non-Hermitian.

The situation described above illustrates how non-Hermiticity can arise in the
energy correction operator. When this occurs, Eq. (missing)9 is not truly an
LME, and the positivity of $\rho$ is no longer guaranteed. Not only that, but
the Hermiticity of $\rho$ is also not preserved and non-negligible imaginary
contributions arise in the populations. As an example, we show in Fig. 3 the
prediction of Eq. (missing)9, with parameters set to $\hbar\omega_{M}=1$ eV,
$\hbar\kappa=0.1$ eV, $g=\kappa$ and $\omega_{i}=\omega_{M}+\delta_{i}$,
choosing $\delta_{2}=-\delta_{1}=\delta$, and $|\psi(0)\rangle=|1\rangle$ as
the initial state. In Fig. 3a, we depict the phase of the “complex population”
of $|1\rangle$. While it starts at 0 for short time scales, it undergoes rapid
changes and displays multiple complete turns. Thus, the Hermiticity of $\rho$
is quickly compromised due to the off-diagonal, non-Hermitian couplings
described above. In the middle and bottom panels, we plot the population of
state $|2\rangle$, showing only its real part, as the imaginary part is
negligible. Nevertheless, Fig. 3b shows positive and negative populations in
red and blue, respectively, and the negative portions are definitely not
negligible. In the bottom plot, a cut at $\delta=\kappa$ accentuates the
negativity, indicated with the blue region.

Despite the observations in Fig. 3, note that these deficiencies do not
necessarily mean that Eq. (missing)9 with a non-Hermitian $\Delta$ is always
inaccurate. Similarly to the BRE, it can remain precise if the underlying
approximations are reasonable [19]. The parameters have indeed been chosen to
highlight the difficulties that can arise from Eq. (missing)9. Still, it
remains true that the attempt at a solution to the BRE problems has lead to a
new, arguably worse issue: populations become complex rather than simply
negative. In the following subsections, we describe and compare several
alternatives designed to eliminate the non-Hermiticity problem.

(a)

(b)

(c)

Figure 3: Breakdown of the $g\Lambda-g\Gamma$ prescription of [22]. (a) Non-
trivial complex phase of the “population” of $|1\rangle$ owing to the loss of
Hermiticity of Eq. (missing)10. (b) Population of $|2\rangle$ obtained with
Eq. (missing)9 with $\omega_{1}=\omega_{M}+\delta_{1}$,
$\omega_{2}=\omega_{M}+\delta_{2}$ and $\delta=-\delta_{1}=\delta_{2}$ (very
close to being real). Blue and red represent positive and negative values,
respectively. (c) Population of $|2\rangle$ at $\delta=\kappa$. The blue
region indicates negative, unphysical population.

###  III.2 Approach 1: Degenerate levels

It was argued in the introduction that the off-diagonal terms are relevant
when the energetic difference between states $|1\rangle$ and $|2\rangle$ in
Fig. 2 is comparable to the energy corrections. For that reason, a possible
solution to the Hermiticity problem would be to ignore the energy difference
between these states while deriving the master equation, and then add it back
as a correction afterwards. In the illustrative system of Fig. 2, this
corresponds to computing the parameters in the master equation as if
$\omega_{1}=\omega_{2}=\bar{\omega}=(\omega_{1}+\omega_{2})/2$, and then
adding the actual energy difference to the effective Hamiltonian. We follow up
on the notation introduced previously and refer to this method as
$d\Lambda-d\Gamma$, as the close-lying energy levels of the system are assumed
to be degenerate from the start in the master equation derivation. With this
procedure, the energy shift is

$\Delta^{d}=\begin{pmatrix}2\lambda(-\bar{\omega})&0&0\\\
0&\lambda(\bar{\omega})&\lambda(\bar{\omega})\\\
0&\lambda(\bar{\omega})&\lambda(\bar{\omega})\end{pmatrix},$ (11)

and, since $\lambda(\bar{\omega})$ is real, $\Delta^{d}$ is automatically
Hermitian. This replacement is numerically accurate when $J(\omega_{2})$ and
$J(\omega_{1})\simeq J(\bar{\omega})$. However, we do not expect it to perform
very well whenever either $J(\omega_{2})$ or $J(\omega_{1})\not\simeq
J(\bar{\omega})$ because, although $\Delta^{d}_{12}$ becomes negligible, the
dissipative terms are also evaluated at $\bar{\omega}$, significantly
different from $\omega_{1}$ and $\omega_{2}$, potentially leading to
inaccurate decay rates.

This procedure is easily applicable to systems exhibiting clearly structured
energy spectra with several narrow sets of levels, because one single
frequency can be reliably assigned to each set. Within each set, the off-
diagonal couplings would resemble the Hermitian expression from Eq.
(missing)11, and although the couplings between different sets would still be
susceptible to the non-Hermiticity issue, these can in principle be eliminated
through secularization. However, for more complicated level structures, where
the energy levels are not distinctly arranged in narrow sets, assigning the
same frequency to several states becomes a less straightforward task.
Exquisite care must then be taken in order to avoid losing pertinent
information about the frequency dependence of $J$ and $\lambda$, which impacts
the precise values of the energy shifts and decay rates. Additionally, this
means that in general the procedure has to be performed manually, or at least
under strict human supervision.

###  III.3 Approach 2: Arithmetic mean

We devise a second approach for addressing the non-Hermitian couplings that
involves replacing the geometric mean discussed in section II.2 with the
arithmetic mean solely in the energy correction, while keeping the geometric
mean in the decay terms:

$\displaystyle\tilde{\Lambda}_{ij}$
$\displaystyle=\frac{\Lambda_{ij}(\omega_{i})+\Lambda_{ij}(\omega_{j})}{2},$
(12a) $\displaystyle\tilde{\Gamma}_{ij}$
$\displaystyle=\sqrt{\Gamma_{ij}(\omega_{i})}\sqrt{\Gamma_{ij}(\omega_{j})}.$
(12b)

Accordingly, we designate this prescription as $a\Lambda-g\Gamma$. This
adjustment allows for the refactoring of the energy shift terms in Eq.
(missing)2 into a commutator. It is worth noting that if the arithmetic mean
were used in the decay terms as well, the Kossakowski matrix for the system in
Fig. 2 would give rise to a negative decay rate, approximately given by

$\pi\left(J_{+}-\sqrt{J_{+}^{2}+J_{-}^{2}}\right)<0,$ (13)

where $J_{\pm}=J(\omega_{1})\pm J(\omega_{2})$. Similar features appear also
in other, more complex systems.

With this option, the energy correction for the three-level system becomes

$\Delta^{a}=\begin{pmatrix}\lambda(-\omega_{1})+\lambda(-\omega_{2})&0&0\\\
0&\lambda(\omega_{1})&\frac{\lambda(\omega_{1})+\lambda(\omega_{2})}{2}\\\
0&\frac{\lambda(\omega_{1})+\lambda(\omega_{2})}{2}&\lambda(\omega_{2})\end{pmatrix}.$
(14)

Unlike the previous fix, the arithmetic mean naturally preserves the
information about the dependence of $J$ and $\lambda$ with $\omega$. Hence,
this solution to the non-Hermiticity problem is more straightforwardly
applicable to more complex systems. We finish by noting that this approach
yields an expression for $\Delta$ similar to the equivalent term in the second
line of Eq. (missing)6, an alternative form of the BRE. However, in that
equation, the $\Lambda_{ij}$ functions appear together with an additional
$i(\Gamma_{ij}(\omega_{i})-\Gamma_{ij}(\omega_{j}))/4$.

###  III.4 Accuracy comparison

(a)

(b)

(c)

(d)

Figure 4: Population of the state $|2\rangle$ calculated with (a) the exact,
(b) $a\Lambda-g\Gamma$ and (c) $d\Lambda-d\Gamma$ methods, as a function of
$\delta$. (d) $d\Lambda-nd\Lambda$, designed to fix the $d\Lambda-d\Gamma$
method by keeping the full information of the system frequencies. The gray
dashed line is proportional to $1/J(\omega_{M}+\delta)$, indicating the
lifetime of $|2\rangle$ when the excited states do not interfere. The
parameters are $\hbar\omega_{M}=1$ eV, $\hbar\kappa=0.1$ eV, $g=0.05\kappa$,
$\omega_{1}=\omega_{M}+\delta_{1}$ and $\omega_{2}=\omega_{M}+\delta_{2}$,
with $-\delta_{1}=\delta_{2}=\delta$. The initial state is
$|\psi(0)\rangle=(|1\rangle-|2\rangle)/\sqrt{2}$.

In this section, we evaluate the dynamics of the three-level system depicted
in Fig. 2 using several methods. We then compare and identify the best option
among the approaches from sections III.2 and III.3, referred to as
$d\Lambda-d\Gamma$ and $a\Lambda-g\Gamma$, respectively. Additionally, we
include an exact method that exploits the well-known equivalence between a
Lorentzian bath and a single bosonic mode, $a$, coupled to a flat (Markovian)
environment, from the perspective of the system [27]. Therefore, the system
plus mode system can be exactly described through

$\dot{\rho}=-\frac{i}{\hbar}[H+\hbar\omega_{M}a^{\dagger}a,\rho]+\frac{\kappa}{2}\Big{(}-\\{a^{\dagger}a,\rho\\}+2a\rho
a^{\dagger}\Big{)},$ (15)

and we may obtain the reduced system density matrix by tracing over the
bosonic degree of freedom. The approximate methods can then be compared to the
exact solution.

For the comparison, we choose the following parameters: $\hbar\omega_{M}=1$
eV, $\hbar\kappa=0.1$ eV, $g=0.05\kappa$, $\omega_{1}=\omega_{M}+\delta_{1}$
and $\omega_{2}=\omega_{M}+\delta_{2}$, with $-\delta_{1}=\delta_{2}=\delta$.
The initial state is $|\psi(0)\rangle=(|1\rangle-|2\rangle)/\sqrt{2}$. We show
the population of state $|2\rangle$ as calculated with the exact,
$a\Lambda-g\Gamma$ and $d\Lambda-d\Gamma$ methods in Fig. 4a, Fig. 4b and Fig.
4c, respectively. When the states $|1\rangle$ and $|2\rangle$ are close to
degeneracy ($\delta/\kappa\approx 0$), both the $a\Lambda-g\Gamma$ and the
$d\Lambda-d\Gamma$ methods agree with the exact result, as expected. However,
as $\delta/\kappa$ increases, the exact and $a\Lambda-g\Gamma$ color maps
reveal a gradual and clear increase in the decay timescale, marked with dashed
lines for ease of comparison between plots. In contrast, this effect is not
present in the $d\Lambda-d\Gamma$ method.

It should be noted that the reason for the disparity of the $d\Lambda-d\Gamma$
method compared to the first two lies in the decay terms, rather than in
$\Delta$. This is because the $d\Lambda-d\Gamma$ implementation approximates
the levels as degenerate during the full derivation of the master equation,
and only afterwards includes the finer structures. Hence, both $\lambda$ and
$J$ are evaluated at $\bar{\omega}$ and, with the parameters used above, the
decay rate far from resonance is still given by $2\pi J(\omega_{M})$ in Fig.
4c, while it becomes smaller for increasing values of $\delta/\kappa$ in Fig.
4a and Fig. 4b. This can be corrected by not replacing $\omega_{1}$ and
$\omega_{2}$ with $\bar{\omega}$ in the decay part, as shown in Fig. 4d, where
the same dependence of the system’s lifetime on $\delta/\kappa$ is observed.
With our notation, we can term this “mixed approach” $d\Lambda-g\Gamma$.
Judging from Fig. 4d, we conclude that, in the situation considered here, both
$a\Lambda-g\Gamma$ and $d\Lambda-g\Gamma$ work well. Still, the energy shift
in $d\Lambda-g\Gamma$ assumes no natural detuning between the excited states,
while in reality such frequency difference is non-zero. Therefore, it could be
the case that some concrete scenarios might compromise the accuracy of
$d\Lambda-g\Gamma$ compared to $a\Lambda-g\Gamma$, where the detuning between
excited states is automatically included with no additional consideration,
thus ensuring a more robust description of the emitter’s dynamics.

Considering these findings, we conclude that the $a\Lambda-g\Gamma$ approach
from section III.3 provides the best solution. It successfully reproduces the
exact solution in the weak coupling regime, addresses the non-Hermiticity
issue associated with the geometric mean prescription, and seamlessly
incorporates the energetic features of the system without additional
supervision. This makes the $a\Lambda-g\Gamma$ approach a robust and effective
method for describing the energy shift of complex quantum systems. We have
also tried several other options, such as manually ensuring that
$\tilde{\Lambda}_{ij}=\tilde{\Lambda}_{ji}^{*}$ by arbitrarily conjugating one
of them. For conciseness, these methods are not discussed here, as they are
worse and more arbitrary than the arithmetic mean procedure.

##  IV Non-positive semidefiniteness

We continue our analysis by clarifying the situations that ensure the positive
semidefiniteness of the Kossakowski matrix. Notably, in the case of a single
system operator $A$ coupling to one bath operator $B$, such that
$H_{\mathrm{sb}}=AB$, the ensuing Kossakowski matrix is positive semidefinite
by construction, as demonstrated in [22]. Indeed, the procedure outlined in
subsection II.2 yields a dissipative term of the form
$\frac{1}{2}\left(-\left\\{\Sigma^{\dagger}\Sigma,\rho\right\\}+2\Sigma\rho\Sigma^{\dagger}\right)$,
with a single “collective” decay operator
$\Sigma=\sum_{j}A_{j}\sqrt{\gamma(\omega_{j})}\sigma_{j}$. Then, the only non-
zero eigenvalue of the Kossakowski matrix is a positive 1.

Nevertheless, the scenario becomes more intricate when allowing for multiple
system-bath interaction terms,
$H_{\mathrm{sb}}=\sum_{\alpha}A_{\alpha}B_{\alpha}$, with the potential
appearance of bath cross-correlations, i.e., off-diagonal elements in Eq.
(missing)4a. In this case, the Kossakowski matrix after the geometric mean
replacement can often be non-positive semidefinite. Yet, specific conditions
can guarantee non-negative eigenvalues. For instance, if the cross-
correlations vanish and each system transition involves only one of the
$A_{\alpha}$ operators, the situation reduces to a sum over several single
system-bath interactions. These requirements can only be met in very
particular configurations, like the one considered in [25], which can be a
significant limitation.

In the following, we perform a structured study of the size and impact of the
negative eigenvalues of the Kossakowski matrix. First, we focus our analysis
on a specific system of small size, and observe that, with the prescription
$a\Lambda-g\Gamma$ from earlier, the negative eigenvalues can be discarded
without significant accuracy loss. Next, we examine a large sample of randomly
generated test systems to confirm the previous conclusion that the negative
eigenvalues can be neglected. At the moment, we do not have a rigorous
justification as to why taking the geometric mean should imply that, but the
statistical study seems to confirm it to be the case.

###  IV.1 Particular case

Here, we generate a particular system and explore the effectiveness of various
approaches aimed at achieving the positivity of the corresponding Kossakowski
matrix. The system is generated according to the following procedure. We
consider a generic $L$-level system coupled to the environment through $M$
interaction terms, such that
$H_{\mathrm{sb}}=\sum_{\alpha=1}^{M}A_{\alpha}B_{\alpha}$. The spectral
density of the bath is, then, represented as a positive semidefinite $M\times
M$ matrix-valued function. To describe the properties of the environment, we
exploit the equivalence between a generic bath and a particular set of $N$
coupled and lossy modes $a_{\beta}$, as discussed in [28]. It can be shown
that the associated spectral density is

$\mathbf{J}(\omega)=\frac{\hbar}{\pi}\mathbf{g}\cdot\textrm{Im}\left\\{(\mathbf{h}-\omega)^{-1}\right\\}\cdot\mathbf{g}^{T},$
(16a) where particular matrices $\mathbf{h}$ ($N\times N$) and $\mathbf{g}$
($M\times N$), along with an integer $N$, can be chosen to accurately fit any
spectral density. The corresponding Lamb shift integral is given by
$\bm{\lambda}(\omega)=-\hbar\mathbf{g}\cdot\textrm{Re}\left\\{(\mathbf{h}-\omega)^{-1}\right\\}\cdot\mathbf{g}^{T}.$
(16b)

More specifically, $h_{\beta\beta}=\Omega_{\beta\beta}-i\kappa_{\beta}/2$
represents the frequencies and losses of the modes,
$h_{\beta\gamma}=\Omega_{\beta\gamma}$ contains the couplings between discrete
modes, and $g_{\alpha\beta}$ determines the coupling strength between the mode
$a_{\beta}$ and the system through $A_{\alpha}$. Additionally, using the
results from [28], the exact dynamics of the system is also available through

$\displaystyle\dot{\rho}_{\mathrm{sm}}={}$
$\displaystyle-\frac{i}{\hbar}\left[H_{\mathrm{sm}},\rho_{\mathrm{sm}}\right]$
$\displaystyle+\sum_{\beta=1}^{N}\frac{\kappa_{\beta}}{2}\left(-\left\\{a_{\beta}^{\dagger}a_{\beta},\rho_{\mathrm{sm}}\right\\}+2a_{\beta}\rho_{\mathrm{sm}}a_{\beta}^{\dagger}\right),$
(17a) $\displaystyle H_{\mathrm{sm}}={}$ $\displaystyle
H+\sum_{\beta\gamma=1}^{N}\hbar\Omega_{\beta\gamma}a_{\beta}^{\dagger}a_{\gamma}$
$\displaystyle+\sum_{\alpha\beta=1}^{M,N}\hbar
A_{\alpha}g_{\alpha\beta}\left(a_{\beta}^{\dagger}+a_{\beta}\right),$ (17b)
$\displaystyle\rho={}$
$\displaystyle\mathrm{Tr}_{\mathrm{m}}\\{\rho_{\mathrm{sm}}\\},$ (17c) where
$\rho_{\mathrm{sm}}$ is the density matrix of the system and the discrete
modes.

Therefore, Eq. (missing)17 provides a natural benchmark for any perturbative
equation we propose. Incidentally, this methodology is a generalization of the
one employed in subsection III.4, where Eq. (missing)15 yielded the exact
dynamics.

(a)

(b)

(c)

(d)

Figure 5: Particular case: (a) Scheme of the levels and interaction operators.
(b) Spectral density (red) and Lamb shift integral (blue). The faint dashed
lines represent the 6 different transition energies present in the system. (c)
Time evolution of the populations. (d) Time evolution of the deviation with
respect to the exact dynamics. The dashed lines have been calculated with the
“full” methods including negative and positive eigenvalues, and the dotted
lines with their “positive” counterparts.

For the analysis of the specific system, we set the number of levels,
interaction terms and discrete modes to $L=3$, $M=2$ and $N=2$; with any
simpler choice, the geometric mean ($g\Gamma$) yields a positive semidefinite
Kossakowski matrix. In Fig. 5a and Fig. 5b, we schematically show the energy
structure of the system and the spectral functions of the environment,
produced with the method described above with random parameters. The vertical
dashed lines in Fig. 5b represent all the possible transition energies in the
system. We then solve the dynamics of the particular system with
$|\psi(0)\rangle=(|1\rangle+\sqrt{2}|2\rangle)/3$ as the initial state. To
that end, we employ seven methods, namely, Exact, BRE, BRE(+),
$a\Lambda-a\Gamma$, $a\Lambda-a\Gamma$(+), $a\Lambda-g\Gamma$ and
$a\Lambda-g\Gamma$(+). The Exact and BRE methods solve Eq. (missing)17 and Eq.
(missing)2, respectively. As for the $a\Lambda-a\Gamma$ and $a\Lambda-g\Gamma$
methods, we take our own advice from subsection III.4 and represent the energy
shift always following the arithmetic mean prescription:
$\tilde{\Lambda}_{ij}=(\Lambda_{ij}(\omega_{j})+\Lambda_{ij}(\omega_{i}))/2$.
The various approaches differ in the procedure to obtain the Kossakowski
matrix, i.e., the dissipative part. For the $a\Lambda-a\Gamma$ we take the
arithmetic mean,
$\tilde{\Gamma}_{ij}=(\Gamma_{ij}(\omega_{j})+\Gamma_{ij}(\omega_{i}))/2$, and
for the $a\Lambda-g\Gamma$ the geometric mean,
$\tilde{\Gamma}_{ij}=\sqrt{\Gamma_{ij}(\omega_{j})}\sqrt{\Gamma_{ij}(\omega_{i})}$.
Last, all 3 approximate methods have a positive, (+), counterpart, in which
the Kossakowski matrix is diagonalized, and the negative eigenvalues are
discarded. As a consequence, the (+) methods are all genuine LMEs, and the
question is which of these, if any, accurately represents the dynamics of the
system.

In Fig. 5c and Fig. 5d, we display the evolution of the populations as
obtained from the methods discussed above and the deviation with respect to
the exact result, respectively. The deviation is obtained through the
Frobenius norm of the corresponding density matrices:
$||\rho_{\mathrm{Exact}}(t)-\rho_{\mathrm{Method}}(t)||_{F}$. First, note that
every “full” method including both positive and negative eigenvalues of the
Kossakowski matrix agrees with the exact result, as expected in a sufficiently
weak coupling regime. Accordingly, the dashed lines lie all on top of each
other. Regarding the positive versions of the methods, there are two sets of
lines distinctly separate from the rest, concretely, the ones corresponding to
the BRE(+) and $a\Lambda-a\Gamma$(+). It can be therefore concluded that
removing the negative eigenvalues of the Kossakowski matrix yields incorrect
predictions for the BRE and the arithmetic mean. However, rather remarkably,
there is barely any difference between the $a\Lambda-g\Gamma$ and
$a\Lambda-g\Gamma$(+) methods. Thus, at least in this particular case, it
seems that the geometric mean greatly reduces the effect of the negative
eigenvalues of the Kossakowski matrix. To be more explicit, we diagonalize the
Kossakowski matrices corresponding to the BRE, $a\Lambda-a\Gamma$ and
$a\Lambda-g\Gamma$ methods and present the most significant eigenvalues in
ascending order (first two and last two) in Table 1. There, the blue and red
columns contain the 2 largest negative and positive eigenvalues, respectively.
We observe a striking difference in the ratio between the largest negative and
positive eigenvalues, which is $\mathcal{O}(-1)$ for both the BRE and
$a\Lambda-a\Gamma$ methods, but only $\mathcal{O}(-10^{-2})$ for the
$a\Lambda-g\Gamma$ approach. The smallness of the negative eigenvalues after
performing the geometric mean procedure in the dissipative part explains why
their removal does not imply a loss of accuracy from $a\Lambda-g\Gamma$ to
$a\Lambda-g\Gamma$(+).

Table 1: Negative (blue) and positive (red) eigenvalues of the corresponding Kossakowski matrices. We show only the 4 largest in absolute value, out of the 9 total eigenvalues for $L=3$, rounded to the first significant digit. Method |  Eigenvalues ($\mathrm{meV}/\hbar$) |  |  |  |
---|---|---|---|---|---
| 1st |  2nd | … | 8th | 9th
BRE | $-5$ |  $-9\cdot 10^{-2}$ |  | $2\cdot 10^{-1}$ | $9$
$a\Lambda-a\Gamma$ | $-2$ |  $-5\cdot 10^{-2}$ |  | $9\cdot 10^{-3}$ | $6$
$a\Lambda-g\Gamma$ | $-5\cdot 10^{-2}$ |  $-5\cdot 10^{-5}$ |  | $1\cdot 10^{-2}$ | $4$

###  IV.2 Statistical study

(a)

(b)

(c)

Figure 6: Distribution of the most significant eigenvalues of the Kossakowski
matrix, obtained with the (a) BRE, (b) $a\Lambda-a\Gamma$ and (c)
$a\Lambda-g\Gamma$ methods. Blue and red histograms represent negative and
positive eigenvalues, respectively. Each eigenvalue is specified by the shade
of the color. The horizontal axes are split along the middle of each plot
(vertical dotted line), and display a negative or positive logarithmic scale.

In this subsection, we check whether the reduction of the negative eigenvalues
after the geometric mean procedure in the Kossakowski matrix generalizes in a
random sample of systems. To that end, we again exploit Eq. (missing)16a and
Eq. (missing)16b to generate a sample of 1000 systems and baths with random
parameters. More details can be found in Appendix A. Then, to assess the
consistency of the size reduction of the negative eigenvalues following the
geometric mean procedure, we diagonalize the resulting Kossakowski matrices
and plot the eigenvalue distributions in Fig. 6, where blue and red histograms
indicate negative and positive eigenvalues, respectively. To homogenize the
systems, we normalize the eigenvalues to the largest positive one (9th) in
these plots, hence the dark red peak at $10^{0}$. Notice that the dark blue
histograms close to $-10^{0}$ in Fig. 6a and Fig. 6b indicate the presence of
large negative eigenvalues for the BRE and $a\Lambda-a\Gamma$ methods.
However, rather remarkably, the histogram corresponding to the first
eigenvalue in Fig. 6c does not cluster close to $-10^{0}$, but around
$10^{-3}$. This is in agreement with the previous observation that the
geometric mean significantly reduces the negative eigenvalues. Therefore, Fig.
6 supports the conclusion that removing the negative eigenvalues after the
geometric mean procedure, and thus reaching an LME, should not impact the
accuracy of the QME.

In order to gain information relative to the breakdown of the equations due to
the onset of the strong coupling regime, we introduce a numerical factor in
the spectral functions to artificially scale up the interaction for each
random system instance, spanning a range of several orders of magnitude. We
then solve the dynamics with each method and evaluate the deviation from the
exact evolution in the whole sample, with
$|\psi(0)\rangle=(|1\rangle+\sqrt{2}|2\rangle)/3$ as the initial state. The
deviation is measured with the time-averaged Frobenius norm of the difference
between the density matrices:

$\Delta_{\mathrm{Method}}=\frac{1}{T}\int_{0}^{T}\mathrm{d}t\
||\rho_{\mathrm{Exact}}(t)-\rho_{\mathrm{Method}}(t)||_{F}.$ (18)

In practice, however, it is more sensible to study the statistical properties
of $\log_{10}\Delta_{\mathrm{Method}}$, to account for the skewness of
$\Delta_{\mathrm{Method}}$, which is due to $\Delta_{\mathrm{Method}}$ begin
small but always positive. The results are condensed in Fig. 7. The first main
piece of information presented there is that the “full” methods (continuous
lines) all converge to the exact result when the coupling strength becomes
smaller, as expected. We show with the dotted lines that removing the negative
eigenvalues of the Kossakowski matrix from the BRE and $a\Lambda-a\Gamma$
approaches results in the stabilization of the average deviation around 0.1,
rather than convergence to the exact solution. Hence, we observe that the
importance of the negative eigenvalues is not really dependent on the coupling
strength in those two methods. Strikingly, notice that the dotted line
corresponding to the $a\Lambda-g\Gamma$(+) method follows closely the
deviation of the “full” methods. Therefore, the effect of the negative
eigenvalues is small enough in the $a\Lambda-g\Gamma$ method that they can
safely be neglected, a procedure that automatically yields a genuine LME. The
lines in Fig. 7 imply that we have found an approximate LME that works
accurately in the same coupling regime as the BRE, with two big advantages.
First, our $a\Lambda-g\Gamma$(+) master equation is always free from
unphysical effects, and, second, it admits a quantum trajectory interpretation
of the dynamics with non-negative jump probabilities. Incidentally, we also
remark here that the removal of the negative eigenvalues is known to yield the
closest symmetric positive semidefinite matrix [29].

(a)

(b)

(c)

(d)

(e)

(f)

Figure 7: Time-averaged deviation of the dynamics calculated with Eq.
(missing)18, as a function of the coupling strength factor. The continuous
lines refer to the “full” methods, keeping the negative eigenvalues, and the
dotted lines represent the “positive” methods. The horizontal dotted line
indicates the point beyond which any method undoubtedly fails to reproduce the
exact dynamics.

We end the discussion with one last comment regarding the statistical
distribution of the deviations. With our large sample size, we can calculate
the standard deviation of the distribution of
$\log_{10}\Delta_{\mathrm{Method}}$. Doing so yields a standard deviation of
approximately a quarter of an order of magnitude quite consistently across
different methods and coupling strengths. This means that the lines in Fig. 7
have an approximate width of half an order of magnitude. We have omitted the
dispersion measurement in Fig. 7, as it did not provide much insight and only
made the figure harder to interpret.

##  V Conclusion

Our work shows that the methodology developed in [22] and extended in [25] to
derive an LME was not free of limitations, which prevented the final equation
from being a desired LME in certain regimes. In particular, the Lamb shift
Hamiltonian was susceptible to non-Hermiticity, leading to the loss of the
Hermiticity of the density matrix and, accordingly, to complex populations.
Additionally, the Kossakowski matrix associated with the derived master
equation exhibited non-positive semidefiniteness, violating one of the
assumptions of Lindblad’s theorem, except under relatively simple
configurations. We solve the Hermiticity concern by replacing the geometric
mean introduced in subsection II.2 with an arithmetic mean in the energy shift
terms. Among various alternatives, we highlight this approach due to its
superior simplicity in implementation. Regarding the non-positive definiteness
of the Kossakowski matrix, we observe that while the geometric mean
replacement did not eliminate the negative eigenvalues, these become
substantially smaller than those in the original BRE, and their impact on the
dynamics is markedly suppressed. Therefore, we adopt a pragmatic perspective
and propose the removal of the negative eigenvalues to satisfy the assumptions
of Lindblad’s theorem, after the geometric mean replacement. This option is
supported by the statistical study performed in this work, and particular
simulations of larger systems not shown here for brevity. The combined
implementation of the solutions to both problems yields an accurate and
general LME that, akin to the BRE, maintains close ties with the physical
environment, while simultaneously overcoming the mathematical deficiencies
inherent to the BRE and the physical inadequacy of the secular approximation.

The codes used to generate the findings of this study are available at
[github page].

## Acknowledgments

This work has been funded by the Spanish Ministry of Science, Innovation and
Universities-Agencia Estatal de Investigación through the FPI contract No.
PRE2019-090589 as well as grants PID2021-125894NB-I00, EUR2023-143478, and
CEX2018-000805-M (through the María de Maeztu program for Units of Excellence
in R&D). We also acknowledge financial support from the Proyecto Sinérgico CAM
2020 Y2020/TCS-6545 (NanoQuCo-CM) of the Community of Madrid.

##  Appendix A Random system generation details

The system and bath parameters were stochastically generated using the Python
library NumPy (version 1.26.2). We employ an initial seed equal to 1234321 for
the random number generator for reproducibility. To align the energy scales of
the system and discrete modes, the system energies $\hbar\omega_{\alpha}$ were
confined to the interval $[0.1~{}\mathrm{eV},5.0~{}\mathrm{eV})$, and the mode
energies $\hbar\Omega_{\beta\beta}$ to
$[0.3~{}\mathrm{eV},2.0~{}\mathrm{eV})$. The mode couplings
$\hbar\Omega_{\beta\gamma}$ are symmetric and ranged from $0.0$ eV to $1.0$
eV, and the rates $\hbar\kappa_{\beta}$ were sampled from
$[0.2~{}\mathrm{eV},0.5~{}\mathrm{eV})$. The transition operators were
constructed as Hermitian matrices with zeros on the diagonal and values
ranging from 0 to 1 elsewhere. Their physical dimensions are absorbed in the
spectral functions, making the transition operators unitless. Concerning the
couplings, $\hbar g_{\alpha\beta}$ were uniformly distributed within the
interval $[0.0~{}\mathrm{meV},1.0~{}\mathrm{meV})$, scaled by the square root
of the strength factor discussed in subsection IV.2. This parametrization
ensures a diverse set of systems for our analysis.

## References

  * Schlosshauer [2007] M. Schlosshauer,  _Decoherence: And the Quantum-To-Classical Transition_, The Frontiers Collection (Springer Berlin Heidelberg, 2007).
  * Breuer and Petruccione [2002] H.-P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford University Press, 2002).
  * Harbola  _et al._ [2006] U. Harbola, M. Esposito, and S. Mukamel, Quantum master equation for electron transport through quantum dots and single molecules, Phys. Rev. B 74, 235309 (2006).
  * Härtle  _et al._ [2013] R. Härtle, G. Cohen, D. R. Reichman, and A. J. Millis, Decoherence and lead-induced interdot coupling in nonequilibrium electron transport through interacting quantum dots: A hierarchical quantum master equation approach, Phys. Rev. B 88, 235426 (2013).
  * Carmichael [1993] H. Carmichael,  _An Open Systems Approach to Quantum Optics: Lectures Presented at the Université Libre de Bruxelles, October 28 to November 4, 1991_ (Springer Berlin Heidelberg, 1993).
  * Lambropoulos  _et al._ [2000] P. Lambropoulos, G. M. Nikolopoulos, T. R. Nielsen, and S. Bay, Fundamental quantum optics in structured reservoirs, Reports on Progress in Physics 63, 455 (2000).
  * Diósi [1990] L. Diósi, Landau’s density matrix in quantum electrodynamics, Foundations of Physics 20, 63 (1990).
  * Colas  _et al._ [2022] T. Colas, J. Grain, and V. Vennin, Benchmarking the cosmological master equations, European Physical Journal C 82, 1085 (2022).
  * Lindblad [1976] G. Lindblad, On the generators of quantum dynamical semigroups, Commun. Math. Phys. 48, 119 (1976).
  * Manzano [2020] D. Manzano, A short introduction to the Lindblad master equation, AIP Adv. 10, 025106 (2020).
  * Carmichael and Walls [1973] H. J. Carmichael and D. F. Walls, Master equation for strongly interacting systems, Journal of Physics A: Mathematical, Nuclear and General 6, 1552 (1973).
  * Beaudoin  _et al._ [2011] F. Beaudoin, J. M. Gambetta, and A. Blais, Dissipation and ultrastrong coupling in circuit qed, Phys. Rev. A 84, 043832 (2011).
  * del Pino  _et al._ [2015] J. del Pino, J. Feist, and F. J. Garcia-Vidal, Quantum theory of collective strong coupling of molecular vibrations with a microcavity mode, New Journal of Physics 17, 053040 (2015).
  * Stassi  _et al._ [2016] R. Stassi, S. Savasta, L. Garziano, B. Spagnolo, and F. Nori, Output field-quadrature measurements and squeezing in ultrastrong cavity-QED, New Journal of Physics 18, 123005 (2016).
  * Settineri  _et al._ [2018] A. Settineri, V. Macrí, A. Ridolfo, O. Di Stefano, A. F. Kockum, F. Nori, and S. Savasta, Dissipation and thermal noise in hybrid quantum systems in the ultrastrong-coupling regime, Phys. Rev. A 98, 053834 (2018).
  * Lednev  _et al._ [2023] M. Lednev, F. J. García-Vidal, and J. Feist, A Lindblad master equation capable of describing hybrid quantum systems in the ultra-strong coupling regime (2023), arXiv:2305.13171 [quant-ph] .
  * Jeske  _et al._ [2015] J. Jeske, D. J. Ing, M. B. Plenio, S. F. Huelga, and J. H. Cole, Bloch-Redfield equations for modeling light-harvesting complexes, J. Chem. Phys. 142, 064104 (2015).
  * Eastham  _et al._ [2016] P. R. Eastham, P. Kirton, H. M. Cammack, B. W. Lovett, and J. Keeling, Bath-induced coherence and the secular approximation, Phys. Rev. A 94, 012110 (2016).
  * Hartmann and Strunz [2020] R. Hartmann and W. T. Strunz, Accuracy assessment of perturbative master equations: Embracing nonpositivity, Phys. Rev. A 101, 012103 (2020).
  * Farina and Giovannetti [2019] D. Farina and V. Giovannetti, Open-quantum-system dynamics: Recovering positivity of the Redfield equation via the partial secular approximation, Phys. Rev. A 100, 012107 (2019).
  * Davidović [2020] D. Davidović, Completely positive, simple, and possibly highly accurate approximation of the Redfield equation, Quantum 4, 326 (2020).
  * McCauley  _et al._ [2020] G. McCauley, B. Cruikshank, D. I. Bondar, and K. Jacobs, Accurate Lindblad-form master equation for weakly damped quantum systems across all regimes, npj Quantum Information 6, 74 (2020).
  * Nathan and Rudner [2020] F. Nathan and M. S. Rudner, Universal Lindblad equation for open quantum systems, Phys. Rev. B 102, 115109 (2020).
  * D’Abbruzzo  _et al._ [2023] A. D’Abbruzzo, V. Cavina, and V. Giovannetti, A time-dependent regularization of the Redfield equation, SciPost Phys. 15, 117 (2023).
  * Fernández de la Pradilla  _et al._ [2023] D. Fernández de la Pradilla, E. Moreno, and J. Feist, Vacuum-field-induced state mixing, SciPost Phys. 15, 252 (2023).
  * Johansson  _et al._ [2013] J. Johansson, P. Nation, and F. Nori, QuTiP 2: A Python framework for the dynamics of open quantum systems, Computer Physics Communications 184, 1234 (2013).
  * Imamoğlu [1994] A. Imamoğlu, Stochastic wave-function approach to non-Markovian systems, Phys. Rev. A 50, 3650 (1994).
  * Sánchez-Barquilla  _et al._ [2022] M. Sánchez-Barquilla, F. J. García-Vidal, A. I. Fernández-Domínguez, and J. Feist, Few-mode field quantization for multiple emitters, Nanophotonics 11, 4363 (2022).
  * Higham [1988] N. J. Higham, Computing a nearest symmetric positive semidefinite matrix, Linear Algebra and its Applications 103, 103 (1988).