Time dynamics of Bethe ansatz solvable models

TL;DR

This paper introduces a method to compute the time evolution of Bethe ansatz solvable models by deriving dynamical Bethe equations, enabling exact solutions for models like the Bose-Hubbard dimer and Tavis-Cummings.

Contribution

It develops a novel approach to determine the dynamics of exactly solvable models using time-dependent Bethe parameters and differential equations.

Findings

Derivation of dynamical Bethe equations for time evolution.

Exact solutions for Bose-Hubbard dimer and Tavis-Cummings models.

Extension beyond the Gaudin class of models.

Abstract

We develop a method for finding the time evolution of exactly solvable models by Bethe ansatz. The dynamical Bethe wavefunction takes the same form as the stationary Bethe wavefunction except for time varying Bethe parameters and a complex phase prefactor. From this, we derive a set of first order nonlinear coupled differential equations for the Bethe parameters, called the dynamical Bethe equations. We find that this gives the exact solution to particular types of exactly solvable models, including the Bose-Hubbard dimer and Tavis-Cummings model. These models go beyond the Gaudin class, and offers an interesting possibility for performing time evolution in exactly solvable models.

Equations

$$\prod _{j=1}^M\mathbf{B}({\lambda }_j)|\text{vac}⟩,$$

$$e^{ip(t)}\prod _{j=1}^M\mathbf{B}({\lambda }_j(t))|\text{vac}⟩,$$

$$i\frac{d}{dt}|\Psi (t)⟩=\tau (\lambda )|\Psi (t)⟩.$$

$$|{\Psi }_M^{\sigma }(t)⟩=e^{i{p}^{\sigma }(t)}\prod _{j=1}^M\mathbf{B}({\lambda }_j^{\sigma }(t))|\text{vac}⟩,$$

$$\sum _{\sigma }|\Psi (\{{\lambda }^{\sigma }(t)\})⟩⟨\Psi (\{{\lambda }^{\sigma }(t)\})|\propto \hat{I},$$

$$⟨\Psi (\{{\lambda }^{{\sigma }^{\prime }}(t)\})|\Psi (\{{\lambda }^{\sigma }(t)\})⟩=0,$$

$$⟨\dot{\Psi }(\{{\lambda }^{{\sigma }^{\prime }}(t)\})|\Psi (\{{\lambda }^{\sigma }(t)\})⟩+⟨\Psi (\{{\lambda }^{{\sigma }^{\prime }}(t)\})|\dot{\Psi }(\{{\lambda }^{\sigma }(t)\})⟩=0,$$

$$\tau (\mu )$$

$$+\sum _{n=1}^M{\varphi }_n(\mu ,\{{\lambda }_j^{\sigma }\})\mathbf{B}(\mu )\prod _{\genfrac{}{}{0ex}{}{j=1}{j\ne n}}^M\mathbf{B}({\lambda }_j^{\sigma })|\text{vac}⟩,$$

$$i(i\frac{dp}{dt}-\Theta (\mu ,\{{\lambda }_j^{\sigma }(t)\})\prod _{j=1}^M\mathbf{B}({\lambda }_j^{\sigma }(t))|\text{vac}⟩$$

$$=-i\frac{d}{dt}\prod _{j=1}^M\mathbf{B}({\lambda }_j^{\sigma }(t))|\text{vac}⟩$$

$$+\sum _{n=1}^M{\varphi }_n(\mu ,\{{\lambda }_j^{\sigma }(t)\})\mathbf{B}(\mu )\prod _{\genfrac{}{}{0ex}{}{j=1}{j\ne n}}^M\mathbf{B}({\lambda }_j^{\sigma }(t))|\text{vac}⟩.$$

$$f(\{{\lambda }_j^{\sigma }(t)\})\prod _{j=1}^M\mathbf{B}({\lambda }_j^{\sigma }(t))|\text{vac}⟩=-i\frac{d}{dt}\prod _{j=1}^M\mathbf{B}({\lambda }_j^{\sigma }(t))|\text{vac}⟩$$

$$+\sum _{n=1}^M{\varphi }_n(\mu ,\{{\lambda }_j^{\sigma }(t)\})\mathbf{B}(\mu )\prod _{\genfrac{}{}{0ex}{}{j=1}{j\ne n}}^M\mathbf{B}({\lambda }_j^{\sigma }(t))|\text{vac}⟩,$$

$$p^{\sigma }(t)=-{\int }_0^t\left[i\Theta (\mu ,\{{\lambda }_j^{\sigma }(t^{\prime })\})+f(\{{\lambda }_j^{\sigma }(t^{\prime })\})\right]d{t}^{\prime }.$$

$$i⟨\Psi (\{{\lambda }^{{\sigma }^{\prime }}(t)\})|\dot{\Psi }(\{{\lambda }^{\sigma }(t)\})⟩=$$

$$+⟨\Psi (\{{\lambda }^{{\sigma }^{\prime }}(t)\})|\sum _{n=1}^M{\varphi }_n(\mu ,\{{\lambda }_j^{\sigma }(t)\})\mathbf{B}(\mu )\prod _{\genfrac{}{}{0ex}{}{j=1}{j\ne n}}^M\mathbf{B}({\lambda }_j^{\sigma }(t))|\text{vac}⟩,$$

$$\hat{\mathcal{H}}=$$

$$+\frac{U}{2}\left(a^{†}{a}^{†}aa+b^{†}{b}^{†}bb\right)+V{a}^{†}a{b}^{†}b,$$

$$\hat{H}=-\frac{1}{J}\left(\hat{\mathcal{H}}-\frac{U}{2}\hat{N}(\hat{N}-1)-ϵ\hat{N}\right),$$

$$\hat{H}=\Delta b^{†}b+a^{†}b+a{b}^{†}+c^2{a}^{†}a{b}^{†}b,$$

$$\mathbf{B}(\lambda )=\lambda b^{†}-\mathbf{X},$$

$$\mathbf{X}=\frac{\Delta }{c}{b}^{†}+c{a}^{†}a{b}^{†}+c^{-1}{a}^{†}.$$

$$|\text{vac}⟩=|0{⟩}_a\otimes |0{⟩}_b.$$

$$i\left(\frac{\dot{\Delta }}{c}-{\dot{\lambda }}_n^{\sigma }(t)\right)={\phi }_n(\{{\lambda }^{\sigma }\}){\lambda }_n^{\sigma }(t),\forall n=1,...,N.$$

$${\phi }_n(\{{\lambda }^{\sigma }\})=$$

$$+\frac{1}{c{\lambda }_n^{\sigma }}\prod _{\begin{array}{c}j=1\\ j\ne n\end{array}}^N\left(1-\frac{c}{{\lambda }_j^{\sigma }-{\lambda }_n^{\sigma }}\right).$$

$$\Delta (t)={\Delta }_0+\cos (t^2),$$

$$\nu (t)=\frac{|⟨a^{†}b⟩|}{N},$$

$$-i\frac{{\dot{\lambda }}_n(t)}{{\lambda }_n(t)}={\phi }_n(\{\lambda \})\forall n=1,...,N.$$

$$\prod _{j=1}^N{\mathbf{B}}^{(c^{\prime },{\Delta }^{\prime })}({\lambda }_j(0))|\text{vac}⟩=\prod _{j=1}^N{\mathbf{B}}^{(c,\Delta )}({\lambda }_j^0)|\text{vac}⟩.$$

$${\lambda }_j(0)={\lambda }_j^0-\frac{{\Delta }^{\prime }-\Delta }{c}.$$

$${\bf R}(\lambda,\mu)=\left(\begin{array}[]{cccc}f(\mu,\lambda)&0&0&0\\ 0&g(\mu,\lambda)&1&0\\ 0&1&g(\mu,\lambda)&0\\ 0&0&0&f(\mu,\lambda)\end{array}\right).$$

$$T(\lambda )=\left(\begin{array}{cc}A(\lambda )& B(\lambda )\\ C(\lambda )& D(\lambda )\end{array}\right),$$

Full text

I Introduction

Exact methods of mathematical physics have substantially pushed our understanding of many paramount nonlinear phenomena. One such method is the Quantum Inverse Method (QIM) which was developed almost 40 years ago by Faddeev, Sklyanin, Takhtadzhyan and others Sklyanin et al. (1979); Sklyanin (1982); Kulish and Sklyanin (1982). QIM together with the algebraic version of Bethe Ansatz Slavnov (2018a); Levkovich-Maslyuk (2016) has been successfully applied to various problems from different areas of physics such as one dimensional BECs Lieb and Liniger (1963); Knap et al. (2014), spin chains Maillet (2007); Kitanine et al. (1999); Kato et al. (2003); Bortz and Göhmann (2005), $(1+1)$ models of quantum field theory Faddeev (1982), $(2+1)$ model of classical statistical physics Thiery and Le Doussal (2016), conformal field theory and string theory Arutyunov et al. (2004), quantum optics Bogoliubov and Kulish (2012), and quantum dots Bortz and Stolze (2007).

Obtaining the time dynamics of quantum many-body systems remains an important but very challenging problem due to the high computational and calculational demands. In the case of the QIM, the dynamics of the system after a quench of one or several parameters has been successfully shown Faribault et al. (2009); Zill et al. (2018). However, in general, QIM without modifications can not be applied to the system with time-dependent parameters. Recently several exact methods for time-dependent Hamiltonians were proposed. In Ref. Sinitsyn et al. (2018), a set of conditions under which the Schrodinger equation can be solved exactly was presented. It was also shown in Ref. Sinitsyn et al. (2018) that among Hamiltonians satisfying these conditions are the multistate Landau-Zener model and the generalized Tavis-Cummings model. Earlier in Ref. Barmettler et al. (2013), Barmettler, Fioretto and Gritsev proposed a generalization of the Bethe wavefunction for the dynamical case and presented its explicit form for the detuning driven Tavis-Cummings model. In Ref. Fioretto et al. (2014), by means of correspondence between the class of Gaudin models and the classical Knizhnik-Zamolodchikov equations some exact solutions for Gaudin-magnets were obtained for special choices of time dependence of the coupling constants. There has also been progress in studies of the exact dynamics of periodically driven systems Gritsev and Polkovnikov (2017). However, to our knowledge a general formulation of how to perform the time evolution of an integrable system has not been shown.

In this paper, we study the generalization of Bethe wavefunction for the time-dependent case. Specifically, consider that we are dealing with an integrable system with Bethe wavefunction

[TABLE]

where $\mathbf{B}(\lambda)$ is an operator which depends on complex parameter $\lambda$ and $|\text{vac}\rangle$ is the pseudo-vacuum reference state, specific to the model being considered. For an initial state that can be represented by the Bethe wavefunction, we show that its time evolution can be described using the dynamical Bethe wavefunction,

[TABLE]

where $p(t)$ is a complex phase. The time dependent wavefunction has exactly the same structure as Bethe wavefunction, but its parameters are functions of time and it has time varying prefactor. One of the most important features of the Bethe vectors is that it allows for the determinant representation for observables Slavnov (1989), which is widely used in calculations of the Bethe ansatz Gamayun et al. (2018); Bulchandani et al. (2018). The fact that the time-dependent wavefunction (2) has the structure of a Bethe vector allows us to transfer all the Bethe ansatz machinery to the time-dependent case.

When a system is exactly solvable by QIM, one can always make (1) an eigenfunction, by choosing special values of the parameters $\lambda_{j}$, which satisfy the Bethe equations. In this paper, we formulate a set of conditions for when the dynamical Bethe wavefunction (2) satisfies the time-dependent Schrodinger equation. The set of conditions is a set of nonlinear coupled differential equations, which we call the dynamical Bethe equations. The time-dependent wavefunction can always be represented by the dynamical Bethe wavefunction (2) for an arbitrary smooth time-dependence of the model parameters if the Hilbert space of the model under consideration is small enough. We provide an explicit example of the dynamical Bethe equations for a detuning driven Bose-Hubbard dimer.

The form of the wavefunction (2) first appeared in Ref. Barmettler et al. (2013) for the Tavis-Cummings model, where the set of dynamical Bethe equations for $\lambda_{j}(t)$ was found, and its connection of trajectories $\lambda_{j}(t)$ with classical motion in a potential was established. So far all the examples of the dynamically integrable models considered in Sinitsyn et al. (2018); Barmettler et al. (2013); Fioretto et al. (2014) belong to the Gaudin class Gaudin (1983) of integrable models or models with a classical R-matrix. Here we show that (2) can be applied to a wider class of integrable models, which goes beyond Gaudin class. The Bose-Hubbard dimer example that we show here belongs to the so-called rational XXX R-matrix class. We note that the set of conditions formulated in Ref. Sinitsyn et al. (2018) does not require that the model should belong to the Gaudin class to be dynamically integrable. Furthermore, models which can be solved by dynamical Bethe wavefunction also do not necessarily satisfy the set of conditions in Ref. Sinitsyn et al. (2018).

The paper is organized as follows. In Sec. II we discuss the general procedure of constructing the dynamical Bethe wavefunction. In Sec. III we derive the dynamical Bethe equations for a Bose-Hubbard dimer with driven detuning and quenching. In Sec. IV we summarize and discuss the future prospects of our method. For more background and details about the Bethe Ansatz technique and our derivations, we refer the reader to the Appendix.

II Dynamical Bethe equations

In this section we discuss the general method of finding the dynamical Bethe wavefunction, without specifying the model. For the reader who is not familiar with Bethe ansatz, we refer them to the Supplemetary Material, or for an extensive review see for example Ref. Slavnov (2018b).

We first assume that the model under consideration can be solved by algebraic Bethe ansatz. We also assume that the set of these Bethe vectors form a complete orthogonal set. This condition should be checked for every specific model separately, but for the vast majority of physically relevant models it is known to be satisfied. Also for simplicity we restrict the considered models to be those with a rational R-matrix and XXX or XXZ-like R-matrices. In practice these three classes cover most physically relevant models.

A central quantity in integrable models is the trace of the monodromy matrix $\tau(\lambda)$ (see Appendix). This operator has many useful algebraic properties provided by the integrability of the model, its specific form should be defined for each model separately. By construction $\tau(\lambda)$ is explicitly connected with the Hamiltonian $\hat{\mathcal{H}}$ of the model under consideration. Usually the Hamiltonian $\hat{\mathcal{H}}$ can be expressed as some elementary function or a residue of $\tau(\lambda)$ at some certain point $\lambda_{0}$. Because of the connection between $\hat{\mathcal{H}}$ and $\tau(\lambda)$, we will see that it is beneficial to consider the following Schrodinger-like equation

[TABLE]

This will allow us to learn the complete information about the time dynamics of the system.

We look for the solution of (3) of the form

[TABLE]

here $M$ is the number of excitations in the system and $\sigma$ enumerates the eigenstates. At $t=0$, the vectors (4) are eigenvectors which form a complete orthogonal set and the set of parameters $\Lambda^{\sigma}_{M}(t)=\{\lambda^{\sigma}_{1}(t),\lambda^{\sigma}_{2}(t),...,\lambda^{\sigma}_{M}(t)\}$ satisfies the stationary Bethe equations for each $\sigma$. We demand time-dependent wavefunctions to also form a complete set

[TABLE]

where we have a proportionality because the wavefunctions are not normalized. The expansion of Bethe vectors (4) over a convenient basis is a difficult problem and in general not solvable, because of the complex structure of (4). For example, $\mathbf{B}(\lambda)$ can be represented as a series of exponential length.

Thus instead of studying the Schrodinger equation (3) directly, we demand that

[TABLE]

for $\sigma<\sigma^{\prime}$, where $\sigma\in[1,L]$ and $L$ is the dimensionality of the Hilbert space under consideration. This states that the Bethe vectors are mutually orthogonal for all $t$. We also demand that

[TABLE]

which must be satisfied for any solution of (3).

We now would like to write (6) and (7) as a set of coupled differential equations. Eq. (6) may re-expressed in this form by writing $|\Psi(\{\lambda^{\sigma}(t)\})\rangle$ in terms of its derivative, which for Bethe vectors always take a special form. To show this, we start by finding the result of operatoring $\tau(\mu)$ on the Bethe wavefunction (4), giving the well-known result

[TABLE]

where $\Theta(\mu,\{\lambda_{j}^{\sigma}\})$ and $\phi_{n}(\mu,\{\lambda_{j}^{\sigma}\})$ are eigenvalues and the off-shell functions defined in (39) and (40) correspondingly. By combining (3), (4), and (8), we obtain

[TABLE]

Demanding that the right hand side is proportional to the left hand side,

[TABLE]

where $f(\{\lambda^{\sigma}_{j}(t)\})$ is a smooth function. If (II) is satisfied we can solve (3) with (4) by choosing special form of phase factor $p^{\sigma}(t)$

[TABLE]

Although it is possible to explicitly find both $f(\{\lambda^{\sigma}_{j}(t)\})$ and $p^{\sigma}(t)$, in practice this is not necessary, because the phase factor $e^{ip^{\sigma}(t)}$ cancels for any observable due to normalization.

After substitution of (II) into (6), the conditions (6) transfers to the set of differential equations:

[TABLE]

Now conditions (7) and (II) are set of $L^{2}-L$ nonlinear differential equations, with $ML$ variables, where $M$ is the number of parameters which parameterize Bethe wavefunction (4). The solution of (7) and (II) is a set of trajectories $\Lambda^{\sigma}(t)=\{\lambda^{\sigma}_{1}(t),...,\lambda^{\sigma}_{M}(t)\}$, for each wavefunction enumerated by $\sigma$. So when $M=L-1$ the number of equations coincides with the number of variables and (6) and (7) always have a solution. So the dynamical Bethe wavefunction can always be constructed if the Hilbert space of the system under consideration is small enough, for arbitrary smooth time dependence of the parameters of the model.

In Bethe ansatz it is typical for the Bethe wavefunction to be parameterized by a number of parameters which is much smaller than size of the Hilbert space. For example, while the Hilbert space of the XXZ Heisenberg magnet has an exponentially large dimension, its Bethe wavefunction is parametrized by a number of parameters linearly proportional to the number of excitations, which provides a great advantage in terms of computational complexity. The equations (7) and (II), however, become overdetermined if the dimensionality of Hilbert space $L>M+1$. Nevertheless, in principle, the existence of solutions for (7) and (II) when it is overdetermined is not prohibited because the equations are nonlinear. A trivial example of such a solution is adiabatic evolution when $\{\lambda^{\sigma}_{1}(t),...,\lambda^{\sigma}_{M}(t)\}$ is a solution of static Bethe equations at every moment.

III Example: Bose-Hubbard dimer

We now illustrate the above method to obtain the time dynamics by applying it to the Bose-Hubbard dimer. This model provides both a simple and non-trivial example of a dynamically integrable model from the XXX class. For this particular case the dynamical Bethe equations can be written in a particularly simple and explicit form.

The Bose-Hubbard dimer in the two mode approximation (equivalent to a two-site Bose-Hubbard model) can be described by the following Hamiltonian Milburn et al. (1997); Ermakov et al. (2018)

[TABLE]

where $a,b$ are bosonic operators for the two sites satisfying $[a,a^{\dagger}]=[b,b^{\dagger}]=1$. The total number operator of particles $\hat{N}=a^{\dagger}a+b^{\dagger}b$ is a conserved quantity, $[\mathcal{\hat{H}},\hat{N}]=0$. For the Bethe ansatz formalism it is convenient rescale and offset the Hamiltonian by defining

[TABLE]

which commutes with (13). Defining the dimensionless detuning $\Delta=\frac{2\epsilon}{J}$ and coupling constant $c^{2}=\frac{U-V}{J}$, the Hamiltonian can be rewritten

[TABLE]

which is the form we shall use. Using the dimensionless Hamiltonian (15) means that all energies are measured in units of $J$ and time is measured in units of $\hbar/J$.

We now assume that the coupling constant $c$ is time-independent whereas the detuning $\Delta(t)$ continuously depends on time. We introduce the generalized creation operator $\mathbf{B}(\lambda)$ which depends on a complex parameter $\lambda$ Bogoliubov (2016)

[TABLE]

where

[TABLE]

The pseudo-vacuum in this case is simply the zero particle Fock state

[TABLE]

The time-dependent Bethe wavefunction can then be written following (4) using the above definitions. We first consider the case where only the detuning is time-dependent. In this case the dynamical Bethe equations can be written in a compact form given by

[TABLE]

Here $\varphi_{n}(\{\lambda^{\sigma}\})$ are the so-called off-shell functions defined as

[TABLE]

When $\varphi_{n}(\{\mathbf{\lambda}\})=0,\;\forall n=1,...,N$, where $N$ is the number of particles, these reduce to the static Bethe equations.

The dynamical Bethe equations are a set of first order coupled ordinary differential equations. As the initial condition for (19) we need to pick a set of parameters $\{\lambda(0)\}=\{\lambda_{1}(0),...,\lambda_{N}(0)\}$, which parametrizes the initial state $|\Psi_{N}(0)\rangle$. For example, if the initial state is an eigenstate, the set $\{\lambda(0)\}$ should satisfy the static Bethe equations.

We numerically solve the set of equations (19) for a detuning with time dependence

[TABLE]

which has a rather non-trivial non-linear and aperiodic dependence. The initial condition was chosen to be the solution of static Bethe equations which corresponds to the ground state of (15). As the observable, we calculate the intersite coherence

[TABLE]

for the details of calculation of $\nu(t)$ see Appendix B. Figure 1(a) shows our results. We find that the method perfectly reproduces time dynamics calculated by exact diagonalization, giving identical curves. The method is computationally efficient the solution requires the evolution of $N$ coupled equations. In Fig. 1(c) we show the stroboscopic maps of the solutions of dynamical Bethe equations (19), point of certain color corresponds to the value of the component of the solution of (19) $\lambda_{j}(t_{k})$ at the moment $t_{k}$. Instead of solving (19) one may solve more general system of equations (7),(II), which is applicable for arbitrary time dependece of both $\Delta$ and $c^{2}$, we checked that solutions of (7),(II) does perfectly coincide with the solutions of (19) when only detuning is driven.

As a second example, we consider the case of a quench, when the parameters are changed suddenly from $c,\Delta$ to $c^{\prime},\Delta^{\prime}$. When all the parameters of the model are constant the set of equations (19) become

[TABLE]

The set of equations (23) describes the evolution of an initial state $|\Psi_{N}(0)\rangle$ with a static Hamiltonian (15). The initial state can be parameterized by a Bethe vector with the set of parameter $\{\lambda^{0}\}$ satisfying (55) for the initial parameters $c,\Delta$. After the quench is performed, the Hamiltonian parameters change to $c^{\prime},\Delta^{\prime}$, hence we need to establish the connection between the old wavefunction expressed in terms of $c,\Delta$, and the new one expressed in terms of $c^{\prime},\Delta^{\prime}$. The initial condition for (23) thus is given by

[TABLE]

For the case that only the detuning is quenched $c^{\prime}=c$, the initial conditions for (23) can be simply found to be

[TABLE]

In Fig. 1(b) we plot an example solution of the intersite coherence (22) from the dynamical Bethe equations (23). We again see that there is perfect agreement of the time dynamics with numerical results obtained from exact diagonalization. In Fig. 1(d) we plot stroboscopic maps for the solution of (23) in the same fashion as we did for (19).

IV Outlook and conclusions

We have described a method for evaluating the time dynamics of systems that are exactly solvable by Bethe ansatz. The method is based on the dynamical Bethe wavefunction (2) which is a straightforward generalization of Bethe ansatz for dynamical case, where the Bethe parameters are time dependent and there is a time varying complex phase. The main advantage of the dynamical Bethe wavefunctions (2) is that they are mathematically manageable thanks to the well-developed Bethe ansatz results which are directly applicable.

The set of differential dynamical Bethe equations (7) and (II) can be applied to any Bethe ansatz solvable model from XXX, XXZ or Gaudin class which has a dimensionality not bigger than $N+1$, where $N$ is the number of parameters in the Bethe wavefunctions. What would be interesting is if the dynamical Bethe wavefunction could describe the diabatic evolution of a non-trivial model with a larger Hilbert space than $N+1$. This would be an example of exact non-ergodic behavior which is a topic of great importance Turner et al. (2018). We have shown that our approach produces exact time dynamics with the Tavis-Cummings model, which possesses formally a Hilbert space of dimension $2^{N}$, but can be restricted by symmetry to a dimension $N+1$. Thus this alone does not demonstrate a completely non-trivial example. Currently, the in terms of computational advantage, the dynamical Bethe equations only provide an equivalent approach to alternative techniques, since both scale as $N$. However, we do not exclude the possibility that there are systems with larger Hilbert space for which overdetermined system (7) and (II) might have a non-trivial solution.

Acknowledgements.

The authors are grateful to O. Lychkovskiy for useful discussions. I. E. is supported by the Russian Science Foundation under the grant No 17-71-20158.

Appendix A Algebraic Bethe ansatz

Here we briefly sketch the main aspects of the algebraic Bethe ansatz technique which are necessary for the understanding of the present paper. For an extensive review of the Bethe ansatz, we refer the reader to Refs. Slavnov (2018b); Korepin et al. (1997); Essler et al. (2005).

The cornerstone of any integrable model is the R-matrix which in this paper always takes the form

[TABLE]

Here the entries $f(\mu,\lambda)$ and $g(\mu,\lambda)$ are specified for each model separately. In general, the R-matrix is the solution of the Yang-Baxter equation Slavnov (2018b), and can take many different forms. The specific form of the R-matrix generates a family of integrable models.

In order to construct an integrable model we need to define the monodromy matrix

[TABLE]

which depends on the complex spectral parameter $\lambda$. Here, $A(\lambda),B(\lambda),C(\lambda)$ and $D(\lambda)$ are operators acting in the Hilbert space of the model under consideration, and their explicit representation depends on the model. The monodromy matrix should also satisfy the Yang-Baxter equation

[TABLE]

To construct the Hamiltonian of a particular integrable model, we define the trace of the monodromy matrix

[TABLE]

The Hamiltonian $\hat{\mathcal{H}}$ of the model may be expressed via the trace of the monodromy matrix $\tau(\lambda)$, or its derivative at some specified $\lambda=\lambda_{0}$. Usually it can be expressed as some elementary function of $\tau(\lambda_{0})$, or as a residue of the $\tau(\lambda)$ at a particular point $\lambda_{0}$.

The pseudovacuum state $|\text{vac}\rangle$ is a state from the Hilbert space of the model, which is annihilated by the operator $C(\lambda)|\text{vac}\rangle=0$. The conjugated operator also satisfies $\langle\text{vac}|B(\lambda)=0$. Usually the pseudovacuum state is an eigenstate of the system, but, in general, it is not required. We also define two eigenvalue functions $a(\lambda)$ and $d(\lambda)$ according to

[TABLE]

The Bethe wavefunction is then defined as

[TABLE]

where $\{\lambda^{\sigma}_{j}\}$ is the set of complex parameters $\{\lambda^{\sigma}_{j}\}=\{\lambda^{\sigma}_{1},\lambda^{\sigma}_{2},...,\lambda^{\sigma}_{M}\}$, and $M$ is the number of excitations in the system, and $\sigma$ labels the wavefunction. The wavefunction (31) is an eigenfunction of the trace of monodromy matrix $\tau(\lambda)$

[TABLE]

if the set $\{\lambda^{\sigma}_{j}\}$ satisfies to the set of Bethe equations

[TABLE]

All the roots within one solution $\{\lambda^{\sigma}_{j}\}$ should be different, otherwise $|\Psi(\{\lambda^{\sigma}_{j}\})\rangle$ can not be an eigenfunction. The Bethe equations (33) are set of coupled nonlinear algebraic equations. It has $M$ equations and $N$ solutions, where $N$ is the size of Hilbert space.

For our purposes it is very important to know what the effect of the transfer matrix $\tau(\lambda)$ acting on the Bethe vector (31). For notational simplicity we omit the index $\sigma$ henceforth, such that $\{\lambda_{j}\}$ denotes the set $\{\lambda_{1}^{\sigma},\lambda_{2}^{\sigma},...,\lambda_{M}^{\sigma}\}$. From Ref. Slavnov (2018b) it is known that

[TABLE]

Here we defined the functions

[TABLE]

For the $D(\lambda)$ operator, we have the similar expressions

[TABLE]

where we defined

[TABLE]

Combining these results we can find the effect of acting $\tau(\lambda)$ on the Bethe wavefunction, given by

[TABLE]

Here we defined

[TABLE]

and the off-shell function as

[TABLE]

If we now demand that the off-shell function (40) is zero, it is evident that the wavefunction (31) is an eigenfunction for $\tau(\lambda)$. The roots of off-shell functions (40) coincide with the roots of Bethe equations (33), but we should distinguish between these since later we will encounter cases where the off-shell function is not zero.

Finally, we mention several important properties of Bethe wavefunctions. The dual Bethe wavefunctions are defined as

[TABLE]

In general, despite the notation, the wavefunction (41) does not coincide with the hermitian conjugate of the function (31), i.e. $\langle\Psi(\{\lambda^{\sigma}_{j}\})|\neq|\Psi(\{\lambda^{\sigma}_{j}\})\rangle^{\dagger}$. Dual vectors like this must be introduced in order to evaluate scalar products and averages of observables. Generally, in the literature devoted to Bethe Ansatz, the left bracket $\langle\Psi|$ implies the dual vector (41).

For most of the integrable models it has been proven that Bethe vectors form a complete set Korepin et al. (1997); Slavnov (2018b)

[TABLE]

where $\hat{I}$ is the identity operator, and $N$ is the size of the Hilbert space. In general Bethe wavefunctions are not normalized.

One of the most important properties of Bethe wavefunctions is that for many models it is possible to evaluate the scalar product of Bethe wavefunctions and averages of the operators by applying Slavnov’s formula Slavnov (1989). This allows one to express scalar product as a determinant. We do not reproduce the general form of the Slavnov’s formula here because of its complexity, and it not very useful to consider it without specifying the model. Application of Slavnov’s formula to the models considered in this paper have been studied in Refs. Bogoliubov (2016); Bogoliubov et al. (2017).

Appendix B Dynamical Bethe equations for the Bose-Hubbard dimer

Here we give more details of the derivation of the dynamical Bethe equations for the detuning driven Bose-Hubbard dimer. A more detailed description regarding the Bethe ansatz solution of this model can be found in Ref. Bogoliubov (2016), we use the same notations as this paper.

The Hamiltonian of the Bose-Hubbard dimer is

[TABLE]

The diagonal elements of the monodromy matrix are in this case

[TABLE]

The Hamiltonian (43) can then be expressed via trace of the monodromy matrix (29) according to

[TABLE]

According to the definitions (A), the eigenvalue functions are then

[TABLE]

The elements of the R-matrix are defined as

[TABLE]

We now wish to look for Bethe eigenfunctions of the form

[TABLE]

where the pseudo-vacuum state is $|\text{vac}\rangle=|0\rangle_{a}\otimes|0\rangle_{b}$. $\sigma$ is the index which labels the energy levels of the system, for the sake of notational simplicity we omit this below. The eigenvector depends on $N$ complex parameters $\{\lambda\}=\{\lambda_{1},\lambda_{2},...,\lambda_{N}\}$. By applying the Bethe ansatz machinery we can evaluate

[TABLE]

where we have defined

[TABLE]

Here $E_{N}(\{\lambda\})$ is the energy and $\varphi_{n}(\{\mathbf{\lambda}\})$ is the off-shell function. From (B) we can see that when set $\{\lambda\}$ satisfies

[TABLE]

the wavefunction (51) becomes an eigenfunction of the Hamiltonian (43). The set of equations (55) are known as the Bethe equations.

We now look for a time-dependent wavefunction of the form

[TABLE]

If only the detuning $\Delta$ is time-dependent, it is easy to see that $[\frac{d}{dt}\mathbf{B},\mathbf{B}]=0$, and the derivative of (56) can be taken easily. Substituting (56) into the time-dependent Schrodinger equation we obtain

[TABLE]

If we demand now that

[TABLE]

the wavefunction (56) will satisfy the time-dependent Schrodinger equation. We call the set of conditions (19) the dynamical Bethe equations. The dynamical Bethe equations are set of first order coupled ordinary differential equations. For the initial condition of (19) we need to pick a set $\{\lambda(0)\}=\{\lambda_{1}(0),...,\lambda_{N}(0)\}$, which parametrizes the initial state $|\Psi_{N}(0)\rangle$. For example if the initial state is an eigenstate, the set $\{\lambda(0)\}$ should satisfy the static Bethe equations (55). The phase factor $p(t)$ is given by

[TABLE]

To evaluate observables one may use the determinant representation as a general approach Slavnov (1989); Bogoliubov (2016). More convenient approach is to use the expansion of Bethe vectors (51) over the Fock space, which was developed in Ermakov et al. (2018):

[TABLE]

where the coefficient $\Gamma_{lmk}$ defined as

[TABLE]

and $D(M,k)$ are coefficients defined by the following recurrence relation

[TABLE]

with the conditions: $D(1,1)=1$ and $D(M,k)=0$ if $k>M$. This coefficient possess the obvious property: $D(M,1)=D(n,n)=1$. The general expression for $D(M,k)$ is given by

[TABLE]