14pt14ptStart-up and Transient Response Characteristics of Natural Draft Direct Dry Cooling Systems at Various Scales

10pt10ptWian Strydom1,*, Johannes Pretorius1, Ryno Laubscher1

10pt10pt1Department of Mechanical and Mechatronic Engineering, University of<EMAIL_ADDRESS>
§ 11PT0PTABSTRACT

Natural draft direct dry cooling systems (NDDDCSs) combine the benefits of two traditional dry cooling systems, air-cooled condensers (ACCs) and indirect natural draft dry cooling systems, to offer advantages such as reduced system complexity, high thermal efficiency, and low parasitic power losses and operational costs. These scalable systems can potentially operate in both large coal-fired power plants (900 MWt) and medium-scale concentrated solar power (CSP) plants (100 MWt). While prior research focused on NDDDCS steady-state performance, this study delves into the transient behaviour and start-up characteristics of this cooling system, crucial for power plants with rapid start-up requirements or load ramps during electricity demand peaks. The paper employs a transient 1-D mathematical model to characterize the NDDDCSs' response to start-up and load ramp conditions under no-wind, presenting results for both large- and medium scales while identifying limiting factors.

11pt0ptKeywords: Natural draft, transient, start-up, condenser, dry cooling

§ 12PT0PT1. INTRODUCTION

Climate change challenges the energy generation industry, through increasing water scarcity [Schellnhuber and Heinke et al., 2014], energy demand, and global population growth [WEC, 2013]. These factors necessitate a shift towards more resource-efficient energy generation technologies. Dry cooling technologies, with reduced water usage compared to wet cooled systems [Macknick and Newmark et al., 2011], are natural candidates.

Two prevalent dry cooling systems, forced draft air-cooled condensers (ACCs) and indirect natural draft dry cooling systems, offer distinct advantages. The ACCs operate using a multitude of large fan drives and have a high thermal efficiency. However, the fans run virtually continuously in harsh conditions and require regular maintenance. These systems also have significant parasitic losses. Indirect natural draft systems, while avoiding fan drives, face challenges in reduced thermal efficiency and high complexity which demands substantial initial capital investments [Conradie and Buys et al., 1998].

The natural draft direct dry cooling system (NDDDCS) combines the strengths of ACCs and indirect systems, offering higher thermal efficiencies and lower operational costs. The absence of fan drives eliminates several failure modes, as passive airflow is driven by a natural draft. Limited research has been performed on NDDDCSs to date, and have focused on steady-state performance. Studies by Kong et al. using 3-D computational fluid dynamics (CFD) simulations highlight the superiority of vertically arranged heat exchanger bundles, with increased performance associated with larger apex angles [Kong and Wang et al., 2018, Kong and Wang et al., 2018]. Normal and radiator type windbreaker walls further enhance performance under crosswind conditions [Kong and Wang et al., 2018, Goodarzi and Keimanesh, 2013].

Following from Kr$\mathrm{\ddot{o}}$ger's work on ACCs and indirect systems [Kr$\mathrm{\ddot{o}}$ger, 2004], 1-D modelling techniques have proven viable for NDDDCS investigations. Strydom et al.'s sensitivity analysis emphasized the impact of geometric variation on NDDDCS performance, while also showcasing its potential as a standalone cooling system for CSP plants [Strydom and Pretorius et al., 2022, Strydom and Pretorius et al., 2023].

This study addresses a gap in the literature by exploring the start-up characteristics of NDDDCSs, crucial for power generation cycles during rapid start-ups or demand peaks. Transient performance is gaining significance as intermittent renewable technologies become more prevalent on the grid. Existing transient investigations have focussed on small-scale indirect natural draft dry cooling towers [Dong and Li et al., 2018, Dong and Li et al., 2019, Dong and Li et al., 2019], leaving a knowledge gap concerning NDDDCS start-up. The NDDDCS relies on the interdependence of the condensation process and system draft, unlike indirect systems which mechanically pump water into the heat exchangers.

This research develops transient 1-D mathematical models to quantify NDDDCS start-up and transient load ramp characteristics at two scales: large-scale coal-fired power plant (900 MWt) and medium-scale CSP plant (100 MWt). Results shed light on NDDDCS response time and potential bottlenecks in start-up and load ramp scenarios.

§ 12PT12PT2. REFERENCE TOWER GEOMETRY AND SIZING

A schematic of the NDDDCS tower and finned tube geometry is shown in Figure <ref>, with the reference tower dimensions used in this study listed in Table <ref>. The tube dimensions, heat transfer- and pressure drop characteristics are identical to those used in previous studies [Strydom and Pretorius et al., 2022, Strydom and Pretorius et al., 2023]. The total tower height is given by $H_5$, the inlet height is represented by $H_4$ and the inlet- and outlet diameters are indicated by $d_3$ and $d_5$ respectively.

Schematic of NDDDCS (a) and effective single-row finned tube (b)

NDDDCS reference dimensions at the two considered scales

Dimension       Symbol    Large-scale     Medium-scale

Total tower height (m)               $H_5$      165.0      65.0

Inlet height (m)                     $H_4$        16.5     6.5

Inlet diameter (m)                  $d_3$       165.0      65.0

Outlet diameter (m)                     $d_5$      82.5     32.5

§ 12PT12PT3. MATHEMATICAL MODEL

The steam- and air-side of the 1-D model are not spatially discretised. The governing equations are solved simultaneously, where required. Transient continuity- and energy equations are derived for the steam-side, while transient continuity-, momentum- and energy equations are derived for the air-side. It is assumed that there is no steam-side pressure drop within the condenser tubes. The density at the next time step on the steam-side is calculated via:

\begin{equation}
{\rho_{si}}^{t+\Delta t} = \frac{\Delta t}{V_c} \left(\dot{m}_{si} - {\dot{m}_{co}}^{~~{t}} \right) + {\rho_{si}}^t
\end{equation}

where ${\rho_{si}}^{t+\Delta t}$ is the steam density at the next time step, $\Delta t$ is the time step, $V_c$ is the total steam-side volume, $\dot{m}_{si}$ the supplied steam mass flow rate boundary condition from the power cycle, ${\dot{m}_{co}}^{~~{t}}$ is the mass flow rate of exiting condensate and ${\rho_{si}}^t$ is the current steam density. The internal energy of the steam is calculated by:

\begin{equation}
{u_{si}}^{t+\Delta t}    = \frac{\Delta t}{{\rho_{si}}^t V_c}\left(\dot{m}_{si} {h_{si}}^t - {\dot{m}_{co}}^{~~{t}} {h_{co}}^t - {\dot{Q}_r}^{~~t} - {u_{si}}^t (\dot{m}_{si} - {\dot{m}_{co}}^{~~{t}} )\right) + {u_{si}}^t
\end{equation}

where ${u_{si}}^{t+\Delta t}$ is the internal energy of steam at the next time step, ${u_{si}}^t$ is the current internal energy, ${h_{si}}^t$ and ${h_{co}}^t$ are the respective enthalpies of entering steam and exiting condensate and ${\dot{Q}_r}^{~~t}$ is the steam-side heat transfer rate. The latter is given by:

\begin{equation}
{\dot{Q}_r}^{~~{t}} = \dot{m}_a^{~~t} {c_{pa25}}^t e^t ({T_w}^t - T_{a2}) = {h_c}^t A_c ({T_{si}}^t - {T_w}^t)
\label{eq:Q_steam_tran}
\end{equation}

where ${T_w}^t$ represents the tube wall temperature, $T_{a2}$ the inlet air temperature, ${T_{si}}^t$ represents the saturated inlet steam temperature, $\dot{m}_a^{~~t}$ is the current air mass flow rate, ${c_{pa25}}^t$ is the average specific heat capacity of air, $e^t$ is the heat exchanger effectiveness, ${h_c}^t$ is the condensation heat transfer coefficient and $A_c$ is the total steam-side wetted area. The tube wall temperature is assumed to be uniform between the inner and outer surfaces of the condenser tubes. Rearranging Equation <ref> to calculate the tube wall temperature yields:

\begin{equation}
{T_w}^t = \frac{{h_c}^t A_c {T_{si}}^t + \dot{m}_a^{~~t} {c_{pa25}}^t e^t T_{a2}}{{h_c}^t A_c + \dot{m}_a^{~~t} {c_{pa25}}^t e^t}
\end{equation}

The heat exchanger effectiveness, $e^t$ is given by:

\begin{equation}
e = 1 - \mathrm{exp}\left( \frac{UA^t}{\dot{m}_{a}^{~~t} {c_{pa25}}^t} \right)
\end{equation}

Due to the inclusion of the condenser tube wall temperature in the energy balance of Equation <ref>, the overall heat conductance, $UA$, is equivalent to the air-side heat conductance. The mass flow rate of condensate leaving the condenser tubes is calculated as:

\begin{equation}
\dot{m}_{co}^{~~t} = \dot{Q}_r^{~~t}/{h_{fg}}^t
\end{equation}

where ${h_{fg}}^t$ is the latent heat of vaporization. The air-side heat transfer, $\dot{Q}_a^{~~t}$, is given by:

\begin{equation}
\dot{Q}_a^{~~t} = \dot{m}_{a}^{~~t} ({h_{a5}}^t - h_{a2})
\end{equation}

where ${h_{a5}}^t$ and $h_{a2}$ are the enthalpies of entering and exiting air, respectively. The air mass flow rate at the next time step is calculated via:

\begin{equation}
\dot{m}_{a}^{~~t+\Delta t} = \frac{A_{35} \Delta t}{L_a} \left(\Delta {P_a}^t - \Delta {P_{L}}^t \right) + \dot{m}_{a}^{~~t}
\label{eq:momentum_discretized_ch8}
\end{equation}

where $\Delta {P_{a}}^t$ is the system driving force, $\Delta {P_{L}}^t$ is the total air-side pressure loss, $A_{35}$ is the average tower flow area and $L_a$ is the air travel distance. The pressure balance in Equation <ref> is equivalent to the draft equation taken from the work of [Kr$\mathrm{\ddot{o}}$ger, 2004]. Thus, the driving force is:

\begin{equation}
\Delta {P_a}^t =  p_{a1}\left[\left(1-0.00975\frac{H_4}{2T_{a1}} \right)^{3.5}  \left(1-0.00975\frac{H_5 - \frac{H_4}{2}}{{T_{a3}}^t} \right)^{3.5} - \left(1-0.00975\frac{H_5}{T_{a1}} \right)^{3.5} \right]
\end{equation}

where $p_{a1}$ is the ground-level atmospheric pressure, $T_{a1}$ is the corresponding temperature and ${T_{a3}}^t$ is the outlet air temperature after the heat exchangers. The flow resistance is given by:

\begin{equation}
\Delta {P_{L}}^t =  \frac{(K_{il} + {K_{he\theta}}^t + {K_{ct}}^t + {K_{ts}}^t)_{he}}{2{\rho_{a25}}^t}
\left(\frac{\dot{m}_{a}^{~~t}}{A_{fr}} \right)^2
\left(1-0.00975\frac{H_5-\frac{H_4}{2}}{{T_{a3}}^t}\right)^{3.5}
+ \frac{({K_{to}}^t + {a_{e5}}^t)}{2{\rho_{a5}}^t}\left(\frac{\dot{m}_{a}^{~~t}}{A_{5}} \right)^2
\end{equation}

where ${K_{il}}$ (constant value of 2.5), ${K_{he\theta}}^t$, ${K_{ct}}^t$, ${K_{ts}}^t$, ${K_{to}}^t$ and ${a_{e5}}^t$ represent the non-dimensional inlet louvre-, heat exchanger-, tower-, tower support-, outlet pressure loss and velocity correction factors. The tower outlet air density is represented by ${\rho_{a5}}^t$, the corresponding area is given by $A_5$ and the frontal area of the heat exchangers is represented by $A_{fr}$. The enthalpy of air leaving the NDDDCS is calculated via:

\begin{equation}
{h_{a5}}^{t+\Delta t}  =  \frac{\Delta t}{{\rho_{a25}}^t V_a} \left( \dot{m}_{a}^{~~t}  (h_{a2} -  {h_{a5}}^t)  + \dot{Q}_{r}^{~~t}  \right) + {h_{a5}}^{t}
\end{equation}

where $V_a$ is half of the internal flow volume. The equations for the loss factors and overall heat conductance are presented in previous papers [Strydom and Pretorius et al., 2022, Strydom and Pretorius et al., 2023].

§ 12PT12PT4. VALIDATION AND SENSITIVITY ANALYSIS

The 1-D steady state NDDDCS model (SS) was previously validated against the work of Kong et al. [Strydom and Pretorius et al., 2022]. Therefore, this study validates the transient 1-D NDDDCS model by comparing the results of the transient model, once it reaches a steady state, to the results of the steady-state model. The large scale tower is used as reference and the same effective boundary conditions are imposed to allow comparison. It is important to note that two transient models were developed: transient steam-side combined with steady-state air-side (TS) and transient steam- and air-side (TT). The SS model was adapted to make use of Equation <ref> and to neglect the steam-side pressure drop for the sake of consistency between the three models. A constant ambient temperature of 293.15 K (20 $^\circ$C) was selected. The steam mass flow rate (444.98 kg/s) attained by the SS model for a constant saturated steam temperature boundary condition of 323.15 K is assigned as a boundary condition for the TS and TT models. The TT model requires an initial steam temperature and air mass flow rate, specified as the ambient temperature and as effectively zero, respectively. The results are summarized in Table <ref>. It is clear from Table <ref> that the TS and TT models converge to the same steady-state condition as the SS model. A sensitivity analysis revealed that a time-step of 0.2 s is suitable for this simulation.

Summary of steady state results for cross validation of TT model

Model       $T_{si}$ (K)       $T_{w}$ (K)     $T_{a5}$ (K)     $v_s$ (m/s)     $\dot{Q}$ (MW)          $\dot{m}_a$ (kg/s)     $\Delta P$ (Pa)

SS         323.15          322.35            315.45             99.58                1059.92         47033.71            137.64

TS         323.15          322.35            315.45             99.58                1059.92         47033.70            137.64

TT         323.17          322.38            315.48             99.48                1059.90         46986.32            137.76

§ 12PT12PT5. RESULTS

All properties are calculated using the CoolProp library within the Python environment. The steam mass flow rate ($\dot{m}_{si}$, step input), ambient temperature ($T_{a1}$, 293.15 K) and pressure ($p_{a1}$, 101325 Pa) at ground level are boundary conditions, while the outlet air temperature ($T_{a3}$, 293.15 K), air mass flow rate ($\dot{m}_a$, $\approx$0 kg/s) and saturated steam temperature ($T_{si}$, 295.15 K) are initial conditions. Figure <ref> serves to illustrate why simulating only the steam-side transiently (TS model) is not adequate to capture the full spectrum of transient start-up effects in NDDDCSs.

The maximum allowable internal pressure is selected as 160 kPa (a), which is based on an approximate rupture disk burst limit, typically installed on the steam ducting to avoid system overpressure. Results indicate that the large-scale NDDDCS can handle 9.6% of its maximum steam admission rate as a step input without exceeding this value. Under this admission rate, the TS and TT models calculate vastly different results. Due to the air-side being calculated steadily for the TS model, the airflow through the NDDDCS can respond instantaneously to changes in the steam-side to provide adequate heat rejection. This is not the case for the TT model, as the response time of the air-side is significantly slower than the steam-side. The initial air mass flow rate is not enough to reject the thermal energy provided by the steam, and an acute build up of steam back pressure is observed until the air mass flow rate is sufficient.

Large-scale system steam back pressure predicted by TS (a) and TT models (b)

Figure <ref> (a, b) shows the transient heat balance and tube wall temperature of the large-scale NDDDCS under the same step input previously mentioned.  The heat addition to the tube wall increases rapidly due to a lack of airflow under start-up, leading to a build up in tube wall temperature. An inflection point is reached in the heat addition curve, where the heat removal from the tube wall increases and the heat addition to the tube wall settles as the imbalance in heat transfer starts to reduce, leading to a reduction in steam back pressure. The steam-side settles about five times faster (30 s) than the air-side (150 s).

Large-scale heat balance (a) and tube wall temperature and air mass flow rate (b) graphs using TT model under maximum steam admission rate step input

The large- (coal-fired) and medium scale (CSP) NDDDCS back pressure and air mass flow rate trends are shown in Figure <ref> for typical start-up steam admission rates, acquired from proprietary sources.

Steam back pressure (a, c) and air mass flow rate (b, d) curves under typical start-up load ramp conditions for coal-fired power plant (a, b) and CSP plant (c, d) using TT model

The large- and medium-scale NDDDCSs experience controlled spikes in steam back pressure (Figure <ref> (a, c)) initially due to airflow in the tower not being established yet, after which the air mass flow rate increases rapidly, as seen in Figure <ref> (b, d). As the air mass flow rate continues to increase, the NDDDCSs can easily supply the required heat rejection rate to match the transient cycle requirements, which controls the back pressure. This indicates that for a typical coal-fired or CSP plant, the NDDDCS will not impose transient operation limits, as the fastest cycle load ramp rates are still slower than the response time of the cooling system. This belays past concerns about the transient response time of natural draft steam condensing systems, which rely on the interplay between the condensing process and draft creation.

§ 12PT12PT6. CONCLUSION

This study investigated the transient start-up characteristics of NDDDCSs at large- and medium-scales under transient start-up and load ramp scenarios, for windless atmospheric conditions. Results indicate that the steam-side responds much faster than the air-side of the system during unsteady conditions. The study finds that for the typical transient scenarios experienced at a coal-fired or CSP plant, the NDDDCS is able to respond rapidly enough to avoid back pressure increases that would limit turbine output during the load ramp.

§ 12PT12PTACKNOWLEDGEMENTS

This work is supported by the National Research Foundation (NRF) of South Africa and the Solar Thermal Energy Research Group (STERG) of the University of Stellenbosch.

§ 12PT12PTNOMENCLATURE

[P01]$A$Area ($\mathrm{m^2}$)
[P02]$a$Major axis or velocity correction factor ($\mathrm{m~or}~-$)
[P03]$b$Minor axis ($\mathrm{m}$)
[P04]$c$Specific heat capacity ($\mathrm{J/kg\cdot K}$)
[P05]$d$Diameter ($\mathrm{m}$)

[P07]$g$Gravitational acceleration ($\mathrm{m/s^{2}}$)
[P08]$H$Height ($\mathrm{m}$)
[P09]$h$ Heat transfer coefficient or latent heat of vaporization ($\mathrm{W/m^{2}\cdot K}$ or $\mathrm{J/kg}$)

[P11]$K$Loss factor ($\mathrm{-}$)
[P12]$k$Thermal conductivity ($\mathrm{W/m\cdot K}$)
[P13]$L$Length ($\mathrm{m}$)

[P15]$\dot{m}$Mass flow rate ($\mathrm{kg/s}$)

[P18]$p$Pressure or pitch ($\mathrm{Pa~or~m}$)
[P19]$Pr$Prandtl number ($\mathrm{-}$)
[P20]$\dot{Q}$Heat transfer rate ($\mathrm{W}$)

[P24]$T$Temperature ($\mathrm{K}$)
[P25]$UA$Overall heat conductance ($\mathrm{W/K}$)

[P28]$v$Velocity ($\mathrm{m/s}$)

[G0]Greek Symbols
[G1]$\Delta$Change or thickness ($\mathrm{-~or~m}$)
[G2]$\delta$Thickness ($\mathrm{m}$)
[G3]$\epsilon$NTU-effectiveness ($\mathrm{-}$)

[G4]$\mu$Dynamic viscosity ($\mathrm{kg/m\cdot s}$)

[G5]$\rho$Density ($\mathrm{kg/m^3}$)

[S1]$a$Air, Atmospheric
[S2]$c$Condensate, Cooling

[S5]$he$Heat exchanger
[S6]$i$Inlet, In

[S9]$o$Outlet, Out

[S11]$s$Steam, Support
[S12]$t$Tube, Tower

[S14]$3$Tower inlet diameter
[S15]$4$Tower inlet height
[S16]$5$Tower outlet diameter

[Z2]ACC     Air-cooled condenser
[Z4]NDDDCSNatural draft direct dry cooling system
[Z1]3-D     3-dimensional
[Z3]CFD     Computational fluid dynamics

[Schellnhuber and Heinke et al., 2014]
H.J. Schellnhuber and J. Heinke et al.
Multimodel assessment of water scarcity under climate change.
Proceedings of the National Academy of Sciences - PNAS, 1110 (9):0 3245–3250, 2014.
ISSN 0027-8424.

[WEC, 2013]
World energy scenarios.
Technical report, World Energy Council, 2013.

[Macknick and Newmark et al., 2011]
J. Macknick and R. Newmark et al.
Review of operational water consumption and withdrawal factors for electricity generating technologies.

[Conradie and Buys et al., 1998]
A.E. Conradie and J.D. Buys et al.
Performance optimization of dry-cooling systems for power plants through sqp methods.
Applied Thermal Engineering, 180 (1):0 25–45, 1998.
ISSN 1359-4311.

[Kong and Wang et al., 2018]
Y. Kong and W. Wang et al.
Thermo-flow performances of natural draft direct dry cooling system at ambient winds.
International Journal of Heat and Mass Transfer, 116:0 173–184, 2018a.
ISSN 0017-9310.

[Kong and Wang et al., 2018]
Y. Kong and W. Wang et al.
Annularly arranged air-cooled condenser to improve cooling efficiency of natural draft direct dry cooling system.
International Journal of Heat and Mass Transfer, 118:0 587–601, 2018b.
ISSN 0017-9310.

[Kong and Wang et al., 2018]
Y. Kong and W. Wang et al.
Wind leading to improve cooling performance of natural draft air-cooled condenser.
Applied Thermal Engineering, 136:0 63–83, 2018c.
ISSN 1359-4311.

[Goodarzi and Keimanesh, 2013]
M. Goodarzi and R. Keimanesh.
Heat rejection enhancement in natural draft cooling tower using radiator-type windbreakers.
Energy Conversion and Management, 71:0 120–125, 2013.

[Kr$\mathrm{\ddot{o}}$ger, 2004]
D.G. Kr$\mathrm{\ddot{o}}$ger.
Air-cooled heat exchangers and cooling towers.
Penwell Corp., Tulsa, Okla, 2004.
ISBN 9780878148967.

[Strydom and Pretorius et al., 2022]
W. Strydom and J.P. Pretorius et al.
Natural draft direct dry cooling system scaling for diverse applications.
16th International Conference on Heat Transfer, Fluid Mechanics and Thermodynamics, 2022.

[Strydom and Pretorius et al., 2023]
W. Strydom and J.P. Pretorius et al.
Sensitivity analysis on the performance of a natural draft direct dry cooling system for a 50mwe csp application.
8th World Congress on Momentum, Heat and Mass Transfer, 2023.

[Dong and Li et al., 2018]
P. Dong and X. Li et al.
The transient start-up process of natural draft dry cooling towers in dispatchable thermal power plants.
International Journal of Heat and Mass Transfer, 123:0 201–212, 2018.
ISSN 0017-9310.

[Dong and Li et al., 2019]
P. Dong and X. Li et al.
Numerical investigation of the influence of local effects on the transient start-up process of natural draft dry cooling towers in dispatchable power plants.
International Journal of Heat and Mass Transfer, 133:0 166–178, 2019a.
ISSN 0017-9310.

[Dong and Li et al., 2019]
P. Dong and X. Li et al.
The crosswind effects on the start-up process of natural draft dry cooling towers in dispatchable power plants.
International Journal of Heat and Mass Transfer, 135:0 950–961, 2019b.
ISSN 0017-9310.


# Nonadiabatic transitions in Landau-Zener grids: integrability and
semiclassical theory

Rajesh K. Malla  Theoretical Division, and the Center for Nonlinear Studies,
Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA  Vladimir Y.
Chernyak  Department of Chemistry, Wayne State University, 5101 Cass Ave,
Detroit, Michigan 48202, USA  Department of Mathematics, Wayne State
University, 656 W. Kirby, Detroit, Michigan 48202, USA  Nikolai A. Sinitsyn
Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico
87545, USA

###### Abstract

We demonstrate that the general model of a linearly time-dependent crossing of
two energy bands is integrable. Namely, the Hamiltonian of this model has a
qaudratically time-dependent commuting operator. We apply this property to
four-state Landau-Zener (LZ) models that have previously been used to describe
the Landau-Stückelberg interferometry experiments with an electron shuttling
between two semiconductor quantum dots. The integrability then leads to simple
but nontrivial exact relations for the transition probabilities. In addition,
the integrability leads to a semiclassical theory that provides analytical
approximation for the transition probabilities in these models for all
parameter values. The results predict a dynamic phase transition, and show
that similarly-looking models belong to different topological classes.

##  I Introduction

The Landau-Zener (LZ) model describes an evolution for amplitudes of two
states with a time-dependent Hamiltonian

$H=\left(\begin{array}[]{cc}b_{1}t&g\\\ g^{*}&b_{2}t\end{array}\right),$ (1)

where $b_{1,2}$ are called slopes of diabatic levels and $g$ is the inter-
level coupling. The basis in which the off-diagonal elements of $H(t)$ are
time-independent is called diabatic basis. The LZ formula provides an exact
analytical expression for the probability to remain in the same diabatic state
after the evolution during time $t\in(-\infty,+\infty)$:

$P_{LZ}=e^{-2\pi|g|^{2}/|b_{1}-b_{2}|}.$ (2)

This formula plays a special role in the theory of nonadiabatic transitions
because it can be used as an approximation when the energy levels are mostly
well separated. The adiabaticity is then broken only in disjoint regions of
time-energy, in which the nonadiabatic dynamics is experienced only by pairs
of states and the parameter time-dependence can be linearized.

For nanoscale systems of modern interest, however, many states may experience
the nonadiabatic transitions simultaneously, even when the linear
approximation of the parameter time-dependence near the nonadiabatic
transitions is still applicable. The state evolution is then described by the
nonstationary Schrödinger equation

$i\frac{d}{dt}|\psi\rangle=H(t)|\psi\rangle,$ (3)

and the time-dependent Hamiltonian of a multistate Landau-Zener (MLZ) process
has generally the form be

$H(t)=Bt+A,$ (4)

where $A$ and $B$ are time-independent matrices, and $B$ is diagonal. Let $E$
be the diagonal part of $A$. The nonzero elements of $Bt+E$ are called
diabatic energies and the corresponding eigenstates are called diabatic
states. As $t\rightarrow\pm\infty$, the diabatic states coincide with the
Hamiltonian eigenstates. The goal of the MLZ theory is to find the amplitudes
$S_{nm}$ and the transition probabilities, $P_{m\rightarrow n}=|S_{nm}|^{2}$,
from the diabatic states $m$ as $t\rightarrow-\infty$ to the states $n$ as
$t\rightarrow+\infty$.

Figure 1: The linearly time-dependent diabatic levels in the model of two
crossing bands are forming a pattern that we call LZ-grid. The first and the
second bands have, respectively, $N$ and $M$ parallel levels. The diabatic
basis states of the same band do not interact with each other directly but any
such a state can be coupled directly to arbitrary diabatic states of the other
band.

Among the MLZ models, there is a class of Hamiltonians that has attracted
special attention previously. It corresponds to the time-dependent crossing of
two bands with parallel diabatic levels, as shown in Fig. 1. Let $N$ and $M$
be the integer numbers of the parallel levels in these bands. Matrices $A$ and
$B$ then have the dimensions $(N+M)\times(N+M)$, and

$A=\left(\begin{array}[]{cccc}E_{1}&G\\\
G^{\dagger}&E_{2}\end{array}\right),\quad
B=\left(\begin{array}[]{cccc}b_{1}{1}_{N}&0\\\
0&b_{2}{1}_{M}\end{array}\right),$ (5)

where $1_{N}$ and $1_{M}$ are the unit, respectively, N$\times$N and
M$\times$M matrices. The diagonal matrices

$E_{1}\equiv{\rm
diag}\\{e_{1}^{1},e_{1}^{2},\ldots,e_{1}^{N}\\},\,\,E_{2}\equiv{\rm
diag}\\{e_{2}^{1},e_{2}^{2},\ldots,e_{2}^{M}\\}$

are responsible for the spacing between the parallel diabatic levels;
$b_{1,2}$ are the slopes of the bands. The lower indices $1$ and $2$ in
$b_{1,2}^{i},\,e_{1,2}^{i}$ refer, respectively, to the $N$-level and
$M$-level bands. $G$ is a $N\times M$ matrix that describes direct coupling
between the two bands. All elements of $G$ can be nonzero and complex-valued.

The band crossing MLZ model (5) was discussed originally in relation to
physics of the Rydberg atoms harmin1 ; harmin2 . Optical realization of this
model was used to create an optical Galton board optics-lzgrid . The early
work sometimes referred to the level crossing pattern in the two-band model as
to LZ-grid of energy levels. Later, LZ-grids attracted attention in relation
to the two-state systems that are coupled to an environment, which splits the
two levels of the LZ model into the two bands of many parallel levels usuki ;
sinitsyn-bath ; hanggie ; garanin-grid ; ashhab-lz , and more recently LZ-
grids emerged in the study of qubits coupled to optical modes, such as in
circuit QED systems osc-grid1 ; osc-grid2 ; osc-grid3 ; osc-grid4 ; recent-
grid1 ; recent-grid2 .

Although certain facts about the LZ-grids have been determined analytically
usuki , such systems remain generally unsolvable. Due to the complex
oscillatory behavior of the transition probabilities as functions of the
parameters, physics of LZ-grids remains poorly studied. Approximations have
been developed but only for limits of either very small ostrovsky-grid  or
very large yurovsky-grid  separations of the parallel levels in the bands.
Such limits are approached very slowly with decreasing/growing $e^{i}_{1,2}$,
so they usually give too crude approximations to realistic choices of the
parameters. There is also one fully solvable LZ-grid model quest-LZ  but it
does not clarify many questions about the general problem.

Our article has two goals. First, we add to the analytical understanding of
the transition probabilities in all LZ-grids by applying the recent
developments on the time-dependent integrability commute ; quadr-LZ . Thus, in
section II, we show that all LZ-grids are formally integrable, in the sense
that there is an analytic expression for the time-quadratic-polynomial
operator that commutes with the Hamiltonian, and satisfies an additional
condition that is needed for the integrability of time-dependent Hamiltonians
as defined in commute . However, unlike many known solvable cases laser-LZ ;
bcs ; gamma-LZ ; yuzbashyan-LZ , LZ-grids generally remain not fully solvable.
Instead, the integrability leads to a simple but nontrivial symmetry for the
transition probabilities.

Our second goal is to apply this symmetry to a specific four-state LZ-grid
model that has attracted attention recently due to experiments with Landau-
Stückelberg interferometry note  in coupled quantum dots petta-4lz1 ;
petta-4lz2 . We will assume that the electrostatic energy of one of the dots
is changing linearly with time and, in section III, show that the
integrability in such models leads to nontrivial exact relations between
different transition probabilities.

Finally, one of the applications of the integrability is the possibility to
reduce the order of the differential equation for amplitudes of LZ-grid
states. This property leads to asymptotically exact expressions for leading
exponents that describe the transition probabilities in the nearly-adiabatic
limit, as it has been recently demonstrated for the general three-state LZ
model quadr-LZ . In section IV we apply this semiclassical approach to the
four-state quantum dot models in order to derive approximate expressions for
the experimentally most relevant probability to remain in the same quantum dot
after a linear sweep of the gate voltage.

##  II Integrability of the general band-crossing model

The integrability conditions for a time-dependent Hamiltonian $H(t)$ are
defined as the possibility to find a parameter combination $\tau$, and an
analytical form of a nontrivial operator $H^{\prime}$ such that commute

$\displaystyle\frac{\partial H}{\partial\tau}-\frac{\partial
H^{\prime}}{\partial t}$ $\displaystyle=$ $\displaystyle 0,$ (6)
$\displaystyle[H,H^{\prime}]$ $\displaystyle=$ $\displaystyle 0.$ (7)

Such conditions are known in the theory of solitons faddeev-book , and the
existence of time-polynomial commuting operators in several MLZ systems was
originally noticed in patra-LZ . Before we discuss the applications, let us
first prove that the conditions (6) and (7) can be always satisfied for the
LZ-grid Hamiltonian with $A$ and $B$ given by (5).

Let us define a continuous family of operators

$H(t,\tau)=B(\tau)t+A(\tau),$ (8)

where $B(\tau)$ and $A(\tau)$ are obtained from the original $B$ and $A$ by
setting

$B(\tau)\equiv B\tau,\quad E_{1}(\tau)\equiv\tau E_{1},\quad G(\tau)\equiv
G\sqrt{\tau},$ (9)

and keeping $E_{2}$ intact. Note that at $\tau=1$, $A(\tau)$ and $B(\tau)$ are
the same as original $A$ and $B$.

Then, the pair of operators, $H(t,\tau)$ and $H^{\prime}(t,\tau)$, where

$H^{\prime}(t,\tau)=\frac{\partial_{\tau}B(\tau)t^{2}}{2}+\partial_{\tau}A(\tau)t-\frac{1}{2(b_{2}-b_{1})\tau^{2}}A^{2}(\tau),$
(10)

satisfy (6) and (7). Indeed, (6) is trivial to verify, whereas (7) leads to
two independent conditions for the terms proportional to, separately, $t^{2}$
and $t$:

$\displaystyle\frac{1}{2}[\partial_{\tau}B(\tau),A(\tau)]+[\partial_{\tau}A(\tau),B(\tau)]=0,$
(12)
$\displaystyle[\partial_{\tau}A(\tau),A(\tau)]-\frac{1}{2(b_{2}-b_{1})\tau^{2}}[A^{2}(\tau),B(\tau)]=0,$

which can be verified by direct substitution of the $\tau$-dependent matrices
(8) and (9).

Due to the satisfied integrability conditions, one can deform the integration
path in the two-time space $(t,\tau)$ without changing the evolution
amplitudes. Namely, let us define the evolution operator

$U=\hat{\cal T}_{\cal P}\exp\left(-i\int_{\cal
P}H(t,\tau)\,dt+H^{\prime}(t,\tau)\,d\tau\right),$ (13)

where $\hat{\cal T}_{\cal P}$ is the path ordering operator along ${\cal P}$
in the two-time space $(t,\tau)$. Equations (6)-(7) mean that the nonabelian
gauge field with components ${\bf A}(t,\tau)=(H,H^{\prime})$ has zero
curvature, so the result of integration in (13) does not change after the
deformations of ${\cal P}$ that keep only the initial and final points of
${\cal P}$ intact commute , and avoid singularities of the $\tau$-dependent
Hamiltonians, as in Fig. 2.

Figure 2: The true time-evolution path ${\cal P}$ (blue arrows) with $\tau=1$
and $t\in(-\infty,+\infty)$ can be deformed into the path ${\cal P}_{\tau}$,
such that the horizontal part of ${\cal P}_{\tau}$ has $\tau={\rm const}\neq
1$ (dashed black arrows). These deformations do not change the evolution
matrix. Vertical legs of ${\cal P}_{\tau}$ have $t=\pm T$ with
$T\rightarrow\infty$, so they contribute only to the trivial adiabatic phases
in the evolution matrix, and do not affect the transition probabilities. For
the three-state LZ model, the path ${\cal P}_{\tau}$ can be chosen so that
$b_{1}\rightarrow\infty$ along the horizontal piece of this path.

Let the physical evolution correspond to the changes of $t$ from $-\infty$ to
$+\infty$ at $\tau=1$. Then, ${\cal P}$ starts at the point
$(t,\tau)=(-\infty,1)$. We fix, initially, $t$ and change $\tau$ from this
point to another value, and only then perform $t$-evolution at fixed new
$\tau$. After this, we bring $\tau$ back to $\tau=1$ at $t=+\infty$ commute ;
faddeev-book .

The $\tau$-evolution at fixed $t=-T\rightarrow-\infty$ or
$t=T\rightarrow+\infty$ is strictly adiabatic due to the quadratic dependence
of the diagonal elements of $H^{\prime}$ on $t$. Therefore, the transition
probabilities in the LZ-grid models that differ only by $\tau$ within the
family (9), which is parametrized by $\tau$, are identical.

This nontrivial invariance can now be used in combination with a trivial
symmetry that is common for all MLZ models. Namely, by rescaling time in the
Schrödinger equation (3),

$t\rightarrow t/\sqrt{\tau},$ (14)

we cannot change the transition probabilities for the evolution in the
interval $t\in(-\infty,\infty)$ at fixed $\tau$. On the other hand, this
rescaling corresponds to the change of the parameters in the original model
(5):

$b_{1,2}\rightarrow b_{1,2}/\tau,\quad E_{1,2}\rightarrow
E_{1,2}/\sqrt{\tau},\quad G\rightarrow G/\sqrt{\tau}.$ (15)

Thus, the transition probabilities are independent of the variable
transformations, simultaneously, (9) and (15). Combining them, we find that
the transition probabilities in model (5) are invariant of a simple
transformation of the diagonal matrices:

$E_{1}\rightarrow E_{1}\sqrt{\tau},\quad E_{2}\rightarrow E_{2}/\sqrt{\tau}.$
(16)

This is the most general exact result of our article. The physical meaning of
this result is illustrated in Fig. 3. For a model with only two levels in each
band in this figure, there are only two independent level splittings

$\Delta e_{1}\equiv e_{1}^{1}-e_{1}^{2},\quad\Delta e_{2}\equiv
e_{2}^{1}-e_{2}^{2},$

so Eq. (16) means that the transition probabilities depend only on the
combination $\Delta e_{1}\Delta e_{2}$ but not on the ratio $\Delta
e_{1}/\Delta e_{2}$. It is easy to verify that

$S=\Delta e_{1}\Delta e_{2}/(b_{1}-b_{2})$

has the physical meaning of an area enclosed by the diabatic levels (the
diamond plaquette in Fig. 3). Hence, we can also formulate (16) as an
invariance of the transition probabilities of the transformations of the LZ-
grid that preserve the areas enclosed by the diabatic levels, as well as the
LZ parameters $|G_{ij}|^{2}/(b_{1}-b_{2})$ if the slopes of the bands are also
allowed to change.

Figure 3: Varying the distances between parallel levels does not change the
state-to-state transition probabilities if the area enclosed by the diabatic
levels (filled by blue color) is conserved.

Unfortunately, the invariance under the transformations (16) is not sufficient
to solve the whole model, i.e., to express the transition probabilities in
terms of the known special functions of the model’s parameters. Nevertheless,
we will show that this symmetry strongly simplifies the analysis, and even
leads to certain further exact relations for the transition probabilities when
a model has additional discrete symmetries.

##  III Quantum dot models

In what follows, we will explore application of the symmetry (16) to two
models, which have been studied to some extent for various reasons previously.
We will refer to these models as to symmetric and antisymmetric. The symmetric
model has the Hamiltonian

$H_{1}(t)=\left(\begin{array}[]{cccc}e_{1}&0&g&g\\\ 0&-e_{1}&g&g\\\
g&g&bt+e_{2}&0\\\ g&g&0&bt-e_{2}\end{array}\right),$ (17)

and the antisymmetric one has the Hamiltonian

$H_{2}(t)=\left(\begin{array}[]{cccc}e_{1}&0&g&-\gamma\\\ 0&-e_{1}&\gamma&g\\\
g&\gamma&bt+e_{2}&0\\\ -\gamma&g&0&bt-e_{2}\end{array}\right),$ (18)

where all parameters are real. The symmetric model has emerged previously in
discussions of nonadiabatic behavior in MLZ systems in the large coupling
limit malla-LZ .

Figure 4: The antisymmetric model describes a single electron with spin
shuttling between two quantum dots with the Hamiltonian (18). Dashed lines
show pairs of energy levels in each dot. Up and down red arrows mark the
orthogonal spin states. The splitting within each pair is the effect of an
external static magnetic field. The relative energies of localized states in
different dots are controlled by time-dependent electric gate voltage. The
couplings $g$ and $\gamma$ describe, respectively, spin-conserving and spin-
flipping electron tunnelings that are permitted by the time-reversal symmetry
of the system in the absence of the magnetic field.

Physically, both models, (17) and (18), can describe a single electron that
jumps between discrete levels of two quantum dots. The linear potential ramp
in one of the dots is then induced by applying a time-dependent electric gate
voltage difference between the dots, as in the experiments petta-4lz1 ;
petta-4lz2 . Each dot has two discrete energy levels. Most naturally, this
happens when electronic spin can flip during the tunneling event due to the
spin orbit coupling, as shown in Fig. 4.

In the antisymmetric case, with the Hamiltonian $H_{2}$, the couplings $g$ and
$\gamma$ describe then the spin preserving and the spin flipping tunneling
between the two dots. The minus sign near $\gamma$ guarantees the time-
reversal invariance of the model. In fact, the Hamiltonian $H_{2}$ at
$e_{1}=e_{2}=0$ is the most general, up to gauge transformations, Hamiltonian
that one can create for a time-reversed four-state system with spin four-LZ ,
and the splittings $e_{1},\,e_{2}\neq 0$ are induced by applying the external
magnetic field. Generally, different quantum dots have different g-factors,
which then means that $e_{1}\neq e_{2}$.

Both models, $H_{1}$ and $H_{2}$, have elementary discrete symmetries. For
example, let

$\Theta=\left(\begin{array}[]{cccc}0&-1&0&0\\\ -1&0&0&0\\\ 0&0&0&1\\\
0&0&1&0\end{array}\right).$

The Hamiltonian $H_{1}$ has an elementary symmetry

$H_{1}(t)=-\Theta H_{1}(-t)\Theta,$

which leads to the symmetry of the evolution operator $U(T|-T)={\cal
T}\exp(-i\int H_{1}(t)\,dt)$:

$U=\Theta U^{\dagger}\Theta,$ (19)

from which follows that the amplitudes $U_{12}$ and $U_{43}$ in the symmetric
model are purely real. There are also relations between different scattering
amplitudes, such as, $U_{13}=-U_{42}^{*}$, $U_{14}=-U_{32}^{*}$, which lead to
the relations between the transition probabilities: $P_{3\rightarrow
1}=P_{2\rightarrow 4}$, $P_{4\rightarrow 1}=P_{2\rightarrow 3}$,
$P_{3\rightarrow 2}=P_{1\rightarrow 4}$, and $P_{4\rightarrow
2}=P_{1\rightarrow 3}$. For a reader interested in examples of discrete
symmetry effects in other MLZ models, we refer to Refs. four-LZ ; li-dynamic .

Similarly, for the antisymmetric model, let

$\Theta_{A}=\left(\begin{array}[]{cccc}0&-i&0&0\\\ i&0&0&0\\\ 0&0&0&i\\\
0&0&-i&0\end{array}\right),$

then

$H_{2}(t)=-\Theta_{A}H_{2}(-t)\Theta_{A},$

from which follows that $U_{12}$ and $U_{43}$ are purely imaginary in this
model, and there is the same set of relations between the transition
probabilities as for the symmetric model.

Apart from this, the transition probability independence of the area
conserving transformations (Fig. 3), leads to less intuitive constraints.
Namely, since different quantum dots generally have different level
splittings, i.e.,

$e_{1}\neq e_{2}$

there is generally an asymmetry of dynamics in respect to the initially chosen
quantum dot. For example, the trivial symmetry (19) does not predict any
relation between the probabilities $P_{2\rightarrow 1}$ and $P_{3\rightarrow
4}$, as we illustrate in Fig. 5(a).

However, the invariance of the transition probabilities of the transformations
(16) allows us to tune $e_{1}=e_{2}$ in both models without affecting the
transition probabilities. Let us also add the gauge transformation
$H_{1,2}\rightarrow H_{1,2}-b1_{4}/2$, which does not change the probabilities
either. At such values of the parameters, the Hamiltonians (17) and (18) have
an additional discrete symmetry. Namely, let

$\theta=\left(\begin{array}[]{cccc}0&0&0&-1\\\ 0&0&-1&0\\\ 0&1&0&0\\\
1&0&0&0\end{array}\right).$

Then, $H_{1,2}^{\prime}=H_{1,2}^{e_{1}=e_{2}}-b1_{4}/2$ satisfy

$H^{\prime}_{1,2}(t)=\theta H^{\prime}_{1,2}(t)\theta,$

from which we obtain additional relations on the transition amplitudes, such
as $U_{12}=U_{43}^{*}$ and consequently $P_{2\rightarrow 1}=P_{3\rightarrow
4}$, e.t.c.. Let us now summarize all the relations among the transition
probabilities that follow from simultaneously the integrability and the
discrete symmetries of the models. Although the latter symmetries are a bit
different for $H_{1}$ and $H_{2}$, they lead to the same relations between the
transition probabilities:

$\displaystyle P_{1\rightarrow 3}$ $\displaystyle=$ $\displaystyle
P_{4\rightarrow 2},$ (20) $\displaystyle P_{3\rightarrow 4}$ $\displaystyle=$
$\displaystyle P_{2\rightarrow 1},$ (21) $\displaystyle P_{3\rightarrow 1}$
$\displaystyle=$ $\displaystyle P_{2\rightarrow 4},$ (22) $\displaystyle
P_{1\rightarrow 4}$ $\displaystyle=$ $\displaystyle P_{2\rightarrow
3}=P_{3\rightarrow 2}=P_{4\rightarrow 1}.$ (23)

Figure 5: Diabatic levels of the quantum dot models. (a) Dashed colored arrows
show the transitions from level 2 to level 1 (violet) and from level 3 to
level 4 (blue) with equal transition probabilities in both the symmetric and
the antisymmetric models. Due to the difference, $e_{1}\neq e_{2}$, there is
no obvious geometric symmetry between such transitions. (b) Two semiclassical
paths in the diabatic level diagram that contribute to the transition
probability from level 2 to level 1 (red and green arrows). Due to the quantum
interference, this probability is expected to show oscillatory dependence on
the area enclosed by the diabatic levels, $S=e_{1}e_{2}/b$.  Figure 6:
Numerically found time-dependence of transition probabilities $P_{2\rightarrow
1}$ (solid curves) and $P_{3\rightarrow 4}$ (dashed curves) in the
antisymmetric model for $\\{e_{1},e_{2}\\}$=$\\{1,3\\}$ (blue), $\\{1,2\\}$
(red). The remaining parameters: $b=2$, $g=2$, and $\gamma=1$. As both
$P_{2\rightarrow 1}$ and $P_{3\rightarrow 4}$ saturate at the same values at
$t=\pm\infty$, their intermediate time dynamics are different.

While for $e_{1}=e_{2}$, the relations (20)-(23) are consequences of trivial
discrete symmetries, for $e_{1}\neq e_{2}$ they are generally the results of
the integrability of the LZ-grid model. We also note that our analysis, and
hence relations (20)-(23) apply only to evolution from $t=-\infty$ to
$t=+\infty$, whereas they do not apply to the probabilities at intermediate
times, except for $e_{1}=e_{2}$, as we illustrate in Fig. 6. Hence, Eqs.
(20)-(23) are our first nontrivial prediction for experimental verification.

The integrability of time-dependent Hamiltonians is a type of quantum
symmetries that have not been studied experimentally previously. The
experiments on Landau-Stückelberg interferometry provide an opportunity to
detect the presence of such unusual quantum symmetries by measuring the state-
to-state transition probabilities in already available solid state and atomic
systems. The deviations from the exact predictions would mean the presence of
the terms beyond the standard Hamiltonian $(\ref{ham-2})$, which may emerge
either due to nonlinear time-dependence of the gate voltage or effects of the
magnetic field on the tunneling amplitudes.

In addition to the integrability conditions, there are six elements of the
transition probability matrix that are known exactly and explicitly due to the
no-go rule and the Brundobler-Elser formula be ; usuki ; nogo-LZ . For the
antisymmetric model they are

$\displaystyle P_{1\rightarrow 2}$ $\displaystyle=$ $\displaystyle
P_{4\rightarrow 3}=0,$ (24) $\displaystyle P_{1\rightarrow 1}$
$\displaystyle=$ $\displaystyle P_{2\rightarrow 2}=P_{3\rightarrow
3}=P_{4\rightarrow 4}=e^{-2\pi(g^{2}+\gamma^{2})/b},$ (25)

and for the symmetric model we should replace $\gamma\rightarrow g$.

In addition to relations (20)-(23), (24), and (25), we can only add the
unitarity of the evolution constraints:

$\sum_{n=1}^{4}P_{m\rightarrow n}=\sum_{n=1}^{4}P_{n\rightarrow m}=1,\,\forall
m\in\\{1,2,3,4\\}.$

It turns out that many of the latter relations are not independent after we
include the already mentioned relations. Thus, even having so many
constraints, the matrix of the transition probabilities has three unknown
independent parameters that have to be calculated separately.

Fortunately, the integrability leads to another simplification of the model’s
analysis. Namely, it was shown in quadr-LZ  that the integrability enables a
semiclassical approach for estimation of the transition probabilities in the
nearly-adiabatic limit. It was also noted in quadr-LZ  that the analytical
formulas that are obtained by this approach often provide a reasonable
approximation for the numerical solutions at arbitrary values of the
parameters.

Hence, in the following sections we apply this semiclassical approach to our
four-state models, $H_{1}$ and $H_{2}$. For simplicity, here we will restrict
ourselves only to the transition probability $P_{2\rightarrow
1}=P_{3\rightarrow 4}$. This transition probability is the most physically
interesting, first, because it is the only one that is needed to estimate the
probability to remain in the same quantum dot after the time-linear sweep of
the gate voltage. Indeed, all the other needed for this probabilities are
given by the exact expressions (24) and (25). Moreover, a simple analysis
shows that the probability $P_{2\rightarrow 1}$ should generally dominate over
$P_{1\rightarrow 1}=P_{2\rightarrow 2}$ in the adiabatic limit because the
latter become nonzero after two, rather than one in the case of
$P_{2\rightarrow 1}$, nonadiabatic overgap transitions.

Second, the transition from level 2 to level 1 is nontrivial even in the limit
of large separation of all level crossings because it is then influenced by
quantum interference of different evolution trajectories, as we show in Fig.
5(b). The dependence of the probabilities of such transitions on the relative
parameter values is understood poorly. Although numerical simulations of few-
state systems are easy, they show complex oscillatory behavior on the
parameters even in simplest and perturbative regimes kiselev-LZ , so it is
hard to see a general pattern for the role of different parameters. Hence, by
developing nonperturbative analytical description of $P_{3\rightarrow 4}$ in
models (17) and (18) we will obtain a useful insight into the interference
effects in the nonadiabatic regime.

##  IV Semiclassical solution for $P_{2\rightarrow 1}=P_{3\rightarrow 4}$ in
quantum dot models

In order to find the transition probability $P_{3\rightarrow 4}$ for models
(17) and (18) in the adiabatic limit, we take advantage of the
$\tau$-independence of the transition probability and make the slopes of the
tilted levels, 3 and 4, infinite by setting $\tau\rightarrow\infty$, as it was
done to study the three-state MLZ model in quadr-LZ . The tilted levels cross
both the levels 1 and 2 then at time moments:

$t_{\pm}=\pm e_{2}/b,$

where “$-$” is for the crossing of level 3 and “$+$” is for the crossing of
level 4.

Let us assume that level 3 is initially populated and we want to find the
probability to end up on level 4. Following quadr-LZ , we can employ the fact
that for high-slope crossing the characteristic time of the nonadiabatic
interactions is vanishing as $\delta t\sim 1/\sqrt{\tau}$, whereas the
nonadiabatic transitions between the diabatic states 3 and 4 take time that is
independent of $\tau$. Hence, we can separately treat the interactions of
levels 3 and 4 with level 1, then with each other, and then with level 2.
Specific details, however, depend on the models, which we will consider
separately.

###  IV.1 Symmetric model

First, we consider the Hamiltonian $H_{1}(t,\tau)$ in the limit
$\tau\rightarrow\infty$, and where $\tau$-dependence is according to (9). We
introduce the symmetric, $|+\rangle$, and the anti-symmetric, $|-\rangle$,
combinations of the diabatic states 1 and 2:

$|\pm\rangle=\frac{1}{\sqrt{2}}\left(|1\rangle\pm|2\rangle\right).$ (26)

As $\tau\rightarrow\infty$, both levels 3 and 4 couple only to the $|+\rangle$
state and the nonadiabatic transitions to and from $|+\rangle$ happen in the
direct vicinity of time moments $t_{-}$ and $t_{+}$. The corresponding
coupling is $g\sqrt{2\tau}$, and the slope difference between the diabatic
levels of $|+\rangle$ and either $|3\rangle$ or $|4\rangle$ is $b\tau$. The
probabilities of fast transitions from $|3\rangle$ to $|+\rangle$ near $t_{-}$
and from $|+\rangle$ to $|4\rangle$ near $t_{+}$ are given by the LZ formula
for two state transitions:

$P_{3\rightarrow+}=P_{+\rightarrow 4}=1-e^{-4\pi g^{2}/b}.$

In addition, during the time interval $t\in(t_{-},t_{+})$ the sates
$|+\rangle$ and $|-\rangle$ interact with each other, in particular, via the
virtual transitions trough $|1\rangle$ and $|2\rangle$, as explained in quadr-
LZ . Hence, the total transition probability from level 3 to level 4 can be
expressed via the product of the probabilities

$P_{3\rightarrow 4}=P_{3\rightarrow+}P^{++}P_{+\rightarrow 4}=\left(1-e^{-4\pi
g^{2}/b}\right)^{2}P^{++},$ (27)

where $P^{++}$ is the probability to remain in $|+\rangle$ after the dynamics
with a $2\times 2$ effective Hamiltonian that acts in the subspace of
$|\pm\rangle$ during the time interval $t\in(t_{-},t_{+})$.

Away from the points $t=t_{\pm}$, as $\tau\rightarrow\infty$, the effect of
virtual transitions to and from levels 1 and 2 can be calculated
perturbatively. Up to the zeroth order in $\tau$, such transitions lead to an
effective Hamiltonian in the subspace of states $|\pm\rangle$:

$h_{v}=\frac{-2g^{2}}{b}\left(\frac{1}{t-t_{-}}+\frac{1}{t-t_{+}}\right)|+\rangle\langle+|,$
(28)

where the index $v$ refers to the interactions appearing due to the virtual
transitions in the second order of perturbation series over $1/\tau^{1/2}$.
The higher order corrections due to such transitions in the limit
$\tau\rightarrow\infty$ vanish.

The other contribution to the effective two-state Hamiltonian arrives from the
splitting $e_{1}$ of levels 1 and 2, which mixes states $|\pm\rangle$:

$h_{e}=e_{1}\left(|-\rangle\langle+|+|+\rangle\langle-|\right).$ (29)

$P^{++}$ can be determined by solving the effective Schrödinger equation with
the Hamiltonian $H_{\rm eff}=h_{v}+h_{e}$. After rescaling time $t\rightarrow
te_{2}/b$, this equation in the basis of states $|\pm\rangle$ is given by

$i\frac{d}{dt}|\psi\rangle=\frac{1}{b}H_{\rm eff}^{s}(t)|\psi\rangle,\quad
t\in(-1,1),$ (30)

where

$H_{\rm
eff}^{s}(t)=\left(\begin{array}[]{cc}\frac{4g^{2}t}{t^{2}-1}&e_{1}e_{2}\\\
e_{1}e_{2}&0\end{array}\right).$ (31)

In order to find $P^{++}$, we should find the amplitude of evolution from
$|+\rangle$ as $t\rightarrow-1$ into $|+\rangle$ as $t\rightarrow+1$, and take
its absolute value squared. Note that the parameter $b$ plays the role in (30)
of an effective Planck constant.

Thus, the integrability reduces the problem to a two-state system with time-
dependent 2$\times$2 Hamiltonian $H_{\rm eff}^{s}$. It is possible to
reformulate equation (30) by changing to the variable $x\in(-\infty,+\infty)$,
where $t=\tanh(x)$, and $dt/(t^{2}-1)=dx$. This makes it standard for
application of the, so-called, Dykhne formula for the overgap transition
probability between two states approx-LZ5 ; approx-LZ0  in the adiabatic
limit, in which overgap transitions are suppressed exponentially. This formula
provides the corresponding slowest decaying exponent and its leading order
prefactor. We will use this formula without switching from $t$ to $x$ because
the direct application of the Dykhne formula to the evolution during
$t\in(-1,1)$ produces the same final result, and can be equally justified for
the evolution (30).

The difference of the eigenvalues of $H_{\rm eff}^{s}$ is found analytically:

$\Delta
E(t)=2\frac{\sqrt{(e_{1}e_{2})^{2}(t^{2}-1)^{2}+4g^{4}t^{2}}}{(1-t^{2})}.$
(32)

By equating this difference to zero we find the branching points:

$t_{1,2}=\left(\frac{-2+r^{2}\pm 2\sqrt{1-r^{2}}}{{}^{2}}\right)^{1/2},\quad
r\equiv\frac{(e_{1}e_{2})}{g^{2}}.$ (33)

where the square root convention in $(\ldots)^{1/2}$ is chosen so that ${\rm
Im}(t_{1,2})>0$, and $t_{1}$ corresponds to the “$+$” sign in (33). Depending
on whether the ratio $r$ is bigger or smaller than $1$, we found two different
types of the behavior.

Phase I: $r<1$, i.e., $e_{1}e_{2}<g^{2}$. In this case, $\sqrt{1-r^{2}}$ is
real, so the branching points $t_{1,2}$ are purely imaginary. The transition
probability $P^{++}$ can be estimated with the standard Dykhne formula as

$P^{++}\large=e^{-(2/b){\rm Im}\left[\int_{0}^{t_{1}}\Delta
E(t)\,dt\right]},\quad r<1,$ (34)

where the final integration point $t_{1}$ is the imaginary root in (33) that
is closer to the real time axis.

Phase II: $r>1$, i.e., $e_{1}e_{2}>g^{2}$. Here, the branching points have
both real and imaginary parts. Moreover, the imaginary parts are equal to each
other, so both the branch cuts are relevant because they correspond to
semiclassical evolution trajectories with comparable amplitudes.

The transition probability is given by a generalized Dykhne formula that sums
the amplitudes of both trajectories and only then takes its absolute value
squared:

$P^{++}=\Big{|}e^{-\frac{i}{b}\int_{0}^{t_{1}}\Delta
E(t)\,dt+i\phi_{g}}+e^{-\frac{i}{b}\int_{0}^{t_{2}}\Delta
E(t)\,dt}\Big{|}^{2},$ (35)

where $\phi_{g}$ is a geometric phase difference between the two trajectories.
It is of subdominant order $O(1)$ in comparison to the integrals in (35), and
for real Hamiltonians can take only discrete values $0$ or $\pi$. In appendix
A, we calculate this phase for both $H_{1}$ and $H_{2}$ and show that for the
symmetric model $\phi_{g}=0$.

The integrals in (34), (35) cannot be simplified anymore. In this form, the
probability $P^{++}$ is already much easier to calculate numerically than by
solving the Schrödinger equation numerically directly. In appendix B, we show
that the integrals can be additionally simplified for some choices of the
parameters. Such cases are useful for developing the intuition about the
magnitude of $P^{++}$ and its dependence on the parameters.

The result (35) is valid for $g^{2}/b\gg 1$. Naturally, at $g=0$ it makes an
unphysical prediction: $P^{++}=2$. In order to adjust Eq. (35), we note that
in our case

$t_{1}=-t_{2}^{*},$

so, we have

${\rm Re}\left[\frac{1}{b}\int_{0}^{t_{1}}\Delta E(t)\,dt\right]=-{\rm
Re}\left[\frac{1}{b}\int_{0}^{t_{2}}\Delta E(t)\,dt\right],$ ${\rm
Im}\left[\frac{1}{b}\int_{0}^{t_{1}}\Delta E(t)\,dt\right]={\rm
Im}\left[\frac{1}{b}\int_{0}^{t_{2}}\Delta E(t)\,dt\right],$

and the desired approximation that makes $P^{++}=1$ at $g=0$ is

$P^{++}\approx\frac{\cos^{2}\left[{\rm
Re}\left(\frac{1}{b}\int_{t_{0}}^{t_{1}}\Delta
E(t)\,dt\right)\right]}{\cosh^{2}\left[{\rm
Im}\left(\frac{1}{b}\int_{t_{0}}^{t_{1}}\Delta E(t)\,dt\right)\right]}.$ (36)

Finally, using (27), we obtain the desired approximation for the four-state
model $H_{1}$:

$P_{3\rightarrow 4}\approx\left(1-e^{-4\pi
g^{2}/b}\right)^{2}\frac{\cos^{2}\left[{\rm
Re}\left(\frac{1}{b}\int_{t_{0}}^{t_{1}}\Delta
E(t)\,dt\right)\right]}{\cosh^{2}\left[{\rm
Im}\left(\frac{1}{b}\int_{t_{0}}^{t_{1}}\Delta E(t)\,dt\right)\right]}.$ (37)

In Fig. 7, the analytical predictions (34) and (36) are compared to the
results obtained by solving the Schrödinger equation (30) numerically for
several values of the parameter combination $e_{1}e_{2}/g^{2}$. The analytical
result for the critical case, $e_{1}e_{2}=g^{2}$ (the red curve), was taken
from appendix Eq. (55). Phases I and II are clearly distinguishable in
numerical simulations: in phase I, $P^{++}$ decays monotonously with
increasing $1/b$, whereas in phase II, $P^{++}$ oscillates as a function of
$1/b$.

Note also that $P^{++}$ increases with increasing coupling strength $g$ in
phase I. For large $g$, the leading contribution to the sum, $P_{3\rightarrow
3}+P_{3\rightarrow 4}$, is dominated by the transition probability
$P_{3\rightarrow 4}$, and consequently $P^{++}$. Hence, our result agrees with
the asymptotic behavior that was found for the symmetric model $H_{1}$ in
malla-LZ,  in the limit $e_{1}e_{2}\ll g^{2}$.

Figure 7: Probability $P^{++}$ calculated for the symmetric model numerically
(plot markers), and the corresponding semiclassical prediction by Eqs. (34)
and (35) for $r=4$ (black), $r=\sqrt{2}$ (blue), $r=1$ (red), $r=4/9$ (green),
and $r=1/4$ (magenta).

###  IV.2 Antisymmetric model with the Hamiltonian $H_{2}$

For the Hamiltonian (18), there are two couplings, $g$ and $\gamma$. The
asymmetry between $g$ and $\gamma$ requires from us to introduce new linear
combinations of the states $|1\rangle$ and $|2\rangle$:

Figure 8: Probability $P^{+-}$ for the antisymmetric model, calculated
numerically (plot marker) for the evolution with the Hamiltonian (39). Solid
curves are the corresponding semiclassical predictions by Eq. (42) at
$e_{1}=e_{2}=2$, $g=2$, and $\gamma=0.8$ (black), $\gamma=1$ (blue),
$\gamma=1.5$ (red).

$|A^{+}\rangle=\frac{g|1\rangle+\gamma|2\rangle}{\sqrt{g^{2}+\gamma^{2}}},\quad|A^{-}\rangle=\frac{-\gamma|1\rangle+g|2\rangle}{\sqrt{g^{2}+\gamma^{2}}},$

such that $|3\rangle$ couples directly to $|A^{+}\rangle$ and $|4\rangle$
couples to $|A^{-}\rangle$. In the limit $\tau\rightarrow\infty$, both level 3
and level 4 cross the diabatic energies of $|1\rangle$ and $|2\rangle$ at time
moments, respectively, $t_{-}$ and $t_{+}$. Hence, the transition probability
$P_{3\rightarrow 4}$ is

$P_{3\rightarrow 4}=P_{3\rightarrow A^{+}}P^{+-}P_{A^{-}\rightarrow 4},$ (38)

where $P^{+-}$ is the probability of the transition from $|A^{+}\rangle$ as
$t\rightarrow t_{-}$ to $|A^{-}\rangle$ as $t\rightarrow t_{+}$, and the fast
nonadiabatic transitions are given by the standard LZ formula:

$P_{3\rightarrow A^{+}}=P_{A^{-}\rightarrow
4}=1-e^{-2\pi(g^{2}+\gamma^{2})/b}.$

In order to find $P^{+-}$, we follow the previous approach: we rescale time
and obtain the analogous to (30) effective Schrödinger equation for evolution
in the subspace $|A^{\pm}\rangle$ during time interval $t\in(-1,1)$ and the
effective 2$\times$2 Hamiltonian

$H^{A}_{\rm
eff}(t)=e_{1}e_{2}\left(\begin{pmatrix}\frac{r_{-}}{r_{+}}-\frac{r_{+}}{t+1}&-\sqrt{1-r_{-}^{2}/r_{+}^{2}}\\\
-\sqrt{1-r_{-}^{2}/r_{+}^{2}}&-\frac{r_{-}}{r_{+}}-\frac{r_{+}}{t-1}\end{pmatrix}\right),$
(39)

where

$r_{\pm}\equiv(g^{2}\pm\gamma^{2})/(e_{1}e_{2}).$

The corresponding distance between the adiabatic levels of $H_{\rm eff}^{A}$
is

$\Delta
E_{A}=\frac{2e_{1}e_{2}}{(1-t^{2})}\sqrt{(t^{2}-1)^{2}+2(t^{2}-1)r_{-}+r_{+}^{2}}.$
(40)

This expression has almost the same structure as Eq. (32) for the symmetric
model. Similarly, the transition from $|A^{+}\rangle$ to $|A^{-}\rangle$
requires a passage through the avoided crossing at $t=0$, and there are two
branching points $t_{1}$ and $t_{2}$ that are obtained by setting $\Delta
E_{A}=0$.

However, we found two qualitative differences of the antisymmetric model from
the symmetric one. First, the antisymmetric model does not have phase I, i.e.,
it always leads to the phase with two equally important branch cuts at time
points

$t_{1,2}=\left(1-r_{-}\pm i\sqrt{r_{+}^{2}-r_{-}^{2}}\right)^{1/2},$ (41)

where the convention for “$(\ldots)^{1/2}$” is to keep only the square roots
in the upper complex plane. Such roots in (41) satisfy the relation
$t_{1}=-t_{2}^{*}$ as in phase II of the symmetric model. Hence, the
antisymmetric model is always in the phase with the oscillatory behavior.
Second, we show in appendix A that the geometric phase for the antisymmetric
model is $\phi_{g}=\pi$. Hence,

$P^{+-}\approx\Big{|}e^{-\frac{i}{b}\int_{0}^{t_{1}}\Delta
E_{A}(t)\,dt}-e^{-\frac{i}{b}\int_{0}^{t_{2}}\Delta E_{A}(t)\,dt}\Big{|}^{2}.$
(42)

In Fig. 8, we provide the numerical check of (42), which confirms the
analytical prediction, including the effect of the topological phase
$\phi_{g}=\pi$.

At $g=0$, this formula predicts correctly that $P^{+-}=0$, and thus does not
need further adjustments. However, we found that a slightly better fit to
numerical simulations in the strongly nonadiabatic regime for the Hamiltonian
$H_{2}$ is given by

$P_{3\rightarrow 4}\approx
2\left(1-e^{-2\pi(g^{2}+\gamma^{2})/b}\right)^{2}\frac{\sin^{2}\left[\frac{1}{b}{\rm
Re}\int_{0}^{t_{1}}\Delta E_{A}(t)\,dt\right]}{\cosh\left[\frac{2}{b}{\rm
Im}\int_{0}^{t_{1}}\Delta E_{A}(t)\,dt\right]}.$ (43)

Formulas (37) and (43) become asymptotically exact in the adiabatic limit
$1/b\rightarrow\infty$. It is instructive to look also at how they perform in
the strongly nonadiabatic regime, for which they cannot be justified
rigorously. In Fig. 9, we show that, although small deviations from the
numerically exact predictions are generally visible, in a broad range of the
parameters both formulas perform quite well. They still correctly predict
multiple oscillations, including phases and amplitudes, even in the
nonadiabatic regime. Thus our semiclassical theory is sufficiently rigorous
for, e.g., planning the future experiments without resorting to exhaustive
numerical simulations.

Figure 9: Transition probabilities $P_{3\rightarrow 4}$ are shown for the
numerically exact simulations (plot markers) and campared to the analytical
approximations (solid lines) from Eqs. (37) and (43) for both the symmetric,
$H_{1}$, and antisymmetric, $H_{2}$, models. In all cases: $e_{1}=e_{2}=2$.
The coupling in the symmetric model is $g=1$ (black) and $g=3$ (blue), and in
the antisymmetric model: $g=2$, and $\gamma=0.8$ (red) and $\gamma=1$ (green).

##  V Discussion

The concept of integrability in MLZ theory originally was introduced to unify
various fully solvable models. Interestingly, there are also large classes of
systems that are formally integrable but cannot be fully solved. Thus, we
showed that all LZ-grid models satisfy the integrability conditions for time-
dependent Hamiltonians commute . As a result, the state-to-state transition
probabilities in all such models are invariant of parameter rescaling (16) but
this is not sufficient to determine these probabilities analytically.

Nevertheless, this symmetry is nontrivial and different from the set of the
previously known exact results that were derived for LZ-grids earlier usuki .
In this sense, the integrability commute  is akin in effect to a conservation
law for a time-independent Hamiltonian.

We are not aware of experimental studies of such symmetries in explicitly
time-dependent Hamiltonians, although these symmetries should be realizable in
a broad class of physical systems li-dynamic . Given that the Landau-
Stückelberg interferometry for a four-state antisymmetric model (18) has been
already demonstrated in quantum dots petta-4lz1 ; petta-4lz2 , the
integrability of LZ-grids can also be tested. Specifically, for the
Hamiltonian (18), we predict nontrivial but simple-looking relations between
the state-to-state transition probabilities (20)-(23) in the dynamics induced
by a linear potential chirp.

The experiment to observe such a nontrivial symmetry requires a modification
of the control and measurement protocol in petta-4lz1 . The standard Landau-
Stückelberg interferometry protocol induces periodic transitions through the
region with nonadiabatic transitions. It leads to spectacular multi-
dimensional plots for the state probability dependence on the potential
sweeping rate and the duration of the periodic control. However, such
measurements are influenced by all parameters of the scattering matrix for a
single passage through the nonadiabatic region. Hence, an analytical
interpretation of such patterns for multi-state quantum systems is usually
problematic due to the large number of relevant parameters.

Here, we propose to explore the multistate Landau-Zener transition
probabilities directly, for example, resetting the state and then measuring it
after each linear potential chirp in a coupled quantum dot system petta-4lz1 .
The advantage of such an experiment would be the direct measurement of the
characteristics, for which either exact or semiclassical analytical
predictions can be established. By comparing the experimental data to the
theory, one can then explore such phenomena as the probability oscillations
due to completely nonadiabatic interference effects that we exposed, as well
as to confirm the exact constraints that follow from the model’s integrablity.

On the side of mathematical physics, we found further confirmations to the
conjecture in quadr-LZ  that the time-dependent Hamiltonian integrability
often leads to a semiclassical approach to calculate the leading contributions
to the overgap transition amplitudes in multistate systems. This approach is
an opportunity to explore the driven quantum models with considerable
complexity. Although generally approximate, it produces reasonable
approximations even in the strongly nonadiabatic regime, in which it cannot be
rigorously justified.

Our semiclassical analysis reveals that even the simplest driven models can
show complex behavior, such as phase transitions, in which the adiabatic limit
plays a similar role to the thermodynamic limit in many-body systems. Thus, in
the symmetric model, the interfering semiclassical trajectories do not
necessarily lead to oscillations of measurable characteristics. Rather there
are two phases: one with oscillatory behavior and another one with a
monotoneous decay of the transition probability. We also found that models
with similar structure and similar oscillatory behavior, can belong to
different topological classes that are characterized by the topological phase
$\phi_{g}$ between the interfering semiclassical trajectories. The
integrability allows us to observe and quantify such phenomena in relatively
complex and strongly driven quantum systems.

##  Appendix A Topological $\pi$-phase

The $\pi$-phase between the amplitudes of different semiclassical trajectories
in complex time has been known in literature, see, e.g., in joye-prefactor .
In practice, however, this phase has been usually inferred numerically, and we
do not know of a published rigorous analytical theory for its calculation.
Therefore, here we provide such a theory, which exposes the topological nature
of the $\pi$-phase and provides a simple path for its calculation. Our models
with the Hamiltonians $H_{\rm eff}^{s}$ and $H_{\rm eff}^{A}$ are particularly
suitable for illustration because all calculations for them can be performed
analytically, and the results differ for the different models.

The geometric phase appears generally in quantum mechanics when evolution is
considered along a closed path in the parameter space. In the context of the
Dykhne formula, such a path can be found if we note that the phase difference
between the amplitudes of the trajectories, ${\bm{C}}^{-}$ and ${\bm{C}}^{+}$
in Fig. 10(a), that go through two different branching points, is the same as
the phase acquired during the evolution along a path
${\bm{C}}=\left({\bm{C}}^{-}\right)^{-1}{\bm{C}}^{+}$, that is, the cyclic
trajectory on the Riemann surface that starts and ends at $t=0$ on the
original real time axis. ${\bm{C}}$ follows ${\bm{C}}^{+}$ along its direction
marked in Fig. 10(a), and then switches to ${\bm{C}}^{-}$ but follows it in
the opposite direction to what is marked in Fig. 10(a). In Fig. 10(b) we show
that this path is a knot that winds around the origin points of the branch
cuts.

Figure 10:  (a) The two trajectories, ${\bm{C}}^{-}$ and ${\bm{C}}^{+}$,
around the corresponding branch points $t_{1}$ and $t_{2}$, that make
comparable contributions to the semiclassical probability amplitude. (b) The
path ${\bm{C}}=({\bm{C}}^{-})^{-1}{\bm{C}}^{+}$ is topologically equivalent to
a knot that winds around the branching points $t_{1,2}$. (c) The integral over
the knot path in Eq. (48) is the sum of contributions from the infinitely
small circular contours that wind in opposite directions around the poles at
$t_{1,2}$.

Any real two-state time-dependent Hamiltonian can be written in the form

$H(t)=f(t){1}_{2}+Z(t)\sigma_{z}+X(t)\sigma_{x},$ (44)

where $f(t)$, $Z(t)$, and $X(t)$ are functions of time. Let us then define

$\sin\theta=\frac{X(t)}{\sqrt{X(t)^{2}+Z(t)^{2}}},\quad\cos\theta=\frac{Z(t)}{\sqrt{X(t)^{2}+Z(t)^{2}}}.\quad$
(45)

Along the cyclic path in complex plane, ${\bm{C}}$, an eigenstate of $H(t)$
can still be parametrized by an angle $\theta$:

$|u(\theta)\rangle=\left(\begin{array}[]{c}-\sin\theta/2\\\
\cos\theta/2\end{array}\right).$ (46)

The topological phase originates from the fact that the Hamiltonian is
periodic function of $\theta$, whereas the eigenstate vector (46) depends on
$\cos(\theta/2)$ and $\sin(\theta/2)$. Hence, if $\theta$ changes from $0$ to
$2\pi$ during an adiabatic evolution, the Hamiltonian returns to its initial
form at the end but the state vector can change sign:

$|u(2\pi)\rangle_{\theta(\bm{C})=2\pi}=-|u(0)\rangle.$ (47)

In other words, the state vector acquires a phase $\pi$. Generally, a periodic
adiabatic evolution of parameters changes the angle $\theta$ by $2\pi n$,
where $n$ is an arbitrary integer. The topological phase is then either $0$,
if $n$ is even, or $\pi$, if $n$ is odd. There cannot be another type of a
geometric phase in this case because $\langle u(\theta)|\frac{\partial
u(\theta)}{\partial\theta}\rangle=0$.

For a closed trajectory ${\bm{C}}$, which can be parametrized by time $t$, the
topological phase is given by

$\phi_{g}=\frac{1}{2}\int_{\bm{C}}\frac{d\theta}{dt}dt.$ (48)

Note that

$\frac{d\theta}{dt}=\frac{1}{\cos\theta}\frac{d\sin\theta}{dt}.$

For the symmetric model (30), we have

$Z=\frac{2g^{2}t}{(t^{2}-1)},\quad X=e_{1}e_{2},$

so,

$\frac{1}{2}\frac{d\theta}{dt}=\frac{r(t^{2}+1)}{4t^{2}+r^{2}(1-t^{2})^{2}}=\frac{(t^{2}+1)}{r(t^{2}-t_{1}^{2})(t^{2}-t_{2}^{2})},$
(49)

where $r=e_{1}e_{2}/g^{2}$, and $t_{1,2}$ are the two roots in the upper
complex half-plane, which are written in (33).

Note that the integrand in (49) has simple poles at $t_{1,2}$ rather than the
branching points. Hence, the integral over ${\bm{C}}$ is given by the
difference of the residues at these poles, as illustrated in Fig. 10(c):

$\phi_{g}=\pi
i\left[Res\left(\frac{d\theta}{dt}\right)_{t_{1}}-Res\left(\frac{d\theta}{dt}\right)_{t_{2}}\right],$
(50)

where $Res(\ldots)_{a}$ is the residue of the expression at a simple pole $a$,
and where the minus sign is because ${\bm{C}}$ winds around $t_{1}$ and
$t_{2}$ in opposite directions. Substituting the roots from (33) to (50), and
making sure that $t_{1}=-t_{2}^{*}$, we find for the effective two-state
system with the Hamiltonian $H_{\rm eff}^{s}$ that

$\phi_{g}=0.$ (51)

Analogously, for the effective Hamiltonian (39) that corresponds to the
antisymmetric model (18), we have

$\frac{d\theta}{dt}=\frac{-2\sqrt{r_{+}^{2}-r_{-}^{2}}t}{(t^{2}-t_{1}^{2})(t^{2}-t_{2}^{2})},$
(52)

where $t_{1,2}$ are given by (41), and $t_{1}=-t_{2}^{*}$. Substituting (52)
and (41) into (50), we find that for the effective Hamiltonian (39) the
topological phase is

$\phi_{g}=\pi.$ (53)

Thus, despite similar oscillatory behavior, the symmetric and antisymmetric
models are characterized by the different values of topological phase
$\phi_{g}$. In this sense the models belong to different topological classes.

##  Appendix B Analytically simple special cases

###  B.1 Symmetric model: the critical point at $e_{1}e_{2}=g^{2}$

The condition $e_{1}e_{2}=g^{2}$ marks the phase transition point between
phase I and phase II. In this special case the two branching points are
degenerate and purely imaginary, i.e., $t_{1}=t_{2}=i$. The adiabatic energy
difference simplifies to

$\Delta E=2g^{2}\frac{t^{2}+1}{1-t^{2}}.$

Substituting $t=ix$ into the integral in the Dykhne formula (34) we find

$2{\rm Im}\left[\int_{0}^{t_{1}}\Delta
E\,dt\right]=2g^{2}\int_{0}^{1}\frac{1-x^{2}}{1+x^{2}}\,dx=2(\pi-2).$

Following joye-prefactor , we should set the exponential prefactor in this
special case to be 2 rather than 1 because the behavior near the zero is
$\Delta E\sim(t-t_{1})$ rather than the typical $\Delta E\sim(t-t_{1})^{1/2}$.
This leads to the transition probability

$P^{++}_{e_{1}e_{2}=g^{2}}=2e^{-2(\pi-2)g^{2}/b}.$ (54)

This formula is valid only asymptotically in the limit $g^{2}/b\gg 1$. Without
affecting the behavior in this domain, we can also write

$P^{++}_{e_{1}e_{2}=g^{2}}\approx\frac{1}{\cosh(2(\pi-2)g^{2}/b)},$ (55)

which is still generally approximate but has right value, $1$, at $g=0$.
Finally, the semiclassical formula for the transition probability in the
original 4-state model for this special case is

$P_{3\rightarrow 4}^{e_{1}e_{2}=g^{2}}\approx\frac{\left(1-e^{-4\pi
g^{2}/b}\right)^{2}}{\cosh(2(\pi-2)g^{2}/b)}.$ (56)

###  B.2 Symmetric model: the case with oscillations at
$e_{1}e_{2}=g^{2}\sqrt{2}$

Substituting $e_{1}e_{2}=g^{2}\sqrt{2}$, to Eq. (32), we find

$\Delta E=2\sqrt{2}g^{2}\frac{\sqrt{t^{4}+1}}{1-t^{2}}.$

The two branch points in the upper plane are then given by

$t_{1}=e^{i3\pi/4},\quad t_{2}=e^{i\pi/4}.$

To calculate the integral, we set

$t=e^{i\pi/4}s,$

which leads to:

$\frac{1}{b}\int_{0}^{t_{2}}dt\Delta
E=\frac{2\sqrt{2}e^{i\pi/4}g^{2}}{b}\int_{0}^{1}ds\,\frac{\sqrt{1-s^{4}}}{1-is^{2}}=\frac{g^{2}(\xi+i\zeta)}{b},$

where the analytic expressions for $\xi$ and $\zeta$ can be written in terms
of the known special functions but the expressions are a bit lengthy. We just
give their numerical values, which can be found to very high precision:

$\xi=1.19814\ldots,\quad\zeta=1.9434\ldots.$

The analogous calculation leads to

$\frac{1}{b}\int_{0}^{t_{1}}dt\Delta E=\frac{g^{2}(-\xi+i\zeta)}{b},$

and then to

$P^{++}\approx\frac{\cos[\xi g^{2}/b]^{2}}{\cosh[\zeta g^{2}/b]^{2}},$

which provides the leading exponent for small $b$ and is $1$ at $g=0$.

###  B.3 Antisymmetric model: the case with $e_{1}e_{2}=g^{2}-\gamma^{2}$

In the antisymmetric model there is also a special case

$e_{1}e_{2}=g^{2}-\gamma^{2},\quad r_{-}=1,$

such that the integrals can be computed explicitly. Substituting $r_{-}=1$ in
(40), we find the eigenvalue difference

$\Delta E_{A}=\frac{2e_{1}e_{2}}{(1-t^{2})}\sqrt{t^{4}-1+r_{+}^{2}},$

and the branch points in the upper half plane are given by

$t_{1}=|r_{+}^{2}-1|^{1/4}e^{i3\pi/4},\quad
t_{2}=|r_{+}^{2}-1|^{1/4}e^{i\pi/4}.$

In order to calculate the integral we set

$t=|r_{+}^{2}-1|^{1/4}e^{i\pi/4}s,$

and find

$\frac{1}{b}\int_{0}^{t_{2}}dt\Delta
E_{A}=\frac{2e_{1}e_{2}|r_{+}^{2}-1|^{3/4}e^{i\pi/4}}{b}\int_{0}^{1}ds\,\frac{\sqrt{1-s^{4}}}{1-i|r_{+}^{2}-1|^{1/2}s^{2}}=\frac{e_{1}e_{2}}{b}\\{F[|1-r_{+}^{2}|]+iG[|1-r_{+}^{2}|]\\},$

where analytical expressions for $F[x]$ and $G[x]$ can be written in terms of
the known special functions, which can be calculated with high precision.
Similarly, we evaluate the integral for the branching point $t_{1}$:

$\frac{1}{b}\int_{0}^{t_{1}}dt\Delta
E_{A}=\frac{e_{1}e_{2}}{b}\\{-F[|1-r_{+}^{2}|]+iG[|1-r_{+}^{2}|]\\}.$

Combining the two contributions, and taking into account the topological phase
$\phi_{g}=\pi$, we finally find:

$P^{+-}\approx\frac{\sin[e_{1}e_{2}F[|1-r_{+}^{2}|]/b]^{2}}{\cosh[e_{1}e_{2}G[|1-r_{+}^{2}|]/b]^{2}}.$

## Acknowledgements

This work was supported by the U.S. Department of Energy, Office of Science,
Basic Energy Sciences, Materials Sciences and Engineering Division, Condensed
Matter Theory Program. R.K.M. was supported by LANL LDRD program under project
number 20190574ECR and Center for Nonlinear Studies under project number
20190495CR.

## References

  * (1) S. Brundobler, and V. Elser. “S-matrix for generalized Landau-Zener problem”, J. Phys. A 26, 1211 (1993).
  * (2) D. A. Harmin and P. N. Price. Incoherent time evolution on a grid of Landau-Zener anticrossings. Phys. Rev. A 49, 1933 (1994).
  * (3) D. A. Harmin. Coherent time evolution on a grid of Landau-Zener anticrossings. Phys. Rev. A 56, 232 (1997).
  * (4) D. Bouwmeester, I. Marzoli, G. P. Karman, W. Schleich, and J. P. Woerdman. Optical Galton board. Phys. Rev. A 61, 013410 (1999).
  * (5) T. Usuki. Theoretical study of Landau-Zener tunneling at the M+N level crossing. Phys. Rev. B 56, 13360 (1997).
  * (6) N. A. Sinitsyn and N. Prokof’ev. Nuclear spin bath effects on Landau-Zener transitions in nanomagnets. Phys. Rev. B 67, 134403 (2003)
  * (7) K. Saito, M. Wubs, S. Kohler, Y. Kayanuma, and P. Hanggi. Dissipative Landau-Zener transitions of a qubit: bath-specific and universal behavior. Phys. Rev. B 75, 214308 (2007).
  * (8) S. Ashhab. Landau-Zener transitions in a two-level system coupled to a finite-temperature harmonic oscillator. Phys. Rev. A 90, 062120 (2014).
  * (9) D. A. Garanin, R. Neb, and R. Schilling. Effect of environmental spins on Landau-Zener transitions. Phys. Rev. B 78, 094405 (2008).
  * (10) M. Werther, F. Grossmann, Z. Huang, and Y Zhao. Davydov-Ansatz for Landau-Zener-Stueckelberg-Majorana transitions in an environment: Tuning the survival probability via number state excitation. J. Chem. Phys. 150, 234109 (2019).
  * (11) F. Bello, N. Kongsuwan, J. F. Donegan, and O. Hess. Controlled Cavity-Free, Single-Photon Emission and Bipartite Entanglement of Near-Field-Excited Quantum Emitters. Nano Lett. 20, 8, 5830 (2020).
  * (12) M.Kervinen, J. E. Ramírez-Munoz, A. Valimaa, and M. A. Sillanpaa. Landau-Zener-Stückelberg interference in a multimode electromechanical system in the quantum regime. Phys. Rev. Lett. 123, 240401 (2019).
  * (13) Z. Sun, J. Ma, X. Wang, and F. Nori. Photon-assisted Landau-Zener transition: Role of coherent superposition states. Phys. Rev. A 86, 012107 (2012).
  * (14) P. Y. Wen, et al.. Landau-Zener-Stückelberg-Majorana interferometry of a superconducting qubit in front of a mirror. Phys. Rev. B 102, 075448 (2020).
  * (15) R. K. Malla and M. E. Raikh. Landau-Zener transition in a two-level system coupled to a single highly excited oscillator. Phys. Rev. B 97, 035428 (2018)
  * (16) V. N. Ostrovsky, M. V. Volkov, J. P. Hansen, and S. Selsto. Four-state nonstationary models in multistate Landau-Zener theory. Phys. Rev. B 75, 014441 (2007)
  * (17) V. A. Yurovsky and A. Ben-Reuven. Curve crossing in linear potential grids: The quasidegeneracy approximation. Phys. Rev. A 63, 043404 – Published 6 March 2001
  * (18) N. A. Sinitsyn, and V. Y. Chernyak. The quest for solvable multistate Landau-Zener models. J. Phys. A: Math. Theor. 50, 255203 (2017).
  * (19) N. A. Sinitsyn, E. A. Yuzbashyan, V. Y. Chernyak, A. Patra, and C. Sun. Integrable time-dependent quantum Hamiltonians. Phys. Rev. Lett. 120, 190402 (2018).
  * (20) V. Y. Chernyak, and N. A. Sinitsyn. Integrability in the multistate Landau-Zener model with time-quadratic commuting operators. arXiv:2006.15144, submitted to J. Phys. A (2021).
  * (21) C. Sun, V. Y. Chernyak, A. Piryatinski, and N. A. Sinitsyn. Cooperative Light Emission in the Presence of Strong Inhomogeneous Broadening. Phys. Rev. Lett. 123, 123605 (2019).
  * (22) F. Li, V. Y. Chernyak, and N. A. Sinitsyn. Quantum annealing and thermalization: insights from integrability. Phys. Rev. Lett. 121, 190601 (2018).
  * (23) V. Y. Chernyak, N. A. Sinitsyn, and C. Sun. Dynamic spin localization and gamma-magnets. Phys. Rev. B 100, 224304 (2019).
  * (24) E. A. Yuzbashyan. Integrable time-dependent Hamiltonians, solvable Landau-Zener models and Gaudin magnets.  Ann. Phys. 392, 323 (2018).
  * (25) Although the LZ formula should be credited to four authors: Landau, Majorana, Zener, and Stückelberg, the interference of semiclassical trajectories, in the context of nonadiabatic transitions, was discussed only by Landau and Stückelberg.
  * (26) F. Ginzel, A. R. Mills, J. R. Petta, G. Burkard. Spin shuttling in a silicon double quantum dot. Phys. Rev. B 102, 195418 (2020).
  * (27) X. Mi, S. Kohler, and J. R. Petta. Landau-Zener interferometry of valley-orbit states in Si/SiGe double quantum dots. Phys. Rev. B 98, 161404(R)(2018).
  * (28) L. D. Faddeev, and L. A. Takhtajan. Hamiltonian Methods in the Theory of Solitons. (Springer, 1987).
  * (29) A. Patra, and E. Yuzbashyan. Quantum integrability in the multistate Landau-Zener problem. J. Phys. A 48, 245303 (2015).
  * (30) R. K. Malla, and M. E. Raikh. Loss of adiabaticity with increasing tunneling gap in nonintegrable multistate Landau-Zener models. Phys. Rev. B 96, 115437 (2017).
  * (31) N. A. Sinitsyn. Solvable four-state Landau-Zener model of two interacting qubits with path interference. Phys. Rev. B 92, 205431 (2015).
  * (32) F. Li, and N. A. Sinitsyn. Dynamic Symmetries and Quantum Nonadiabatic Transitions. Chem. Phys. 481, 28 (2016).
  * (33) N. A. Sinitsyn. Counterintuitive transitions in the multistate Landau-Zener problem with linear level crossings. J. Phys. A: Math. Gen. 37 10691 (2004).
  * (34) M. N. Kiselev, K. Kikoin, and M. B. Kenmoe. SU(3) Landau-Zener interferometry. EPL 104 57004 (2013).
  * (35) K. Fukushima, and T. Shimazaki. Lefschetz-thimble inspired analysis of the Dykhne-Davis-Pechukas method and an application for the Schwinger Mechanism. Ann. Phys. 415, 168111 (2020).
  * (36) J.‐T. Hwang, and P. Pechukas. The adiabatic theorem in the complex plane and the semiclassical calculation of nonadiabatic transition amplitudes. J. Chem. Phys. 67, 4640 (1977).
  * (37) A. Joye. Non-trivial prefactors in adiabatic transition probabilities induced by high-order complex degeneracies. J. Phys. A.: Math. Gen. 26 6517 (1993).