Bifurcation and chaos in nonlinear Lindblad equations

TL;DR

This paper investigates nonlinear Lindblad equations in open quantum systems, revealing that they can exhibit bifurcations and chaos, which are relevant for certain experimental scenarios involving mean field theories or classical control.

Contribution

It introduces the analysis of nonlinear Lindblad equations using classical dynamical systems techniques, highlighting the emergence of bifurcations and chaos in quantum dynamics.

Findings

Nonlinear Lindblad equations can exhibit bifurcations.

Chaotic dynamics are possible in these systems.

Relevance to experimental quantum systems is discussed.

Abstract

The Lindblad equation describes the dissipative time evolution of a density matrix that characterizes an open quantum system in contact with its environment. The widespread ensemble interpretation of a density matrix requires its time evolution to be linear. However, when the density matrix is obtained from a mean field theory of interacting quantum systems or from a top-down control by a changing classical environment, the ensemble interpretation is inappropriate and nonlinear dynamics arise naturally. We therefore study the dynamical behavior of nonlinear Lindblad equations using the example of a two-level system. By using techniques developed for classical dynamical systems we show that various types of bifurcations and even chaotic dynamics can occur. We also discuss experimental situations for which our results could be relevant.

Equations

$${\partial }_t\mathbf{\rho }=\underset{\text{von Neumann term}}{\underbrace{-\frac{i}{\hslash }[H,\mathbf{\rho }]}}+\underset{\text{dissipator}}{\underbrace{\sum _k{\gamma }_k\left(A_k\mathbf{\rho }{A}_k^{†}-\frac{1}{2}\{A_k^{†}{A}_k,\mathbf{\rho }\}\right)}},$$

$$L_1:=\left(\begin{array}{cc}0& 0\\ 1& 0\end{array}\right),L_2:=\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right),L_3:=\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right).$$

$$\dot{\mathbf{\rho }}=-\frac{i}{\hslash }[H,\mathbf{\rho }]+\sum _{i,j=1}^3{h}_{ij}\left(L_i\mathbf{\rho }{L}_j^{†}-\frac{1}{2}\{L_j^{†}{L}_i,\mathbf{\rho }\}\right),$$

$$\mathbf{\rho }=\frac{1}{2}\left(\text{mathds}1+\mathbf{x}\cdot \mathbf{\sigma }\right)=\frac{1}{2}\left(\begin{array}{cc}1+z& x-iy\\ x+iy& 1-z\end{array}\right),$$

$$U:=\left\{(x,y,z)\in {\mathbb{R}}^3,|x^2+y^2+z^2|\le 1\right\}.$$

$$\dot{z}$$

$$\dot{x}$$

$$\dot{y}$$

$$\dot{s}=-\Gamma s$$

$$\dot{z}=\left(h_{11}-h_{22}\right)-\left(h_{11}+h_{22}\right)z,$$

$$\underset{t\to \infty }{\lim }\mathbf{\rho }(t)=\frac{1}{h_{11}+h_{22}}\left(\begin{array}{cc}{h}_{11}& 0\\ 0& h_{22}\end{array}\right).$$

$$\dot{z}=-z\cdot \left(t+z^2\right),$$

$$h_{11}(z)=$$

$$h_{22}(z)=$$

$$H=-\frac{1}{2}\sum _{ij}{J}_{ij}{\mathbf{\sigma }}_i\cdot {\mathbf{\sigma }}_j=-\frac{1}{2}\sum _i{\mathbf{\sigma }}_i\cdot \sum _j{J}_{ij}{\mathbf{\sigma }}_j\approx -\frac{1}{2}\sum _i{\mathbf{\sigma }}_i\cdot {⟨\sum _j{J}_{ij}{\mathbf{\sigma }}_j⟩}_i=-\frac{1}{2}\sum _i\alpha \mathbf{m}\cdot {\mathbf{\sigma }}_i\equiv \frac{1}{2}\sum _i{\hat{h}}_i,$$

$$h_{11}={\gamma }_1=\gamma \frac{e^{-\beta \alpha z}}{e^{-\beta \alpha z}+e^{\beta \alpha z}}\text{ and }{h}_{22}={\gamma }_2=\gamma \frac{e^{\beta \alpha z}}{e^{-\beta \alpha z}+e^{\beta \alpha z}}$$

$$\dot{z}=\gamma \tanh (\beta \alpha z)-\gamma z,$$

$$\dot{z}=-z\cdot \left(t+z^2\right)+b.$$

$$h_{11}(z)=$$

$$h_{22}(z)=$$

$$\alpha$$

$$|t-2\alpha |$$

$$|b|$$

$$\hat{h}=-\alpha \mathbf{m}\cdot \mathbf{\sigma }-\mu \mathbf{B}\cdot \mathbf{\sigma }.$$

$$\dot{s}=-\Gamma s-i\frac{2\mu }{\hslash }Bs,$$

$$\dot{z}=\gamma \tanh (\beta (\alpha z+\mu B))-\gamma z,$$

$$\dot{x}$$

$$\dot{y}$$

$$\dot{x}=$$

$$\dot{y}=$$

$$\dot{r}=$$

$$\dot{\varphi }=$$

$$\text{Re}[h_{23}]=-\text{Re}[h_{13}]=\frac{(ϵ+a-\kappa )x-by}{4}$$

$$\text{Im}[h_{23}]=\text{Im}[h_{13}]=\frac{(ϵ+a+\kappa )y+bx}{4}$$

$$h_{11}\ge$$

Full text

K****eywords Lindblad equation  $\cdot$ quantum dynamics  $\cdot$ nonlinear dynamics  $\cdot$ bifurcations

1 Introduction

The description of a quantum system by a Schrödinger equation is based on the assumption that the system is sufficiently isolated from the rest of the world that no uncontrollable influences affect its time evolution [1]. In such a situation, the external environment needs to be taken into account only via boundary conditions, electromagnetic potentials or other potential energy terms. However, when environmental influences that cannot be controlled in detail affect the quantum system, the Schrödinger equation is not appropriate any more [2, 3]. As is familiar from statistical mechanics, an environment can induce random transitions between the states of a quantum system so that it does not undergo a unitary time evolution any more. A frequently used equation that captures this effect of the environment is the Lindblad equation [4]. It is a Master Equation for the density matrix $\boldsymbol{\rho}$ of the quantum system that contains transitions between states in addition to the von Neumann term that describes unitary time evolution. There are two general ways to derive the Lindblad equation: the first one starts from a quantum description of the system and the bath and takes the trace over the states of the bath in order to obtain a time evolution equation for the reduced density matrix that describes only the system [2]. This derivation requires a couple of ad-hoc assumptions that cannot really be justified, in particular that the combined state of system and bath can be written as a product state. The second way to derive the Lindblad equation consists in writing down the most general equation that satisfies the requirement that the density matrix remains a density matrix, i.e., that it remains Hermitian, positive semi-definite, and with trace 1. This means that the time evolution must be completely positive and trace preserving [5]. Furthermore, in order to preserve the ensemble interpretation of the density matrix, the time evolution is required to be linear in $\boldsymbol{\rho}$. It can be proven mathematically that the most general equation that satisfies all these criteria is the Lindblad equation. This second derivation makes no assumption about the nature of the environment, except that its influence on the quantum system depends only on the present state of the quantum system and not on its past, as no memory terms are included.

Due to the dissipative terms unravelings of Lindblad Equations are non-linear in the wave function and hence they violate superposition principle. The Lindblad equation is only one of several examples where the theoretical description of a quantum system does not follow unitary time evolution. In particular in condensed matter theory, methods that violate unitary time evolution are widely employed [6]. One well-known example are the Hartree and Hartree-Fock theory for many-fermion systems, where the influence of the other particles on a given particle is taken into account via the potential generated by the other particles, which in turn is proportional to their charge density [7]. These theories are nonlinear in the wave functions. Applying them to bosons gives the Gross-Pitaevskii equation used for Bose condensates. All these equations are mean-field equations that replace the explicit interaction terms that should occur in a many-particle Schrödinger equation by a classical interaction via a potential.

When such a nonlinear version of the Schrödinger equation is translated into the corresponding von Neumann equation, the quantum-state dependent terms of the Hamilton operator $H$ cause the commutator $[H,\rho]$ to become nonlinear in the density matrix $\rho$ [8]. This shows that widely used quantum theories can involve time evolutions that are nonlinear not only in the wave function, but also in the density matrix. In fact, Breuer and Petruccione (see Section 3.7 in [8]) list two additional classes of Master equations that are nonlinear in the density matrix: (i) nonlinear Boltzmann equations, which describe the time evolution of a one-particle density matrix due to collisions with other particles, (ii) mean-field Master Equations for $n$ identical interacting quantum systems, where the average influence of the other systems on the time evolution of the density matrix of one system is included in a similar way as in the Hartree method.

A discrete version of nonlinear dynamics of the density matrix is obtained via carefully designed quantum state transformations, as for instance described by Bechmann et al. [9]. In such transformations, several spins are prepared in specified states and are then subjected to an interaction that causes them to become entangled (for example via a C-NOT gate). When then a measurement with an additional post-selection is performed, the resulting quantum state transformation is nonlinear. By applying such a quantum state transformation iteratively, one can obtain complex dynamics and even chaos [10, 11, 12]. By chaos, the authors mean exponential sensitivity of the dynamics of the density matrix with respect to the initial condition.

By comparing this quantum state transformation system that shows a discrete dynamics of the density matrix with the above-mentioned three types of systems that show a continuous nonlinear dynamics of the density matrix, we can identify two shared features: The first feature is an interaction between identical particles that are part of the quantum system of interest. Since the density matrix for all these models is that of one quantum particle, this interaction makes the interpretation of the density matrix as an ensemble of independent systems invalid. The second feature is the interaction with a classical external world that destroys unitary time evolution. In the simplest cases, the role of this external world consists merely in imposing a temperature on the system and causing it to relax to thermal equilibrium (for instance for systems described by a Boltzmann equation). In the quantum information example, this external world consists in a device for state preparation and measurement/postselection and thus establishes a nonequilibrium situation, which in turn is a prerequisite for obtaining dynamical trajectories that do not relax towards fixed points.

Another way of establishing a nonequilibrium situation consists in periodic driving. Indeed, a periodically driven quantum many-particle system was found to show periodic attractors, bifurcations, and chaos when a mean-field calculation of the dynamics was performed [13].

This widespread occurrence of time evolutions that are nonlinear in the density matrix suggests that bifurcations and nontrivial attractors can occur generically in open quantum systems if they are driven away from equilibrium, although the quantum state transformation example suggests that this requires a careful design of the system. In this paper, we therefore want to explore in a more general way the occurrence of various types of bifurcations, and in particular of limit cycles and chaos within the framework of the Lindblad equation. We will focus on the simplest possible quantum system, namely a two-level system, or, equivalently, a spin-1/2. We will present two different types of approaches: The first one will be a generic demonstration that various well-known types of bifurcations can be obtained by making the transition rates in the Lindblad equation dependent on the density matrix and ensuring that the density matrix remains positive semi-definite with trace one under the resulting time evolution. This includes the occurrence of limit cycles and chaos. The second one will be a more microscopic approach, where a series of explicit physical model will be constructed in order to illustrate that one can indeed conceive of real-world systems that show these nonlinear dynamics phenomena. The ingredients of these explicit models are the following: First, we include an interaction between the individual quantum systems that constitute the ensemble, having mostly the spins of a solid from condensed matter physics in mind. The interaction between these spins, which will be taken into account on the level of a mean-field approximation, leads to nonlinearities, making the energy eigenstates of the spins dependent on the density matrix. By applying a finite temperature, the dynamics of the Lindblad equation consists in spin flips between and projections onto energy eigenstates of the spins. In order to create an ongoing nonequilibrium situation, we couple the spins to an environment that measures the magnetization and responds by applying a magnetic field the strength and direction of which depend on the measured magnetization, which in turn is described by the density matrix.

2 The Lindblad equation

The Lindblad equation describes the time evolution of the density matrix $\boldsymbol{\rho}=\boldsymbol{\rho}(t)$ of an open quantum system in the limit where no memory terms are required, [8]

[TABLE]

with $H$ being the Hamilton operator, $[\cdot,\cdot]$ and $\{\cdot,\cdot\}$ denoting the commutator and anticommutator respectively, and the $A_{k}$ being the Lindblad operators that introduce dissipation into the quantum system. They can be chosen such that they form an orthonormal basis, that is $\text{trace}(A_{k}^{\dagger}A_{j})=\delta_{kj}$. Below, we will only consider examples where the $A_{k}$ are trace-free. This is no real restriction as a model with nonvanishing traces of the Lindblad operator can be mapped onto one with vanishing traces by modifying the Hamilton operator [8]. We will consider a two-state model. This can represent any two-level system, but the examples discussed below are best suited for a spin-1/2 system. In the following, we will use the non-diagonal form of the Lindblad equation that is obtained by expressing the Lindblad operators in a given basis of trace-free matrices,

[TABLE]

By writing $A_{k}=\sum\limits_{i=1}^{3}a_{i}^{(k)}L_{i}$, the Lindblad equation becomes

[TABLE]

with $h_{ij}=\gamma_{k}\,a_{i}^{(k)}a_{j}^{(k)^{*}}=h_{ji}^{*}$ being positive semi-definite.

In the standard Lindblad equation, the transition rates $\gamma_{k}$ and the Lindblad operators $A_{k}$ are fixed model ingredients that capture the effect of the environment on the quantum system. We will, however, consider the possibility that the transition rates and Lindblad operators depend on the density matrix $\boldsymbol{\rho}$, as explained briefly in the Introduction, and as made explicit by the specific models given further below. In the representation (3), the dependence of the $A_{k}$ on the density matrix goes into the coefficients $a_{i}^{(k)}$, as does the dependence of the $\gamma_{k}$ on the density matrix.

It is useful to write $\boldsymbol{\rho}$ in terms of three parameters,

[TABLE]

where $\boldsymbol{\sigma}$ is the vector of Pauli spin matrices, and $x$, $y$ and $z$ are real-valued. The only constraint on their values is that $(x,y,z)$ must stay inside the admissible region $U$, which assures that $\boldsymbol{\rho}$ stays positive semi-definite,

[TABLE]

By defining $\Gamma:=\frac{1}{2}\left(h_{11}+h_{22}+4\,h_{33}\right)$, we can write equation (3) in the following form,

[TABLE]

The last terms are due to the von Neumann term, with $H_{ij}$ being the matrix elements of the Hamiltonian with respect to the two basis states $\left|0\right>$ and $\left|1\right>$. As we assume that the transition rates and the Hamilton operator can depend on the density matrix, we have in general $h_{ij}=h_{ij}(x,y,z)$ and $H=H(x,y,z)$. We can ensure that the matrix $h_{ij}$ remains positive semi-definite at all times either by calculating explicitly the relevant constraints on the parameters of the model, or by using the diagonal form (1) with positive transition rates $\gamma$.

Below, we will consider only special cases where part of the terms vanish or have a simpler form. We will start with the case that the system of equations (6) is one-dimensional, where we can already investigate pitchfork bifurcations and saddle-node bifurcations. When we go to the two-dimensional case, we will additionally obtain Hopf bifurcations and limit cycles, and the three-dimensional case will furthermore yield strange attractors.

3 One-dimensional case

The system of equations (6) becomes one-dimensional when the transition matrices $A_{k}$ make transitions between and measurements of the two eigenstates of the Hamilton operator. If we choose as basis states the two eigenstates of the Hamilton operator, the matrices $A_{k}$ become identical to the operators $L_{1}$ to $L_{3}$, and the coefficient matrix $h$ becomes diagonal. Furthermore, a diagonal density matrix $\boldsymbol{\rho}$ remains diagonal under time evolution, and the von Neumann term vanishes. In this situation, equations (6) reduce to the equation

[TABLE]

for the off-diagonal elements $s=x+i\,y$ and

[TABLE]

which determines the time evolution of the diagonal elements. As the off-diagonal entries of the density matrix decrease to zero, it is sufficient to investigate the time evolution of $z$.

When the transition rates are independent of $\boldsymbol{\rho}$, the time evolution goes to

[TABLE]

If the transition rates depend on $z$, several fixed points can occur, as illustrated in Figure 1(b), and a parameter change can lead to various bifurcations as fixed points are created or destroyed.

In the following, we will perform a generic study of the pitchfork bifurcation and the saddle-node bifurcation, and we will build specific physical systems that will display these bifurcations.

3.1 Pitchfork bifurcation

When the system has an intrinsic symmetry around $z=0$, equation (8) contains only odd powers when written as a polynomial of $z$. Such a symmetry occurs for instance for a two-level spin system when none of the two spin orientations is preferred by the environment. The minimum model showing this pitchfork bifurcation is given by a polynomial of the order 3. Through an appropriate choice of the time scale, this polynomial can be brought to the form

[TABLE]

which is the normal form of the pitchfork bifurcation and depends only on the parameter $t$, with a stable fixed point at $z=0$ for $t>0$ and two stable fixed points at $z=\pm\sqrt{-t}$ for $t<0$. An explicit form for the transition rates is in this case

[TABLE]

This ansatz respects the symmetry between the states $\left|0\right>$ and $\left|1\right>$ and their transition rates, namely $h_{11}(z)=h_{22}(-z)$. A sufficient criterion for the transition rates $h_{ii}(z)$ being non-negative is that the parameters $t$ and $\alpha$ satisfy $\alpha\in(0,2)$, $|t-2\,\alpha|\leq\sqrt{8\,\alpha}$. The normal form of the pitchfork bifurcation (expression (10)) can only be a good approximation in a sufficiently small region around this bifurcation. We will see this in the following example.

An explicit, microscopic model is given by an ensemble of spins that are associated with a magnetic moment, which are embedded in a medium, for instance a crystal. They are coupled to each other and experience flips due to the interaction with the phonon heat bath of temperature $T$ of the medium. We interpret now the density matrix $\boldsymbol{\rho}{\,=\,\frac{1}{2}\left(\mathds{1}+\boldsymbol{r}\cdot\boldsymbol{\sigma}\right)}$ as describing the ensemble of these spins, assuming that their number is so large that stochastic fluctuations of the state of this ensemble can be neglected. The vector ${\boldsymbol{r}\,=\,}(x,y,z)$ is in this case proportional to the magnetization $\boldsymbol{m}$ of the spin system, and simultaneously it is the averaged magnetic moment of a single spin that flips between its two orientations. The Hamilton operator reads

[TABLE]

where we have made the mean-field approximation that each spin sees the same local field $\alpha\,\boldsymbol{m}$ caused by the other spins. We choose the normalization of $\boldsymbol{m}$ such that its components are identical to the entries of the density matrix representing the ensemble of spins, i.e., $\boldsymbol{m}=(x,y,z)$.

An effectively one-dimensional model is obtained if the embedding medium imposes an easy axis on the system, as is the case in uniaxial magnetic systems. Orientation of the spins along this axis is energetically favored in such systems. We implement in the following the easy axis by the choice of Lindblad operators. By choosing $L_{1}$ to $L_{3}$, see (2), the easy axis is the $z$ axis. We have already shown above that in this case the off-diagonal elements of the density matrix relax to zero (see (7)), and that the dynamics of the diagonal elements, i.e., the dynamics of $z$, decouples from the other two equations and is described by (8). We consider therefore in the following only the $z$ dynamics.

By setting $\boldsymbol{m}=(0,0,z)$, the two eigenvalues of the single-spin hamiltonian $\hat{h}$ are $\pm\,\alpha z$. The interaction with the heat bath causes the spins to flip between these two energy eigenstates, which means that the Lindblad equation contains the two Lindblad operators $L_{1}$ and $L_{2}$. The operator $L_{1}$ makes a spin flip from the spin orientation in positive $z$ direction to that in negative $z$ direction, while $L_{2}$ makes the reverse flip. The flip rates must be chosen such that we obtain Boltzmann weights in thermal equilibrium, leading to

[TABLE]

with $\beta=1/k_{B}T$. Since the spins flip between their energy eigenstates, the von Neumann term, which is obtained by using the single-spin hamiltonian (see equation (12)) $\hat{h}=-\alpha\boldsymbol{\sigma}\cdot\boldsymbol{r}=-\alpha z\sigma_{z}$, vanishes.

The dynamical equation (8) for $z$ then becomes

[TABLE]

For $\alpha\,\beta\leq 1$ , equation (14) has only one fixed point at $z=0$, which is stable. For $\alpha\,\beta>1$, the fixed point at $z=0$ is unstable, and there are two additional fixed points $\pm\,z^{*}$ that are given by $\tanh(\alpha\,\beta\,z^{*})=z^{*}$. These fixed points are stable. We thus have a pitchfork bifurcation at $\beta\alpha=1$.

Not surprisingly, our calculation has resulted in a phase transition to ferromagnetism, as it has used the ideas behind the mean-field theory by Weiss. When $\beta\alpha$ is close to 1, $z$ is small, we have $\tanh(x)=x-x^{3}/3+\mathcal{O}(x^{5})$, and the normal form (10) of the bifurcation is a good approximation, but for larger values of $z$ higher-order terms of the Taylor expansion become important.

3.2 Saddle-node bifurcation

When the system contains no symmetry around $z=0$, we can add a constant term to the right-hand side of equation (10), resulting in

[TABLE]

With the choice $t>0$, this system has 3 fixed points for $b=0$. When the parameter $b$ crosses a critical value $b_{c}=\pm 2\left(\frac{|t|}{3}\right)^{\frac{3}{2}}b_{c}$, a saddle-node bifurcation happens as the unstable fixed point collides with one of the stable fixed points, with the outcome that both of them vanish.

The bifurcation equation (15) can be expressed in the form of equation (8) with non-negative functions $h_{ii}(q)$ for instance by setting

[TABLE]

The parameters $\alpha$, $t$ and $b$ must satisfy

[TABLE]

Figure 3 shows the stability diagram of this system.

An explicit microscopic model that shows such a saddle-node bifurcation is obtained when model (12) is supplemented by an external magnetic field $\boldsymbol{B}$, resulting in the single-spin hamiltonian

[TABLE]

We choose $\boldsymbol{B}=(0,0,B)$ so that it points along the easy axis in positive $z$ direction. As before, we implement the easy axis in the simplest possible way by choosing $L_{1}$ and $L_{2}$, see (2) as Lindblad operators. Due to the von Neumann term, see equation (6), (7) is now extended to

[TABLE]

which relaxes again to $s=0$. The $z$ dynamics, given by (8), now takes the form

[TABLE]

which shows saddle-node bifurcations as $B$ or $T$ are changed. For parameter values such that $z$ is small at the fixed points, we can again set $\tanh(x)=x-x^{3}/3+\mathcal{O}(x^{5})$ and regain the form (15).

When the medium in which the spins are embedded is isotropic, there is no easy axis with respect to which we can choose the Lindblad operators. Nevertheless, the magnetic field field $\boldsymbol{B}$ now imposes a preferred direction, which we can define to be the $z$ direction. If initially $\boldsymbol{m}$ is not parallel to $\boldsymbol{B}$, the effective field that a spin sees is the superposition of $\boldsymbol{B}$ and the mean field $\boldsymbol{m}$, and the Lindblad operators must be chosen according to the orientation of that effective field. The dynamical equation of $\boldsymbol{m}$ in this more general case will be treated further below in Section 4. However, as $\boldsymbol{m}$ will become parallel to $\boldsymbol{B}$ after a relaxation time, the dynamical equation (20) describes what happens in the subspace relevant for the bifurcation.

4 Two-dimensional case

When in equation (6) $z=0$ is an attracting value for the $z$ dynamics irrespective of the values of $x$ and $y$, we obtain an effectively two-dimensional system. A sufficient criterion for this to happen is that $h_{11}=h_{22}$ and $h_{23}=-h_{13}^{*}$ and $H=0$. Then the system of equations (6) becomes (if we set $z=0$)

[TABLE]

From a physical point of view, a decaying value for $z$ means that there is no equilibrium polarization along the $z$ axis. The probabilities for measuring spin-up in the $z$ direction becomes identical to that for measuring spin-down, namely 1/2. In the geometric view of the Bloch sphere, this means that the dynamics becomes restricted to a circle of radius one, namely the intersection of the Bloch sphere with the $z=0$ plane.

Now, there are two qualitatively different possible dynamical scenarios in the vicinity of a fixed point: If the Jacobian matrix corresponding to equations (21) has two real eigenvalues $\lambda_{1}<\lambda_{2}<0$, we can call the eigenspaces corresponding to $\lambda_{1}$ and $\lambda_{2}$ the fast and slowly decaying directions. After some time, the dynamics along the fast direction has relaxed, and we have an effectively one-dimensional system along the slow direction. This is similar to the situation given in equations (7) above. When a bifurcation occurs, the larger of the two eigenvalues goes through zero, which means that the bifurcation occurs within the one-dimensional subspace given by the slow direction, and all results obtained in the previous section can be applied also to this two-dimensional situation. In the degenerate case where the two eigenvalues coincide (i.e., when the fixed point becomes a "star" [14]), the dynamics follows a straight trajectory with the distance to the fixed point fully characterizing the dynamics, i.e., the dynamics is again effectively one-dimensional [14].

The second scenario is the one where the two eigenvalues of the Jacobian matrix are complex conjugate. In this situation, the dynamical system can exhibit a Hopf bifurcation, where a stable spiral becomes unstable and a limit cycle (a closed, isolated trajectory) is generated. The normal form of the Hopf bifurcation is

[TABLE]

or, in polar coordinates,

[TABLE]

If we impose the condition $|\epsilon|<1$, the value of $r$ remains in the interval $[0,1]$. The bifurcation occurs at $\epsilon=0$. There are no restrictions on the rate of angular increase $b$.

In the following we show that this normal form can be obtained by choosing suitable expressions for the $h_{ij}$ that ensure furthermore that the matrix $h_{ij}$ is positive semi-definite. For the diagonal elements, we set $\Gamma=a+r^{2}$ and $h_{11}=h_{22}=2h_{33}=\frac{\Gamma}{2}$. We will see that for sufficiently large $a$ the matrix $h_{ij}$ is positive semi-definite. For the nondiagonal elements, we write $\text{Re}[h_{12}]=\kappa$ and set $\text{Im}[h_{12}]=0$. Together with the choice

[TABLE]

and

[TABLE]

equations (21) take the normal form (4).

The ansatz for the diagonal elements $h_{ii}$ ensures that the coefficient matrix $h$ is positive semi-definite for all $x,y\in[-1,1]$ if $a$ is large enough compared to $|\kappa|,|\epsilon|,|b|$. A sufficient criterion is that all principal minors of $h$ are non-negative [15]. Hence, the following inequalities must be satisfied for all $x,y$:

[TABLE]

By considering the leading terms in powers of $a$, it can be verified that these conditions are indeed fulfilled when $a$ is large enough. Here, the evaluation of $\text{Det}(h)$ is the most complicated one as not only the highest order in $a$, which is $a^{3}$, but also the next order, which is $a^{2}$ must be considered if the conditions shall be satisfied for all $r\in[0,1]$: We have

[TABLE]

For $r=1$ the first term vanishes, and we must require that the second term is positive. This is satisfied for all $r\in[0,1]$ if $\epsilon<1$, for sufficiently large $a$. the system one-dimensional. In the above example for a limit cycle, the rotating frame has an angular velocity of $b$. Within this rotating frame dynamics is one-dimensional with the relevant variable being $r$, which approaches the fixed point $r^{*}=\sqrt{\epsilon}$.

Next, we construct again a microscopic model that undergoes such a bifurcation. To this purpose, we let the system of coupled spins interacts with an environment that produces a magnetic field in response to the magnetization of the spin system, with the field pointing in a direction that is rotated with respect to that of the magnetization. We make again the mean-field approximation. The Hamiltonian for one spin is thus again (compare (18))

[TABLE]

where we have defined an effective field $\boldsymbol{B}_{\text{eff}}=\boldsymbol{B}+\frac{\alpha}{\mu}\boldsymbol{m}$. Without loss of generality we choose the $x$-$y$ plane to be the plane spanned by the initial state of $\boldsymbol{m}$ and $\boldsymbol{B}$. The strength of the field $\boldsymbol{B}$ depends on $r$ and should vanish for $r=0$. We therefore set $|\boldsymbol{B}|=\text{const}\cdot|\boldsymbol{m}|+o(|\boldsymbol{m}|)$.

Its direction is given by $\boldsymbol{e}_{\phi+\delta_{1}}$, with $\phi$ being the direction in which $\boldsymbol{m}$ points, and $\delta_{1}$ being the angle by which $\boldsymbol{B}$ is rotated with respect to $\boldsymbol{m}$. The effective field $\boldsymbol{B}_{\text{eff}}$ is therefore also rotated with respect to $\boldsymbol{m}$, but by a smaller angle that we call $\delta$. Such a sustained rotation of the magnetic field constitutes an active driving of the system, and it therefore requires an energy input into the device that measures $\boldsymbol{m}$ and generates the field $\boldsymbol{B}$.

We assume again that the medium in which the spins are embedded has a finite temperature $T$, and that the interaction with the phonons of the medium causes spin flips between the two energy eigenstates of the spin. We denote the two rates by $\gamma_{f-}$ for the flip from the lower-energy state to the higher-energy state and $\gamma_{f+}$ for the reverse flip,

[TABLE]

The sum of the two spin-flip rates is $\gamma_{f}$. In addition to the spin flips, the interaction with the phonons can also result in a measurement of the energy eigenstate, and we label the rate of this process with $\gamma_{m}$, where the letter $m$ stands for "measurement". The Hamiltonian ((25) does not only affect the flip rates, but it enters also the von Neumann term. In contrast to the previous examples, the density matrix is not diagonal in the eigenbasis of $\hat{h}$, which means that the von Neumann term does not vanish but causes a precession of the spin around the direction of $\boldsymbol{B}$ . This makes the system three-dimensional. For the following analytical calculation, we therefore make the assumption that the rates $\gamma_{m}$ and $\gamma_{f}$ are so large that the unitary dynamics according to the von Neumann term can be neglected. In this case the vectors $\boldsymbol{m}$ and $\boldsymbol{B}$ stay in the $x$-$y$ plane. But we will additionally show below the result of a numerical integration of the full Lindblad equation with the von Neumann term in order to illustrate that the dynamics still converges to a limit cycle, but with its center shifted vertically. The density matrix of the system has then the form

[TABLE]

with $r$ and $\phi$ being the representation of $\boldsymbol{m}$ in polar coordinates in the $x$-$y$ plane.

In order to obtain the Lindblad equation that describes the dynamics of the density matrix $\rho$, we need to calculate the two operators that flip the spin between the + and - orientation with respect to an angular direction $\boldsymbol{e}_{\tilde{\phi}}$ in the $x$-$y$ plane, and the operator that performs a projection on an energy eigenstate. The latter one was not needed in the one-dimensional model as the density matrix was diagonal in the basis of the Hamiltonian. But when a spin is not in an eigenstate of $\hat{h}$, the interaction with the heat bath (phonons) will bring it into a state with definite energy, in addition to inducing flips between the two possible energy eigenstates, depending on whether an interaction event with the heat bath leads to the absorption or emission of a phonon, or not. The required rotation operator is given by (with $\boldsymbol{n}$ being the unit vector parellel to the axis of rotation)

[TABLE]

For a spin that flips between the $+x$ and $-x$ direction, the flip matrices are

[TABLE]

The measurement matrix is

[TABLE]

These are the matrices with respect to the direction ${\tilde{\phi}}=0$ in the $x$-$y$ plane. Performing a rotation by an angle ${\tilde{\phi}}$ gives the general measurement matrix with respect to the direction of $\boldsymbol{B}_{\text{eff}}$,

[TABLE]

and the flip matrix

[TABLE]

and the associated matrix $\textbf{F}^{\dagger}$.

The Lindblad equation with the three Lindblad operators M, F and $\textbf{F}^{\dagger}$ is

[TABLE]

with

[TABLE]

When we set ${\tilde{\phi}}=\phi+\delta$ and consider the nonzero matrix elements of the Lindblad equation, we obtain

[TABLE]

Separating the real and imaginary part on the left- and right-hand side and setting $\gamma_{f+}-\gamma_{f-}=\Delta\gamma_{f}$ and $\gamma_{f+}+\gamma_{f-}=\gamma_{f}$, this gives the two equations

[TABLE]

For $\delta=0$, the radial equation reduces to (14) for the one-dimensional system, as then $r=z$. Since $r$ is positive, we obtain only the positive branch of (14).

Since we have now assumed that $\boldsymbol{B}$ depends on $\boldsymbol{m}$ and that it vanishes when $\boldsymbol{m}$ vanishes, $\Delta\gamma_{f}$ must be proportional to $r$ to leading order in $r$ ($\Delta\gamma_{f}=\gamma_{f+}-\gamma_{f+}=\gamma_{f}\,\tanh\left[\beta|\alpha\,\boldsymbol{m}+\mu\,\boldsymbol{B}|\right]=\gamma_{f}\,\tanh\left(\beta\,c\,r\right)=\gamma_{f}\,\beta\,c\,r+\mathcal{O}(r^{3})$ ), and $r=0$ is always a fixed point.

If the fixed point at $r=0$ is unstable, we have a stable fixed point $r^{*}\in(0,1)$, which is determined by the equation

[TABLE]

In this case we obtain

[TABLE]

The angle $\phi$ increases at a constant rate. For the temperature $T$ (or, equivalently, the value of $\beta$) at which the stable stationary solution for $r$ becomes nonzero, a Hopf bifurcation occurs. As a side result, we have shown that the change of the direction of $\boldsymbol{m}$ is towards $\boldsymbol{B}$. If the magnetic field was constant in time, $\boldsymbol{m}$ would relax towards it. This completes the remark made at the end of section 3.

Figure 6 shows the trajectory obtained by a numerical integration of the model for a value of $\beta$ large enough to obtain a limit cycle. When performing this simulation, we additionally included the von Neumann term in the Lindblad equation. We see a limit cycle that is shifted in the $z$ direction. Without the von Neumann term, the limit cycle would be exactly in the $x$-$y$ plane.

5 Three-dimensional case

We now look at the full system of equations (6) with no constraint forcing the dynamics onto a lower-dimensional manifold. A nonlinear three-dimensional dynamical system can exhibit the same types of bifurcations as the lower-dimensional systems, but in addition it can make the transition to chaos and show a strange attractor [14]. In general, we can expect that equations of the type (6) display also strange attractors. In order to demonstrate this explicitly, we construct a specific model that is similar in spirit to the previous ones. We use again a system of interacting spins in mean-field approximation coupled to a heat bath and subjected to a magnetic field the direction and strength of which depends on the magnetization $\boldsymbol{m}$. As in the previous section, the interaction with the phonons causes the spins to make transitions to and between their energy eigenstates.

The main idea now is to choose the function $\boldsymbol{B}(\boldsymbol{m})$ such that it makes $\boldsymbol{m}$ move along a trajectory that bears similarity to that of a classical strange attractor.Similarly to the previous example, such an ongoing change of $\boldsymbol{B}$ in response to $\boldsymbol{m}$ leads to the dissipation of energy in the system and can be sustained only by providing an energy supply. Since $|\boldsymbol{m}|$ is confined to the interval $[0,1]$, it is important to choose parameters such that $|\boldsymbol{m}|$ is not close to 1 at all times since it will then hardly respond to a change of the strength $B$ of the magnetic field and can then not show a chaotic trajectory. We choose the cyclically symmetric Thomas model [16] as inspirational source and set

[TABLE]

At any moment in time, we can define a two-dimensional plane spanned by the vectors $\boldsymbol{m}$ and $\boldsymbol{B}$. During an infinitesimal time interval $dt$ the dissipative dynamics according to the flip and measurement matrices cause a change $d\boldsymbol{m}$ of $\boldsymbol{m}$ towards $\boldsymbol{B}$, with the angular change of $\boldsymbol{m}$ in this plane and the radial change being given by the equations (4). Since the plane within which this motion takes place changes with time, the equations of motion (4) must now be transformed to a form that is independent of the coordinate system. We have $r=\sqrt{x^{2}+y^{2}+z^{2}}=|\boldsymbol{m}|$, and the direction of angular change is that of the vector $(\boldsymbol{m}\times\boldsymbol{B})\times\boldsymbol{m}$. We therefore write (4) in the form

[TABLE]

Here, $\phi$ is no longer an angle in the $x$-$y$ plane, but in the plane spanned by $\boldsymbol{m}$ and $\boldsymbol{B}$. With this definition of $\phi$, the expressions for $\dot{r}$ and $\dot{\phi}$ are given by (4) and the flip rates by (26), as before.

We also want to consider the case that the von Neumann term cannot be neglected. We have already given a component-wise version of the von Neumann term as the last terms in equations (6). Evaluating these terms with the Hamiltonian (25) gives the final version of Lindblad equation (remember that $|\boldsymbol{m}|=r=\sqrt{x^{2}+y^{2}+z^{2}}$ and $\boldsymbol{B}_{\text{eff}}=\boldsymbol{B}+\frac{\alpha}{\mu}\boldsymbol{m}$)

[TABLE]

This dynamics of the density matrix with the choice (5) for $\boldsymbol{B}(\boldsymbol{m})$ was evaluated by numerical integration.

Figure 7 shows the trajectory of $\boldsymbol{m}$ obtained with a parameter set that leads to chaos when the von Neumann term is sufficiently small compared to the dissipative terms in the left graph, and a limit cycle when the von Neumann term is larger.

6 Discussion and Conclusion

We have studied the dynamics of Lindblad equations where the transition rates are nonlinear in the density matrix and have constructed explicit models that show such a dynamics. This approach raises the interesting issue of how the density matrix should be interpreted. In the conventional ensemble interpretation, the Lindblad equation must be linear in the density matrix. Often, the density matrix is interpreted as a device to reflect our limited knowledge of the precise state a quantum system, which is assumed to be a pure state. This applies in particular to the ensemble description of statistical mechanics systems, where the probabilistic description is often viewed as consequence of our ignorance of the exact quantum state of the system. In this case, the time evolution must be linear in the density matrix since the different pure states of the ensemble cannot interact with one another.

Decoherence theory, which is often taken as the background theory for deriving the Lindblad equation, is based on a similar assumption: here, the system is considered as part of a larger system that includes an environment with which it interacts, and the dynamics of this larger system is taken to be unitary. Due to the interaction, the system becomes entangled with the bath. Assuming a sufficient amount of randomness or uncorrelatedness of the environmental degrees of freedom, the environmental state that co-occur with the difference system states, can be argued to be orthogonal to each other. This is the central assumption of decoherence theory. The full quantum state of system (S) and environment (E) thus has the form $\sum_{n}c_{n}\left|\psi\right>^{(\text{S})}_{n}\left|\Psi\right>^{(\text{E})}_{n}$ with orthogonal $\{\left|\Psi\right>_{n}^{(\text{E})}\}$. The reduced density matrix of the system, which is obtained from the full density matrix by taking the trace over the environmental states, is then that of the ensemble of the $\left|\psi\right>^{(\text{S})}_{n}$, and therefore the time evolution of this reduced density matrix is linear. This type of consideration is also applied in approaches to the quantum measurement problem [17] and in the field of quantum foundations of statistical mechanics [18].

Although decoherence theory is by some authors considered as an explanation of the quantum-classical transition and in particular of the measurement problem [17], others emphasize that it cannot explain why only one of the different possible outcomes is realized in a one-time run of the measurement experiment [19, 20], since the full state of system and environment is assumed to be an entangled state that contains all the $\left|\psi\right>^{(\text{S})}_{n}$ at the same time. Furthermore, it is argued that a unitary time evolution of macroscopic systems is incompatible with the actual calculations performed in statistical mechanics and condensed matter theory [6], and that the quantum-classical transition indicates that there are limits of validity to unitary time evolution [7].

It is therefore apparent that the assumptions that lead to a time evolution that is linear in the density matrix are not applicable in all situations. In particular the idea that the system, possibly taken together with its environment, is in a pure quantum state, can be challenged, since in most experimental situations, the preparation of a quantum mechanical state cannot be controlled in all detail. One reason for this is the nonzero temperature of the preparation device. In this situation, it is impossible even in principle to prepare a pure state that can be verified experimentally, for instance by quantum state tomography [21]. When considering this imprecision of a quantum state as fundamental, the density matrix does not only reflect our limited knowledge but is the best description of the system we can possibly give. This also means that the density matrix now can be understood as describing a single system and not an ensemble of systems. Interpreting the density matrix as describing a single system is also advocated by other authors [22].

Furthermore, when not just the preparation, but also the subsequent time evolution of the quantum system is subject to uncontrollable stochastic influences, this density matrix does not evolve according to the von Neumann equation but requires a Lindblad equation (or a non-Markovian equation, but we do not consider this situation here). Since now the density matrix represents a single system and not an ensemble, the Lindblad equation is no longer required to be linear in the density matrix, and nonlinearities can arise under suitable circumstances.

In this paper, we have obtained nonlinearities in the Lindblad equation by applying a mean-field approximation to the interactions between spins. When performing a mean-field approximation on a quantum system, one usually replaces quantum mechanical operators by their expectation value, with the result that the influence of the other quantum particles take the shape of an external potential (which is a feature of classical physics). The nonlinearities arose in our model because this mean-field approximation made the one-particle Hamiltonian dependent on the density matrix, and this Hamiltonian in turn enters the von Neumann term as well as the temperature-dependent transition rates. We think that the reason why mean.-field like approximations are so widespread and successful is that they capture features of reality, which is quantum only up to certain (temperature-dependent) length scales. In our calculations, this approximation included furthermore the neglect of spatial correlations, making the one-dimensional versions of our model formally equivalent to the mean-field theory of ferromagnetism: According to the Landau theory of phase transitions, the free energy of a uniaxial ferromagnet in the presence of a magnetic field $h$ is $F=\,r\,m^{2}+u\,m^{4}-hm$, with $m$ being the order parameter (magnetisation) and with the parameters satisfying $u>0$ and $r\propto(T-T_{C})$. Relaxation dynamics toward equilibrium takes the form

[TABLE]

which corresponds to (14) for $h=0$ and to (20) for $h\neq 0$. Mean-field theories are a useful tool to capture by an analytical calculations the phenomenology of the overall behavior. Since spatial correlations are neglected, they have their limitations, as they cannot describe phenomena such as the formation of domains, of spatiotemporal patterns or spatiotemporal chaos, or how fluctuations modify the dynamics. In addition to the pitchfork and saddle-node bifuractions, one-dimensional systems also show a transcritical bifurcation, where two fixed points exchange their stability, as one of them moves into the physically feasible coordinate region. An example for a simple model that shows this bifurcation is the Laser model by Haken (see for instance [14]). As it is different in nature from the finite-temperature mean-field spin models on which we focused here, we did not include the transcritical bifurcation in Section 3.

In order to obtain limit cycles and chaos, we had to introduce a feedback between the magnetic field $\boldsymbol{B}$ and the entries $(x,y,z)$ of the density matrix. If a system of coupled spins in a magnetic field is left to itself and its interaction with a heat bath, an equilibrium state with detailed balance must be reached, and this is not compatible with a sustained change of the magnetization of the coupled spin system. Sustained changes require an ongoing driving of the system away from equilibrium, and this implies that the feedback mechanism requires an ongoing energy supply. This is a comparable situation to the quantum state transformation experiments described in the Introduction. There, the apparatus that prepares and measures or selects the quantum states is also a classical device that is connected to a power supply.

One can also think of other ways to drive the quantum system to a periodic of chaotic trajectory. For example, one could apply fields that do not change in time, and add an active feedback process that prevents relaxation away from an unstable trajectory. Such a stabilization can be done in a minimally invasive way [23]. Another and somewhat trivial way to drive a quantum system on a periodic or chaotic trajectory would be to impose this trajectory by top-down control. If we place a quantum two-level system in an environment that has a magnetic field and a temperature that changes on a time scale that is much slower than the thermal equilibration of the quantum system, we obtain,

[TABLE]

In such a situation, the density matrix is well defined at every moment in time, as is also argued by authors working in the field of NMR [24]. A system that leads to a chaotic equation of motion for the density matrix of the spin can be built in the following way: First, one constructs a electromechanical device that leads to a chaotic equation of motion in spherical coordinates $r$, $\theta$, $\phi$ of the end point of an arm of variable length and orientation. Then one mounts a magnet on the end of this arm such that the field points inwards in radial direction, and one places the quantum spin system in the center at $r=0$. The direction and strength of the magnetic field at the location of the spin are then fixed by the angle $(\theta,\phi)$ and radius $r$ respectively. This means that the equation of motion of the electromechanical device translates into an equation of motion of the magnetic field at the center, and this in turn translates into an equation of motion of the type (6) for the density matrix of the spin.

The last example has demonstrated most clearly that the density matrix is determined by the environment of the quantum system. In this example, the environment is given by a classical device that determines the magnetic field and the temperature to which the quantum system is exposed. There is no other way in which a quantum system can be controlled or influenced apart from changing classical, macroscopic control parameters, which then can cascade down to the quantum system. When the environment is simple, for instance a heat bath that does not change in time, this top-down determination of the density matrix is often overlooked, although temperature is a thermodynamic variable that is imposed by the environment. In the examples discussed in Sections 4 and 5, there was additionally a bottom-up influence from the quantum system to the classical environment, as the classical environment responded to the magnetization of the system. In quantum-state transformation protocols, there is also such a feedback from the quantum system to the classical system, as the postselection depends on the measurement outcome. The difference to our examples is that in our mean-field model the measurement is done on all interacting spins of the system, which leads to a classical measurement result (the vector $\boldsymbol{m}$), and not to the binary outcome of the measurement on a spin-1/2 system.

To conclude, the interaction of a quantum system with a classical system can lead in various ways to a violation of the unitary dynamics of wave functions as well as to the linear dynamics of the density matrix, both of which are only valid under restricted conditions and assumptions. Since our world is full of nonlinear dynamical phenomena, linear dynamics can only be a special case, even though a lot of research focuses on this case. It is our belief and that of various other authors [25, 26, 27] that the top-down influence from the classical world on the quantum world is an irreducible feature of nature that must be taken into account when one wants to find a solution to the puzzles surrounding the interpretation of quantum mechanics.