# Analytical derivation and expansion of the anti-Kibble–Zurek scaling in the
transverse field Ising model

Kaito Iwamura  Department of Physics, Waseda University, Tokyo 169-8555, Japan
Takayuki Suzuki  National Institute of Information and Communications
Technology, Nukui-Kitamachi 4-2-1, Koganei, Tokyo 184-8795, Japan

###### Abstract

A defect density which quantifies the deviation from the spin ground state
characterizes non-equilibrium dynamics during phase transitions. The widely
recognized Kibble-Zurek scaling predicts how the defect density evolves during
phase transitions. However, it can be perturbed by noise, leading to anti-
Kibble-Zurek scaling. In this research, we analytically investigate the effect
of Gaussian white noise on the transition probabilities of the Landau–Zener
model. We apply this model to the one-dimensional transverse field Ising model
and derive an analytical approximation for the defect density. Our analysis
reveals that under small noise conditions, the model follows an anti-
Kibble–Zurek scaling. As the noise increases, a new scaling behavior emerges,
showing higher accuracy than previously reported. Furthermore, we identify the
parameters that optimize the defect density based on the new scaling. This
allows for the refinement of optimized parameters with greater precision and
provides further validations of previously established scaling.

##  I Introduction

In modern quantum physics, understanding the dynamics of phase transitions is
essential for quantum information and quantum technology. Notably, the Quantum
Kibble–Zurek mechanism (QKZM) has attracted considerable attention Zurek  _et
al._ (2005); Damski (2005); Dziarmaga (2005); Polkovnikov (2005). QKZM serves
as a theoretical framework that provides scaling laws for quantum phase
transitions in the adiabatic limit. QKZM is an application of the Kibble–Zurek
mechanism to quantum systems, which was originally proposed within cosmology
Kibble (1976, 1980) and later applied to condensed matter systems Zurek (1985,
1993, 1996).

The transverse field XY-chain model and the transverse field Ising model are
systems where the spins are excited by an external magnetic field. The QKZM
suggests that quantum phase transitions follow particular scaling laws in the
adiabatic regime Zurek  _et al._ (2005); Damski (2005); Dziarmaga (2005);
Polkovnikov (2005); Mukherjee  _et al._ (2007). It allows us to predict the
defect density, which represents deviations from the ground state.

In the context of quantum control, it is desirable to minimize defects.
According to the QKZM law, as the driving speed $v$ decreases, the number of
defects decreases as well. The limit as $v$ approaches zero is referred to as
the adiabatic limit. However, recent study have observed an interesting
phenomenon called the anti-Kibble–Zurek mechanism (anti-KZM) Griffin  _et al._
(2012); Gao  _et al._ (2017); Dutta  _et al._ (2016); Singh and Gangadharaiah
(2021); Puebla  _et al._ (2020), where the scaling breaks down in the
adiabatic region when external noise is introduced into the system. Anti-KZM
may be effective in some experiments deriving scaling with the transverse
Ising chain Polkovnikov  _et al._ (2011); Cui  _et al._ (2016); Ai  _et al._
(2021); Wang  _et al._ (2014); Gong  _et al._ (2016) or with other systems Ulm
_et al._ (2013); Clark  _et al._ (2016); Chen  _et al._ (2011); Pyka  _et al._
(2013); Del Campo  _et al._ (2013); Braun  _et al._ (2015); Anquez  _et al._
(2016). The scaling law proposed by anti-KZM argues that defects become larger
when the driving speed $v$ is small, contrary to the QKZM. This suggests that
there exists a specific non-zero optimized driving speed, denoted by
$v_{\mathrm{opt}}$, which minimizes defects. Therefore, clarifying the scaling
based on anti-KZM is important and it is necessary to conduct theoretical
research.

In previous research simulating the defect density in the transverse field XY-
chain model with a Gaussian white noise Gao  _et al._ (2017), the Hamiltonian
has been considered

$\displaystyle
H_{XY}(t)=-\sum^{N}_{j=1}\quantity(J_{x}\sigma^{x}_{j}\sigma^{x}_{j+1}+J_{y}\sigma^{y}_{j}\sigma^{y}_{j+1}+h\sigma^{z}_{j}),$
(1) $\displaystyle=$
$\displaystyle-\frac{K}{2}\sum^{N}_{j=1}\quantity((1+\phi)\sigma^{x}_{j}\sigma^{x}_{j+1}+(1-\phi)\sigma^{y}_{j}\sigma^{y}_{j+1})-h\sum^{N}_{j=1}\sigma^{z}_{j}.$
(2)

Here, $K=J_{x}+J_{y}$, $\phi=\frac{J_{x}-J_{y}}{K}$. The number of spins $N$
is assumed to be sufficiently large, implying a thermodynamic limit. The
parameter $t$ takes from an infinite past to an infinite future. The system
evolves from the initial state, which is the ground state. At the final time,
the defect density which is the deviation from the ground state is calculated
numerically. The time-dependence is introduced through either
$h(t)=vt+\gamma(t)$, $J_{x}(t)=vt+\gamma(t)$, or $\phi(t)=vt+\gamma(t)$. Here,
$v$ is the driving speed, and $\gamma(t)$ denotes Gaussian noise with zero
mean and the second moment that satisfies $\overline{\gamma(t)}=0,\
\overline{\gamma(t)\gamma(t^{\prime})}=W^{2}\delta(t-t^{\prime})$. The scaling
of the defect density $n$ for small $W$ is represented as Gao  _et al._ (2017)

$\displaystyle n\sim cv^{\beta}+rv^{-1}.$ (3)

Here, $c$, $r$, and $\beta$ are coefficients determined by the choice of the
system and the parameters. The first term represents the well-known term
identified in QKZM, while the second term which is proper to $v^{-1}$ emerges
due to the introduction of noise. The mechanism is referred to as anti-KZM,
where $n$ increases as $v$ decreases. This is a numerically established result
and are not analytically demonstrated.

The defect density is derived by transforming the system into a time-dependent
quantum two-level system through the Jordan–Wigner transformation. Models such
as the one-dimensional transverse field XY-chain model and the one-dimensional
transverse field Ising model, which incorporate noise, are transformed into
noise-induced Landau–Zener models. When introducing the Gaussian white noise
$\gamma(t)$, it is established that the density matrix follows a master
equation Dutta  _et al._ (2016); Pichler  _et al._ (2013); Rahmani (2015). The
density matrix derives the transition probability and the defect density.

In a numerical research on the noise-induced transverse field Ising model with
Gaussian white noise Dutta  _et al._ (2016), corresponding to the case of
$J_{x}(t)=vt+\gamma(t),\ J_{y}=0,\ h(t)=J-vt-\gamma(t),\ t\in[0,J/v]$, the
similar anti-KZM scaling as in Eq. (3) has been confirmed. Unfortunately, the
scaling of $v^{-1}$ in Eq. (3) leads to divergence for small $v$. However, as
stated in Dutta  _et al._ (2016), applying the Kayanuma formula Kayanuma
(1984) provides an asymptotic expression in the small $v$ regime

$\displaystyle n^{\text{Kayanuma}}\sim\frac{1}{2}-\frac{\sqrt{v}}{4\pi J},$
(4)

which does not diverge. This formula is derived under the approximation that
the noise is sufficiently large. This result converges asymptotically to 1/2
as $v$ approaches zero, which is consistent with the numerical simulation.
However, in the other region this formula is not valid. Therefore, for small
$v$, a more precise approximate solution is required.

The analytical derivation of the transition probability in the noise-induced
Landau–Zener model has been done in previous studies Kayanuma (1984); Kenmoe
_et al._ (2013). However, the noise considered in these analyses is a
different types of noise, expressed as $\tilde{\gamma}(t)$. Here,
$\tilde{\gamma}(t)$ represents Gaussian noise with
$\overline{\tilde{\gamma}(t)}=0$ and
$\overline{\tilde{\gamma}(t)\tilde{\gamma}(t^{\prime})}=\tilde{W}^{2}e^{-\Gamma|t-t^{\prime}|}$.
It may not be feasible to directly apply Gaussian white noise to derive these
approximate solutions. This is because the study Kayanuma (1984) assumes a
large noise in the last and the system considered in Kenmoe  _et al._ (2013)
does not cover all cases. The research Singh and Gangadharaiah (2021)
analytically derives the defect density introduced the noise
$\tilde{\gamma}(t)$ in the nearest-neighbor interaction term by using the
result from Kenmoe  _et al._ (2013). In addition, the work Puebla  _et al._
(2020) has established anti-KZM scaling in fully connected models. Other
studies have explored systems under various perturbations, such as a recycling
term Dóra  _et al._ (2019); Gulácsi and Dóra (2021), a time-dependent simple
harmonic oscillator Mukherjee and Dutta (2009); Mullen  _et al._ (1989);
Kayanuma and Mizumoto (2000); Suzuki and Iwamura (2023), and a dissipative
thermal bath Patane  _et al._ (2008, 2009); Yin  _et al._ (2014); Nalbach  _et
al._ (2015); Nigro  _et al._ (2019); Rossini and Vicari (2019, 2020).

This research aims to analytically derive the defect density and to reproduce
the anti-KZM in the one-dimensional transverse field Ising model with Gaussian
white noise in the transverse field term. In addition, we derive a more
precise scaling than the previously proposed anti-KZM scaling in previously
unexplored regions. To this end, we first perform an analytical derivation of
the transition probabilities from the ground state in the Landau–Zener model
with Gaussian white noise. Second, by applying approximate solutions of the
Landau–Zener model to the transverse field Ising model, we find the anti-KZM
scaling. Finally, we examine the correlation between the optimized driving
speed $v_{\mathrm{opt}}$ with respect to the parameter $W^{2}$, which
represents the magnitude of the noise. While the scaling of $v_{\mathrm{opt}}$
with $W^{2}$ is already known to be $v_{\mathrm{opt}}\propto W^{\frac{4}{3}}$
in the small $W^{2}$ regime from numerical calculations Gao  _et al._ (2017);
Dutta  _et al._ (2016), our research extends this analysis to the unexplored
regime of large $W^{2}$.

In Sec. II, we calculate the transition probability of the Landau–Zener model
with Gaussian white noise. In Sec. III, we apply the formula discussed in Sec.
II to the transverse field ising chain with Gaussian white noise in the
adiabatic limit to determine the scaling of the defect density. Section IV is
dedicated to a discussion on the parameters that minimize the defect density.
Section V is the conclusion.

##  II The Landau–Zener model with Gaussian white noise

In the natural units, consider the following two-level Hamiltonian system

$\displaystyle H_{\text{LZ+noise}}(t)=H_{0}(t)+H_{\text{noise}}(t),$ (5)
$\displaystyle H_{0}(t)=\frac{vt}{2}\sigma^{z}+{J}\sigma^{x},\quad
H_{\text{noise}}(t)=\gamma(t)\sigma^{z}.$ (6)

Here, $\gamma(t)$ represents Gaussian white noise satisfying
$\overline{\gamma(t)}=0$ and
$\overline{\gamma(t)\gamma(t^{\prime})}=W^{2}\delta(t-t^{\prime})$, where
$W^{2}$ and ${J}$ are constants with an energy dimension. The quantities with
an overline indicate that they are ensemble-averaged. The driving speed $v$
has dimensions of the square of the energy. The term $H_{0}(t)$ corresponds to
a Landau–Zener-type Hamiltonian. It has already been well analyzed, and the
time-evolution operator is obtained with the parabolic cylinder functions.

The expectation value of a time-independent and noise-independent observable
$P$ at an arbitrary time $t$ can be expressed as

$\displaystyle\expectationvalue{\overline{P(t)}}=\Tr[\overline{\rho(t)}P]=\Tr[\overline{U(t,t_{i})\rho(t_{i})U^{\dagger}(t,t_{i})}P],$
(7)

where the density matrix $\rho(t)$ and the time-evolution operator
$U(t,t_{i})$ include contributions of noise. $t_{i}$ is an initial time and
$U(t_{i},t_{i})=1$. The noise-averaged density matrix $\overline{\rho(t)}$
follows the master equation

$\displaystyle\frac{d}{dt}\overline{\rho(t)}=-i\quantity[H_{0}(t),\overline{\rho(t)}]+\frac{1}{2}W^{2}\quantity[\quantity[\sigma^{z},\overline{\rho(t)}],\sigma^{z}].$
(8)

This equation is verified in Dutta  _et al._ (2016) by Novikov’s theorem
Novikov (1965) and by the other derivations Pichler  _et al._ (2013); Rahmani
(2015). Employing the dimensionless parameters: the time parameter
$\tau=\sqrt{v}t$, the adiabatic parameter $\kappa={J}^{2}/v$, and the noise
strength parameter $\lambda=W^{2}/J$, the master equation is expressed as

$\displaystyle\frac{d}{d\tau}\overline{\rho(\tau)}$
$\displaystyle=-i\quantity[H_{0}(\tau),\overline{\rho(\tau)}]+\frac{\lambda}{2}\sqrt{\kappa}\quantity[\quantity[\sigma^{z},\overline{\rho(\tau)}],\sigma^{z}],$
(9) $\displaystyle H_{0}(\tau)$
$\displaystyle=\frac{\tau}{2}\sigma^{z}+\sqrt{\kappa}\sigma^{x}.$ (10)

At an initial time $\tau_{i}$ chosen to be the infinite past, we set the
initial condition as the ground state
$\rho(\tau_{i})=\ket{\uparrow}\bra{\uparrow}$. Here,
$\sigma_{z}\ket{\uparrow}=\ket{\uparrow}$. To determine the probability of
maintaining $\ket{\uparrow}\bra{\uparrow}$, which corresponds to the
transition probability, we define $P=\ket{\uparrow}\bra{\uparrow}$ in Eq. (7).
We solve the time evolution equation for $\overline{\rho(\tau)}$ up to the
final time $\tau_{f}$ by the master equation Eq. (9) and calculate the
transition probability $\expectationvalue{\overline{P(\tau_{f})}}$. The first
term on the right-hand side of Eq. (9) represents the unitary term. If only
this term is considered, analytical solutions for the density matrix
$\overline{\rho(\tau)}$ at any time are obtained with the parabolic cylinder
functions. It is known that
$\expectationvalue{\overline{P(\tau_{f})}}=e^{-2\pi\kappa}$. As the parameter
$\kappa$ increases, the transition probability decreases. The second term in
Eq. (9), the dissipation term, accounts for the interaction with noise. The
presence of this term leads to an increase in the transition probability as
$\kappa$ increases. The master equation Eq. (9) exhibits explicit time
dependence solely through $H_{0}$. If the initial time is the infinite past
and the final time is the infinite future, the shift of the time $\tau$ of
$H_{0}$ does not depend on the transition probability at the final time.

In this study, we consider the case where the noise is sufficiently small,
i.e., $\lambda\ll 1$. We divide $\kappa$ into three distinct regions: (i) the
small region ($\kappa\ll 1$), (ii) the medium region
($1\ll\kappa\ll\frac{1}{\lambda}$), and (iii) the large region
($\frac{1}{\lambda}\ll\kappa$). In the regions (i) and (ii), we employ
perturbation approximation up to the first order with respect to $\lambda$ in
Sec. II.1. On the other hand, in the regions (ii) and (iii), we employ the
adiabatic approximation in Sec. II.2. By connecting the approximations, we
derive the globally effective approximation in Sec. II.3.

Figure 1: A schematic diagram illustrating the domains of $\kappa$: (i), (ii),
and (iii). Over each region, the corresponding approximation methods are
indicated.

###  II.1 First-order perturbation for noise

In this subsection, we consider (i) the small region ($\kappa\ll 1$) and (ii)
the medium region ($1\ll\kappa\ll 1/\lambda$). In these regions, we
analytically approximate the solution up to the first order with respect to
$\lambda$.

The Hamiltonian $H_{0}(\tau)$ is referred to as the Landau–Zener model, and
the time-evolution operator for $H_{0}(\tau)$ is expressed as

$\displaystyle U_{0}(\tau,\tau_{i})=$
$\displaystyle\begin{pmatrix}f(\tau,\tau_{i})&-g^{\ast}(\tau,\tau_{i})\\\
g(\tau,\tau_{i})&f^{\ast}(\tau,\tau_{i})\end{pmatrix},$ (11) $\displaystyle
f\left(\tau,\tau_{i}\right)=$ $\displaystyle
e^{-\pi\kappa/2}D_{i\kappa}\left(e^{-\frac{\pi}{4}i}\tau_{i}\right)D_{-i\kappa}\left(e^{\frac{\pi}{4}i}\tau\right)$
(12)
$\displaystyle{}+e^{-\pi\kappa/2}\sqrt{\kappa}D_{-i\kappa-1}\left(e^{\frac{\pi}{4}i}\tau_{i}\right)\sqrt{\kappa}D_{i\kappa-1}\left(e^{-\frac{\pi}{4}i}\tau\right),$
(13) $\displaystyle g\left(\tau,\tau_{i}\right)=$ $\displaystyle
e^{-\pi\kappa/2}e^{\frac{\pi}{4}i}D_{i\kappa}\left(e^{-\frac{\pi}{4}i}\tau_{i}\right)\sqrt{\kappa}D_{-i\kappa-1}\left(e^{\frac{\pi}{4}i}\tau\right)$
(14)
$\displaystyle{}-e^{-\pi\kappa/2}e^{\frac{\pi}{4}i}\sqrt{\kappa}D_{-i\kappa-1}\left(e^{\frac{\pi}{4}i}\tau_{i}\right)D_{i\kappa}\left(e^{-\frac{\pi}{4}i}\tau\right).$
(15)

We set the initial condition as $U(\tau_{i},\tau_{i})=1$. Moving to the
interaction picture

$\displaystyle\tilde{\rho}(\tau)=$ $\displaystyle
U^{\dagger}_{0}(\tau,\tau_{i})\rho(\tau)U_{0}(\tau,\tau_{i}),$ (16)
$\displaystyle\tilde{\sigma}_{z}(\tau,\tau_{i})=$ $\displaystyle
U^{\dagger}_{0}(\tau,\tau_{i})\sigma_{z}U_{0}(\tau,\tau_{i}),$ (17)

the transition probability $\expectationvalue{\overline{P(\tau)}}$ and the
time-evolution equation for $\overline{\tilde{\rho}(\tau)}$ are given by

$\displaystyle\expectationvalue{\overline{P(\tau)}}=\Tr\quantity[U_{0}(\tau,\tau_{i})\overline{\tilde{\rho}(\tau)}U_{0}^{\dagger}(\tau,\tau_{i})P],$
(18)
$\displaystyle\frac{d}{d\tau}\overline{\tilde{\rho}(\tau)}=\frac{1}{2}\lambda\sqrt{\kappa}\quantity[\tilde{\sigma}_{z}(\tau,\tau_{i}),\quantity[\overline{\tilde{\rho}(\tau)},\tilde{\sigma}_{z}(\tau,\tau_{i})]].$
(19)

By integrating Eq. (19) with respect to $\tau$, the transition probability up
to the first order in $\lambda$ is given by

$\displaystyle\expectationvalue{\overline{P(\tau_{f})}}=$
$\displaystyle\bra{\uparrow}U_{0}(\tau_{f},\tau_{i})\overline{\tilde{\rho}(\tau_{f})}U_{0}^{\dagger}(\tau_{f},\tau_{i})\ket{\uparrow}$
(20) $\displaystyle=$
$\displaystyle\bra{\uparrow}U_{0}(\tau_{f},\tau_{i})\rho(\tau_{i})U_{0}^{\dagger}(\tau_{f},\tau_{i})\ket{\uparrow}+\int^{\tau_{f}}_{\tau_{i}}d\tau\bra{\uparrow}U_{0}(\tau_{f},\tau_{i})\frac{\lambda\sqrt{\kappa}}{2}[\tilde{\sigma}_{z}(\tau,\tau_{i}),[\rho(\tau_{i}),\tilde{\sigma}_{z}(\tau,\tau_{i})]]U_{0}^{\dagger}(\tau_{f},\tau_{i})\ket{\uparrow}+\mathcal{O}\quantity(\lambda^{2}).$
(21)

Substituting the specific form of $U_{0}$ and taking the limits as
$\tau_{i}\to-\infty$ and $\tau_{f}\to\infty$, the expression can be
represented as

$\displaystyle\expectationvalue{\overline{P(\tau_{f})}}=e^{-2\pi\kappa}+$
$\displaystyle 4\lambda\sqrt{\kappa}\int^{\infty}_{-\infty}d\tau\
\quantity(|X(\tau)|^{2}+\frac{2e^{-2\pi\kappa}}{1-e^{-2\pi\kappa}}\quantity(\operatorname{Re}(X\quantity(\tau)))^{2}+\frac{2e^{-\pi\kappa}}{1-e^{-2\pi\kappa}}\operatorname{Re}\quantity(X\quantity(\tau))Y(\tau))+\mathcal{O}\quantity(\lambda^{2}),$
(22) $\displaystyle X(\tau)=$ $\displaystyle\kappa
e^{-\frac{\pi\kappa}{2}}D_{-i\kappa-1}(e^{\frac{i\pi}{4}}\tau)D_{i\kappa-1}(-e^{-\frac{i\pi}{4}}\tau),$
(23) $\displaystyle Y(\tau)=$
$\displaystyle\frac{\kappa}{2}e^{-\frac{\pi\kappa}{2}}\quantity(|D_{-i\kappa-1}(e^{\frac{i\pi}{4}}\tau)|^{2}+|D_{i\kappa-1}(-e^{-\frac{i\pi}{4}}\tau)|^{2}).$
(24)

The formulas used in the derivation are given in Appendix A. Both $|X(\tau)|$
and $|Y(\tau)|$ do not exceed 1. We derive $\overline{\tilde{\rho}(\tau)}$ by
evolving using Eq. (19) from the initial time $\tau_{i}$, and obtain
$\expectationvalue{\overline{P(\tau_{f})}}$, which includes contributions up
to the infinite order of $\lambda$. Eq. (19) does not have an analytical
solution and is therefore computed numerically. On the other hand, the
approximate solution up to the first order of $\lambda$ is derived by
calculating Eq. (22). While the analytical solutions for transition
probability at arbitrary times are not derived, this subsection focuses on
deriving an approximate solution at the final time. These two numerical plots
are illustrated in Fig. 2. According to Fig. 2, the transition probability
approaches the results of numerical simulations and first-order approximations
when $\kappa$ is small. As $\kappa$ increases beyond this point, the accuracy
of the approximation deteriorates due to the failure of the condition
$\lambda\kappa\ll 1$.

Figure 2: Comparison between numerical results obtained from solving the time
evolution through the master equation Eq. (19) and the first-order
approximation of $\lambda$ Eq. (22). The solid lines (Numerical) represent the
numerical results. The dashed lines (first order) represent the first-order
approximation. Both simulations cover a time interval from $\tau_{i}=-200$ to
$\tau_{f}=200$, with $\lambda$ set to $10^{-3}$ and $\kappa$ set to 100, 10,
and 1.

We seek the approximate solution for Eq. (22) in the regions (i) and (ii).
First, we begin the region (i): small $\kappa$. When $\kappa=0$, both
$X(\tau)$ and $Y(\tau)$ equal zero, leading to the dominance of the first term
of the right-hand side in Eq. (22). Consequently, we can approximate the
expression as

$\displaystyle\expectationvalue{\overline{P(\tau_{f})}}\stackrel{{\scriptstyle\mathrm{(i)}}}{{\sim}}e^{-2\pi\kappa}.$
(25)

This is strongly suppressed where $\kappa\gg 1$, and it only holds
significance in regions where $\kappa$ is small.

Second, consider the region (ii): medium $\kappa$. Given the condition
$\kappa\gg 1$, dropping terms proportional to $e^{-\pi\kappa}$ is an effective
approximation. Under this approximation, only the term with
$\quantity|X(\tau)|^{2}$ in Eq. (22) contributes. Thus, the transition
probability at the final time can be approximated as

$\displaystyle\expectationvalue{\overline{P(\tau_{f})}}\sim$ $\displaystyle
4\lambda\sqrt{\kappa}e^{-\pi\kappa}\kappa^{2}$ (26)
$\displaystyle\times\int^{\infty}_{-\infty}d\tau\
\quantity|D_{-i\kappa-1}(e^{\frac{i\pi}{4}}\tau)D_{-i\kappa-1}(-e^{\frac{i\pi}{4}}\tau)|^{2}.$
(27)

Ignoring the term of $\mathcal{O}(\kappa^{-2})$ and
$\mathcal{O}(e^{-\pi\kappa})$, $\expectationvalue{\overline{P(\tau_{f})}}$ is
approximated as

$\displaystyle\expectationvalue{\overline{P(\tau_{f})}}\stackrel{{\scriptstyle\mathrm{(ii)}}}{{\sim}}2\pi\lambda\kappa.$
(28)

The detailed derivation is conducted in Appendix B and C. It is evident that
Eq. (28) becomes zero when $\kappa=0$.

Therefore, effective solutions in the region where $\kappa$ is small are
dominant only within the region, and conversely, effective solutions in the
region where $\kappa$ is large are dominant only within the region. Hence, the
transition probability in both regions (i) and (ii) can be approximated as the
sum of effective solutions,

$\displaystyle\expectationvalue{\overline{P(\tau_{f})}}\stackrel{{\scriptstyle\mathrm{(i)},\mathrm{(ii)}}}{{\sim}}e^{-2\pi\kappa}+2\pi\lambda\kappa.$
(29)

The first term on the right-hand side represents the transition probability of
the Landau–Zener model in the absence of noise. The second term corresponds to
the influence of the noise. Despite being a transition probability, the right-
hand side diverges as $\kappa\to\infty$ due to the noise effect. This
approximate solution is valid only in the region where $\lambda\kappa\ll 1$,
and it is not allowed to take the limit of $\kappa$ approaching infinity. This
is confirmed by the comparison between numerical and approximate solutions in
Fig. 2. As $\kappa$ increases, the influence of higher-order contributions of
$\lambda$ becomes significant, leading to a deterioration in the accuracy of
the first-order approximation.

###  II.2 Adiabatic approximation

In the regions (ii) and (iii), that is $\kappa\gg 1$, the adiabatic
approximation is employed as Gulácsi and Dóra (2021). The equations of the
components of the density matrix $\overline{\rho(\tau)}$, derived from the
master equation Eq. (9), are

$\displaystyle\overline{\rho(\tau)}=\begin{pmatrix}\rho_{11}(\tau)&\rho_{12}(\tau)\\\
\rho_{21}(\tau)&\rho_{22}(\tau)\end{pmatrix},$ (30)
$\displaystyle\frac{d}{d\tau}|\rho(\tau)\rangle\\!\rangle=\frac{d}{d\tau}\quantity(\rho_{11}(\tau)\
\rho_{12}(\tau)\ \rho_{21}(\tau)\ \rho_{22}(\tau))^{T}$ (31) $\displaystyle=$
$\displaystyle\begin{pmatrix}0&i\sqrt{\kappa}&-i\sqrt{\kappa}&0\\\
i\sqrt{\kappa}&-i\frac{\tau}{2}-\lambda\sqrt{\kappa}&0&-i\sqrt{\kappa}\\\
-i\sqrt{\kappa}&0&i\frac{\tau}{2}-\lambda\sqrt{\kappa}&i\sqrt{\kappa}\\\
0&-i\sqrt{\kappa}&i\sqrt{\kappa}&0\end{pmatrix}\begin{pmatrix}\rho_{11}(\tau)\\\
\rho_{12}(\tau)\\\ \rho_{21}(\tau)\\\ \rho_{22}(\tau)\end{pmatrix}.$ (32)

The term $\mathcal{O}(\lambda^{2})$ is neglected in the following analysis in
this subsection. Let $z(\tau)=\frac{\tau}{4\sqrt{\kappa}}$. The time evolution
with this Liouvillian is calculated by using the adiabatic approximation
Sarandy and Lidar (2005). The specific calculations are detailed in Appendix
D. The eigenvalues of the Liouvillian are

$\displaystyle\chi_{1}(\tau)=$ $\displaystyle 0,$ (33)
$\displaystyle\chi_{2}(\tau)=$
$\displaystyle-\frac{\lambda}{z^{2}(\tau)+1}\sqrt{\kappa},$ (34)
$\displaystyle\chi_{3}(\tau)=$
$\displaystyle-\frac{2z^{2}(\tau)+1}{2z^{2}(\tau)+2}\lambda\sqrt{\kappa}-2i\sqrt{z^{2}(\tau)+1}\sqrt{\kappa},$
(35) $\displaystyle\chi_{4}(\tau)=$
$\displaystyle-\frac{2z^{2}(\tau)+1}{2z^{2}(\tau)+2}\lambda\sqrt{\kappa}+2i\sqrt{z^{2}(\tau)+1}\sqrt{\kappa}.$
(36)

The Liouvillian is diagonalized with $4\times 4$ matrix $S(\tau)$ as

$\displaystyle S(\tau)\begin{pmatrix}\chi_{1}(\tau)&0&0&0\\\
0&\chi_{2}(\tau)&0&0\\\ 0&0&\chi_{3}(\tau)&0\\\
0&0&0&\chi_{4}(\tau)\end{pmatrix}S^{-1}(\tau).$ (37)

The adiabatic basis is defined as
$|\tilde{\rho}(\tau)\rangle\\!\rangle=S^{-1}(\tau)|\rho(\tau)\rangle\\!\rangle$.
In this adiabatic basis, the master equation takes

$\displaystyle\frac{d}{d\tau}|\tilde{\rho}(\tau)\rangle\\!\rangle\sim$
$\displaystyle\begin{pmatrix}\chi_{1}(\tau)&0&0&0\\\ 0&\chi_{2}(\tau)&0&0\\\
0&0&\chi_{3}(\tau)&0\\\
0&0&0&\chi_{4}(\tau)\end{pmatrix}|\tilde{\rho}(\tau)\rangle\\!\rangle.$ (38)

The approximation is referred adiabatic approximation and is effective in the
region where $\kappa\gg 1$. Using the initial conditions
$\rho_{11}(\tau_{i})=1,\
\rho_{12}(\tau_{i})=\rho_{21}(\tau_{i})=\rho_{22}(\tau_{i})=0$ and the
conditions $\tau_{i}\to-\infty$ and $\tau_{f}\to\infty$, transition
probability is approximated as

$\displaystyle\expectationvalue{\overline{P(\tau_{f})}}\stackrel{{\scriptstyle\mathrm{(ii)},\mathrm{(iii)}}}{{\sim}}\frac{1}{2}\quantity(1-\exp\quantity(-4\pi\lambda\kappa)).$
(39)

This is the approximate solution when $\kappa\gg 1$.

###  II.3 Approximate solution for the transition probability

Connecting the approximation Eq. (29) derived in the regions (i) and (ii) with
the approximation Eq. (39) derived in the regions (ii) and (iii), the
effective transition probability in the regions (i), (ii), and (iii) is
expressed as

$\displaystyle\expectationvalue{\overline{P(\tau_{f})}}\stackrel{{\scriptstyle\text{all}}}{{\sim}}p^{\text{non-
ad}}(\kappa)+p^{\text{ad}}(\lambda,\kappa),$ (40) $\displaystyle p^{\text{non-
ad}}(\kappa)=e^{-2\pi\kappa},$ (41) $\displaystyle
p^{\text{ad}}(\lambda,\kappa)=\frac{1}{2}\quantity(1-\exp\quantity(-4\pi\lambda\kappa)).$
(42)

This coincides with Eq. (29) to the first-order perturbation of $\lambda$.
Fig. 3 displays the numerical result obtained from solving the master equation
Eq. (19) and the approximation solutions Eq. (40).

Figure 3: Plots of the transition probabilities at $\tau_{f}$ against
$\kappa$. The blue line (Numerical) represents the numerical simulation of Eq.
(19) over the time span from $\tau_{i}=-200$ to $\tau_{f}=200$. The yellow
dashed line (Non-adiabatic) represents $p^{\text{non-ad}}$, while the red
dashed line (Adiabatic) represents $p^{\text{ad}}$, with $\lambda=10^{-3}$.
The gray dot line represents $1/2$.

It becomes apparent that the non-adiabatic approximation solution is effective
in the small $\kappa$ region, while the adiabatic approximation solution is
effective in the large $\kappa$ region. We illustrate the comparison between
the numerically derived transition probability and $p^{\text{non-
ad}}+p^{\text{ad}}$ in the $\kappa\sim 1$ domain in Fig. 4.

Figure 4: Plots of the transition probabilities at $\tau_{f}$ against
$\kappa$, focusing around $\kappa=1$. The blue line (Numerical) represents the
numerical simulation of Eq. (19) over the time interval from $\tau_{i}=-200$
to $\tau_{f}=200$. The green dashed line (Non-adiabatic+Adiabatic) represents
$p^{\text{non-ad}}+p^{\text{ad}}$, with $\lambda=10^{-3}$.

While we have not found a precise solution near $\kappa\sim 1$, it is
confirmed that Eq. (40) does not significantly deviate from the numerical
results. The behavior of the transition probability in the regime where
$\kappa$ is sufficiently large is provided in Appendix E.

##  III TRANSVERSE FIELD ISING CHAIN with Gaussian white noise

In Sec. II, We have considered a system involving the two-level Landau–Zener
model with Gaussian white noise to analyze transition probability. In this
section, we extend this theoretical framework to apply it to the transverse
field Ising model following the approach Dziarmaga (2005). The Hamiltonian is
described as

$\displaystyle
H_{\text{IC}}(t)=-\frac{1}{2}\sum^{N}_{j=1}\quantity(\quantity(\frac{vt}{2}+\gamma(t))\sigma^{z}_{j}+J\sigma^{x}_{j}\sigma^{x}_{j+1}).$
(43)

Here, a periodic boundary condition $\sigma^{a}_{N+1}=\sigma^{a}_{1}$ is
imposed. We consider a sufficiently large number of spins $N$, corresponding
to a thermodynamic limit. For simplicity, we assume $N$ is an even number. In
the adiabatic regime, where $\kappa$ is large, we seek the defect density and
verify its scaling. Taking the initial state as the ground state of this
Hamiltonian, the state evolves under this Hamiltonian $H_{\text{IC}}$. The
initial state at initial time $\tau_{i}\sim-\infty$ consists of spin-down
particles $\ket{\psi(\tau_{i})}=\ket{\downarrow}_{j}^{\otimes N}$, where
$\sigma^{z}_{j}\ket{\downarrow}_{j}=-\ket{\downarrow}_{j}$. The ground state
at the final time $\tau_{f}\sim\infty$ consists of spin-up particles
$\ket{\uparrow}_{j}^{\otimes N}$. The defect density corresponds to the
density of spin-down particles at the final time. The aim is to derive
analytically the defect density at the final time. Here, we introduce the
spinless fermion operator $\hat{c}_{j}$ and perform the Jordan-Wigner
transformation as

$\displaystyle\sigma^{z}_{j}=1-2\hat{c}^{\dagger}_{j}\hat{c}_{j},\quad\sigma^{x}_{j}=\quantity(\hat{c}^{\dagger}_{j}+\hat{c}_{j})\prod_{l<j}\quantity(-\sigma^{z}_{l}).$
(44)

The operators $\hat{c}_{j}$ satisfy the anticommutation relations
$\\{\hat{c}_{j},\hat{c}_{k}^{\dagger}\\}=\delta_{jk}$ and
$\\{\hat{c}_{j},\hat{c}_{k}\\}=\\{\hat{c}_{j}^{\dagger},\hat{c}_{k}^{\dagger}\\}=0$.
When operating $\prod_{j=1}^{N}\sigma^{z}_{j}$ on the initial state, the
eigenvalue is 1. It indicates $\hat{c}_{j}$ follows the anti-periodic boundary
conditions $\hat{c}_{N+1}=-\hat{c}_{1}$ Dziarmaga (2005). Therefore,
performing the Fourier transform of the operator $\hat{c}_{j}$, we obtain

$\displaystyle\hat{c}_{j}=\frac{1}{\sqrt{N}}e^{-i\frac{\pi}{4}}\sum_{q}e^{-iqj}\hat{c}_{q},$
(45)

where $q=\pm(2n-1)\pi/N,\ n\in{1,\cdots,N/2}$. Using the operators
$\hat{c}_{q}$, the Hamiltonian Eq. (43) can be expressed as

$\displaystyle H_{\text{IC}}(t)=\sum_{q>0}(\hat{c}^{\dagger}_{q}\
\hat{c}_{-q})\quantity(h_{q}(t)+h_{\text{noise}}(t))\matrixquantity(\hat{c}_{q}\\\
\hat{c}^{\dagger}_{-q}),$ (46)

where

$\displaystyle h_{q}(t)=\quantity(\frac{vt}{2}+J\cos q)\sigma^{z}+{J}\sin
q\sigma^{x},$ (47) $\displaystyle\quad
h_{\text{noise}}(t)=\gamma(t)\sigma^{z}.$ (48)

We define the states as
$\ket{\uparrow}_{q}=\hat{c}_{q}^{\dagger}\hat{c}_{-q}^{\dagger}\ket{0}_{q}$
and $\ket{\downarrow}_{q}=\ket{0}_{q}$, where $\ket{0}_{q}$ represents the
normalized vacuum state for mode $q$, satisfying
$\hat{c}_{q}\ket{0}_{q}=\hat{c}_{-q}\ket{0}_{q}=0$. The ground state at
infinite past for mode $q$ is obtained by $\ket{\uparrow}_{q}^{\otimes{q>0}}$
from Eq. (46). Since the Hamiltonian is separable into sums for each mode $q$
and the total sum of $q$ remains conserved, the density operator at any time
$t$ can be decomposed as

$\displaystyle\hat{\rho}(t)=\bigotimes_{q>0}\hat{\rho}_{q}(t).$ (49)

Here, $\hat{\rho}_{q}$ is a density operator constructed only from
$\ket{\uparrow}_{q}$ and $\ket{\downarrow}_{q}$. For the initial conditions in
mode $q$, we select the pure state of the ground state
$\ket{\uparrow}_{q}\bra{\uparrow}$. The density operator $\hat{\rho}_{q}(t)$
is represented in matrix form as

$\displaystyle\rho_{q}(t)=\begin{pmatrix}{}_{q}\bra{\uparrow}\hat{\rho}_{q}(t)\ket{\uparrow}_{q}&{}_{q}\bra{\uparrow}\hat{\rho}_{q}(t)\ket{\downarrow}_{q}\\\
{}_{q}\bra{\downarrow}\hat{\rho}_{q}(t)\ket{\uparrow}_{q}&{}_{q}\bra{\downarrow}\hat{\rho}_{q}(t)\ket{\downarrow}_{q}\end{pmatrix}.$
(50)

The density matrix $\rho_{q}(t)$ evolves under the Hamiltonian
$h_{q}(t)+h_{\text{noise}}(t)$. The transition probability in mode $q$ at time
$t$ is expressed using the time evolution operator $U_{q}(t,t_{i})$ and the
noise-averaged density matrix $\overline{\rho_{q}(t)}$ as

$\displaystyle\expectationvalue{\overline{P_{q}(t)}}=$
$\displaystyle\Tr\quantity[\overline{\rho_{q}(t)}P_{q}]$ (51) $\displaystyle=$
$\displaystyle\Tr[\overline{U_{q}(t,t_{i})\rho_{q}(t_{i})U_{q}^{\dagger}(t,t_{i})}P_{q}],$
(52)

where $P_{q}=\text{diag}(1,0)$, which corresponds to
$\ket{\uparrow}_{q}\bra{\uparrow}$. Similar to Sec. II, the noise-averaged
density matrix $\overline{\rho_{q}(t)}$ satisfies the master equation

$\displaystyle\frac{d}{dt}\overline{\rho_{q}(t)}=-i\quantity[h_{q}(t),\overline{\rho_{q}(t)}]+\frac{1}{2}W^{2}\quantity[\quantity[\sigma^{z},\overline{\rho_{q}(t)}],\sigma^{z}].$
(53)

The differences with Eq. (5) involve a time shift in $h_{q}(t)$ compared to
$H_{0}(t)$, and the modification of the coupling term from $J$ in $H_{0}(t)$
to $J\sin q$ in $h_{q}(t)$. Considering the initial and final times extending
to the infinite past and future, the time shift in $h_{q}(t)$ can be
neglected. The dimensionless master equation is expressed as

$\displaystyle\frac{d}{d\tau}\overline{\rho_{q}(\tau)}$
$\displaystyle=-i\quantity[\tilde{h}_{q}(\tau),\overline{\rho_{q}(\tau)}]+\frac{\lambda}{2}\sqrt{\kappa}\quantity[\quantity[\sigma^{z},\overline{\rho_{q}(\tau)}],\sigma^{z}],$
(54) $\displaystyle\tilde{h}_{q}(\tau)$
$\displaystyle=\frac{\tau}{2}\sigma^{z}+\sqrt{\kappa_{q}}\sigma^{x},$ (55)
$\displaystyle\kappa_{q}$ $\displaystyle=\kappa\sin^{2}q.$ (56)

This master equation is similar to Eq. (9). The parameter $\kappa_{q}$
satisfies the conditions $\kappa_{\frac{\pi}{2}}=\kappa$ and
$\kappa_{q}\leq\kappa=J^{2}/v$. In the vicinity of $q\sim 0,\ \pi$, the energy
gap diminishes, indicating the non-adiabatic regime. Whereas elsewhere, a
significant energy gap prevails, indicating the adiabatic domain.

Subsequently, we consider the adiabatic regime. Using the transition
probability of $q$ mode, we calculate the defect density. As the state
$\ket{\uparrow}_{q}$ denotes a two-particle state, the expected defects number
is twice the transition probability. Consequently, the expected defects number
at the final time $\tau_{f}$ can be expressed as

$\displaystyle\expectationvalue{\mathcal{N}(\tau_{f})}=2\sum_{q>0}\quantity|\expectationvalue{\overline{P_{q}(\tau_{f})}}|^{2}.$
(57)

In the thermodynamic limit, the normalized expectation value of defects can be
approximated as

$\displaystyle n(\tau_{f})$
$\displaystyle=\frac{\expectationvalue{\mathcal{N}(\tau_{f})}}{N}\to\int^{\pi}_{0}\frac{dq}{\pi}\quantity|\expectationvalue{\overline{P_{q}(\tau_{f})}}|^{2},$
(58)

where $\sum_{q>0}\to\int^{\pi}_{0}\frac{dq}{2\pi}$ is used. By using the
result of the Landau-Zener model, we obtain the transition probability at
final time $\tau_{f}\sim\infty$ as

$\displaystyle\quantity|\expectationvalue{\overline{P_{q}(\tau_{f})}}|^{2}$
$\displaystyle\sim l^{\text{non-ad}}(q)+l^{\text{ad}}(q),$ (59) $\displaystyle
l^{\text{non-ad}}(q)$ $\displaystyle=e^{-2\pi\kappa_{q}},$ (60) $\displaystyle
l^{\text{ad}}(q)$
$\displaystyle=\frac{1}{2}\quantity(1-\exp\quantity(-4\pi\lambda\sqrt{\kappa}\sqrt{\kappa_{q}})).$
(61)

We present plots showing both the numerical results and the approximate
solution Eq. (59) in Fig. 5. The Numerical result converges to
$l^{\mathrm{non-ad}}(q)$ when $q$ is close to 0 or $\pi$, while it converges
to $l^{ad}(q)$ elsewhere. Therefore, Eq. (59) provides an approximate solution
and its integration yields an approximation of the defect.

Figure 5: Plots of the transition probability
$\quantity|\expectationvalue{\overline{P_{q}(\tau_{f})}}|^{2}$ against $q$.
The blue line (Numerical) represents the numerical simulation of Eq. (19) over
the time span from $\tau_{i}=-200$ to $\tau_{f}=200$. The yellow dashed line
(Non-adiabatic) represents $l^{\text{non-ad}}(q)$, while the red dashed line
(Adiabatic) represents $l^{\text{ad}}(q)$, with $\lambda=10^{-3}$.

Using Eq. (59), the defect density $n(\lambda,\kappa)$ can be approximated as

$\displaystyle n^{\infty\text{-order}}(\lambda,\kappa)=$ $\displaystyle
n^{\text{KZM}}(\kappa)+n^{\text{noise}}(\lambda,\kappa),$ (62) $\displaystyle
n^{\text{KZM}}(\kappa)$
$\displaystyle=\int^{\pi}_{0}\frac{dq}{\pi}l^{\text{non-
ad}}(q)=\int^{\pi}_{0}\frac{dq}{\pi}e^{-2\pi\kappa\sin^{2}q}$ (63)
$\displaystyle\sim 2\int^{\infty}_{0}\frac{dq}{\pi}e^{-2\pi\kappa
q^{2}}=\frac{1}{\pi\sqrt{2\kappa}},$ (64) $\displaystyle
n^{\text{noise}}(\lambda,\kappa)$
$\displaystyle=\int^{\pi}_{0}\frac{dq}{\pi}l^{\text{ad}}(q)$ (65)
$\displaystyle=\int^{\pi}_{0}\frac{dq}{\pi}\quantity(\frac{1}{2}-\frac{1}{2}\exp\quantity(-4\pi\lambda\kappa|\sin
q|))$ (66)
$\displaystyle=\frac{1}{2}-\frac{1}{2}I_{0}\quantity(4\pi\lambda\kappa)+\frac{1}{2}L_{0}\quantity(4\pi\lambda\kappa)$
(67)

Here, the condition $\kappa\gg 1$ is used in the calculation of Eq. (64).
Consequently, the approximate solution $n^{\infty\text{-order}}$ is valid only
in the regime $\kappa\gg 1$. $I_{0}(z)$ is the modified Bessel function of the
first kind and $L_{0}(z)$ is the modified Struve function Gradshteyn and
Ryzhik (2014). $n^{\text{KZM}}(\kappa)$ represents the defect density in the
absence of noise and corresponds with the result from the QKZM.
$n^{\text{noise}}$ denotes the contribution due to the influence of Gaussian
white noise. $n^{\infty\text{-order}}$ includes contributions up to infinite
order with respect to $\lambda$. By using the asymptotic conditions of
$I_{0}(z)$ and $L_{0}(z)$, $n^{\infty\text{-order}}\to 1/2$ can be verified
when $\lambda\kappa\to\infty$. This implies that the system asymptotically
approaches the infinite-temperature state as the noise becomes large. In
contrast, when $\lambda\kappa\ll 1$, using the properties
$I_{0}(z)=1+\frac{z^{2}}{4}+\mathcal{O}(z^{4})$ and
$L_{0}(z)=\frac{2z}{\pi}+\mathcal{O}(z^{3})$, we can further refine the
expression

$\displaystyle
n^{\infty\text{-order}}(\lambda,\kappa)\sim\frac{1}{\pi\sqrt{2\kappa}}+4\lambda\kappa=:n^{\text{1st}}(\lambda,\kappa).$
(68)

The second term in $n^{\text{1st}}$ is proportional to $W^{2}\propto\lambda$,
which is expected from computational researches Gao  _et al._ (2017); Dutta
_et al._ (2016). Eq. (68) demonstrates the successful reproduction of the
anti-KZM scaling. In the absence of Gaussian white noise, a larger $\kappa$,
i.e., a smaller $v$, leads to a smaller defect density as expected by QKZM.
However, when Gaussian white noise is induced, the defect density increases in
the region where $\kappa$ is large. The second term diverges as
$\kappa\to\infty$, which is due to neglecting higher-order contributions of
$\lambda\kappa$. Incorporating these higher-order $\lambda$ contributions
improves the accuracy of the approximation. The approximate solution including
contributions up to the second order with respect to $\lambda$ is given by

$\displaystyle n^{\infty\text{-order}}(\lambda,\kappa)$
$\displaystyle\sim\frac{1}{\pi\sqrt{2\kappa}}+4\lambda\kappa-2\pi^{2}\lambda^{2}\kappa^{2}$
(69) $\displaystyle=:n^{\text{2nd}}(\lambda,\kappa).$ (70)

The results of the defect density $n$ obtained numerically, along with
$n^{\infty\text{-order}}$, $n^{\text{1st}}$, and $n^{\text{2nd}}$ are compared
in Fig. 6. The spin number $N$ is chosen to be a large, specifically $N=100$.
$\lambda$ is smaller in the upper figure.

Figure 6: Plots of the defect density against $\kappa$. The upper (a)
corresponds to $\lambda=10^{-3}$, while the lower (b) corresponds to
$\lambda=5\times 10^{-3}$. The blue solid lines (Numerical) represent the
numerical simulation of Eq. (19) over the time span $\tau_{i}=-200$ to
$\tau_{f}=200$, plotting the discrete sum in Eq. (58) with $N=100$. The green
dashed lines ($\infty$-order) represent $n^{\infty\text{-order}}$, the yellow
dashed lines (First-order) represent $n^{\text{1st}}$, and the red dot lines
(Second-order) represent $n^{\text{2nd}}$. The pink solid lines (KZ) represent
$n^{\text{KZM}}$.

It has been verified that $n^{\infty\text{-order}}$ closely matches the defect
density $n$ in the region where $\kappa$ is large, which is adiabatic region.
In the region where $\kappa$ is small, the approximation used to perform the
integration for deriving $n^{\text{KZM}}$ becomes less effective, and
$n^{\infty\text{-order}}$ is not a suitable approximation. Both
$n^{\text{1st}}$ and $n^{\text{2nd}}$ serve as approximate solutions to
$n^{\infty\text{-order}}$ when $\kappa$ is small. It is observed that the
region where the approximation is effective for $n^{\text{2nd}}$ is larger
than that for $n^{\text{1st}}$. Hence, we have successfully obtained a new
scaling $n^{\text{2nd}}$. As $\lambda$ increases, the $\kappa$ which
represents the location of the minimum becomes smaller and the validity of
$n^{\text{1st}}$ in its neighborhood becomes worse. The behavior of the defect
density in the regime where $\kappa$ is sufficiently large is provided in
Appendix E.

##  IV The derivation of the optimization parameter

Let $v_{\mathrm{opt}}$ denote the value at which the defect density $n$
attains its minimum. The previous studies Gao  _et al._ (2017); Dutta  _et
al._ (2016) numerically seek the minimum of the defect density $n$ and obtain
the scaling of $v_{\mathrm{opt}}$ with respect to $W^{2}\propto\lambda$. We
seek the approximation of $v_{\mathrm{opt}}$ through the critical point. As
defined by
$\kappa^{\text{1st}}_{\min}=(2^{7}\pi^{2}\lambda^{2})^{-\frac{1}{3}}$, the
derivative is obtained by

$\displaystyle\left.\frac{d}{d\kappa}n^{\text{1st}}\right|_{\kappa=\kappa^{\text{1st}}_{\min}}=-\frac{1}{\pi(2\kappa^{\text{1st}}_{\min})^{\frac{3}{2}}}+4\lambda=0.$
(71)

The critical point of $n^{\text{1st}}$ is $\kappa=\kappa^{\text{1st}}_{\min}$,
where $n^{\text{1st}}$ attains its minimum value. Since $\kappa=J^{2}/v$,
$v_{\mathrm{opt}}$ for $n^{\text{1st}}$ is given by

$\displaystyle
v^{\text{1st}}_{\text{opt}}/J^{2}=(2^{7}\pi^{2}\lambda^{2})^{\frac{1}{3}}.$
(72)

Therefore, similar to the previous study Dutta  _et al._ (2016), we have
demonstrated $v_{\mathrm{opt}}^{\text{1st}}\propto\lambda^{\frac{2}{3}}\propto
W^{\frac{4}{3}}$. Next, we choose $n^{\text{2nd}}$ as an approximate solution
for the defect density $n$. The critical point $\kappa^{\text{2nd}}$ cannot be
determined analytically. We make an assumption that $\kappa^{\text{2nd}}$ is
close to $\kappa^{\text{1st}}$. We determine $\kappa^{\text{2nd}}$ as the
perturbation from $\kappa^{\text{1st}}_{\min}$. Setting
$\kappa^{\text{2nd}}_{\min}=\kappa^{\text{1st}}_{\min}\quantity(1+\varepsilon)$,
the derivative is expressed as

$\displaystyle\left.\frac{d}{d\kappa}n^{\text{2nd}}\right|_{\kappa=\kappa^{\text{2nd}}_{\min}}$
(73) $\displaystyle=$
$\displaystyle-\frac{1}{\pi(2\kappa^{\text{2nd}}_{\min})^{\frac{3}{2}}}+4\lambda-4\pi^{2}\lambda^{2}\kappa^{\text{2nd}}_{\min}$
(74) $\displaystyle=$
$\displaystyle-\frac{1}{\pi(2\kappa^{\text{1st}}_{\min})^{\frac{3}{2}}}(1+\varepsilon)^{-\frac{3}{2}}+4\lambda-4\pi^{2}\lambda^{2}\kappa^{\text{1st}}_{\min}(1+\varepsilon).$
(75)

Assuming $\varepsilon$ to be sufficiently small, incorporating perturbations
up to the first order in $\varepsilon$ yields

$\displaystyle\left.\frac{d}{d\kappa}n^{\text{2nd}}\right|_{\kappa=\kappa^{\text{2nd}}_{\min}}\sim\frac{3\varepsilon}{2\pi(2\kappa^{\text{1st}}_{\min})^{\frac{3}{2}}}-4\pi^{2}\lambda^{2}\kappa^{\text{1st}}_{\min}(1+\varepsilon).$
(76)

Here, Eq. (71) is utilized. For the right-hand side to be zero, $\varepsilon$
is represented as

$\displaystyle\varepsilon=\frac{1}{3}\quantity(\frac{\pi^{4}\lambda}{2^{4}})^{\frac{1}{3}}+\mathcal{O}\quantity(\lambda^{\frac{2}{3}}).$
(77)

Thus, ignoring $\mathcal{O}\quantity(\lambda^{\frac{4}{3}})$,
$v_{\mathrm{opt}}$ for $n^{\text{2nd}}$ is expressed as

$\displaystyle
v^{\text{2nd}}_{\mathrm{opt}}/J^{2}=(2^{7}\pi^{2}\lambda^{2})^{\frac{1}{3}}-\frac{2}{3}\pi^{2}\lambda.$
(78)

If selected $n^{\infty\text{-order}}$ as an approximate solution for the
defect density, $v_{\mathrm{opt}}$ is obtained by finding numerically the
minimum value as $v$ varies. This is denoted as
$v_{\mathrm{opt}}^{\infty\text{-order}}$. By solving Eq. (19) and calculating
the discrete sum in Eq. (58), the minimum values are obtained numerically for
$v_{\mathrm{opt}}$. $v_{\mathrm{opt}}$,
$v_{\mathrm{opt}}^{\infty\text{-order}}$, $v^{\text{1st}}_{\mathrm{opt}}$, and
$v^{\text{2nd}}_{\mathrm{opt}}$ against $\lambda$ are plotted respectively in
Fig. 7.

Figure 7: Plots of $v_{\mathrm{opt}}$ against $\lambda$. The blue line
(Numerical) represents the numerical simulation of Eq. (19) over the time span
$\tau_{i}=-200$ to $\tau_{f}=200$, plotting the discrete sum in Eq. (58), with
$N=100$. The green dashed line ($\infty$-order) represents
$v^{\infty\text{-order}}_{\mathrm{opt}}$, the yellow dashed line (First-order)
represents $v^{\text{1st}}_{\mathrm{opt}}$, and the red dot line (Second-
order) represents $v^{\text{2nd}}_{\mathrm{opt}}$.

It can be confirmed that $v_{\mathrm{opt}}^{\infty\text{-order}}$ serves as an
approximation of $v_{\mathrm{opt}}$. Since $v^{\text{2nd}}_{\mathrm{opt}}$
approaches the numerical results more closely than
$v^{\text{1st}}_{\mathrm{opt}}$, we suggest that
$v^{\text{2nd}}_{\mathrm{opt}}$ is more effective. The relative deviation of
$v^{\text{2nd}}_{\mathrm{opt}}$ over $v^{\text{1st}}_{\mathrm{opt}}$ is
expressed as

$\displaystyle\frac{v^{\text{2nd}}_{\text{opt}}}{v^{\text{1st}}_{\text{opt}}}\sim
1-\left(\frac{\pi^{4}\lambda}{2^{4}\times
3^{3}}\right)^{\frac{1}{3}}=:1-\zeta(\lambda).$ (79)

Let the defect density in each optimized $v$ be $n_{\text{opt}}$, which
corresponds to the minimum value of $n$. The relative deviation is expressed
as

$\displaystyle\frac{n^{\text{2nd}}_{\text{opt}}}{n^{\text{1st}}_{\text{opt}}}\sim
1-\quantity(\frac{\pi^{4}\lambda}{2^{10}\times
3^{3}})^{\frac{1}{3}}=:1-\xi(\lambda).$ (80)

The magnitude of $\zeta$ and $\xi$ indicates the effectiveness of the
approximation of $v^{\text{1st}}_{\text{opt}}$ and
$n^{\text{1st}}_{\text{opt}}$ respectively. If both parameters are
sufficiently small compared to 1, $v^{\text{1st}}_{\text{opt}}$ and
$n^{\text{1st}}_{\text{opt}}$ can be considered effective. For example, when
$\lambda=10^{-3}$, $\zeta\approx 0.061$ and $\xi\approx 0.015$. When
$\lambda=5\times 10^{-3}$, $\zeta\approx 0.10$ and $\xi\approx 0.026$. We see
that as $\lambda$ increases, $\zeta$ and $\xi$ increase. This leads to a loss
of validity of the first-order approximation in the regime of large $\lambda$.
However, if $\lambda\gtrsim 0.01$, approximation formulas break down. This is
because, as $\lambda$ increases, the critical point of $\kappa$ decreases and
the approximation utilized to perform the integration in Eq. (64) becomes less
effective near the critical point.

##  V Conclusion

Under the approximation of a small noise, $\lambda\ll 1$, we analytically
derived the transition probability of the Landau–Zener model with Gaussian
white noise. We considered three $\kappa$ regions: (i) small, (ii) medium, and
(iii) large. For the regions (i) and (ii), we solved a master equation under a
first-order approximation for $\lambda$. This involved performing integrals of
products of special functions, and the feasibility of these integrals depends
strongly on the system. Hence, it is anticipated that this method may not be
applied to other systems. However, by selecting parameters that asymptotically
approach those of the system we addressed, we suggested that other systems may
converge to the solutions obtained in this research. For the regions (ii) and
(iii), we employed an adiabatic approximation method, which is considered to
be useful in various other systems.

Applying the analysis of the Landau–Zener model, we analytically derived the
defect density in the transverse field Ising model with Gaussian white noise.
The results reveal that in the regime where $\kappa\ll 1/\lambda$, the anti-
KZM scaling previously known through numerical computations Dutta  _et al._
(2016); Gao  _et al._ (2017) is effective. Furthermore, the optimal parameter
$v_{\mathrm{opt}}$ follows the same scaling as the anti-KZM previously known.
On the other hand, in the region where $1/\lambda\ll\kappa$, a new scaling was
found. It can be argued that to minimize the defect density, it is desirable
to choose $v_{\mathrm{opt}}$ smaller than the value obtained from anti-KZM
previously known in the noisy region. Furthermore, we have achieved the
capability to discuss the efficacy of first-order perturbation.

While the effectiveness of the approximate solutions for $\kappa\sim 1$ has
not been formally proven, they are promising approximations by compared to
numerical results. Therefore, further investigation is necessary to confirm
the reason. Additionally, we aim at extending the methodology employed in this
study to other systems such as XY chain models.

As an application to another system, we consider the Schwinger mechanism
Heisenberg and Euler (1936); Weisskopf (1936); Schwinger (1951), where
electron-positron pairs are expected to be generated when an electric field is
applied to a vacuum. Since the time evolution of this system is related to the
Landau–Zener system, we propose applying the methods of this research to
derive the number density of electron-positron pairs generated by an electric
field, accounting for the influence of noise. It is anticipated that, due to
the induced noise, the number of generated electrons significantly increases
as the dynamically assited Schwinger mechanism Schützhold  _et al._ (2008).

## Acknowledgement

We thank H. Nakazato for helpful discussions. T.S. was supported by Japan’s
MEXT Quantum Leap Flagship Program Grant No. JPMXS0120319794.

##  Appendix A The derivation of the transition probability with perturbation

We derive the first-order approximation solution Eq. (22) for the
perturbation. $f(\tau_{f},\tau_{i})$ and $g(\tau_{f},\tau_{i})$ are components
of the time evolution operator $U_{0}(\tau,\tau_{i})$. These are obtained
analytically and are expressed as linear combinations of parabolic cylinder
functions Gradshteyn and Ryzhik (2014) as Eq. (11). The transition probability
is expressed up to first order in $\lambda$ as

$\displaystyle\expectationvalue{\overline{P(\tau_{f})}}=$
$\displaystyle|f(\tau_{f},\tau_{i})|^{2}-4\lambda\sqrt{\kappa}(|f(\tau_{f},\tau_{i})|^{2}-|g(\tau_{f},\tau_{i})|^{2})\int^{\tau_{f}}_{\tau_{i}}d\tau\
|f(\tau,\tau_{i})|^{2}|g(\tau,\tau_{i})|^{2}$ (81)
$\displaystyle+4\lambda\sqrt{\kappa}\operatorname{Re}\quantity(f(\tau_{f},\tau_{i})g(\tau_{f},\tau_{i})\int^{\infty}_{-\infty}d\tau\
f^{*}(\tau,\tau_{i})g^{*}(\tau,\tau_{i})(|f(\tau,\tau_{i})|^{2}-|g(\tau,\tau_{i})|^{2})).$
(82)

Using the equations for a positive value of $\tau$

$\displaystyle D_{i\kappa}\quantity(e^{\frac{3\pi}{4}i}\tau)$
$\displaystyle\sim
e^{\frac{i}{4}\tau^{2}}e^{-\frac{3\pi}{4}\kappa}\tau^{i\kappa}+\mathcal{O}(\tau^{-1}),$
(83) $\displaystyle D_{-i\kappa}\quantity(e^{\frac{\pi}{4}i}\tau)$
$\displaystyle\sim
e^{-\frac{i}{4}\tau^{2}}e^{\frac{\pi}{4}\kappa}\tau^{-i\kappa}+\mathcal{O}(\tau^{-1}),$
(84) $\displaystyle D_{-i\kappa-1}\quantity(e^{\frac{5\pi}{4}i}\tau)$
$\displaystyle\sim\frac{\sqrt{2\pi}}{\Gamma(i\kappa+1)}e^{\frac{i}{4}\tau^{2}}e^{-\frac{\pi}{4}\kappa}\tau^{i\kappa}+\mathcal{O}(\tau^{-1}),$
(85) $\displaystyle D_{i\kappa-1}\quantity(e^{-\frac{\pi}{4}i}\tau)$
$\displaystyle\sim\mathcal{O}(\tau^{-1}),$ (86)

$|\Gamma(1-ix)|^{2}=\frac{\pi x}{\sinh(\pi x)}$ and limiting
$\tau_{i}\to-\infty$, $\tau_{f}\to\infty$, the transition probability is
obtained as

$\displaystyle\tilde{X}(\tau)=\frac{\sqrt{2\pi\kappa}}{\Gamma(1+i\kappa)}e^{-2\pi\kappa}D_{i\kappa}(e^{-\frac{i\pi}{4}}\tau)\sqrt{\kappa}D_{i\kappa-1}(e^{-\frac{i\pi}{4}}\tau),\quad\tilde{Y}(\tau)=e^{-\frac{\pi\kappa}{2}}\quantity(|D_{i\kappa}(e^{-i\frac{\pi}{4}}\tau)|^{2}-|\sqrt{\kappa}D_{i\kappa-1}(e^{-i\frac{\pi}{4}}\tau)|^{2}),$
(87)
$\displaystyle\expectationvalue{\overline{P(\tau_{f})}}=e^{-2\pi\kappa}+4\lambda\sqrt{\kappa}\int^{\infty}_{-\infty}d\tau\
\quantity(\frac{1-2e^{-2\pi\kappa}}{1-e^{-2\pi\kappa}}e^{2\pi\kappa}|\tilde{X}(\tau)|^{2}+\operatorname{Re}(\tilde{X}(\tau))\tilde{Y}(\tau)).$
(88)

Using the formula
$\frac{\sqrt{2\pi}}{\Gamma(\nu+1)}D_{\nu}(z)=e^{-\frac{i\pi\nu}{2}}D_{-\nu-1}(-iz)+e^{\frac{i\pi\nu}{2}}D_{-\nu-1}(iz)$,
we obtain Eq. (22).

##  Appendix B The derivation of approximation with perturbation

We demonstrate Eq. (28) in the region (ii) ($1\ll\kappa\ll\frac{1}{\lambda}$).
The approximate solution for the transition probability is approximated as

$\displaystyle\expectationvalue{\overline{P(\tau_{f})}}\sim$ $\displaystyle
4\lambda\sqrt{\kappa}e^{-\pi\kappa}\kappa^{2}$ (89)
$\displaystyle\times\int^{\infty}_{-\infty}d\tau\
\quantity|D_{-i\kappa-1}(e^{\frac{i\pi}{4}}\tau)D_{-i\kappa-1}(-e^{\frac{i\pi}{4}}\tau)|^{2}.$
(90)

To perform this calculation, we utilize the following formula Nasri (2016)

$\displaystyle D_{-\nu}(z)D_{-\nu}(-z)$ (91) $\displaystyle=$
$\displaystyle\frac{2e^{\frac{z^{2}}{2}}}{\Gamma\quantity(\nu+\frac{1}{2})}\int^{\infty}_{0}dy\
\cos(2zy)e^{-2y^{2}}{}_{1}F_{1}\quantity(\frac{1}{2},\nu+\frac{1}{2},y^{2}).$
(92)

Here ${}_{1}F_{1}$ denotes the confluent hypergeometric function. By
substituting $z=e^{i\frac{\pi}{4}}\tau$ and $\nu=i\kappa+1$, and changing the
integration path to $t=e^{i\frac{\pi}{4}}y$, we get

$\displaystyle
D_{-i\kappa-1}(e^{i\frac{\pi}{4}}\tau)D_{-i\kappa-1}(-e^{i\frac{\pi}{4}}\tau)$
(93) $\displaystyle=$
$\displaystyle\frac{2e^{-i\frac{\pi}{4}}e^{i\frac{\tau^{2}}{2}}}{\Gamma\quantity(i\kappa+\frac{3}{2})}\int^{\infty}_{0}dt\
\cos(2\tau
t)e^{2it^{2}}{}_{1}F_{1}\quantity(\frac{1}{2},i\kappa+\frac{3}{2},-it^{2}).$
(94)

Therefore, $\expectationvalue{\overline{P(\tau_{f})}}$ is approximated as

$\displaystyle\expectationvalue{\overline{P(\tau_{f})}}\sim$ $\displaystyle
4\lambda\sqrt{\kappa}e^{-\pi\kappa}\kappa^{2}\int^{\infty}_{-\infty}d\tau\
\quantity|D_{-i\kappa-1}(e^{\frac{i\pi}{4}}\tau)D_{-i\kappa-1}(-e^{\frac{i\pi}{4}}\tau)|^{2}$
(95) $\displaystyle=$ $\displaystyle
4\lambda\sqrt{\kappa}e^{-\pi\kappa}\kappa^{2}\int^{\infty}_{-\infty}d\tau\
\frac{1}{\quantity|\Gamma\quantity(i\kappa+\frac{3}{2})|^{2}}\int^{\infty}_{0}dr\
(e^{2i\tau r}+e^{-2i\tau
r})e^{2ir^{2}}{}_{1}F_{1}\quantity(\frac{1}{2},i\kappa+\frac{3}{2},-ir^{2})$
(96) $\displaystyle\times\int^{\infty}_{0}ds\ (e^{2i\tau s}+e^{-2i\tau
s})e^{-2is^{2}}{}_{1}F_{1}\quantity(\frac{1}{2},-i\kappa+\frac{3}{2},is^{2}).$
(97)

We change the order of integration, performing the integration over $\tau$
first. Utilizing the conditions that $r$ and $s$ are positive, we employ the
formulas $\int^{\infty}_{-\infty}d\tau\ e^{2i\tau(r-s)}=\pi\delta(r-s)$ and
$\int^{\infty}_{-\infty}d\tau\ e^{2i\tau(r+s)}=0$. Consequently, we get

$\displaystyle\expectationvalue{\overline{P(\tau_{f})}}\sim$ $\displaystyle
4\lambda\sqrt{\kappa}\frac{2\pi}{\quantity|\Gamma\quantity(i\kappa+\frac{3}{2})|^{2}}e^{-\pi\kappa}\kappa^{2}$
(98) $\displaystyle\times\int^{\infty}_{0}dr\
\quantity|{}_{1}F_{1}\quantity(\frac{1}{2},i\kappa+\frac{3}{2},-ir^{2})|^{2}.$
(99)

Note that this derivation is based on the assumption of infinite time
integration. Therefore, this relationship cannot be applied in systems where
the time span is finite. Furthermore, using the following formula

$\displaystyle\frac{2\pi}{\quantity|\Gamma\quantity(i\kappa+\frac{3}{2})|^{2}}e^{-\pi\kappa}\kappa^{2}=\frac{\kappa^{2}}{\kappa^{2}+\frac{1}{4}}\quantity(1+e^{-2\pi\kappa})$
(100)

and the approximations $e^{-2\pi\kappa}\ll 1$ and $\kappa^{2}\gg\frac{1}{4}$,
the transition probability is expressed as

$\displaystyle\expectationvalue{\overline{P(\tau_{f})}}\sim
4\lambda\sqrt{\kappa}\int^{\infty}_{0}dr\
\quantity|{}_{1}F_{1}\quantity(\frac{1}{2},i\kappa+\frac{3}{2},-ir^{2})|^{2}.$
(101)

By employing the variable transformation $r=\sqrt{\kappa}x$, the expression
can be rewritten as

$\displaystyle\expectationvalue{\overline{P(\tau_{f})}}\sim
4\lambda\kappa\int^{\infty}_{0}dx\
\quantity|{}_{1}F_{1}\quantity(\frac{1}{2},i\kappa+\frac{3}{2},-i\kappa
x^{2})|^{2}.$ (102)

Furthermore, when $\kappa$ is large,

$\displaystyle{}_{1}F_{1}\quantity(\frac{1}{2},i\kappa+\frac{3}{2},-i\kappa
x^{2})\sim\frac{1}{\sqrt{1+x^{2}}}.$ (103)

We confirm in Appendix C that the approximate formula Eq. (103) is effective
when $\kappa$ is large. Consequently, by utilizing

$\displaystyle\int^{\infty}_{0}dt\frac{1}{1+x^{2}}=\frac{\pi}{2},$ (104)

we can demonstrate the validity of Eq. (28).

##  Appendix C The approximation of the special function

As stated in Appendix B, the approximation formula Eq. (103) is shown to be
effective when $\kappa$ is large. By the properties of a confluent
hypergeometric function Gradshteyn and Ryzhik (2014)

$\displaystyle\frac{d}{dz}{}_{1}F_{1}\quantity(a,b,z)=\frac{a}{b}{}_{1}F_{1}\quantity(a+1,b+1,z)$
(105)

and

$\displaystyle{}_{1}F_{1}\quantity(a,b,z)=$
$\displaystyle\frac{b-a-z-1}{b}{}_{1}F_{1}\quantity(a+1,b+1,z)$ (106)
$\displaystyle+\frac{a+1}{b}{}_{1}F_{1}\quantity(a+2,b+1,z),$ (107)

we get

$\displaystyle\frac{d}{dx}\quantity(\sqrt{1+x}{}_{1}F_{1}\quantity(\frac{1}{2},i\kappa+\frac{3}{2},-i\kappa
x))$ (108) $\displaystyle=$
$\displaystyle\frac{3}{4\sqrt{1+x}(i\kappa+\frac{3}{2})}{}_{1}F_{1}\quantity(\frac{5}{2},i\kappa+\frac{5}{2},-i\kappa
x).$ (109)

Integrating both sides with respect to $x$ from 0 and using
${}_{1}F_{1}\quantity(a,b,0)=1$ yield

$\displaystyle{}_{1}F_{1}\quantity(\frac{1}{2},i\kappa+\frac{3}{2},-i\kappa
x)=\frac{1}{\sqrt{1+x}}\quantity(1+E_{\kappa}(x)),$ (110) $\displaystyle
E_{\kappa}(x)=\frac{3}{4(i\kappa+\frac{3}{2})}\int^{x}_{0}\frac{dy}{\sqrt{1+y}}{}_{1}F_{1}\quantity(\frac{5}{2},i\kappa+\frac{5}{2},-i\kappa
y).$ (111)

Therefore, the approximation Eq. (103) is valid in the regime $E(x)\ll 1$.
Fig. 8 plots the $x$-dependence of $\quantity|E_{\kappa}(x)|$. It can be
observed $\quantity|E_{\kappa}(x)|\ll 1$ for large $\kappa$, confirming the
validity of the approximation Eq. (103).

Figure 8: Plots of $\quantity|E_{\kappa}(x)|$ against $x$. As $\kappa$
increases, it is observed that this value becomes significantly smaller than 1
for any $x$.

##  Appendix D Adiabatic approximation

In the regions (ii) and (iii) ($1\ll\kappa$), the effectiveness of the
adiabatic approximation is demonstrated. Let
$|\rho(\tau)\rangle\\!\rangle=\quantity(\rho_{11}(\tau)\ \rho_{12}(\tau)\
\rho_{21}(\tau)\ \rho_{22}(\tau))^{T}$. The master equation Eq. (9), is
expressed as

$\displaystyle|\dot{\rho}(\tau)\rangle\\!\rangle=\mathcal{L}(\tau)|\rho\rangle\\!\rangle.$
(112)

Here, $\mathcal{L}(\tau)$ denotes the Liouvillian at time $\tau$, which is
non-Hermitian. Diagonalizing $\mathcal{L}$ yields

$\displaystyle\tilde{\mathcal{L}}(\tau)=S^{-1}(\tau)\mathcal{L}(\tau)S(\tau)$
(113)

Here,
$\tilde{\mathcal{L}}(\tau)=\operatorname{diag}(\chi_{1}(\tau),\chi_{2}(\tau),\chi_{3}(\tau),\chi_{4}(\tau))$.
$|\chi_{\alpha}(\tau)\rangle\\!\rangle$ and
$\langle\\!\langle\hat{\chi}_{\alpha}(\tau)|$ are the normalized right
eigenvectors and the normalized left eigenvectors of $\mathcal{L}(\tau)$,
satisfying

$\displaystyle\mathcal{L}(\tau)|\chi_{\alpha}(\tau)\rangle\\!\rangle=$
$\displaystyle\chi_{\alpha}(\tau)|\chi_{\alpha}(\tau)\rangle\\!\rangle,$ (114)
$\displaystyle\langle\\!\langle\hat{\chi}_{\alpha}(\tau)|\mathcal{L}(\tau)=$
$\displaystyle\langle\\!\langle\hat{\chi}_{\alpha}(\tau)|\chi_{\alpha}(\tau),$
(115)

with $\alpha=1,2,3,4$. When there is no degeneracy in the eigenvalues of the
Liouvillian, meaning $\chi_{\alpha}\neq\chi_{\beta}$, the condition
$\langle\\!\langle\hat{\chi}_{\alpha}(\tau)|\chi_{\beta}(\tau)\rangle\\!\rangle=0$
is satisfied. Normalization is performed such that
$\langle\\!\langle\hat{\chi}_{\alpha}(\tau)|\chi_{\beta}(\tau)\rangle\\!\rangle=\delta_{\alpha\beta}$.
Expanding the state in terms of the right eigenvectors as

$\displaystyle|\rho(\tau)\rangle\\!\rangle=\sum_{\beta=1}^{4}c_{\beta}(\tau)|\chi_{\beta}(\tau)\rangle\\!\rangle,$
(116)

the master equation Eq. (112) takes

$\displaystyle\sum_{\beta=1}^{4}\quantity(\dot{c}_{\beta}(\tau)|\chi_{\beta}(\tau)\rangle\\!\rangle+c_{\beta}(\tau)|\dot{\chi}_{\beta}(\tau)\rangle\\!\rangle)$
(117) $\displaystyle=$
$\displaystyle\sum_{\beta=1}^{4}c_{\beta}(\tau)\chi_{\beta}(\tau)|\chi_{\beta}(\tau)\rangle\\!\rangle.$
(118)

Applying the left eigenvector $\langle\\!\langle\hat{\chi}_{\alpha}(\tau)|$
yields

$\displaystyle\dot{c}_{\alpha}(\tau)=\chi_{\alpha}(\tau)c_{\alpha}(\tau)-\sum_{\beta=1}^{4}c_{\alpha}(\tau)\langle\\!\langle\chi_{\alpha}(\tau)|\dot{\hat{\chi}}_{\beta}(\tau)\rangle\\!\rangle.$
(119)

Therefore, if
$\langle\\!\langle\hat{\chi}_{\alpha}(\tau)|\dot{\chi}_{\beta}(\tau)\rangle\\!\rangle$
is sufficiently small, the equation becomes closed for $\alpha$, resulting in
adiabatic time evolution.
$\langle\\!\langle\chi_{\alpha}(\tau)|\dot{\chi}_{\beta}(\tau)\rangle\\!\rangle$
has the dimension of frequency. The relevant physical frequency scale is given
by the energy gap
$\omega_{\alpha\beta}(\tau)=\chi_{\alpha}(\tau)-\chi_{\beta}(\tau)$.
Therefore, the adiabatic condition is given by Sarandy and Lidar (2005)

$\displaystyle l_{\alpha\beta}\ll r_{\alpha\beta},$ (120) $\displaystyle
l_{\alpha\beta}=\max_{\tau}\quantity|\langle\\!\langle\hat{\chi}_{\alpha}(\tau)|\dot{\chi}_{\beta}(\tau)\rangle\\!\rangle|,\
r_{\alpha\beta}=\min_{\tau}\quantity|\omega_{\alpha\beta}(\tau)|.$ (121)

Here, $r_{\alpha\beta}=r_{\beta\alpha}$. The Liouvillian under consideration
is given as

$\displaystyle\mathcal{L}(\tau)=\begin{pmatrix}0&i\sqrt{\kappa}&-i\sqrt{\kappa}&0\\\
i\sqrt{\kappa}&-i\frac{\tau}{2}-\lambda\sqrt{\kappa}&0&-i\sqrt{\kappa}\\\
-i\sqrt{\kappa}&0&i\frac{\tau}{2}-\lambda\sqrt{\kappa}&i\sqrt{\kappa}\\\
0&-i\sqrt{\kappa}&i\sqrt{\kappa}&0\end{pmatrix}.$ (122)

We calculate the eigenvalues and eigenvectors up to the first order of
perturbation with respect to $\lambda$. Let
$z(\tau)=\frac{\tau}{4\sqrt{\kappa}}$. The term $\mathcal{O}(\lambda^{2})$ is
neglected in the following analysis in this subsection. The eigenvalues and
right eigenvectors of the Liouvillian are

$\displaystyle\chi_{1}=$ $\displaystyle 0,$ (123) $\displaystyle\chi_{2}=$
$\displaystyle-\frac{\lambda}{z^{2}+1}\sqrt{\kappa},$ (124)
$\displaystyle\chi_{3}=$
$\displaystyle-\frac{2z^{2}+1}{2z^{2}+2}\lambda\sqrt{\kappa}-2i\sqrt{z^{2}+1}\sqrt{\kappa},$
(125) $\displaystyle\chi_{4}=$
$\displaystyle-\frac{2z^{2}+1}{2z^{2}+2}\lambda\sqrt{\kappa}+2i\sqrt{z^{2}+1}\sqrt{\kappa},$
(126)
$\displaystyle|\chi_{1}(\tau)\rangle\\!\rangle=\begin{pmatrix}\frac{1}{\sqrt{2}}\\\
0\\\ 0\\\ \frac{1}{\sqrt{2}}\end{pmatrix},$ (127)
$\displaystyle|\chi_{2}(\tau)\rangle\\!\rangle=\begin{pmatrix}\frac{z}{\sqrt{2(z^{2}+1)}}\\\
\frac{1}{\sqrt{2(z^{2}+1)}}+i\lambda\frac{z\sqrt{z^{2}+1}}{2\sqrt{2}(z^{2}+1)^{2}}\\\
\frac{1}{\sqrt{2(z^{2}+1)}}-i\lambda\frac{z\sqrt{z^{2}+1}}{2\sqrt{2}(z^{2}+1)^{2}}\\\
-\frac{z}{\sqrt{2(z^{2}+1)}}\end{pmatrix},$ (128)
$\displaystyle|\chi_{3}(\tau)\rangle\\!\rangle=\begin{pmatrix}\frac{1}{2\sqrt{z^{2}+1}}+i\lambda\frac{4z^{2}+1}{16(z^{2}+1)^{2}}\\\
-\frac{1}{2}\quantity(1+\frac{z}{\sqrt{z^{2}+1}})+i\lambda\frac{3z+\sqrt{z^{2}+1}}{16(z^{2}+1)^{2}}\\\
\frac{1}{2}\quantity(1-\frac{z}{\sqrt{z^{2}+1}})+i\lambda\frac{3z-\sqrt{z^{2}+1}}{16(z^{2}+1)^{2}}\\\
-\frac{1}{2\sqrt{z^{2}+1}}-i\lambda\frac{4z^{2}+1}{16(z^{2}+1)^{2}}\end{pmatrix},$
(129)
$\displaystyle|\chi_{4}(\tau)\rangle\\!\rangle=\begin{pmatrix}\frac{1}{2\sqrt{z^{2}+1}}-i\lambda\frac{4z^{2}+1}{16(z^{2}+1)^{2}}\\\
\frac{1}{2}\quantity(1-\frac{z}{\sqrt{z^{2}+1}})-i\lambda\frac{3z-\sqrt{z^{2}+1}}{16(z^{2}+1)^{2}}\\\
-\frac{1}{2}\quantity(1+\frac{z}{\sqrt{z^{2}+1}})-i\lambda\frac{3z+\sqrt{z^{2}+1}}{16(z^{2}+1)^{2}}\\\
-\frac{1}{2\sqrt{z^{2}+1}}+i\lambda\frac{4z^{2}+1}{16(z^{2}+1)^{2}}\end{pmatrix}.$
(130)

Here, the tight eigenvectors are chosen such that they satisfy
$l_{\alpha\beta}=l_{\beta\alpha}$ up to the first order of $\lambda$. The
eigenvector corresponding to the eigenvalue $\chi_{1}=0$ is independent of the
time $\tau$ and is considered to be the infinite-temperature state. This state
implies that both the probability of transitioning and the probability of not
transitioning are ${1}/{2}$. The $4\times 4$ matrix $S(\tau)$ is defined as

$\displaystyle S(\tau)=\quantity(|\chi_{1}(\tau)\rangle\\!\rangle\
|\chi_{2}(\tau)\rangle\\!\rangle\ |\chi_{3}(\tau)\rangle\\!\rangle\
|\chi_{4}(\tau)\rangle\\!\rangle).$ (131)

Based on this, the inverse matrix $S^{-1}(\tau)$ can also be obtained as

$\displaystyle S^{-1}(\tau)=\begin{pmatrix}\langle\\!\langle\chi_{1}(\tau)|\\\
\langle\\!\langle\chi_{2}(\tau)|\\\ \langle\\!\langle\chi_{3}(\tau)|\\\
\langle\\!\langle\chi_{4}(\tau)|\end{pmatrix}.$ (132)

First, We verify the adiabatic condition Eq. (120). From
$l_{12}(\tau)=l_{13}(\tau)=l_{14}(\tau)=0$, the state $\ket{\chi_{1}(\tau)}$
is adiabatic. By calculating other states, the following results are obtained.
Up to the first order of $\lambda$,
$l_{23}(\tau)=l_{24}(\tau)=\frac{1}{4\sqrt{2\kappa}}$ and
$r_{23}(\tau)=r_{24}(\tau)=2\sqrt{\kappa}$, as well as
$l_{34}(\tau)=\frac{\lambda}{50\sqrt{5\kappa}}$ and
$r_{34}(\tau)=4\sqrt{\kappa}$. Therefore, the adiabatic condition is
$\kappa\gg 1$.

Second, we calculate the transition probability $\tilde{\rho}(\tau_{f})_{11}$
under the adiabatic approximation condition. From Eq. (116), the adiabatic
basis is defined as

$\displaystyle|\tilde{\rho}(\tau)\rangle\\!\rangle=S^{-1}(\tau)|\rho(\tau)\rangle\\!\rangle=\begin{pmatrix}c_{1}(\tau)\
c_{2}(\tau)\ c_{3}(\tau)\ c_{4}(\tau)\end{pmatrix}^{T},$ (133)

Under the adiabatic approximation, the differential equation is expressed as

$\displaystyle\dot{c}_{\alpha}(\tau)\sim\chi_{\alpha}(\tau)c_{\alpha}(\tau).$
(134)

Therefore, $c_{\alpha}(\tau)$ is given by

$\displaystyle
c_{\alpha}(\tau)=c_{\alpha}(\tau_{i})\exp\quantity(\int^{\tau}_{\tau_{i}}d\tau^{\prime}\chi_{\alpha}(\tau^{\prime})),$
(135)

The initial condition is $|\rho(\tau_{i})\rangle\\!\rangle=(1\ 0\ 0\ 0)^{T}$.
In the limit $\tau_{i}\to-\infty$, the initial conditions in the adiabatic
basis are given by $c_{1}(\tau_{i})=\frac{1}{\sqrt{2}},\
c_{2}(\tau_{i})=-\frac{1}{\sqrt{2}},\ c_{3}(\tau_{i})=c_{4}(\tau_{i})=0$. When
evolved up to time $\tau$, the state can be approximated in the original basis
as

$\displaystyle\ket{\rho(\tau)}\sim$ $\displaystyle
S(\tau)\begin{pmatrix}\frac{1}{\sqrt{2}}\\\
-\frac{1}{\sqrt{2}}\exp\quantity(-\int^{\tau}_{\tau_{i}}d\tau^{\prime}\frac{\lambda\sqrt{\kappa}}{z^{2}(\tau^{\prime})+1})\\\
0\\\ 0\end{pmatrix}$ (136) $\displaystyle=$
$\displaystyle\frac{1}{2}\begin{pmatrix}1-\frac{z(\tau)}{\sqrt{z^{2}(\tau)+1}}\exp\quantity(-\int^{\tau}_{\tau_{i}}d\tau^{\prime}\frac{\lambda\sqrt{\kappa}}{z^{2}(\tau^{\prime})+1})\\\
-\quantity(\frac{1}{\sqrt{z^{2}(\tau)+1}}+i\lambda\frac{z(\tau)}{2(z^{2}(\tau)+1)^{\frac{3}{2}}})\exp\quantity(-\int^{\tau}_{\tau_{i}}d\tau^{\prime}\frac{\lambda\sqrt{\kappa}}{z^{2}(\tau^{\prime})+1})\\\
-\quantity(\frac{1}{\sqrt{z^{2}(\tau)+1}}-i\lambda\frac{z(\tau)}{2(z^{2}(\tau)+1)^{\frac{3}{2}}})\exp\quantity(-\int^{\tau}_{\tau_{i}}d\tau^{\prime}\frac{\lambda\sqrt{\kappa}}{z^{2}(\tau^{\prime})+1})\\\
1+\frac{z(\tau)}{\sqrt{z^{2}(\tau)+1}}\exp\quantity(-\int^{\tau}_{\tau_{i}}d\tau^{\prime}\frac{\lambda\sqrt{\kappa}}{z^{2}(\tau^{\prime})+1}).\end{pmatrix}$
(137)

The normalization and Hermiticity conditions of the density matrix are
verified. Thus, the state at the final time is determined by substituting
$\tau_{i}\sim-\infty$, $\tau_{f}\sim\infty$, resulting in

$\displaystyle\ket{\rho(\tau_{f})}\sim\begin{pmatrix}\frac{1}{2}\quantity(1-\exp\quantity(-4\pi\lambda\kappa))\\\
0\\\ 0\\\
\frac{1}{2}\quantity(1+\exp\quantity(-4\pi\lambda\kappa))\end{pmatrix}.$ (138)

Therefore, under the adiabatic approximation, $\rho(\tau_{f})_{11}$ is given
by

$\displaystyle\rho_{11}(\tau_{f})\sim\frac{1}{2}\left(1-\exp\left(-4\pi\lambda\kappa\right)\right).$
(139)

##  Appendix E The comparisons with the Kayanuma formula

We first compute the transition probability in the regime where $\kappa$ is
sufficiently large ($\lambda\kappa\gg 1$). Next, we verify the defect density
within the transverse-field Ising model under the condition of large $\kappa$.
Our approximate solutions of the transition probability and the defect density
are compared with those derived from the Kayanuma formula, which is effective
in large $\lambda$.

First, we verify the behavior of the transition probability in the regime
where $\kappa$ is sufficiently large. the approximate solution Eq. (40)
approaches $\frac{1}{2}$ as $\kappa\to\infty$, similar to the Kayanuma formula
Kayanuma (1984)

$\displaystyle\expectationvalue{\overline{P(\tau_{f})}}\sim
p^{\text{Kayanuma}}(\kappa)=\frac{1}{2}\quantity(1-\exp\quantity(-4\pi\kappa)).$
(140)

The comparison between $p^{\text{Kayanuma}}$ and $p^{\text{non-
ad}}+p^{\text{ad}}$ is shown in Fig. 9. In the regime where $\kappa$ is large,
all cases asymptotically approach $1/2$. However, since we are considering
regions where $\lambda$ is small, the $p^{\text{Kayanuma}}$ failed to
reproduce the behavior of the numerical results.

Figure 9: Plots of the transition probabilities at $\tau_{f}$ against
$\kappa$. The blue line (Numerical) represents the numerical simulation of Eq.
(19) over the time span from $\tau_{i}=-200$ to $\tau_{f}=200$. The green
dashed line (Non-adiabatic+Adiabatic) represents $p^{\text{non-
ad}}+p^{\text{ad}}$, while the red dashed line (Kayanuma formula) represents
$p^{\text{Kayanuma}}$, with $\lambda=10^{-3}$. The gray dashed line represents
$1/2$.

Next, We verify the defect density. The approximate solutions of Eq. (68) and
Eq. (70) break down in the region where $\kappa$ is large. In previous
research Dutta  _et al._ (2016), the following scaling for $n$ in the regime
of sufficiently large $\kappa$ was proposed.

$\displaystyle
n^{\text{Kayanuma}}(\kappa)\sim\frac{1}{2}-\frac{1}{4\pi\sqrt{\kappa}}$ (141)

On the other hand, by employing the asymptotic formulas of the functions
$I_{0}(z)$ and $L_{0}(z)$, the expression can be reduced to the following by
neglecting $\mathcal{O}(4\pi\lambda\kappa)^{-2}$

$\displaystyle n^{\text{reciprocal}}(\lambda,\kappa)\sim
n^{\text{KZM}}(\kappa)+\frac{1}{2}-\frac{1}{4\pi^{2}\lambda\kappa}.$ (142)

We were able to confirm through Fig. 10 that this scaling yields the closer
result to the numerical one than that of Kayanuma formula.

Figure 10: Plots of the defect density against $\kappa$, focusing around
$\kappa\gg 1$. The blue solid line (Numerical) represents the numerical
simulation of Eq. (19) over the time span $\tau_{i}=-200$ to $\tau_{f}=200$,
plotting the discrete sum in Eq. (58), with $N=100$. The green dashed line
($\infty$-order) represents $n^{\infty\text{-order}}$, the yellow dashed line
(Reciprocal) represents $n^{\text{reciprocal}}$, and the pink dashed line
(Kayanuma) represents $n^{\text{Kayanuma}}$, with $\lambda=5\times 10^{-3}$.
The gray dot line represents $1/2$.

## References

  * Zurek  _et al._ (2005) W. H. Zurek, U. Dorner, and P. Zoller, Phys. Rev. lett. 95, 105701 (2005).
  * Damski (2005) B. Damski, Phys. Rev. lett. 95, 035701 (2005).
  * Dziarmaga (2005) J. Dziarmaga, Phys. Rev. lett. 95, 245701 (2005).
  * Polkovnikov (2005) A. Polkovnikov, Phys. Rev. B 72, 161201 (2005).
  * Kibble (1976) T. W. Kibble, J. Phys. A: Math. Gen. 9, 1387 (1976).
  * Kibble (1980) T. W. Kibble, Phys. Rep. 67, 183 (1980).
  * Zurek (1985) W. H. Zurek, Nature 317, 505 (1985).
  * Zurek (1993) W. H. Zurek, Acta Phys. Pol. B 24, 1301 (1993).
  * Zurek (1996) W. H. Zurek, Phys. Rep. 276, 177 (1996).
  * Mukherjee  _et al._ (2007) V. Mukherjee, U. Divakaran, A. Dutta, and D. Sen, Phys. Rev. B 76, 174303 (2007).
  * Griffin  _et al._ (2012) S. M. Griffin, M. Lilienblum, K. T. Delaney, Y. Kumagai, M. Fiebig, and N. A. Spaldin, Phys. Rev. X 2, 041022 (2012).
  * Gao  _et al._ (2017) Z.-P. Gao, D.-W. Zhang, Y. Yu, and S.-L. Zhu, Phys. Rev. B 95, 224303 (2017).
  * Dutta  _et al._ (2016) A. Dutta, A. Rahmani, and A. Del Campo, Phys. rev. lett. 117, 080402 (2016).
  * Singh and Gangadharaiah (2021) M. Singh and S. Gangadharaiah, Phys. Rev. B 104, 064313 (2021).
  * Puebla  _et al._ (2020) R. Puebla, A. Smirne, S. F. Huelga, and M. B. Plenio, Phys. Rev. Lett. 124, 230602 (2020).
  * Polkovnikov  _et al._ (2011) A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Rev. Mod. Phys. 83, 863 (2011).
  * Cui  _et al._ (2016) J.-M. Cui, Y.-F. Huang, Z. Wang, D.-Y. Cao, J. Wang, W.-M. Lv, L. Luo, A. Del Campo, Y.-J. Han, C.-F. Li, _et al._ , Sci. rep. 6, 33381 (2016).
  * Ai  _et al._ (2021) M.-Z. Ai, J.-M. Cui, R. He, Z.-H. Qian, X.-X. Gao, Y.-F. Huang, C.-F. Li, and G.-C. Guo, Phys. Rev. A 103, 012608 (2021).
  * Wang  _et al._ (2014) L. Wang, C. Zhou, T. Tu, H.-W. Jiang, G.-P. Guo, and G.-C. Guo, Phys. Rev. A 89, 022337 (2014).
  * Gong  _et al._ (2016) M. Gong, X. Wen, G. Sun, D.-W. Zhang, D. Lan, Y. Zhou, Y. Fan, Y. Liu, X. Tan, H. Yu, _et al._ , Sci. Rep. 6, 22667 (2016).
  * Ulm  _et al._ (2013) S. Ulm, J. Roßnagel, G. Jacob, C. Degünther, S. Dawkins, U. Poschinger, R. Nigmatullin, A. Retzker, M. Plenio, F. Schmidt-Kaler, _et al._ , Nat. commun. 4, 2290 (2013).
  * Clark  _et al._ (2016) L. W. Clark, L. Feng, and C. Chin, Science 354, 606 (2016).
  * Chen  _et al._ (2011) D. Chen, M. White, C. Borries, and B. DeMarco, Phys. Rev. Lett. 106, 235304 (2011).
  * Pyka  _et al._ (2013) K. Pyka, J. Keller, H. Partner, R. Nigmatullin, T. Burgermeister, D. Meier, K. Kuhlmann, A. Retzker, M. B. Plenio, W. Zurek, _et al._ , Nat. commun. 4, 2291 (2013).
  * Del Campo  _et al._ (2013) A. Del Campo, T. Kibble, and W. Zurek, J. Phys. Condens. Matter 25, 404210 (2013).
  * Braun  _et al._ (2015) S. Braun, M. Friesdorf, S. S. Hodgman, M. Schreiber, J. P. Ronzheimer, A. Riera, M. Del Rey, I. Bloch, J. Eisert, and U. Schneider, Proc. Natl. Acad. Sci. USA 112, 3641 (2015).
  * Anquez  _et al._ (2016) M. Anquez, B. Robbins, H. Bharath, M. Boguslawski, T. Hoang, and M. Chapman, Phys. rev. lett. 116, 155301 (2016).
  * Pichler  _et al._ (2013) H. Pichler, J. Schachenmayer, A. J. Daley, and P. Zoller, Phys. Rev. A 87, 033606 (2013).
  * Rahmani (2015) A. Rahmani, Phys. Rev. A 92, 042110 (2015).
  * Kayanuma (1984) Y. Kayanuma, J. Phys. Soc. Jpn. 53, 108 (1984).
  * Kenmoe  _et al._ (2013) M. Kenmoe, H. Phien, M. Kiselev, and L. Fai, Phys. Rev. B 87, 224301 (2013).
  * Dóra  _et al._ (2019) B. Dóra, M. Heyl, and R. Moessner, Nat. commun. 10, 2254 (2019).
  * Gulácsi and Dóra (2021) B. Gulácsi and B. Dóra, Phys. Rev. B 103, 205153 (2021).
  * Mukherjee and Dutta (2009) V. Mukherjee and A. Dutta, J. Stat. Mech. Theory Exp. 2009, P05005 (2009).
  * Mullen  _et al._ (1989) K. Mullen, E. Ben-Jacob, Y. Gefen, and Z. Schuss, Phys. rev. lett. 62, 2543 (1989).
  * Kayanuma and Mizumoto (2000) Y. Kayanuma and Y. Mizumoto, Phys. Rev. A 62, 061401 (2000).
  * Suzuki and Iwamura (2023) T. Suzuki and K. Iwamura, Phys. Rev. B 108, 014110 (2023).
  * Patane  _et al._ (2008) D. Patane, A. Silva, L. Amico, R. Fazio, and G. E. Santoro, Phys. rev. lett. 101, 175701 (2008).
  * Patane  _et al._ (2009) D. Patane, L. Amico, A. Silva, R. Fazio, and G. E. Santoro, Phys. Rev. B 80, 024302 (2009).
  * Yin  _et al._ (2014) S. Yin, P. Mai, and F. Zhong, Phys. Rev. B 89, 094108 (2014).
  * Nalbach  _et al._ (2015) P. Nalbach, S. Vishveshwara, and A. A. Clerk, Phys. Rev. B 92, 014306 (2015).
  * Nigro  _et al._ (2019) D. Nigro, D. Rossini, and E. Vicari, Phys. Rev. A 100, 052108 (2019).
  * Rossini and Vicari (2019) D. Rossini and E. Vicari, Phys. Rev. B 100, 174303 (2019).
  * Rossini and Vicari (2020) D. Rossini and E. Vicari, Phys. Rev. Res. 2, 023211 (2020).
  * Novikov (1965) E. A. Novikov, Sov. Phys. JETP 20, 1290 (1965).
  * Sarandy and Lidar (2005) M. Sarandy and D. Lidar, Phys. Rev. A 71, 012331 (2005).
  * Gradshteyn and Ryzhik (2014) I. S. Gradshteyn and I. M. Ryzhik,  _Table of integrals, series, and products_ (Academic press, 2014).
  * Heisenberg and Euler (1936) W. Heisenberg and H. Euler, Z. Phys. 98, 714 (1936).
  * Weisskopf (1936) V. Weisskopf, Mat. Fys. Medd. 14, 6 (1936).
  * Schwinger (1951) J. Schwinger, Phys. Rev. 82, 664 (1951).
  * Schützhold  _et al._ (2008) R. Schützhold, H. Gies, and G. Dunne, Phys. Rev. Lett. 101, 130404 (2008).
  * Nasri (2016) R. Nasri, Integral Transforms Spec. Funct. 27, 181 (2016).