∎ e1e-mail<EMAIL_ADDRESS>11institutetext: Department of Physics,
Jahangirnagar University, Savar, Dhaka-1342, Bangladesh

# Self-gravitating envelope solitons in a degenerate quantum plasma system

N. Ahmede1  N.A. Chowdhury  A. Mannan  A.A. Mamun

(Received: date / Accepted: date)

###### Abstract

The existence and the basic features of ion-acoustic (IA) envelope solitons in
a self-gravitating degenerate quantum plasma system (SG-DQPS), containing
inertial non-relativistically degenerate light and heavy ion species as well
as inertialess non-relativistically degenerate positron and electron species,
have been theoretically investigated by deriving the nonlinear Schrödinger
(NLS) equation. The NLS equation, which governs the dynamics of the IA waves,
has disclosed the modulationally stable and unstable regions for the IA waves.
The unstable region allows to generate bright envelope solitons which are
modulationaly stable. It is found that the stability and the growth rate
dependent on the plasma parameters (like, mass and number density of the
plasma species). The implications of our results in astronomical compact
object (viz. white dwarfs, neutron stars, and black holes, etc.) are briefly
discussed.

###### Keywords:

Envelope solitons Modulatioanal Insatability Reductive Perturbation Method

††journal: Eur. Phys. J. D

##  1 Introduction

The field of self-gravitating degenerate quantum plasma (DQP) physics is one
of the current interesting research field among the plasma physicists because
of the painstaking observational evidence confirms the existence of such kind
of extreme plasma in astronomical compact objects (viz. white dwarfs, neutron
stars, and black holes, etc. Chandrasekhar1931 ; Koester1990 ; Fowler1994 ;
Koester2002 ; Zaman2017 ) and potential applications in modern technology
(viz. metallic and semiconductor nano-structures, quantum x-ray free-electron
lasers, nano-plasmonic devices Atwater2007 ; Stockmann2011 , metallic nano-
particles, spintronics Wolf2001 , thin metal films, nano-tubes, quantum dots,
and quantum well Manfredi2007 , etc.). The number density of the plasma
species is extremely high in self-gravitating DQP system (SG-DQPS) (order of
$10^{30}~{}cm^{-3}$ in white dwarfs Koester1990 ; Shapiro1983  and order of
$10^{36}~{}cm^{-3}$ or even more in neutron stars Koester1990 ; Shapiro1983 )
which leads to generate a strong gravitational field inside the plasma medium.
Basically, the SG-DQPS contains degenerate inertial light (viz.
${\rm~{}^{1}_{1}H}$ Fletcher2006 ; Killian2006  or ${\rm~{}^{4}_{2}He}$
Chandrasekhar1931 ; Fowler1994  or ${\rm~{}^{12}_{~{}6}C}$ Koester1990 ;
Koester2002 ) and heavy (viz. ${\rm~{}^{56}_{26}Fe}$ Vanderburg2015  or
${\rm~{}^{85}_{37}Rb}$ Witze2014  or ${\rm~{}^{96}_{42}Mo}$ Witze2014 ) ion
species and inertialess degenerate electron and positron species. Heisenberg’s
uncertainty principle established the relationship between the uncertainty to
determine the position and momentum of a particle simultaneously, and
mathematically it can be expressed as, $\Delta x\Delta p\geq\hbar/2$ (where
$\Delta x$ is the uncertainty in position of the particle and $\Delta p$ is
the uncertainty in momentum of the same particle, and $\hbar$ is the reduced
Planck constant). This indicates that the position of the plasma species are
very certain (because of highly dense and compressed plasma species) inside
the plasma system but the momenta of the plasma species are extremely
uncertain. Therefore, these plasma species are uncertain (certain) in momentum
(position) give rise to a very high pressure known as “degenerate pressure”.
The expression for the degenerate pressure $P_{j}$ (degenerate plasma particle
species $j$) as a function of number density ($N_{j}$) is given by
Chandrasekhar1931 ; Fowler1994 ; Mamun2011

$\displaystyle
P_{j}=K_{j}N_{j}^{\gamma},~{}~{}~{}~{}\gamma=\frac{5}{3},~{}~{}~{}K_{j}\simeq\frac{3}{5}\frac{\pi\hbar^{2}}{m_{j}},$
(1)

where $j=e$ ($p$) for the electron (positron) species, and $j=l$ ($h$) for the
light (heavy) ion species, respectively. The $\gamma$ is the relativistic
factor ($\gamma=5/3$ stands for non-relativistic case and $\gamma=4/3$ stands
for ultra-relativistic case) and $m_{j}$ is the mass. It is clear from
equation (1) that the degenerate pressure $P_{j}$ is independent on thermal
temperature but depends on degenerate particle number density $N_{j}$ and mass
$m_{j}$. Finally, the strong gravitational field (degenerate pressure) of the
SG-DQPS wants to squeeze (stretch) the plasma system but they are counter-
balanced to each other.

During the last few years, a large number of authors have studied the
propagation of nonlinear waves in DQP by considering self-gravitational or
without self-gravitational field. Asaduzzaman et al. Asaduzzaman2017  have
investigated the linear and nonlinear propagation of self-gravitational
perturbation mode in a SG-DQPS and found that self-gravitational perturbation
mode becomes unstable when the wavelength of the perturbation mode is minimum.
Mamun Mamun2017  examined the self-gravito shock structures in a SG-DQPS.
Chowdhury et al. Chowdhury2018  have studied the modulational instability (MI)
of nucleus-acoustic waves in a DQP system and found that the bright and dark
envelope solitons are modulationally stable. But to the best of our knowledge,
no attempt has been made to study MI of the ion-acoustic waves (IAWs) by
deriving a nonlinear Schrödinger equation (NLS) and formation of the envelope
solitons in any kind of SG-DQPS. Therefore, in the present work, a SG-DQPS
(containing inertialess degenerate electron and positron species, inertial
degenerate light as well as heavy ion species) has been considered to obtain
the conditions of MI of the IAWs and the formation of the envelope solitons,
and also to identify their basic features.

The manuscript is organized as follows. The basic governing equations for the
dynamics of the SG-DQPS are descried in Sec. 2. The derivation of the NLS
equation is provided in Sec. 3. The stability of the IAWs and envelope
solitons are examined in Sec. 4. A brief discussion is finally presented in
Sec. 5.

##  2 Governing equation

We consider a SG-DQPS containing inertialess degenerate electrons (mass
$m_{e}$; number density $N_{e}$), positrons (mass $m_{p}$; number density
$N_{p}$), inertial degenerate light ions (mass $m_{l}$; number density
$N_{l}$), and heavy ions (mass $m_{h}$; number density $N_{h}$). The detail
information about the light and heavy nuclei is provided in Table 1. The
nonlinear dynamics of the SG-DQPS is described by

$\displaystyle\frac{\partial P_{e}}{\partial
X}=-m_{e}N_{e}\frac{\partial\tilde{\phi}}{\partial X},$ (2)
$\displaystyle\frac{\partial P_{p}}{\partial
X}=m_{p}N_{p}\frac{\partial\tilde{\phi}}{\partial X},$ (3)
$\displaystyle\frac{\partial N_{l}}{\partial T}+\frac{\partial}{\partial
X}(N_{l}U_{l})=0,$ (4) $\displaystyle\frac{\partial U_{l}}{\partial
T}+U_{l}\frac{\partial U_{l}}{\partial
X}=-\frac{\partial\tilde{\phi}}{\partial X}-\frac{1}{m_{l}N_{l}}\frac{\partial
P_{l}}{\partial X},$ (5)
$\displaystyle\frac{\partial^{2}\tilde{\phi}}{\partial X^{2}}=4\pi
G(m_{l}N_{l}+m_{h}N_{h}+m_{e}N_{e}+m_{p}N_{p}),\ $ (6)

where $P_{e}$, $P_{p}$, and $P_{l}$ are the degenerate pressure of the
degenerate electrons, positrons, and light ions, respectively; $X$ ($T$) is
the space (time) variable; $U_{l}$ is the light ion fluid speed;
$\tilde{\phi}$ is the self-gravitational potential; $G$ is the universal
gravitational constant. Now, the charge neutrality condition for the
electrostatic wave potential is

$\displaystyle N_{e}=N_{p}+Z_{l}N_{l}+Z_{h}N_{h},\ $ (7)

where $Z_{l}$ and $Z_{h}$ are the charge state of light and heavy ions,
respectively. Here, it may be noted that the effect of the electrostatic wave
potential has been neglected. Now, we consider normalized variables, namely,
$x=X/L_{q}$, $t=T\omega_{jl}$, $n_{l}=N_{l}/n_{l0}$, $u_{l}=U_{l}/C_{q}$,
$C_{q}=\sqrt{\pi}\hbar n_{e0}^{1/3}/m_{l}$, $\phi=\tilde{\phi}/C_{q}^{2}$,
$\omega_{jl}^{-1}=(4\pi Gm_{l}n_{l0})^{-1/2}$ (where $n_{l0}$ and $n_{e0}$ are
the equilibrium number densities of the light ion and electron species,
respectively). After normalization, Eqs. ($2$)-($6$) can be taken in the
following form

$\displaystyle\frac{\partial\phi}{\partial
x}=-\frac{3}{2}\alpha^{2}\frac{\partial n_{e}^{2/3}}{\partial x},$ (8)
$\displaystyle\frac{\partial\phi}{\partial
x}=\frac{3}{2}\sigma_{1}^{2}\sigma_{2}^{\frac{2}{3}}\frac{\partial
n_{p}^{2/3}}{\partial x},$ (9) $\displaystyle\frac{\partial n_{l}}{\partial
t}+\frac{\partial}{\partial x}(n_{l}u_{l})=0,$ (10)
$\displaystyle\frac{\partial u_{l}}{\partial t}+u_{l}\frac{\partial
u_{l}}{\partial x}=-\frac{\partial\phi}{\partial x}-\beta\frac{\partial
n_{l}^{2/3}}{\partial x},$ (11) $\displaystyle\frac{\partial^{2}\phi}{\partial
x^{2}}=\gamma_{e}(n_{e}-1)-\gamma_{l}(n_{l}-1)+\gamma_{p}(n_{p}-1),\ $ (12)

where $\alpha=m_{l}/m_{e}$, $\sigma_{1}=m_{l}/m_{p}$,
$\sigma_{2}=n_{p0}/n_{e0}$, $\mu=n_{e0}/n_{l0}$, $\beta=(3/2)\mu^{-2/3}$,
$\lambda=n_{p0}/n_{l0}$, $\gamma=Z_{l}m_{h}/Z_{h}m_{l}$ (which is greater than
$1$ for any set of heavy and light ion species),
$\gamma_{e}=\mu(1/\alpha+\gamma/Z_{l})$ (here, $1/\alpha\ll\gamma/Z_{l}$,
where $1/\alpha$ varies from $\sim 10^{-4}$ to $\sim 10^{-3}$, and
$\gamma/Z_{l}$ varies from $\sim 0.1$ to $2.0$, and this means that
$\gamma_{e}\simeq\mu\gamma/Z_{l}$), $\gamma_{l}=\gamma-1$,
$\gamma_{p}=\lambda(1/\sigma_{1}-\gamma/Z_{l})$. For inertialess degenerate
electron and positron, the number densities can be expressed as

$\displaystyle n_{e}=\left(1-\frac{2\phi}{3\alpha^{2}}\right)^{\frac{3}{2}},$
(13) $\displaystyle
n_{p}=\left(1+\frac{2\phi}{3\sigma_{1}^{2}\sigma_{2}^{\frac{2}{3}}}\right)^{\frac{3}{2}}.$
(14)

Now, we substitute Eqs. (13) and (14) into Eq. (12) and extend the resulting
equation up to third order, we get

$\displaystyle\frac{\partial^{2}\phi}{\partial
x^{2}}-\gamma_{l}+\gamma_{l}n_{l}=\gamma_{1}\phi+\gamma_{2}\phi^{2}+\gamma_{3}\phi^{3}+\ldots,\
$ (15)

where

$\displaystyle\gamma_{1}=\left(\frac{\gamma_{p}}{\sigma_{1}^{2}\sigma_{2}^{2/3}}-\frac{\gamma_{e}}{\alpha^{2}}\right),$
$\displaystyle\gamma_{2}=\left(\frac{\gamma_{e}}{6\alpha^{4}}+\frac{\gamma_{p}}{6\sigma_{1}^{4}\sigma_{2}^{4/3}}\right),$
$\displaystyle\gamma_{3}=\left(\frac{\gamma_{e}}{54\alpha^{6}}-\frac{\gamma_{p}}{54\sigma_{1}^{6}\sigma_{2}^{2}}\right).\
$

Table 1: The values of $\gamma$ when ${\rm~{}^{1}_{1}H}$ Killian2006 ; Fletcher2006 , ${\rm~{}^{4}_{2}He}$ Chandrasekhar1931 , and ${\rm~{}^{12}_{~{}6}C}$ Koester1990 ; Koester2002  are considered as the light ion species, and ${\rm~{}^{56}_{26}Fe}$ Vanderburg2015 , ${\rm~{}^{85}_{37}Rb}$ Witze2014 , and ${\rm~{}^{96}_{42}Mo}$ Witze2014  are considered as the heavy ion species. Light ion species |  Heavy ion species |  ${\bf\gamma}$
---|---|---
|  ${\rm~{}^{56}_{26}Fe}$ Vanderburg2015  |  2.16
$~{}~{}~{}~{}{\rm~{}^{1}_{1}H}$ Killian2006 ; Fletcher2006  |  ${\rm~{}^{85}_{37}Rb}$ Witze2014  |  2.30
|  ${\rm~{}^{96}_{42}Mo}$ Witze2014  |  2.28
|  ${\rm~{}^{56}_{26}Fe}$ Vanderburg2015  |  1.08
$~{}~{}~{}~{}{\rm~{}^{4}_{2}He}$ Chandrasekhar1931  |  ${\rm~{}^{85}_{37}Rb}$ Witze2014  |  1.15
|  ${\rm~{}^{96}_{42}Mo}$ Witze2014  |  1.14
|  ${\rm~{}^{56}_{26}Fe}$ Vanderburg2015  |  1.08
$~{}~{}~{}~{}{\rm~{}^{12}_{~{}6}C}$ Koester1990 ; Koester2002  |  ${\rm~{}^{85}_{37}Rb}$ Witze2014  |  1.15
|  ${\rm~{}^{96}_{42}Mo}$ Witze2014  |  1.14

We note that the terms on the right hand side of Eq. (15) are the contribution
of electron and positron species. Thus, Eqs. (10), (11), and (15) describe the
dynamics of the gravitational envelop solitons in the SG-DQPS under
consideration.

##  3 Derivation of NLS equation

To investigate the MI of the IA waves in SG-DQPS, we will derive the NLS
equation by employing the reductive perturbation method Taniuti1969 ;
Chowdhury2017a . So, we first introduce the stretched co-ordinates for
independent variables $x$ and $t$ in terms of $\xi$ and $\tau$ as follows:

$\left.\begin{array}[]{l}\xi=\epsilon(x-v_{g}t),\\\
\tau=\epsilon^{2}t,\end{array}\right\\}$ (16)

where $v_{g}$ is the envelope group velocity and $\epsilon$ is a small
dimensionless expansion parameter. Then we can expand all dependent physical
variables $n_{l}$, $u_{l}$, and $\phi$ in power series of $\epsilon$ as

$\displaystyle
n_{l}=1+\sum_{m=1}^{\infty}\epsilon^{(m)}\sum_{l^{\prime}=-\infty}^{\infty}n_{ll^{\prime}}^{(m)}(\xi,\tau)~{}\mbox{exp}(il^{\prime}\Upsilon),$
(17) $\displaystyle
u_{l}=\sum_{m=1}^{\infty}\epsilon^{(m)}\sum_{l^{\prime}=-\infty}^{\infty}u_{ll^{\prime}}^{(m)}(\xi,\tau)~{}\mbox{exp}(il^{\prime}\Upsilon),$
(18)
$\displaystyle\phi=\sum_{m=1}^{\infty}\epsilon^{(m)}\sum_{l^{\prime}=-\infty}^{\infty}\phi_{l^{\prime}}^{(m)}(\xi,\tau)~{}\mbox{exp}(il^{\prime}\Upsilon),\
$ (19)

where $\Upsilon=kx-wt$ and $k$ ($\omega$) is the real variable representing
the fundamental carrier wave number (frequency). The derivative operators in
Eqs. (10), (11), and (15) are regarded as

$\displaystyle\frac{\partial}{\partial t}\rightarrow\frac{\partial}{\partial
t}-\epsilon
v_{g}\frac{\partial}{\partial\xi}+\epsilon^{2}\frac{\partial}{\partial\tau},$
(20) $\displaystyle\frac{\partial}{\partial
x}\rightarrow\frac{\partial}{\partial
x}+\epsilon\frac{\partial}{\partial\xi}.\ $ (21)

Now, by substituting Eqs. (17)-(21) into Eqs. (10), (11), and(15) and
collecting the different powers of $\epsilon$. Now, the first order ($m=1$)
reduced equations with $l^{\prime}=1$ can be expressed as

$\displaystyle n_{l1}^{(1)}=\frac{k^{2}}{S}\phi_{1}^{(1)},$ (22)
$\displaystyle u_{l1}^{(1)}=\frac{k\omega}{S}\phi_{1}^{(1)},\ $ (23)

where $S=\omega^{2}-\beta_{1}k^{2}$ and $\beta_{1}=2\beta/3$. The
compatibility of the system leads to the linear dispersion relation as

$\displaystyle\omega^{2}=\frac{\gamma_{l}k^{2}}{\gamma_{1}+k^{2}}+\beta_{1}k^{2}.$
(24)

Figure 1: The variation of $\omega$ with $k$ for different values of $\mu$;
along with $\alpha=3.67\times 10^{3}$, $\gamma=2.16$, $\gamma/Z_{l}=0.5$,
$\sigma_{1}=3.68\times 10^{3}$, $\sigma_{2}=0.3$, and $\lambda=0.2$.

Figure 2: The variation of $P/Q$ with $k$ for different values of $\mu$; along
with $\alpha=3.67\times 10^{3}$, $\gamma=2.16$, $\gamma/Z_{l}=0.5$,
$\sigma_{1}=3.68\times 10^{3}$, $\sigma_{2}=0.3$, and $\lambda=0.2$.

The dispersion characteristics of the wave are depicted in Fig. 2 [obtained
from Eq. (24)], which indicates that (a) the angular wave frequency ($\omega$)
of the IAWs exponentially decreases with the increase of $k$; (b) the value of
$\omega$ increases with the increase of $n_{e0}$ for the fixed value of
$n_{l0}$ (via $\mu=n_{e0}/n_{l0}$). The second order ($m=2$) reduced equations
with $l^{\prime}=1$ are given by,

$\displaystyle
n_{l1}^{(2)}=\frac{k^{2}}{S}\phi_{1}^{(2)}+\frac{2ik\omega(kv_{g}-\omega)}{S^{2}}\frac{\partial\phi_{1}^{(1)}}{\partial\xi},$
(25) $\displaystyle
u_{l1}^{(2)}=\frac{k\omega}{S}\phi_{1}^{(2)}+\frac{i(\omega^{2}+\beta_{1}k^{2})(kv_{g}-\omega)}{S^{2}}\frac{\partial\phi_{1}^{(1)}}{\partial\xi},\
$ (26)

thus, the expression for $v_{g}$ is obtained as

$\displaystyle v_{g}=\frac{\partial\omega}{\partial
k}=\frac{\gamma_{l}\omega^{2}-(\omega^{2}-\beta_{1}k^{2})^{2}}{k\omega\gamma_{l}}.\
$ (27)

The amplitude of the second-order harmonics is found to be proportional to
$|\phi_{1}^{(1)}|^{2}$ and these are expressed as

$\left.\begin{array}[]{l}n_{l2}^{(2)}=C_{1}|\phi_{1}^{(1)}|^{2},\\\
u_{l2}^{(2)}=C_{2}|\phi_{1}^{(1)}|^{2},\\\
\phi_{2}^{(2)}=C_{3}|\phi_{1}^{(1)}|^{2},\\\
n_{l0}^{(2)}=C_{4}|\phi_{1}^{(1)}|^{2},\\\
u_{l0}^{(2)}=C_{5}|\phi_{1}^{(1)}|^{2},\\\
\phi_{0}^{(2)}=C_{6}|\phi_{1}^{(1)}|^{2},\\\ \end{array}\right\\}$ (28)

where the coefficients are

$\displaystyle C_{1}=\frac{2C_{3}k^{2}S^{2}+3\omega^{2}k^{4}}{2S^{3}},$
$\displaystyle C_{2}=\frac{C_{1}\omega S^{2}-\omega k^{4}}{kS^{2}},$
$\displaystyle
C_{3}=\frac{3\gamma_{l}\omega^{2}k^{4}-2\gamma_{2}S^{3}}{2S^{3}(\gamma_{1}+4k^{2})-2\gamma_{l}k^{2}S^{2}},$
$\displaystyle C_{4}=\frac{C_{6}S^{2}+k^{2}\omega^{2}+2\omega
v_{g}k^{3}-\beta_{2}k^{4}}{S^{2}(v_{g}^{2}-\beta_{1})},$ $\displaystyle
C_{5}=\frac{C_{4}v_{g}S^{2}-2\omega
k^{3}}{S^{2}},~{}~{}~{}~{}\beta_{2}=\beta/9,$ $\displaystyle
C_{6}=\frac{(k^{2}\omega^{2}+2\omega
v_{g}k^{3}-\beta_{2}k^{4})\gamma_{l}-2\gamma_{2}S^{2}(v_{g}^{2}-\beta_{1})}{\gamma_{1}S^{2}(v_{g}^{2}-\beta_{1})-\gamma_{l}S^{2}}.\
$

Finally, by substituting all the Eqs. (22)$-$(28) into the third order part
($m=3$) and $l^{\prime}=1$ and simplifying them, we can obtain the following
NLS equation:

$\displaystyle
i\frac{\partial\Phi}{\partial\tau}+P\frac{\partial^{2}\Phi}{\partial\xi^{2}}+Q|\Phi|^{2}\Phi=0,\
$ (29)

where $\Phi=\phi_{1}^{(1)}$ for simplicity. The coefficient of dispersion and
nonlinear terms $P$ & $Q$ are given by

$\displaystyle
P=\frac{F_{1}-4\omega\beta_{1}^{2}k^{4}-3kv_{g}\omega^{4}}{2\gamma_{l}k^{2}\omega^{2}},$
(30) $\displaystyle
Q=\frac{S^{2}[3\gamma_{3}+2\gamma_{2}(C_{3}+C_{6})-F_{2}]}{2\gamma_{l}\omega
k^{2}},\ $ (31)

where
$F_{1}=4\beta_{1}k^{2}\omega^{3}+2\beta_{1}v_{g}\omega^{2}k^{3}+v_{g}\beta_{1}^{2}k^{5}$,
$F_{2}=k^{2}/S^{2}[2\omega
k\gamma_{l}(C_{2}+C_{5})+\gamma_{l}\omega^{2}(C_{1}+C_{4})+(\gamma_{l}\beta_{3}k^{6}/S^{2})]$,
and $\beta_{3}=4\beta/81$.

##  4 Stability analysis and envelope solitons

Let us now analysis the MI of IAWs by considering the linear solution of the
NLS equation (29) in the form
$\Phi=\hat{\Phi}~{}e^{iQ|\hat{\Phi}|^{2}\tau}+c.~{}c$ (c. c denotes the
complex conjugate), where $\hat{\Phi}=\hat{\Phi}_{0}+\epsilon\hat{\Phi}_{1}$
and $\hat{\Phi}_{1}=\hat{\Phi}_{1,0}e^{i(\tilde{k}\xi-\tilde{\omega}\tau)}$\+
c. c. Now, by substituting these values into Eq. (29), one readily obtains the
following nonlinear dispersion relation Chowdhury2018 ; Sultana2011 ;
Schamel2002 ; Kourakis2005 ; Fedele2002 ; Chowdhury2017b

$\displaystyle\tilde{\omega}^{2}=P^{2}\tilde{k}^{2}\left(\tilde{k}^{2}-\frac{2|\hat{\Phi}_{0}|^{2}}{P/Q}\right).\
$ (32)

Here, the perturbed wave number $\tilde{k}$ and the perturbed frequency
$\tilde{\omega}$ are different from the carrier wave number $k$ and frequency
$\omega$. It is observed from Eq. (32) that the IAWs will be modulationally
stable (unstable) in SG-DQPS for that range of values of $\tilde{k}$ in which
$P/Q$ is negative (positive), i.e., $P/Q<0$ ($P/Q>0$). When
$P/Q\rightarrow\pm\infty$, the corresponding value of $k$ ($=k_{c}$) is known
as the critical or threshold wave number ($k_{c}$) for the onset of MI. The
variation of $P/Q$ with $k$ for $\mu$ is shown in Fig. 2 and which clearly
indicates that (a) the IAWs are modulatonally stable (unstable) in SG-DQPS for
small (long) wavelength; (b) the $k_{c}$ increases with the increase of
$n_{e0}$ for constant value of $n_{l0}$ (via $\mu=n_{e0}/n_{l0}$). In the
modulationally unstable ($P/Q>0$) region and under this condition
$\tilde{k}<\tilde{k}_{c}=\sqrt{2|\hat{\Phi}_{0}|^{2}(Q/P)}$, the MI growth
rate can be written [from Eq. (32)] as

$\displaystyle\Gamma=|P|\tilde{k}^{2}\sqrt{\frac{\tilde{k}^{2}_{c}}{\tilde{k}^{2}}-1}.\
$ (33)

Figure 3: The variation of $\Gamma$ with $k$ for different values of
$\sigma_{2}$; along with $\mu=1.5$, $\alpha=3.67\times 10^{3}$, $\gamma=2.16$,
$\gamma/Z_{l}=0.5$, $\sigma_{1}=3.68\times 10^{3}$, $\lambda=0.2$, $k=0.0004$,
and $\phi=0.8$.

Figure 4: The variation of $\Gamma$ with $k$ for different values of $\gamma$;
along with $\mu=1.5$, $\alpha=3.67\times 10^{3}$, $\gamma/Z_{l}=0.5$,
$\sigma_{1}=3.68\times 10^{3}$, $\sigma_{2}=0.3$, $\lambda=0.2$, $k=0.0004$,
and $\phi=0.8$.

The effect of $\sigma_{2}$ and $\gamma$ on the growth rate are presented in
Figs. 4 and 4, where $\Gamma$ is plotted against $\tilde{k}$ and it is
observed that (a) the growth rate ($\Gamma$) increases with the increase in
the value of positron number density $n_{p0}$, but decreases with increase of
the electron number density $n_{e0}$ (via $\sigma_{2}=n_{p0}/n_{e0}$); (b) the
maximum value of $\Gamma$ increases (decreases) with the decrease of $m_{h}$
($m_{l}$) for the fixed value of $Z_{l}$ and $Z_{h}$ (via
$\gamma=Z_{l}m_{h}/Z_{h}m_{l}$); (c) on the other hand, the maximum value of
$\Gamma$ increases (decreases) with the decrease of $Z_{l}$ ($Z_{h}$) for the
fixed value of $m_{h}$ and $m_{l}$ (via $\gamma=Z_{l}m_{h}/Z_{h}m_{l}$). So,
the charge state and mass of the light and heavy ion plays an opposite role to
manifest the $\Gamma$ in SG-DQPS. The physics of this result is that the
nonlinearity of the SG-DQPS increases (decreases) with the increase of the
value of $m_{l}$ or $Z_{h}$ ($m_{h}$ or $Z_{l}$) which enhance (suppress) the
maximum value of the $\Gamma$.

The self-gravitating bright envelop solitons are generated in the
modulationally unstable region (when $P/Q>0$) and the solitonic solution of
Eq. (29) for the self-gravitating bright envelope solitons can be written as
Chowdhury2018 ; Sultana2011 ; Schamel2002 ; Kourakis2005 ; Fedele2002

$\displaystyle\Phi(\xi,\tau)=\left[\psi_{0}~{}\mbox{sech}^{2}\left(\frac{\xi-U\tau}{W}\right)\right]^{1/2}$
$\displaystyle\hskip
28.45274pt\times\exp\left[\frac{i}{2P}\left\\{U\xi+\left(\Omega_{0}-\frac{U^{2}}{2}\right)\tau\right\\}\right],$
(34)

Figure 5: The variation of Re($\Phi$) with $\xi$ for bright envelope solitons;
along with $\mu=1.5$, $\alpha=3.67\times 10^{3}$, $\gamma=2.16$,
$\gamma/Z_{l}=0.5$, $\sigma_{1}=3.68\times 10^{3}$, $\sigma_{2}=0.3$,
$\lambda=0.2$, $U=0.001$, $k=0.0004$, $\psi_{0}=0.005$, $\tau=0$,
$\Omega_{0}=0.04$.

Figure 6: The variation of the $|\Phi|$ with $\xi$ and $\tau$ for bright
envelope solitons; along with $\mu=1.5$, $\alpha=3.67\times 10^{3}$,
$\gamma=2.16$, $\gamma/Z_{l}=0.5$, $\sigma_{1}=3.68\times 10^{3}$,
$\sigma_{2}=0.3$, $\lambda=0.2$, $U=0.001$, $k=0.0004$, $\psi_{0}=0.005$,
$\tau=0$, $\Omega_{0}=0.04$.

where $U$ is the propagation speed of the localized pulse, $W$ is the pulse
width which can be written as $W=\sqrt{2|P/Q|/\psi_{0}}$ ($\psi_{0}$ is the
constant amplitude), and $\Omega_{0}$ is the oscillating frequency for $U=0$.
The self-gravitating bright envelop solitons which are obtained from the
numerical analysis of Eq. (34), are depicted in Figs. 6 and 6. The bright
envelop solitons remain same as time($\tau$) passes, i.e., the self-
gravitating bright envelop solitons are modulationally stable (please see Fig.
6).

##  5 Discussion

In our above analysis, we have considered an unmagnetized realistic laboratory
or astrophysical SG-DQPS consisting of inertialess non-relativistically
degenerate electron and positron species, inertial non-relativistically
degenerate light ion species as well as heavy ion species. The NLS equation
has been derived by employing the well-known reductive perturbation method,
which governs the evolution of nonlinear IAWs. The notable informations that
have been found from our theoretical investigation, can be pin-pointed as
follows:

  1. 1.

The angular wave frequency ($\omega$) of the IAWs exponentially decreases with
the increase of $k$. On the other hand, the value of $\omega$ increases with
the increase of $n_{e0}$ for the fixed value of $n_{l0}$ (via
$\mu=n_{e0}/n_{l0}$).

  2. 2.

The IAWs will be modulationally stable (unstable) for that range of values of
$k$ in which $P/Q$ is negative (positive), i.e., $P/Q<0$ ($P/Q>0$).

  3. 3.

The growth rate ($\Gamma$) increases with the increase in the value of
positron number density $n_{p0}$, but decreases with increase of the electron
number density $n_{e0}$ (via $\sigma_{2}=n_{p0}/n_{e0}$). On the other hand,
the maximum value of $\Gamma$ increases (decreases) with the decrease of
$m_{h}$ ($m_{l}$) for the fixed value of $Z_{l}$ and $Z_{h}$ (via
$\gamma=Z_{l}m_{h}/Z_{h}m_{l}$). Furthermore, the maximum value of $\Gamma$
increases (decreases) with the decrease of $Z_{l}$ ($Z_{h}$) for the fixed
value of $m_{h}$ and $m_{l}$ (via $\gamma=Z_{l}m_{h}/Z_{h}m_{l}$).

  4. 4.

The self-gravitating bright envelop solitons remain same (modulationally
stable) as time passes.

The findings of this theoretical investigation may be useful for understanding
the nonlinear structure (bright envelope solitons) of a SG-DQPS in space (viz.
neutron stars and white dwarf Chandrasekhar1931 ; Koester1990 ; Fowler1994 ;
Koester2002 ; Zaman2017 ).

## References

  * (1) S. Chandrasekhar, Astrophys. J. 74, 81 (1931)
  * (2) D. Koester, G. Chanmugam, Rep. Prog. Phys. 53, 837 (1990)
  * (3) R.H. Fowler, J. Astrophys. Astr. 15, 115 (1994)
  * (4) D. Koester, Astron. Astrophys. Rev. 11, 33 (2002)
  * (5) D.M.S. Zaman, M. Amina, P.R. Dip, A.A. Mamun 2017 Eur. Phys. J. Plus 132, 457.
  * (6) H.A. Atwater, Sci. Am. 296, 56 (2007)
  * (7) M.I. Stockmann, Phys. Today 64, 39 (2011)
  * (8) S.A. Wolf, D. Awschalom, R.A. Buhrman, Science 294, 1488 (2001)
  * (9) G. Manfredi, P.A. Hervieux, Appl. Phys. Lett. 91, 061108 (2007)
  * (10) S.L. Shapiro, S.A. Teukolsky, Black Holes, White Dwarfs and Neutron Stars: the Physics of Compact Objects (John Wiley & Sons, New York, 1983)
  * (11) R.S. Fletcher, X.L. Zhang, S.L. Rolston, Phys. Rev. Lett. 96, 105003 (2006)
  * (12) T.C. Killian, Nature (London) 441, 297 (2006)
  * (13) A. Vanderburg et al., Nature (London) 526, 546 (2015)
  * (14) A. Witze, Nature 510, 196 (2014)
  * (15) A.A. Mamun, P.K. Shukla, Europhys. Lett. 94, 65002 (2011)
  * (16) M. Asaduzzaman, A. Mannan, A.A. Mamun, Phys. Plasmas 24, 052102 (2017)
  * (17) A.A. Mamun, Phys. Plasmas 24, 102306 (2017)
  * (18) N.A. Chowdhury, M.M. Hasan, A. Mannan, A.A. Manun, Vacuum 147, 31 (2018)
  * (19) T. Taniuti, N. Yajima, J. Math. Phys. 10, 1369 (1969)
  * (20) N.A. Chowdhury, A. Mannan, M.M. Hasan, A.A. Mamun, Chaos 27, (2017) 093105.
  * (21) S. Sultana, I. Kourakis, Plasma Phys. Control. Fusion 53, 045003 (2011)
  * (22) R. Fedele, H. Schamel, Eur. Phys. J. B 27, 313 (2002)
  * (23) I. Kourakis, P.K. Shukla, Nonlinear Proc. Geophys. 12, 407 (2005)
  * (24) R. Fedele, H. Schamel, Eur. Phys. J. B 27, 313 (2002)
  * (25) N.A. Chowdhury, A. Mannan, A.A. Mamun, Phys. Plasmas 24, (2017) 113701.


# Large-scale Epidemiological modeling: Scanning for Mosquito-Borne Diseases
Spatio-temporal Patterns in Brazil

Eduardo C. Araujo1, Claudia T. Codeço2, Sandro Loch1,
Luã B. Vacaro1, Laís P. Freitas3,4, Raquel M. Lana5,
Leonardo S. Bastos2, Iasmim F. de Almeida6, Fernanda Valente7,
Luiz M. Carvalho1 and Flávio C. Coelho1
1EMAp - FGV, Brazil, 2PROCC-FIOCRUZ, Brazil, 3ESPUM, Canada,
4CReSP, Canada, 5BSC, Spain, 6ENSP-FIOCRUZ, Brazil, 7Observatório de
Bioeconomia, FGV, Brazil

###### Abstract

The influence of climate on mosquito-borne diseases like dengue and
chikungunya is well-established, but comprehensively tracking long-term
spatial and temporal trends across large areas has been hindered by fragmented
data and limited analysis tools. This study presents an unprecedented
analysis, in terms of breadth, estimating the SIR transmission parameters from
incidence data in all 5,570 municipalities in Brazil over 14 years (2010-2023)
for both dengue and chikungunya. We describe the Episcanner computational
pipeline, developed to estimate these parameters, producing a reusable dataset
describing all dengue and chikungunya epidemics that have taken place in this
period, in Brazil. The analysis reveals new insights into the climate-epidemic
nexus: We identify distinct geographical and temporal patterns of arbovirus
disease incidence across Brazil, highlighting how climatic factors like
temperature and precipitation influence the timing and intensity of dengue and
chikungunya epidemics. The innovative Episcanner tool empowers researchers and
public health officials to explore these patterns in detail, facilitating
targeted interventions and risk assessments. This research offers a new
perspective on the long-term dynamics of climate-driven mosquito-borne
diseases and their geographical specificities linked to the effects of global
temperature fluctuations such as those captured by the ENSO index.

Keywords: episcanner, mosquito-borne diseases, climate-epidemic nexus,
geographical and temporal patterns

##  1 Introduction

Dengue and chikungunya are among the most significant mosquito-borne diseases
in terms of disease burden and potential for global expansion. Worldwide, the
estimated number of dengue cases has increased from 30 million, in 1990, to 56
million in 2019 (Yang et al.,, 2021), reaching more than 100 countries.
Chikungunya, primarily transmitted by the same mosquitoes Aedes aegypti and
Aedes albopictus, has been spreading worldwide since the early 2000s and is
now present in all continents with tropical zones (Manzoor et al.,, 2022).

The temporal dynamics of these vector-borne diseases is characterized by
seasonal and multi-year cycles that can vary from place to place. Within
seasons, they can further differ in their reproductive number, peak size,
timing and total attack rate (de Almeida et al.,, 2022). This spatial-temporal
heterogeneity is strongly associated with the climate and meteorological
conditions that affect the mosquito vectorial capacity and viral transmission,
in combination with the history of previous population exposure that affects
the level of collective immunity. The urban environment, characterized by
crowding, inadequate garbage and water services, can further increase the
population vulnerability and exposure, while effective vector control measures
can potentially flatten the curve. As a result, all of these local factors
contribute to generate very diverse spatially-temporal epidemic dynamics
locally.

Dengue has maintained a continuous presence in the tropical regions of Brazil
since 1986 (Rodriguez-Barraquer et al.,, 2011), but with patterns varying from
transient transmission to epidemic and endemic persistence (de Almeida et
al.,, 2022). Chikungunya fever was introduced in Brazil in 2014 and persisted,
with outbreaks increasing in frequency (de Souza et al.,, 2023). In the last
years, probably due to the gradual warming of the southern states, outbreaks
of both diseases have begun to occur in areas with subtropical and temperate
climates (Codeco et al.,, 2022; de Almeida et al.,, 2023). Understanding how
climate can affect the reproductive number and timing of dengue and
chikungunya epidemics locally is useful for guiding decisions at the national
level.

This study is centred on estimating and examining the characteristics of all
the dengue and chikungunya epidemics that occurred in Brazilian municipalities
from 2010 to 2023. The key descriptors are the reproductive number, size and
timing, estimated by the Episcanner pipeline, a computational pipeline
developed for this end.

##  2 Methods

###  2.1 Data Sources

Brazil comprises 5570 municipalities, divided into 26 states, plus the Federal
District, the country’s capital. The disease data for all the municipalities
was obtained from the Infodengue project (Codeco et al.,, 2018), available at
info.dengue.mat.br, for the period between 2010 and 2023. The Infodengue
project gathers notification data for dengue, chikungunya, and Zika from the
Brazilian Ministry of Health, sanitizes and re-publishes it, along with
epidemiological analyses, such as estimates for the effective reproduction
number, ${\cal R}_{t}$, by week, for all municipalities. The code for
obtaining the Infodengue data is available on the GitHub repository associated
with this paper. The climate data were obtained from the Copernicus ERA5
reanalysis dataset (Muñoz Sabater and Service,, 2019) and daily averaged by
the municipality. The climate variables included land-surface temperature,
relative humidity, precipitation and atmospheric pressure.

###  2.2 The Episcanner Pipeline

If all municipalities experienced a dengue epidemic per year, in 10 years,
there would be 55,700 epidemics, varying in velocity, timing, and magnitude.
This is a large number of models to be fitted individually. Episcanner is a
computational pipeline that efficiently scans all these potential epidemics
and estimates the epidemiological descriptors whenever an epidemic is
detected. The epidemic descriptors are the basic reproduction number (${\cal
R}_{0}$), peak week, the total outbreak size, as well as timing and duration
of the epidemic.

The core of the Episcanner pipeline is the Richards logistic growth function,
a sigmoid curve that has a one-to-one relationship with the SIR
epidemiological model (Wang et al.,, 2012). Based on this relationship, the
epidemiologically meaningful parameters of the SIR (Susceptible-Infectious-
Recovered) transmission model can be efficiently derived from the fitted
Richard model, which serves as a curve template for identifying an epidemic.
This simplification is pivotal for curtailing the computational cost
associated with estimating the epidemic descriptors of interest for all
municipalities and epidemic years considered.

The Episcanner pipeline is depicted in Figure 1, and can be summarized in the
following steps:

  1. 1.

Infodengue’s incidence time series are obtained from their API.

  2. 2.

Filter years with a potential epidemic, defined as an annual time series
having a minimum of 3 weeks with at least 0,9 probability of ${\cal R}_{t}>1$
and more than fifty cumulative cases. Any year with less than that is
discarded as a short outbreak. This is done for all municipalities.

  3. 3.

Split all-time series at week 45 to obtain a collection of 52 week-long time
series, one for each year and each municipality. The cut point defines the
typical beginning of the dengue season in Brazil.

  4. 4.

Richards model is fitted to the observed curves, using an optimization
algorithm, generating the SIR parameter estimates (see below for details).

  5. 5.

Epidemiological parameters dataset is created and made available through the
Mosqlimate API.

Figure 1: The computational pipeline of Episcanner inferential engine is the
dark blue block. Data from Infodengue, after being filtered to identify
candidates for epidemic years, is then fitted to the Richards model. The
resulting estimated parameter dataset is made available through an API as a
JSON object (https://api.mosqlimate.org/datastore/). This data can be
visualized in the Episcanner dashboard, which is part of the Infodengue
website (https://info.dengue.mat.br/epi-scanner).

The complete set of parameters produced from all Brazilian municipalities can
be browsed through the online dashboard (Coelho,, 2024) called Episcanner,
within the Infodengue website. This dashboard is updated weekly and allows
anyone to browse the estimated epidemic characteristics of dengue and
chikungunya, as described in this paper.

###  2.3 Estimation of the Epidemiological Parameters

The Richards model is defined by the following equation:

$\displaystyle\begin{split}J(t)=L-L[1+\alpha
e^{b(t-t_{j})}]^{-1/\alpha},\end{split}$ (1)

where $J(t)$ is interpreted as the accumulated number of cases at week $t$,
$L>0$ is the estimated total number of cases at the end of the epidemic,
$t_{j}\in[1,52]$ is the week of the inflexion point of the sigmoid curve,
which corresponds to the peak of the epidemic. The remaining parameters
$\alpha$ and $b$ are additional parameters to be estimated.

When mapping to the SIR model, the function $J(t)$ corresponds to the number
of recovered individuals at time $t$, which in the SIR model is denoted
$R(t)$. The equivalent SIR model’s parameters can be derived from the Richards
equation’s parameters (1) through the equations given below (refer to Wang et
al., (2012) for derivation):

$\displaystyle\begin{split}\beta=\cfrac{b}{\alpha},\end{split}$ (2)

$\displaystyle\begin{split}\gamma=b\left(\cfrac{1}{\alpha}-1\right),\end{split}$
(3)

$\displaystyle\begin{split}R_{0}=\cfrac{\beta}{\gamma}=\cfrac{1}{1-\alpha}.\end{split}$
(4)

Besides the reproductive number ($R_{0}$), epidemic size ($L$) and peak week
($t_{j}$), the Episcanner dataset also includes the onset of the epidemic
($w_{s}$) that is defined as the week in which the new cases crossed the 5%
percentile, and the final week of the epidemic ($w_{e}$) defined as the week
when the cases dropped below this threshold (See Figure 2). From these points,
the epidemic duration ($w_{e}-w_{s}$) in weeks is computed and also added to
the dataset.

From now on, we will refer to these derived parameters, collectively, as
epidemic parameters or epidemic descriptors.

Figure 2: Example of the outcome of estimating the start and end weeks of an
epidemic for Rio de Janeiro municipality in all epidemic years between 2010
and 2023. The red area represents the cumulative number of cases. Vertical
dashed lines indicate the start and end of the epidemic. The blue curve is the
fitted Richards curve, and the orange one at the bottom is the absolute error
of the fit every week.

###  2.4 Fitting the model to data

For every selected city and epidemic year, as defined in section 2.1, we
fitted the Richards model (eq. 1) to the data by solving the following
optimization problem: determining the values of the parameters set
$\xi:=\\{L,\alpha,b,t_{j}\mid L,\alpha,b,t_{j}\\}\in\mathbb{R}^{+}$ that
minimizes the error function below:

$\displaystyle\begin{split}\operatorname{argmin}_{\xi}\sum_{t=1}^{52}\frac{(C(t)-J(t,\xi))^{2}}{52},\end{split}$
(5)

where $C(t)$ is the observed cumulative number of cases at week $t$. The
following ranges restricted the possible values of the parameters of $\xi$:
$\gamma\in[0.95,1.05]$ (1-week infectious period), $\alpha\in[0.001,1]$,
$b\in[10^{-6},1]$, $t_{j}\in[5,35]$ (based on the overall distribution of peak
weeks). Additionally, we restricted the possible values of $\alpha$ to always
satisfy $\alpha=\frac{b}{\gamma+b}$ (Wang et al.,, 2012).

This fitting process is computationally more efficient than to fit the SIR
ordinary differential equations, since it avoids numerical integration steps.
The optimization was implemented with the _lmfit_ Python
package111https://lmfit.github.io/lmfit-py/intro.html, using a global
optimization approach based on the differential evolution algorithm (Storn and
Price,, 1997). Using the estimated epidemiological parameters, we explored
their association with climate variables across time and space.

###  2.5 Predicting the week of epidemic peak

One question of interest is how early the epidemic peak will occur in a given
year, in a given place. To test if the dengue epidemic descriptors, estimated
by Episcanner, show a correlation to epidemiological, demographic and climate
variables, we built a Histogram Gradient Boosting Regression (HGBR) model. The
outcome of the model is the peak week of epidemics, $W^{peak}_{i,y}$for every
epidemic year in every location. Predictors included lagged climate variables,
case counts and population sizes, as well as lagged epidemic descriptors. The
regression model proposed can be written as below

$\displaystyle\begin{split}W^{peak}_{i,y}\sim
C_{i,y-1}+D_{i,y-1}+E_{i,y-1}+\varepsilon,\end{split}$ (6)

where $W^{peak}_{i,y}$ is the peak week of the epidemic at municipality $i$ on
year $y$. $C_{i,y-1}$, $D_{i,y-1}$ and $E_{i,y-1}$ are respective sets of
climate, demographic, and epidemiological predictors. A comprehensive
description of all model features is provided in Table S1 of the supplementary
material. We used the Scikit-learn Python package to fit the model
(https://scikit-learn.org).

The model was fitted using data encompassing all Brazilian cities for the
entire period, filtered as described above for epidemic parameter estimation.
A single model was fitted to each geographical region of Brazil. Brazil has 5
geographical regions, namely: north (N), northeast (NE), midwest (MW),
southeast (SE) and south (S). This model was fitted only for dengue due to its
broader geographic coverage. We used Shapley Additive Explanations (SHAP
values) framework (Lundberg and Lee,, 2017) to analyze the importance of each
feature (predictor) of the model.

##  3 Results

###  3.1 Spatio-temporal patterns of dengue and chikungunya epidemics

Out of the 5570 Brazilian cities, 4096 (73%) had at least one epidemic year
for dengue and 1263 cities (22%) for chikungunya between 2010 and 2023. In
total, there were 18567 dengue epidemic events and 2199 chikungunya epidemic
events. Table 1 contains summary statistics for some of the parameters
estimated. The full set of parameters also includes the SIR model’s $\beta$
and $\gamma$ parameters, the $\alpha$ parameter from the Richards equation,
and the starting and ending weeks of the epidemic.

Table 1: Summary statistics for selected parameters estimated for dengue and chikungunya across all epidemics. dengue | ${\cal R}_{0}$ | duration | peak week | total cases | incidence (per 100k)
---|---|---|---|---|---
mean | 1.97 | 26.62 | 22.04 | 992.09 | 2146.22
std | 0.52 | 10.38 | 5.80 | 4262.09 | 2841.39
25% | 1.60 | 19.00 | 18.40 | 113.75 | 523.45
50% | 1.86 | 25.00 | 22.11 | 240.40 | 1203.41
75% | 2.20 | 33.00 | 25.24 | 624.70 | 2626.13
chikungunya |  |  |  |  |
mean | 2.12 | 25.47 | 24.69 | 624.43 | 1156.76
std | 0.70 | 12.17 | 7.29 | 2341.62 | 1620.65
25% | 1.60 | 16.0 | 19.78 | 98.65 | 213.80
50% | 1.94 | 24.0 | 24.68 | 195.33 | 594.58
75% | 2.45 | 34.0 | 29.51 | 461.99 | 1436.29

The spatial distribution of epidemic parameters (${\cal R}_{0}$, peak week,
and epidemic duration) exhibit discernible patterns across different regions
of Brazil. Notably, the average duration of dengue epidemics exhibits an
ascending gradient from the southeast (SE) to the northwest (NW) regions, see
Figure 3. The reproductive numbers of both dengue and chikungunya exhibited
spatial and temporal variability across the country (see Figure S1 of the
supplementary material). This variability is evident not only between
different geographic regions but also over time.

Figure 3: Median Epidemic duration (in weeks) estimated for Dengue and
Chikungunya cases. Note the increase in the duration of dengue epidemics as
one moves from the southern to the northern states. This trend is less marked
but also visible for chikungunya epidemics. In the smaller map in the middle,
we have the geographic regions of Brazil: North(N), Northeast (NE), Midwest
(MW), Southeast (SE), and South (S).

The estimated duration of dengue epidemics varied from short ($<10$ weeks) to
year-long. Longer epidemics were seen mostly in the North region (Figure 3).
The peak week tended to occur earlier for dengue than chikungunya and in the
west earlier than the east side of the country (see Figure 4). Looking at the
country-wide time series for dengue, we can see that higher incidence peaks
are associated with shorter epidemics with higher ${\cal R}_{0}$ (Figure 5).

Figure 4: Median peak weeks estimated for Dengue and Chikungunya epidemics
between 2010 and 2023 by municipality. The inset map shows Brazil’s geographic
regions: North(N), Northeast (NE), Midwest (MW), Southeast (SE), and South
(S). Figure 5: The boxplots, from top to bottom, show the estimated peak week,
the estimated basic reproduction number, and the epidemic duration in weeks,
for dengue epidemics by year. Above, we find the multivariate ENSO index v2
series, and the last panel is the time series of dengue cases in Brazil
between November 2010 and November 2023.

A nonlinear relationship was seen between epidemics duration and $R_{0}$ (see
Figure S2). Epidemics with high $R_{0}$ tended to be shorter, while outbreaks
with low $R_{0}$ tended to last longer. This trend was also discernible,
albeit less prominently, in the context of chikungunya epidemics, see Figure
S2 of the supplementary material.

Figure 6: Estimated duration of the dengue epidemic as a function of estimated
${\cal R}_{0}$.Longer epidemics are associated with lower estimates for R0.
Dots are individual cities.

###  3.2 Earlier epidemics and climate

The Figure 7 depicts the HGBR model’s performance in predicting the dengue
epidemic peak week as a function of climate and epidemiological covariates
measured one year before. The model provided more accurate predictions in the
South, Southeast and Midwest, where seasonality is stronger.

Figure 7: Performance of the Histogram Gradient Boosting regression model
trained with the entire dataset. Each point represents the predicted vs.
observed weeks for a city. The dashed black line represents the identity.

When we look at the importance of the model’s features (Figure 8), in all
regions, the features related to the peak week, population, and cases in the
previous year were always among the top ten most important features. Also, all
regions had at least two climate predictors in the top ten, with the average
temperature in January of the current year being the most common among them.

Figure 8: The top 10 most important predictors in each region, encoded as SHAP
values. The SHAP scale represents how much each feature affects the model’s
output (peak week), with positive values associated with an increase in peak
week and negative with a decrease. Each dot represents one epidemic. Their
colour is mapped to the value of each feature, with blue standing for lower
values and pink for higher. The 0 in the SHAP scale is aligned with the
expected value of each feature above.

Ocean surface temperature oscillations, as measured by the El niño/La niña
index (ENSO), also stood out as an important feature to explain variations in
the timing of epidemics. We looked specifically at the value of the ENSO index
in the last quarter of the previous year. Positive average ENSO values (El
Niño) in that quarter were associated with earlier epidemics in the following
year (Figure 9), while negative 4th quarter ENSO index led to later epidemics.

Figure 9: Estimated median peak week of dengue epidemics by region for El Niño
and La Niña years. La Niña years, in this analysis, are years when the mean
ENSO index in the last quarter of the previous year was below 0, and El niño
years are the ones where the mean ENSO index in the last quarter of the year
was above 0.

##  4 Discussion

This study introduces a computational pipeline (Figure 1) capable of
processing incidence data from all municipalities and estimating key epidemic
statistics in a standardized format. To achieve the required computational
scalability, we employed optimization techniques to fit the Richards model,
which has a one-to-one relationship with the archetypical SIR mathematical
model. Even though the SIR model is an approximation to the transmission
mechanisms taking place in the field, we believe that this simplified
description of the epidemics is a small price to pay in exchange for
developing a consistent set of epidemic parameters over a large spatio-
temporal range. Using more detailed transmission models would be undesirable,
not only because of the added computational cost to fit but also for the added
uncertainty associated with the estimation of a larger number of parameters
subject to identifiability issues (Gallo et al.,, 2022; Kao and Eisenberg,,
2018).

Following parameter estimation, a regression analysis was conducted to
elucidate the influence of local and global climate patterns on epidemic
parameters. Although not exhaustive, these exploratory analyses revealed a
discernible association between the El Niño Southern Oscillation (ENSO) and
epidemic timing (Figure 9), represented by the peak week of the epidemic. The
regression analysis results show that various epidemiological, demographic,
and climatic factors determine the shape of future epidemics. This association
has been reported before for Singapore (Earnest et al.,, 2012), Thailand
(Tipayamongkholgul et al.,, 2009), and other regions of the world (Johansson
et al.,, 2009; Moraes et al.,, 2019; Banu et al.,, 2015; Do et al.,, 2014).
However, our study differentiates itself by the larger spatio-temporal scale
of the analysis and for generating an open dataset that enables further
explorations and will be continuously updated.

The observed association between ${\cal R}_{0}$ and epidemic duration, shown
in figure S2, exhibits remarkable consistency, partly reflecting the
interdependence between these parameters within the transmission model. Given
the environmental influences driving epidemics, for example, the well-
described correlation between transmission and temperature, it is noteworthy
that these two parameters remain highly correlated, as per theoretical
expectations, that assume constant transmission rates and endogenous
determination of duration (by depletion of susceptibles). The small size of
the epidemics (relative to population size) suggests that such depletion is
not to be expected, except in the hypothesis of a massive underreporting of
cases in most cities, which could mask larger attack ratios.

The extensive dataset of epidemiological parameters generated by Episcanner
has demonstrated utility in developing predictive models, particularly when
integrated with complementary datasets such as climate and demographic series.
Although predictive modelling was only superficially explored in this study,
its potential for future research is considerable. To our knowledge, this
dataset represents a unique large-scale compilation of epidemic parameters
estimated using a consistent methodology, enabling robust comparisons. We
anticipate that statistical modellers will leverage this open dataset for
future investigations.

##  5 Conclusion

The methodological framework introduced in this paper facilitates the
comprehensive characterization of epidemic dynamics across both temporal and
spatial dimensions. This approach offers invaluable insights into the
transmission patterns of arbovirus diseases, enabling public health
authorities to devise timely interventions to mitigate disease spread
effectively. Furthermore, the identification of spatial variations in epidemic
parameters underscores the imperative for region-specific strategies in
disease surveillance and control initiatives.

The predictive potential for this parameter set has been demonstrated and can
be further elaborated to reveal more specific associations at smaller spatial
scales. One limitation of such models is the lack of data on the immunological
status of the population, i.e., the proportion of the population previously
exposed to each dengue serotype. Data about the proportion of cases infected
with each DENV type per week or at least per year would contribute
significantly to understanding these epidemics’ long-term dynamics. We believe
that making such data available would be the most important governmental
investment in controlling arbovirus diseases.

The epidemiological parameters used in this paper are updated weekly and are
made available through Episcanner, a publicly available web dashboard
(Coelho,, 2024).

Continuous monitoring of the epidemiological characteristics (described in
section 2.3) of dengue and chikungunya can be a key tool for efficient control
of their incidence and morbidity. When it comes to accurate epidemiological
assessments, methodological consistency over time and geographical space is
indispensable. In particular, for robust comparisons across time and space,
employing the same analytical methods is crucial to ensure that methodological
differences are not confounded with data variations.

This study elucidates longer-term association patterns between dengue
epidemics’ intensity and global climate variations, as exemplified by the
multivariate ENSO index. These findings reinforce the longstanding hypothesis
(Reiter,, 2001) that the ongoing escalation in global average temperatures
will inevitably exacerbate the burden of mosquito-borne diseases.

## Acknowledgment

This work was funded by a grant from the Wellcome Trust (226088/Z/22/Z).

## References

  * Banu et al., (2015) Banu, S., Guo, Y., Hu, W., Dale, P., Mackenzie, J. S., Mengersen, K., and Tong, S. (2015).  Impacts of el niño southern oscillation and indian ocean dipole on dengue incidence in bangladesh.  Scientific reports, 5(1):16105.
  * Codeco et al., (2018) Codeco, C., Coelho, F., Cruz, O., Oliveira, S., Castro, T., and Bastos, L. (2018).  Infodengue: A nowcasting system for the surveillance of arboviruses in Brazil.  Revue d’Épidémiologie et de Santé Publique, 66:S386.
  * Codeco et al., (2022) Codeco, C. T., Oliveira, S. S., Ferreira, D. A., Riback, T. I., Bastos, L. S., Lana, R. M., Almeida, I., Godinho, V. B., Cruz, O. G., and Coelho, F. C. (2022).  Fast expansion of dengue in Brazil.  The Lancet Regional Health–Americas, 12.
  * Coelho, (2024) Coelho, F. C. (2024).
  * de Almeida et al., (2023) de Almeida, I. F., Codeço, C. T., Lana, R. M., Bastos, L. S., de Souza Oliveira, S., da Cruz Ferreira, D. A., Godinho, V. B., Riback, T. I. S., Cruz, O. G., and Coelho, F. C. (2023).  The expansion of chikungunya in brazil.  The Lancet Regional Health–Americas, 25.
  * de Almeida et al., (2022) de Almeida, I. F., Lana, R. M., and Codeço, C. T. (2022).  How heterogeneous is the dengue transmission profile in brazil? a study in six brazilian states.  PLOS Neglected Tropical Diseases, 16(9):1–20.
  * de Souza et al., (2023) de Souza, W. M., de Lima, S. T., Mello, L. M. S., Candido, D. S., Buss, L., Whittaker, C., Claro, I. M., Chandradeva, N., Granja, F., de Jesus, R., et al. (2023).  Spatiotemporal dynamics and recurrence of chikungunya virus in brazil: an epidemiological study.  The Lancet Microbe, 4(5):e319–e329.
  * Do et al., (2014) Do, T. T. T., Martens, P., Luu, N. H., Wright, P., and Choisy, M. (2014).  Climatic-driven seasonality of emerging dengue fever in Hanoi, vietnam.  BMC public health, 14:1–10.
  * Earnest et al., (2012) Earnest, A., Tan, S., and Wilder-Smith, A. (2012).  Meteorological factors and el nino southern oscillation are independently associated with dengue infections.  Epidemiology & Infection, 140(7):1244–1251.
  * Gallo et al., (2022) Gallo, L., Frasca, M., Latora, V., and Russo, G. (2022).  Lack of practical identifiability may hamper reliable predictions in covid-19 epidemic models.  Science advances, 8(3):eabg5234.
  * Johansson et al., (2009) Johansson, M. A., Cummings, D. A., and Glass, G. E. (2009).  Multiyear climate variability and dengue—el nino southern oscillation, weather, and dengue incidence in puerto rico, mexico, and thailand: a longitudinal data analysis.  PLoS medicine, 6(11):e1000168.
  * Kao and Eisenberg, (2018) Kao, Y.-H. and Eisenberg, M. C. (2018).  Practical unidentifiability of a simple vector-borne disease model: Implications for parameter estimation and intervention assessment.  Epidemics, 25:89–100.
  * Lundberg and Lee, (2017) Lundberg, S. M. and Lee, S.-I. (2017).  A unified approach to interpreting model predictions.  Advances in neural information processing systems, 30.
  * Manzoor et al., (2022) Manzoor, K. N., Javed, F., Ejaz, M., Ali, M., Mujaddadi, N., Khan, A. A., Khattak, A. A., Zaib, A., Ahmad, I., Saeed, W. K., and Manzoor, S. (2022).  The global emergence of chikungunya infection: An integrated view.  Reviews in Medical Virology, 32(3):e2287.
  * Moraes et al., (2019) Moraes, B. C. d., Souza, E. B. d., Sodré, G. R. C., Ferreira, D. B. d. S., and Ribeiro, J. B. M. (2019).  Sazonalidade nas notificações de dengue das capitais da amazônia e os impactos do el niño/la niña.  Cadernos de Saúde Pública, 35:e00123417.
  * Muñoz Sabater and Service, (2019) Muñoz Sabater, J. and Service, C. C. C. (2019).  Era5-land hourly data from 1950 to present.
  * Reiter, (2001) Reiter, P. (2001).  Climate change and mosquito-borne disease.  Environmental health perspectives, 109(suppl 1):141–161.
  * Rodriguez-Barraquer et al., (2011) Rodriguez-Barraquer, I., Cordeiro, M. T., Braga, C., de Souza, W. V., Marques, E. T., and Cummings, D. A. (2011).  From re-emergence to hyperendemicity: the natural history of the dengue epidemic in brazil.  PLoS neglected tropical diseases, 5(1):e935.
  * Storn and Price, (1997) Storn, R. and Price, K. (1997).  Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces.  Journal of global optimization, 11(4):341–359.
  * Tipayamongkholgul et al., (2009) Tipayamongkholgul, M., Fang, C.-T., Klinchan, S., Liu, C.-M., and King, C.-C. (2009).  Effects of the el niño-southern oscillation on dengue epidemics in thailand, 1996-2005.  BMC public health, 9:1–15.
  * Wang et al., (2012) Wang, X.-S., Wu, J., and Yang, Y. (2012).  Richards model revisited: Validation by and application to infection dynamics.  Journal of theoretical biology, 313:12–19.
  * Yang et al., (2021) Yang, X., Quam, M. B. M., Zhang, T., and Sang, S. (2021).  Global burden for dengue and the evolving pattern in the past 30 years.  Journal of Travel Medicine, 28(8):taab146.

Table S1: Description of the features used in the Histogram Gradient Boosting Regressor model. The climatic data are obtained daily, and the number of cases is obtained weekly. In the features identified with long tails as cases and rainy days, a log transformation was applied to improve the model performance. Feature description |  Type
---|---
year: year of the peak week that are being predicted. |  temporal
casos_01: sum of cases in the January of the year whose peak week are being predicted. |  epidemiological
casos_1_3: sum of cases on the third quarter of the previous year. |  epidemiological
casos_1_4: sum of cases on the fourth quarter of the previous year. |  epidemiological
populacao_1: population in the previous year. |  demographic
peak_week_1: peak week estimated on the previous year. |  epidemiological
R0_1: reproduction number estimated on the previous year. |  epidemiological
ep_dur_1: epidemic duration estimated on the previous year. |  epidemiological
dummy_ep: 1 if the previous year were an epidemic identified by Episcanner and 0 otherwise. |  epidemiological
temp_med_4: average of the average temperature over the fourth quarter of the previous year. |  climatic
temp_amp_4: average of the temperature amplitude over the fourth quarter of the previous year. |  climatic
temp_max_4: average of the maximum temperature over the fourth quarter of the previous year. |  climatic
temp_min_4: average of the minimum temperature over the fourth quarter of the previous year. |  climatic
umid_min_4: average of the minimum humidity over the fourth quarter of the previous year. |  climatic
umid_max_4: average of the maximum humidity over the fourth quarter of the previous year. |  climatic
umid_amp_4:average of the humidity amplitude over the fourth quarter of the previous year. |  climatic
enso_4: average of the multivariate ENSO (El Niño-Southern Oscillation) in the fourth quarter of the previous year. |  climatic
precip_tot_4: sum of the total precipitation over the fourth quarter of the previous year. |  climatic
rainy_day_4: sum of days with rain (precipitation above zero) over the fourth quarter of the previous year. |  climatic
thr_temp_min_4: sum of days with minimum temperature below 15-celsius degrees |  climatic
thr_temp_amp_4: sum of days with temperature amplitude above celsius degrees. |  climatic
thr_umid_med_4: sum of days with average humidity above 0.8. |  climatic
temp_med_1_current: average of the average temperature over the january of the year whose peak week is being predicted. |  climatic
temp_amp_1_current: average of the temperature amplitude over the january of the year whose peak week is being predicted. |  climatic
temp_max_1_current: average of the maximum temperature over the january of the year whose peak week is being predicted. |  climatic
temp_min_1_current: average of the minimum temperature over the january of the year whose peak week is being predicted. |  climatic
precip_tot_1_current: sum of total precipitation over the january of the year whose peak week is being predicted. |  climatic
rainy_day_1_current: sum of days with precipitation over the january of the year whose peak week is being predicted. |  climatic
enso_1_current: average of the ENSO (El Niño-Southern Oscillation) over the January of the year whose peak week is being predicted. |  climatic
latitude: latitude of the city center. |  spatial
longitude: longitude of the city center. |  spatial
Figure S1: Median of the reproduction number ($R_{0}$) estimated for Dengue
and Chikungunya cases. In the middle map, N refers to the North, NE to the
North East, MW to the Midwest, SE to the Southeast, and S to the South. Figure
S2: Duration of the chikungunya epidemic as a function of R0. Longer epidemics
are associated with lower estimates for R0.


# Photo-induced Ferromagnetic and Superconducting Orders in Multi-orbital
Hubbard Models

Sujay Ray  Philipp Werner  Department of Physics, University of Fribourg, 1700
Fribourg, Switzerland

###### Abstract

The search for hidden orders in photoexcited lattice systems is an active
research field driven by experimental reports of light-induced or light-
stabilized phases. In this study, we investigate hidden electronic orders in
strongly correlated two-orbital Hubbard models with orbital-dependent
bandwidths. In equilibrium, the half-filled systems are antiferromagnetically
ordered. Using non-equilibrium dynamical mean field theory we demonstrate the
appearance of nonthermal ferromagnetic order in the photo-doped state, if the
two bandwidths are sufficiently different, and its coexistence with spin-
singlet $\eta$-superconductivity in the high photo-doping region. Spin-triplet
$\eta$-superconducting order appears instead if the two bandwidths are
comparable. The rich nonequilibrium phasediagram uncovered in this work shows
that Mott insulating multi-orbital systems provide an interesting platform for
the realization of nonthermal electronic orders.

Introduction: The exploration of metastable or hidden phases, which cannot be
accessed under equilibrium conditions, is an exciting research direction in
condensed matter physics lat_struc ; mang_hid ; Fausti2011 ; Cavalleri2018 ;
Torre2021 ; ph_noneq . Due to the interplay between different active degrees
of freedom Dagotto2005 , strongly correlated electron systems are particularly
interesting candidates for the realization of such hidden phases yuta_etasc ;
yuta_etasc1 ; opt_magord ; ind_sc_cuprate ; phind_sc . Among them, Mott
insulators have been shown to host different types of hidden phases opt_exint
; eta_jiajun ; eta_spintrip ; tri_chiral ; exciton_insulator ; Kugel_Khomskii
when excited by a laser pulse. Recent theoretical studies of photo-doped Mott
insulators identified nonthermal staggered superconducting (SC) orders in
single-band and two-band Hubbard models yuta_etasc ; yuta_etasc1 ; kaneko1 ;
kaneko2 ; eta_jiajun ; eta_spintrip  and chiral superconducting orders in
geometrically frustrated systems tri_chiral . Two-orbital Hubbard models have
also been shown to exhibit nonthermal excitonic order exciton_insulator  and
Kugel-Khomoskii spin-orbital order Kugel_Khomskii . Because of the long
lifetime of charge excitations in large-gap Mott insulators Sensarma2010 ;
Eckstein2011 ; Strohmaier2010 ; Lenarcic2014 , a partial thermalization of
charge carriers may occur in each Hubbard band, and these photo-doped
metastable states can be susceptible to new types of ordering instabilities.
Due to the rapid advancement in laser technology, it becomes possible to
induce and detect such metastable states VO2_dynamics ;
ferro_chrgtrnsfr_dynamics , which also provides interesting perspectives for
technological applications.

Several experimental and theoretical studies considered magnetic orders in
photoexcited systems mag_switch ; opt_magord ; manganites_1 ; VO2_dynamics ;
ferro_chrgtrnsfr_dynamics ; Qint_Mott ; Ishihara_2007 ; Ishihara_2011 ;
Ishihara_DE1 ; Ishihara_DE2 ; Qdyn_spin_struc ; opt_magcont . The phenomena
range from ultrafast demagnetization ultrafast_demag1 ; ultrafast_demag2  to
magnetic switching and the appearance of nonthermal magnetic orders mag_switch
; Kugel_Khomskii . The double exchange model plays a prominent role in studies
of photo-induced magnetic phases. In the context of perovskite manganites, a
series of numerical investigations based on exact diagonalization and density
matrix renormalization group calculations reported nonthermal magnetic phases
manganites_1 ; manganites_2 ; manganites_3 ; Ishihara_DE1 ; Ishihara_DE2 ;
Ishihara_DE3 . However, these studies of closed systems lack a systematic
control over the density and effective temperature of the photo-carriers.
While the nonthermal melting of Kugel-Khomskii spin-orbital order in a photo-
doped system has been investigated with non-equilibrium dynamical mean field
theory (DMFT) Kugel_Khomskii , the emergence of magnetic order is difficult to
study in real-time simulations, due to heating effects and intrinsically slow
dynamics.

Here we use the Nonequilibrium Steady State (NESS) approach NESS_martin ,
which provides a high level of control over the effective temperature and
density of photo-carriers, to study the two-orbital Hubbard model with
orbital-dependent bandwidths. The interplay between more localized and more
itinerant electrons leads to nontrivial couplings between the spin, orbital,
and charge degrees of freedom and results in a multitude of hidden electronic
orders.

Model and Method: We consider the two-orbital Hubbard model with Hamiltonian

$\displaystyle
H=-\sum_{\left<ij\right>,\sigma}\sum_{\alpha=1,2}t^{\text{hop}}_{\alpha}c^{\dagger}_{i,\alpha\sigma}c_{j,\alpha\sigma}$
$\displaystyle+U\sum_{i}\sum_{\alpha=1,2}n_{i,\alpha\uparrow}n_{i,\alpha\downarrow}-\mu\sum_{i}\sum_{\alpha=1,2}\left(n_{i,\alpha\uparrow}+n_{i,\alpha\downarrow}\right)$
$\displaystyle+(U-2J)\sum_{i,\sigma}n_{i,1\sigma}n_{i,2\bar{\sigma}}+(U-3J)\sum_{i,\sigma}n_{i,1\sigma}n_{i,2\sigma},$
(1)

where $c^{\dagger}_{i,\alpha\sigma}$ $(c_{i,\alpha\sigma})$ is the creation
(annihilation) operator at site $i$, $\alpha$ and $\sigma$ denote the orbital
and spin indices, respectively, $t^{\text{hop}}_{\alpha}$ is the orbital
dependent nearest neighbor hopping amplitude between sites $i$ and $j$, $U$ is
the intra-orbital Hubbard repulsion, $J$ the Hund coupling, and $\mu$ the
chemical potential. In general, we assume different hopping parameters for the
two orbitals and hence different bandwidths. We are interested in the strongly
correlated regime of the half filled Hubbard model, where both bands are Mott
insulating, see the upper panel of Fig. 1. When a laser pulse is applied to
such a system, the wider band gets more strongly photo-doped.

In the NESS simulations NESS_martin , the system is weakly coupled to cold
Fermion baths at each site. These baths induce charge carriers (triply and
singly occupied sites) into our half-filled system, thus creating an effective
photoexcited state. The Fermion baths are aligned with the lower and the upper
Hubbard bands. By choosing the chemical potential of the Fermion baths
$\mu_{b}$, we can control the density of photo-carriers, while changing the
temperature of the Fermion baths $T_{b}$ allows to control the effective
temperature of the photo-doped system. We can thus map out the nonequilibrium
phase diagram of the photo-excited system as a function of photo-doping
concentration and effective temperature.

The DMFT simulations of the NESS are implemented for an infinitely connected
Bethe lattice using the non-crossing approximation (NCA) as impurity solver
Keiter_1970 ; Eckstein_2010 . The DMFT self-consistency equation in the Nambu
Keldysh formalism takes the form

$\displaystyle\Delta(t,t^{\prime})=t_{0}^{2}\gamma
G(t,t^{\prime})\gamma+\sum_{b}D_{b}(t,t^{\prime}).$ (2)

Here $\Delta$ is the hybridization function, which has two contributions. The
first term represents the lattice self-consistency for the Bethe lattice, with
$t_{0}=W/2$, where $W$ is the half-bandwidth. The second term comes from the
Fermionic baths and is defined in frequency space by $D_{b}(\omega)=$
$g^{2}\rho_{b}(\omega)=\Gamma\sqrt{W_{b}^{2}-(\omega-\omega_{b})^{2}}$, where
$\Gamma=g^{2}/W_{b}^{2}$ is a dimensionless coupling constant, $\omega_{b}$
indicates the center of the energy spectrum and $W_{b}$ indicates the half-
bandwidth of the bath $b$. In Eq. (2) we use the spinor basis
$\psi^{\dagger}=(c^{\dagger}_{1\uparrow}\hskip
4.26773ptc^{\dagger}_{2\uparrow}\hskip 4.26773ptc_{1\downarrow}\hskip
4.26773ptc_{2\downarrow})$ to study spin-singlet SC order and
$\psi^{\dagger}=(c^{\dagger}_{1\uparrow}\hskip 4.26773ptc_{2\uparrow}\hskip
4.26773ptc^{\dagger}_{1\downarrow}\hskip 4.26773ptc_{2\downarrow})$ for spin-
triplet SC order. The choice of the spinor $\psi^{\dagger}$ and the $\gamma$
matrix allows us to study different magnetic and SC phases, see Supplemental
Material (SM). Because of the steady state assumption, the hybridization
functions and Green’s functions depend only on the time difference
$t-t^{\prime}$.

Figure 1:  Schematic diagram showing the photo-excitation in two-orbital Mott
insulators with different bandwidths. After the applied laser pulse, the wider
band is photo-doped more strongly than the narrower band. $t^{\text{hop}}_{1}$
and $t^{\text{hop}}_{2}$ are the hopping parameters for orbital $1$ and $2$,
respectively. The bottom panels show the dominant hopping processes at
holon/triplon density 0.25, which stabilize the FM state.

In all our calculations we use four fermion baths with $\omega_{b}=\pm U/2,\pm
3U/2$ and $W_{b}=3.0$ and $\Gamma=0.028$ for efficient photo-doping and we use
the interaction parameters $U=20$ and $J=1$. The hopping parameter for orbital
$1$ is kept fixed at $t^{\text{hop}}_{1}=1$, while $t^{\text{hop}}_{2}$ for
orbital $2$ is varied from $0.1$ to $1.0$. We measure the photo-doping
concentration by the triplon density, which we define as
$\left<n_{1\uparrow}n_{1\downarrow}n_{2}+n_{1}n_{2\uparrow}n_{2\downarrow}-4n_{1\uparrow}n_{1\downarrow}n_{2\uparrow}n_{2\downarrow}\right>$,
with $n_{\alpha}=n_{\alpha\uparrow}+n_{\alpha\downarrow}$ the occupation of
orbital $\alpha$. In addition, we define the magnetic orders via the local
magnetization
$m_{i}=\sum_{\alpha}n_{i,\alpha\uparrow}-n_{i,\alpha\downarrow}$, and two
types of SC order - spin-singlet orbital-triplet order
$p_{i,oy}=\frac{1}{4}\sum_{\alpha}(c^{\dagger}_{i,\alpha\uparrow}c^{\dagger}_{i,\alpha\downarrow}+h.c.)$,
and spin-triplet orbital-singlet SC order
$p_{i,sy}=\frac{1}{4}\sum_{\sigma}(c^{\dagger}_{i,1\sigma}c^{\dagger}_{i,2\sigma}+h.c.)$.
In a staggered SC state ($\eta$ order) or staggered antiferromagnetic state
(AFM order) these order parameters have opposite signs on the two sublattices
$A$ and $B$,

$\displaystyle\eta^{x}_{i,oy(sy)}$ $\displaystyle=$
$\displaystyle\delta_{i}^{A/B}p_{i,oy(sy)},$ (3) $\displaystyle M_{i,AFM}$
$\displaystyle=$ $\displaystyle\delta_{i}^{A/B}m_{i},$ (4)

with $\delta_{i}^{A(B)}=1$ $(-1)$ for $i$ on the $A$ $(B)$ sublattice.

In order to study different phases we apply small seed fields which couple to
the order parameters of interest. In a symmetry-broken state the order
parameter grows to a high value and becomes almost independent of the value of
the seed field. In the DMFT calculations, we switch off the seed field after
some iterations, so that the symmetry-broken phase can be unambiguously
identified by a non-zero order-parameter. We cannot simultaneously realize
$\eta^{x}_{i,oy}$ and $\eta^{x}_{i,sy}$ SC orders, since they correspond to
different $\psi^{\dagger}$ and $\gamma$, while FM and $\eta^{x}_{i,oy}$ order
can be simultaneously stabilized.

Figure 2: (a) Nonequilibrium phase diagram of the photo-doped half-filled two-
orbital Hubbard model with different bandwidths ($U=20$, $J=1$) as a function
of the triplon density and hopping parameter $t^{\text{hop}}_{2}$ of the
narrow band ($t^{\text{hop}}_{1}=1$). AFM order (an extension of the
equilibrium phase) appears for small triplon density, whereas a FM state is
stabilized for intermediate photo-doping and low $t^{\text{hop}}_{2}$. For
higher photo-doping a spin-singlet $\eta$-SC order is stabilized for smaller
values of $t^{\text{hop}}_{2}$ and spin-triplet $\eta$-SC appears for very
strong photo-doping when $t^{\text{hop}}_{2}\sim t^{\text{hop}}_{1}$. (b) FM
order parameter as a function of Hund coupling $J$ at triplon density $0.27$
for $t^{\text{hop}}_{2}=0.3$. The vanishing magnetization at small $J$ shows
the role of Hund coupling in stabilizing the FM state. (c) Magnetization and
$\eta^{x}_{oy}$ order in band $1$ and band $2$ for $t^{\text{hop}}_{2}=0.3$.
The magnetization is dominated by band $2$, while the $\eta$-SC order is
hosted by band $1$. The effective temperature in all calculations is kept at
$\beta_{\text{eff}}\sim 55$.

Photo-induced hidden orders: We first study the nonequilibrium phase diagram
as a function of photo-doping (triplon density) and the hopping parameter of
the narrow band ($t^{\text{hop}}_{2}$), see Fig. 2. To probe the ferromagnetic
(FM) and spin-singlet $\eta$-SC states, we add the seed terms
$H_{\text{seed}}=B_{\text{seed}}\sum_{i}m_{i}+P_{\text{seed}}\sum_{i}(c^{\dagger}_{i,1\uparrow}c^{\dagger}_{i,1\downarrow}+c^{\dagger}_{i,2\uparrow}c^{\dagger}_{i,2\downarrow}+h.c.)$
to our Hamiltonian (1), where $B_{\text{seed}}$ and $P_{\text{seed}}$ are,
respectively, small magnetic and superconducting seed fields. Due to the
partial thermalization of the electrons in each Hubbard band, the photodoped
state can be represented by an equilibrium state with fixed number of
multiplets (triplons, doublons and singlons), described by the effective low
energy model SWT ; eta_spintrip

$\displaystyle H_{\text{eff}}$ $\displaystyle=$ $\displaystyle
H_{\text{hop}}(P_{in}P_{jn-1})+H_{s/o}(P_{in}P_{jn})+H_{\eta}(P_{in}P_{jn-2}),$

where $P_{in}$ is the projection operator at site $i$ on states with particle
number $n$. The first term in Eq. (2) is the hopping term
$-\sum_{\left<ij\right>n}\sum_{\alpha\sigma}t^{\text{hop}}_{\alpha}c^{\dagger}_{i,\alpha\sigma}c_{j,\alpha\sigma}$
which acts on sites which differ in particle number by one. The second term is
the spin-orbital term
$2(t^{\text{hop}}_{\alpha})^{2}/U\sum_{\left<ij\right>}\left(\frac{1}{2}+2\bm{s_{i}}.\bm{s_{j}}\right)\left(\frac{1}{2}+2\bm{\tau_{i}}.\bm{\tau_{j}}\right)$
and acts on sites with the same particle numbers, where $\bm{s_{i}}$ and
$\bm{\tau_{i}}$ are spin and orbital pseudospin operators at site $i$. Hund
coupling $J>0$ ($J<0$) favors spin (orbital) order. The last term is the
$\eta$-pseudospin term given by
$-4(t^{\text{hop}}_{\alpha})^{2}/U\sum_{\left<ij\right>,\nu}\left(\eta_{i\nu}^{+}\eta_{j\nu}^{-}+\eta_{i\nu}^{-}\eta_{j\nu}^{+}\right)$
with $\nu$ an index for the six different SC orders (three spin-singlet
orbital-triplet and three spin-triplet orbital-singlet orders, see SM). This
term acts on sites which differ in particle number by two eta_spintrip .

Near zero photo-doping, there are negligible charge excitations and the system
is governed by the $H_{s/o}$ term, which stabilizes the antiferromagnetic
(AFM) phase with a transition temperature ($T_{c}$) of the order of
$(t^{\text{hop}}_{\alpha})^{2}/U$. As the photo-doping is increased, the
diffusing charge carriers melt the AFM phase, which results in a nonthermal PM
metal. As the bandwidth of the narrower band ($t^{\text{hop}}_{2}$) is
decreased, the $T_{c}$ also decreases, since we are on the strongly correlated
side of the AFM dome. This shifts the boundary between the AFM and PM phases
towards lower photo-doping. At intermediate photo-doping, the hopping term
$H_{\text{hop}}$ starts to dominate, because of the substantial density of
singly occupied (singlon), doubly occupied (doublon) and triply occupied
(triplon) sites. If we decrease $t^{\text{hop}}_{2}$ in this regime to
localize the electrons in the narrower band, a nonthermal FM state is
stabilized by the interplay of Hund coupling and kinetic energy. A positive
Hund coupling $J$ favors high-spin doublon states. On the other hand, the
hopping processes in the wide band from doublon to singlon and from triplon to
doublon sites, illustrated in the bottom panel of Fig. 1, stabilize the FM
state since this state allows to gain kinetic energy while keeping the local
spin arrangement unchanged. Any other hopping processes would create a
configuration with higher energy. Thus a nonthermal FM phase appears at lower
$t^{\text{hop}}_{2}$ values and near triplon/singlon density 0.25 for $J>0$.
In the absence of Hund coupling, all the local doublon states are degenerate
in energy and the hopping term is not expected to favor the FM state. Indeed,
by decreasing the Hund coupling, the FM magnetization disappears below a
critical value of $J$, as shown in Fig. 2(b).

Figure 3: NESS phase diagram of the two-orbital Hubbard model with $U=20$,
$J=1$, $t^{\text{hop}}_{1}=1$ and $t^{\text{hop}}_{2}=0.4$ in the space of
$\beta_{\text{eff}}$ and triplon density. The AFM phase appears in the undoped
and weakly photo-doped region, while the $\eta$ SC state is stabilized in the
high photo-doping region. At intermediate photo-doping we have the FM state.
In these simulations, $T_{b}$ is varied to achieve different
$\beta_{\text{eff}}$. The inset of the figure shows the FM susceptibility
$\chi_{\text{FM}}$, estimated from the real-time calculations, as a function
of doping. The optimal doping condition for FM order is consistent in both
calculations.

When the photo-doping is increased further, more charge excitations are
created and for $t^{\text{hop}}_{1}=t^{\text{hop}}_{2}=1$, in the high doping
limit (triplon density $\gtrsim 0.46$) spin-triplet orbital-singlet $\eta$-SC
order ($\eta^{x}_{sy}$) develops eta_spintrip . As shown in Fig. 2(a), this
order disappears when $t^{\text{hop}}_{2}$ is decreased. For small
$t^{\text{hop}}_{2}$, the charge excitations are mostly concentrated in the
wider band, while the narrower band remains close to half-filled (e. g., for
$t^{\text{hop}}_{1}=1$, $t^{\text{hop}}_{2}=0.3$, triplon density $\sim 0.28$
and $\beta_{\text{eff}}\sim 55$, the double occupancies in the wider band and
the narrower band are 0.25 and 0.03, respectively). At sufficiently low
effective temperature, the doublons and singlons in the wide band condense to
give rise to spin-singlet orbital-triplet $\eta_{oy}$ order. Since the
dominant contribution to the $\eta_{oy}$ order comes from the wider band
($c^{\dagger}_{i,1\uparrow}c^{\dagger}_{i,1\downarrow}$), we have an orbital
selective $\eta_{oy}$ order (see Fig. 2(c)).

Interestingly, in the photo-doping region with triplon density between $0.3$
to $0.4$, both the first and last term in Eq. (2) contribute and we find a
coexistence of FM and $\eta$-SC order. At first sight, this is surprising,
since magnetic order is usually detrimental to singlet pairing. Here, such a
co-existence is possible because the two photo-induced order parameters are
supported by different orbitals – the wider band hosts the SC order, while the
FM order is dominated by the spin-aligned localized electrons in the narrower
band, as shown in Fig. 2(c).

Next, we fix the $t^{\text{hop}}_{1}=1,t^{\text{hop}}_{2}=0.4$ and investigate
the phase diagram in the space of triplon density and inverse effective
temperature $\beta_{\text{eff}}$, see Fig. 3. At low $T_{\text{eff}}$ (high
$\beta_{\text{eff}}\sim 50$), we have AFM order at low photo-doping. Around
triplon density 0.25 the FM phase appears, while at large photo-doping, we
find the $\eta_{oy}$-SC phase, which partly coexists with the FM phase. With
increasing temperature (decreasing $\beta_{\text{eff}}$) the ordered regions
shrink and they disappear at similar $T_{c}$ values. The highest effective
$T_{c}$ for the nonthermal FM order is around triplon density $\sim 0.29$,
which is consistent with the optimal doping deduced from the phase diagram in
Fig. 2(a). The $\eta_{oy}$ phase shows a similar temperature dependence as in
the one band Hubbard model eta_jiajun , consistent with our observation that
the $\eta_{oy}$-SC state is essentially decoupled from the narrow band and has
a single-band characteristics.

We see from Fig. 3 that both the nonthermal FM and $\eta$-SC orders appear at
$\beta_{\text{eff}}>30$. While this is a high temperature for realistic
parameters, it is difficult to realize these effective temperatures in real-
time simulations of a photoexcited multi-orbital system, even using the
entropy cooling technique entropy_cooling . However, near the optimal photo-
doping density, we still expect to find an enhanced FM susceptibility at the
accessible higher $T_{\text{eff}}$.

Figure 4: (a) Illustration of the entropy cooling technique, where electrons
are injected into the upper Hubbard bands from a fully occupied narrow Fermion
bath and electrons are ejected from the lower Hubbard bands to an empty
Fermion bath. (b) Schematic diagram of the four-impurity DMFT calculation
mimicking the dynamics in two bilayer systems – an AFM system (bottom) and an
AFM stacked FM system (top). (c) Magnetic susceptibility
$\chi_{\text{FM/AFM}}$ in the FM and AFM state (red curves) and time-evolution
of various local observables – triplon density, double occupancy in orbital
$1$ ($d_{1}$) and orbital $2$ ($d_{2}$). The yellow curve shows the hopping
pulse $v_{sb}(t)$ (not to scale) used to produce the photo-doping.

Real-time simulations: To investigate the dynamics and measure the magnetic
suceptibilities, we employ the entropy cooling technique entropy_cooling ;
entropy_cooling1  to minimize the heating effect in nonequilibrium photo-doped
Mott insulators. In this setup, we start with a half-filled Mott insulating
state and add two additional narrow noninteracting Fermion baths (one
completely full and one completely empty) as shown in Fig. 4(a). These baths
are coupled to the Mott insulating system by a hopping term $v_{sb}$. To
facilitate an efficient photo-doping, while keeping the effective temperature
of the system low, the hopping term is modulated in time as
$v_{sb}(t)=A_{0}f_{\text{en}}(t)\sin(\Omega(t)t)$, where $A_{0}=1$ denotes the
amplitude of the hopping, $f_{\text{en}}(t)$ is an envelope function
determining the switch-on and switch-off times of the hopping and $\Omega(t)$
is a time-dependent (chirped) frequency. In our calculations with $U=20$,
$J=1$ model, we chose $f_{\text{en}}(t)=\exp(-((t-t_{0})/2\sigma)^{2})$ as a
Gaussian with $t_{0}=12$ and $\sigma=2.8$. $\Omega(t)$ is increased linearly
from an initial value $\Omega_{\text{ini}}=27.0$ to a final value
$\Omega_{\text{fin}}=29.8$. The hopping modulation pulse is shown in Fig. 4(c)
(yellow curve). In this set-up, the DMFT self-consistency condition becomes

$\displaystyle\Delta(t,t^{\prime})=t_{0}(t)^{2}\gamma
G(t,t^{\prime})\gamma+\sum_{\alpha}v_{sb}(t)^{2}G_{\text{bath},\alpha}(t,t^{\prime})\hskip
14.22636pt$ (6)

where $t_{0}(t)=t_{0}$ is kept fixed, $G(t,t^{\prime})$ is the Green’s
function of the system and $G_{\text{bath},\alpha}(t,t^{\prime})$ is the
Green’s function of the noninteracting bath with flavor $\alpha$ (full,empty).

Because of the conservation of total spin in our Hamiltonian, it is not
possible to simulate the switching to a purely FM state within this framework.
Thus, in order to study the photo-induced FM state, we implemented a four-
impurity calculation, which mimics the bi-layer system shown in Fig. 4(b).
Here, the impurity models $1$ and $2$ represent two neighboring sites in layer
1, whereas the impurity models $3$ and $4$ represent the corresponding sites
in layer 2. In this set-up, an AFM state is identified by AFM spin alignments
in all four pairs of neighboring impurities, while a FM state is indicated by
a layer-dependent magnetization, e. g. with spin-up polarization on impurities
$1$ and $2$ and spin-down polarization on impurities $3$ and $4$ (see Fig.
4(b)).

In Fig. 4(c) we show the result of the real-time DMFT calculation. As soon as
the coupling to the narrow baths is switched on, the system gets photo-doped,
as indicated by the increasing triplon density during the pulse. As a result
of the different bandwidths, the double occupancy in orbital 1 ($d_{1}$) grows
larger than in orbital 2 ($d_{2}$), which is consistent with the results of
the NESS calculations. The magnetic susceptibility is measured by dividing the
magnetization (magnetic order parameter) by the corresponding seed field as
$\chi=\left<m\right>/B_{\text{seed}}$. In Fig. 4(c), we show the
$\chi_{\text{FM}}$ and $\chi_{\text{AFM}}$ estimates for FM and AFM order. In
the FM case, we choose $B_{\text{seed}}>0$ for impurity 1 and 2,
$B_{\text{seed}}<0$ for impurity 3 and 4, while in the AFM case
$B_{\text{seed}}>0$ for impurity 1 and 4 and $B_{\text{seed}}<0$ for impurity
2 and 3. The FM susceptibility grows larger than the AFM one and keeps
increasing at the longest simulation times, indicating a dominant tendency
towards FM ordering. We can also control the photo-doping concentration by
controlling the final frequency $\Omega_{\text{fin}}$ of the chirped pulse and
extract $\chi_{\text{FM}}$ as a function of the triplon density, as shown in
the inset of Fig. 3. We find that $\chi_{\text{FM}}$ reaches a maximum value
near triplon density $\sim 0.3$, which agrees with the optimal photodoping
concentration for the FM phase in the NESS simulations. Due to numerical
limitations, we cannot drive the system slow enough to reduce $T_{\text{eff}}$
further or drive the system for long enough time to demonstrate an actual
symmetry breaking. But the close correlation between the FM phase boundary in
the NESS approach and the FM susceptibility in the real-time simulations
supports the appearance of a nonthermal FM phase in photo-doped two-orbital
Mott insulators.

In conclusion, we demonstrated that the photo-doped two-band Hubbard model
with orbital-dependent bandwidths hosts nonthermal FM and staggered SC order,
and even a co-existence of both orders. Manganites, which are often described
by the double exchange model, and NiO, with two Mott insulating half-filled
$e_{g}$ bands, would be possible candidates to search for these nonthermal
electronic orders.

Acknowledgements: We thank Y. Murakami for helpful comments on the manuscript
and the Swiss National Science Foundation for the funding (Grant No.
200021-196966). The nonequilibrium DMFT calculations are based on the NESSi
library NESSi  and the NESS simulations on a code originally developed by J.
Li and M. Eckstein. The calculations were run on the beo06 cluster at the
University of Fribourg.

## References

  * (1) H. Ichikawa, S. Nozawa, T. Sato, A. Tomita, K. Ichiyanagi, M. Chollet, L. Guerin, N. Dean, A. Cavalleri, S.-i. Adachi, et al., Nat. Mater. 10, 101 (2011).
  * (2) L. Stojchevska, I. Vaskivskyi, T. Mertelj, P. Kusar, D. Svetin, S. Brazovskii, and D. Mihailovic, Science 344, 177 (2014).
  * (3) D. Fausti, R. I. Tobey, N. Dean, S. Kaiser, A. Dienst, M. C. Hoffmann, S. Pyon, T. Takayama, H. Takagi, and A. Cavalleri, Science 331, 189 (2011).
  * (4) A. Cavalleri, Contemp. Phys. 59, 31 (2018).
  * (5) A. de la Torre, D. M. Kennes, M. Claassen, S. Gerber, J. W. McIver, and M. A. Sentef, Rev. Mod. Phys. 93, 041002 (2021).
  * (6) Y. Murakami, D. Golež, M. Eckstein, P. Werner, arXiv:2310.05201 (2023).
  * (7) E. Dagotto, Science 309, 257 (2005).
  * (8) Y. Murakami, S. Takayoshi, T. Kaneko, A. M. Läuchli, and P. Werner, Phys. Rev. Lett. 130, 106501 (2023).
  * (9) Y. Murakami, S. Takayoshi, T. Kaneko, Z. Sun, D. Golež, A. J. Millis, and P. Werner, Comm. Phys. 5, 23 (2022).
  * (10) A. Kirilyuk, A. V. Kimel, and T. Rasing, Rev. Mod. Phys. 82, 2731 (2010).
  * (11) D. Fausti, R. I. Tobey, N. Dean, S. Kaiser, A. Dienst, M. C. Hoffmann, S. Pyon, T. Takayama, H. Takagi, and A. Cavalleri, Science 331, 189 (2011).
  * (12) M. Mitrano, A. Cantaluppi, D. Nicoletti, S. Kaiser, A. Perucchi, S. Lupi, P. Di Pietro, D. Pontiroli, M. Ricco, S. R. Clark, et al., Nature 530, 461 (2016).
  * (13) P. Werner, and Y. Murakami, Phys. Rev. B 102, 241103(R) (2020).
  * (14) R. Mikhaylovskiy, E. Hendry, A. Secchi, J. H. Mentink, M. Eckstein, A. Wu, R. Pisarev, V. Kruglyak, M. Katsnel-son, T. Rasing, et al., Nat. Commun. 6, 8190 (2015).
  * (15) J. Li, M. Müller, A. Kim, A. Läeuchli, and P. Werner, Phys. Rev. B 107, 205115 (2023).
  * (16) J. Li, H. U. R. Strand, P. Werner and M. Eckstein , Nat. Com. 9, 4581 (2018).
  * (17) J. Li, D. Golez, P. Werner, and M. Eckstein, Phys. Rev. B 102, 165136 (2020)
  * (18) S. Ray, Y. Murakami, and P. Werner, Phys. Rev. B 108, 174515 (2023).
  * (19) S. Ejima, T. Kaneko, F. Lange, S. Yunoki, and H. Fehske, Phys. Rev. Research 2, 032008(R) (2020).
  * (20) T. Kaneko, S. Yunoki, and A. J. Millis, Phys. Rev. Research 2, 032027(R) (2020).
  * (21) R. Sensarma, D. Pekker, E. Altman, E. Demler, N. Strohmaier, D. Greif, R. Jördens, L. Tarruell, H. Moritz, and T. Esslinger, Phys. Rev. B 82, 224302 (2010).
  * (22) M. Eckstein and P. Werner, Phys. Rev. B 84, 035122 (2011).
  * (23) N. Strohmaier, D. Greif, R. Jördens, L. Tarruell, H. Moritz, T. Esslinger, R. Sensarma, D. Pekker, E. Alt- man, and E. Demler, Phys. Rev. Lett. 104, 080401 (2010).
  * (24) Z. Lenarčič and P. Prelovšek, Phys. Rev. B 90, 235136 (2014).
  * (25) A. Cavalleri, C. Tóth, C. W. Siders, J. A. Squier, F. Ráksi, P. Forget, and J. C. Kieffer, Phys. Rev. Lett. 87, 237401 (2001).
  * (26) E. Collet, M. H. Lemée-Cailleau, M. Le Cointe, H. Cailleau, M. Wulff, T. Luty, S. Koshihara, M. Meyer, L. Toupet, P. Rabiller, and S. Techert, Science 300, 612 (2003).
  * (27) T. Li, A. Patz, L. Mouchliadis, J. Yan, T. A. Lograsso, I. E. Perakis, and J. Wang, Nature 496, 69 (2013).
  * (28) S. Wall, D. Brida, S. R. Clark, H. P. Ehrke, D. Jaksch, A. Ardavan, S. Bonora, H. Uemura, Y. Takahashi, T. Hasegawa, H. Okamoto, G. Cerullo, and A. Cavalleri Nature Phys. 7, 114 (2011).
  * (29) J. H. Mentink, J. Phys.: Condens. Matter 29, 453001 (2017)
  * (30) H. Matsueda and S. Ishihara, J. Phys. Soc. Jpn. 76, 083703 (2007).
  * (31) Y. Kanamori, H. Matsueda, and S. Ishihara, PRL 107, 167403 (2011)
  * (32) Y. Kanamori, H. Matsueda, and S. Ishihara, Phys. Rev. Lett. 103, 267401 (2009).
  * (33) Y. Kanamori, H. Matsueda, and S. Ishihara, Phys. Rev. B 82, 115101 (2010).
  * (34) K. Miyano, T. Tanaka, Y. Tomioka, and Y. Tokura, Phys. Rev. Lett. 78, 4257 (1997).
  * (35) W. Koshibae, N. Furukawa, and N. Nagaosa, Phys. Rev. Lett. 103, 266402 (2009).
  * (36) E. Beaurepaire, J. C. Merle, A. Daunois, and J. Y. Bigot, Phys. Rev. Lett. 76, 4250 (1996).
  * (37) B. Koopmans, G. Malinowski, F. Dalla Longa, D. Steiauf, M. Fähnle, T. Roth, M. Cinchetti, and M. Aeschlimann, Nat. Mater. 9, 259 (2010).
  * (38) M. Fiebig, K. Miyano, Y. Tomioka, Y. Tokura, Science 280, 1925 (1998)
  * (39) R. D. Averitt, A. I. Lobad, C. Kwon, S. A. Trugman, V. K. Thorsmølle, and A. J. Taylor, Phys. Rev. Lett. 87, 017401 (2001).
  * (40) A. Ono, S. Ishihara, Phys. Rev. B 98, 214408 (2018).
  * (41) J. Li and M. Eckstein, Phys. Rev. B 103, 045133 (2021).
  * (42) H. Keiter and J. C. Kimball, Phys. Rev. Lett. 25, 672 (1970).
  * (43) M. Eckstein and P. Werner, Phys. Rev. B 82, 115115 (2010).
  * (44) S.-S. B. Lee, J. v. Delft, and A. Weichselbaum, Phys. Rev. B 96, 245106 (2017).
  * (45) P. Werner, M. Eckstein, M. Müller, and G. Refael, Nat. Comm. 10, 5556 (2019).
  * (46) P. Werner, J. Li, D. Golez, and M. Eckstein, Phys. Rev. B 100, 155130 (2019).
  * (47) M. Schüler, D. Golež, Y. Murakami, N. Bittner, A. Herrmann, H. U. R. Strand, P. Werner, and M. Eckstein, Comput. Phys. Commun. 257, 107484 (2020).