# Optimal quantum control via genetic algorithms

Jonathon Brown  Centre for Theoretical Atomic, Molecular, and Optical Physics,
School of Mathematics and Physics, Queens University, Belfast BT7 1NN, United
Kingdom  Mauro Paternostro  Centre for Theoretical Atomic, Molecular, and
Optical Physics, School of Mathematics and Physics, Queens University, Belfast
BT7 1NN, United Kingdom  Alessandro Ferraro  Centre for Theoretical Atomic,
Molecular, and Optical Physics, School of Mathematics and Physics, Queens
University, Belfast BT7 1NN, United Kingdom  Dipartimento di Fisica “Aldo
Pontremoli,” Università degli Studi di Milano, I-20133 Milano, Italy

###### Abstract

We devise a machine learning-enabled approach to quantum state engineering
based on evolutionary algorithms. In particular, we focus on superconducting
platforms and consider a network of qubits – encoded in the states of
artificial atoms with no direct coupling – interacting via a common single-
mode driven microwave resonator. The qubit-resonator couplings are assumed to
be in the resonant regime and tunable in time. A genetic algorithm is designed
in order to find the functional time-dependence of the couplings that optimise
the fidelity between the evolved state and a variety of targets, including
three-qubit GHZ and Dicke states and four-qubit graph states. We observe high
quantum fidelities (above 0.96 in the worst case setting of a system of
effective dimension 96) and resilience to noise, despite the algorithm being
trained in the ideal noise-free setting. These results show that the genetic
algorithms represent an effective approach to control quantum systems of large
dimensions.

##  I Introduction

Quantum state engineering is an essential enabling step for a variety of
quantum information tasks, including the initialization of quantum simulators
Georgescu  _et al._ (2014), the loading of classical data for quantum-enhanced
analysis Aaronson (2015), or the generation of resourceful states in quantum
communication networks Muralidharan  _et al._ (2016). In particular, quantum
entangled states, which embody a stark departure from classicality, often
provide the main resources towards quantum advantage Horodecki  _et al._
(2009). Therefore the preparation of entangled states in multi-node quantum
networks presents a key challenge for realising quantum protocols on near term
devices.

A successful approach towards state and resource generation consists of two
steps: (i) engineer a suitable time-dependent Hamiltonian with tunable
parameters, as allowed by current experimental capabilities in a given
physical platform; (ii) find and implement proper temporal dependencies (pulse
shapes) for these parameters by invoking optimal quantum control techniques
D’Alessandro (2007); Glaser  _et al._ (2015). This results in tailor-made
control schemes pertinent to the specific platform at hand.

One platform at the forefront of engineering flexible multi-node quantum
networks, to which such approach has been successfully applied, is that of
superconducting quantum circuits Clarke and Wilhelm (2008); Kjaergaard  _et
al._ (2020); Huang  _et al._ (2020); Koch  _et al._ (2007). So far
superconducting architectures have been employed to realise two-qubit gates
using frequency-tunable DiCarlo  _et al._ (2009); Barends  _et al._ (2014); Li
_et al._ (2019); Barends  _et al._ (2019) and microwave-driven Chow  _et al._
(2011); Córcoles  _et al._ (2013) artificial atoms. In addition, demonstration
of the coupling between artificial atoms and microwave resonators Majer  _et
al._ (2007) opened the door for resonator mediated two-qubit gates Cross and
Gambetta (2015); Puri and Blais (2016); Paik  _et al._ (2016) and provided an
alternative platform to study cavity quantum electrodynamics Blais  _et al._
(2007) leading to the field of Circuit QED Blais  _et al._ (2021). Whilst
extremely flexible in their design, it has been shown however that operating
these superconducting systems with a reduced level control is not only
desirable, but necessary in some cases Falci  _et al._ (2017); Di Stefano  _et
al._ (2015). Thus finding optimal control protocols that utilize a limited but
effective level of control is of practical interest.

Optimal control of quantum systems has yielded a range of new methods
inspired, in part, by the development of modern machine learning methods.
Specifically, neural-network-based reinforcement learning methods Sutton and
Barto (2018); Mnih  _et al._ (2013) have been recognised as useful tools to
study quantum systems Giannelli  _et al._ (2022) in a variety of contexts
including state transfer Paparelle  _et al._ (2020); Porotti  _et al._ (2019);
Sgroi  _et al._ (2021); Brown  _et al._ (2021), quantum thermodynamics Sgroi
_et al._ (2021), circuit architecture search Kuo  _et al._ (2021) and control
of dissipative systems An  _et al._ (2021). Reinforcement learning techniques
have proven particularly suitable for control problems of increasing dimension
when compared to more standard techniques Zhang  _et al._ (2019). However, for
the most part these techniques cannot be used as closed-loop optimization
schemes and therefore are of limited use for optimization on physical quantum
systems 111Noticeable exceptions are discussed in Refs. Sgroi  _et al._
(2021); Brown  _et al._ (2021) and have relatively poor convergence guarantees
which necessitate, often expensive, hyper parameter tuning steps.
Interestingly, evolutionary strategies have been proposed as a scalable
substitute to reinforcement learning methods Salimans  _et al._ (2017) and
have already been used as an alternative to gradient based parameter updates
in both deep reinforcement learning Such  _et al._ (2017) and quantum
reservoir computing Ghosh  _et al._ (2021). Their possible applications to
state and resource generation is however still vastly unexplored. The main aim
of this work is to investigate the potential of such an approach when applied
to multiple qubits in quantum networks.

Here we consider a register of qubits, coupled using a reduced-control
architecture, and employ a genetic algorithm to find optimal pulse sequences
to drive their dynamics. In order to illustrate our approach, we present
efficient control schemes for preparing entangled three- and four-qubit
states, including GHZ, Dicke, and graph states, and assess performance against
relevant decoherence sources finding the thresholds that limit the quality of
our results. The generation of high-quality states is thus demonstrated while
identifying fully the sequence of driving pulses to use. Our approach holds
the potential to establish a hybrid take to quantum control that mixes machine
learning and optimal control as a viable route to the engineering of crucial
resources for quantum information processing.

The remainder of the manuscript is organised as follows. In Sec. II we present
the specifics of the system considered and formalise the Hamiltonian. In Sec.
III, we present an overview of the Continuous Genetic Algorithm employed then
follow by formalising the algorithm for the problem of quantum control in
section IV. Sec. V focuses on the presentation of the resulting control
schemes for the preparation of three-qubit GHZ and Dicke states as well as a
specific instance of a four-qubit graph states. In Sec. V.3, the effect of
decoherence is investigated. Sec. VI offers our closing remarks and a forward
look.

##  II System

We consider a system composed of $N$ identical and non-interacting qubits
coupled via a common single-mode driven resonator. We model such system with
the Hamiltonian $H=H_{0}+H_{int}+H_{d}$, where

$\displaystyle
H_{0}=\omega_{c}a^{\dagger}a+\sum_{j=1}^{N}\omega_{j}\sigma_{j}^{+}\sigma_{j}^{-},$
(1) $\displaystyle
H_{int}=\sum_{j=1}^{N}g_{j}(a^{\dagger}\sigma_{j}^{-}+a\sigma_{j}^{+}),$ (2)
$\displaystyle H_{d}=\xi(ae^{i\omega_{d}t}+a^{\dagger}e^{-i\omega_{d}t}).$ (3)

In writing these Hamiltonians, we have assumed the rotating wave approximation
and used units such that $\hbar=1$. Here $a$ denotes the annihilation operator
of the resonator mode, whereas $\sigma_{j}^{+}=|1\rangle\langle 0|$ and
$\sigma_{j}^{-}=|0\rangle\langle 1|$ are the raising and lowering operators of
the $j$-th qubit. Moreover, $\omega_{c}$ denotes the resonator frequency,
$\omega_{j}$ the transition frequency of the $j$-th qubit and $g_{j}$ its
coupling strength with the resonator. The driving amplitude (assumed to be
real) and the carried frequency of the drive are indicated with $\xi$ and
$\omega_{d}$ respectively.

The use of a resonator to mediate two-qubit gate interactions is well studied.
For example, in the context of superconducting systems the Resonator-Induced
Phase (RIP) gate Puri and Blais (2016); Cross and Gambetta (2015); Paik  _et
al._ (2016) utilizes two qubits dispersively coupled to a common microwave
resonator. In such a dispersive regime, the coherent qubit-resonator
interaction term becomes negligible leading to an effective qubit-qubit
interaction term. This is beneficial for protection against decoherence at the
expense of gate operation time. On the other end of the spectrum, resonant
regimes have been considered to selectively tune qubits in and out of coupling
with a common resonator Haack  _et al._ (2010). In addition, external
microwave drives are also commonly used such as with the the cross-resonance
gate Chow  _et al._ (2011); Córcoles  _et al._ (2013). In contrast to such
previous works, we assume in $H$ a fully resonant regime and adopt a limited
control scheme, in that the drive amplitude of the resonator mode and the
qubit-resonator couplings are tunable.

Taking an interaction picture to the frequency of the harmonic mode, using
$\tilde{H}=e^{iRt}He^{-iRt}-R$, where
$R=\omega_{c}\left(a^{\dagger}a+\sum_{j=1}^{N}\sigma_{j}^{+}\sigma_{j}^{-}\right)$,
we get

$\tilde{H}=\sum_{j=1}^{N}\left[\delta_{j}\sigma_{j}^{+}\sigma_{j}^{-}+g_{j}(a^{\dagger}\sigma_{j}^{-}+a\sigma_{j}^{+})\right]+\xi(ae^{i\delta_{d}t}+a^{\dagger}e^{-i\delta_{d}t}),$
(4)

where $\delta_{j}=\omega_{j}-\omega_{c}$ and
$\delta_{d}=\omega_{d}-\omega_{c}$ are the detunings between the
$j^{\text{th}}$ qubit and the harmonic mode, and between the drive frequency
and the harmonic mode respectively. The usefulness of this transformation
becomes apparent if one assume full resonance between qubit transition
frequencies, drive frequency and resonator frequency – namely, setting
$\omega_{c}=\omega_{d}=\omega_{j}$, $\forall j=1,...,N$. The Hamiltonian thus
takes the simpler form

$\tilde{H}=\sum_{j=1}^{N}g_{j}(t)(a^{\dagger}\sigma_{j}^{-}+a\sigma_{j}^{+})+\xi(t)(a+a^{\dagger}),$
(5)

which comprises $N+1$ terms embodying to an equal number of controls, where we
assume each of the qubit-cavity couplings $g_{j}(t)$ and the drive amplitude
$\xi(t)$ to be time-dependent controllable parameters.

The above Hamiltonian governs the dynamics of the system via the time-
dependent Schrödinger equation. Thus if the system is initially in some state
$|\psi(t_{0})\rangle$, then the time-evolved state of the system at any future
time $t>t_{0}$ is given by 222An interesting consequence of the particular
form of $\tilde{H}$ in Eq. (5) is that, if the system is initialized in a pure
state with real coefficients, then the time-evolved reduced state of the qubit
will always be described by real coefficients. This is of course a restriction
if one wishes to prepare more general states, however it simplifies the
optimization when considering completely real target states. One can ease this
constraint if necessary by introducing non-zero detunings in Eq. (4)

$|\psi(t)\rangle=\mathcal{T}e^{i\int_{t_{0}}^{t}\tilde{H}(t^{\prime})dt^{\prime}}|\psi(t_{0})\rangle.$
(6)

The goal here is to find optimal functional forms for $g_{j}(t)$ and $\xi(t)$,
such that the system is dynamically steered in some desired way. Specifically,
we are interested in state preparation within the qubit subspace, so we first
determine the reduced state of the qubit network

$\displaystyle\rho_{Q}(t)=\text{Tr}_{c}\Big{(}|\psi(t)\rangle\langle\psi(t)|\Big{)}=\sum_{i=0}^{\infty}{}_{c}\langle{i}|\psi(t)\rangle\langle\psi(t)|{i}\rangle_{c},$
(7)

and work to find $g_{j}(t)$ and $\xi(t)$ such that the fidelity

$\mathcal{F}(\rho_{Q},\sigma)=\langle\psi_{\sigma}|\rho_{Q}|\psi_{\sigma}\rangle$
(8)

is maximized, where $|\psi_{\sigma}\rangle$ is the target state of interest.
It is worth stressing that our approach would work equally – mutatis mutandis
– with mixed target states.

##  III The Continuous Genetic Algorithm

Evolutionary strategies are a class of direct search optimization techniques,
drawing inspiration from Darwinian evolution, that have recently been proposed
as a viable substitute for gradient based parameter optimization in Neural
Networks Such  _et al._ (2017) and quantum reservoirs Ghosh  _et al._ (2021),
as well as a scalable alternative to Reinforcement learning techniques
Salimans  _et al._ (2017). Of particular interest to continuous-control
problems is the so-called “Continuous Genetic Algorithm” (CGA) Haupt and Haupt
(2004), which generalises the more traditional Discrete Genetic algorithm to
allow for continuous parameter values.

Consider an optimization problem with $\mathcal{N}_{var}$ parameter variables
$p_{i}$. We call a specific instance of these parameters,
$\mathcal{C}=[p_{1},p_{2},...,p_{\mathcal{N}_{var}}]$, a Chromosome. This
embodies one proposed solution to the optimization problem. One then defines a
Fitness function, $f(\mathcal{C})=f({p_{1},p_{2},...,p_{N_{var}}})$:
$\mathbb{R}^{N_{var}}\to\mathbb{R}$. This function will be determined by the
optimization task under consideration and will assign a numerical score to
each proposed solution (chromosome) depending on its usefulness (fitness) with
respect to the specified task at hand. In the name of generality we assert
that $p_{i}\in\left[-1,1\right]$, where the parameter values are suitably
scaled within the fitness function. We call each new iteration of the
algorithm a Generation where the algorithm proposes several chromosomes to
make up a Population. The algorithm can be qualitatively summarised using the
following steps:

  1. 1.

Initialization. Define the fitness function (determined by the optimization
task) and fix the hyper-parameters for the genetic algorithm. One then
generates an initial population with $\mathcal{N}_{pop}$ chromosomes,
typically randomly, which acts as the zeroth generation of the algorithm.

  2. 2.

Natural Selection. The fitness of each chromosome is assessed with a call to
the fitness function. The population is then ranked based on these fitness
values and the $\mathcal{N}_{survive}$ highest scoring chromosomes are chosen
to survive, while the rest are discarded to be replaced by new offspring
chromosomes in the next generation.

  3. 3.

Pairing. From the $\mathcal{N}_{survive}$ surviving chromosomes chose
$\mathcal{N}_{parents}$ pairs of parent chromosomes in order to produce
$(\mathcal{N}_{pop}-\mathcal{N}_{survive})$ offspring to repopulate. Parent
chromosomes are sampled probabilistically based on their relative fitness,
such that the fittest chromosomes reproduce more frequently, where cloning
(using the same chromosome for both parents) is disallowed.

  4. 4.

Mating. Here, each pair of parent chromosomes is combined in some manner to
produce enough offspring to replenish the population. In the simplest case,
one can consider cutting the parent chromosomes in half and producing two
offspring, made up of the two possible combinations of these halves.
Alternatively, one could pick indices along the chromosomes at random, and
generate two offspring by first copying the parent chromosomes, then swapping
the parameter values corresponding to these indices between the two copies.
This direct transfer of parameter values to generate offspring is the simplest
and most obvious way of mating two parent chromosomes. However this approach
simply provides offspring chromosomes made up entirely of parameter values
that were already present in the previous generation and, as such, introduces
no new “genetic material” in the form of novel parameter values for a given
parameter. To tackle this we can implement a combination of the form

$p_{i,\textit{ offspring}}=\beta p_{i,\textit{ parent
1}}+(1-\beta)p_{i,\textit{ parent 2}}.$ (9)

This allows us to introduce an element of exploration by allowing not only
parameter values present in the previous generation but also any continuous
value in between, modulated by random variable $\beta\in[0,1]$.

  5. 5.

Mutation. This final step introduces further exploration into the search by
choosing, at random, a number of elements within each chromosome to be
replaced with a new random value. The rate of this mutation is set by a
parameter $0<\alpha<1$, which determines the number of indices to be targeted
relative to $\mathcal{N}_{var}$. This mutation step is applied to the entire
population except for the single chromosome with the highest fitness, which is
known as Elitism. This is important to ensure theoretical guarantees of
convergence. Specifically it ensures that the maximum achieved fitness is
always at least maintained in new generations.

The algorithm repeats steps 2-5 until convergence or an acceptable level of
fitness has been achieved. The aforementioned hyper-parameters associated with
CGA then are: population size $\mathcal{N}_{pop}$, number of survivors
$\mathcal{N}_{survive}$ to keep in each iteration, number of parental pairs to
mate $\mathcal{N}_{parents}$ and the mutation rate $\alpha$. Also, as
discussed above, we have some freedom in how we implement the mating
procedure.

##  IV Continuous Genetic Algorithms for Optimal Quantum Control

In order to apply the CGA we first must formulate the optimal control problem
in a suitable manner. As is common practice, we discretize the evolution time
into $T$ time intervals of equal duration $\tau$ — which is manually chosen —
and assume the functional form for each control to be defined by its values at
the $T+1$ times $t=0,\tau,2\tau,...,T\tau$ (connecting the latter with a
simple tanh function, see Appendix A for specifics). Therefore, given that
each control function is completely defined by these $T+1$ values, and there
are $N+1$ controls in total as per Eq. (5), then the total number of
parameters that we have to optimize over is $\mathcal{N}_{var}=(N+1)(T+1)$
(which as said is the length of the chromosomes). We then need to outline how
we asses the fitness of each chromosome and in doing so define the
optimization task.

Figure 1: Schematic representation of how the variables within each chromosome
(visually represented by the sequence of connected squares) are assigned to
which control pulse, depicted along the top of the diagram. Along the bottom
shows sample control pulses generated using random parameter values and the
function construction method outlined in Appendix A.

###  IV.1 The Fitness Function

As said, each chromosome is a sequence of $\mathcal{N}_{var}$ parameter
values, and there are $N+1$ control pulses. Therefore we first break each
chromosome up into $N+1$ sequences where the first $T+1$ parameters determine
the first control pulse and so on, as in Fig. 1. Then the parameter values
corresponding to each control are scaled accordingly where each of the qubit-
resonator couplings assume a common range $g_{j}\in[-g_{0},g_{0}]$, and the
drive amplitude is in the range $\xi\in[-\xi_{0},\xi_{0}]$. The scaled
parameter values are then used to build the control pulses as in Appendix A,
which fully define the time dependent Hamiltonian $\tilde{H}(t)$ in Eq. 5.
After defining an initial state $|\psi(t=0)\rangle$ the system evolves
unitarily according to $\tilde{H}(t)$ (as in Eq. 6), where we keep track of
the state of the system at every intermediate time between $t=0$ and $t=\tau
T$. For each state in this history of states, we trace out the cavity system
to obtain a state history for the qubit subspace, $\rho_{Q}(t)$. Next we
calculate $\mathcal{F}(\rho_{Q}(t),\sigma)$, which is the time dependent
fidelity induced by the specific control pulses. In order to assign a
numerical fitness value to the chromosome we simply take the maximum fidelity
achieved in the qubit subspace throughout the induced dynamics, i.e

$\textit{Chromosome
fitness}=\max_{t}\left[\mathcal{F}(\rho_{Q}(t),\sigma)\right].$ (10)

In fact, the function actually used is a slight variation of that presented
above necessitated by the specific details of the simulation and the ability
to “extract” the state with the highest fidelity (see Appendix C for details).

##  V Results

###  V.1 Analysis of performance: case studies

Below we outline the optimal control schemes found when applying the CGA
approach to prepare genuinely entangled 3 and 4 -qubit states from completely
separable initial states. In doing so this acts as a proof of principle — for
both the Hamiltonian $\tilde{H}$, and the optimization method — with respect
to entanglement generation in general and state preparation specifically.
Below we consider the physically reasonable maximum coupling and driving
$g_{0},\xi_{0}=2\pi~{}200$MHz, and total time-scales of the order of $\tau
T\approx 10$ns (see appendix B for full details). We specifically consider 3
states: GHZ state, three-qubit Dicke states with two excitation, and the 4
qubit “Box” Cluster state.

#### GHZ state.

GHZ states are relevant genuinely tripartite entangled states Dür  _et al._
(2000) defined by

$|GHZ\rangle=\frac{|000\rangle+|111\rangle}{\sqrt{2}},$ (11)

We assume the system to initially be in the state
$|\psi_{0}\rangle=|010\rangle_{Q}|0\rangle_{cav}$. We use $|010\rangle$ as the
initial state of the qubit network in this specific instance, as opposed to
simply the global vacuum, since the latter has non-zero fidelity to the target
GHZ state and thus starts the optimization in at undesirable local maxima.
Such initial preparation is not so restrictive given that single qubit
rotations are easily implemented in many quantum systems when compared with
multi-qubit operations. We use values of $\tau=1ns$ for the duration of each
time interval and a total number of intervals afforded to the optimizing agent
of $T=10$. The results are shown in Fig. 2 where a highest fidelity of
$\mathcal{F}(\rho_{Q},\sigma)=0.9746$ is achieved in $\approx 8$ns.

Figure 2: Results for three-qubit GHZ state. We have the optimal control
pulses [Top], the fidelity in the qubit subspace during the dynamics [Middle]
and the matrix histogram for the target state (left) and the state with the
highest fidelity during the dynamics (right) [Bottom]. Hyper-parameters
$\tau=1$ns and $T=10$ are used here, with maximal fidelity achieved, and
maintained, within $\approx 8$ns. On the fidelity plot, the horizontal green
line is drawn at $c_{n}=1/2$ and all fidelities above this exhibit GME (green
region). The blue region $c_{n}=3/4$ on the other hand shows those fidelities
that exhibit both GME and GHZ-class entanglement in particular

#### Three-qubit Dicke state.

Dicke states embody another class of genuinely tripartite entangled states,
inequivalent to GHZ states Dür  _et al._ (2000), which have been
experimentally realised, projected onto lower dimensional entangled states and
employed in open destination teleportation and telecloning Prevedel  _et al._
(2009); Chiuri  _et al._ (2012). Specifically we consider the three-qubit two-
excitation Dicke state defined as

$|{D}^{(2)}_{3}\rangle=\frac{\left(|011\rangle+|101\rangle+|011\rangle\right)}{\sqrt{3}}.$
(12)

We assume initial state preparation of $|000\rangle_{Q}|0\rangle_{cav}$, and
allow more fine control this time with $\tau=0.5ns$ and $T=20$ The maximum
fidelity achieved in this instance is
$\mathcal{F(\rho_{Q},\sigma)}_{max}=0.9896$ within $\approx 7$ns, as shown in
figure 3.

Figure 3: Results for three-qubit Dicke state. We have the optimal control
pulses [Top], the fidelity in the qubit subspace during the dynamics [Middle]
and the matrix histogram for the target state (left) and the state with the
highest fidelity during the dynamics (right) [Bottom]. The hyper-parameters
$\tau=0.5$ns and $T=20$ have been used for these calculations. The
corresponding maximal fidelity was achieved within $\approx 7$ns. The green
region on the fidelity plot highlights fidelities for which genuine
multipartite entanglement (GME) is detected based on the fidelity based
witness Eq. (15) with $c_{n}=2/3$. The inset shows
$\text{Tr}(\rho\mathcal{\hat{W}})$ at each instant of time, and the grey
region on both the inset plot and the fidelity plot show the temporal region
where GME is detected via the collective-spin witness in Eq. (16) with
$b_{s}=3.12$ Tóth (2007).

#### Box Cluster State.

Cluster states are a class of graph states useful for measurement based
quantum computing Nielsen (2006); Bartlett and Rudolph (2006); Walther  _et
al._ (2005); Briegel  _et al._ (2009). Cluster states are represented by
graphs of interconnected vertices, where if one starts with each qubit in the
ground state, they are defined procedurally by applying a Hadamard gate to
each qubit (represented by the vertices of the graph) then applying
conditional-phase gates between qubits whose vertices are connected by edges.
Here we consider a so-called Box cluster state which is defined by four
vertices connected by four edges to form a square. This can be written as,

$|\psi_{Box}\rangle=CZ_{41}\left(\Pi^{3}_{j=1}CZ_{i,i+1}\right)\left(\Pi^{4}_{j=1}{\cal
H}_{j}\right)|0000\rangle_{1234},$ (13)

where $CZ_{i,j}$ is the controlled-Z gate between qubits $i$ and $j$ and
${\cal H}_{i}$ is the Hadamard gate on qubit $i$. Explicitly, the state reads:

$|\psi_{Box}\rangle=\sum_{i,j,k,l=0,1}(-1)^{x_{i}+x_{j}+x_{k}+x_{l}}|x_{i}x_{j}x_{k}x_{l}\rangle.$
(14)

As before, assuming $\tau=0.5ns$ and $T=20$ and starting from the initial
vacuum state a maximum fidelity of $\mathcal{F}(\rho_{Q},\sigma)=0.9642$ was
achieved (cf. Fig. 4), this time requiring the full $10ns$ to achieve and
maintain maximal fidelity.

Figure 4: Results for a four-qubit Box Cluster state. We have plotted the
optimal control pulses [Top panel], the fidelity in the qubit subspace during
the dynamics [Middle panel] and the matrix histogram for both the target state
(left-most figure) and the state with the highest fidelity during the dynamics
(right-most one) [Bottom panel]. We have used the hyper-parameters
$\tau=0.5$ns and $T=20$ and the maximum fidelity was achieved within $\approx
10$ns. The values of fidelity for which GME is detected via the fidelity based
witness Eq. (15) with $c_{n}=1/2$ Gühne and Tóth (2009) are shown in the green
region.

###  V.2 Entanglement Detection

A central question in any state-engineering scheme is the characterization of
the features of the state that has been synthesized. Within the context of our
investigation, the core aspect to address is the quantification of
multipartite entanglement.

The task of accurately determining the amount of entanglement in a given
quantum state, which is challenging in general, is made even more difficult in
multipartite settings due to the hierarchical structure of entanglement in
many-body systems Horodecki  _et al._ (2009); Plenio and Virmani (2014);
Christandl and Winter (2004); Guo and Zhang (2020); Szalay (2015) and the need
to construct convex-roof extensions of any pure-state quantifier, when dealing
with mixed states Gao  _et al._ (2008); Cao  _et al._ (2010); Röthlisberger
_et al._ (2009).

A significant tool in these endeavours is embodied by entanglement witnesses,
which often offer experimentally viable ways of detecting (or even quantifying
Brandão (2005)) entanglement Chruściński and Sarbicki (2014); Gühne and Tóth
(2009) that are suitable for mixed states and have been used to detect genuine
multipartite entanglement (GME) already for GHZ class, W-class and graph
states Tóth (2007); Campbell  _et al._ (2009); Jungnitsch  _et al._ (2011).
Whilst high fidelity with a maximally entangled state is a good indicator that
we have generated entanglement close to the right structure, it is useful to
quantitatively check for GME in the system. A natural approach is to use so-
called “fidelity based” entanglement witnesses. These are of the general form

$\mathcal{\hat{W}}_{F}=c_{n}\mathds{1}-|\psi\rangle\langle\psi|,$ (15)

where $|\psi\rangle$ is the state of interest and $c_{n}$ is the maximal
overlap between $|\psi\rangle$ and all bi-separable states. Thus, any state
for which Tr$(\rho\hat{\mathcal{W}}_{F})\geq 0$ is bi-separable and
consequently Tr$(\rho\hat{\mathcal{W}}_{F})<0$ indicates genuine multipartite
entanglement. The overlap values $c_{n}$ have already been calculated for
3-qubit GHZ and W states, and 4-qubit linear cluster states Tóth (2007); Gühne
and Tóth (2009), which up to local unitaries and swaps coincide exactly with
the three cases considered above. In light of the analysis reported in Figs.
2, 3, and 4 these witnesses are readily implemented. By considering the
quantity $-\text{Tr}(\rho\mathcal{\hat{W}}_{F})+c_{n}$, we can highlight a
threshold value of fidelity above which GME can be detected. The region where
this happens is highlighted in green in each of the figures.

One may also be interested in not only ensuring that the state has GME but
also, in the GHZ case, if we have generated GHZ-class entanglement Acín  _et
al._ (2001). In this case, $c_{n}$ will be the maximal overlap between
$|GHZ\rangle$ and all W-class entangled states. This region is highlighted in
blue in Fig. 2.

Finally, other constructions of witnesses are useful for specific states; for
example, witnesses based on collective spin operators have been used to more
efficiently detect GME in symmetric Dicke states as in Tóth (2007). Here the
witness takes the form

$\hat{\mathcal{W}}=b_{s}\mathds{1}-(\hat{J}_{x}^{2}+\hat{J}_{y}^{2}),$ (16)

where $\hat{J}_{k}$ is the collective $k=x,y$ spin operator Campbell  _et al._
(2009). It is not so straightforward here to place a delineation at a given
fidelity as with the fidelity based witnesses, so we first plot
$\text{Tr}(\rho\mathcal{\hat{W}})$ as an inset within figure 3 and highlight
the region for which $\text{Tr}(\rho\mathcal{\hat{W}})<0$. This region is then
shown as the grey hatched area on the larger fidelity plot. It can be seen
that both the fidelity based witness and the collective spin based witness
detect GME at strikingly similar times.

###  V.3 Decoherence

So far, we have exclusively considered closed system dynamics and as such the
optimality of the control schemes presented is limited to the noiseless case.
It is of interest then to assess how these control schemes perform in the
presence of decoherence. Specifically we can write the following Lindblad
master equation to model the effect of decay and dephasing acting on each of
the constituent subsystems Breuer  _et al._ (2002); Manzano (2020)

$\dot{\rho}=-i\left[\tilde{H},\rho\right]+\kappa\mathcal{D}\left[a\right]\rho+\sum_{j}^{N}\left(\gamma_{j}\mathcal{D}\left[\sigma^{-}_{j}\right]\rho+2\gamma_{\phi,j}\mathcal{D}\left[\sigma^{+}_{j}\sigma^{-}_{j}\right]\rho\right),$
(17)

where $\kappa$ is the cavity decay rate, $\gamma_{\phi,j}$ and $\gamma_{j}$
are the dephasing and decay rate for the $j^{th}$ qubit, respectively, and we
have introduced the superoperators

${\cal D}[Q]\rho=Q\rho
Q^{\dagger}-\frac{1}{2}\left\\{Q^{\dagger}Q,\rho\right\\}$ (18)

for an arbitrary operator $Q$. If we adopt physically reasonable values for
these rates we can obtain an estimate of the performance of a physical system.
For example, considering superconducting systems we set $\kappa=2\pi\times
5\text{ kHz}$ for the cavity damping and the typical values of $2\pi\times
300\text{ KHz}$ and $2\pi\times 5\text{ MHz}$ for the dephasing and decay rate
for each of the qubits Majer  _et al._ (2007); Blais  _et al._ (2021). In
these systems, the use of high-Q cavities and low temperatures leads to a
reduced cavity decay and dephasing rate. The qubit decay is thus the main
source of decoherence. In these conditions, the effects of noise is reported
in Fig. 5. We can clearly see that the control protocols are almost completely
insensitive to cavity decay and qubit dephasing, whilst still reasonably
robust against qubit decay.

Figure 5: A comparison between the fidelity achieved for the (Top) GHZ state
preparation, (Middle) Dicke state preparation and (Bottom) Box Cluster state
preparation schemes in the presence of cavity decay and qubit dephasing only
(Grey) and qubit decay only (Red). Each plot shows the maximum achieved
fidelity in the latter case. The horizontal blue line shows the maximum
fidelity achieved in the ideal, noiseless case.

##  VI Conclusions

We have investigated the use of evolutionary algorithms, which are well-
established classical machine learning strategies, for quantum state
engineering and multipartite entanglement generation within the framework of
quantum control. Specifically, we considered the effective Hamiltonian of a
network of non-interacting qubits jointly addressed by a common driving bus
and applied a genetic algorithm to identify the set of optimal pulses to drive
the evolution of the qubit register. This has allowed us to successfully put
forward robust protocols for the engineering of three- and four-qubit states
that play a crucial role in quantum metrology and computing, including Dicke
and cluster states. The protocols, which offer significant robustness to the
most common and crucial sources of imperfection, provide evidence of the
benefit of a hybrid approach to quantum control that puts together the insight
provided by machine learning strategies to well-established schemes for
optimal control. The extension of our approach to larger registers and non-
unitary dynamics will pave the way to quantum process engineering enhanced by
machine learning and optimised by quantum control methods.

###### Acknowledgements.

We acknowledge support from the European Union’s Horizon 2020 FET-Open project
TEQ (766900), the Leverhulme Trust Research Project Grant UltraQuTe (grant
RGP-2018-266), the Royal Society Wolfson Fellowship (RSWF/R3/183013), the UK
EPSRC (grant EP/T028424/1) and the Department for the Economy Northern Ireland
under the US-Ireland R&D Partnership Programme.

## References

  * Georgescu  _et al._ (2014) Iulia M Georgescu, Sahel Ashhab, and Franco Nori, “Quantum simulation,” Rev. Mod. Phys. 86, 153 (2014).
  * Aaronson (2015) Scott Aaronson, “Read the fine print,” Nat. Phys. 11, 291–293 (2015).
  * Muralidharan  _et al._ (2016) Sreraman Muralidharan, Linshu Li, Jungsang Kim, Norbert Lütkenhaus, Mikhail D Lukin, and Liang Jiang, “Optimal architectures for long distance quantum communication,” Sci. Rep. 6, 1–10 (2016).
  * Horodecki  _et al._ (2009) Ryszard Horodecki, Paweł Horodecki, Michał Horodecki, and Karol Horodecki, “Quantum entanglement,” Rev. Mod. Phys. 81, 865 (2009).
  * D’Alessandro (2007) D. D’Alessandro,  _Introduction to Quantum Control and Dynamics_, Chapman & Hall/CRC Applied Mathematics & Nonlinear Science (Taylor & Francis, 2007).
  * Glaser  _et al._ (2015) Steffen J. Glaser, Ugo Boscain, Tommaso Calarco, Christiane P. Koch, Walter Köckenberger, Ronnie Kosloff, Ilya Kuprov, Burkhard Luy, Sophie Schirmer, Thomas Schulte-Herbrüggen, Dominique Sugny, and Frank.K. Wilhelm, “Training schrödinger’s cat: quantum optimal control,” Eur. Phys. J. D 69, 279 (2015).
  * Clarke and Wilhelm (2008) John Clarke and Frank Wilhelm, “Superconducting quantum bits,” Nature 453, 1031–42 (2008).
  * Kjaergaard  _et al._ (2020) Morten Kjaergaard, Mollie E Schwartz, Jochen Braumüller, Philip Krantz, Joel I-J Wang, Simon Gustavsson, and William D Oliver, “Superconducting qubits: Current state of play,” Annu. Rev. Condens. Matter Phys. 11, 369–395 (2020).
  * Huang  _et al._ (2020) He-Liang Huang, Dachao Wu, Daojin Fan, and Xiaobo Zhu, “Superconducting quantum computing: a review,” Science China Information Sciences 63, 1–32 (2020).
  * Koch  _et al._ (2007) Jens Koch, Terri M. Yu, Jay Gambetta, A. A. Houck, D. I. Schuster, J. Majer, Alexandre Blais, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, “Charge-insensitive qubit design derived from the cooper pair box,” Phys. Rev. A 76, 042319 (2007).
  * DiCarlo  _et al._ (2009) Leonardo DiCarlo, Jerry M Chow, Jay M Gambetta, Lev S Bishop, Blake R Johnson, DI Schuster, J Majer, Alexandre Blais, Luigi Frunzio, SM Girvin, _et al._ , “Demonstration of two-qubit algorithms with a superconducting quantum processor,” Nature 460, 240–244 (2009).
  * Barends  _et al._ (2014) Rami Barends, J Kelly, A Megrant, A Veitia, D Sank, E Jeffrey, TC White, J Mutus, AG Fowler, B Campbell, _et al._ , “Logic gates at the surface code threshold: Superconducting qubits poised for fault-tolerant quantum computing,” Nature 508, 500 (2014).
  * Li  _et al._ (2019) Shaowei Li, Anthony D. Castellano, Shiyu Wang, Yulin Wu, Ming Gong, Zhiguang Yan, Hao Rong, Hui Deng, Chen Zha, Cheng Guo, Lihua Sun, Chengzhi Peng, Xiaobo Zhu, and Jian-Wei Pan, “Realisation of high-fidelity nonadiabatic CZ gates with superconducting qubits,” npj Quant. Inf. 5, 1–7 (2019).
  * Barends  _et al._ (2019) Rami Barends, CM Quintana, AG Petukhov, Yu Chen, Dvir Kafri, Kostyantyn Kechedzhi, Roberto Collins, Ofer Naaman, Sergio Boixo, F Arute, _et al._ , “Diabatic gates for frequency-tunable superconducting qubits,” Phys. Rev. Lett. 123, 210501 (2019).
  * Chow  _et al._ (2011) Jerry M Chow, Antonio D Córcoles, Jay M Gambetta, Chad Rigetti, Blake R Johnson, John A Smolin, Jim R Rozen, George A Keefe, Mary B Rothwell, Mark B Ketchen, _et al._ , “Simple all-microwave entangling gate for fixed-frequency superconducting qubits,” Phys. Rev. Lett. 107, 080502 (2011).
  * Córcoles  _et al._ (2013) Antonio D Córcoles, Jay M Gambetta, Jerry M Chow, John A Smolin, Matthew Ware, Joel Strand, Britton LT Plourde, and Matthias Steffen, “Process verification of two-qubit quantum gates by randomized benchmarking,” Phys. Rev. A 87, 030301(R) (2013).
  * Majer  _et al._ (2007) Johannes Majer, J M Chow, J M Gambetta, Jens Koch, B R Johnson, J A Schreier, L Frunzio, D I Schuster, A A Houck, A Wallraff, _et al._ , “Coupling superconducting qubits via a cavity bus,” Nature 449, 443–447 (2007).
  * Cross and Gambetta (2015) Andrew W Cross and Jay M Gambetta, “Optimized pulse shapes for a resonator-induced phase gate,” Phys. Rev. A 91, 032325 (2015).
  * Puri and Blais (2016) Shruti Puri and Alexandre Blais, “High-fidelity resonator-induced phase gate with single-mode squeezing,” Phys. Rev. Lett. 116, 180501 (2016).
  * Paik  _et al._ (2016) Hanhee Paik, A Mezzacapo, Martin Sandberg, DT McClure, B Abdo, AD Córcoles, O Dial, DF Bogorin, BLT Plourde, M Steffen, _et al._ , “Experimental demonstration of a resonator-induced phase gate in a multiqubit circuit-qed system,” Phys. Rev. Lett. 117, 250502 (2016).
  * Blais  _et al._ (2007) Alexandre Blais, Jay Gambetta, Andreas Wallraff, David I Schuster, Steven M Girvin, Michel H Devoret, and Robert J Schoelkopf, “Quantum-information processing with circuit quantum electrodynamics,” Phys. Rev. A 75, 032329 (2007).
  * Blais  _et al._ (2021) Alexandre Blais, Arne L. Grimsmo, S. M. Girvin, and Andreas Wallraff, “Circuit quantum electrodynamics,” Rev. Mod. Phys. 93, 025005 (2021).
  * Falci  _et al._ (2017) G Falci, P G Di Stefano, A Ridolfo, A D’Arrigo, G S Paraoanu, and E Paladino, “Advances in quantum control of three-level superconducting circuit architectures,” Fortschr. Phys. 65, 1600077 (2017).
  * Di Stefano  _et al._ (2015) P. G. Di Stefano, E. Paladino, A. D’Arrigo, and G. Falci, “Population transfer in a lambda system induced by detunings,” Phys. Rev. B 91 (2015).
  * Sutton and Barto (2018) Richard S. Sutton and Andrew G. Barto,  _Reinforcement Learning: An Introduction_, 2nd ed. (The MIT Press, 2018).
  * Mnih  _et al._ (2013) Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Alex Graves, Ioannis Antonoglou, Daan Wierstra, and Martin Riedmiller, “Playing atari with deep reinforcement learning,” (2013), arXiv:1312.5602 [cs.LG] .
  * Giannelli  _et al._ (2022) Luigi Giannelli, Pierpaolo Sgroi, Jonathon Brown, Gheorghe Sorin Paraoanu, Mauro Paternostro, Elisabetta Paladino, and Giuseppe Falci, “A tutorial on optimal control and reinforcement learning methods for quantum technologies,” Phys. Lett. A , 128054 (2022).
  * Paparelle  _et al._ (2020) Iris Paparelle, Lorenzo Moro, and Enrico Prati, “Digitally stimulated Raman passage by deep reinforcement learning,” Phys. Lett. A 384, 126266 (2020).
  * Porotti  _et al._ (2019) Riccardo Porotti, Dario Tamascelli, Marcello Restelli, and Enrico Prati, “Coherent transport of quantum states by deep reinforcement learning,” Commun. Phys. 2, 61 (2019).
  * Sgroi  _et al._ (2021) Pierpaolo Sgroi, G. Massimo Palma, and Mauro Paternostro, “Reinforcement learning approach to non-equilibrium quantum thermodynamics,” Phys. Rev. Lett. 126, 020601 (2021).
  * Brown  _et al._ (2021) Jonathon Brown, Pierpaolo Sgroi, Luigi Giannelli, Gheorghe Sorin Paraoanu, Elisabetta Paladino, Giuseppe Falci, Mauro Paternostro, and Alessandro Ferraro, “Reinforcement learning-enhanced protocols for coherent population-transfer in three-level quantum systems,” New J. Phys. 23, 093035 (2021).
  * Kuo  _et al._ (2021) En-Jui Kuo, Yao-Lung L Fang, and Samuel Yen-Chi Chen, “Quantum architecture search via deep reinforcement learning,” arXiv:2104.07715 (2021).
  * An  _et al._ (2021) Zheng An, Hai-Jing Song, Qi-Kai He, and DL Zhou, “Quantum optimal control of multilevel dissipative quantum systems with reinforcement learning,” Phys. Rev. A 103, 012404 (2021).
  * Zhang  _et al._ (2019) Xiao-Ming Zhang, Zezhu Wei, Raza Asad, Xu-Chen Yang, and Xin Wang, “When does reinforcement learning stand out in quantum control? A comparative study on state preparation,” npj Quant. Inf. 5, 1–7 (2019).
  * Note (1) Noticeable exceptions are discussed in Refs. Sgroi  _et al._ (2021); Brown  _et al._ (2021).
  * Salimans  _et al._ (2017) Tim Salimans, Jonathan Ho, Xi Chen, Szymon Sidor, and Ilya Sutskever, “Evolution strategies as a scalable alternative to reinforcement learning,” arXiv:1703.03864 (2017).
  * Such  _et al._ (2017) Felipe Petroski Such, Vashisht Madhavan, Edoardo Conti, Joel Lehman, Kenneth O Stanley, and Jeff Clune, “Deep neuroevolution: Genetic algorithms are a competitive alternative for training deep neural networks for reinforcement learning,” arXiv:1712.06567 (2017).
  * Ghosh  _et al._ (2021) Sanjib Ghosh, Tanjung Krisnanda, Tomasz Paterek, and Timothy CH Liew, “Realising and compressing quantum circuits with quantum reservoir computing,” Commun. Phys. 4, 1–7 (2021).
  * Haack  _et al._ (2010) Géraldine Haack, Ferdinand Helmer, Matteo Mariantoni, Florian Marquardt, and Enrique Solano, “Resonant quantum gates in circuit quantum electrodynamics,” Phys. Rev. B 82, 024514 (2010).
  * Note (2) An interesting consequence of the particular form of $\tilde{H}$ in Eq. (5) is that, if the system is initialized in a pure state with real coefficients, then the time-evolved reduced state of the qubit will always be described by real coefficients. This is of course a restriction if one wishes to prepare more general states, however it simplifies the optimization when considering completely real target states. One can ease this constraint if necessary by introducing non-zero detunings in Eq. (4).
  * Haupt and Haupt (2004) Randy L Haupt and Sue Ellen Haupt,  _Practical genetic algorithms_ (John Wiley & Sons, 2004).
  * Dür  _et al._ (2000) Wolfgang Dür, Guifre Vidal, and J Ignacio Cirac, “Three qubits can be entangled in two inequivalent ways,” Phys. Rev. A 62, 062314 (2000).
  * Prevedel  _et al._ (2009) Robert Prevedel, Gunther Cronenberg, Mark S Tame, Mauro Paternostro, Philip Walther, Mu-Seong Kim, and Anton Zeilinger, “Experimental realization of Dicke states of up to six qubits for multiparty quantum networking,” Phys. Rev. Lett. 103, 020503 (2009).
  * Chiuri  _et al._ (2012) A Chiuri, C Greganti, M Paternostro, G Vallone, and P Mataloni, “Experimental quantum networking protocols via four-qubit hyperentangled Dicke states,” Phys. Rev. Lett. 109, 173604 (2012).
  * Tóth (2007) Géza Tóth, “Detection of multipartite entanglement in the vicinity of symmetric dicke states,” JOSA B 24, 275–282 (2007).
  * Nielsen (2006) Michael A Nielsen, “Cluster-state quantum computation,” Rep. Math. Phys. 57, 147–161 (2006).
  * Bartlett and Rudolph (2006) Stephen D Bartlett and Terry Rudolph, “Simple nearest-neighbor two-body hamiltonian system for which the ground state is a universal resource for quantum computation,” Phys. Rev. A 74, 040302 (2006).
  * Walther  _et al._ (2005) Philip Walther, Kevin J Resch, Terry Rudolph, Emmanuel Schenck, Harald Weinfurter, Vlatko Vedral, Markus Aspelmeyer, and Anton Zeilinger, “Experimental one-way quantum computing,” Nature 434, 169–176 (2005).
  * Briegel  _et al._ (2009) Hans J Briegel, David E Browne, Wolfgang Dür, Robert Raussendorf, and Maarten Van den Nest, “Measurement-based quantum computation,” Nat. Phys. 5, 19–26 (2009).
  * Gühne and Tóth (2009) Otfried Gühne and Géza Tóth, “Entanglement detection,” Physics Reports 474, 1–75 (2009).
  * Plenio and Virmani (2014) Martin B Plenio and Shashank S Virmani, “An introduction to entanglement theory,” in  _Quantum information and coherence_ (Springer, 2014) pp. 173–209.
  * Christandl and Winter (2004) Matthias Christandl and Andreas Winter, ““Squashed entanglement”: an additive entanglement measure,” J. Math. Phys. 45, 829–840 (2004).
  * Guo and Zhang (2020) Yu Guo and Lin Zhang, “Multipartite entanglement measure and complete monogamy relation,” Phys. Rev. A 101, 032301 (2020).
  * Szalay (2015) Szilárd Szalay, “Multipartite entanglement measures,” Phys. Rev. A 92, 042329 (2015).
  * Gao  _et al._ (2008) Xiuhong Gao, Sergio Albeverio, Kai Chen, Shao-Ming Fei, and Xianqing Li-Jost, “Entanglement of formation and concurrence for mixed states,” Frontiers of Computer Science in China 2, 114–128 (2008).
  * Cao  _et al._ (2010) Kun Cao, Zheng-Wei Zhou, Guang-Can Guo, and Lixin He, “Efficient numerical method to calculate the three-tangle of mixed states,” Phys. Rev. A 81, 034302 (2010).
  * Röthlisberger  _et al._ (2009) Beat Röthlisberger, Jörg Lehmann, and Daniel Loss, “Numerical evaluation of convex-roof entanglement measures with applications to spin rings,” Phys. Rev. A 80, 042301 (2009).
  * Brandão (2005) Fernando G. S. L. Brandão, “Quantifying entanglement with witness operators,” Phys. Rev. A 72, 022310 (2005).
  * Chruściński and Sarbicki (2014) Dariusz Chruściński and Gniewomir Sarbicki, “Entanglement witnesses: construction, analysis and classification,” J. Phys. A: Math. Theor. 47, 483001 (2014).
  * Campbell  _et al._ (2009) Steve Campbell, M S Tame, and Mauro Paternostro, “Characterizing multipartite symmetric Dicke states under the effects of noise,” New J. Phys. 11, 073039 (2009).
  * Jungnitsch  _et al._ (2011) Bastian Jungnitsch, Tobias Moroder, and Otfried Gühne, “Entanglement witnesses for graph states: General theory and examples,” Phys, Rev. A 84, 032310 (2011).
  * Acín  _et al._ (2001) Antonio Acín, Dagmar Bruß, Maciej Lewenstein, and Anna Sanpera, “Classification of mixed three-qubit states,” Phys. Rev. Lett. 87, 040401 (2001).
  * Breuer  _et al._ (2002) H.P. Breuer, P.I.H.P. Breuer, F. Petruccione, and S.P.A.P.F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2002).
  * Manzano (2020) Daniel Manzano, “A short introduction to the lindblad master equation,” AIP Advances 10, 025106 (2020).
  * Johansson  _et al._ (2012) J Robert Johansson, Paul D Nation, and Franco Nori, “QuTiP: An open-source python framework for the dynamics of open quantum systems,” Comput. Phys. Commun. 183, 1760–1772 (2012).

##  Appendix A Functional Form for controls

A common approach to optimal control problems is to break the dynamics up into
$T$ time intervals of equal duration $\tau$, then at each interval assign a
constant value to each control parameter. This results in piece-wise constant
(PWC) control “pulses”, $u(t)\to\\{u_{i}\\}_{i=0,...T-1}$, where $u$ is
routinely used to denote a generic control function. Such discretization is
useful for application of Reinforcement Learning techniques as in Porotti  _et
al._ (2019); Paparelle  _et al._ (2020); Sgroi  _et al._ (2021); Brown  _et
al._ (2021); Giannelli  _et al._ (2022). Whilst useful, this type of
functional form includes discontinuities, often in the form of instantaneous
jumps in the control, which is experimentally unfeasible and subsequently
requires some method of smoothing post-optimization often to the detriment of
performance. Here we use a generalised form of this discretisation, where
instead of allowing the optimization to chose the $T$ constant values
corresponding to the control value during each time interval, we allow it to
chose the $T+1$ values corresponding to the start(/end) of each interval, then
connect each value with a smooth, time dependent function during the interval.
Namely, the functional time dependence of each control pulse is made to be a
clipped $\tanh$ function centred in the middle of the time interval. For
example consider a simple $\tanh$ function, shown in figure A. If we clip the
$\tanh$ function within windows of different widths centred around zero we can
make a “step-like” like function with varying severity. This is tantamount to
scaling the $\tanh$ function along the y-axis and clipping it at $\pm 1$. So
the time dependence within time interval $i$, $t_{i}\leq t\leq t_{i+1}$ can be
written as,

$f_{i}(t)=\begin{cases}0&\mbox{if }t<t_{i}\\\
(u_{i+1}-u_{i})\left(\frac{S\tanh(Wt-(t_{i}+\frac{\tau}{2}))+1}{2}\right)+u_{i}&\mbox{if
}t_{i}\leq t<t_{i+1}\\\ 0&\mbox{if }t\geq t_{i+1}\end{cases}$ (19)

where $u_{i}$, $u_{i+1}$ is the value of the control at times $t_{i}$,
$t_{i+1}$ respectively, $W$ determines the severity of the step, and $S$ is a
scaling factor introduced to deal with the error $\varepsilon$, as in fig 6.

Figure 6: Comparison of how the clipping window affects the shape of the
interconnecting tanh functions used to construct the control functions. In (a)
the tanh function is shown with over-laid dashed lines representing different
widths of this clipping window. From (a) we can see how the clipping window
width determines the severity of the time dependence. Namely, taking the
largest width (grey) leads to a more “step-like” dependence, as evidenced in
(b) where each of the clipped functions are scaled into a time interval of
equal duration $\tau$ and re-scaled to account for the error $\varepsilon$.
(a) shows how the error $\varepsilon$ increases as the clipping window
narrows, i.e the narrowest (yellow) window has the largest error thus
requiring the most re-scaling.

Therefore the complete functional for a general control, under this scheme, is
given by

$f(t)=\sum_{i=0}^{T}f_{i}(t),$ (20)

where, again, $T$ is the number of time intervals. Full functional forms can
be clearly seen in figures 2, 3, and 4.

##  Appendix B Algorithm Implementation Details.

In each case the population size for the algorithm was $\mathcal{N}_{pop}=48$
and was determined by the number of CPUs available, since each chromosome in
the population was evaluated in parallel to significantly speed up computation
time. The number of survivors was fixed at $\mathcal{N}_{pop}/2=24$, from
which $\mathcal{N}_{pop}/4=12$ pairs of parents were selected according to a
probability distribution determined by their relative fitness. Each pair of
parent chromosomes produce 2 offspring chromosome to re-populate. The mating
procedure employed involved 2 steps, after making a one copy of each parent
chromosomes: 1) Each separate section of $T+1$ elements - corresponding to the
different control pulses - were completely swapped between the two copied
chromosomes with $\approx 50\%$ probability, otherwise they were left
unchanged. 2) random indices were then selected and the combination given by
equation 9 was applied, where $\beta$ was randomly sampled from [0,1] for each
separate index, in an inverse manner to the two copies. This results in 2 new
offspring chromosomes that are completely complementary to one and other.
Mutation was then applied by selecting chromosomes (apart from the fittest) at
random and random indices within these chromosomes to replace with completely
random values. The rate $\alpha$ determined the total number of parameter
values within the entire generation that were flipped and generally assumed
values of $\alpha\approx 0.2$.

##  Appendix C Simulation Details

Simulations were carried out using the Numerical Schrödinger Equation solver
within the QuTiP package in python Johansson  _et al._ (2012). Thus, for the
simulations we had to approximate the harmonic mode by a $d$-dimensional
harmonic system, however care must be taken. Since the cavity is assumed to be
resonantly driven, then if we use a low value of $d$ and/or too strong a
drive, then the $d$-dimensional harmonic system will quickly become saturated,
and in fact early optimization used this fact to produce extremely good
results with low dimensional cavities. Of course when one then simulated these
controls with higher levels in the cavity the performance was destroyed and so
the results were neither realistic nor practically useful. One can avoid this
in two ways: 1) By encoding enough redundancy in the cavity by using very high
values of $d$, which increases computational cost; or 2) Including a term in
the numercal fitness function that punishes population of the higher energy
states of the cavity approximation. Here we apply both by chosing $d\in[5,6]$
as well as including the term

$\phi_{1}=-\frac{\nu}{\tau T}\int_{0}^{\tau T}\langle n|\rho_{cav}(t)|n\rangle
dt,$ (21)

where $|n\rangle$ is the highest excited state withing the $d$-dimensional
approximation and $\rho_{cav}$ is the reduced state of the cavity. $\nu$
determines the strength of punishment and was set at $\nu=0.1$. (Practically,
since the dynamics is solved numerically the integral was actually a
summation).

Another issue one encounters with this type of simulation is that, if we
assess fitness based on the outright maximum value of the fidelity during the
induced dynamics, then it is possible to observe successful control schemes
that induce sharp spikes in the fidelity landscape. If the control scheme
induces such spikes on a time scales shorter than that required to completely
“switch off” all of the controls then it is impossible to extract the state of
maximum fidelity and the control sequences again cease to be practically
useful. We can combat this by including an additional term in the fitness
function that rewards control schemes that briefly maintain near maximal
fidelity for a short time, allowing us to selectively uncouple the cavity
whilst maintaining the state of maximal fidelity in the qubit subspace. The
term used was

$\phi_{2}=+\frac{\mu}{m\tau}\int_{t_{max}}^{t_{max}+m\tau}\mathcal{F}(\rho_{Q}(t),\sigma)dt,$
(22)

where $m$ is the number of time intervals of length $\tau$ to include in the
numerical bonus beyond which we no longer care if the fidelity deteriorates.
$\mu$ is again a variable that determines the relative importance of
maintaining fidelity after maximum and was set to $\mu=0.5$. Thus the actual
reward function employed was

$\displaystyle\textit{Chromosome fitness}=$
$\displaystyle\left[\mathcal{F}(\rho_{Q}(t_{max}),\sigma)\right]+\phi_{1}+\phi_{2}$
(23)

where $t_{max}$ is the time at which maximum fidelity is achieved during the
induced dynamics. This leads to to control schemes that maximise and briefly
maintain fidelity allowing us to selective switch of the couplings, whilst
also only exclusively utilising lower lying levels of the harmonic mode. Thus
in principle the resulting controls could be practically implemented and yield
identical performance. Clearly, these considerations are necessitated by the
use of simulation and the case is much simpler if one wishes to use the
algorithm on a physical system, however in this case the ability to
parallelise the computational steps, one major advantage of the Genetic
Algorithm, is suppressed.