Multistate Landau-Zener models with all levels crossing at one point

TL;DR

This paper explores the properties and integrability of multistate Landau-Zener models where all levels cross at a single point, providing analytical solutions and applications to quantum systems.

Contribution

It introduces dual models for previously solved MLZ models, analyzes their integrability, and applies the results to Bose condensate conversion and specific quantum models.

Findings

Analytical scattering matrices for dual MLZ models

Purely algebraic solutions for bowtie and Tavis-Cummings models

Application to Bose condensate conversion during Feshbach resonance

Abstract

We discuss common properties and reasons for integrability in the class of multistate Landau-Zener (MLZ) models with all diabatic levels crossing at one point. Exploring the Stokes phenomenon, we show that each previously solved model has a dual one, whose scattering matrix can be also obtained analytically. For applications, we demonstrate how our results can be used to study conversion of molecular into atomic Bose condensates during passage through the Feshbach resonance, and provide purely algebraic solutions of the bowtie and special cases of the driven Tavis-Cummings model (DTCM).

Figures

Equations

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

$$\displaystyle\hat{H}(t)=\left(\begin{array}[]{cccc}\beta_{1}t&xg_{1}&g_{2}&0\\ xg_{1}&\beta t&0&-yg_{2}\\ g_{2}&0&-\beta_{2}t&g_{1}\\ 0&-yg_{2}&g_{1}&-\beta_{1}t\end{array}\right),$$

$$x=\sqrt{\frac{{\beta }_1-\beta }{{\beta }_1-{\beta }_2}},y=\sqrt{\frac{{\beta }_1+\beta }{{\beta }_1+{\beta }_2}}.$$

$$i{\dot{a}}_0={\beta }_{0}t{a}_0+\sum _{k=1}^{N-1}{g}_k{a}_k;i{\dot{a}}_k={\beta }_{k}t+g_k{a}_0,k\ne 0.$$

$$i{\dot{a}}_n={\beta }_{n}t{a}_n+g_n{a}_{n+1}+g_{n-1}{a}_{n-1},n=1,\dots ,N,$$

$$g_n=g\sqrt{N_B+n}\sqrt{S(S+1)-(S-n+1)(S-n)}.$$

$$\hat{H}=\left(\begin{array}[]{cccccc}\beta_{1}t&0&0&g_{13}&g_{12}&0\\ 0&\beta_{2}t&0&g_{23}&0&-g_{12}\\ 0&0&\beta_{3}t&0&-g_{23}&-g_{13}\\ g_{13}&g_{23}&0&-\beta_{3}t&0&0\\ g_{12}&0&-g_{23}&0&-\beta_{2}t&0\\ 0&-g_{12}&-g_{13}&0&0&-\beta_{1}t\end{array}\right).$$

$${\Theta }_{kk}=(-1{)}^{f_k},{\Theta }_{ij}=0,i\ne j.$$

$$\hat{\Theta }\hat{H}(-t)\hat{\Theta }=-\hat{H}(t),$$

$$\hat{\Theta }{e}^{-i\hat{H}(-t)dt}\hat{\Theta }=e^{i\hat{H}(t)dt}.$$

$$\hat{U}(T,-T)=\underset{dt\to 0}{\lim }\prod _{n=-T/dt}^{T/dt}{e}^{-i\hat{H}(t_n)dt},$$

$$\hat{\Theta }\hat{U}(T,-T)\hat{\Theta }=\prod _{n=-T/dt}^{T/dt}\hat{\Theta }{e}^{-i\hat{H}(t_n)dt}\hat{\Theta }=$$

$$=\prod _{n=-T/dt}^{T/dt}{e}^{i\hat{H}(-t_n)dt}={\hat{U}}^{†}(T,-T).$$

$$S_{nm}=(-1{)}^{f_n+f_m}{S}_{mn}^{\ast },P_{nm}=P_{mn}.$$

$$\mathrm{Tr}\left[\hat{U}(T,-T)\hat{\Theta }\right]=\mathrm{Tr}\left[\hat{U}(0,-T)\hat{\Theta }\hat{U}(T,0)\right].$$

$$e^{-i\hat{H}(T)dt}\hat{\Theta }{e}^{i\hat{H}(-T)dt}=e^{-i\hat{H}(T)dt}\hat{\Theta }{e}^{-i\hat{H}(-T)dt}\hat{\Theta }\hat{\Theta }=$$

$$=e^{-i\hat{H}(T)dt}{e}^{i\hat{H}(T)dt}\hat{\Theta }=\hat{\Theta }.$$

$$\mathrm{Tr}\left[\hat{S}\hat{\Theta }\right]=\mathrm{Tr}\left[\hat{\Theta }\right],$$

$$\sum _{n\in G_2}{S}_{nn}-\sum _{m\in G_1}{S}_{mm}=N-2M,$$

$$\sum _{n\in G_2}\sqrt{P_{nn}}-\sum _{m\in G_1}\sqrt{P_{mm}}=N-2M.$$

$$\text{noindent}{S}_{11}$$

$$\displaystyle\noindent{\rm Det}\left(\begin{array}[]{cc}S_{11}&S_{12}\\ S_{21}&S_{22}\end{array}\right)$$

$$\displaystyle{\rm Det}\left(\begin{array}[]{cccc}S_{11}&S_{12}&\cdots&S_{1M}\\ S_{21}&S_{22}&\cdots&S_{2M}\\ \vdots&\cdots&\ddots&\vdots\\ S_{M1}&\cdots&\cdots&S_{MM}\end{array}\right)=e^{-\pi\sum\limits_{k=M+1}^{N}\sum\limits_{r=1}^{M}\frac{|A_{kr}|^{2}}{(\beta_{r}-\beta_{k})}},\quad M=1,\ldots,N-1.$$

$$S_{11}{S}_{22}-(-1{)}^{f_1+f_2}|S_{12}{|}^2=e^{-\pi \left(\sum _{k\in \bar{G}(1)}|A_{1k}{|}^2/({\beta }_1-{\beta }_k)+\sum _{k\in \bar{G}(2)}|A_{2k}{|}^2/({\beta }_2-{\beta }_k)\right)},$$

$$e^{-2\pi \sum _{k\in \bar{G}(1)}|A_{1k}{|}^2/({\beta }_1-{\beta }_k)}{P}_{22}={\left((-1{)}^{f_1+f_2}{P}_{12}+e^{-\pi \left(\sum _{k\in \bar{G}(1)}|A_{1k}{|}^2/({\beta }_1-{\beta }_k)+\sum _{k\in \bar{G}(2)}|A_{2k}{|}^2/({\beta }_2-{\beta }_k)\right)}\right)}^2.$$

$$S_{mn}={\alpha }_{mn}{i}^{f_n+f_m}{e}^{i({\varphi }_m-{\varphi }_n)},$$

$$c_3\equiv S_{13}{S}_{32}{S}_{21},c_4\equiv S_{12}{S}_{24}{S}_{43}{S}_{31}.$$

$$\mathrm{Im}[c_3]=\mathrm{Im}[c_4]=0.$$

$$\displaystyle\hat{H}=\left(\begin{array}[]{cccc}\beta_{1}t&g_{1}&0&0\\ g_{1}&\beta_{2}t&g_{2}&0\\ 0&g_{2}&\beta_{3}t&g_{3}\\ 0&0&g_{3}&\beta_{4}t\end{array}\right).$$

$${\alpha }_{mn}={\alpha }_{nm},$$

$${\hat{\alpha }}^2=\hat{1}.$$

$${\alpha }_{nn}=(-1{)}^{f_n},\mathrm{for}\hat{A}=0.$$

Full text

I Introduction

This article continues the series of publications multiparticle ; six-LZ ; cQED-LZ ; chen-LZ ; constraints about exact results in the MLZ theory be . This theory deals with explicitly time-dependent Schrödinger equations of the form

[TABLE]

where $\hat{A}$ and $\hat{B}$ are Hermitian time-independent $N\times N$ matrices. An MLZ model is called solvable if one can determine the scattering matrix $\hat{S}$ or transition probability matrix $\hat{P}$, $P_{nm}\equiv|S_{nm}|^{2}$, for evolution from infinite negative to infinite positive time values. Recent discovery of integrability conditions six-LZ has led to numerous fully solvable models of this type. Complexity of these models varies from the elementary two-state case LZ ; book-LZ and few-state models six-LZ to systems of spins and fermions with nonlinear interactions and combinatorially large phase space cQED-LZ ; chen-LZ .

There are two reasons making the MLZ theory attractive. First, scattering matrices in solvable MLZ models are obtained explicitly in terms of the known special functions. This is in contrast to time-independent systems that are considered solvable by the algebraic Bethe ansatz. For finite size models, Bethe ansatz generally produces only implicit solutions that depend on parameters satisfying complex nonlinear algebraic equations tv-bethe ; qd-goden , which in turn are often almost as hard to study numerically as to diagonalize the full Hamiltonian by generic algorithms. For this reason, there are still debates about what are the criteria for calling a finite size quantum model integrable yuzbashyan-LZ . Within the MLZ theory, the notion of quantum integrability acquires conventional meaning of completely analytical description of the scattering process.

Second, MLZ theory opens unusual opportunities for studies of mesoscopically complex quantum systems without approximations. Many solved MLZ models describe situations that have previously been discussed in experimental context, for example, they describe interacting spin systems in linearly time-dependent magnetic fields li-16LZ . The same models have been used previously to explain experiments with molecular nanomagnets app-spin , although analytical solutions of these models had not been known. Similarly, the DTCM cQED-LZ , which belongs to the class (1), had been extensively studied by variety of approximate and numerical methods before its solution was found gurarie-LZ . Originally, this model attracted attention due to applications to experiments with molecular Bose condensates app-bose . More generally, multistate Landau-Zener transitions appear in numerous modern fields of research, including quantum coherence coher , dynamic phase transitions dziarmaga , Landau-Zener interferometry LZ-interferometry , metrology metrology-LZ , quantum control qcontrol , superconductors sc-LZ , ultracold atoms ultracold-LZ , and quantum dots dot-lz-exp1 .

In MLZ theory, eigenvalues of matrix $\hat{B}t$ are called diabatic levels. It is usually convenient to visualize parameters of the Hamiltonian by plotting diabatic levels in a time-energy diagram, as shown in Fig. 1(a). Models in which all diabatic levels are crossing at one point, as shown in Fig. 1(b), are obtained by setting diagonal elements of $\hat{A}$ to zero. Such models are particularly important in MLZ theory because they are used as elementary blocks to construct more complex solvable systems using integrability conditions Quest-LZ . Unfortunately, these conditions do not apply to cases with all crossings happening at a single point. Therefore, such models need specific research on their own.

In this article, we present several findings that reveal symmetries and some of the reasons for integrability in MLZ theory. We also develop methods and mathematical concepts that should be handy in future studies of this topic. The structure of our article is as follows. In section II we collect known solvable MLZ models featuring all levels crossing in one point, and then discuss their bipartite property. In section III we point to another symmetry of solvable systems, and then use these symmetries to design purely algebraic derivation of transition probabilities in the bowtie and five-state DTCM. Section IV discusses Stokes phenomenon in MLZ systems. We show that bipartite property leads to a new class of solvable models with applications to physics of molecular Bose condensates. We summarize our findings in section V.

II Bipartite Models

When all levels are crossing at a single point, it is convenient to illustrate parameters of the Hamiltonian as a connectivity graph, which is different from the diabatic level diagram in Fig. 1. Let us think about the diabatic states as the nodes of the connectivity graph. If direct coupling between two such states, i.e., the corresponding off-diagonal element of the matrix $\hat{A}$, is nonzero, we will connect corresponding nodes with a link. One can then mark links by corresponding couplings and nodes by corresponding level slopes, as shown in Fig. 2(a). This figure contains all information about the Hamiltonian:

[TABLE]

where $g_{1},\,g_{2},\,x,\,y,\,\beta,\,\beta_{1},\,\beta_{2}$ are constant parameters. It is known that this model is solvable Quest-LZ under condition

[TABLE]

Interestingly, the procedure that leads to solution of this model remains mathematically unjustified, although numerous careful numerical checks have always confirmed its validity. Only at special choice of parameters, shown in Fig. 2(b), has this solution been proved both analytically and exactly constraints .

Let us now show, in Fig. 3, connectivity graphs of other known solvable models. Figure 3(a) corresponds to the bowtie model bow-tie , in which an arbitrary number of crossing levels are coupled directly to only one of them. Corresponding Schrödinger equation for amplitudes $a_{n}$, $n=0,1,\ldots,N-1$ reads

[TABLE]

Figure 3(b) corresponds to the LZ-chain model constraints , in which each state is coupled only to states with nearby indexes. Diabatic state amplitudes of this model evolve according to equations

[TABLE]

where we assume that $g_{0}=g_{N}=0$.

Not all models of the form (9) have scattering matrices that are known. The most general solvable LZ-chain model is the DTCM. This model has equidistant level slopes, $\beta_{n}=\beta n$, and couplings given by,

[TABLE]

where $g$ and $N_{B}>-1$ are constant parameters, and $2S+1=N$. In the limit $N_{B}\rightarrow\infty$ but $g^{2}N_{B}={\rm const}$, this model describes a quantum spin of size $\hat{S}$ in a linearly time-dependent field.

General solution of the DTCM is currently known only in the form of an algorithm that allows derivation of transition probabilities in any invariant sector of this model recursively from solutions for lower dimensional sectors cQED-LZ . Complexity of this procedure is growing very quickly. Therefore, only special cases are well understood, such as the four-state sector ($S=3/2$), and the case of arbitrary $S$ but specific initial conditions at $t=-\infty$: with either $|a_{n}|^{2}=\delta_{n1}$ or $|a_{n}|^{2}=\delta_{nN}$.

Any solvable MLZ model can be used to construct more complex ones by “populating” the original model with noninteracting fermions or bosons multiparticle ; dziarmaga ; jipdi . An example of such a composite model is shown in Fig. 3(c). Its Hamiltonian reads

[TABLE]

We provide details of its solution in appendix A.

Finally, Fig. 3(d) shows the connectivity graph of the previously unknown five-state model whose solution we present in appendix B. This model was found and solved by considering a special case of a more general five-state system that can be solved using integrability conditions, as it was done for the four-state model (6) in Quest-LZ . There is no doubt that this class of models can be extended to higher dimensions, although their solutions quickly become very complex.

Let us now observe the most obvious common property of graphs in Figs. 2 and 3: all of the known solvable models are represented by bipartite graphs, i.e. graphs whose nodes can be partitioned into two groups so that nodes of the same group are connected only to nodes of the other group. It is natural to expect then that new solvable MLZ systems can be found among similar bipartite systems.

In what follows, we will refer to all MLZ models with a bipartite connectivity graph and with all diabatic level crossings at a single point in the diabatic level diagram as simply bipartite models. So, let us explore this property and its effect on the scattering matrix.

As there are two interacting groups of states, $G_{1}$ and $G_{2}$, we give the index 1 to the smaller group of nodes and the index 2 to the bigger group in the connectivity graph (or arbitrarily if these groups have equal numbers of elements), and then introduce operator $\hat{\Theta}$, which is represented by a diagonal $N\times N$ matrix with

[TABLE]

Here and in what follows, $f_{k}$ is the index (either 1 or 2) of the group to which the $k$-th diabatic state belongs. It is straightforward to verify that Hamiltonians of bipartite MLZ models satisfy the following identity:

[TABLE]

from which it follows also that

[TABLE]

Consider the formal expression for the evolution matrix

[TABLE]

where $t_{n}=ndt$ and $n$ is increasing in the product from the right to the left. We can then insert resolution of identity $\hat{1}=\hat{\Theta}\hat{\Theta}$ after each factor of the product in (15) to obtain that

[TABLE]

Using definition of the scattering matrix constraints , $\hat{S}=\lim\limits_{T\rightarrow\infty}U(T,-T)$, and writing (16) in components, we find relation between scattering matrix elements and between corresponding transition probabilities:

[TABLE]

Thus, the transition probability matrix of any bipartite model is symmetric. Moreover, diagonal elements $S_{nn}$ of the scattering matrix are purely real.

Symmetries (17) have been known in bowtie and LZ-chain models bow-tie ; constraints . However, algebraic derivation along (13)-(16) can reveal a new and less obvious than (17) constraint: let us take the trace of the evolution operator times $\hat{\Theta}$. Using the fact that under the trace we can perform cyclic permutations, we can rewrite

[TABLE]

Identity (14) then leads to simplifications:

[TABLE]

We can repeat this step until we get rid of all exponential factors. The result is a new identity for the scattering matrix elements:

[TABLE]

or in components:

[TABLE]

where $G_{1}$ and $G_{2}$ are, respectively, the first and the second groups of diabatic states; $N$ is the size of the Hilbert space and $M$ is the number of diabatic states in the first group.

Although constraint (20) is linear in scattering matrix elements it leads to a nonlinear constraint on transition probabilities. For example, for sufficiently small couplings we can rewrite (20) as

[TABLE]

At larger couplings, we already should care about signs of scattering matrix elements, so that generally we can only say that either $S_{nn}=\sqrt{P_{nn}}$ or $S_{nn}=-\sqrt{P_{nn}}$. The fact that simple symmetries, such as Eq. (13), can lead to nonlinear constraints on transition probabilities, such as Eq. (21), is the first important finding of this article.

We conclude this section by pointing to one more consequence of symmetry (17). It is known that elements of the scattering matrix of any MLZ model satisfy the so-called hierarchy constraints (HCs) constraints . Let us assume that indexes of diabatic states are chosen so that $\beta_{1}>\beta_{2}>\ldots>\beta_{N}$. Then the two lowest order HCs read:

[TABLE]

and generally:

[TABLE]

In an arbitrary MLZ model, the first constraint (22) can be converted into expression for the probability to remain on the same level if the system is initially on the level with the extremal slope. The second constraint (25) already depends on phases of amplitudes $S_{22}$, $S_{12}$, and $S_{21}$, so Eq. (25) cannot generally lead to a constraint that is written only in terms of transition probabilities. However, for bipartite models we know that $S_{22}$ is real, and from (17) we find

[TABLE]

where $\bar{G}(k)$ means the set of elements of the group of states connected to the node $k$ directly. Using that $|S_{12}|^{2}\equiv P_{12}$ and (22), we find relation between probabilities $P_{22}$ and $P_{12}$:

[TABLE]

Finally, we note that HCs can be reformulated so that they start not with the highest but rather with the lowest slope level. Such HCs lead to equations analogous to (22)-(32), in which indexes $1,2,\ldots,N$ are replaced by, respectively, $N,N-1,\ldots,1$. This is followed then by change of slope signs, e.g., $\beta_{1}-\beta_{2}\rightarrow\beta_{N-1}-\beta_{N}$, e.t.c.. So, in any bipartite model, not only survival probabilities for two extremal levels are known but also there are two equations of the form (32) that impose coupling-dependent relation between $P_{12}$ and $P_{22}$, as well as between $P_{N,N-1}$ and $P_{N-1,N-1}$.

III Pure gauge phase ansatz

The class of bipartite models is considerably larger than the class of solved models shown in Figs. 2 and 3. What else is special to solvable cases?

It has been noticed previously six-LZ ; Quest-LZ ; gbow-tie that many MLZ models can be solved by a semiclassical ansatz in which direct transition to another level brings a simple phase $i$ to the transition amplitude. This suggests the following ansatz to satisfy constraints (17):

[TABLE]

where $\phi_{m}$ are phases that may depend on parameters of the model, possibly nontrivially, and where $\alpha_{mn}$ are real parameters. We will think of $\alpha_{mn}$ as elements of a new matrix $\hat{\alpha}$. This ansatz reduces the number of unknown phases from $N(N-1)/2$ to just $N$. Symmetry (33) can be verified analytically in several cases: the bowtie model, the models in Figs. 2(b) and 3(c), and the model of arbitrary spin in a linearly time-dependent magnetic field, because their scattering matrices are known bow-tie ; jipdi ; four-LZ .

However, in the DTCM and the general case of the model (6), only transition probabilities have been found so far, while phases of scattering matrix elements remain unknown. So, we cannot verify (33) in them analytically. Instead, we tested for applicability of the ansatz (33) numerically. Since phases $\phi_{m}$ are unknown, we did it indirectly, namely, Eq. (33) predicts that any cyclic product of amplitudes should not depend on all phases $\phi_{m}$. So, if we define

[TABLE]

and set indexes so that $1,3\in G_{1}$ and $2,4\in G_{2}$, then

[TABLE]

We found that, indeed, scattering matrix elements in all previously solved bipartite models satisfy (35).

There is difference from Eq. (17), however, in that Eq. (33) generally does not work if evolution operator $\hat{U}(T,-T)$ at finite $T$ is substituted instead of the scattering matrix. So, if we define $c_{3}(T)$ and $c_{4}(T)$ by replacing elements of $\hat{S}$ by elements of $\hat{U}(T,-T)$ and plot the ratios $r_{3,4}={\rm Im}\left(c_{3,4}(T)\right)/{\rm Re}\left(c_{3,4}(T)\right)$ as functions of $T$, we find that these ratios are initially noticeably nonzero but decrease to zero at increasing $T$, as we show in Fig. 4. Hence, the ansatz (35) works only asymptotically. We found this behavior also in the bowtie and driven Tavis-Cummings models. There are also exceptions, where (33) holds for evolution operator $c_{4}(T)$ at arbitrary $T$, e.g., in the model in Fig. 2(b). We attribute this to presence of additional symmetry in this model, which was discussed in four-LZ .

It is expected that known solvable models do not exhaust the solvable class of bipartite systems. If so, conditions (33) are likely common or even determining this class. Their understanding may shed new light on the phenomenon of integrability in MLZ theory and provide the way to generate new solvable cases. Therefore, below we explore some of the consequences of this symmetry.

III.1 Numerical search for models with symmetry (33)

If MLZ-integrability is indeed related to the symmetry (33), then conditions (35) provide a strategy to search for new solvable models. Consider for example a four-state LZ-chain:

[TABLE]

The most general known solvable model of the type (40) is the four-state DTCM that has equally distant slopes $\beta_{n}=\beta n$ and couplings constants specified by Eq. (10).

Let us consider some case of Eq (40) with a different, from DTCM, choice of slopes and also with fixed values of two out of the three couplings. We can then simulate evolution with the Hamiltonian (40) and check the validity of (35) at different values of the variable coupling. We show results of such a test in Fig. 5. Although constant parameters are chosen differently from the DTCM, this figure shows that there is a coupling value $g_{3}\approx 0.47$, at which imaginary parts of both $c_{3}$ and $c_{4}$ change their signs simultaneously. So, at this point, symmetry (33) is valid and the model is likely solvable although its solution is currently not known.

Several similar tests correctly reproduced coupling values in a few already solved models (not shown). We also found that such points are rather common in bipartite models with four interacting states. However, in models with higher dimensions, scanning of a single parameter is usually insufficient to observe such points. We also found that some symmetric choices of parameters, for example LZ-chains with equal couplings and equidistant level slopes, do not generally satisfy (35) for $N\geq 4$. So, our numerical investigation suggests that there are numerous solvable but still unknown bipartite models with asymmetric choices of parameters.

III.2 Parametrization by orthogonal symmetric matrix

Let us now note that $N$ phases $\phi_{n}$ in (33) can be removed by mere change of phases in definition of $N$ diabatic states followed by redefinition of phases of coupling constants. These phases drop from all expressions for HCs because determinants depend only on cyclic products such as in (34). Therefore, each of $N-1$ independent HCs leads to an equation on elements of $\hat{\alpha}$, and consequently on transition probabilities. This is in contrast to the general bipartite case in which only four such constraints are allowed.

Second, from (17) it follows that $\hat{\alpha}$ has to be symmetric:

[TABLE]

and from unitarity, $\hat{S}\hat{S}^{{\dagger}}=\hat{1}$, we find that $\hat{\alpha}$ is orthogonal. Together, Eqs. (41) and (42) mean that

[TABLE]

An immediate consequence is that eigenvalues of $\hat{\alpha}$ are all $\pm 1$. Moreover, we can find multiplicities of such eigenvalues. At zero couplings we have $S_{nm}=\delta_{nm}$ and eigenvalues of $\hat{\alpha}$ coinciding with its diagonal elements:

[TABLE]

Switching on couplings continuously cannot change discrete eigenvalues. So, we find that the number of $-1$ eigenvalues of $\hat{\alpha}$ equals the number $M$ of nodes in the smaller group, $G_{1}$.

Symmetric orthogonal matrix with $M\leq N/2$ negative eigenvalues can be parametrized with components of $M$ orthonormal vectors:

[TABLE]

Indeed, ${\bf v_{k}}$ is the eigenvector of $\hat{\alpha}$ with eigenvalue $-1$. All real vectors orthogonal to ${\bf v}_{1}\ldots{\bf v}_{M}$ belong to the space with eigenvalues $+1$.

Vector ${\bf v_{1}}$ depends on $N-1$ independent parameters (the last one is fixed by requirement of its unit norm). Vector ${\bf v_{2}}$ depends on $N-2$ unknown parameters because, in addition to normalization, it must be orthogonal to vector ${\bf v}_{1}$. We can continue this counting until we reach ${\bf v_{M}}$ that depends on $N-M$ independent parameters. This corresponds totally to $M(2N-M-1)/2$ unknown parameters of the matrix $\hat{\alpha}$. On the other hand, we have only $N-1$ HCs that express combinations of these parameters via the elements of the model’s Hamiltonian.

It is interesting to know whether knowledge of the symmetry (33) supplemented by HCs, and possibly the constraint (20), is sufficient to derive all transition probabilities in some models. For example, the bowtie model has $M=1$. The number of independent parameters in $\hat{\alpha}$ is then $N-1$, which is the same as the number of independent HCs. The next in complexity case with $M=2$ shows that HCs alone are not sufficient to reconstruct all parameters of $\hat{\alpha}$. For example, for $N=4$ we have then five parameters for $\hat{\alpha}$ and only three independent HCs. In the case of the driven Tavis-Cummings model, however, we have additional information because all transition probabilities from the extremal slope levels are known when $N$ is arbitrary cQED-LZ . In any case, HCs are highly nonlinear, so it is not obvious that they can be disentangled to reconstruct elements of the scattering matrix. In the following two subsections we will show that in some cases this can be done.

III.3 Algebraic solution of the bowtie model

Previously, solution of the bowtie model has been obtained by methods of complex analysis. It is expected, however, that there is a hidden algebraic reason behind any exact solution. Below we claim that, in the bowtie model, the reason is in combination of HCs (30) and the property (33).

For a general case of the bowtie model, it is convenient to choose state indexes so that $\beta_{i}>\beta_{j}$ if $i>j$, and the special level has the slope $\beta_{0}$. So, we allow negative integer indexes. Let there be $m$ diabatic level slopes lower than $\beta_{0}$ and $n$ higher slopes ($m+n+1=N$). Nonzero elements of the Hamiltonian are:

[TABLE]

According to (44), matrix $\hat{\alpha}$ can be expressed via only one real vector

[TABLE]

with a constraint $\sum_{j}v_{j}^{2}=1$. Moreover, it is easy to verify that transition probabilities in this case must depend only on combinations $2v_{k}^{2}$, which we are going to determine now.

Let us consider the first $m$ equations of HCs (30) for the Hamiltonian (45). In the left hand sides, they will have the determinant

[TABLE]

where, for $(j+m+1)$-st hierarchy equation, vector $\mathbf{w}$ is obtained by keeping only components of $\mathbf{v}$ up to $j$-th index (note that $j<0$). A determinant written as (46) can be simplified using Sylvester’s identity sylvester :

[TABLE]

where $\hat{R}$ is $m\times n$ and $\hat{Q}$ is $n\times m$ matrices. Applying this formula to (46) we find

[TABLE]

Thus, first $m$ HCs can be written as

[TABLE]

where we defined

[TABLE]

Solving Eqs. (48) recursively starting from $j=-m$ we find

[TABLE]

Similarly, we can consider the HCs for $n\times n$ part in the right-lower corner of the scattering matrix. They lead to $n$ constraints:

[TABLE]

The component $v_{0}$ of the vector $\mathbf{v}$ is obtained by requiring the unit norm of $\mathbf{v}$:

[TABLE]

Substituting (48)-(52) into (44) we find the scattering matrix up to $N$ phases $\phi_{k}$, which can be gauged away. Taking the absolute value squared of components of this matrix we obtain transition probabilities, which indeed coincide with the ones known for the bowtie model. For example,

[TABLE]

Such a solution of the bowtie model is considerably simpler than the original one in bow-tie .

III.4 Driven Tavis-Cummings model

DTCM is a LZ-chain model with equidistant slopes $\beta_{n}=\beta n$ and couplings given in (10). Alternatively, its Hamiltonian describes interaction of a bosonic mode with an arbitrary size spin:

[TABLE]

where we used the fact that we can always set $\beta=1$ by rescaling time, as discussed in cQED-LZ . Implicitly, the model (54) depends on the additional parameter $N_{B}$, which is the number of bosons in the state with the fully up-polarized spin. It is known how to derive transition probabilities in this model by taking some limit of the known solution of a much more complex model cQED-LZ . The latter, however, has combinatorially large phase space, so this process quickly becomes too involved. The exception is for conditions when the spin is initially fully polarized either down or up, i.e., when we start on the level with, respectively, highest or lowest slopes. Then, according to chen-LZ , transition probabilities to other states are described by q-deformed binomial distributions. For example, let us introduce parameters

[TABLE]

We know then from cQED-LZ that

[TABLE]

[TABLE]

Let us rewrite these elements in terms of $x$ and $N_{B}$:

[TABLE]

[TABLE]

In what follows, to simplify expressions, we will consider only the case with $N_{B}=0$. Then

[TABLE]

At $N_{B}=0$ we have couplings

[TABLE]

Equations (31) and (III.4) allow us to derive $S_{22}$ and $S_{N-1,N-1}$ explicitly for any DTCM chain:

[TABLE]

For $N=4$, i.e., $S=3/2$, this is sufficient to reconstruct all elements of the transition probability matrix using the fact that this matrix is doubly stochastic and that for bipartite models it is symmetric:

[TABLE]

[TABLE]

We checked that this solution coincides with the one derived in cQED-LZ .

Let us now look at previously unexplored case with $N=5$, i.e., $S=2$. The following elements of the transition probability matrix are known from cQED-LZ :

[TABLE]

Using the 2nd order HCs, we have

[TABLE]

[TABLE]

The trace formula (20) gives

[TABLE]

from where we find

[TABLE]

The remaining three elements $P_{23}$, $P_{24}$, $P_{34}$ can be again obtained using the fact that summation of elements of each column or each row of the transition probability matrix is equal to one. We found:

[TABLE]

Figure 6 shows perfect agreement of Eqs. (III.4)-(III.4) with results of numerical simulations.

The ansatz (33) was not needed to solve $N=4$ and $N=5$ sectors. Rather we used Eqs. (17), (20), and (31) that are common to all bipartite models. For the higher dimensional DTCM sector with $N=6$, however, counting free parameters shows that even with additional input about phases and transition probabilities from fully polarized states, we cannot reconstruct all elements of the transition probability matrix. Apparently, unlike the bowtie model, simple discrete symmetries plus HCs are already not sufficient to obtain complete analytical understanding of DTCM. So, it seems that validity of the ansatz (33) is important but this is not the last piece of the puzzle about the origin of LZ-integrability.

IV Stokes phenomenon in MLZ models

IV.1 General case

Evolution equation with an MLZ Hamiltonian corresponds to a system of ordinary differential equations with a single irregular point in a complex plane at $t\rightarrow\infty$. It has been discussed previously no-go ; shytov ; joye-LZ that since the MLZ Hamiltonian has no singular points at finite complex time values, it is possible to deform the time-evolution contour so that it connects real $t\rightarrow\pm\infty$ limits with a contour that goes along a semicircle with radius $R\rightarrow\infty$. Deformations of time contour in both upper and lower halves of the complex plane are allowed.

The advantage of this deformation is that if the integration contour always has $|t|\rightarrow\infty$ then the diabatic levels are always well separated, and the eigenvectors of the Hamiltonian always coincide with diabatic states. Adiabatic approximation predicts the amplitudes of solutions

[TABLE]

where $\epsilon_{n}$ is the adiabatic energy of the $n$-th eigenstate of the Hamiltonian. Moreover, one does not need to know exact expressions for eigenvalues and eigenvectors because only nonvanishing in $|t|\rightarrow\infty$ limit terms contribute substantially to Eq. (73). Thus, eigenvectors coincide with diabatic states and, for MLZ models with all levels crossing at one point, we can approximate

[TABLE]

For definiteness, we will assume that $t_{0}$ in (73) is taken at infinite negative real values.

Asymptotic functions (73), however, cannot describe the complete solution along the contour that goes near the irregular point (at infinity in our case), which is the essence of the Stokes phenomenon math ; aoki-LZ . For example, the true solution must be a single-valued function of $t$ but, upon integrating along the full circle ${\bm{C}}$, Eq. (73) acquires the factor: $\psi_{n}(te^{2\pi i})=e^{2\pi\eta_{n}}\psi_{n}(t),$ $n=1,\ldots,N$, where we introduced

[TABLE]

So, solution cannot behave as in Eq. (73) everywhere along the contour ${\bm{C}}$. Let us define anti-Stokes lines note along which solutions (73) do not grow or decay exponentially in the limit $|t|\rightarrow\infty$. Stokes sectors are defined then as parts of the complex plane that contain one and only one anti-Stokes line, as we additionally explain in Fig. 7. It is an established mathematical fact that it is always possible to find a set of $N$ solutions with asymptotic behavior (73) everywhere within such a sector deform-book .

The problem is that knowledge of the leading term of the solution in one Stokes sector is not sufficient to find asymptotic behavior of the solution in the full neighboring sector. What can be shown, however, is that such a continuation is described by transition matrices associated with the Stokes lines note . Along these lines, the functions (73) have the highest or lowest rates of the amplitude growth when approaching the singular point at infinity.

Figure 7 shows that there are four Stokes lines in MLZ models given by $t=Re^{i(-3\pi/4+(k-1)\pi/2)}$, where $R\in(0,+\infty)$, and $k=1,2,3,4$, and there are four anti-Stokes lines that are negative and positive parts of the real and imaginary axes. Thus, the covering of the complex plane by Stokes sectors is relatively simple: there are only four of them, marked by $I-IV$. Note that neighboring sectors overlap along Stokes lines.

Let us now assume that at $t=-\infty$ on the real time axis the state vector is a sum of solutions (73) taken with some amplitudes $a_{n}^{I}$. In sector-I this solution can be analytically continued, and along the contour ${\bm{C}}$ with infinite radius we have

[TABLE]

where $\psi^{I}_{n}(t)$ are solutions of the full Schrödinger equation that have asymptotics (73) in the whole sector-I. Our goal is to find the scattering matrix for evolution from $t=-\infty$ to $t=+\infty$ on the real time axis, which is essentially the same as to find how the function $\psi(t)$ defined in (76) behaves near the point of crossing of the positive real axis with the contour ${\bm{C}}$, which is inside sector-III. To get there along the contour ${\bm{C}}$, we should first find solution in the intermediate sector-II.

Solutions with asymptotic behavior (73) that can be analytically continued throughout the whole sector-II have different subleading contributions from functions $\psi^{I}_{n}(t)$ at the same time values in the region where the two sectors overlap. Therefore, the vector amplitudes in overlapping sectors must be related:

[TABLE]

and in sector-II, the state vector is given by

[TABLE]

where $\psi^{II}_{n}(t)$ are solutions of the Schrödinger equation that have leading asymptotic values (73) everywhere inside sector-II.

More generally, after the time-integration contour crosses a $k$-th Stokes line, we can switch to the basis of states that have leading behavior (73) everywhere inside $(k+1)$-st Stokes sector. Amplitudes in the previous and new bases become connected by some transition matrix $\hat{S}_{k}$. So, $\psi(t)$ defined in (76) can be rewritten as

[TABLE]

In (79) and later we will identify indexes $I,II,III,...$ with $1,2,3,...$.

Matrices $\hat{S}_{k}$ are nontrivial, i.e., there is no known general way to write them in a closed form for an arbitrary ordinary differential equation with time-dependent coefficients, including for high order equations that describe MLZ models. However, several properties of these matrices can be rigorously proved.

First, from the fact that solutions $\psi^{k}_{n}(t)$ and $\psi^{(k+1)}_{n}(t)$ in the neighboring sectors can be simultaneously analytically continued in the overlapping region between sectors $k$ and $(k+1)$, it follows that matrix $\hat{S}_{k}$ is time independent. Second, since such two solutions with the same leading exponents along the Stokes line are different only by asymptotically subdominant contributions in the overlapping region, diagonal components of $\hat{S}_{k}$ are all units. Third, one can show berry-stokes that inside the $k$-th Stokes sector if one moves along the contour ${\bm{C}}$ then new sub-leading exponents in solutions $\psi_{n}^{k}(t)$ appear sharply in the vicinity of the $k$-th Stokes line. Therefore, further along this contour, amplitudes of terms with smaller ramp rate along the Stokes line will be influenced by amplitudes of terms with the higher rates but not vice versa. Hence, matrices $\hat{S}_{k}$ have a triangular form, i.e., they have all zero entries either above or below the main diagonal when state indexes follow the order of the level slopes. More rigorous proof of this fact can be found in deform-book .

In our MLZ models, along $L_{1}$, we have leading asymptotics of diabatic states

[TABLE]

Since we assume that $\beta_{1}>\ldots>\beta_{N}$, we have that, along $L_{1}$, the state with index $n$ dominates over states with higher indexes $n^{\prime}>n$, so the corresponding transition matrix has the following triangular form:

[TABLE]

with some constant parameters $x_{ij}$, $i,j=1,\ldots,N$. Along Stokes line $L_{2}$, basis states have opposite ramp rate order, i.e.,

[TABLE]

so the transition matrix is again triangular but has nonzero elements in the upper angle:

[TABLE]

Analogously, we have:

[TABLE]

with some parameters $y_{ij}$, $i,j=1,\ldots,N$.

Note also that if we go around the closed contour ${\bm{C}}$, we should obtain the original state, so $\psi(t)_{t=-\infty}=\sum_{n=1}^{N}a_{n}^{V}\psi_{n}^{V}(t)=\sum_{n=1}^{N}a_{n}^{I}\psi_{n}^{I}(t)$. Using Eq. (79) we find

[TABLE]

where the last factor on the right hand side follows from the integral in (73) over the full cycle ${\bm{C}}$.

IV.2 Scattering matrix

Similarly to Eq. (104), the desired scattering matrix for transition from large negative to large positive real time values is given by

[TABLE]

From Eq. (105), we find that the number of unknown parameters on the right hand side is lower than the number of entries of the scattering matrix $\hat{S}$. So, if the scattering matrix is known, we can reconstruct Stokes matrices $\hat{S}_{k}$. For example, comparing the last columns and rows in both sides of Eq. (105) we have

[TABLE]

where $j=1,\ldots,N-1$. Comparing the next to the last columns and rows and using (106) we can find parameters $x_{N-1,j}$ and $x_{j,N-1}$, and so on. At each step of such calculations we would have to solve only linear equations although final dependence of parameters $x_{ij}$ on $S_{ij}$ would be nonlinear. So, the first feature of Stokes phenomenon in MLZ systems is that knowledge of the scattering matrix determines all independent elements of the Stokes matrices. This is not typical for systems of ordinary differential equations, and follows from relative simplicity of covering of the complex plane by Stokes sectors in MLZ models.

The second interesting property follows from the fact that the number of independent parameters $x_{ij}$ is even lower than the number of elements of the scattering matrix. There are $N$ remaining equations for diagonal elements of matrices in both sides of Eq. (105). These equations produce nontrivial constraints on scattering amplitudes, i.e., on the physically important characteristics. For example, the first pair of such constraints reads

[TABLE]

Comparing Eqs. (107) and (30), we find that Eqs. (107) are nothing but the HCs, which were derived for MLZ theory in constraints by a different approach (expression for $S_{NN}$ has been known for much longer time and is often called the Brundobler-Elser formula be ).

Although we merely reproduced already known Eq. (30), its derivation in this section by exploring the Stokes phenomenon has useful consequences. First, so far, our discussion of transition matrices $\hat{S}_{k}$ has not used unitarity of evolution along the real time axis. Therefore, formulas such as (107) are valid even if the Hamiltonian is not Hermitian. Second, unlike the derivation in constraints , it is now clear that there are no other similar nontrivial constraints in MLZ theory, unless we impose additional symmetries on the model.

IV.3 Case of bipartite crossing at one point

Let us now consider the unitary evolution of a bipartite system of levels crossing at one point. From (17) and (104)-(105) we get

[TABLE]

or

[TABLE]

As we discussed in the previous subsection, there is a unique way to relate matrices $\hat{S}_{3}$, $\hat{S}_{4}$ with $\hat{S}$, and consequently with $\hat{S}_{1}$ and $\hat{S}_{2}$. Comparing (109) with (105) we find

[TABLE]

or in components

[TABLE]

As an example, let us look at the 3-state bowtie model with $\beta_{1}>\beta_{2}>0$:

[TABLE]

Evolution with this Hamiltonian is described by the scattering matrix hioe ; bow-tie

[TABLE]

where we denoted

[TABLE]

In the first equality in (122) we used the symmetry (17) of bipartite models, and in the second equality we wrote the known explicit solution, in which we assumed that the gauge for diabatic states is chosen to remove nontrivial phases of scattering amplitudes. Matrix $\hat{\eta}$ has components

[TABLE]

Substituting (122) and (124) into (105) we find equations for $x_{ij}$, which are straightforward to solve:

[TABLE]

and from (110) we have

[TABLE]

IV.4 Dual semi-infinite MLZ models

As an application of the above formalism, we will show now that with any $N$-state bipartite model we can identify a dual to it MLZ model that describes interactions of $N-1$ semi-infinite Landau-Zener chains. If one of these models is solvable then we automatically know solution of its dual one.

First, we note that the scattering matrix between oscillatory non-decaying asymptotic functions is well defined in MLZ models because both negative and positive pieces of the real time axis coincide with anti-Stokes lines. What is also special about MLZ models whose diabatic levels cross in one point is that there is another, namely the imaginary time, axis with the same property. Therefore, we can similarly define the scattering matrix $\hat{S}^{\prime}$ for evolution from $t=-i\infty$ to $t=i\infty$ along the imaginary time.

Consider evolution (1) along the imaginary time axis. It is convenient then to return to real time by replacing $t=i\tau$, where $\tau\in(-\infty,+\infty)$ is real. We find then

[TABLE]

where we again assume that $\hat{B}={\rm diag}(\beta_{1},\ldots,\beta_{N})$ is nondegenerate diagonal and $\hat{A}$ is a Hermitian matrix with zero diagonal elements.

Imagine that at $\tau\rightarrow-\infty$ the state vector $\psi(\tau)$ is normalized, $|\psi(\tau)_{\tau\rightarrow-\infty}|=1$, and that its oscillating part behaves as $\psi(\tau)_{\tau\rightarrow-\infty}\sim e^{i\beta_{n}t^{2}}$. The element of the scattering matrix, $S_{mn}^{\prime}$, is then defined to be the time-independent prefactor of the amplitude of the solution term that oscillates as $\sim e^{i\beta_{m}\tau^{2}}$ at $\tau\rightarrow+\infty$.

Let us assume that matrix $\hat{A}$ describes bipartite couplings, i.e., that diabatic states can be split into two groups, $G_{1}$ and $G_{2}$, so that states of one group are coupled directly only to states of the other group. Then, according to previous discussion in this section, we can determine transition matrices through Stokes lines in terms of scattering matrices of the evolution with MLZ Hamiltonian $\hat{H}=\hat{B}t+\hat{A}$.

In terms of transition matrices, scattering matrix for evolution (141) is given by

[TABLE]

To see this, we note that if $t_{0}$ in (73) were chosen to be on the negative imaginary axis, we would be able to repeat previous arguments leading to the formula analogous to (105). However, changing $t_{0}$ to be real and negative results in additional factors: $\hat{S}_{k}\rightarrow e^{\pi\hat{\eta}/2}\hat{S}_{k}e^{-\pi\hat{\eta}/2}$, which account for the difference between (142) and (105).

Thus, we conclude that knowledge of the solution of some bipartite MLZ model automatically results in knowledge of the scattering matrix of the model (141). The latter, however, describes non-Hermitian evolution, which is formally not an MLZ model. Nevertheless, consider the following secondary quantized quadratic Hamiltonian of $N$ interacting bosonic modes:

[TABLE]

This Hamiltonian is Hermitian, and in the Fock space of population states $|n_{1},\ldots,n_{N}\rangle$ it has the form (1). Unitary evolution with $\hat{H}_{B}$ conserves the difference between the total number of bosons between the two groups. So the phase space splits into invariant subspaces that are different by the value of this conserved number. States in each subspace can be marked by $N-1$ discrete indexes.

Let us introduce the following vector of operators

[TABLE]

and consider its evolution with the Hamiltonian (IV.4) in the Heisenberg picture:

[TABLE]

This equation is equivalent to

[TABLE]

where $\hat{H}^{\prime}$ has the same form as defined in Eq. (141). This observation means that if the scattering matrix $\hat{S}^{\prime}$ is known for evolution equation (141), we automatically know asymptotic values of operators in the Heisenberg picture of the model (IV.4). For example,

[TABLE]

where we identify $\hat{a}_{k}$ and $\hat{b}^{\dagger}_{s}$ with, respectively, $\hat{a}_{k}(-\infty)$ and $\hat{b}^{\dagger}_{s}(-\infty)$. From such expressions we can find arbitrary correlators of operators in the final state. For example, let us assume that the initial state is characterized by well defined number of bosons of each atomic type. The final population of the bosonic mode $n\in G_{1}$ is then given by

[TABLE]

where the averaging is over the initial state at $t\rightarrow-\infty$, so cross-correlators average to zero. Knowledge of such correlators can be then used to derive transition probabilities between any pair of diabatic states of the Hamiltonian (IV.4), as it was done in gurarie-LZ for the simplest case with $N=2$.

Summarizing, any bipartite MLZ model of the form (1), has a dual counterpart with the Hamiltonian (IV.4). The latter is also an MLZ model but with the infinite number of interacting states. Scattering matrices in dual models are related by equations (105) and (142). Namely, assuming knowledge of the scattering matrix in any one of these equations, we can derive triangular matrices $\hat{S}_{k}$, $k=1,2,3,4$ given by (85), (91) and (97), and then obtain the scattering matrix from the other equation.

IV.5 Coherent conversions between molecular and atomic condensates

Models of the type (IV.4) have actually been encountered in physics in relation to atomic condensate creation, which is induced by applying a linearly time-dependent magnetic field, from a molecular condensate during the passage through the Feshbach resonance gurarie-LZ ; yuzbashyan-LZ ; kayali-LZ . The number of molecules in the molecular condensate is assumed macroscopic so that this number is absorbed in the definition of coupling constants in the mean field approximation. Coupling terms in (IV.4) correspond then to quantum coherent chemical reaction channels:

[TABLE]

meaning that molecules $[AB]$ can split into pairs of atoms $A$ and $B$. Single atoms may have internal states marked by indexes $k$ and $s$. When the external field is far from the resonance, this reaction is suppressed. By sweeping the magnetic field linearly throughout the resonance, we can control the number of single atoms that we produce under condition that this sweep is much faster than the lifetime of the molecular condensate. Internal states can be different localized modes in a trap that interact differently with the linearly time-dependent external magnetic field. This leads to differences in slope parameters in (IV.4). So far, only simplest cases without internal states of the same atoms have been studied yuzbashyan-LZ ; kayali-LZ . Hence, below we will use duality to explore the effect of competition between different channels of this quantum coherent chemical reaction.

For demonstration, we consider the case when the first atomic field has no internal states interacting with the external field while the second atom has two such states. The Hamiltonian of our model is

[TABLE]

where we will assume that $\beta_{1}>\beta_{2}>0$.

The corresponding dual model is the three-state bowtie model, whose scattering matrix is given in (122) and transition Stokes matrices are given in Eqs. (128)-(140). Substituting them in (142) we find

[TABLE]

Let us now look at the simplest and experimentally most relevant case when all atomic modes are initially unpopulated. Final average numbers of atoms in each of the three atomic types then read:

[TABLE]

These equations show that new reaction channels make multiplicative effect on the number of produced single atoms and that relative importance of different reaction channels is exponentially sensitive to couplings and the chemical potential crossing rates.

V Discussion

We considered the special but important subclass of MLZ models with all diabatic levels crossing in one point. Our goal was to broaden the scope of exact results by exploring specific features that are common to all known solvable models. We identified several such features.

First, it is the bipartite structure of all solvable models, i.e., all of them describe interactions between two groups of states, such that direct couplings between diabatic states within the same group are zero. We found that this simple symmetry leads to Eqs. (17), (20), and (31) that are not encountered in other MLZ models and constrain elements of the scattering matrix considerably.

Our second finding is the additional symmetry (33), which is rather the ansatz with phases of scattering amplitudes chosen to be removable by gauge transformations. Both analytical and numerical studies showed that this ansatz worked for all known solvable models.

One consequence of the validity of this ansatz is that the scattering matrix can be parametrized by a set of orthonormal vectors with real entries. In known solvable models, the latter are always polynomials of some parameter combinations. So, we can speculate that this property relates MLZ-integrability to the theory of orthogonal polynomials, whose importance for mathematical physics has been well recognized aborodin . Moreover, all independent hierarchy constraints (30), which generally tell very little about the transition probability matrix, now lead to direct constraints on it. In the bowtie model, these constraints are sufficient to derive the complete transition probability matrix. However, in the case of the driven Tavis-Cummings model we achieved only limited success by solving two low-dimensional sectors. Hence, there should be other symmetries that are to be discovered in order to understand the integrability of the DTCM.

Finally, we showed that the Stokes phenomenon in all MLZ models has several simplifying properties that allow determination of all nontrivial Stokes multipliers from the scattering matrix. For bipartite models this leads to duality with a new solvable class of MLZ systems. Models of this class describe interactions of molecular and atomic Bose condensates in time-dependent fields.

Our results definitely do not exhaust all interesting properties of bipartite models. It seems to us that future progress will be achieved by understanding how and why various discrete and dynamic symmetries constrain Stokes multipliers in MLZ systems. Achieving this goal will enable discoveries of new classes of explicitly time-dependent solvable systems.

Appendix A Composite solvable model with particle-hole symmetry

First, we note that (11), is related to the Hamiltonian that describes interaction of three fermionic modes:

[TABLE]

We will assume that $\beta_{1}>\beta_{2}>\beta_{3}>0$. For evolution with $\hat{H}^{f}$, Heisenberg equations for annihilation and creation operators read:

[TABLE]

where $\eta_{nm}=1$ if $m<n$ and $\eta_{nm}=0$ if $m>n$. Equations (A) for operators are the same as equations for amplitudes in the model (11), so if we know the scattering matrix $\hat{S}$ for the model (11) we can connect also values of operators at $t=\pm\infty$:

[TABLE]

where we identify $\hat{c}_{i}\equiv\hat{c}_{i}(-\infty)$, e.t.c.. Due to the particle-hole symmetry, evolution equations (A) satisfy conditions of the No-Scattering theorem in li-16LZ , according to which elements of the scattering matrix $\hat{S}$ in (175) along the second main diagonal are all zero:

[TABLE]

Next, consider evolution with the Hamiltonian (A) in the sector that contains the vacuum state. Other orthogonal states in this sector include:

[TABLE]

In the basis (176) the Hamiltonian (A) reads:

[TABLE]

It describes the four-state bowtie model, whose expression for the transition probability matrix is known exactly bow-tie :

[TABLE]

where

[TABLE]

We can now express transition probabilities in the model (177) via transition probabilities in the model (11). For example, let initial state be $|0\rangle$ and we consider the probability to find the fermion $\hat{c}_{1}$ at the end of evolution:

[TABLE]

On the other hand, this is the probability of that, in the four-state bowtie model, the final state is either $|12\rangle$ or $|13\rangle$:

[TABLE]

where the upper index “$bt$” means that this is the element of the transition probability matrix of the bowtie model (177). Similarly, we find

[TABLE]

[TABLE]

From Eqs. (180)-(183), we derive

[TABLE]

Considering other initial conditions, we can reconstruct all other elements of the probability matrix in the six-state model with the Hamiltonian (11):

[TABLE]

Appendix B Solvable five-state bipartite model

The general bipartite five-state model has the Hamiltonian

[TABLE]

where we used the fact that one can fix slopes of any two levels by gauge transformation and time-rescaling. To find parameters at which this model is solvable, we follow the same steps as in solution of the four-state model in Eq. (6), which was developed in Quest-LZ . First, we add different constant parameters to the diagonal of the Hamiltonian, so that the pattern of diabatic level crossings looks like Fig. 1(a). Next, imposing integrability conditions Quest-LZ , we determine the constraints on parameters that make this model solvable. We then apply the semiclassical ansatz, as it is explained in Quest-LZ in order to derive transition probabilities. Finally, we set constant parameters at the main diagonal of the Hamiltonian to zero and obtain the model of the type (186) with parameters that correspond to the already known transition probability matrix.

This whole procedure is straightforward but quite lengthy, so we will only list the final results. We found that solvable models, up to redefinition of indexes, correspond to level slopes $b_{3}>b_{4}>b>0>-b>b_{5}$, and the couplings satisfying the following constraints:

[TABLE]

where $\tau_{i}$ is the sign of $g_{2i}/g_{1i}$ which can be taken as either 1 or $-1$. Depending on the signs of $\tau_{i}$’s, there are three phases in this model that correspond to different forms of transition probability matrices. Let

[TABLE]

For $\tau_{3}=\tau_{4}=\tau_{5}=1$, we found

[TABLE]

for $\tau_{3}=\tau_{4}=-\tau_{5}=1$:

[TABLE]

and for $\tau_{3}=-\tau_{4}=-\tau_{5}=1$:

[TABLE]

Figure 8 compares some of these predictions to results of our numerical simulations. The agreement is excellent.

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. and F.L. also thank the support from the LDRD program at LANL.

F. Li and C. Sun made equal contributions to this article.