Speed of Evolution and Correlations in Multi-Mode Bosonic Systems

We employ an exact solution of the thermal bath Lindblad master equation with the Liouvillian superoperator that takes into account both dynamic and environment-induced intermode couplings to study the speed of evolution and quantum speed limit (QSL) times of a open multi-mode bosonic system. The time-dependent QSL times are defined from quantum speed limits, giving upper bounds on the rate of change of two different measures of distinguishability: the fidelity of evolution and the Hilbert–Schmidt distance. For Gaussian states, we derive explicit expressions for the evolution speed and the QSL times. General analytical results are applied to the special case of a two-mode system where the intermode couplings can be characterized by two intermode coupling vectors: the frequency vector and the relaxation rate vector. For the system initially prepared in a two-mode squeezed state, dynamical regimes are generally determined by the intermode coupling vectors, the squeezing parameter and temperature. When the vectors are parallel, different regimes may be associated with the disentanglement time, which is found to be an increasing (a decreasing) function of the length of the relaxation vector when the squeezing parameter is below (above) its temperature-dependent critical value. Alternatively, we study dynamical regimes related to the long-time asymptotic behavior of the QSL times, which is characterized by linear time dependence with the proportionality coefficients defined as the long-time asymptotic ratios. These coefficients are evaluated as a function of the squeezing parameter at varying temperatures and relaxation vector lengths. We also discuss how the magnitude and orientation of the intermode coupling vectors influence the maximum speed of evolution and dynamics of the entropy and the mutual information.


Introduction
Among a wide range of multifaceted fundamental problems dealing with different aspects of dynamics of quantum systems, the importance of the problems related to the speed of evolution of a quantum state for quantum technologies such as quantum communications and quantum computations cannot be overestimated. In processes underlying these technologies, this speed generally determines how fast a given task can be processed.
The maximum speed of the evolution can be quantified using the quantum speed limit (QSL) time defined as the minimal time needed for a system to perform transition between predetermined initial and final (target) states. An important point is that quantum mechanics dictates QSL bounds for the minimal evolution time which is inversely related to the evolution speed. For the case of unitary evolution, the well-known Mandelstam-Tamm [1,2] and Margolus-Levitin [3][4][5] inequalities provide general limits on the speed of dynamical evolution and set a dynamical intrinsic time scale. Interesting alternative geometric derivations for both inequalities are obtained from the statistical distance between quantum states in [6].
In [7], the geometric approach was used to generalize the inequalities to the case of nonunitary dynamics of open quantum systems. Following a similar general approach, where quantum speed limits are derived as upper bounds on the rate of change of a geometric measure of distinguishability, different expressions for generalized QSL bounds were also obtained in [8,9] using quantum Fisher information and the relative purity. The QSL inequality for open quantum systems using the trace norm was derived and interpreted in [10]. More details on QSLs and their applications in the fields such as quantum metrology and optimal control theory can be found in, e.g., Ref. [11] (see also references cited therein).
Recently, in Ref. [12], the QSL approach was applied to solve the optimal control problem dealing with the minimum energetic cost required to implement a quantum gate. Speed limits on the informational measures such as the von Neumann entropy, maximal information, and coherence of quantum systems were introduced in [13] where the speed limit on the maximal information is used to obtain the minimum time required to erase the information of quantum systems via some quantum processes of interest.
Quantum navigation is another important class of controlled quantum dynamics whereby the objective is to transport one quantum state into another in the shortest possible time under the influence of uncontrollable perturbations. Problems of this kind can be regarded as the quantum counterpart of the classical Zermelo navigation problem of finding the time-optimal control that takes a ship from one location to another, under the influence of external wind or currents [14][15][16][17]. Analysis of the evolution speed is a necessary step toward solving the quantum Zermelo problem.
For the last half of the decade, the evolution speed of two-level open quantum systems interacting with both Markovian and non-Markovian environments has been the subject of intense studies [16,[18][19][20][21]. A systematic analysis of the most common QSL bounds in the damped Jaynes-Cummings model was performed in [22]. The problem of quantum metrology in the context of non-Markovian quantum evolution of a two-qubit system was explored in [23].
As compared to the case of finite dimensional quantum systems, for bosonic multimode quantum systems representing a family of open continuous-variable systems, the problem of evolution speed has not yet received a proper attention. In this paper, we intend to fill the gap.
In a recent paper [24], the evolution speed and various QSL times were computed for the special case of a single quantized cavity mode that undergoes Markovian dynamics. Obviously, in order to study the effects determined by dynamics of intermode correlations, we shall need to go beyond the scope of single-mode models. These effects are governed by the intermode interactions arising from both dynamic and environment mediated couplings.
For instance, when a two-mode boson system represents the polarization degrees of freedom of of a quantized radiation mode [25], these couplings may have a profound effect on dynamics of the averages of the quantum Stokes operators that underlies the depolarization of quantum light propagating in a fiber [26][27][28]. Note that a similar two-boson system was formulated in Ref. [29] to model two cavity modes interacting with a single three-level atom.
In this paper, we explore how the evolution speed of multi-mode bosonic systems is influenced by the dynamical regimes dictated by the intermode couplings. The plan of this paper is as follows.
Section 2 begins with a brief discussion of QSL. We will introduce the fidelity of evolution and the Hilbert-Schmidt (HS) distance as measures of the distinguishability of quantum states applicable to the special case of a pure initial state [24]. For the rates of change of both of these measures, the speed of evolution appears as the upper bound that defines the time-dependent QSL times. In Section 2.1, our starting point is the Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) master equation (the general form of quantum kinetic equations which preserves the complete positivity of the dynamics [52][53][54]) for the multi-mode bosonic system interacting with the thermal bath and we, following Refs. [43,55], present its exact solution derived by solving the dynamical equation for the normally ordered characteristic function. In Section 2.2, we deduce general formulas that characterize the temporal evolution of Gaussian states using the generalized two-point fidelity, the evolution speed and the covariance matrix. In Section 3, general theoretical results are applied to the special case of a two-mode bosonic system which can be regarded as a model describing the propagation of quantized polarization modes in an optical fiber [27]. For this system, the intermode interactions can be conveniently treated using a geometric representation based on the frequency and the relaxation rate vectors.
In Section 4, we present a number of closed-form analytical relations for the cases where these vectors are either parallel or mutually orthogonal. For the case where the system is initially prepared in the two-mode squeezed state, different dynamical regimes are studied depending on the squeezing parameter and temperature. Analytical results are used to compute the disentanglement time and dynamics of the purity, the fidelity and the HS distance. Squeezing the parameter dependence of the coefficients characterizing long-time asymptotics of the QSL times is evaluated at different temperatures and relaxation rates. We also study how the orientation of the frequency and the relaxation rate vectors influence the maximum speed of evolution and dynamics of the entropy and the mutual information.
Finally, in Section 5, we draw the results together and make some concluding remarks. Technical details on algebraic identities related to the Liouvillian superoperator and its conjugate are relegated to Appendix A. In Appendix B, a closed-form formula for the normally ordered characteristic function is used to obtain the expression for the kernel of the Green's function describing the evolution of the Glauber-Sudarshan P function.
This distance satisfy the inequalities (see, e.g., the book [56]) that involve the Uhlmann fidelity F U (ρ(0),ρ(t)) = || ρ(0) ρ(t)|| 2 1 . Clearly, the inequalities (1) can also be regarded as relations linking the temporal evolution of the trace distance and dynamics of the Uhlmann fidelity that defines the geometric characteristics of density matrix operators such as the Bures distance and the Bures angle.
An important point is that in our case, the initial stateρ(0) is pure state and the Uhlmann fidelity can be simplified to the form: which is identical to the so-called fidelity of evolution, F (t). A two-point generalization of this fidelity reduces to the purity and the fidelity of evolution in the limiting cases with t 1 = t 2 = t and t 1 = 0, t 2 = t, respectively. It will also be used to compute the speed of evolution: This speed determines the upper bound for the rate of change of the fidelity (2) where a dot stands for the derivative with respect to time, which is an immediate consequence of the well-known Cauchy-Schwarz inequality for the inner product of Hilbert-Schmidt operators: where Â 2 = Tr |Â| 2 is the Hilbert-Schmidt norm.
Interestingly, the Hilbert-Schmidt (HS) distance, ρ(t) −ρ(0) 2 = 1 − 2F (t) + P (t), can also be used as a measure of distinguishability, and it is not difficult to show that the speed of evolution (5) provides the upper bound for the magnitude of its rate of change Since F (t) ≤ P (t) and 1 − F (t) ≤ ρ(t) −ρ(0) 2 , we have the time-dependent quantum speed limit (QSL) times which are associated with Equations (8) and (6) and lead to the QSL inequalities: From Equation (11), it might be concluded that the quantum speed limit imposed by t HS is tighter than the fidelity-based one (see Ref. [24] for discussion). Clearly, all the above QSL times given by Equations (9) and (10) depend on the speed and the fidelity of evolution (see Equations (5) and (2)) along with the purity (see Equation (4)).
In the remaining part of this section, our task is to provide general analytical expressions for these dynamical characteristics. More specifically, in Section 2.1, we, following our previous studies [43,55], present an exactly solvable model of Lindblad dynamics and its solution in the form suitable for our purposes. Then, in Section 2.2, we apply the results to deduce general relations for Gaussian states.

Master Equation and Characteristic Functions
The general form of the master equation ∂ρ ∂t = Lρ (12) suggests that the temporal evolution of the density matrixρ representing the state of an open quantum system is governed by the Liouvillian superoperator, L. In this paper, we consider the thermal bath version of the Liouvillian superoperator describing the Lindblad dynamics of an N-mode bosonic system in the Markovian approximation. This superoperator is expressed in terms of two superoperators: the commutator CÂB given by and the the dissipator DÂB given by where the dagger denotes Hermitian conjugation;â † n (â n ) is the creation (annihilation) operator of the nth mode; [Â,B] =ÂB −BÂ and {Â,B} =ÂB +BÂ stand for the commutator and the anticommutator, respectively; Ω nm (K nm ) is the element of the frequency (dissipation) matrix, Ω (K); z T =¯h Ω 0 k B T is the dimensionless inverse temperature parameter, where Ω 0 is the bare frequency,h is the reduced Planck constant, k B is the Boltzmann constant and T is the temperature of the environment. The frequency and dissipation matrices are both Hermitian: Ω = Ω † and K = K † . The dissipation matrix is also positive definite: K > 0.
It is well known that averages of the displacement operators,D(α) andD N (α), given byD where and an asterisk will indicate complex conjugation, define the symmetrically and normally ordered characteristic functions, χ(α) and χ N (α), as follows where Ô = Tr{Ôρ}.
The well-known Fourier-Weyl formula (see, e.g., the books [57,58]) in combination with the orthogonality relations can be used to obtain the identity linking the HS scalar product of the density operators,ρ 1 andρ 2 , and the inner product of their characteristic functions, χ 1 (α) and χ 2 (α). In what follows, this identity will be employed to evaluate the purity and the evolution fidelity needed to characterize the speed of evolution. The dynamics of the normally ordered characteristic function, χ N , is governed by the evolution equationχ with the differential operator,L † , determined by action of the conjugate of the Liouvillian superoperator, L † , on the displacement operator, :D : (see Appendix A for details). This operator is given byL whereD n T = (e z T − 1) −1 is the mean number of thermal photons and Γ = (1 − e −z T )K is the relaxation matrix giving the rates of thermalization. By using the method of characteristics to solve the initial value problem for the dynamical Equation (24) with the initial condition we can derive the following expression for the normally ordered characteristic function [43,55]: where I N is the identity N × N matrix. In Appendix B, this result is used to deduce the expression for the kernel of the Green's function that determines time dependence of the Glauber-Sudarshan P function. From Equation (28), it is not difficult to obtain the relations describing the temporal evolution of the symmetrically ordered characteristic function, χ(α). This formula in combination with Equation (23) will be utilized in our subsequent calculations.

Evolution Speed of Gaussian States
In this section, we concentrate on the special case where the initial state is Gaussian and its characteristic function is of the following form By using Equation (30), the time dependence of the characteristic function is expressed in terms of the two matrices where A T is the transpose of A. Similar to Q 0 and P 0 , the matrices Q and P are Hermitian and symmetric, respectively. The exponent on the right-hand side of Equation (33) can be conveniently cast into the real-valued form where α is replaced by s = (x 1 , . . . , (37), the subscripts "R" and "I" indicate the real and the imaginary parts of a matrix, respectively: A R = Re A and A I = Im A.
The 2 × 2 block-structure of the matrix (37) with rearranged elements reads where T is the permutation matrix, and it is closely related to the block structure of the covariance matrix, Σ, which is defined in terms of the quadrature operators (see, e.g., Refs. [45,57] The relation linking the block matrices of Σ andG reads where σ 2 is one of the Pauli matrices given by Now, we pass on to computing of the generalized two-point fidelity given by Equation (3). To this end, we substitute the characteristic functions of the form (36) into the relation (23). After performing the resulting Gaussian integral, we derive the identity where d(s) = ∏ N i=1 dx i dp i , giving the expressions for the purity (4) and for the fidelity of evolution (2) We are also in a position to evaluate the squared speed of evolution by substituting derivatives of the determinant that enters Equation (42) with respect t 1 and t 2 into the relation An alternative way to compute v 2 (t) is to evaluate the integral for the squared time derivative of the characteristic function It leads to the following formula for the squared speed of evolution: We conclude this section with the remark that the algebraic identities valid for a symmetric matrix G where G (c) = (det G)G −1 is the matrix of cofactors, which can be used to show that differentiating the two-point fidelity (42) and computing traces that enter the last expression on the right-hand side of Equation (47) yield identical results for the evolution speed. With the help of Equation (43) giving the purity, the latter can be further simplified as follows:

Evolution of Two-Mode Systems
In this section, general results of the previous section will be applied to the special case with N = 2. Thus, we focus our attention on the evolution of two-mode bosonic systems. As it was discussed above, such systems play an important part in the quantum optics as they represent polarization degrees of freedom of a quantized photonic mode.

Purity, Fidelity and Speed of Evolution
For the two-mode system, the frequency and relaxation matrices can be written as linear combinations of the Pauli matrices where σ ≡ (σ 1 , σ 2 , σ 3 ); ω ≡ (ω 1 , ω 2 , ω 3 ) and γ ≡ (γ 1 , γ 2 , γ 3 ) is the frequency and the relaxation rate vector, respectively. Thus, the dynamical interactions of the modes are described by the frequency vector ω, whereas the environment-mediated intermode couplings are determined by the relaxation rate vector, γ. These vectors that might be called the intermode coupling vectors govern the regime of dissipative dynamics and will be parameterized using the angular representation: It is rather straightforward to obtain the matrix exponential for A(t) (see Equation (29)) in the following explicit form: where where n R = Re n, n I = Im n, |n R | 2 − |n I | 2 = 1, n R ⊥ n I . We assume that the system is initially prepared in the two-mode squeezed state where r is the squeezing parameter, andK + =â † 1â † 2 (K − =â 1â2 ) is the raising (lowering) generator of the su(1, 1) Lie algebra. For this state, the characteristic function is given by Equation (32) with the matrices For the evolved quantum state with the density matrixρ(t), the explicit expression for the matrix Q(t) (see Equation (34)) that enters the characteristic function (33) is whereas the result for the matrix P(t) (see Equation (35)) is given by where n = n ⊥ + n 3 e 3 , n ⊥ = (n 1 , n 2 , 0).
Our calculations can be divided into the three following steps: (a) formulas (57)-(60) are substituted into Equation (37) to obtain the matrix G; (b) by using this matrix, we compute the purity (43), the fidelity of evolution (44) and the speed of evolution (47); (c) the results are used to evaluate the QSL time for the HS distance (9) and the fidelity-based QSL time (10).

Mutual Information and Disentanglement Time
The covariance matrix Σ can also be readily computed from the matrix G. According to Equation (38), the matrix gives the following block structure of the covariance matrix: For Gaussian states, the symplectic eigenvalues of the covariance matrix (63) are given by [45,57]: where ∆ ± = det Σ 11 + det Σ 22 ± 2 det Σ 12 , determine the entropy and the mutual information where µ 1 and µ 2 are the one-mode symplectic eigenvalues for the reduced one-mode density matrices,ρ 1 = Tr 2ρ andρ 2 = Tr 1ρ , respectively. Note that the mutual information can be defined as the relative entropy, D(ρ||ρ 1 ⊗ρ 2 ), of the two-mode density matrixρ to the product of the reduced one-mode density matriceŝ ρ 1 ⊗ρ 2 , which is an important entropic measure of the degree of correlation between the modes. An important point is that it contains the classical part of correlations and thus cannot be employed as an indicator of non-local quantum correlations such as entanglement.
For two-mode Gaussian states, different measures of entanglement have been proposed. These include the entanglement of formation, the Bures distance, and the Gaussian measures of entanglement [59][60][61]. For our purposes, similar to [43], it is convenient to deploy a quantifier of bipartite entanglement in Gaussian states which is expressed in terms of the lowest symplectic eigenvalue, λ PT , of the partial transpose of the two-mode density matrix,ρ PT , which is given by One of such quantifiers is the logarithmic negativity [57] E N (ρ) = log 2 ||ρ PT || 1 = max{0, − log 2 λ PT }.
It implies that the entangled states meet the condition: λ PT < 1 and, in the regime of finite-time disentanglement, the solution of equation gives the disentanglement time t d .

Results
As we have seen, for the two-mode systems, dynamical and environment-mediated intermode interactions appear to be encoded into the frequency and the relaxation vectors, ω and γ. Thus, we can expect that temporal evolution will crucially depend on the orientation of these vectors. According to Refs. [28,43], the angle between ω and γ is one of the important factors determining dynamical regimes of depolarization and disentanglement. In order to understand the orientation-dependent effects, we begin with the simplest limiting case where γ ω and iω − γ = iΩ − Γ.

Dynamical Regimes at γ ω and γ ⊥ ω
To be more specific, we assume that the vectors are oriented along the e 3 axis, iω − γ = (iΩ − Γ)e 3 , so that the matrices iΩ − Γ and A(t) are both diagonal. After evaluating the matrices Q and P (see Equations (57) and (59)) we obtain the matrices G (see Equation (37)) andG (see Equation (62)). The latter is described by the 2 × 2 block matrices that determine the block structure of the covariance matrix (63). The two-point fidelity (42) can now be derived as an explicit function of t 1 and t 2 of the following form: Then, the relations give the purity (43) and the fidelity of evolution (44).
It is now rather straightforward to evaluate the derivatives (45) of the fidelity (72) with respect to t 1 and t 2 and deduce the following closed-form formula for the speed of evolution: Since det Σ ij = detG ij and ∆ ± = q 2 1 + q 2 2 ∓ 2p 2 , it is not difficult to obtain the symplectic eigenvalues (64) and the lowest eigenvalue for the partial transpose of the density matrix (67) In the long-time limit, all these eigenvalues approach the limiting value that, similar to the steady-state density matrix:ρ(∞) =ρ eq ∝ exp[−z T (â † 1â 1 +â † 2â 2 )], is determined by the inverse temperature parameter z T . Hence, at non-vanishing temperatures with n T > 0, the regime of finite-time disentanglement takes place, and the disentanglement time can be found by solving Equation (69).
At Γ = γ 3 = 0, we have q 1 = q 2 = q and the simplified expression for λ PT : λ PT = |p − q|, gives the solution in the closed-form as follows This result was originally reported in [33]. Figure 1 illustrates the behavior of the disentanglement time, t d , as the relaxation rate, Γ, increases approaching its limiting value γ 0 . It can be seen that when the squeezing parameter r is sufficiently high, the disentanglement time decreases with Γ, whereas an increase in Γ will reduce the time t d provided the parameter r is below its critical value, r c .
The critical value of the squeezing parameter, r c , is plotted against the mean number of thermal photons, n T , in Figure 2. From the r c vs. n T curve depicted in Figure 2, it can be inferred that the critical squeezing parameter grows with temperature.  As is shown in Figure 3a,b, the time dependencies of the purity (43) and the fidelity of evolution (44) (see also Equation (73)) reveal non-monotonic behavior when the squeezing parameter exceeds the threshold r c . By contrast, Figure 3c shows that the HS distance is a monotonically increasing function of time that in the long-time limit with t → ∞, approaches its asymptotic value: where F (∞) and P (∞) are the long-time asymptotics for the fidelity and the purity given by The long-time behavior of the HS distance and the fidelity-based QSL times, t HS and t F , (see Equations (9) and (10)), is determined by the coefficients that might be called the QSL long-time asymptotic ratios. Typically, the evolution speed decays rapidly with time, and thus, the nonlinear corrections for the time dependence of QSL times are small. So, we concentrate on the QSL time ratios, κ HS and κ F .  (a) (b) Figure 5. Dependence of long-time asymptotic ratios for (a) fidelity-based and (b) HS distance-based QSL times, t F (see Equation (10)) and t HS (see Equation (9)), on the squeezing parameter computed at different values of the relaxation rate ratio Γ/γ 0 for n T = 0.1 and ω 0 = 0.
(a) (b) Figure 6. Dependence of long-time asymptotic ratios for (a) fidelity-based and (b) HS distance-based QSL times, t F (see Equation (10)) and t HS (see Equation (9)), on the squeezing parameter computed at different values of the relaxation rate ratio Γ/γ 0 for n T = 0.45 and ω 0 = 0.
Referring to Figure 4, the low-temperature curves for the ratio of the fidelity-based QSL time plotted in Figure 4a demonstrate non-monotonic behavior with a pronounced minimum. By contrast, the curves for the long-time ratio κ HS shown in Figure 4b, after reaching a local maximum in the low squeezing region, decline as the squeezing parameter r increases.
In Figures 4-6, an increase in the relaxation rate Γ leads to reduction of both the ratios κ HS and κ F . An additional effect of the relaxation rate is to enhance the non-monotonic behavior of κ HS in the low squeezing region (see Figure 5b).
A comparison between the results depicted in Figures 4a and 5a shows that the minimum is suppressed at elevated temperature. Interestingly, as is shown in Figure 6, in the vicinity of the high-temperature limit with n T = 0.5, the curves for κ HS and κ F become remarkably similar. When the relaxation rate Γ is close to γ 0 , the shape of the squeezing parameter dependence of the long-time ratios reveals two dynamical regimes that take place in the well-separated low-squeezing and high-squeezing regions.

Dynamical Regimes at γ ⊥ ω
We conclude our discussion of the dynamical regimes with remarks about another limiting case where the relaxation vector γ is reoriented along the normal to the frequency vector ω as follows In this case, the matrix P (see Equation (59)) can be cast into the form: where, similar to the matrix Q given by Equation (57), the time-dependent part is solely determined by the matrix product A † (t)A(t). Therefore, it might be expected that the time dependence of A † (t)A(t) dictates the dynamical regimes. According to Equation (82), similar to the simple P T -symmetric twolevel system studied in [16], the value of γ − iω is either real or imaginary depending on the sign of the difference between Γ and |Ω|.
When Γ > Ω and ω = 0 (γ = By contrast, in the opposite case with Γ < Ω and γ = 0 (ω = describes an alternative regime associated with the presence of oscillations. Finally, the critical regime that occurs when Γ = |Ω| and ω = γ = 0 is represented by the matrix This expression implies the dissipative dynamics will start to slow down when the oscillatory regime is changed to the critical one.

Effects of Intermode Couplings
As it was mentioned above, the evolution speed typically evolves monotonically approaching zero and the speed evaluated at initial instant of time, t = 0, gives its maximum value. This value is one of the factors that influence the QSL times and dynamical regimes. Figures 7 and 8 show how the maximum value of the evolution speed depends on the orientation of the frequency and relaxation vectors in the e 1 -e 3 plane. The surfaces plotted in Figures 7 and 8 represent angular distributions for v(0) computed as a function of the polar angles θ Γ and θ Ω at the two values of the relaxation rate ratio: Γ/γ 0 = 0.5 and Γ/γ 0 = 0.9, respectively. For both values of the relaxation rate, we have evaluated the low-temperature (n T = 0.1) and high-temperature (n T = 0.45) angular distributions.
At Γ/γ 0 = 0.5, referring to Figure 7a, it can be seen that the depth of modulation of the speed v(0) induced variations in the angle θ Γ is noticeably smaller than that induced by the angle θ Ω . As is shown in Figure 7b, an increase in temperature will raise the speed and make the θ Γ -induced modulation of v(0) more pronounced. When the relaxation rate is increased up to Γ/γ 0 = 0.9 (see Figure 8), we arrive at similar conclusions. In this case, both the temperature-induced speed up of the evolution and the θ Γ -induced modulation of v(0) are noticeably stronger as compared to the results presented in Figure 7.
The results for the entropy and the mutual information are depicted in Figure 9a,b, respectively. As is shown in Figure 9b, the mutual information monotonically decays, whereas the time dependence of the entropy (see Figure 9a) is non-monotonic. It starts from zero and reaches its maximum value, which is well above its equilibrium value S eq . Referring to Figure 10, by contrast to the entropy curves presented in Figure 10a, the time dependence of the mutual information shown in Figure 10b

Conclusions
In this paper, we have adopted a QSL-based approach to dynamics of a multi-mode bosonic system governed by the thermal bath Lindblad master equation. Within this approach, we have focused our attention on the speed of evolution and related timedependent QSL times. It is found that the evolution speed gives the upper bound for the rate of change of both the fidelity and the Hilbert-Schmidt distance (see Equations (6) and (8)). So, the HS distance and the fidelity-based QSL times, t HS and t F , (see Equations (9) and (10)) are both inversely proportional to the time average of the evolution speed.
General analytical results are applied to the special case of the two-mode bosonic system. In this case, the dynamics is determined by the intermode couplings that enter the dynamical and relaxation parts of the Lindblad superoperator L (see Equation (13)) and can be conveniently described in terms of the frequency and the relaxation rate vectors (see Equation (51)), ω and γ, that might be referred to as the intermode coupling vectors.
For the system initially prepared in the two-mode squeezed state, we have analyzed dynamical regimes assuming that the intermode coupling vectors are parallel, ω γ e 3 . In general, these regimes depend on a number of factors such as the squeezing parameter, the intermode interactions and temperature.
It turned out that different regimes can be related to the disentanglement time, t d , (see Equations (69) and (76)) which is found to either increase or decrease with the relaxation rate, Γ, equal to the magnitude of the relaxation vector depending on the value of the squeezing parameter (see Figure 1). The critical value of the squeezing parameter, r c , separating the low-squeezing region, where t d grows with Γ, from the high-squeezing one is found to be an increasing function of temperature (see Figure 2). (In this paper, we use the temperature-dependent mean value of thermal photons n T instead of temperature.) Figure 3 shows that in the region above the threshold r c , the time dependence of the purity and the fidelity becomes non-monotonic, whereas this is not the case for the HS distance.
The starting point of our study of the regimes associated with the QSL times was the long-time asymptotic behavior of the time-dependent QSL times characterized by the proportionality coefficients, κ HS and κ F , that enter the linear relations linking the QSL times, t HS and t F , and the actual time t (see Equation (81)). These coefficients called the large-time asymptotic ratios are plotted against the squeezing parameter at varying temperature and the relaxation rate Γ in Figures 4-6. Interestingly, the curves for κ HS with Γ close to the overall relaxation rate γ 0 reveal a well-pronounced local maximum in the vicinity of the above squeezing threshold r c . Note also that at high temperature, the curves for κ HS and κ F are similar in shape (see Figure 6). Referring to Figure 4, this is clearly no longer the case when the temperature is low.
For mutually orthogonal intermode coupling vectors, according to Equations (84)-(86), the dynamical regimes are dictated by the sign of the difference between the vector magnitudes: Γ − Ω. In a more general case, where the vectors are arbitrarily oriented in the e 1 -e 3 plane, we have evaluated the maximum speed of evolution in relation to the polar angles θ Γ and θ Ω (angular parametrization of the vectors is given by Equation (52)). The surfaces representing angular distributions of this speed are plotted in Figures 7-8. These distributions demonstrate that the orientation of the intermode coupling vectors has a profound effect on the evolution speed.
Similar to the evolution speed, from Figures 9 and 10, it can be inferred that the temporal evolution of the entropy and the mutual information is sensitive to variations in the polar angles. Note that the dynamics of the mutual information, I(1; 2), qualitatively resembles the dynamical behavior of the logarithmic negativity, E N , previously reported in [43].
The difference is that the decay rate of I(1; 2) is typically noticeably slower as compared to the decay rate of E N . In addition, E N vanishes at t > t d , whereas the mutual information falls off approaching zero in the long-time limit, t → ∞. These differences can be attributed to both the classical part of intermode correlations quantified by I(1; 2) and the nonclassical correlations that differ from entanglement (such correlations can be characterized using the quantum discord that describes the nonlocal fraction of I(1; 2)).
We conclude with the remark that our analytical approach provides tools for a more comprehensive analysis of the evolution speed and QSL-related problems for continuous variable systems in a realistic setup involving open environments and interacting modes. We hope that our work will stimulate further progress in the field.

Conflicts of Interest:
The authors declare no conflict of interest. that follows from Equations (A6) and (A7), we obtain the differential operator (25) that governs the dynamical equation for the normally ordered characteristic function χ N .
Our concluding remark is that in the relation the images of the superoperatorL and of its conjugateL † can be derived as an immediate consequence of the above algebraic identities.
On the other hand, the relation gives the characteristic function of the initial state χ (ini) N expressed in terms of P(β, 0). The Green's function is defined as follows P(α, t) = dµ(β)G(α, β, t)P(β, 0), where G(α, β, t) is the kernel of the superoperator.