Approaching a fully-polarized state of nuclear spins in a solid

Magnetic noise of atomic nuclear spins is a major source of decoherence in solid-state spin qubits. In theory, near-unity nuclear spin polarization can eliminate decoherence of the electron spin qubit, while turning the nuclei into a useful quantum information resource. However, achieving sufficiently high nuclear polarizations has remained an evasive goal. Here we implement a nuclear spin polarization protocol which combines strong optical pumping and fast electron tunneling. Nuclear polarizations well above 95% are generated in GaAs semiconductor quantum dots on a timescale of 1 minute. The technique is compatible with standard quantum dot device designs, where highly-polarized nuclear spins can simplify implementations of qubits and quantum memories, as well as offer a testbed for studies of many-body quantum dynamics and magnetism.

95% are generated in GaAs semiconductor quantum dots on a timescale of 1 minute. The technique is compatible with standard quantum dot device designs, where highly-polarized nuclear spins can simplify implementations of quantum bits and memories, as well as offer a testbed for studies of many-body quantum dynamics and magnetism.
The capability of initializing a quantum system into a well-defined eigenstate is one of the fundamental requirements in quantum science and technology. This has been demonstrated for individual and dilute nuclear spins in the solid state [1, 2], but remains a long-standing challenge for dense three-dimensional lattices of nuclear spins. The spin-ensemble ground state is characterized by a polarization degree P N = ±100%, which is equivalent to absolute zero spin temperature.
Very high polarizations, P N ≈ 95 − 99%, have been demonstrated in bulk materials through brute-force cooling of the lattice, but the cooling cycle may take hours and even days [3,4].
More scalable approaches seek to use individual or dilute electron spins to polarize the dense nuclear ensembles. Microwave pumping of paramagnetic impurities in bulk solids [5,6] provides polarizations up to P N ≈ 80 − 90%. In semiconductor nanostructures, P N ≈ 50 − 80% is achieved either through electronic transport [7] or optical excitation [8]. At such polarization degrees, nuclear spin fluctuations are still similar to their thermal-equilibrium maximum, with 1 − P 2 N being the figure of merit in reducing the nuclear spin noise [9]. Therefore, different techniques are needed to approach a fully-polarized nuclear state. c Typical photoluminescence spectra of a negatively charged trion X − in an individual QD. The spectral splitting ∆E PL depends both on B z and the helicity of the optical pumping (σ ± ) due to the buildup of the nuclear spin polarization. d Experimental cycle consisting of nuclear spin optical pumping, nuclear magnetic resonance (NMR) excitation, and optical probing of the photoluminescence spectrum. V Gate is varied to switch the QD state between electron-charged (1e) and neutral (0e) or is set to an arbitrary value V Pump during pumping.
Extensive theoretical studies have been conducted to understand what limits nuclear spin pumping in a central-spin scenario, where the electron can be polarized on demand, while the ensemble of N nuclei can only be accessed through hyperfine (magnetic) coupling with that central electron ( Fig. 1a). The formation of coherent "dark" states [10] has been shown to suppress the transfer of spin from the electron to nuclei [11]. Thus an open question remains -is it possible, even in principle, to reach a fully-polarized nuclear state in a real central-spin system?
Here we work with epitaxial GaAs/AlGaAs quantum dots (QDs) and use optical techniques to polarize nuclear spins. While the optical method is well known, its bottleneck is the slow, nanosecond-scale, recombination of the photo-generated electrons. We resolve this issue by in-troducing charge transport -optical recombination is replaced with fast sub-picosecond electron tunneling. Moreover, no "dark"-state limitation occurs, which we also attribute to the extremely short lifetime of the electron spin. As a result, we achieve nuclear polarization degrees well above P N > 95%. The maximum polarizations vary between individual QDs, which we ascribe to random QD anisotropies. For the best dots we derive P N 99%, limited only by the accuracy of the existing measurement techniques. These high polarizations surpass the predicted P N 90% threshold for achieving non-trivial regimes, manifested in extended electron spin qubit coherence [12,13], quantum memory operation [13,14], superradiant electron-nuclear spin dynamics [15,16], as well as magnetic-ordering phase transition [17,18].
The semiconductor device, sketched in Fig. 1b, is a p − i − n diode with epitaxial GaAs QDs embedded into the AlGaAs barrier layers (see Supplementary Note 1). By changing the gate bias V Gate it is possible to charge the QD with individual resident electrons [19,20] and apply electric field. Each individual QD contains N ≈ 10 5 nuclei, with the three abundant isotopes 75 As, 69 Ga and 71 Ga, all possessing spin momentum I = 3/2. The sample is cooled to ≈ 4.25 K and placed in a magnetic field B z parallel to electric field and sample growth direction (see Supplementary Note 3). Thanks to the selection rules [21], optical laser excitation can create spin-polarized electronhole pairs: photons with +1 (−1) angular momentum (in units of ), corresponding to a σ + (σ − ) polarized beam, generate electrons with spin projection s z = −1/2 (s z = +1/2). Owing to the electron-nuclear hyperfine interaction (Fig. 1a), a polarized electron can transfer its spin to one of the nuclei and, through repeated optical pumping, induce a substantial polarization |P N |.
Conversely, the energy of the photon emitted from electron-hole recombination depends on the mutual alignment of s z and the total magnetic field, which is a sum of B z and the effective field of the polarized nuclear spins. The resulting optical spectrum is a doublet (Fig. 1c), whose splitting ∆E PL is used as a sensitive probe of the nuclear spin polarization state. We define the exciton hyperfine shift E hf = −(∆E PL − ∆E PL,0 ), where ∆E PL,0 is the splitting measured for depolarized nuclei (P N ≈ 0).
The high resolution optical spectra (Fig. 1c), required for accurate measurement of E hf , can only be observed for a narrow range of sample biases and optical excitation powers. In order to cover a wide range of pumping parameters we use a pump-probe technique (Fig. 1d). We maximize the hyperfine shift |E hf | by optimizing the following four parameters: the elliptical polarization of the optical pump, its power P Pump , photon energy E Pump and the sample bias V Pump during pumping.
The corresponding results are interpreted with reference to the broad range photoluminescence spectra shown in Figs. 2a, b. Fig. 2a shows low power spectra, which reveal a well-known bias- is completely quenched. Moreover, the optimal pump laser power P Pump = 1.5 mW is five orders of magnitude higher than the s-shell saturation power. Based on these observations, the nuclear spin pumping effect can be understood as a cyclic process sketched in Fig. 2d. First, circularlypolarized resonant optical excitation creates a spin-polarized electron-hole pair in the quantum dot. Then, the electron has a small but finite probability to undergo a flip-flop with one of the nuclei, increasing the ensemble polarization |P N |. Finally, in order to proceed to the next cycle, the electron is removed through tunneling. The tunneling time, estimated from bias-dependent photoluminescence in Supplementary Note 4, is 0.1 ps, much shorter than the ≈ 300 ps radiative recombination time [24]. The combination of high-power optical pumping and fast tunnel escape results in rapid cycling. This in turn leads to a high rate of nuclear spin pumping, which helps to outpace the inevitable nuclear spin relaxation. The cycling time is also much shorter than the period of coherent electron precession 20 ps, ensuring the spin-flipped electrons are removed before they can undergo a reverse flip-flop [25]. The ultimate result is a large steady-state |E hf |.
Although the hyperfine shift E hf scales linearly with polarization degree P N , its absolute value depends on the QD structure. The electron wavefunction leaks into the barriers where the fraction of Ga atoms replaced with Al atoms is not known precisely. A more reliable measurement of the P N is achieved through nuclear magnetic resonance (NMR) spin thermometry (see Supplementary Note 5 for details). The method rests on the assumption that the probability p m for each nucleus to occupy a state with spin projections m follows the Boltzmann distribution p m ∝ e mβ , where β = hν L /k b T N is the dimensionless inverse spin temperature, expressed in terms of nuclear spin Larmor frequency ν L and spin temperature T N (h is the Planck's constant and k b is the Boltzmann constant). For spin I=1/2 where m = ±1/2, any statistical distribution has the Boltzmann form.
By contrast, for I >1/2 Boltzmann distribution expresses the non-trivial nuclear spin temperature hypothesis [26], which was verified for epitaxial GaAs quantum dots previously [8].
In order to perform spin thermometry, we first measure the single-QD NMR spectra [27], as exemplified in Fig. 3a for 69 Ga spins. The three magnetic-dipole transitions of the 3/2 spins are well resolved thanks to the quadrupolar shifts ν Q , which originate from the lattice mismatch of GaAs and AlGaAs. On the other hand, these quadrupolar effects are too small to impede nuclear spin cooling -a significant advantage over the highly-strained Stranski-Krastanov QDs, where quadrupolar shifts are large and disordered [27]. The resolved NMR triplet is essential, as it allows β to be derived from the Boltzmann exponent, which then relates to P N through the standard in QD1 at B z = 10 T using "inverse NMR" signal enhancement technique [27]. Inset shows energy levels of a spin-3/2 nucleus. The resonant frequency of the central transition between m = ±1/2 is ν L , whereas the satellite transitions involving m = ±3/2 are split off by the quadrupolar shifts ±ν Q . b Low-resolution spectrum of the same QD1, but measured using the saturation technique in order to reveal the population probabilities of the nuclear spin levels. c Populations of spin levels with different projection m = ±1/2, ±3/2, sketched for the two different levels of optical nuclear spin polarization corresponding to the data in (b).
d Hyperfine shift variation arising from selective manipulation of the 69 Ga nuclear spin plotted against the photoluminescence spectral splitting ∆E PL in measurements where the degree of optical nuclear spin pumping is varied. Squares show the total 69 Ga hyperfine shift measured by broadband saturation of the entire NMR triplet, which equalizes populations p m for all m. Circles and triangles show the selective signals of the ±1/2 ↔ ±3/2 resonances measured via frequency-swept adiabatic inversion. Lines show fitting, from which nuclear spin polarization degree is derived and plotted in the top horizontal scale (see Supplementary Note 5). e Maximum positive and minimum negative nuclear spin polarization degrees P N derived for 69 Ga (triangles) and 75 As (circles) in individual dots QD1 -QD3 at B z = 10 T (solid symbols and QD numbers) and B z = 4 T (open symbols). The nonlinear scale ∝ 1 − P 2 N is used to highlight the areas around |P N | ≈ 1. Error bars are 95% confidence intervals. f Maximum positive and minimum negative hyperfine shifts measured on individual dots QD1 -QD12 at B z = 10 T. Fig. 3b with simple saturation NMR spectroscopy [28]. At moderate polarization P N ≈ −0.6 (dashed lines) all three magnetic-dipole transitions m ↔ m + 1 are observed, and their amplitudes are proportional to the differences |p m+1 − p m | (Fig. 3c). At the maximum positive polarization (solid line) a single NMR peak +1/2 ↔ +3/2 is observed, indicating that nearly all spins have been cooled to the m = +3/2 state. By changing the helicity of the optical pump it is possible to cool the nuclei towards the

Brillouin function. Qualitatively this is demonstrated in
For quantitative spin thermometry we measure the peak areas of the −3/2 ↔ −1/2 and +1/2 ↔ +3/2 NMR transitions at different initial polarizations of the nuclei, quantified by ∆E PL . The results are shown in Fig. 3d (circles and triangles), together with the total signal obtained by saturating all three NMR transitions (squares). Fitting with Boltzmann model is shown by the lines, together with the derived polarization degree P N in the top axis. The model reproduces well both the linear dependence of the total NMR signal and the non-linear dependencies of the selective ±1/2 ↔ ±3/2 signals, revealing a close approach to P N ≈ −1. Qualitatively, at P N = −1 the m = +1/2, +3/2 states must be depopulated, resulting in a vanishing +1/2 ↔ +3/2 signal, as indeed observed experimentally. Moreover, at P N = −1 the −3/2 ↔ −1/2 signal must be 2/3 of the total NMR signal, also in good agreement with experiment.
The largest positive and negative P N derived from spin thermometry on individual dots QD1 -QD3, chosen for their highest |E hf |, are shown in Fig. 3e. At the highest static field B z = 10 T the best fit estimates for 69 Ga are around |P N | ≈ 0.99, with somewhat lower |P N | ≈ 0.98 for 75 As. Spin thermometry conducted at B z = 4 T for one of the QDs also yields high polarizations, although the measurement accuracy is reduced due to the less efficient optical probing.
A simpler measurement of the largest positive and negative E hf is shown in Fig. 3f for 12 randomly chosen dots. For some QDs, nuclear polarization is reduced to P N ≈ 0.9. We also observe for all studied QDs that the optimal optical polarization of the pump is not circular, having a randomly-oriented linearly-polarized contribution ranging between 0 and 0.4. This points to in-plane anisotropy of QDs. From a control measurement, with magnetic field tilted by ≈ 12 • away from the growth axis, we find a reduction in maximum |E hf |. A reduction is also found in a piece of the same QD structure subject to a uniaxial stress in the sample plane. Therefore, the reduced P N in some QDs is attribute to the random anisotropy of the confining potential or strain. Low-symmetry confinement is known to result in heavy-light hole mixing [29,30], which may explain why optimal electron and nuclear spin pumping requires elliptically-polarized light.
The buildup dynamics, measured under optimal nuclear spin pumping, are shown in Fig. 4a. b Nuclear spin relaxation dynamics in the dark measured in a neutral state following σ + optical pumping (squares). The same relaxation dynamics are also measured with partial NMR saturation after the optical pump, which reduces the initial P N (triangles and stars). Lines show fitting used to derive the nuclear spin half-lifetimes T 1,N . Inset shows the same data, but normalized by the hyperfine shift at short dark times.
c Nuclear spin relaxation times T 1,N as a function of the initial hyperfine shift E hf . The corresponding approximate initial P N is shown on the top axis. Error bars are 95% confidence intervals. d Density of states calculated for N = 6 dipolar-coupled I = 3/2 nuclei (without the electron). Each band, broadened by dipolar couplings hν dd ∝ max |b j,k |, corresponds to a well-defined total spin projection M . The adjacent bands are split by the Zeeman energy hν L . e Population probability of the eigenstates, calculated for the spectrum in (d) and for two types of mixed states: Boltzmann distribution of Zeeman energies with high polarization P N ≈ 0.99 (red) and a narrowed Gaussian distribution with P N = 0 (blue).
The approach to the steady state is non-exponential since the nuclei that are further away from the center of the QD are less coupled to the electron and take longer to polarize. It takes on the order of ≈ 60 s to reach the steady-state P N within the measurement accuracy. Once optical pumping is switched off, nuclear spins depolarize in the dark (squares in Fig. 4b) on a timescale of minutes, mainly through spin diffusion [31]. Such long lifetimes mean that a highly-polarized nuclear spin state can be prepared and used for the subsequent fast (nanosecond) control of the electron spin qubit. We further examine the effect of the initial P N on the relaxation dynamics by augmenting the optically-pumped nuclear state with a short partially-depolarizing NMR pulse (triangles and stars in Fig. 4b). When normalized by the initial polarization, the plot reveals accelerated nuclear spin relaxation under reduced initial polarization (inset in Fig. 4b). This is quantified in Fig. 4c, where at high polarization the nuclear spin relaxation half-lifetime T 1,N is seen to be a factor of ≈ 2 − 3 larger than in case of low initial polarization (lowest initial polarization is limited by the accuracy of the T 1,N measurement). This is a non-trivial result: the spin diffusion model, as well as non-diffusion relaxation mechanisms are linear, so that scaling of initial P N should not change In order to explain the non-linear relaxation, we consider the eigenspectrum of a nuclear spin ensemble, with an example shown in Fig. 4d for N = 6 spins I = 3/2. The adjacent bands, separated by the Zeeman energy, typically hν L ≈ 1 − 100 MHz, correspond to a flip of a single nucleus, which changes the total ensemble spin projection M by ±1. Each band consists of all possible superpositions with a given M , with degeneracy lifted by the nuclear-nuclear dipolar magnetic interaction. For M ≈ 0 (i.e. P N ≈ 0) the broadening of each band is maximal, characterized by the dipole-dipole energy hν dd ≈ 1 kHz. In the opposite limit, there are only two non-degenerate fully-polarized states with M = ±N I (i.e. P N = ±1). Thus at |P N | → 1, the distribution of the available dipolar energies is narrower than at P N ≈ 0. The dipolar reservoir can act as a source or sink of energy for a flip-flop spin exchange between two nuclei whose energy gaps are slightly different (for example due to the inhomogeneity of the quadrupolar shifts ν Q ). Therefore, the slowdown of nuclear spin diffusion, which proceeds through pairwise nuclear flip-flops, is interpreted as a witness of dipolar reservoir narrowing at high |P N |.
The aforementioned narrowing of the nuclear dipolar reservoir is conceptually similar to the state-narrowing technique, which aims to reduce the statistical dispersion of the nuclear Zeeman energies ∝ M in order to enhance the coherence of the electron spin qubit. An example of a narrowed mixed state is sketched in Fig. 4e for P N ≈ 0, but with uncertainty in M reduced down to a few units, as demonstrated experimentally previously [32,33]. The fundamental advantage of a polarized state (also sketched in Fig. 4e), is that it both narrows the uncertainty in M by a factor ∝ 1 − P 2 N and reduces the dipolar broadening. The ultimate limit of P N = ±1 is the only case where electron spin qubit coherence is predicted to be essentially non-decaying [12,13].
By contrast, even if the dispersion of M is reduced to zero, the dipolar energy uncertainty of a depolarized ensemble may still cause dynamics on the timescales of 1/ν dd ≈ 1 ms, leading in turn to electron spin qubit decoherence. Evaluation of electron spin coherence in a highly-polarized nuclear spin environment is an interesting subject for future work and may also provide a more sensitive tool for nuclear spin thermometry near |P N | ≈ 1. Alternatively, more accurate measurement of P N can be sought through nuclear-nuclear interactions and the "trigger" detection method designed for dilute spins [26] but applied to the few abundant nuclei occupying the thermally excited spin states.
The nuclear spin cooling method reported here is applicable to a standard p − i − n diode structure, fully compatible with high-quality electron spin qubit operation, as demonstrated recently in the same semiconductor structure [34]. The technique is simple to implement and robust -once optical pumping parameters are optimized for a certain QD, they do not require any correction over months of experiments. Even larger nuclear polarizations can be sought by combining QDs of high in-plane symmetry with biaxial strain in order to reduce the heavy-light hole mixing. Our nuclear spin cooling method uses the purity of the optical pump polarization as the final heat sink, ultimately limiting the achievable P N even for a perfectly symmetric QD. This is different from the resonant "dragging" schemes [35][36][37] where the ultimate heat sink is the photon number in the optical mode, offering in principle a much closer approach to |P N | ≈ 1, provided the dark-state bottleneck could be avoided. Combining the advantages of the two approaches in a two-stage cooling cycle can be a route towards the ultimate goal of initializing a nuclear spin ensemble into its fully-polarized quantum ground state. This would be a prerequisite for turning the enormously large Hilbert space of the N ≈ 10 5 QD nuclei into a high-capacity quantum information resource.   [38][39][40]. The n-type doped layer is followed by the electron tunnel barrier layers: first a 5 nm thick Al 0.15 Ga 0.85 As layer is grown at a reduced temperature of 560 • C to suppress Si segregation, followed by a 10 nm thick Al 0.15 Ga 0.85 As and then a 15 nm thick Al 0.33 Ga 0.67 As layer grown at 600 • C. Aluminium droplets are grown on the surface of the Al 0.33 Ga 0.67 As layer and are used to etch the nanoholes [41,42]. Atomic force microscopy shows that typical nanoholes have a depth of ≈ 6.5 nm and are ≈ 70 nm in diameter.
Next, a 2.1 nm thick layer of GaAs is grown to form QDs by infilling the nanoholes as well as to form the quantum well (QW) layer. Thus, the maximum height of the QDs in the growth z direction is ≈ 9 nm. The GaAs layer is followed by a 268 nm thick Al 0.33 Ga 0.67 As barrier layer.
Finally, the p-type contact layers doped with C are grown: a 65 nm thick layer of Al 0.15 Ga 0.85 As with a 5 × 10 18 cm −3 doping concentration, followed by a 5 nm thick layer of Al 0.15 Ga 0.85 As with a 9 × 10 18 cm −3 concentration, and a 10 nm thick layer of GaAs with a 9 × 10 18 cm −3 concentration. The band structure of the electrons and holes in a GaAs QD is sketched in Supplementary Fig. 6 (see for example Ref. [21] for a review). The electron conduction band in GaAs has spin s = 1/2, with two possible spin projections s z = ±1/2 along the quantizing magnetic field. The valence band is four-fold degenerate at the center of the Brillouin zone in bulk GaAs. The confinement along the z growth axis is sufficient to split the valence band into the heavy hole and light hole subbands with total momentum projections j z = ±3/2 and j z = ±1/2, respectively. The typical heavy-light hole splitting is ∆E hh−lh ≈ 10 − 15 meV in GaAs/AlGaAs quantum wells [43,44]. The   Apart from the quantum-well type of confinement along the z axis, carriers in a QD are also confined in the orthogonal xy plane. In a real semiconductor structure there is always some breaking of the symmetry in the xy plane. Such in-plane anisotropy can mix the heavy and light holes, so that the eigenstates are no longer described by pure j z = ±3/2 or j z = ±1/2 projections.
As a result of such mixing the selection rules change, and the optical transitions in general become elliptically polarized.
In the pump-probe experiments we use photoluminescence of a negatively charged trion X − , where the electron-hole recombination occurs in presence of another resident electron. The X − spectra, such as shown in Supplementary Fig. 6b tend to have the narrowest linewidths and their Zeeman splittings are free from the non-linearity which is presented for neutral excitons X 0 due to the fine structure splitting. The energies of the states involved in X − photoluminescence are shown in Supplementary Fig. 6c. The energy of the ground state resident electron is (µ B g e B z + E e hf )s z , whereas the energy of the optically excited trion is E Gap +(µ B g h B z +E h hf )j z /3. Taking the differences and substituting the momentum projections allowed by the selection rules s z + j z = ±1, we find the photon energies of the two optically-allowed transitions is the electron (heavy hole) hyperfine shift and E Gap is the difference between the conduction band and valence band ground states in a QD. The splitting of the spectral doublet is then (S1) Next we eliminate the Zeeman contribution and define the excitonic hyperfine shift: where ∆E PL,0 is the photoluminescence doublet splitting at zero nuclear spin polarization. Valence band hole hyperfine interaction is of the order of 10% of the electron hyperfine interaction [45].
Consequently, the excitonic hyperfine shift E hf is dominated by the electronic contribution E e hf . The Hamiltonian describing the nuclear spin system alone includes the Zeeman, the quadrupolar and the dipole-dipole terms. The Zeeman term accounts for the coupling of the QD nuclear spins I j to the static magnetic field B z directed along the z axis: where the summation goes over all individual nuclei 1 ≤ j ≤ N , = h/(2π) is the reduced Planck's constant, γ j is the gyromagnetic ratio of the j-th nuclear spin andÎ j is a vector of spin operators with Cartesian components (Î x,j ,Î y,j ,Î z,j ). The result of the Zeeman term alone is a spectrum of equidistant single-spin eigenenergies m γ j B z , corresponding to 2I + 1 eigenstates witĥ The interaction of the nuclear electric quadrupolar moment with the electric field gradients is described by the term (Ch. 10 in Ref. [46]): where q j and η j describe the magnitude and asymmetry of the electric field gradient tensor, whose principal axes are x y z . The strain is inhomogeneous within the QD volume, so that q j and η j vary between the individual nuclei. The axes x y z are different for each nucleus and generally do not coincide with crystallographic axes or magnetic field direction. In lattice-matched GaAs/AlGaAs QDs the electric field gradients at the nuclear sites do not exceed q j /h ≈ 100 kHz, as witnessed via NMR spectroscopy. At sufficiently strong magnetic fields | γ j B z | |q j |, quadrupolar effects can be treated perturbatively -the main effect is the anharmonicity of the nuclear spin energies and the resulting quadrupolar NMR multiplet of 2I magnetic-dipole transitions, split by ν Q ≈ q j /h.
The m = ±1/2 states of a half-integer nuclear spin are influenced by quadrupolar effects only in the second order, resulting in a smaller inhomogeneous broadening, which scales as ∝ ν 2 Q /ν L with nuclear spin Larmor frequency ν L = γB z /(2π).
Direct interaction between the nuclei is described by the dipole-dipole Hamiltonian: Here, µ 0 = 4π×10 −7 NA −2 is the magnetic constant and r j,k denotes the length of the vector, which forms an angle θ with the z axis and connects the two spins j and k. The interaction of the conduction band electron spin s with the ensemble of the QD nuclear spins is dominated by the contact (Fermi) hyperfine interaction, with the following Hamiltonian: where the hyperfine constant of an individual nucleus j is a j = A (j) |ψ(r j )| 2 v . Unlike a j , the A (j) hyperfine constant is a parameter describing only the material and the isotope type to which nucleus j belongs, |ψ(r j )| 2 is the density of the electron envelope wavefunction at the nuclear site r j of the crystal lattice, and v is the crystal volume per one cation or one anion. The definitions of the hyperfine constants differ between different sources. With the definition adopted here, a fully polarized isotope with spin I, hyperfine constant A and a 100% abundance (e.g. 75 As), would shift the energies of the electron spin states s z = ±1/2 by ±AI/2, irrespective of the shape of |ψ(r j )| 2 .
For valence band holes the contact (Fermi) contribution vanishes, leaving the weaker dipoledipole terms to dominate the hyperfine interaction. Compared to the valence band electrons, the coupling has a more complicated non-Ising form [45]. The effect of the net nuclear polarization on the heavy-hole spin splitting can be captured by a simplified expression: whereĵ z is the z component of the hole spin momentum operator. The valence band hyperfine material constants C (j) are sensitive to heavy-light hole mixing and both their signs and magnitudes depend on the material [45].
Owing to the flip-flop term ∝ (ŝ xÎx,j +ŝ yÎy,j ) of the hyperfine Hamiltonian (Supplementary Eq. S6) the eigenstates of the electron-nuclear central spin system are in general entangled, i.e. they cannot be written as a direct product of the electron spin single-particle state and the nuclear spin ensemble state. Consequently, when such product state is generated through optical injection of a spin-polarized electron into the quantum dot, the wavefunction of the central spin system starts evolving. We estimate the rate of coherent evolution using the Rabi frequency ∝ j a 2 [25] for the limit of vanishing electron spin splitting. For a fully polarized nuclear spin ensemble coupled to an electron spin polarized in the opposite direction, this Rabi frequency describes the exact solution of periodic spin exchange between the electron and the collective nuclear spin state. Therefore, in order for dynamic nuclear spin polarization to be efficient, the polarized electron spins need to be removed and injected much faster than the hyperfine-induced Rabi rotations (otherwise the electron spin will periodically polarize and depolarize the nuclei, without any net spin transfer). For a typical GaAs QD with N ≈ 10 5 nuclear spins we have h √ N /A ≈ 25 ns. When electron spin splitting is not zero, there is an increase in the frequency of coherent oscillations that follow initialization into a product electron-nuclear state. In the limit of large electron spin splitting this is approximately the electron spin resonance frequency. For experimental conditions used in our work the maximum sum of the net hyperfine shift and the electron Zeeman splitting at B z = 10 T is within 180 µeV, which corresponds to electron Larmor period of 20 ps. From these basic derivations, we arrive to a rough estimate that electron spin recycling must occur on a sub-picosecond timescale in order to achieve near-unity nuclear spin polarization.

Supplementary Section 3. EXPERIMENTAL METHODS AND TECHNIQUES
All measurements are performed in a liquid helium bath cryostat. The sample is placed in an insert tube filled with a low-pressure heat-exchange helium gas. The base temperature is ≈ 4.25 K.
We use confocal microscopy configuration where QD photoluminescence (PL) is excited by a laser beam focused by a cryo-compatible apochromatic objective with a focal length of 2.89 mm and a numerical aperture of 0.81. The excitation spot diameter is ≈ 1 µm. Both the optical excitation and a static magnetic field B z up to 10 T are applied along the sample growth axis z (Faraday geometry).
Quantum dot photoluminescence is collected and collimated by the same cryo-compatible objective.
The PL signal is dispersed in a two-stage grating spectrometer, followed by a pair of achromatic doublets, which transfers the spectral image onto a charge-coupled device (CCD) detector with a magnification of 3.75. The orientation of the semiconductor sample is verified by reflecting a collimated laser off the sample surface -the small unintentional tilt of the sample is found to be ≈ 0.7 • . The laser used for optical pumping of the nuclear spins is a ring-cavity tunable titanium sapphire (Ti:Sa) laser, operating in a single-mode continuous-wave regime. This laser is coupled with a wavelength meter (30 MHz accuracy) for precise tuning and stabilization of the optical pumping wavelength. The sample gate bias is connected by a combination of a twisted pair (inside the cryostat) and a 50 Ω coaxial cable (outside the cryostat) to an arbitrary function generator through a low-pass LC filter with a 1.9 MHz cut-off frequency. Selective manipulation of the nuclear spins is achieved with a resonant radiofrequency oscillating magnetic field, generated by a small copper wire coil. This coil is placed to have its axis within the top surface of the semiconductor sample and perpendicular to the static magnetic field. A 50 Ω cryogenic coaxial cable is used to connect the coil to a radiofrequency amplifier with a maximum rated power of 100 W. Supplementary Fig. 7c shows a cycle used in the nuclear spin relaxation measurements. The cycles start with radiofrequency depolarization that is sufficiently long to eliminate any effect of the nuclear polarization left over from the previous measurement cycle. This is followed by optical pumping for T Pump = 60 s. After the pump, the sample is kept in the dark for a time x p e r i m e n t c y c l e P u m p P r o b e T Dark under gate bias V Dark . After that, a probe pulse is applied to measure the fraction of the nuclear spin polarization that decayed during T Dark . In this type of experiments the number of the pump-probe cycles used to collect the probe photoluminescence signal varies between 1 and 10 -single-shot probing is required when T Dark exceeds a few hundred seconds. Furthermore, we perform measurements using the cycle of Supplementary Fig. 7c, but with a second radiofrequency pulse added after the pump and before the dark interval. The duration of this second Rf pulse T Rf,2 is varied between 0 and 5 s to control the degree of the initial nuclear spin polarization. In principle, there are multiple ways to control the degree of the initial nuclear spin polarization in the quantum dot, such as the duration T Pump of the pump or its power and wavelength. However, any such changes in the optical pumping also affect the rate of nuclear spin diffusion into the barriers around the quantum dot [31]. The degree of nuclear polarization in the barriers then affects the rate of nuclear spin relaxation in the subsequent dark interval. The advantage of the second radiofrequency pulse is that it depolarizes the nuclei at the same rate in the entire sample.
Therefore, the spatial profile of the nuclear spin polarization P N after the second pulse is simply a scaled profile of the P N profile produced by the optical pulse. Spin diffusion is described by a linear differential equation, so that proportional reduction of P N in the entire sample should not affect the timescales of the subsequent nuclear spin diffusion and relaxation in the dark. Consequently, any dependence of the relaxation time on the degree of the initial nuclear spin polarization (left after the second radiofrequency pulse) is ascribed purely to the reduction (narrowing) of the energy that the dipole-dipole reservoir can supply or absorb during the nuclear flip-flop events of the spin diffusion process. We note that the second radiofrequency pulse has minimal effect on the measurement of the subsequent relaxation dynamics, since its duration is no more than a factor of 0.1 of the shortest measured nuclear spin relaxation time T 1,N .
Supplementary Fig. 7d shows the timing of the experiment used to study the dependence of the steady-state nuclear spin polarization on the optical pumping parameters such as sample bias V Pump , pump power and wavelength. The experiment cycle starts with a radiofrequency erase that eliminates any leftover nuclear polarization. Then the pump and probe pulses start, but the acquisition (CCD detector exposure) of the probe photoluminescence begins only after a delay T Buildup = 50 s. This delay allows nuclear spin polarization to build up closer towards its steady state so that relatively short pump pulses T Pump = 5 s can be used for faster acquisition of the photoluminescence signal. The probe pulses are kept short, in order to produce minimal nuclear spin depolarization during each pump-probe cycle (see details in Supplementary Section 3 C). The probe pulses are also much shorter than the pump T Probe /T Pump < 0.003, to ensure minimal effect on the steady-state nuclear spin polarization.

B. Optical pumping of quantum dot nuclear spins
The steady-state nuclear spin polarization depends on the wavelength of the pump laser. Maximum hyperfine shifts |E hf | are found to occur when the laser is resonant with a certain optical Inset shows a zoomed in dependence for the s-shell peak measured at slightly different transition of the quantum dot. Calibration of the optimal pumping parameters starts with a measurement of a broad-range wavelength dependence -an example is shown in Supplementary Fig. 8.
Once the individual spectral features, such as s−, p−, d− and f −shell peaks, are identified, a more detailed optimization is performed. We focus on the s−shell pumping peak and measure more detailed dependencies on the wavelength (or equivalently the pump photon energy) at different values of pump power P Pump and sample gate bias V Pump . The insert in Supplementary Fig. 8 shows an example of such a detailed dependence at the optimum P Pump = 2.7 mW, V Pump = −2.7 V. It can be seen that the pump laser needs to be tuned to within a narrow margin of ≈ 0.2 meV in order to achieve the highest possible |E hf |.
In the experimental setup the collimated pump laser beam first passes through a linear polarizer and then a λ/2 waveplate installed in a motorized rotation mount. This way it is possible to create arbitrary orientation of the linearly polarized beam, which is then directed to a cube beamsplitter, followed by a λ/4 waveplate installed in another motorized rotation mount. By placing the λ/4 waveplate last, it is possible to compensate for any polarization imperfections of the nominally nonpolarizing beamsplitter and obtain a beam with high degree of circular polarization. This beam is then directed through a quartz window of the cryostat insert and the cryogenic objective, which focuses it on the surface of the QD semiconductor sample. In order to account for any polarization imperfections in the optical path we perform calibration measurements where both the λ/2 and λ/4 waveplate orientations are scanned and the resulting hyperfine shifts E hf are measured. The results shown in Supplementary Fig. 9, indicate that the waveplate must be set within ±2 • in order to attain the highest nuclear spin polarization degree. The optimal orientations of the λ/2 and λ/4 waveplates are different for the minimum negative (triangles) and the maximum positive (squares) Once the orientations of the λ/2 and λ/4 waveplates are optimised, we examine the polarization state of the pump beam directed to the cryostat. To this end, we place a linear polarizer (analyzer) after the λ/4 waveplate, followed by a power meter. The linear polarizer is rotated to find the minimum (I min ) and the maximum (I max ) intensities of the transmitted beam. For a perfect circularly polarized beam I min = I max , whereas for a linearly polarized beam I min = 0. We characterise the optimized beams using the degree of linear polarization ρ lin = (I max − I min )/I max and the analyzer orientation angle α max where the maximum intensity is achieved. These results are summarized in the polar plot of Supplementary Fig. 10. The optimal degree ρ lin and orientation α max of the linearly polarized components vary between individual quantum dots. Moreover, optimal polarization parameters are different for σ + and σ − pumping and even depend on magnetic field for the same QD1. Such variability, as well as the large values of ρ lin 0.4 suggest that polarization imperfections in the optical elements (e.g. mechanical stress in the cryo-objective) are not the major contribution. The large deviation of the optimal optical pumping from pure circular polarization is therefore attributed to the properties of the individual quantum dots. These may include anisotropy of the QD shape and inhomogeneous microstrains that make semiconductor material around the QD act as an optical waveplate and give rise to heavy-light hole mixing [47] that causes optical selection rules to depart from those of the bulk GaAs. Indeed, previous studies have shown that a sufficiently large uniaxial strain 0.5 % can flip the valence band hole quantization axis into the sample growth plane [48]. Another measure of the QD anisotropy is the fine structure splitting (FSS) of a neutral exciton at zero magnetic field.
While we have not conducted systematic correlation studies, selective measurement on QD1, where very large nuclear spin polarization was achieved, revealed a FSS of ≈ 28 µeV. This is considerably larger than the few-µeV FSS observed in symmetric GaAs QDs [49]. This comparison suggests  For brevity, throughout this work we use the term "σ + pumping" ("σ − pumping") to describe the optimized elliptically-polarized optical pumping with σ + (σ − ) character of the circularly polarized component.

C. Optical probing of quantum dot nuclear spins
For optical probing of the nuclear spin polarization we use a diode laser emitting at 690 nm.
Sample forward bias, typically +0.7 V, and the probe power are chosen to maximize (nearly saturate) PL intensity of the ground state X − trion. The difference between the spectral splitting ∆E PL of the X − trion doublet and the same splitting ∆E PL,0 measured for depolarized nuclei reveals the hyperfine shifts E hf = −(∆E PL − ∆E PL,0 ). Illumination with a probe laser inevitably acts back on the nuclear spin polarization. In order to quantify such back-action we perform calibration measurements with examples shown in Supplementary Fig. 11. In these experiments the QD is first pumped with a σ + or σ − polarized laser with power and bias set to maximize the steady state nuclear polarization. Then the pump is switched off and the probe laser pulse is applied. The hyperfine shift E hf is measured from PL spectroscopy at the end of this probe.
It can be seen that the probe induces decay of the nuclear spin polarization. For QD1 we use the same probe power of P Probe = 30 nW, but the unwanted probe-induced depolarization is faster at B z = 4 T (solid symbols) compared to B z = 10 T (open symbols). For selective-NMR measurements of the nuclear spin polarization (spin thermometry) in QD1 we use T Probe = 15 ms at B z = 10 T, so that the resulting depolarization is negligible (< 1%). PL intensity of the same QD1 is weaker at B z = 4 T so we use a longer T Probe = 24 ms in order to obtain a sufficiently strong probe PL signal. However, this leads to a larger depolarization of ≈ 4% under σ − pumping.
Depolarization itself is not an issue, since it would simply rescale all the measured E hf , which would not affect the differential NMR spin thermometry. In practice, the probe-induced depolarization also depends on the instantaneous E hf -such nonlinearity is what causes the distortion, resulting in larger uncertainties of the nuclear spin polarization measured at B z = 4 T. For QD2 we use a lower probe power P Probe = 7 nW (crossed symbols in Supplementary Fig. 11), which leads to an even slower probe-induced depolarization than for QD1. This allows to have a longer probe pulse (T Probe = 40 ms) for QD2, while keeping parasitic depolarization small (< 1%).

D. Radiofrequency control of nuclear spins
The radiofrequency oscillating magnetic field B x ⊥ z is produced by a coil placed at a distance of ≈ 0.5 mm from the QD sample. The coil is made of 10 turns of a 0.1 mm diameter enameled copper wire wound on a ≈ 0.4 mm diameter spool in 5 layers, with 2 turns in each layer. Two main types of radiofrequency signals are used in this work. The first type is a frequency-swept monochromatic excitation which is used for adiabatic inversion of the nuclear spin population.
The amplitude of the radiofrequency field is constant and the frequency is swept linearly in time.
Radiofrequency sweeps are discussed further in Supplementary Section 5 D.
The second type is the broadband radiofrequency excitation which is required to saturate inhomogeneously broadened quadrupolar resonances. The typical width of the resonances that needs to be saturated is tens to hundreds of kHz (further details are given in Supplementary Section 5 B), which is significantly larger than the typical homogeneous NMR linewidth (< 1 kHz). Therefore monochromatic radiofrequency excitation cannot provide a sufficiently uniform saturation of the entire inhomogeneously broadened resonance. This necessitates the use of a broadband radiofrequency excitation. Ideally, one wants a signal with a rectangular spectral profile, that has a constant spectral density in the required frequency interval, and a zero intensity outside that interval. In practice, when implementing the radiofrequency waveforms on a digital generator, it is convenient to approximate the required rectangular spectral band with a frequency comb.
In spectral domain, the comb consists of periodically spaced monochromatic modes of constant amplitude, covering the desired interval of frequencies. The mode spacing of 120 Hz is chosen to be smaller than the homogeneous NMR linewidth. Under these conditions, by using a sufficiently small amplitude of each mode we achieve exponential depolarization (i.e. without nuclear spin Rabi oscillations) of the nuclear spin ensemble [50] with a typical time constant of τ ≈ 30 ms. The saturation of a chosen NMR resonance is achieved by applying a frequency comb excitation for a period of ≈ 5τ . When subject to such excitation, the nuclear spins undergo slow Rabi rotation, transitioning between the spin states parallel and antiparallel to the external magnetic field [28].
Due to the nuclear-nuclear dipole-dipole interactions each nuclear spin is subject to a local field.
The randomness of these local fields results in dephasing between Rabi precessions of the individual nuclei. Consequently, the nuclear spin ensemble becomes depolarized (i.e. each nucleus is randomly polarized) after a long saturation pulse.
For the saturation NMR spectra, shown in Fig. 3b of the main text, we use a frequency comb with a total width of 6 kHz. The central frequency of the comb is scanned to obtain the spectra -this frequency is the horizontal axis of the spectral plots. For the high-resolution NMR spectra, shown in Fig. 3a of the main text, we employ the "inverse" NMR technique [27] which enhances the NMR signal and allows the spectra to be measured even on those nuclear spin transitions that are depopulated at high polarization degrees. In this approach the radiofrequency excitation spectrum is a broadband frequency comb with a narrow gap. The central frequency of the gap is scanned and is used for the horizontal axis of the "inverse" NMR spectra. The width of the gap controls the balance between the NMR signal amplitude and the spectral resolution. For the spectra of

Supplementary Section 4. ADDITIONAL EXPERIMENTAL DATA
A. Extended data from nuclear spin pumping measurements Our approach to maximizing the nuclear spin polarization is through line-search optimization of the optical pumping parameters, such as pump photon energy E Pump , pump power P Pump , optical polarization and the sample bias V Pump . In order to understand the physics of the nuclear spin pumping process we also measure a systematic parametric dependence. Given the typical timescales of the nuclear spin process (1 − 100 s) it is not possible to explore the entire parameter space within a reasonable experimental time. Therefore, we measure various one-and two-dimensional sections in the multidimensional parameter space. An example is shown in Supplementary Fig. 8, where a one-dimensional dependence on E Pump is shown. Interpretation of the nuclear spin polarization data is conducted by correlating it with photoluminescence spectroscopy.
Supplementary Fig. 12a is a bias-dependent photoluminescence spectroscopy map measured under non-resonant optical excitation (HeNe laser emitting at 632.8 nm). At low excitation power P Exc = 0.5 nW we observe multiple spectral features, labelled accordingly. The emission of the bulk GaAs free exciton is observed at ≈ 1.518 eV, accompanied by a low energy band at ≈ 1.50 eV arising from doping and impurities. The emission of the Si-doped AlGaAs layer is observed as broad spectral features between ≈ 1.60 − 1.70 eV, and the peak at ≈ 1.71 eV originates from the GaAs quantum well (QW). Single-QD emission is observed between ≈ 1.56 − 1.59 eV as a series of Zeeman doublets that switch over as the gate bias V Gate is changed. The higher-resolution spectra of the QD excitons are shown in Supplementary Fig. 13. At V Gate ≈ +0.5 V photoluminescence is dominated by the neutral exciton X 0 , identified from its fine structure splitting at B z = 0 T. At more negative biases the emission of positively charged excitons dominates, since electrons rapidly tunnel out of the dot, leaving excess photogenerated (non-equilibrium) holes. At more positive biases the emission of X 0 is superseded by the negatively charged trion X − , which becomes dominant when QD confines a resident (equilibrium) electron. At even more positive biases the dot is charged with multiple resident electrons. The spectral features originating from doubly (X 2− ) and triply (X 3− ) charged excitons can be distinguished, while the photoluminescence peaks at even higher charge numbers tend to overlap.
When the power is increased (P Exc = 20 µW in Supplementary Fig. 12b)  intensity of the QW increases with power without saturation.
Supplementary Fig. 12c shows the two-dimensional map of the hyperfine shift E hf measured as a function of the pump power P Pump and photon energy E Pump at a fixed V Pump = +0.5 V, which roughly corresponds to a bias where the QD equilibrium state switches from 0 to 1 electron.
Nuclear spin pumping evidently takes place at powers as low as P Pump ≈ 10 nW (which is close to saturation power of the s-shell QD excitons) as long as the pump laser is tuned above the s-shell exciton transition, so that the QD can absorb the pump photons. However, the resulting nuclear spin polarization degree is low, characterized by E hf ≈ −30 µeV. Spin pumping efficiency increases when the pump power is increased to hundreds of µW, which is well above the ground state exciton saturation. The lowest negative E hf ≈ −90 µeV is achieved at P Pump ≈ 1 mW. At this high power, a series of spectral peaks is observed. Their periodicity matches the periodicity observed in the high-power photoluminescence spectra ( Supplementary Fig. 12b), which allows us to identify the peaks as originating from different excitonic shells (up to six visible). The mechanism of nuclear spin pumping can then be understood to arise from resonant absorption of the circularly polarized pump photons, which generate spin-polarized electrons and holes in the excited orbital states. It also follows from Supplementary Fig. 12c that at V Pump = +0.5 V the steady-state |E hf | produced by pumping via the higher d and f shells is larger than via the p and ground s shells. Excitation via higher shells means that excitons can relax towards the ground state before recombination. Such energy relaxation provides a route for a simultaneous exchange of spin with the nuclei [51], since it helps to absorb or supply a small amount of energy required to compensate the mismatch of the electron and nuclear Zeeman energies. Without the coupling to external energy reservoirs the electron-nuclear spin flip-flop would be energetically forbidden. At high powers P Pump 100 µW, nuclear spins can be polarized even via optical excitation well below the ground state QD exciton transition, indicating that it's a distinct spin pumping mechanism which we further discuss below.
Supplementary Fig. 12d shows the same dependence of E hf on P It is also worth noting that inverted E hf is observed under certain pumping conditions. For example, in Supplementary Fig. 12d we observe E hf > 0 around P Pump ≈ 10 µW and E Pump ≈ 1.580 eV, which is ≈ 16 meV above the energy of the s-shell nuclear spin pumping peak. This feature at E Pump ≈ 1.580 eV, labelled "LH", is ascribed to a light-hole exciton. Optical excitation of a heavy-hole exciton transition with a σ + polarized light (with photons carrying a +1 momentum) generates a hole with momentum projection j z = +3/2 and an electron with a s z = −1/2 spin projection. The s z = −1/2 electrons then lead to nuclear spin pumping with a negative E hf < 0, as indeed observed in Supplementary Fig. 12d for a wide range of the pump parameters. However, when resonant with a light-hole exciton transition, the same σ + photon generates a hole with j z = +1/2 and an electron with s z = +1/2 (see Supplementary Fig. 6a), which then leads to an inverted E hf > 0. The same argument applies to σ − optical excitation, and manifests in experiments as E hf < 0, observed under resonant excitation of the light-hole exciton transition.
In order to investigate the nuclear spin pumping mechanisms further, we fix the pump power at P Exc = 1.5 mW and measure the dependence of E hf on E Pump and V Pump . The results are presented in Supplementary Fig. 14c, which shows an extension of the data from Fig. 2c of the main text. The spectral peaks, ascribed to individual excitonic shells, are seen to Stark-shift with the applied bias ( Supplementary Fig. 14c). The largest hyperfine shift |E hf | is again observed for the s-shell exciton at reverse bias, which varies between −2.7 and −2.1 V for different individual QDs. The higher-shell excitons (p, d, etc.) differ from the s shell in that their excitation can be followed by relaxation into the lower energy shell(s). The slightly lower |E hf | can then be ascribed to such relaxation between shells, which may involve flipping of the electron without spin transfer to the nuclei. In other words, relaxation between shells may result in a reduced electron spin polarization, which in turn leads to reduction of the maximum achievable |E hf |. Resonant pumping into the s-shell can only generate up to two electrons and two holes in the QD. This leaves only four excitonic complexes that can take part in dynamic nuclear spin polarization: neutral exciton X 0 , neutral biexciton XX 0 , negatively charged trion X − and positively charged trion X + . We have performed the same measurements as in Supplementary Figs. 8 and 14c but with σ − pump polarization (producing positive E hf ). We find that the spectral positions E Pump of the optimal nuclear spin pumping peaks under σ + and σ − pumping are split by ≈ 900 µeV at B z = 10 T, matching excitonic spectral splitting observed in photoluminescence (Supplementary Fig. 13). However, it is not possible to determine directly which excitonic feature is responsible for nuclear spin pumping with maximum |E hf |. The biexciton XX 0 is unlikely to play a role -it consists of an electron spin singlet (two electrons, one with spin projection s z = −1/2 and one with s z = +1/2) and a hole spin singlet (two holes, one with momentum projection j z = −3/2 and one with j z = +3/2) which do not couple to nuclear spins. The X − trion is also unlikely to cause efficient nuclear spin pumping, because the two electrons are in a singlet state. Moreover, spin-selective optical excitation of X − requires prior injection of another spin-flipped electron, for which there is no sufficiently fast process that could compete with rapid tunneling. The X + trion is more likely to contribute to nuclear spin pumping, since it contains only one (spin-polarized) electron. However, to form a hole spin singlet, X + excitation would still need to be accompanied by hole spin flipping, which would create a bottleneck and slow down the cyclic nuclear spin pumping process. Therefore, we argue that resonant optical excitation of X 0 is the most likely route for efficient nuclear spin pumping, as it enables fast optical reexcitation upon tunneling of the previously-excited electron-hole pair out of the QD. In addition to the broad resonance that gives the most efficient nuclear spin pumping (E Pump ≈ 1.564 eV at V Pump = −2.3 V in Supplementary Fig. 14c), there are narrower and less efficient Stark-shifting resonances observed at lower E Pump (intersecting around E Pump ≈ 1.552 eV and E Pump ≈ 1.546 eV for the same V Pump = −2.3 V). These narrow resonances may correspond to optical excitation of X + and X − . This would be consistent with photoluminescence spectra ( Supplementary Fig. 14a), which show that all charged exciton transitions appear on the low-energy side of X 0 . Further investigation, both experimental and theoretical, would be needed to establish with certainty which excitonic transition is responsible for high-efficiency nuclear spin pumping in the regime of fast tunneling.
Supplementary Fig. 14c shows that nuclear spins can be polarized at photon energies down to E Pump ≈ 1.25 eV, which is well below the ground state QD exciton energy and bulk GaAs bandgap.
This mechanism leads to negative hyperfine shifts E hf ≈ −30 µeV for both σ + and σ − pumping, suggesting that optical excitation plays a different role, possibly related to activation of charge traps or Auger effect. The trapped-charge hypothesis is further supported by the bias dependence, which shows that the sub-bandgap nuclear spin pumping disappears for V Pump < −0.5 V and V Pump > +0.9 V. The buildup time of the nuclear spin polarization is found to be around ≈ 5 s, which is approximately an order of magnitude slower than nuclear spin pumping via resonant excitation of the QD excitons (see Supplementary Section 4 C). The exact mechanism of subbandgap optical nuclear spin pumping is currently unclear and would require a separate systematic investigation.
Focusing on the s-shell pumping, we plot the minimum E hf (i.e. maximum |E hf |) as a function of bias V Pump in Supplementary Fig. 15a. The corresponding photoluminescence intensity of the s-shell exciton under above-gap excitation is shown in Supplementary Fig. 15b. At positive V Gate the electric field in the sample is small and the band structure is close to flat-band (right sketch in Supplementary Fig. 15c). As a result, optical recombination is the only way the photo-generated carriers can leave the QD. The typical radiative lifetimes for the studied type of GaAs QDs are ≈ 300 ps, which creates a bottleneck for how quickly the spin-polarized electrons can be injected into the QD, limiting in turn the rate of the nuclear spin pumping. When the sample gate is tuned towards larger reverse (negative) bias, the maximum hyperfine shift |E hf | increases. At the same time, photoluminescence intensity gradually decreases when V Gate < −1 V, indicating that electrons and holes tunnel out of the QD (left sketch in Supplementary Fig. 15c) faster than they can recombine optically. We attribute this correlation to the key role that the tunneling plays in dynamic nuclear spin polarization. Fast tunneling allows to overcome the radiative-recombination bottleneck, so that high-power optical excitation can be used to inject spin polarized electrons at a high rate. The efficiency of nuclear spin pumping peaks at V Pump = −2.3 V. For even larger reverse bias (i.e. more negative V Pump ) the maximum hyperfine shift |E hf | is seen to reduce slightly. For V Pump < −2.3 V tunneling becomes even faster, which would require an even higher pump power P Pump 10 mW to maintain steady-state occupation of the QD with spin-polarized electrons.
However, when focused into a diffraction-limited spot, such high-power optical excitation causes heating of the crystal lattice, which may result in accelerated nuclear spin relaxation, explaining why the highest achievable nuclear spin polarization is reduced at very large reverse biases.

B. Estimate of the electron tunnel rate
In order to quantify the optical nuclear spin pumping process, we estimate the electron tunneling rate using the photoluminescence intensity data. We employ a rate equation approach by considering the probability p eh that the QD s-shell is occupied by an electron-hole pair. The (1 − p eh )Γ exc − p eh (Γ R + Γ Tun ) = 0 ⇒ The intensity of photoluminescence I PL is proportional to p eh . We further assume that the maximum observed PL intensity I PL,max corresponds to negligible tunneling rate Γ Tun → 0. With this assumption we eliminate the unknown PL intensity that would be observed at p eh = 1 and find for the tunneling rate at an arbitrary bias: It then follows that the reduction of photoluminescence intensity under reverse bias signifies that the tunneling rate exceeds the sum (Γ R + Γ exc ) of the radiative recombination and optical excitation rates. The radiative recombination time is 1/Γ R ≈ 300 ps (Ref. [24]) for the studied type of QDs (Γ R ≈ 3.3 × 10 9 s −1 ). Since photoluminescence intensity of the s shell exciton is saturated, we assume that the excitation rate exceeds the recombination rate Γ Exc > Γ R . This gives the lower bound estimate for the tunneling rate. The photoluminescence measurement shown in Sup- In order to estimate Γ Tun below V Gate < −2 V we consider the well-known WKB approximation of the tunneling rate through a triangular barrier (see e.g. Ref. [52]). Up to a constant factor, we have: where m * e is the effective electron mass, e is the electron ionization energy and F z is the electric field in the growth direction. The total thickness of the structure between the doped layers is ≈ 300 nm, so the electric field is estimated as We can now independently estimate the optical excitation rate that leads to optimal nuclear spin pumping. Resonance fluorescence intensity, measured on InGaAs/GaAs QDs in the same setup and under similar experimental conditions, saturates at P Exc ≈ 5 nW. We assume that Γ Exc ≈ Γ R at saturation, where the radiative rate is Γ R ≈ 10 9 s −1 for InGaAs QDs. Assuming that the excitation rate scales linearly with optical power, we find that the resonant pumping power of P Pump ≈ 1.5 mW, used for optimal nuclear spin pumping, corresponds to Γ Exc ≈ 1.5 mW 5 nW 10 9 s −1 ≈ 3 × 10 14 s −1 . This corresponds to optical reexcitation time of ≈ 0.0033 ps. These estimates yield Γ Exc that are comparable to or somewhat higher than the above-calculated Γ Tun , as would be Apart from the fast cycling of the optically-generated electrons, we expect that fast tunneling also facilitates the nuclear spin pumping by disrupting the formation of coherent nuclear "dark" states, which are otherwise predicted to prevent the approach to a near-unity nuclear spin polarization [10,11]. Moreover, the short (tunneling-limited) lifetime of the electron spin can also be interpreted as spectral broadening of the electron spin levels. Such spectral broadening can facilitate nuclear spin pumping by compensating the energy mismatch of the electron and nuclear spin Zeeman energies, which otherwise inhibits the electron-nuclear spin flip-flops. The role played by tunneling is then similar to the effect that elevated lattice temperatures have on nuclear spin polarization, as studied previously in InGaAs QDs [51]. Future theoretical work may explain the details of the nuclear spin pumping process by treating optical excitation, tunneling and electron-nuclear spin interactions in a unified framework.

C. Nuclear spin buildup dynamics
The nuclear spin buildup dynamics under optical pumping are non-exponential (Fig. 4a of the main text). Therefore we fit the data by a sum of two stretched exponentials: The best fits are shown by the solid lines in Fig. 4a of the main text. The fitting parameters for the optimal steady-state nuclear spin pumping of QD2 at B z = 10 T are as follows: Parameter σ + pumping σ − pumping The fit is empirical in nature, so its parameters should be treated as estimates. Nevertheless, we can establish that the fast initial buildup occurs on a 0.1 − 0.3 s timescale, slowing down to 0.5 − 1.5 s when the nuclear spin polarization approaches closer to its steady state. Electronnuclear spin dynamics become nonlinear when electron spin Zeeman splitting is cancelled by the hyperfine shift. This is manifested in a kink in the nuclear spin buildup dynamics, observed at Fig. 4a of the main text. From the zero-splitting condition E hf ≈ −µ B g e,z B z we can estimate the electron g-factor g e,z ≈ −0.09, in agreement with previous measurements on the same structure [31].
It is also interesting to estimate the rate of the electron-nuclear spin flip-flops. Starting from a where β = hν L /k b T N is the dimensionless inverse temperature, expressed in terms of the nuclear spin Larmor frequency ν L and the spin temperature T N (h is the Planck's constant and k b is the Boltzmann constant). In this definition we assume ν L > 0 and B z > 0, so that m = +I is the ground state for nuclei with γ > 0, in agreement with Supplementary Eq. S3. For spin I=1/2 where m = ±1/2 any statistical distribution is described by Supplementary Eq. S13 with some T N . By contrast, for I >1/2 Supplementary Eq. S13 states the non-trivial nuclear spin temperature hypothesis [26] -previous experimental studies on low-strain epitaxial quantum dots [8] have shown its validity for the state induced by optical dynamical nuclear polarization. Nuclear spin polarization degree is defined as For the Boltzmann distribution of Supplementary Eq. S13 the polarization degree is given by the Brillouin function: It is worth noting that the polarization degree P N and the dimensionless inverse temperature β provide more relevant description of the nuclear spin state than the temperature T N . Indeed, when the Larmor frequency ν L is changed (by varying the external magnetic field) P N and β are preserved, whereas T N is not constant for a given optically-pumped nuclear spin state. The temperature T N only gains physical meaning at low magnetic fields, comparable to the local nuclear dipolar fields.
P N and β are also related to entropy [4,54] -the minimum in entropy is achieved only for P N = ±1 The hyperfine shift experienced by the quantum dot exciton is linearly proportional to the nuclear spin polarization degree P N : where the sum is over individual isotope species. Although the hyperfine shift E hf can be measured accurately from the photoluminescence spectra, the proportionality factor F (j) depends not only on the material's hyperfine constants A (j) , but also on the leakage of the electron wavefunction into the AlGaAs barrier. Since it is difficult to estimate this leakage independently, the measurement of E hf alone is not suitable for accurate derivation of P N . The unknown proportionality factor between E hf and P N can be eliminated for I > 1/2 nuclei, provided that it is possible to address One can then substitute Supplementary Eq. S14 into Supplementary Eq. S16 to calculate the change in the optically detected hyperfine shift ∆E hf resulting from selective saturation of a single NMR transition m ↔ m + 1. For example, for m = +1/2 and m + 1 = +3/2 of the j-th isotope, we calculate ∆E . This result has a simple interpretation that the hyperfine shift variation ∆E hf depends only on the difference in the initial populations of the states that are selectively saturated with Rf.

B. Corrections for the nuclei with small or inverted quadrupolar shifts
For a fully resolved NMR triplet, Supplementary Eq. S17 is sufficient to extract the inverse temperatures β j and derive the polarization degree of the spin-3/2 nuclei. In a real semiconductor system the separation of the quadrupolar NMR components is not perfect. Here we examine the role that the nuclei with small or inverted quadrupolar shift ν Q have on the derivation of nuclear spin polarization from experimental data. The case of an experiment where a radiofrequency comb is used to saturate two out of three NMR transitions is considered in Supplementary Fig. 16. For unstrained GaAs, nuclear quadrupolar effects are absent (ν Q = 0) and all NMR transitions of the spin-3/2 nucleus appear at the same Larmor frequency ν L . Strain induces quadrupolar effects which are characterized to first order by the shift ν Q . In all our experiments ν Q ν L , so that first-order approximation is valid. The central transition −1/2 ↔ +1/2 between nuclear states with spin projections m = ±1/2 is unaffected by quadrupolar shifts in the first order, hence its NMR frequency is ν L for all nuclei (vertical solid line in Supplementary Fig. 16). The satellite transitions are affected by quadrupolar shifts: the NMR frequency of the −3/2 ↔ −1/2 transition is ν L + ν Q , whereas the NMR frequency of the +1/2 ↔ +3/2 transition is ν L − ν Q (solid lines in Supplementary Fig. 16 with slopes +1 and −1, respectively).
The strain varies within the QD volume, so there is a statistical distribution of ν Q values within the ensemble of the nuclei (sketched in the left part of Supplementary Fig. 16). The majority of the 75 As nuclei have a positive quadrupolar shift ν Q > 0 (for Ga nuclei the shift is predominantly negative ν Q < 0). Therefore, if we want to saturate simultaneously the two NMR transitions −3/2 ↔ −1/2 and −1/2 ↔ +1/2 (labelled −3/2 ↔ +1/2 for brevity) we choose a radiofrequency comb band sketched by the shaded area in Supplementary Fig. 16. As can be seen in Supplementary Fig. 16 the nominal −3/2 ↔ +1/2 comb saturates the desired transitions for the majority of nuclei, which have −ν Offs,L < ν Q < ν Offs,H (note that the lower bound −ν Offs,L is positive for a negative ν Offs,L < 0). These nuclei give a correct contribution to the Rf-induced hyperfine shifts. But there are also several cases, where nuclei give contributions that differ from those intended. For a fraction of nuclei with small quadrupolar shifts ν Offs,L < ν Q < −ν Offs,L all three NMR transitions are excited by the comb, resulting in full depolarization of such nuclei. Furthermore, for those 75 As nuclei where quadrupolar shift is inverted −ν Offs,H < ν Q < ν Offs,L , the −3/2 ↔ −1/2 transition will be out of resonance with the Rf band, while the +1/2 ↔ +3/2 transition will be saturated. Such nuclei will produce hyperfine shifts that would correspond to the −1/2 ↔ +3/2 two-transition saturation, rather than the intended −3/2 ↔ +1/2 saturation. Finally, for a small fraction of nuclei with very large absolute quadrupolar shifts |ν Q | > |ν Offs,H |, only the central −1/2 ↔ +1/2 transition will be saturated. Thus, introducing the empirical coefficients c, the observed hyperfine shift ∆E −3/2↔+1/2 hf,Obs in the −3/2 ↔ +1/2 two-transition saturation experiment can be written as: where we have dropped the isotope index. If the quadrupolar NMR triplet is fully resolved, then c Sat,Ideal = 1 with all other c coefficient equal to zero. In a real quantum dot c Sat,Ideal < 1 and the remaining coefficients are non-zero, so that the observed hyperfine shift ∆E The main difference of the single-transition selective excitation is that the nuclei with small quadrupolar shifts ν Q are eliminated from the measured hyperfine shifts. This is preferred over the two-transition saturation measurement, where such nuclei are fully depolarized, resulting in a parasitic hyperfine shift characterised by the c Sat,Full coefficient in Supplementary Eq. S18.
C. NMR spectra of the QD nuclei Supplementary Fig. 18a,b shows typical nuclear magnetic resonance spectra measured on 75 As and 69 Ga nuclei in QD1. The three magnetic dipole transitions of each of the spin-3/2 isotopes are split due to the natural elastic strain within the quantum dot volume, arising most likely from the residual lattice mismatch between the GaAs QD and the AlGaAs barriers. Although the NMR triplet is well resolved, there is a few-percent overlap between the spectral components.
When quantifying nuclear spin polarization degrees close to unity, such overlap must be taken into consideration. In order to quantify the spectral overlap, we study a piece of the same QD sample but subject to a uniaxial stress along the [110] crystallographic direction (i.e. the strain is applied perpendicular to the sample growth direction). Nuclear quadrupolar shifts induced by the external stress significantly exceed the intrinsic quadrupolar shifts. As a result the NMR triplet is fully resolved, as can be seen in the inverse NMR spectra of Supplementary Figs. 19a,b where we focus on the −1/2 ↔ +1/2 and −3/2 ↔ −1/2 transitions.
The spectral shapes of the −3/2 ↔ −1/2 satellites measured in a stressed QD sample are similar to those in the unstressed sample ( Supplementary Fig. 18). The satellite lineshape consists mainly of an asymmetric peak, but also shows evidence of spectral wings that are broad enough to overlap with the −1/2 ↔ +1/2 central transition in an unstressed sample. In principle, the overlap can be derived by integrating the relevant part of the satellite lineshape, measured with inverse NMR and shown in Supplementary Figs. 19a,b. However, this approach is vulnerable to noiseit is more efficient to incorporate integration into the NMR spectroscopy method [34,55]. Such an integral NMR measurement is performed by selectively saturating all nuclear spin transitions within a certain spectral band. The high-frequency edge of the saturating band (implemented as a frequency comb) is kept fixed. The low-frequency edge is scanned and the resulting change in the hyperfine shift ∆E hf is measured. This dependence of ∆E hf reveals the fraction of the nuclei covered by the saturating Rf band, and therefore provides a scaled definite integral of the NMR lineshape. For 75 As nuclei the fixed-frequency edge of the Rf band is detuned by +500 kHz from the central transition to ensure that the entire −3/2 ↔ −1/2 transition can be covered. For 69 Ga nuclei the fixed-frequency edge of the Rf band is detuned by +230 kHz from the central transition, so that the entire −1/2 ↔ +1/2 and +1/2 ↔ +3/2 transitions are also included in the band, in order to amplify the integral NMR signal of the −3/2 ↔ −1/2 satellite. Supplementary Fig. 19c shows the integral NMR spectrum of the −3/2 ↔ −1/2 transition of 75 As nuclei in a stressed sample. The steep rise in the integral signal matches the position of the sharp peak in the inverse NMR spectrum of Supplementary Fig. 19a. However, we also observe the slopes that stretch as far as ≈ ±100 kHz from the satellite peak maximum, indicating the contribution of a broad NMR signal. Broad spectral features are also observed in Supplementary   Fig. 19d for 69 Ga nuclei, though in a narrower spectral range and with an overall smaller contribution. This broad background can be ascribed to the NMR signal from AlGaAs barriers or any Al atoms diffusing into the GaAs QD layer [34]. The Al atoms that randomly replace the Ga atoms distort the tetrahedral symmetry of the four nearest neighbours surrounding each As atom.
The resulting unit-cell-scale strain results in pronounced quadrupolar shifts. By contrast, all Ga atoms have four identical As atoms as nearest neighbours. Therefore, Ga atoms are affected by Al/Ga random alloying only through next-nearest neighbours, explaining why the broad nuclear quadrupolar wings are smaller than for As nuclei.
Integral NMR spectra are processed in order to derive the correction coefficients. The experimental data is first smoothed (lines in Supplementary Fig. 19c,d) by fitting with a sum of three skew normal distribution peaks. The integral lineshapes are then normalized and shifted along the frequency scale to have the −3/2 ↔ −1/2 satellite NMR peaks in the stressed sample (Supplementary Fig. 19a,b) match the peak positions in the unstressed sample ( Supplementary Fig. 18a,b). It can be seen that the contributions of the ideal signals are higher for the 69 Ga nuclei due to their smaller quadrupolar broadening. As a result, nuclear spin polarization measurements are more accurate for 69 Ga than for 75 As. It is worth noting that the nuclear spin thermometry data measured on QDs in an unstrained sample is corrected with the c coefficients measured on a different individual QD (in a stressed sample). However, measurements conducted on several individual QDs from the same sample reveal NMR spectra very similar to those shown Supplementary Fig. 18a,b.
Thus, while there is always some uncertainty arising from dot-to-dot variation, its effect is expected to be smaller than the actual correction introduced through the c coefficients. magnitude of the hyperfine shift ∆E hf in the slow-sweep limit decreases, indicating that population transfer becomes non-adiabatic. This non-adiabaticity is a result of demagnetization in the rotating frame, where Zeeman energy is transferred into the nuclear dipole-dipole interaction reservoir [46,56,57]. For all ν 1 the sweep also becomes non-adiabatic in the large-rate limit.
Supplementary Fig. 20b shows sweep rate dependence for the range starting from −5 kHz to −50 kHz, which selectively covers the −3/2 ↔ −1/2 satellite NMR transition (the starting points of the sweeps are shown by the dashed lines in Supplementary Fig. 18b). selective Rf excitation of a certain NMR transition for a chosen isotope, whereas hyperfine shifts arising from other transitions and isotopes remain unaffected. In the experiment we avoid a certain range of positive initial nuclear spin polarizations (characterised by 850 µeV < ∆E PL < 900 µeV for QD1 at B z = 10 T) where electron spin energy splitting is close to zero due to the hyperfine shift and the Zeeman effect cancelling each other. Such cancellation is characterised by accelerated nuclear spin dynamics (see Fig. 4a of the main text), making it difficult to perform non-perturbing optical probing. Therefore, in order to discuss the spin thermometry fitting, we focus on the negative nuclear polarizations, where most of the datapoints are collected (Supplementary Fig. 22).
The splitting in the photoluminescence spectrum of a negatively charged trion X − (see Supple- Data for the non-selective saturation (same data as in (a), squares) and selective adiabatic frequency sweeps over NMR satellites (circles for −3/2 ↔ −1/2, triangles for +1/2 ↔ +3/2). c,d Same as (a,b) but for 75 As nuclei.
mentary Fig. 6b) can be written as (see Supplementary Eq. S2): where ∆E PL,0 is the trion splitting corresponding to depolarized nuclei and the summation goes over all isotopes with their individual polarization degrees P N,j . The individual proportionality constants can be written as F (j) = k j (A (j) − C (j) ), where both the electron (A (j) ) and the hole (C (j) ) hyperfine material constants are included since the photoluminescence of the trion is only observed for recombination of an electron and a hole with the opposite spin z projections (see Supplementary Section 2). The factors 0 < k j ≤ 1 account for the Ga nuclei atoms replaced by Al, resulting in a reduced hyperfine shift experienced by the electron spin. If all I = 3/2 isotopes have the same polarization degree, Supplementary Eq. S21 simplifies to Resolving this for P N and substituting into the last of Supplementary Eq. S17, we find that the change in the hyperfine shift arising from non-selective saturation of all 3 NMR transitions of the j-th isotope is a linear function of the trion spectral splitting: where w j is the weight coefficient of the j-th isotope in the total hyperfine shift E hf , and F tot is the total proportionality factor. The measured ∆E −I↔+I hf,j are shown by the squares in Supplementary   Fig. 22a and all isotopes to form the total χ 2 functional.
As a last step, we include in our model the possibility that different isotopes have different polarization degrees P N,j . For arbitrary P N,j it is not possible to resolve the trion PL splitting ∆E PL as a function of P N,j , requiring some explicit assumptions. As a simplest approximation, we assume that polarization degrees of the three abundant spin-3/2 isotopes ( 75 As, 69 Ga and 71 Ga) are linearly interdependent. Mathematically, this is equivalent to allowing the total scaling factor F tot to deviate for the different measured isotopes. The introduction of F tot and w j as model parameters is also convenient in that the data does not have to be collected on all isotopes, in particular on 71 Ga, which has not be studied in this work.
The χ 2 is a function of only five fitting parameters: the zero-polarization trion splitting ∆E PL,0 , the total scaling factors F tot of 75 As and 69 Ga as well as the weight coefficients w j of 75 As and 69 Ga. In case of QD1, where the data was measured at two different magnetic fields, we fit these datasets independently in order to account for the different degree to which the probe laser pulse introduces parasitic depolarization in the optically measured hyperfine shifts. The best fits obtained by minimizing χ 2 functional are plotted by the solid lines in Supplementary Fig. 22 and show a good match to the measured data. The best fit parameters are listed below together with the total number of experimental datapoints N data and the root-mean-square (RMS) residual R min = χ 2 min /N data derived from the minimized functional value χ 2 min . In addition, we quote the residual obtained from a separate linear fit where only the full-saturation hyperfine shift ∆E −I↔+I hf,j is considered: As discussed above, the full-saturation experiment is the most robust against the errors arising from NMR spectral overlaps. Therefore, the RMS residual R min obtained from linear fitting of ∆E −I↔+I hf,j alone characterizes the random measurement errors. These errors in the opticallydetected hyperfine shifts originate mainly from the noise of the CCD detector used to collect the optical photoluminescence spectra of a single quantum dot. Any excess of R min obtained from the nonlinear fit of the selective-NMR data is an indicator of systematic deviation from the fitting model. According to Supplementary Eq. S24 such excess is small, confirming the validity of the Boltzmann distribution model (Supplementary Eq. S13). The spread in the isotope-specific weights w i , and the scaling factor F tot is on the order of 1% for the B z = 10 T data collected from three individual quantum dots, affirming the systematic nature of these results. The fit of the B z = 4 T data shows the most deviation, which is explained by the need for shorter probe pulses T Probe , resulting in more noisy photoluminescence spectra as well as larger systematic deviations arising from the probe-induced nuclear spin depolarization.
In order to derive the nuclear spin polarization degree, we use the highest and the lowest trion spectral splitting ∆E PL detected without any radiofrequency manipulation. For this measurement we use the timing diagram of Supplementary Fig. 7d, where we allow the nuclear spin polarization to build up over T Buildup > 100 s, giving a closer approach to the steady state than what can be achieved in the NMR thermometry measurements ( Supplementary Fig. 22), where the pumping time T Pump 30 s is limited by the maximum duration of the CCD detector exposure. It is worth noting that this approach of using the separately measured steady-state ∆E PL is the reason why we build our fitting model on relating the polarization degree P N to spectral splitting ∆E PL via Supplementary Eq. S21. Otherwise, ∆E PL can be eliminated and polarization degree can be derived purely from the Rf-induced hyperfine shifts of Supplementary Eq. S17. The best fit value of ∆E PL,0 is subtracted from the steady state ∆E PL to derive the lowest negative and the highest positive hyperfine shifts. For QD1 from at B z = 10 T we find E hf = −109.6 µeV and E hf = +112.3 µeV.
The latter number exceeds by ≈ 1% the best-fit product IF tot ( 69 Ga). By definition, the IF tot product is the maximum |E hf | corresponding to full polarization P N = ±1. The discrepancy with the measured E hf reveals the scale of errors in the derived polarization degrees P N , both due to the random noise in the raw data and any systematic inaccuracy of the fitting model.

F. Error analysis in model fitting of the nuclear spin polarization data
In order to systematically analyze the fitting errors we construct a multidimensional confidence region (Chapter 9 in Ref. [58]) defined as a collection of all points in the fitting parameter space for which χ 2 < (1 + Q(γ, n)/N data )χ 2 min , where we have approximated the standard error in the experimental data by the RMS fit residual R min = χ 2 min /N data . (Note that here we define χ 2 and χ min as sums that are not normalized by the standard error.) We define Q(γ, n) as a quantile of the χ 2 -distribution with n parameters corresponding to the confidence level 1 − γ. We use 1 − γ = 0.95 where the relevant quantile is Q(1 − 0.95, 5) ≈ 11.07. We implement a Monte-Carlo calculation, where the χ 2 sum is computed for a large number of random sets of the fitting parameters around the best-fit point. For each trial point that satisfies Supplementary Eq. S25 we calculate the polarization degrees P N from the maximum and minimum steady-state spectral splitting ∆E PL measured with long T Buildup > 100 s.
This systematic evaluation agrees with the rough estimates above, confirming that the accuracy of our P N estimates is on the other of a few percent. The fit returns similar values for polarization degrees of 75 As and 69 Ga -this is expected for a spin pumping mechanism [8] where the inverse temperature β of each individual nucleus is independently equilibrated with the β of a spin-polarized electron. The somewhat broader confidence intervals of 75 As could be simply due to the larger overlaps of the NMR spectral components, making the fit less sensitive to P N and more dependent on the accuracy of the c coefficients tabulated in Supplementary Eq. S20.
In order to further evaluate the error estimates we approach the problem of data modelling from the opposite direction. Namely, we start with a hypothesis that the maximum absolute polarization degree |P N | is no more than a certain value < 1, and then evaluate how well our experimental data can be matched to this hypothesis. Constraining the fit to |P N | < 0.9 we search for the fitting parameter combination that minimizes the χ 2 functional for the same model as the one used to derive the best fit. The resulting constrained best-fit is shown for QD1 by the dashed lines in Supplementary Fig. 22. There are visible systematic deviations from the measured data, already suggesting that the |P N | < 0.9 hypothesis is inadequate, and the actual absolute polarization degree is well above 0.9. Quantitatively, the RMS residual from the fit of the B z = 10 T datasets constrained to |P N | < 0.9 is ≈ 1.73 µeV for QD1, ≈ 2.91 µeV for QD2 and ≈ 2.17 µeV for QD3.
These residuals are a factor of 2 higher than the best-fit values tabulated in Supplementary Eq. S24 -statistically, such deviations are improbable for our datasets containing hundreds of datapoints.
Next we fit the same experimental data but without correcting for the overlaps of the NMR spectral components. This is equivalent to setting c Sat,Ideal = c Swp,Ideal = 1 with the remaining c coefficients set to 0. Such a fit can be seen as a lower bound estimate for the absolute polarization degree |P N |. Without correction, the RMS fit residual slightly increases from 0.749 µeV to 0.819 µeV. The resulting uncorrected P N are shown in Supplementary Fig. 23. The uncorrected polarization degrees for 75 As (P N ≈ 0.88) are lower than for 69 Ga (P N ≈ 0.94), contradicting the expectation of equal β across different isotopes. Moreover, the maximum uncorrected |P N | are very close to the corresponding c Sat,Ideal coefficients. Additional computations confirm that this is to be expected -a naive uncorrected fit of the data affected by NMR spectral overlap returns for fully polarized nuclei a reduced polarization |P N | which roughly equals the fraction of the "ideal" nuclei that are not affected by the overlap. Nevertheless, even without the corrections, a high polarization degree is derived for 69 Ga nuclei, since they are less prone to NMR spectral overlaps than 75 As.
We now consider the different possible sources of systematic errors. Since electron localization in a GaAs quantum dot is not infinitely strong, the electron wavefunction leaks into the AlGaAs barriers where it gradually decays with the increasing distance from the dot. This means that both the nuclear spin pumping efficiency and the sensitivity of the electron hyperfine shift to nuclear spin polarization are spatially inhomogeneous. The resulting steady-state nuclear spin polarization is also spatially inhomogeneous, if only because the QD layer is sandwiched between the two doped semiconductor layers, where free charge carriers result in P N ≈ 0. On the other hand, when considering spins as a quantum resource, a significant role is played only by the nuclei within the QD electron wavefunction. Recent studies of nuclear spin relaxation in the same sample have shown that spin diffusion is the dominant mechanism of nuclear spin decay in a QD [31].
The nuclear spins at the center of the QD are quickly polarized by the optically-pumped electron spin and then transfer their polarization to more distant nuclei via nuclear spin flip-flops. For a sufficiently long pumping the nuclei in the AlGaAs barriers around the dot become gradually polarized. This manifests in a slow-down of the subsequent relaxation without the pump. The relaxation times of the nuclear spins are on the order of hundreds of seconds, much longer than the nuclear spin buildup times, which are less than one second. Such a large ratio of the timescales suggests that competition between spin pumping and polarization leakage would not be a limiting factor for achieving |P N | up to ≈ 0.99. Moreover, relaxation much slower than pumping means that spin diffusion creates a smooth spatial profile of the nuclear spin polarization -the extent of the polarized volume is larger than the volume of the electron wavefunction. Therefore, we expect that the electron probes a volume with a nearly uniform nuclear spin polarization degree P N . In other words, the existence of spin diffusion combined with long nuclear spin pumping means that there is no realistic mechanism that would result in abrupt spatial variations of P N . Under these conditions the selective-NMR thermometry measurements return the average of nuclear polarization, weighted by the electron envelope wavefunction density. Observation of a near-unity average polarization itself implies that polarization is very homogeneous for all nuclei within the electron wavefunction volume. For example, if we take the typical leakage of the electron wavefunction into the AlGaAs barriers at ≈ 0.1 (estimated previously in Ref. [8]) and assume full polarization within the GaAs layer (|P N | = 1), then observation of a weighted average of |P N | = 0.99 implies that polarization within the AlGaAs barriers cannot be much smaller than |P N | = 0.9. To summarize, although our present technique is not capable of revealing the spatial profile of the nuclear spin polarization, the most plausible hypothesis is that nuclear spin polarization achieved under steady-state optical pumping is nearly uniform within the volume of the QD electron wavefunction. In practice, this implies the ability to polarize nearly all nuclei whose coupling to the electron is strong enough to have any relevance to electron-nuclear coherent spin dynamics [25]. This also justifies our model, which assumes P N to be constant within the quantum dot volume and its surrounding.
Although 75 As, 69 Ga and 71 Ga are the three abundant isotopes, the inevitable penetration of the electron into the AlGaAs barriers implies some hyperfine interaction with the spin-5/2 27 Al nuclei. And yet it turns out that 27 Al hyperfine shift is to small to be studied quantitatively. As it has been shown previously [8], this is a combined effect of several factors. The small fraction of the wavefunction overlapping with AlGaAs (≈ 0.1), the small fraction of Al atoms (0.33 in our sample) and the small hyperfine constant (≈ 0.3 of that of As and Ga) mean that the 27 Al relative contribution to the total electron hyperfine shift is within ≈ 1%. In addition to that, the lack of Al at the center of the QD, where the overlap with the electron is the strongest, suggests inhibition of the pumping-through-diffusion mechanism discussed above. The 27 Al spins can only be polarized through direct (and weak) contact with the spin-polarized electron, meaning that aluminium polarization can be reduced. In the context of the present work where we focus on polarization of As and Ga, these observations mean that any systematic errors arising from 27 Al are small (within ≈ 1%). Investigation of 27 Al spin polarization would be an interesting subject for future work -this would require more sensitive experimental techniques, such as trigger detection via abundant isotopes [26].
Summarising this analysis, we see that there is a handful of potential error sources, both random and systematic, but all on the order of 1%. Taking a conservative approach we conclude with confidence that nuclear spin polarization degrees well above 0.95 are achieved. In reality, the polarization is likely to be higher, with rigorous confidence-interval analysis returning polarization degrees as high as |P N | ≈ 0.99 (for 69 Ga at high magnetic field B z = 10 T in all three selected individual quantum dots). Achieving even higher polarizations would depend critically on development of more sensitive thermometry techniques. One possibility is to use the dephasing dynamics of the electron spin qubit, since this would gain sensitivity at high polarizations as ∝ 1/ 1 − P 2 N