# Training variational quantum algorithms with random gate activation

Shuo Liu  Institute for Advanced Study, Tsinghua University, Beijing 100084,
China  Shi-Xin Zhang<EMAIL_ADDRESS>Tencent Quantum Laboratory,
Tencent, Shenzhen, Guangdong 518057, China  Shao-Kai Jian<EMAIL_ADDRESS>Department of Physics & Engineering Physics, Tulane University, New Orleans,
LA 70118, USA  Hong Yao<EMAIL_ADDRESS>Institute for Advanced Study,
Tsinghua University, Beijing 100084, China

###### Abstract

Variational quantum algorithms (VQAs) hold great potentials for near-term
applications and are promising to achieve quantum advantage on practical
tasks. However, VQAs suffer from severe barren plateau problem as well as have
a large probability of being trapped in local minima. In this Letter, we
propose a novel training algorithm with random quantum gate activation for
VQAs to efficiently address these two issues. This new algorithm processes
effectively much fewer training parameters than the conventional plain
optimization strategy, which efficiently mitigates barren plateaus with the
same expressive capability. Additionally, by randomly adding two-qubit gates
to the circuit ansatz, the optimization trajectories can escape from local
minima and reach the global minimum more frequently due to more sources of
randomness. In real quantum experiments, the new training algorithm can also
reduce the quantum computational resources required and be more quantum noise
resilient. We apply our training algorithm to solve variational quantum
simulation problems for ground states and present convincing results that
showcase the advantages of our novel strategy where better performance is
achieved by the combination of mitigating barren plateaus, escaping from local
minima, and reducing the effect of quantum noises. We further propose that the
entanglement phase transition could be one underlying reason why our RA
training is so effective.

Introduction: Recently, various quantum-classical hybrid variational
algorithms, such as variational quantum eigensolver (VQE) [1, 2], quantum
approximate optimization algorithm (QAOA) [3], and quantum neural networks
[4], have been proposed with the vision of establishing valuable quantum
killer-apps in Noisy Intermediate-Scale Quantum (NISQ) era [5]. In the
standard VQA setting, we minimize the expectation value of the objective
function $O$ designed for the target problem with respect to a variational
state (ansatz) $|\psi(\vec{\theta})\rangle$ prepared by a parameterized
quantum circuit (PQC). This hybrid scheme is accomplished with a feedback loop
between classical computers and quantum devices: the quantum devices
repeatedly prepare the variational state $|\psi(\vec{\theta})\rangle$ for
estimation of the expectation value of the objective function
$\langle\psi(\vec{\theta})|O|\psi(\vec{\theta})\rangle$, while the classical
computers are utilized to optimize the parameters $\vec{\theta}$ based on
classical optimization strategies such as gradient descent.

To guarantee that the variational quantum algorithm (VQA) solution is close
enough to the exact solution for the target problem, the high expressibility
of quantum ansatz $U$ is required. In principle, we can increase the depth of
the PQC with more quantum gates and training parameters or apply neural
network post-processing modules [6, 7] to achieve a higher expressibility. For
VQAs on real quantum hardware, the trainability is also an important factor.
Unfortunately, there is a trade-off between expressibility and trainability
[8] in VQA, and the performance of VQAs is severely limited by optimization
issues such as barren plateaus [9], which can be induced by the entanglement
in the quantum circuit [10, 11, 12] and the noise on NISQ devices [13]. The
gradients vanish exponentially with the depth of PQC. Therefore, exponential
computational resources are required to accurately estimate the gradients
value and a large number of iterations is required for the VQA solution to
converge. Even the gradient-free optimization approaches are also suppressed
by the barren plateaus [14]. In order to mitigate the barren plateaus and
achieve better performance from VQA, a series of strategies have been
proposed, including parameters initialization methods [15, 16, 17, 18, 19, 20,
21, 22, 23], local objective function [24, 25] and special quantum circuit
ansatz [26, 27, 28, 29, 30, 31].

Apart from the barren plateaus, the non-convexity of energy landscapes and the
existence of a large number of local minima also strongly limit the
trainability of VQAs [32] and render training VQAs unscalable. The VQA
solutions provided by gradient descent can be easily trapped in local minima
which are correlated with the initialization position and far away from the
global minimum [33]. To escape from the local minima, exponential trials with
random initialization of parameters, i.e. exponential optimization
trajectories, are required in the general case [34]. To avoid the local
minima, strategies utilize classical neural network [35], unitary block
optimization [36] and quantum dropout [37] have been investigated. To better
improve the trainability of VQAs, more effective and universal strategies are
still required.

In this Letter, we propose a training algorithm for VQAs with trainable gates
activated randomly and progressively to overcome the above challenges. We use
random activation (RA) to represent this strategy. In our new approach, the
number of trainable quantum gates is effectively much lower, and thus the
barren plateaus can be efficiently mitigated with the expressibility
unchanged. Additionally, the randomness of activating quantum gates increases
the probability of escaping from the local minima. Moreover, the quantum
computational resources required as well as the negative effect of quantum
noise can be both greatly reduced when experimentally realizing VQA on real
quantum hardware.

Without loss of generality, we focus on VQE task in this work, while our
approach is applicable to other types of VQAs. The VQE has been exploited in a
variety of contexts from quantum chemistry [38, 39, 40, 41, 42, 43, 44], many-
body physics [45, 46, 47, 48, 49, 50, 51, 52], to lattice gauge theories [53,
54] and is specified by a triplet ($O$, $|\psi_{i}\rangle$, $U(\vec{\theta})$)
including an objective function $O$, an initial quantum state
$|\psi_{i}\rangle$, and a PQC ansatz $U(\vec{\theta})$. The objective function
is usually chosen as the expectation or the variance [55, 56, 57] of the
system Hamiltonian operator $H$, whose solution gives the ground state or the
excited state, respectively. There are many choices for PQC structures such as
hardware-efficient ansatz (HEA) [58], Hamiltonian variational ansatz (HVA)
[59, 60] and the unitary coupled cluster (UCC) ansatz [61, 62, 63, 64, 65, 66,
67, 68].

Hamiltonian and circuit ansatz: In this Letter, we utilize the new proposed
training algorithm to solve the ground state energy problem of the one-
dimensional (1D) antiferromagnetic XXZ model, a representative lattice spin
model in quantum many-body physics. With convincing numerical results and
detailed ablation studies, we demonstrate the improved optimization results
using our approach and pin down the contributing factors as the mitigation of
the barren plateaus and a better chance to reach the global minimum.

The Hamiltonian of 1D antiferromagnetic XXZ model with periodic boundary
condition reads:

$\displaystyle H=\sum_{i}X_{i}X_{i+1}+Y_{i}Y_{i+1}+J_{z}Z_{i}Z_{i+1},$ (1)

where $X,Y,Z$ are Pauli matrices and $J_{z}$ is the $zz$ interaction strength.
We choose HVA for the PQC ansatz in this Letter. Specifically, a depth-$l$ HVA
for XXZ model is shown in Fig. 1 and the corresponding circuit unitary is:

$\displaystyle
U(\vec{\theta})=\prod_{i=1}^{l}R^{i,\mathrm{odd}}_{zz}R^{i,\mathrm{odd}}_{yy}R^{i,\mathrm{odd}}_{xx}R^{i,\mathrm{even}}_{zz}R^{i,\mathrm{even}}_{yy}R^{i,\mathrm{even}}_{xx},$
(2)

where
$R^{i,\mathrm{odd}}_{\sigma\sigma}=\prod_{j}\exp(i\theta_{i,2j+1,\sigma}\sigma_{2j+1}\sigma_{2j+2})$
and
$R^{i,\mathrm{even}}_{\sigma\sigma}=\prod_{j}\exp(i\theta_{i,2j,\sigma}\sigma_{2j}\sigma_{2j+1})$
$(\sigma=x,y,z)$. The parameters $\theta$ of quantum gates are independent
tuning parameters.

Figure 1: The HVA PQC for XXZ model with $4$ qubits: red rectangles represent
the activated two-qubit gates in which the parameters are updated in VQE
iterative optimization; white rectangles represent the unactivated two-qubit
gates which can be regarded as identity and no need to implement in real
experiments. For the training procedure, firstly, we randomly activate $10\%$
two-qubit gates with random initial parameters $\theta_{i}\in[0,2\pi]$, while
other unactivated two-qubit gates’ parameters remain as $0$. During training,
more and more two-qubit gates are randomly activated and their parameters will
unfreeze and update. $l$ is the PQC depth, i.e. how many times the Hamiltonian
block repeats. Figure 2: The converged VQE energy in $N=12$ system for the RA
training (green) and plain strategy (blue) with different PQC depths $l$: (a),
$J_{z}=0.5$; (b), $J_{z}=1.0$; (c), $J_{z}=2.0$. The black star represents the
averaged VQE energy and the red line represents the median VQE energy across
$500$ independent optimization trials. The inset is the zoom-in of VQE results
with deeper PQC. The outliers beyond the caps, which are much larger than the
ground truth are not shown (see the Supplemental Material for more details
about the standard box plots). The performance of RA training is much better:
the mean and median of VQE energy from RA is substantially lower than that
obtained from plain training, as more trials are trapped in the local minima
for plain training.

Applying the PQC $U(\vec{\theta})$ to the initial prepared state
$|\psi_{i}\rangle$ ($|\Psi^{-}\rangle$ Bell state, i.e.
$\bigotimes_{i=0}^{\frac{N}{2}-1}\frac{1}{\sqrt{2}}(|01\rangle-|10\rangle)_{2i,2i+1}$),
we obtain the output state
$|\psi(\vec{\theta})\rangle=U(\vec{\theta})|\psi_{i}\rangle$ and the VQE
energy $E(\vec{\theta})=\langle\psi(\vec{\theta})|H|\psi(\vec{\theta})\rangle$
as the objective function. For plain training, all two-qubit gates are
activated from the beginning of optimization with random initial parameters
$\vec{\theta}$ sampled uniformly from $[0,2\pi]$. Then the parameters
$\vec{\theta}$ are optimized by gradient descent to reach the lowest energy
$\langle\psi(\vec{\theta}^{*})|H|\psi(\vec{\theta}^{*})\rangle$, where
$\vec{\theta}^{*}$ are the optimal parameters. In principle, a deeper PQC has
higher expressibility and a lower converged VQE energy. However, barren
plateau problems become more severe with large PQC depth, posing a fundamental
challenge to identify the optimal parameters $\vec{\theta}^{*}$ as discussed
above.

Training with random gate activation: To overcome the challenges presented in
the VQA training procedure, we introduce the training algorithm with
incrementally random gate activation. Different from the plain training
method, only a small fraction of two-qubit gates, for example, $10\%$, are
activated for the initial optimization iterations (i.e. the two-qubit gates
with the red color shown in Fig. 1). Meanwhile, other two-qubit gates stay as
identity gates, i.e. $\theta=0$ for these unactivated gates. During the
optimization, the parameters $\vec{\theta}_{10\%}$ of the initial $10\%$
activated gates will be updated to approach the optimal parameters
$\vec{\theta}^{*}_{10\%}$ corresponding to the lowest VQE energy estimation
$\langle\psi(\vec{\theta}^{*}_{10\%})|H|\psi(\vec{\theta}^{*}_{10\%})\rangle\leq\langle\psi(\vec{\theta}_{10\%})|H|\psi(\vec{\theta}_{10\%})\rangle$.
Then other random $10\%$ unactivated two-qubit gates are activated. Since the
parameters for these newly activated gates are set to $0$ previously, the VQE
energy has no sudden change for the incremental activation procedure. Now the
parameters of the $20\%$ activated gates are optimized together, and this
procedure repeats until all parameterized gates are activated, and the optimal
parameters $\vec{\theta}_{100\%}^{*}$ are obtained. This approach has
effectively fewer parameterized gates and can hopefully mitigate the barren
plateaus with better VQE performance.

VQE results comparison: We provide numerical analyses on the performance of
different training strategies, with numerical implementation based on
TensorCircuit package [69]. The results of the plain training and the RA
training are shown in Fig. 2 with varying PQC depths and Hamiltonian
parameters. The results consist of $500$ independent optimization trials on a
$12$-qubit system. We use Adam [70] optimizer with hyperparameters
learning_rate=0.01, decay_rate=0.9, decay_steps=100, and the percentage of
two-qubit gates activated each time is $10\%$ (results are similar for other
reasonable hyperparameter choices and details can be found in the Supplemental
Material). We keep the maximal number of iterations maxiter=5000 which is
large enough to guarantee that the two different training approaches achieve
their best performance. This is different from the usual approach where the
converge speed instead of the final converged position is the focus [27, 30,
31]. We believe our comparison is more suitable to exploit the potential
quantum advantage of VQAs.

For the RA training, there are much fewer trials trapped in the bad local
minima and the averaged VQE energy is much lower than that obtained from the
plain training as shown in Fig. 2. There are two sources for the advantageous
performance for the RA training approach, and we will investigate the two
factors in detail below: escaping from the local minima more frequently and
mitigating the barren plateaus.

The improvement of averaged VQE energy exists for different PQC depths and is
particularly significant for shallower quantum circuits with $l=1,2$. Note
that the barren plateaus are not severe in such a shallow setup. Considering
an optimization trajectory of the plain training with random initial
parameters $\vec{\theta}$, the converged local minimum of the trajectory is
soled determined by the initialization [33] for plain training. And as long as
the trajectory is trapped in a local minimum, the VQE performance can be
rather bad, providing local minima far away from the exact ground state.
Instead, if we randomly activate the gate during optimization, for the fixed
initial parameters, there are multiple different random optimization
trajectories caused by the randomness of how to activate two-qubit gates. The
extra randomness sources greatly enlarge the parameter space that can be
explored for the optimization trajectories, and efficiently help the
optimization trajectories escape from local minima more frequently.

Since the effective circuit depth of our algorithm with only a fraction number
of gates is much lower than that of a plain training algorithm, our algorithm
also helps to mitigate barren plateaus. To observe this quantitatively, the
averaged variance of the VQE energy gradients with respect to the parameters
of the activated gates are shown in Fig. 3. With the number of activated gates
and the system size increase, the averaged variance decreases exponentially.
The case of $100\%$ density corresponds to the plain VQE training where all
gates are activated. And such a plain training strategy suffers from the most
severe barren plateau problems. As shown in Fig. 3, the gradient variance can
be efficiently increased during the RA training procedure. The mitigation of
barren plateaus can also be confirmed via the better median VQE energy of the
training. The median energy is not affected by the outlier local minima trials
and the improvement on the median is more likely a result from mitigating
barren plateaus instead of escaping from a local minimum. In sum, the
contribution from escaping local minimum can be found via the mean of VQE
energy where the distribution of converged energy matters, and the
contribution from mitigating barren plateaus can be identified via the median
of VQE energy.

Figure 3: The averaged variance of energy gradients with respect to the
parameters of the activated two-qubit gates. The depth of PQC is $l=7$ and the
Hamiltonian parameter $J_{z}=1.0$. The barren plateau problems become worse
with more activated gates.

Layerwise activation: We carry out further ablation studies to confirm our
analysis on the advantages of RA training. A straightforward question is what
will happen if there is no randomness in the gate activation procedure, e.g.,
activating the two-qubit gates in a layerwise fashion. The barren plateaus can
also be mitigated due to the effective lower depth for layerwise activation
idea. But the extra randomness source for escaping from local minima now
disappears. We consider two classes of layerwise optimization: one is the
training with layerwise appending activation (LAA), namely, we append the new
activated identity layer to the end of the activated layers; the other is the
training with layerwise prepending activation (LPA), namely, we prepend the
new identity layer to the activated ones.

The distribution of the converged VQE results for different training
strategies (Plain, RA, LAA, LPA) are shown in Fig. 4. When the quantum circuit
is shallow, almost all the trajectories of RA training successfully escape
from the local minima, while many trials of the other three training
strategies are trapped in local minima and lead to a larger fluctuation in the
results due to the absence of the randomness in the activation. Therefore, the
randomness of activating two-qubit gates is vital and relevant for better
escaping from local minima. When the quantum circuit is deep, the VQE
performance of RA training is still much better than that of the plain
training as discussed above. We have also observed that the VQE performance of
RA and LPA are very close. This result indicates that the improvement of VQE
performance with deep circuits is mainly due to the mitigation of barren
plateaus where RA and LPA are both capable of.

Figure 4: The accumulated distribution of $500$ converged energies from
independent optimization trials using four training strategies based on the
Adam optimizer with decay_steps=100 in $N=12$ system: (a), $J_{z}=0.5$, $l=2$;
(b), $J_{z}=0.5$, $l=7$; (c), $J_{z}=1.0$, $l=2$; (d), $J_{z}=1.0$, $l=7$;
(e), $J_{z}=2.0$, $l=2$; (f), $J_{z}=2.0$, $l=7$. When the quantum circuit is
shallow, the VQE performance of RA training is the best because this strategy
can escape from local minima more easily. When the quantum circuit is deep,
the VQE performance of RA and LPA is very close for $J_{z}=0.5$ and
$J_{z}=1.0$ and the VQE performance of RA is the best over $4$ different
training strategies for $J_{z}=2.0$.

Entanglement phase transition in RA training: A possible reason why RA
training works so well even if we consider a fraction of the quantum gates can
be demystified from entanglement phase transition perspective. We take the HVA
circuit structure utilized in this work as the model and the ratio of
activated gates $p$ as the driven parameter. We can then identify an
entanglement phase transition for the output state from an area law to a
volume law, similar to the measurement induced entanglement phase transition
[71, 72, 73, 74, 75, 76, 77]. In our setup, the activation ratio $p$ is
relevant, namely $p_{c}=0$. This fact indicates that even with only 10% gate
activated in this circuit, the output state still lies in the same phase of
the full circuit output. This entanglement phase transition perspective may
partially explain why the RA training is effective. See more theoretical and
numerical details in the Supplemental Material.

Discussions: We have demonstrated that the RA training can efficiently
mitigate the barren plateau problem and help the optimization trajectories
escape from local minima based on numerical results and theoretical analysis.
This new algorithm can be naturally applied to other practical tasks and other
flavors of VQAs due to its universal form. Additionally, it can be integrated
with other strategies [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28,
29, 30, 35, 36, 78, 79] to further improve performance and solution quality of
VQAs.

From a practical point of view, the new training algorithm brings more
benefits when considering implementation on real quantum hardware in the NISQ
era. In experiments, the circuit gradients are evaluated via parameter shift
rule [80, 81, 82, 83]. For a PQC with $p$ training parameters, we have to
evaluate the observable expectation on $2p$ sets of parameters. In RA training
scheme, the effective number of training parameters during the optimization is
much smaller than the plain training, which greatly reduces the required
number of measurement shots. Moreover, the required number of measurement
shots for each observable is also lowered due to the mitigated barren
plateaus. Meanwhile, the number of two-qubit gates required to implement on
the hardware in RA training is also lower than the plain training case.
Therefore, the efficiency in terms of the practical quantum computational
resources required is a big advantage for RA training approach. And for the
experiments we conducted in this Letter, roughly 500 times less quantum gates
are used for RA training (see the Supplemental Material for more detail about
the quantum resource efficiency analysis).

In terms of hardware implementation, RA training also has better noise
resilience. On the one hand, we can stop the activation iteration earlier to
achieve a better trade-off between the circuit expressive power and the
accumulated quantum noise. On the other hand, the effectively shallower
circuits also suppress barren plateaus induced by the quantum noise [13]. We
show experimental results on noisy VQE settings, and the obtained energies are
substantially improved when using RA training (the numerical results and data
can be found in the Supplemental Material). Therefore, RA training is a
practical and standalone training approach that we recommend trying on noisy
quantum devices.

Acknowledgements: We thank Zhou-Quan Wan for helpful discussions. This work is
supported in part by the NSFC under Grant No. 11825404 (SXZ, SL, and HY) and
by the MOSTC under the MOSTC Grants No. 2021YFA1400100 and No. 2018YFA0305604
(HY). S.-K.J is supported by a startup fund at Tulane University.

## References

  * Peruzzo  _et al._ [2014a] 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, Nature Communications 5, 4213 (2014a).
  * McClean  _et al._ [2016a] J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, The theory of variational hybrid quantum-classical algorithms, New Journal of Physics 18, 023023 (2016a).
  * Farhi  _et al._ [2014] E. Farhi, J. Goldstone, and S. Gutmann, A Quantum Approximate Optimization Algorithm, arXiv:1411.4028 (2014).
  * Farhi and Neven [2018] E. Farhi and H. Neven, Classification with quantum neural networks on near term processors, arXiv:1802.06002 (2018).
  * Preskill [2018] J. Preskill, Quantum Computing in the NISQ era and beyond, Quantum 2, 79 (2018).
  * Zhang  _et al._ [2022a] S.-X. Zhang, Z.-Q. Wan, C.-K. Lee, C.-Y. Hsieh, S. Zhang, and H. Yao, Variational quantum-neural hybrid eigensolver, Phys. Rev. Lett. 128, 120502 (2022a).
  * Zhang  _et al._ [2021a] S.-X. Zhang, Z.-Q. Wan, C.-Y. Hsieh, H. Yao, and S. Zhang, Variational Quantum-Neural Hybrid Error Mitigation, arXiv:2112.10380 (2021a).
  * Holmes  _et al._ [2022] Z. Holmes, K. Sharma, M. Cerezo, and P. J. Coles, Connecting ansatz expressibility to gradient magnitudes and barren plateaus, PRX Quantum 3, 010313 (2022).
  * 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, Nature Communications 9, 4812 (2018).
  * Patti  _et al._ [2021] T. L. Patti, K. Najafi, X. Gao, and S. F. Yelin, Entanglement devised barren plateau mitigation, Phys. Rev. Research 3, 033090 (2021).
  * Sharma  _et al._ [2022] K. Sharma, M. Cerezo, L. Cincio, and P. J. Coles, Trainability of dissipative perceptron-based quantum neural networks, Phys. Rev. Lett. 128, 180505 (2022).
  * 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._ [2021] 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, 6961 (2021).
  * Arrasmith  _et al._ [2021] A. Arrasmith, M. Cerezo, P. Czarnik, L. Cincio, and P. J. Coles, Effect of barren plateaus on gradient-free optimization, Quantum 5, 558 (2021).
  * Grant  _et al._ [2019] E. Grant, L. Wossnig, M. Ostaszewski, and M. Benedetti, An initialization strategy for addressing barren plateaus in parametrized quantum circuits, Quantum 3, 214 (2019).
  * Verdon  _et al._ [2019] G. Verdon, M. Broughton, J. R. McClean, K. J. Sung, R. Babbush, Z. Jiang, H. Neven, and M. Mohseni, Learning to learn with quantum neural networks via classical neural networks, arXiv:1907.05415 (2019).
  * Liu  _et al._ [2021] H.-Y. Liu, T.-P. Sun, Y.-C. Wu, Y.-J. Han, and G.-P. Guo, A parameter initialization method for variational quantum algorithms to mitigate barren plateaus based on transfer learning, arXiv:2112.10952 (2021).
  * Rad  _et al._ [2022] A. Rad, A. Seif, and N. M. Linke, Surviving the barren plateau in variational quantum circuits with bayesian learning initialization, arXiv:2203.02464 (2022).
  * Kulshrestha and Safro [2022] A. Kulshrestha and I. Safro, Beinit: Avoiding barren plateaus in variational quantum algorithms, arXiv:2204.13751 (2022).
  * Sauvage  _et al._ [2021] F. Sauvage, S. Sim, A. A. Kunitsa, W. A. Simon, M. Mauri, and A. Perdomo-Ortiz, Flip: A flexible initializer for arbitrarily-sized parametrized quantum circuits, arXiv:2103.08572 (2021).
  * Sack  _et al._ [2022] S. H. Sack, R. A. Medina, A. A. Michailidis, R. Kueng, and M. Serbyn, Avoiding barren plateaus using classical shadows, arXiv:2201.08194 (2022).
  * Grimsley  _et al._ [2022] H. R. Grimsley, G. S. Barron, E. Barnes, S. E. Economou, and N. J. Mayhall, Adapt-vqe is insensitive to rough parameter landscapes and barren plateaus, arXiv:2204.07179 (2022).
  * Zhang  _et al._ [2022b] K. Zhang, M.-H. Hsieh, L. Liu, and D. Tao, Gaussian initializations help deep variational quantum circuits escape from the barren plateau, arXiv:2203.09376 (2022b).
  * Cerezo  _et al._ [2021] M. Cerezo, A. Sone, T. Volkoff, L. Cincio, and P. J. Coles, Cost function dependent barren plateaus in shallow parametrized quantum circuits, Nature Communications 12, 1791 (2021).
  * Uvarov and Biamonte [2021] A. V. Uvarov and J. D. Biamonte, On barren plateaus and cost function locality in variational quantum algorithms, Journal of Physics A: Mathematical and Theoretical 54, 245301 (2021).
  * Pesah  _et al._ [2021] A. Pesah, M. Cerezo, S. Wang, T. Volkoff, A. T. Sornborger, and P. J. Coles, Absence of barren plateaus in quantum convolutional neural networks, Phys. Rev. X 11, 041011 (2021).
  * Liu  _et al._ [2022a] X. Liu, G. Liu, J. Huang, and X. Wang, Mitigating barren plateaus of variational quantum eigensolvers, arXiv:2205.13539 (2022a).
  * Bilkis  _et al._ [2021] M. Bilkis, M. Cerezo, G. Verdon, P. J. Coles, and L. Cincio, A semi-agnostic ansatz with variable structure for quantum machine learning, arXiv:2103.06712 (2021).
  * Du  _et al._ [2022] Y. Du, T. Huang, S. You, M.-H. Hsieh, and D. Tao, Quantum circuit architecture search for variational quantum algorithms, npj Quantum Information 8, 62 (2022).
  * Zhang  _et al._ [2021b] K. Zhang, M.-H. Hsieh, L. Liu, and D. Tao, Toward Trainability of Deep Quantum Neural Networks, arXiv:2112.15002 (2021b).
  * Anand  _et al._ [2022] A. Anand, L. B. Kristensen, F. Frohnert, S. Sim, and A. Aspuru-Guzik, Information flow in parameterized quantum circuits, arXiv:2207.05149 (2022).
  * Bittel and Kliesch [2021] L. Bittel and M. Kliesch, Training variational quantum algorithms is np-hard, Phys. Rev. Lett. 127, 120502 (2021).
  * Anschuetz [2022] E. R. Anschuetz, Critical Points in Quantum Generative Models, arXiv:2109.06957 (2022).
  * Zhou  _et al._ [2020] L. Zhou, S.-T. Wang, S. Choi, H. Pichler, and M. D. Lukin, Quantum approximate optimization algorithm: Performance, mechanism, and implementation on near-term devices, Phys. Rev. X 10, 021067 (2020).
  * Rivera-Dean  _et al._ [2021] J. Rivera-Dean, P. Huembeli, A. Acín, and J. Bowles, Avoiding local minima in Variational Quantum Algorithms with Neural Networks, arXiv:2104.02955 (2021).
  * Slattery  _et al._ [2022] L. Slattery, B. Villalonga, and B. K. Clark, Unitary block optimization for variational quantum algorithms, Phys. Rev. Research 4, 023072 (2022).
  * Wang  _et al._ [2022] Z.-D. Wang, P.-L. Zheng, B. Wu, and Y. Zhang, Quantum Dropout for Efficient Quantum Approximate Optimization Algorithm on Combinatorial Optimization Problems, arXiv:2203.10101 (2022).
  * O’Malley  _et al._ [2016] P. J. J. O’Malley, R. Babbush, I. D. Kivlichan, J. Romero, J. R. McClean, R. Barends, J. Kelly, P. Roushan, A. Tranter, N. Ding, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. G. Fowler, E. Jeffrey, E. Lucero, A. Megrant, J. Y. Mutus, M. Neeley, C. Neill, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. C. White, P. V. Coveney, P. J. Love, H. Neven, A. Aspuru-Guzik, and J. M. Martinis, Scalable quantum simulation of molecular energies, Phys. Rev. X 6, 031007 (2016).
  * Kandala  _et al._ [2017a] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets, Nature 549, 242 (2017a).
  * Colless  _et al._ [2018] J. I. Colless, V. V. Ramasesh, D. Dahlen, M. S. Blok, M. E. Kimchi-Schwartz, J. R. McClean, J. Carter, W. A. de Jong, and I. Siddiqi, Computation of molecular spectra on a quantum processor with an error-resilient algorithm, Phys. Rev. X 8, 011021 (2018).
  * 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, Chemical reviews 119, 10856 (2019).
  * 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, Nature Communications 10, 3007 (2019).
  * 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).
  * Tilly  _et al._ [2020] J. Tilly, G. Jones, H. Chen, L. Wossnig, and E. Grant, Computation of molecular excited states on ibm quantum computers using a discriminative variational quantum eigensolver, Phys. Rev. A 102, 062425 (2020).
  * Liu  _et al._ [2019] J.-G. Liu, Y.-H. Zhang, Y. Wan, and L. Wang, Variational Quantum Eigensolver with Fewer Qubits, Physical Review Research 1, 023025 (2019).
  * Cai [2020] Z. Cai, Resource Estimation for Quantum Variational Simulations of the Hubbard Model, Physical Review Applied 14, 014059 (2020).
  * McClean  _et al._ [2016b] J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, The theory of variational hybrid quantum-classical algorithms, New Journal of Physics 18, 023023 (2016b).
  * Wecker  _et al._ [2015] D. Wecker, M. B. Hastings, and M. Troyer, Progress towards practical quantum variational algorithms, Physical Review A 92, 042303 (2015).
  * Ho and Hsieh [2019] W. W. Ho and T. H. Hsieh, Efficient variational simulation of non-trivial quantum states, SciPost Phys. 6, 29 (2019).
  * Dallaire-Demers  _et al._ [2018] P.-L. Dallaire-Demers, J. Romero, L. Veis, S. Sim, and A. Aspuru-Guzik, Low-depth circuit ansatz for preparing correlated fermionic states on a quantum computer, arXiv:1801.01053 (2018).
  * Tamiya  _et al._ [2021] S. Tamiya, S. Koh, and Y. O. Nakagawa, Calculating nonadiabatic couplings and Berry’s phase by variational quantum eigensolvers, Physical Review Research 3, 023244 (2021).
  * Uvarov  _et al._ [2020] A. Uvarov, J. Biamonte, and D. Yudin, Variational Quantum Eigensolver for Frustrated Quantum Systems, Physical Review B 102, 075104 (2020).
  * Peruzzo  _et al._ [2014b] 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, Nature Communications 5, 4213 (2014b).
  * Klco  _et al._ [2018] N. Klco, E. F. Dumitrescu, A. J. McCaskey, T. D. Morris, R. C. Pooser, M. Sanz, E. Solano, P. Lougovski, and M. J. Savage, Quantum-classical computation of schwinger model dynamics using quantum computers, Phys. Rev. A 98, 032331 (2018).
  * Zhang  _et al._ [2020] D.-B. Zhang, Z.-H. Yuan, and T. Yin, Variational quantum eigensolvers by variance minimization, arXiv:2006.15781 (2020).
  * Zhang  _et al._ [2021c] F. Zhang, N. Gomes, Y. Yao, P. P. Orth, and T. Iadecola, Adaptive variational quantum eigensolvers for highly excited states, Phys. Rev. B 104, 075159 (2021c).
  * Liu  _et al._ [2023] S. Liu, S.-X. Zhang, C.-Y. Hsieh, S. Zhang, and H. Yao, Probing many-body localization by excited-state variational quantum eigensolver, Phys. Rev. B 107, 024204 (2023).
  * Kandala  _et al._ [2017b] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets, Nature 549, 242 (2017b).
  * Wiersema  _et al._ [2020] R. Wiersema, C. Zhou, Y. de Sereville, J. F. Carrasquilla, Y. B. Kim, and H. Yuen, Exploring entanglement and optimization within the hamiltonian variational ansatz, PRX Quantum 1, 020319 (2020).
  * Li  _et al._ [2021] A. C. Li, M. S. Alam, T. Iadecola, A. Jahin, D. M. Kurkcuoglu, R. Li, P. P. Orth, A. B. Özgüler, G. N. Perdue, and N. M. Tubman, Benchmarking variational quantum eigensolvers for the square-octagon-lattice kitaev model, arXiv:2108.13375 (2021).
  * Kutzelnigg [1982] W. Kutzelnigg, Quantum chemistry in Fock space. I. The universal wave and energy operators, The Journal of Chemical Physics 77, 3081 (1982).
  * Kutzelnigg and Koch [1983] W. Kutzelnigg and S. Koch, Quantum chemistry in Fock space. II. Effective Hamiltonians in Fock space, The Journal of Chemical Physics 79, 4315 (1983).
  * Kutzelnigg [1985] W. Kutzelnigg, Quantum chemistry in Fock space. IV. The treatment of permutational symmetry. Spin‐free diagrams with symmetrized vertices, The Journal of Chemical Physics 82, 4166 (1985).
  * Bartlett  _et al._ [1989] R. J. Bartlett, S. A. Kucharski, and J. Noga, Alternative coupled-cluster ansätze II. The unitary coupled-cluster method, Chemical Physics Letters 155, 133 (1989).
  * Taube and Bartlett [2006] A. G. Taube and R. J. Bartlett, New perspectives on unitary coupled-cluster theory, International Journal of Quantum Chemistry 106, 3393 (2006).
  * Yanai  _et al._ [2009] T. Yanai, Y. Kurashige, D. Ghosh, and G. K.-L. Chan, Accelerating convergence in iterative solution for large-scale complete active space self-consistent-field calculations, International Journal of Quantum Chemistry 109, 2178 (2009).
  * Yanai and Chan [2006] T. Yanai and G. K.-L. Chan, Canonical transformation theory for multireference problems, The Journal of Chemical Physics 124, 194106 (2006).
  * Harsha  _et al._ [2018] G. Harsha, T. Shiozaki, and G. E. Scuseria, On the difference between variational and unitary coupled cluster theories, The Journal of Chemical Physics 148, 044107 (2018).
  * Zhang  _et al._ [2023] S.-X. Zhang, J. Allcock, Z.-Q. Wan, S. Liu, J. Sun, H. Yu, X.-H. Yang, J. Qiu, Z. Ye, Y.-Q. Chen, C.-K. Lee, Y.-C. Zheng, S.-K. Jian, H. Yao, C.-Y. Hsieh, and S. Zhang, TensorCircuit: a Quantum Software Framework for the NISQ Era, Quantum 7, 912 (2023).
  * Kingma and Ba [2017] D. P. Kingma and J. Ba, Adam: A Method for Stochastic Optimization, arXiv:1412.6980 (2017).
  * Li  _et al._ [2018] Y. Li, X. Chen, and M. P. A. Fisher, Quantum zeno effect and the many-body entanglement transition, Phys. Rev. B 98, 205136 (2018).
  * Skinner  _et al._ [2019] B. Skinner, J. Ruhman, and A. Nahum, Measurement-induced phase transitions in the dynamics of entanglement, Phys. Rev. X 9, 031009 (2019).
  * Li  _et al._ [2019] Y. Li, X. Chen, and M. P. A. Fisher, Measurement-driven entanglement transition in hybrid quantum circuits, Phys. Rev. B 100, 134306 (2019).
  * Choi  _et al._ [2020] S. Choi, Y. Bao, X.-L. Qi, and E. Altman, Quantum error correction in scrambling dynamics and measurement-induced phase transition, Phys. Rev. Lett. 125, 030505 (2020).
  * Gullans and Huse [2020] M. J. Gullans and D. A. Huse, Dynamical purification phase transition induced by quantum measurements, Phys. Rev. X 10, 041020 (2020).
  * Jian  _et al._ [2021] S.-K. Jian, C. Liu, X. Chen, B. Swingle, and P. Zhang, Measurement-induced phase transition in the monitored sachdev-ye-kitaev model, Phys. Rev. Lett. 127, 140601 (2021).
  * Liu  _et al._ [2022b] S. Liu, M.-R. Li, S.-X. Zhang, S.-K. Jian, and H. Yao, Universal KPZ scaling in noisy hybrid quantum circuits, arXiv:2212.03901 (2022b).
  * Zhang  _et al._ [2021d] S.-X. Zhang, C.-Y. Hsieh, S. Zhang, and H. Yao, Neural predictor based quantum architecture search, Machine Learning: Science and Technology 2, 045027 (2021d).
  * Zhang  _et al._ [2022c] S.-X. Zhang, C.-Y. Hsieh, S. Zhang, and H. Yao, Differentiable quantum architecture search, Quantum Science and Technology 7, 045023 (2022c).
  * 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).
  * Banchi and Crooks [2021] L. Banchi and G. E. Crooks, Measuring Analytic Gradients of General Quantum Evolution with the Stochastic Parameter Shift Rule, Quantum 5, 386 (2021).
  * Bottou [2012] L. Bottou, Stochastic Gradient Descent Tricks (Springer Berlin Heidelberg, Berlin, Heidelberg, 2012).

## Supplemental Materials

###  S1 Box plots

In this work, we use box plots extensively to show the distribution of
converged energies from different optimization trails. Box plots are used to
visualize summary statistics of a dataset, displaying attributes of the
distribution like the data’s range and distribution.

In our case, the dataset consists of the converged VQE energies for $N_{t}$
($N_{t}=500$ for noiseless case and $N_{t}=100$ for noisy case) independent
optimization trials. As shown in Fig. S1, the black star represents the
averaged VQE energy and the red line represents the median VQE energy across
$N_{t}$ independent optimization trials. The left edge and the right edge of
the box correspond to the first quartile ($Q_{1}$) and the third quartile
($Q_{3}$) respectively, where $Q_{1}$ is defined as the median VQE energy
between the best VQE energy (minimum) and the median VQE energy, and $Q_{3}$
is the median VQE energy between the median VQE energy and the worst VQE
energy (maximum). $Q_{3}-Q_{1}$ is known as the interquartile range (IQR).
After determining the first and third quartiles and the interquartile range as
outlined above, then fences are calculated using the following formula:

$\displaystyle\text{Lower fence}=Q_{1}-1.5\text{IQR},$ (S1)
$\displaystyle\text{Upper fence}=Q_{3}+1.5\text{IQR}.$ (S2)

Any data lying outside these defined bounds can be considered as an outlier.
If the best (worst) VQE energy is higher (lower) than the Lower (Upper) fence,
the left (right) bar is choosen as the best (worst) VQE energy; otherwise, the
left (right) bar is chosen as the best (worst) VQE energy inside the fences.
For Plain training, many trials are trapped in the local minima which are
higher than the Upper fence and thus not shown in the box plots, namely, the
right bar in the box plots corresponds to the worst VQE energy inside the
fences instead of the real worst VQE energy. The VQE energy distribution
including the worst results trapped in the local minima can be identified in
the plots of the accumulated distribution instead, e.g., Fig. 4. On the
contrary, due to the large IQR and the ground state energy lower bound, the
lower fence always represents the best VQE energy obtained in the experiments.

Figure S1: Different parts of a boxplot: the black star represents the
averaged VQE energy and the red line represents the median VQE energy across
$N_{t}$ independent optimization trials.

###  S2 Barren plateaus

The averaged gradient variances for different percentages of randomly
activated two-qubit gates with $J_{z}=0.5$ and $J_{z}=2.0$ are shown in Fig.
S2. The barren plateau problems become worse with more parameterized quantum
gates activated consistent with the results shown in Fig. 3 with $J_{z}=1.0$.

Figure S2: Demonstration of barren plateaus in the VQE problem. The averaged
variance of energy gradients with respect to the parameters of the activated
two-qubit gates are shown: (a), $J_{z}=0.5$; (b), $J_{z}=2.0$. The PQC depth
is $l=7$. The barren plateau problems become worse with more activated gates.

###  S3 Best VQE energy analysis

We also observe that the best VQE energy over $500$ independent trials is
similar for two different types of training when the quantum circuit depth is
shallow ($l=1,2$), while the best VQE energy of RA training is much lower than
that of the plain training when the quantum circuit depth is deep ($l=5,6,7$)
as shown in Fig. 2. Theoretically, the best VQE energy should be the same for
different strategies since the PQC ansatz is the same. However, due to the
local minima issue, the theoretical global minimum of the ansatz can only be
identified with exponential ($O(2^{l})$) parameters initialization trials for
plain training. When the quantum circuit is shallow, $500$ independent trials
are enough to reach this global minimum though many of the trials are trapped
in local minima. So the best VQE energy of two training algorithms is close to
each other as the theoretical minimum is reached in both cases. However, when
the quantum circuit is deep, $500$ independent trials are not enough for plain
training to explore the full energy landscape and reach the global minimum.
Even if there is a lucky trajectory that does not meet any local minima, the
global minimum can not be reached via a large number of iterations due to the
vanishing gradients. The improvement of the best VQE energy thus owes to both
of the contribution factors of our algorithm.

###  S4 Adam optimizer and different decap_steps

In the main text, we have fixed the optimizer type and the optimizer’s
hyperparameters. It is worth noting that the conclusions drawn in the main
text don’t rely on the special choice of optimizer and its hyperparameters.
The results for different decay_steps and fixed learning_rate=0.01,
decay_rate=0.9 are shown in Fig. S3 (the quantum circuit depth is $l=7$). The
learning rate on $t$ iterations is given by
learning_rate($t$)=learning_rate$*$decay_rate${}^{\frac{t}{\text{decay}\\_\text{steps}}}$.
The choices of decay_steps influence the decay speed of the learning rate and
thus the VQE optimization performance. However, RA training can consistently
achieve better performance than that of the plain training for different
decay_steps: the best VQE energy, averaged VQE energy and median VQE energy
are all much lower and closer to the ground state energy as shown in Fig. S3.
The results for different quantum circuit depth $l$ and fixed decay_steps=80
are shown in Fig. S4, which also gives similar results as discussed in the
main text with decay_steps=100. With fixed decay_steps=80, the distribution of
the results for different independent trials is shown in Fig. S5.

Figure S3: The VQE energy obtained from plain training (blue) and RA training
(green) on a $12$-qubits system with different decay_steps for Adam optimizer
and fixed PQC depth $l=7$: (a), $J_{z}=0.5$; (b), $J_{z}=1.0$; (c),
$J_{z}=2.0$. The black star represents the averaged VQE energy and the red
line represents the median VQE energy over $500$ independent trials. The
performance of RA training is much better than that of the plain training with
lower best VQE energy, lower averaged VQE energy, and lower median VQE energy.
When $J_{z}=1.0$ and decay_steps=60, the VQE performance of RA training is
worse, which may be caused by the fast-decaying learning rate of the
optimizer. Figure S4: The converged VQE energy for RA training (green) and
plain training (blue) with decay_steps=80 and different PQC depths $l$: (a),
$J_{z}=0.5$; (b), $J_{z}=1.0$; (c), $J_{z}=2.0$. The black star represents the
averaged VQE energy and the red line represents the median VQE energy across
$500$ independent optimization trials. The inset is the zoom-in of VQE results
with deeper PQC. The outliers beyond the caps, which are much larger than the
ground truth are not shown. The performance of RA training is much better: the
averaged VQE energy and median VQE energy from RA training are substantially
lower than that obtained from plain training. More trials are trapped in the
local minima for plain training. The results are similar to those shown in
Fig. 2 with decay_steps=100. Figure S5: The accumulated distribution of $500$
converged energies from independent optimization trials with four training
strategies with decay_steps=80: (a), $J_{z}=0.5$, $l=2$; (b), $J_{z}=0.5$,
$l=7$; (c), $J_{z}=1.0$, $l=2$; (d), $J_{z}=1.0$, $l=7$; (e), $J_{z}=2.0$,
$l=2$; (f), $J_{z}=2.0$, $l=7$. When the quantum circuit is shallow, the VQE
performance of RA training is the best because this strategy can mitigate
barren plateaus and escape from local minima more easily at the same time.
When the quantum circuit is deep, the VQE performance of RA and LPA is very
close.

###  S5 SGD optimizer

We have also tried different types of gradient descent optimizers such as
stochastic gradient descent (SGD) [84] optimizer to verify the availability of
RA training. The results for different $J_{z}$ and different quantum circuit
depths $l$ are shown in Fig. S6. The distribution of the results for different
independent trials is shown in Fig. S7. The conclusions remains similar to the
case with Adam optimizer discussed before.

Figure S6: The VQE energy for the plain training (blue) and the RA training
(green) with different quantum circuit depths $l$ and fixed learning_rate=0.01
based on SGD optimizer: (a), $J_{z}=0.5$; (b), $J_{z}=1.0$; (c), $J_{z}=2.0$.
The black star represents the averaged VQE energy and the red line represents
the median VQE energy. The inset is the zoom-in of VQE results with deeper
PQC. The performance of RA training is much better, similar to the case shown
in Fig. 2. Figure S7: The accumulated distribution of $500$ converged energies
from independent optimization trials with four training strategies with
learning_rate=0.01 based on SGD optimizer: (a), $J_{z}=0.5$, $l=2$; (b),
$J_{z}=0.5$, $l=7$; (c), $J_{z}=1.0$, $l=2$; (d), $J_{z}=1.0$, $l=7$; (e),
$J_{z}=2.0$, $l=2$; (f), $J_{z}=2.0$, $l=7$. When the quantum circuit is
shallow, the VQE performance of RA training is the best because this strategy
can mitigate barren plateaus and escape from local minima more easily at the
same time; when the quantum circuit is deep, it should be noted that the VQE
performance of RA is slightly worse than the performance of the LAA, but still
better than that of Plain and LPA.

###  S6 Different rates of random gate activation

The rate to randomly add two-qubit gates each time is $10\%$ in the main text.
In this section, we have also tried different rates for RA training. The rate
$1$ corresponds to the plain training while the rate $1/m$ indicates that we
activate $1/m$ gates each time. The results with fixed PQC depth $l=2$ are
shown in Fig. S8 and the results with fixed PQC depth $l=7$ are shown in Fig.
S9.

When the quantum circuit is shallow ($l=2$) and the barren plateau problems
are not severe, with the decrease in the rate of adding two-qubit gates, the
averaged VQE energy becomes lower and fewer trials are trapped in the local
minima. The RA training can help the optimization trajectories escape from the
local minima. The best VQE energy and median VQE energy is close with
different activation rates because the barren plateau problems are not severe.
When this rate is too large, for example, $1/2$, there is no obvious VQE
performance improvement between the plain training and the RA training.

When the quantum circuit is deep ($l=7$) and the barren plateau problems
become severe, with the decrease in the rate of adding two-qubit gates, the
averaged VQE energy also becomes lower though not as substantial as the
shallow circuit case. The fluctuation from different optimization trials is
still smaller than that of the plain training. However, the better VQE
performance is mainly determined by the mitigation of barren plateaus. The
best VQE energy and the median VQE energy of the RA training is closer to the
ground state energy with the decrease rate of adding two-qubit gates.

Figure S8: The VQE performance with different rates of adding two-qubit gates
each time and fixed PQC depth $l=2$. Figure S9: The VQE performance with
different rates of adding two-qubit gates each time and fixed PQC depth $l=7$.

In practice, the numerical implementation of RA training regards the gate-
adding rate in the sense of statistical average. We use a structure factor to
control the activation of each parameterized two-qubit gate. The initial
structure factors $\vec{g}$ are sampled uniformly from $(0,1)$. If the
structure factor is positive, the two-qubit gate stays as identity gate, i.e.
$\theta=0$ for the unactivated gate, and if the structure factor is negative,
the two-qubit gate is activated whose parameter $\theta$ will update. Assuming
the rate of random gate activation is $1/m$ ($m=10$ in the main text), the
structure factors $\vec{g}\mapsto\vec{g}^{\prime}=\vec{g}-\frac{1}{m}\vec{I}$
after one round of random gate activation. Equivalently, $\frac{1}{m}$
fraction of random two-qubit gates activate on average. It is also worth
noting that the vvag method in TensorCircuit package [69] admits independent
but vectorized optimization loops for batched VQE with different initial
parameters $\vec{\theta}$ and structure factors $\vec{g}$ running
simultaneously, which unlocks infinite possibilities to explore.

###  S7 Comparsion with the initialization strategy

As discussed in the Refs. [15, 16, 17, 18, 19, 20, 21, 22, 23], better
parameters initialization methods can help VQA reach better performance. In
this section, we investigate the plain training where the initial parameters
of $10\%$ two-qubit gates are sampled uniformly from $[0,2\pi]$ and the
remaining initial parameters are set to $0$. We use Plain∗ to represent this
strategy. This ablation study is to demonstrate that the gain of RA approach
is indeed from the randomness of adding two-qubit gates instead of solely from
the better initialization strategy where $90\%$ parameters are set to zero.

The results are shown in Fig. S10. When the circuit is shallow, more trials of
Plain∗ are trapped in the local minima due to the absence of the randomness of
how to add two-qubit gates into the circuit. However, the improvement of the
VQE performance is not as pronounced as the case of Plain training with all
initial parameters sampled uniformly from $[0,2\pi]$. Namely, the gain that RA
training obtains over plain training consists of both contributions from
better initialization strategy and the randomness of gate activation.

When considering real hardware implementation, compared with Plain∗, RA can
reach a slightly better performance with effectively fewer parameterized gates
required, i.e. fewer gradients to be estimated as we discussed in the main
text.

Figure S10: The converged VQE energy ($N=12$, Adam optimizer, decay_steps=100)
for RA training (green) and plain training with better initialization strategy
(blue) with different PQC depths $l$: (a), $J_{z}=0.5$; (b), $J_{z}=1.0$; (c),
$J_{z}=2.0$. The black star represents the averaged VQE energy and the red
line represents the median VQE energy across $500$ independent optimization
trials. The results from RA training are slightly better than Plain∗.

###  S8 Robustness against the quantum noise

To verify the robustness of this new training algorithm against the quantum
noise and the barren plateaus induced by the quantum noise in real hardware,
we carry out VQE optimization in the presence of a quantum depolarizing noise
with noise strength $p=10^{-3}$ after each two-qubit gate. The depolarizing
channel is commonly used in numerical simulation to approximate the noise in a
NISQ device. In these experiments, we also use a total iteration step $1000$
instead of $5000$ above, demonstrating that the approach can be much more
iteration efficient. The results are shown in Fig. S11. The quantum noise
makes the VQE results much worse. However, RA training can still achieve much
better performance than the plain training, in particular with deep PQC in
which the plain training totally fails since the noise accumulates with the
depth. RA training together with early stop technique indeed provides a
feasible and reasonable way to train VQA on noisy devices.

Figure S11: The converged VQE energy ($N=6$, SGD optimizer,
learning$\\_$rate=0.01) for RA training (green) and plain training (blue) with
different PQC depths $l$: (a), $J_{z}=0.5$; (b), $J_{z}=1.0$; (c),
$J_{z}=2.0$. The black star represents the averaged VQE energy and the red
line represents the median VQE energy across $100$ independent optimization
trials. The performance of RA training is much better especially when the
circuit depth is large: the averaged VQE energy and median VQE energy from RA
training is substantially lower than that obtained from plain training.

###  S9 Efficiency comparison for Plain and RA strategies

In this section, we will give a rough estimation of the number of two-qubit
quantum gates required for Plain and RA strategies in total. And this quantity
can be well served as a measure for the total quantum computational resources
consumed.

Assuming $p$ parameters to be optimized in the ansatz (here $p=3Nl$ in our VQE
setup), we have to implement the ansatz with $p$ two qubit gates (we omit the
constant factor brought by the decomposition of
$e^{i\theta\sigma_{i}\sigma_{i+1}}$ gates). To calculate the circuit gradients
via the parameter shift rule [80, 81, 82, 83], it requires to evaluate the
expectation twice for each trainable parameters. And the accuracy of the
gradient for the parameter $\theta$ obtained from $n$ measurement shots is
proportional to $\epsilon\sim\frac{1}{\sqrt{N_{m}}}$ where $N_{m}$ is the
number of measurement shots. Due to the exponentially ($O(e^{-p})$) vanishing
variance induced by barren plateaus, we need $O(e^{2p})$ measurement shots to
achieve an acceptable accuracy. Recall each measurement shot requires to run a
circuit with $p$ quantum gates. Therefore, the number of total quantum gates
required in Plain training is:

$\displaystyle N_{g}=2p^{2}e^{2p},$ (S3)

where the factor $2$ corresponds to the two expectations and we have neglected
the unimportant constants.

For RA training, assuming the rate of adding two-qubit gates is $1/m$ ($m=10$
in the main text), the number of required quantum gates is on average:

$\displaystyle
N^{\prime}_{g}=\frac{\sum_{k=1}^{m}(\frac{k}{m})^{2}2p^{2}e^{\frac{k}{m}2p}}{m}.$
(S4)

In the $m\longrightarrow\infty$ limit, i.e., we only randomly activate one
two-qubit gate each time during the optimization,

$\displaystyle N^{\prime}_{g}$ $\displaystyle=$
$\displaystyle\int_{0}^{1}2p^{2}x^{2}e^{2xp}dx$ $\displaystyle=$
$\displaystyle\frac{-1+(2p^{2}-2p+1)e^{2p}}{2p}.$

Then,

$\displaystyle\frac{N^{\prime}_{g}}{N_{g}}=\frac{1}{2p},$ (S6)

in the leading order. Therefore, RA training can have a $2p$ times
computational resources efficiency improvement compared with the Plain
training (including the Plain∗ training). In the practical experiments we
consider above, $p=252(N=12,l=7)$, there is a roughly $500$ times
computational resources efficiency improvement. And the computational
resources save is more evident when the system gets larger.

###  S10 Entanglement Phase Transition

We focus on the HVA circuit in this work. The circuit consists of blocks of
two-qubit gates. We consider the case when only $0\leq p\leq 1$ ratio of the
two qubit gates are activated (selected randomly) and other gates are just
identities. These activated two-qubit gates can be regarded as drawn from Haar
ensemble or random Clifford ensemble. By measuring the half-chain entanglement
entropy of the final output state, we conclude that there is an entanglement
phase transition on the random circuit setup driven by the filling ratio $p$.

Fig. S12 shows the numerical results on the entanglement phase transition
which is numerically simulated using stabilizer circuit, namely, each
activated two-qubit gate in the original HVA circuit is replaced with one
random two-qubit Clifford gate. And the entanglement of the output state is
further averaged over different choice of activated gate and different
Clifford implementation of each activated gate. We use $8L$ blocks for the HVA
circuit such that the entanglement of the final state is saturated to the
infinite time limit. From the results, we clearly see a entanglement phase
transition induced by the filling ratio of two-qubit gates $p$, which turns
out to be a relevant perturbation. Therefore, we believe that the final state
during the RA training is always in the same phase of the final state of plain
training. This fact contributes to the effectiveness and correctness of the RA
training algorithm.

Figure S12: When $p>p_{c}$, the entanglement entropy obeys “volume law”. The
inset is the finite size data collapse with the scaling form:
$S_{A}(p)-S_{A}(p_{c})=F((p-p_{c})L^{1/\nu})$, where $S_{A}(p_{c})=0$.
$p_{c}=0$, and the critical exponent is fit to $\nu=0.9$.


# Electron spin spectroscopy of single TEMPO dimers coupled via interfering
tunneling currents

Yishay Manassen1,2, Moamen Jbara1,3, Michael Averbukh1, Zion Hazan1, Carsten
Henkel4 and Baruch Horovitz1,2 1Department of Physics, Ben Gurion University
of the Negev, Beer Sheva 84105, Israel
2 The Ilse Katz Institute for Nanoscale Science and Technology, Ben-Gurion
University of the Negev, Beer Sheva 84105, Israel
3 Department of Chemistry, Technion, Haifa, Israel 32000
4 University of Potsdam, Institute of Physics and Astronomy, 14476 Potsdam,
Germany

###### Abstract

We report the detection of electron spin resonance (ESR) in individual dimers
of the stable free radical 2,2,6,6-tetramethyl-piperidine-1-oxyl (TEMPO). ESR
is measured by the current fluctuations in a scanning tunnelling microscope
(ESR-STM method). The multi-peak power spectra, distinct from macroscopic
data, are assigned to dimers having exchange and Dzyaloshinskii-Moriya
interactions in presence of spin-orbit coupling. These interactions are
generated in our model by interfering electronic tunneling pathways from tip
to sample via the dimer’s two molecules. This is the first demonstration that
tunneling via two spins is a valid mechanism of the ESR-STM method.

The attempt to detect and manipulate a single spin in individual molecules is
a fundamental challenge wrachtrup ; koehler ; rugar ; elzerman ; xiao . A
promising tool to monitor the electron spin resonance (ESR) on the nm-scale is
based on a scanning tunneling microscope (STM) that measures current-current
correlations in a static magnetic field (ESR-STM) manassen1 ; manassen2 ;
manassen3 ; durkan1 ; durkan2 ; komeda ; saaino , rather than using external
radiofrequency fields. The experiments so far resulted in a signal at the
Larmor frequency, a signal that is sharp even at room temperature and whose
frequency varies linearly with the applied magnetic field durkan1 ; komeda .
Such a current spectrum has been observed in several spin systems, including
dangling bonds manassen1 ; manassen2 , metal impurities in silicon manassen3
and adsorbed paramagnetic molecular radicals. durkan1 ; durkan2  On the
Si(111)$7\times 7$ surface, two peaks show up and relate to defects that
differ in STM images komeda ; saaino . ESR-STM has been used to detect the
hyperfine spectrum of a single spin in SiC balatsky2 ; manassen5 . Similarity
to macroscopic ESR was demonstrated in the spectrum of silicon vacancy
manassen5 , showing hyperfine contributions from 29Si nuclei.

Recently, a different type of ESR-STM was observed at low temperatures, using
a spin polarized tip and rf irradiation mulleger ; baumann ; willke ; seifert
. Furthermore, single spin ENDOR (electron nuclear double resonance) was
performed manassen6  by applying an rf field at frequencies of the nuclear
transitions and monitoring the intensity of the hyperfine line observed by
ESR-STM; this facilitated measurements of the hyperfine coupling, the
quadrupole coupling and the nuclear g-factors.

Several theoretical models for ESR-STM have been put forward caso ; golub ;
balatsky1 ; manassen4 ; horovitz1 , and it was shown that tunneling via a
single spin cannot explain the observations caso ; golub . Instead there must
be at least two tunneling channels whose interference generates ESR-STM caso ;
golub ; horovitz1 . The two channels are most likely due to two distinct spin
sites: one is the target scanned by the STM probe, while the other is possibly
located on the tip itself. The lack of direct evidence for the second spin has
made elusive the interpretation of the experiments.

In the present work we consider a system of 2,2,6,6-tetramethyl-
piperidine-1-oxyl (TEMPO) molecules that agglomerate and are likely to form
dimers, see Fig. 1. The TEMPO molecule is a stable free radical carrying a
spin 1/2, hence a dimer would be an ideal setup for forming parallel tunneling
routes that lead to the ESR-STM phenomenon. In our experiment, neither tip nor
substrate are spin polarized and rf radiation is not applied. We model the
data with a two-spin model that allows for exchange and Dzyaloshinskii-Moriya
interactions horovitz1 ; horovitz2 . We find reasonable agreement of the
theory with many experimental spectra, thus providing a complete theoretical
interpretation. We note also that the dimer scenario is of much interest to
quantum information science, since by tuning parameters, a long-lived dark
state is available and quantum entanglement can be achieved horovitz2 .

In the ESR-STM experiments we performed, the molecules are deposited on gold
films of thickness 100nm on Mica buchholz . TEMPO is dissolved in toluene and
drop casted on the surface at a concentration of $0.041\,{\rm g}/25\,{\rm
ml}$, corresponding to one monolayer. After drop casting the sample is put in
a UHV chamber (pressure range $10^{-10}\,{\rm torr}$), and tunneling probe
data are taken at room temperature. The clean Au surface shows flat terraces
of variable shape, some of them triangular (Fig. 1a). After deposition, the
molecules disperse as individual entities fairly uniformly on the surface
(Fig. 1b), afterwards there is a slow agglomeration process (several days)
where some parts of the surface end up clean gold (as in Fig. 1a), while in
other parts there is a dense monolayer of TEMPO molecules (Fig. 1c). STM
studies of TEMPO adsorbed on Si(111)$7\times 7$ have shown that the molecule
adsorbs with its NO axis (Fig. 1d) normal to the surface pitters . More
relevant here are studies of TEMPO with Au spheres, showing the disappearance
of the ESR signal when the NO group is close to the Au surface zhang .

Figure 1: STM image of Au(111) surface before deposition of TEMPO. (b) Au(111)
after deposition with single molecules. (c) TEMPO adsorbed on Au(111) when
agglomeration is reached. All images are $30\times 30\,{\rm nm}^{2}$.
Tunneling conditions: (a) $-1\,{\rm V},\,0.4\,{\rm nA}$; (b) $-0.1\,{\rm
V},\,0.5\,{\rm nA}$; (c) $-0.5\,{\rm V},\,0.5\,{\rm nA}$. (d) The molecular
structure of TEMPO. The unpaired electron occupies a $p_{z}$ orbital on the
N$-$O bond kobayashi .

The g-tensor and hyperfine coupling of TEMPO are well documented kobayashi ,
showing an almost isotropic g-factor $g\approx 2.007$ (within .5%). The
eigenvalues of the hyperfine interaction with the nuclear ${}^{14}{\rm N}$
spin are 17, 15, and 94 MHz for the $x$, $y$ and $z$ directions, respectively
(the singly-occupied p-orbital of N defines the $z$-axis). The dominant
hyperfine coupling $a$ in the direction parallel to the magnetic field depends
on the molecular orientation and possibly on averaging due to rotational
motion. Hence $a$ is a fitted parameter that turns out to be $50{-}80\,{\rm
MHz}$ between various molecular sites. The nominal external magnetic field,
perpendicular to the surface, is $230\,{\rm G}$, which with $g=2$ would
correspond to a Larmor frequency of $\nu=644\,{\rm MHz}$, yet we take $\nu$ as
a fitting parameter, allowing for uncertainties in the actual field.
Experimental spectra of the tunnel current are taken in 9800 channels covering
the range of 580-780 MHz, each channel requiring about $50\,\mu{\rm s}$
acquisition time. This is short compared to the spin lifetime of the 14N
nucleus, typically in the range of $0.5{-}1\,{\rm ms}$ [Ref. odell, ], so that
we may assume fixed spin states. The spectrum analyzer averages 200 spectra,
however, so that eventually nuclear spin flips occur, and an ensemble of all
levels is probed. A whole spectrum of a single site takes 90 s.

We first outline our model and results. Each TEMPO molecule has three
${}^{14}{\rm N}$ nuclear spin states whose hyperfine coupling to the unpaired
electron spin splits the ESR into $\nu,\,\nu\pm a$. A dimer with two spins has
then nine hyperfine states with nuclear spin projections $m,\,m^{\prime}=0,\pm
1$. When electrons transit between the electrodes (tip and substrate) through
two TEMPO molecules, interfering exchange tunneling events occur and generate
interactions between the two molecular spins as well as dissipation (linewidth
$\Gamma$) horovitz1 ; horovitz2 . Both interactions and dissipation depend
sensitively on the energy levels of the two molecules being of two types,
either degenerate ($m=m^{\prime}$) or non-degenerate ($m\neq m^{\prime}$). We
assume that $a\gg\Gamma$ (confirmed by our analysis) so that the distinction
between the two types is well defined. There are three dimer states with
$m=m^{\prime}$ that are degenerate and therefore have more (secular) terms to
be kept when deriving the master equation; these generate exchange as well as
Dzyaloshinskii-Moriya interactions horovitz2 . In addition there are six non-
degenerate pairs $m\neq m^{\prime}$ with a weaker exchange interaction that
splits the three ESR transitions into six peaks. Although the degenerate pairs
give more peaks, their statistical weight is smaller, and only the two extreme
ones are usually visible, leading to overall eight peaks. It is interesting to
note that the behaviour in the more prominent dips between the peaks is most
informative: their positions determine the bare hyperfine parameter, and their
shapes are sensitive to couplings among degenerate pairs.

We proceed to describe our theoretical model for TEMPO dimers in more detail.
Compared to its component $a$ parallel to the magnetic field, the transverse
hyperfine splitting $b$ can be neglected because it affects the spectra only
in second order $\sim(b/\nu)^{2}$. Such terms may cause a small difference
between the two hyperfine splittings and even shift $\nu$. In addition to
Larmor and hyperfine terms, the two-spin Hamiltonian contains effective
interactions due to tunneling to either tip or substrate. They were derived in
Refs. horovitz1, ; horovitz2,  within a Born-Markov master equation and turn
out fairly large, as they depend on the electron bandwidth $\Delta$ of the
electrodes. The degenerate cases (nuclear levels $m=m^{\prime}$) involve
exchange coupling ($J_{\rm ex}$) and Dzyaloshinskii-Moriya coupling ($J_{\rm
DM}$). The latter appears only when spin-orbit coupling in the tunneling
junction is taken into account which is actually essential to the presence of
ESR-STM horovitz1 . It is parametrized by two angles $\theta,\phi$ of an SU(2)
spin rotation matrix acting on the tunneling electron spin. The interaction
strengths are

$\displaystyle\left.\begin{array}[]{r}J_{\rm ex}\\\ J_{\rm
DM}\end{array}\right\\}$ $\displaystyle=$ $\displaystyle
4J_{1}J_{2}N^{2}(0)\Delta\cos\mbox{\small$\frac{1}{2}$}\theta\left\\{\begin{array}[]{r}\cos\phi\\\
\sin\phi\end{array}\right.$ (5)

where $J_{1},\,J_{2}$ are the exchange tunneling elements via the two
localised spins, respectively, and $N(0)$ is the density of states of either
electrode at its Fermi level. The Hamiltonian of a degenerate dimer, for
$m=-1,\,0,\,1$, is

$\displaystyle{\cal H}_{\rm deg}$ $\displaystyle=$
$\displaystyle\mbox{\small$\frac{1}{2}$}(\nu+am)[\tau_{z}\otimes\mathbbm{1}+\mathbbm{1}\otimes\tau_{z}]-J_{\rm
ex}\bm{\tau}\otimes\bm{\tau}$ (6) $\displaystyle{}+J_{\rm
DM}[\tau_{x}\otimes\tau_{y}-\tau_{y}\otimes\tau_{x}]$

where the tensor products describe the two localized spins. The exchange
tunneling between the tip and substrate, including the spin-orbit coupling,
also leads to a finite linewidth for each spin, $\Gamma_{1},\,\Gamma_{2}$. We
note that in this model the dimer molecules need not to be close to each
other, since the interactions are generated via tunneling to tip or substrate.
Indeed, the separation of molecules in Fig. 1b is 2-4 nm, a scale on which the
direct dipole-dipole interaction can be neglected, justifying our model
Hamiltonian.

The non-degenerate dimer states are decoupled from the degenerate ones in the
master equation by the secular approximation, valid in our experiment where
$a\gg\Gamma_{1,2}$. Each pair $m\neq m^{\prime}$ has an anisotropic exchange
coupling (being the only secular term horovitz2 ), and the Hamiltonian is

$\displaystyle{\cal H}_{\rm nondeg}$ $\displaystyle=$
$\displaystyle\mbox{\small$\frac{1}{2}$}(\nu+am)\tau_{z}\otimes\mathbbm{1}+\mbox{\small$\frac{1}{2}$}(\nu+am^{\prime})\mathbbm{1}\otimes\tau_{z}$
(7) $\displaystyle{}-J_{\rm ex}\tau_{z}\otimes\tau_{z}$

It is useful to list the eigenstates, levels and transition frequencies, see
Table I. We note that the hyperfine transitions $\nu+am,\,m=-1,\,0,\,1$ are
split by the interactions. In particular $J_{\rm DM}$ is responsible for the
splitting of the $T1\rightarrow T2,\,T2\rightarrow T3$ transitions.

Table 1: Spectra and transitions of degenerate dimers (first 4 lines where $\tan\psi=J_{\rm DM}/J_{\rm ex}$, $J=\sqrt{J_{\rm ex}^{2}+J_{\rm DM}^{2}}$ and nuclear levels $m=m^{\prime}=-1,\,0,\,1$) and non-degenerate dimers (last line where $s,s^{\prime}=\pm 1$ are electron spin states). The $S,\,T$ notation refers to singlet or triplet states that decouple exactly in the limit $J_{\rm DM}=0$. Transitions are given only near the Larmor frequency $\nu$.  eigenstates | energy levels | transition frequencies
---|---|---
$S:|{\uparrow\downarrow}\,mm\rangle-\mbox{e}^{i\psi}|{\downarrow\uparrow}\,mm\rangle$ | $2J+J_{\rm ex}$ | $S\rightarrow T1,\,T3:$
$T1:|{\uparrow\uparrow}\,mm\rangle$ | $\nu-J_{\rm ex}+am$ | $(\nu+am)\pm 2(J+J_{\rm ex})$
$T2:|{\uparrow\downarrow}\,mm\rangle+\mbox{e}^{i\psi}|{\downarrow\uparrow}\,mm\rangle$ | $-2J+J_{\rm ex}$ | $T1\rightarrow T2,\,T2\rightarrow T3:$
$T3:|{\downarrow\downarrow}\,mm\rangle$ | $-\nu-J_{\rm ex}-am$ | $\nu+am\pm 2(J-J_{\rm ex})$
$|ss^{\prime}m\neq m^{\prime}\rangle$ | $\mbox{\small$\frac{1}{2}$}s(\nu+am)+\mbox{\small$\frac{1}{2}$}s^{\prime}(\nu+am^{\prime})-ss^{\prime}J_{\rm ex}$ | $\nu+am\pm 2J_{\rm ex}$

The fluctuations in the tunneling current arise from spin flips. We model
these by time correlations of two spin operators, one carries a current with
spin flip $\tau_{\pm}$, the other carrying a current without spin flip
$\tau_{z}$. (There are additional combinations from flips of both spins
horovitz2 , but outside our frequency range and not considered here.) Hence
the current correlations are proportional to

$\displaystyle C_{2}(\nu_{1},\nu_{2},\omega)$ $\displaystyle=$
$\displaystyle\langle(\tau_{-}\otimes\tau_{z})_{t}(\tau_{+}\otimes\tau_{z})_{0}\rangle_{\omega}$
$\displaystyle{}+\langle(\tau_{z}\otimes\tau_{-})_{t}(\tau_{z}\otimes\tau_{+})_{0}\rangle_{\omega}+(+\leftrightarrow-)$

where $\langle\ldots\rangle_{\omega}$ is a Fourier transform. The frequencies
$\nu_{1,2}=\nu-a,\nu,\nu+a$ arise from the nuclear quantum numbers
$m,m^{\prime}$ near the first and second spin. Note that the interchange
$\nu_{1},\nu_{2}$ yields distinct results, e.g. if $\Gamma_{1}\neq\Gamma_{2}$.
For all nuclear configurations of a dimer, the non-degenerate (six pairs) and
degenerate (three pairs) contributions are

$\displaystyle C_{\rm nondeg}(\omega)$ $\displaystyle=$
$\displaystyle\sum_{\pm}\big{[}C_{2}(\nu\pm a,\nu,\omega)+C_{2}(\nu\pm
a,\nu\mp a,\omega)$ $\displaystyle\qquad{}+C_{2}(\nu,\nu\pm a,\omega)\big{]}$
$\displaystyle C_{\rm deg}(\omega)$ $\displaystyle=$ $\displaystyle
C_{2}(\nu-a,\nu-a,\omega)+C_{2}(\nu,\nu,\omega)$ (9)
$\displaystyle\qquad{}+C_{2}(\nu+a,\nu+a,\omega)$

The total observable is $C_{\rm nondeg}(\omega)+C_{\rm deg}(\omega)$.

We evaluate the spin-spin correlations via the master equation for both the
degenerate and non-degenerate cases and fit parameters to the experimental
data. We consider first the general features of the theoretical results, an
example is in Fig. 2(left). The dominant terms are the six non-degenerate
terms (higher red lines) that are split by $\pm 2J_{\rm ex}$. Between these
six peaks are strong dips at the bare (non-interacting) hyperfine transitions
$\nu,\,\nu\pm a$, thus readily determining these parameters. The degenerate
terms (lower black lines) split from the bare transitions by both a strong and
a weak splitting $\pm 2(J\pm J_{\rm ex})$ (see Table I). The weak splitting
$2(J-J_{\rm ex})>0$ arises from $J_{\rm DM}$, if the latter vanishes, some
lines do not split and there would be a (small) peak in the bare locations.
Since this was not seen in any of our data we conclude that $J_{\rm DM}$ is
significant, i.e. at least of order $J_{\rm ex}$. The weak splittings produce
side shoulders on the main peaks, seen in some of our data. The strong
splitting is larger than that of the non-degenerate terms ($2(J+J_{\rm
ex})>2J_{\rm ex}$), hence, although the overall intensity of the degenerate
terms is weaker, the split peaks at the edges of the spectrum become dominant
and lead in general to an apparent 8-peak structure.

Before fitting our data, we note that macroscopic ESR for dimers with exchange
is totally different from the above. In the macroscopic case a homogeneous rf
field is applied so that the measured correlations are symmetric in both
spins,
$\langle(\tau_{-}\otimes\tau_{z}+\tau_{z}\otimes\tau_{-})_{t}(\tau_{+}\otimes\tau_{z}+\tau_{z}\otimes\tau_{+})_{0}\rangle_{\omega}$.
The permutation symmetry forbids in particular $S\rightarrow T$ transitions
for $m=m^{\prime}$ and $J_{\rm DM}=0$. Furthermore, in the STM case the
exchange is generated via the electrodes and is different for degenerate and
non-degenerate pairs [compare Eqs. (6, 7)], while in the macroscopic case,
there is a direct exchange that applies equally to all $m,\,m^{\prime}$.
Finally, the spin-orbit interaction responsible for the coupling $J_{\rm DM}$
is most likely generated by the heavy metal atoms in the tip, hence absent in
a macroscopic setup.

The spectra are limited by the detection sensitivity of the impedance matching
circuits to the frequencies $580{-}780\,{\rm MHz}$, thus some of the expected
peaks are probably outside this range. We have included in our analysis only
data that have at least seven peaks (except for Fig. 3d, see below); we have
found 15 such spectra. Fig. 2(right) shows the theory spectrum (black smooth
line), the sum of the two curves in Fig. 2(left), compared to experimental
data, the magenta jagged curve. The parameters of the fit are given in the
caption, in addition we use the linewidths $\Gamma_{1},\Gamma_{2}\approx
4\,{\rm MHz}$. (We note that in the absence of splitting due to $J_{\rm ex}$
and $J_{\rm DM}$, the linewidth becomes $\Gamma_{1}+\Gamma_{2}$.)

Figure 2: Left: Spectra of nondegenerate correlations (upper red line) and of
degenerate ones (lower black line) vs. $f=\omega/2\pi$. Right: Experimental
(magenta) and theory (black) spectra. Fitting parameters in MHz: $\nu:665$,
$a:50,55$, $J_{\rm ex}:8$, $J_{\rm DM}:9$; hyperfine $a$ is somewhat different
between left and right.

Fig. 3(a-c) shows additional fitting curves with fairly similar parameters,
obtained for different spots on the sample. In all cases, reasonable agreement
is found with realistic fit parameters, note in particular that the couplings
$J_{\rm ex}$, $J_{\rm DM}$ are comparable in magnitude, which points, on the
basis of the model, towards a sizable rotation angle $\tan\phi\sim 1$ in the
spin-orbit interaction. The fitted values of $a$ are in between the hyperfine
eigenvalues for the $z$ and $x$, $y$ axes kobayashi . This suggests some tilt
or rotational averaging of the molecular $z$ axis relative to the magnetic
field (normal to the surface).

An exception is Fig. 3d that appears to have only 3 peaks, as if $J_{\rm
ex}=J_{\rm DM}=0$. Fitting the data this way yields line widths
$\Gamma_{1},\,\Gamma_{2}$ that are much larger than all other cases, however.
We believe that it is more likely that a weak, but finite exchange applies in
this case (caption of Fig. 3d), causing shoulders and an apparent increase in
width of these lines.

In conclusion, we have found a large set of ESR-STM spectra that that fits
well to a theory of two spins located on a molecular dimer and coupled via
electrons that tunnel between tip and substrate. The fitting parameters give
effective exchange and Dzyaloshinskii-Moriya couplings that are comparable and
fairly strong (comparable to the hyperfine splitting, using the conventional
definitions $4J_{\rm ex},\,4J_{\rm DM}$ of Ref. wertz, ). Our analysis of
these dimers opens a route for studying hyperfine interactions and g-factors
in molecules and determining their parameters. It also paves the road to
measure spin-orbit coupling for tunneling electrons.

Figure 3: Experimental ESR-STM spectra (jagged magenta lines) and fitted
theoretical curves using the following parameters (in MHz): (a) $\nu:666$,
$a:71,\,69$, $J_{\rm ex}:7$, $J_{\rm DM}:12$. (b) $\nu:683$, $a:67$, $J_{\rm
ex}:7$, $J_{\rm DM}:12$. (c) $\nu:663$, $a:60,\,59$, $J_{\rm ex}:7$, $J_{\rm
DM}:10$. (d) $\nu:677$, $a:80,\,73$, $J_{\rm ex}:1$, $J_{\rm DM}:8$. The small
differences in $a$ correspond to different nuclear $m$ levels and may arise
from second-order shifts due to sub-leading hyperfine couplings.

###### Acknowledgements.

This work was funded by Attract – Developing breakthrough technologies for
science and society, grant NMR(1). Additional funding was provided by the ISF
collaboration with China, grant 2811/19, and the German DFG project number
282981637. C.H. and B.H. gratefully acknowledge funding by the German DFG
through the DIP programme (Fo 703/2-1).

## References

  * (1) J. Wrachtrup, C. von Borczyskowski, J. Bernard, M. Orrit, and R. Brown, Optical detection of magnetic resonance in a single molecule, Nature 363, 244 (1993).; Optically detected spin coherence of single molecule, Phys. Rev. Lett. 71, 3565 (1993).
  * (2) J. Koehler, J. A. J. M. Disselhorst, M. C. J. M. Donckers, E. J. J. Groenen, J. Schmidt, and W. E. Moerner, Magnetic resonance of a single molecular spin, Nature 363, 242 (1993).
  * (3) D. Rugar, R. Budakian, H. J. Mamin, and B. W. Chui, Single spin detection by magnetic resonance force microscope, Nature 430, 329 (2004).
  * (4) J. M. Elzerman, P. Hanson, L. H. Willems van Beveren, B. Witkamp, L. M. K. Vandersypen, and L. P. Kouwenhoven, Single shot read-out of an individual electron spin in a quantum dot, Nature 430, 431 (2004).
  * (5) M. Xiao, I. Martin, E. Yablonovitch, and H. W. Jiang, Electrical detection of the spin resonance of a single electron in a silicon-field effect transistor, Nature 430, 435 (2004).
  * (6) Y. Manassen, E. Ter-Ovanesyan, D. Shachal, and S. Richter, Electron-Spin Resonance-Scanning Tunneling Microscopy Experiments on Thermally Oxidized Si(111), Phys. Rev. B 48, 4887 (1993).
  * (7) Y. Manassen, Real-time response and phase-sensitive detection to demonstrate the validity of ESR-STM results, J. Magn. Reson. 126, 133 (1997).
  * (8) Y. Manassen, I. Mukhopadhyay, and N. Ramesh Rao, Electron-spin-resonance STM on iron atoms in silicon, Phys. Rev. B 61, 16223 (2000).
  * (9) C. Durkan and M. E. Welland, Electronic spin detection in molecules using scanning-tunneling microscopy assisted electron-spin resonance, Appl. Phys. Lett. 80, 458 (2002).
  * (10) C. Durkan, Detection of single electronic spins by scanning tunneling microscopy, Contemp. Phys. 45, 1 (2004).
  * (11) T. Komeda and Y. Manassen, Distribution of frequencies of a single precessing spin detected by scanning tunneling microscope, Appl. Phys. Lett. 92, 212506 (2008).
  * (12) Y. Saaino, H. Isshiki, S.M.F. Shahed, T. Takaoka, and T. Komeda, Atomically resolved Larmor frequency detection on Si(111)-$7\times 7$ oxide surface, Appl. Phys. Lett. 95, 082504 (2009).
  * (13) A V. Balatsky, M. Nishijima and Y. Manassen, Electron spin resonance scanning tunneling microscopy, Adv. Phys. 61, 117 (2012).
  * (14) Y. Manassen, M. Averbukh and M. Morgenstern, Analyzing multiple encounter as a possible origin of electron spin resonance signals in scanning tunneling microscopy on Si(111) featuring C and O defects, Surf. Sci. 623, 47 (2014).
  * (15) S. Müllegger, S. Tebi, A. K. Das, W. Schöfberger, F. Faschinger, and R. Koch, Radio Frequency Scanning Tunneling Spectroscopy for Single-Molecule Spin Resonance, Phys. Rev. Lett. 113, 133001 (2014).
  * (16) S. Baumann, W. Paul, T. Choi, C. P. Lutz, A. Ardavan, and A. J. Heinrich, Electron paramagnetic resonance of individual atoms on a surface, Science 350, 417 (2015).
  * (17) P. Willke et al., Probing quantum coherence in single-atom electron spin resonance, Science Adv. 4, eaaq1543 (2018).
  * (18) T. S. Seifert, S. Kovarik, C. Nistor, L. Persichetti, S. Stepanow, P. Gambardella, Single-atom electron paramagnetic resonance in a scanning tunneling microscope driven by a radiofrequency antenna at 4K, Phys. Rev. Research 2, 013032 (2020)
  * (19) Y. Manassen, M. Averbukh, M. Jbara, B. Siebenhofer, A. Shnirman, B. Horovitz, Fingerprints of single nuclear spin energy levels using STM-ENDOR, J. Magn. Reson. 289, 107 (2018).
  * (20) A. Caso, B. Horovitz and L. Arrachea, Model for electron spin resonance in STM noise, Phys. Rev. B 89, 075412 (2014).
  * (21) A. Golub and B. Horovitz, nanoscopic interferometer model for spin resonance in current noise, Phys. Rev. B 88, 115423 (2013).
  * (22) A. V. Balatsky, Y. Manassen, and R. Salem, Exchange-based noise spectroscopy of a single precessing spin with scanning tunnelling microscopy, Phil. Mag. B 82, 1291 (2002); ESR-STM of a single precessing spin: Detection of exchange-based spin noise, Phys. Rev. B 66, 195416 (2002).
  * (23) Y. Manassen and A. V. Balatsky, $1/f$ Spin Noise and a Single Spin Detection with STM, Israel J. Chem. 44, 401 (2004), special issue on single molecule spectroscopy, arXiv: cond-mat/0402460.
  * (24) B. Horovitz and A. Golub, Double Quantum Dot scenario for spin resonance in current noise, Phys. Rev. B 99, 241407(R) (2019).
  * (25) B. Horovitz and C. Henkel, Spin entanglement via scanning tunneling microscope current, Phys. Rev. B 104, L081405 (2021)
  * (26) S. Buchholz, H. Fuchs, and J. P. Rabe, Surface structure of thin metallic films on mica as seen by scanning tunneling microscopy, scanning electron microscopy, and low energy electron diffraction, J. Vac. Sci. Technol. B 9, 857 (1991).
  * (27) J. L. Pitters, P. G. Piva, X. Tong, and R. A. Wolkow, Reversible Passivation of Silicon Dangling Bonds with the Stable Radical TEMPO, Nano Lett. 3, 1431 (2003).
  * (28) Z. Zhang, A. Berg, H. Levanon, R. W Fessenden and D. Meisel, On the interactions of free radicals with gold nanoparticles, J. Am. Chem. Soc. 125, 7960 (2003).
  * (29) H. Kobayashi, T. Ueda, K. Miyakubo, T. Eguchi and A. Tani, ESR study of molecular dynamics and orientation of TEMPO included in organic 1-D nanochannel, Phys. Chem. Chem. Phys. 10, 1263 (2008).
  * (30) Luke A. O’Dell and Robert W. Schurko, Static solid-state 14N NMR and computational studies of nitrogen EFG tensors in some crystalline amino acids, Phys. Chem. Chem. Phys. 11, 7069 (2009).
  * (31) J. E. Wertz and J. R. Bolton, Electron Spin Resonance, Elementary Theory and Practical Applications (McGraw-Hill, New York, 1986), Chap. 10, Figs. 10-15b.