The Quest for Solvable Multistate Landau-Zener Models

TL;DR

This paper demonstrates that integrability conditions facilitate the discovery of new solvable multistate Landau-Zener models, revealing their abundance and extending their applicability to complex many-body systems with multiple level crossings.

Contribution

The authors refine integrability conditions and develop a computer-assisted method to find new solvable MLZ models, including those with multiple diabatic level crossings.

Findings

Integrability conditions enable efficient search for new solvable MLZ models.

Solvable models are numerous, spanning simple to complex many-body systems.

Extended the class of solvable models to include multiple level crossing points.

Abstract

Recently, integrability conditions (ICs) in mutistate Landau-Zener (MLZ) theory were proposed [1]. They describe common properties of all known solved systems with linearly time-dependent Hamiltonians. Here we show that ICs enable efficient computer assisted search for new solvable MLZ models that span complexity range from several interacting states to mesoscopic systems with many-body dynamics and combinatorially large phase space. This diversity suggests that nontrivial solvable MLZ models are numerous. In addition, we refine the formulation of ICs and extend the class of solvable systems to models with points of multiple diabatic level crossing.

Equations

$$i\frac{d\Psi }{dt}=\hat{H}(t)\Psi ,\hat{H}(t)=\hat{A}+\hat{B}t.$$

$$\hat{H}_{\rm DO}=\left(\begin{array}[]{ccc}t/2&g_{1}&g_{2}\\ g_{1}&e_{1}-t/2&0\\ g_{2}&0&e_{2}-t/2\end{array}\right),\quad e_{1}<e_{2},$$

$$\hat{H}_{\rm 2B}=\left(\begin{array}[]{cccccc}2e_{2}-t&0&0&\sqrt{2}g_{2}&0&0\\ 0&e_{1}+e_{2}-t&0&g_{1}&g_{2}&0\\ 0&0&2e_{1}-t&0&\sqrt{2}g_{1}&0\\ \sqrt{2}g_{2}&g_{1}&0&e_{2}&0&\sqrt{2}g_{2}\\ 0&g_{2}&\sqrt{2}g_{1}&0&e_{1}&\sqrt{2}g_{1}\\ 0&0&0&\sqrt{2}g_{2}&\sqrt{2}g_{1}&t\end{array}\right),\quad e_{1}<e_{2}.$$

$$\hat{H}=\frac{1}{2}t({\hat{a}}^{†}\hat{a}-{\hat{b}}^{†}\hat{b}-{\hat{c}}^{†}\hat{c})+e_1{\hat{b}}^{†}\hat{b}+e_2{\hat{c}}^{†}\hat{c}+\left(g_1{\hat{a}}^{†}\hat{b}+g_2{\hat{a}}^{†}\hat{c}+h.c.\right),$$

$$|1⟩$$

$$|4⟩$$

$$i\frac{d}{dt}\left(\begin{array}[]{c}\hat{a}\\ \hat{b}\\ \hat{c}\end{array}\right)=\hat{H}_{\rm DO}\left(\begin{array}[]{c}\hat{a}\\ \hat{b}\\ \hat{c}\end{array}\right),$$

$$\hat{H}_{\rm LZ}=\left(\begin{array}[]{cc}\beta_{1}t&g\\ g&\beta_{2}t\\ \end{array}\right).$$

$$\hat{S}_{\rm LZ}=\left(\begin{array}[]{cc}\sqrt{p}&\pm i\sqrt{1-p}\\ \pm i\sqrt{1-p}&\sqrt{p}\\ \end{array}\right),\quad p=e^{-2\pi g^{2}/|\beta_{1}-\beta_{2}|},$$

$$p_{1,2}=e^{-2\pi g_{1,2}^2},q_{1,2}=1-p_{1,2}.$$

$$\hat{P}_{1}=\left(\begin{array}[]{ccc}p_{1}&q_{1}&0\\ q_{1}&p_{1}&0\\ 0&0&1\\ \end{array}\right),\quad\hat{P}_{2}=\left(\begin{array}[]{ccc}p_{2}&0&q_{2}\\ 0&1&0\\ q_{2}&0&p_{2}\\ \end{array}\right).$$

$${\hat{P}}_{\mathrm{DO}}={\hat{P}}_2{\hat{P}}_1.$$

$${\hat{P}}_{2B}={\hat{P}}_4{\hat{P}}_3{\hat{P}}_2{\hat{P}}_1.$$

$$\hat{H}_{1}=\left(\begin{array}[]{ccc}-t&\sqrt{2}g_{1}&0\\ \sqrt{2}g_{1}&0&\sqrt{2}g_{1}\\ 0&\sqrt{2}g_{1}&t\\ \end{array}\right).$$

$$\hat{P}_{1}=\left(\begin{array}[]{ccc}p_{1}^{2}&2p_{1}q_{1}&q_{1}^{2}\\ 2p_{1}q_{1}&(1-2p_{1})^{2}&2p_{1}q_{1}\\ q_{1}^{2}&2p_{1}q_{1}&p_{1}^{2}\\ \end{array}\right),$$

$$\hat{P}_{1}=\left(\begin{array}[]{cccccc}1&0&0&0&0&0\\ 0&1&0&0&0&0\\ 0&0&p_{1}^{2}&0&2p_{1}q_{1}&q_{1}^{2}\\ 0&0&0&1&0&0\\ 0&0&2p_{1}q_{1}&0&(1-2p_{1})^{2}&2p_{1}q_{1}\\ 0&0&q_{1}^{2}&0&2p_{1}q_{1}&p_{1}^{2}\end{array}\right).$$

$$\hat{P}_{2}=\left(\begin{array}[]{cccccc}1&0&0&0&0&0\\ 0&p_{1}&0&q_{1}&0&0\\ 0&0&1&0&0&0\\ 0&q_{1}&0&p_{1}&0&0\\ 0&0&0&0&1&0\\ 0&0&0&0&0&1\end{array}\right).$$

$$\hat{P}_{3}=\left(\begin{array}[]{cccccc}1&0&0&0&0&0\\ 0&p_{2}&0&0&q_{2}&0\\ 0&0&1&0&0&0\\ 0&0&0&1&0&0\\ 0&q_{2}&0&0&p_{2}&0\\ 0&0&0&0&0&1\end{array}\right),\quad\hat{P}_{4}=\left(\begin{array}[]{cccccc}p_{2}^{2}&0&0&2p_{2}q_{2}&0&q_{2}^{2}\\ 0&1&0&0&0&0\\ 0&0&1&0&0&0\\ 2p_{2}q_{2}&0&0&(1-2p_{2})^{2}&0&2p_{2}q_{2}\\ 0&0&0&0&1&0\\ q_{2}^{2}&0&0&2p_{2}q_{2}&0&p_{2}^{2}\end{array}\right).$$

$$\hat{H}(t)=\left(\begin{array}[]{cccc}-b_{1}t+e&0&g&-\gamma\\ 0&b_{1}t+e&\gamma&g\\ g&\gamma&-b_{2}t&0\\ -\gamma&g&0&b_{2}t\\ \end{array}\right),$$

$$\hat{H}(t)=\left(\begin{array}[]{cccc}-b_{1}t+e+e_{2}&0&gy&-\gamma x\\ 0&b_{1}t+e-e_{2}&\gamma&g\\ gy&\gamma&bt+e_{1}&0\\ -\gamma x&g&0&b_{2}t\\ \end{array}\right),$$

$$b_1>b_2,\mathrm{and}b<b_2.$$

$$t_{12}=0,t_{34}=e_1/(b_2-b),$$

$$t_{13}=\frac{e-e_1}{b_1+b},t_{14}=\frac{e}{b_1+b_2},t_{23}=-\frac{e-e_1}{b_1-b},t_{24}=-\frac{e}{b_1-b_2}.$$

$$A_1={\int }_{t_{24}}^{t_{23}}(b_{1}t+e)dt+{\int }_{t_{23}}^{t_{13}}(bt+e_1)dt.$$

$$A_1=\frac{b_1}{2}(t_{23}^2-t_{24}^2)+\frac{b}{2}(t_{13}^2-t_{23}^2)+e(t_{23}-t_{24})+e_1(t_{13}-t_{23}),$$

$$A_2=\frac{b_2}{2}(t_{14}^2-t_{24}^2)-\frac{b_1}{2}(t_{13}^2-t_{14}^2)+e(t_{13}-t_{14}).$$

$$e_1^{\pm }=e\pm |e|\sqrt{\frac{b_1^2-b^2}{b_1^2-b_2^2}}.$$

$$E_{13}\equiv E_1-E_3=e-e_1,E_{14}\equiv E_1-E_4=e.$$

$$E_{32}^{\prime }\equiv E_3^{\prime }-E_2^{\prime }=-e-\frac{e_1(b_1-b_2)}{b_2-b},E_{31}^{\prime }\equiv E_3^{\prime }-E_1^{\prime }=-e+\frac{e_1(b_1+b_2)}{b_2-b}.$$

$$-\frac{\gamma xg}{E_{14}}+\frac{\gamma gy}{E_{13}}=0.$$

$$xy=E_{31}^{\prime }/E_{32}^{\prime },$$

Full text

I Introduction

The future of quantum technology will depend crucially on our ability to manipulate complex quantum systems by time-dependent fields. Achieving sufficient control of driven quantum matter is thus an important challenge, including to mathematical physics. For evolution of more than two interacting states, one of the difficulties here is the scarcity of solvable models that are described by the nonstationary Schrödinger equation with practically interesting explicitly time-dependent Hamiltonians. Most of known solvable driven models are either trivially reducible to independent harmonic oscillators and two-state systems maj ; multiparticle ; reducible ; yuzbashyan-LZ ; galitski or correspond to almost impossible for implementation driving protocols such as multistate shortcuts to adiabaticity shortcuts . This lack of nontrivial useful results strongly restricts our understanding of quantum dynamics in explicitly time-dependent fields, especially at the level of mesoscopic interacting systems with a combinatorially large accessible phase space.

Considerable progress on resolving this problem has been achieved within the MLZ theory, whose goal is to find the scattering $N\times N$ matrix $\hat{S}$ that describes evolution with a linearly time-dependent Hamiltonian be :

[TABLE]

Here, $\Psi$ is the state vector in the space of $N$ states; $\hat{A}$ and $\hat{B}$ are constant Hermitian $N\times N$ matrices. Linear time-dependence of parameters is relatively easy to implement in experiment, so MLZ models find numerous practical applications apps .

One can always choose the, so-called, diabatic basis in which the matrix $\hat{B}$ is diagonal, and if any pair of its elements are degenerate then the corresponding off-diagonal element of the matrix $\hat{A}$ is zero. In the diabatic basis, elements of the matrix $\hat{B}$ are called the slopes of diabatic levels, nonzero off-diagonal elements of the matrix $\hat{A}$ are called the coupling constants, and diagonal elements of the Hamiltonian are called diabatic energies.

An element $S_{nn^{\prime}}$ of the scattering matrix is the amplitude of the state $n$ at $t\rightarrow+\infty$, given that at $t\rightarrow-\infty$ the system was in the $n^{\prime}$-th diabatic state. In many applications, only the matrix $\hat{P}$, with elements $P_{nn^{\prime}}\equiv P_{n^{\prime}\rightarrow n}\equiv|S_{nn^{\prime}}|^{2}$ called transition probabilities, is needed. More precise definition of the scattering matrix and discussion of its general properties can be found in previous reviews be ; constraints . Historically, the first study of a model type (1) beyond $N=2$ was done by Majorana in maj .

We will call a particular choice of Eq. (1) solvable or integrable if one can write a compact analytical expression for elements of the transition probability matrix in terms of the well understood special functions of model parameters. Although the general solution of MLZ problem remains unknown, the number of exactly solved realizations of Eq. (1) has been growing quickly recently due to the discovery of integrability conditions (ICs) that identify models whose scattering matrices can be found exactly in the form of a product of specific elementary matrices six-LZ ; four-LZ ; cQED-LZ . Unfortunately, there is still no known algorithm to derive systems that satisfy ICs. Therefore the search for solvable models has required many trial-and-error efforts, so far.

In this article, we describe tricks that we found useful to alleviate the latter problem. While there is still no algorithm for detailed classification of solvable MLZ models, our goal is to show that there are systematic ways to search for new models that satisfy ICs. For illustrations, we derive and investigate several new solvable cases with different patterns of energy level crossings and different sizes of the phase space. We suggest analytical expressions for elements of the transition probability matrices in these models and provide results of rigorous numerical tests of our solutions.

We would like to stress that ICs in MLZ theory have not been proved mathematically rigorously yet. However, computational algorithms, described in Ref. cQED-LZ , can be used to test our predictions with high precision. For example, transition probabilities in models of the type (1) with up to $N=10$ states can be found numerically with up to three significant digits precision within a few minutes using a standard PC. Numerical algorithms allow “high throughput” tests of our predictions for a broad range of model parameters and different initial conditions. We considered such a test as confirming the theory only when all tested multiple choices of parameters produced numerical results that were indistinguishable from all our analytical predictions for a given model. We tested usually more than 100 different choices of constant parameters in the nonperturbative regime per model. Thus, although this article is about mathematical questions, it is experimental in spirit. The role of an experimental setup is played here by a computer and results of our numerical experiments are obtained in order to provide an input for more mathematically rigorous studies of MLZ theory in the future.

The plan of our article is as follows. In section II we summarize basic knowledge about integrability in MLZ theory and suggest refined version of ICs. In section III, we demonstrate our approach using a simple four-state MLZ Hamiltonian. In section IV, we show how one can generalize this approach to a higher dimensional case using a six-state model as an example, and in section V we show how similar intuition can be used to derive a new combinatorially complex but solvable MLZ model. We then summarize our findings in Discussion.

II MLZ models and integrability

In this section, we review basic established facts about integrability in MLZ theory using two previously solved models as examples. We will frequently use this section for references later. The additional goal is to refine definition of ICs.

II.1 Diabatic level diagram

It is convenient to illustrate parameters of any MLZ model in a graph with time-energy axes, as shown in Fig. 1(a) and Fig. 1(b). Lines of this graph are time-dependent diabatic levels, i.e. diagonal elements of the Hamiltonian in the diabatic basis as functions of time. Small black filled circles mark intersections of levels with nonzero pair-vise couplings. Integers on the left side of diabatic levels mark level indexes. On the right, levels are sometimes marked by analytic expressions for diabatic energies. If intersection of two levels is not specially marked then corresponding coupling between crossing diabatic levels is assumed to be zero, as the coupling between levels 2 and 6 in Fig. 1(b). If more than two levels intersect at one point, then nonzero couplings between different levels are additionally marked, as between levels 3 and 5 or 5 and 6 in Fig. 1(b).

Direct couplings between parallel levels are always considered zero, such as between levels 2 and 3 in Fig. 1(a); this is included in the definition of the diabatic basis. Using these rules, we can read off the Hamiltonian from such a diagram. For example, for Fig. 1(a) the Hamiltonian is a 3$\times$3 matrix

[TABLE]

and the Hamiltonian of the model in Fig. 1(b) is a 6$\times$6 matrix:

[TABLE]

Model (2) is known as the 3-state Demkov-Osherov (DO) model be ; constraints ; do . It is one of the first MLZ models that were solved by methods of complex analysis do . Model (3) is also solvable as it is a bosonic extension of the three-state DO model in the sense that it is realized by a quadratic Hamiltonian of three interacting bosonic modes:

[TABLE]

where $\hat{a}$, $\hat{b}$ and $\hat{c}$ are bosonic annihilation operators. This Hamiltonian conserves the number of bosons, so if we look at the two-boson sector of this model then in the basis

[TABLE]

where $|0\rangle$ is the vacuum state, the Hamiltonian (4) has the matrix form (3). On the other hand, if we switch to the Heisenberg picture, then the quadratic Hamiltonian (4) leads to evolution of operators

[TABLE]

where $\hat{H}_{\rm DO}$ is the same as in (2). Let $S_{ij}$, $i,j\in 1,2,3$, be known scattering matrix elements of the DO model. Then we also have relations among operators, such as $\hat{a}(+\infty)=S_{11}\hat{a}(-\infty)+S_{12}\hat{b}(-\infty)+S_{13}\hat{c}(-\infty)$ e.t.c.. One can use them to connect correlators of an arbitrary number of bosonic operators at $t=+\infty$ with correlators at $t=-\infty$, from which one can derive transition probabilities in the model (3). This procedure was designed in multiparticle . It was applied to solve the model (3) in jipdi .

II.2 Integrability Conditions

The search for new solvable MLZ models by methods of complex analysis and multiparticle generalizations, as described above, has achieved only limited success. Recently, a much bigger class of solvable MLZ systems was discovered with a different approach based on ICs six-LZ ; four-LZ ; cQED-LZ .

The simplest way to introduce ICs is by using the diabatic level diagram. Let us say that a path in this diagram is closed if it goes along diabatic levels to produce a closed loop such that switching levels along this path is allowed only at level crossings with nonzero couplings. An example of such a path in Fig. 1(b) is the one that starts at the crossing of levels 3, 5 and 6, then goes along level 6 forward in time, switches to level 4 at crossing with levels 4 and 1, then goes backward in time to level 2; then switches to level 5, and then goes backward in time to return to the original crossing point with levels 3 and 6.

We define the area inside such a closed path as the sum of areas of enclosed plaquettes in the diabatic level diagram counting clockwise and counterclockwise enclosed areas with opposite signs. We also introduce the notion of a Hamiltonian projected to a subset of levels. This Hamiltonian is obtained by setting couplings of given levels with all other states to zero and restricting the considered phase space only to the resulting evolution within the given subset of levels. ICs then say that for the full integrability of a MLZ model the following two conditions should be satisfied:

(i) All closed paths of the graph in the diabatic level diagram should enclose zero areas.

(ii) If two or several diabatic levels intersect at one point, and if there is no coupling path that connects two diabatic states within the Hamiltonian projected to this subspace of levels, then there is a corresponding pair of adiabatic energies of the full Hamiltonian that experiences exact crossing at all sufficiently small but finite values of coupling parameters.

For example, the DO model satisfies condition (i) trivially because the graph in Fig. 1(a) is a simple tree so any closed path on it has zero area naturally. In this model, there are also no diabatic level crossing points without direct coupling between levels, so condition (ii) is also satisfied trivially.

The case of the model in Fig. 1(b) is less trivial. It is easy to check that the triangle enclosed by levels 2, 4 and 6 is equal to the triangle enclosed by levels 2, 5 and 6, and that there is a single independent closed path in this level diagram that has to enclose one of these triangles clockwise and another one counterclockwise. So, the area of this path is zero, which proves condition (i).

We will not prove validity of (ii) completely analytically because it is generally a complex task while the lack of minigaps at eigenvalue crossing points can be visualized and then tested numerically very quickly with an approach that was described in Ref. cQED-LZ . For example, to make the presence of exact level crossings visually obvious, we diagonalize the Hamiltonian (3) and plot its eigenvalues as function of time in Fig. 1(c). This figure shows three pairs of adiabatic level crossings that correspond to intersections of diabatic levels in Fig. 1(b) with zero pairwise couplings, i.e. the pairs (3,4), (2,6) and (1,5), which proves property (ii) for this model.

There are two more pairs of levels without direct couplings, namely (3,6) and (1,6) that, however, do not lead to exact crossings of adiabatic levels. This is because these levels cross simultaneously with other levels, respectively 5 and 4, to which they are coupled directly. So, if the Hamiltonian is projected to three such simultaneously crossing levels the phase space does not split into disjoined subspaces. Condition (ii) allows such crossings to happen without emergence of exact adiabatic level crossings.

Finally, we note that our definition of ICs (i-ii) is slightly different from the one that was originally suggested in six-LZ . Condition (i) here is more restrictive than just assumption of equality between dynamic phases of interfering semiclassical trajectories. This would make no difference for all previously found solvable models. However, in appendix A, we discuss an example of a model that breaks the condition (i) but does not break its counterpart in Ref. six-LZ . Our numerical tests showed that, at least by the standard semiclassical ansatz, this model is no longer solvable. So the new version of condition (i) should be preferred. In appendix B we also point that the new definition of the first integrability condition raises interesting analogies with properties of WKB quantized classically integrable multidimensional systems. On the other hand, our definition of condition (ii) is less restrictive than in Ref. six-LZ because now we allow level diagrams with more than two diabatic levels to cross at one point, as it is the case in Fig. 1(b).

II.3 Solution in the form of the semiclassical ansatz

If ICs are satisfied then one can obtain scattering amplitudes of the model in the form of a product of simple matrices. This form of the solution is called the matrix product ansatz or simply the semiclassical ansatz. The name follows from the fact that this ansatz coincides with prediction of the semiclassical approximation in which nonadiabatic processes near different crossing points happen independently of each other. Assuming such real values of all couplings, semiclassical ansatz is generated as follows six-LZ :

1. We list all level crossing points of the diabatic level diagram in chronological order of their appearance.

2. For each crossing point, we identify the projected Hamiltonian that is obtained by keeping only diabatic levels that cross at this point and assuming that couplings to all other levels are zero.

3. For each projected Hamiltonian, we determine the projected scattering matrix, i.e. the scattering matrix of MLZ model described by the projected Hamiltonian. We augment this scattering matrix then to the full phase space of the original model by assuming that other states simply do not experience interactions.

4. The scattering matrix of the full model is then the chronologically ordered product of projected scattering matrices at all encountered crossing points.

5. Final transition probabilities are obtained by taking the absolute value squared of elements of this scattering matrix.

An equivalent formulation of the semiclassical ansatz can be done in terms of the semiclassical trajectories six-LZ that are the paths that go forward in time along the diabatic levels, so that they start and end at, respectively, initially and finally populated levels. Each time such a path goes through a crossing point, one should prescribe a specific amplitude factor that this point contributes to the full trajectory amplitude as described in six-LZ . The final amplitude of the transition is then the sum of amplitudes of all semiclassical trajectories connecting the given pair of states.

Here we note that all parameters of scattering matrices at individual crossing points are often not needed to be known. If there are no pairs of semiclassical trajectories that interfere with each other in the full level diagram, then one can skip the step 4) and write the final transition probability matrix as the product of individual transition probability matrices at each crossing point.

It was also found that, in the case when interference of trajectories matters, only trivial phase factors of elements of projected scattering matrices should be included in order to determine final transition probabilities six-LZ . For example, if only two levels intersect at a given crossing point then the projected Hamiltonian at this crossing is the one of the two-state Landau-Zener model book-LZ with some coupling $g$ and crossing level slopes $\beta_{1}$ and $\beta_{2}$:

[TABLE]

In all solvable models with a finite number of states, nonzero coupling constants can always be chosen real but not always positive. The simplified scattering matrix that corresponds to the projected two-state Hamiltonian (7) is then known six-LZ . It depends on the sign of $g$:

[TABLE]

where sign ($\pm$) is the same as the sign of the coupling constant $g$. Note that the matrix (8) does not include dynamic coupling-dependent phases of the full scattering matrix of the Landau-Zener model. The reason is that in all exactly solved models such phases factorize in the final result and do not influence final transition probabilities. Therefore, in what follows, we will use only such reduced versions of scattering matrices in our intermediate calculations, and verify validity of this choice at the end numerically.

II.4 Examples of semiclassical ansatz solutions

Let us consider the DO model (2) as an example. Since there is no interference between different paths connecting initial and final states, we can work directly with transition probability matrices that describe each crossing point. Let us define parameters

[TABLE]

There are two crossing points with nonzero couplings. In each case, the Hamiltonians projected to the crossing points are just the two state ones of the form (7), so we can introduce matrices of transition probabilities for evolution with such projected Hamiltonians:

[TABLE]

The full transition probability of the three-state DO model is then the product of these matrices taken in chronological order:

[TABLE]

Let us now turn to the model (3). As in any solvable model, its scattering matrix factorizes in chronologically ordered contributions from all diabatic level crossing points with nonzero couplings. There is again no path interference in Fig. 1(b) of the model (3), so the full transition probability matrix can be directly written as a product of four transition probability matrices describing four crossing points:

[TABLE]

The matrix $\hat{P}_{1}$ here describes transitions at simultaneous crossing of levels 3, 5, and 6. Projected to the subspace of these three levels Hamiltonian is a 3$\times$3 matrix. After the shift of the crossing point to $t=0$, this Hamiltonian reads

[TABLE]

It describes the model of a spin-1 in a linearly time-dependent magnetic field. The transition probability matrix for the model (13) is known since the time of Majorana maj :

[TABLE]

or in the basis of all six states of the original model

[TABLE]

Similarly, $\hat{P}_{2}$ is a unit matrix except the 2$\times$2 block that corresponds to transitions between levels 2 and 4. Corresponding probabilities are given by the standard two-level Landau-Zener formula:

[TABLE]

Matrices $\hat{P}_{3}$ and $\hat{P}_{4}$ are constructed analogously:

[TABLE]

Such an approach to obtain the full transition probability matrix (12) of the model (3) is much more straightforward than dealing with correlators of bosonic operators and relating them to transition probabilities.

III General 4-state system satisfying integrability conditions

III.1 Distorting previously solved model

In previous section, we showed that if a model satisfies ICs then derivation of its transition probability matrix is straightforward. However, the main challenge usually is to identify systems satisfying these conditions. Here we propose an approach that is based on our new conjecture stating that already solved MLZ models are actually special limits of some more complex solvable models.

For example, we know that solution for spin-1 in a linearly time-dependent field maj was later found to be a special case of the three-state bow-tie model hioe , which in turn was generalized later to the solvable multistate version bow-tie , and eventually to the solvable generalized bow-tie model gbow-tie . So, it is reasonable to expect that other known exact results in MLZ theory have similar generalizations.

The utility of this conjecture is that an already solved system provides the graph of diabatic levels with required topology. Here we say that graphs of diabatic levels, or corresponding models, have the same topology if the graphs have arbitrarily different parameters (couplings, level slopes) but the same set of level crossing points and types of level crossings. The latter means that two conditions are satisfied. First, levels with the same indexes in both graphs cross in the same points. Second, choices of parameters of corresponding Hamiltonians should be such that points leading to exact eigenvalue crossings should be the same in both graphs. One consequence of the latter rule is that, for a pairwise level crossing, if the direct coupling between corresponding two levels is zero in one model it should be also zero in the topologically equivalent model.

We can use the graph of the previously solved MLZ model to design a more complex solvable model by distorting this graph, which means assuming a more general choice of level slopes and coupling parameters while preserving the model’s topology. Generally, assuming arbitrary values of level slopes and couplings would violate ICs, however, for an available graph, we can reimpose such conditions in the form of relations within the new set of parameters. We are going to show that the resulting model is generally richer than the original one in the sense of the larger number of independent parameters. Moreover, it usually contains the original model as a special case.

As the first example, let us consider a simple four-state model that was shown to satisfy ICs and solved in four-LZ . It has the Hamiltonian

[TABLE]

and its level diagram is plotted in Fig. 2(a). Exact validity of the semiclassical ansatz for this model was rigrously proved in constraints .

Let us now deform this graph as shown in Fig. 2(b) so that the new graph still has two crossings of diabatic levels without direct couplings. All nonzero couplings can now be considered different, which leads to the following Hamiltonian:

[TABLE]

with some constants $b_{1},b_{2},b,e,e_{1},e_{2},x,y,g,\gamma$. Here we used the gauge freedom be so that slopes of levels 1 and 2 differ only by sign, and we set constant diagonal element of level 4 to zero.

Parameter $e_{2}$ can be set to zero by time shift $t\rightarrow t+e_{2}/b_{1}$ that does not affect transition probabilities for evolution from $-\infty$ to $+\infty$. After this, we can make gauge transformation of diabatic state amplitudes $a_{i}\rightarrow e^{-itb_{2}e_{2}/b_{1}}a_{i}$, $i=1,\ldots,4$, which results in $e_{2}=0$ at redefined parameters $e\rightarrow e-b_{2}e_{2}/b_{1}$ and $e_{1}\rightarrow e_{1}+(b-b_{2})e_{2}/b_{1}$. We also note that the gauge freedom to choose zero energy point and the arbitrariness of choosing indexes of states allow us to assume that

[TABLE]

Otherwise, we can just rename indexes, followed by proper renaming of constant parameters and obtain the same Hamiltonian at desired relations (20). Therefore, in what follows, we will assume that $e_{2}=0$ with constraints in (20), unless specially stated.

Let $t_{ij}$ be the time moment of diabatic crossing between levels $i$ and $j$. For crossings without direct couplings we have:

[TABLE]

and for time moments of avoided crossings (i.e. moments of crossing of corresponding diabatic energy levels):

[TABLE]

Let us now consider the condition on that the area inside the closed path in the graph in Fig. 2(b) is zero. This area is the difference between areas swept over the time axis by two trajectories that connect crossing points at time moments $t_{24}$ and $t_{13}$. One such a trajectory starts at $t_{24}$, then stays on level 2 till time $t_{23}$, and then goes along level 3 till time $t_{13}$. The area that this trajectory sweeps over the time axis is

[TABLE]

The second trajectory starts at $t_{24}$, then stays on level 4 till time $t_{14}$, and then goes along level 1 till time $t_{13}$. Explicitly, for two trajectories we have then

[TABLE]

[TABLE]

Setting $A_{1}-A_{2}=0$ we find equation that parameters of the model satisfy in order to realize IC (i). In terms of $e_{1}$ this equation is quadratic with two roots:

[TABLE]

Equation (25) shows that only values of $b$ in the range $-b_{1}<b<b_{1}$ lead to physical real valued $e_{1}^{\pm}$.

Our parameters do not satisfy IC (ii) yet. Generally, an analytical proof that some crossing of diabatic levels corresponds to an exact crossing point of adiabatic energy levels is a complex problem. However, there are perturbative conditions that must be satisfied by parameters if there is an exact crossing point. We will derive such constraints analytically and then check numerically whether imposing them is sufficient to satisfy condition (ii).

Let $E_{k}$, $k=1,\ldots,4$ be diabatic energies of levels, respectively, $1,\ldots,4$ at time $t_{12}$. We also define energy gaps

[TABLE]

Similarly, let $E_{k}^{\prime}$ be energies of levels at time $t_{34}$. We have then

[TABLE]

Let us now treat all couplings as small. The necessary condition for the crossing point between levels 1 and 2 to be exact is that coupling between these diabatic states remains zero at 2nd order perturbation in $\gamma$ and $g$. This leads to the constraint

[TABLE]

Similarly, coupling between levels 3 and 4 should be zero at 2nd order perturbation. This gives the constraint

[TABLE]

which leads to

[TABLE]

Thus, in order to satisfy ICs, there are three constraints in (25) and (28) on parameters so that only parameter $b$ remains independent among the newly introduced ones. The original model in Fig. 2(a) is recovered at $b=-b_{2}$ but other values of $b$ produce models that have not encountered previously in MLZ theory.

Conditions (28) are necessary but not sufficient in order to guarantee the presence of exact crossing points near times $t_{12}$ and $t_{34}$. Therefore, we resort to numerical simulations for additional proof. Figures 3(a-b) and Figs. 4(a-b) show results of explicit numerical calculation of the adiabatic energy levels (i.e., eigevnalues of the Hamiltonian as functions of time $t$) at different values of parameters. It is visually clear that there are, indeed, pairs of exact eigenvalue crossings in each figure when conditions (25) and (28) are satisfied.

Depending on the sign in Eq. (25) and parameter $e$, one can generate different graph types. Two types correspond to models in which different semiclassical trajectories can interfere, while other two graph types do not feature path interference.

III.2 Phase with path interference

In Figs. 3(a-b), we show adiabatic energies of the Hamiltonain in cases when either (a) $e_{1}=e_{1}^{-}$ with $e>0$ or (b) $e_{1}=e_{1}^{+}$ with $e<0$. Geometries of these graphs are different by the relative positions of the exact crossing points. The common feature of both graphs is the possibility of more than one semiclassical trajectories that connect one diabatic state at $t=-\infty$ with another diabatic state at $t=+\infty$. In both cases we find

[TABLE]

with $x$ and $y$ having the same sign. Let us introduce

[TABLE]

Semiclassical ansatz calculations with parameters (29) are analogous to calculations in Ref. four-LZ that were done for the original model from Fig. 2(a). They lead to the following simple form of the transition probability matrix:

[TABLE]

[TABLE]

Figures 3(c-d) compare solution (31)-(32) with results of direct numerical simulations at different values of parameter $b$. The agreement between analytical predictions and numerics is perfect for all values of parameters. This leaves no doubt that semiclasscial ansatz prediction here is not an approximation but rather the exact result although, certainly, such an agreement with numerics cannot be accepted as a mathematically rigorous proof. Behavior of parameters $x$ and $y$ is provided in Fig. 3(e) in order to show that, in our numerical tests, these parameters varied substantially. For generality, we assumed in simulations that $e_{2}\neq 0$, so we also show Fig. 3(f) confirming that varying $e_{2}$ does not influence transition probabilities, as long as we adjust $e_{1}$ properly to satisfy integrability conditions.

Transition probabilities in Figs. 3(c-d) show discontinuous behavior at $b=b_{2}$. This is not surprising because at $b=b_{1}$ we have $e_{1}=0$, so diabatic levels 3 and 2 become degenerate not at a single time moment but during all the time. Transition probabilities in MLZ models are known to behave generally discontinuously when parameters pass through values that create such permanent degeneracies do .

Interestingly, solution (31)-(32) for the model in Fig. 2(b) has exactly the same form as the solution for the original model in Fig. 2(a) with only minor generalization that allows free choice of parameter $b$. This hints to possible symmetry that relates arbitrary values of $b$ to the case with $b=-b_{2}$, which was rigorously proved in constraints .

III.3 Phase without path interference

Let us now turn to phases such that (a) $e_{1}=e_{1}^{-}$ with $e<0$ or (b) $e_{1}=e_{1}^{+}$ with $e>0$ (see Eq. (25)). Repeating perturbative calculations, we find that in this case parameters $x$ and $y$ have opposite signs:

[TABLE]

i.e, negative couplings appear in pairs.

In Figs. 4(a-b) we show possible patterns of adiabatic energy levels and demonstrate two exact crossing points. Apparently, there is no path interference in this case. In both cases we find a solution

[TABLE]

Results of numerical tests of this prediction are shown in Figs. 4(c-d). Again, we find perfect agreement between the semiclassical ansatz and numerical predictions. In Figs. 4(e,f) we show, respectively, adiabatic levels at $b=b_{2}$ and dependence of parameters $x$, $y$, $e_{1}$ on $b$ that we used to satisfy ICs.

Several properties of solution (34) are to be mentioned. First, up to exchange of positions of some raws and columns, this transition probability matrix is a direct product of 2$\times$2 matrices

[TABLE]

which are transition probability matrices of two two-state Landau-Zener models. In fact, the special case with $b=-b_{2}$ corresponds to the solvable model from Ref. multiparticle . This model is, indeed, constructed as a direct product of two two-state Hamiltonians. This fact also suggests that the case with arbitrary $b$ may be related to the direct product of two decoupled two-state Landau-Zener systems.

Second interesting property is the absence of discontinuity at $b=b_{2}$. This follows from the fact that in this case we have $e_{1}=2e\neq 0$, i.e., there is no permanent degeneracy of two levels. At this point, we have $x=-y=\sqrt{(b_{1}+b_{2})/(b_{1}-b_{2})}$, and two levels are parallel to each other, as illustrated in Figs. 4(e). This case, however, is different from another known solvable model, called the four-state bow-tie model bow-tie , that has the same geometry of level crossings but different constraints on coupling constants. We find it surprising that the four-state bow-tie model is not a special case of our model. This observation means that, when some of the diabatic levels are parallel to each other, one has more options to satisfy ICs. Determining constraints on parameters in such cases is, however, more difficult because nontrivial constraints emerge then in higher order perturbation theory. We consider such a model in the following section.

IV Distorting Bosonic model

Let us now consider the bosonic model represented by Fig. 1(b) and look for possibility to distort it in order to generate a new solvable model.

IV.1 Model with new couplings only

For simplicity, first, we will not change the slopes of diabatic levels at all and try only to find a different choice of coupling parameters that would satisfy condition (ii). According to the strategy that we described in previous section, we augment the set of parameters to generate the model shown in Fig. 5(a) with six unknown coupling parameters $g_{1}$-$g_{6}$. It is easy to check numerically that a random choice of these couplings leads to permanent gaps between all adiabatic energy levels, so this model is generally not integrable.

To narrow the search for proper couplings, we treat coupling constants as small. Diabatic levels 2 and 6 cross at $t_{cr}=(e_{1}+e_{2})/2$. Let $E_{k}$, $k=1,\ldots,6$ be diabatic energies of states at this time moment. We will need

[TABLE]

At second order in strength of parameters $g_{i}$, coupling between levels 2 and 6 is then given by

[TABLE]

In order not to open the minigap between levels 2 and 6, this coupling should be zero. Since $E_{2}-E_{4}=-(E_{2}-E_{5})$, this gives a constraint on couplings:

[TABLE]

Let us test now the “no-coupling” condition for levels 3 and 4. These diabatic levels cross at $t_{cr}=2e_{1}-e_{2}$, and at this moment

[TABLE]

Generally, a minigap between levels 3 and 4 would open as the third order correction:

[TABLE]

By equating $v_{34}$ to zero we find

[TABLE]

Taking ratio of Eqs. (36) and (38) we find a constraint $g_{3}=\pm\sqrt{2}g_{1}$. Similarly, by considering crossing of levels 1 and 5, we find constraint $g_{6}=\pm\sqrt{2}g_{2}$.

Such constraints are insufficient to determine values of parameters $g_{4}$ and $g_{5}$. This is expected because when a model has parallel diabatic levels, such levels do not give rise to adiabatic crossing between them, which reduces the number of independent constraints on coupling parameters. So, for such cases, one should develop higher order perturbative tests for existence of exact level crossings. Here, we will not go for this because the already found constraints are sufficiently helpful. We just recall that presence of exact crossing points usually requires some symmetric choices of couplings. After trying several candidate sets, we found that the choice shown in Fig. 5(b) does lead to the desired exact eigenvalue crossings, as we prove in Fig. 5(c). The Hamiltonian of this model reads:

[TABLE]

with arbitrary real parameters $e_{1}$, $e_{2}$, $g$, and $\gamma$.

According to Sec. II.3, as in the original bosonic model, the matrix of transition probabilities factorizes as a product of four transition probability matrices describing all crossing points with nonzero couplings. Consider first the crossing point of diabatic levels 3, 5, and 6. After the time shift that places this crossing point at $t=0$, corresponding projected to the phase space of these levels Hamiltonian becomes a 3$\times$3 matrix:

[TABLE]

Transition probabilities for evolution with this Hamiltonian are known because it is a special case of the Hamiltonian of the solvable bow-tie model. Let us define parameters

[TABLE]

The transition probability matrix restricted to the crossing levels is then given by hioe

[TABLE]

or in the basis of all six states

[TABLE]

Similarly, $\hat{P}_{2}$ and $\hat{P}_{3}$ are unit matrices except the 2$\times$2 blocks that describe transitions near crossings of level pairs (2,4) and (2,5):

[TABLE]

and transition probabilities at the crossing of levels 1, 4, and 6 are again obtained from the solution of the three-state bow-tie model:

[TABLE]

According to Sec. II.3, the total transition probability matrix is given by $\hat{P}_{\rm fin}=\hat{P}_{4}\hat{P}_{3}\hat{P}_{2}\hat{P}_{1}$, or explicitly

[TABLE]

In Fig. (6) we test predictions of Eq. (46) and find that agreement with numerical simulations is again excellent. Note that one of the signatures of integrability, which can be tested numerically, is that all transition probabilities are independent of rescaling the distances between parallel diabatic levels, as we show in Fig. 6(a). This property of the semiclassical ansatz is found in all solved models with parallel diabatic levels cQED-LZ ; do ; gbow-tie ; no-go .

IV.2 Deformation leading to path interference

Further extensions of the found solution can be searched among models with different level slopes but the same positions of nonzero couplings and keeping parallel levels to be parallel after the deformation. As physical meaning of such models is obscure, we will restrict ourselves to a particularly symmetric case in which we exchange the slope of level 6 with that of levels 4 and 5, as shown in Fig. 7(a). We adjust distances between parallel diabatic levels to satisfy IC (i) without changing the number and types of the crossing points. Couplings are found using IC (ii) by using the third order perturbation theory, just as in the previous subsection. The resulting model has the Hamiltonian

[TABLE]

which is parametrized by three independent parameters $e$, $g$, and $\gamma$. Figure 7(b) shows numerical proof of the existence of three exact level crossings, as required by IC (ii) in this case. The fact that some of the couplings have different signs and that all three exact crossings happen simultaneously at $t=0$ indicates that this model can describe an odd-spin system four-LZ .

Examining Fig. 7(a) we find that there is a possibility for interference between semiclassical trajectories, so we have to construct semiclassical ansatz from scattering rather than probability matrices of each crossing point. Consider first the crossing point of levels 3, 5, and 6. The effective Hamiltonian near this point is a 3$\times$3 matrix

[TABLE]

This is again the Hamiltonian of the three-state bow-tie model hioe but with different from Eq. (40) order of level slopes. Let

[TABLE]

The scattering matrix for the Hamiltonian (48) can be obtained by taking square roots of known transition probabilities (given explicitly in Ref. hioe ), and adding factors $i$ taken to the power of the number of direct transitions that connect the levels:

[TABLE]

Note that this scattering matrix is truncated so that we drop its dependence on coupling-dependent phase factors which do not change final transition probabilities in all known exactly solvable models. Lifting $\hat{S}_{1}$ to the phase space of all six states, we get

[TABLE]

Chronologically, the next two crossing points correspond to pairwise transitions between levels (2,5) and then (2,4). Corresponding scattering matrices have the form (8) with one subtlety: due to different signs of couplings near those intersections, off-diagonal elements of corresponding scattering matrices have opposite signs:

[TABLE]

The last crossing point is between levels 1, 4 and 6, which is again described by the three-state bow-tie model:

[TABLE]

The final scattering matrix of the model (47) is given by

[TABLE]

Let us introduce parameters

[TABLE]

Multiplying scattering matrices in (54) and taking the absolute values squared of all elements of $\hat{S}^{\prime}$ we obtain the final matrix of transition probabilities:

[TABLE]

Figures 7(c-d) show perfect agreement of this prediction with results of our direct numerical simulations. This model example demonstrates basic nontrivial features of the semiclassical ansatz such as the possibility to describe simultaneous multiple level intersections and interference of semiclassical trajectories that connect complex crossing points.

V Model of Interacting Fermions

Now, we turn to the question whether one can use our strategy in order to generate physically meaningful solvable MLZ models with combinatorial complexity. At least one such a model has been already identified previously cQED-LZ . It is important, however, to know how common such models are and how different they can be. Here we will present quite a different from cQED-LZ model, which describes nonlinear interactions of fermions.

V.1 Six-state sector

The logic that brought us to this model was the following. We started with observation that there are two already known six-state solvable MLZ models with very similar diabatic level diagrams:

[TABLE]

and

[TABLE]

The matrix of transition probabilities in the model (56) was derived in multiparticle as a two-fermion extension of the DO model. The model (57) was found and solved only recently using ICs six-LZ .

Similarity between Hamiltonians (56) and (57) suggests that they can be special instances of a more general model. The simplest candidate for such a generalization is the following 6$\times$6 Hamiltonian matrix:

[TABLE]

where $x$ is a new parameter. At $x=0$ and $x=2$ the model (58) transfers into the models (56) and (57), respectively.

It is straightforward to check that IC (i) is always satisfied at $x\neq 1$, however, when $x$ crosses the value $x=1$ chronological order of diabatic level crossigs changes. We tested validity of IC (ii) for both $x<1$ and $x>1$ numerically and found that three crossing points of diabatic levels with zero direct coupling between them always resulted in exact crossings of adiabatic energy levels, as illustrated in Figs. 8(a-b).

Since ICs are satisfied, we can apply the semiclassical ansatz to derive transition probability matrices. However, long calculations can be avoided because this analysis is completely analogous to the one in Ref. six-LZ that was used to solve the model (57). Namely, in the case with $x>1$, topology of the diabatic level graph of the model (58) and chronological order of level crossings is the same as for the graph of the model (57), while coupling constants are also the same. So the transition probability matrix does not depend on $x$, as far as $x>1$, and solution derived in Ref. six-LZ for the model (57) is equally applicable to any case with $x>1$. Similarly, at any $x<1$, topology of the graph of diabatic levels and chronological order of level crossings is the same as in the model (56). Hence, the transition probability matrix that is given by the semiclassical ansatz at any $x<1$ is the same as the one at $x=0$, which was found analytically in Ref. multiparticle .

Let us introduce parameters:

[TABLE]

Final results for the two sectors of the model (58) read:

[TABLE]

and

[TABLE]

In Figs. 8(c-d) we show comparison of theoretical predictions of the semiclassical ansatz, Eqs. (60)-(61), with results of direct numerical simulations. Transition probabilities are plotted versus parameter $x$. As predicted by the semiclassical ansatz, there is no dependence of such probabilities on $x$ except at the critical point $x=1$, where sharp change of behavior happens. Very small deviations from theoretical predictions can be noticed for some numerical points near $x=1$, however, this is explained by the fact that it takes more time for transition probabilities to saturate in the vicinity of $x=1$ than the time interval considered in our numerical simulations. We also note that parameters in Figs. 8(c-d) were intentionally chosen mainly from the nonperturbative domain $|e_{i}-e_{j}|\sim|g_{ij}|$ for $i<j=1,2,3$ in order to minimize the possibility that agreement with theory is accidental.

V.2 Multiparticle generalization

To search for further extension of the model (58) we recall that the model (56) was obtained in Ref. multiparticle by populating the four-state DO model with two noninteracting fermions. So, the Hamiltonian of the model (56) can be written as a quadradic Hamiltonian of four interacting fermionic modes multiparticle . As the model (58) is different from (56) only by constant terms on the main diagonal, we searched for the secondary quantized Hamiltonian of the model (58) by adding additional quartic terms to the fermionic representation of the model (56). We found the Hamiltonian

[TABLE]

where $\hat{c}_{k}$ and $\hat{d}$ are fermionic annihilation operators and $e_{k}$, $g_{k}$, $x$ are constants. The matrix form of the Hamiltonian (62) is the same as in Eq. (58) if we choose $N=4$ in the two-particle sector and the basis

[TABLE]

The Hamiltonian (62) conserves the number of fermions, $N_{F}$, which we can treat as a free integer parameter of the model. We can also allow $N$ to be an arbitrary positive integer. Generally, we can interpret $N$ as the number of different “quantum dots” in which spinless fermions of the model reside. There is a single energy level in each quantum dot so, due to the Pauli principle, maximum one fermion can reside in each dot. One of such quantum dots, which is described by the operator $\hat{d}$, is special: fermions can produce transitions between different states only via this dot, which is coupled to all other dots by tunneling amplitudes $g_{k}$, $k=1,\ldots,N-1$, as shown in Fig. 9(a) for the case $N=4$. Moreover, we assume that only the energy of the fermion in this dot is time-dependent, e.g., due to a linearly time-dependent gate voltage. Parameter $x$ describes the nonlinear effect of coupling of the special dot to other dots. This coupling is activated every time the special dot is occupied. Physical meaning of parameter $e_{k}$ is merely the energy of a fermion when it resides in the quantum dot with index $k\in 1,\ldots,N-1$.

Does the Hamiltonian (62) satisfy ICs? In this article, we do not provide the definite answer but we do find strong evidence in support of the conjecture that this model is, indeed, solvable. For $N_{F}>1$, there are $n=N!/[N_{F}!(N-N_{F})!]$ different states in which one can place $N_{F}$ fermions among $N$ sites. For example, at $N=4$ and $N_{F}=2$, there are six states available for dynamics, and we already provided the evidence that this model satisfies ICs and its solution coincides with prediction of the semiclassical ansatz.

Let us next consider the case $N_{F}=1$, and look at evolution in the basis of Fock states: $|k\rangle\equiv\hat{c}_{k}^{\dagger}|0\rangle$, $|N\rangle\equiv\hat{d}^{\dagger}|0\rangle$, where $|0\rangle$ is the vacuum state of fermions. The Hamiltonian in this basis has the form

[TABLE]

where the upper index in $\hat{H}^{(N,1)}$ means that the Hamiltonian (62) is restricted to one fermion that can be in $N$ different states. The Hamiltonian (64) is precisely the $N$-state DO model, which is exactly solvable do and for which ICs are trivially satisfied.

V.3 Ten-state sector

The next in complexity sector of the model (62) has two fermions that can occupy any of five quantum dots, as shown in Fig. 9(b). In the basis

[TABLE]

the Hamiltonian (62) is a 10$\times$10 matrix, which we will not show here explicitly. Diabatic levels of this model are shown in Fig. 10. The first IC is trivial but tedious to verify, so we will skip this step here. In Fig. 11(a), we show adiabatic energies of this Hamiltonian that we calculated numerically at $x=2.5$. One can visually see that there are twelve exact crossing points in this figure, as it is needed to satisfy IC (ii).

Here we note that not all exact crossings in Fig. 11(a) are robust in the sense that by increasing couplings some of the crossings can annihilate with each other. However, as it was discussed in cQED-LZ , IC (ii) requires presence of a specific number of exact crossing points only at finite but sufficiently small values of all couplings.

V.3.1 Transition probabilities at $x<1$

If we restrict attention only to values $x<1$, topology of the graph of diabatic states is the same as at $x=0$; hence we can safely set $x=0$ in order to determine the transition probability matrix. The point $x=0$ corresponds to the noninteracting case with a quadratic fermionic Hamiltonian. In the Heisenberg’s picture, evolution of operators has then the same form as evolution of amplitudes in the DO model multiparticle . Consider the DO model with $N$ levels marked so that the level with index $N$ has the positive slope and levels with indexes $i=1,\ldots,N-1$ belong to the band of parallel levels with diabatic energies $e_{1}<e_{2}<\ldots<e_{N-1}$. Solution of the Heisenberg equation for annihilation operators can then be written in the form

[TABLE]

where we identify $\hat{d}\equiv\hat{c}_{N}$. Transition amplitudes in the DO model are known be ; do ; multiparticle :

[TABLE]

where we dropped the coupling-dependent phases of amplitudes except the imaginary unit factors that are acquired at every turning to another level. It was discussed in four-LZ that all other phase factors cancel in semiclassical calculations of transition probabilities in all known exactly solvable models, and we assume this is happening here too. We do not pursue the mathematical proof of this fact because we check the final result numerically.

Consider a sector of such a model with $M$ fermions. Let indexes $\gamma_{1},\ldots,\gamma_{M}$ correspond to quantum dots that are initially populated with fermions at $t\rightarrow-\infty$, and let $\alpha_{1},\ldots,\alpha_{M}$ be indexes of the dots populated with fermions at $t\rightarrow+\infty$ ($\alpha_{1}<\alpha_{2}<\ldots<\alpha_{M}$ etc.). The transition probability between such two states is the same as the product of average populations of final dots:

[TABLE]

where

[TABLE]

and relation between combined indexes, such as $\gamma_{1},\ldots,\gamma_{M}$, and numbering of diabatic levels is explained in Fig. 9 for sectors with six and ten interacting states. Simple algebra leads to the final result:

[TABLE]

where

[TABLE]

Equations (74)-(75) with (67)-(71) provide a simple formal algebraic solution of the problem at $x=0$ and consequently at any $x<1$.

Physically more useful characteristic can be the average population of the site (quantum dot) $\alpha$ irrespectively of occupations of other states. For independent fermions, this characteristic was previously obtained in multiparticle :

[TABLE]

So, we can only add here that Eq. (76) is equally valid for the case with quartic interactions at arbitrary $x<1$.

Let us now restrict to the case $N=5$ and $N_{F}=2$, which is illustrated in Fig. 10(a). The initially occupied state 3 in Fig. 10(a) corresponds to initially occupied sites $\gamma_{1}=3$ and $\gamma_{2}=5$ in Fig. 9(b). Finally occupied level 2 in Fig. 10(a), for example, corresponds to finally occupied quantum dots $\alpha_{1}=2$ and $\alpha_{2}=5$ in Fig. 9(b).

For initial state 3 (see Fig. 9(b)) and arbitrary final state $n=1,\ldots,10$, Eq. (74) then predicts the following vector of transition probabilities:

[TABLE]

Similarly, if initial level is 2 then transition probability vector is

[TABLE]

V.3.2 Transition probabilities for $x>1$ in ten-state sector

We do not have a simple formula for $x>1$. Behavior of transition probabilities in this case is more complex because of the possibility of interference between different semiclassical trajectories. For example, there are two paths that connect initial level 2 and final level 10 in Fig. 10(b), which we mark by red and blue arrows. So, we have to apply the semiclassical ansatz, which is equivalent to deriving all semiclassical paths that connect pairs of initial and final states in the diabatic level diagram and summing amplitudes of all these paths assuming that contributions from all level crossings to each path amplitude factorize.

According to the semiclassical rules, each time a trajectory turns to another level at a crossing point with some pairwise coupling $g_{k}$ it receives the factor $\pm i\sqrt{q_{k}}$, where ($\pm$) is the sign of the coupling. The amplitude of the trajectory marked by red arrows in Fig. 10(b) is then

[TABLE]

Similarly, the trajectory marked by blue arrows has the amplitude

[TABLE]

Comparing them we find that amplitudes of red and blue trajectories have the same absolute magnitude but different signs. Therefore, the total amplitude is $S_{10,2}=S_{\rm red}+S_{\rm blue}=0$. Similarly, we can find other transition probabilities starting with initial level 2:

[TABLE]

If initial state is 3 then there is no path interference and transition probability vector is

[TABLE]

We compare predictions of Eqs. (77)-(80) with nuremical results in Fig. 11(b-c) and again find perfect agreement. Thus, we verified that ICs are valid for, at least, three simplest sectors of the model (62), and that numerically calculated transition probabilities in these sectors coincide with analytical predictions of the semiclassical ansatz within three significant digit precision for all considered parameter values. Complexity of the phase space of the model (62) is growing very quickly with increasing numbers of fermions and quantum dots, which makes further numerical tests problematic but already available results do support our conjecture that the full model (62) is integrable.

VI Discussion

We demonstrated the strategy to search for new solvable MLZ models. The idea is to start with a graph of some already known solvable model and then distort it keeping only topology of the original graph and allowing a broader choice of parameters. ICs can then be converted to a set of algebraic equations on new parameters. Thus, IC (i) of zero area inside any closed path of a graph leads to constraints on slopes of diabatic levels. Treating couplings as small, one can then derive the necessary conditions that such couplings should satisfy in order to agree with IC (ii), which is achieved with the quantum perturbation theory. Sufficiency of these conditions is tested numerically. Following this path, we constructed several new solvable models of different complexity, and we provided detailed numerical check of our predictions for transition probabilities in these models.

Based on our findings, it seems that solvable MLZ models are much more common than it has been previously expected. By creating symmetrized products of already solved systems or populating these systems with some number of noninteracting bosons and fermions, as it is described in multiparticle , we can generate more complex graphs that, in turn, can be distorted, leading to new integrable MLZ models and so on. Moreover, special cases when some of the diabatic levels become parallel produce additional opportunities to distort parameters without breaking ICs. Such opportunities to find new solvable models using distortions can be combined also with previously developed methods to derive systems with smaller phase space sizes from the bigger models. This can be done either by identifying the bigger model as a composite one constraints or by considering limits with level degeneracies and decoupling some of the states in these limits gbow-tie ; cQED-LZ .

While the origin of the MLZ integrability remains puzzling, discoveries of new solvable models can have important consequences. For example, we found that transition probability matrices in newly solved models had the same symmetries as in models that had been rigorously solved by other means. Understanding such relations may lead to the rigorous proof of ICs and validity of the semiclassical ansatz. We also find it surprising that, in all promising candidate models, perturbatively derived constraints on coupling constants have lead to truly exact eigenvalue crossings. Formally, perturbatively derived constraints are only necessary conditions that do not have to be sufficient. It seems that, somehow, imposing IC (i) facilitates the emergence of exact crossing points so that a simple perturbative test, which is performed up to some finite order of perturbation theory, provides not only necessary but also sufficient conditions for existence of such points. This observation is one of already many surprisingly fortunate but puzzling properties of MLZ models.

Exact and approximate methods to connect asymptotic behavior of ordinary differential equations are of general interest to mathematics and mathematical physics joye-LZ ; aoki-LZ . For example, the most physically important special functions are solutions of differential equations with well understood asymptotics. Future studies should answer the question whether solvable MLZ models represent a new class of physically important solvable differential equations or such models can be understood in terms of the already well studied special functions, such as the generalized hypergeometric function.

Finally, some of the solutions that have been obtained with ICs are important for practical reasons. For example, our fermionic model in Sec. V describes the process of a quantum level interaction with a fermionic bath when this level is driven by a time-dependent gate voltage. Such models have been of interest previously due to the possibility for control of quantum states with minimal dissipation multiparticle ; levitov-lz1 . So, with MLZ solutions, it becomes possible to explore complex explicitly driven quantum mechanical systems analytically and exactly. We hope that many other physically interesting processes will be eventually understood with such solvable models.

Appendix A Counterexample to previous version of the first IC

Our definition of the first integrability condition in Sec. II.2 is different from the original definition in six-LZ . Previously, condition (i) requested that if, on the diabatic level diagram, there are more than one semiclassical paths that connect one initial level at $t=-\infty$ with another level at $t=+\infty$ then dynamic phases, i.e. integrals of diabatic energies over time, of such trajectories should be the same. The latter property is actually the consequence of IC (i) defined in the present article because the dynamic phase is the area that a trajectory in the level diagram sweeps over the time axis. Equality of such phases means that the path that goes forward along one of these trajectories and then returns to the initial point along the other trajectory encloses zero area in this diagram.

The new definition, however, is more restrictive because it includes the case when different trajectories do not interfere in the physical sense but the graph still has closed loops. We show an example of such a model in Fig. 12(a). One can verify numerically that IC (ii) is satisfied when

[TABLE]

i.e., at such relations among parameters there is an exact crossing point of adiabatic levels near the intersection of levels 1 and 2. As for the first IC (i), it is not satisfied at $e\neq 0$ because the graph has a closed path that encloses a finite area. So, semiclassical ansatz should fail.

On the other hand, each semiclassically allowed trajectory in Fig. 12(a) is unique in the sense that there is no other path that respects causality and connects the same pair of levels, i.e., there are no physically interfering trajectories. So, if definition in Ref. six-LZ is generally true then the semiclassical ansatz should apply.

Figure 12(b) compares prediction of the semiclassical ansatz for this model with results of numerical simulations. Although the semiclassical ansatz gives a reasonable fit to the data, there is noticable disagreement. This makes us to conclude that the more restrictive definition of condition (i) in Sec. II.2 of the present article should be preferred. So far, we are not aware of any MLZ model that would satisfy new version of ICs (i)-(ii) but would not be solvable by the semiclassical ansatz.

Appendix B Analogy between first integrability condition and WKB description of classically integrable systems

Here we would like to point that the geometric interpretation of the first integrability condition, presented in Sec. II.2, may provide a link between the integrable MLZ problems and semiclassical WKB description of classically integrable systems. Consider the graph $G=(G_{0},G_{1})$ associated with the diabatic level diagram of a given $n\times n$ MLZ problem. The vertices in $G_{0}$ of this graph are given by those intersection points for which the coupling between the intersecting diabatic levels is non-zero, and the edges in $G_{1}$ represent the segments of the straight lines between two consecutive scattering events. Let $X$ be the geometric realization of the graph $G$, i.e., $X$ stands for the graph that is viewed as a set of segments $s_{e}$, associated with its edges $e\in G_{1}$, connected to each other at the vertices, rather than just a collection $G=(G_{0},G_{1})$ of edges and vertices. By definition, every point of $X$ defines diabatic energy $\omega(t)$ of a diabatic state represented by this point. Hence, the geometric space $X$ of our MLZ problem has a natural immersion $f:X\looparrowright\mathbb{R}^{2}$, which is what is shown in the diabatic level diagram.

In the $(t,\omega)$ plane we have a rank $1$ form $\alpha=\omega dt$, which, being restricted (strictly speaking pulled-back) to $X$ gives rise to a form $f^{*}\alpha$ on $X$. The first integrability condition (area rule) is then equivalent to the requirement that for any cycle $C$ in $X$ we have

[TABLE]

which, due to the quasi-one-dimensional nature of $X$ is equivalent to existence of a global function ${\cal S}:X\rightarrow\mathbb{R}$ with $d{\cal S}=f^{*}\alpha$, or, stated differently, ${\cal S}=\int f^{*}\alpha$. We will refer to ${\cal S}$ as the action function.

For $1$-dimensional quantum systems, the WKB approximation is based on the action function that can be represented as ${\cal S}=\int\alpha$, with $\alpha=pdx$. The WKB form $\psi(x)\sim\exp\left(i\hbar^{-1}{\cal S}(x)\right)$ of an eigenstate may not be generalizable to the multi-dimensional case, however if an $n$-dimensional system is integrable, i.e, it has $n$ mutually commuting integrals $I_{j}$ of motion, a multidimensional counterpart of a WKB wave function can be obtained. Indeed, switching to the classical limit and further fixing the values of all integrals of motion $I_{j}(p,x)=E_{j}$, where $E_{j}$ are some constants for $j=1,\ldots,n$, we can solve these equations with respect to $p$, resulting in the set of functions $p_{j}(x)$, which give rise to a $1$-form $\alpha=\sum_{j}p_{j}dx^{j}$. It is possible to show that the integrability condition $\{I_{j},I_{k}\}=0$, where $\{\ldots\}$ stands for Poisson brackets for all $j$ and $k$, is equivalent to the closeness condition $d\alpha=0$, the latter meaning that one can define a global action function ${\cal S}=\int\alpha$ that determines the zero-order approximation for a semiclassical WKB function that has exactly the same form as in $1$-dimensional case.

Summarizing the analogy: in both cases, namely MLZ problems and quantized integrable systems, the integrability condition, at least in the lowest (principle) order of the semiclassical expansion, boils down to the existence of a globally defined action function ${\cal S}$. We would like to admit explicitly that at this point we do not see how exactly the above analogy can be applied. It may mean, for example, that integrable MLZ systems describe some ballistic, i.e. large kinetic energy, limits of integrable many-body $1$-dimensional quantum systems. Finally, we note that this analogy resembles the property, found at least in some MLZ models, that solvable MLZ Hamiltonians have nontrivial commuting operators that depend polynomially on time yuzbashyan-LZ .

Acknowledgment. This work was carried out under the auspices of the National Nuclear Security Administration of the U.S. Department of Energy at Los Alamos National Laboratory under Contract No. DE-AC52-06NA25396. V.Y.C. was supported by the National Science Foundation under Grant No. CHE-1111350. N.A.S. also thanks the support from the LDRD program at LANL.