A NUMERICAL APPROACH FOR THE STOCHASTIC CONTROL OF A TWO-LEVEL QUANTUM SYSTEM

The aim of this paper is to study the problem of controlling the stochastic evolution of a two-level quantum system in the presence of two randomly fluctuating electromagnetic fields, given by a Wiener process. The system is modeled by the stochastic Schrödinger equation dependent on time. We set up the quantum optimal control problem by choosing a cost functional type Bolza. By applying the Pontryagin Maximum Stochastic Principle to an extended Hamiltonian, we express the stochastic optimal controls in terms of the co-state of the system. To solve numerically the resulting stochastic differential equations we propose an iterative algorithm using the Euler-Maruyama method. Finally, we obtain the optimal trajectories on the Bloch sphere.


Introduction
Studying the motion of a 1 2 -spin particle in a two-level system and the interaction with its surronding, sometimes is necessary to introduce a stochastic model. We consider the Brownian motion of the spin represented by a randomly fluctuating electromagnetic field [Kubo and Hashitsume, 1970]. The efficient computation of numerical solutions of quantum stochastic diffusion equations, computing quantum trajectories in heterodyne measurement [Campagne-Ibarcq, P. et al., 2016], tracking optimal quantum trajectories in Bloch sphere [Weber, S.J. et al., 2014] and controlling the stochastic evolution of quantum system in presence of random perturbations or fluctuations [De Vega, 2005] are the motivations for the present work.
Stochastic simulation algorithms are used to study many phenomena with stochastic noise present, as the dynamics of stock prices or the motions of a atomic particle subject to external random fluctuations. Many of the models correspondig to these phenomena, which are described by stochastic differential equations, don't have analytic solutions.
There is no analytical solution to the stochastic model studied in this paper, so numerical approximations are required. There are several different approaches to numerically solve this kind of problems, namely, the gradient approach, the general Monte Carlo simulation, the Runge-Kutta approximations, the Markov chains approximations and the Euler-Maruyama approximations, among others. Using time discrete approximations, it is important to considerate the efficiency of a numerical scheme, by identifying the type of convergence that make practical sense to the model. In deterministic optimal control problems within of Nuclear Magnetic Resonance, several algorithms have been used, based on gradient approach, which have shown high performance in problems related to efficiency of coherence in the transfer spin states. In this work, we have extended to case stochastic, the algorithm used in the deterministic model studied in [Romero, 2017]. So, we don't have considerate the possibility of use this gradient approach. The Monte Carlo simulation is very general and accurate, but it is highly sophisticated and it requires large sample sizes. Given that our stochastic model was con-sidered as a perturbed system in the diffusion term of the corresponding deterministic system [Borkar, 2005], we have used the Euler-Maruyama scheme. The Euler-Maruyama approximation is the simplest method which can be implemented and which converges to an appropriate step size, although, it is sometimes inefficient and gives poor stability properties. In our case, this numerical scheme, by its simplicity, allowed us to develop a simple program in MATLAB to simulate sample paths close to those of the solution of the deterministic system. The simplicity of the program developed allows fast convergence of the trajectories. We consider a two-level open quantum system interacting with an environment, althought an electromagnetic field. When a sample is placed in a constant and longitudinal static magnetic field with magnitude B 0 in the direction of the Z axis, the magnetic moments of this sample are aligned. Then, if it is exposed to a randomly time varying radio frequency (r.f.) electromagnetic field with magnitude u 1 (t) along the X, and a r.f electromagnetic field with magnitude u 2 (t) along the Y axes, varying also in a random way, the sample absorbs the energy through a sequence of transverse magnetic pulses. The total magnetic field to which the sample is subjected is B(t) = u 1 (t) ı + u 2 (t)  + B 0 k. When the magnetic moment vector of the system is transferred to the XY plane, the sequence of transverse magnetic pulses is stopped, causing the magnetic moment vector to precess. Repetitions of this process produce fluctuations in B 0 and eventually, decoherence. The pulse sequence should be as short as possible to minimize the effects of relaxation, to optimize the sensitivity to the experiment and the contrast of the obtained image. This is achieved by controlling the sequence of pulses that create a unitary transformation in the shortest possible time. For Control Theory the minimization in time of a sequence of pulses equals the minimization of lengths of trajectories of vector states (in homogeneous spaces).
We study the simplest control system of a -1 2 spin particle interacting with a electromagnetic field, neglecting other interactions with the system. This quantum control system describes the dynamics of a system like a 2-level quantum system, governed by the Schrödinger equation for a pure state (we set = 1) where the state ψ : [0, T ] → C 2 is a vector representing the unitary ket |ψ , T ∈ R is the duration of the process, u : [0, T ] → R, u(t) = u 1 (t) + u 2 (t) is the control and the energy of the system is represented by the Hamiltonian H(t) that, in our case, is the interaction of the total spin angular momentum, including the interaction with the external electromagnetic fields B(t), which plays the role of the control u(t). The Hamiltonian satisfies where S = S 1 ı + S 2  + S 3 k is the spin angular momentum operator and γ is the gyromagnetic ratio of the system or constant of proportionality between the magnetic momentum and the angular momentum. The state vector in (1) is written as |ψ(t) = α|+ + β|− , where |+ and |− are the orthonormal eigenvectors corresponding to eigenvalues + 2 and − 2 , respectively, of S 3 . So, in the {|+ , |− } basis, the matrix representing S 3 is S 3 = −i 0 0 i . In the same way, we trices {S 1 , S 2 , S 3 } are the Pauli matrices, they form a basis for su(2), Lie{S z , S x , S y } = su(2), making sense of the controllability problem for the system (1), [D'Alessandro, 2001]. Unlike our previous paper [Romero, 2017], we introduce disturbance in the system (1), considering the components B x (t) = u 1 (t) and B y (t) = u 2 (t) of the electromagnetic field varies randomly in time, and B 0 (t) is the constant B 0 = 1. So, using the 1 2 S 1 , 1 2 S 2 , 1 2 S 3 Pauli matrices, representing the operators S 1 , S 2 , S 3 , we write We can separate the Hamiltonian H(u(t)) into H 0 (u(t)), corresponding to and H 1 (u(t)), corresponding to Therefore H 0 (u(t)) = −γ(S 1 u 1 (t) + S 2 u 2 (t) + S 3 B 0 ) (6) H(t) contains deterministic and stochastic interactions so, we can separate them as where H 0 (t) is the free system Hamiltonian corresponding to the external electromagnetic field (4) and the random Hamiltonian H 1 (t) of the system under the influence of one stochastic perturbation (S 1 u 1 (t) + S 2 u 2 (t))s(t), with s(t) assumed a stationary and Gaussian process, acting by means of the field (5) Rescaling time, considering B 0 = 1 and denoting S 1 = − 1 2 γS 1 , S 2 = − 1 2 γS 2 and S 3 = − 1 2 γS 3 , we have where s(t) can be represented as a Brownian motion process or a Wiener process W (t) by dW (t) = s(t)dt. Therefore, the stochastic linear Schrödinger equation for pure states corresponding to equation (1) is: ) For a better visual description of the evolution of the trajectory of the controlled states ψ = (ψ 1 , ψ 2 ) ∈ C 2 we will use the surface of the Bloch sphere S B ⊂ R 3 . An adequate projection of the space of the pure states over the Bloch sphere is given by the Hopf projection Π : with ψ * 1 the conjugate of ψ 1 in C, [Boscain, 2006].

The Quantum Optimal Control Problem
We consider the model describing a spin 1 2 -particle in an electromagnetic field fluctuating in the X, Ydirections by white noise. The optimal control problem for this quantum system is the following: to the final state ψ(T ) in C 2 and minimizes over 0 ≤ t ≤ T the energy cost functional of Bolza type, following: where E(f ) denotes the conditional expectation with respect to f , T = π √ 2 and O is the observable operator with the target information O = ψ(T ) ψ (T ), which will allow an optimal evolution of the system.
We split into real and imaginary parts of ψ = ( ψ 1 , ψ 2 ) and we set x = (Re ψ 1 , Re ψ 2 , Im ψ 1 , Im ψ 2 ) , Re S 3 Im S 3 Im S 3 Re S 3 and also for S 1 , S 2 . We obtain the whole dimensional extended system: where, again S 1 , S 2 and S 3 are denoted by the real matrices:

Adjoint Process
To find optimal controls in the quantum optimal stochastic problem (14), subject to a cost functional (13), we will use the stochastic Pontryagin maximum principle, which is formulated following the pseudo-Hamiltonian or Stochastic Control Hamiltonian associated with the stochastic general system and his adjoint stochastic linear system is called the stochastic adjoint system associated to (16) if the scalar product x(t) · λ(t) is a global stochastic first integral of (16), namely, if d( x(t) · λ(t)) = 0 The adjoint process has the following characterization: Theorem 1. Associated to the time-dependent Schrödinger linear equation in the stochastic control problem (14) and the corresponding variational stochastic system with ξ(t) ∈ T ψ R n , is the following adjoint backward stochastic linear system where λ(t) ∈ T * x R n is the adjoint vector of x(t) and λ T = (0, 0, 0, 1) Proof In general, we write the adjoint backward stochastic linear system Applying the transversality conditions of the Pontryagin Maximum Principle and the Itô product formula we obtain from (18) and (20) the following and then, proving equation (19) We use Theorem 1 to get the following expression for the Hamiltonian: The necessary conditions of optimality for the controls, according to the Pontryagin Maximum Principle, so, if u 1 (t) = 0 and u 2 (t) = 0, the condition reads and the Hamiltonian results Definition 2. An admissible controlū(t) is called a singular control on the control region V, if V ⊂ U is nonempty and for a.e., t ∈ [0, 1] and ∀v ∈ V, a.s., we have H(t, x(t,ū(t)),ū(t), λ(t,ū(t)), q(t,ū(t))) = H(t, x(t;ū(t)), v, p(t,ū(t)), q(t,ū(t))) So, if the Hamiltonian does not depend on the control, this is a singular control. In our case, any control u(t) = 0 is a singular control on U . Looking for optimal controls, when the Hamiltonian is invariant with respect to the control variables on some admissible control region, further necessary optimality conditions are needed aside from the Pontryagin Maximum Principle, namely the Second-Order Maximum Principle [Tang, 2010], [Yong-Zhou, 1999].

The Second-Order Stochastic Maximum Principle
Let U be a nonempty Borel subset of R m and K 0 a positive constant. We assume: are Borel measurable and continuous in x and for K 0 , ∀i = 1, . . . m A2) The first-order derivatives above are continuous in u onŪ . The functionsf , σ, l y h have continuous second-order derivatives in x. The second-order derivatives are Borel measurable with respect to (t, x, u), are bounded by the constant K 0 , that is, In this context we define: and now, we have the fundamental result: Theorem 2 (Second-Order Maximum Principle (Tang)). Assume that (A1) and (A2) are satisfied. Let (y(·), u(·) be an optimal pair and u(·) be singular on the control region V . Then, there is a subset I 0 ⊂ [0, 1] which is of full measure, such that at each t ∈ I 0 , (y(·), u(·) satisfy, in addition to the first-order maximum condition (22), the following second-order maximum condition ∀v ∈ V, a.s.: To apply the second-order maximum principle, we consider the second order adjoint equation, [Yong-Zhou, 1999]: From (26) and (32) we have and ∆b x (t, x, u)P (t) ∆b x (t, x, u) = 0 (35) so, the second order necessary minimum condition (30), means that any control u(t) = 0 satisfies the secondorder maximum principle, that is to say, that each singular control u(t) = 0 is a candidate for optimal control. Otherwise, the second-order maximum condition would not be satisfied and the constant control u(t) would necessarily not be optimal. We propose the following design for the control fields, based on the last considerations and according to [Romero, 2017]: Replacing the equations (37) into equations (14) and (19) we obtain a coupled two-point boundary value problem of stochastic differential equations, whose analytic solutions can not be found, so numerical solutions are needed.

Numerical Approximation
Several facts have motivated the increasing interest in the development of computational techniques related with optimal quantum control and, recently, with the stochastic optimal quantum control in the context of the magnetic resonance theory. Among them is the fact that in very few cases analytical solutions of the Schrödinger equation can be found, so numerical solutions are needed. In particular, the evolution of the Stochastic Schrödinger Equation numerically can be a useful tool by understanding the quantum systems interacting with the environment, namely, the open quantum systems. In this section we describe the proposed algorithm and prove its monotonic convergence. We focus in this work the monotonic iterative scheme.
The Euler scheme is the simplest effective computational method used in deterministic differential equations and the Euler-Maruyama scheme is it's analogue for stochastic differential equations. We have used this method to develop the next algorithm to determine the optimal controls to steer a two-level quantum system from an ground state, represented by the north pole of the Bloch sphere to an excited state, represented by the south pole of the Bloch sphere.

Algorithm outline
We extend the ideas in [Maday, 2003] and [Romero, 2017] to the stochastic case.
Given four constants δ 1 , δ 2 , η 1 , η 2 , ∈ (0, 2), λ 0 (t), v 0 (t), w 0 (t) real valued functions and a natural number k ≥ 1, we define the following iterative functions: and the following forward stochastic differential system: the backward stochastic differential system: The Euler-Maruyama method is obtained by truncating Itô's formula of the stochastic Taylor series after the first terms. It computes approximations x i ≈ x(i∆t). We will apply this method to solve the former Itô equations.
To solve the forward and backward processes (42) and (43) we select a grid of [0, T ]: So, we have, for i = 1, . . . N : where x i = x(t i ), λ i = λ(t i ) and we can consider with s i a random real number in [0, 1].

Continue until convergence
This algorithm produces ensembles of functions {λ (k) (t), x (k) (t), u (k) (t)} k≥1 that eventually converge to the solution of the system: (46) The efficiency of this algorithm depends on the adjustment of the initial parameters. We start by choosing δ 1 = 0.01, δ 2 = 0.40, η 1 = 0.01, η 2 = 0.40, v (0) = 0.589 sin(6 t + 0.41), w (0) = 0.589 sin(6 t + 0.41), λ 0 (t) = (0, 0, 0, 4.53). The process was convergent at k = 15 for these choices. Figure 1 shows the stochastic optimal controls u 1 and u 2 . Figure 1 shows the corresponding stochastic trajectory x = (x 1 , x 2 , x 3 , x 4 ) to the optimal controls. Figure  1 shows the projection of the optimal trajectories on the Bloch sphere and Figure 4 shows the comparative trajectory for the deterministic case.    (x1, x2, x3, x4) to steer |0 to |1 , using two controls found with the proposed algorithm. Figure 2. The stochastic optimal controls found u2(t) and u1(t) to steer |0 to |1 , using the proposed algorithm . 5 Controlling the system only with one stochastic control along the X-axis Let us consider a more detailed problem of controlling a nuclear spin subjected to a static magnetic field along the Z-axis and a single time varying radio frequency magnetic field along the X-axis. According to the previous equations, in this case we have (47) subject to the constraint: Now, the Hamiltonian of the system (47) is: and the necessary conditions for optimality of the controls get: Applying the algorithm adapted to this case, with δ 2 = 0.40, η 2 = 0.40, w (0) = 0.59 sin(6.1 t + 0.43) and λ 0 (t) = (0, 0, 0, 4.53), the process was convergent at k = 10 for this choice. Figure 5 shows state transfer limit trajectory (x 1 , x 2 , x 3 , x 4 ) and the stochastic optimal control found u 1 (t), that steers |0 to |1 , using the proposed algorithm. Figure 7 shows four iterations of the algorithm of the state transfer trajectories on Bloch sphere and the state transfer limit trajectory.
5.1 Controlling the system only with one stochastic control along the Y -axis Let us consider now the problem of controlling a nuclear spin subjected to a static magnetic field along the to |1 on the Bloch sphere, using the stochastic optimal control u1(t) for the algorithm, when u2(t) = 0. on the Bloch sphere, using the stochastic optimal control u1(t) for the algorithm, when u2(t) = 0. Figure 9. State transfer trajectory (x1, x2, x3, x4) to steer |0 to |1 , using the proposed algorithm. Figure 10. The stochastic optimal control found u2(t) to steer |0 to |1 , using the proposed algorithm. Figure 5. State transfer limit trajectory (x1, x2, x3, x4) that steers |0 to |1 , using the proposed algorithm. Figure 6. The stochastic optimal control found u1(t) that steers |0 to |1 , using the proposed algorithm.
Z-axis and a single time varying radio frequency magnetic field along the Y -axis. In this case we have the Figure 11. State transfer trajectory from |0 to |1 on Bloch sphere, using the stochastic optimal control u2(t) for the algorithm, when (51) subject to constraint: The Hamiltonian of the system (51) is: and the necessary conditions of optimality for the controls get: In Figure 10 appears the state transfer limit trajectory (x 1 , x 2 , x 3 , x 4 ) and the stochastic optimal control found u 2 (t) that steers |0 to |1 , using the proposed algorithm. In Figure 11 appears the state transfer limit trajectory on the Bloch sphere.

Discussion and Conclusion
In the proposed algorithm, the initial conditions play a very important role. On the one hand, we have the initial condition of the system of stochastic differential equations and, on the other hand, the initial value of the parameters δ, ν, and the initial functions λ 0 (t), v 0 (t) and w 0 (t). In the first case, small perturbations of the initial condition (ground state (1, 0, 0, 0)), do not affect the controls found or the corresponding optimal trajectory on the Bloch sphere, namely, by changing the initial condition slightly, the target (0, 0, 0, 1) remains still in a small neighborhood and we have also the fast convergence of the algorithm. On the other hand, conversely, small perturbations of the values δ and ν greatly affect the optimal trajectory, arriving, for some values, to compromise the convergence of the algorithm. The initial functions, λ 0 (t), v 0 (t) and w 0 (t), can eventually affect the optimality of the optimal trajectory. This algorithm can be used to steer any two points on the sphere of Bloch, since their convergence does not depend on the initial point or end point, see Fig. 12. We have discussed an Krotov-like algorithm, implemented with Euler-Maruyama method, to generate the stochastic optimal control of a state-to-state transition that maximize the cost functional Bolza type given, by equation (13). The advantage of this algorithm over other approaches is their simplicity to reach the target in few iterations. In the deterministic case, this algorithm is called algorithm of fast convergence. This algorithm is more related to the wave function formalism than the density-matrix approach, from the point of view of the optimal control time and the cost functional cost Bolza type [Tremblay, 2011]. The disadvantages of this algorithm are their high sensitive to choose of restricted parameters, the restriction of the step size and the need for a good reference adjoint function and the pseudo controls λ 0 (t), v 0 (t) and w 0 (t).
We are also interested in a kind of exponential and robust stability for the numerical scheme proposed in this paper and in his long-time asymptotic behavior , as an extension of stability defined in [Mao, 2015] and [Mora, 2017]. In a next paper, which is in preparation, we establish the conditions that must satisfy the algorithm to have these types of stability.
In this paper ,we have addressed the problem of finding one o two stochastic controls that enable to steer a quantum system with two levels from an initial state given to a final given state, minimizing the energy represented by a functional of Bolza type. Through the application of the Pontryagin Maximum Principle and a Stochastic Maximum Principle of the Second Order, we obtained in both cases that the proposed controls are optimal. These optimal stochastic controls allowed to steer the system of the initial position to the final position of the system, represented by the north and south pole of the Bloch sphere (Figure 3), respectively. The trajectory obtained is similar to the found in the deterministic case ( Figure 4).