Nonadiabatic transitions during a passage near a critical point

For finite $k$, the exact analytical analysis becomes impossible, so we consider the limit of a very slow transition: $g^{2}/\beta\gg 1$. Near the critical point and for $N\gg 1$, the all-to-all interactions justify the applicability of the mean field approximation, so that the effective Hamiltonian for the $i$-th spin is given by

where $S_{z}(t)$ is the mean polarization field of all spins:

The evolution starts as $t\rightarrow-\infty$ in the ground state with $\langle\sigma_{i}^{z}\rangle=1$ for all $i$. For the adiabatic evolution, as $t\rightarrow+\infty$, all spins must flip to the states with $\langle\sigma_{i}^{z}\rangle=-1$. For intermediate time, the adiabatic evolution makes the spin always point against the direction of the net magnetic field that has components $(B_{x},B_{y},B_{z})=(g,0,\beta t-kS_{z}(t))$. This gives us the value of $\langle\sigma_{i}^{z}\rangle$ in the strict adiabatic limit:

Substituting this into (5) we find the equation that defines $S_{z}(t)$ in the quasi-adiabatic case:

We consider the regime in which the nonadiabatic excitations are strongly suppressed, so we disregard their effect on $S_{z}(t)$.

For each spin, the probability of the nonadiabatic excitation can be estimated now from a spin-1/2 Hamiltonian (4), in which $S_{z}(t)$ is defined implicitly by Eq. (6). Let us now change the time variable in (6) to

so that

By differentiating, we find

Using this, we change the variables, from $t$ to $x$, in the time-dependent Schrödinger equation

where $H_{i}$ is defined in (4) and write this equation in the matrix form

where $a$ and $b$ are the amplitudes of, respectively, the spin up and down states, and

Equation (8) is very similar to the Landau-Zener model but it has an additional $x$-dependent factor, $\alpha(x)$, in front of the 2$\times$2 Landau-Zener matrix Hamiltonian. This difference, however, is crucially important.

For $k>g$ there is a range of $x$ for which $\alpha(x)<0$. From Eq. (7) this means that the function $x(t)$ is no longer monotonically growing, and hence the change of the time-variable is ill-defined. This problem can be traced to the fact that for $k>g$ Eq. (6) that defines $S_{z}(t)$ has more than one solution for some $t$.

This means that the original model (2) for $k>g$ describes the passage through a 1st order phase transition at some time moment. The nonadiabatic excitations in the nLZ model in that case are abundant and have been already well studied for the nLZ model niu-nlz . Hence, here we consider only the situation with $k\leq g$, which avoids this phase transition, and for which our mean-field approximations leading to Eq. (8) in the quasi-adiabatic limit are justified. Our primary interest, however, will be in the passage near the phase transition, which corresponds to $g-k\ll g$. A special value, $k=k_{c}=g$, corresponds to the control protocol that merely touches the critical point during the time evolution. In this case, Eq. (7) still has a unique solution for $t\in(-\infty,+\infty)$ but $S_{z}(t)$ changes so quickly near $t=0$ that the adiabaticity conditions cannot be satisfied even for $\beta\rightarrow 0$.

The Hamiltonian in Eq. (8) has eigenvalues

In the Dykhne’s approach, the adiabatic evolution is extended to the complex time so that we pass through the exact eigenenergy degeneracy point, at which $\varepsilon_{+}=\varepsilon_{-}$. If $\alpha(x)$ in (8) were a constant, as in the standard Landau-Zener model ($k=0$), this point would be the branch cut at $x=ig$. However, for $k>0$ at this point $\alpha(x)$ diverges. Instead, the energy degeneracy is found now when the prefactor $\alpha(x)$ is zero. Let us introduce a parameter

that characterizes the relative strength of the spin interactions. Then, $\alpha(x_{0})=0$ is satisfied at

For the moderate ferromagnetic interactions, $1>\kappa>0$, this zero is closer to the real axis than the point of the branch cut of the linear Landau-Zener model, at $x=ig$. This is the sign for the enhanced nonadadiabatic effects. At $\kappa=1$, which is the critical value, the degeneracy point touches the real time axis.

Since near $x=x_{0}$ the adiabaticity is broken, we expect that the probability of the nonadiabatic excitation, which in our case is to remain in the spin up state as $t\rightarrow+\infty$, is given by the adiabatic evolution along the time contour that passes through the degeneracy point:

Equation (13) is known as the Dykhne formula dykhne . Dykhne originally considered the typical situation with real parameters and the energy difference in (13) having a branch cut at $x=x_{0}$, as in the standard Landau-Zener Hamiltonian ($k=0$) for $x\approx x_{0}$: $\varepsilon_{+}-\varepsilon_{-}\propto\sqrt{x-x_{0}}$. For such cases, not only the leading exponent was found. Dykhne proved that the exponential prefactor in Eq. (13) at the leading order in small $\beta$ is generally $1$, as confirmed e.g. by the solution of the Landau-Zener model.

In our case, there is a complication. Although $x=x_{0}$ is the point of the energy degeneracy, this degeneracy is not a branch cut point. Rather the entire Hamiltonian

is a simple zero at this point. Namely, for $x\rightarrow x_{0}$,

The energy difference for $H(x)$ at $x_{0}$ has also a similar simple zero. By investigating Eq. (10), with $\alpha(x)$ given in Eq. (9), near $x=x_{0}$ given in Eq. (12), we find

The effect of more complex degeneracy points on the Dykhne formula was explored by Joye in joye-prj . His theory says that if near the energy degeneracy point the energy difference behaves as $\varepsilon_{+}-\varepsilon_{-}\propto(x-x_{0})^{n/2}$, and the entire Hamiltonian behaves as $H(x)\propto(x-x_{0})^{m}$, then the exponent in the Dykhne formula has a prefactor

Given Eq. (15), we have $n=2$ and $m=1$, so the Joye’s formula predicts the zero prefactor, which would mean no nonadiabatic transitions. However, Joye interpreted this fact in joye-prj , so that a further theory was needed for this particular case. Hence, we are to develop this theory in order to derive the probability of the nonadiabatic excitations for the evolution with the Hamiltonian (14).

We will need the Schrödinger equation with $H(x)$ in the adiabatic basis. Let $|+\rangle$ and $|-\rangle$ be the eigenstates that correspond to the $\varepsilon_{\pm}$ eigenvalues. This basis is $x$-dependent, so after switching the basis, the adiabatic eigenstates of $H(x)$ are coupled by the time derivative: $i\langle+|\partial_{x}|-\rangle=-ig/[2(x^{2}+g^{2})]$, so the Schrödinger equation is given by

where

We assume that we start as $x\rightarrow-\infty$ in the ground state, so in this basis

and as $x\rightarrow+\infty$ the state behaves as

where $B_{\pm}$ are still unknown amplitudes and $\varepsilon_{\pm}=\pm\varepsilon/\beta$. The excitation probability is identified with $|B_{+}|^{2}$.

Note that the coefficient $1/\beta$ is large because we are interested in the quasi-adiabatic transition, which corresponds to a small transition rate $\beta$. This makes the diagonal part of the Hamiltonian in (17) generally large. The region with the nonadiabatic transitions then is crossed during a short time, which becomes shorter with decreasing $\beta$. The transition probability between the adiabatic states can then be, naively, estimated using the standard lowest nonzero order Born approximation. According to it, we treat the off-diagonal part of the matrix (17) as a small perturbation and calculate the amplitude of the initially empty state up to the first order in this perturbation. The excitation probability is the square of this amplitude:

It is well known that such a Born approximation to $P_{ex}$ generally fails to predict the excitation probability when the latter is suppressed exponentially, e.g., as $P_{ex}\sim e^{-f/\beta}$. The reason is that the standard Born series is intended to produce the contributions to $P_{ex}$ as powers of $\beta$. Therefore, formally different order terms in the Born series have to be suppressed additionally and thus can produce comparable contributions to $P_{ex}$.

The break down of the lowest Born approximation can be illustrated within the standard Landau-Zener model ($k=0$) for which the exponential prefactor in (24) at the saddle point of the integral, at $x=ig$, diverges. By calculating such an integral, one would find a prefactor to the exponent in (3) different from the exactly known value, $1$. However, in the context of the properties of the point $x_{0}$ that we encounter here, there is no such a singularity, so this formula is justified when it is applied to a piece of the integration time-path that goes through $x_{0}$.

Consider first the situation with touching the critical point, which takes place for $\kappa=1$ at $x_{0}=0$. The nonadiabatic transitions are mainly produced near the energy degeneracy point, so we can safely approximate near this point $\varepsilon(x)\approx\frac{3x^{2}}{2g}+O(x^{3}),$ and $g/(x^{2}+g^{2})\approx 1/g+O(x^{2})$. The Schrödinger equation then simplifies to

The lowest Born approximation then predicts

where $\Gamma(\ldots)$ is the Euler’s Gamma-function. The lowest order approximation leads in (26) to the power law, $P_{ex}\sim\beta^{2/3}$, rather than the exponential suppression of $P_{ex}$. Therefore, this approximation was well justified, and the higher order terms in the Born series would produce only the corrections of higher power in $\beta$. Therefore, the numerical coefficient

can be trusted.

In Fig. 1, we compare the numerically obtained solution for $P_{ex}$ (black dots) with our analytical predictions. In particular, in Fig. 1(b), the solid curve for $\kappa=1$ represents the prediction of the formula (26). The agreement with the numerical results is excellent.

Finally, we note that in the same regime Ref. niu-nlz predicted $P_{ex}\sim\beta^{3/4}$. The value $3/4$ is close to our $2/3$ but different. This discrepancy may be the result of the difference in the way we developed the mean field approximation.

For the differential equation (8), the piece of the imaginary time axis that goes from zero to the branch cut at $x=ig$ is the Stokes line, along which the solution amplitudes are changing exponentially quickly. The strict adiabatic approximation generally fails near this line. One possibility to derive the transition amplitude between the adiabatic states then is to identify the region at the Stokes line such that the behavior of the time-evolution across it simplifies. Inside this region, the evolution can be described analytically beyond the adiabatic approximation, and then it can be smoothly connected to the adiabatic states from both sides of the Stokes line joye-prj . In our case, this region is found near the energy degeneracy point at $x=x_{0}$.

Equation (8) has a single analytical solution on the complex plane with branch cuts starting at $x=\pm ig$ and going to the infinity along the imaginary time axis. Therefore, it is safe to deform the integration path from the ${\rm Re}(x)$-axis to the path shown by thick arrows in Fig. 2. The points $x_{\pm}$ along this path have the same imaginary value as the point $x_{0}$ in (12). The integration path turns from the real axis to reach the point $x_{-}$ along the vertical arrow. Then it goes through the point $x_{0}$ towards $x_{+}$, and then returns to the real time axis.

The vertical pieces are chosen to be sufficiently far from the point $x_{0}$, so that one can apply the adiabatic approximation. For the points that lie on the path before $x_{-}$ this leads us to an adiabatic amplitude $|\psi(x)\rangle\sim e^{-i\int^{x}_{-\infty}\varepsilon_{-}(\tau)\,d\tau}$. After passing through the degeneracy point, another asymptotic solution that behaves as $\sim e^{-i\int^{x}_{x_{+}}\varepsilon_{+}(\tau)\,d\tau}$ acquires a finite amplitude, and its magnitude at the real axis is what determines the final excitation probability.

Along the arrow that connects $x_{-}$ and $x_{+}$, the adiabaticity is broken, so our goal now is to find the absolute value squared of the amplitude of the new asymptotic solution that emerges by reaching the point $x_{+}$. By linearizing the time-dependence of the Hamiltonian parameters in Eq. (17) near $x_{0}$ (see the first term in Eq. (15)) in the adiabatic basis and setting $\tau\equiv x-x_{0}$ we find the effective Hamiltonian along the arrow that points to $x_{+}$:

Due to large values of $\beta$, the adiabaticity is violated in a very short region near $x_{0}$, so it is safe to assume that along the horizontal arrow in Fig. 2 the time $\tau$ changes in the interval $\tau\in(-\infty,+\infty)$. Although $H(\tau)$ is non-Hermitian, the perturbative analysis applies to it equally well, as for the Hermitian evolution. This Hamiltonian has the same property as the one from Eq. (25), for the case with $\kappa=1$. Namely, along the integration path, at $\tau=0$, the energy gap almost closes, while the off-diagonal terms are non-singular at $x_{0}$. Hence, the transition amplitude from the initial adiabatic state at the point $x_{-}$ to the new state with the growing phase asymptotic $\sim e^{-i\int^{x}\varepsilon_{+}(\tau)\,d\tau}$ can be estimated safely, up to an unimportant unitary phase factor, by applying the lowest Born approximation:

Finally, along the down-pointing arrow in Fig. 2, this new state gains an additional adiabatic phase factor $e^{-i\int_{x_{+}}^{x}d\tau\,\varepsilon_{+}(\tau)}$.

Due to the large values of $1/\beta$, the adiabaticity is broken only in the vanishing with $\beta$ vicinity of $x_{0}$. Hence, $x_{\pm}$ in this limit can be placed arbitrarily close to $x_{0}$. Then, for the adiabatic evolution along the vertical arrows, it is safe to disregard the real components of $x_{+}$ and $x_{-}$, so that the contribution of the vertical parts of the integration contour produce the standard Dykhne exponent (13) factor in the excitation probability. The horizontal path that goes through $x_{0}$ contributes with an additional $|A|^{2}$ factor to this probability. By collecting the contributions from all pieces of the integration contour we find

All integrals in (27) can be calculated analytically, so we finally find

where

and where $\kappa$ is defined in terms of the parameters of the nLZ model (2) by Eq. (11). Equations (28) and (29) are our main results. They show explicitly that the excitation probability has the functional form $P_{ex}\propto\beta e^{-f/\beta}$, i.e., it has a nontrivial prefactor that decays to zero as $\beta^{\nu}$. Interestingly, the exponent $\nu=1$ in the prefactor in (28) is different from the exponent $\nu_{c}=2/3$ for the case of touching the critical point at $\kappa=1$. Mathematically, this follows from Eq. (15). Namely, for $1>\kappa>0$ the leading term in the energy difference near $x=x_{0}$ is growing linearly with $x-x_{0}$ but for $\kappa=1$ this term is identically zero, so the next order in $(x-x_{0})$ term controls the energy splitting.

The function $F(\kappa)$ is shown in Fig. 3. As $\kappa\rightarrow 0$, it approaches the limit

which reproduces the exponent in the Landau-Zener formula (3). Near the critical point, $1-\kappa\ll 1$, it behaves as a power law

As expected, $F(1)=0$, which means a transition to the regime with only a power law suppression of the nonadiabatic excitations. Our expression for $F(k)$ differs quantitatively from the one that was derived in Ref. niu-nlz for the nLZ model. Thus, instead of the exponent $\mu=3/2$ in (30), Ref. niu-nlz predicted $\mu=2$. Nevertheless, qualitatively our Fig. 3 is very similar to its analog, Fig. 5 in niu-nlz . Within our mean field approximation, however, we achieved an excellent agreement with exact numerical simulations.

For comparison with numerical results, we noted that Eq. (28) was derived for the limit $\beta/g^{2}\ll 1$, so it made unphysical predictions for the strongly nonadiabatic regime $\beta/g^{2}\gg 1$. To avoid this complication, we replaced the prefactor

by the function

This did not spoil the asymptotic behavior of the prefactor as $\beta\rightarrow 0$, but interpolated it to a known value $P_{ex}(\beta\rightarrow\infty)=1$ in order to avoid unphysical predictions with $P_{ex}>1$.

In Fig. 1, we compare such adjusted predictions of Eq. (28) to the result of the exact numerical integration of the evolution equation (17), which is equivalent to the original Schrödinger equation (8). The agreement is excellent even for the small values of $\kappa$, which are far from the critical point $\kappa_{c}=1$. It seems that for $0<\kappa<1$ the only condition for the validity of Eq. (28) is the smallness of the prefactor: $|A|^{2}\ll 1$, which is guaranteed at sufficiently small ratio, $\beta/g^{2}$. However, for $\kappa\ll 1$, this is achieved for relatively small $P_{ex}$. Hence, for the passage at a larger distance from the critical point the role of the prefactor is smaller.

Finally, we note that our analysis is straightforward to extend to the antiferromagnetic interactions with $\kappa<0$ but the analytical formulas become much more complex. For $\kappa<0$ there are two rather than one zeros of the function $\alpha(x)$ at the same distance from the real time axis. The vicinity of each zero can be treated similarly to the zero at $x_{0}$ for $\kappa>0$. However, in addition, we then must take into account the interference between the contribution of each zero to the transition amplitude. This leads to oscillations in the exponential prefactor. As this regime is far from the critical point $\kappa_{c}=1$, we leave it without further discussion.

Consider a model of a linear chain of sites with quantum coherent hopping of a charged particle. The chain is placed in an external constant electric field. Let $|n\rangle$, where $n=-\infty,\ldots,\infty$, be orthogonal states of the particle, with spatially separated average positions. Let $a$ be the distance between the sites $n$ and $n\pm 1$, as shown in Fig. 4(a). The direct hopping between the neighboring states has the energy amplitude $\Delta$.

The electric field creates a linear potential for the particle with electric charge $e$ along the $x$-axis of the chain: $U(x)=e{\cal E}x$. In the tight-binding model, this corresponds to the energy $e{\cal E}an$ for the particle in the $n$-th node of the chain. Thus, the entire tight-binding Hamiltonian is given by

For any finite electric field, and in the absence of decoherence/relaxation, the eigenstates of the Hamiltonian (31) are the localized Stark states sinitsyn-prb2002 , which have a spectrum with a fixed energy gap between the nearest energy levels. Such a spectrum is called the Stark ladder. Nevertheless, at ${\cal E}=0$, the eigenstates are the delocalized Bloch waves (see Fig. 4(b)). Hence, $b=0$ in (31) is the critical point. If the electric field changes sufficiently slowly without crossing the value ${\cal E}=0$, the predictions of the adiabatic theorem are satisfied and we expect to find an exponential suppression of the excitations. However, if we cross the point ${\cal E}=0$ we expect an abundant production of the excitations even for a quasi-adiabatic evolution.

Let the time-dependent electric field start and end at large absolute values: $|b|_{t\rightarrow\pm\infty}\gg\Delta$ in (31). Then, the Stark states are asymptotically the basis states $|n\rangle$. Hence, the number of the excitations can be identified with the change of the index $n$. Let $n_{0}$ be the initial node of the particle as $t\rightarrow-\infty$ and let as $t\rightarrow+\infty$ there is a probability distribution $P(n)$ to find the particle at the $n$-th node. On the Stark ladder, between the energies of the Stark states $|n\rangle$ and $|n_{0}\rangle$ there are $(n-n_{0})$ elementary energy gaps. Hence, it is natural to define the average number of the nonadiabatic excitations as

where the average is over the final distribution $P(n)$.

The linear dynamic transition through the critical point in the model (31), namely the case with $b=\beta t$, was previously studied in sinitsyn-prb2002 . As for other critical phenomena, the exact solution of this model produced a power law for the number of the excitations:

Let us now consider a different protocol. In what follows, we are interested in the Hamiltonian (31) with

Here, the electric field becomes infinitely large as $t\rightarrow\pm\infty$. However, it is never zero and reaches the minimum $b=g$ at $t=0$. The rate of the passage through this minimum is controlled by the parameter $\varepsilon$. Our goal is to find the scaling of $n_{ex}(\varepsilon)$ for $\varepsilon\rightarrow 0$ after the time evolution during $t\in(-\infty,+\infty)$.

We will search for the solution to the Schrödinger equation of the form

and introduce the amplitude generating function (AGF)

from which we can find all state amplitudes in the original basis by the inverse Fourier transformation:

The AGF is also convenient for finding the moments of the position index $n$. Thus, if we normalize

then the average of the $m$-th power of $n$ is given by

The time-dependent Schrödinger equation for the state amplitudes is

Changing the variables $a_{n}\rightarrow a_{n}e^{-in\left(\frac{\varepsilon^{2}t^{3}}{3}+gt\right)}$, then multiplying the $n$-th equation in (39) by $e^{i\phi n}$, and summing all these equations, we find the first order differential equation for the AGF, $u(t)\equiv u(t;\phi)$:

where $\phi\in[0,2\pi)$ should be treated as a constant parameter.

Let as $t\rightarrow-\infty$ the system start with the particle at the node with $n=0$. The initial condition for Eq. (40) is then $u(\--\infty)=1$. The solution for the AGF as $t\rightarrow+\infty$ then is given by

The integral in (41) can be expressed via the Airy function:

in terms of which

Using Eq. (38), we get

Thus, the number of the nonadiabatic excitations, $n_{ex}$, up to a constant prefactor, is described by the function $x{\rm Ai}[x]$, where $x=g\varepsilon^{-2/3}$. Note that $\varepsilon^{-2/3}$ has the dimension of time. Large values of this time correspond to the quasi-adiabatic regime. Figure 5, confirms the validity of Eq. (43) by comparing the analytical results with the numerical solution of the Schrödinger equation.

In the quasi-adiabatic limit the argument of the Airy function in (43) is large, so we can replace ${\rm Ai}[x]$ by its known asymptotic value

Finally, we find for the quasi-adiabatic evolution

This result shows explicitly that for $g>0$, i.e., as far as there is no crossing of the critical point by the external field, the number of the excitations is suppressed in the adiabatic limit exponentially but with a nontrivial prefactor of the exponent. This prefactor scales as the power law, $\varepsilon^{\nu}$, where $\nu=-1/2$, on the transition rate $\varepsilon$. Equation (44) also reveals the “exponent-in-exponent” behavior. Namely, in our case the parameter $g$ controls the minimal distance to the critical point that is encountered during the process. The factor in the exponent in (44) scales as $\propto g^{\mu}$, where $\mu=3/2$.

We also comment on the case of merely touching the critical point. This corresponds in our protocol to $g=0$. Equation (44) at $g=0$ diverges, which suggests a possible new power law behavior of $n_{ex}$. Indeed, by setting $g=0$ in Eq. (41) we find

This is the power law, $n_{ex}\propto\varepsilon^{\nu_{c}}$, with $\nu_{c}=-2/3$, which is different from $\nu=-1/2$ in the prefactor of Eq. (44).