$F^{\prime}:={\left.\kern-1.2ptf\vphantom{\big{|}}\right|_{l^{\prime}}}\rightarrow\mathbb{R}$.
Up to a reordering of the simplices of
$S=\left\\{\sigma_{1},\dots,\sigma_{N}\right\\}$, we assume without loss of
generality that $F(\sigma_{1})\leq\dots\leq F(\sigma_{N})$. Let
$\mathfrak{S}_{N}^{S}\subseteq\mathfrak{S}_{N}$ be the subset of permutations
that satisfies
$\tau\in\mathfrak{S}_{N}^{S}\Longrightarrow\left(C^{\tau}_{k}\right)_{k\in\llbracket
1,N\rrbracket}=\left(\left\\{\sigma_{\tau(1)},\dots,\sigma_{\tau(k)}\right\\}\right)_{k\in\llbracket
1,N\rrbracket}$ is a filtration of $S$. Finally, let
$N=(N_{\tau,k})_{(\tau,k)\in\mathfrak{S}_{N}^{S}\times\llbracket
1,N\rrbracket}$, where $N_{\tau,k}=H_{*}(C^{\tau}_{k})$ be the
multipersistence module indexed over $\mathfrak{S}_{N}^{S}\times\llbracket
1,N\rrbracket$, with arrows induced by inclusion. Note that the one-
dimensional persistence module $N_{\rm id}:=(N_{{\rm id},k})_{k\in\llbracket
1,N\rrbracket}$ is isomorphic to
${\left.\kern-1.2ptM\vphantom{\big{|}}\right|_{l}}$.

Now, let $\tau$ be the simplex permutation that matches the simplex ordering
induced by $F$ to the one induced by $F^{\prime}$. By definition, one has
$\tau\in\mathfrak{S}_{N}^{S}$. Let $i$ be the index of the first transposition
in the Coxeter decomposition of $\tau$. Then, the matching ${\rm vine}$
associated to the transposition $(i\,i+1)$ induces a matching between the bars
of ${\mathcal{B}}({\left.\kern-1.2ptM\vphantom{\big{|}}\right|_{l}})$ and
${\mathcal{B}}({\left.\kern-1.2ptM\vphantom{\big{|}}\right|_{l^{\prime}}})$,
or equivalently, between the bars of ${\mathcal{B}}(N_{\rm id})$ and
${\mathcal{B}}(N_{(i\,i+1)})$, as well as a morphism between $N_{\rm id}$ and
$N_{(i\,i+1)}$. Let us now consider the diamond:

${N_{(i\,i+1),i}}$${\cdots}$${N_{\text{id},i-1}=N_{(i\,i+1),i-1}}$${N_{\text{id},i+1}=N_{(i\,i+1),i+1}}$${\cdots}$${N_{\text{id},i}}$
(9)

First, note that all maps in that diamond are either injective of corank 1 or
surjective of nullity 1 (since they are all induced by adding a positive or
negative chain of the corresponding reduced bases). Moreover, by the Mayer-
Vietoris theorem, the following sequence is exact:

$N_{\text{id},i-1}=N_{(i\,i+1),i-1}=H_{*}(C^{\text{id}}_{i-1})\rightarrow
N_{(i\,i+1),i}\oplus N_{\text{id},i}=H_{*}(C^{(i\,i+1)}_{i})\oplus
H_{*}(C^{\text{id}}_{i})\xrightarrow{(x,y)\mapsto
y-x}H_{*}(C^{\text{id}}_{i+1})=N_{\text{id},i+1}=N_{(i\,i+1),i+1}.$

Such diamonds are called transposition diamonds, and it has been shown in
[MO15, Theorem 2.4] that the morphism induced by $\mathrm{vine}$ between the
lower and upper parts of the diamond matches bars with same representative
positive chains together. Hence, bars in ${\mathcal{B}}(N_{\rm id})$ and
${\mathcal{B}}(N_{(i\,i+1)})$ that are matched under $\mathrm{vine}$ are
associated to the same summand of $N$. Moreover, by repeating this
argumentation with the other transpositions in the decomposition of $\tau$,
one has that the same is true for bars in
${\mathcal{B}}({\left.\kern-1.2ptM\vphantom{\big{|}}\right|_{l}})$ and
${\mathcal{B}}({\left.\kern-1.2ptM\vphantom{\big{|}}\right|_{l^{\prime}}})$.

Now, if $N$ is interval decomposable, then, since the transition maps of $M$
can be seen as transition maps of $N$, it follows that interval summands of
$M$ correspond to interval summands of $N$. In that case, using a
dimensionality argument with Theorem 6.6, one can show that two bars of
${\mathcal{B}}({\left.\kern-1.2ptM\vphantom{\big{|}}\right|_{l}})$ and
${\mathcal{B}}({\left.\kern-1.2ptM\vphantom{\big{|}}\right|_{l^{\prime}}})$
that are matched under $\mathrm{vine}$ through arrows of $N$ also belong to
the same interval summand of $M$. Unfortunately, even though $N$ and $M$ are
constructed from the same chain complex, $N$ contains much more arrows that
$M$, and we cannot guarantee in the general case that $N$ is interval
decomposable.

We will thus conclude with Conjecture 6.8. By assumption, there exists a
sequence of diagonal lines between $l$ and $l^{\prime}$ such that the simplex
orderings of two consecutive lines differ by at most a single transposition.
The vineyard matching $\mathrm{vine}$ associated to these transpositions
induce morphisms of $M$ that can be seen as morphisms of $N$, thus ensuring
that it is exact on both $M$ and $N$ between $l$ and $l^{\prime}$. Hence, by
Conjecture 6.8, the vineyard matching is unique and exact.

∎

######  Remark 6.10.

Note that the assumption of Proposition 6.9 is always satisfied when $\delta$
becomes smaller than the smallest distance (in filtration function values)
between critical points of the module.

## 7\. Experiments

In this section, we showcase the performances of Algorithm 1 on various data
sets. More precisely, we evaluate the running times and approximation errors
of our approximation scheme on both synthetic and real data sets, and we
measure the empirical dependencies on the number of simplices, on the number
of lines that are used, and on the dimension $n$. We also compare our approach
to Rivet when $n=2$. All experiments were done on a laptop with AMD Ryzen 4800
CPU. Our code is publicly available at
https://gitlab.inria.fr/dloiseau/multipers, and is implemented in C++, with
Python interface.

### 7.1. Simple examples

In this section, we first provide examples of multipersistence module
approximation when the underlying multipersistence module is manually crafted
and known. In Figures 15 and 16, we provide two examples of pairs of distinct
interval decomposable multipersistence modules that have the same pointwise
Betti numbers and rank invariants. In both examples, our approximation scheme
manages to recover the correct decompositions. In Figure 17, we provide a
multipersistence module that is not interval decomposable, and our (fake)
candidate decomposition.

${k}$${k^{2}}$${k}$$\scriptstyle{\begin{pmatrix}1\\\
0\end{pmatrix}}$$\scriptstyle{\begin{pmatrix}1\\\ 0\end{pmatrix}}$

${k}$${k^{2}}$${k}$$\scriptstyle{\begin{pmatrix}1\\\
0\end{pmatrix}}$$\scriptstyle{\begin{pmatrix}0\\\ 1\end{pmatrix}}$

Figure 15:  (Top) Two distincts interval decomposable modules having the same
rank invariant. (Bottom) Output of Algorithm 1; each color corresponds to a
different summand.

${k}$${k}$${k^{2}}$${k}$${k}$${k^{2}}$${k}$${k}$$\scriptstyle{1}$$\scriptstyle{(1\quad
0)}$$\scriptstyle{(1\quad
0)}$$\scriptstyle{1}$$\scriptstyle{\begin{pmatrix}1\\\
0\end{pmatrix}}$$\scriptstyle{\begin{pmatrix}1\\\
1\end{pmatrix}}$$\scriptstyle{\begin{pmatrix}1\\\
0\end{pmatrix}}$$\scriptstyle{1}$$\scriptstyle{1}$

${k}$${k}$${k^{2}}$${k}$${k}$${k^{2}}$${k}$${k}$$\scriptstyle{1}$$\scriptstyle{(0\quad
1)}$$\scriptstyle{(1\quad
0)}$$\scriptstyle{1}$$\scriptstyle{\begin{pmatrix}1\\\
0\end{pmatrix}}$$\scriptstyle{\begin{pmatrix}1\\\
1\end{pmatrix}}$$\scriptstyle{\begin{pmatrix}0\\\
1\end{pmatrix}}$$\scriptstyle{1}$$\scriptstyle{1}$

Figure 16:  (Top) Two distincts interval decomposable modules having the same
rank invariant. (Bottom) Output of Algorithm 1; each color corresponds to a
different summand.

${k}$${k^{2}}$${k}$${0}$${k}$${k}$$\scriptstyle{\begin{pmatrix}1\\\
1\end{pmatrix}}$$\scriptstyle{(0\quad
1)}$$\scriptstyle{0}$$\scriptstyle{0}$$\scriptstyle{\begin{pmatrix}0\\\
1\end{pmatrix}}$$\scriptstyle{1}$$\scriptstyle{1}$

Figure 17:  (Top left) Filtered simplicial chain complex that leads to the
(Top right) indecomposable multipersistence module. (Bottom) Output of
Algorithm 1; each color corresponds to a different summand. The decomposition
is not real, although our output still preserves the rank invariant.

### 7.2. Convergence

In this section, we check the empirical convergence of Algorithm 1 on real
data sets. We look at a noisy circle with $1,000$ points ($60\%$ of those are
on the annulus, and the remaining $40\%$ are outliers in the square) in Figure
18, as well as three time series from the UCR archive [Che15] that were
embedded in $\mathbb{R}^{3}$ using time delay embedding in Figures 19, 20 and
21. Our approximation was computed with $n=2$ filtrations, the Alpha complex
filtration and a log-density estimation. Since we do not know the underlying
multipersistence module, we use the Euclidean norm between (a) the
multipersistence image [CB20] of our approximation and (b) a limit
multipersistence image computed on our approximation with a limit precision
(i.e., distance between consecutive lines) of $\delta=10^{-4}$ as a proxy to
measure the distance between our approximation and the true underlying module.
In all cases, one can see that the error curve decreases fast and smoothly, as
predicted by our approximation result Proposition 5.5.

Figure 18:  (Left) Noisy annulus data set colored by log density. (Middle)
Multipersistence image in dimensions 0 and 1. (Right) Error plot showing
convergence.

Figure 19:  (Left) Multipersistence image in dimensions 0 (left), 1 (up right)
and 2 (bottom right). (Right) Error plot showing convergence for the first
time series of the Coffee data set.

Figure 20:  (Left) Multipersistence image in dimensions 0 (left), 1 (up right)
and 2 (bottom right). (Right) Error plot showing convergence for the first
time series of the Worms data set.

Figure 21:  (Left) Multipersistence image in dimensions 0 (left), 1 (up right)
and 2 (bottom right). (Right) Error plot showing convergence for the first
time series of the Wine data set.

Figure 22:  (Left) Multipersistence image in dimensions 0 (left), 1 (up right)
and 2 (bottom right). (Right) Error plot showing convergence for the first
time series of the Ham data set.

### 7.3. Performance

In this section, we empirically check the dependencies between running time
and numbers of lines, simplices and dimensions.

#### 7.3.1. Synthetic data with $n=2$

We first focus on two synthetic data sets: (a) the noisy annulus of Section
7.2, and (b) a random point cloud in the unit square $[0,1]^{2}$ with one
filtration being the usual Alpha complex filtration, and the other being the
lower star filtration of a random function on the points. We show the
influence of the number of lines and simplices on the running times for both
data sets in Figure 23. As expected, the running time is linear w.r.t. the
number of lines. Furthermore, one can see that the linear coefficient depends
on the complexity of the data set, as we can see that the random filtration on
the points of the square yields longer running times than the noisy annulus.
As for the number of simplices, we empirically noted a quadratic dependency
with the running times.

Figure 23:  Running times for the noisy annulus (red) and the random point
cloud in the unit square (blue). (Top) Running times with respect to the
number of lines of the two datasets; the number of points is here fixed at
$1,000$ points. (Bottom left) Running times with respect to the number of
simplices; the number of lines is here fixed at $1,000$ lines. We also show
the curve in log-log scale in (Bottom right).

#### 7.3.2. Higher dimension

To illustrate the fact that Algorithm 1 can run with more than two
filtrations, we now focus on a synthetic data set. We uniformly sample $300$
points in the unit square $[0,1]^{2}$, and then compute its Alpha complex.
Finally, we assign to each vertex a random function value in $[0,1]^{n}$, and
compute the approximate multipersistence module (with $\delta=0.7$) induced by
the lower-star filtration of this random function. We repeated this experiment
for several numbers of dimensions $n$, and show the result in Figure 24. As
expected, there is an exponential scaling with respect to the dimension when
$\delta$ is fixed.

Figure 24: Running time w.r.t. the number of dimensions. Note that since the
precision $\delta$ is fixed, the number of lines grows exponentially with the
number of dimensions.

#### 7.3.3. Comparison with Rivet when $n=2$

As mentioned in the previous sections, Rivet [LW15] is a tool for computing
minimal presentations of 2-multipersistence modules. We provide performance
comparisons between Rivet and Algorithm 1 in the tables below.

In this first table, we focus on the noisy annulus of Section 7.2, with $80\%$
of the points uniformly sampled on the annulus, and $20\%$ are outliers in the
square; and the same filtrations than in Section 7.2.

$n$ | #simplices | Rivet |  Rivet, peak RAM | Alg. 1, $\delta=0.01$ | Alg. 1, $\delta=0.001$ | Alg. 1, peak RAM
---|---|---|---|---|---|---
$100$ | 563 | 0.02s | 9MB | 0.004s | 0.03s | 140MB
$1~{}000$ | $5~{}943$ | 0.49s | 220MB | 0.18s | 0.69s | 150MB
$5~{}000$ | $29~{}907$ | 22.13s | 5.41GB | 4.15s | 8.6s | 180MB
$7~{}000$ | $41~{}879$ | 59s | 10.29GB | 8.00s | 14.39s | 187MB
$10~{}000$ | $59~{}887$ | OOM | 12.8GB | 16s | 26s | 204MB
$20~{}000$ | $119~{}831$ | - | - | 71s | 85s | 237MB

One can see that Algorithm 1 significantly outperforms Rivet both in terms of
RAM usage and running time.

In our second table, we focus on the first time series from the data set
Coffee, and processed it as in Section 7.2. Then, we used the Vietoris-Rips
filtration as the first filtration, and density estimation as the second one.

threshold | #simplices | Rivet |  Rivet, peak RAM | Alg. 1, $\delta=0.01$ | Alg. 1, peak RAM
---|---|---|---|---|---
0.1 | $9~{}961$ | 0.28s | 38MB | 0.96s | 170MB
0.2 | $35~{}620$ | 0.79s | 80MB | 12s | 201MB
0.3 | $71~{}230$ | 1.45s | 122MB | 48s | 281MB
0.4 | $114~{}144$ | 2.6s | 166MB | 124s | 396MB
0.5 | $168~{}513$ | 5.1s | 219MB | 263s | 576MB

One can see that while Rivet works remarkably well on flag complexes such as
Vietoris-Rips, Algorithm 1 is still able to run in a reasonable amount of
time.

## 8\. Conclusion

In this article, we presented an algorithm for approximating any
$n$-multipersistence module, whose complexity, running time, and approximation
error can be controlled by user-defined parameters. We then showcased the
performances of our method on synthetic and real data sets, and provided our
code in an open-source package available at
https://gitlab.inria.fr/dloiseau/multipers.

Several questions remain open for future work. While we proved that our
candidate has bounded approximation error when approximating interval
decomposable modules, can we prove that it is optimal (in some way) among the
family of interval decomposable modules when the input is _not_ interval
decomposable? What are the properties of this output in the general case, and
can it be used in practice instead of the original (non interval decomposable)
module? Finally, can we use the stability results that are now becoming
available for specific multipersistence modules [BL21] to infer confidence
regions and convergence rates for our candidate decompositions?

## Acknowledgments

The authors would like to thank Hannah Schreiber for her tremendous help in
implementing our approximation and running our numerical experiments.

## References

  * [ABE+21] Hideto Asashiba, Mickaël Buchet, Emerson Escolar, Ken Nakashima, and Michio Yoshiwaki.  On interval decomposability of 2D persistence modules.  arXiv:1812.05261 [math], May 2021.
  * [AEK+17] Henry Adams, Tegan Emerson, Michael Kirby, Rachel Neville, Chris Peterson, Patrick Shipman, Sofya Chepushtanova, Eric Hanson, Francis Motta, and Lori Ziegelmeier.  Persistence images: a stable vector representation of persistent homology.  Journal of Machine Learning Research, 18(8):1–35, 2017.
  * [AENY19] Hideto Asashiba, Emerson Escolar, Ken Nakashima, and Michio Yoshiwaki.  On approximation of 2D persistence modules by interval-decomposables.  arXiv:1911.01637 [math], November 2019.
  * [Bje20] Håvard Bjerkevik.  Stability of higher-dimensional interval decomposable persistence modules.  arXiv:1609.02086 [cs, math], January 2020.
  * [BL18] Magnus Botnan and Michael Lesnick.  Algebraic stability of zigzag persistence modules.  Algebraic & Geometric Topology, 18(6):3133–3204, October 2018\.
  * [BL21] Andrew Blumberg and Michael Lesnick.  Stability of 2-parameter persistent homology.  arXiv:2010.09628 [cs, math], February 2021.
  * [BLO20] Magnus Botnan, Vadim Lebovici, and Steve Oudot.  Local characterizations for decomposability of 2-parameter persistence modules.  arXiv:2008.02345 [math], August 2020.
  * [BLO22] Magnus Botnan, Vadim Lebovici, and Steve Oudot.  On rectangle-decomposable 2-parameter persistence modules.  Discrete & Computational Geometry, 2022.
  * [BOO21] Magnus Botnan, Steffen Oppermann, and Steve Oudot.  Signed barcodes for multi-parameter persistence via rank decompositions and rank-exact resolutions.  arXiv:2107.06800 [cs, math], July 2021.
  * [Bub15] Peter Bubenik.  Statistical topological data analysis using persistence landscapes.  Journal of Machine Learning Research, 16(3):77–102, 2015.
  * [Car09] Gunnar Carlsson.  Topology and data.  Bulletin of The American Mathematical Society, 46:255–308, 2009\.
  * [CB20] Mathieu Carrière and Andrew Blumberg.  Multiparameter persistence image for topological machine learning.  In Advances in Neural Information Processing Systems 34 (NeurIPS 2020), pages 22432–22444. Curran Associates, Inc., 2020.
  * [CCI+20] Mathieu Carrière, Frédéric Chazal, Yuichi Ike, Théo Lacombe, Martin Royer, and Yuhei Umeda.  PersLay: a neural network layer for persistence diagrams and new graph topological signatures.  In 23rd International Conference on Artificial Intelligence and Statistics (AISTATS 2020), pages 2786–2796. PMLR, 2020.
  * [CCO17] Mathieu Carrière, Marco Cuturi, and Steve Oudot.  Sliced Wasserstein kernel for persistence diagrams.  In 34th International Conference on Machine Learning (ICML 2017), volume 70, pages 664–673. PMLR, 2017.
  * [CdGO16] Frédéric Chazal, Vin de Silva, Marc Glisse, and Steve Oudot.  The Structure and Stability of Persistence Modules.  SpringerBriefs in Mathematics. Springer International Publishing, 2016.
  * [CFK+19] René Corbet, Ulderico Fugacci, Michael Kerber, Claudia Landi, and Bei Wang.  A kernel for multi-parameter persistent homology.  Computers & Graphics: X, 2:100005, 2019.
  * [Che15] Chen, Yanping and Keogh, Eamonn and Hu, Bing and Begum, Nurjahan and Bagnall, Anthony and Mueen, Abdullah and Batista, Gustavo.  The UCR time series classification archive, 2015.
  * [CO19] Jérémy Cochoy and Steve Oudot.  Decomposition of exact pfd persistence bimodules.  Discrete & Computational Geometry, pages 1–39, 2019.
  * [CSEM06] David Cohen-Steiner, Herbert Edelsbrunner, and Dmitriy Morozov.  Vines and vineyards by updating persistence in linear time.  In 22nd Annual Symposium on Computational Geometry (SoCG 2006), pages 119–126. Association for Computing Machinery, 2006.
  * [dMV11] Vin de Silva, Dmitriy Morozov, and Mikael Vejdemo-Johansson.  Dualities in persistent (co)homology.  Inverse Problems, 27(12):124003, 2011.
  * [DW22] Tamal Dey and Yusu Wang.  Computational Topology for Data Analysis.  Cambridge University Press, 2022.
  * [DX21] Tamal Dey and Cheng Xin.  Rectangular approximation and stability of 2-parameter persistence modules.  arXiv:2108.07429 [cs, math], August 2021.
  * [DX22] Tamal Dey and Cheng Xin.  Generalized persistence algorithm for decomposing multiparameter persistence modules.  Journal of Applied and Computational Topology, 2022.
  * [EH10] Herbert Edelsbrunner and John Harer.  Computational Topology: An Introduction.  American Mathematical Society, 2010.
  * [KLO19] Michael Kerber, Michael Lesnick, and Steve Oudot.  Exact computation of the matching distance on 2-parameter persistence modules.  In 35th International Symposium on Computational Geometry (SoCG 2019), volume 129, pages 46:1–46:15. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik, 2019.
  * [KN20] Michael Kerber and Arnur Nigmetov.  Efficient approximation of the matching distance for 2-parameter persistence.  arXiv:1912.05826 [cs], March 2020.
  * [Lan18] Claudia Landi.  The rank invariant stability via interleavings.  In Research in Computational Topology, pages 1–10. Springer, 2018\.
  * [Les15] Michael Lesnick.  The theory of the interleaving distance on multidimensional persistence modules.  Foundations of Computational Mathematics, 15(3):613–650, 2015.
  * [LW15] Michael Lesnick and Matthew Wright.  Interactive visualization of 2-D persistence modules.  arXiv:1512.00180 [cs, math], December 2015.
  * [Mil20] Ezra Miller.  Homological algebra of modules over posets.  arXiv:2008.00063 [math], August 2020.
  * [MO15] Clément Maria and Steve Oudot.  Zigzag persistence via reflections and transpositions.  In 26th Annual ACM-SIAM Symposium on Discrete Algorithms (SODA 2015), pages 181–199. Society for Industrial and Applied Mathematics, 2015.
  * [Mun84] James Munkres.  Elements of algebraic topology.  CRC Press, 1984.
  * [Oud15] Steve Oudot.  Persistence Theory: From Quiver Representations to Data Analysis, volume 209 of Mathematical Surveys and Monographs.  American Mathematical Society, 2015.
  * [RHBK15] Jan Reininghaus, Stefan Huber, Ulrich Bauer, and Roland Kwitt.  A stable multi-scale kernel for topological machine learning.  In 28th IEEE Conference on Computer Vision and Pattern Recognition (CVPR 2015), pages 4741–4748. IEEE Computer Society, 2015.
  * [Vip20a] Oliver Vipond.  Local equivalence of metrics for multiparameter persistence modules.  arXiv:2004.11926 [math], April 2020.
  * [Vip20b] Oliver Vipond.  Multiparameter persistence landscapes.  Journal of Machine Learning Research, 21(61):1–38, 2020.


# Optimization of a quantum control sequence for initializing an NV spin
register

T. Chakraborty1,2, J. Zhang1 and D. Suter1,3 1Fakultät Physik, Technische
Universität Dortmund, D-44221 Dortmund, Germany  2QuTech and Kavli Institute
of Nanoscience, Delft University of Technology, 2600 GA Delft, The Netherlands
3Department of Physics and NMR Research Center, Indian Institute of Science
Education and Research, Pune 411008, India

###### Abstract

Implementation of many quantum information protocols require an efficient
initialization of the quantum register. In the present report, we optimize a
population trapping protocol for initializing a hybrid spin register
associated a single nitrogen vacancy (NV) center in diamond. We initialize the
quantum register by polarizing the electronic and the nuclear spins of the NV
with a sequence of microwave, radio-frequency and optical pulses. We use a
rate equation model to explain the distribution of population under the effect
of the optical pulses. The model is compared to the experimental data obtained
by performing partial quantum state tomography. To further increase the spin
polarisation, we propose a recursive protocol with optimized optical pulses.

## I. Introduction

A nitrogen atom replacing one carbon atom of a diamond lattice and an adjacent
carbon vacancy site form a nitrogen vacancy (NV) center in diamond, which is
an optically active atomic defect center [1, 2]. The NV center, having
excellent room temperature (RT) quantum properties like long spin coherence
time [3], offers a stable spin-photon interface [4, 5], coherent optical
transitions and so on, which have been successfully applied in a number of
aspects of the emerging field of quantum technology [6, 7, 8]. For instance,
NV centers were used to demonstrate spin-photon entanglement with scalability
[9], quantum teleportation [10], efficient implementation of quantum
algorithms [11], entanglement distribution over a multinode quantum network
[12], photonic quantum repeater [13], quantum sensing platforms with
remarkable sensitivity [14, 15] and so on.

Having a large band gap ($5.4$ eV), diamond has an unoccupied conduction band
which does not allow interaction of free electrons with the NV centers [16].
Because of its high Debye temperature ($\approx 1800$ 0C), the interaction
with phonons is weak at RT [17]. In addition, naturally abundant diamond
crystals contain magnetic impurities of very low concentration ($1.1\%$ 13C)
and growth of ultrapure single crystal diamond contain NVs with slowly
dephasing spins [18]. Thus diamond, as a host, ensures that NV centers can be
considered as isolated quantum systems, where the long lived electronic and
nuclear spins associated with NV can act as a quantum register [19]. Several
experiments demonstrated local control of the NV spins with high fidelity
through initialization, coherent manipulation and read out using optical,
microwave (MW), and radio frequency (RF) pulses [19, 20, 21, 22].

Often for quantum operations, hybrid spin register associated with NVs are
exploited [19, 23, 22, 24]. Commonly, these registers include vacancy
electronic spins being coupled to the nuclear spins from 14N or neighboring
13C atoms [2]. For efficient implementation of several NV-based quantum
information protocols, it is necessary that the initial state of the spin
register has high purity. Hence, one needs to optimize methods which
simultaneously polarize the electronic and nuclear spins of the NV center.
Although it is straightforward to achieve almost complete polarization of the
electronic spins through optical pumping, it does not allow polarization of
the nuclear spins. Due to having an ultra long spin coherence time, nuclear
spins are important candidates for storing and processing quantum information.
Hence, to use nuclear spin as a resource in quantum protocols, it is necessary
to perform the nuclear polarization deterministically. In this context a
number of techniques are reported. Polarization of 15N nuclear spin is
performed by accessing the excited state level anti-crossing for a single NV
at RT [25]. By means of transferring polarization from electronic spins,
hyperpolarization of 13C nuclear spin ensembles is achieved using a sequence
of laser and microwave pulses, and by exploiting the fact that the nuclear
spin eigenstates do not have same quantization axes in different electron spin
manifolds, which occurs because NV electron spin S=1.[26]. Moreover, bulk
nuclear polarization for NV centers are achieved by applying laser induced
dynamic nuclear polarization (DNP) schemes [27, 28].

In the present case, we initialize the spin register through polarizing the
electronic and nuclear spins by applying a sequence of laser, MW and RF
pulses. We consider a hybrid two qutrit system where the vacancy electron is
coupled to the 14N nuclear spin, and perform the spin manipulation in a
subspace of the spins. Such a protocol for initialization of the NV spins has
been demonstrated earlier [23, 22, 24]. However, In the present case we report
a detailed optimization of this sequence, specifically the optical pumping
procedure. Employing optically detected magnetic resonance (ODMR) spectroscopy
techniques we perform partial state tomography and investigate the population
dynamics of the relevant quantum states using a rate equation model. We have
performed an additional simulation to test if the present sequence, when
executed in an iterative way with optimized laser pulses, can enhance the
polarization further. Thus, our optimized sequence with known optical pumping
parameter can perform NV spin manipulation in a deterministic way.

The paper is organized as the following. First, we describe the system and the
initialization sequence. Then we present the experimental results and
demonstrate its comparison with the solution of the formulated rate equation
model. Finally, we present the results of the multi-cycle simulation.

## II. System and polarization procedure

Figure 1: Sequence of microwave (MW), radio frequency (RF) and $532$ nm laser
pulses used for initialization. The sequence contains three parts: segment $1$
and $2$ for addressing the nuclear spins $m_{I}=-1$ and $+1$, and segment $3$
measures free induction decay of the electronic spin. The energy levels with
the relevant transitions are shown where the population of the states are
represented by filled circles.

The experiments were performed in a home-built confocal microscope setup
integrated with an electronic circuit which generates the necessary MW and RF
signals for controlling the electronic and nuclear spins, as described in ref.
[22]. We chose the vacancy electronic spin $S=1$ coupled with the 14N nuclear
spin $I=1$ of an isolated NV center embedded in a $99.998\%$ 12C enriched
diamond as the quantum register for our experiments. The system Hamiltonian,
when the field is aligned along the NV axis, is

$\displaystyle H$ $\displaystyle=$ $\displaystyle
DS_{z}^{2}-{\gamma}_{e}BS_{z}+PI_{z}^{2}-{\gamma}_{n}BI_{z}+AS_{z}I_{z}.$ (1)

$S_{z}$ and $I_{z}$ are the $z$-components of the electronic and nuclear
spins, the zero field splitting is $D=2.87$ GHz, electronic and nuclear
gyromagnetic ratios are ${\gamma}_{e}=-28$ GHzT-1 and ${\gamma}_{n}=-3.1$
MHzT-1, respectively, the nuclear quadrapole coupling $P=4.5$ MHz and the
hyperfine coupling $A=-2.16$ MHz [22]. We performed the experiments in a
magnetic field $B=6.1$ mT along the symmetry axis of the NV.

Our initialization procedure relies on a population trapping (PT) protocol
where by applying a number of $\pi$ pulses we sequentially accumulate the
population into the $m_{I}=0$ state. For this purpose, we perform the spin
manipulation in a subspace spanned by $m_{s}=0,-1$ and $m_{I}=0,\pm 1$. We
write the states in the notation $|m_{s},m_{I}\rangle$ in this paper. A
coherent population trapping (CPT) protocol couples an initial and a target
state to an intermediate excited state by applying coherent electromagnetic
fields, which trap the population in the superposition of the initial and the
target states [29, 30]. However, in our case, instead of applying coherent
driving fields simultaneously, we sequentially apply MW, RF and optical pulses
to drive the population between the initial state $|0,\pm 1\rangle$ and the
target state $|0,0\rangle$ through certain intermediate states. Therefore, we
term this method a “population trapping” protocol. Fig.1 schematically
represents the PT pulse sequence. We illustrate the PT protocol later in this
section. The relevant energy eigenstates of the six level system are shown in
Fig.1. The present sequence manipulates the spins associated with the NV-
charge state. The transition frequencies and Rabi frequencies for the
electronic and nuclear transitions, which we address in the sequence, are
shown in table 1. In the beginning, a $5$$\mu$s long laser pulse initializes
the system into the bright ($m_{s}=0$) state of the electron spin while the
nuclear spin is fully depolarized. Now we aim to optimally polarize both the
electron and the nuclear spins by applying a combination of MW, RF and laser
pulses. We perform this task in two segments: through segment(seg)1 and seg2,
we address the population associated with the $m_{I}=-1$ and $+1$ states in
turn. Seg3 contains the pulses which we use for analysing the resulting state.
Seg1 includes a MW and an RF $\pi$ pulse which drive the transitions
$|0,-1\rangle$ $\leftrightarrow$ $|-1,-1\rangle$ and $|-1,-1\rangle$
$\leftrightarrow$ $|-1,0\rangle$. Although these pulses polarize the nuclear
spin, the electron spin becomes depolarized. To repolarize the electronic
spin, we apply a second laser pulse. The effect of this laser pulse is
primarily a polarisation of the electron spin, but it also partly depolarizes
the 14N nuclear spin [22]. We therefore adjust the duration of the laser pulse
to maximize the electronic and the nuclear spin polarization. This task is
performed by measuring the influence of the laser pulse duration on the
population transfer dynamics and by interpreting the results using a rate
equation model, which is discussed in the subsequent section.

Table 1:  Relevant transition frequencies and Rabi frequencies, that we use in our initialization sequence Transitions | Transition frequency | Rabi frequency
---|---|---
$|0,-1\rangle$ $\leftrightarrow$ $|-1,-1\rangle$ |  $2.696$ GHz |  $8.3$ MHz
$|0,+1\rangle$ $\leftrightarrow$ $|-1,+1\rangle$ |  $2.694$ GHz |  $8.3$ MHz
$|-1,-1\rangle$ $\leftrightarrow$ $|-1,0\rangle$ |  $2.801$ MHz |  $3.87$ kHz
$|-1,+1\rangle$ $\leftrightarrow$ $|-1,0\rangle$ |  $7.095$ MHz |  $3.55$ kHz

Seg2 transfers the population from $|0,+1\rangle$ to $|-1,0\rangle$ using a MW
and an RF $\pi$ pulse which act selectively on the transitions $|0,+1\rangle$
$\leftrightarrow$ $|-1,+1\rangle$ and $|-1,+1\rangle$ $\leftrightarrow$
$|-1,0\rangle$, and swap the populations between them. Next, in close analogy
to seg1, we repolarize the electronic spin with a laser pulse and optimize its
duration. A proper optimization of the control pulses maximises the population
of the $|00\rangle$ state. The optimization of the laser pulses and the
experimental analysis of the final populations of the relevant states are
described in the next section.

## III. Experimental results and comparison with the rate equation model

Figure 2: Electron spin resonance spectra obtained by Fourier transforming the
free induction decay data for the electronic spin. The four spectra
corresponds to four different durations of the laser pulse L1 in seg1. The
amplitudes of the three lines in the spectra are proportional to the
populations of the three nuclear spin states as marked in the figure.  Figure
3: The variation of spectral amplitudes that correspond to the populations of
the nuclear spin states $m_{I}=0,\pm 1$ and their sum as a function of the
duration of the laser pulse L1 in seg1. The solid lines represent the
simulation using the rate equation model. We used $1/k_{S}=0.27$ $\mu$s and
$1/k_{I}=4.76$ $\mu$s [22].

We implement the proposed initialization protocol and perform the optimization
of the L1 and L2 laser pulses in two steps. After polarizing the electronic
spin using the $5$$\mu$s laser pulse, we test and optimize seg1, and analyze
the resulting populations of the states through the measurement protocol of
seg3. During this process, we do not implement seg2. To investigate the
redistribution of population under the influence of the laser pulse L1, we run
seg1 for different values of the laser pulse duration, and perform partial
state tomography. This is done by capturing Ramsey type free induction decays
(FIDs) of the electronic spin using the control sequence in seg3. In seg3, the
FID sequence contains two ${\pi}/2$ MW pulses which are separated by a free
evolution time $\tau$ and are resonant between the $m_{s}=0$ and $-1$ levels.
The first ${\pi}/2$ pulse creates a coherent superposition of the
$|m_{s}=0\rangle$ and $|m_{s}=-1\rangle$ states, which evolves for a time
$\tau$. The coherence is then converted into population difference by the
second ${\pi}/2$ pulse and read out using a laser pulse of $300$ ns duration.
We apply hard (high power) ${\pi}/2$ pulses such that all the three nuclear
spin-conserving transitions between $m_{s}=0$ and $-1$ are excited. A Fourier
transform of the time domain FID data allows us to capture spectral amplitudes
associated with the $|m_{I}=0,\pm 1\rangle$ states in the frequency domain.
Thus the spectra contain three lines whose amplitudes $A_{m_{I}}$ are
proportional to the difference of populations between the $m_{S}=0$ and $-1$
states: $A_{m_{I}}=P_{|0,m_{I}\rangle}-P_{|-1,m_{I}\rangle}$.

To explain the laser induced population transfer dynamics among the $6$ energy
levels of the spin register, we formulate a rate equation model. We define a
column matrix representing the populations of the six states and call it the
population vector $\vec{P}$. We express the rate equation in the following
form:

$\displaystyle\frac{d\vec{P}}{dt}$ $\displaystyle=$ $\displaystyle
M(k_{s},k_{I})\vec{P}.$ (2)

Since the spin manipulation is performed in the subspace spanned by the
$|m_{S}=0,-1\rangle\otimes|m_{I}=0,\pm 1\rangle$ states, $m_{S}=+1$ state
remains unoccupied throughout the spin control sequence. The populations of
the states in $\vec{P}$ follow the order
$(m_{s},m_{I})=(0,-1;0,+1;0,0;-1,-1;-1,+1;-1,0)$ which is the same order for
the energy levels shown in Fig.1. The rate matrix $M$ describes the transfer
of population between the states due to the effect of the laser pulse and is
given in the appendix. We assume that the electronic spin changes from
$m_{S}=-1$ to $0$ with a rate of $k_{s}$ and the nuclear spin changes between
any of the states $|m_{I}=0,\pm 1\rangle$ at a rate $k_{I}$. Assuming that the
sum of the populations is unity, the spectral amplitudes for the 14N nuclear
spins $m_{I}=0$ and $\pm 1$ are normalized with reference to the initial state
amplitudes experimentally measured after applying the $5$$\mu$s long laser
pulse, which equally populates all nuclear spin sublevels. Thus, the
normalized population vector for this state is
$\vec{P}=\frac{1}{3}(1,1,1,0,0,0)$. The MW and RF pulses of seg1 ideally
convert this into $\vec{P}=\frac{1}{3}(0,1,1,0,0,1)$ , which is the initial
condition for the dynamics during the laser pulse. The solution of Eq.2 for
this initial condition is given in the appendix.

Fig.2 exhibits $4$ experimental spectra measured after the implementation of
seg1 with durations of $50$, $500$, $1700$ and $4000$ ns of the laser pulse
L1. The spectra consist of three lines associated with the nuclear spins
states $m_{I}=0$ and $\pm 1$ which are separated by the NV hyperfine
interaction. The MW and RF pulses in seg1 do not control the population of the
$|0,+1\rangle$ state. Thus it maintains its initial population of $1/3$ after
seg1. For the shortest duration of the L1 pulse ($50$ ns), the electron spin
does not repolarize significantly. Accordingly, the $m_{I}=0$ states in the
$m_{s}=0$ and $-1$ subspaces have almost equal population, resulting in a very
small amplitude $A_{0}$ of the $m_{I}=0$ line in the spectrum. However, with
increasing the duration of the L1 pulse, the population of the
$|m_{I}=0,m_{s}=0\rangle$ state increases, as the nuclear spin is mostly
unaffected while the electronic spin is transferred from $m_{S}=-1$ to
$m_{S}=0$. If the laser pulse duration is increased further, the nuclear spin
becomes depolarized. Fig.2 shows that when the duration of L1 is $4$$\mu$s,
the three lines have equal amplitudes.

The evolution of the amplitudes $A_{0}$ and $A_{\pm 1}$ and the total
population of the three nuclear spin states in the $m_{S}=0$ subspace are
shown in Fig.3. The amplitudes increase with increasing pulse duration.
$A_{0}$ reaches a maximum when the duration of the L1 pulse is $500$ ns. We
optimize the laser pulse duration such that we obtain the maximum population
of the $|0,0\rangle$ state. This optimized pulse allows us to obtain the
maximum polarization of the electron spin with the minimum loss of the nuclear
spin polarization. When the duration of the L1 pulse reaches $4$$\mu$s, the
total population of the $m_{S}=0$ state reaches unity which signifies complete
electronic polarization. $A_{0}$ and $A_{-1}$ converge to $1/3$ which
indicates complete nuclear spin depolarization. $A_{+1}$ maintains a constant
value of $1/3$ for all durations of the laser pulse. The amplitudes simulated
using our rate equation model are in excellent agreement with the experimental
data. Under our experimental conditions, $k_{s}$ and $k_{I}$ are similar to
the values reported in ref. [22].

Next, we set the duration of the L1 pulse to $500$ ns and proceed with the
experiments using the sequence of pulses contained in seg2. In the beginning
of seg2, we start with the populations $\vec{P}=(0.07,0.33,0.55,0,0,0.05)$.
Since we address the $m_{I}=+1$ state in this stage, the MW and RF pulses
create the population vector $\vec{P}=(0.07,0,0.55,0,0.05,0.33)$. Using this
as the initial condition we solve the rate equation model. The solution gives
analytical expression of the population of the six levels as a function of the
laser pulse duration which is shown in the appendix.

Four spectra measured for $60$, $460$, $1410$ and $4000$ ns duration of the
laser pulse L2 are shown in Fig.4. Similar to the earlier case, due to
electron spin polarization, an increase in pulse duration increases the
population difference between the states $|0,0\rangle$ and $|-1,0\rangle$.
$A_{0}$ reaches its maximum for a pulse duration of $460$ ns. The evolution of
the amplitudes with the laser pulse duration can better be observed in Fig. 5.
The nuclear spin starts to depolarize, populating the $m_{I}$ $=\pm 1$ states
as the laser pulse duration becomes longer and eventually the nuclear spin
depolarizes completely . The rate equation model is well consistent with
experimental data as it can be seen in Fig. 5.

Figure 4: Fourier transform of free induction decay of the electronic spin for
four different durations of L2. The spectra show three peaks for the three
nuclear spin states $m_{I}=0$ and $\pm 1$ of 14N.  Figure 5: Amplitudes of the
three resonance lines corresponding to the populations of $m_{I}=0,\pm 1$ as a
function of the laser pulse (L2) duration. The simulation performed using the
rate equation model (with $1/k_{S}=0.27$ $\mu$s and $1/k_{I}=4.76$ $\mu$s
[22]) are shown by the solid curves.

## IV. Iteration to enhance purity

Earlier experimental studies of population transfer dynamics [22, 24]have
shown that laser induced depolarization of the nuclear spin limits the overall
purity. An important contribution to loosing nuclear spin polarization is the
charge state conversion dynamics of the NV center. The two charge states have
different hyperfine interactions. Thus, switching from one charge state to
another can result in nuclear depolarization [22]. This prediction is
supported by recent results showing that the NV- to NV0 conversion rate is
reduced if light $594$ nm light is used instead of 532 nm [24]. In this
section we show that it is possible to enhance the polarization of the hybrid
spin register beyond the values obtained for a single cycle in section III,
when MW, RF and properly optimized $532$ nm laser pulses are applied in an
iterative way. We rely on the idea that when we apply the initialization
sequence $[$seg1$+$seg2$]$ as shown in Fig.1 for multiple cycles, each
subsequent cycle uses the final populations of a given cycle as the initial
state. Thus, when a cycle starts with an improved purity of the target
$|0,0\rangle$ state, the purity of the state is enhanced further at the end of
the cycle. We experimentally have implemented the sequence $[$seg1 $+$
seg2$]$$\times N$ for $N=1$ to $5$ and have analyzed the resulting state using
seg3. However, in this way the resulting gain in the purity of the
$|0,0\rangle$ state is very small. An additional improvement can be obtained
by adjusting the durations of the laser pulses independently, using the rate
equation discussed in section III for each pulse depending on the relevant
initial states. Below we explain a protocol and demonstrate the results of the
simulation we have performed.

Table 2: Simulated values of the laser pulse duration and population of the $|0,0\rangle$ state after each of the three cycles for seg1 and seg2 Number of cycle | Duration of L1(ns) |  $P_{|0,0\rangle}$ after seg1 (%) |  Duration of L2(ns) |  $P_{|0,0\rangle}$ after seg2 (%)
---|---|---|---|---
$1$ |  $500$ |  $55$ |  $460$ |  $70.9$
$2$ |  $156$ |  $72.3$ |  $140$ |  $73.5$
$3$ |  $20$ |  $74$ |  $0$ |  $74$

We simulate the implementation of seg1 and seg2 in an iterative way: we
simulate the population redistribution for a sequence $[$seg1 $+$
seg2$]$$\times N$ where the duration of the pulses L1 and L2 are optimized for
each cycle. For the $1^{st}$ cycle we choose the duration of the L1 and L2
pulses as $500$ and $460$ ns as the resonance amplitude for the nuclear spin
$m_{I}=0$ is maximum for this value [Fig.3]. For the $2^{nd}$ cycle, we
consider the experimentally obtained final population after $1^{st}$ cycle as
the initial state and solve the rate equation model for L1and L2 pulses. In
this cycle, we observe that maximum transfer of population to $|0,0\rangle$
state occurs when the duration of the L1and L2 pulses are $156$ and $140$ ns.
In this similar way, we solve the rate equation for the $3^{rd}$ cycle. The
optimized duration of the laser pulses and the obtained purity
$P_{|0,0\rangle}$ for three cycles are shown in Fig.6(a) and (b). Table 2
summarizes these values for both the segments. Moreover, we also have
simulated the population redistribution for a sequence $[$(seg1)$\times N$ $+$
(seg2)$\times N]$ where we first optimize the duration of the laser pulses for
three consecutive cycles of seg1 and then we do the similar task for seg2.
However, this method could not enhance the purity than what we obtained using
the sequence $[$seg1 $+$ seg2$]$$\times N$.

Figure 6: (a) Optimized duration of the laser pulses for three cycles for seg1
and (b) the resulting population of the $|0,0\rangle$ state after each cycle.
(c) Optimized laser pulse duration for seg2 and (d) purity after each cycle.

## V. Discussions and conclusion

The present paper demonstrates a detailed optimization procedure for
initializing a hybrid quantum spin register by performing the spin
manipulation in a reduced subspace of the spins. We apply a series of MW and
RF pulses to swap the electronic and nuclear spin polarization, and laser
pulses for optically pumping the electronic spin. The duration of the laser
pulses are optimized in such a way that we maximize the electronic and nuclear
spin polarization. For this purpose, we formulate a rate equation model and
investigate the polarization dynamics as a function of the duration of the
laser pulse. By performing partial state tomography we experimentally
determine the population of the relevant eigenstates.

We follow an approach developed in ref. [22] to optimize the laser pulse
duration. However, in this case we have applied a different polarization
method which is based on a population trapping protocol. With this protocol,
we access the subspace spanned by $|m_{S}=0,-1\rangle\otimes|m_{I}=0,\pm
1\rangle$ states and sequentially populate the target state $|0,0\rangle$. For
this purpose, a number of MW, RF and optical pulses are employed to drive the
population from the initial state $|0,\pm 1\rangle$ to the target state
through certain intermediate states in $m_{S}=-1$ subspace. We have optimized
the optical pumping processes to maximize the spin polarization. To describe
the present experimental results, we have used a rate equation model which
explains laser-induced redistribution of the population among the six spin-
levels. Furthermore, to reduce the nuclear depolarization, we have formulated
a recursive method where in each cycle the laser pulses perform the optical
pumping in an optimum way. Our simulation shows that when multiple cycles of
the control sequence, with each cycle having an optimum laser pulse duration,
are implemented, an enhancement of the polarization can be achieved.

Initializing a quantum register to an eigenstate is important for
implementation of quantum information tasks. In these protocols, optical
control of the spins plays a significant role. The present results of the
optimization of the optical control pulses can help developing spin
manipulation protocols in deterministic way. One advantage of our method is
that the spin control can be performed in an arbitrary magnetic field. We have
developed the present method for an NV spin register. However, this spin
manipulation approach is general. Hence, for other atomic defect centers in
diamond which has optically addressable spin transitions, a similar method can
be formulated.

###### Acknowledgements.

This work was supported by the DFG through grants 192/19-2, 192/34-1 and by
the MERCUR foundation through grant Pr-2013-0003.

## appendix

During the optical pulse, the redistribution of the populations $\vec{P}$ can
be described as

$\displaystyle\frac{d\vec{P}}{dt}$ $\displaystyle=$ $\displaystyle
M(k_{s},k_{I})\vec{P},$

where the matrix $M(k_{S},k_{I})$ of rate constants is

$M(k_{S},k_{I})=\left[\begin{array}[]{cccccc}-2k_{I}&k_{I}&k_{I}&k_{S}&0&0\\\
k_{I}&-2k_{I}&k_{I}&0&k_{S}&0\\\ k_{I}&k_{I}&-2k_{I}&0&0&k_{S}\\\
0&0&0&-k_{S}&0&0\\\ 0&0&0&0&-k_{S}&0\\\ 0&0&0&0&0&-k_{S}\end{array}\right]$

As discussed in the main text, the basis states are
$(m_{s},m_{I})=(0,-1;0,+1;0,0;-1,-1;-1,+1;-1,0)$. The constants $k_{I}$ and
$k_{S}$ describe the rates at which the electronic and nuclear spin
polarisation are changed. Under our experimental conditions, $1/k_{S}=0.27$
$\mu$s and $1/k_{I}=4.76$ $\mu$s [22].

In Seg1 the solution obtained for the initial condition
$\vec{P}(0)=\frac{1}{3}(1,1,1,0,0,0)$ is

$\displaystyle\vec{P}_{seg1}(t)$
$\displaystyle=\frac{1}{3}[1-\frac{k_{I}(e^{-k_{S}t}-e^{-3k_{I}t})}{(3k_{I}-k_{S})},$
$\displaystyle
1-\frac{(2k_{I}-k_{S})e^{-3k_{I}t}+k_{I}e^{-k_{S}t}}{(3k_{I}-k_{S})},$
$\displaystyle
1-\frac{(k_{I}-k_{S})e^{-k_{S}t}+(k_{s}-k_{I})e^{-3k_{I}t}}{(3k_{I}-k_{S})},$
$\displaystyle 0,0,e^{-k_{S}t}].$

For initial condition $\vec{P}=(0.07,0,0.55,0,0.05,0.33)$ in seg2, the
solution is

$\displaystyle\vec{P}_{seg2}(t)=$
$\displaystyle[0.34+\frac{e^{-3k_{I}t}(0.26k_{S}-0.4k_{I})-0.38k_{I}e^{-k_{S}t}}{3k_{I}-k_{S}},$
$\displaystyle
0.34-\frac{e^{-k_{S}t}(0.38k_{I}-0.05k_{S})-e^{-3k_{I}t}(0.63k_{I}-0.29k_{S})}{3k_{I}-k_{S}},$
$\displaystyle
0.34+\frac{e^{-3k_{I}t}(1.03k_{I}-0.55k_{S})-e^{-k_{S}t}(0.38k_{I}-0.33k_{S})}{3k_{I}-k_{S}}]$
$\displaystyle 0$ $\displaystyle,0.05e^{-k_{S}t},0.33e^{-k_{S}t}].$

## References

  * Jelezko and Wrachtrup [2006] F. Jelezko and J. Wrachtrup, Phys. Stat. Sol. (a) 13, 3207 (2006).
  * Suter and Jelezko [2017] D. Suter and F. Jelezko, Progress in nuclear magnetic resonance spectroscopy 98, 50 (2017).
  * Balasubramanian  _et al._ [2009] G. Balasubramanian, P. Neumann, D. Twitchen, M. Markham, R. Kolesov, N. Mizuochi, J. Isoya, J. Achard, J. Beck, J. Tissler, _et al._ , Nature materials 8, 383 (2009).
  * Kubanek  _et al._ [2012] A. Kubanek, A. Sipahigil, E. Togan, M. Goldman, Y. Chu, N. de Leon, A. Zibrov, and M. Lukin, in  _APS March Meeting Abstracts_ , Vol. 2012 (2012) pp. D29–003.
  * Bernien  _et al._ [2012] H. Bernien, L. Childress, L. Robledo, M. Markham, D. Twitchen, and R. Hanson, Physical Review Letters 108, 043604 (2012).
  * Pezzagna and Meijer [2021] S. Pezzagna and J. Meijer, Applied Physics Reviews 8, 011308 (2021).
  * Childress and Hanson [2013] L. Childress and R. Hanson, MRS bulletin 38, 134 (2013).
  * Fabbri and Hernández-Gómez [2020] N. Fabbri and S. Hernández-Gómez, Frontiers in Physics 8, 652 (2020).
  * Vasconcelos  _et al._ [2020] R. Vasconcelos, S. Reisenbauer, C. Salter, G. Wachter, D. Wirtitsch, J. Schmiedmayer, P. Walther, and M. Trupke, npj Quantum Information 6, 1 (2020).
  * Pfaff  _et al._ [2014] W. Pfaff, B. J. Hensen, H. Bernien, S. B. van Dam, M. S. Blok, T. H. Taminiau, M. J. Tiggelman, R. N. Schouten, M. Markham, D. J. Twitchen, _et al._ , Science 345, 532 (2014).
  * Zhang  _et al._ [2020] J. Zhang, S. S. Hegde, and D. Suter, Physical Review Letters 125, 030501 (2020).
  * Pompili  _et al._ [2021] M. Pompili, S. L. Hermans, S. Baier, H. K. Beukers, P. C. Humphreys, R. N. Schouten, R. F. Vermeulen, M. J. Tiggelman, L. d. S. Martins, B. Dirkse, _et al._ , arXiv preprint arXiv:2102.04471 (2021).
  * Nemoto  _et al._ [2016] K. Nemoto, M. Trupke, S. J. Devitt, B. Scharfenberger, K. Buczak, J. Schmiedmayer, and W. J. Munro, Scientific reports 6, 1 (2016).
  * Degen  _et al._ [2017] C. L. Degen, F. Reinhard, and P. Cappellaro, Reviews of modern physics 89, 035002 (2017).
  * Schirhagl  _et al._ [2014] R. Schirhagl, K. Chang, M. Loretz, and C. L. Degen, Annual review of physical chemistry 65, 83 (2014).
  * Larsson and Delaney [2008] J. Larsson and P. Delaney, Physical Review B 77, 165201 (2008).
  * Doherty  _et al._ [2012] M. Doherty, F. Dolde, H. Fedder, F. Jelezko, J. Wrachtrup, N. Manson, and L. Hollenberg, Physical Review B 85, 205203 (2012).
  * Chakraborty  _et al._ [2019] T. Chakraborty, F. Lehmann, J. Zhang, S. Borgsdorf, N. Wöhrl, R. Remfort, V. Buck, U. Köhler, and D. Suter, Physical Review Materials 3, 065205 (2019).
  * Dutt  _et al._ [2007] M. V. G. Dutt, L. Childress, L. Jiang, E. Togan, J. Maze, F. Jelezko, A. S. Zibrov, P. R. Hemmer, and M. D. Lukin, Phys. Rev. Lett. 316, 1312 (2007).
  * Childress  _et al._ [2006] L. Childress, M. G. Dutt, J. Taylor, A. Zibrov, F. Jelezko, J. Wrachtrup, P. Hemmer, and M. Lukin, Science 314, 281 (2006).
  * Gaebel  _et al._ [2006] T. Gaebel, M. Domhan, I. Popa, C. Wittmann, P. Neumann, F. Jelezko, J. R. Rabeau, N. Stavrias, A. D. Greentree, S. Prawer, _et al._ , Nature Physics 2, 408 (2006).
  * Chakraborty  _et al._ [2017] T. Chakraborty, J. Zhang, and D. Suter, New Journal of Physics 19, 073030 (2017).
  * Pagliero  _et al._ [2014] D. Pagliero, A. Laraoui, J. D. Henshaw, and C. A. Meriles, Applied Physics Letters 105, 242402 (2014).
  * Rao  _et al._ [2020] K. R. K. Rao, Y. Wang, J. Zhang, and D. Suter, Physical Review A 101, 013835 (2020).
  * Jacques  _et al._ [2009] V. Jacques, P. Neumann, J. Beck, M. Markham, D. Twitchen, J. Meijer, F. Kaiser, G. Balasubramanian, F. Jelezko, and J. Wrachtrup, Physical review letters 102, 057403 (2009).
  * Álvarez  _et al._ [2015] G. A. Álvarez, C. O. Bretschneider, R. Fischer, P. London, H. Kanda, S. Onoda, J. Isoya, D. Gershoni, and L. Frydman, Nature communications 6, 1 (2015).
  * Scheuer  _et al._ [2016] J. Scheuer, I. Schwartz, Q. Chen, D. Schulze-Sünninghausen, P. Carl, P. Höfer, A. Retzker, H. Sumiya, J. Isoya, B. Luy, _et al._ , New journal of Physics 18, 013040 (2016).
  * King  _et al._ [2015] J. P. King, K. Jeong, C. C. Vassiliou, C. S. Shin, R. H. Page, C. E. Avalos, H.-J. Wang, and A. Pines, Nature communications 6, 1 (2015).
  * Jamonneau  _et al._ [2016] P. Jamonneau, G. Hétet, A. Dréau, J.-F. Roch, and V. Jacques, Physical review letters 116, 043603 (2016).
  * Nicolas  _et al._ [2018] L. Nicolas, T. Delord, P. Jamonneau, R. Coto, J. Maze, V. Jacques, and G. Hétet, New Journal of Physics 20, 033007 (2018).


# Randomized adaptive quantum state preparation

Alicia B. Magann  Quantum Algorithms and Applications Collaboratory, Sandia
National Laboratories, Albuquerque, New Mexico 87185, USA  Sophia E. Economou
Department of Physics, Virginia Tech, Blacksburg, Virginia 24061, USA
Christian Arenz  School of Electrical, Computer, and Energy Engineering,
Arizona State University, Tempe, Arizona 85287, USA

(August 28, 2024)

###### Abstract

We develop an adaptive method for quantum state preparation that utilizes
randomness as an essential component and that does not require classical
optimization. Instead, a cost function is minimized to prepare a desired
quantum state through an adaptively constructed quantum circuit, where each
adaptive step is informed by feedback from gradient measurements in which the
associated tangent space directions are randomized. We provide theoretical
arguments and numerical evidence that convergence to the target state can be
achieved for almost all initial states. We investigate different randomization
procedures and develop lower bounds on the expected cost function change,
which allows for drawing connections to barren plateaus and for assessing the
applicability of the algorithm to large-scale problems.

_Introduction_.— Methods for preparing quantum states are an integral
component of any quantum technology. For example, the preparation of states
that encode ground, excited, and thermal states of many-body systems is a key
element in quantum simulation [1, 2]. Ground state preparation can also be
leveraged to solve combinatorial optimization problems, with a variety of
applications including in routing and scheduling [3]. While a variety of
quantum state preparation algorithms exist, the growing availability of noisy,
intermediate-scale quantum [4] devices has inspired immense interest in
variational methods for preparing desired quantum states [5, 6, 7, 8]. These
variational quantum algorithms (VQAs) are heuristics that function by
classically optimizing over a set of parameters that enter into a quantum
circuit whose structure is typically fixed. The parameterized quantum circuit
is executed on a quantum device and serves as an ansatz to minimize a cost
function, $J$, whose global minimum is achieved for the desired target state.

Even in the absence of noise, a variety of challenges are present in VQAs. On
the quantum device, for example, one must select a quantum circuit ansatz and
associated initial state out of a formidably large design space. Meanwhile,
the difficulty of the cost function minimization means that the challenges on
the classical side can be even more significant [9], often involving the
navigation of optimization landscapes that contain barren plateaus and
suboptimal local minima.

Figure 1: Schematic representation of randomized adaptive quantum state
preparation. Each adaptive step $k$ involves first estimating the gradient
$\frac{dJ(0)}{d\theta_{k}}$ at $\theta_{k}=0$ of a cost function $J$ with
respect to a parameter $\theta_{k}$, e.g., a rotation angle, whose
corresponding direction is randomized through conjugation with a random
unitary transformation $V_{k}$ (brown). The parameter $\theta_{k}$ is then
updated in the negative direction of the gradient (blue arrow), which moves
the system closer to the target state, and the quantum circuit $U_{k}$ is
extended accordingly.

To overcome these challenges, approaches have been proposed that utilize
feedback from qubit measurements to _adaptively_ construct a quantum circuit
to minimize $J$. The first such algorithm was the Adaptive Derivative-
Assembled Problem-Tailored Variational Quantum Eigensolver (ADAPT-VQE) [10,
11], which was also adapted to combinatorial optimization problems [12].
Instead of relying on a predefined ansatz, ADAPT-VQE grows it in a layer-wise
manner tailored to the problem, in tandem with classical optimization over the
quantum circuit parameters.

Other methods have considered defining the structure of the ansatz _a priori_
and then performing layer-wise optimization to adaptively set the circuit
parameter values [13, 14, 15, 16]. Adaptive methods that do not require any
classical optimization have also been developed, including the feedback-based
algorithm for quantum optimization [17, 18] and methods based on Riemannian
gradient flows [19, 20]. In methods where each adaptive step is efficiently
implementable, e.g., [10, 11, 12, 13, 14, 15, 16, 17, 18] and the approximate
method in [19], the circuit growth can get stuck when the gradient vanishes.
This means that these methods can face similar issues as conventional VQAs
that are prone to converge to suboptimal solutions.

Here, we propose _randomized_ adaptive quantum state preparation as a generic,
adaptive quantum algorithm that minimizes up-front design choices, does not
require classical optimization, and allows for preparing target quantum states
from arbitrary (random) initial states. As a consequence of the latter point,
the algorithm, or its key subroutine depicted in Fig. 1, can be readily
combined with other state preparation methods to improve convergence and state
preparation fidelities. We provide theoretical arguments and numerical
evidence that substantiate our claims about convergence, explore different
methods for achieving randomization in practice, and develop lower bounds on
the expected change in $J$ at each adaptive step. These bounds give guarantees
for how much the cost function value changes when randomization is used,
thereby allowing us to relate the efficiency of this randomized approach to
the existence of barren plateaus [21]. We go on to discuss how this approach
can be applied to cooling in open quantum systems and mixed state preparation
in general.

_Randomized adaptive quantum algorithms._ — We consider minimizing cost
functions of the form

$\displaystyle
J_{k}=\mathinner{\langle{\psi_{k}}|}H_{p}\mathinner{|{\psi_{k}}\rangle},$ (1)

by creating states
$\mathinner{|{\psi_{k}}\rangle}=U_{k}\mathinner{|{\psi_{0}}\rangle}$ via an
adaptively constructed quantum circuit $U_{k}$, $k=0,1,\cdots$. The goal is to
apply this circuit to a fixed initial state $\mathinner{|{\psi_{0}}\rangle}$
to achieve $J_{k+1}\leq J_{k}$ in each adaptive step $k$. Here, $H_{p}$ is a
Hermitian operator whose ground state $\mathinner{|{E_{\text{min}}}\rangle}$,
with corresponding eigenvalue $E_{\text{min}}$, is taken to be the target
state. We note that in general, knowledge of the initial and target states
$\mathinner{|{\psi_{0}}\rangle}$ and $\mathinner{|{E_{\text{min}}}\rangle}$ is
not required. However, for the preparation of an arbitrary, known target state
$\mathinner{|{\psi_{T}}\rangle}$, $H_{p}$ can be defined as
$H_{p}=\mathds{1}-\mathinner{|{\psi_{T}}\rangle}\mathinner{\langle{\psi_{T}}|}$.

We consider a quantum circuit that is adaptively created according to

$\displaystyle U_{k+1}=e^{-i\theta_{k}H_{k}}U_{k},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
U_{0}=\mathds{1},$ (2)

where in each step $k$ we move into the negative direction of the gradient of
$J$ with respect to $\theta_{k}$ by setting
$\theta_{k}\leftarrow-\gamma\frac{dJ(0)}{d\theta_{k}}$, where
$\frac{dJ(0)}{d\theta_{k}}=i\mathinner{\langle{\psi_{k}}|}[H_{k},H_{p}]\mathinner{|{\psi_{k}}\rangle}$
is the gradient of
$J(\theta_{k})=\mathinner{\langle{\psi_{k}}|}e^{i\theta_{k}H_{k}}H_{p}e^{-i\theta_{k}H_{k}}\mathinner{|{\psi_{k}}\rangle}$
evaluated at $\theta_{k}=0$. For sufficiently small learning rates $\gamma$,
this ensures that $J_{k+1}\leq J_{k}$. References [10, 12, 17, 18, 19]
consider adaptive procedures similar to Eq. (2), and in cases where the
circuit growth gets stuck in suboptimal solutions, i.e., when the gradient
vanishes, the utility of incorporating randomness into the circuit structure
to escape these suboptima has been observed numerically in [19].

In this work, we utilize randomness to overcome challenges associated with
convergence through the introduction of an _intrinsically randomized_
framework for quantum state preparation in which the $H_{k}$’s are selected at
random. In the following, we provide theoretical arguments that this
randomization enables convergence from almost all initial states to arbitrary
target states.

We first note that the cost function gradient can be expressed as

$\displaystyle\frac{dJ(0)}{d\theta_{k}}=\langle\text{grad}J[U_{k}],iH_{k}\rangle,$
(3)

where $\langle\cdot,\cdot\rangle$ denotes the Hilbert-Schmidt inner product
and
$\text{grad}J[U_{k}]=[\mathinner{|{\psi_{k}}\rangle}\mathinner{\langle{\psi_{k}}|},H_{p}]$
is the (Riemannian) gradient (up to multiplication with $U_{k}$ from the
right) of $J[U_{k}]=\langle\psi_{0}|U_{k}^{\dagger}H_{p}U_{k}|\psi_{0}\rangle$
with respect to the unitary transformation $U_{k}$ [20, 22, 23]. Both $iH_{k}$
and $\text{grad}J[U_{k}]$ belong to the special unitary algebra
$\mathfrak{su}(2^{n})$ consisting of all traceless and anti-Hermitian
$2^{n}\times 2^{n}$ matrices where $n$ is the number of qubits. From this
geometric perspective we can now deduce two different cases (i) and (ii) for
when $\frac{dJ(0)}{d\theta_{k}}$ vanishes 111We neglect the trivial third case
that $iH_{k}=0$.. In case (i), the gradient vanishes when $iH_{k}$ is
orthogonal to $\text{grad}J[U_{k}]\neq 0$, while in case (ii)
$\text{grad}J[U_{k}]=0$, which happens when $H_{p}$ commutes with
$\mathinner{|{\psi_{k}}\rangle}\mathinner{\langle{\psi_{k}}|}$.

We first discuss case (i). Typically, no assumptions can be made on whether
$\text{grad}J[U_{k}]$ moves into a lower dimensional subspace of
$\mathfrak{su}(2^{n})$ when growing the circuit. As such, the situation that
$iH_{k}$ has no overlap with $\text{grad}J[U_{k}]$ can occur, e.g., when
$iH_{k}$ is an element of a subspace over which $\text{grad}J[U_{k}]$ has no
support. To overcome this issue, we propose to create each $H_{k}$ at random.
This can be achieved by conjugating a traceless Hermitian operator $H$ by a
Haar random unitary transformation $V_{k}$ in each adaptive step, i.e., such
that $H_{k}=V_{k}^{\dagger}HV_{k}$. Since $H_{k}$ is created uniformly
randomly according to the Haar measure, the probability that $iH_{k}$ is
orthogonal to $\text{grad}J[U_{k}]$ is zero. That is, for almost all $iH_{k}$,
but a set of measure zero, case (i) does not occur. While Haar random
unitaries are not efficiently implementable, below we discuss the efficient
implementation via approximate unitary 2-designs [25].

We now focus on case (ii). For cost functions of the form (1), the set of
critical points where $\text{grad}J[U_{k}]$ vanishes consists of global optima
and saddle points only [26, 27]. Under mild assumptions on the nature of the
saddle points (strict saddles), relevant works from the classical machine
learning and optimization literature have found that saddle points are avoided
for almost all initial conditions [28, 29, 30]. We thus expect that randomized
adaptive quantum state preparation will almost surely converge to the ground
state of $H_{p}$ for almost all initial states
$\mathinner{|{\psi_{0}}\rangle}$. We remark that the convergence result cannot
hold for _all_ initial states, as we immediately see that for eigenstates of
$H_{p}$, $\text{grad}J[U_{0}]=0$. However, the situation that
$\text{grad}J[U_{0}]=0$ can be avoided with probability one when the initial
state is randomized too.

Each step $k$ of randomized adaptive quantum state preparation can now be
summarized as follows. First, the unitary transformation
$V^{\dagger}_{k}e^{-i\theta_{k}H}V_{k}$, whose generator $H$ is randomized
through conjugation with a random unitary $V_{k}$, is applied to the state
$\mathinner{|{\psi_{k}}\rangle}$. Second, the gradient,
$\frac{dJ(0)}{d\theta_{k}}$, is estimated, e.g., using the parameter shift
rule [31, 32, 33, 34, 35]. Third, the parameter $\theta_{k}$ is updated in the
negative direction of the gradient, i.e.,
$\theta_{k}\leftarrow-\gamma\frac{dJ(0)}{d\theta_{k}}$.

_Efficiency and relation to barren plateaus_.—Randomized adaptive quantum
state preparation is not expected to be efficient in general, as ground state
preparation is QMA complete [36, 37]. Here, we show that for a given problem,
the efficiency can be related to the scaling of the gradient (3) with the
system size, and therefore, to the existence of barren plateaus [21], i.e.,
exponentially flat regions in the optimization landscape where the variance of
the gradient vanishes exponentially in the number of qubits $n$.

In the Supplemental Material, we show that if we select
$\gamma=1/(4\|H_{p}\|_{2})$ and assume $\|H_{k}\|_{2}=1$, where
$\|\cdot\|_{2}$ denotes the spectral norm, then the cost function change
$\Delta J_{k}=J_{k}-J_{k+1}$ is lower bounded by

$\displaystyle\Delta
J_{k}\geq\frac{1}{8\|H_{p}\|_{2}}\left(\frac{dJ(0)}{d\theta_{k}}\right)^{2}.$
(4)

If we assume that it takes $M$ steps to create the ground state up to an error
$\epsilon$, i.e., $J_{M}=E_{\text{min}}+\epsilon$, we find from (4) that $M$
is upper bounded by,

$\displaystyle
M\leq\frac{C_{\epsilon}}{\min\limits_{k\in[0,M]}\langle\text{grad}J[U_{k}],iH_{k}\rangle^{2}},$
(5)

where the constant in the numerator is given by
$C_{\epsilon}=8\|H_{p}\|_{2}\big{(}J_{0}-(E_{\text{min}}+\epsilon)\big{)}$.
Thus, if $\langle\text{grad}J_{k}[U_{k}],iH_{k}\rangle^{2}$ does not vanish
faster than $1/\text{poly}(n)$, the ground state of $H_{p}$ can be prepared up
to precision $\epsilon$ (in the corresponding eigenvalue) in polynomially many
steps. We remark that the minimum in the denominator of Eq. (5) explicitly
depends on $M$, and thereby on the random path taken. It is interesting to
note that similar expressions are obtained in adiabatic state preparation
[38], where the scaling of the adiabatic state preparation time $T$ is
determined by the smallest value of the spectral gap $\Delta(t)$, taken over
all times, i.e., $\min_{t\in[0,T]}\Delta(t)$ [39].

We proceed by investigating the efficiency of different randomization
strategies. At each step $k$, the expected cost function change
$\mathbb{E}_{H_{k}}\Delta J_{k}$ is lower bounded by the variance
$\mathbb{E}_{H_{k}}\langle\text{grad}J[U_{k}],iH_{k}\rangle^{2}$ of the
gradient, up to the prefactor in (4), assuming that the expectation vanishes.
This suggests that sampling $V_{k}$ from unitary 2-designs suffices to obtain
convergence to the ground state. This observation is further substantiated by
the numerical simulations in Fig. 2(a). However, the appearance of the
variance in the bound (4) means that we additionally expect the efficiency and
practical utility of this method to be closely related to the existence of
barren plateaus. While barren plateaus pose a major challenge to the
scalability of VQAs [21, 40, 41], examples have been found where the variance
of the gradient does not vanish faster than $1/\text{poly}(n)$ [42].
Leveraging these instances for efficient realizations of randomized adaptive
quantum state preparation will be the subject of future studies.

We now consider lower bounds for $\mathbb{E}_{H_{k}}\Delta J_{k}$ to derive
guarantees for how much $J_{k}$ can be improved when randomization is used.
When $H_{k}$ is created through conjugation by a unitary $V_{k}$ sampled from
a unitary 2-design, we show in the Supplemental Material that
$\mathbb{E}_{V_{k}}\Delta J_{k}$ is lower bounded by

$\displaystyle\mathbb{E}_{V_{k}}\Delta
J_{k}\geq\frac{\text{Tr}\\{H^{2}\\}}{4\|H_{p}\|_{2}}\frac{\text{Var}_{\psi_{k}}(H_{p})}{2^{2n}-1},$
(6)

where
$\text{Var}_{\psi_{k}}(H_{p})=\mathinner{\langle{\psi_{k}}|}H_{p}^{2}\mathinner{|{\psi_{k}}\rangle}-\mathinner{\langle{\psi_{k}}|}H_{p}\mathinner{|{\psi_{k}}\rangle}^{2}$
is the variance of $H_{p}$ with respect to the state
$\mathinner{|{\psi_{k}}\rangle}$. Since $\text{Var}_{\psi_{k}}(H_{p})$ is
bounded from above by a constant that is independent of the system dimension,
we see that the bound in (6) vanishes exponentially in $n$.

Another way of creating random $H_{k}$’s is by sampling uniformly from an
operator pool $\mathcal{A}$ whose size we denote by $|\mathcal{A}|$. While in
this case, situation (i) can occur when an operator is selected that is
orthogonal to $\text{grad}J_{k}[U_{k}]$, on average, we have that

$\displaystyle\mathbb{E}_{H_{k}\in\mathcal{A}}\Delta
J_{k}\geq\frac{1}{8\|H_{p}\|_{2}|\mathcal{A}|}\sum_{H_{k}\in\mathcal{A}}\langle\text{grad}J[U_{k}],iH_{k}\rangle^{2},$
(7)

which suggests that situation (i) can be avoided on average for sufficiently
large operator pools, e.g., when
$\text{span}\\{\mathcal{A}\\}=\mathfrak{su}(2^{n})$. Note that in this case,
convergence to the ground state can also be obtained by simply continuing to
sample from $\mathcal{A}$ when (i) occurs. However, in this situation we
expect an exponential runtime, as $M$ scales as $1/|\mathcal{A}|$. For
$|\mathcal{A}|=\text{poly}(n)$, on the other hand, convergence to the ground
state is no longer guaranteed, as the adaptive procedure can get stuck in
suboptima where the gradient vanishes, due to (i).

Figure 2: Performance of randomized adaptive quantum state preparation when
$H_{p}$ is an Ising Hamiltonian with all-to-all couplings between $n$ spins.
Each data point corresponds to the average taken over $100$ different
algorithm realizations and initial states. In (a), the approximation ratio,
$\alpha$, is plotted as a function of the number of adaptive steps, $M$, shown
on a logarithmic scale. In each step, the randomization of $H_{k}$ is
implemented by conjugating the Pauli operator $X\otimes\mathds{1}$ with a Haar
random unitary transformation (right triangles) and with a unitary
transformation sampled from an approximate unitary 2-design (left triangles),
as depicted in Fig. 1, created using the sequence in [43] with $\ell=1$, and
also by sampling $H_{k}$ uniformly from an operator pool $\mathcal{A}$
containing all $2^{2n}-1$ Pauli operators (circles). In (b), $\alpha$ is
plotted as a function of the number of pool elements, $|\mathcal{A}|$, shown
on a logarithmic scale, for fixed $M$. Insets show semilog plots of the
scaling with respect to $n$ for both cases to achieve $\alpha>0.99$ in (a) and
$\alpha>0.9$ in (b).

In Fig. 2, we numerically investigate this tradeoff and study the convergence
of randomized adaptive quantum state preparation for the problem of preparing
the ground state of an Ising Hamiltonian with couplings present between all
$n$ spins, which is equivalent to solving the combinatorial optimization
problem MaxCut on a complete graph [44]. We consider the approximation ratio
$\alpha=J_{k}/E_{\text{min}}$ as our figure of merit. In Fig. 2(a), we plot
$\alpha$ as a function of the number of adaptive steps $M$ for the different
randomization strategies described above. In Fig. 2(b), we consider fixed
$M=35000$ and investigate $\alpha$ as a function of the pool size
$|\mathcal{A}|$. The pool size is increased by adding successively heigher
weight Pauli operators to $\mathcal{A}$, until the full pool used in Fig. 2(a)
containing all $2^{2n}-1$ terms is formed. We observe that the curves in Fig.
2(a) are almost identical, suggesting that the three different randomization
strategies converge in the same manner to the ground state for which
$\alpha=1$. Furthermore, the curves in Fig. 2(b) suggest that a full Pauli
operator pool is not needed to obtain convergence to the ground state.
However, the inset semilog plots in 2(a) and 2(b) do suggest an exponential
scaling of the number of adaptive steps, $M$, and the number of operators in
the pool, $|\mathcal{A}|$, with respect to $n$.

_Mixed states_.— Thus far, we have discussed the preparation of pure (ground)
states from initial pure states. Here, we explore generalizations to preparing
a target pure state from an initially mixed state, and vice versa. Since it is
not possible to create pure quantum states from mixed quantum states in a
closed quantum system through unitary transformations, we consider an
extended, “dilated” space by coupling a set of system qubits, $S$, to a set of
auxiliary qubits, $A$. We assume that the combined system is initially in a
separable state $\rho_{0}^{SA}=\rho_{0}^{S}\otimes\rho_{0}^{A}$, where we
denote the initial states of $S$ and $A$ by $\rho_{0}^{S}$ and $\rho_{0}^{A}$,
respectively. We then consider growing a quantum circuit over the full
composite system according to (2), in order to adaptively create the state
$\rho_{k}^{SA}=U_{k}\rho_{0}^{SA}U_{k}^{\dagger}$. The state $\rho_{k}^{S}$ at
the $k$-th adaptive step is then obtained by tracing over the degrees of
freedom of the auxiliary qubits in subsystem $A$, i.e.,
$\rho_{k}^{S}=\text{Tr}_{A}\\{\rho_{k}^{SA}\\}$.

In this setting, we now consider the task of preparing the system qubits in a
target pure state $\mathinner{|{\psi_{T}}\rangle}$. The cost function (1)
becomes
$J_{k}=1-\mathinner{\langle{\psi_{T}}|}\rho_{k}^{S}\mathinner{|{\psi_{T}}\rangle}$,
and an adaptive change is described by
$J(\theta_{k})=1-\mathinner{\langle{\psi_{T}}|}\text{Tr}_{A}\\{e^{i\theta_{k}H_{k}}\rho_{k}^{SA}e^{-i\theta_{k}H_{k}}\\}\mathinner{|{\psi_{T}}\rangle}$.
This yields the gradient

$\displaystyle\frac{dJ(0)}{d\theta_{k}}=-\langle[\rho_{k}^{SA},\mathinner{|{\psi_{T}}\rangle}\mathinner{\langle{\psi_{T}}|}\otimes\mathds{1}_{A}],iH_{k}\rangle,$
(8)

where we have made use of the fact that the cost function above can be
obtained from (1) by setting
$H_{p}=\mathds{1}-\mathinner{|{\psi_{T}}\rangle}\mathinner{\langle{\psi_{T}}|}\otimes\mathds{1}_{A}$
where $\mathds{1}_{A}$ denotes the identity operator on subsystem $A$.
Decreasing $J_{k}$ in each step can now be achieved by moving into the
negative direction of the gradient given by Eq. (8).

If the auxiliary qubits are initially in a pure state, then due to the
Stinespring dilation [45], there exists a unitary transformation over the
composite system that allows for creating every state for the system qubits in
subsystem $S$ as long as subsystem $A$ has dimension at most $d_{S}^{2}$,
where $d_{S}$ is the dimension of subsystem $S$. For uniformly randomized
$H_{k}$, the gradient can only vanish (with probability 1) when
$\text{grad}J[U_{k}]=[\rho_{k}^{SA},\mathinner{|{\psi_{T}}\rangle}\mathinner{\langle{\psi_{T}}|}\otimes\mathds{1}_{A}]=0$,
i.e., at critical points that are given by saddle points and global optima, as
in the closed system case [27]. Thus, when the initial system state
$\rho_{0}^{S}$ does not commute with the target system state
$\mathinner{|{\psi_{T}}\rangle}\mathinner{\langle{\psi_{T}}|}$, we expect to
obtain convergence almost surely to the target state. Although the fully mixed
state $\rho_{0}^{S}=\mathds{1}/d_{A}$ trivially commutes with every target
state, if we additionally apply a random unitary transformation to
$\rho_{0}^{SA}$, we can almost surely obtain convergence to a generic pure
state. Thus, randomized adaptive quantum state preparation allows for
“cooling” the system from a state of infinite temperature, i.e., the fully
mixed state, to a pure state of zero temperature by adaptively “dumping”
entropy into the auxiliary system.

_Conclusions_.— We have introduced an algorithm for preparing quantum states
that has favorable convergence properties and is applicable to almost all
initial states. Knowledge of the initial and target quantum states is not
required. The algorithm leverages randomization as the primary innovation, and
operates by minimizing a cost function through an adaptively constructed
quantum circuit. Each adaptive update step is informed by gradient
measurements in which the associated tangent space directions are randomized.
We have presented lower bounds on the average improvement that is obtained in
each step, and have numerically studied the behavior of the algorithm for
different randomization methods. We additionally discussed a generalization to
mixed states that could be leveraged for thermal state preparation, an
application area which is currently receiving significant interest [46, 47,
48, 49, 50, 51, 52].

A tradeoff is that, on one hand, the consideration of a uniformly random
tangent space direction in each step allows for convergence to the target
state, while on the other hand, the gradients may become exponentially small
in the system size, thereby causing the convergence time to diverge and making
the algorithm impractical for large-scale problems. This should not be
surprising, as selecting directions at random is far from being optimal. The
largest gradient, and thereby the largest (guaranteed) cost function change,
is obtained when $iH_{k}=\text{grad}J[U_{k}]$. This situation corresponds to
the full gradient flow, which is, in general, not efficiently implementable.
It is an interesting question how gradient flows can be efficiently
approximated [19] while maintaining convergence. To that end, the algorithm we
present could be modified to project into random subspaces, rather than into a
single random direction.

Furthermore, we emphasize that the algorithm possesses significant design
flexibility that can be harnessed to increase its practicality, for example,
by developing application-specific adaptations that tailor the randomness to
the problem instance, e.g., by taking into account the symmetries of the
target state [53, 54, 55, 56], and by studying how much the rate of
convergence can be improved by incorporating classical optimization in
different ways [10]. More fundamentally, we view the randomized adaptive
circuit update, which can be applied to arbitrary input states and is depicted
in Fig. 1, as a _subroutine_ that could be incorporated into or appended onto
other algorithms [38, 57], e.g., those utilizing problem-specific ansätze that
aren’t random, to improve convergence and state preparation fidelities.

We conclude by noting that we have focused this work on randomized adaptive
quantum state preparation in the circuit model of quantum computing. However,
randomness and adaptive constructions can also be created outside the circuit
model, e.g. through random fields [58, 59]. Extensions to this latter setting
could open up new approaches for state preparation in both closed and open
quantum systems that can be driven by an applied field, including qubit
systems, analog quantum simulators, molecules, and materials.

###### Acknowledgements.

_Acknowledgements_.— We gratefully acknowledge feedback from Mohan Sarovar and
Andrew Baczewski on this work. We also would like to thank Thomas Schulte-
Herbrüggen, Gunther Dirr, and Emanuel Malvetti for fruitful discussions.
A.B.M. acknowledges support from Sandia National Laboratories’ Laboratory
Directed Research and Development Program under the Truman Fellowship. Sandia
National Laboratories is a multimission laboratory managed and operated by
National Technology & Engineering Solutions of Sandia, LLC, a wholly owned
subsidiary of Honeywell International Inc., for the U.S. Department of
Energy’s National Nuclear Security Administration under contract DE-NA0003525.
This paper describes objective technical results and analysis. Any subjective
views or opinions that might be expressed in the paper do not necessarily
represent the views of the U.S. Department of Energy or the United States
Government. C. A. and S. E. acknowledge support from the National Science
Foundation (Grant No. 2231328).

## References

  * McArdle  _et al._ [2020] S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, Quantum computational chemistry, Rev. Mod. Phys. 92, 015003 (2020).
  * Cao  _et al._ [2019] Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. Sawaya, _et al._ , Quantum chemistry in the age of quantum computing, Chem. Rev. 119, 10856 (2019).
  * Papadimitriou and Steiglitz [1998] C. Papadimitriou and K. Steiglitz,  _Combinatorial Optimization: Algorithms and Complexity_ (Dover Publications, 1998).
  * Preskill [2018] J. Preskill, Quantum Computing in the NISQ era and beyond, Quantum 2, 79 (2018).
  * Peruzzo  _et al._ [2014] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’brien, A variational eigenvalue solver on a photonic quantum processor, Nat. Commun. 5, 1 (2014).
  * Cerezo  _et al._ [2021a] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, _et al._ , Variational quantum algorithms, Nat. Rev. Phys. 3, 625 (2021a).
  * Farhi  _et al._ [2014] E. Farhi, J. Goldstone, and S. Gutmann, A Quantum Approximate Optimization Algorithm (2014), arXiv:1411.4028 [quant-ph] .
  * Magann  _et al._ [2021] A. B. Magann, C. Arenz, M. D. Grace, T.-S. Ho, R. L. Kosut, J. R. McClean, H. A. Rabitz, and M. Sarovar, From pulses to circuits and back again: A quantum optimal control perspective on variational quantum algorithms, PRX Quantum 2, 010101 (2021).
  * Bittel and Kliesch [2021] L. Bittel and M. Kliesch, Training variational quantum algorithms is np-hard, Phys. Rev. Lett. 127, 120502 (2021).
  * Grimsley  _et al._ [2019] H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, An adaptive variational algorithm for exact molecular simulations on a quantum computer, Nat. Commun. 10, 1 (2019).
  * Tang  _et al._ [2021] H. L. Tang, V. Shkolnikov, G. S. Barron, H. R. Grimsley, N. J. Mayhall, E. Barnes, and S. E. Economou, Qubit-adapt-vqe: An adaptive algorithm for constructing hardware-efficient ansätze on a quantum processor, PRX Quantum 2, 020310 (2021).
  * Zhu  _et al._ [2022] L. Zhu, H. L. Tang, G. S. Barron, F. A. Calderon-Vargas, N. J. Mayhall, E. Barnes, and S. E. Economou, Adaptive quantum approximate optimization algorithm for solving combinatorial problems on a quantum computer, Phys. Rev. Research 4, 033029 (2022).
  * Carolan  _et al._ [2020] J. Carolan, M. Mohseni, J. P. Olson, M. Prabhu, C. Chen, D. Bunandar, M. Y. Niu, N. C. Harris, F. N. Wong, M. Hochberg, _et al._ , Variational quantum unsampling on a quantum photonic processor, Nat. Phys. 16, 322 (2020).
  * Skolik  _et al._ [2021] A. Skolik, J. R. McClean, M. Mohseni, P. van der Smagt, and M. Leib, Layerwise learning for quantum neural networks, Quantum Mach. Intell. 3, 1 (2021).
  * Campos  _et al._ [2021a] E. Campos, A. Nasrallah, and J. Biamonte, Abrupt transitions in variational quantum circuit training, Phys. Rev. A 103, 032607 (2021a).
  * Campos  _et al._ [2021b] E. Campos, D. Rabinovich, V. Akshay, and J. Biamonte, Training saturation in layerwise quantum approximate optimization, Phys. Rev. A 104, L030401 (2021b).
  * Magann  _et al._ [2022a] A. B. Magann, K. M. Rudinger, M. D. Grace, and M. Sarovar, Feedback-based quantum optimization, Phys. Rev. Lett. 129, 250502 (2022a).
  * Magann  _et al._ [2022b] A. B. Magann, K. M. Rudinger, M. D. Grace, and M. Sarovar, Lyapunov-control-inspired strategies for quantum combinatorial optimization, Phys. Rev. A 106, 062414 (2022b).
  * Wiersema and Killoran [2022] R. Wiersema and N. Killoran, Optimizing quantum circuits with riemannian gradient flow (2022), arXiv:2202.06976 [quant-ph] .
  * Schulte-Herbrüggen  _et al._ [2010] T. Schulte-Herbrüggen, S. J. Glaser, G. Dirr, and U. Helmke, Gradient flows for optimization in quantum information and quantum dynamics: Foundations and applications, Rev. Math. Phys. 22, 597 (2010).
  * McClean  _et al._ [2018] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, Barren plateaus in quantum neural network training landscapes, Nat. Commun. 9, 1 (2018).
  * Helmke and Moore [2012] U. Helmke and J. B. Moore,  _Optimization and dynamical systems_ (Springer Science & Business Media, 2012).
  * Absil  _et al._ [2009] P.-A. Absil, R. Mahony, and R. Sepulchre, Optimization algorithms on matrix manifolds, in  _Optimization Algorithms on Matrix Manifolds_ (Princeton University Press, 2009).
  * Note [1] We neglect the trivial third case that $iH_{k}=0$.
  * Dankert  _et al._ [2009] C. Dankert, R. Cleve, J. Emerson, and E. Livine, Exact and approximate unitary 2-designs and their application to fidelity estimation, Phys. Rev. A 80, 012304 (2009).
  * Brockett [1989] R. Brockett, Least squares matching problems, Linear Algebra and its Applications 122-124, 761 (1989), special Issue on Linear Systems and Control.
  * Wu  _et al._ [2007] R. Wu, H. Rabitz, and M. Hsieh, Characterization of the critical submanifolds in quantum ensemble control landscapes, J. Phys. A 41, 015006 (2007).
  * Ge  _et al._ [2015] R. Ge, F. Huang, C. Jin, and Y. Yuan, Escaping from saddle points — online stochastic gradient for tensor decomposition, in  _Proceedings of The 28th Conference on Learning Theory_, Proceedings of Machine Learning Research, Vol. 40, edited by P. Grünwald, E. Hazan, and S. Kale (PMLR, Paris, France, 2015) pp. 797–842.
  * Jin  _et al._ [2017] C. Jin, R. Ge, P. Netrapalli, S. M. Kakade, and M. I. Jordan, How to escape saddle points efficiently, in  _Proceedings of the 34th International Conference on Machine Learning_, Proceedings of Machine Learning Research, Vol. 70, edited by D. Precup and Y. W. Teh (PMLR, International Convention Centre, Sydney, Australia, 2017) pp. 1724–1732.
  * Lee  _et al._ [2019] J. D. Lee, I. Panageas, G. Piliouras, M. Simchowitz, M. I. Jordan, and B. Recht, First-order methods almost always avoid strict saddle points, Math. Program 176, 311 (2019).
  * Li  _et al._ [2017] J. Li, X. Yang, X. Peng, and C.-P. Sun, Hybrid quantum-classical approach to quantum optimal control, Phys. Rev. Lett. 118, 150503 (2017).
  * Mitarai  _et al._ [2018] K. Mitarai, M. Negoro, M. Kitagawa, and K. Fujii, Quantum circuit learning, Phys. Rev. A 98, 032309 (2018).
  * Schuld  _et al._ [2019] M. Schuld, V. Bergholm, C. Gogolin, J. Izaac, and N. Killoran, Evaluating analytic gradients on quantum hardware, Phys. Rev. A 99, 032331 (2019).
  * Wierichs  _et al._ [2022] D. Wierichs, J. Izaac, C. Wang, and C. Y.-Y. Lin, General parameter-shift rules for quantum gradients, Quantum 6, 677 (2022).
  * Kyriienko and Elfving [2021] O. Kyriienko and V. E. Elfving, Generalized quantum circuit differentiation rules, Phys. Rev. A 104, 052417 (2021).
  * Kitaev  _et al._ [2002] A. Y. Kitaev, A. Shen, M. N. Vyalyi, and M. N. Vyalyi,  _Classical and quantum computation_, 47 (American Mathematical Soc., 2002).
  * Kempe  _et al._ [2006] J. Kempe, A. Kitaev, and O. Regev, The complexity of the local hamiltonian problem, SIAM J. Comput. 35, 1070 (2006).
  * Farhi  _et al._ [2000] E. Farhi, J. Goldstone, S. Gutmann, and M. Sipser, Quantum computation by adiabatic evolution, arXiv:quant-ph/0001106 (2000).
  * Albash and Lidar [2018] T. Albash and D. A. Lidar, Adiabatic quantum computation, Rev. Mod. Phys. 90, 015002 (2018).
  * Ortiz Marrero  _et al._ [2021] C. Ortiz Marrero, M. Kieferová, and N. Wiebe, Entanglement-induced barren plateaus, PRX Quantum 2, 040316 (2021).
  * Wang  _et al._ [2021a] S. Wang, E. Fontana, M. Cerezo, K. Sharma, A. Sone, L. Cincio, and P. J. Coles, Noise-induced barren plateaus in variational quantum algorithms, Nature communications 12, 1 (2021a).
  * Cerezo  _et al._ [2021b] M. Cerezo, A. Sone, T. Volkoff, L. Cincio, and P. J. Coles, Cost function dependent barren plateaus in shallow parametrized quantum circuits, Nat. Commun. 12, 1791 (2021b).
  * Nakata  _et al._ [2017] Y. Nakata, C. Hirche, C. Morgan, and A. Winter, Unitary 2-designs from random x- and z-diagonal unitaries, J. Math. Phys 58, 052203 (2017).
  * Barahona  _et al._ [1988] F. Barahona, M. Grötschel, M. Jünger, and G. Reinelt, An application of combinatorial optimization to statistical physics and circuit layout design, Oper. Res 36, 493 (1988).
  * Nielsen and Chuang [2002] M. A. Nielsen and I. Chuang, Quantum computation and quantum information (2002).
  * Wu and Hsieh [2019] J. Wu and T. H. Hsieh, Variational thermal quantum simulation via thermofield double states, Phys. Rev. Lett. 123, 220502 (2019).
  * Moussa [2019] J. E. Moussa, Low-depth quantum metropolis algorithm (2019), arXiv:1903.01451 [quant-ph] .
  * Motta  _et al._ [2020] M. Motta, C. Sun, A. T. K. Tan, M. J. O’Rourke, E. Ye, A. J. Minnich, F. G. S. L. Brandão, and G. K.-L. Chan, Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution, Nat. Phys. 16, 205 (2020).
  * Wang  _et al._ [2021b] Y. Wang, G. Li, and X. Wang, Variational quantum gibbs state preparation with a truncated taylor series, Phys. Rev. Applied 16, 054035 (2021b).
  * Warren  _et al._ [2022] A. Warren, L. Zhu, N. J. Mayhall, E. Barnes, and S. E. Economou, Adaptive variational algorithms for quantum gibbs state preparation (2022), arXiv:2203.12757 [quant-ph] .
  * Chowdhury  _et al._ [2020] A. N. Chowdhury, G. H. Low, and N. Wiebe, A variational quantum algorithm for preparing quantum gibbs states (2020), arXiv:2002.00055 [quant-ph] .
  * Sewell  _et al._ [2022] T. J. Sewell, C. D. White, and B. Swingle, Thermal multi-scale entanglement renormalization ansatz for variational gibbs state preparation (2022), arXiv:2210.16419 [quant-ph] .
  * Gard  _et al._ [2020] B. T. Gard, L. Zhu, G. S. Barron, N. J. Mayhall, S. E. Economou, and E. Barnes, Efficient symmetry-preserving state preparation circuits for the variational quantum eigensolver algorithm, npj Quantum Inf. 6, 1 (2020).
  * Shkolnikov  _et al._ [2021] V. O. Shkolnikov, N. J. Mayhall, S. E. Economou, and E. Barnes, Avoiding symmetry roadblocks and minimizing the measurement overhead of adaptive variational quantum eigensolvers (2021), arXiv:2109.05340 [quant-ph] .
  * Yen  _et al._ [2019] T.-C. Yen, R. A. Lang, and A. F. Izmaylov, Exact and approximate symmetry projectors for the electronic structure problem on a quantum computer, J. Chem. Phys 151, 164111 (2019), https://doi.org/10.1063/1.5110682 .
  * Anselmetti  _et al._ [2021] G.-L. R. Anselmetti, D. Wierichs, C. Gogolin, and R. M. Parrish, Local, expressive, quantum-number-preserving vqe ansätze for fermionic systems, New J. Phys. 23, 113010 (2021).
  * Poulin and Wocjan [2009] D. Poulin and P. Wocjan, Preparing ground states of quantum many-body systems on a quantum computer, Phys. Rev. Lett. 102, 130503 (2009).
  * Banchi  _et al._ [2017] L. Banchi, D. Burgarth, and M. J. Kastoryano, Driven quantum dynamics: Will it blend?, Phys. Rev. X 7, 041015 (2017).
  * Yang  _et al._ [2020] P. Yang, M. Yu, R. Betzholz, C. Arenz, and J. Cai, Complete quantum-state tomography with a local random field, Phys. Rev. Lett. 124, 010405 (2020).
  * Beck [2017] A. Beck,  _First-order methods in optimization_ (SIAM, 2017).

##  Appendix A Lower bounding the cost function change

We consider

$\displaystyle
J(x)=\mathinner{\langle{\psi_{k}}|}e^{ixH_{k}}H_{p}e^{-ixH_{k}}\mathinner{|{\psi_{k}}\rangle},$
(9)

and note that
$J(x=0)=\mathinner{\langle{\psi_{k}}|}H_{p}\mathinner{|{\psi_{k}}\rangle}=J_{k}$
and $J\left(x=-\gamma\frac{dJ(0)}{d\theta}\right)=J_{k+1}$.

##  A.0.1 Recalling results from classical optimization

In order to lower bound the cost function change $\Delta J_{k}=J_{k}-J_{k+1}$,
we first recall some standard tools from classical optimization, particularly
the gradient descent algorithm. For a continuously differentiable cost
function $J(x)$, $x\in\mathbb{R}$, with Lipschitz continuous gradient
$\frac{\partial J(x)}{\partial x}$, a direct application of Taylor’s theorem
[60] yields,

$\displaystyle J(y)\leq J(x)+(y-x)\frac{\partial J(x)}{\partial
x}+\frac{L_{f}^{2}}{2}(y-x)^{2},$ (10)

where $L_{f}$ is the Lipschitz constant satisfying,

$\displaystyle\frac{\left|\frac{\partial J(x)}{\partial x}-\frac{\partial
J(y)}{\partial y}\right|}{|x-y|}\leq L_{f}.$ (11)

Setting $x=0$ and $y=-\gamma\frac{\partial J(x)}{\partial x}|_{x=0}$ in (10)
and picking $\gamma=\frac{1}{L_{f}}$ gives

$\displaystyle J_{k}-J_{k+1}\geq\frac{1}{2L_{f}}\left(\frac{\partial
J(0)}{\partial x}\right)^{2},$ (12)

which is a well known result for gradient algorithms [60]. Thus, in order to
lower bound the cost function change $\Delta J_{k}$, it remains to determine
the Lipschitz constant $L_{f}$ for our setting.

##  A.0.2 Determination of the Lipschitz constant

We first note that $\frac{\partial J(x)}{\partial
x}=\mathinner{\langle{\psi_{k}}|}e^{ixH_{k}}[iH_{k},H_{p}]e^{-ixH_{k}}\mathinner{|{\psi_{k}}\rangle}\equiv\mathinner{\langle{\psi(x)}|}A\mathinner{|{\psi(x)}\rangle}$,
where we defined
$\mathinner{|{\psi(x)}\rangle}=e^{-ixH_{k}}\mathinner{|{\psi_{k}}\rangle}$ and
$A=[iH_{k},H_{p}]$. We can than upper bound $\left|\frac{\partial
J(x)}{\partial x}-\frac{\partial J(y)}{\partial y}\right|$ by

$\displaystyle\left|\frac{\partial J(x)}{\partial x}-\frac{\partial
J(y)}{\partial y}\right|$
$\displaystyle=|\mathinner{\langle{\psi(x)}|}A\mathinner{|{\psi(x)}\rangle}-\mathinner{\langle{\psi(y)}|}A\mathinner{|{\psi(y)}\rangle}|$
(13)
$\displaystyle=\frac{1}{2}\left|\big{(}\mathinner{\langle{\psi(x)}|}+\mathinner{\langle{\psi(y)}|}\big{)}A\big{(}\mathinner{|{\psi(x)}\rangle}-\mathinner{|{\psi(y)}\rangle}\big{)}+\big{(}\mathinner{\langle{\psi(x)}|}-\mathinner{\langle{\psi(y)}|}\big{)}A\big{(}\mathinner{|{\psi(x)}\rangle}+\mathinner{|{\psi(y)}\rangle}\big{)}\right|$
(14) $\displaystyle\leq
2\|A\big{(}\mathinner{|{\psi(x)}\rangle}-\mathinner{|{\psi(y)}\rangle}\big{)}\|_{2}\leq
2\|A\|_{2}\|\mathinner{|{\psi(x)}\rangle}-\mathinner{|{\psi(y)}\rangle}\|_{2}$
(15)
$\displaystyle=2\|A\|_{2}\|(U_{x}-U_{y})\mathinner{|{\psi_{k}}\rangle}\|_{2}=2\|A\|_{2}\|(U_{y}^{\dagger}U_{x}-\mathds{1})\mathinner{|{\psi_{k}}\rangle}\|_{2},$
(16)

where we used unitary invariance of the vector norm, $U_{x}(s)=e^{-ixsH_{k}}$
and $U_{y}(s)=e^{-iysH_{k}}$, using the short-hand notation $U_{x}(1)=U_{x}$,
$U_{y}(1)=U_{y}$, and $\|A\|_{2}$ and $\|\mathinner{|{\psi}\rangle}\|_{2}$
denote the spectral norm of $A$ and the Euclidean vector norm of
$\mathinner{|{\psi}\rangle}$, respectively. Since

$\displaystyle
U_{y}^{\dagger}U_{x}-\mathds{1}=i\int_{0}^{1}ds\,U_{y}^{\dagger}(s)(y-x)H_{k}U_{x}(s),$
(17)

we have

$\displaystyle\|(U_{y}^{\dagger}U_{x}-\mathds{1})\mathinner{|{\psi_{k}}\rangle}\|_{2}$
$\displaystyle\leq\int_{0}^{1}ds\,\|(U_{y}^{\dagger}(s)(y-x)H_{k}U_{x}(s))\mathinner{|{\psi_{k}}\rangle}\|_{2}$
(18)
$\displaystyle=|x-y|\int_{0}^{1}ds\,\|H_{k}\|_{2}\|\mathinner{|{\psi(x)}\rangle}\|_{2}=|x-y|\|H_{k}\|_{2}=|x-y|\|H\|_{2}.$
(19)

In the last step we used unitary invariance of the spectral norm, assuming
that $H_{k}=V_{k}^{\dagger}HV_{k}$ where $V_{k}$ is unitary and $H$ is a
Hermitian matrix. Since

$\displaystyle\|A\|_{2}=\|[H_{k},H_{p}]\|_{2}\leq
2\|H_{k}\|_{2}\|H_{p}\|_{2}=2\|H_{p}\|_{2}\|H\|_{2},$ (20)

we arrive at

$\displaystyle\frac{\left|\frac{\partial J}{\partial x}-\frac{\partial
J}{\partial y}\right|}{|x-y|}\leq 4\|H_{p}\|_{2}\|H\|_{2}^{2},$ (21)

from which we conclude that the Lipschitz constant is given by
$L_{f}=4\|H_{p}\|_{2}\|H\|_{2}^{2}$. Further assuming $\|H\|_{2}=1$, we arrive
at the lower bound (4) given in the paper.

##  Appendix B Expected cost function change when sampling from unitary
2-designs

We first note that the expected value of the gradient
$\mathbb{E}_{V_{k}}\langle\text{grad}J[U_{k}],iH_{k}\rangle$ is given by

$\displaystyle\mathbb{E}_{V_{k}}\langle\text{grad}J[U_{k}],iV_{k}^{\dagger}HV_{k}\rangle=\langle\text{grad}J[U_{k}],i\mathds{1}\rangle\frac{\text{Tr}\\{H\\}}{d}=0,$
(22)

when $V_{k}$ is sampled according to the Haar measure on $\text{SU}(d)$ or
from a unitary 2-design. Consequently, in this case
$\mathbb{E}_{V_{k}}\langle\text{grad}J[U_{k}],iH_{k}\rangle^{2}$ is equal to
the variance.

In order to lower bound the expected cost function change
$\mathbb{E}_{V_{k}}\Delta J_{k}$ when $V_{k}$ is drawn according to the Haar
measure on $\text{SU}(d)$ or sampled from a unitary 2-design, we first note
that from the lower bound (4) we have $\mathbb{E}_{V_{k}}\Delta
J_{k}\geq\frac{1}{8\|H_{p}\|_{2}}\mathbb{E}_{V_{k}}\langle\text{grad}J[U_{k}],iH_{k}\rangle^{2}=\frac{1}{8\|H_{p}\|_{2}}\mathbb{E}_{V_{k}}\mathinner{\langle{\psi_{k}}|}[iV_{k}^{\dagger}HV_{k},H_{p}]\mathinner{|{\psi_{k}}\rangle}^{2}$.To
calculate the expectation on the right-hand side, we make use of the relation
[42],

$\displaystyle\int_{V\in\text{SU}(d)}\,dV\,\text{Tr}\\{V^{\dagger}BVCV^{\dagger}DVA\\}$
$\displaystyle=\frac{1}{d^{2}-1}\left(\text{Tr}\\{A\\}\text{Tr}\\{C\\}\text{Tr}\\{BD\\}+\text{Tr}\\{AC\\}\text{Tr}\\{B\\}\text{Tr}\\{D\\}\right)$
$\displaystyle-\frac{1}{d(d^{2}-1)}\left(\text{Tr}\\{AC\\}\text{Tr}\\{BD\\}+\text{Tr}\\{A\\}\text{Tr}\\{B\\}\text{Tr}\\{C\\}\text{Tr}\\{D\\}\right).$
(23)

We find,

$\displaystyle\mathbb{E}_{V_{k}}\mathinner{\langle{\psi_{k}}|}[iV_{k}^{\dagger}HV_{k},H_{p}]\mathinner{|{\psi_{k}}\rangle}^{2}$
$\displaystyle=\mathbb{E}_{V_{k}}\text{Tr}\\{[iV_{k}^{\dagger}HV_{k},H_{p}]\mathinner{|{\psi_{k}}\rangle}\mathinner{\langle{\psi_{k}}|}[iV_{k}^{\dagger}HV_{k},H_{p}]\mathinner{|{\psi_{k}}\rangle}\mathinner{\langle{\psi_{k}}|}\\},$
$\displaystyle=-2\text{Tr}\\{H^{2}\\}\left(\frac{1}{d^{2}-1}\mathinner{\langle{\psi_{k}}|}H_{p}\mathinner{|{\psi_{k}}\rangle}^{2}-\frac{1}{d(d^{2}-1)}\mathinner{\langle{\psi_{k}}|}H_{p}\mathinner{|{\psi_{k}}\rangle}^{2}\right.$
$\displaystyle\left.-\frac{1}{d^{2}-1}\mathinner{\langle{\psi_{k}}|}H_{p}^{2}\mathinner{|{\psi_{k}}\rangle}+\frac{1}{d(d^{2}-1)}\mathinner{\langle{\psi_{k}}|}H_{p}\mathinner{|{\psi_{k}}\rangle}^{2}\right)$
$\displaystyle=\frac{2\text{Tr}\\{H^{2}\\}}{d^{2}-1}\left(\mathinner{\langle{\psi_{k}}|}H_{p}^{2}\mathinner{|{\psi_{k}}\rangle}-\mathinner{\langle{\psi_{k}}|}H_{p}\mathinner{|{\psi_{k}}\rangle}^{2}\right),$
(24)

which yields the desired result given in (6) in the paper.