To show that this event hold with vanishing probability we first write
$\sum_{i\in S_{+}}z_{i}\omega_{i}=\xi^{T}Ae_{s}$ where $e_{s}$ is the vector
whose $j^{th}$ entry is $\omega_{j}\mathds{1}\\{j\in S_{+}\\}$. We observe
that $(Ae_{s})_{i}$ is subGaussian with subGaussian factor $Cc^{2}s/n$ where
$s=|S_{+}|$, since $A=X/\sqrt{n}$ and $X$ has isotropic subGaussian rows.
Hence conditioned on $\xi$ the random variable $\xi^{T}Ae_{s}$ has subGaussian
factor $Cc^{2}\|\xi\|^{2}s/n$ since $\xi$ is independent of $A$. We define the
event $F:=\\{\|\xi\|\leq C\sqrt{n}\\}$. Since $\|\xi\|\leq C\sqrt{n}$ with
probability $1-n^{-c}$, for a given subset $S$ of size s we have

$P(\\{c^{2}\lambda^{2}s/2\leq 2\lambda|\sum_{i\in
S_{+}}z_{i}\omega_{i}|\\}\cap F)\leq\exp(-c^{4}\lambda^{2}(s/4)^{2}/2Cs)$

where we have simplified the constants. Finally we obtain,

$\displaystyle P(\\{|S_{+}|\geq c_{0}n\\}\cap E\cap F)$ $\displaystyle\leq
P(\\{\lambda c^{2}|S_{+}|/4\leq|\sum_{i\in S_{+}}z_{i}\omega_{i}|\\}\cap F).$
$\displaystyle\leq\sum_{s=c_{0}n}^{p}\sum_{S:|S|=s}\exp(-c^{4}\lambda^{2}(s/4)^{2}/2Cs)$
$\displaystyle\leq C\exp(-Cc_{0}\lambda^{2}n),$

where we have simplified all the constants except $c_{0}$. We can remove the
conditioning on $E$ and $F$ since both holds with probability at least
$1-n^{-c}$.

Similar bounds hold for $S_{-}$ and we have shown (redefining $c_{0}$) that
for $\lambda\geq\lambda_{0}$, the Lasso solution is $c_{0}n$ sparse with
probability at least $1-n^{-c}$. ∎

###### Proof of Lemma A.3.

We start with the Fenchel dual for the function $h$,

$h^{*}(x)=\sup_{z}\langle z,x\rangle-\frac{1}{2}\|A_{S}z\|^{2}-\langle
A_{S}z,y\rangle.$

It can be easily verified that the maximum is attained at
$z=(A_{S}^{T}A_{S})^{-1}(x-A_{S}^{T}y)$ and plugging in the value we obtain

$h^{*}(x)=\frac{1}{2}(x-A^{T}_{S}y)^{T}(A^{T}_{S}A_{S})^{-1}(x-A^{T}_{S}y).$

To compute $b^{*}(x)$ we start with defining the function
$\bar{b}(x)=\lambda\omega|x+y|$ for scalars $x$ and $y$. Thus we have,

$\bar{b}^{*}(x)=\sup_{z}xz-\lambda\omega|z+y|.$

It can be easily observed that $\bar{b}^{*}(x)=\infty$ unless
$|x|\leq\lambda\omega$.

Further define the function

$\alpha(x,z)=xz-\lambda\omega|z+y|.$

If $z\geq-y$, the function $\alpha(x,z)=z(x-\lambda\omega)-\lambda y$ attains
maximum value at $z=-y$ as $x-\lambda\omega\leq 0$. Similarly if $z\leq-y$,
$\alpha(x,z)=z(x+\lambda\omega)+\lambda y$ attains maximum value at $z=-y$ as
$x+\lambda\omega\geq 0$. Hence

$\bar{b}^{*}(x)=\begin{cases}\sup_{z}\alpha(x,z)=-xy,\quad\text{ if
}|x|\leq\lambda\omega\\\ \infty,\text{ otherwise }.\end{cases}$

Now we consider the function

$b^{*}(x)=\sup_{z}\langle
z,x\rangle-\lambda\|\tilde{\Omega}_{S}(z+\mu_{0,S})\|_{1}$

and note that the optimization problem in the RHS is separable in its
arguments. Hence we have

$b^{*}(x)=\begin{cases}-\langle\mu_{0,S},x\rangle,\quad\text{ if
}|x_{i}/\omega_{i}|\leq\lambda\text{ for all }i\in S\\\ \infty,\text{
otherwise },\end{cases}$

which was to be shown. ∎

###  A.3 Universality of the Optimizer

###### Proof of Lemma 7.3.

We only prove the result for the ridge case, since the Lasso case follows by
suitably modifying (to account for the term $\tilde{\Omega}$) the proof of
Theorem 3.8 of [HS23] together with the statements from [MM21]. Recall that
$G$ is the design matrix with i.i.d. Gaussian entries. Denoting
$\bar{G}=G/\sqrt{n}$ problem, we consider the following functions:

$\displaystyle r(w,u)$
$\displaystyle=\frac{1}{n}u^{T}\bar{G}w-\frac{1}{n}u^{T}\xi-\frac{1}{2n}\|u\|^{2}+\frac{\lambda}{2n}(\|\tilde{\Omega}(w+\mu_{0})\|^{2}-\|\tilde{\Omega}\mu_{0}\|^{2}),$
$\displaystyle\ell(w,u)$
$\displaystyle=-\frac{1}{n^{-3/2}}\|u\|g^{T}w+\frac{\|w\|h^{T}u}{n^{3/2}}-\frac{1}{n}u^{T}\xi-\frac{1}{2n}\|u\|^{2}+\frac{\lambda}{2n}(\|\tilde{\Omega}(w+\mu_{0})\|^{2}-\|\tilde{\Omega}\mu_{0}\|^{2}),$

where $u\in\mathbb{R}^{n}$ and $h$ and $g$ are independent random vectors with
distributions $h\sim N(0,I_{n})$, $g\sim N(0,I_{p})$ respectively. We denote
by $R(w)=\max_{u}h(w,u)$ and $L(w)=\max_{u}\ell(w,u)$ the ridge cost and the
associated Gordon cost respectively.

Note that using $w=\sqrt{n}\theta$, $R(w)=R^{K}(\theta,G)$ for $K=R$ and we
work with the notation $R(w)$ since the Gaussian design and the ridge setup is
clear from the context. Our goal is to demonstrate that the optimal value of
$R(w)$ i.e. $\min_{w}R(w)$ is close to $\psi^{R}$. We will use CGMT to show
that $\min_{w}R(w)$ is well approximated by $\min_{w}L(w)$ and the latter is
close to the desired constant $\psi^{R}$ (a constant that we will define
shortly). The proof proceeds by showing that $\min_{w}L(w)$ is well
approximated by the value at saddle point of a function
$\psi_{n}(\beta,\gamma)$, which is further approximated by the value at saddle
point of its non-stochastic version $\psi(\beta,\gamma)$. Denoting the saddle
point of the latter by $(\beta_{*},\gamma_{*})$, $\psi^{R}$ will be defined as
$\psi(\beta_{*},\gamma_{*})$. We start by defining $\psi_{n}$ and $\psi$.

Define the function $\psi_{n}(\beta,\gamma)$ as follows:

$\psi_{n}(\beta,\gamma)=\left(\frac{\sigma^{2}}{\gamma}+\gamma\right)\frac{\beta}{2}-\frac{\beta^{2}}{2}-\frac{1}{n/p}\frac{1}{p}\sum_{i=1}^{p}\frac{\gamma\lambda^{2}\omega^{4}_{i}}{2(\beta+\lambda\omega^{2}_{i}\gamma)}\left(\mu_{0,i}-\frac{\beta
g_{i}}{\lambda\omega^{2}_{i}}\right)^{2}.$ (A.9)

The above equation can be written compactly as

$\psi_{n}(\beta,\gamma)=\left(\frac{\sigma^{2}}{\gamma}+\gamma\right)\frac{\beta}{2}-\frac{\beta^{2}}{2}-\frac{1}{n/p}E_{P_{n}}\left[\frac{\gamma\lambda^{2}\Omega^{4}}{2(\beta+\lambda\Omega^{2}\gamma)}\left(M-\frac{\beta\Gamma}{\lambda\Omega^{2}}\right)^{2}\right],$
(A.10)

where the expectation is taken with respect to the law $P_{n}$ of the
$3$-tuple of random variables
$(M,\Omega,\Gamma)\sim\frac{1}{p}\sum_{i=1}^{p}\delta_{\mu_{i},\omega_{i},g_{i}}$.

Also define the function $\psi(\beta,\gamma)$ as follows:

$\psi(\beta,\gamma)=\left(\frac{\sigma^{2}}{\gamma}+\gamma\right)\frac{\beta}{2}-\frac{\beta^{2}}{2}-\frac{1}{n/p}E_{Q_{n}}\left[\frac{\gamma\lambda^{2}\Omega^{4}}{2(\beta+\lambda\Omega^{2}\gamma)}\left(M-\frac{\beta\Gamma}{\lambda\Omega^{2}}\right)^{2}\right],$
(A.11)

where the expectation is taken with respect to the law $Q_{n}$ of the
$3$-tuple of random variables
$(M,\Omega,\Gamma)\sim(\frac{1}{p}\sum_{i=1}^{p}\delta_{\mu_{0,i},\omega_{i}})\otimes
N(0,1)$. Taking the derivates of the function $\psi_{n}$ with respect to
$\beta$ and $\gamma$ one obtains the following fixed point equations for the
saddle points:

$\displaystyle\gamma^{2}_{*,n}$
$\displaystyle=\sigma^{2}+\frac{1}{n/p}E_{P_{n}}\left[\eta_{2}\left(M+\gamma_{*,n}\Gamma;\frac{\gamma_{*,n}\lambda\Omega^{2}}{\beta_{*,n}}\right)-M\right]^{2},$
(A.12) $\displaystyle\beta_{*,n}$
$\displaystyle=\gamma_{*,n}-\frac{1}{n/p}E_{P_{n}}\left[\Gamma.\left(\eta_{2}\left(M+\gamma_{*,n}\Gamma;\frac{\gamma_{*,n}\lambda\Omega^{2}}{\beta_{*,n}}\right)-M\right)\right].$
(A.13)

Similarly from equation (A.11) we get the fixed point equations for the
limiting case as follows:

$\displaystyle\gamma^{2}_{*}$
$\displaystyle=\sigma^{2}+\frac{1}{n/p}E_{Q_{n}}\left[\eta_{2}\left(M+\gamma_{*}\Gamma;\frac{\gamma_{*}\lambda\Omega^{2}}{\beta_{*}}\right)-M\right]^{2},$
(A.14) $\displaystyle\beta_{*}$
$\displaystyle=\gamma_{*}-\frac{1}{n/p}E_{Q_{n}}\left[\Gamma.\left(\eta_{2}\left(M+\gamma_{*}\Gamma;\frac{\gamma_{*}\lambda\Omega^{2}}{\beta_{*}}\right)-M\right)\right].$
(A.15)

It can be shown following the steps of Proposition 5.2 of [HS23] that the
solutions $\beta_{*,n}$ and $\gamma_{*,n}$ lie in some compact subset of
$\mathbb{R}^{+}$ bounded away from $0$ with probability at least $1-n^{-c}$
and the same can be shown for $\beta_{*}$ and $\gamma_{*})$. Next we show that
$|\beta_{*,n}-\beta_{*}|$ and $|\gamma_{*,n}-\gamma_{*}|$ converges to $0$
with probability at least $1-n^{-c}$. However because of the matrix $\Omega$
this convergence does not follow easily from [HS23]. We use different
techniques and prove the statement in a separate lemma whose proof is
technical and hence we defer for the moment

######  Lemma A.4.

With probability at least $1-n^{-c}$, we have

$|\beta_{*,n}-\beta_{*}|=O(\epsilon_{n})\quad|\gamma_{*,n}-\gamma_{*}|=O(\epsilon_{n}).$

for some $\epsilon_{n}\rightarrow 0$.

The rest of the proof follows in the lines of [HS23]. Since the modification
involving the matrix $\Omega$ is routine we only sketch the outline and omit
detailed proof. One shows (see [HS23, Lemma 5.2]) that with probability at
least $1-n^{-c}$ the following statements hold:

  1. 1.

$|\psi_{n}(\beta_{*,n},\gamma_{*,n})-\psi(\beta_{*},\gamma_{*})|=O(\epsilon_{n})$,

  2. 2.

$|\min_{w}L(w)-\max_{\beta>0}\min_{\gamma>0}\psi_{n}(\beta,\gamma)|=O(\epsilon_{n})$,

  3. 3.

$|\min_{w}H(w)-\psi(\beta_{*},\gamma_{*})|=O(\epsilon_{n}).$

Now for any Lipshitz function $\varrho$ define the following set

$S_{n,\epsilon}(\varrho)=\\{w\in\mathbb{R}^{p}:|\varrho(w/\sqrt{n})-E\varrho(w_{*}/\sqrt{n})|\geq\epsilon^{-1/2}\\}.$

Then it can be shown (see [HS23, Theorem 3.4] and its proof) that the
following statement holds:

$\displaystyle
P\left(\min_{w}R(w)\geq\psi(\beta_{*},\gamma_{*})+\epsilon\right)$
$\displaystyle\leq Cn^{-c},$ $\displaystyle P\left(\min_{w\in
S_{n,\epsilon}(\varrho)}R(w)\leq\psi(\beta_{*},\gamma_{*})+K^{-1}\epsilon\right)$
$\displaystyle\leq Cn^{-c}$

for some $K>0$. Recall that $\theta=w/\sqrt{n}$,
$\psi^{R}=\psi(\beta_{*},\gamma_{*})$. Then letting $\delta_{n}=Cn^{-c}$ and
choosing $z$ suitably the above inequalities imply that the statement of the
Lemma 7.3 is true if we can prove Lemma A.4. ∎

###### Proof of lemma A.4.

First we adopt the method used in [MRSY19] to show that with probability at
least $1-n^{-c}$

$\|(\beta_{*,n},\gamma_{*,n})-(\beta^{\dagger},\gamma^{\dagger})\|=O(\epsilon_{n})$
(A.16)

holds for some $\epsilon_{n}\rightarrow 0$ and for some constants
$\beta^{\dagger},\gamma^{\dagger}$. Let us denote the joint law of the random
variables $(\bar{M},\bar{\Omega},Z)$ as $P$, where $(\bar{M},\bar{\Omega})$
has law $\mu$ (see Assumption 4) and $Z$ is independently distributed as
$N(0,1)$. Note that by Assumption 4 $P_{n}\xRightarrow{W_{2}}P$ and
$Q_{n}\xRightarrow{W_{2}}P$. For notational convenience we will define
functions of a general 3-touple of random variables $(M,\Omega,\Gamma)$ and
depending on the context they will be distributed as $P_{n},Q_{n}$ or $P$.

Step:1 We write the optimizers as solutions of some functions i.e. define

$\displaystyle F_{n}(\beta,\gamma)$
$\displaystyle=\frac{1}{n/p}\left[\eta_{2}\left(M+\gamma\Gamma;\frac{\gamma\lambda\Omega^{2}}{\beta}\right)-M\right]^{2}-\sigma^{2}-\gamma^{2},$
(A.17) $\displaystyle G_{n}(\beta,\gamma)$
$\displaystyle=\frac{1}{n/p}\left[\Gamma.\left(\eta_{2}\left(M+\gamma\Gamma;\frac{\gamma\lambda\Omega^{2}}{\beta}\right)-M\right)\right]+\beta-\gamma,$
(A.18) $\displaystyle F(\beta,\gamma)$
$\displaystyle=\kappa\left[\eta_{2}\left(M+\gamma\Gamma;\frac{\gamma\lambda\Omega^{2}}{\beta}\right)-M\right]^{2}-\sigma^{2}-\gamma^{2},$
(A.19) $\displaystyle G(\beta,\gamma)$
$\displaystyle=\kappa\left[\Gamma.\left(\eta_{2}\left(M+\gamma\Gamma;\frac{\gamma\lambda\Omega^{2}}{\beta}\right)-M\right)\right]+\beta-\gamma.$
(A.20)

Note that $(\beta_{*,n},\gamma_{*,n})$ are solutions of the equations
$E_{P_{n}}F_{n}(\beta,\gamma)$ and $E_{P_{n}}G_{n}(\beta,\gamma)$. We define
$(\beta^{\dagger},\gamma^{\dagger})$ as solutions of the equation
$E_{P}F(\beta,\gamma)$ and $E_{P}G(\beta,\gamma)$. It can be shown along the
lines of Lemma 5.2 of [HS23] that both solutions are unique.

Step 2: We will show that the following uniform convergence holds in some
compact set $\mathcal{C}$ of $\mathbb{R}^{+}\times\mathbb{R}^{+}$bounded away
from $(0,0)$:

$\displaystyle\lim_{n\rightarrow\infty}\sup_{(\gamma,\beta)\in\mathcal{C}}|E_{P_{n}}F_{n}(\beta,\gamma)-E_{P}F(\beta,\gamma)|$
$\displaystyle=0,$ (A.21)
$\displaystyle\lim_{n\rightarrow\infty}\sup_{(\gamma,\beta)\in\mathcal{C}}|E_{P_{n}}G_{n}(\beta,\gamma)-E_{P}G(\beta,\gamma)|$
$\displaystyle=0.$ (A.22)

We defer the proof for the moment.

Step 3: Since the optimizers in both the finite sample case
($(\beta_{*,n},\gamma_{*,n})$) and the limiting case
($(\beta^{\dagger},\gamma^{\dagger})$) lie in compact sets with probability at
least $1-n^{-c}$ (see Lemma 5.2 of [HS23]), we can choose $\mathcal{C}$, so
that the optimizers belong to $\mathcal{C}$ with probability at least
$1-n^{-c}$. Now let $(\tilde{\beta},\tilde{\gamma})$ be a limit point of the
sequence $(\beta_{*,n},\gamma_{*,n})$ then we prove that
$(\tilde{\beta},\tilde{\gamma})=(\beta^{\dagger},\gamma^{\dagger})$.

Indeed,

$|E_{P}(F(\tilde{\beta},\tilde{\gamma}))|\leq|E_{P}(F(\tilde{\beta},\tilde{\gamma}))-E_{P}(F(\beta_{*,n},\gamma_{*,n}))|+|E_{P_{n}}(F(\beta_{*,n},\gamma_{*,n}))-E_{P}(F(\beta_{*,n},\gamma_{*,n}))|,$

since $E_{P_{n}}(F(\beta_{*,n},\gamma_{*,n}))=0$. Now from the continuity of
$E_{P}(F(.,.))$ we have

$|E_{P}(F(\tilde{\beta},\tilde{\gamma}))-E_{P}(F(\beta_{*,n},\gamma_{*,n}))|\rightarrow
0,$

while from equation (A.21) we conclude that

$|E_{P_{n}}(F(\beta_{*,n},\gamma_{*,n}))-E_{P}(F(\beta_{*,n},\gamma_{*,n}))|\rightarrow
0.$

Thus $E_{P}(F(\tilde{\beta},\tilde{\gamma}))=0$ and similarly we can prove
$E_{P}(G(\tilde{\beta},\tilde{\gamma}))=0$. Since
$(\beta^{\dagger},\gamma^{\dagger})$ is the unique solution of the set of
equations $E_{P}(F(\beta,\gamma))=0$ and $E_{P}(G(\beta,\gamma))=0$, we obtain
$(\tilde{\beta},\tilde{\gamma})=(\beta^{\dagger},\gamma^{\dagger})$. Noting
that the optimizers $(\beta_{*n},\gamma_{*n})$ and
$(\beta^{\dagger},\gamma^{\dagger})$ lie in the compact set $\mathcal{C}$ with
probability at least $1-n^{-c}$ this proves that for some $\epsilon_{n}$, we
have
$\|(\beta_{*,n},\gamma_{*,n})-(\beta^{\dagger},\gamma^{\dagger})\|=O(\epsilon_{n})$
with probability at least $1-n^{-c}$.

Thus, to prove the lemma it is sufficient to prove (A.21) and (A.22). Since
$p/n\rightarrow\kappa$ it is enough to prove

$\lim_{n\rightarrow\infty}\sup_{(\gamma,\beta)\in\mathcal{C}}|E_{P_{n}}F(\beta,\gamma)-E_{P}F(\beta,\gamma)|=0.$

By Arzela-Ascoli theorem the proof follows if we can show the following
statements:

  1. 1.

Pointwise convergence: For a given $(\gamma,\beta)$ in $\mathcal{C}$,
$\lim_{n\rightarrow\infty}|E_{P_{n}}F(\beta,\gamma)-E_{P}F(\beta,\gamma)|=0.$

  2. 2.

Equicontinuity: $\forall\epsilon,\exists\delta$ such that if
$\|(\beta,\gamma)-(\beta^{\prime},\gamma^{\prime})\|<\delta$,

$\sup_{n}|E_{P_{n}}F(\beta,\gamma)-E_{P_{n}}F(\beta^{\prime},\gamma^{\prime})|<\epsilon.$

Writing $F(\beta,\gamma)$ as $f(M,\Gamma,\Omega)$, a function of
$M,\Gamma,\Omega$, one observes that to show (1), the following result can be
used: $P_{n}\xRightarrow{W_{2}}P$ implies
$|E_{P_{n}}F(\beta,\gamma)-E_{P}F(\beta,\gamma)|=0$, if one can show
$\sup_{(M,\Gamma,\Omega)}\frac{f(M,\Gamma,\Omega)}{\|M,\Gamma,\Omega\|^{2}}<\infty$.
Noting the form of $F(\beta,\gamma)$ it is enough to show that,

$\sup_{(M,\Gamma,\Omega)}\frac{\left[\eta_{2}\left(M+\gamma\Gamma;\frac{\gamma\lambda\Omega^{2}}{\beta}\right)-M\right]^{2}}{\|M,\Gamma,\Omega\|^{2}}\leq\sup_{(M,\Gamma,\Omega)}\frac{\left[(-\lambda\Omega^{2}M+\beta\Gamma)/(\beta\gamma^{-1}+\lambda\Omega^{2})\right]^{2}}{\|M,\Gamma\|^{2}}<\infty.$

Since $\Omega$ is bounded, the numerator is bounded by a quadratic function of
$M$ and $G$ and hence the claim is proved.

For the function $G$, we similarly have to show

$\sup_{(M,\Gamma,\Omega)}\frac{\Gamma(\eta_{2}\left(M+\gamma\Gamma;\frac{\gamma^{2}\lambda\Omega^{2}}{\beta}\right)-M)}{\|M,\Gamma,\Omega\|^{2}}\leq\sup_{(M,\Gamma,\Omega)}\frac{\Gamma(-\lambda\Omega^{2}M+\beta\Gamma)/(\beta\gamma^{-1}+\lambda\Omega^{2})}{\|M,\Gamma\|^{2}}<\infty.$

Again, the numerator is upper bounded by a quadratic function of $\Gamma$ and
$M$ is bounded; hence, the claim is proved.

To show (2) we observe that,

$\displaystyle F(\beta,\gamma)-F(\beta^{\prime},\gamma^{\prime})$
$\displaystyle=\left[(-\lambda\Omega^{2}M+\beta\Gamma)/(\beta\gamma^{-1}+\lambda\Omega^{2})\right]^{2}-\left[(-\lambda\Omega^{2}M+\beta^{\prime}\Gamma)/(\beta^{\prime}\gamma^{\prime-1}+\lambda\Omega^{2})\right]^{2}+\gamma^{{}^{\prime}2}-\gamma^{2}$
$\displaystyle=\frac{(\beta^{\prime}\gamma^{\prime-1}+\lambda\Omega^{2})^{2}(-\lambda\Omega^{2}M+\beta\Gamma)^{2}-(\beta\gamma^{-1}+\lambda\Omega^{2})^{2}(-\lambda\Omega^{2}M+\beta^{\prime}\Gamma)^{2}}{(\beta\gamma^{-1}+\lambda\Omega^{2})^{2}(\beta^{\prime}\gamma^{\prime-1}+\lambda\Omega^{2})^{2}}+\gamma^{{}^{\prime}2}-\gamma^{2}.$

Noting that $\beta,\beta^{\prime},\gamma,\gamma^{\prime}$ lies in compact sets
bounded away from $0$, we need to show that
$\|(\beta,\gamma)-(\beta^{\prime},\gamma^{\prime})\|<\delta$ implies

$|E_{P_{n}}[(\beta^{\prime}\gamma^{\prime-1}+\lambda\Omega^{2})^{2}(-\lambda\Omega^{2}M+\beta\Gamma)^{2}-(\beta\gamma^{-1}+\lambda\Omega^{2})^{2}(-\lambda\Omega^{2}M+\beta^{\prime}\Gamma)^{2}]|<\epsilon.$
(A.23)

First, we show that

$|E_{P_{n}}[(\beta^{\prime}\gamma^{\prime-1}+\lambda\Omega^{2})^{2}((-\lambda\Omega^{2}M+\beta\Gamma)^{2}-(-\lambda\Omega^{2}M+\beta^{\prime}\Gamma)^{2})]|<\epsilon/2.$
(A.24)

Indeed since $(\beta^{\prime}\gamma^{\prime-1}+\lambda\Omega^{2})^{2}$ is
bounded by some constant $C$, we consider the expression the following set of
inequalities:

$\displaystyle|E_{P_{n}}[((-\lambda\Omega^{2}M+\beta\Gamma)^{2}-(-\lambda\Omega^{2}M+\beta^{\prime}\Gamma)^{2})]|$
$\displaystyle\leq
E_{P_{n}}[|\beta^{{}^{\prime}2}-\beta^{2}|\Gamma^{2}]+E_{P_{n}}[2|\beta-\beta^{\prime}|\lambda\Omega^{2}M\Gamma]|$
$\displaystyle\leq
E_{P_{n}}[|\beta^{{}^{\prime}2}-\beta^{2}|\Gamma^{2}]+CE_{P_{n}}[|\beta-\beta^{\prime}|\lambda(M^{2}+\Gamma^{2})]|.$

Since both $x^{2}$ and $x$ are uniformly continuous functions in the compact
domain and $P_{n}\xRightarrow{W_{2}}P$ implies that $E_{P_{n}}\Gamma^{2}$ and
$E_{P_{n}}(M^{2}+\Gamma^{2})$ is uniformly bounded we can choose $\delta$
small enough so that equation (A.24) holds.

Next we show that

$|E_{P_{n}}[((\beta^{\prime}\gamma^{\prime-1}+\lambda\Omega^{2})^{2}-(\beta\gamma^{-1}+\lambda\Omega^{2})^{2})(-\lambda\Omega^{2}M+\beta^{\prime}\Gamma)^{2})]|<\epsilon/2.$
(A.25)

Using the facts that $\Omega$ is bounded, $\beta\gamma^{-1}$ and
$\beta^{2}\gamma^{-2}$ are uniformly continuous functions in the compact
domain (the domain also being bounded away from $0$), we have
$((\beta^{\prime}\gamma^{\prime-1}+\lambda\Omega^{2})^{2}-(\beta\gamma^{-1}+\lambda\Omega^{2})^{2})<\varepsilon$
for sufficiently small $\delta$. Also, $P_{n}\xRightarrow{W_{2}}P$ implies
that $E_{P_{n}}(-\lambda\Omega^{2}M+\beta^{\prime}\Gamma)^{2}$ is uniformly
bounded in $n$, so that choosing $\varepsilon$ small enough, equation (A.25)
holds. Finally we use triangle inequality to show (A.23) from (A.24) and
(A.25). The equicontinuity of $G$ can be shown along similar lines; hence we
omit the proof.

Step 4: A calculation analogous to the one shown above proves that the
following statements hold:

$\displaystyle\lim_{n}\sup_{(\gamma,\beta)\in\mathcal{C}}|E_{Q_{n}}F_{n}(\beta,\gamma)-E_{P}F(\beta,\gamma)|$
$\displaystyle=0,$ (A.26)
$\displaystyle\lim_{n}\sup_{(\gamma,\beta)\in\mathcal{C}}|E_{Q_{n}}G_{n}(\beta,\gamma)-E_{P}G(\beta,\gamma)|$
$\displaystyle=0.$ (A.27)

An argument similar to the one used in Step 3 can be used to conclude that
with probability at least $1-n^{-c}$ the following holds

$\|(\beta_{*},\gamma_{*})-(\beta^{\dagger},\gamma^{\dagger})\|=O(\epsilon_{n})$
(A.28)

Using equations (A.16) and (A.28) we conclude that the desired result holds. ∎


# Supplementary Materials for:
Ideal refocusing of an optically active spin qubit under strong hyperfine
interactions

Leon Zaporski1,$\dagger$ Noah Shofer1,* Jonathan H. Bodey1,* Santanu Manna2,*
George Gillard3 Daniel M. Jackson1 Martin Hayhurst Appel1 Christian Schimpf2
Saimon Covre da Silva2 John Jarman1 Geoffroy Delamare1 Gunhee Park1 Urs
Haeusler1 Evgeny A. Chekhovich3 Armando Rastelli2 Dorian A. Gangloff1,4 Mete
Atatüre1,$\dagger$ Claire Le Gall1,$\dagger$ Cavendish Laboratory, University
of Cambridge, JJ Thomson Avenue, Cambridge, CB3 0HE, United Kingdom  Institute
of Semiconductor and Solid State Physics, Johannes Kepler University,
Altenbergerstraße 69, Linz 4040, Austria  Department of Physics and Astronomy,
University of Sheffield, Sheffield, S3 7RH, United Kingdom  Department of
Engineering Science, University of Oxford, Parks Road, Oxford, OX1 3PJ

* These authors contributed equally to this work.
$\dagger$ Correspondence should be addressed to<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Contents

  1. I Sample structure
  2. II Experimental Set-up for spin-control
  3. III Supplementary Measurements
    1. III.1 Saturation curve
    2. III.2 Electron spin $T_{1}$
    3. III.3 Gate Fidelity
  4. IV Free Induction Decay: data analysis
    1. IV.1 Experimental details
      1. IV.1.1 Fit function and parameters
    2. IV.2 Constraining the total number of nuclei
  5. V Reconstructing distribution of quadrupolar shifts from NMR data
    1. V.1 Inverse NMR
    2. V.2 Integral NMR
    3. V.3 Fits to the data
  6. VI Modelling of a dynamically-decoupled electron spin
    1. VI.1 Hamiltonian of the electron-nuclear system
    2. VI.2 Hamiltonian of the dynamically-decoupled electron
    3. VI.3 Overview of existing models of dynamical decoupling
    4. VI.4 Conventions
      1. VI.4.1 Ensemble averaging
      2. VI.4.2 Keldysh contour
    5. VI.5 Ring diagram formalism
    6. VI.6 Generalised $\mathbf{T}$-matrix
  7. VII Numerical simulation
    1. VII.1 Simultaneous fit to the data
  8. VIII Material constants and constrained parameters

##  I Sample structure

Figure 1: QD device sample structure design (Linz Ref. SA0553)

Fig. 1 shows the heterostructure containing the quantum dots used in this
work. The structure is grown via molecular beam epitaxy (MBE). Quantum dots
are grown via nanohole etching with $\mathrm{Al}$-droplets and infilling of
the $\mathrm{Al_{0.33}Ga_{0.67}As}$ layer with GaAs da Silva  _et al._ (2021);
Heyn  _et al._ (2009); Atkinson  _et al._ (2012). During growth, the quantum
dot containing layer is embedded between a heavily n-doped (n+) and heavily
p-doped (p+, p++) layers of $\mathrm{Al_{0.15}Ga_{0.85}As}$ and
$\mathrm{Al_{0.15}Ga_{0.85}As}/\mathrm{GaAs}$, respectively. To make the n-i-
p-diode, contacts are fabricated onto the n+ and the top GaAs p++ layers. The
back n-contact is fabricated by first etching down roughly $300$ nm using a
sulfuric acid/hydrogen peroxide solution (1 $\mathrm{H}_{2}\mathrm{SO}_{4}$ :
8 $\mathrm{H}_{2}\mathrm{O}_{2}$ : 80 $\mathrm{H}_{2}\mathrm{O}$). About $100$
nm of AuGeNi alloy is then evaporated onto the surface of the sample to form
the n-contact. The sample is then annealed at $220{}^{\circ}\mathrm{C}$ for
$30$ seconds, followed by another step at $430{}^{\circ}\mathrm{C}$ for $160$
seconds. The top p-contact is fabricated by evaporating $10$ nm of Ti onto the
top GaAs layer, followed by $100$ nm of Au. Prior to metal evaporation the
surface of the sample is dipped in a $\sim 10-50\%$ HCl solution to remove any
native surface oxides that negatively affect the conductance of the top
contact. The contacts are then wirebonded to a chip carrier to provide
electrical connections.

Figure 2: Experimental set-up: AOM $=$ acousto optic modulator; EOM $=$
electro-optic modulator; APD $=$ avalanche photodiode; AWG $=$ arbitrary wave
generator. Colors of the laser beams correspond to the driven processes
(readout, depolarization and spin control) outlined in the energy level
diagram.

##  II Experimental Set-up for spin-control

The experimental setup is outlined in Fig. 2, together with the energy level
diagram showing three processes: readout, depolarization and spin control -
each driven by one of the three, color-coded lasers.

Spin control is done via a Ti:Sapphire pulsed laser (Mira Optima 900-P) in a
pico-second mode detuned by $500-600$ $\mathrm{GHz}$ from the
resonanceBerezovsky  _et al._ (2008); Press  _et al._ (2008). Pulses are split
into two optical paths, one of which contains a delay stage which allows to
set an arbitrary time-offset between pulses in both paths. Readout and
depolarization are effectuated via two continuous-wave diode lasers.

All the lasers are fed through acousto-optic modulators (AOMs) to enable pulse
picking essential for programming the pulse sequences. To improve the
suppression of a powerful depolarization laser, it is also passed through an
electro-optical modulator (EOM).

The modulators are locked to $76$-$\mathrm{MHz}$ repetition rate of the pulsed
laser and programmed using a combination of delay generators with $8$
$\mathrm{ps}$ jitter as well as an arbitrary wave generator (AWG).

The optical paths are combined in a confocal microscope and sent down the bath
cryostat to the sample held at $T=4$ $\mathrm{K}$. A $B$-field perpendicular
to the growth axis is generated by a superconducting magnet.

To suppress the resonant laser background we use cross-polarized linear
polarizers with a quarter wave-plate in between to convert all laser
polarizations arriving to the sample from linear to circular. Following the
polarization suppresion, emission from the QD is spectrally filtered from the
spin control pulses on a diffraction grating with a $20$ $\mathrm{GHz}$
bandwidth, and then detected by an avalanche photo diode (APD).

##  III Supplementary Measurements

###  III.1 Saturation curve

In order to express the laser powers (stated in the units of the voltage
measured on a photodiode in the confocal microscope) via system-specific
parameters, we measure the saturation curve of counts scattered when
resonantly driving the transition between the crystal ground state
($\ket{\mathrm{c.g.s.}}$) and $\ket{\Uparrow\downarrow}$ state at $B=0$
$\mathrm{T}$. We fit the two-level system model:

$I(P)=\frac{I_{\mathrm{max}}}{1+(P/P_{\mathrm{sat}})^{-1}}$ (1)

to the measured countrates $I$ for all of the considered laser powers to find
$P_{\mathrm{sat}}=2.0(2)\times 10^{-7}$ $\mathrm{V}$ and
$I_{\mathrm{max}}=1.03(3)$ $\mathrm{MCounts/s}$.

Figure 3: Saturation curve of counts scattered when driving the
$\ket{\mathrm{c.g.s.}}\to\ket{\Uparrow\downarrow}$ transition resonantly. The
red points are the countrates measured for a range of fixed laser powers,
whereas the solid black curve is a fit to the data.

###  III.2 Electron spin $T_{1}$

To place a lower bound on the electron spin $T_{1}$ we perform a population-
relaxation measurement using two resonant lasers set to frequencies
$\nu_{1}=380767$ $\mathrm{GHz}$ and $\nu_{2}=380793$ $\mathrm{GHz}$ (Fig. 1e
of the main text), at the magnetic field of $6.5$ $\mathrm{T}$ (electron
Zeeman splitting $\nu_{\mathrm{e}}=4.453(6)$ $\mathrm{GHz}$). At the beginning
of the experimental sequence, the first laser initialises the electron spin in
the $\ket{\downarrow}$ state, whereas the second laser re-initialises it in
$\ket{\uparrow}$ simultaneously revealing the (reference) population of
$\ket{\downarrow}$. Following the delays of $1$, $199$ and $399$ $\mu s$, the
second laser is pulsed again to measure the population of $\ket{\downarrow}$
again; photon counts recorded during the experiments with said delays are
plotted in the panels a, b and c of the Fig. 4, respectively.

Small electron spin relaxation is observed over $399$ $\mu\mathrm{s}$,
however, measured pick-up in $\ket{\downarrow}$-population following this
delay is negligibly small to the reference. This allows us to place a weak
lower bound $T_{1}>400$ $\mu\mathrm{s}$.

Figure 4: a Measurement of electron spin relaxation after delay of $1$
$\mu\mathrm{s}$; insets correspond to transitions driven during subsequent
parts of the experiment (consistent across all panels). b Measurement of
electron spin relaxation after delay of $199$ $\mu\mathrm{s}$. c Measurement
of electron spin relaxation after delay of $399$ $\mu\mathrm{s}$.

###  III.3 Gate Fidelity

Increasing the number of $\pi$-pulses - $N_{\pi}$ \- in the CPMG sequence
exposes the electron spin to increased amount of laser-induced dephasingBodey
_et al._ (2019). During each pulse the Bloch vector’s magnitude scales by a
fraction $0\%\leq\mathcal{F}\leq 100\%$, i.e. _gate fidelity_. In this simple
model, visibility recorded at the beginning of the CPMG sequence
($v_{\mathrm{max}}=v(\tau\to 0)$ \- see Eq. 2 of the main text) with $N_{\pi}$
pulses is given simply by $\mathcal{F}^{N_{\pi}}$.

To measure the $\pi$-gate fidelity $\mathcal{F}$ we run a short-time CPMG
experiment with $N_{\pi}=162$ pulses and a total time of $\tau=4.26$
$\mu\mathrm{s}$. Keeping the preparation-to-readout delay fixed to $6.50$
$\mu\mathrm{s}$, i.e. much shorter than the electron spin $T_{1}$, we find
$v_{\mathrm{max}}=0.32(3)$.

To do this we first integrate the photon counts within $41$ $\mathrm{ns}$ wide
windows of the histograms presented in panels a and b of the Fig. 5 (shaded
areas). Integrated counts of the light-shaded areas are a measure of
background, and they are subtracted from the integrated counts of the
corresponding dark-shaded areas in order to isolate the effect of dephasing
from that of the imperfect state preparation. The signal in panel b has been
obtained after CPMG sequence preceded by an additional $\pi$-pulse.
Consequently, our measurement reveals the visibility as defined in Eq. 2 of
the main text. The errors on the integrated counts are assumed to follow shot
noise statistics, and are subsequently propagated into the expression for
visibility. This analysis allows to quantify the best-achieved $\pi$-gate
fidelity:

$\mathcal{F}=99.30(5)\%$ (2)

This fidelity exceeds the best fidelities obtained in InGaAs/GaAs QDs
structures Bodey  _et al._ (2019). The QD device structure used here does not
contain a blocking-barrier – an interface Houel  _et al._ (2012) which may be
responsible for limiting gate fidelities in other devices Stockill  _et al._
(2016); Berezovsky  _et al._ (2008); Press  _et al._ (2008); Bodey  _et al._
(2019).

Measuring the coherence of dynamically decoupled electron spin requires hours
of integration time. Mechanical stability of our current set-up on hour-long
timescales bounds the achieved fidelity to $\bar{\mathcal{F}}<\mathcal{F}$,
which we now constrain, for completeness. In this analysis, for every single
recorded CPMG data set (plotted, following normalisation, in the Fig. 10) we
find $v_{\mathrm{max}}$ from the stretched exponential fit
($f(\tau)=v_{\mathrm{max}}\times\exp{-(\tau/T_{2})^{\alpha}}+c$). We plot the
fitted values of $v_{\mathrm{max}}$ as a function of $N_{\pi}$ in the Fig. 6,
and find the average gate fidelity as $\bar{\mathcal{F}}=97.81(5)\%$ from a
fit of $v_{\mathrm{max}}(N_{\pi})=\bar{\mathcal{F}}^{N_{\pi}}$ to the data
(except for $N_{\pi}=1,3$ and $9$, due to experimental errors or lack of a
visibility measurement).

Figure 5: (a) Readout signal following CPMG sequence with $N_{\pi}=162$.
Integrated counts under dark-shaded area were equal to 3062, and 1309 under
light-shaded area (background). (b) Readout signal following the same pulse
sequence with an additional $\pi$-pulse added after preparation. Integrated
counts under dark-shaded area were equal to 2206, and 1289 under light-shaded
area (background).  Figure 6: Maximal visibility as a function of the number
of $\pi$-pulses in the CPMG sequence. Datapoints were obtained from stretched-
exponential fits
($f(\tau)=v_{\mathrm{max}}\times\exp{-(\tau/T_{2})^{\alpha}}+c$) to the
unnormalized dynamical decoupling data; these are plotted in Fig. 10 following
normalization.

##  IV Free Induction Decay: data analysis

###  IV.1 Experimental details

The experiment was performed at $B_{z}=3.5$ $\mathrm{T}$ and applied gate
voltage of $517$ $\mathrm{mV}$ (c.f. Fig. 1e of the main text). The readout
laser frequency was set to that of the highest-frequency transition -
$380754.5$ $\mathrm{GHz}$ \- and its power was stabilised around
$0.9P_{\mathrm{sat}}$. The depolarisation laser was set to drive the two
highest-energy transitions (frequency of $380751$ $\mathrm{GHz}$) with a power
of $7P_{\mathrm{sat}}$. To compensate for the small offset of depolarisation
laser’s polarization from circular, its frequency did not coincide with the
arithmetic mean of the two frequencies, exactly. Visibility constructed
according to Eq. 2 of the main text (which coincided with normalized spin
projection $S_{y}/S_{\mathrm{max}}$ in this instance - see the main text) was
taken as a signal in further parts of data analysis. Preparation of the
opposite spin states relied on insertion of the additional $\pi$-pulse
following the first resonant laser pulse.

Figure 7: FID signal (solid red curve with faint grey error bars), overlaid
with the damped sinusoidal fit
$f(\tau)=A\sin(2\pi\nu_{\mathrm{L}}\tau+\phi)e^{-(\tau/T_{2}^{*})^{\alpha}}+C$
(solid black curve). Dotted curves represent the fit envelopes,
$Ae^{-(\tau/T_{2}^{*})^{\alpha}}+C$ .

####  IV.1.1 Fit function and parameters

The signal was fitted with
$f(\tau)=A\sin(2\pi\nu_{\mathrm{L}}\tau+\phi)e^{-(\tau/T_{2}^{*})^{\alpha}}+C$,
where:

  * •

$A=0.93(2)$

  * •

$\nu_{\mathrm{L}}=2.398(3)$ $\mathrm{GHz}$

  * •

$\phi=1.65(2)$

  * •

$T_{2}^{*}=2.55(5)$ $\mathrm{ns}$

  * •

$\alpha=2.3(1)$

  * •

$C=-0.0234(3)$.

The envelope and fit (offset by the mean normalized countrate) are presented
in Fig. 7, overlaid on top of the signal. Fitted frequency of Larmor
precession $\nu_{\mathrm{L}}$ corresponds to the $g$-factor of magnitude
$|g_{\mathrm{e}}|=0.04895(6)$ \- in a good agreement with $|g_{e}|=0.0495(9)$
extracted from the RF map of $X^{-}$ in co-tunnelling regime (c.f. Fig. 1e of
the main text).

###  IV.2 Constraining the total number of nuclei

The drop of coherence in Free Induction Decay (FID) by a factor $e$ is
observed over a timescale $T_{2}^{*}$. This characterizes the magnitude of
shot-to-shot fluctuations of the Overhauser field $B_{\mathrm{OH}}$, and the
relationship is expressed via:

$T_{2}^{*}=\frac{\sqrt{6}\hbar}{g_{\mathrm{e}}\mu_{\mathrm{B}}\sqrt{\Delta^{2}B_{\mathrm{OH}}}}$
(3)

Calculating the variance of the Overhauser field from the infinite temperature
distribution of nuclear spin projections yields:

$\Delta^{2}B_{\mathrm{OH}}=\frac{2}{(g_{\mathrm{e}}\mu_{\mathrm{B}}\sqrt{N})^{2}}\sum_{k}\eta_{k}A_{k}^{2}I_{k}(I_{k}+1)$
(4)

where the sum runs over distinct nuclear species with hyperfine constants,
$A_{k}$, and concentrations, $\eta_{k}$, equal to those stated in the Table 1.
This allows us to express the total number of nuclei through:

$N=\frac{5}{4\hbar^{2}}\sum_{k}\eta_{k}A_{k}^{2}(T_{2}^{*})^{2}$ (5)

Using the measured value of $T_{2}^{*}=2.55(5)\,\mathrm{ns}$, we constrain the
total number of nuclei to $N=6.5(3)\times 10^{4}$.

##  V Reconstructing distribution of quadrupolar shifts from NMR data

Figure 8: (a) Inverse NMR spectrum showing a narrow central transition and
broadened satellite transitions. The x-axis states the frequency offset of
applied $\nu_{\mathrm{gap}}$ from central transition. The horizontal bars with
arrows indicate the positions of the fixed edges and directions of the
integral NMR scans shown in the panels below. (b) Integral NMR signal from the
Central Transition. (c) Integral NMR signal from the higher-frequency
Satellite Transition taken with a long NMR pumping time of $95$ $\mathrm{s}$.
(d) Integral NMR signal from the higher-frequency Satellite Transition taken
with a short NMR pumping time of $2$ $\mathrm{s}$

In order to reconstruct the distribution of quadrupolar shifts,
$P_{\mathrm{As}}(\nu)$, we perform NMR measured optically via
photoluminescence detection on another piece of the same wafer (i.e. a
different QD device).

Following the optical preparation of a low-temperature nuclear spin state
Chekhovich  _et al._ (2017), a burst of radiofrequency oscillating magnetic
field is sent to the QD for a time sufficiently long to equilibrate the
populations of addressed subsets of nuclear spin states. Regardless of whether
we apply ‘Inverse NMR’ or ‘Integral NMR’ excitation schemes, the direct
observables are the hyperfine splittings of $X^{-}$ transitions that reveal
$\sum_{i}A^{\mathrm{(e)}}_{i}\langle\hat{I}^{i}_{z}\rangle$ and
$\sum_{i}A^{(\mathrm{h})}_{i}\langle\hat{I}^{i}_{z}\rangle$, where hyperfine
constants $A^{(\mathrm{e})}_{i}$ and $A^{(\mathrm{h})}_{i}$ are those of the
electron and the hole, respectively. Under radiofrequency-induced relaxation
the splittings change, informing about the spectral overlap of nuclear spin
transitions with the radio-frequency comb. Such measurement is therefore
sensitive to Zeeman splitting of different nuclear species, together with the
quadrupolar shifts.

In all these experiments, a magnetic field is applied parallel to the growth
axis (Faraday geometry). External strain causing a constant quadrupolar shift
of $\sim 250$ $\mathrm{kHz}$ (for Arsenic nuclei) is applied along the $[110]$
crystal axis, and perpendicular to the growth axis in order to detect the
broad spectral features of the $+\frac{1}{2}\to+\frac{3}{2}$ satellite
transition, unobstructed by the other two NMR transitions of the spin-$3/2$
nuclei.

###  V.1 Inverse NMR

In the Inverse NMR excitation scheme, a broad frequency comb with a $6$
$\mathrm{kHz}$-wide gap centered around the variable frequency
$\nu_{\mathrm{gap}}$ is sent to the QD Bulutay  _et al._ (2014). The pulse
drives all the allowed NMR transitions except for those in the frequency gap.

Panel a of Fig. 8 shows the complete Inverse NMR spectrum of the Arsenic-75
nuclei, featuring a sharp central transition and the two broadened
transitions.

###  V.2 Integral NMR

In the Integral NMR excitation scheme, a frequency comb of increasing width
(from 0 $\mathrm{kHz}$ to 600 $\mathrm{kHz}$) is sent to the QDRagunathan
(2019). Throughout the measurement either the lower frequency edge ($38.45$
$\mathrm{MHz}$) or higher frequency edge ($38.92$ $\mathrm{MHz}$) of the comb
is fixed (see Fig. 8a). Increasing the width of the comb (i.e. sweeping the
frequency of a moving edge: $\nu_{\mathrm{edge}}$) relaxes the populations of
an increasing number of nuclei giving rise to a signal propotional to the
integral of the spectral distribution of NMR transitions. Such measurement is
more sensitive to broad, weak features which are much harder to resolve under
the Inverse NMR scheme.

Panel b of Fig. 8 shows the integral NMR signal obtained via sweeping the comb
edge over the central transition. Panels c and d of Fig. 8 depict the
corresponding signal for a high-frequency (quadrupolar-shifted) satellite
transition with long and short NMR pumping times, respectively. Regardless of
the pumping time and direction of sweep the integral NMR signal picks up a
$200$ $\mathrm{kHz}$-broad spectral feature.

###  V.3 Fits to the data

Figure 9: (a) Inverse NMR signal (black dots) alongside the skew normal fit
(solid red curve). (b) Integral NMR signal (black dots) alongside scaled
$\mathrm{CDF}_{(\mathrm{A})}(\nu)$ (solid red curve), and fitted
$\mathrm{CDF}(\nu)$ (solid blue curve) which accommodates a residual Gaussian
feature.

First, we fit the spectral lineshape of the high-frequency satellite
transition addressed via Inverse NMR with a scaled skew normal distribution:

$S_{(\mathrm{A})}(\nu)=A\frac{2}{\sigma_{(\mathrm{A})}}\phi\bigg{(}\frac{\nu-\mu_{(\mathrm{A})}}{\sigma_{(\mathrm{A})}}\bigg{)}\Phi\bigg{(}\zeta_{(\mathrm{A})}\bigg{(}\frac{\nu-\mu_{(\mathrm{A})}}{\sigma_{(\mathrm{A})}}\bigg{)}\bigg{)}+C$
(6)

where:

$\begin{split}\phi(x)&=\frac{1}{\sqrt{2\pi}}e^{-x^{2}/2}\\\
\Phi(x)&=\int_{-\infty}^{x}dx\,\phi(x)=\frac{1}{2}\bigg{(}1+\erf\bigg(\frac{x}{\sqrt{2}}\bigg{missing})\bigg{)}\end{split}$
(7)

The best fit (shown in Fig. 9a alongside the data) gives:

  * •

$A=158(1)$ $\mu\mathrm{eV}\times\mathrm{kHz}$

  * •

$\mu_{(\mathrm{A})}=252.21(4)$ $\mathrm{kHz}$

  * •

$\sigma_{(\mathrm{A})}=9.74(8)$ $\mathrm{kHz}$

  * •

$\zeta_{(\mathrm{A})}=3.22(8)$

  * •

$C=0.72(3)$ $\mu\mathrm{eV}$

The cumulative distribution function corresponding to the fitted skew normal
distribution is:

$\mathrm{CDF}_{(\mathrm{A})}(\nu)=\Phi\bigg{(}\frac{\nu-\mu_{(\mathrm{A})}}{\sigma_{(\mathrm{A})}}\bigg{)}-2T\bigg{(}\frac{\nu-\mu_{(\mathrm{A})}}{\sigma_{(\mathrm{A})}},\zeta_{(\mathrm{A})}\bigg{)}$
(8)

where $T(x,y)$ is an Owen’s $T$ function.

We juxtapose the $\mathrm{CDF}_{(\mathrm{A})}(\nu)$ (the red curve in Fig. 9b)
with the normalized Integral NMR data to arrive at a clear mismatch,
suggestive of the presence of two sub-ensembles ($\mathrm{A}$) and
($\mathrm{B}$) featuring lower and higher degrees of broadening, respectively.

This inconsistency is best resolved by fitting the integral NMR data with:

$\mathrm{CDF}(\nu)=\beta\mathrm{CDF}_{(\mathrm{B})}(\nu)+(1-\beta)\mathrm{CDF}_{(A)}(\nu)$
(9)

where $\mathrm{CDF}_{(\mathrm{A})}(\nu)$ is fixed by the inverse NMR fit and
the added weighted residual $\mathrm{CDF}_{(\mathrm{B})}(\nu)$ corresponds to
a broad Gaussian feature:

$\mathrm{CDF}_{(\mathrm{B})}(\nu)=\Phi\bigg{(}\frac{\nu-\mu_{(\mathrm{B})}}{\sigma_{(\mathrm{B})}}\bigg{)}$
(10)

Best fit is found for the following set of parameters:

  * •

$\mu_{(\mathrm{B})}=246(2)$ $\mathrm{kHz}$

  * •

$\sigma_{(\mathrm{B})}=73(4)$ $\mathrm{kHz}$

  * •

$\beta=0.244(7)$

The (scaled) fit is presented in the Fig. 9b as the solid blue curve.

Next, by differentiating the fitted $\mathrm{CDF}(\nu)$ we reconstruct the
probability distribution of Arsenic quadrupolar shifts:

$\begin{split}P_{\mathrm{As}}(\nu)&=\beta\frac{1}{\sigma_{(\mathrm{B})}}\phi\bigg{(}\frac{\nu-\mu_{(\mathrm{B})}}{\sigma_{(\mathrm{B})}}\bigg{)}\\\
&+(1-\beta)\frac{2}{\sigma_{(\mathrm{A})}}\phi\bigg{(}\frac{\nu-\mu_{(\mathrm{A})}}{\sigma_{(\mathrm{A})}}\bigg{)}\Phi\bigg{(}\zeta_{(\mathrm{A})}\bigg{(}\frac{\nu-\mu_{(\mathrm{A})}}{\sigma_{(\mathrm{A})}}\bigg{)}\bigg{)}\end{split}$
(11)

Finally, to use this distribution in our CPMG model (see section VI), we must
account for the fact that the NMR experiments were realized in Faraday
geometry ($B$-field parallel to the growth axis) and under external applied
strain, while the all-optical readout and control of the electron spin
required operation in Voigt geometry ($B$-field perpendicular to the growth
axis) and were performed on unstrained piece of the same wafer. We assume that
the broadening of the sub-ensemble ($\mathrm{B}$) results from random alloying
in $\mathrm{Al}_{x}\mathrm{Ga}_{1-x}\mathrm{As}$, which leads to the same
contribution to $P_{\mathrm{As}}(\nu)$ regardless of which geometry is
considered, due to the rotation symmetry of zincblende structure. To factor
out the effect of the externally applied strain, we further centralize the
entire distribution via translating it by $\mu_{(\mathrm{B})}$. The sub-
ensemble $(\mathrm{A})$ is assumed to feature quadrupolar broadening coming
purely from the residual strain, which is pinned to the growth axis. This
necessitates a transformation of the quadrupolar transition frequencies as
$\nu\to-\nu/2$ when transforming the distribution from Faraday to Voigt
geometry Chekhovich  _et al._ (2018). The above considerations lead us to:

$\begin{split}P_{\mathrm{As}}(\nu)=\beta\frac{1}{\sigma_{(\mathrm{B})}}\phi\bigg{(}\frac{\nu}{\sigma_{(\mathrm{B})}}\bigg{)}+(1-\beta)\frac{4}{\sigma_{(\mathrm{A})}}\phi\bigg{(}\frac{\nu-(\mu_{(\mathrm{B})}-\mu_{(\mathrm{A})})}{\sigma_{(\mathrm{A})}/2}\bigg{)}\Phi\bigg{(}-\zeta_{(\mathrm{A})}\bigg{(}\frac{\nu-(\mu_{(\mathrm{B})}-\mu_{(\mathrm{A})})}{\sigma_{(\mathrm{A})}/2}\bigg{)}\bigg{)}\end{split}$
(12)

In the process of fitting we made use of the SciPy implementation of the error
function $\erf(x)$ and the Owen’s $T$ function $T(x,y)$ Virtanen  _et al._
(2020).

##  VI Modelling of a dynamically-decoupled electron spin

###  VI.1 Hamiltonian of the electron-nuclear system

In the absence of a magnetic field the Hamiltonian for electron-nuclear system
is given by:

$\hat{H}(B=0)=\underbrace{\sum_{i}A_{i}\hat{\mathbf{S}}\cdot\hat{\mathbf{I}}_{i}}_{\hat{H}_{\mathrm{hf}}}+\underbrace{\sum_{i}\omega^{i}_{\mathrm{Q}}(\hat{\mathbf{n}}_{i}\cdot\hat{\mathbf{I}}_{i})^{2}}_{\hat{H}_{\mathrm{Q}}}$
(13)

The term denoted by $\hat{H}_{\mathrm{hf}}$ stands for the hyperfine
interaction between the electron and the nuclei, whereas
$\hat{H}_{\mathrm{Q}}$ captures the quadrupolar interaction of high-spin
($I_{i}=\tfrac{3}{2}$) nuclei originating from strain-induced and alloying-
induced electric field gradients. The unit vector $\hat{\mathbf{n}}_{i}$
defines the unique principal axis of an electric field gradient tensor Abragam
(1961).

In the regime of high magnetic fields ($B>1$ T) aligned with the $z$-axis, the
Zeeman interaction ($\hat{H}_{\mathrm{z}}$) splits the electron spin states,
as well as the nuclear spin states, modifying the system Hamiltonian in the
following way:

$\hat{H}(B\neq
0)=\hat{H}(B=0)+\underbrace{\omega_{\mathrm{e}}\hat{S}_{z}+\sum_{i}\omega^{i}_{\mathrm{n}}\hat{I}^{i}_{z}}_{\hat{H}_{\mathrm{z}}}$
(14)

The energy scales involved in the dynamics are:

$\begin{split}\omega_{\mathrm{e}}&\sim
1\,\mathrm{GHz},\quad\omega^{i}_{\mathrm{n}}\sim 10\,\mathrm{MHz}\\\
A_{i}&\sim 100\,\mathrm{kHz},\quad\omega^{i}_{\mathrm{Q}}\sim
10\,\mathrm{kHz}\end{split}$ (15)

which dictates the direction of quantisation axis along the $z$-axis.
Consequently, the terms in the Hamiltonian that are non-diagonal in the
$\ket{S_{z},\\{m_{i}\\}_{i=1,..,N}}$ basis should be treated perturbatively.

Projecting the electron-nuclear flip-flop terms
$\propto\hat{S}_{+}\hat{I}^{i}_{-}+\hat{S}_{-}\hat{I}^{i}_{+}$ present in
$\hat{H}_{\mathrm{hf}}$ onto a low-energy subspace (and ignoring small
corrections to the electron and nuclear Zeeman interaction) via a canonical
transformation one arrives at Cywiński  _et al._ (2009):

$\hat{H}_{\mathrm{hf}}\approx\hat{S}_{z}\sum_{i}A_{i}\hat{I}^{i}_{z}+\hat{S}_{z}\sum_{i\neq
j}\frac{A_{i}A_{j}}{2\omega_{\mathrm{e}}}\hat{I}^{i}_{+}\hat{I}^{j}_{-}$ (16)

where the second term introduces the electron-mediated nuclear spin flip-
flops.

Turning the attention to the quadrupolar Hamiltonian $\hat{H}_{\mathrm{Q}}$
one can isolate the diagonal part:

$\hat{H}^{\mathrm{D}}_{\mathrm{Q}}=\sum_{i}\sum_{m_{i}}\omega^{i}_{\mathrm{Q}}\bra{m_{i}}(\hat{\mathbf{n}}_{i}\cdot\hat{\mathbf{I}}_{i})^{2}\ket{m_{i}}\outerproduct{m_{i}}{m_{i}}$
(17)

which simply changes the nuclear Zeeman transition frequencies to
$\omega^{i}_{\mathrm{n}}+\Delta^{i}_{\mathrm{Q}}$, $\omega^{i}_{\mathrm{n}}$
and $\omega^{i}_{\mathrm{n}}-\Delta^{i}_{\mathrm{Q}}$, giving rise to two
satellite transitions on top of a central transition.

In the particular case of $\hat{\mathbf{n}}_{i}\cdot\hat{\mathbf{z}}=1$ one
arrives at $\Delta^{i}_{\mathrm{Q}}=2\omega^{i}_{\mathrm{Q}}$, whereas for
$\hat{\mathbf{n}}_{i}\cdot\hat{\mathbf{z}}=0$ the quadrupolar shift is
$\Delta^{i}_{\mathrm{Q}}=-\omega^{i}_{\mathrm{Q}}$.

The non-diagonal part of quadrupolar Hamiltonian $\hat{H}_{\mathrm{Q}}$ leads
to a non-collinear interaction captured by:

$\hat{H}_{\mathrm{nc}}\approx\sum_{i}A_{\mathrm{nc}}^{i}\hat{S}_{z}\hat{I}^{i}_{x},$
(18)

where
$A^{i}_{\mathrm{nc}}\propto\omega^{i}_{\mathrm{Q}}A_{i}/\omega^{i}_{\mathrm{n}}$.
This term leads to a number of rich phenomena such as dynamic nuclear
polarization studied in Ref. Högele  _et al._ (2012), or coherence modulation
observed in Ref. Botzem  _et al._ (2016).

In our system, the rate of the processes induced via the non-collinear
interaction is $\sqrt{N}A^{i}_{\mathrm{nc}}\sim 10\,\mathrm{kHz}$. In
contrast, the hyperfine interaction leads to dynamics at rates of
$\sqrt{N}A_{i}\sim 10\,\mathrm{MHz}$, for the collinear
($\propto\hat{I}_{z}^{i}$) part, and
$(\sqrt{N})^{2}\frac{A_{i}Aj}{2\omega_{\mathrm{e}}}\sim 1\,\mathrm{MHz}$ for
the electron-mediated nuclear spin flip-flops. Consequently, the non-collinear
interaction can be ignored when modelling the electron spin decoherence.

###  VI.2 Hamiltonian of the dynamically-decoupled electron

A particularly useful picture of the electron-nuclear interaction comes from
transforming the system Hamiltonian into the frame evolving due to
$\hat{H}_{z}+\hat{H}^{\mathrm{D}}_{\mathrm{Q}}$ (c.f. ‘interaction picture’).
Following this transformation, the Hamiltonian can be recast as:

$\hat{H}_{I}(t)=g\mu_{\mathrm{B}}\Big{(}\underbrace{\sum_{i}\frac{A_{i}}{g\mu_{\mathrm{B}}}\hat{I}_{z}^{i}(t)}_{B_{\mathrm{OH}}^{\parallel}(t)}+\underbrace{\sum_{i\neq
j}\frac{A_{i}A_{j}}{2g\mu_{\mathrm{B}}\omega_{\mathrm{e}}}\hat{I}_{+}^{i}(t)\hat{I}_{-}^{j}(t)}_{(B_{\mathrm{OH}}^{\perp}(t))^{2}/2B_{\mathrm{ext}}}\Big{)}\hat{S}_{z}$
(19)

where $B_{\mathrm{OH}}^{\parallel}$ and $B_{\mathrm{OH}}^{\perp}$ are the
time-dependent longitudinal and transverse components of the nuclear
Overhauser field, respectively. In this picture, species-specific transverse
components of the Overhauser field precess under the combined effect of the
external $B$-field and the quadrupolar shifts (Eq. 17). While the external
$B$-field is uniform across the nuclei, their quadrupolar shifts follow a
$10-100$ $\mathrm{kHz}$ wide distribution (see section V) which results in the
electronic dephasing.

Since the spin-control pulses are much shorter than any system-specific
timescale, including the period of the electronic Larmor precession
$2\pi/\omega_{\mathrm{e}}$ and the inhomogeneous dephasing time $T_{2}^{*}$,
they are assumed to be instantaneous, which allows to write the interaction
Hamiltonian as:

$\hat{H}_{\mathrm{I}}(t)=f(\tau;t)g\mu_{\mathrm{B}}\Big{(}B_{\mathrm{OH}}^{\parallel}(t)+\frac{(B_{\mathrm{OH}}^{\perp}(t))^{2}}{2B_{\mathrm{ext}}}\Big{)}\hat{S}_{z}$
(20)

where the filter function $f(\tau;t)$ encodes the action of ultra-fast
$\pi$-pulses that flip the electron spin during the CPMG sequence of total
duration $\tau$. The filter function is a piecewise constant function
$f:\mathbb{R}\to\\{-1,0,+1\\}$ that changes sign during pulse sequence every
time the $\pi$-pulse is applied, and stays zero outside of the pulse sequence.

Importantly, at the time $t=0$ the state of the system is given by:

$\hat{\rho}(t=0)=\hat{\rho}_{\mathrm{e}}(t=0)\otimes\hat{\rho}_{\mathrm{n}}$
(21)

with $\hat{\rho}_{\mathrm{e}}(t=0)=\tfrac{1}{2}(1+\sigma_{x})$ and
$\hat{\rho}_{\mathrm{n}}\propto\mathds{1}$, since at $T=4$ $\mathrm{K}$ the
nuclear energy scales are much smaller that the magnitude of a thermal
fluctuation ($\hbar\omega_{\mathrm{N}}^{i}\ll k_{\mathrm{B}}T$).

To capture the electronic coherence that remains after the CPMG pulse
sequence, it is convenient to use the decoherence function defined as follows:

$W(\tau)=\frac{\bra{\uparrow}\Tr_{\mathrm{n}}\hat{\rho}(\tau)\ket{\downarrow}}{\bra{\uparrow}\Tr_{\mathrm{n}}\hat{\rho}(0)\ket{\downarrow}}$
(22)

The quasi-static noise component of $B_{\mathrm{OH}}^{\parallel}$ originating
from the shot-to-shot fluctuation of the total nuclear polarization is
filtered out via successive electron spin-inversions, and is only relevant to
the decoherence process in the absence of decoupling (see section IV). The
remaining spectral-diffusion processes that contribute to the evolution of
longitudinal component of the nuclear noise (like dipolar nuclear interactions
Cywiński  _et al._ (2009)) are very slow and estimated to proceed on tens of
microsecond timescales Neder  _et al._ (2011).To a good approximation the
longitudinal noise component $B_{\mathrm{OH}}^{\parallel}$ is uncorrelated
with the transverse noise component $B_{\mathrm{OH}}^{\perp}$, leading to
Neder  _et al._ (2011):

$W(\tau)\approx W_{\perp}(\tau)\underbrace{W_{\parallel}(\tau)}_{\sim 1}$ (23)

The interaction picture Hamiltonian is further simplified in this limit. By
transforming it into the frame evolving due to:

$\hat{H}_{0}(t)=\hat{H}_{z}+\hat{H}^{\mathrm{D}}_{\mathrm{Q}}+f(\tau;t)\hat{S}_{z}\sum_{i}A_{i}\hat{I}^{i}_{z},$
(24)

finally, one arrives at:

$\hat{H}_{\mathrm{I}}(t)\approx
f(\tau;t)g\mu_{\mathrm{B}}\frac{(B_{\mathrm{OH}}^{\perp}(t))^{2}}{2B_{\mathrm{ext}}}\hat{S}_{z}$
(25)

###  VI.3 Overview of existing models of dynamical decoupling

Due to the strong hyperfine interaction of the electron with the nuclear spins
($\omega_{\mathrm{e}}<\mathcal{A}$), as well as high nuclear homogeneity, the
environmental noise felt by the electron spin qubit is non-Gaussian.
Therefore, modeling the electron dephasing requires a treatment more advanced
than the usual Filter Function formalismCywiński  _et al._ (2008); Cywiński
(2014); Stockill  _et al._ (2016). Instead, it should be based on casting the
decoherence of the electron as a fully unitary process involving interaction
with nuclear ensemble with infinitely long bath correlation time and at an
infinite temperature. In such a model, the exact quantum evolution of the
system is calculated and averaged over all equiprobable initial conditions of
the ensemble (as also done in this work).

Models applicable here, both of which lead to equivalent results, have been
developed in Ref. Cywiński  _et al._ (2009) (using ring diagram formalism,
briefly explained in subsection VI.5) and Ref. Neder  _et al._ (2011) (using
semi-classical model). The key conclusion presented in both references was
that the decoherence function in presence of dynamical decoupling could be
straightforwardly related to a sequence-specific $\mathbf{T}$-matrix via:

$W_{\perp}(\tau)=\frac{1}{\det(1+i\mathbf{T}(\tau))}$ (26)

In particular, for a CPMG sequence with $n$ $\pi$-pulses, a
$\mathbf{T}$-matrix was found as Malinowski  _et al._ (2017):

$\mathbf{T}_{i,j}(\tau)=\frac{5A_{i}A_{j}}{2g_{\mathrm{e}}\mu_{B}B}\frac{\omega_{i,j}}{\omega^{2}_{i,j}-A^{2}_{i,j}}\Bigg{[}1-\frac{\cos{\frac{A_{i,j}\tau}{2n}}}{\cos{\frac{\omega_{i,j}\tau}{2n}}}\Bigg{]}\sin{\frac{\omega_{i,j}\tau+n\pi}{2}}e^{i\frac{\omega_{i,j}\tau+n\pi}{2}}$
(27)

where indices $i$ and $j$ run over all nuclei, and:

$\omega_{i,j}=\omega^{i}_{\mathrm{n}}-\omega^{j}_{\mathrm{n}},\quad
A_{i,j}=A_{i}-A_{j}$ (28)

Since the ensemble consists of $\sim 10^{5}$ nuclei, the size of the
$\mathbf{T}$-matrix is prohibitively large to perform any simulation, which
motivates coarse graining the nuclei into $K$ groups per nuclear species, and
multiplying $(k,l)$-th component of the matrix by factor $\sqrt{N_{k}N_{l}}$
\- where $N_{k}$ is a number of the nuclei in the $k$-th group - in order to
account for collective enhancement of dynamics of a resulting ‘macrospin’.

It is readily seen that at Tesla-strong $B$-fields where
$\omega_{\mathrm{n}}^{i}\gg A_{i}$, the dominant terms in the
$\mathbf{T}$-matrix come from the nuclei of similar Zeeman frequency. This has
motivated a phenomenological treatment of quadrupolar effects in the Ref.
Cywiński  _et al._ (2009) and Ref. Neder  _et al._ (2011), where their
presence has been simply emulated by the broadening of single-species Larmor
precession frequencies $\omega_{\mathrm{n}}^{i}$. Non-trivial broadening of
the nuclear spectra measured in our QD devices (see section V) resulting in
emergence of resolvable satellite transitions in the NMR, invites a more
careful treatment. We develop an extension of the model which treats the
quadrupolar effects from first principles, and therefore allows us to
constrain the simulation with complimentary NMR data that reveals the exact
distribution of quadrupolar shifts $\omega_{\mathrm{Q}}^{i}$. A posteriori, we
verify that whilst the effect of quadrupolar Hamiltonian
$\hat{H}_{\mathrm{Q}}$ on the state-mixing (i.e. non-collinear interactions)
can be ignored, the quadrupolar shifts to the nuclear spin transitions are
essential in capturing the decoherence of the electron spin.

Development of our theory extension follows closely the reasoning from Ref.
Cywiński  _et al._ (2009). For that reason, we stick to the notation and
conventions from that work, and we briefly outline them in the next two
subsections.

###  VI.4 Conventions

####  VI.4.1 Ensemble averaging

In the sections to follow, we write

$\langle...\rangle=\Tr_{\mathrm{n}}...$ (29)

as $\hat{\rho}_{\mathrm{n}}(t=0)\propto\mathds{1}$.

####  VI.4.2 Keldysh contour

When the evolution of the system is governed by a time-dependent Hamiltonian
$\hat{H}(t)$, the time evolution of pure quantum states is calculated via the
action of a unitary operator $\hat{U}(t)$:

$\ket{\psi(t)}=\underbrace{\mathcal{T}e^{-i\int^{t}_{0}\hat{H}(t^{\prime})dt^{\prime}}}_{\hat{U}(t)}\ket{\psi(0)}$
(30)

where $\mathcal{T}$ is the time ordering operator.

Correspondingly, action of the propagator on the density operator looks as
follows:

$\hat{\rho}(t)=\hat{U}^{\dagger}(t)\hat{\rho}(0)\hat{U}(t)$ (31)

where the hermitian adjoint of the unitary operator $\hat{U}^{\dagger}(t)$
involves the reverse time ordering operator $\tilde{\mathcal{T}}$ and the
flipped sign of the interaction:

$\hat{U}^{\dagger}(t)=\tilde{\mathcal{T}}e^{i\int^{t}_{0}\hat{H}(t^{\prime})dt^{\prime}}$
(32)

Naturally, one expects that in such case, any calculation that involves
propagating the density operator $\hat{\rho}(t)$ in time is a tedious task.

A compact shortcut is introduced via Keldysh contour path-ordering Cywiński
_et al._ (2009) imposed by the operator $\mathcal{T}_{C}$. All the operators
following $\mathcal{T}_{C}$ have to follow a time contour $0\to\infty\to 0$
(also called a Keldysh contour). In addition, an integral of an arbitrary
function along the Keldysh contour is defined via:

$\int_{C}F(t,c)dt_{c}=\int_{0}^{\infty}F(t,1)dt+\int_{\infty}^{0}F(t,-1)dt$
(33)

where $c$ has been introduced as a sign of interaction, which changes from
positive along $0\to\infty$, to negative along $\infty\to 0$.

###  VI.5 Ring diagram formalism

In solving the problem of decoherence of the electron interacting with a
nuclear ensemble, the decoherence function (defined in the Eq. 22) has been
shown to satisfy Cywiński  _et al._ (2009):

$W_{\perp}(\tau)=\Big{\langle}\mathcal{T}_{C}\exp(-i\int_{C}\hat{\mathcal{V}}(t_{c})dt_{c})\Big{\rangle}$
(34)

i.e. to be an expectation value of the Keldysh-contour ordered propagator,
generated by the following nuclear operator:

$\hat{\mathcal{V}}(t_{c})=c\bra{\uparrow}\hat{H}_{\mathrm{I}}(t)\ket{\uparrow}$
(35)

where $\hat{H}_{\mathrm{I}}(t)$ is the interaction picture Hamiltonian derived
at the end of the section VI.2 (Eq. 25) and $c$ is the sign of interaction
(see Eq. 33).

Explicitly:

$\hat{\mathcal{V}}(t_{c})=cf(\tau;t)\hat{U}_{0}^{\dagger}(t_{c})\Big{(}\sum_{i\neq
j}\frac{A_{i}A_{j}}{4\omega_{\mathrm{e}}}\hat{I}_{i}^{+}\hat{I}_{j}^{-}\Big{)}\hat{U}_{0}(t_{c})$
(36)

where:

$\hat{U}_{0}(t_{c})=\exp(-i\sum_{j}\omega^{j}_{\mathrm{n}}\hat{I}_{z}^{j}t-ic\sum_{j}\frac{A_{j}}{2}\hat{I}_{z}^{j}\int_{0}^{t}f(\tau;t^{\prime})dt^{\prime}-i\hat{H}^{\mathrm{D}}_{\mathrm{Q}}t)$
(37)

is simply a unitary transformation generated by $\hat{H}_{0}$ from Eq. 24.

By expanding the Keldysh-contour ordered exponential from Eq. 34 into Dyson
series, one arrives at
$W_{\perp}(\tau)=1+\sum_{k=1}^{\infty}W_{\perp}^{(2k)}(\tau)$ where $2k$-th
order terms in the expansion are related to Keldysh-contour ordered $2k$-point
auto-correlators of $\hat{\mathcal{V}}(t_{c})$, i.e.:

$W_{\perp}^{(2k)}(\tau)\sim\int_{C}\prod^{2k}_{i=1}dt_{c_{i}}\bigg{\langle}\mathcal{T}_{C}\prod^{2k}_{i=1}\hat{\mathcal{V}}(t_{c_{i}})\bigg{\rangle}$
(38)

Authors of Ref. Cywiński  _et al._ (2009) recognise an algebraic structure
behind these correlators which allows expressing $W_{\perp}^{(2k)}(\tau)$, for
arbitrary $k$, via weighted sums of _ring diagrams_ :

$R_{2k}(\tau)=\Tr_{\mathrm{n}}\mathbf{T}^{2k}\Big{(}1+\mathcal{O}\Big{(}\frac{1}{N}\Big{)}\Big{)},$
(39)

where for the case of $\hat{H}^{\mathrm{D}}_{\mathrm{Q}}=0$ the
$\mathbf{T}$-matrix is the same as that from the equation 27.

Central result of the Ref. Cywiński  _et al._ (2009) shows that resummation of
all the ring diagrams in the expansion of $W_{\perp}(\tau)$ results in:

$W_{\perp}(\tau)\approx\exp(\sum_{k=1}^{\infty}\frac{(-i)^{2k}}{2k}R_{2k}(\tau))\approx\exp(\sum_{k=1}^{\infty}\frac{(-i)^{2k}}{2k}\Tr_{\mathrm{n}}\mathbf{T}^{2k}(\tau))$
(40)

which, following a straight-forward calculation involving the diagonalization
of $\mathbf{T}$-matrix, results in a compact relationship:

$W_{\perp}(\tau)\approx\frac{1}{\det(1+i\mathbf{T})}$ (41)

The approximate nature of the expression from Eq. 40 comes from
$\mathcal{O}(N^{-1})$ relative error in the relation from Eq. 39, and allowing
the sum in the exponent to run up to $\infty$, rather than $N$. Both
approximations can be safely invoked in the regime of $N\sim 10^{5-6}$.

###  VI.6 Generalised $\mathbf{T}$-matrix

The Ring diagram formalism applied to the case
$\hat{H}^{\mathrm{D}}_{\mathrm{Q}}=0$ has been clearly laid out in a great
detail in Ref. Cywiński  _et al._ (2009). Here we generalise this result to
the case of $\hat{H}^{\mathrm{D}}_{\mathrm{Q}}\neq 0$ and we show that the
algebraic Ring diagram structure of the problem survives, albeit with a
generalised expression for a $\mathbf{T}$-matrix.

As in the Ref. Cywiński  _et al._ (2009), we begin by calculating the lowest
order contribution to the decoherence function ($k=1$ term from the Eq. 38):

$W^{(2)}_{\perp}(\tau)=-\frac{1}{2}\int_{C}dt_{c_{1}}\int_{C}dt_{c_{2}}\langle\mathcal{T}_{C}\hat{\mathcal{V}}(t_{c_{1}})\hat{\mathcal{V}}(t_{c_{2}})\rangle$
(42)

where $\hat{\mathcal{V}}(t_{c})$ is given by Eq. 36. The calculation starts
from expanding the $+$-rotating component of the $\mathbf{B}_{\mathrm{OH}}(t)$
in the nuclear product basis:

$\begin{split}\hat{U}_{0}^{\dagger}(t_{c})\Big{(}\sum_{j}A_{j}\hat{I}^{j}_{+}\Big{)}\hat{U}_{0}(t_{c})&=\sum_{j}A_{j}\Big{(}\sum_{m_{j}}\sqrt{I_{j}(I_{j}+1)-m_{j}(m_{j}+1)}(\hat{U}_{0}^{j})^{\dagger}(t_{c})\outerproduct{m_{j}+1}{m_{j}}\hat{U}_{0}^{j}(t_{c})\Big{)}\end{split}$
(43)

The action of unitaries on $\outerproduct{m_{j}+1}{m_{j}}$ operators simply
brings state-dependent oscillating factors:

$\exp(i(\omega_{m_{j}+1}-\omega_{m_{j}})t+ic\frac{A_{j}}{2}\int_{0}^{t}f(\tau;t^{\prime})dt^{\prime})\equiv
e_{m_{j}}(t_{c})$ (44)

Where $\omega_{m_{j}+1}-\omega_{m_{j}}$ are equal to
$\ket{m_{j}+1}\to\ket{m_{j}}$ nuclear transition frequencies; here
$\hat{H}_{\mathrm{Q}}\neq 0$ introduces quadrupolar shifts, which split these
transition frequencies into $\omega^{j}_{\mathrm{n}}+\Delta^{j}_{\mathrm{Q}}$,
$\omega^{j}_{\mathrm{n}}$ and
$\omega^{j}_{\mathrm{n}}-\Delta^{j}_{\mathrm{Q}}$ for $m_{j}=1/2$,
$m_{j}=-1/2$, and $m_{j}=-3/2$, respectively.

For convenience, we define the set of operators:

$\hat{\mathcal{O}}^{j}_{m_{j}}\equiv\sqrt{2I_{j}+1}\outerproduct{m_{j}+1}{m_{j}}$
(45)

and turn our attention to the following correlator:

$\begin{split}\langle\mathcal{T}_{C}\hat{\mathcal{V}}(t_{c_{1}})\hat{\mathcal{V}}(t_{c_{2}})\rangle&=\sum_{k\neq
l}\sum_{p\neq
q}\sum_{m_{k},m_{l},m^{\prime}_{p},m^{\prime}_{q}}\mathcal{V}_{k,l}^{m_{k},m_{l}}(t_{c_{1}})\mathcal{V}_{p,q}^{m^{\prime}_{p},m^{\prime}_{q}}(t_{c_{2}})\langle\mathcal{T}_{C}\hat{\mathcal{O}}^{k}_{m_{k}}\hat{\mathcal{O}}^{l,\dagger}_{m_{l}}\hat{\mathcal{O}}^{p}_{m^{\prime}_{p}}\hat{\mathcal{O}}^{q,\dagger}_{m^{\prime}_{q}}\rangle\end{split}$
(46)

With:

$\mathcal{V}_{k,l}^{m_{k},m_{l}}(t_{c})=cf(\tau;t)e_{m_{k}}(t_{c})e^{*}_{m_{l}}(t_{c})\underbrace{\sqrt{\frac{I_{k}(I_{k}+1)-m_{k}(m_{k}+1)}{2I_{k}+1}}}_{\equiv\sqrt{a(I_{k},m_{k})}}\underbrace{\sqrt{\frac{I_{l}(I_{l}+1)-m_{l}(m_{l}+1)}{2I_{l}+1}}}_{\equiv\sqrt{a(I_{l},m_{l})}}\frac{A_{k}A_{l}}{4\omega_{\mathrm{e}}}$
(47)

The remaining correlator is further simplified by considering the only
possible non-zero contractions of nuclear indices:

$\begin{split}\langle\mathcal{T}_{C}\hat{\mathcal{O}}^{k}_{m_{k}}\hat{\mathcal{O}}^{l,\dagger}_{m_{l}}\hat{\mathcal{O}}^{p}_{m^{\prime}_{p}}\hat{\mathcal{O}}^{q,\dagger}_{m^{\prime}_{q}}\rangle&=\delta_{k,q}\delta_{l,p}\langle\mathcal{T}_{C}\hat{\mathcal{O}}^{k}_{m_{k}}\hat{\mathcal{O}}^{k,\dagger}_{m^{\prime}_{k}}\rangle\langle\mathcal{T}_{C}\hat{\mathcal{O}}^{l,\dagger}_{m_{l}}\hat{\mathcal{O}}^{l}_{m_{l}^{\prime}}\rangle\\\
&=\delta_{k,q}\delta_{l,p}\delta_{m_{k},m_{k}^{\prime}}\delta_{m_{l},m_{l}^{\prime}}\langle\mathcal{T}_{C}\hat{\mathcal{O}}^{k}_{m_{k}}\hat{\mathcal{O}}^{k,\dagger}_{m_{k}}\rangle\langle\mathcal{T}_{C}\hat{\mathcal{O}}^{l,\dagger}_{m_{l}}\hat{\mathcal{O}}^{l}_{m_{l}}\rangle\end{split}$
(48)

The first line follows from the fact that $k\neq l$ and $p\neq q$, as well as
the lack of inter-nuclear coherences in a thermal state. Similarly, the second
line results from the lack of inter-state coherences for a single nucleus in a
thermal state.

At this stage it becomes clear that the Keldysh-contour ordering operator can
be removed. For a $k$-th nucleus, permutation of the two operators will bring
either $(2I_{k}+1)\langle\outerproduct{m_{k}+1}{m_{k}+1}\rangle$ or
$(2I_{k}+1)\langle\outerproduct{m_{k}}{m_{k}}\rangle$; for an infinite
temperature bath the populations of all nuclear spin states are equal:

$\langle\outerproduct{m_{k}}{m_{k}}\rangle=1/(2I_{k}+1)$ (49)

rendering the two results identical.

Finally:

$\begin{split}\langle\mathcal{T}_{C}\hat{\mathcal{O}}^{k}_{m_{k}}\hat{\mathcal{O}}^{l,\dagger}_{m_{l}}\hat{\mathcal{O}}^{p}_{m^{\prime}_{p}}\hat{\mathcal{O}}^{q,\dagger}_{m^{\prime}_{q}}\rangle&=\delta_{k,q}\delta_{l,p}\delta_{m_{k},m_{k}^{\prime}}\delta_{m_{l},m_{l}^{\prime}}\end{split}$
(50)

so the correlator from the equation 46 simplifies significantly:

$\langle\mathcal{T}_{C}\hat{\mathcal{V}}(t_{c_{1}})\hat{\mathcal{V}}(t_{c_{2}})\rangle=\sum_{k\neq
l}\sum_{m_{k},m_{l}}\mathcal{V}_{k,l}^{m_{k},m_{l}}(t_{c_{1}})\mathcal{V}_{l,k}^{m_{l},m_{k}}(t_{c_{2}})$
(51)

For convenience we free the sums to run over all indices, and we re-arrange
the $4$ dimensional tensors into matrices by introducing combined indices
$\alpha=(k,m_{k})$ and $\beta=(l,m_{l})$, such that:

$\langle\mathcal{T}_{C}\hat{\mathcal{V}}(t_{c_{1}})\hat{\mathcal{V}}(t_{c_{2}})\rangle=\sum_{\alpha,\beta}(1-\delta_{k,l})\mathcal{V}_{\alpha,\beta}(t_{c_{1}})(1-\delta_{l,k})\mathcal{V}_{\beta,\alpha}(t_{c_{2}})$
(52)

Finally, we write:

$W^{(2)}_{\perp}(\tau)=-\frac{1}{2}\sum_{\alpha,\beta}\mathbf{T}_{\alpha,\beta}\mathbf{T}_{\beta,\alpha}=-\frac{1}{2}\Tr_{\mathrm{n}}\mathbf{T}^{2}$
(53)

having arrived at a generalised $\mathbf{T}$-matrix of the following form:

$\mathbf{T}_{\alpha,\beta}=(1-\delta_{k,l})\int_{C}dt_{c}\mathcal{V}^{m_{k},m_{l}}_{k,l}(t_{c})$
(54)

Where $\alpha$ and $\beta$ are combined indices.

Higher-order contributions ($k>1$) to the decoherence function will also
follow exactly the ring diagram structure of the $\hat{H}_{\mathrm{Q}}=0$
problem (see Ref. Cywiński  _et al._ (2009)), since the combinatorics of
contractions of the $k$-point correlators here is analogous, as it also relies
on:

  * •

lack of inter-nuclear correlations,

  * •

lack of inter-state correlations,

  * •

equal populations of nuclear spin-states,

all resulting from the bath having an infinite temperature. Again, this leads
to the identity from Eq. 41, where now the $\mathbf{T}$-matrix has a
generalized form which we identified.

For spin-$\frac{3}{2}$ species our theory extension requires working with
$12K\times 12K$ matrices rather than $3K\times 3K$ matrices, but it straight-
forwardly captures the effect of quadrupolar broadening from first principles.
Performing the remaining integrals can be done exactly like in Ref. Cywiński
_et al._ (2009) and Malinowski  _et al._ (2017). In particular, for a CPMG
sequence with $n$ $\pi$-pulses we recover:

$\begin{split}\mathbf{T}_{\alpha,\beta}=&(1-\delta_{k,l})\sqrt{a(I_{k},m_{k})a(I_{l},m_{l})}\frac{A_{k}A_{l}}{\omega_{\mathrm{E}}}\frac{\omega^{m_{k},m_{l}}_{k,l}}{(\omega^{m_{k},m_{l}}_{k,l})^{2}-A^{2}_{k,l}}\\\
&\times\Bigg{\\{}1-\frac{\cos(\tfrac{A_{k,l}\tau}{2n})}{\cos(\tfrac{\omega^{m_{k},m_{l}}_{k,l}\tau}{2n})}\Bigg{\\}}\sin(\frac{\omega^{m_{k},m_{l}}_{k,l}\tau+n\pi}{2})\exp(i\frac{\omega^{m_{k},m_{l}}_{k,l}\tau+n\pi}{2})\end{split}$
(55)

where:

$\begin{split}A_{k,l}&\equiv A_{k}-A_{l}\\\
\omega^{m_{k},m_{l}}_{k,l}&\equiv(\omega_{m_{k}+1}-\omega_{m_{k}})-(\omega_{m_{l}+1}-\omega_{m_{l}})\end{split}$
(56)

are differences of single-nucleus hyperfine constants and transition
frequencies. Coarse-graining the $\mathbf{T}$-matrix to describe the
interactions between groups of nuclei again results in the appearance of
additional factors $\sqrt{N_{k}N_{l}}$ where $N_{k}$ and $N_{l}$ are numbers
of nuclei in $k$-th and $l$-th groups, respectively.

##  VII Numerical simulation

We use a simple python code to assemble $\mathbf{T}$-matrix and calculate the
determinant of $1+i\mathbf{T}$ for delays $\tau$ between two $\pi/2$-puses in
CPMG sequence. This is straightforwardly related to the electronic coherence
according to Eq. 41.

We set species-dependent hyperfine constants, concentrations and Zeeman
splittings present in the model to the material constants presented in the
table 1. The dot-specific parameters, i.e. the total number of nuclei and the
g-factor, are constrained through independent measurements, and summarised in
the table 2.

The fit of the model to the data relies on the distribution of quadrupolar
broadening, $P_{\mathrm{As}}(\nu)$, which is reconstructed in the section V
for another QD device made from the same wafer. In order to account for the
possible dot-to-dot variations, our model has two free parameters: the
electron’s wavefunction overlap with an inhomogeneous nuclear sub-ensemble
affected by random alloying, $\beta$, and a ‘scaling factor’, $\kappa$, that
accounts for the variation in the magnitude of quadrupolar shifts. The
‘scaling factor’ simply scales the width of distribution via
$P_{\mathrm{As}}(\nu)\to\kappa^{-1}P_{\mathrm{As}}(\nu/\kappa)$
transformation.

In order to sample the distribution $P_{\mathrm{As}}(\nu)$ reconstructed via
the Integral and Inverse NMR spectra (see section V), we partition nuclear
spins into $K=200$ groups of unequal size, each containing $N_{k}$ nuclear
spins with a quadrupolar shift $\nu_{k}$. We consider $K$ equally spaced
quadrupolar shifts $\nu_{k}$ within
$[-5\sigma_{(\mathrm{B})},5\sigma_{(\mathrm{B})}]$ interval (where
$\sigma_{(\mathrm{B})}$ is the width of the detected broad Gaussian
distribution - see section V), and we assign each shift to the following
number of nuclei:

$N_{k}=N\frac{c_{\mathrm{As}}}{2}\int_{\nu_{k}-\Delta\nu/2}^{\nu_{k}+\Delta\nu/2}d\nu\,P_{\mathrm{As}}(\nu)$
(57)

where $c_{\mathrm{As}}$ is the concentration of Arsenic, $N$ is a total number
of nuclei and the bin width $\Delta\nu=10\sigma_{(\mathrm{B})}/K$.

In order to capture the residual hetero-nuclear effects we include both
Gallium isotopes (${}^{71}\mathrm{Ga}$ and ${}^{69}\mathrm{Ga}$) in the
simulation, under the assumption that the effect of random alloying in
$\mathrm{AlGaAs}$ on their spectra of quadrupolar inhomogeneities is
negligibleChekhovich  _et al._ (2015). For each of the Gallium isotopes we
constrain the distribution of quadrupolar shifts to:

$P_{\mathrm{Ga}}(\nu)\propto P^{(\mathrm{A})}_{\mathrm{As}}(-2\nu)$ (58)

where $P^{(\mathrm{A})}_{\mathrm{As}}(\nu)$ is the narrow mode of the
Arsenic’s distribution (see section V). The relative scaling of width follows
from the measurement on ${}^{75}\mathrm{As}$ and ${}^{69}\mathrm{Ga}$ nuclei
performed in Ref. Chekhovich  _et al._ (2018), according to which:

$\frac{Q^{\mathrm{As}}S^{\mathrm{As}}_{11}}{Q^{\mathrm{Ga}}S^{\mathrm{Ga}}_{11}}\approx-2$
(59)

where $Q^{\alpha}$ is the nuclear quadrupolar moment of species $\alpha$, and
$S^{\alpha}_{11}$ is a diagonal component of its gradient-elastic tensor.

To represent the distribution of quadrupolar inhomogeneities for both Gallium
isotopes - $P_{\mathrm{Ga}}(\nu)$ \- we arrange the Gallium nuclei into bins
of the same width as those of Arsenic ($\Delta\nu=10\sigma_{(\mathrm{B})}/K$),
and retain only the bins that have no fewer nuclei than in any of Arsenic’s
bins. Such cut-off prevents the simulation from crashing due to the presence
of empty bins resulting from rounding errors.

Figure 10: (a)-(i) CPMG datasets (black dots) for $N_{\pi}=1,2,3,4,5,8,9,27$
and $81$, plotted together with best simultaneous fit of the theory curves
(red solid curve). Figure 11: Weighted residual sum of squares averaged over
all the CPMG datasets, as a function of two free parameters in a numerical
simulation: $\beta$ and scaling factor $\kappa$. Best fit is found for
$\beta=0.35$ and $\kappa=1.45$.

###  VII.1 Simultaneous fit to the data

We use our model with two free parameters, $\beta$ and $\kappa$, to
simultaneously fit the coherence decay in all the CPMG datasets presented in
Fig. 10. The goodness of the simultaneous fit is characterized by the weighted
residual sum of squares per data point ($\mathrm{WRSS}$) averaged over fits to
distinct $N_{\pi}$ datasets. For each CPMG dataset containing $M$ datapoints
$y_{i}$ with errors $\Delta_{i}$ we calculate:

$\mathrm{WRSS(N_{\pi})}=\Big{(}\sum^{M}_{i}(y_{i}-\hat{y}_{i})^{2}\tfrac{1}{\Delta_{i}^{2}}\Big{)}/\Big{(}M\sum^{M}_{i}\tfrac{1}{\Delta_{i}^{2}}\Big{)},$
(60)

where $\hat{y}_{i}$ are the corresponding model values, and we subsequently
average the expression over different $N_{\pi}$ datasets to arrive at a robust
measure: $\overline{\mathrm{WRSS}}$.

Due to the size of the $\mathbf{T}$-matrix, the simulation time is too long to
allow for an automated least-squares optimization. We thus perform a raster
scan of the $\overline{\mathrm{WRSS}}$ over models with
$(\beta,\kappa)\in\\{0.05,0.15,0.25,0.35,0.45\\}\times\\{1.000,1.225,1.450,1.675\\}$.
The calculated $\overline{\mathrm{WRSS}}(\beta,\kappa)$ has a clear local
minimum on that set, corresponding to $\beta=0.35$ and $\kappa=1.450$ (see
Fig. 11). For our best fit $\overline{\mathrm{WRSS}}=3.9\times 10^{-4}$.
Theory curves of the best simultaneous fit are plotted on top of the data in
Fig. 10.

##  VIII Material constants and constrained parameters

Species $\alpha$ | ${}^{75}\mathrm{As}$ | ${}^{71}\mathrm{Ga}$ | ${}^{69}\mathrm{Ga}$
---|---|---|---
Hyperfine interaction constant, $A_{\alpha}$ ($\mathrm{GHz}$) | $65.3$ | $69.9$ | $54.7$
Concentration, $c_{\alpha}$ | $1$ | $0.396$ | $0.604$
Zeeman splitting, $\omega_{\alpha}/B$ $(\mathrm{MHz}/\mathrm{T})$ | $7.22$ | $12.98$ | $10.22$
Table 1: Material constants consistent with Ref. Malinowski  _et al._ (2017). External magnetic field, $B$ | $6.5~{}\mathrm{T}$
---|---
Total number of nuclei, N | $6.5(3)\times 10^{4}$
Electron $g$-factor, $|g_{\mathrm{e}}|$ | $0.04895(6)$
$\pi$-gate fidelity, $\mathcal{F}$ | $99.30(5)\%$
Table 2: Parameters constrained through independent measurement or set
experimentally.

## References

  * da Silva  _et al._ (2021) S. F. C. da Silva, G. Undeutsch, B. Lehner, S. Manna, T. M. Krieger, M. Reindl, C. Schimpf, R. Trotta, and A. Rastelli, Applied Physics Letters 119, 120502 (2021), https://doi.org/10.1063/5.0057070 .
  * Heyn  _et al._ (2009) C. Heyn, A. Stemmann, T. Köppen, C. Strelow, T. Kipp, M. Grave, S. Mendach, and W. Hansen, Applied Physics Letters 94, 183113 (2009).
  * Atkinson  _et al._ (2012) P. Atkinson, E. Zallo, and O. G. Schmidt, Journal of Applied Physics 112, 054303 (2012).
  * Berezovsky  _et al._ (2008) J. Berezovsky, M. H. Mikkelsen, N. G. Stoltz, L. A. Coldren, and D. D. Awschalom, Science 320, 349 (2008), https://www.science.org/doi/pdf/10.1126/science.1154798 .
  * Press  _et al._ (2008) D. Press, T. D. Ladd, B. Zhang, and Y. Yamamoto, Nature 456, 218 (2008).
  * Bodey  _et al._ (2019) J. H. Bodey, R. Stockill, E. V. Denning, D. A. Gangloff, G. Éthier-Majcher, D. M. Jackson, E. Clarke, M. Hugues, C. L. Gall, and M. Atatüre, npj Quantum Information 5, 95 (2019).
  * Houel  _et al._ (2012) J. Houel, A. V. Kuhlmann, L. Greuter, F. Xue, M. Poggio, B. D. Gerardot, P. A. Dalgarno, A. Badolato, P. M. Petroff, A. Ludwig, D. Reuter, A. D. Wieck, and R. J. Warburton, Phys. Rev. Lett. 108, 107401 (2012).
  * Stockill  _et al._ (2016) R. Stockill, C. Le Gall, C. Matthiesen, L. Huthmacher, E. Clarke, M. Hugues, and M. Atatüre, Nature Communications 7, 12745 (2016).
  * Chekhovich  _et al._ (2017) E. A. Chekhovich, A. Ulhaq, E. Zallo, F. Ding, O. G. Schmidt, and M. S. Skolnick, Nature Materials 16, 982 (2017).
  * Bulutay  _et al._ (2014) C. Bulutay, E. A. Chekhovich, and A. I. Tartakovskii, Phys. Rev. B 90, 205425 (2014).
  * Ragunathan (2019) G. Ragunathan, Ph.D. thesis (2019).
  * Chekhovich  _et al._ (2018) E. A. Chekhovich, I. M. Griffiths, M. S. Skolnick, H. Huang, S. F. C. da Silva, X. Yuan, and A. Rastelli, Phys. Rev. B 97, 235311 (2018).
  * Virtanen  _et al._ (2020) P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, Nature Methods 17, 261 (2020).
  * Abragam (1961) A. Abragam,  _The Principles of Nuclear Magnetism_, International series of monographs on physics (Clarendon Press, 1961).
  * Cywiński  _et al._ (2009) L. Cywiński, W. M. Witzel, and S. Das Sarma, Phys. Rev. B 79, 245314 (2009).
  * Högele  _et al._ (2012) A. Högele, M. Kroner, C. Latta, M. Claassen, I. Carusotto, C. Bulutay, and A. Imamoglu, Phys. Rev. Lett. 108, 197403 (2012).
  * Botzem  _et al._ (2016) T. Botzem, R. P. G. McNeil, J.-M. Mol, D. Schuh, D. Bougeard, and H. Bluhm, Nature Communications 7, 11170 (2016).
  * Neder  _et al._ (2011) I. Neder, M. S. Rudner, H. Bluhm, S. Foletti, B. I. Halperin, and A. Yacoby, Phys. Rev. B 84, 035441 (2011).
  * Cywiński  _et al._ (2008) L. Cywiński, R. M. Lutchyn, C. P. Nave, and S. Das Sarma, Phys. Rev. B 77, 174509 (2008).
  * Cywiński (2014) L. Cywiński, Phys. Rev. A 90, 042307 (2014).
  * Malinowski  _et al._ (2017) F. K. Malinowski, F. Martins, P. D. Nissen, E. Barnes, Ł. Cywiński, M. S. Rudner, S. Fallahi, G. C. Gardner, M. J. Manfra, C. M. Marcus, and F. Kuemmeth, Nature Nanotechnology 12, 16 (2017).
  * Chekhovich  _et al._ (2015) E. A. Chekhovich, M. Hopkinson, M. S. Skolnick, and A. I. Tartakovskii, Nature Communications 6, 6348 (2015).