# Probabilistic Pathways: New Frontiers in Quantum Ensemble Control

Randa Herzallah
Warwick Mathematics Institute, The University of Warwick, Coventry, CV4 7AL,
UK.

Abdessamad Belfakir
School of Physics and Astronomy, University of Nottingham, Nottingham, NG7
2RD, UK randa.herzallah<EMAIL_ADDRESS>
###### Abstract

In this paper, we propose a novel probabilistic control framework for
efficiently controlling an ensemble of quantum systems that can also
compensate for the interaction of the systems with the external environment.
The main challenge in this problem is to simultaneously steer an ensemble of
systems with variation in their internal parameters from an initial state to a
desired final state. The minimisation of the discrepancy between the
probabilistic description of the dynamics of a quantum ensemble and a
predefined desired probabilistic description is the key step in the proposed
framework. With this objective, the derived solution will not only allow the
transitioning of the ensemble from one state to another, but will generally
allow steering an initial distribution of the ensemble to a final
distribution. Numerical results are presented, demonstrating the effectiveness
of the proposed probabilistic control framework.

> Keywords: Probabilistic control, Kullback-Leibler Divergence (KLD), quantum
> control.

##  1 Introduction

Designing and implementing control solutions to simultaneously steer a large
quantum ensemble between states of interest is a critical aspect of many
advanced quantum technologies. These technologies include laser cooling,
nuclear magnetic resonance spectroscopy, magnetic resonance imaging, quantum
computation [1], and long-distance quantum communication [2]. Quantum
ensembles consist of a vast number of quantum systems (e.g., spin systems [3,
4, 5, 6]) with variations in their parameters. The primary challenge in
controlling ensembles is to derive a common external control field capable of
concurrently evolving each system in the ensemble from an initial state to a
desired final state, by adjusting global parameters of the overall system,
rather than controlling individual members.

Control strategies for quantum ensembles have been successfully applied to
molecular ensembles [7] and have proven effective in quantum information
processing [8]. Recent studies have focused on controlling inhomogeneous spin
ensembles and stabilizing spin system ensembles using Lyapunov control
methodology [9]. Ensemble controllability has been discussed in [10],
emphasizing the importance of Lie brackets and noncommutativity in designing
compensating controls. Further research has investigated the development of
unitary control in homogeneous quantum ensembles for maximizing signal
intensity in coherent spectroscopy [5] and employing the sampling-based
learning control (SLC) approach [11, 12].

Environmental decoherence, a crucial aspect of quantum mechanics, arises from
the unavoidable interaction between a system and its environment. While
extensively studied in controlling and manipulating quantum systems [13, 14,
15, 16, 17, 18], the issue of environmental decoherence has not yet been
adequately addressed within the context of quantum ensembles. Decoherence
poses a significant challenge for quantum control development, as it demands
precise control of system dynamics and results from the interaction between a
quantum system and its noisy environment [19, 20, 21]. Protecting quantum
system interferences forms the foundation of quantum information applications
[22]. Consequently, the ability to control and manage environmentally induced
decoherence would provide substantial benefits for controlling numerous
physical and chemical phenomena. Prior work on quantum optimal control in the
presence of decoherence has utilized adaptive feedback control methods [23]
and control optimization algorithms to maximize population transfer fidelity
or quantum gates. Genetic algorithms have been proposed to control and
suppress decoherence effects [24], and experimental control of decoherence in
molecular vibrational wave packets has been conducted in the laboratory [25].

Given the current state of the art in controlling quantum ensembles, there is
an urgent need to develop more accurate control methods that enable the
simultaneous steering of ensemble members to the same desired state in the
presence of uncertainty or field noise introduced by the system’s Hamiltonian
or its interaction with the environment. One highly efficient control method
for uncertain stochastic classical control systems is the fully probabilistic
control approach [26, 27, 28, 29]. This approach characterizes the generative
models of classical system dynamics using probabilistic descriptions. A
control law is then derived to minimize the discrepancy between the
probabilistic description of the joint distribution of the system dynamics and
its controller and a predefined desired joint distribution [27, 28]. This
method has been demonstrated on various classical systems and has proven
effective in deriving accurate control laws under high levels of uncertainty
and stochasticity [30, 31, 32]. Consequently, extending and adapting the fully
probabilistic control method to control quantum systems is an appealing
prospect.

In this paper, we focus on controlling ensembles of finite-dimensional quantum
systems that interact with external environments using the fully probabilistic
approach. Our goal is to develop a probabilistic control framework, building
upon its classical counterpart, capable of accurately steering the members of
a quantum ensemble to a desired state value while considering the interaction
of the quantum ensemble with the external environment in the optimal control
law design.

Extending the fully probabilistic control approach to the quantum mechanical
regime presents numerous challenges. A significant challenge is the coupling
of a quantum system with its external environment [33, 34, 35]. To address
this issue, we employ the Lindblad master equation formalism. This formalism
models the dynamics of open quantum systems, providing a framework that
incorporates the effects of environmental decoherence and dissipation into the
system’s evolution [33, 36]. Another hurdle is the complex nature of quantum
system state variables and operators. We tackle this by characterizing
generative probabilistic models based on complex probability density
functions, crucial for accurately capturing the unique features of quantum
systems and ensuring the effectiveness of our probabilistic control framework.
A key measure of success in our proposed control methodology is the ability to
guide each member of a quantum ensemble to a predetermined target state.
Achieving this necessitates the measurement of the final state of the quantum
systems in the ensemble and its comparison with the desired state. However,
measuring a quantum state is not as straightforward as in classical systems.
Quantum systems, defined by their wave functions or density operators, are
generally not directly observable. Instead, one must perform a series of
measurements on many identically prepared quantum systems to estimate their
state, a process known as quantum state tomography. This process aids in
estimating the system’s density matrix. Intricately tied to this measurement
process is the phenomenon of measurement backaction, another significant
challenge. This unique quantum trait, where the act of measurement can alter
the state of the system, could influence the post-control state of the system,
potentially affecting the control methodology’s performance. While intriguing,
we opt not to consider measurement backaction in this paper to maintain the
simplicity of our presentation of the proposed fully probabilistic approach.
However, it’s worth noting that a similar approach to environmental
interaction could be employed to account for measurement backaction [37, 38].

Consequently, our primary focus is on the mathematical formulation of fully
probabilistic control within quantum settings, as well as the analysis of a
class of bilinear quantum control systems. We aim to establish a rigorous
probabilistic control framework for quantum ensembles that takes into account
environmental decoherence, uncertainty, and the complex nature of quantum
state variables. This framework will facilitate the development of more robust
and accurate control methods for a variety of quantum technologies. By
presenting our approach in a simplified manner and systematically addressing
the key challenges, we provide a solid foundation for future work that could
extend our proposed approach to include the effects of measurement backaction
and other complex phenomena.

This paper is organised as follows: In Section (2) we briefly recall the
evolution of open quantum systems and develop the corresponding state space
model using the vectorization of their corresponding density operator. In
Section (3), we introduce a general theory to fully control quantum systems
affected by state dependent noisy environment using a probabilistic approach
and demonstrate its general solution. Then, we apply the developed approach in
Section (4) to systems described by Gaussian pdfs. In Section (5) we apply the
method to control an ensemble of particular systems and show the effectiveness
and the applicability of the method. Finally, we conclude in Section (6).

##  2 Problem formulation and Quantum systems models

There are numerous dynamical models describing quantum systems interacting
with external environments, depending on the type of system-environment
coupling [33, 34, 35]. An alternative approach to describe the time evolution
of the open quantum system is through master equations, the form of which
depends on the character of the system-environment interaction. Under the
Markovian assumption, an open quantum system can be described by the following
master equation:

$i\dfrac{d\rho(t)}{dt}=[H,\rho(t)]+\mathcal{L}(\rho(t)),\hskip
5.69046pt\rho(0)=\rho_{0},$ (1)

where $H$ is the Hamiltonian of the system and $\mathcal{L}[\rho(t)]$ is the
open system super-operator. In this paper, we set $\hbar=1$ and consider
Hamiltonians of the form:

$H=H_{0}+u(t)H_{1},$ (2)

where $H_{0}$ is the Hamiltonian of the isolated free system, and the term
$u(t)H_{1}$ is associated with the time-dependent control Hamiltonian defined
through the interaction of the system with the external field $u(t)$. In the
Lindblad approach, the open system Liouvillian is:

$\mathcal{L}(\rho(t))=i\sum_{s}\big{(}L_{s}\rho(t)L_{s}^{\dagger}-\dfrac{1}{2}\\{L_{s}^{\dagger}L_{s},\rho(t)\\}\big{)},$
(3)

where $L_{s}$ are Lindblad operators and $L_{s}^{\dagger}$ is the adjoint
operator of $L_{s}$. The time evolution of a physical property described by a
Hermitian operator $\hat{o}$ is given by:

${o}(t)=\braket{\hat{o}}=\operatorname{Tr}(\rho(t)\hat{o}).$ (4)

Given the definition of the master equation (1), we will now demonstrate that
the open quantum system defined by this master equation can be equivalently
described using a time-dependent equation that involves a state vector
representing the vectorization of the density operator. For this purpose, the
free Hamiltonian $H_{0}$ can be expanded element-wise as follows:

$H_{0}=\sum_{k=0}^{l-1}E_{k}\ket{k}\bra{k},$ (5)

where $E_{k}$ is the $k$th eigenvalue of $H_{0}$ associated with the
eigenvector $\ket{k}$, with $k\in{0,1,\dots,l-1}$. Similarly, the Lindblad
operators appearing in (3) can be written element-wise as follows:

$L_{s}=L_{j,k}=\sqrt{\Theta_{k\to j}}\ket{j}\bra{k},$ (6)

where $\Theta_{k\to j}$ is the dissipative transition rate from the eigenstate
$\ket{k}$ to the eigenstate $\ket{j}$. Following this element-wise
presentation, the master equation (1) can be equivalently written as [34]:

$\dfrac{d\rho_{n,q}(t)}{dt}=(-i{(E_{n}-E_{q})}-\gamma_{n,q})\rho_{n,q}(t)+\sum_{k=0}^{l-1}\Theta_{k\to
n}\rho_{k,k}(t)\delta_{n,q}+i{u(t)}\sum_{k=0}^{l-1}\bigg{(}\rho_{n,k}(t)\bra{k}H_{1}\ket{q}-\bra{n}H_{1}\ket{k}\rho_{k,q}(t)\bigg{)},$
(7)

where,

$\gamma_{n,q}:=\dfrac{1}{2}\sum_{j=0}^{l-1}(\Theta_{n\to j}+\Theta_{q\to j}),$
(8)

for $n,q={0,1\dots,l-1}$. Using (7), the vectorization $\tilde{x}$ of the
density operator $\rho(t)$ provided in Appendix (A) is shown to satisfy the
following differential equation:

$\dfrac{d\tilde{x}(t)}{dt}=(\tilde{A}+iu(t)\tilde{N})\tilde{x}(t),\hskip
14.22636pt\tilde{x}(0)=\tilde{x}_{0},$ (9)

where the matrix elements of the operators
$\tilde{A}\in\mathbf{C}^{l^{2}\times l^{2}}$ and
$\tilde{N}\in\mathbf{C}^{l^{2}\times l^{2}}$ can be easily found from (7). The
vectorization of the initial density operator $\rho_{0}$ is given by
$\tilde{x}_{0}$ [39, 40]. Considering the following transformation:

$x(t)\to\tilde{x}(t)-x_{e},$ (10)

the time evolution equation (9) can be rewritten as follows:

$\dfrac{dx({t})}{dt}=\tilde{A}x({t})+\tilde{B}(x(t))u(t),\hskip
14.22636ptx(0)=\tilde{x}_{0}-x_{e},$ (11)

where $x_{e}$ is an eigenvector of $\tilde{A}$ defined by $\tilde{A}x_{e}=0$
and $\tilde{B}(x(t))=i\tilde{N}(x(t)+x_{e})$. This equation, known as the
input equation, describes the dependence of the system’s dynamics, represented
by $x(t)$, on the input electric field ${u(t)}$ [40]. Thus, by solving the
ordinary differential equations (11), we can accurately determine the
evolution of the open quantum system. The solution to (11) can be expressed
as:

$\displaystyle x(t+1)=Ax(t)+B(x(t))u(t),$ (12)

where the matrices $A$ and $B(x({t}))$ are defined as:

$A=e^{\tilde{A}\Delta t},\hskip 14.22636pt\text{and}\hskip
14.22636ptB(x({t}))=\bigg{(}\int_{0}^{\Delta
t}e^{\tilde{A}\lambda}\tilde{B}(x(\lambda))\mathrm{d}\lambda\bigg{)},$ (13)

with $\lambda=\Delta t-t$ and $\Delta t$ representing the sampling period.
Equivalently, relation (12) can be rewritten as:

$\displaystyle x({t})=Ax({t-1})+B(x({t-1}))u(t-1).$ (14)

In practice, control of quantum physical systems can be achieved through the
use of multiple electric fields in the control Hamiltonian within the von
Neumann equation [11]. By denoting $x_{t}\equiv{x}(t)$ and $u(t)\equiv u_{t}$,
the relation (14) can be reformulated as:

$\displaystyle x_{t}=Ax_{t-1}+B(x_{t-1})u_{t-1}.$ (15)

In this work, we aim to control an ensemble of inhomogeneous quantum systems,
each characterized by a unique Hamiltonian due to the dispersion in the
parameters that describe them. This inhomogeneity or dispersion results in
each system responding slightly differently to the same control signal, a
challenge that needs to be addressed in control problems. To account for this
parameter dispersion and its impact on the state evolution of the systems in
the ensemble, we introduce multiplicative noise into our model. Unlike
additive noise, multiplicative noise varies with the state of the system,
making it a suitable choice for modeling uncertainties or variations that
scale with the system’s state. Therefore, we add the multiplicative noise
term, $\eta(x_{t-1})$, to the discretized equation (15). The process or
measurement noise levels in $\eta(x_{t-1})$ depend on the system state vector,
resulting in a state-dependent noise term. The final form of our model is
given by:

$\displaystyle x_{t}=Ax_{t-1}+B(x_{t-1})u_{t-1}+\eta(x_{t-1}).$ (16)

where,

$\eta(x_{t-1})=\zeta_{t}Ax_{t-1},$ (17)

and $\zeta_{t}$ is a scalar noise. This modified equation (16) incorporates
the effects of multiplicative noise and thus accommodates the parameter
dispersion inherent in our ensemble of inhomogeneous quantum systems. By
employing a vectorized description of the density matrix, we can now
reformulate the time evolution of the observable $\hat{o}$, given in (4), in
terms of the vectorized density operator. These adaptations to the model, both
in terms of multiplicative noise and the vectorized description of the density
matrix, allow us to better capture the unique dynamics of each member of the
quantum ensemble, enhancing the effectiveness of our control methodology. To
further improve the precision of our model, we consider additional sources of
noise and uncertainties that could affect the observed values $o_{t}$. This
consideration leads to the inclusion of the noise term $\sigma_{t}$ in our
output equation:

$o_{t}=Dx_{t}+\sigma_{t}.$ (18)

Here, $D=(\text{vec}(\hat{o}^{T}))^{T}$ represents the matrix obtained by
transposing the vectorized version of the observable operator $\hat{o}$, with
"vec" denoting the vectorization operator and $\hat{o}^{T}$ being the
transpose of $\hat{o}$. The term $\sigma_{t}$ is a multivariate Gaussian
noise, introduced to model various sources of noise and uncertainties, such as
measurement error, environmental disturbances, and process noise [38]. The
integration of this term helps enhance the realism and robustness of our
model, thereby improving its predictive accuracy and control performance.

By defining the quantum system’s dynamics through both the input (16) and
output (18) equations, we establish a bilinear state-space model that
effectively characterizes the behavior of the quantum systems in our ensemble.
This comprehensive model provides a strong foundation for the subsequent
design of effective control strategies.

A key innovation introduced in this paper lies in the extension of the
bilinear state-space model, originally formulated in previous works [38, 39,
40], to incorporate both multiplicative and additive Gaussian noise. By
transforming the quantum open system’s evolution, typically described by the
von Neumann-Lindblad equation, into this enhanced state-space model, we
provide a novel representation that not only offers a more comprehensive
picture of the quantum system’s dynamics but also explicitly accounts for
various sources of uncertainties such as manufacturing differences, different
environmental conditions, the inherent uncertainty in quantum mechanics,
sampling, parameter, and functional uncertainties. These are captured by the
inclusion of noise terms in both the input and output equations. This
innovative formulation serves as the foundation of the fully probabilistic
control framework for quantum systems proposed in this paper, setting the
stage for the introduction of advanced control strategies in the following
sections.

In the sections to follow, we will demonstrate how utilizing this bilinear
state-space model with noise terms enables the development of efficient
control strategies within a fully probabilistic framework. This approach
provides a robust characterization of the system’s time evolution,
specifically tailored for quantum systems, and allows for enhanced control and
management of their dynamics. Given the precision and robustness of this
methodology, it holds significant potential for a wide range of applications
in quantum information processing, quantum communication, and quantum sensing.
In these domains, precise control and estimation of quantum states are of
paramount importance [41, 42]..

##  3 Fully Probabilistic Control for Quantum Systems

Building upon the stochastic bilinear state-space model formulated in the
previous section, we now present the control objective and introduce the
general framework for the fully probabilistic control of ensembles of quantum
systems. This innovative framework is designed to account for various sources
of uncertainty and stochasticity affecting the quantum systems in a
constructive manner, ultimately yielding more accurate control laws under
these challenging conditions.

###  3.1 Objectives for the Fully Probabilistic Control of Quantum Systems

As previously discussed, this paper focuses on controlling ensembles of
quantum systems, where individual elements are subject to dispersion in their
parameters. These uncertainties are incorporated into the stochastic bilinear
state-space model represented by equations (16) and (18). Consequently, the
state of the system at each time step can be represented by the following
conditional probability density function (pdf):

$s(x_{t}\left|{x_{t-1}},u_{t-1}\right.).$ (19)

This representation provides a comprehensive description of the current state
of the system $x_{t}$, as influenced by the previous states of the system
$x_{t-1}$ and electric field $u_{t-1}$. Similarly, due to the presence of the
noise term $\sigma_{t}$ in equation (18), the state of the measurement $o_{t}$
can be characterized by the following pdf:

$\displaystyle s({o_{t}}\left|{{x_{t}}}\right.).$ (20)

Armed with these conditional pdfs, the objective of the fully probabilistic
control can be stated as designing a controller’s pdf, $c(u_{t}|x_{t})$, that
minimizes the Kullback-Leibler divergence (KLD) between the joint pdf of the
closed-loop system dynamics, $f(\mathcal{Z}(t,\mathcal{H}))$, and a pre-
specified ideal joint pdf, ${}^{I}f(\mathcal{Z}(t,\mathcal{H}))$:

$\mathcal{D}(f||^{I}f)=\int
f(\mathcal{Z}(t,\mathcal{H}))\ln\big{(}\dfrac{f(\mathcal{Z}(t,\mathcal{H}))}{{}^{I}f(\mathcal{Z}(t,\mathcal{H}))}\big{)}d\mathcal{Z}(t,\mathcal{H}),$
(21)

where

$f(\mathcal{Z}(t,\mathcal{H}))=\prod_{t=1}^{\mathcal{H}}s(x_{t}|x_{t-1},u_{t-1})s(o_{t}|x_{t})c(u_{t-1}|x_{t-1}),$
(22)

and

$^{I}f(\mathcal{Z}(t,\mathcal{H}))=\prod_{t=1}^{\mathcal{H}}{{}^{I}s(x_{t}|x_{t-1},u_{t-1})}{{}^{I}s(o_{t}|x_{t})}{{}^{I}c(u_{t-1}|x_{t-1})}.$
(23)

In these equations, ${}^{I}s(x_{t}|x_{t-1},u_{t-1})$ denotes the ideal pdf of
the state vector $x_{t}$, ${}^{I}s(o_{t}|x_{t})$ denotes the ideal pdf of the
measurement $o_{t}$,
$\mathcal{Z}(t,\mathcal{H})=\\{x_{t},\dots,x_{\mathcal{H}},o_{t},\dots,o_{\mathcal{H}},u_{t-1},\dots,u_{\mathcal{H}}\\}$
is the closed-loop observed data sequence and $\mathcal{H}<\infty$ is a given
control horizon, and ${}^{I}c(u_{t-1}|x_{t-1})$ represents the ideal pdf of
the controller. Given that the state $x_{t}$ is not directly observable and
that the output pdf is conditioned on the state $x_{t}$ of the quantum system,
it is reasonable to allow the state of the quantum system to evolve naturally.
Consequently, we assume that
${}^{I}s(x_{t}|x_{t-1},u_{t-1})=s(x_{t}|x_{t-1},u_{t-1})$. This assumption
simplifies the problem and focuses the controller design on accounting for
uncertainties in the measurements and control inputs, rather than attempting
to manipulate the natural evolution of the quantum system. By doing so, the
proposed fully probabilistic control framework remains compatible with the
inherent characteristics of quantum systems while effectively addressing the
challenges associated with uncertainties and parameter dispersion in ensembles
of quantum systems.

Previous works, such as those cited in [26, 27, 28], have demonstrated that
the minimization of the KLD equation (21) with respect to the control input is
achieved by initially defining $-\ln(\gamma(x_{t-1}))$ as the expected cost-
to-go function. In the context of this work, we express
$-\ln(\gamma(x_{t-1}))$ as follows:

$\displaystyle-\ln(\gamma(x_{t-1}))$
$\displaystyle=\underset{c(u_{t-1}|x_{t-1})}{\text{min}}\sum_{\tau=t}^{\mathcal{H}}\int
f(\mathcal{Z}_{t},\dots,\mathcal{Z}_{\mathcal{H}}|x_{t-1})\times\ln\left(\dfrac{s(o_{\tau}|x_{\tau})c(u_{\tau-1}|x_{\tau-1})}{{}^{I}s(o_{\tau}|x_{\tau})^{I}c(u_{\tau-1}|x_{\tau-1})}\right)d(\mathcal{Z}_{t},\dots,\mathcal{Z}_{\mathcal{H}}),$
(24)

for any arbitrary $\tau\in{1,...,\mathcal{H}}$. In this equation,
$\mathcal{Z}_{t}$ denotes the set ${x_{t},o_{t},u_{t-1}}$. Importantly,
observe that the argument of the $\ln$ function does not include the
conditional density of the quantum system state, $s(x_{t}|x_{t-1},u_{t-1})$,
nor its corresponding ideal distribution, ${}^{I}s(x_{t}|x_{t-1},u_{t-1})$.
This omission is not an oversight, but instead reflects a fundamental quantum
property: the state of the dynamic system, being unobservable, evolves
naturally according to the laws of quantum evolution. This quantum
characteristic is manifested by equating the actual distribution of the
quantum system state to its corresponding ideal distribution, i.e.,
${}^{I}s(x_{t}|x_{t-1},u_{t-1})=s(x_{t}|x_{t-1},u_{t-1})$. This equality
underpins the unobserved natural evolution of the quantum state, a principle
fundamental to quantum dynamics and a cornerstone of the proposed control
methodology.

Minimization of the cost to go function (24) can then be performed recursively
to give the following recurrence equation analogical to the dynamic
programming solution:

$\displaystyle-\ln(\gamma(x_{t-1}))$
$\displaystyle=\underset{c(u_{t-1}|x_{t-1})}{\text{min}}\int\bigg{[}s(x_{t}|x_{t-1},u_{t-1})s(o_{t}|x_{t})c(u_{t-1}|x_{t-1})\bigg{(}\ln\left(\dfrac{s(o_{t}|x_{t})c(u_{t-1}|x_{t-1})}{{}^{I}s(o_{t}|x_{t})^{I}c(u_{t-1}|x_{t-1})}\right)$
$\displaystyle-\ln(\gamma(x_{t}))\bigg{)}\bigg{]}d(x_{t},o_{t},u_{t}).$ (25)

This innovative approach for devising a controller pdf provides an elegant
solution for addressing the inherent uncertainties in quantum systems. By
minimizing the KLD between the actual and ideal joint pdfs, our method
facilitates more accurate control laws under various conditions of
uncertainty, thereby enhancing the overall performance and reliability of
quantum systems. The fully probabilistic control framework enables us to
account for various sources of uncertainty and stochasticity affecting the
quantum systems, and to construct efficient control strategies tailored
specifically for quantum systems. Moreover, the recursive nature of the
controller pdf design allows for efficient online implementation, making it
adaptable to real-time applications in quantum information processing, quantum
communication, and quantum sensing. The proposed approach addresses the
challenges associated with controlling ensembles of quantum systems and offers
a novel, comprehensive framework for the robust management of their dynamics.

###  3.2 General Solution to the Fully Probabilistic Quantum Control Problem

In this section, we derive the general solution of the optimal randomized
controller that minimizes the KLD given in equation (21) for systems defined
by arbitrary pdfs. This solution is not restricted to specific forms of the
pdfs and can be applied to any quantum system described by arbitrary
conditional pdfs.

###### Theorem 1

The pdf of the optimal control law $c(u_{t-1}|x_{t-1})$ that minimizes the
cost-to-go function (3.1) between the joint pdf of the closed loop description
of a quantum system (22) and a pre-specified ideal one (23) is given by:

$c(u_{t-1}|x_{t-1})=\frac{{}^{I}c(u_{t-1}|x_{t-1})\exp[-\beta(u_{t-1},x_{t-1})]}{\gamma(x_{t-1})},$
(26)

where,

$\gamma(x_{t-1})=\int{{}^{I}c(u_{t-1}|x_{t-1})\exp[{-\beta(u_{t-1},x_{t-1})}]\mathrm{d}u_{t-1}},$
(27)

and,

$\beta(u_{t-1},x_{t-1})=\int
s(x_{t}|u_{t-1},x_{t-1})s(o_{t}|x_{t})\times\ln\bigg{(}\frac{s(o_{t}|x_{t})}{{}^{I}s(o_{t}|x_{t})}\frac{1}{{\gamma}(x_{t})}\bigg{)}\mathrm{d}x_{t}\mathrm{d}o_{t}.$
(28)

###### Proof 1

This theorem can be easily proven by adapting the proof of Proposition 2 in
Ref [29].

The above theorem presents a general solution for the fully probabilistic
control problem of quantum systems, without being restricted to a specific
form of the generative probabilistic model describing the system dynamics.
This solution is applicable to any quantum system described by an arbitrary
pdf.

##  4 Solution to Quantum Control Problems with Gaussian pdfs

Building upon the general solution of the fully probabilistic control approach
presented in Section (3.2), we now focus on quantum systems characterized by
Gaussian probability density functions (pdfs). The linear Gaussian setting
simplifies the problem and provides an opportunity to demonstrate the efficacy
of the method.

###  4.1 System Description and Assumptions

This section considers quantum systems described by a bilinear state space
model of the form given in (16) and (18) when the noises $\eta_{t}$ in (17)
and $\sigma_{t}$ are generated using a Gaussian process. Under these
conditions, the pdf of the quantum system state (19) can be described by
linear Gaussian pdfs as follows:

$\displaystyle s(x_{t}\left|{x_{t-1}},u_{t-1}\right.)$
$\displaystyle\sim\mathcal{N}_{\mathcal{C}}(\mu_{t},{\Gamma_{t}}),$ (29)

where $\mu_{t}$ and $\Gamma_{t}$ are the mean and covariance matrices,
respectively:

$\displaystyle\mu_{t}$ $\displaystyle=E(x_{t})=Ax_{t-1}+Bu_{t-1},$ (30)
$\displaystyle\Gamma_{t}$
$\displaystyle=E((x_{t}-\mu_{t})(x_{t}-\mu_{t})^{\dagger}),$
$\displaystyle=E(\zeta_{t}Ax_{t-1}\zeta_{t}x_{t-1}^{\dagger}A^{\dagger}),$
$\displaystyle=Ax_{t-1}\Sigma x_{t-1}^{\dagger}A^{\dagger},$ (31)

where we used equations (16) and (17) in order to evaluate the covariance
matrix $\Gamma_{t}$. In this case, the matrices $A$ and $B$ represent the
state and control matrices, respectively, appearing in (16). The functional
$E(.)$ denotes the expectation value, $x_{t}^{\dagger}$ is the conjugate
transpose of $x_{t}$, $\Sigma=E(\zeta_{t}^{2})$ and
$\mathcal{N}_{\mathcal{C}}$ denotes a complex Gaussian distribution. The forms
of the complex Gaussian distribution and the normal Gaussian distribution are
recalled in Appendix (B).

Furthermore, with the noise term $\sigma_{t}$ in (18) generated by a Gaussian
distribution, the pdf of the measurement $o_{t}$ in (20) can be characterized
by Gaussian pdf as follows:

$\displaystyle s({o_{t}}\left|{{x_{t}}}\right.)$
$\displaystyle\sim\mathcal{N}(o_{d},G),$ (32)

where $o_{d}=Dx_{t}$ and
$G=E(\sigma_{t}\sigma_{t}^{T})=E((o_{t}-Dx_{t})(o_{t}-Dx_{t})^{T})$ are the
mean and covariance matrices, respectively of this pdf.

Under the aforementioned conditions, the quantum system state $x_{t}$,
measurement $o_{t}$, and the electric field $u_{t}$ at each time $t$ can be
fully described by the following joint pdf:

$f(x_{t},o_{t},u_{t-1}|x_{t-1})=c(u_{t-1}|x_{t-1})\mathcal{N}(Dx_{t},G)\mathcal{N}_{\mathcal{C}}(\mu{t},{\Gamma_{t}}).$
(33)

Taking this into account, we postulate that the ideal joint pdf governing the
entire system– including the system state, measurement, and controller– is
given by:

$^{I}f(x_{t},o_{t},u_{t-1}|x_{t-1})={}^{I}s(x_{t}|x_{t-1},{u_{t-1}}){}^{I}s({o_{t}}\left|{{x_{t}}}\right.){}^{I}c({u_{t-1}}\left|{x_{t-1}}\right),$
(34)

where the ideal distributions are given by:

$\displaystyle{}^{I}s(x_{t}|x_{t-1},{u_{t-1}})$
$\displaystyle=s(x_{t}|x_{t-1},{u_{t-1}})\sim\mathcal{N}_{\mathcal{C}}(\mu_{t},{\Gamma_{t}}),$
(35) $\displaystyle{}^{I}s({o_{t}}\left|{{x_{t}}}\right.)$
$\displaystyle\sim\mathcal{N}({o}_{d},G_{r}),$ (36)
$\displaystyle{}^{I}c({u_{t-1}}\left|{x_{t-1}}\right.)$
$\displaystyle\sim\mathcal{N}(u_{r},\Omega),$ (37)

and where $o_{d}$ and $G_{r}$ are the mean and covariance matrices of the
ideal pdf of the measurement respectively, and $u_{r}$ and $\Omega$ are the
mean and the covariance matrices respectively of the ideal pdf of the
controller. Furthermore, as illustrated in Equation (35), the ideal
distribution of the quantum state, ${}^{I}s(x_{t}|x_{t-1},{u_{t-1}})$, is set
to be identical to the actual distribution, $s(x_{t}|x_{t-1},{u_{t-1}})$,
which is characterized by the master equation (1). This equivalence signifies
the undetected natural evolution of the quantum state, an intrinsic
characteristic of quantum dynamics, which forms a pivotal component of the
proposed control methodology.

###  4.2 Solution to the linear Gaussian Quantum Control Problem

Given the defined joint pdf of the quantum system (33) and the ideal joint pdf
(34), we can now apply the general solution of the fully probabilistic control
approach presented in Section (3.2) to solve the quantum control problem for
linear Gaussian systems. Specifically, by incorporating the components of the
joint pdf of the closed-loop description of the quantum system (33) and the
ideal joint pdf (34) into equations (26)-(28), we can derive analytic
solutions for the optimal cost-to-go function and the randomized controller.
We will soon present the form of the pdf for the optimal control, but first,
let us outline the structure of the optimal cost-to-go function.

###### Theorem 2

By substituting the ideal distribution of the measurements (36), the ideal
distribution of the controller (37), the real distribution of the measurement
(32), and the real distribution of system dynamics (29) into (27), the optimal
cost to go function can be shown to be given by,

$\displaystyle-\ln\left({\gamma\left({{x_{t-1}}}\right)}\right)$
$\displaystyle=0.5{x}_{t-1}^{\dagger}{M_{t-1}}{x_{t-1}}+0.5P_{t-1}{{x}_{t-1}}+{0.5\omega_{t-1}},$
(38)

where,

$\displaystyle{M_{t-1}}$
$\displaystyle=A^{\dagger}(D^{\dagger}G_{r}^{-1}D+M_{t})A-A^{\dagger}(D^{\dagger}G_{r}^{-1}D+M_{t})^{\dagger}B\big{(}\Omega^{-1}+B^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t})B\big{)}^{-1}$
$\displaystyle
B^{\dagger}(D^{\dagger}G_{r}^{-1}D+M_{t})A+A^{\dagger}(D^{\dagger}G_{r}^{-1}D+M_{t}){\Sigma}A,$
(39)

$\displaystyle P_{t-1}$
$\displaystyle=(P_{t}-2{o}_{d}^{\dagger}G_{r}^{-1}D)A+2(\Omega^{-1}u_{r}-0.5B^{\dagger}(P_{t}^{\dagger}-2D^{\dagger}G_{r}^{-1}{o}_{d}))^{\dagger}\big{(}\Omega^{-1}+B^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t})B\big{)}^{-1}$
$\displaystyle B^{\dagger}(D^{\dagger}G_{r}^{-1}D+M_{t})A,$ (40)

and,

$\displaystyle\omega_{t-1}$
$\displaystyle={\omega_{t}}+{o}_{d}^{\dagger}G_{r}^{-1}{{o}_{d}}+\ln\bigg{(}\dfrac{|G_{r}|}{|G|}\bigg{)}-\operatorname{Tr}\big{(}G(G^{-1}-G_{r}^{-1})\big{)}+u_{r}^{\dagger}\Omega^{-1}u_{r}$
$\displaystyle-\big{(}\Omega^{-1}u_{r}-0.5B^{\dagger}(P_{t}^{\dagger}-2{D}^{\dagger}G_{r}^{-1}{{o}_{d}})\big{)}^{\dagger}\big{(}\Omega^{-1}+B^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t})B\big{)}^{-1}\big{(}\Omega^{-1}u_{r}-0.5B^{\dagger}(P_{t}^{\dagger}-2{D}^{\dagger}G_{r}^{-1}{{o}_{d}})\big{)}$
$\displaystyle-2\ln({|\Omega|^{-1/2}}|\Omega^{-1}+B^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t})B|^{-1/2}),$
(41)

where $|G|$ stands for the determinant of the matrix $G$.

###### Proof 2

The proof of this theorem is given in Appendix (C).

Having established the form of the optimal cost-to-go function (38), we can
utilize it in equation (26) to determine the distribution of the optimal
controller $c(u_{t-1}|x_{t-1})$. This optimizes the pdf of the entire system
(33), essentially bringing the pdf of the system as close as possible to the
ideal one. The form of this optimal controller is presented in the subsequent
theorem.

###### Theorem 3

The controller’s distribution, which minimizes the Kullback-Leibler Divergence
(KLD) between the joint pdf of the closed-loop system dynamics described in
(33) and the ideal joint pdf presented in (34), is given by:

$\displaystyle c({u_{t-1}}\left|{x_{t-1}}\right.)$
$\displaystyle\sim\mathcal{N}(v_{t-1},R_{t}),$ (42)

where,

$\displaystyle v_{t-1}$
$\displaystyle=\bigg{(}{\Omega^{-1}}+{B^{\dagger}}({D^{\dagger}}G_{r}^{-1}D+{M_{t}})B\bigg{)}^{-1}\bigg{(}{\Omega^{-1}}{u_{r}}-{B^{\dagger}}({D^{\dagger}}G_{r}^{-1}D+{M_{t}})A{x_{t-1}}-0.5{B^{\dagger}}(P_{t}^{\dagger}-2{D^{\dagger}}G_{r}^{-1}{{o}_{d}})\bigg{)},$
(43)

and,

$\displaystyle R_{t}$
$\displaystyle=\bigg{(}{\Omega^{-1}}+{B^{\dagger}}({D^{\dagger}}G_{r}^{-1}D+{M_{t}})B\bigg{)}^{-1}.$
(44)

###### Proof 3

The proof is given in Appendix (D).

The randomised optimal controller derived in Theorem (3) represents a novel
approach to controlling an ensemble of quantum systems. By introducing a
multiplicative noise, $\eta(x_{t-1})$ and an additive noise, $\sigma_{t}$, to
the state vector $x_{t}$ and the observations $o_{t}$ respectively, this
innovative linear form of the controller effectively addresses the challenges
of controlling randomly selected members of quantum ensembles. This
breakthrough in quantum control theory enhances the efficacy and adaptability
of control strategies for quantum systems.

###  4.3 Implementation of the Gaussian controller

The implementation of the Gaussian controller, as described in Theorem (3),
can be broken down into a step-by-step process, which is summarized in the
following algorithm. The implementation consists of two primary phases:
Optimization and Testing.

During the Optimization phase, the objective is to characterize the
probability density function (pdf) of the randomized controller by optimizing
its parameters. The optimization process utilizes historical data, simulated
data, or real-time data obtained online to achieve the desired performance of
the controlled quantum system. By optimizing these parameters, the randomized
controller will be able to drive the pdf of the joint distribution of the
closed-loop description of the quantum system to a pre-specified ideal joint
pdf.

The Testing phase involves using the Gaussian controller, constructed during
the Optimization phase, to drive members of an ensemble to the same target
state. The ensemble members are randomly generated, allowing for an assessment
of the controller’s performance under a range of initial conditions. This
benchmarking process demonstrates the robustness and efficiency of the
Gaussian controller when applied to real-world quantum systems.

Algorithm 1 Fully probabilistic control of quantum systems

1:Optimization:

2:Evaluate the operator $D$ associated with the target operator $\hat{o}$;

3:Compute the matrices $\tilde{A}$ and $\tilde{N}$ from equation (7);

4:Determine the predefined desired value ${o}_{d}$;

5:Specify the initial state $x_{0}$, the shifting state $x_{e}$, and then
calculate the initial value of the measurement state $o_{0}\leftarrow Dx_{0}$;

6:Provide the covariances $\Sigma,G,G_{r}$, and $\Omega$ of the multiplicative
noise, $\eta(x_{t-1})$, the measurement vector, $o_{t}$, the ideal
distribution of $o_{t}$, and the controller, $u_{t}$ respectively;

7:Initialize: $t\leftarrow 0$, $M_{0}\leftarrow\text{rand}$,
$P_{0}\leftarrow\text{rand}$;

8:while $t\neq\mathcal{H}$ do

9: Set: $\zeta_{t}=$ sample_from_Gaussian($0$, $\Sigma$);

10: Evaluate $A$ and $B$ using equation (13);

11: Calculate the steady state solutions of $M_{t}$ and $P_{t}$ following the
formulas provided in equations (2) and (2), respectively;

12: Use $M_{t}$ and $P_{t}$ to compute the optimal control input, $v_{t-1}$
following equation (43) given in Theorem (3);

13: Set: $u_{t-1}\leftarrow v_{t-1}$;

14: Use the obtained control input from the previous step to evaluate $x_{t}$
according to equation (16), $x_{t}\leftarrow
Ax_{t-1}+B(x_{t-1})u_{t-1}+\zeta_{t}Ax_{t-1}$;

15: Set: $\sigma_{t}=$ sample_from_Gaussian($0$, $G$);

16: Following equation (18), evaluate $o_{t}$ to find the measurement state at
time instant $t$, $o_{t}\leftarrow Dx_{t}+\sigma_{t}$;

17: Update time: $t\leftarrow t+1$;

18:end while

19:Testing:

20:Set $N$: the number of members of the ensemble;

21:Initialize: $k\leftarrow 1$;

22:while $k\leq N$ do

23: while $t\neq\mathcal{H}$ do

24: Set: $\zeta^{k}_{t}=$ sample_from_Gaussian($0$, $\Sigma$);

25: Evaluate $A^{k}$ and $B^{k}$ using equation (13);

26: Using the obtained control input from the optimization phase, evaluate
$x_{t}^{k}$ according to equation (16), $x_{t}^{k}\leftarrow
A^{k}x^{k}_{t-1}+B^{k}(x^{k}_{t-1})u_{t-1}+\zeta_{t}^{k}A^{k}{x}_{t-1}$;

27: Set: $\sigma^{k}_{t}=$ sample_from_Gaussian($0$, $G$);

28: Following equation (18), evaluate $o^{k}_{t}$ to find the measurement
state at time instant $t$, $o^{k}_{t}\leftarrow Dx_{t}^{k}+\sigma_{t}^{k}$;

29: end while

30: Increment ensemble member: $k\leftarrow k+1$;

31:end while

Algorithm (1) will be applied in the next section to control an ensemble of
randomly generated systems to the same desired state.

##  5 Results and discussions

###  5.1 Control of quantum ensemble

In this section, we demonstrate the effectiveness of our algorithm when
applied to an ensemble of quantum systems. We consider a simplified
Hamiltonian $H$ as described in (2), given by:

$H=H_{0}+{u}_{t}H_{1},$ (45)

where we assume that the system interacts with only one electric field,
represented by ${u}_{t}$. Our control objective is to transfer an ensemble of
quantum systems, initially prepared in a state $\rho_{i}$, to a predefined
target state $\rho_{d}$. This transfer is achieved through the interaction
between the ensemble members and the optimal randomized electric field. The
parameters of this electric field are optimized using Theorem (3) to achieve a
pre-specified desired joint density function for the joint distribution of the
ensemble dynamics. We will demonstrate that the optimized randomized
controller is highly effective in controlling all members of the ensemble.

The control problem we aim to address is to maximize the fidelity between the
actual state $\rho_{t}$ and the target state $\rho_{d}$ for all ensemble
members. Here, $\rho_{t}$ represents the density operator describing the
systems at time $t$. By applying our algorithm (1) to this ensemble of quantum
systems, we showcase its efficiency and versatility in controlling complex
quantum ensembles. This highlights the potential of our algorithm to tackle
real-world quantum control challenges.

In the following, we present the simulation results for the control of the
quantum ensemble using our proposed Gaussian controller. The simulations
demonstrate the robustness and effectiveness of the controller in driving the
ensemble members towards the desired target state. Furthermore, we will
discuss the implications of our results and the potential applications of our
algorithm in quantum control and quantum information processing tasks.

####  5.1.1 spin-1/2 ensemble

Consider a spin-$1/2$ system interacting with an electric field $u_{t}$ [11].
The Hamiltonian describing this interaction is given by:

$H=H_{0}+{u}_{t}H_{1}=\dfrac{1}{2}\sigma_{3}+\dfrac{{u}_{t}}{2}(\sigma_{1}+\sigma_{2}),$
(46)

where $\sigma_{1}$, $\sigma_{2}$, and $\sigma_{3}$ are the Pauli operators
given in the basis ${\ket{1},\ket{0}}$ by the Pauli matrices,

$\sigma_{1}=\left(\begin{array}[]{cc}0&1\\\ 1&0\end{array}\right),\hskip
11.38092pt\sigma_{2}=\left(\begin{array}[]{cc}0&-i\\\
i&0\end{array}\right),\hskip
11.38092pt\sigma_{3}=\left(\begin{array}[]{cc}1&0\\\ 0&-1\end{array}\right).$
(47)

Furthermore, we consider that the system interacts with an environment, e.g.,
a surrounding vacuum state. The master equation describing this interaction is
given by:

$\dfrac{d\rho(t)}{dt}=-i[H,\rho(t)]+\Theta\bigg{(}\sigma_{-}\rho(t)\sigma_{+}-\dfrac{1}{2}{\sigma_{+}\sigma_{-},\rho(t)}\bigg{)},\hskip
5.69046pt\rho(0)=\rho_{0}=\ket{0}\bra{0},$ (48)

where $\sigma_{\pm}=\dfrac{\sigma_{1}\pm i\sigma_{2}}{2}$ and $\Theta$
describes the coupling between the system and the environment.

In the optimization phase, the control objective is to transform the system
state from its initial value $x_{0}=(0,1,0,0)^{T}$ which corresponds to the
first eigenstate of the free Hamiltonian, $\ket{0}=(0,1)^{T}$, to the desired
state $x_{d}=(1,0,0,0)^{T}$ which corresponds to the second eigenstate of the
free Hamiltonian $\ket{1}=(1,0)^{T}$, using the proposed probabilistic quantum
control approach. Consequently, the target operator is taken to be $D=[1\hskip
5.69046pt0\hskip 5.69046pt0\hskip 5.69046pt0]$, which is equivalent to
$D=(\text{vec}(\Pi_{1}))^{T}$, where $\Pi_{1}=\ket{1}\bra{1}$. The matrix
elements of the $\tilde{A}$ and $\tilde{N}$ matrices are provided in Appendix
(E). In the optimization phase of Algorithm (1), we consider
$x_{e}=(0,0,0,0)^{T}$, $\Delta{t}=2.5\times 10^{-6}$, $G_{r}=0.00001$,
$\Omega=10$, and $\Theta=0.1$, and the system is evolved until a sufficient
fidelity between the actual state and the target state is achieved.

Following the optimization phase, the obtained control signal is used to
control one of the members of the spin-1/2 system. Figure (1(a)) illustrates
the obtained time evolution of the populations of the two states
$\ket{0}\bra{0}$ and $\ket{1}\bra{1}$ of the selected member, and Figure
(1(b)) displays the time evolution of the optimal control signal responsible
for achieving the control objective. We can clearly see that the population of
the state $\ket{1}\bra{1}$ has reached its desired value, i.e., $o_{d}=1$, and
it is maintained at that value for extended periods. This indicates that the
optimal electric field is efficient in preventing the environmentally induced
relaxation of the system to the ground state $\ket{0}\bra{0}$, signifying that
the optimal control can counteract the natural effect of the environment,
i.e., relaxation.

Furthermore, the optimized randomized controller obtained from the
optimization phase is applied to a sample of ensemble members, each
corresponding to different realizations of the state time evolution generated
through the noise term $\eta(x_{t-1})$. This experiment allows us to assess
the performance of the controller on a diverse set of members, each
representing a unique noise-driven evolution. Figure (2) displays the
fidelities between the target state $\ket{1}\bra{1}$ and the states of $1000$
ensemble members interacting with the optimal control electric field. From
this figure, we can clearly see that the fidelities for the state transition
of all tested 1000 members lie in the interval of $[0.9932,1]$ with a mean
value equal to $0.9945$. This demonstrates that the algorithm exhibits
exceptional performance in this specific example. It is worth mentioning that
we have considered members which are interacting with their external
environment. This implies that the optimal electric field is capable of
driving the ensemble members to the desired state and preventing their
relaxation, even in the presence of environmental effects. The results
underscore the robustness and effectiveness of our algorithm in addressing
real-world quantum control challenges that involve noise and environmental
interactions.

(a)

(b)

Figure 1: (1(a) ): Time evolution of the populations of the $\ket{0}\bra{0}$
and $\ket{1}\bra{1}$ states. (1(b)): Time evolution of the control signal
$u_{t}$ responsible of achieving the control objective. Figure 2:  The testing
performance of the optimal control for spin-$1/2$ ensemble. The samples are
randomly generated through the noise term $\eta(x_{t-1})$.

####  5.1.2 $\Lambda$-type atomic ensemle

(a)

(b)

Figure 3: (3(a) ): Time evolution of the populations the $\ket{0}\bra{0}$,
$\ket{1}\bra{1}$ and $\ket{2}\bra{2}$ states. (3(b)): Time evolution of the
control signal, $u_{t}$ responsible of achieving the control objective. Figure
4:  The testing performance of the optimal control for three level systems
ensemble. The samples are randomly generated through the noise term
$\eta(x_{t-1})$.

We now investigate a $\Lambda$-type atomic system [11] prepared in an initial
state and interacts with an electric field ${u}_{t}$. The Hamiltonian
describing this interaction is represented by matrices in the basis
${\ket{2},\ket{1},\ket{0}}$ as follows:

$H=H_{0}+{u}_{t}H_{1}=\left(\begin{array}[]{ccc}\dfrac{3}{2}&0&0\\\ 0&1&0\\\
0&0&0\\\ \end{array}\right)+\left(\begin{array}[]{ccc}0&0&1\\\ 0&0&1\\\
1&1&0\\\ \end{array}\right){u}_{t}.$ (49)

Moreover, we assume that the system interacts with an environment, with the
system-environment coupling being described by a single Lindblad operator,
$L_{02}=\sqrt{\Theta}\ket{0}\bra{2}$, where $\Theta\equiv\Theta_{2\to 0}$
represents the dissipative transition rate from Hamiltonian eigenstate
$\ket{2}$ to the Hamiltonian ground state $\ket{0}$.

The control objective aims to transfer the system from state $\ket{0}\bra{0}$
to the desired state $\rho_{d}=\ket{\psi_{d}}\bra{\psi_{d}}$, with
$\ket{\psi_{d}}=\dfrac{1}{\sqrt{2}}(\ket{1}+\ket{2})$. To achieve this, we
must maximize the fidelity between the actual state and the target state
$\rho_{d}$. The target operator $D$ in (18) now takes the form
$D=(\text{vec}(\Pi_{d})^{T})^{T}$, where $\Pi_{d}=\ket{\psi_{d}}\hskip
1.42271pt\bra{\psi_{d}}$ represents the projector onto the target state
$\ket{\psi_{d}}$. By focusing on this control objective, we can ensure a high
degree of precision in the transfer of the system between states, even when
accounting for the environmental interactions. This approach highlights the
versatility of the probabilistic quantum control method in effectively
managing different quantum systems and control objectives.

The elements of the $\tilde{A}$ and $\tilde{N}$ matrices are provided in
Appendix (F). These matrices are then employed to compute the time-evolution
of matrices $A$ and $B$ using (13) at each time step. With these matrices and
setting $x_{e}=(0,0,0,0,0,0,0,0,0)^{T}$, $\Delta{t}=0.0000000001$,
$G_{r}=0.0000001$, $\Omega=10000$, and $\Theta=0.9$, we evaluate the matrices
$M_{t}$ and $P_{t}$ defined in equations (2) and (2), respectively, at each
instant of time, as discussed in the optimization phase of Algorithm (1). The
evolution of $u_{t}$ is then obtained by substituting the calculated $M_{t}$
and $P_{t}$ matrices into (43). We repeat these steps until the fidelity
between the actual state and the target state reaches its maximum value of
unity, signifying that the system’s state has reached the desired target
state.

Following the optimization phase, the obtained optimal randomized controller
is used to control one of the members of the $\Lambda-$type system. Figure
(3(a)) illustrates the time evolution of the populations of states
$\ket{0}\bra{0}$, $\ket{1}\bra{1}$, and $\ket{2}\bra{2}$ of the selected
ensemble interacting with the calculated electric field. The time evolution of
the obtained optimal control signal is depicted in Figure (3(b)). It is
evident that the optimal control objective is achieved after only a few steps.
As a result, the optimal control is capable of preserving quantum coherence,
i.e., the quantum superposition, under the influence of decoherence.

To further demonstrate the capability of the designed randomised controller to
effectively control the members of the $\Lambda-$type system, the optimal
controller is applied to other randomly selected ensemble members generated
through the noise term $\eta(x_{t-1})$. These ensemble members represent
different realizations of the state time evolution influenced by the noise
term $\eta(x_{t-1})$. Figure (4) displays the fidelities between the target
state and the final states of $1000$ ensemble members interacting with the
optimal controller. We observe that the fidelities are almost equal to 1,
implying that all members successfully reached the target state
$\ket{\psi_{d}}$. This result highlights the potential effectiveness of the
proposed method in controlling $\Lambda$-type atomic systems under
environmental effects.

##  6 Final comments

In this paper, we have presented a novel approach for deriving optimal
randomized controllers for quantum ensembles. Our proposed method consists of
two main phases: optimization and testing. The optimization phase involves
constructing a control for a generalized system of the ensemble using a
probabilistic approach, which is based on describing the system dynamics with
probability density functions and minimizing the distance between the actual
closed-loop joint pdf and the desired joint pdf. We first established a
general control solution for quantum systems described by arbitrary pdfs. The
solution is then demonstrated on a specific case involving quantum systems
affected by Gaussian noise and described by Gaussian pdfs, which yielded an
analytical form for the randomized controller.

In the testing phase, the controller obtained during the optimization phase is
applied to a large number of randomly selected ensemble members. Our results
demonstrated that, with the optimized controller, all members of the quantum
ensemble that were also affected by environmental noises and interactions were
successfully steered to the same target state, illustrating the effectiveness
of our proposed methodology. In practical applications, the verification of
successful steering to the target state would involve performing quantum state
tomography on the final state of the system, and calculating the fidelity
between the actual state and the desired state. Moreover, it’s important to
note a crucial aspect of quantum mechanics: the act of measurement itself can
influence the state of the quantum system, a phenomenon known as measurement
backaction. While this effect was not explicitly accounted for in the current
work, it’s an important consideration as it can potentially impact the
effectiveness of the control field. Future work could explore integrating the
effect of measurement backaction into the control methodology to ensure robust
control even in the presence of frequent measurements.

In practical settings, the approach entails deriving the control field using
historical, simulated, or real-time data obtained from quantum systems with
the introduced method, and subsequently applying the electric field to the
quantum ensemble.

The main advantages of the proposed methodology include:

  * •

It provides a general form for the time evolution of the controller for
systems affected by arbitrary noise, making it suitable for a wide range of
quantum systems.

  * •

It is easy to implement numerically, as explained in Algorithm (1), and does
not require significant computational resources.

  * •

The method is highly adaptable and can be applied to various quantum systems
and control objectives, making it appealing to physicists and experimentalists
working with quantum systems.

Our approach allows for the efficient control and manipulation of quantum
systems, even in the presence of environmental noises and interactions. This
can be particularly useful for preserving quantum coherence and achieving
specific control objectives in various experimental settings.

Possible extensions of this work include considering quantum systems with
arbitrary forms of noise affecting their evolutions. This would further expand
the applicability and versatility of our proposed method, enabling it to
address a broader range of quantum systems and control challenges.

## Acknowledgements:

This work was supported by the EPSRC grant EP/V048074/1.

## Appendices

##  Appendix A Vectorization of the density operator

In equation (1), we consider a density matrix $\rho(t)\in\mathbf{C}^{l\times
l}$ given by:

$\rho(t)=(\rho(t))^{\dagger}=\left(\begin{array}[]{ccccccccc}\rho_{0,0}(t)&\rho_{0,1}(t)&\dots&\rho_{0,l-1}(t)\\\
\rho_{1,0}(t)&\rho_{1,1}(t)&\dots&\rho_{1,l-1}(t)\\\
\vdots&\vdots&\ddots&\vdots\\\
\rho_{l-1,0}(t)&\rho_{l-1,1}(t)&\dots&\rho_{l-1,l-1}(t)\\\
\end{array}\right).$ (A.1)

The vectorization of $\rho(t)$ used in (9) is:

$\displaystyle\tilde{x}(t)$ $\displaystyle=\text{vec}(\rho(t))$
$\displaystyle=\begin{bmatrix}\rho_{0,0}(t)&\rho_{1,1}(t)\dots\rho_{l-1,l-1}(t)&\rho_{0,1}(t)\dots\rho_{0,l-1}(t)&\rho_{1,0}(t)\dots\rho_{l-1,0}(t)\dots\dots\dots\rho_{l-2,l-1}(t)&\rho_{l-1,l-2}(t)\end{bmatrix}^{T},$
(A.2)

where $T$ stands for the transpose operation. Thus, for $l=3$ the
vectorization of the density matrix
$\rho(t)=\left(\begin{array}[]{ccc}\rho_{00}(t)&\rho_{01}(t)&\rho_{02}(t)\\\
\rho_{10}(t)&\rho_{11}(t)&\rho_{12}(t)\\\
\rho_{20}(t)&\rho_{21}(t)&\rho_{22}(t)\end{array}\right)$ is trivially given
by:

$\tilde{x}(t)=\tilde{x}_{t}=\begin{bmatrix}\rho_{00}(t)&\rho_{11}(t)&\rho_{22}(t)&\rho_{01}(t)&\rho_{02}(t)&\rho_{10}(t)&\rho_{20}(t)&\rho_{12}(t)&\rho_{21}(t)\end{bmatrix}^{T}.$
(A.3)

##  Appendix B Complex Gaussian and Gaussian distributions

For nonsingular covariance matrix $\Gamma_{t}$, the complex Gaussian
distribution of a complex random variable $x_{t}\in\mathbf{C}^{n}$ is given
by:

$\displaystyle\mathcal{N}_{\mathcal{C}}(\mu_{t},{\Gamma_{t}})=\dfrac{1}{\pi^{n}|\Gamma_{t}|}\exp\bigg{[}-{({{x}_{t}}-{{\mu}_{t}})^{\dagger}}{\Gamma_{t}}^{-1}({x_{t}}-{\mu_{t}})\bigg{]},$
(B.1)

where $\mu_{t}=E(x_{t})$, and
$\Gamma_{t}=\text{E}((x_{t}-\mu_{t})(x_{t}-\mu_{t})^{\dagger})$, with
$\text{E}(x_{t})$ being the expected value of $x_{t}$, $|\Gamma_{t}|$ is the
determinant of $\Gamma_{t}$, and $\dagger$ stands for the complex conjugate
transpose operation.
The Gaussian distribution for a real random variable $o_{t}\in\mathbf{R}^{n}$
is given by:

$\displaystyle\mathcal{N}(o_{m},{G})$
$\displaystyle=\dfrac{1}{\sqrt{(2\pi)^{n}|G|}}\exp\bigg{[}-{0.5({{o}_{t}}-{{o}_{m}})^{T}}{G}^{-1}({o_{t}}-{o_{m}})\bigg{]}.$
(B.2)

where $o_{m}=E(o_{t})$ and $G=\text{E}((o_{t}-o_{m})(o_{t}-o_{m})^{T})$ are
respectively the mean and covariance matrices, and $T$ stands for the
transpose operation. Since $o_{t}$ belongs to the set of real numbers
$\mathbf{R}^{n}$, the Gaussian distribution (B.2) can be equivalently
restructured as:

$\displaystyle\mathcal{N}(o_{m},{G})$
$\displaystyle=\dfrac{1}{\sqrt{(2\pi)^{n}|G|}}\exp\bigg{[}-{0.5({{o}_{t}}-{{o}_{m}})^{\dagger}}{G}^{-1}({o_{t}}-{o_{m}})\bigg{]}.$
(B.3)

This recasting of the real Gaussian distribution into the form of a complex
Gaussian distribution provides a consistency in the derivation of the
randomized controller. It not only simplifies the presentation but also aligns
well with the inherent nature of quantum systems, where the state variable is
complex. Thus, the adapted representation facilitates a more coherent
theoretical framework for our probabilistic quantum control approach.

##  Appendix C Evaluation of the optimal cost-to-go, $\gamma(x_{t-1})$

In this appendix we show how to obtain the optimal cost-to-go function (38)
stated in Theorem (2). For this purpose, we first evaluate the coefficient
$\beta(u_{t-1},x_{t-1})$ defined in equation (28), repeated here:

$\displaystyle\beta(u_{t-1},x_{t-1})=\int
s(x_{t}|u_{t-1},x_{t-1})s(o_{t}|x_{t})\times\ln\bigg{(}\frac{s(o_{t}|x_{t})}{{}^{I}s(o_{t}|x_{t})}\dfrac{1}{{\gamma}(x_{t})}\bigg{)}\mathrm{d}x_{t}\mathrm{d}o_{t}.$
(C.1)

Using equations (32), (36) and (38) we evaluate:

$\displaystyle\ln\bigg{(}\frac{s(o_{t}|x_{t})}{{}^{I}s(o_{t}|x_{t})}\dfrac{1}{{\gamma}(x_{t})}\bigg{)}$
$\displaystyle=-0.5{({o_{t}}-D{x_{t}})^{\dagger}}{G^{-1}}({o_{t}}-D{x_{t}})+0.5{({o_{t}}-{{o}_{d}})^{\dagger}}G_{r}^{-1}({o_{t}}-{{o}_{d}})+0.5x_{t}^{\dagger}{M_{t}}{x_{t}}+0.5P_{t}{x_{t}}+0.5{\omega_{t}}+0.5\ln\bigg{(}\dfrac{|G_{r}|}{|G|}\bigg{)},$
(C.2)

which implies that:

$\displaystyle\ln\bigg{(}\frac{s(o_{t}|x_{t})}{{}^{I}s(o_{t}|x_{t})}\dfrac{1}{{\gamma}(x_{t})}\bigg{)}$
$\displaystyle=-0.5o_{t}^{\dagger}{G^{-1}}{o_{t}}-0.5x_{t}^{\dagger}{D^{\dagger}}{G^{-1}}D{x_{t}}+o_{t}^{\dagger}{G^{-1}}D{x_{t}}+0.5o_{t}^{\dagger}G_{r}^{-1}{o_{t}}-o_{t}^{\dagger}G_{r}^{-1}{{o}_{d}}+0.5x_{t}^{\dagger}{M_{t}}{x_{t}}+0.5P_{t}{x_{t}}+0.5{\omega_{t}}$
$\displaystyle+0.5{o}_{d}^{\dagger}G_{r}^{-1}{{o}_{d}}+0.5\ln\bigg{(}\dfrac{|G_{r}|}{|G|}\bigg{)}$
$\displaystyle=-0.5o_{t}^{\dagger}({G^{-1}}-G_{r}^{-1}){o_{t}}+o_{t}^{\dagger}({G^{-1}}D{x_{t}}-G_{r}^{-1}{{o}_{d}})-0.5x_{t}^{\dagger}D^{\dagger}{G^{-1}}D{x_{t}}+0.5x_{t}^{\dagger}{M_{t}}{x_{t}}+0.5P_{t}{x_{t}}+0.5{\omega_{t}}$
$\displaystyle+0.5{o}_{d}^{\dagger}G_{r}^{-1}{{o}_{d}}+0.5\ln\bigg{(}\dfrac{|G_{r}|}{|G|}\bigg{)}.$
(C.3)

By substituting (C) into (C.1), and integrating over $o_{t}$ we get:

$\displaystyle\beta(u_{t-1},x_{t-1})=\int
s(x_{t}|u_{t-1},x_{t-1})s(o_{t}|x_{t})\bigg{(}-0.5o_{t}^{\dagger}({G^{-1}}-G_{r}^{-1}){o_{t}}+o_{t}^{\dagger}({G^{-1}}D{x_{t}}-G_{r}^{-1}{{o}_{d}})-0.5x_{t}^{\dagger}D^{\dagger}{G^{-1}}D{x_{t}}$
$\displaystyle+0.5x_{t}^{\dagger}{M_{t}}{x_{t}}+0.5P_{t}{x_{t}}+0.5{\omega_{t}}+0.5{o}_{d}^{\dagger}G_{r}^{-1}{{o}_{d}}+0.5\ln\bigg{(}\dfrac{|G_{r}|}{|G|})\bigg{)}\mathrm{d}x_{t}\mathrm{d}o_{t}$
$\displaystyle=\int{s\left({{x_{t}}|{u_{t-1}},{\rm{}}{x_{t-1}}}\right){\rm{}}}\bigg{(}0.5x_{t}^{\dagger}{D^{\dagger}}G_{r}^{-1}D{x_{t}}-x_{t}^{\dagger}{D^{\dagger}}G_{r}^{-1}{{o}_{d}}+0.5x_{t}^{\dagger}{M_{t}}{x_{t}}+0.5P_{t}{x_{t}}+0.5{\omega_{t}+}0.5{o}_{d}^{\dagger}G_{r}^{-1}{{o}_{d}}$
$\displaystyle-0.5\operatorname{Tr}((G^{-1}-G_{r}^{-1})G)+0.5\ln\bigg{(}\dfrac{|G_{r}|}{|G|}\bigg{)}\bigg{)}d{x_{t}}$
$\displaystyle=\int{s\left({{x_{t}}|{u_{t-1}},{\rm{}}{x_{t-1}}}\right){\rm{}}}\bigg{(}0.5x_{t}^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t}){x_{t}}+0.5(P_{t}-2{o}_{d}^{\dagger}G_{r}^{-1}{D})x_{t}+0.5{\omega_{t}}+0.5{o}_{d}^{\dagger}G_{r}^{-1}{{o}_{d}}+0.5\ln\bigg{(}\dfrac{|G_{r}|}{|G|}\bigg{)}$
$\displaystyle-0.5\operatorname{Tr}((G^{-1}-G_{r}^{-1})G)\bigg{)}d{x_{t}}.$
(C.4)

Integrating over $x_{t}$ yields the following form of
$\beta(u_{t-1},x_{t-1})$:

$\displaystyle\beta(u_{t-1},x_{t-1})=0.5\mu_{t}^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t}){\mu_{t}}+0.5(P_{t}-2{o}_{d}^{\dagger}G_{r}^{-1}{D})\mu_{t}+0.5{\omega_{t}}+0.5{o}_{d}^{\dagger}G_{r}^{-1}{{o}_{d}}+0.5\ln\bigg{(}\dfrac{|G_{r}|}{|G|}\bigg{)}$
$\displaystyle-0.5\operatorname{Tr}((G^{-1}-G_{r}^{-1})G)+0.5\operatorname{Tr}(({D^{\dagger}}G_{r}^{-1}D+M_{t})\Gamma_{t}).$
(C.5)

Using the definition of $\Gamma_{t}$ given in equation (31), the last term in
the above equation can be evaluated to give:

$\displaystyle\operatorname{Tr}(({D^{\dagger}}G_{r}^{-1}D+M_{t})\Gamma_{t})$
$\displaystyle=\operatorname{Tr}(({D^{\dagger}}G_{r}^{-1}D+M_{t})A{x}_{t-1}\Sigma
x_{t-1}^{\dagger}{A}^{\dagger}),$ $\displaystyle=<\text{cyclic permutation}>$
$\displaystyle\operatorname{Tr}({A}^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t})A{x}_{t-1}\Sigma
x_{t-1}^{\dagger}),$ $\displaystyle=<\text{identity
$\operatorname{Tr}(Qzz^{\dagger})=z^{\dagger}Qz$}>$ $\displaystyle
x_{t-1}^{\dagger}{A}^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t})A{x}_{t-1}\Sigma.$
(C.6)

Remembering that $\mu_{t}=Ax_{t-1}+Bu_{t-1}$, the relation (C) becomes:

$\displaystyle\beta(u_{t-1},x_{t-1})=0.5(Ax_{t-1}+Bu_{t-1})^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t}){(Ax_{t-1}+Bu_{t-1})}+0.5(P_{t}-2{{o}_{d}^{\dagger}}G_{r}^{-1}D)(Ax_{t-1}+Bu_{t-1})$
$\displaystyle+0.5{\omega_{t}}+0.5{o}_{d}^{\dagger}G_{r}^{-1}{{o}_{d}}+0.5\ln\bigg{(}\dfrac{|G_{r}|}{|G|}\bigg{)}-0.5\operatorname{Tr}((G^{-1}-G_{r}^{-1})G)+0.5x_{t-1}^{\dagger}{A}^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t}){\Sigma}A{x}_{t-1}$
$\displaystyle=0.5x_{t-1}^{\dagger}A^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t}){Ax_{t-1}}+0.5(P_{t}-2{{o}_{d}^{\dagger}}G_{r}^{-1}D)Ax_{t-1}+0.5u_{t-1}^{\dagger}B^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t}){Bu_{t-1}}$
$\displaystyle+u_{t-1}^{\dagger}B^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t}){Ax_{t-1}}+0.5(P_{t}-2{{o}_{d}^{\dagger}}G_{r}^{-1}D)Bu_{t-1}+0.5{\omega_{t}}+0.5{o}_{d}^{\dagger}G_{r}^{-1}{{o}_{d}}+0.5\ln\bigg{(}\dfrac{|G_{r}|}{|G|}\bigg{)}$
$\displaystyle-0.5\operatorname{Tr}((G^{-1}-G_{r}^{-1})G)+0.5x_{t-1}^{\dagger}{A}^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t}){\Sigma}A{x}_{t-1}$
(C.7)

By substituting the form of $\beta$ found in (C) and the ideal distribution of
the controller provided in equation (37), into the definition of
$\gamma(x_{t-1})$ given in (27) we find that:

$\displaystyle\gamma(x_{t-1})=\int{{}^{I}c(u_{t-1}|x_{t-1})\exp[{-\beta(u_{t-1},x_{t-1})}]\mathrm{d}u_{t-1}},$
$\displaystyle=(2\pi)^{-1/2}|\Omega|^{-1/2}\int{\exp}\bigg{[}-0.5{({u_{t-1}}-{u_{r}})^{\dagger}}{\Omega^{-1}}({u_{t-1}}-{u_{r}})-0.5x_{t-1}^{\dagger}A^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t}){Ax_{t-1}}$
$\displaystyle-0.5(P_{t}-2{{o}_{d}^{\dagger}}G_{r}^{-1}D)Ax_{t-1}-0.5u_{t-1}^{\dagger}B^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t}){Bu_{t-1}}-u_{t-1}^{\dagger}B^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t}){Ax_{t-1}}$
$\displaystyle-0.5(P_{t}-2{{o}_{d}^{\dagger}}G_{r}^{-1}D)Bu_{t-1}-0.5{\omega_{t}}-0.5{o}_{d}^{\dagger}G_{r}^{-1}{{o}_{d}}-0.5\ln\bigg{(}\dfrac{|G_{r}|}{|G|}\bigg{)}+0.5\operatorname{Tr}((G^{-1}-G_{r}^{-1})G)$
$\displaystyle-0.5x_{t-1}^{\dagger}{A}^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t}){\Sigma}A{x}_{t-1}\bigg{]}{\rm{d}}{u_{t-1}}$
$\displaystyle=(2\pi)^{-1/2}|\Omega|^{-1/2}{\exp}\bigg{(}-0.5x_{t-1}^{\dagger}A^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t}){Ax_{t-1}}-0.5(P_{t}-2{{o}_{d}^{\dagger}}G_{r}^{-1}D)Ax_{t-1}-0.5{\omega_{t}}-0.5{o}_{d}^{\dagger}G_{r}^{-1}{{o}_{d}}$
$\displaystyle-0.5\ln\bigg{(}\dfrac{|G_{r}|}{|G|}\bigg{)}+0.5\operatorname{Tr}((G^{-1}-G_{r}^{-1})G)-0.5x_{t-1}^{\dagger}{A}^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t}){\Sigma}A{x}_{t-1}-0.5u_{r}^{\dagger}\Omega^{-1}u_{r}\bigg{)}$
$\displaystyle\int{\exp}\bigg{(}-0.5u_{t-1}^{\dagger}(\Omega^{-1}+B^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t})B){u_{t-1}}+u_{t-1}^{\dagger}\big{(}\Omega^{-1}u_{r}-B^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t}){Ax_{t-1}}$
$\displaystyle-0.5B^{\dagger}(P_{t}-2{{o}_{d}^{\dagger}}G_{r}^{-1}D)^{\dagger}\big{)}\bigg{)}{\rm{d}}{u_{t-1}}$
(C.8)

By completing the square with respect to $u_{t-1}$ in the above equation, and
using the solution of the general multiple integral given in Theorem (10.5.1)
in Ref [43], it follows that:

$\displaystyle\gamma(x_{t-1})={\exp}\bigg{(}-0.5x_{t-1}^{\dagger}A^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t}){Ax_{t-1}}-0.5(P_{t}-2{{o}_{d}^{\dagger}}G_{r}^{-1}D)Ax_{t-1}-0.5{\omega_{t}}-0.5{o}_{d}^{\dagger}G_{r}^{-1}{{o}_{d}}$
$\displaystyle-0.5\ln\bigg{(}\dfrac{|G_{r}|}{|G|}\bigg{)}+0.5\operatorname{Tr}((G^{-1}-G_{r}^{-1}))-0.5x_{t-1}^{\dagger}{A}^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t}){\Sigma}A{x}_{t-1}-0.5u_{r}^{\dagger}\Omega^{-1}u_{r}\bigg{)}$
$\displaystyle\times{\exp}\bigg{(}0.5\big{(}\Omega^{-1}u_{r}-B^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t}){Ax_{t-1}}-0.5B^{\dagger}(P_{t}-2{{o}_{d}^{\dagger}}G_{r}^{-1}D)^{\dagger}\big{)}^{\dagger}$
$\displaystyle\big{(}\Omega^{-1}+B^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t})B\big{)}^{-1}\big{(}\Omega^{-1}u_{r}-B^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t}){Ax_{t-1}}-0.5B^{\dagger}(P_{t}-2{{o}_{d}^{\dagger}}G_{r}^{-1}D)^{\dagger}\big{)}$
$\displaystyle\bigg{)}\times{|\Omega|^{-1/2}}|\Omega^{-1}+B^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t})B|^{-1/2}$
(C.9)

Finally the form provided in (38) in Theorem (2) is found:

$\displaystyle-\ln{(\gamma(x_{t-1}))}$
$\displaystyle=0.5x_{t-1}^{\dagger}\bigg{(}A^{\dagger}(D^{\dagger}G_{r}^{-1}D+M_{t})A-A^{\dagger}(D^{\dagger}G_{r}^{-1}D+M_{t})^{\dagger}B\big{(}\Omega^{-1}+B^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t})B\big{)}^{-1}$
$\displaystyle
B^{\dagger}(D^{\dagger}G_{r}^{-1}D+M_{t})A+{A}^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t}){\Sigma}A\bigg{)}x_{t-1}$
$\displaystyle+0.5\bigg{(}(P_{t}-2{o}_{d}^{\dagger}G_{r}^{-1}D)A+2(\Omega^{-1}u_{r}-0.5B^{\dagger}(P_{t}^{\dagger}-2D^{\dagger}G_{r}^{-1}{o}_{d}))^{\dagger}\big{(}\Omega^{-1}+B^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t})B\big{)}^{-1}$
$\displaystyle B^{\dagger}(D^{\dagger}G_{r}^{-1}D+M_{t})A\bigg{)}x_{t-1}$
$\displaystyle+0.5\bigg{(}{\omega_{t}}+{o}_{d}^{\dagger}G_{r}^{-1}{{o}_{d}}+\ln\bigg{(}\dfrac{|G_{r}|}{|G|}\bigg{)}-\operatorname{Tr}\big{(}(G^{-1}-G_{r}^{-1})G\big{)}+u_{r}^{\dagger}\Omega^{-1}u_{r}$
$\displaystyle-\big{(}\Omega^{-1}u_{r}-0.5B^{\dagger}(P_{t}^{\dagger}-2{D}^{\dagger}G_{r}^{-1}{{o}_{d}})\big{)}^{\dagger}\big{(}\Omega^{-1}+B^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t})B\big{)}^{-1}\big{(}\Omega^{-1}u_{r}-0.5B^{\dagger}(P_{t}^{\dagger}-2{D}^{\dagger}G_{r}^{-1}{{o}_{d}})\big{)}$
$\displaystyle-2\ln({|\Omega|^{-1/2}}|\Omega^{-1}+B^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t})B|^{-1/2})\bigg{)}.$
(C.10)

##  Appendix D Calculating of the control distribution function

By substituting (37), (C) and (C) into (26), it follows that the optimal
control distribution is given by:

$\displaystyle
c(u_{t-1}|x_{t-1})=(2\pi)^{-1/2}|\Omega^{-1}+B^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t})B|^{1/2}\exp\bigg{[}-0.5u_{t-1}^{\dagger}\bigg{(}\Omega^{-1}+B^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t})B\bigg{)}{u_{t-1}}$
$\displaystyle+u_{t-1}^{\dagger}\bigg{(}\Omega^{-1}u_{r}-B^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t}){Ax_{t-1}}-0.5B^{\dagger}(P_{t}-2{{o}_{d}^{\dagger}}G_{r}^{-1}D)^{\dagger}\bigg{)}-0.5\bigg{(}{\Omega^{-1}}{u_{r}}-{B^{\dagger}}({D^{\dagger}}G_{r}^{-1}D+{M_{t}})A{x_{t-1}}$
$\displaystyle-0.5{B^{\dagger}}(P_{t}^{\dagger}-2{D^{\dagger}}G_{r}^{-1}{{o}_{d}})\bigg{)}^{\dagger}\bigg{(}{\Omega^{-1}}+{B^{\dagger}}({D^{\dagger}}G_{r}^{-1}D+{M_{t}})B\bigg{)}^{-1}\bigg{(}{\Omega^{-1}}{u_{r}}-{B^{\dagger}}({D^{\dagger}}G_{r}^{-1}D+{M_{t}})A{x_{t-1}}$
$\displaystyle-0.5{B^{\dagger}}(P_{t}^{\dagger}-2{D^{\dagger}}G_{r}^{-1}{{o}_{d}})\bigg{)}\bigg{]}.$
(D.1)

which can be written as:

$\displaystyle c(u_{t-1}|x_{t-1})$
$\displaystyle=(2\pi)^{-1/2}|\Omega^{-1}+B^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t})B|^{1/2}$
(D.2)
$\displaystyle\times\exp\bigg{[}-{0.5({{u}_{t-1}}-{{v}_{t-1}})^{\dagger}}\bigg{(}\Omega^{-1}+B^{\dagger}({D^{\dagger}}G_{r}^{-1}D+M_{t})B\bigg{)}({u_{t-1}}-{v_{t-1}})\bigg{]},$

where,

$\displaystyle v_{t-1}$
$\displaystyle=\bigg{(}{\Omega^{-1}}+{B^{\dagger}}({D^{\dagger}}G_{r}^{-1}D+{M_{t}})B\bigg{)}^{-1}\bigg{(}{\Omega^{-1}}{u_{r}}-{B^{\dagger}}({D^{\dagger}}G_{r}^{-1}D+{M_{t}})A{x_{t-1}}-0.5{B^{\dagger}}(P_{t}^{\dagger}-2{D^{\dagger}}G_{r}^{-1}{{o}_{d}})\bigg{)},$
(D.3)

This means that the pdf of the controller follows a normal distribution, as
given by:

$c(u_{t-1}|x_{t-1})\sim\mathcal{N}(v_{t-1},R_{t}),$ (D.4)

where

$\displaystyle R_{t}$
$\displaystyle=\bigg{(}{\Omega^{-1}}+{B^{\dagger}}({D^{\dagger}}G_{r}^{-1}D+{M_{t}})B\bigg{)}^{-1}.$
(D.5)

##  Appendix E State space model for spin $\frac{1}{2}$

The evolution of a spin-1/2 system in interaction with an external environment
can be described by the following master equation:

$\dfrac{d\rho(t)}{dt}=-i[H,\rho(t)]+\Theta\bigg{(}\sigma_{-}\rho(t)\sigma_{+}-\dfrac{1}{2}\\{\sigma_{+}\sigma_{-},\rho(t)\\}\bigg{)},$
(E.1)

where, as given in equation (46),
$H=\dfrac{1}{2}\sigma_{3}+\dfrac{1}{2}(\sigma_{1}+\sigma_{2}){u}_{t}$ such
that ${u}_{t}$ is an external electric field, and
$\sigma_{1},\sigma_{2},\sigma_{3}$ are the relevant Pauli matrices and
$\sigma_{\pm}=\dfrac{\sigma_{1}\pm i\sigma_{2}}{2}$. The Hamiltonian $H$ can
be written in matrix form as follows:

$H=\dfrac{1}{2}\left(\begin{array}[]{cc}1&{u}_{t}(1-i)\\\ {u}_{t}(1+i)&-1\\\
\end{array}\right).$ (E.2)

Writing the density operator in terms of its elements, the von-Neumann
equation (E.1) can be re-written as:

$\displaystyle\dfrac{d}{dt}\left(\begin{array}[]{cc}\rho_{00}(t)&\rho_{01}(t)\\\
\rho^{*}_{01}(t)&\rho_{11}(t)\end{array}\right)=\dfrac{-i}{2}\left(\begin{array}[]{cc}{u}_{t}\big{(}(1-i)\rho^{*}_{01}(t)-(1+i)\rho_{01}(t)\big{)}&2\rho_{01}(t)+{u}_{t}(1-i)(\rho_{11}(t)-\rho_{00}(t))\\\
-2\rho^{*}_{01}(t)+{u}_{t}(1+i)(\rho_{00}(t)-\rho_{11}(t))&{u}_{t}\big{(}(1+i)\rho_{01}(t)-(1-i)\rho^{*}_{01}(t)\big{)}\\\
\end{array}\right)$ (E.7)
$\displaystyle-\left(\begin{array}[]{cc}\Theta\rho_{00}(t)&\dfrac{\Theta}{2}\rho_{01}(t)\\\
\dfrac{\Theta}{2}\rho^{*}_{01}(t)&-\Theta\rho_{00}(t)\end{array}\right)$
(E.10)

where the elements $\rho_{00}(t),\rho_{01}(t),\rho^{*}_{01}(t)$ and
$\rho_{11}(t)$ are the elements of the density operator $\rho(t)$. Vectorizing
the above equation yields:

$\displaystyle\dfrac{d}{dt}\underbrace{\left(\begin{array}[]{cc}\rho_{00}(t)\\\
\rho_{11}(t)\\\ \rho_{01}(t)\\\
\rho^{*}_{01}(t)\end{array}\right)}_{x(t)}=\underbrace{\left(\begin{array}[]{cccc}-\Theta&0&0&0\\\
\Theta&0&0&0\\\ 0&0&-i-\dfrac{\Theta}{2}&0\\\ 0&0&0&i-\dfrac{\Theta}{2}\\\
\end{array}\right)}_{\tilde{A}}\underbrace{\left(\begin{array}[]{cc}\rho_{00}(t)\\\
\rho_{11}(t)\\\ \rho_{01}(t)\\\ \rho^{*}_{01}(t)\end{array}\right)}_{x(t)}$
(E.23)
$\displaystyle+i\underbrace{\dfrac{1}{2}\left(\begin{array}[]{cccc}0&0&(1+i)&-(1-i)\\\
0&0&-(1+i)&(1-i)\\\ (1-i)&-(1-i)&0&0\\\ -(1+i)&(1+i)&0&0\\\
\end{array}\right)}_{\tilde{N}}\underbrace{\left(\begin{array}[]{cc}\rho_{00}(t)\\\
\rho_{11}(t)\\\ \rho_{01}(t)\\\
\rho^{*}_{01}(t)\end{array}\right)}_{x(t)}{u}_{t},$ (E.32)

from which we find the state equation for the spin-1/2 system to be given by:

$\dfrac{d{x}(t)}{dt}=(\tilde{A}+i\tilde{N}{u}_{t}){x}(t).$ (E.33)

##  Appendix F State space model for a three level system

The Hamiltonian describing the interaction between a $\Lambda$-type atomic
system, examined in Section (5.1.2), and an electric field ${u}_{t}$ can be
written in matrix form as follows:

$H=H_{0}+{u}_{t}H_{1}=\left(\begin{array}[]{ccc}\dfrac{3}{2}&0&0\\\ 0&1&0\\\
0&0&0\\\ \end{array}\right)+\left(\begin{array}[]{ccc}0&0&1\\\ 0&0&1\\\
1&1&0\\\
\end{array}\right){u}_{t}=\left(\begin{array}[]{ccc}\dfrac{3}{2}&0&{u}_{t}\\\
0&1&{u}_{t}\\\ {u}_{t}&{u}_{t}&0\\\ \end{array}\right)$ (F.1)

By assuming that the system interacts with an external environment described
by one Linblad operator, $L_{02}=\sqrt{\Theta}\ket{0}\bra{2}$, and expanding
the density operator in terms of its matrix elements the von-Neumann equation
given in (1) can be written as:

$\displaystyle
i\dfrac{d}{dt}\left(\begin{array}[]{ccc}\rho_{00}(t)&\rho_{01}(t)&\rho_{02}(t)\\\
\rho^{*}_{01}(t)&\rho_{11}(t)&\rho_{12}(t)\\\
\rho^{*}_{02}(t)&\rho^{*}_{12}(t)&\rho_{22}(t)\end{array}\right)=$ (F.5)
$\displaystyle\left(\begin{array}[]{ccc}u_{t}\big{(}\rho^{*}_{02}(t)-\rho_{02}(t)\big{)}&\dfrac{1}{2}\rho_{01}(t)+u_{t}\big{(}\rho^{*}_{12}(t)-\rho_{02}(t)\big{)}&\dfrac{3}{2}\rho_{02}(t)+u_{t}\big{(}\rho_{22}(t)-\rho_{00}(t)-\rho_{01}(t)\big{)}\\\
\\\
-\dfrac{1}{2}\rho^{*}_{01}(t)+{u}_{t}(\rho^{*}_{02}(t)-\rho_{12}(t))&{u}_{t}\big{(}\rho_{12}^{*}(t)-\rho_{12}(t)\big{)}&\rho_{12}(t)+{u}_{t}\big{(}\rho_{22}(t)-\rho_{01}^{*}(t)-\rho_{11}(t)\big{)}\\\
\\\
-\dfrac{3}{2}\rho_{02}^{*}(t)+{u}_{t}\big{(}\rho_{00}+\rho_{01}^{*}(t)-\rho_{22}(t)\big{)}&-\rho_{12}^{*}(t)+{u}_{t}(\rho_{01}(t)+\rho_{11}(t)-\rho_{22}(t))&{u}_{t}\big{(}\rho_{02}(t)+\rho_{12}(t)-\rho_{02}^{*}(t)-\rho_{12}^{*}(t)\big{)}\\\
\end{array}\right)$ (F.11)
$\displaystyle+\dfrac{i}{2}\Theta\left(\begin{array}[]{ccc}-2\rho_{00}(t)&-\rho_{01}(t)&-\rho_{02}(t)\\\
-\rho_{01}^{*}(t)&0&0\\\ -\rho^{*}_{02}(t)&0&2\rho_{00}(t)\end{array}\right).$
(F.15)

which when using the vectorisation defined in Eq.(A), gives,

$\displaystyle\dfrac{d}{dt}\underbrace{\left(\begin{array}[]{ccccccccc}\rho_{00}(t)\\\
\rho_{11}(t)\\\ \rho_{22}(t)\\\ \rho_{01}(t)\\\ \rho_{02}(t)\\\
\rho_{01}^{*}(t)\\\ \rho_{02}^{*}(t)\\\ \rho_{12}(t)\\\ \rho_{12}^{*}(t)\\\
\end{array}\right)}_{x(t)}=\underbrace{\left(\begin{array}[]{ccccccccc}-\Theta&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0\\\ \Theta&0&0&0&0&0&0&0&0\\\
0&0&0&-\dfrac{i}{2}-\dfrac{\Theta}{2}&0&0&0&0&0\\\
0&0&0&0&-\dfrac{3i}{2}-\dfrac{\Theta}{2}&0&0&0&0\\\
0&0&0&0&0&\dfrac{i}{2}-\dfrac{\Theta}{2}&0&0&0\\\
0&0&0&0&0&0&\dfrac{3i}{2}-\dfrac{\Theta}{2}&0&0\\\ 0&0&0&0&0&0&0&-i&0\\\
0&0&0&0&0&0&0&0&i\\\
\end{array}\right)}_{\tilde{A}}\underbrace{\left(\begin{array}[]{ccccccccc}\rho_{00}(t)\\\
\rho_{11}(t)\\\ \rho_{22}(t)\\\ \rho_{01}(t)\\\ \rho_{02}(t)\\\
\rho_{01}^{*}(t)\\\ \rho_{02}^{*}(t)\\\ \rho_{12}(t)\\\ \rho_{12}^{*}(t)\\\
\end{array}\right)}_{x(t)}$ (F.43)
$\displaystyle+i\underbrace{\left(\begin{array}[]{ccccccccc}0&0&0&0&1&0&-1&0&0\\\
0&0&0&0&0&0&0&1&-1\\\ 0&0&0&0&-1&0&1&-1&1\\\ 0&0&0&0&1&0&0&0&-1\\\
1&0&-1&1&0&0&0&0&0\\\ 0&0&0&0&0&0&-1&1&0\\\ -1&0&1&0&0&-1&0&0&0\\\
0&1&-1&0&0&1&0&0&0\\\ 0&-1&1&-1&0&0&0&0&0\\\
\end{array}\right)}_{\tilde{N}}\underbrace{\left(\begin{array}[]{ccccccccc}\rho_{00}(t)\\\
\rho_{11}(t)\\\ \rho_{22}(t)\\\ \rho_{01}(t)\\\ \rho_{02}(t)\\\
\rho_{01}^{*}(t)\\\ \rho_{02}^{*}(t)\\\ \rho_{12}(t)\\\ \rho_{12}^{*}(t)\\\
\end{array}\right)}_{x(t)}{u}_{t}.$ (F.62)

Hence, we find the form of the state equation given in Eq.(9) as follows,

$\dfrac{d{x}(t)}{dt}=(\tilde{A}+i\tilde{N}u(t)){x}(t).$ (F.63)

## References

  * [1] D. G. Cory, A. F. Fahmy, and T. F. Havel,  _Proc. Natl. Acad. Sci. USA_ , 94, 1634 (1997).
  * [2] L. M. Duan, M. D. Lukin, J. I. Cirac, and P. Zoller, _Nature (London)_ 414, 413 (2001)
  * [3] G. Bensky, D. Petrosyan, J. Majer, J. Schmiedmayer, and G. Kurizki, _Phys. Rev. A_ , 86, 012310 (2012).
  * [4] J. S. Li, J. Ruths, T. Y. Yu, H. Arthanari, and G. Wagner,  _Proc. Natl. Acad. Sci. USA_ , 108, 1879 (2011).
  * [5] S. J. Glaser, T. Schulte-Herbrüggen, M. Sieveking, O. Schedletzky, N. C. Nielsen, O. W. Sorensen, and C. Griesinger, _Science_ 280, 421 (1998).
  * [6] C. Altafini, _IEEE Trans. Automat. Control_ 52, 2019 (2007).
  * [7] G. Turinici and H. Rabitz,  _Phys. Rev. A_ 70, 063412 (2004).
  * [8] J. S. Li and N. Khaneja, _Phys. Rev. A_ 73, 030302 (R) (2006)
  * [9] K. Beauchard, P.S. Pereira da Silva and P. Rouchon, _Automatica_. 48, 68 (2012).
  * [10] J. S. Li and N. Khaneja, _IEEE Trans. Automat. Control_ 54, 528 (2009).
  * [11] Chunlin Chen, Daoyi Dong, Ruixing Long, Ian R. Petersen, and Herschel A. Rabitz, _Phys. Rev. A._ 89, 023402 (2014).
  * [12] Chunlin Chen, Daoyi Dong, Bo Qi, Ian R. Petersen, and Herschel A. Rabitz, _IEEE Transactions on Neural Networks and Learning Systems._ 28, 1345 (2017).
  * [13] H. Jirari H and W. Pötz, _Phys. Rev. A_ 72, 013409 (2005).
  * [14] W. Cui, Z. R. Xi and Y. Pan, _Phys. Rev. A_ 77, 032117 (2008).
  * [15] U. Hohenester and G. Stadler, _Phys. Rev. Lett._ 92, 196801 (2004).
  * [16] G. Gordon, _J. Phys. B: At. Mol. Opt. Phys._ 42, 223001 (2009).
  * [17] M. D. Grace, J. Dominy, R. L. Kosut, C. Brif and H. Rabitz  _New J. Phys._ 12, 015001 (2010).
  * [18] I. A. Grigorenko and D. V. Khveshchenko  _Phys. Rev. Lett._ 95, 110501 (2005).
  * [19] W. Zurek,  _Rev. Mod. Phys._ 75, 715 (2003).
  * [20] E. Joos. et al. , “Decoherence and the Appearance of a Classical World in Quantum Theory",_Springer, Berlin, second edn._ (2003).
  * [21] M. Schlosshauer, “ Decoherence and the Quantum-to-Classical Transition",_Springer, Berlin_ (2007).
  * [22] M. A. Nielsen and I. L. Chuang , “Quantum Computation and Quantum Information",_Cambridge, UK: Cambridge University Press_ (2000).
  * [23] C. Brif, H. Rabitz, S. Wallentowitz, and I. A. Walmsley,  _Phys. Rev. A._ 63, 063404 (2001).
  * [24] W. Zhu and H. Rabitz,  _J. Chem. Phys._ 118, 6751 (2003).
  * [25] M. P. A. Branderhorst, P. Londero, P. Wasylczyk, C. Brif, R. L. Kosut, H. Rabitz and I. A. Walmsley,  _Science_ 320, 638 (2008).
  * [26] M. Kárný  _Automatica_. 32, 1719 (1996).
  * [27] R. Herzallah and M. Kárný, _Neural Networks_. 24, 1128 (2011).
  * [28] R. Herzallah, _Neural Networks_. 42, 74 (2013).
  * [29] M. Kárný ad T. V. Guy  _Syst. Control Lett_. 55, 259 (2006).
  * [30] R. Herzallah, _Neural Networks_. 63, 199 (2015).
  * [31] R. Herzallah, _Asian Journal of Control_. 20, 2065 (2018).
  * [32] Y. Zhou and R. Herzallah, _AIMS Electronics and Electrical Engineering._ 4, 216 (2020).
  * [33] K. Kraus, A. Böhm, J. D. Dollard, W. H. Wootters, “ Effects and Operations: Fundamental Notions of Quantum Theory (Lecture Notes in Physics vol 190)",_Berlin: Springer_ (1983).
  * [34] H. P. Breuer and F. Petruccione, “The Theory of Open Quantum Systems",_New York: Oxford University Press_ (2007).
  * [35] C. W. Gardiner and P. Zoller, “ Quantum Noise: A Handbook of Markovian and Non-Markovian Quantum Stochastic Methods with Applications to Quantum Optics",_Berlin: Springer_ (1999).
  * [36] J. von Neumann, “Mathematical Foundations of Quantum Mechanics", _Princeton University Press_ , (1996).
  * [37] Todd A. Brun  _Phys. Rev. A_. 61, 042107 (2000).
  * [38] A. C. Doherty and K. Jacobs, _Phys. Rev. A_. 60, 2700 (1999).
  * [39] B. Schmidt and C. Hartmann, _Computer Physics Communications_. 228, 229 (2018).
  * [40] Boris Schäfer-Bung, Carsten Hartmann, Burkhard Schmidt, and Christof Schütte,_J. Chem. Phys._ 135, 014112 (2011).
  * [41] C. A. Ryan, C. Negrevergne, M. Laforest, E. Knill and R. Laflamme  _Phys. Rev. A_. 78, 012328 (2008).
  * [42] N. Khaneja, T. Reiss, C. Kehlet, E. Knill, T. Schulte-Herbrüggen and S. J. Glaser  _J. Magn. Reson_. 172, 296 (2005).
  * [43] F. A. Graybill, “Matrices with Applications in Statistics",_Belmont, Calif. : Wadsworth International Group_ (1983).