The impact of stochastic incorporation on atomic-precision Si:P arrays

Scanning tunneling microscope lithography can be used to create nanoelectronic devices in which dopant atoms are precisely positioned in a Si lattice within $\sim$1 nm of a target position. This exquisite precision is promising for realizing various quantum technologies. However, a potentially impactful form of disorder is due to incorporation kinetics, in which the number of P atoms that incorporate into a single lithographic window is manifestly uncertain. We present experimental results indicating that the likelihood of incorporating into an ideally written three-dimer single-donor window is $63 \pm 10\%$ for room-temperature dosing, and corroborate these results with a model for the incorporation kinetics. Nevertheless, further analysis of this model suggests conditions that might raise the incorporation rate to near-deterministic levels. We simulate bias spectroscopy on a chain of comparable dimensions to the array in our yield study, indicating that such an experiment may help confirm the inferred incorporation rate.


INTRODUCTION
Atomic precision (AP) placement of individual dopant atoms in Si nanoelectronic devices is a promising avenue for realizing a variety of technologies ranging from analog quantum simulators [1][2][3][4][5][6], to qubits [7][8][9][10][11][12][13], to digital electronics [14,15]. This paper considers a particular limitation of AP donor placement with scanning tunneling microscope (STM) lithography [16][17][18] (see Fig. 1). STM lithography allows for the fabrication of devices in which single donor atoms are positioned to within ∼1 nm of a target lattice site [19]. This weak placement disorder has previously been predicted not to be of concern when it comes to the prospects for realizing analog quantum simulators [5] or nuclear spin qubits [20] with this approach. However, these models do not account for the failure of a precisely placed donor to incorporate. While STM lithography achieves AP placement of individual donor atoms, we will show that a single precisely placed donor will incorporate with probability less than one for room-temperature dosing. While this is an indirect observation and quantitatively similar to previously reported results [21], our inference is corroborated by a kinetic model. That same model suggests conditions that might lead to deterministic incorporation.
Analog quantum simulation of the Fermi-Hubbard model using donor arrays provides a useful example that both introduces the physics of interest and the potential impacts of stochastic incorporation. The physics that is most relevant to analog quantum simulation is the hydrogenic nature of shallow donors [22]. In its neutral ground state, a single bulk-like P donor will have one weakly bound electron 45.59 meV from the bulk conduction band edge [23][24][25]. The affiliated orbital will have an anisotropic effective Bohr radius that exceeds that of a hydrogen atom in vacuum [26]. This renormalization is due to the bulk dielectric constant of the host crystal ( Si = 11.7) and the anisotropic conduction band effective mass (m ⊥ = 0.19m e and m = 0.92m e ). This artificial hydrogen atom can also be positively charged when it is stripped of its single weakly bound electron, or negatively charge when it is doubly occupied.
The low-energy effective theory that governs the behavior of the weakly bound electrons on a chain or array of such artificial hydrogen atoms is an extended Fermi-Hubbard model [4,5]. The tunnel coupling between sites is determined by the physical distance between donor atoms and it is feasible to realize spacings that are within an order of magnitude of the effective Bohr radius. This can then realize instances of the extended Fermi-Hubbard model that are "hard" in the sense that there is competition between itinerance and localization of the weakly bound electrons.
Atomically precise arrays of more than ∼50 donors are needed to realize a system that is large enough to represent an instance of the extended Fermi-Hubbard model that exceeds the scale for which exact simulation on a classical computer is feasible [27]. Prior theoretical analyses [4,5] have shown that the placement disorder realized in STM lithography is sufficiently weak to realize large arrays that exhibit the desired physics, assuming that each lithographic window yields exactly one P donor arXiv:2105.12074v1 [cond-mat.mes-hall] 25 May 2021 (1) Here P L (w) is the probability of patterning a lithographic window consisting of w dimers within a single row and P I (n|w) is the probability of incorporating n donors into such a window. Implicit in this formula is the independence of success probabilities from site to site and the neglect of failure modes outside of lithography and incorporation chemistry [28]. Both assumptions make our formula an optimistic estimate. P L (w) is ultimately a function of the STM hardware and control software whereas P I (n|w) is a function of the surface chemistry of window w. To reliably fabricate arrays with N >> 1 sites, i.e., P (N ) ≈ 1, it is evident that a window must be identified such that P I (n = 1|w) N ≈ 1. The data in Fig. 2 suggest that even for the best choice of w, a window consisting of 3 dimers, P I (n = 1|w = 3) = 63 ± 10%. Even with perfect lithography and processing the probability of successfully fabricating a ∼50-donor array for analog quantum simulation is ∼10 −10 .
In what follows, we substantiate this dismal observation with a model of the incorporation chemistry. However, we also suggest a path forward. While it may at first seem as if only P L (w) can be driven arbitrarily close to 1 through improvements to the STM hardware and control software, our incorporation model indicates that dosing at lower pressures for longer times and/or slightly higher temperatures can lead to P I (n|w) → 1. Intuitively, the w = 3 windows need sufficient time and space for a single PH 3 molecule to land in the window without interference and then follow the correct pathway to a configuration in which incorporation occurs. Our model suggests that the stochastic incorporation that we have observed is potentially a feature of room-temperature dosing. Because our incorporation statistics are based on inference, we briefly consider more direct measurable signatures of missing donors. We propose that bias spectroscopy measurements on a chain of dimensions comparable to the array fabricated in this work can yield insights into whether stochastic incorporation is occurring in a given process.

STM experiments
To generate a statistical sample sufficient to estimate P I (n|w) for different values of n and w, 50 windows were patterned using standard STM lithography techniques (see Fig. 1). Our aim was to pattern an array consisting entirely of w = 3 dimer windows into which the incorporation of n = 1 P atom should be maximized. While feedback lithography would result in near-perfect patterning of any target window size [29][30][31], we did not make use of it. Thus P L (w = 3) is less than one but this is irrelevant for our analysis as subsequent scans of the depassivated windows were taken to establish w for each window. After dosing and incorporation, STM scans of the areas containing the depassivated windows were taken again to determine n. Estimates for P I (n|w) are the relative frequencies of each observed n for the measured w.
We note that STM scans of the pre-dose and postincorporation surfaces were adjusted to bring the two scans into registration. The post-incorporation scan is aligned to the pre-dose scan using intrinsic defects in the reconstructed Si(100)-2×1 surface as markers (see Appendix B). After alignment, each window in both the pre-dose and post-incorporation STM scans are extracted and gridded with respect to each Si lattice site (lattice constant 5.4Å) for analysis. We expect that this alignment procedure results in a near-perfect overlap between the scans, where, in principle, only features present within a given window would be considered for evidence of incorporation. We allow for a single exception to this rule, where a Si:P heterodimer, as defined below, is allowed to be a single 2 × 1 dimer away from the window (Appendix B), as Si:P heterodimers represent direct evidence of P incorporation. It should be noted that windows that demonstrated no change between STM scans of the post-dosing and post-incorporation surfaces are removed from the sample as incorporation cannot be established for such windows (total of 3).
The value of w was determined for a given lithographic window by counting the number of contiguous depassivated Si atoms within a given dimer row. Depassivated atoms found in adjacent dimer rows or located non-contiguously to the primary pattern were excluded. If no more than one depassivated Si atom satisfied the criteria set above, then the pattern was considered disordered and removed from our sample. The variability in the number of depassivated Si atoms for the 47 windows is depicted in Fig. 2a. Here and elsewhere we specify w, the number of depassivated dimers rather than atoms, resulting in half-integer values for dimers that are partially depassivated.
For our target w = 3 pattern, only 20 windows were determined to meet the above criteria, i.e., P L (w = 3) = 43%. However, based on chemical considerations [32,33] the w = 3.5 pattern should exhibit similar incorporation chemistry, and including these windows in our sample increases the number of windows that realize the target pattern to 27, i.e., P L (w = 3) = 57%. Thus the estimates of P I (n|w) for w = 2, 3 in Fig. 2 include the corresponding half-dimer results (i.e., w → w , see Appendix A for results separating the half-integer values of w).
To determine the number of donors in a given window from the post-incorporation scans, it was necessary to categorize features according to their height relative to the surrounding H-resist. The distribution of feature heights in the w = 2, 3 windows are shown in Fig. 2. For a single donor to have incorporated successfully, two features should be present [33,34]-an incorporated P atom in the form of a Si:P heterodimer and an ejected Si adatom. Based on methodology used previously in the literature [21], we assign the first feature height grouping to Si:P heterodimers, with a relative height of 0.3-0.6 A [34][35][36] and the second to Si adatoms, with a relative height of 1.2-1.8Å [21].
However, we also observe two feature categories that aren't associated with successful incorporation. The first is a set of features with a relative height of 1Å that was previously assigned to depassivated Si atoms [37], consistent with height measurements of the pre-dose surface. The second set of features have relative heights greater than 1.8Å, which we ascribe to either background contamination in the UHV chamber or ejected Si adatoms adsorbed on H [38,39]. We presume that these possible Si adatoms are purely adventitious as they never correspond to windows with zero features, and therefore neglect them in further analysis.
Having assigned values of w to each window and categorized the features in the post-incorporation scans, we evaluated P I (n|w) from the relative frequencies of particular values of n. This required us to assign values of n to a particular window based upon the presence of particular features, specifically the Si:P heterodimer and Si adatom features. We note that both features were counted individually in Fig. 2b. Due to the different diffusion barriers for the Si:P heterodimer [35,40] and the Si adatom [41,42], we do not expect that both features are necessarily present within a particular lithographic window. Thus, if either feature was found within a particular window it was counted toward an incorporation event. This brings us to one of our key results, that P I (n = 1|w = 3) = 63 ± 10%, consistent with previously reported statistics for similarly sized windows [21].
We note that we expect P I (n > 1|w = 3) = 0% due to a paucity of dangling bonds on which 2 PH 3 molecules might shed 4 H atoms to arrive at a bridging P-H configuration in which the P can incorporate. 2-donor incorporation requires at least w = 4 depassivated dimers. However, we observe a small but non-zero probability for P I (n = 2|w = 3). We hypothesize that this is a consequence of features diffusing into the window, rather than 2 donors incorporating within the same w = 3 window. We suggest that a similar mechanism might also be biasing our estimate of P I (n = 1|w = 2) upwards. We note that this bias is not accounted for in our quoted uncertainties.
Returning to Eq. 1, our results indicate that P (N ) = (40% ± 22%) N for our particular process [43]. However, if P L (w = 3) → 100% using, e.g., feedback lithography, we see that P (N ) → (63% ± 10%) N . In both cases, the prospects for incorporating a large array for analog quantum simulation are dismal. To better understand these statistics we corroborate our results with an atomistic model that accounts for stochastic single-donor incorporation. P likely incorporates via two different pathways: a one-step process, or a two-step process with a thermodynamically unfavorable middle step. We indicate the lower end of the dimer row by underlining the Si label in red. Labels here refer to the notation used by Warschkow et al. [32]. (c) Incorporated P calculated using an instance of our KMC model in which the barrier of the B3 to C1 reaction is lowered by 0.1 eV. Incorporation is generally increased, matching 3-dimer incorporation to the experimentally measured probability

Incorporation model
We developed a Kinetic Monte Carlo (KMC) model to estimate P I (n|w). Our model includes reaction barriers from Warschkow et al. [32] and Wilson et al. [33] with additional features described in the supplementary information. Focusing particularly on the windows for which we expect single-donor incorporation, we compute P I (n|w) for w = 2, 3, and 4. The predictions of our model are compared to our experimental results in Fig. 3a. We match the ratio of half dimers included in each width bin to the measured ratio from this work.
For n = 1 and w = 2, 3 our KMC model matches the experimental incorporation rate within error bars. In the three dimer wide window, single-donor case, in particular, we predict P I (n = 1|w = 3) = 59% ± 2% compared to the measured value of 63% ± 10%. For each reaction within our model chemistry, the rate is governed by an attempt frequency and a reaction barrier according to the Arrhenius equation. As the rates are exponentially sensitive to the reaction barriers and errors ∼0.1 eV are typical of the semilocal DFT calculations that predicted them, we examine whether we can better reproduce experiment by comparably sized adjustments to our barrier heights. Rather than assessing the sensitivity of incorporation rates to every barrier independently, we first identify particularly sensitive steps in the reaction pathway.
The two reaction pathways that dominate the incorporation of a P atom are depicted in Fig. 3b. In the first pathway, outlined in the supporting information, a PH 2 molecule in a lower dimer end position (B2, using the naming conventions of Warschkow et al. [32]) moves into a bridging position between the two dimer atoms while simultaneously losing a hydrogen to the nearby dimer (C1). Alternatively, the PH 2 molecule can first move to a nearby dimer, but on the raised end of the dimer (B3), before finally moving to a bridging PH position (C1). This latter pathway involves a two-step process to incorporation, where the middle step is not thermodynamically favorable. Since the barrier for moving from B2 to C1 and the barrier for moving from B2 to B3 are within 0.05 eV of each other, a system that starts in a B2 configuration has a roughly equal chance of moving directly to C1 and incorporating or moving to B3. From B3, the system again has a roughly equal chance of incorporating to C1 or transitioning back to B2.
This balance of options in the two-step pathway means that even slight changes that increase the likelihood of moving from B3 to C1 instead of from B3 back to B2 can have significant impacts on the final incorporation. In Fig. 3c, we examine the impact of lowering the B3 to C1 barrier by 0.1 eV while keeping all other reaction barriers constant. We find that lowering the barrier for the B3 to C1 reaction leads to P I (n = 1|w = 3) = 67% ± 2%. Further decreases to this reaction barrier eventually saturate before deterministic incorporation is achieved: lowering the B3 to C1 barrier by 0.2 eV, for instance, only results in P I (n = 1|w = 3) = 70% ± 2%.
Accordingly, it is worth examining our model to see whether we can increase P I (n = 1|w = 3) to 100%. We also consider whether we can corroborate the number of donors in an incorporated device through measurements other than the number of ejected Si adatoms and Si:P heterodimers.

Potential for deterministic incorporation
Dose time and pressure can also strongly affect incorporation and we present predictions for the incorporation rate as a function of both in Fig. 4. Notably the incorporation varies even at constant exposure, indicated as diagonal lines. A system at room temperature exposed to a 0.18 Langmuir dose of PH 3 at a pressure of 3 × 10 −8 Torr for 6 seconds has P I (n = 1|w = 3) = 52%±2%, while the same exposure at a pressure of 3 × 10 −11 Torr for 6,000 seconds has P I (n = 1|w = 3) = 64% ± 2%. In general, moving toward conditions of lower pressure and longer dose time leads to increased incorporation.
In fact, for extreme conditions of low pressure and long time, we predict near-deterministic incorporation. At 3 × 10 −12 Torr for 6 × 10 5 s (∼6.9 days), P I (n = 1|w = 3) = 99% ± 1%. Currently these conditions may be prohibitively difficult to achieve, thanks to the need for extreme ultrahigh vacuum and exceptionally pure PH 3 dosing over several days, but it is instructive to understand the factors that lead our model to predict this.
Given enough space and time for a single PH 3 molecule in a window to dissociate before another PH 3 molecule adsorbs in the same space, incorporation will occur. The thermodynamic pathway for PH 3 dissociation is entirely downhill. There are no dead ends in the path that would prevent it from eventually incorporating, only wrong turns that can be reversed given sufficient attempts. As pressure decreases, the likelihood of a single PH 3 molecule adsorbing into a w = 3 window fully dissociating before additional PH 3 adsorption increases. At 3 × 10 −12 Torr, the average time for an incorporation event is 3.4 × 10 5 s. If only one PH 3 molecule adsorbs onto the w = 3 window during the 6 × 10 5 s dose, that molecule will be able to to shed its hydrogen and incorporate without any competition from neighboring molecules. In contrast, for a dose pressure of 3 × 10 −8 Torr, increasing the dose time can decrease the likelihood of incorporation due to additional phosphine molecules adsorbing into the window before the molecule has fully dissociated. Once PH 3 evolves into a bridging PH configuration, however, the barriers for returning to a PH 2 state are so high that any recombination is unlikely, even in the extreme time scales considered here. Moving from a C1 configuration to a B3 configuration requires overcoming a barrier of 1.86 eV, giving an average reaction time of 1.8 × 10 19 s at room temperature and 6.4 × 10 3 s at the anneal temperature of 320°C.
Incorporation rates can therefore be increased by ensuring that the first PH 3 molecule adsorbed on the surface dissociates before another PH 3 molecule is adsorbed, blocking its pathway to incorporation. By increasing the sample temperature during dosing, the typical time for the dissociation reaction can be decreased significantly, as demonstrated in Fig. 4b-c. At temperatures as low as 75°C, near-deterministic incorporation is possible within a more typical range of pressures and temperatures, with conditions of 3 × 10 −10 Torr for 6000 s (∼1.6 hours), P I (n = 1|w = 3) = 98.8% ± 1%. This trend continues to expand at higher temperatures with all exposures > 1 L producing deterministic incorporation at 150°C. Given the likelihood of hydrogen migration ruining the fidelity of the lithography window at higher temperatures, we expect that finding a feasible combination of pressure, time, and moderate temperature will produce the highest chance of success.
Our simulations predict that stochastic incorporation of a P donor in a w = 3 window is not due to any inherent limitations of the PH 3 chemistry. It is instead due to the limitations of current room-temperature dosing techniques. This suggests that heating the sample during dosing and/or the development of ultra-pure chemical precursors and the ability to dose at pressures at or below 3 × 10 −12 Torr may make deterministic donor incorporation with atomic precision feasible using a PH 3 chemistry. Our results also suggest that with current technology and dose temperatures, incorporation can still be increased by moving to lower pressures and longer dose times. We additionally consider the benefits of lowtemperature (∼100 K) dosing for deterministic incorporation in the supplementary information.

Corroborating incorporation statistics through transport measurements
In assessing whether stochastic incorporation occurs in a given fabrication process, it is critical to consider other data that could support or refute incorporation statistics derived from observed Si adatoms and Si:P heterodimers. We have already noted that our inferred value of P I (n > 1|w = 3) > 0% is suspect and it is desirable to have a direct measurement of the electrically active substituted donors, more so because these are the technologically relevant elements. Bias spectroscopy is thus a natural choice and a powerful characterization technique that can be used to directly measure the spectrum of current-carrying many-body states in a donor array. We consider the properties of a small donor array that might be probed with this technique to corroborate our incorporation statistics.
A bias spectroscopy measurement would require fabricating such a donor array between two δ-doped nanowire source and drain leads with either another in-plane lead or top gate [44] to adjust the on-site potential across the array. To reduce the likelihood of confounding the measurement with sites consisting of n > 1 donors, we propose making an array comprised exclusively of w = 3 windows. Then it is likely the case that each site will consist of n = 0 or n = 1 donor. A one-dimensional chain with a short spacing between sites will be ideal. This spacing should be short enough that we can still expect conduction through the chain in case of missing sites. The number of sites in the chain should be both short enough to maintain good stability in the patterning and long enough to have a reasonable expectation of missing donors. A 5-site chain strikes a balance between being long enough to have 2 missing sites (in expectation) and short enough to expect good tip stability in patterning. Indeed, a single row of the array in Fig. 1 would suffice. For such a chain, the spacing between sites would correspond to a nearest-neighbor tunnel coupling of ≈ 10 meV, diminishing to ≈ 1 meV in the event of one missing site or ≈ 0.1 meV for two consecutive missing sites [45]. Going to half this spacing is still viable and further increases the minimum tunnel coupling realized in the event of two missing sites by over an order of magnitude. We consider simulations of such a chain in Fig. 5.
Our simulated bias spectra are intended to be qualitatively illustrative of the impacts of missing sites. We present an exhaustive enumeration of all possible missing donors, as well as other spacings and orientations that may be of experimental interest in the supplementary information. Here we simply compare the spectrum of the ideal chain to the most likely configurations consistent with P I (n = 1|w = 3) = 63 ± 10% (i.e., two missing donors). We expect the spectrum for the chain in which every other donor is missing (1 of the 10 possible configurations for a five donor chain) to be qualitatively similar to that of the unbroken chain in this regime, with four prominent diamonds. But the majority of configurations' bias spectra (9 of 10) will more closely resemble the case in which there are two consecutive donors and two missing sites, with three prominent diamonds, as demonstrated in the supporting information. We thus expect that even in conditions favorable to tunnel coupling through missing donor sites, missing donors consistent with an incorporation rate of ∼60% can be clearly detected with bias spectroscopy in 90% of possible configurations for a five donor chain. A larger spacing will yield even greater sensitivity to the precise locations of the missing donors.
A quantitatively accurate prediction of the current In all cases, the faint transitions at higher occupancies (i.e., on the left of the plots) are diminished due to the increased impact of the intersite Coulomb interaction. We note that our calculations do not include current through scattering states to clearly demarcate the transport through the eigenstates of the chain and thus our Coulomb diamonds do not have signatures of transport "above" or "below" them, as would be expected in a real experiment.
through a few-donor chain would require some improvements to our model that are discussed in the supplementary information. However, we expect the basic principle illustrated in Fig. 5 to hold-for a sufficiently short chain there will still be measurable Coulomb diamonds even in the presence of missing donors and deviations from the transport signatures expected in the ideal chain might herald a particular pattern of missing donors. We note that the transport signatures of missing donors are stronger for calculations using the traditional Fermi-Hubbard model (i.e., without the intersite Coulomb in-teraction). In that limit, the incipient Hubbard bands above and below half filling give indications that clearly reflect the longest unbroken chain. Because the intersite Coulomb interaction in the extended Fermi-Hubbard model suppresses the formation of these incipient bands in small chains and breaks the symmetry between the bands above and below half-filling, we expect these signatures to be a bit messier. Accordingly, we also expect that transport measurements of a chain of this dimension would provide experimental confirmation that the intersite term is a critical feature to include in effective models of Si:P arrays.

Conclusion
We demonstrated that even for perfectly patterned lithographic windows, single-donor incorporation is stochastic with a 63 ± 10% likelihood of success at roomtemperature dosing conditions, which we corroborate with a kinetic model. We thus conclude that building classically intractable single donor arrays is overwhelmingly unlikely under these conditions. Our model nevertheless suggests that near-deterministic incorporation may be possible above room temperature or at lowpressure, long-time dosing conditions. We also considered bias spectroscopy as a method to directly measure the number of incorporated donors in an array, demonstrating that irregularities in the pattern of Coulomb diamonds can give clear indications of missing donors.
It is also worth remarking on the prospects of using AP technologies to realize analog quantum simulators with single acceptors [46,47], in addition to their prospective use in qubit technologies [48,49]. The Fermi-Hubbard model parameters for a pair of substitutional B dopants have been measured in a sample prepared without AP placement [50]. With the recent demonstration of B 2 H 6 as a precursor for delta doping [51], it is clear that AP placement of individual B atoms might also be possible. An incorporation model for B 2 H 6 , similar to the one in this paper, suggests that the tendency for dimerized B to incorporate as electrically inactive might limit the prospects for achieving deterministic single B atom placement for analog quantum simulation [52]. This is qualitatively different from the case presented here for PH 3 , in which we predict that changing the pressure and temperature of the PH 3 dosing conditions can ultimately lead to deterministic AP placement. Nevertheless, the model incorporation analysis also suggested that non-dimer acceptor precursors and halogen resists might provide a path forward [52].
Finally, it is worth emphasizing that stochastic incorporation is not necessarily a no-go for the application of STM lithography to realizing analog quantum simulators, qubits, or digital electronics. A simple workaround is to target the fabrication of multi-donor clusters [7,53] to realize quantum dots with precise dimensions. While our results suggest that there will be uncertainty in the number of donors in such a cluster, it may be that there is a "sweet spot" in target cluster dimensions for which the resulting variance in the resulting application-relevant properties (e.g., tunnel couplings, charging energies, etc.) is sufficiently low. In this context, the single-donor limit is somewhat pathological as there is a significant probability of not fabricating a cluster at all. The demonstration of a process that reliably realizes deterministic incorporation of a single donor into a three-dimer window, such as suggested by our modeling, would obviate some of the need for these considerations.
Samples were terminated with a monohydride resist for STM lithography by heating the sample to 350 • C while exposing to atomic H supplied by cracking H 2 (2 × 10 −6 Torr, 10 mins) with a W filament positioned 1 cm away from the sample. All STM operations utilized PtIr STM tips from NaugaNeedles (NN-USPtIr-W250, radius of curvature 25-50 nm), degassed at ∼ 175 • C for 30 minutes and further heated by electron bombardment using a commercial Omicron tip preparation tool. All AP patterns were generated using +3.5V bias, 6.0 nA tunneling current, 6.0 mC/cm dose, and 10 nm/sec tip speed.
PH 3 is introduced into the UHV chamber via a precision leak valve, with PH 3 coverage calculated using the partial pressure of the PH 3 introduced during dosing measured by an RGA (specifically the 34 amu fragment). In a typical process, a PH 3 partial pressure of 3.0×10 −10 Torr is introduced for 10 minutes at room temperature, resulting in a coverage of 0.15 Langmuir, which has been previously demonstrated to provide sufficient coverage for atomically precise donor structures [54,55]. While the partial pressure reported through the PH 3 RGA value may represent a lower bound on the actual partial pressure introduced, the reported dosing conditions have also been demonstrated to provide the expected 0.25 ML coverage on a clean Si(100) surface for P incorporated at 310 • C [56], further demonstrating that the dosing conditions utilized will result in sufficient coverage in the depassivated windows. A subsequent 310 • C anneal is used to incorporate P into Si lattice sites, demonstrated to be sufficient for full donor incorporation [57].
All STM scans shown were taken under -2.75V, 300 pA imaging conditions.

Theory and modeling
We use a Kinetic Monte Carlo model [58,59] as implemented in the KMClib package [60] to determine P I (n|w) as a function of n and w. Each KMC calculation is repeated 200 times with different random seeds, and the sample mean of the results is reported. We calculate error bars by assuming a binomial distribution of measured counts and using the standard error based on sample size. In all calculations unless otherwise stated, we use the experimental dose and anneal temperatures, pressures, and times.
Our transport calculations are based on extended Fermi-Hubbard model parameters extracted from the multi-valley effective mass theory first presented in Ref. 20. The Meir-Wingreen formula [61] is used to compute the current through interacting donor chains and the differential conductance is evaluated from the current using a finite difference formula.

AUTHOR CONTRIBUTIONS
JAI, QC, JCK contributed equally to this work. EB, SM, and MSC coordinated the experimental effort. ADB coordinated the theoretical effort and the writing of the manuscript. JCK planned and conducted the experiments generating the STM data. DRW made critical advances to fabrication. JAI, JCK, QC, and AMM analyzed the experimental data. QC, PAS, and RPM developed the model for donor incorporation. MIB carried out the transport calculations.
Supplementary information: The impact of stochastic incorporation on atomic-precision Si:P arrays

CONTENTS
The supplementary information elaborates on details of some of the central results in the main body of the paper.
• Appendix A describes experimental details regarding lithographic outcomes deviating from perfect two or three dimer windows and post-lithography/post-dose images.
• Appendix B describes the image registration procedure used to compare the pre-dose and post-incorporation STM scans.
• Appendix C describes details of our incorporation model, including an analysis of a low-temperature process described in a recent patent (Ref. 62).
• Appendix D describes the extended Fermi-Hubbard model and our simulated bias spectra on short donor chains.

APPENDIX A: EXPERIMENTAL DETAILS
While using STM lithography to create three silicon dimer wide windows (w = 3), we often create alternative outcomes such as windows of two (w = 2) or four dimer (w = 4) width, the inclusion of half dimers, and disordered arrays. We display a selection of common lithographic outcomes in Fig. A1. After dosing the lithographic patterns with PH 3 gas, we examine the adsorption configurations of PH 3 fragments in the sites. We observe three general outcomes, as illustrated in Figure A2: (a) saturation of the lithographic template with PH 3 gas fragments, (b) less than saturation coverage of the lithographic template with PH 3 gas fragments, and infrequently (c) complete repassivation of the litho site with hydrogen. These observations of sites saturated with PH 3 fragments and sites with less than saturation coverage with PH 3 fragments occur in spite of dosing conditions which lead to 0.25 ML P dopant density for large area patterns.
In the main text, we have compressed results for two (w = 2) and two and a half dimers (w = 2.5) into two dimers, and three (w = 3) and three and a half dimers (w = 3.5) into three dimers. In Fig. A3, we separate these contributions out. We find that adding a half dimer to the lithographic window tends to decrease the overall rate of incorporation, although the total sample size, particularly for two dimer wide windows, remains too low to draw significant conclusions.  (20) 3.5 dimers (7) FIG. A3. The statistics of how many donors are incorporated expanded to include half dimer windows.

APPENDIX B: MINESWEEPER ANALYSIS
In order to determine P I (n|w), the pre-dosing and post-incorporation datasets had to be "atomically registered" to each other such that each lattice site corresponds to the same pixel(s) in both datasets. This was accomplished by analyzing and transforming data sets with OpenCV followed by labeling with a custom labeling graphical user interface as detailed in the following.
The first step is to ensure that the pre-dosing to post-incorporation data sets had the same pixel-to-atomic site scaling. This was accomplished by finding the lattice periodicity via Fourier transform and scaling the post-incorporation data to match the scaling of the pre-dosing data sets. We then used the lattice defects as alignment markers between the two data sets by thresholding both data sets and nulling out all data above the approximate height of the defects. A correlation function was applied to the datasets to find the best overlap. Occasionally this produced an insufficient matching and an affine transform was produced by matching three specific defects across the images. Both methods produced a transform function that effectively "atomically registered" the pre-dosing and post-incorporation datasets to one another.
After the images had been "atomically registered", the datasets were segmented into small regions around the depassivation window and corresponding post-incorporation area. These sub-images were loaded into a custom labeling graphical user interface where the authors shifted each sub-image to adjust the image so that each lattice site corresponded to an 8x8 set of pixels. Each lattice site is labeled as XXX, YYY through analysis of the contrast relatively to the neighboring lattice sites in the STM topographical image. In order to account of uncertainty in both the registration and feature height assignment procedures, both Si ad atom and Si:P heterodimer features within a lithographic window are labeled as well as Si:P heterodimer features found within a single 2 × 1 dimer away from a given window. Such a Si:P heterodimer is included in the analysis only if neither a Si ad atom or Si:P heterodimer is found within the lithographic window, and is considered conclusive evidence of P incorporation as Si:P heterodimers can only result from P incorporation. Several features which would correspond to Si ad atoms were rarely found to be located close to a patterned windows (within several 2 × 1 dimers) where no features were assigned, however as no Si:P heterodimers were found within those respective windows or surrounding them, we do not conclude that these ad atoms indicate the presence of P incorporation. After labeling, each labeled site is located within the overall postincorporation scan, with the corresponding feature height measured relative to the surrounding interdimer H-resist.
In both the pre-dosing and post-incorporation scans, care was taken to avoid incorrect height assignments resulting from "halo" that surrounds dangling bonds and other absorbates on the H-terminated Si(100)-2×1 surface while using empty-state imaging conditions [63].

A. Novel Density Functional Theory Configurations
Our Kinetic Monte Carlo model uses reactions and rates that are aggregated from results in the extensive literature on PH 3 dissociation on this particular surface [32,33,40,[64][65][66]. In addition to the phosphine configurations and reactions, we add an initial adsorption condition and a concerted reaction.
Our initial condition (which is implicit in Ref. 32, but we are explicitly stating here), relates to the up-down tilt of silicon dimers a bare Si(100)-2×1. We predict phosphine to have an adsorption energy of -0.91 eV at the lower side of the silicon dimer as opposed to an adsorption energy of -0.56 eV on the higher side. In our model, we therefore assume that a phosphine molecule will preferentially adsorb on the lower end of the silicon.
We also discover an additional concerted reaction (labeled as B2 to C1 in Fig. 3b within the main text). In this reaction, a PH 2 in the favorable lowered dimer position loses a hydrogen atom to a nearby dimer and moves to a bridging PH position at the same time. This reaction is exothermic with a thermodynamic gain of -1.03 eV and a barrier of 0.96 eV. This contrasts with the similar reaction pathway (B3 to C1) from Warschkow et al. as the PH 2 molecule starts in the thermodynamically favored lower end of the silicon dimer (B2) instead of the higher end (B3), which requires an endothermic reaction to move the PH2 from the favored lowered end (B2) to the higher end (B3) with a thermodynamic cost of 0.24 eV and a barrier of 0.92 eV. (While B3 to C1 is technically a two step reaction within the Warschkow et al. framework, since the second step is thermodynamically downhill and barrierless, we have consolidated it as one reaction within this work.) Electronic structure theory calculations for these additions are performed with the Gaussian-basis, local-orbital pseudopotential code SeqQuest. [67] We use a norm-conserving Perdew-Burke-Ernzerhof (PBE) functional. [68] A seven layer thick 4×4 supercell slab of the Si(100)-2×1 is constructed and the dangling bonds at the unreconstructed bottom surface were terminated with selenium atoms, as these were found to minimize strain on the slab. We use a 20Å vacuum region to ensure the surfaces at either end of the slab are sufficiently isolated from each other. We sample the Brillioun zone with a 2×2×1 grid. To achieve total energies relaxed within ∼0.01 eV of the structural minimum, atomic positions are relaxed until the forces are reduced to less than 0.0003 Ry/Bohr. The barriers described above are then used as inputs in an Arrhenius model [69] for the transition rate Γ = A exp ∆/k B T , where A is the attempt frequency, ∆ is the reaction barrier from DFT, k B is Boltzmann's constant, and T is the temperature. We consider A = 1 × 10 12 s −1 for all reactions, which is the order of magnitude calculated for these reactions by Warschkow et al. [32].
We calculate the effusive flow rate of molecules landing on any particular silicon dimer as Φ ef f usion = P A/ √ 2πmk B T , where P is the pressure of the incoming precursor gas, A is the area of impingement, taken here as a single silicon dimer, m is the mass of the precursor gas, k B is the Boltzmann constant, and T is the temperature.
We repeat each KMC calculation 200 times to obtain a meaningful statistical sampling of likely outcomes and report the average outcomes, along with standard deviations as applicable. For ease of reproducibility, we have placed our Kinetic Monte Carlo code on GitHub [70]. We calculate the standard error by assuming a binomial distribution of measured counts and using the standard error based on sample size.

B. Low-Temperature Dosing Results
We also use our kinetic Monte Carlo model to examine low-temperature ( < 100 K) dosing techniques to achieve saturation before annealing such as in a recent patent. [62] We simulate a process where the initial dosing step is done at a temperature of 100 K, and subsequent annealing at a temperature of 500°C for 15 s. As shown in Fig. C1 in saturated dosing conditions (i.e. > 1 L coverage), near-deterministic doping can be achieved. Dosing at low temperatures allows saturation coverage of the exposed window to occur. At 100 K, enough energy is provided for an initially adsorbed phosphine to shed one hydrogen, but not for further reactions. The subsequent anneal at 500°C provides enough energy for the phosphine fragments to desorb in a reasonable time frame. This provides enough room for the remaining phosphine fragments to fully dissociate into bridging PH fragments. For a 2 dimer wide window, this process is nearly deterministic with a predicted incorporation rate of P I (n = 1|w = 2) = 96% for a dose pressure of 3 × 10 −8 Torr and a dose time of 600 s, as shown in Fig. C1c. However, for 3 dimer window, this process leads to two incorporations with a non-negligible probability; P I (n = 2|w = 3) = 13% for a dose pressure of 3 × 10 −8 Torr and a dose time of 600 s.
While this process can lead to near deterministic incorporation, it is also highly likely to lead to hydrogen desorption in the surrounding resist. This is due to the barriers for H 2 desorption and PH 2 + H desorption being essentially the same. After a low-temperature dose, the entire window will be saturated with PH 2 + H fragments and not full PH 3 molecules, making the PH 2 desorption barrier the relevant rate limiting step for full dissociation. Desorption of PH 2 + H has a barrier of approximately 2.0 eV, while H 2 desorption has a barrier of 2.1 eV. [33] Any anneal sufficiently hot to desorb PH 2 in a reasonable time frame would therefore almost certainly also desorb hydrogen from the surrounding silicon and destroy the integrity of the surrounding resist. While the small difference in barriers between PH 2 + H and H 2 desorption can theoretically be exploited to find an anneal schedule that maximizes PH 2 + H dissociation while leaving H 2 , we find that this would lead to non-deterministic incorporation. Using a 100 K dose with a 10 minute 400°C anneal had the best effect for maximally encouraging PH 2 + H dissociation with minimal hydrogen resist degradation, but the maximum incorporation rate was only P I (n = 1|w = 3) = 79%.
This analysis predicts that saturation dosing at low temperature and then high temperature annealing techniques are a viable route to deterministic incorporation, albeit with the potentially significant drawback of losing hydrogen resist integrity and thus introducing more uncertainty into the precise location of incorporation. We are careful to note the obvious fact that experiment is the ultimate arbiter of whether the resist integrity is degraded in practice. Combined with the findings in the main text, this indicates that the experimentally stochastic incorporation rates may, in fact, be an artifact of room temperature dosing. Other regimes of both low and moderate temperature dosing with the correct combination of dose pressure and time are predicted to be more feasible routes toward atomic precision deterministic single donor doping. The four sums comprising this Hamiltonian describe the following physics: 1. The energy ( j,σ ) of a single-particle orbital localized on site j with spin σ. This includes the effects of externally applied electric fields. The particle number operator for site j and spin σ isn j,σ =ĉ † j,σĉ j,σ , whereĉ † j,σ (ĉ j,σ ) is the associated creation (annihilation) operator associated with that fermionic spin-orbital.
2. The spin-independent tunnel coupling (t j,k ) between sites j and k. The brackets in the summation indicate that it is restricted to nearest-neighbors, though we note that a more detailed examination of this assumption for arrays with short pitches will be the topic of future work.
3. The Coulomb repulsion (U ) associated with double occupation of site j by electrons with opposite spins.
4. The Coulomb repulsion (V j,k ) between electrons occupying sites j and k. The density operator associated with site j,ρ j =n j,↑ +n j,↓ . We note that the conventional version of this model is limited to nearest-neighbor intersite repulsion.
We rely on an implementation of multi-valley effective mass theory (MV-EMT) [20] first developed at Sandia National Laboratories for the parameters of this model. This particular approach to MV-EMT due to Shindo and Nara [71] goes beyond the more common Kohn-Luttinger theory [26] and includes a non-perturbative and consistent treatment of both the central-cell correction and valley-orbit coupling that captures the fine structure of the valleyorbit split 1s-like states for a single donor. Recent extensions of this implementation of MV-EMT that account for electron-electron interactions in the modeling of exchange-coupled donors indicate that this theory can be adjusted to capture the charging energy of a single doubly-occupied donor, as well [72]. The t j,k parameters are taken from an exhaustive set of tunnel couplings between individual pairs of donors in bulk silicon that were made available in the arXiv version of Ref. 20 [73]. The onsite repulsion is U = 43.8 meV to reproduce the charging energy of a single bulk-like D − state. V is simply the static screened Coulomb interaction between charge centers localized on donor sites with Si = 11.7.
We expect the accuracy of this particular parameterization will be lacking relative to an experiment for a number of reasons. The tunnel couplings were computed between pairs of bulk-like donors, the U term does not account for the fact that the doubly occupied orbitals will have a different character in a chain with a short pitch (i.e., we expect U to be lower), and overall we expect that the single-particle orbitals will look quite different from our base MV-EMT. All of these issues are compounded by the fact that we are studying non-uniform chains in which some sites have failed to incorporate. For the purposes of the qualitative analysis in this paper, we argue that these limitations are acceptable and that experimental design and interpretation will require a more careful parameterization. We note that ongoing work is aimed at developing such a parameterization using an mesh-based MV-EMT capability that makes use of a discontinuous Galerkin discretization of the Shindo-Nara equations and device-specific self-consistent electrostatics. While donor chains are particularly challenging to model with this level of detail, results from this solver (adapted to the Luttinger-Kohn Hamiltonian) have been presented for the analysis of relatively simple Ge hole quantum dots in Ref. [74].
To simulate bias spectroscopy, we use the Meir-Wingreen formula [61] to compute the current through a donor chain modeled as an instance of the extended Fermi-Hubbard Hamiltonian in Eq. D1. We use exact diagonalization to compute the necessary Green's functions, which is computationally tractable for this simplistic model in which each site is comprised of two spin orbitals. A two-point finite difference formula is used to compute the differential conductance from the current, computed at 501 points for each instance.
We present simulated bias spectra for the exhaustive set of one-dimensional chains with two pitches, corresponding to ≈ 2.5 nm and ≈ 5 nm spacings along [100] in Fig. D1 and D2 (respectively) and [110] in Fig. D3 and Fig. D4. The spacing and length total length of the chain are chosen to be sufficiently short that it is reasonable to expect that the tunneling current can be measured even in the absence of missing sites. We note that our calculations do not include current through scattering states to clearly demarcate the transport through the eigenstates of the chain and thus our Coulomb diamonds do not have signatures of transport "above" or "below" them, as would be expected in a real experiment.
In each of the following figures, the individual rows correspond to a particular number of single-donor sites missing. For a 5-donor chain, there are thus 20 distinct non-trivial instances (i.e., 5 instances with one donor missing, 10 with two, 10 with three, and 5 with four). These are indicated on top of each bias spectrum with "x" indicating a missing site and "o" indicating a donor. The scale of the color bars is allowed to vary to capture the relevant details in each plot and only the maximum value of the differential conductance is specified in each panel.