# Hybrid Exchange Measurement-Based Qubit Operations in Semiconductor Double
Quantum Dot Qubits

Matthew Brooks<EMAIL_ADDRESS>Charles Tahan<EMAIL_ADDRESS>Laboratory for Physical Sciences, 8050 Greenmead Dr., College Park, MD 20740,
USA

###### Abstract

Measurement-based quantum computing (MBQC) promises natural compatibility with
quantum error correcting codes at the cost of a polynomial increase in
physical qubits. MBQC proposals have largely focused on photonic systems,
where 2-qubit gates are difficult. Semiconductor spin qubits in quantum dots,
on the other hand, offer fast 2-qubit gates via the exchange interaction. In
exchange-based quantum computing, as with other solid-state qubits, leakage to
higher states is a serious problem that must be mitigated. Here, two hybrid
measurement-exchange schemes are proposed which quantify the benefits of MBQC
on quantum dot-based quantum computing. Measurement of double quantum dot
encoded qubits in the singlet-triplet basis, along with inter- and intra-qubit
exchange interaction, are used to perform one and two qubit operations. Both
schemes suppress individual qubit spin-state leakage errors, offer fast gate
times, and require only controllable exchange couplings, up to known phase and
Pauli corrections.

##  I Introduction

Measurement-based quantum computing (MBQC) performs unitary rotations of
encoded quantum states in a quantum processor by taking measurements,
projective or partial, of qubits within entangled systemsRaussendorf and
Briegel (2000); Verstraete and Cirac (2004); Raussendorf  _et al._ (2002);
Nielsen (2001); Raussendorf (2009); Jozsa (2006); Browne and Briegel (2019).
The amount of calculations that can be performed on an MBQC processor is
limited by the amount of entangled qubits in the initial state known as the
resource state. This allows for parallelization of some gates on such a
processor by simultaneously performing commuting measurements of qubits within
the resource stateRaussendorf and Briegel (2000); Verstraete and Cirac (2004);
Raussendorf  _et al._ (2002); Nielsen (2001); Raussendorf (2009); Jozsa
(2006); Browne and Briegel (2019); Raussendorf and Harrington (2007); Van den
Nest  _et al._ (2006); Plato  _et al._ (2008); Shettell and Markham (2020).
This is unlike more conventional quantum circuit-based models realized in
solid state systems where unitaries are achieved by driving and interacting
the qubitsBarends  _et al._ (2014); Arute  _et al._ (2019); Preskill (2018);
Wang  _et al._ (2018). In such systems like superconducting circuitsArute  _et
al._ (2019); Wang  _et al._ (2018), trapped ionsZhang  _et al._ (2017);
Ballance  _et al._ (2016) or semiconductor quantum dots (QDs)Hanson  _et al._
(2007); Russ and Burkard (2017); Mi  _et al._ (2018); Veldhorst  _et al._
(2017), qubit control typically stems from interactions with externally-driven
electric, magnetic or laser fields.

Semiconductor spin qubits offer long coherence times and precise control of
inter-QD energy detuning and tunnel couplingsTakeda  _et al._ (2020); Yang
_et al._ (2019); Zajac  _et al._ (2018); Russ  _et al._ (2018) within both
commercially available and isotopically purifiable semiconductors, making them
excellent candidates for quantum processors. Exchange interactions allow for
fast spin-spin entangling operations susceptible to errors due to charge noise
in the controlling gatesRuss and Burkard (2017); Reed  _et al._ (2016); Yoneda
_et al._ (2018). Encoded qubits in the larger Hilbert space of two or more
spinsDiVincenzo  _et al._ (2000); Bacon  _et al._ (2000); Levy (2002); Johnson
_et al._ (2005); Laird  _et al._ (2010); Russ and Burkard (2017) are employed
to allow operations-via-exchange, at the cost of more qubits and increased
susceptibility to leakage errors. Leakage errors occur when the qubit couples
to spin states outside the encoded qubit subspace and are generally caused by
inter-qubit exchange interactionWardrop and Doherty (2014) and charge noise or
hyperfine interactionsAndrews  _et al._ (2019), depending on the encoded
subspace. In this work we are interested in exploring MBQC with encoded spin
qubits to avoid these issues by demonstrating full qubit control, within
approaches that suppress qubit leakage and are inherently compatible with
quantum error correction at the cost of a polynomial increase in the required
qubitsRaussendorf (2009); Raussendorf  _et al._ (2006).

Proof-of-concept experimental realizations of MBQC exist in photonic
systemsBell  _et al._ (2014a); Russo  _et al._ (2018); Bodiya and Duan (2006),
and also in trapped ion systemsLanyon  _et al._ (2013). The challenges of
photonic-based QC (photons do not naturally interact) traditionally motivated
MBQC, while the assumed large physical qubit resource overhead has discouraged
exploration of solid-state MBQC architectures. A notable exception in solid-
state implementations of quantum processors are charge parity measurement-
based anyonic braiding operations in Majorana fermion qubits coupled to
QDsKarzig  _et al._ (2017); Xue (2013); Bonderson and Lutchyn (2011); Grimsmo
and Smith (2019). This scheme however is functionally equivalent to a typical
gate-based approach as there is no larger resource state employed, thus there
is less of a qubit overhead problem bar the ancila QDs coupling the logical
Majorana qubits.

Here we discuss two protocols for performing MBQC with semiconductor spin
qubits assumed to be given by Si double-quantum-dot (DQD) qubitsZwanenburg
_et al._ (2013); Zajac  _et al._ (2018). In the first protocol, 1 and 2D
cluster states in spin qubits are realized. This scheme offers efficient
arbitrary single qubit state preparation, as well as full single and two qubit
control and compatibility with stabilizer codes by design of static qubit
graph geometries. The second protocol employs entangling measurements to
achieve a hybrid scheme more akin to conventional gate-based quantum
computation, similar to the Majorana-qubit measurement gates schemes. Both
schemes are robust against spin state leakage and offer fast gate times.

This paper is organized as follows, firstly the ingredients of cluster state
MBQC exploited by one of our proposals are introduced in Sec. II. Then, our
first hybrid exchange-MBQC control scheme is given in Sec. III, followed by
the second protocol employing entangling measurements in Sec. IV. In Sec. V
the methods of qubit measurement are discussed, then lastly the results are
summarized in Sec.VI.

##  II Cluster States

The foundation of MBQC are quantum processors made up of large entangled
states, such as a cluster stateRaussendorf (2009); Raussendorf  _et al._
(2006). A cluster state $\ket{\Psi}$ is given by

$\ket{\Psi}=\prod_{j\in N(i)}CZ\ket{+_{i}}\otimes\ket{+_{j}}$ (1)

where $N(i)$ gives all adjacent interacting qubits to qubit $i$, $CZ$ gives
the entangling two qubit control-phase ($CZ$) gate diag$\\{1,1,1,-1\\}$ and
$\ket{\pm}=(\ket{0}\pm\ket{1})/\sqrt{2}$ are the eigenstates of the Pauli-$x$
($\sigma_{x}$) matrix. In the simplest geometry, a 1D array of qubits each
entangled by nearest neighbor, the first of which is in some prepared state
$\ket{\psi}$ (see Fig. 1(a)). Single-qubit rotations on $\ket{\psi}$ are
achieved by measuring each qubit along the array sequentially in the
$\ket{M_{\pm}(\theta)}=(\ket{0}\pm e^{-i\theta}\ket{1})/\sqrt{2}$ basis, each
performing the following rotation

$\ket{\tilde{\psi}}=\sigma_{x}^{s}HP(\theta)\ket{\psi}=\sigma_{x}^{s}W(\theta)\ket{\psi}$
(2)

onto its neighboring qubit, where $s$ is the measurement outcome ($s=0,1$ for
$\ket{M_{\pm}(\theta)}$), $H$ is the Hadamard gate and $P(\theta)$ is a phase
gate, the phase of which depends on the angle of the rotated measurement
basis. The $\sigma_{x}^{s}$ term is a Pauli correction dependent on the
measurement outcome. To perform a controlled arbitrary single qubit rotation,
a maximum of four sequential measurements accounting for the previous Pauli
correction are needed since $H$ and $\sigma_{x}$ do not commute.
Experimentally, performing an $\ket{M_{\pm}(\theta)}$ measurement is
equivalent to a Pauli-$x$ measurement ($\ket{M_{\pm}(0)}$) between phase
gates, i.e. $P(\theta)\ket{M_{\pm}(0)}\bra{M_{\pm}(0)}P(-\theta)$.

Using the same tools to achieve single qubit control, two qubit rotations can
be done in 2D arrays of cluster states. For example, a control-not ($CX$) gate
requires only a resource state of four qubits arranged by three in a line
containing the target qubit and the fourth entangled to the middle qubit as
the control qubit. By two simultaneous $\ket{M_{\pm}(0)}$ measurements on the
wire, the target is teleported along the wire whilst the control qubit is
stationary. Therefore if the necessary ingredients for full single qubit
control can be given for a cluster state scheme, two qubit control can be
assumed to be achievable if a 2D graph can be initialized.

Figure 1: (a) Grid cluster state, consisting of qubits (black dots) all
initialized in the $\ket{+}$ state, entangled to their neighboring qubits by a
$CZ$ drawn as a black line. Two logical qubit states $\ket{\psi_{1}}$ (red)
and $\ket{\psi_{2}}$ (blue) start on the left hand side of the graph, and
follow a path of measurements through the graph given by their respective
color. The qubits are then entangled by a two qubit gate measurement sequence
outputting the state $\ket{\tilde{\Psi}}$, and all unused qubits in the graph
can be disentangled from the output state by Pauli-$z$ measurements. (b)
Linear two semiconductor double quantum dot qubits diagram. Each qubit
consists of two dots (red) is charged with one electron (blue) in a $s-t_{0}$
charge configuration (blue dotted line). The two adjacent dots from each qubit
are coupled by en exchange interaction $J_{23}$ (magenta dotted line) given by
the depicted potential shape (magenta line). The magnetic field gradients for
each qubit as well as the whole system are given (black).

##  III Hybrid Exchange-Measurement Graph Protocol

###  III.1 Qubit Model: Magnetic gradients

In the first protocol, the qubits considered are assumed to be Si DQD spin
qubits charged with two electrons, one in each dot in a charge-0 singlet-
triplet
$\ket{s}(\ket{t_{0}})=(\ket{\uparrow\downarrow}\pm\ket{\downarrow\uparrow})/\sqrt{2}$
configuration isolated by an external magnetic $B_{0}$ field operating at its
symmetric operation point (SOP)Reed  _et al._ (2016). We assume that in the
idle state there is no exchange interaction between the two dots.
Additionally, the qubits are assumed to experience static magnetic field
gradients, implemented by micromagnetsMcNeil  _et al._ (2010), or engineered
$g$-factorsRuskov  _et al._ (2018); Harvey-Collard  _et al._ (2019), across
the two dots. These allow for exchange-only $CZ$ gates to be implemented
between two tunnel-coupled qubits in the limit of weak tunneling $t_{ij}\ll
E_{c}$ (where $t_{ij}$ is the tunnel coupling between dots $i$ and $j$,
$E_{c}$ is the charging energy of the dots) and the exchange-interaction
$J_{ij}$ between the dots satisfies $J_{ij}\ll B_{0}$. In this limit the
exchange interaction between dots is given by $J_{ij}=4t_{ij}^{2}/E_{c}$. This
is described for two qubits coupled linearly by the following exchange
effective Hamiltonian (assuming $\hbar=1$ throughout)Wardrop and Doherty
(2014)

$H_{\text{eff}}=(\mu\Delta_{12}+\tilde{B})\sigma_{z}^{1}+(\mu\Delta_{34}+\tilde{B})\sigma_{z}^{2}+\frac{1}{4}J_{23}\left(\sigma_{z}^{1}\sigma_{z}^{2}-\mathbb{I}\right)$
(3)

where $\mu$ is the Bohr magneton, $\Delta_{ij}=B_{i}-B_{j}$ are the intra-
qubit magnetic field gradients where $B_{i}$ is the magnetic field experienced
by the $i^{\text{th}}$ dot, and $\tilde{B}$ is given as

$\begin{split}\tilde{B}=\frac{1}{2}\left[-\frac{1}{2}\left(\mu\Delta_{12}+\mu\Delta_{34}\right)+\mu\Delta
B-\right.\\\
\left.-\sqrt{\left[\mu\Delta_{B}-\frac{1}{2}\left(\mu\Delta_{12}+\mu\Delta_{34}\right)\right]^{2}+\frac{1}{4}J_{23}^{2}}\right]\end{split}$
(4)

where $\Delta B$ is the inter-qubit magnetic field difference, i.e. $\Delta
B=(B_{1}+B_{2}-B_{3}-B_{4})/2$, and $\sigma_{z}^{i}$ is the Pauli-$z$ matrix
of the $i^{\text{th}}$ qubit in the $x$-basis of a singlet-triplet qubit, i.e.
$\sigma_{z}=\ket{1}\bra{1}-\ket{0}\bra{0}$ where
$\ket{0}=\left(\ket{t_{0}}+\ket{s}\right)/\sqrt{2}=\ket{\uparrow\downarrow}$
and
$\ket{1}=\left(\ket{t_{0}}-\ket{s}\right)/\sqrt{2}=\ket{\downarrow\uparrow}$.
Assuming $\Delta_{B}=\Delta_{12}=\Delta_{34}=\Delta$, a $CZ$ is achieved
between the two qubits at times

$\tau=\frac{n_{1}\pi}{\mu\Delta}$ (5)

when the inter-qubit exchange satisfies

$J_{23}=\frac{\mu\Delta(4n_{1}-2n_{2}-1)}{n_{1}}$ (6)

where $n_{1}\in\mathbb{Z},\neq 0$ and $n_{2}\in\mathbb{Z}$. Leakage of the
qubits into spin states outside of the computational basis (i.e.
$\ket{\uparrow\uparrow}$ and $\ket{\downarrow\downarrow}$) is suppressed if
$J_{23}\ll\mu\Delta$Wardrop and Doherty (2014), which is assumed here. A two
qubit DQD set up like this is depicted in Fig. 1(b). Additionally it is
assumed that the qubits in both proposed protocols are operated at their SOP.
The SOP is the point at which $\partial E_{s/t}/\partial\epsilon=0$, where
$E_{s/t}$ are the eigenvalues of the singlet/triplet spin-states and
$\epsilon$ is the detuning between the two dots of the qubit. Operating here
protects against first-order charge noise.

Armed with a $CZ$ gate, the necessary ingredients for initializing a cluster
state are given. Firstly, the qubits are initialized into a spin-triplet state
(a $\ket{+}$ in the qubit space) by charging two electrons into one dot of
each qubit and tuning the intra-qubit tunnel coupling to the $(1,1)$ charge
regime, i.e. one electron per dot. Then the graph is initialized by turning on
exchange interactions between the desired qubits. This encoding of the qubit-
space is key to exploiting fast QND measurement methods measuring in the
singlet-triplet basis for measurement-based qubit control of the cluster
state.

###  III.2 Exchange Mediated Rotated Measurements

Consider a simple two DQD spin-qubit example, the first qubit in some
arbitrary qubit state
$\ket{\psi}=\\{\cos(\theta/2),e^{i\phi}\sin(\theta/2)\\}$, the second in a
triplet state with exchange interaction as described by Eq. (3), again
assuming $\Delta B=\Delta_{12}=\Delta_{34}=\Delta$. At some time $\tau$ the
qubit initially containing $\ket{\psi}$ is measured projectively in the
singlet-triplet basis, using the schemes discussed in Sec. (V) with outcome
$s=0,1$ for a triplet or singlet measurement respectively. The resulting
unnormalized state on the second qubit is

$\ket{\tilde{\psi}}=\frac{1}{\sqrt{2}}\left(\begin{matrix}e^{-2i\tau\mu\Delta}&(-1)^{s}\\\
1&(-1)^{s}e^{-i\tau(J_{23}-2\mu\Delta)}\end{matrix}\right)\ket{\psi}.$ (7)

Using an exchange $J_{23}=(2n+1)\pi/\tau$ where $n\in\mathbb{Z}$, this
simplifies to

$\ket{\tilde{\psi}}=\sigma_{x}^{s}P\left((-1)^{s}\theta[\tau]\right)W\left(\theta[\tau]\right)\ket{\psi}$
(8)

where

$\theta[\tau]=2\tau\mu\Delta.$ (9)

The resulting projected unitary is the desired exchange controlled single
qubit $W(\theta)$ gate, repetitions of which allow for full single qubit
control, with additional measurement dependent phase
$P\left((-1)^{s}\theta\right)$ and Pauli $\sigma_{x}^{s}$ corrections. This
additional phase may be corrected after-the-fact by intra-qubit phase gates
due to magnetic gradient precession. Importantly however, the $W(\theta)$
control gate is not dependent on the measurement outcome, only the phase
correction. The gate time of this projective measurement-based gate depends on
the angle $\theta$ of the gate. As a benchmark, to project a Hadamard within
this scheme the gate time is $\sim 20\,\mathrm{ns}$Ruskov and Tahan (2019);
Wardrop and Doherty (2014). It is important to note here that after the
application of this protocol, the qubit that initially contained the encoded
state is projected onto either a triplet or singlet state that can be
immediately reused in the same protocol. This notion of re-entangling spent
qubits of the cluster state helps mediate the overall physical qubit cost of
building a processor based on this protocol. This is explored further in Sec.
III.5.

###  III.3 Efficient Arbitrary State Preparation

Here the exchange-measurement hybrid protocols outlined in Sec. III.2 shall be
demonstrated to be useful for single measurement arbitrary state preparation.
This is achieved with two qubits, assumed to both start in the spin triplet
state ($\ket{+}$ in the qubit basis). The resulting unnormalized teleported
state is

$\begin{split}&\ket{\psi}=\\\ &\left(\begin{matrix}\cos\left[\frac{1}{2}(\pi
s+\theta)\right]\\\
e^{-\frac{1}{2}i(J_{23}\tau-2(2\mu\Delta\tau+\phi_{\Delta}))}\cos\left[\frac{1}{2}(-J_{23}\tau+\pi
s+\theta)\right]\end{matrix}\right)\end{split}$ (10)

where $\theta=2\mu\Delta\tau+\phi_{\Delta}$. At times $\tau=\pi(2n+1)/J_{23}$
where $n\in\mathbb{Z}$, this simplifies to the normalized state

$\ket{\psi}=\left(\begin{matrix}\cos\left[\frac{1}{2}(\pi s+\theta)\right]\\\
-ie^{i\theta}\sin\left[\frac{1}{2}(\pi s+\theta)\right]\end{matrix}\right).$
(11)

which can be made to be any arbitrary single qubit state by application of a
phase gate.

With this protocol any arbitrary single qubit state can be written with a
measurement-based protocol with a single measurement. However, the state is
dependent upon the measurement outcome $s$. The two possible states generated
by the protocol are always antipodal to one another on the Bloch sphere. To
recover the $s=0$ state from a $s=1$ measurement requires a Pauli and phase
correction $\sigma_{x}$ and appropriate phase.

###  III.4 Two-qubit Gates

$\Qcircuit@C=1em@R=1em{\raisebox{-50.00008pt}{$(a)$}&}$

$\Qcircuit@C=1em@R=1em{\raisebox{-50.00008pt}{$(b)$}&}$

$\Qcircuit<EMAIL_ADDRESS>\lstick{\ket{\psi_{2}}}\targ\gate{H}\gate{X^{s_{2}}}\qw}$

Figure 2: (a) Four-qubit square geometry to perform an entangling two qubit
gate in exchange-measurement hybrid cluster state protocol. The initial two
qubits start in the top right and bottom left qubits, all qubits are then
entangled by exchange to their two neighbors, and singlet-triplet
measurements, equivalent to projective Pauli-$x$ measurements, are performed,
resulting in an entangling two qubit state in the remaining two qubits. (b)
The equivalent gate unitary, including Pauli corrections of the exchange-
measurement hybrid protocol given in (a). Here the bottom right qubit in (a)
is the first logical qubit and the top left is the second. The measurement
outcomes $s_{1}$ and $s_{2}$ corresponds to the measurements performed on the
physical qubits that initially encoded $\ket{\psi_{1}}$ and $\ket{\psi_{2}}$
respectively.

In a cluster state, to perform two qubit gates, 2D graphs are necessary. With
spin qubits, 2D geometries are limited in scope due to considerations such as
plunger-gate placements and cross-talk. However, one of the simplest
geometries, four qubits arranged in a square, can be used with the hybrid
exchange-measurement graph protocol described here to perform entangling two
qubit gates. An example of this is given in Fig. 2. In this protocol, all the
qubits are identical, with magnetic field gradients $\Delta$. Initially the
two logical qubits start in the top right and bottom left physical qubits of
the square, and the remaining qubits in the $\ket{+}$ state. Inter-qubit
interaction is then turned on between all neighboring qubits with time $\tau$
at which point the qubits that previously housed the logical states are
measured. The resulting states encoded in the two qubits that weren’t measured
depend on the values of interaction strength used for each of the four inter-
qubit exchange gates. For example, if the value of the interaction strength in
Fig. 2(a) used is $J=2(2n_{1}-1)\mu\Delta$ for times
$\tau=(2n_{2}-1)\pi/(2\mu\Delta)$ where $\\{n_{1},n_{2}\\}\in\mathbb{Z},\neq
0$, the resulting two qubit gate performed by this scheme is given in Fig.
2(b), up to two Pauli-$x$ corrections.

###  III.5 Re-entangling Measured Qubits

Figure 3: Example of a simple 2-qubit recycled cluster state with accumulated
Pauli corrections. Initially the left hand (LH) qubit is prepared in some
state $\ket{\psi}$ and the right hand (RH) qubit is in either a singlet or
triplet state. The protocol given in (8) is then applied to both qubits
measuring the LH qubit leaving it in either a singlet or triplet state and the
RH qubit in the state given. The protocol is then applied again measuring the
RH qubit, adjusting the rotation $\theta^{\prime}$ of the second application
of the protocol based on the phase correction $\theta$ of the previous
iteration of the protocol.

Projective measurements, like that assumed in this work, leave behind a qubit
prepared in a singlet or triplet state which may be re-entangled to the
resource state. By again engaging the exchange interaction $J_{23}$ the qubit
may be re-entangled with the information carrying qubit to allow for
repetition of a measurement-based gate on the same geometry. This circumvents
the issues of resource overhead for solid state implementations of MBQC for
non-commuting measurements. The adjustments needed to do so depend upon the
previous measurement result. As such the projected gate given by the protocol
(8) becomes

$\ket{\tilde{\psi}}=\sigma_{z}^{s_{1}}\sigma_{x}^{s_{2}}P\left((-1)^{s_{2}}\theta\right)W\left(\theta\right)\ket{\psi}$
(12)

where $s_{1}$ is the result of the previous measurement performed on the qubit
the state is being projected to and $s_{2}$ is the result of the current
measurement projecting the state back onto the first qubit. Performing
multiple gates with this protocol, the overall number of post gate
corrections, i.e. up to two Pauli corrections and one phase correction,
remains constant regardless of the number of iterations if the previous known
anomalous phase is considered when performing the next gate. For example,
assuming an initial two qubit state $\ket{\psi}\otimes\ket{(-1)^{s_{1}}+}$, if
the protocol (8) is repeated twice for two different angles $\theta$ and
$\theta^{\prime}$ respectively, the applied unitary is

$U=\sigma_{x}^{s_{1}+s_{3}}P\left((-1)^{s_{1}+s_{3}}\theta^{\prime}\right)W\left((-1)^{s_{2}}\left(\theta^{\prime}+\theta\right)\right)W\left(\theta\right).$
(13)

This shows that subsequent rotations $\theta^{\prime}$ may be adjusted to
account for the measurement outcomes of the previous iteration of the protocol
to apply the desired series of unitaries without performing a phase correction
at each step. Note that if the protocol is applied an odd number of times,
such as in (13), a maximum of one Pauli-$x$ correction may be needed, whilst
if the protocol is applied an even number of times, a maximum of Pauli-$x$ and
an Pauli-$z$ correction may be needed.

Figure 4: (a) Linear cluster state geometry with measurements (magenta) to
perform an arbitrary qubit rotation $U$ on state $\ket{\psi}$ initially
encoded in qubit $Q1$. The final state $\ket{\tilde{\psi}}=U\ket{\psi}$ is
teleported to qubit $Q5$. The black dashed lines denote exchange interaction
to allow for necessary measurements. (b) Graph state geometry to encode the
state $\ket{\tilde{\psi}}$ to the $[[4,1,2]]$ stabilizer state given by the
solid black lines in (c), by measuring $Q3$. (d) Protocol for decoding the
stabilizer state given in (c), projected $\ket{\tilde{\psi}}$ to $Q1$.

With an architecture in which each qubit can be measured and re-entangled to
the graph on demand, a processor which can be used both for qubit operations
and quantum error correction encoding can be designed. Graph states are
naturally compatible with stabilizer codes as the stabilizers $K_{i}$ of a
cluster state may be inferred from its geometry as

$K_{i}=(\sigma_{x})_{i}\otimes_{j\in N(i)}(\sigma_{z})_{j}.$ (14)

For example, encoding and decoding a qubit with a $[[4,1,2]]$ stabilizer state
has been demonstrated experimentally with a five qubit cluster state in an
optical systemBell  _et al._ (2014b). The graph geometry used for this code is
not feasibly compatible with a spin-qubit implementation. However, a similar,
linear code can be used, at the cost of distinguishing some single qubit Pauli
errors. This is achieved with an architecture of five qubits in a line. The
same architecture is described in Sec. II to achieve arbitrary single qubit
rotations of the form

$\ket{\tilde{\psi}}=W(0)W(\theta_{3})W(\theta_{2})W(\theta_{1})\ket{\psi}$
(15)

with the rotated state ending in the fifth physical qubit. This state may then
be transferred back to the middle (third) qubit by two $W(0)$ measurements.
Then, by entangling all adjacent qubits with a $CZ$ interaction and measuring
the third qubit in the Pauli-$x$ basis, $\ket{\tilde{\psi}}$ is then encoded
in a $[[4,1,2]]$ stabilizer state across the remaining four qubits. The qubit
can then be recovered to the first qubit form the stabilizer to perform
additional rotations by measuring the second, fourth and fifth qubit in the
Pauli-$x$. This protocol, depicted in Fig. 4, has stabilizers
$\\{S_{1},S_{2},S_{3}\\}=\\{\sigma_{x}\otimes\sigma_{z}\otimes\mathbb{I}\otimes\mathbb{I},\sigma_{z}\otimes\sigma_{x}\otimes\sigma_{x}\otimes\sigma_{z},\mathbb{I}\otimes\mathbb{I}\otimes\sigma_{z}\otimes\sigma_{x}\\}$,
and can detect any arbitrary single qubit Pauli error, though cannot
distinguish between all such errors, see Fig. 5.

Figure 5: Table of the expectation values of the stabilizers
$\\{S_{1},S_{2},S_{3}\\}$ and products of the linear stabilizer code given in
Fig. 4 under every single qubit Pauli on the $i$th error
$\\{X_{i},Y_{i},Z_{i}\\}$. The occurrence of any single qubit error is
detectable, but some are indistinguishable from a single measurement of the
syndromes.

##  IV Entangling Measurements Protocol

The methods presented up to this point are a form of MBQC whereby the
measurements provide unitary operations on the encoded qubits, at the cost of
entanglement between qubits generated by exchange gates. This allows for a
powerful method of processing and encoding quantum information with an at
worse polynomial cost in physical qubits, which can be mediated by re-
entangling measured qubits. The following protocol shown here, however,
implements measurements that entangle qubits, as well as measurements that
rotate qubits at the cost of entanglement.

###  IV.1 Qubit Model: No Magnetic Gradients

In this scheme qubits are encoded in the Pauli-$z$ basis of the singlet-
triplet spin space
($\ket{0}=\left(\ket{\uparrow\downarrow}-\ket{\downarrow\uparrow}\right)/\sqrt{2}$
and
$\ket{1}=\left(\ket{\uparrow\downarrow}+\ket{\downarrow\uparrow}\right)/\sqrt{2}$)
and require only controllable exchange coupling between dots within the same
qubit. No magnetic field gradients or inter-qubit exchange coupling are
required. Therefore, the effective Hamiltonian for the inter- and intra-qubit
exchange interaction of two qubits comprised of four dots (as in Fig. 1(b))
can be written asLi  _et al._ (2012)

$H_{\text{eff}}=\frac{J_{12}}{2}\sigma_{z}^{1}+\frac{J_{34}}{2}\sigma_{z}^{2}-\frac{J_{23}}{4}\sigma_{x}^{1}\sigma_{x}^{2}.$
(16)

###  IV.2 Inter-Qubit Measurement

Alternatively to measuring in the basis of a logical qubit, an entangling
operation can be achieved by measuring a dot from each adjacent encoded qubit,
i.e. measuring in the singlet-triplet basis of the second and third QDs in
Fig. 1. Forgoing any exchange interaction, the resulting 2-qubit projector in
the qubit basis used in Sec. III is

$\mathcal{M}=\ket{00}\bra{00}+\ket{11}\bra{11}.$ (17)

However, this measurement also entangles the encoded qubits to spin-states
outside of the qubit space ($\ket{\uparrow\uparrow}$ and
$\ket{\downarrow\downarrow}$), i.e. it allows for leakage of the two qubits.
The measurement cannot however change the net-spin state of the two qubits it
acts on and, as such, measuring either of the individual qubits in the
singlet-triplet basis corrects for the leakage, at the cost of the
entanglement of the qubits.

A simple example of a use of the projector in Eq. (17), in this qubit basis,
is as a quantum bus. This is achieved by encoding qubit one in some arbitrary
state $\ket{\psi}$ and qubit two in a $\ket{+}$ or singlet state. By first
performing the entangling measurement between the two qubits, then measuring
qubit one, the resulting state given to the second qubit is
$\sigma_{z}^{s}\ket{\psi}$ where $s$ is the result of the second measurement
and the leakage is corrected for. With the exchange interaction and phase
gates considered in Sec. III, the projected state after both measurements
includes a phase gate
$\ket{\tilde{\psi}}=\sigma_{z}^{s}P\left(-J_{23}\tau+4\mu\Delta
t+2\phi_{\Delta}\right)\ket{\psi}$. In this basis however, full qubit control
cannot be performed with measurements of this type.

###  IV.3 Qubit Control with Entangling Measurements

The entangling measurement described can be used along with intra-qubit
exchange-interaction to implement an alternative non-cluster state approach to
MBQC. In this regime it can be assumed that phase gates of the formLi  _et
al._ (2012)

$P(J,\tau)=\left(\begin{matrix}1&0\\\ 0&e^{iJ\tau}\end{matrix}\right)$ (18)

where the qubit with internal exchange-interaction $J$ is evolved for a period
of time $\tau$. The final tool needed for single and two qubit control is
access to an ancila qubit initialized in the
$\ket{+}=\ket{\uparrow\downarrow}$ state for each logical qubit (polynomial
increase in physical qubits) and the ability to measure each qubit in the
singlet-triplet basis.

$\Qcircuit@C=1em@R=1em{\lstick{\ket{\psi}}&\gate{P_{\frac{\pi}{2}}}\multigate{1}{\mathcal{M}}\gate{P_{\frac{\pi}{2}}}\gate{Y^{s}}\qw\rstick{H\ket{\psi}}\\\
\lstick{\ket{+}}\gate{P^{\dagger}_{\frac{\pi}{2}}}\ghost{\mathcal{M}}\qw\meter\rstick{s}}$
Figure 6: Circuit diagram of the phase gates and measurement sequence needed
to achieve a single qubit Hadamard gate up to a Pauli-$y$ correction.

With these tools, single qubit unitaries can be produced as follows. Suppose
initially there are two qubits, the first initialized to some arbitrary state
$\ket{\psi}$ and the second in the $\ket{+}$ state. Then both qubits
experience an exchange phase gate imparting phases $\phi_{1}$ and $\phi_{2}$
to the respective qubits, before which an entangling measurement is performed
on the qubits and after which another set of phase gates $\phi_{3}$ and
$\phi_{4}$ are performed before the second qubit is measured with outcome $s$.
This results in the following unitary acting on the input state in the first
qubit

$\ket{\tilde{\psi}}=\frac{1}{\sqrt{2}}\left(\begin{matrix}1&e^{i(\phi_{1}+(-1)^{s}\phi_{2})}\\\
e^{i((-1)^{s}\phi_{1}+\phi_{2})}&e^{i(\phi_{2}+\phi_{3})}\end{matrix}\right)\ket{\psi}.$
(19)

Then by selecting appropriate phases, useful gates up to some known correction
may be implemented. For example, Fig. 6 shows the sequence of phase gates and
measurements needed to perform a Hadamard gate on the logical qubit
$\ket{\psi}$ up to a Pauli-$y$ correction. The only measurement outcomes that
determine the correction is the ancila qubit measurement $s$ at the end of the
sequence; the entangling measurement result does not contribute. The same
tools can be used to implement a two qubit entangling gate with two ancila
qubits. The measurement and phase gate sequence needed is given in Fig. 7(a).
The resulting two qubit unitary of the sequence given in Fig. 7(a) is shown in
Fig. 7(b) including the single Pauli-$x$ correction given by the measurement
results $s_{1}+s_{2}$ of the measurement of the ancila. This is a
demonstration of implementing a Clifford gate group with the proposed scheme.
This is done for simplicity to demonstrate that a universal gate set can be
implemented, however, non-Clifford two qubit gates are possible with the same
tools.

For the single qubit gate scheme the entangling measurement with the
measurement on the ancila qubit does not allow for leakage of the logical
qubit, regardless of the phase gates used and the or ancila qubit initial
state. The two qubit gate example given also surpreses leakage, however, this
is due to the selection of the phase gates, initial ancila states, and
measurement sequence and is therefore a significant factor to consider when
designing gates within this scheme. The gate times of both the single and two
qubit gates are limited by the single qubit phase gates. Assuming an exchange
coupling of $J=160\,\mathrm{MHz}$Reed  _et al._ (2016) gate times of $\sim
80\,\mathrm{ns}$ for the Hadamard and $\sim 140\,\mathrm{ns}$ for the rotated
$CZ$ gates shown in Fig. 6 and Fig. 7 respectively.

This scheme offers a method of hybrid exchange measurement-based approach to
quantum computation similar to a traditional gate-based approach without the
need for magnetic field gradients or driving fields and can be made to be
robust to leakage as well as operated at the SOP of the physical qubits. The
main limitation is the need for ancila qubits initialized to the $\ket{+}$
state as this method offers no way of preparing such states from a $\ket{0}$
or $\ket{1}$ state. Such states can be initialized by adiabatically turning
off the exchange interactionFoletti  _et al._ (2009); Fogarty  _et al._ (2018)
of a qubit in the triplet state, which is within the control parameters of the
scheme. Unfortunately however, since this scheme is more akin to traditional
gate-based approaches to quantum computations, the inherent capability of the
scheme with stabilizer codes described in Sec. III.5 does not apply to this
scheme.

$\Qcircuit@C=0.4em@R=1em{\lstick{(a)}\\\
\lstick{\ket{+}}&\gate{P^{\dagger}_{\frac{\pi}{2}}}\qw\qw\qw\qw\multigate{1}{\mathcal{M}_{12}}\qw\qw\qw\raisebox{-35.00005pt}{$\ket{\tilde{\Psi}}$}\\\
\lstick{\ket{\psi_{1}}}\qw\qw\qw\multigate{1}{\mathcal{M}_{23}}\gate{P^{\dagger}_{\frac{\pi}{2}}}\ghost{\mathcal{M}_{12}}\gate{Z}\qw\qw\\\
\lstick{\ket{+}}\gate{P_{\frac{\pi}{2}}}\multigate{1}{\mathcal{M}_{34}}\qw\ghost{\mathcal{M}_{23}}\qw\qw\gate{Z}\meter\rstick{s_{1}}\\\
\lstick{\ket{\psi_{2}}}\gate{Z}\ghost{\mathcal{M}_{34}}\gate{Z}\qw\qw\qw\qw\meter\rstick{s_{2}}}$

$\Qcircuit@C=1em@R=1em{\lstick{(b)\ket{\psi_{1}}}&\gate{H}\ctrl{1}\gate{H}\qw\qw\raisebox{-30.00005pt}{$\ket{\tilde{\Psi}}$}\\\
\lstick{\ket{\psi_{2}}}\gate{H}\gate{Z}\gate{H}\gate{X^{s_{1}+s_{2}}}\qw}$

Figure 7: (a) Circuit diagram of the phase gates and measurement sequence
needed to achieve a two qubit entangling unitary with two ancila qubits. (b)
Circuit diagram of the two qubit unitary performed by the measurement and
phase sequence given in (a) with single Pauli-$x$ correction.

##  V Measurement Methods

In this work two schemes for MBQC with singlet-triplet spin qubits are given
assuming only measurement of the singlet-triplet basis. This can be achieved
in a number of different ways, however here we shall consider two such
methods. The conventional spin-to-charge method and the QND curvature method.

The conventional spin-to-charge method of measuring in the singlet-triplet in
achieved by adiabatically tuning the tunnel coupling between the two dots
being measured. Due to Pauli-exclusion, the singlet state allows for hopping
from the $(1,1)$ to the $(0,2)$ while the triplet state does not couple to the
doubly occupied stateHanson  _et al._ (2007). If the singlet state is
measured, it can then be re-initialized back to the $(1,1)$ singlet state by
again adiabatically tuning the tunnel coupling.

Alternatively, the singlet-triplet basis can be measured by coupling to a
superconducting (SC) resonatorRuskov and Tahan (2019). The resonator is driven
with frequencies $\omega_{m}\approx\omega_{r}$ and
$\omega_{s}\approx\omega_{q}$ with strength $\epsilon_{s}$, where $\omega_{r}$
is the frequency of the resonator, $\omega_{m}$ is the frequency of an
external gate modulation near the resonator and $\omega_{q}$ is the singlet-
triplet qubit frequency. This scheme exploits the quantum curvature of the
qubits and is described as follows in the rotating frame and qubit-space of
the qubit-qubit interaction Hamiltonian (3)

$\begin{split}H_{\text{meas}}=(\omega_{r}-\omega_{m})a^{\dagger}a+\frac{(\omega_{q}-\omega_{s})}{2}\sigma_{z}+&\\\
(\chi+\delta\omega)a^{\dagger}a\sigma_{z}+\frac{1}{2}(g_{\parallel}\sigma_{z}+g_{0}\sigma_{0})(a+a^{\dagger})&\\\
+\epsilon_{d}(ae^{-i\omega_{m}t}+a^{\dagger}e^{i\omega_{m}t})&\end{split}$
(20)

where $a(a^{\dagger})$ are the annihilation (creation) operators of the
resonator, $g_{\parallel}$ is the longitudinal coupling between resonator and
qubit and $g_{0}$ the static coupling. Here the important term for measuring
the state of the qubit is the $\Delta\omega$ term. This term is the change in
the SC resonator frequency due to the state of the qubit. This term is given
by the quantum curvature of the eigen-energies of the qubit and is defined as

$\delta\omega=\hbar\omega_{r}^{2}\left(\frac{\eta}{e_{0}}\right)^{2}\frac{\partial^{2}G(V_{m}^{0})}{\partial
V_{m}^{2}}$ (21)

where $\omega_{r}$ is the resonator frequency, $\eta$ is the dimensionless
parameter for resonator inductance and $G(V_{m}^{0})$ is defined as follows,
assuming the singlet-triplet qubit Hamiltonian as

$\begin{split}H&=E_{s}\ket{s}\bra{s}+E_{t_{0}}\ket{t_{0}}\bra{t_{0}}\\\
&=G(V_{m}^{0})\mathbb{I}+C(V_{m}^{0})\sigma_{z}.\end{split}$ (22)

Measuring $\delta\omega$ is assumed to quickly perform a projective
measurement in the singlet-triplet spin basis when the cavity is driven.

Both methods of measurement can be used for the two MBQC protocols given here.
The primary difference between the two is the additional step needed by the
spin-to-charge method to re-initialize the $(1,1)$ charge state after a
singlet state is measured. This can be limiting on the speed of processes that
reuse physical qubits, such as the stabilizer state encoding, and the
entangling measurement protocol. Both methods are viable, however the
curvature method is optimal. All gate times given in this paper are calculated
assuming the curvature method.

##  VI Discussion

Two protocols for full measurement-based qubit control with Si DQD qubits
using two hybrid measurement-exchange schemes are given. Tuning of the
exchange interaction between two adjacent qubits allows for a Pauli-$x$
measurement in the qubit space to project an equivalent unitary, as is
achieved by an arbitrary $M(\theta)$ measurement in a cluster state, up to a
Pauli and phase correction. This allows for full single and two-quit control
as a measurement-based scheme with phase corrections by magnetic precession of
the individual qubits. This additionally allows for the subsequent QDs of the
qubit to be set up such that $\Delta B\gg J_{23}$, thus suppressing leakage of
the individual qubits to spin $\pm 1$ states. This scheme is shown to also be
useful for efficient arbitrary qubit state preparation. The nature of the
measurements allows for the limiting of the required qubit resource overhead
by recycling measured qubits back into the resource state with additional
known required Pauli corrections. This allows for novel protocols for the
arbitrary manipulation as well as encoding and decoding of a qubit state into
stabilizer state with five physical qubits in a 1D geometry.

A second protocol with entangling measurements is also shown by measuring in
the singlet-triplet basis between one dot from one qubit with one dot from a
different qubit. This measurement-based protocol is closer to a traditional
gate-based protocol, and can be made to be robust against leakage and can be
operated at the SOP of all the physical qubits, requiring only measurements
and phase gates.

In the first protocol, a single single qubit gate consists of one exchange
gate and one measurement, whilst a two qubit control consists of four exchange
gates and two measurements, all performed simultaneously. In the second
protocol single qubit gates consist of three exchange gates and two
measurements whilst two qubit gates consist of seven exchange gates, five
measurements and $n$ ancilla qubits per $n$ data qubits (total $2n$ physical
qubits). The first proposed protocol is an improvement compared to
conventional gate-based leakage resistant schemes in singlet-triplet encoded
qubits with a magnetic field gradient, where single qubit leakage errors can
be corrected with a minimum of five exchange gates and one ancilla qubitMehl
_et al._ (2015), and two qubit gates are assumed to be given by long range
cavity coupling, unlike in the work presented. Additionally, the second
proposed protocol is more efficient than leakage protection and correction in
similar systems without magnetic field gradients, three-spin exchange-only
qubits. In such systems, active leakage reduction requires 30 exchange
gatesFong and Wandzura (2011) and protected two qubit gates require 16
exchange gatesSetiawan  _et al._ (2014). A direct comparison of these schemes
is made in Table 1. Equally, both schemes are superior when compared to
superconducting systems that can also experience leakage errors, requiring 16
single and two qubits gates and two measurementsGhosh and Fowler (2015).

| Protocol 1 | $S/T_{0}$ | Protocol 2 | 3-Spin
---|---|---|---|---
Single Qubit Gate
Gates | 1 | 1 | 3 | 4
Measurements | 1 | - | 2 | -
Ancillae | 1 | 0 | 1 | 0
Leakage Protected | Yes | No | Yes | No
Correction Gates | - | 7 | - | 30
Correction Ancillae | - | 1 | - | 1
Two Qubit Gate
Gates | 4 | 1 | 7 | 16
Measurements | 2 | - | 5 | -
Ancillae | 2 | 0 | 2 | 0
Leakage Protected | Yes | No | Yes | Yes
Correction Gates | - | 14 | - | -
Correction Ancillae | - | 2 | - | -
Table 1: Table comparing the number of exchange gates, timesteps, measurements
and ancillae between the first proposed cluster state protocol, and leakage
protected operations in singlet-triplet qubits with magnetic field gradients
from literatureMehl  _et al._ (2015) and the second entangling measurements
protocol and gradient free exchange only three spin-encoded qubits from
literatureFong and Wandzura (2011); Setiawan  _et al._ (2014). Single qubits
operations from literature aren’t robust to leakage, so given numbers include
leakage correction circuits after a rotation. There are no two qubit leakage
robust sequences in the literature, only leakage suppressed, so for comparison
the leakage correction circuit is considered twice. Both proposed protocols
offer significant improvements over their similar respective counterparts.

This work promises alternative hybrid exchange-measurement-based quantum
computation schemes. The estimated gate times for all the schemes proposed
here are $<150\,\mathrm{ns}$ and so are all well within expected relaxation
times of Si DQD spin qubits of tens of microsecondsSigillito  _et al._ (2019).
This allows for semiconductor qubits to take advantage of MBQC without the
resource overhead or necessity for arbitrary rotated measurement schemes.

##  VII Acknowledgements

We acknowledge helpful discussions with M. Benito, T. Ladd, R. Ruskov, M.
Russ, YP. Shim, V. Shkolnikov and Y. Yanay.

## References

  * Raussendorf and Briegel (2000) R. Raussendorf and H. J. Briegel, arXiv preprint quant-ph/0010033 (2000).
  * Verstraete and Cirac (2004) F. Verstraete and J. I. Cirac, Phys. Rev. A 70, 060302 (2004).
  * Raussendorf  _et al._ (2002) R. Raussendorf, D. Browne, and H. Briegel, J. Mod. Opt. 49, 1299 (2002).
  * Nielsen (2001) M. Nielsen, arXiv preprint quant-ph/0108020 (2001).
  * Raussendorf (2009) R. Raussendorf, International Journal of Quantum Information 7, 1053 (2009).
  * Jozsa (2006) R. Jozsa, NATO Science Series, III: Computer and Systems Sciences. Quantum Information Processing-From Theory to Experiment 199, 137 (2006).
  * Browne and Briegel (2019) D. Browne and H. Briegel, “One-way quantum computation,” in  _Quantum Information_ (John Wiley & Sons, Ltd, 2019) Chap. 21, pp. 449–473.
  * Raussendorf and Harrington (2007) R. Raussendorf and J. Harrington, Phys. Rev. Letts. 98, 190504 (2007).
  * Van den Nest  _et al._ (2006) M. Van den Nest, A. Miyake, W. Dür, and H. J. Briegel, Phys. Rev. Letts. 97, 150504 (2006).
  * Plato  _et al._ (2008) A. D. K. Plato, O. C. Dahlsten, and M. B. Plenio, Phys. Rev. A 78, 042332 (2008).
  * Shettell and Markham (2020) N. Shettell and D. Markham, Phys. Rev. Letts. 124, 110502 (2020).
  * Barends  _et al._ (2014) R. Barends, J. Kelly, A. Megrant, A. Veitia, D. Sank, E. Jeffrey, T. C. White, J. Mutus, A. G. Fowler, B. Campbell, _et al._ , Nature 508, 500 (2014).
  * Arute  _et al._ (2019) F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brandao, D. A. Buell, _et al._ , Nature 574, 505 (2019).
  * Preskill (2018) J. Preskill, Quantum 2, 79 (2018).
  * Wang  _et al._ (2018) Y. Wang, Y. Li, Z.-q. Yin, and B. Zeng, npj Quantum Information 4, 1 (2018).
  * Zhang  _et al._ (2017) J. Zhang, G. Pagano, P. W. Hess, A. Kyprianidis, P. Becker, H. Kaplan, A. V. Gorshkov, Z.-X. Gong, and C. Monroe, Nature 551, 601 (2017).
  * Ballance  _et al._ (2016) C. Ballance, T. Harty, N. Linke, M. Sepiol, and D. Lucas, Phys. Rev. Letts. 117, 060504 (2016).
  * Hanson  _et al._ (2007) R. Hanson, L. P. Kouwenhoven, J. R. Petta, S. Tarucha, and L. M. Vandersypen, Rev. Mod. Phys. 79, 1217 (2007).
  * Russ and Burkard (2017) M. Russ and G. Burkard, Journal of Physics: Condensed Matter 29, 393001 (2017).
  * Mi  _et al._ (2018) X. Mi, M. Benito, S. Putz, D. M. Zajac, J. M. Taylor, G. Burkard, and J. R. Petta, Nature 555, 599 (2018).
  * Veldhorst  _et al._ (2017) M. Veldhorst, H. Eenink, C. Yang, and A. S. Dzurak, Nat. Comms. 8, 1 (2017).
  * Takeda  _et al._ (2020) K. Takeda, A. Noiri, J. Yoneda, T. Nakajima, and S. Tarucha, Phys. Rev. Letts. 124, 117701 (2020).
  * Yang  _et al._ (2019) C. Yang, K. Chan, R. Harper, W. Huang, T. Evans, J. Hwang, B. Hensen, A. Laucht, T. Tanttu, F. Hudson, _et al._ , Nat. Electron. 2, 151 (2019).
  * Zajac  _et al._ (2018) D. M. Zajac, A. J. Sigillito, M. Russ, F. Borjans, J. M. Taylor, G. Burkard, and J. R. Petta, Science 359, 439 (2018).
  * Russ  _et al._ (2018) M. Russ, D. M. Zajac, A. J. Sigillito, F. Borjans, J. M. Taylor, J. R. Petta, and G. Burkard, Phys. Rev. B 97, 085421 (2018).
  * Reed  _et al._ (2016) M. Reed, B. Maune, R. Andrews, M. Borselli, K. Eng, M. Jura, A. Kiselev, T. Ladd, S. Merkel, I. Milosavljevic, _et al._ , Phys. Rev. Letts. 116, 110402 (2016).
  * Yoneda  _et al._ (2018) J. Yoneda, K. Takeda, T. Otsuka, T. Nakajima, M. R. Delbecq, G. Allison, T. Honda, T. Kodera, S. Oda, Y. Hoshi, _et al._ , Nat. Nanotechnol. 13, 102 (2018).
  * DiVincenzo  _et al._ (2000) D. P. DiVincenzo, D. Bacon, J. Kempe, G. Burkard, and K. B. Whaley, Nature 408, 339 (2000).
  * Bacon  _et al._ (2000) D. Bacon, J. Kempe, D. A. Lidar, and K. B. Whaley, Phys. Rev. Letts. 85, 1758 (2000).
  * Levy (2002) J. Levy, Phys. Rev. Letts. 89, 147902 (2002).
  * Johnson  _et al._ (2005) A. Johnson, J. Petta, J. Taylor, A. Yacoby, M. Lukin, C. Marcus, M. Hanson, and A. Gossard, Nature 435, 925 (2005).
  * Laird  _et al._ (2010) E. A. Laird, J. M. Taylor, D. P. DiVincenzo, C. M. Marcus, M. P. Hanson, and A. C. Gossard, Phys. Rev. B 82, 075403 (2010).
  * Wardrop and Doherty (2014) M. P. Wardrop and A. C. Doherty, Phys. Rev. B 90, 045418 (2014).
  * Andrews  _et al._ (2019) R. W. Andrews, C. Jones, M. D. Reed, A. M. Jones, S. D. Ha, M. P. Jura, J. Kerckhoff, M. Levendorf, S. Meenehan, S. T. Merkel, _et al._ , Nat. Nanotechnol. 14, 747 (2019).
  * Raussendorf  _et al._ (2006) R. Raussendorf, J. Harrington, and K. Goyal, Annals of physics 321, 2242 (2006).
  * Bell  _et al._ (2014a) B. Bell, D. Herrera-Martí, M. Tame, D. Markham, W. Wadsworth, and J. Rarity, Nat. Comms. 5, 1 (2014a).
  * Russo  _et al._ (2018) A. Russo, E. Barnes, and S. E. Economou, Phys. Rev. B 98, 085303 (2018).
  * Bodiya and Duan (2006) T. Bodiya and L.-M. Duan, Phys. Rev. Letts. 97, 143601 (2006).
  * Lanyon  _et al._ (2013) B. Lanyon, P. Jurcevic, M. Zwerger, C. Hempel, E. Martinez, W. Dür, H. Briegel, R. Blatt, and C. F. Roos, Phys. Rev. Letts. 111, 210501 (2013).
  * Karzig  _et al._ (2017) T. Karzig, C. Knapp, R. M. Lutchyn, P. Bonderson, M. B. Hastings, C. Nayak, J. Alicea, K. Flensberg, S. Plugge, Y. Oreg, _et al._ , Phys. Rev. B 95, 235305 (2017).
  * Xue (2013) Z.-Y. Xue, The European Physical Journal D 67, 89 (2013).
  * Bonderson and Lutchyn (2011) P. Bonderson and R. M. Lutchyn, Phys. Rev. Letts. 106, 130505 (2011).
  * Grimsmo and Smith (2019) A. L. Grimsmo and T. B. Smith, Phys. Rev. B 99, 235420 (2019).
  * Zwanenburg  _et al._ (2013) F. A. Zwanenburg, A. S. Dzurak, A. Morello, M. Y. Simmons, L. C. Hollenberg, G. Klimeck, S. Rogge, S. N. Coppersmith, and M. A. Eriksson, Rev. Mod. Phys. 85, 961 (2013).
  * McNeil  _et al._ (2010) R. P. McNeil, R. J. Schneble, M. Kataoka, C. J. Ford, T. Kasama, R. E. Dunin-Borkowski, J. M. Feinberg, R. J. Harrison, C. H. Barnes, D. H. Tse, _et al._ , Nano Letts. 10, 1549 (2010).
  * Ruskov  _et al._ (2018) R. Ruskov, M. Veldhorst, A. S. Dzurak, and C. Tahan, Phys. Rev. B 98, 245424 (2018).
  * Harvey-Collard  _et al._ (2019) P. Harvey-Collard, N. T. Jacobson, C. Bureau-Oxton, R. M. Jock, V. Srinivasa, A. M. Mounce, D. R. Ward, J. M. Anderson, R. P. Manginell, J. R. Wendt, _et al._ , Phys. Rev. Letts. 122, 217702 (2019).
  * Ruskov and Tahan (2019) R. Ruskov and C. Tahan, Phys. Rev. B 99, 245306 (2019).
  * Bell  _et al._ (2014b) B. A. Bell, D. A. Herrera-Martí, M. S. Tame, D. Markham, W. J. Wadsworth, and J. G. Rarity, Nat. Comms. 5 (2014b), 10.1038/ncomms4658, arXiv:1404.5498 .
  * Li  _et al._ (2012) R. Li, X. Hu, and J. You, Phys. Rev. B 86, 205306 (2012).
  * Foletti  _et al._ (2009) S. Foletti, H. Bluhm, D. Mahalu, V. Umansky, and A. Yacoby, Nat. Phys. 5, 903 (2009).
  * Fogarty  _et al._ (2018) M. Fogarty, K. Chan, B. Hensen, W. Huang, T. Tanttu, C. Yang, A. Laucht, M. Veldhorst, F. Hudson, K. M. Itoh, _et al._ , Nat. Comms. 9, 1 (2018).
  * Mehl  _et al._ (2015) S. Mehl, H. Bluhm, and D. P. DiVincenzo, Phys. Rev. B 91, 085419 (2015).
  * Fong and Wandzura (2011) B. H. Fong and S. M. Wandzura, arXiv preprint arXiv:1102.2909 (2011).
  * Setiawan  _et al._ (2014) F. Setiawan, H.-Y. Hui, J. Kestner, X. Wang, and S. D. Sarma, Phys. Rev. B 89, 085314 (2014).
  * Ghosh and Fowler (2015) J. Ghosh and A. G. Fowler, Phys. Rev. A 91, 020302 (2015).
  * Sigillito  _et al._ (2019) A. Sigillito, J. Loy, D. Zajac, M. Gullans, L. Edge, and J. Petta, Phys. Rev. Appl. 11, 061006 (2019).


Instead, when a process is executed only by one subcommunicator, MASKE
synchronizes all the other processors to recover the same system before
progressing with the loop.

When the chosen process is instead of KMC-free, discrete type, the $krun$
object first updates the $ratesT$ arrays. These arrays are similar to the
$rates$ array, in that each subcommunicator has one stored for $delete$ events
and one stored for $nucleate$ events in its submaster processor: see Fig. 13.
However, unlike the $rates$ arrays, the generic $i^{th}$ entry in the $ratesT$
arrays is obtained by taking its current value and adding the product between
$r_{k,i}$ and $\Delta t_{k}$, where $r_{k,i}$ is the corresponding rate value
stored in the $rates$ array. In other words, the $ratesT$ arrays store the
$\sum_{j=1}^{n_{cp}}{\Delta t^{*}_{j}r_{k,i,j}}$ terms in Eq. 22, and now the
latest contributions to them, $\Delta t_{k}r_{k,n_{cp}+1}$ in Eq. 22, is being
added. As a result, $krun$ can now select one discrete event to carry out with
probability $P_{i}$ as per Eq. 22.

The execution of a discrete event entails either the deletion or the creation
of a particle, in the current version of MASKE. This is carried out by all
subcommunicators at the same time. Following the stoichiometry specified in
the $chemDB$ file, the discrete event adds or removes molecules from the
implicit solution (see the $reax$ keyword in Section LABEL:secChemDB, and
specifically the $bkg$ attribute in it). Therefore the composition of the
implicit solution is modified accordingly, recomputing the new concentrations
of all the solvated molecules based on the new number of molecules in solution
and on their corresponding apparent volumes, which change the total volume of
the solution itself. Also this operation is carried out by all
subcommunicators. After having executed a discrete, KMC-free event, the
$ratesT$ arrays are zeroed across all subcommunicators, to restart the sums in
Eqs. 21 and 22 (note that the term $\sum_{j=1}^{n_{cp}}{\Delta
t^{*}_{j}R_{k,j}}$ in Eq. 21 is just the cumulative of all $ratesT$ arrays).

After any process is executed, be it continuous or discrete, the trial
particles associated to the various nucleation events (of KMC-free type) are
brought back to their initial positions. Then processes of type $Every$ are
run, provided that the correct, user-specified number of loops have been
completed since their last execution. Outputs are then printed out every so
many (user-specified) number of loops; the outputs include MASKE’s $thermo$
and $dump$ files described in C. This triggers the beginning of the next loop,
and thus the next step in MASKE. The loop ends when the simulated time $t^{*}$
is greater than the end time specified by the user via the input command
$Krun$: seeC.

##  Appendix F Sampling deletion and nucleation rates

In the current version of MASKE, the sampling of discrete KMC-free events is
limited to the possibility of deleting or creating (nucleating) a full
particle. When prescribing a deletion or nucleation process in the main input
script of MASKE, the user specifies a mechanism in the $ChemDB$ file to attach
to the process, which in turn invokes one of the user-specified chemical
reactions. A reaction involves a number of molecules in solution (called
’background’, in the $bkg$ attribute of the reaction) and a number of
molecules in the solid phase (called ’foreground’, via the $fgd$ attribute).
The number of reactions to nucleate a generic particle $j$ is the closest
integer to:

$n_{r,j,nuc}=\frac{V_{p,j}}{\sum_{i}\nu{{}_{i}}v_{m,i,fgd}}$ (24)

where $V_{p,j}$ is the volume of particle $j$ and $v_{m,i,fgd}$ is the volume
of the $i^{th}$ foreground molecule involved in the reaction, with $\nu_{i}$
its stoichiometric coefficient. For particle dissolution, the expression is
the same as in Eq. 24, except $-\nu_{i}$ should replace $\nu_{i}$.

For now, MASKE only features one mechanism, called $allpar$, whereby the
nucleation or deletion of an entire particle, irrespective of its size, is
fully determined by one single chemical reaction, assumed to take place at the
same time for all the solid molecules in the particle. This means that the
time to fully dissolve a particle, $\Delta t_{p,j}$ equals the time for one
single chemical reaction to take place, $r^{-1}_{1}$, where $r_{1}$ is the
rate of that chemical reaction. This mechanism is correct when a particle is
only as large as the sum of solid molecules involved in one reaction, i.e.
$V_{p,j}=\sum\limits_{i}\nu{{}_{i}}v_{m,i,fgd}$ leading to $n_{r,j,nuc}=1$ in
Eq. 24. For larger particles, the other extreme case would be to consider all
reactions taking place in series, in which case the particle nucleation or
dissolution time would be $\Delta
t_{p,j}=\sum\limits_{s=1}^{n_{r,j}}r_{s}^{-1}$, where $r_{s}$ is the rate of
the $s^{th}$ reaction in the series. More realistic mechanisms for larger
particles would involve a combination of reactions taking place in series and
in parallel, such as in classical and non-classical nucleation and
growth/dissolution theories. The interested reader can find an example of such
a mechanism in Shvab et al. (2017).

The $allpar$ mechanism in MASKE uses a modified version of Transition State
Theory (TST) for the rates of single dissolution/nucleation reactions, leading
to the unit rate expressions in Eqs. 3 and 4 in the main body of this article.
There it is pointed out how the $c^{o{\ddagger}}$ and $\Delta G^{o{\ddagger}}$
terms (standard concentration and activation free energy of the activated
complex) are related to the reaction rate constant $k$:

$k=\kappa\frac{k_{B}T}{h}\frac{c^{o{\ddagger}}}{\gamma^{\ddagger}_{d}}\exp\left(\frac{-\Delta
G^{o{\ddagger}}}{k_{B}T}\right)$ (25)

Eq. 25 indicates that the arbitrary choice of $c^{o{\ddagger}}$ has an effect
on the value of $\Delta G^{o{\ddagger}}$ to be employed. Consider a set of
experiments of dissolution or precipitation conducted at different
temperatures $T$. The Gibbs free energy can be expressed from its enthalpy $H$
and entropy $S$ terms: $\Delta G^{o{\ddagger}}=\Delta H^{o{\ddagger}}-T\Delta
S^{o{\ddagger}}$. One can thus rearrange Eq. 25 to obtain:

$\ln\left(\frac{k}{c^{0{\ddagger}}}\cdot\frac{h\gamma^{\ddagger}}{\kappa
k_{B}T}\right)=\frac{\Delta S^{o{\ddagger}}}{k_{B}}-\frac{\Delta
H^{o{\ddagger}}}{k_{B}}\frac{1}{T}$ (26)

If the $\exp\left(\frac{-\Delta H^{o{\ddagger}}}{k_{B}T}\right)$ term is
assumed to lead the temperature dependence, a set of experimental data
obtained at different $T$ would yield a linear plot of
$\ln\left(\frac{k}{c^{0{\ddagger}}}\cdot\frac{h\gamma^{\ddagger}}{\kappa
k_{B}T}\right)$ as a function of $T^{-1}$, with $\frac{\Delta
S^{o{\ddagger}}}{k_{B}}$ as the intercept when $T^{-1}\rightarrow 0$ and
$\frac{\Delta H^{o{\ddagger}}}{k_{B}}$ as the slope. Therefore, for a given
choice of $c^{o{\ddagger}}$, the linear plot provides the corresponding
$\Delta G^{o{\ddagger}}$. If $c^{o{\ddagger}}_{new}$ is a new, different
choice of $c^{o{\ddagger}}$, the numerical value of $k$ expressed in this new
units will change to $k_{new}$ (note from Eq. 5 that the dimensionality of
$c^{o{\ddagger}}$ determines that of $k$; for example, if $c^{o{\ddagger}}$ is
1 mol/m2 then $k$ will be in mol/m2 too). This will lead to a new linear plot
yielding different values of enthalpy and entropy, i.e. :

$\ln\left(\frac{k_{new}}{c^{0{\ddagger}}_{new}}\cdot\frac{h\gamma^{\ddagger}}{\kappa
k_{B}T}\right)=\frac{\Delta S^{o{\ddagger}}_{new}}{k_{B}}-\frac{\Delta
H^{o{\ddagger}}_{new}}{k_{B}}\frac{1}{T}$ (27)

By writing
$\frac{k_{new}}{c^{0{\ddagger}}_{new}}=\frac{k_{new}}{c^{0{\ddagger}}}\cdot\frac{c^{0{\ddagger}}}{c^{0{\ddagger}}_{new}}$,
one finds that the following relationships hold: $\Delta
H^{o{\ddagger}}_{new}=\Delta H^{o{\ddagger}}$ and $\Delta
S^{o{\ddagger}}_{new}=\Delta
S^{o{\ddagger}}+k_{B}\ln\left(\frac{c^{0{\ddagger}}}{c^{0{\ddagger}}_{new}}\right)$.
All this shows that the user must be careful to provide, as inputs, consistent
values for $c^{o{\ddagger}}$ and $\Delta G^{o{\ddagger}}$. The conversion of
$\Delta G^{o{\ddagger}}$ values from its original $c^{o{\ddagger}}$ to a new
$c^{0{\ddagger}}_{new}$ requires an additional entropy term due to the
difference in standard state concentrations:

$\Delta G^{o{\ddagger}}_{new}=\Delta
G^{o{\ddagger}}+k_{B}\ln\left(\frac{c^{0{\ddagger}}}{c^{0{\ddagger}}_{new}}\right)$
(28)

Eq. 28 confirms the intuition that more concentrated standard states have
higher entropy and thus lower Gibbs free energy. Furthermore, substituting Eq.
28 into Eq. 5 confirms that the value of the rate constant $k$, and thus also
the rates in Eqs. 3 and 4, do not depend on the arbitrary choice of standard
state concentration for the activated complex.

Going back to the rates in Eqs. 3 and 4, the activity products of the
reactants in dissolution and precipitation reactions, $Q_{r,d}$ and
$Q_{r,pr}$, can be defined more strictly than in the main body of this
article. Since solids have unit activity, only the molecules in solution
contribute to $Q\neq 1$, and these molecules are assigned to the ’background’
in MASKE. More precisely, for the generic chemical reaction associated to the
$allpar$ deletion or nucleation process, all the background species $i$
featuring in the reaction with negative stoichiometric coefficient (i.e.
consumed by the reaction) contribute to the activity product:

$Q_{r}=\prod_{i\in bkg\;,\;\nu_{i}<0}(\gamma_{i}c_{i})^{-\nu_{i}}$ (29)

where $c_{i}$ is the concentration in solution of the molecular species,
$\nu_{i}$ its stoichiometric coefficient in the reaction, and $\gamma_{i}$ its
activity coefficient. If a molecule is labelled as $solvent$ in the
$sol\\_start$ command, the molecule is not included in the calculation of
$Q_{r}$.

MASKE computes $\gamma_{i}$ using the WATEQ Debye-Hückel formula discussed in
the main body of this article, in Eq. 7. If the user does not know the value
of $a_{i}^{o}$ required there, by inputting a negative value in its place
MASKE will use the Davies equation instead:

$log_{10}(\gamma_{i})=-z_{i}^{2}A\left(\frac{\sqrt{I}}{1+\sqrt{I}}-0.3I\right)$
(30)

For uncharged species, with $z_{i}=0$, Eq. 7 reduces to Setchenow equation
$log_{10}(\gamma_{i})=b_{i}I$ (Langmuir, 1997) ; if the user does not know the
value of $b_{i}$, inputting a negative value in the $ChemDB$ file signals
MASKE to employ $b_{i}=0.1$ by default.

The contribution of mechanical interactions to the reaction rates in Eqs. 3
and 4 come from their excess enthalpy terms, and specifically by the changes
in total interaction energy $\Delta U_{d}$ and $\Delta U_{pr}$ caused by the
dissolution or nucleation of the particle being sampled. Fig. 14 shows how,
for pairwise interactions, these energy changes can be split into two
contributions: a bond term $\varepsilon_{0}$ which sets the local minimum
between two neighbouring particles at their equilibrium distance $r_{0}$, and
a strain energy term $U_{s}$ that emerges when the local minimum resulting
from multiple interactions leaves the sampled particle away from its
equilibrium distance $r_{0}$ from its neighbours: see Fig. 14c. For the
example in Fig. 14, therefore, if the blue particle is being sampled for
deletion, its $\Delta U_{d}$ would be $2|\varepsilon_{0}|-U_{s}$ for the
scenario in Fig. 14c, whereas for nucleation $\Delta
U_{pr}=-2|\varepsilon_{0}|+U_{s}$. As expected, a strong bond energy term
(viz. a large, negative $\varepsilon_{0}$) leads to a large negative $\Delta
U_{pr}$, which increases the precipitation rate in Eq. 4, and to a large
positive $\Delta U_{d}$, which decreases the dissolution rate in Eq. 4. The
opposite effect on rates comes from a large strain energy $U_{s}$, due to the
local particle arrangement, which leads to a less negative $\Delta U_{pr}$ and
to a less positive $\Delta U_{d}$.

Figure 14: Changes in interaction energy accompanying the deletion or
nucleation of a new particle, in blue, within a 1D array of other solid
particles. (a) A typical pairwise potential of interaction between solid
particles, here harmonic with cutoff. In the 1D array in (b) and (c), the
change in interaction energy accompanying the deletion or nucleation of the
blue particle comes from two contributions: $U_{L}$ for the interaction with
the black particle on the left, and $U_{R}$ with the black particle on the
right. (b) Energy change from particle deletion or nucleation when the sampled
particle perfectly fits the 1D array, and (c) when it does not fit the array
perfectly hence it experiences a strain energy $U_{s}$.

Fig. 14 assumes that the sampled particle, in blue, sits in a local energy
minimum resulting from the interactions with all its neighbours. In reality,
if small nanoparticles are considered, thermal fluctuations will make the
particle oscillate around its local energy minimum. For example, Fig. 15 shows
the energy landscape for particle $j$ depending on its location $r$ within a
simple system consisting of a 1D array of other particles and going from
$r=-\infty$ to $r=r_{A}$. Various local minima exist for particle $j$, the
rightmost being the most energetically favourable in this case. This rightmost
basin extends from $r_{A}$ to $r_{B}=r_{A}+r_{c}$, where $r_{c}$ is the
interaction cutoff distance from Fig. 14a. To account for thermal oscillations
within a local energy basin, the rates in Eqs. 3 and 4 should feature
interaction energies that are averaged over such oscillations (i.e. over all
the possible configurations where the particle sits within the local energy
basin, weighted by the corresponding Boltzmann factors).

Figure 15: Energy change associated to the deletion or nucleation of particle
$j$, in blue, depending on its position $r$ in the semi-infinite 1D array
(i.e. with no other particles on the right).

The $allpar$ mechanism in MASKE uses the simplifying assumption that the local
energy minimum dominates the average interaction, which is reasonable for
strong bond energies $\varepsilon_{0}$. This, however, means that the
particles sampled for deletion or nucleation should occupy such local energy
minima before computing the rates in Eqs. 3 and 4. When particle deletion is
sampled, this requirement is not automatically satisfied by MASKE, but the
user can easily ensure it by specifying a process of type $Every$ that
performs an energy minimization in LAMMPS every time a continuous or discrete
process is carried out. When instead nucleation is sampled, MASKE
automatically searches for local energy minima by creating $M$ trial particles
on a user-specified lattice and by running an energy minimization for them.
This energy minimization requires each trial particle to individually find a
local minimum: none of the current minimization styles in LAMMPS allows for
this convergence condition to be imposed, hence an additional, bespoke LAMMPS
minimization style called $quickmaske$ is included in the MASKE distribution
(it uses the same algorithm as the $quickmin$ style in LAMMPS ). During this
minimization, all the particles in MASKE are fixed to their positions except
the $M$ trial particles being sampled for nucleation. Furthermore, the only
active interactions are those between the $M$ sampled trial particles and the
other particles that already exist in MASKE, with $real\\_types$ in the input
script; all the interaction between trial particles are instead turned off. As
a result, the $M$ trial particles effectively sample the energy landscape of a
new particle nucleating at any generic position within a user-specified
region, which may well be the entire simulation box. The resolution of this
sampling depends on how small the cell of the user-specified lattice is;
$\Delta V$ in the precipitation rate in Eq. 4 is the volume of such lattice
cell. For the example in Fig. 15, assume that a 1D lattice with $M$ positions
defining small cells of size $\Delta r$ is used to sample possible nucleation
sites between $r=0$ and $r_{B}$. A number $m$ of such positions will fall
within the $A-B$ basin, i.e. between $r_{A}$ and $r_{B}$. When running energy
minimization on the trial particles, these $m$ particles will move to the
local energy minimum, where the blue particle $j$ sits in Fig. 15. As a
result, the nucleation of a particle in that energy minimum carries a weight
equal to $m\Delta r$ or, for a more general 3D case, $m\Delta V=V_{t}$. This
tributary volume $V_{t}$ is the volume of the energy basin, viz. the integral
of all the locations $r$ pertaining to the basin (if sampled through
infinitely small $\Delta V$). For the 1D case in Fig. 15, $V_{t}$ would be
simply $r_{B}-r_{A}$.

The tributary volume of the local energy basin being sampled, $V_{t}$, appears
in both the dissolution and the precipitation rates in Eqs. 3 and 4. In the
dissolution rate, $V_{t}$ is raised to the power of $\alpha/3$, where $\alpha$
is the spatial dimensionality of $c^{o{\ddagger}}_{d}$, which in turn matches
the spatial dimensionality of the rate constants measured during the
experiments that are used to obtain $\Delta G^{o{\ddagger}}_{d}$. The
$\alpha/3$ exponent ensures that the event-specific rate $r_{1,d}$ is overall
in units of number of events per unit time. Dissolution and precipitation rate
constants are typically expressed per unit surface, in which case $\alpha=2$.
The precipitation rate in Eq. 4 features a $\Delta V\cdot V_{t}^{\alpha/3-1}$
term instead. However, when all the $m$ trial particles pertaining to the same
local energy basin are energy-minimised to the same local minimum, the
cumulative rate of all these identical nucleation events (and thus the
probability of selecting that site for nucleation, from Eq. 22) ends up
scaling as $m\Delta V\cdot V_{t}^{\alpha/3-1}=V_{t}\cdot
V_{t}^{\alpha/3-1}=V_{t}^{\alpha/3}$, which brings back the same factor as in
the dissolution rate.

There is no easy way to compute $V_{t}$, which usually changes from one local
basin to another. For example, the energetically favourable, rightmost basin
in Fig. 15 has $V_{t}=r_{c}$ (in 1D), the less energetically favourable basins
to its left have $V_{t}=r_{0}<r_{c}$, whereas the basins associated to the
blue particle in Fig. 14 have $V_{t}=2r_{0}$ and $V_{t}=2(r_{0}+\delta r)$
respectively in the stress-free system in Fig. 14b and in the stressed system
in in Fig. 14c. Amorphous systems in 3D may also feature rugged energy
landscapes, with narrow energy basins with $V_{t}<r_{0}$. The $allpar$
mechanism in MASKE works under the approximation that one same representative
value of $V_{t}$ can be attributed to all basins; the user provides this value
through scaling parameters attached to each $molecule$ defined in the $ChemDB$
file. As a result, for a particle whose deletion or nucleation features a
given set of molecules in the foreground of its corresponding chemical
reaction, $V_{t}$ is obtained as the sum of the scaled volumes of all these
foreground molecules in the particle. For example, for the rightmost basin in
Fig. 15, the user may provide a scaling factor of $r_{c}/r_{0}$ to the
molecular volumes, since $r_{0}$ determines the volume that a particle in the
array occupies, and $r_{c}$ is the size of the basin. In a typical
dissolution/precipitation process, if the interaction potentials are single-
well and relatively short-range, such as in Fig. 14a, the most important
events for the overall kinetics are particle deletions and nucleations on
existing crystal surfaces, for which the scenario of the rightmost basin in
Fig. 15, and thus the $r_{c}/r_{0}$ correction, are still relevant. This is
because the volume that a particle occupies in a 3D crystal structure is still
controlled by $r_{0}$, except that the interaction potential extends to
$r_{c}$ in the direction perpendicular to the surface.

For the $\Delta V\cdot V_{t}^{\alpha/3-1}$ term in Eq. 4 to return the desired
$V_{t}^{\alpha/3}$ weight factor after summing it over trial particles, the
user must employ a sufficiently fine lattice, with small $\Delta V$, which may
require a large number of trial particles $M$ to sample the nucleation region.
For crystalline systems, the size of typical energy basins is in the range of
a particle volume (e.g. in Fig. 15) and the results in the main body of this
article have shown that $\Delta V\approx V_{p}/8$ already provide a
satisfactory sampling of nucleation sites.

For the excess enthalpy in the $allpar$ mechanism, Eq. 16 in the the main body
of this articles has shown that a relationship must hold between the
interaction energy of two bonded particles in stress-free conditions,
$\varepsilon_{0}$, and the solid-solution interfacial energy of a particle
$\gamma\Omega$. Namely:

$\varepsilon_{0}=\frac{2\gamma\Omega}{n_{b}}$ (31)

where $n_{b}$ is the number of neighbours for a bulk particle. The user can
manipulate the excess enthalpy term through their input of $\gamma$ in the
$ChemDB$ file, and through their choice and parametrization of the interaction
potentials in LAMMPS (via the MASKE input script). This can lead to a wide
range of different system behaviours. For example, by setting $\gamma=0$ and
an interaction potential where the bond energy $\varepsilon_{0}$ is set to
zero as well, the user would create a simulation where particles would have
the same probability of being deleted or nucleated irrespective of whether
they sit next to other existing particles, or they are isolated in the
solution, or even they completely overlap with other existing particles. This
scenario is unphysical, hence the user must be careful in setting a meaningful
relationship between $\gamma$ and interaction potentials when preparing a
simulation. Eq. 31 is valid only for particles interacting through single-well
potentials and that can form an orderly lattice. Modifications of the equation
or of the definition of excess enthalpy altogether are needed, for example
when dealing with amorphous agglomerates of solid particles that feature
eigenstresses Alex et al. (2023), or with larger particles that cannot be
properly treated with the $allpar$ mechanism.

The last term discussed in relation to the the rate equations, Eqs. 3 and 4,
is the fraction $\chi$, of excess enthalpy affecting the precipitation rate.
If $\chi$ is set to zero, the excess enthalpy from solid-solid and solid-
solution interactions only affects the dissolution rate. This is the usual
assumption in Transition State Theory, where one assumes that the excess
enthalpy raises the free energy level of the solid phases involved in the
reaction and it is already fully lost in the activated complex state. As a
result, with $\chi=0$ the excess enthalpy has no effect on the free energy
barrier for precipitation, which involves solvated molecules as reactants and
the activated complex. As noted in Shvab et al. (2017), however, $\chi=0$
entails that particle nucleation becomes position independent (under the
current solution model in MASKE, where ion concentrations and temperature are
assumed uniform everywhere in a simulation box). This leads to the
computationally inefficient scenario whereby many nucleation events may be
carried out in energetically very unfavourable locations, e.g. overlapping
with other existing particles, just for the opposite deletion event to be
carried out immediately afterwards due to the high excess enthalpy of the
recently nucleated particle. Setting $\chi>0$ coincides with assuming that a
fraction of the excess enthalpy is still present in the activated complex. If
a dissolution event is being sampled and the excess enthalpy change is
positive, i.e. the event is enthalpy-wise unfavourable (e.g. dissolving a
well-connected particle in the bulk of a crystal), using $\chi>0$ reduces the
activation energy for dissolution, thus increasing the dissolution rate. On
the other hand, the opposite event of nucleating the same particle at the same
site will also feature a lower activation energy compared to the $\chi=0$
case, because the activated complex now has a lower free energy due to its
portion of favourable, i.e. negative, free energy on the way towards the solid
state. Hence also the nucleation rate back to the same position will increase.
The increase of both dissolution and precipitation rates produces a
compensation effect on the probability of a particle actually existing in that
site, with more likelihood of seeing a particle dissolving and reforming there
compared to the $\chi=0$ case. The case of a deletion event with negative
change in enthalpy instead, i.e. an event that is enthalpy-wise favourable
(e.g. dissolving an particle under stress), will lead to a higher activation
energy for dissolution and for precipitation back in the same site. As a
result both dissolution and nucleation rates will decrease, once again
producing a compensation effect on the probability of a particle existing in
that site, but with less likelihood of seeing a particle dissolving and
reforming there compared to the $\chi=0$ case. This mitigates the
aforementioned issue of nucleating particles in sites that are so
energetically unfavourable that a subsequent deletion event would immediately
remove them.

##  Appendix G Net rates equations in $allpar$, and recovering Transition
State Theory

In the main body of this articles, the expressions of reaction rates for
dissolution and precipitation, Eqs. 3 and 4 (as per $allpar$ mechanism in
MASTKE), have been used to also derive net reaction rates, in Eqs. 9 and 10. A
few considerations underpin the derivation of the net rate expressions.
Starting with the net dissolution rate in Eq. 9, a first consideration is that
in the backward reaction, viz. straight precipitation, all the trial particles
sampling precipitation into the same energy basin will have identical straight
rate $r_{1,pr}$. This leads to the $V_{t}^{\alpha/3}$ term in Eq. 9, which
comes from summing all the $\Delta V_{t}^{\alpha/3-1}$ terms from $r_{1,pr}$
pertaining to the same energy basin. Another consideration is the assumption
that the activated state is identical for the forward (here, dissolution) and
backward (here, precipitation) reactions, hence
$c^{o{\ddagger}}_{d}=c^{o{\ddagger}}_{pr}$ and
$\gamma^{\ddagger}_{d}=\gamma^{\ddagger}_{pr}$. Since precipitation in this
case is simply the reverse reaction of dissolution, $n_{r,j,diss}=n_{r,j,nuc}$
and $Q_{r,pr}$, which is the activity product of the reactants in the
precipitation reaction, has been re-labelled as $Q_{p,d}$, i.e. the activity
product of the products of the dissolution reaction. Lastly, Eq. 9 features
the equilibrium constant of the dissolution reaction,
$K_{eq,d}=\left(\frac{Q_{p,d}}{Q_{r,d}}\right)_{eq}$, where the $eq$ subscript
indicates that the activity products are measured at equilibrium. $K_{eq,d}$
appears in the rate equation because of its relationship to the standard
change in free energy for the dissolution reaction, $\Delta G^{o}_{d}$, i.e.
$K_{eq,d}=\exp\left(-\frac{\Delta G^{o}_{d}}{k_{B}T}\right)$. This $\Delta
G^{o}_{d}$ is the difference in standard free energy between products and
reactants in the dissolution reaction, hence it determines the difference
between forward and backward standard activation energies: $\Delta
G^{o{\ddagger}}_{pr}+\Delta G^{o}_{d}=\Delta G^{o{\ddagger}}_{d}$.

All the terms in the net dissolution rate expression, in Eq. 9, ultimately
relate only to quantities that are defined by the user when inputting the
parameters in the $ChemDB$ file, which that inform the chemical reactions
associated to the dissolution itself. In evaluating the net dissolution rate,
MASKE does not explicit sampling any nucleation event. Therefore, a simulation
may sample exclusively particle deletion events, knowing that Eq. 9 will
reduce the dissolution rate by the rate of the corresponding but reverse
precipitation reactions. A limitation of Eq. 9 is that, if the contribution
from the precipitation rate is greater than from the dissolution rate, then
$r_{1,d}^{net}$ will be zero and the dissolution event will never take place.
This excludes the possibility of any energetically unfavourable fluctuation,
which limits the ability to simulate some mechanisms such as crystal
nucleation. This limitation has been addressed in the main body of this
article.

Similar considerations as above apply to the derivation of the net
precipitation rate, in Eq. 10. There the backward reaction, which is a
straight dissolution, refers to a corresponding precipitation event that
covers only a portion $\Delta V$ of an energy basin with volume $V_{t}$. Hence
$r_{1,d}$ carries the same $\Delta V\cdot V_{t}^{\alpha/3-1}$ as the
precipitation event which is actually being sampled. Indeed, all term in Eq.
10 ultimately only refer to a precipitation reaction being sampled, similar to
how all terms in Eq. 9 only referred to the dissolution reaction. $Q_{p,pr}$
is the activity product of the products of the precipitation reaction, which
replaces the activity product of the reactants of the dissolution reaction
that, here, is just the reverse of the precipitation one. The equilibrium
constant of the precipitation reaction,
$K_{eq,pr}=\left(\frac{Q_{p,pr}}{Q_{r,p}}\right)_{eq}$ emerges in a similar
way as $K_{eq,d}$ in Eq. 9. Namely, $\Delta G^{o}_{pr}$ is the standard change
in free energy for the precipitation reaction,
$K_{eq,pr}=\exp\left(-\frac{\Delta G^{o}_{pr}}{k_{B}T}\right)$ and $\Delta
G^{o{\ddagger}}_{d}+\Delta G^{o}_{pr}=\Delta G^{o{\ddagger}}_{pr}$.

The net rate equations developed in this article are extensions of Transition
State Theory (TST). TST ultimately leads to simple equations of net rate of
the form Lasaga (2014):

$\mathcal{R}_{1,d}^{TST}=k(1-\beta)$ (32)

$\mathcal{R}$ is the net rate of reaction per unit area or per unit volume,
$k$ is the rate constant for the dissolution reaction, and
$\beta=\frac{Q_{p,d}/Q_{r,d}}{K_{eq,d}}$ is the saturation index of the
dissolution reaction. For solid dissolution, since a stress-free solid is
commonly assumed to be in its standard state, $Q_{r,d}=1$ and thus $\beta$
simplifies to $\frac{Q_{p,d}}{K_{eq,d}}$ or, that is the same if one considers
the precipitation reaction opposite to the dissolution one,
$\beta=\frac{Q_{r,pr}}{K_{eq,d}}$. The TST rate equation for the opposite
precipitation reaction is obtained simply by detailed balance:

$\mathcal{R}_{1,pr}^{TST}=k(\beta-1)$ (33)

where $k$ and $\beta$ are still those for the dissolution reaction. Under a
set of circumstances, the TST rates in Eqs. 32 and 33 are recovered from the
$allpar$ rates in MASKE, viz. from Eqs. 9 and 10. First of all, take $\chi=0$,
meaning that any excess enthalpy from solid-solid and solid-solution
interactions only exists in the solid and is fully lost already in the
activated complex state. Let us then define:

$k=\kappa\frac{k_{B}T}{h}\frac{c^{o{\ddagger}}_{d}}{\gamma^{\ddagger}_{d}}\exp\left(\frac{-\Delta
G^{o{\ddagger}}_{d}}{k_{B}T}\right)$ (34)

which has units of number of events per unit time, per unit area or volume
(depending on the dimensionality $\alpha$ of $c^{\ddagger}$). With this one
can already rewrite the non-zero case in Eq. 9 as:

$\mathcal{R}_{1,d}=\frac{r_{1,d}^{net}}{V_{t}^{\alpha/3}}=k\left\\{\exp\left[-\frac{\Delta
U_{d}-\gamma\Omega}{k_{B}T}\cdot\frac{1}{n_{r,j,diss}}\right]Q_{r,d}-\frac{Q_{p,d}}{K_{eq,d}}\right\\}$
(35)

If the physically meaningful relationship between bond strength and surface
energy in Eq. 16 is imposed, the excess enthalpy $\Delta U_{d}-\gamma\Omega$
in Eq. 35 is zero for a stress-free particle occupying a kink position in its
crystal lattice. Noting then that $Q_{r,d}=1$ for a dissolving solid and that
$\frac{Q_{p,d}}{K_{eq,d}}$ is indeed the definition of $\beta$, one recovers
indeed the TST net rate equation in Eq. 32. This means that the rate equation
in MASKE is a more general form of Eq. 32, which reduces to the usual TST form
for dissolution processes that are controlled by stress-free kink particles,
if one assumes that excess enthalpy can only be present in the non-activated
solid state. At the macroscopic level reaction rates are often expressed
following Eq. 32; despite tempting, it would be incorrect to define a
macroscopic rate constant $k$ as per Eq. 34. Indeed, a macroscopic rate
constant implicitly accounts for morphological effect, such as for example the
number of kink sites per unit area in a kink-controlled dissolution process,
or for complex series of reactions all subsumed into a single, effective one
(Cefis and Comi, 2017). Eq. 34 instead applies only to a one-step reaction and
does not subsume any morphological effect, as in the $allpar$ mechanism in
MASKE they are accounted for through the excess enthalpy term in Eq. 9.

The TST precipitation rate in Eq. 33 is also a special case of MASKE’s rate
from Eq. 10. To show this, let us first recall the common assumptions that
$c^{o{\ddagger}}_{pr}=c^{o{\ddagger}}_{pr}$ and
$\gamma^{\ddagger}_{pr}={\ddagger}_{d}$, and the following relationships:
$\Delta G^{o{\ddagger}}_{pr}+\Delta G^{o}_{d}=\Delta G^{o{\ddagger}}_{d}$ and
$K_{eq,d}=\exp\left(-\frac{\Delta G^{o}_{d}}{k_{B}T}\right)$. These lets us
reformulate the prefactor in Eq. 10 (which one can identify as a precipitation
rate constant $k_{pr}$) as:

$k_{pr}=\kappa\frac{k_{B}T}{h}\frac{c^{o{\ddagger}}_{pr}}{\gamma^{\ddagger}_{pr}}\exp\left(\frac{-\Delta
G^{o{\ddagger}}_{pr}}{k_{B}T}\right)=k\exp\left(\frac{-\Delta
G^{o}_{d}}{k_{B}T}\right)=\frac{k}{K_{eq,d}}$ (36)

where $k$ is the dissolution rate constant in Eq. 34. In fact, Eq. 36 simply
states the well-known equilibrium condition whereby
$\frac{k}{k_{pr}}=K_{eq,d}$. If one assumes again $\chi=0$ and stress-free
kink particles with no excess enthalpy (from Eq. 16), Eq. 10 becomes:

$\mathcal{R}_{1,pr}=\frac{r_{1,pr}^{net}}{\Delta V\cdot
V_{t}^{\alpha/3-1}}=\frac{k}{K_{eq,d}}\left(Q_{r,pr}-\frac{Q_{p,pr}}{K_{eq,pr}}\right)$
(37)

Noting finally that $\frac{Q_{r,pr}}{K_{eq,d}}$ is the saturation index
$\beta$ for the dissolution reaction, that $Q_{p,pr}=1$ when the products are
solids, and that $K_{eq,pr}=K_{eq,d}^{-1}$, one recovers the usual TST rate in
Eq. 33.