A Computational Workflow for Designing Silicon Donor Qubits

Developing devices that can reliably and accurately demonstrate the principles of superposition and entanglement is an on-going challenge for the quantum computing community. Modeling and simulation offer attractive means of testing early device designs and establishing expectations for operational performance. However, the complex integrated material systems required by quantum device designs are not captured by any single existing computational modeling method. We examine the development and analysis of a multi-staged computational workflow that can be used to design and characterize silicon donor qubit systems with modeling and simulation. Our approach integrates quantum computational chemistry calculations with electrostatic field solvers to perform detailed simulations of a phosphorus dopant in silicon. We show how atomistic details can be synthesized into an operational model for the logical gates that define quantum computation in this particular technology. The resulting computational workflow realizes a design tool for silicon donor qubits that can help verify and validate current and near-term experimental devices.


Introduction
Theory predicts quantum computing will speedup solutions to important problems in science, engineering, and security [1]. Realizing these speedups will require a machine that is capable of executing computational gates on arrays of physical qubits over timescales that are much shorter than the influence of environmental noise and decoherence [2]. The fault-tolerant operation of quantum computers will require manipulating millions of qubits over billions of instruction cycles [3], while the architectures for these large-scale systems will require sophisticated execution and run-time control systems [4]. However, the basic operating principles of quantum computation can be demonstrated using a more modest number of qubits. This includes demonstrating the non-local correlations inherent to entangled particles and showing the computational steps needed by few-qubit algorithms [5,6].
Recent breakthroughs in the realization of donor qubits in silicon with high gate fidelities and long qubit coherence times have increased interest in these quantum computing architectures [7,8,9]. However, the current proof-of-concept experimental demonstrations must be viewed as point solutions for the larger problem of designing robust quantum computing devices. Key qualities of robust device operation include stability, reproducibility, and reliability, which are currently lacking in existing experimental systems. This is partly due to the uncertainty in the physical processes used for material fabrication, which hinders the experimental effort to realize specific physical designs. As insights into fabrication are gathered from early experiments, we may also expect modeling and simulation to provide guidance on how to design stable and reproducible qubit devices. Modeling and simulation have proven especially useful for the development of conventional CMOS (complementary metaloxide semiconductor) processors, where TCAD (technology computer aided design) tools are standard for verifying designs generated from well characterized process models. Near-term progress in the development of robust quantum devices is poised to benefit from similar TCAD tools that capture the details of process variation, physical model uncertainty, and operational details in the computing context.
Currently, there is no single, standalone tool capable of capturing the quantummechanical, electrodynamical, and quantum computational features needed for the verification of silicon donor qubit device designs. This is because the multi-physics and multi-scale models required for describing such devices span across the simulation domains of existing state-of-the-art tool sets. Consequently, there is an outstanding need to integrate the physical and logical models for qubit devices into a single computational workflow that can provide assessments for expected performance against known limitations. In this contribution, we address the development of a computational workflow that can be used for modeling and simulation to verify silicon quantum computing devices. A qubit design tool is necessary to support efforts in fabricating quantum physical systems that express individual qubits, controlling their programmed interactions, and carrying out these interactions in a computationally meaningful way. Our approach is based on a workflow framework that addresses the performance of small-scale doped silicon nanostructures, where models of physical layout and material design must be combined in order to evaluate the operating principles of an encoded qubit. This computational workflow exposes the ability to investigate key qualities of a device, e.g., gate fidelity, by tuning the design parameters. The computational models used here provide a natural bridge between efforts to fabricate and characterize silicon donor systems and efforts to program future multi-qubit systems [10,11].
Our computational workflow complements state of the art experimental capabilities for investigating the development of silicon donor systems under different material and environmental conditions [12]. In addition, our approach to modeling and simulation mimics some of the long-standing methods used in CMOS device fabrication. This includes conventional TCAD tools, like Sentaurus Device from Synopsys, that provide robust methods for simulating the electrical, thermal, and optical properties of CMOS-based devices. These conventional TCAD tools are typically part of a larger simulation suite that intends to capture all aspects of conventional semiconductor chip development. However, the physical assumptions underlying these models are often limited to classical or semi-classical physics and, therefore, are not suitable for modeling strictly quantum effects. A recent example of a quantum TCAD tool was put forward by Gao et al., who developed QCAD as a TCAD device simulator that accounts for quantum effects in silicon quantum dot devices [13,14]. The work of Gao et al. not only demonstrates the need for specialized TCAD tools that address quantum computing devices, but also serves to clarify the diversity of quantum technologies available requiring such tools.
The paper is organized as follows: We first present an overview of the silicon donor device model for quantum computing in Sec. 2 followed by the development of a computational workflow to design silicon donor qubits in Sec. 3. We then describe a complete implementation of the design workflow in Sec. 4 with details about the electrostatic and quantum chemistry solvers in Sec. 4.1 and Sec. 4.2, respectively, and details about the gate operation model in Sec. 4.3. Our final remarks are presented in Sec. 5.

Silicon Donor Qubit Models
Among the many different approaches to qubit fabrication, donors in silicon represent a promising path for encoding and manipulating qubits of information [15]. As originally proposed by Kane, the idea is to encode quantum information into the nuclear spin state of a dopant atom like phosphorous ( 31 P) that is shallowly embedded in a silicon lattice composed of 28 Si [16]. Because this silicon isotope has zero nuclear spin, the spin-1/2 phosphorous atom is ideally isolated from other spin contaminants [17]. Several works that followed the Kane proposal have also outlined how to realize qubits using the unpaired donor-bound electron, opening up the potential for easier addressability and faster gate operation [10,18,19,20]. Recent experiments using isotopically purified 28 Si have demonstrated single-qubit coherence times (T 2 ) over 30 seconds (nuclear spin) and 0.5 seconds (electron spin) at cryogenic temperatures [9].
A silicon donor qubit can be electrically manipulated by addressing the donor electron wave function with external gate electrodes. The donor electron wave function serves as a 'handle' for addressing the spin states through the hyperfine interaction, which is directly proportional to the spin density at the donor site [21]. The spin states of neighboring donors are coupled via the electron exchange interaction, which can be controlled by an external potential. Recent experiments have realized basic controls for phosphorus qubits in silicon (Si:P), however, a multi-qubit testbed is yet to be realized [7,8,9,22,23,24,25,26].
A key requirement for donor-based quantum computing architectures is the manipulation of individual donor electron wave functions over well-defined regions of space and time in the presence of material variations, finite temperatures, environmental noise, and electronic fluctuations [15]. The donor electron wave function oscillates several times across the interaction region between qubits, due to interference effects arising from the degeneracy of the 6-valleys in the conduction band of silicon. The exchange coupling (J) between donors is sensitive to these oscillations, which occur with a periodicity on the order of the silicon lattice constant [27]. Uncertainty in the position of the donor atoms therefore leads to unpredictable exchange couplings. Given the variability in the placement of the donors [28], the external potential must be tuned sufficiently to obtain the desired value of J. In addition, the donor position needs to be well-defined relative to the surface electrodes that lie along an atomically rough boundary. These process errors complicate control of the resulting qubit. Moreover, they are in addition to the decoherence caused by external noise, defects, spectral diffusion, and charge buildup at the electrode interfaces, which undermines the desired dynamics [29].
Device sensitivity to the electronic structure details imposes a formidable challenge to fabricating silicon donor qubits. But this challenge is not insurmountable as evidenced by increasing precision in sample fabrication [30], reduced decoherence from isotopic purification of the silicon substrate, and recent efforts to directly visualize the donor electronic state [31]. Very recently, Laucht et al. realized a single-qubit device from a lone 31 P atom embedded in isotopically pure 28 Si [25]. In particular, the electron and nuclear spin states of the P donor were manipulated using nanoscale gates to perform single-qubit logic gates and measure experimental signatures of coherence. Accurate characterization of the device in Ref. [25] required precise knowledge about the donor, including parameters such as the exact donor position, material strain and electrostatics in the vicinity of the donor. This demonstration underscores the need to better understand silicon donor physics at both a material and design level in order to improve device behavior.

Computational Workflow
The sensitivity of silicon donor physics to atomistic details, including single-atom defects, charge noise, and inter-donor spacing, emphasizes the significance of using microscopic models to guide device development. However, these models challenge existing TCAD tools because they require material models and simulation techniques that incorporate quantum physics. In addition, device models for silicon donor qubits are sufficiently distinct from existing CMOS-based transistors and integrated circuits as to warrant a dedicated design framework. In this section, we develop a computational work flow sufficient to express the behavior and function of state of the art silicon donor devices.
Our computational workflow for designing silicon donor qubits makes use of a multi-stage model that synthesizes intermediate solutions into a complete device representation. An overview of this workflow is shown in Fig. 1, in which electromagnetic and material input models are first simulated to recover the electrostatic field and electronic structure, respectively. These results are synthesized into a field-dependent model of the donor electron wave function. The synthesized model represents a silicon donor qubit design parameterized by the value of the electrostatic field. In the next stage, we use this model to simulate changes in the donor electron wave function induced by varying the applied electric field. The results of these simulations provide the necessary atomic details to describe the donor electron wave function, the hyperfine splitting, and electron exchange interaction. The final state of the workflow uses these atomic details to construct a gate operation model that describes the behavior of the silicon qubits driven by specific voltage (field) sequences. The output of the gate operation simulation provides the time-dependent electronic and nuclear spin states induced by the applied electric field. In particular, these spin states describe the expected gate fidelities and estimate the system decoherence time.
The required input to the workflow is a physical device layout defining the electronic circuit geometry, voltages, sources and drains as well as a partitioning of different material regions in the device. These material regions define the physical properties needed to calculate the electric field that exists at the donor location. An example physical layout is shown in Fig. 2, which is a model of the device used in recent experimental demonstrations by Laucht et al. [25]. Material properties for the surfaces and layers of this device must be specified to facilitate subsequent simulation of the electric field near the donor system. The specification also includes the operating temperature for the device, which may be on the order of sub-Kelvin temperatures for silicon donor qubits.
In addition to the electromagnetic properties of the layout, the workflow also requires a microscopic description of the silicon lattice in the region of the donor atom. This is necessary to simulate the electronic state of the donor atom. Simulations of the electron wave function may make use of different levels of theory, and these choices will impact the format by which the material model should be presented. A common format is the crystal geometry of the sample materials including the location of the donor atom. Any input must also account for the relative position of the donor with respect to the electrode layout including its material layers. For example, the position of the donor atom relative to the surface of the device determines the amount of strain present in the silicon due to the lattice mismatch of the substrate with the silicon dioxide and metal electrode layers.
An important consideration for current qubit development efforts is understanding how variations in material properties, environmental conditions, and electrode design impact the utility of silicon donor qubits for quantum logic operations. In particular, donor species and placement as well as applied fields influence the behavior of the hyperfine coupling, the spin states and the resulting gate fidelities. Consequently, the input models to our computational workflow may have a probabilisitic component that defines a priori distributions for the fluctuating physical parameters. For example, the relative position of the donor atom may be known with limited accuracy, as found in recent metrological studies [32]. The variability of positions can be accounted for by considering a statistical distribution, and the results may then be statistically combined to obtain an average fidelity.

Qubit Design
As shown in Fig. 1, the electric and material input models pass through initial simulation stages that generate, respectively, the electric field defined across the entire device and the electron orbitals for the unperturbed donor atom in the presence of material defects. The methods that are available for these simulations steps are varied, and we describe the expectations for the electromagnetic and material simulations in Sec. 4.1 and Sec. 4.2, respectively. In our workflow, the output from these models are synthesized together by first extracting the electric potential near the donor atom in the material sample. In Sec. 4.3, we describe how this model can be used to simulate gate operation with atomistic detail.

Electrostatic Field Simulations
The first stage of the workflow requires a specification of the physical layout of the device. An example of this specification is shown in Fig. 2, where a physical layout based on the experimental device tested by Laucht et al. is presented [25]. This particular device consists of multiple DC and AC electrodes used to tailor the electromagnetic environment of the donor atom embedded in the underlying silicon lattice. In general, the layout must account for the bulk material properties of the aluminum electrodes, the SiO 2 interface layer and any intrinsic charge regions. This requires assigning material properties to each region of the layout, including the silicon, silicon oxide, aluminium, and vacuum regions. In addition, the input specification must also assign voltages to each (DC) electrode. Although the source and drain regions are not visible in Fig. 2, they are included in the model as well as additional electrodes are connected to ground the substrate.
Simulating the electromagnetic field corresponding to the physical layout requires solving Poisson's equation over the entire device model. This is typically accomplished by first meshing the entire device with a resolution sufficient to accurately represent Figure 2. A physical layout representing the model electrodes and underlying substrate. The model presented here is based on the recent device used in Ref. [25] changes in the underlying potential. As shown in Fig. 3, meshing defines the points at which the potential is calculated and these points must be of sufficient density such that the boundary conditions between different materials satisfy the electromagnetic field equations. The process of finding the optimal mesh depends on the geometry of the device, although many implementations of electrostatic solvers provide options for automated mesh generation. For the physical layout shown in Fig. 2, we obtained converged solutions for the field near the location of the donor using an extremely fine mesh created with mesh spacing restricted to 0.1 nm for a 3x3x3 nm 3 block.
We have tested several commercial solvers for solving the electrostatic field of the prototype device given in Fig. 3. We have found the commercial electrostatic solver packages from COMSOL to be capable of modeling and simulating the electrostatics of these complex nanostructure geometries. An example of the resulting electrostatic field produced is presented in Fig. 4, which demonstrates the degree of expected field variability across a mesoscopic slice of the device. Regions in close proximity to the addressing electrodes can have relatively large electric field derivatives. These regions are of specific interest for controlling the electron wave function of the donor atom.
However, we anticipate that a challenge for these types of electrostatic field simulations is the degree of meshing needed. Whereas the solution to the electrostatic fields depends on satisfying boundary conditions over the complete nanoscale structure, there is also a requirement to provide sub-nanometer spatial resolution in the region of the donor. This is because the best basis functions for calculating the donor electron wave function will typically represent well-localized atomic orbitals. Because these orbitals are defined on much smaller length scales than the rest of the device, high density meshes will be needed in regions close to the donor location. Adaptive meshing alleviate some of the concern from increasing size, but when the location of the donor is variable, repeated simulations with different meshing will become a significant computation. However, sampling over the distribution of donor location will be an important measure of donor sensitivity to electrode design.

Donor Electron Wave Function Simulations
Quantum mechanical modeling of the silicon lattice and the embedded donor require solving for the electronic structure of the combined system. In particular, the electronic structure of the donor electron plays a bridging role that defines the effect of the applied electromagnetic fields on the time-dependent dynamics of the donor spin states. We now discuss the challenges related to electronic structure simulation for the donor electron wave function, especially with respect to the computational complexity and desired accuracy needed for subsequently describing gate operation.
A main challenge in the electronic structure simulation of the electron donor state is the computational complexity of the underlying methods. For example, density functional theory (DFT) is a workhorse of material science modeling that typically operates within the 100 to 1,000 atom regime on standalone workstations while implementations leveraging massively parallel computing resources can reach This plot shows the value of the field in the direction vertical to the plane of the device measured in V/µm. We use values of the simulated field taken from the highlighted square located under the right electrode. This corresponds with is the estimated location of the phosphorus atom and the area where increased resolution meshing is desired. up to 10,000 atoms. A 5 nm wide silicon crystal cube consists of 10,000 atoms and 140,000 electrons. This volume is relatively small compared to the physical layout in Fig. 3 but represents the high end of what is possible for DFT calculations. The root cause behind this practical limitation is that the computational cost for DFT typically scales cubically with system size due to the required linear algebra matrixmatrix operations.
An example of the scaling expected for these DFT calculations is shown in Table  1. Our example uses the NWChem suite [33] to calculate the electronic structure for Si:P nanocrystals up to 3 nm in size. Using Becke's three-parameter hybrid functional in combination with Lee-Yang-Parr correlation functional (B3LYP) [34,35], we performed tests of DFT scaling with various basis sets and degrees of computational parallelism. We typically used between 10 and 100 basis functions per atom, where the total number of basis functions determines the size of the Hamiltonian matrix. The electronic structure of the donor is then solved by iteration in a self-consistent field (SCF) process through diagonalization of the Hamiltonian. Typically up to 50 iterations (diagonalizations) were required. Our test used parallelized computations with up to 72 CPU nodes, where each CPU node has 2 Haswell processors (thus 24 cores per node). As shown in Table 1, over 1,728 CPU cores were used for several hours to calculate the electronic structure energy for the largest Si:) nanocrystal considered here (diameter 3 nm).
As expected, DFT calculations are not tractable with increasing system size. The development of models with reduced complexity for electronic structure calculations is therefore necessary for applying computational chemistry to models of silicon donor devices. One approach for reducing the computational complexity is to simplify the theory by neglecting core electrons that are typically not involved in chemical bonding. Another simplification is to leverage the locality of interactions by ignoring integrals between well-separated orbitals. Tight-binding (TB) [36,37] and densityfunctional tight binding (DFTB) [38] methods that make use of these approximations are computationally attractive alternatives for electronic structure modeling of silicon donor qubit devices. Whereas TB methods have been used extensively for modeling silicon donor systems [31,39,40,41,42,43], DFTB methods have not yet been extensively applied.
The relative error acquired by adopting different theoretical models for solving the electronic structure of the donor atom play an important role in simulating silicon qubit device physics. For example, a central physical quantity in the model of the donor atom as a qubit is the value of the hyperfine splitting (HFS). In principle, the HFS is the leading term in the nuclear spin-electronic interaction for the Si:P qubit. The isotropic HFS, also known as the Fermi contact interaction, is proportional to the spin density on the donor atom, where the spin density is defined as the difference between the densities for the spin up (ρ α (r)) and spin down (ρ β (r)) electronic states. However, this difference is often approximated by the density of the unpaired electron |φ(r)| 2 that is bound to the donor atom. Therefore, the isotropic HFS is given as where µ n is the nuclear magnetic moment, µ e is the electron magnetic moment, and r 0 is the donor position.
As an example, Fig. 5(a) shows the highest occupied molecular orbital |φ(r)| 2 for a Si:P nanocrystal obtained using DFT simulations and Fig. 5(b) shows the influence of an applied electric field on the differential distribution. An analysis of the HOMO density indicates that the radial distribution of the electron orbital closely matches that of a hydrogenic S type function that can be represented by a Slatertype orbital. The observed agreement is in contrast to the Gaussian type orbital expected for a quantum dot bound by a harmonic potential. In particular, the HOMO density decays very rapidly with respect to the distance from the donor, as shown in Fig. 5(c). These results suggest that although the shape of the orbital for the unpaired electron is complicated, it may be acceptable to use approximate the spin density by a much simpler, analytic form. This approximation can greatly simplify the electronic structure calculations for the donor atom in the presence of variable applied fields. For example, using a Slater-type orbital as an approximation to the actual donor electron orbital provides a convenient method for calculating how the HFS in Eq. (1) scales with applied electric field. In lieu of recomputing the exact electron orbital, which may be time consuming, a perturbation to the approximate analytical orbital can be used to estimate the corresponding HFS. The validity of this perturbative approach depends on the material environment and symmetry of the applied field. This approximation is most likely to be valid when the environment is highly symmetric. We now address how the theoretical difference between the spin density and the HOMO density may be handled and what trade-off may be gained between the accuracy of adopting an approximate orbital. The answers to these questions are obtained through the comparison of theoretically calculations and experimental results. As a specific example, we compare calculations of the HFS for Si:P nanocrystals using various theoretical methods against experimental measurements of the same quantity, as shown in Fig. 6. The experimental HFS value for P in bulk Si is known to be 42 G, while experimental values for Si:P nanocrystals embedded in insulating phosphosilicate glass matrices were previoulsy measured by Fujii et al. using electronic spin resonance (ESR) [44]. In those experiments, the Si:P nanocrystals have a diameter d in the range of 4.4-6.4 nm with a standard deviation of ∼1 nm. They were found to have a small number of deep dangling bond defects at the Si-SiO 2 interfaces which are filled by P doping as indicated by the increasing photoluminescence intensity at a low P concentration range and by the infrared optical absorption beyond a low P concentration threshold.
We make our comparison using several theoretical models that differ in accuracy and computational cost, as shown in Fig. 6. Hhybrid DFT calculations with B3LYP are very expensive, cf. Table 1, especially when the relativistic zeroth-order regular approximation (ZORA) is included [45,46]. DFT calculations within the local density approximation (LDA) performed by Melnikov and Chelikowsky [47] are expected to be computational more efficient due to the lack of exact Hartree-Fock exchange and the employment of pseudopotentials. These ab initio DFT methods can easily be afforded for the smallest nanocrystals with d in the range of 1-2 nm, with the nanocrystals modeled by spherical Si nanoparticles with H terminations on the surface to passivate the dangling bonds. However, for even larger diameters, the computational costs of DFT were too demanding. However, we found that the semi-empirical TB models implemented within the NanoElectronic Modeling in three dimensions (NEMO 3-D) code were capable of addressing these larger sized particles [36,37]. Using NEMO-3D, the nanocrystals were modeled by a cubic box with a zinc blende lattice of Si and a cube length of a in the range of 1-30 nm. The surface was not H-terminated but rather passivated by shifting the energy of the dangling bonds by 30 eV [48]. The Si atom in the center of the particles was replaced with a dopant P atom giving Si:P nanocrystals of different sizes. The geometries were not relaxed because it has been shown in the literature that there is no significant relaxation of the surrounding Si atoms when P serves as the dopant in a Si lattice [47]. However, doping with other group V elements such as As, Sb, or Bi may render a large geometrical distortion [49].
For the ab initio DFT calculations, the isotropic HFS values were calculated by the Fermi contact interaction between the electron spin and the nuclear spin at the dopant P using Eq. (1). For the semi-empirical TB method, the HFS was not explicitly calculated but instead found using a scaling method given a known HFS value for bulk Si as where A iso (bulk) = 42 G and we assume the bulk configuration corresponds with a = 30 nm. This scaling method is also used in the literature to study the Stark shift of the HFS in Si:P devices [39,50] A iso (E) where E is the electric field. As can be seen from Fig. 6, the experimental HFS values for Si:P nanocrystals with d in the range of 4.4-6.4 nm are 1.5-2.5 times larger than the bulk value and display a trend of increasing HFS with decreasing d. This strong particle size dependence was attributed to the quantum confinement of the unpaired electron spin in the nanocrystal. This trend is further verified by the ab initio DFT calculations of Melnikov and Chelikowsky using the LDA functional [47]. Our own B3LYP calculations for d=1-2 nm nanocrystals using Eq. (1) (1) is used to replace the scaling method, the HFS values would be lower than the bulk HFS value of 42 G beyond a particle size of 3 nm, due to the underestimated spin densities at the P nucleus. In addition, in the ab initio DFT calculations, spin-unrestricted open-shell wave functions were employed by giving spin-up and spin-down electrons different spatial orbitals. Consequently, the electron density |φ(r 0 )| 2 for the HOMO at the P nucleus are smaller than the spin density (ρ α (r 0 ) − ρ β (r 0 )) at the P nucleus by 10%. By comparison, the TB calculations were performed using spin-restricted open-shell wave functions and therefore the HOMO densities are exactly the same as the spin densities at the P nucleus. We expect that the HOMO densities in the tight binding calculations are good approximations of the real spin densities due to the small changes from the spin densities as can be seen from the ab initio DFT calculations. In short, the HFS values were found to depend on the level of theory, the size of the basis sets, and the spin restriction. Quantitative theoretical predictions require rigorous validations of theory which will be the focus of our future work.
A key challenge for fabricating these devices is controlling the location of the donor atom. In addition, there are no characterization methods that can pinpoint the location of the donor after device fabrication. State-of-the-art efforts reduce the uncertainty in the donor location using indirect position measurements, which effectively perform triangulation to narrow down the possible locations of the donor [32]. However, there are recent results that suggest the HFS at these locations may not be necessarily uncorrelated. In particular, the electron orbital wave function and HFS may have a weak dependence on donor translations along a plane parallel to the Si/SiO 2 interface [42]. By contrast, translations in the depth of the donor show an approximately exponential dependence for the HFS [42]. These relatively simple models for the dependence of the donor wave function and the associated hyperfine splitting on donor position suggests that sampling over the distribution of positions can be made efficient, thereby reducing the computational cost considerably.

Gate Operation Model
We now describe how the electrostatic and material models for silicon donor devices can be synthesized to generate a gate model describing how the nuclear and electron spin states of the donor system evolve in the presence of time-dependent changes to applied bias. In particular, we intend to define how this new model accounts for the instantaneous state of the donor and its changes in time. The instantaneous state is a fundamental description of the qubits represented by the electron and nuclear spin states that, within the context of quantum physics, provides all knowable information about the system. Therefore, we present a framework by which the complete information of the donor can be derived from the earlier workflow stages.
Our use case for this stage of the workflow is to simulate the behavior of a device whose physical layout and composition served as input to the previous stages. Gates within this single-donor model correspond to well-characterized (continuous) sequences of the applied bias that tune the hyperfine and Zeeman level spacings. An externally applied AC magnetic field then induces transitions between these levels that lead to changes in the joint electron-nuclear state. A gate operation model defines the fidelity with which the prepared state compares to the intended results. Notably, the gate operation model is synthesized from the electrostatic fields and electron orbitals generated by the earlier stages of the computational workflow. This requires the model to account for the time varying magnetic field and the varying hyperfine coupling that drives the dynamics of the electron and nuclear spins. Its ultimate purpose is to provide insight and feedback into the control parameters that are necessary to demonstrate gate operation at a desired fidelity.
A single-donor qubit gate operation model is defined by the interactions between the electron and nuclear spin states and the applied fields described by the timedependent Hamiltonian (4) H(t) = A iso (t)S · I + γ e B 0 S z − γ n B 0 I z + γ e B ac cos(ωt)S x − γ n B ac cos(ωt)I x where S is the electron spin operator and S x,z is its x, z-component, I the nuclear spin operator and I x,z is its x, z−component, B 0 is the applied DC magnetic field, and B ac is the amplitude of the oscillating magnetic field used to drive the spin qubits at frequency ω. The constants γ e = 28.025 GHz/T and γ n = 17.235 MHz/T are the gyromagnetic ratios of the electron and nucleus respectively. In the regime where (γ e + γ n )B 0 A, the eigenstates for the combined electron-nuclear spin state are well represented by the four-level energy diagram shown in Fig. 7. For example, in the device demonstrated previously by Laucht et al. [25], the time-dependence of the Hamiltonian is reflected in (a) the time-varying hyperfine coupling A(t), which is controlled through changes to the applied bias, and (b) the oscillating magnetic field used to drive the spins. Transitions between the nuclear and electron spin states can be induced when A is tuned such that the resulting energy gap matches the frequency ω of the applied magnetic field B ac .
The spin dynamics of the system can be obtained by solving the Schrodinger equation, or its noisy equivalent, the Master equation. For either the Schrodinger equation or its noisy equivalent, the state of the donor system is recovered by integrating over the time-dependent hyperfine coupling and magnetic field captured by Eq. 4. However, it is apparent from the gyromagnetic ratios that the time scale . Donor electron and nuclear spins states with their relevant energy separations not drawn to scale. Transitions between the electron (|↑ , |↓ ) or nuclear spin states (|⇑ , |⇓ ) are induced when A iso is tuned such that the energy gap between the spin states matches the frequency ω of an externally applied oscillating magnetic field. We assume γnB 0 < A iso /2.
for electron dynamics is much shorter than the time scale for nuclear spin dynamics, i.e., approximately three orders of magnitude. Consequently, it is essential to solve the equations of motion with sufficient temporal resolution. In particular, the much shorter time scale characterizing the electron spin must be adequately sampled to accurately simulate the complete joint spin state. This situation represents a simplified example of the multi-scale physics that arises in the simulation of heterogeneous coupled quantum systems. For our discussion of a single-qubit gate operation, the four-level system represented by Fig. 7 is relatively straightforward to solve using direct diagonalization of the time-dependent Hamiltonian. Consequently, the most significant computational complexity stems from integrating the underlying Schrodinger equation, cf. discussion above. Other methods can be applied to solving the differential equations defining changes to the qubit state. Finite element methods as well as finite differencing can be used to alleviate some of the problems with multiscale resolution. There are a variety of numerical packages available for performing any of these calculations, including general commercial solvers like MATLAB and Mathematica as well as open source solutions like the QuTiP Python module.
An example of how the gate operation model can be used to simulate the performance of a silicon donor qubits is shown in Fig. 8. This figure demonstrates preparation of a Bell entangled state using a sequence of NMR and ESR transitions for the joint electronic and nuclear spin state. In particular, we transform the joint spin state as We realize this operation by applying a sequence of controlled rotations for the encoded Figure 8. Transient fidelity during simulated operation of creating an entangled Bell state between the electronic and nuclear spin states of a single donor atom. The system is initialized into the joint spin state |↓, ⇓ . A conditional π/2-rotation is applied to the nuclear spin state using NMR with frequency ω = A iso /2 + γnB 0 = 477.52 Mrad/s. This is followed by a conditional π-rotation of the electron with ESR using ω = γeB 0 + A iso /2. The split x-axis changes resolution due to the distinct time scales required by the nuclear and electronic transitions. The simulated fidelities assume that the phase between the states are tracked and corrected accordingly. Prior to the rotation sequence described, A iso must be tuned to 117.6 MHz by the gate electrodes in the nanostructure such that ω is on resonance with the transition frequencies shown in Fig. 7.
qubits. We begin by assuming the electron-nuclear spin system is initialized as |↓⇓ , which can be achieved using the techniques mentioned in Refs. [7,8]. We then apply a π/2-rotation to the nuclear spin state conditioned on the electron being in the state |↓ . This rotation is implemented using the I x interaction found in Eq. (4) with B ac = 0.1 mT and ω = A iso /2 + γ n B 0 . The rotation takes approximately 65 µ s. The electron spin state is then conditionally rotated using the S x interaction and a second applied magnetic field B ac = 0.1 mT and ω = γ e B 0 + A iso /2. This rotation requires 360 ns to complete. Figure 8 illustrates that perfect fidelity is obtained when using these control conditions for an exactly characterized four-level system. A more significant burden arises when modeling realistic devices that are marked by uncertainty in the device fabrication or operation. In particular, the statistical variability in a diverse range of parameters such as the donor position and other material imperfections as well as electromagnetic field noise and uncertainties in device fabrication must be taken into account. Whereas classical averaging of the quantum state representation can be performed by integrating over a known noise distribution, this may be approximated numerically by drawing samples for Monte Carlo analysis. However, the influence of these noise sources on the hyperfine coupling and, more generally, the electron orbital, are not directly evident. Instead, these noise models require averaging over the output of the electrostatic and materials models described in Secs. 4.1 and 4.2. The need to repeat these calculations several times serves to greatly increase the computational complexity required for accurately simulating gate operations. As noted in Sec. 4.2, the dependence of the hyperfine splitting (HFS) with location may be estimated based on expectations for the behavior of the electron wave function. Given a model for the distribution of donor location, the predicted insensitivity to lateral displacements and the near exponential decrease with depth can provide a more computationally efficient model for the distribution of possible HFS values. Accordingly, this analytic distribution for HFS can be used in place of repeated electronic structure calculations.
Our model spin Hamiltonian in Eq. (4) has the electron and nuclear spin states coupled via the hyperfine interaction A iso (t). The hyperfine coupling is time dependent since it has to be tuned appropriately for gate operations, such that the frequency ω of the oscillating magnetic field matches the relevant energy separation between the different eigenstates in Fig. 7. We highlight that the variation of A iso (t) should be adiabatic, such that the excited donor orbital states are not populated during the process. Our computational chemistry methods described in the previous section do not permit for solving the time-dependent variation of A iso (t) to quantify the adiabaticity.
To maintain adiabaticty, the required timscale for varying A iso (t) should be much longer than h/E, where h is Planck's constant and E is the energy separation between the donor ground and first excited states. In realistic nanostructures, E takes a value of ∼ 10 meV [32]. This implies that the hyperfine coupling should be varied at timescales much longer than ∼ 0.4 ps to ensure negligible population of the excited donor states. A recent experiment indeed showed that the hyperfine coupling could be tuned over several microseconds to demonstrate single-qubit operations [25]. While high-frequency noise in the gate voltages can also potentially cause non-adiabatic transitions, we emphasize that these experiments typically use attenuation filters (beyond MHz frequencies) on the gate electrodes that limit high frequency noise. Similarly, the noise that arises from intrinsic charge fluctuations in nanostructures typically have an inverse dependence on frequency [9], such that any contributions at terahertz frequencies or greater would be substantially smaller. As a result, we expect that it is reasonable to assume the electron orbital quickly equilibrates during the rise time of the applied voltage used to control A iso (t).
There is some uncertainty as to whether the dynamics arising from electronphonon coupling and fluctuations arisign from thermal motion of the nuclei in the nanostructure may impact A iso (t). While the low temperature ∼ 100 mK operation of such devices may rule out this possibility, it is possible to further investigate this situation by solving for the dynamics of electron orbital wave function using the timedependent Schroedinger equation (TDSE). Generally, DFT and TB approaches allow for directly simulating excitation of the donor wave function in response to timedependent fields/perturbation and thermal motion of nuclei. For silicon donors, the dynamics of the donor electron is especially significant since it serves as the handle by which logical operations are carried out, and the simulation of real-time electron dynamics using DFT is a promising tool. This requires solving the TDSE for electrons in the von Neumann density matrix representation as where ρ e (t) and H e (t) stand, respectively, for the time-dependent electronic density and electronic Hamiltonian, and the right hand side is the commutator [H e , ρ e ] = H e ρ e − ρ e H e . The Hamiltonian H e (t) now describes mean field interactions of the donor electron with the nuclei, other electrons and external fluctuations. The interaction of the materials with external fields from the electrodes or other sources is included by adding time-dependent external potentials to the Hamiltonian. The solution to Eq. (6) is usually obtained from diagonalization of the Hamiltonian H(t). This may be achieved by the formation of complex-valued exponential time-evolution operator, from the eigenstates of the H e (t )-matrix, where t ∈ [t, t ]. The evolution operator U (t, t ) is then directly applied to the density matrix ρ(t) and the new (propagated) density matrix at time time t is given by ρ e (t ) = U e (t, t ) † ρ e (t)U e (t, t ), where U (t, t ) † stands for a Hermitian conjugation of U (t, t ). It is possible to implement the time-propagation using a first-order Magnus expansion and iterative diagonalization [51,52,53,54]. This procedure for the time dependent electron dynamics in conjunction with tight-binding DFT has been proved to be very efficient for molecular systems consisting thousands of atoms. The computational scaling is comparable with standard ground state calculations while allowing inclusion of nuclei motion and beyond linear response effect. The set of coupled equations given by Eqs. (1), (4) and (6) capture the timedependent dynamics of the electron and nuclear spin states. Whereas the timedependence of the electron wave function can be modeled explicitly using Eq. (6), the corresponding HFS can be calculated and substituted into the spin Hamiltonian given by Eq. (4). The complete dynamics of the donor systems can be simulated. However, the large discrepancy in time scales for the electron and the spin dynamics relevant to logical computation make it clear that the granularity for these simulations must interpolate between these two extremes. Foremost, it is only necessary to simulate the electron dynamics on those timescales that induce changes. This includes the rise time for the applied voltage and the correlation time for the electronic noise. Outside of these timescales, we may approximate the HFS value by a constant.
Ultimately, our goal of modeling gate operation in realistic silicon donor devices has emphasized the coupled nature of electromagnetic and material modeling. These simulations produce numerical representations of the instantaneous quantum state that can be used to compare against the expected state dynamics and experimental observables. The resulting fidelity, as measured by the overlap between the observed and expected state values, is a useful quantity for assessing the quality of a given design, material sample, or control sequence. In addition, the gate operation model can eventually be validated against experimental samples by making comparison between measured fidelities and those calculated under the model. However, the importance of high-fidelity models and accurately dynamical simulation must be stressed, since these models can then provide feedback for the improved design of new devices.

Conclusions
Given recent advances in the demonstration of prototype quantum computing device, new tools are needed for modeling and simulation of the fabricated qubits. Silicon donor qubits offer an example of how a multi-staged modeling workflow is needed to synthesize electrical and quantum physical properties into a single consistent model. Our multi-physics approach integrates models of the electrostatics, quantum chemistry, and device operation into a common framework. We must caution that the scaling of direct simulation is very unfavorable for quantum mechanical systems. In particular, we do not expect DFT calculations to remain tractable for more than two donor atoms separated by more than a few nanometers. However, alternative approaches like TB and DFTB will be useful for simulating larger systems. We believe that the modeling and simulation workflow presented here for silicon donor devices will provide valuable insights into the ongoing development of quantum computing devices that require atomistic control. We expect the silicon donor TCAD tool to be essential for maturing silicon donor systems into functioning qubits, including integration into the device fabrication process.