{'fraction_non_alphanumeric': 0.04263915049272514, 'fraction_numerical': 0.016842737422542515, 'mean_word_length': 4.568323453412373, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimal probabilistic quantum control theory

There is a fundamental limit to what is knowable about atomic and molecular scale systems. This fuzziness is not always due to the act of measurement. Other contributing factors include system parameter uncertainty, functional uncertainty that originates from input functions, and sensors noises to mention a few. This indeterminism has led to major challenges in the development of accurate control methods for atomic scale systems. To address the probabilistic and uncertain nature of these systems, this work proposes a novel control framework that considers the representation of the system quantum states and the quantification of its physical properties following a probabilistic approach. Our framework is fully probabilistic. It uses the Shannon relative entropy from information theory to design optimal randomised controllers that can achieve a desired outcome of an atomic scale system. Several experiments are carried out to illustrate the applicability and effectiveness of the proposed approach.


Introduction
In recent years, quantum control theory has been acknowledged as an enabling tool for the development of new quantum technology and quantum information theory applications [1][2][3][4][5]. One of the main objectives in quantum control theory is to develop methods that can manipulate and control quantum systems. This has been achieved for the first time at the end of the last century [6][7][8][9][10][11][12] following the advent of laser femtosecond pulses. Consequent research studies investigated the development of optimal laser pulses as a control strategy to accomplish the required control objective [13]. These studies were based on the variation of the laser pulse shape to optimize the outcome of an experiment, i.e. the result of a predefined reaction product. This variation is repeated until the desired result is achieved [14][15][16][17][18][19].
For the design of control strategy and methods, optimal control theory [20,21], Lyapunov control approaches [22][23][24], learning control algorithms [25] and robust control methods [26][27][28][29][30] have been developed for the manipulation of quantum systems and the achievement of various control objectives. Among the aforementioned control design approaches, quantum optimal control is recognized as a powerful method for many complex quantum control tasks and has been successfully implemented for finding a control strategy for controlling molecules. In quantum optimal control, the control objective is usually casted as the optimisation of optimal cost function which is specified as the expected value of a target operator such as the projector onto a certain bound state or any other arbitrary operators [31][32][33][34][35][36]. This optimisation of the cost function is commonly subject to some constraints including a penalty term on the energy radiation and the satisfaction of the Schrödinger equation by the wave function of the system. The solution of this constrained optimisation problem results into solving coupled non-linear Schrödinger equations [35,36].
Hence, several numerical methods have been introduced to optimise the cost functional including the conjugate gradient method [37] and the Krotov iteration method [38]. However, it has been shown that most of these iterative methods are unreliable and computationally inefficient [35]. Consequently, Zhu, Botina and Rabitz have proposed the rapid monotonically convergent iteration method that solves the optimal control equations [35,36]. It was proven that this algorithm recovers the Krotov iterative method as a particular case which was also tested for quantum optimal control of population and for quantum optimal control over the expectation value of a positive definite operator [35,36]. In [35], the positivity of the target operator is shown to be mandatory for the assurance of the convergence of the cost functional. Since then, several extensions to this algorithm have been constructed and practically applied for the control of chemical reactions such as the feedback and the resonant excitation strategies [39,40].
Despite the success made so far in the area of quantum control, advances in this area are limited due to the high level of uncertainty, the effect of dissipation, and the probabilistic nature of atomic-scale physical systems.
A coherent control framework that can effectively address the aforementioned challenges is still lacking. This is mainly due to the fact that most of the existing advances in the field of quantum control are mostly based on designing deterministic controllers, thus they overlook the aforementioned challenges [41,42]. As such, this paper proposes a fundamental new probabilistic feedback control framework for quantum systems. The proposed framework is fully probabilistic. It uses the Shannon relative entropy from information theory to design optimal randomised controllers [43][44][45][46][47][48] that can achieve a desired closed loop performance of quantum systems under high level of uncertainty and stochasticity. To reemphasise, the derived control strategy under the proposed framework is fully probabilistic and it is based on the fact that the dynamics of a quantum system can be estimated using probability density functions (pdfs). We start by developing the general solution for the state space models of quantum systems whose state and output equations can be described by arbitrary probabilistic models. Using the vectorisation of the density operator obtained from the Liouville-von Neumann equation as the state equation and the measurements of the corresponding physical properties as the output equation, the solution to the proposed probabilistic control problem is then obtained and discussed in detail.
Here, it will be shown that the distribution of the vectorised state of the density operator and output equations can be described using complex Gaussian distributions thus facilitating the analytic evaluation of the optimal randomised controller. For these quantum systems described by the Liouville-von Neumann equation, the optimal control solution is obtained as a feedback control law, which is linear in the state and whose gain matrix satisfies a Riccati equation.
The proposed fully probabilistic control framework is different and more general to what has already been established in the literature. Its main characteristics are highlighted as follows: Firstly, it is a unified probabilistic control framework, all the components within this framework including the quantum controller, and quantum systems models are modelled using probabilistic models. This probabilistic characterisation of the individual components required to control quantum systems provides complete descriptions of their behaviours and depicts the inherent uncertainty in their dynamics. Secondly, it addresses the high level of uncertainty and the inherent probabilistic nature of atomic scale systems. Thirdly, although an analytic solution of the proposed method is only feasible for a vectorised quantum state whose time evolution is governed by linear and Gaussian distribution, the solution can be obtained in a closed form for any quantum system that can be described by arbitrary pdfs.
Contrary to existing deterministic control laws for quantum systems, the randomised controller obtained from the proposed probabilistic control framework is more explorative and provides complete information in the decision-making process, therefore is the natural solution to stochastic and uncertain atomic scale systems.
It gives an optimal control solution in a feedback form that can be implemented either in real time by taking measurements from the controlled quantum system, or offline by feeding back estimated state values from the estimated probabilistic model of the time evolution of the system. For the online implementation though, quantum measurement backaction needs to be taken into consideration or otherwise a weak measurement procedure where only partial information on the measured observable is obtained will need to be followed.
While it is interesting, this aspect of the study is beyond the scope of the current paper and thus further development will only focus on the establishment of the proposed probabilistic control framework.
The rest of this paper is organised as follows: Section (2) briefly recalls some preliminaries on the evolution of quantum systems and shows that their dynamics can be modelled using pdfs. It also states the objective of the considered control problem and provides its general solution. The proposed probabilistic control framework is then demonstrated in Section (3) on the bilinear representation of the vectorised time evolution of a quantum system described by the Liouville-von Neumann equation. In Section (4), we particularly apply the proposed probabilistic control framework to the Lithium hybrid molecule 7 Li 2 H and to particular spin systems interacting with external electric fields. Finally, some conclusions are provided in Section (5).

Quantum System Description
The time evolution of a quantum open system interacting with its environment can be described by the Liouville von-Neumann equation, where H 0 is the system's free Hamiltonian, andh is the reduced Planck's constant. The system can be controlled using the electric field u(t) where the interaction of the electric field with the quantum system can be described by the term µu(t) and where µ is an operator related to the system (e.g., electric dipole moment, polarizability).
In addition, ρ(t) ∈ C l×l is the density operator which is a positive hermitian operator with Tr(ρ(t)) = 1, i.e.
the eigenvalues of the density operator are interpreted as probabilities [49,50]. The first term on the right hand side of Eq.(1) is the part that describes a closed system and it is called the Hamiltonian part. The second term is associated with the coupling of the system with the environment that is responsible of the dissipation. In the Lindblad approach, the coupling term is provided by, where {.} stands for the anti-commutation operator, L s are the Lindblad operators and s runs over all dissipation channels. The operators L s are defined in terms of the dissipative transition rates Γ k→j from the free Hamiltonian eigenstate |k to the eigenstate |j as, In this equation, it is assumed that k takes a finite number of values, i.e., {|k , k = 0, ..., l − 1} with l being the number of the eigenvectors of the free Hamiltonian H 0 . In Appendix (A), we show that the Liouville-von Neumann equation (1) can be written as, where δ n,m is the Kronecker symbol, {n, m} = {0, 1, . . . , l − 1}, µ k,n := k| µ |n are the matrix elements of the operator µ, ω n,m := E n − E m h are the Bohr frequencies, with E n being the energy eigenvalue of the free Hamiltonian H 0 associated with the eigenvector |n , and γ n,m is the total dephasing rate defined by, γ n,m := 1 2 The solution of the Liouville von-Neumann equation provides the information on the evolution of physical properties of quantum open systems. The average value of an observableô at a time instant t is given by, ô(t) = Tr(ôρ(t)).

Control Objectives of the Quantum Control Problem
Although open quantum systems can be characterised by the time evolution of the Liouville von-Neumann equation, this evolution may not be completely known and may be subject to uncertainty. This could be due to uncertainties from uncontrollable experimental parameter variations, and functional uncertainty that originates from input functions such as the variations in the manufacturing process for example. Considering the definition given in Appendix (B) for the vectorised state vector, x t of the density matrix ρ(t), the development in this paper will be based on the characterisation of its time evolution by a probability density function, where as previously defined, u t is the electric field at time instant t. Note that because of the previously discussed uncertainties, the probabilistic description of the vectorised state vector of the density matrix as given in Eq. (7) provides a complete specification of the present state x t as a function of the previous state, x t−1 and previous control, u t−1 .
Similarly, observations of the observable,ô are subject to different sources of uncertainties including sensors noises, and measurement uncertainties. Thus, the probability density function of the observations provides the most complete specification of their values, The probabilistic description of the time evolution of the quantum system as given in Eq. (7), and the observations as given in Eq.(8) is general and can be characterised by continuously monitoring its underlying stochastic evolution. It is not constrained by the linearity or Gaussian assumption of the stochastic evolution of the system. This characterisation is taken as a ready methodology in this paper, thus is not discussed further.
Interested readers are referred to [51,52] on some of the available methodologies.
Following this formulation, the objective of the quantum control problem can be stated as follows: design a randomised controller, c(u t−1 |x t−1 ) that minimises the Kullback-Leibler divergence between the joint pdf of the closed-loop description of the quantum system, f (Z(t, H)), and a predefined ideal joint pdf, I f (Z(t, H)), where Z(t, H) = {x t , . . . , x H , o t , . . . , o H , u t−1 , . . . , u H } is the closed-loop observed data sequence and H ≤ ∞ is a given control horizon. The joint pdf of the closed-loop description of the system dynamics is the most complete probabilistic description of its behaviour. For the vectorised representation of the density matrix given in Eq.(7), it can be evaluated using the chain rule [53] as follows, where the pdf s(x t |x t−1 , u t−1 ) describes the dynamics of the vectorised density matrix, s(o t |x t ) represents the dynamics of the observations o t , and c(u t−1 |x t−1 ) represents the pdf of the required randomised controller as mentioned earlier.
Similarly, the ideal joint pdf of the closed-loop data can be factorised as follows, where the pdf I s(x t |x t−1 , u t−1 ) describes the ideal distribution of the state vector of the vectorised density matrix, x t , I s(o t |x t ) is the ideal distribution of the observations, and I c(u t−1 |x t−1 ) represents the ideal pdf of the randomised controller. The definition of these ideal distributions is analogous to the target values in [54] that the system is required to achieve.
Given the definitions of the joint pdf of the closed-loop system and its ideal joint pdf as specified by Eqs. (10) and (11), respectively, minimisation of Eq.(9) can then be obtained recursively by introducing the following definition, for arbitrary τ ∈ {1, ..., H}. In Eq.(12), − ln(γ(x t−1 )) specifies the expected minimum cost-to-go function, Following the same procedure of classical physical systems [43][44][45][46], using the definition of the expected cost-to-go function given in Eq. (12), minimisation is then performed recursively to give the following recurrence functional equation,

General Control Solution to the Quantum Control Problem
Following the representation given in Section (2.2) for the probabilistic state space models of open quantum systems described by the Liouville von-Neumann equation, the general solution for the optimal randomised controller that minimises the recurrence functional equation defined in Eq.(13) is given in the following proposition.

Proposition 1
The pdf of the optimal control law, c(u t−1 |x t−1 ), that minimises the cost-to-go function (13) can be shown to be given by, where and Proof 1 The derivation of the above result can be obtained by evaluating the optimal cost-to-go function specified in Eq. (13).
To reemphasise, the randomised control solution given in Proposition (1) provides a general solution for quantum systems that are affected by various sources of uncertainties as explained earlier. Furthermore, it is not restricted by the Gaussian assumption of the pdfs of the quantum system states and observations or their ideal distributions. It provides the general solution for any arbitrary pdf. However, as will be seen in the following section, if all of the generative probabilistic models of the system dynamics, controller and ideal outcomes are Gaussian pdfs, an analytic form for the randomised controller can be obtained.

Solution of the Quantum Control Problem for Gaussian pdfs
This section will demonstrate the bilinear representation of the vectorised time evolution of a quantum system described by the Liouville-von Neumann equation (1) and discuss its characterisation with probabilistic models.
The theory developed in the previous section will then be applied here to derive the analytic solution of this bilinear state space model.

Bi-Linear State Space Model of the density matrix and measurements
For quantum systems governed by the Liouville-von Neumann equation (1), that are also driven by external control field, u t , the state vector of the vectorised density matrix is shown in Appendix (B) to be given by the following bilinear equation [55], where x t ∈ C n is the vectorisation of the density operator ρ(t) ∈ C l×l with n = l 2 ,x is the state of the system at time t = 0, which is the initial condition of the differential equation (17), [55]. The elements of the matricesÃ,Ñ and q can be determined from the Liouville-von Neumann equation as explained in Appendix (B). In contrast to linear systems, the coupling of quantum systems with the external electric field induces a bilinear termB(x t ) that depends on the system's state x t . Discretising the state space equation (17), the discrete time state space representation can be obtained as follows, where, and where λ = ∆t − t with ∆t being the sampling period. We also have added the fictitious noise termw t+1 to account for any uncertainty associated in the time evaluation of the density matrix as per the discussion in Section (2.2). Equation (18) can be equivalently written as, Similarly, the time evolution of the measurements of the operatorô is given by, where D = vec(ô T ) T ∈ C (1×n) andṽ t is a fictitious noise that accounts for uncertainties caused by the measurement and environment.
Assuming that the noisew t affecting Eq. (21) is a Gaussian noise, the pdf of the system state, x t can be considered as a complex normal pdf, where, Here E(.) stands for the expected value, x † t =x T is the conjugate transpose of x t , and x T t is the matrix transpose of x t . The matrices µ t , and Γ are respectively the mean and the covariance matrices. Also, note that the characterisation of the system state, x t by complex normal distribution is due to the fact that the system state vector is complex. The form of the complex normal distribution is recalled in Appendix (C).
Similarly, assuming that the noiseṽ t affecting Eq.(23) is a Gaussian noise, the pdf associated with the measurement o t can be described by a standard normal distribution, where o mt = Dx t is the mean matrix and where, is the covariance matrix of the Gaussian distribution of the measurement vector, o t . It is worth mentioning that the variable o t is real and for this reason its distribution should be taken as a standard normal distribution.

Optimal Control Law
In Section (3.1), we have shown that the state vector of the vectorised density matrix described by Eq. (21) can be characterised by a complex normal distribution and the measurements described by Eq.(23) can be characterised by a standard normal distribution, as stated in Eqs. (24) and (26) respectively. Thus, the instantaneous joint pdf of the quantum system state, measurements and control signal can be written as, where c(u t−1 |x t−1 ) is the randomised optimal control to be derived. Accordingly, the ideal pdfs of the system state, the measurements and the controller are taken to be Gaussians as follows, Thus, the instantaneous ideal joint pdf can now be written as, The objective here is to find the distribution of the controller c(u t−1 |x t−1 ) provided in Eq.(28) that minimises the Kullback-Leibler distance between the joint distribution (28) and the predefined ideal one (32). The minimisation of the Kullback-Leibler distance as proposed in the current paper is the analogous objective functional to the deterministic one in [54] for obtaining optimal control laws for stochastic systems described by probabilistic models as given in Eqs. (24) and (26) .
The derivation of the controller distribution appearing in Eq.(28) that minimises the Kullback-Leibler divergence between the pdf (28) and the ideal one (32) is based on the evaluation of the optimal performance index, − ln(γ(x t )). The form of this controller will be stated shortly, but first the form of the optimal performance index is provided in the following theorem.
Theorem 1 By substituting the ideal distribution of the system dynamics (29), (30), the ideal distribution of the controller (31), and the real distribution of the system state (24) and measurements (26) into Eq. (15), the performance index can be shown to be given by, where, and where, In additions, in Eqs. (34)-(36) the following definitions were used, Proof 2 The proof of this theorem can be obtained by using backward induction as shown in Appendix (D).
We start by computing the coefficient β(u t−1 , x t−1 ) using the definition (16), which will then be followed by the computation of γ(x t−1 ) using the definition (15). This will yield the assumed form in the theorem.
Remark 1 From Eq.(33), it can be clearly seen that the optimal cost-to-go function, − ln(x t ) maps the complex vector x t ∈ C n to a real number, since M is an hermitian operator. This preserves the important necessity that cost functions are designed to be ordered and that their values can be compared in a consistent ordering manner.
Following the form of γ(x t−1 ) specified in Theorem (1), it is now straight forward to derive the pdf of the optimal control law that minimises the Kullback Leibler divergence between the actual pdf (28) and the ideal one (32). It can be obtained by substituting Eqs. (15) and (16) in Eq. (14) yielding the following theorem.
Theorem 2 The distribution of the optimal control law that minimises the Kullback-Leibler distance between (28) and (32) is Gaussian distribution given by, where with F t , C and J being provided in Eq. (37). In addition, v t−1 is the mean of the Gaussian distribution of the optimal control law and F −1 t is its variance.
Proof 3 The proof of this theorem is given in Appendix (E).
Remark 2 Although the derived probabilistic controller (38) for the assumed probability distributions maintains the standard form of linear quadratic controllers, it is more exploratory due to its probabilistic nature. For the implementation of this controller to real world systems, either the mean of the randomised controller can be taken to be the optimal control input, or control inputs should be ideally sampled from the obtained pdf of the randomised controller. Sampling control inputs from the obtained pdf however results in slightly worse control quality, but randomisation makes the controller more explorative.

Implementation Procedure of the Randomised Controller
The step by step implementation of the fully probabilistic control problem of quantum systems where the characterising pdfs are assumed to be Gaussians is provided as pseudocode in Algorithm (1). The execution of this algorithm facilitates the evaluation of the optimal electric field that drives the quantum system from its initial state to a predefined desired one. Use M t , and P t to compute the mean of the optimal control input, v t−1 following Eq.(39) given in Theorem (2);

12:
Using the obtained control input from the previous step, evaluate x t according to Eq.(21), Following Eq.(23), evaluate o t to find the measurement state at time instant t, o t ← Dx t +ṽ t ; 14: t ← t + 1;

15: end while
Algorithm (1) will be applied in the next section to control certain quantum physical systems.
4 Application to particular physical systems 4

.1 Generalities on Morse potential
Following [56,57], the Morse potential is a solvable model that very well describes the vibrations inside diatomic molecules [58,59]. For a diatomic molecule, the Morse potential is governed by, where r is the distance between the two atoms and r eq is the corresponding equilibrium position [56,57], while D 0 and α are related to the depth and the width of the potential, respectively. The parameter D 0 is related to the molecule properties by ν = 8mrD 0 α 2h2 , where ν is a parameter corresponding to the spectroscopic constants of the molecule and m r is the reduced mass of the oscillating atoms. The Morse potential is a solvable model and its energy eigenvalues are given by, where n = {0, 1, . . . This means that the Linblad operators appearing in Eq.(1) are now L s = 0 for any s. The interaction between the system and the electric field is described within the semi-classical approximation through the electric dipole µ ≡ µ(r) = µ 0 re −r/r * where µ 0 and r * are parameters related to the molecule [60]. The matrix elements of µ can be computed by, where {n, m} = {0, . . . , [p]} and ψ ν n (r) are the energy eigenfunctions associated with E n . They are given by, ψ ν n (r) = N n e − y 2 y j L 2j n (y), where we have used the change of variable y = νe −α(r−req) , and where L 2j n (y) are the Laguerre polynomials. In addition, 2j = ν − 2n − 1, and N n is the normalization factor governed by, where Γ is the gamma function [56].

Application to the Lithium hybrid molecule
Here we consider the Lithium hybrid molecule whose Morse parameters can be found in [60]. Their values are given by D 0 = 2.45090 eV, r eq = 2.379Å and ν ≈ 6.1346. Hence, in this case p = 2.5673, which means that the quantum number n, appearing in Eq.(41) of the Morse potential associated with the Lithium hybrid molecule takes three values 0, 1, 2. As explained in Section (1), the interaction between the electric field and the physical system is realised through the electric dipole operator µ which is now given by µ = µ(r) = µ 0 re −r/r * where µ 0 and r * are parameters related to the molecule [60]. Their values for Lithium hybrid molecule are given by µ 0 ≈ 5.8677 Debye and r * ≈ 1.595Å. By omitting the effect of the environment, the time evolution of the density operator (4) can be governed by, where in this example, {n, m} = {0, 1, 2}.
Although any target operator can be considered, here we consider the following Gaussian operator given in the position representation {|r } as follows,ô where γ 0 = 25 and r = 2.4871Å. The matrix representation of the target operator (46) can then be obtained as follows, Therefore, following the vectorisation method described in Eq.(B.2), the operator D that appears in the measurement equation (23) whereÃ is a diagonal 9 × 9-matrix. The shifting vector x e appearing in Eq.(B.4) is x e = x 0 with x 0 being the vectorisation of the density matrix at t = 0, considered here to be the ground state. Furthermore, Using Eqs. (19) and (20), we can easily find the matrix elements of A and B operators determining the evolution of the vectorised state x t at each time step t, as governed in Eq. (18). Taking ∆t = 0.0167, G r = 0.000000014, Γ r = 10 and Ω = 0.28950, we evaluate the matrices M t and P t defined in Eqs. (34) and (35), respectively at each instant of time as discussed in Algorithm (1). Then, we use these matrices to evaluate the control signal, u t−1 following Eq. (39). The obtained electric field is used in Eq. (21) to find the evolution of the state vector x t and the corresponding observation can be determined using Eq. (23). We repeat these steps until the measurement output, o t becomes as close as possible to the predefined desired value o r , which is taken to be equal to 1 in this example, i.e., o r = 1.

Spin systems
The introduced probabilistic control strategy will be implemented here to manipulate the transition probabilities of a number of examples of spin-j systems. A spin-j system can be described by, where c m are complex coefficients that satisfy j m=−j |c m | 2 = 1, and the ensemble {|j, m ≡ |m , m = −j, . . . , j} are the common eigenstates of the operators J 2 = J 2 1 + J 2 2 + J 2 3 and J 3 . In addition, J k with (k = 1, 2, 3) are some observables satisfying the angular momentum commutation relations [J 1 , J 2 ] = ih 123 J 3 (with 123 being the Levi-civita symbol). We consider that the spin-j system interacts with an external electric field and its evolution can be totally described by the following Hamiltonian, where H 0 is the system's free Hamiltonian and H u (t) is the operator describing the interaction between the system and the electric field u t ≡ u(t). The control objective here is to transfer the system from an initial state, |ψ i to a predefined final state |ψ f through its interaction with the optimised electric field. Hence, the problem can be seen as a maximisation of the fidelity between the initial and the final states; |ψ i and |ψ f , respectively.
This means that the target operator here is nothing but the projector Π f = |ψ f ψ f |. is o r = 1. Indeed, the problem of maximisation of the fidelity between the initial and the final states has been solved using different approaches such as the rapidly convergent algorithm in [36] and the method investigated in [61].

spin-1 2 system
The probabilistic control method discussed in Section(2.2) will be applied here to a spin-1 2 system interacting with an electric field u t . According to [61], the dimensionless Hamiltonian describing the interaction between the system and the electric filed is given by, where σ 1 , σ 2 and σ 3 are the Pauli matrices given in the basis We consider that the system is initially prepared in the state |ψ i = |− , and the objective is to transfer it to the state |ψ f = |+ through its interaction with the optimised electric field u t that is designed to achieve that objective. The calculation of this electric field at each time step is explained in Algorithm (1). In particular, using Eq.(23), the target operator is evaluated to be D = [0 1 0 0], which is nothing but D = (vect(Π + )) T , where Π + = |+ +|. Next, the values of the matricesÃ andÑ are computed as explained in Appendix (F) followed by the calculation of the operators A and B at each time instant using Eqs. (19) and (20)  considered spin-1 2 system interacting with the electric field. We can clearly see that the population ρ 11 has reached the predefined target value o r in few time steps. This demonstrates the effectiveness of the proposed control method. Figure (2b) shows the time evolution of the optimal electric field responsible of transferring the system from its initial state to the desired one.

Spin-1
To further demonstrate the effectiveness of the probabilistic approach introduced in this work, it is applied to a spin-1 system interacting with an electric field, u t . We consider that the interaction can be described by the following dimensionless Hamiltonian, The objective here is to transfer the system from an initial state to the desired state |1, 1 , which means that the target operator is D = (vect(Π 1 )) T where Π 1 = |1, 1 1, 1|. Appendix (G) demonstrates how to evaluate the matricesÃ andÑ associated with the spin-1 system whose evolution is described by Eq.(55), and gives their values. The first objective is to transfer the system from the ground state |1, −1 to state |1, 1 and the second objective is to transfer the system from the state |1, 0 to the state |1, 1 . Taking

Final comments
In this paper we have introduced a new probabilistic strategy to control physical systems at the atomic and molecular scales. The main particularity of our approach compared to the existing ones is that it is fully probabilistic and considers the probabilistic nature of the dynamics of quantum systems due to different sources of uncertainty and stochasity. Due to these uncertainties, we have shown that the dynamics of quantum systems can be characterised by probability density functions. the state of the system and the measurement of its physical properties. Furthermore, using the minimisation of the Kullback-Leibler divergence between the joint pdf of the actual state, measurements and controller and a predefined ideal joint pdf, we have provided the form of the control law that transfers the system from its initial state to the desired one. Moreover, we have applied the proposed approach to the Lithium hybrid and spin systems in interaction with an external electric field and shown that the state of the outcome associated with particular target operators can be controlled using the derived optimal control signal.

Appendices A Evolution of quantum open systems and state space model
The time evolution of a quantum open system interacting with its environment can be described by the Liouville von-Neumann equation (1). In Dirac notations, the free Hamiltonian H 0 can be written as, where E n are the energy eigenvalues of H 0 . In the vector space spanned by the energy eigenvectors, the density matrix can be written as, where ρ kj are the matrix elements of ρ(t). Hence, it follows that, Similarly, the commutator between the operator µ and the density matrix ρ(t) is given by,

B Vectorisation of the density operator
In this section we provide the vectorisation of the density operator introduced in Section (3). Although any vectorisation of the density operator can be used, we adopt the tetradic notation introduced in [62]. The density matrix ρ(t) described in Eq.(1) can be written as follows, By using the following vectorisation adopted from [62], the differential equations (4) can be equivalently written as, whereÃ ∈ C l 2 ×l 2 ,Ñ ∈ C l 2 ×l 2 are some matrices whose elements can be found from Eq.(4), andx 0 is the vectorisation of the initial density operator [63]. In addition T in Eq.(B.2) stands for the transpose operation.
Let x e be an eigenvector ofÃ with zero eigenvalue,Ãx e = 0 and setting x(t) =x(t) − x e and q = N x e , it follows that, whereB(x(t)) = N (x(t) + x e ) = N x(t) + q.
To give an example, using Eq.(4) the matricesÃ and N appearing in Eq.(B.4) for a two dimensional system, i. e. l = 0, 1 can be easily constructed. They are given by, The results can be generalised to any l > 1 in a straightforward manner.

C Complex normal distribution
For nonsingular covariance matrix Γ, the complex normal distribution for a complex random variable x t ∈ C n is given by, where |Γ| denotes the determinant of Γ.

D Calculation of the performance index γ(x t−1 )
In this section, we aim to calculate the form of the performance index given in Eq. (33). Let us first evaluate the coefficient β(u t−1 , x t−1 ) defined in Eq. (16) and repeated here, Let us first calculate ln s(xt|u t−1 ,x t−1 )s(ot|xt) Simplifying the above equation by collecting quadratic, linear and constant terms with respect to the control signal, u t−1 together, yields, Introduce the following definitions, It thus follows that, The first exponential function in Eq.(D.11) is nothing but the complex normal distribution recalled in Eq.(C.1) which when integrated over u t−1 yields a normalisation constant. Hence, Taking − ln of both sides of the above equation yields, This completes the proof of the Theorem (1).

E Distribution of the optimal control
In this section we provide the form of the distribution of the optimal control that minimises the Kullback-Leibler distance between (28) and (32). It is given in Eq. (14), repeated here, where v t and F t are defined in Eq.(D.10). This means that, F State space model for spin 1 2 In this section we show the forms of the matricesÃ andÑ appearing in Eq.(17) for a spin 1 2 interacting only with an external electric field. We assume that this interaction can be described by the following master equation, where as already stated in Eq.(53), H = 1 2 σ 3 + 1 2 (σ 1 + σ 2 )u(t), is defined in terms of the electric field u(t) and the Pauli matrices σ 1 , σ 2 , σ 3 are as given in Eq.(54), repeated here, The matrix representation of the Hamiltonian, H can be then easily evaluated. It is given by, By considering, −2ρ * 01 (t) + u(t)(1 + i)(ρ 00 (t) − ρ 11 (t)) u(t) (1 + i)ρ 01 (t) − (1 − i)ρ * 01 (t)   . (F.5) Using the vectorisation defined in Eq.(B.2), it follows that, yielding the state equation for the spin-1/2 system in the form given in Eq.(B.3), repeated here, dx(t) dt = (Ã + iÑ u(t))x(t).