On the rate of convergence of an exponential scheme for the non-linear stochastic Schrödinger equation with finite-dimensional state space

We address the numerical solution of the finite-dimensional non-linear stochastic Schrödinger equation, which is a locally Lipschitz stochastic differential equation modeling, for instance, quantum measurement processes. We study the rate of weak convergence of an exponential scheme that reproduces the norm of the desired solution by projecting onto the unit sphere. This justifies the use of the Talay-Tubaro extrapolation procedure in the numerical simulation of open quantum systems. In particular, we prove that an Euler-Exponential scheme converges with weak-order one, and we obtain the leading order term of its weak error expansion with respect to the step-size. Then, applying the Talay-Tubaro extrapolation procedure to the Euler-Exponential scheme under consideration we get a second-order method for computing the mean values of smooth functions of the solution of the non-linear stochastic Schrödinger equation. We also prove that the exponential scheme under study has order of strong convergence 1/2, which gives theoretical support to the use of the multilevel Monte Carlo method in simulating open quantum systems. We present a numerical experiment with a quantized electromagnetic field in interaction with a reservoir that illustrates the good performance of the weak second-order method, and the multilevel Monte Carlo method.


Introduction
This paper develops the numerical solution of the stochastic Schrödinger equation describing the evolution of open quantum systems, which is also called the quantum state diffusion model and the non-linear Belavkin equation. We consider a small quantum system with state space  (see, e.g., [1][2][3] for derivations of (1)). The complex SDE (1) has a unique strong solution (see, e.g., [1,4]), which satisfies ]. The evolution of quantum measurement processes is described by the operator X X t t ñá | | , which is written using Dirac notation. Namely, X X t t ñá | |represents the density operator conditioned on the measurement outcomes (see, e.g., [1][2][3] for some derivations). In the homodyne or heterodyne measurement of the observable L t k ( ), for instance, the integral from 0 to t of the photocurrent is proportional to W tr L s X X ds W X L s X ds 2 2 , .
On the other hand, the Gorini-Kossakowski-Lindblad-Sudarshan quantum master equation describes the evolution of the density operator ρ t , which is the unit-trace positive operator acting on h h that represents the quantum state of the physical system. We have that X X t t t r = ñá (| |) (see, e.g., [1,5] for details, and [6] for a recent development). That's why it is said that the quantum density operator ρ t is unraveled in the quantum trajectories X X t t ñá | | . If d is not small enough (in the range of tens), then the calculation of the mean value of the quantum observable A (i.e., tr A t r ( )) by computing X AX , t t á ñ ( )via Monte Carlo methods has great advantages over the numerical solution of (2) (see, e.g., [2,3,[7][8][9][10]), which involves d 2 unknown complex functions. Hence, (1) is widely used for the numerical solution of (2), together with the quantum-jump version of (1) (see, e.g., [2,3])).
This paper deals with the computation of the mean value of X t j ( ), when j is smooth. We focus on the following Euler-Exponential scheme developed by [4,11] that solves accurately (1) with low computational cost (see, e.g., [2,7,11,12] for alternative schemes). First,we obtain that the rate of convergence of Y N j ( (ˆ)) to X T j ( ( )) with respect to the maximum step-size of n n N 0, , is equal to 1 for any function j whose partial Wirtinger derivatives up to order 4 have polynomial growth. An example of such j is x x Ax , ,which appears in the computation of both the variance of 〈X t , AX t 〉 and the average of the quantum variance of the observable A at time t,i.e., X A X X AX , , ). Previously, [4] proved that Y A Y , x + has compact support. This excludes to take W W n k k k 1 , a choice used, for instance, in the multilevel Monte Carlo methods (see, e.g., [13,14]). Here, we drop the assumption that the support of n k 1 x + is compact. At the same time, we establish that ( (ˆ))| is less than or equal to a constant multiplied by max n N n n 0, , )for a general class of functions j. This allows us to select adequate step-sizes τ n+1 − τ n for achieving a desired global error.
Second, we study for the first time, to the best of our knowledge, the application of the Talay-Tubaro extrapolation procedure to the numerical solution of (1). In [15] the Romberg-Richardson extrapolation method for ordinary differential equations is extended to solve globally Lipschitz real SDEs by weak secondorder schemes based on the Euler and Milstein schemes (see, e.g., [16][17][18]). In this paper, we obtain the leading term of the error , where Y N follows scheme 1. Then, we prove that the rate of convergence of obtained by halving the step-sizes of τ n . This yields a weak second-order approximation of X T j  ( )that inherits the good dynamical properties of scheme 1 and has a low computational cost.
We face two obstacles in the error analysis. Since the coefficients of (1) are not globally Lipschitz, the techniques currently used to estimate the rate of weak convergence of the numerical schemes for locally Lipschitz SDEs (see, e.g., [19][20][21][22][23]) do not apply to (1), at least no directly. In section 5.3.2 we overcome the first obstacle by introducing a localization procedure that allows us to treat (1) by using a globally Lipschitz SDE whose drift and diffusion coefficients coincide with those of (1) in the unit ball. To this end, we exploit the fact that that has a singularity at 0. In order to treat the singularity appearing in the projection onto the unit sphere, we derive a short-time asymptotic expansion of Y Y n n 1 1 + +  with respect to the step-size. Third, we study the strong convergence of Scheme 1, which ensures the convergence of the trajectories of ) . Namely, we prove that the rate of strong convergence of Scheme 1 is equal to 1/2. This, together with the linear order of weak convergence of Scheme 1, provides a theoretical basis for the numerical solution of (1) by the multilevel Monte Carlo method (see, e.g., [13]) with Scheme 1 as integrator. Other studies on the strong convergence of the numerical schemes for locally Lipschitz SDEs are, e.g., [24,25].
We organize the paper as follows. Section 2 sets up notation and basic definitions. Section 3 addresses the rates of strong and weak convergence of Scheme 1. Section 4 provides a numerical experiment with a quantum oscillator in the interaction representation. All proofs are deferred to section 5. Section 6 presents conclusions.

Notation
The scalar product , á ñ · · is anti-linear in the first variable. The partial Wirtinger derivatives are, by definition, , and the complex-valued function ) has partial derivatives of first order (see, e.g., [26,27]). For each n Î  we consider the multi-index set n , , : 0 , 1 , a n d 4 Here and subsequently, z z z , = á ñ   , and we use the same symbols K (·), K, and q for different nonnegative increasing functions, real numbers, and natural numbers, respectively. A family of functions f , and θ ä Θ.

Rates of strong and weak convergence
From lemma 3.1, given below, we obtain that we can project Y n+1 onto the unit sphere, and so Scheme 1 is welldefined. Proof. Deferred to section 5.2. , Next, we establish that the rate of weak convergence of Scheme 1 is equal to 1 provided that j and the partial Wirtinger derivatives of j up to order 4 have polynomial growth. From [20] we have that the Euler scheme could converge weakly with a rate slower than 1 when applied to SDEs with non-globally Lipschitz smooth coefficients. And worse still, the Euler scheme diverges, in the weak sense with fixed time horizon, in some cases where the SDE coefficients grow superlinearly (see, e.g., [21]). The property Y 1 n 1 = +   allow us to use a regularized Kolmogorov equation, together with X 1 t =   , to overcome the difficulties caused by the fact that the coefficients of (1) are not globally Lipschitz. On the other hand, we obtain Y 1 n 1 = +   by projecting Y n+1 onto the unit ball, and so the map Y Y n n 1 1 + + ˆhas a singularity at 0, which adds complexity to the convergence analysis.
] be continuously differentiable. Suppose that X is the solution to (1) with ), then there exists an increasing function K (·) depending on j such that for any time discretization n n N 0, , Here, K (·) does not depend on the time discretization n t ( ) and the natural number N.
Proof. The desired result follows from proposition 5.2 given in section 5.3. , In theorem 3.1 we use the partial Wirtinger derivatives -rather than the classical complex derivative-to describe the regularity of j. This enables us to consider a wide range of functions j. For example, is a non-null matrix, has continuous partial Wirtinger derivatives of any order, while x a 〈x, Ax〉 is not a holomorphic function. Moreover, using the Wirtinger derivatives we re-write the Kolmogorov equation corresponding to an SDE on d  as two Kolmogorov equations on d 2  , and so its regularity properties follow from known results on the Kolmogorov equation with real variables. Next, we modify the proof of theorem 3.1 to obtain the leading term of the error To this end, we adapt the error analysis carried out in [15].
given by Scheme 1 with Y X ), then there exists a function T 0, , and an increasing function K (·) depending on j such that for any time discretization n n N 0, , Here, K (·) and Ψ do not depend on n n N 0, , Proof. Deferred to section 5.4. , From theorem 3.2 we deduce that the extrapolation procedure introduced by [15] applies to Scheme 1. That is, in the notation of corollary 3.1, we have that ( (ˆ( ))) is a second-order approximation of X T j ( ( )), which requires one and one-half times the amount of evaluations of Scheme 1.
]is obtained by halving the step-sizes of τ. Let Y n t ( ) and Y n t (¯) be given by Scheme 1 with discretization τ and t , respectively.
, then for any time discretization n n N 0, , Here, K (·) does not depend on the time discretization n t ( ) and the natural number N.
Using the short-time asymptotic expansion of Y n 1 + obtained in the proof of theorem 3.1 we now establish that Scheme 1 converges strongly with order 1/2. Theorem 3.3. Let X be solution of (1) and Ŷ given by Scheme 1 with ) . Assume that Proof. Deferred to section 5.5. , Remark 3.1. The implementation of Scheme 1 requires the computation of the matrix exponential G exp n n t D ( ( ) ) times a vector (see, e.g., [28] for a review). This calculation can be done by using Krylov subspace methods whenever G n t ( )is a sparse matrix, which is common in many physical systems. If d is less than a few thousand, we compute G exp n t ( ( ))by means of Padé approximations combined with a scaling and squaring strategy, and so we calculate only one matrix exponential when we take a sample from each integration step of Scheme 1.
where a † and a are the creation and annihilation operators. We recall that the domain of a † and a is x l k e x : , = -for any k > 0, and a e 0 = 0. The Number operator is equal to a † a. Thus, the evolution of the quantum system is described by the mean values of functions of X X ¥ ¥ | |, where we use the Dirac notation and X t ¥ is the solution of (1) with H =  and L k k =  for any k = 1,K,4. The existence and weak uniqueness of the regular solution to (1) with an infinite-dimensional state-space is proved in, e.g., [29,30].
The model under consideration describes, for instance, a single mode of a quantized electromagnetic field in interaction with a heat bath. In the Hamiltonian , the term a † a is the Hamiltonian of the quantum harmonic oscillator, where the energy origin has been redefined to eliminate the vacuum fluctuation energy, and a a i - represents a linear pumping produced by an electric field (see, e.g., [8]). Moreover, 1  and 2  add the damping caused by photon emissions, 3  induces pure dephasing, and 4  models photon gain (see, e..g., [2,31,32]). We approximate X t ¥ by the solution X t of (1) with is the projection on the linear manifold spanned by e k d : The numerical experiment of section 6 of [4] examines the computation of X a aX , t t á ñ  † by using Scheme 1, the explicit Euler scheme and a version of the implicit Euler method, where the two later ones are projecting onto the unit sphere in each time integration step. In [4], the performance of Scheme 1 outshone that of the other two numerical methods.
As in [4] we compute X a aX , We take d = 50. This selection allows us to get X a aX , ) of the adjoint quantum master equation In figure 1 the solid line shows the interpolate values of e e , t 6 6 á ñ  , which are the reference values for X a aX , It is worth pointing out that the numerical solution of (6) has serious drawbacks if d is large (see, e.g., [2,8,9]), and according to proposition 6.1 of [33] we have that for any p Î , reaches a peak value at approximately T 10 = 9.4. We run Scheme 1 with two constant step-sizes τ n+1 − τ n = 2Δ, and n n 1 ( (ˆ( ))). We get more accurate estimations of X a aX , when we use smaller step-sizes Δ. In figure 1 we see that both Scheme 1 and the Talay-Tubaro extrapolation method reproduce well the oscillatory behavior of the quantum mean value of the number operator, even we use large step-sizes. Moreover, figure 1 highlights the very good accuracy of the extrapolation method described by corollary 3.1. Figure 2 and table 1 provide the errors corresponding to the Monte Carlo simulation of Scheme 1 and the Talay-Tubaro extrapolation method, respectively. Here, the sample size is equal to 10 8 (i.e., M = 10 8 ), and Y Y ,..., (see, e.g., section 2.2.2 of [34]). Similarly, the endpoints of a confidence interval for the Talay According to figure 2 and table 1 we have that the extrapolation method, given by corollary 3.1, improves the accuracy of Scheme 1. Applying a non-linear regression analysis we find that the estimated convergence rates of the errors ò 1 and ò 2 are 0.88054 and 2.073, respectively, in good agreement with theorem 3.1 and corollary 3.1. This motivates to apply the Talay-Tubaro extrapolation procedure in the numerical simulation of open quantum systems.

Multilevel Monte Carlo method
We calculate X a aX ,   (7) and (8)     are independent random variables distributed according to the laws of Y j , respectively. We adjust automatically the parameters L, M 0 ,K,M L by using the MLMC described in section 3 of [13], with the Euler-Maruyama scheme replaced by Scheme 1. We define the cost function at level ℓwith ℓ 1 (resp. ℓ = 0) to be the number of steps taken by Y , ℓ · and Y , 1 Table 2 provides the error 3 ¯defined by

The QuTip toolbox
In this section we compute X a aX , using the current version of the open-source software Quantum Toolbox in Python (see, e.g., [35]). That is, we use the version 4.7.1 of the QuTiP library, which can be found at https://qutip.org. We solve (1) by means of the QuTiP routine ssesolve with the solver option 'platen', as we present in Code 1. Thus, we run an implementation of the stochastic Runge-Kutta method stated in the equation (7.47) of [2], which was designed by Platen (see, e.g., [17] for a version of this scheme that achieves the second weak order of convergence for globally Lipschitz real SDEs).
In table 3 we compare the Talay-Tubaro extrapolation of Scheme 1 with the QuTiP function ssesolve. We restrict the sample-size to 4 · 10 5 and the step-size to 9.4/1280, because the intensive memory usage of Code 1 results in the Pyton reboot in more computationally demanding cases (we use a computer with an Intel Xeon E5-2699 v4 (44) processor and 64 GB of RAM). In table 3, 4 D  ( )is defined by (7) with , where E 2 D ( )is given by (9). According to table 3 we have that the method given by corollary 3.1 outperforms the QuTiP function ssesolve in terms of accuracy. In this example, table 3 suggests that the stochastic Runge-Kutta method designed by Platen shows numerical instabilities when Δ 0.0294.  5. Proofs of the results of section 3

Preliminaries
For completeness, we present the extension of two classical theorems to the framework of the Wirtinger calculus. The following Taylor's formula is used, for instance, in complex-valued signal processing.
Proof. Using the chain rule for Wirtinger's derivatives (see, e.g., [26,27]) we get . As in the derivation of the classical Taylor's formula, using induction and applying integration by parts formula we obtain the theorem.
The following Itô formula for complex-valued semimartingales has been used in the study of conformal martingales.   [36,37]). Replacing X t ℓ and Y t ℓ by Z , respectively, we obtain the desire formula after basic algebraic manipulations involving (3) and (10). , , n n n n n n k m k n n n k 1 1 1 n n n n n n n k m n k n n 1 1 . ,

Proof of theorem 3.1 5.3.1. Asymptotic behavior of Scheme 1
In Scheme 1 the preliminary approximation Y n+1 is projected onto the unit sphere. The following lemma allows us to treat theoretically the singularity at 0 of the function z z z   .
and the lemma follows. , Next, we study the short-time asymptotic behavior of Y n 1 +  .
Notation 1. By 0 G we mean the σ-algebra generated by Y 0 . For a given discretization n n N 0, , G denote the conditional expectation given n G . We write n n for different collections of N random variables such that n  is n G -measurable, and K T n N 0, , 1, is a nonnegative increasing function that does not depends on the partition k k N 0, , x ʼs, and the supremum on where the sum over j k ¹ means the sum over all j k m , 1, , , .
More precisely, here the random variable n 1 , , .
, collecting terms we get the lemma.
In (12) we collect the terms with coefficients n Since the distribution of the n k x ʼs is symmetric about 0, 0

Localization procedure
We will associate (1) with a globally Lipschitz SDE whose coefficients coincide with the drift and diffusion coefficients of (1) in the (closed) unit ball. According to the smooth Urysohn lemma we have that there exists a function As a smooth cutoff function r  we can take, for instance, is the characteristic function of the ball centered at 0 with radius b 1 2 < < , and r r h where g, G, σ k are as in (1) and theorem 3.1. The solution of (1), with X ].Therefore, using the uniqueness of the solution to (14) we deduce that Z X where Z T s z , is given by (14). Then, u T 0, , , and u s z Z , ). Replacing f z ¶ a ( ) by the right-hand side of (3) yields Hence, and (17) we obtain (16), and (iii) completes the proof.
, .  = ¥ , d = 50, and T j = 0.94 j with j = 0, 1,K,10. We use the Talay-Tubaro extrapolation of Scheme 1 ( 2 D  ( )), and the QuTiP library ( 4 D  ( )).   (ˆ) | )| ( )( ) G use a localization procedure. As a preparatory step, we prove that the Euler-Exponential scheme converges with weak-order 1, which is important in its own right. In this direction, we improve the numerical analysis presented in [4] with the purpose of considering general smooth functionals of the solution, and random variables without compact support that simulate the Brownian motion increments. Moreover, we proved that the Euler-Exponential scheme converges with strong-order 1/2, which lead us to implement the multilevel Monte Carlo method to solve numerically the non-linear stochastic Schrödinger equation with finite state space. A numerical experiment, with a quantized electromagnetic field in interaction with a reservoir, illustrates the good performance of the weak second-order method based on extrapolation, and that the multilevel Monte Carlo method and the Euler-Exponential scheme working together produce satisfactory results.