Optically Controlled Entangling Gates in Randomly Doped Silicon

TL;DR

This paper investigates the density and arrangement of dopants in randomly doped silicon to optimize entangling gates for quantum computing, combining theoretical models, simulations, and experimental predictions.

Contribution

It introduces a comprehensive approach combining Poisson point process theory, Monte Carlo simulations, and cluster expansion to optimize dopant densities and structures for quantum gates in silicon.

Findings

Calculated doping densities for maximal entangling gate densities.

Validated theoretical results with Monte Carlo simulations.

Predicted control of one dopant species via orbital excitation of another.

Abstract

Randomly-doped silicon has many competitive advantages for quantum computation; not only is it fast to fabricate but it could naturally contain high numbers of qubits and logic gates as a function of doping densities. We determine the densities of entangling gates in randomly doped silicon comprising two different dopant species. First, we define conditions and plot maps of the relative locations of the dopants necessary for them to form exchange interaction mediated entangling gates. Second, using nearest neighbour Poisson point process theory, we calculate the doping densities necessary for maximal densities of single and dual-species gates. We find agreement of our results with a Monte Carlo simulation, for which we present the algorithms, which handles multiple donor structures and scales optimally with the number of dopants and use it to extract donor structures not captured by our Poisson point process theory. Third, using the moving average cluster expansion technique, we make predictions for a proof of principle experiment demonstrating the control of one species by the orbital excitation of another. These combined approaches to density optimization in random distributions may be useful for other condensed matter systems as well as applications outside physics.

Figures

Tables

Table 1. Table 1: Values for the minimum/maximum radius between dopants occurring when their interaction exchange energy is equal to J dec subscript 𝐽 dec J_{\mathrm{dec}} , defined in Eq. 3, plotted in Fig. 2 for the operation of the optical entangling gate within the decay time. *Bilayer case - calculated for readouts (Si:As) in a plane 13.2nm from and parallel to the control (Si:P) plane.

Polarisation	Interaction =$J_{\mathrm{dec}}$	Name	Distance
axis of light			$\left(nm\right)$
In plane	J(As1s-As1s)	$R_{rr}$	11.0
	J(P1s-As1s)	$R_{min}$	11.4
	J(P2p-As1s)	$R_{max}$	17.9
	J(P2p-P2p)	$R_{cc}$	42.2
	J(P1s-P1s)	$R_{min}^{\prime }$	11.8
Out of plane	J*(P2p-As${}_{1s}{}^{})$	$R_{max}$*	10.2
	J(P2p-P2p)	$R_{cc}$	28.5

Table 2. Table 2: Total doping densities giving maximum number of dopants in 3D, and in 2D isolated by R. *radius corresponding to laser light polarisation out of plane.

	3D (${10}^{15}$/cm3)	2D (${10}^{10}$/cm2)
	R = 42.2 nm		R* = 28.5 nm
Total density	3.18	1.79	3.91
Isolated density	1.17	0.66	1.44
Isolated fraction	1/e

Table 3. Table 3: Analytical (matching to simulation) results for total densities of control (P) and readout (As) dopants and densities of readouts which are part of gates (D qubits ). Densities are in dopants per cm 2 (in 2D) or cm 3 (in 3D).

•	•	$D_{qubits}$	%	$D_{As}$	$D_P$
2D	P-As-As	$6\times {10}^8$		$1.5\times {10}^{11}$	$1.8\times {10}^{10}$
•	P-As	$1.2\times {10}^9$	1.5	$8\times {10}^{10}$	$1.8\times {10}^{10}$
•	P-P	$5.4\times {10}^9$	20	0	$2.9\times {10}^{10}$
3D	SFG	$2\times {10}^{14}$		$7\times {10}^{16}$	$3.2\times {10}^{15}$
•	P-As	$3\times {10}^{14}$	0.8	$3.5\times {10}^{16}$	$3.2\times {10}^{15}$
bilayer	P-As-As	$6\times {10}^8$		$3\times {10}^{11}$	$3.9\times {10}^{10}$
•	P-As	$3.5\times {10}^9$	1.8	$2\times {10}^{11}$	$3.9\times {10}^{10}$

Equations

$$\hat{\textit{u}}_{phase}=\hat{r_{2}^{z}}\big{(}-\frac{\pi}{2}\big{)}\hat{r_{1}^{z}}\big{(}\frac{\pi}{2}\big{)}\hat{e}_{12}\big{(}\frac{\pi}{2}\big{)}\hat{r_{1}^{z}}\big{(}\pi\big{)}\hat{e}_{12}\big{(}\frac{\pi}{2}\big{)},$$

$$J=\frac{h}{4t}\mathrm{or}t=\frac{h}{4J}$$

$$J_{\mathrm{dec}}=\frac{h}{4T_{\mathrm{dec}}}.$$

$$D_i\equiv D_t{P}_0(R)=D_t{e}^{-V(R)D_t},$$

$$D_t(D_{imax})=\frac{1}{V(R)}.$$

$$D_{imax}=D_t{e}^{-1}.$$

$${\int }_{R_{min}}^{R_{max}}d{r}_1{\int }_{r_1}^{R_{max}}d{r}_{2}S(r_1)S(r_2)D_r^2{e}^{-V(R_{max})D_r}.$$

$$e^{-\delta D_r}.$$

$$D_{\mathrm{sfg}}(D_c,D_r)=D_c{e}^{-V(R_{cc})D_c}\times$$

$${\int }_{R_{min}}^{R_{max}}d{r}_1{\int }_{r_1}^{R_{max}}d{r}_{2}S(r_1)\tilde{S}(r_2)D_r^2{e}^{-(V(R_{max})+{\delta }_{SFG})D_r}$$

$$D_{Heis.\mathrm{ex}-\mathrm{gd}}(D_c,D_r)=D_c{e}^{-V(R_{cc})D_c}\times$$

$${\int }_{R_{min}}^{R_{max}}drS(r)D_r^2{e}^{-(V(R_{max})+\delta )D_r};$$

$$D_{Heis.\mathrm{ex}-\mathrm{ex}}(D)$$

$$|S_i⟩=\{\begin{array}{ll}|\uparrow ⟩& \mathrm{if}i=\mathrm{readout}\\ |\downarrow ⟩& \mathrm{if}i=\mathrm{control}\end{array}$$

$$|{\psi }_0⟩=|S_1{S}_2...S_n⟩.$$

$$\mathcal{H}=$$

$$+\frac{1}{2}\sum _{i\ne j}{J}_{ij}^{PP}{\text{S}}_i^P{\text{S}}_j^P,$$

$$⟨S^{z\mathrm{As}}(t)⟩=⟨{\psi }_0⟩\frac{1}{N}\sum _i{S}_i^{z\mathrm{As}}(t){\psi }_0.$$

$$P_{\mathrm{sf}}^{\mathrm{As}}(t)=$$

$$=$$

$$P_{\mathrm{sf}}^{\mathrm{As}}(t)\to \{\begin{array}{ll}\approx 0& \mathrm{for}\mathrm{control}\mathrm{in}1s\\ \ne 0& \mathrm{for}\mathrm{control}\mathrm{in}2p\end{array},\mathrm{for}\mathrm{late}\mathrm{times}\mathrm{t}.$$

$$J=\int d{\mathbf{r}}_{1}d{\mathbf{r}}_2$$

$$\times {\psi }_2({\mathbf{r}}_1-\mathbf{R}){\psi }_1({\mathbf{r}}_2),$$

$${\psi }^{(j)}(\mathbf{r})=\sum _{\mu }{F}_{\mu }^{(j)}(\mathbf{r}){\varphi }_{\mu }(\mathbf{r}),$$

$$F_{+z}^{(1)}(\mathbf{r})=\frac{\exp \left[-{\left(\frac{x^2+y^2}{(\alpha a_1{)}^2}+\frac{z^2}{(\alpha b_1{)}^2}\right)}^{1/2}\right]}{\sqrt{6\pi (\alpha a_1{)}^2(\alpha b_1)}},$$

$$F_{+z}^{(2)}(\mathbf{r}$$

$$F_{+z}^{(3)}(\mathbf{r}$$

$$J=2\sum _{\mu ,\nu }{j}_{\mu ,\nu }\cos ({\mathbf{k}}_{\mu }.\mathbf{R})\cos ({\mathbf{k}}_{\nu }.\mathbf{R}),$$

$$j_{\mu ,\nu }=\int d{\mathbf{r}}_{1}d{\mathbf{r}}_2$$

$$\times F_{2,\nu }({\mathbf{r}}_1-\mathbf{R})F_{1,\mu }({\mathbf{r}}_1).$$

Full text

I Introduction

The ability to perform coherent quantum operations on large collections of quantum bits (qubits) could lead not only superior scaling to classical machines on date-based algorithms such as factoring Shor (1994), unstructured search Grover (1996) and machine learning  Biamonte et al. (2017) but may also lead to a better understanding of strongly correlated systems through Feynman’s idea of quantum simulation Feynman (1982); King et al. (2018); Braumüller et al. (2017); Lamata et al. (2018). Initial successes in the realization of these applications in the last decades, mainly on AMO platforms such as trapped ions Häffner et al. (2008); Zhang et al. (2017) and ultracold atoms in optical lattices Greiner et al. (2002); Gross and Bloch (2017); Mazurenko et al. (2017); Bernien et al. (2017), have also shown how hard it is to reach a regime where quantum computers are large enough to outperform their classical counterparts either due to limitations in the qubit addressability, interactions or sheer quantity. This scalability challenge could be more easily overcome in solid-state realizations, where a wide range of qubits have been proposed including Majorana fermions in nanowires Sau (2017); Lahtinen and Pachos (2017), superconducting qubits Devoret et al. (2004); Barends et al. (2014); Kelly et al. (2015), nitrogen-vacancy centers in diamond Childress and Hanson (2013), quantum dots Mi et al. (2017), and donor impurities in silicon Hill et al. (2015); Morley (2015).

A realization in doped silicon, the most important material in the electronics industry today, could point the way forward for the widespread introduction of quantum computation.

Donor impurities have several spin 1/2 degrees of freedom which provide natural Hilbert spaces to encode qubits DiVincenzo and Loss (1998). The nuclear spins of 31P Kane (1998); Morton et al. (2008) and the electron spin of the valence electron of several species have some of the longest decoherence times of any qubit realizations ranging from hours to several days Pla et al. (2013) especially when benefiting from isotopically pure 28Si and/ or field-insensitive “clock transitions” Saeedi et al. (2013); Saeedi (2014). Initialisation and readout Morello et al. (2010) and single-qubit operations have been demonstrated on single electron spins  Pla et al. (2013), via metallic microwave strip lines Fuechsle et al. (2012) and donor-bound exciton spectroscopy Lo et al. (2015); Saeedi et al. (2015). However, in addition to initialisation, readout and single-qubit operations, a two-qubit entangling gate is required for universal quantum computation Lloyd (1996).

Several proposals to realise entangling gates in silicon exploit the exchange interactions between ground-state electronic states; these interactions are strongly oscillatory because of interference between the different valley components of donor wavefunction Koiller et al. (2001) and typically require atomic-scale precision in the positioning of the donors Kane (1998). This needs specialised lithography based on selective removal of hydrogen atoms in the scanning tunnelling microscope Hersam et al. (2000) and subsequent implantation of donor species  O’Brien et al. (2002), a delicate and resource-intensive techniqueO’Brien et al. (2002); Fuechsle et al. (2012); Broome et al. (2018). Schemes based on electric-dipole interactions have larger donor spacings and tolerate less stringent fabrication requirements Tosi et al. (2017). However, a scheme capable of inducing entangling interactions with truly random placement of donors would allow relatively seamless integration of the fabrication into the standard widespread silicon processing industry.

In this paper, we discuss a family of entangling gates, sketched in Fig. 1, that make use of physically localised or extended orbital excited states to control the interactions between donor qubits in a disordered ensemble Stoneham et al. (2003). These orbital excited states can be produced by excitation with terahertz radiation Greenland et al. (2010). The random locations of the donors provides inhomogeneous broadening due to differing local magnetic field provided by hyperfine spins which could enable selective addressing of sub-ensembles with a specific frequency of terahertz radiation Stoneham (2009), and the spread of frequencies could be further increased by placing the sample in a field gradient, leading to each donor having a different transition frequency Stoneham (2009). The excited states produced (for example, 2p orbitals) have a spatial structure that depends on the axis of polarisation of the laser, thus providing another tool for spatial selectivity within local configurations. Important steps have been made towards the realization of exchange interaction-based gates in donors in silicon. The disentanglement of the control particle from the two entangled qubits was shown to be feasible after the gate operation Rodriquez et al. (2004). Coherent control of the valence electron orbital states of dopants such as phosphorous, bismuth, antimony and arsenic has been demonstrated with terahertz light, using a Free Electron Laser tuned to the low tens of meV Greenland et al. (2010). Furthermore, initialisation and readout can be done optically, using donor-bound exciton spectroscopy in the terahertz frequencies which was demonstrated experimentally with phosphorous spins in silicon Morse et al. (2018). Until now it has remained unclear what densities of ‘viable’ entangling gates could be reached and even what the requirements are for a configuration to be ‘viable’.

The present paper deals only with the special case of homogeneous distributions, either throughout a 3D sample or on a 2D plane. These correspond to the ideal limits of a uniformly doped bulk sample or a perfect delta-layer respectively; we also consider the case where two parallel ideal delta-layers are implanted with different species. The general case of a spatially varying density (hence allowing both for graded doping of a bulk material and for the inevitable broadening of real delta-layers) is treated in a companion paper Stockbridge et al. ; in that paper a heuristic is also given which enables fast estimates of the densities of viable clusters over a wide density range.

This paper is organized as follows. First, we estimate appropriate conditions on the spacing of the dopants required for them to form a viable entangling gate. Second, optimal doping densities to produce the highest densities of entangling gates are calculated using a Monte-Carlo simulation and Poisson point process theory. Third, we propose a proof-of-principle experiment to show control over the exchange interaction between donors at the optimal doping densities and predict its results with the Moving Average Cluster Expansion (MACE) method Hazzard et al. (2014). Our numerical results apply to the case of Group V donors in silicon; however, our methods apply to any material where particular configurations of multiple impurities are required, and the conditions on the configurations needed can be expressed in terms of distances between different species.

II Geometric conditions on optically controlled entangling gate configurations

First, we define conditions on the distances between dopants to form an SFG or Heisenberg entangling gate by comparing energy scales corresponding to the lifetime of the excited state Vinh et al. (2008); Litvinenko et al. (2015), the strength of the entangling exchange interaction and the requirement of qubits to be isolated from each other when not part of the entangling operation. Based on these conditions we then calculate maps and line scans of the interaction exchange energy from a theory including multi-valley interference Wu et al. (2018) to delimit areas around donors in which entangling gate operation is possible.

II.1 Conditions on the exchange interaction energy

A phase gate (which can be combined with Hadamard gates to result in a Controlled NOT gate) Levy (2002) can be implemented by the sequence

[TABLE]

where the single qubit operations are given by the Zeeman rotations $\hat{r}_{i}^{\alpha}=e^{-i\theta\hat{H}_{i}^{\alpha}/gB^{\alpha}}$. The two qubit operation is given by $\hat{e}_{12}(\frac{\pi}{2})=e^{-i\frac{\pi}{2}\hat{H}_{12}/J}$, where $\hat{H}_{12}=J\hat{\mathbf{S}}_{1}\hat{\mathbf{S}}_{2}$ is the Heisenberg Hamiltonian between spins $\mathbf{S}_{1}$ and $\mathbf{S}_{2}$. Equating this with the time evolution operator $\hat{U}=e^{-\frac{i}{\hbar}t\hat{H}_{12}}$ leads to the condition

[TABLE]

for a successful phase gate operation.

The exchange interaction constant $J$ arises from overlap of spratially separated wavefunctions. The value of $J$ (or $t$) must be controlled such that Eq.(2) is fulfilled, without perturbing interactions between the gate dopants and other donors. In silicon this is especially challenging as the exchange interaction between valence dopant electrons not only decays exponentially with distance but also oscillates with a period of the order of the lattice spacing due to intervalley-interference Keyes (2005). The lifetimes $T_{\mathrm{dec}}$ of the donor excited states are limited by phonon-mediated decay, and for phosphorous and arsenic they have been measured as $T_{\mathrm{dec}}=200\,\mathrm{ps}$Vinh et al. (2008), which limits the gate operation time to $t<T_{\mathrm{dec}}$; hence, we define an energy scale

[TABLE]

which sets the minimum exchange energy needed for being able to perform a successful 2-qubit gate.

In the rest of this paper we will consider two group-V donor species: the shallow donor phosphorous (Si:P) as the controls (also referred to as c in the general case), whose excited state wavefunctions (P2p+-) mediate the interactions between the deeper arsenic (Si:As) ground state (As1s) spins as the readout qubits (also referred to as r if not specifically considering an implementation relying on using this particular species).

II.1.1 SFG entangling gate

We assume that only the control donors are excited to a higher-lying orbital state which we will choose to be the 2p excited state; the readout donors remain in the ground state. The following conditions then need to be met for successful gate operation:

Readouts

Interaction between readouts should be smaller than the decoherence interaction exchange energy defined in Eq.(3), ie. $J$(r1s-r${}_{1s})(r)<J_{\mathrm{dec}}$. We will denote the distance at which these two scales are equal as $R_{rr}$.

Readout-Control

If the control is in the ground (1s) state, the interactions between readout and control should also be small: $J$(c1s-r${}_{1s})(r)<J_{\mathrm{dec}}$, with equality at $R_{min}$. Furthermore, if the control is in the excited state, control and readout should interact, therefore $J$(c2p-r${}_{1s})(r)>J_{\mathrm{dec}}$ with equality at $R_{max}$.

Controls

If two neighbouring controls are excited to the $2p$ their interaction should not dominate the process, i.e. $J$(c2p-c${}_{2p})(r)<J_{\mathrm{dec}}$ with equality at $R_{cc}$.

The distances introduced above are shown in Fig. 2.

Please note that in the above considerations we did not need to specify the number of readouts per control. In fact, although originally devised for two readouts Stoneham et al. (2003), we will leave this number general.

II.1.2 Heisenberg entangling gate

There are two ways of implementing this gate. In the first, which we refer to as ‘excited-ground Heisenberg gate’, one donor is in the excited and the other in the ground state; this configuration is the same as the SFG gate with just one readout. In the second case, the ‘excited-excited Heisenberg gate’, both donors are excited. Therefore this can be a single-species gate composed of two Si:P electrons, for example.

In the excited-ground case, the conditions are the same as for the SFG gate, except that only one readout should interact strongly with a control, and not two.

For the excited-excited Heisenberg gate the dominant process should be the 2p+- state of both controls interacting with each other: J(c2p-c${}_{2p})(r)>J_{\mathrm{dec}}$ with equality at $R_{max}=R_{cc}$ and the entire gate should also be isolated from other donors by this distance. The interaction of ground-state controls in $1s$ should be small: J(c1s-c${}_{1s})(r)<J_{\mathrm{dec}}$ with equality at $R^{\prime}_{min}$.

Having defined the conditions that viable entangling gate configurations need to fulfil, we can now find the corresponding distances by calculating maps of the exchange interaction between Si:P and Si:As.

II.2 Constraints on dopant separations from conditions on exchange interaction

The exchange interaction energy calculated using the Heitler-London approximation Koiller et al. (2001); Wu and Fisher (2008) taking into account the multi-valley coupled wavefunction and central cell correction can be found in Appendix A. The maps in Fig. 2 show the regions in which $J>J_{\mathrm{dec}}$, for $J$ connecting different types of dopant species in various states.

For information to not remain localised in the quantum computer, it is essential that qubits can communicate with two or more entangling gates. This has the implication that the maximum separation between controls and qubits is larger than half of the distance separating two controls / gates. This is the case, as can be seen in Fig. 2.

To simplify the subsequent treatment of the dopant distributions, we approximate the interaction zones outlined by $J_{\mathrm{dec}}$ as circles or spheres; to do this, we determine discs whose areas are equal to those of the interaction zones. The disc radii corresponding to the monolayer doping in plane laser polarisation setup and the bilayer out of plane laser polarisation described in the next paragraph are displayed in Table. 1.

II.2.1 Bilayers

For two-species gates, we also consider separating the controls and the readouts into two separate parallel planes, sufficiently distant that the ground states do not interact on the timescale of the gate operation but sufficiently close that the qubits can interact with the excited state of the controls. For this to be optimised in a randomly doped sample, the readout layer must be at a distance from the control layer which maximises exposure to the area in which $J$(P2p-As${}_{1s})<$ $J_{\mathrm{dec}}$ and minimum exposure to the area in which $J$(P1s-As${}_{1s})>$ $J_{\mathrm{dec}}$, dark blue on Fig. 2. The small periodic oscillation put apart, the ground state interaction exchange energies are equal $J_{\mathrm{dec}}$ in a spherically symmetric way, and can be modelled as a sphere. As is clear from Fig. 4b), the optimal distance $d$ separating the two layers should be equal to $R_{min}$. For the 1s states of Si:As and Si:P this is 13.2 nm. Note that the areas which the qubit can occupy in its own layer correspond to a disc of inner radius 0 and outer radius $\sqrt{R_{max}^{2}\>-\>d^{2}}$, as can be seen in Fig.  4b). By polarising the light out of plane the controls can pack closer together and thus make more qubits viable thereby increasing the density of entangling gates.

III Maximising the entangling gate density

In a randomly doped sample, the only free parameters are the doping densities. In this section, we establish the optimal doping densities of both species such that the entangling gate density is maximal. Phosphorous atoms have a random spatial arrangement on the silicon surface resulting from exposure to phosphine gas Trappmann et al. (1997). The atoms cannot land in exactly the same location, thus their arrangement would correspond at best to a hard sphere Poisson point process, however we make the approximation that they correspond to a pure Poisson point process. Many methods, including analytical and simulation-based approaches, are available to treat Poisson point processes: Monte-Carlo simulations can treat problems with complex geometries and have the advantage of being flexible, but must be coded efficiently to scale well. Analytical approaches can supply results in closed form but rely on particular assumptions and are valid only if these apply. Analytical nearest-neighbour methods have been used to study distributions of $m^{th}$ nearest neighbours Torquato et al. (1990), isolated pairs of points Pickard (1982) and the probability of occurrence of reflexive nearest neighbours Cox (1981)Dacey (2010).

Here we employ a Monte Carlo simulation to count viable configurations (defined in the previous section). It scales linearly with the number of dopants, i.e. as $\mathcal{O}(n_{c}+n_{r})$, where $n_{c}$ is the number of control dopants and $n_{r}$ is the number of readout dopants. This allows large samples (up to a few million dopants per run) to be used routinely, leading to lower statistical errors. The algorithm used is described in Appendix B. We find agreement of these results with an analytical solution for a Poisson point process in both two and three spatial dimensions, using nearest-neighbour methods to determine the probability of occurrence of a nearest event within a given radius from a point chosen at random Bahcall J. N. (1981). We treat the two dopant species as separate independent sets of events occupying the same volume, with different densities.

First, we study the control dopant distribution and calculate the control density which gives the highest density of points separated by at least $R_{cc}$; we call the total density of such controls the viable control density.

Second, we find the cumulative probability distribution function (CPDF) describing the probability that a control chosen at random (we make the approximation that the isolated control dopants are randomly distributed) is surrounded by a configuration of readouts fulfilling the distance conditions in Table. 1; we seek the optimal readout density where the greatest number of controls have both 1st and 2nd nearest readout neighbours that are situated in the ‘viability shell’ between $R_{min}$ and $R_{max}$, are further than $R_{rr}$ from each other, and are also further than $R_{rr}$ from any other neighbour. We obtain a CPDF which depends on $R_{min}$ and $R_{max}$ and $R_{rr}$.

Third, we find the optimal density of viable entangling gate configurations by combining the two previous results, by making the approximation that both probabilities relative to the control and readout are independent.

We find slightly higher entangling gate densities for the bilayer configuration described in Fig. 4, for which each species is contained in a 2D layer separated by a small distance. In the Monte Carlo simulation, we are also able to identify configurations with more than two readouts within the viability shell surrounding the same control. These could be useful for implementing multi-qubit entangling gates such as the Toffoli gate Monz et al. (2009).

III.1 Control doping density: maximising the density of events isolated by a fixed radius

We seek to maximise the density of events isolated by at least a fixed radius $R_{cc}$. The probability of finding an event in a spherical shell of thickness dr at a distance $r$ from a randomly chosen event in the distribution of density $D_{t}$ is $\frac{\mathrm{d}V(r)}{\mathrm{d}r}D_{t}\,dr$, where $V(r)$ is the volume of the sphere of radius $r$ Scott and Tout (1989); Moltchanov (2012). The probability that there are no events closer to the randomly chosen event than $r$ is then $P_{0}(r)=e^{-V(r)D_{t}}$, from Bahcall J. N. (1981)Scott and Tout (1989). Optimising this probability leads to setting the total density of events to 1 (such that there is no chance of finding any other event within $r$), as can be seen on Fig. 3. However, optimising the density of points $D_{i}$ isolated by exactly $R$ from their neighbours does present a maximum (see Fig. 3).

[TABLE]

which is maximum when

[TABLE]

At this density, the fraction of points isolated by $R$ is $1/e$:

[TABLE]

The values for the total and isolated control densities corresponding to the control radii (Table. 1) are displayed in Table. 2.

III.2 Readout doping density

Let $S(r)=\frac{\mathrm{d}V(r)}{\mathrm{d}r}$ be the surface area of the $n$-dimensional sphere of radius $r$ (i.e. $S(r)=2\pi r$ in 2D and $S(r)=4\pi r^{2}$ in 3D).

If the readout density is $D_{r}$, the probability of finding the first and second nearest neighbours between $R_{min}$ and $R_{max}$ is

[TABLE]

However, the number of viable configurations is reduced by the additional requirement that the second nearest neighbour must be at least $R_{rr}$ from the first. If $r_{2}<(r_{1}+R_{rr})$, the sphere of radius $r_{2}$ defining the viable positions for the second nearest neighbour must therefore have a spherical cap removed from it (see Fig. 4) subtending an angle \alpha=\cos^{-1}\big{(}\frac{r_{1}^{2}+r_{2}^{2}-R_{rr}^{2}}{2r_{1}r_{2}}\big{)}. This gives rise to a new surface of smaller area $\tilde{S}(r_{2})=2(\pi-\alpha)r_{2}$ in 2D and $\tilde{S}(r_{2})=4(\pi-\alpha)r_{2}^{2}$ in 3D. In the bilayer case, $R_{min}=0$ which further constrains $\alpha$ to be $\pi$ if $r_{1}<(R_{rr}-R_{max})$, ie. if $\alpha$ is complex.

Finally, cases where further readout donors lie outside the sphere of radius $R_{max}$ but within a volume $\delta$ defined as being within radius $R_{rr}$ of the first or second nearest neighbour must be excluded; this corresponds to reducing the probability by a factor

[TABLE]

In the SFG case, if $\delta_{1}(r_{1})$ and $\delta_{2}(r_{2})$ are the volumes of the spheres of radius $R_{rr}$ centred on the first and second nearest neighbours which lie outside the sphere of volume $R_{max}$, and $\delta_{ov}(r_{1},r_{2},\theta)$ is the overlap between these two volumes where $\theta$ is the angle between the readouts, $\delta_{SFG}=\delta_{1}(r_{1})+\delta_{2}(r_{2})-\delta_{ov}(r_{1},r_{2},\theta)$. In both other gate types, $\delta=\delta_{1}(r_{1})$. This three circle overlap configuration and area calculation has been made into a Wolfram Demonstrations project Crane and the functions corresponding to the overlap of two and three circles case was taken from Fewell (2006).

We can now construct the total densities of viable configurations for our three types of entangling gate, for homogeneous doping with a control density $D_{c}$ and a readout density $D_{r}$.

• (i)

The density of SFG gate configurations is the product of the viable control density with the probability for a control to be surrounded by a successful configuration of readouts (the product of the CPDFs defined in equations 7 and 8):

[TABLE]

• (ii)

The density of two-species excited-ground Heisenberg entangling gates is:

[TABLE]

• (iii)

The density of single species excited-excited Heisenberg entangling gates is

[TABLE]

III.3 Results

The results of Eq. 9, 10, 11 in 2D for the distances corresponding to Si:P and Si:As in Table. 1 are plotted in Fig. 5. The plotted quantity is the density of readouts which are actively participating in ‘successful’ entangling gate configurations, ie. twice the density of entangling gates in the SFG case and the same density as that of entangling gates for the Heisenberg single readout gate cases. The percentage of active readouts is simply the density of readouts which are a part of an entangling gate divided by the total density of readouts.

In all cases, the total density of controls is maximised by choosing $D_{c}=\frac{1}{V(R_{cc})}$ from Eq. 5, yielding the viable control density $\frac{1}{e\times\>V(R_{cc})}$ from Eq. 6, as can be seen on Fig. 3. The optimal control densities calculated with the distances of Table. 1 can be seen in Table. 3.

In 2D, the maximum density of readout dopants (Si:As) which are a part of ‘successful’ SFG gates is 6 $\times 10^{8}$ dopants per cm2, which corresponds to a total readout doping density of 1.5 $\times 10^{11}$ dopants per cm2. The bilayer case shows a negligible increase in the density of readouts which are part of SFG gates, for double the total density (3 $\times 10^{11}$ dopants per cm2). The percentages of readouts involved in ‘successful’ SFG configurations in both the monolayer and bilayer cases remain below 0.5%.

The maximum density of readout dopants (Si:As) which are a part of ‘successful’ excited-readout Heisenberg gates is 1.2 $\times$ 109 dopants per cm2, corresponding to a total density of 8 $\times$ 1010 dopants per cm2 and 1.5% active readouts. A very clear increase can be seen in the bilayer case. The maximum density of readout dopants part of ‘successful’ excited-readout Heisenberg gates in bilayers separated by 13.2nm is 3.5 $\times$ 109 dopants per cm2, corresponding to 2 $\times$ 1011 dopants per cm2 and 1.8% active readouts.

The control (Si:P) doping density that achieves maximum density of active readouts is 1.79 $\times$ 1010 dopants per cm2 in the monolayer case and 3.91 $\times$ 1010 dopants per cm2 in the bilayer case. It is possible to achieve up to 3.9% active readouts in the monolayer case and 4.7% active readouts in the bilayer case for low doping densities. Thus, using equal total densities of Si:P and Si:As leads to some of the highest active percentages of Si:As contributing to excited-ground Heisenberg gates.

The 3D bulk doped results can be seen in Fig. 6. The maximum density of readout dopants (Si:As) active in excited-ground Heisenberg type gates is 3 $\times$ 1014 dopants per cm3 corresponding to a total doping density of readouts of 3.5 $\times$ 1016 dopants per cm3 and 0.8% active readouts. Similar to the 2D case, the maximum active readout percentages of 2.1% can be reached for total readout doping densities comparable to that of the controls (3.2 $\times$ 1015 dopants per cm3). The maximum density of readout dopants (Si:As) active in SFG type gates is 2 $\times$ 1014 dopants per cm3 corresponding to a total doping density of readouts of 7 $\times$ 1016 dopants per cm3. It was possible to gain insights about the densities of gates containing three or more readouts from the Monte-Carlo simulation, which may be of interest for different types of quantum gate, e.g. Toffoli gates. Unsurprisingly, they peak at far higher total readout doping densities, such as 7.5 $\times$ 1016 dopants per cm3 for the three readout case (see Fig. 6), but provide lower active readout densities, such as 7 $\times$ 1013 dopants per cm3 for the same case.

Finally, the single species excited-excited Heisenberg gate in 2D yields the highest densities and percentages of active dopants. Active dopant densities of 5.4 $\times$ 109 dopants per cm2 can be reached for total doping densities of 2.9 $\times$ 1010 dopants per cm2, correponding to 20% of dopants being involved in ‘successful’ excited-excited Heisenberg gates! The active percentage maximum is at 27%, which corresponds to an active density of 3.8 $\times$ 109 dopants per cm2 and a total doping density of 1.4 $\times$ 1010 dopants per cm2.

IV Proof of principle experiment: Control of magnetization dynamics by orbital excitation

As a first step towards the implementation of the optically excited exchange entangling gates we propose a proof of principle experiment to show the control of the magnetization dynamics of the readout (Si:As) donor electrons by the control (Si:P) electron’s orbital state, due to enhancement of exchange interactions from orbital excitations. This is the most important building block of the SFG entangling gate. We calculate the quantum many-body magnetization dynamics in the thermodynamic limit using the diagonalization method Moving Average Cluster Expansion (MACE) Hazzard et al. (2014). Because deep donors in silicon act according to the Heisenberg Hamiltonian, this experiment can also be interpreted as a quantum simulation of the two-species S=1/2 Heisenberg anti-ferromagnet111This has been shown to be an effective low energy description of the half-filled Hubbard model, in the case where the interaction energy between the spins of the lattice is far greater than the hopping strength Auerbach (2012); Kubo and Uchinami (1975); Polatsek and Becker (1996). with quenched disorder.

IV.1 Experimental proposal

Starting from a spin polarized state we show below that when all donors are in the orbital ground state, the magnetization of Si:As stays constant on observable timescales, while it changes to a vastly different value when the control species is excited to the 2p state. This implies that while exchange interactions are negligible in the orbital ground state, they are large enough to exhibit non-trivial many body dynamics when in the excited state. Realizing this in experiment would at least partially prove the necessary control over two-body interactions needed for the quantum logic gates discussed in the previous sections.

In experiment, readout of the magnetisation of the Si:As is done via donor-bound exciton spectroscopy (D0X spectroscopy)Saeedi et al. (2015); Morley (2015). This requires electrical detection, which can be achieved in a dilute 2D layer of P and As donors using STM hydrogen lithography to pattern highly conductive metallic-doped phosphorous pads into the same plane and overgrowing a protective thin-film of crystalline silicon with Molecular beam epitaxy. The impurity sheet’s metallically doped pads are electrically contacted using electron beam lithography coupled to reactive ion etching to create features which are filled with aluminium using a metal evaporator. The metallic pads are in turn electrically contacted to obtain, in conjunction with terahertz radiation from a Free Electron Laser (FEL) ($\lambda=31.6\mu$m for the Si:P 1s to 2p+- transition) and D0X spectroscopy, an electrical signal from the 2D Si:As sheet which is a response to the coherent and non-linear excitations of the Si:P electrons. The sample fabrication and electrical detection technique briefly described above and which we have in mind for the experiment we propose here are described in Crane et al. (2018). This detection technique enables the precision condensed matter samples to remain intact after exposure to a FEL pulse.

The contacted impurity layer in silicon is mounted on the bore of a water-cooled Bitter magnet to Zeeman split the ground state impurity electron spin energies, leading to six pairs of dipole-allowed transitions ($\Delta$m = 0, $\pm$1). Electrical detection of electron spin resonance using D0X spectroscopy has been demonstrated for magnetic fields of around 0.35 Tesla Lo et al. (2015). A donor-bound exciton can be formed by a direct 1.15 eV photon (the silicon indirect bandgap is of 1.17 eV). The photon excites an electron from the valence band, leaving behind a hole. When the electron-hole pair recombine via an Auger recombination process, their energy ejects another electron from the donor site, leaving behind a positively charged donor ionPortis et al. (1953). To relax all the electron spin states of both donor species to the lower energy spin state, the sample needs to be cooled down in a dilution fridge to milli-Kelvin temperatures. D0X optical pumping then initialises all the readouts to the opposite spin state from the controls (and later to read out the occupation of one of the Si:As spin states). Each electron has the initial spin state $\ket{S_{i}}$ with

[TABLE]

and the whole system is initially in a product state $\ket{\psi_{0}}$ such that:

[TABLE]

When the terahertz frequency FEL is switched on, it illuminates the entire sample and excites the control (Si:P) electrons to the 2p+- state. The experimentally demonstrated decay time of 3D bulk doped Si:P in the 2p+- state in non-isotopic silicon is 200ps Litvinenko et al. (2015), which gives the timescale available to our experiment.

The dopants evolve in time according to the Heisenberg Hamiltonian

[TABLE]

where we made the species dependence of the interactions explicit. The strength of the interaction $J^{PAs},J^{PP}$ depends explicitly on the orbital state of the control species (Si:P) (see Appendix A). If the densities are high enough, the Si:P wavefunctions will overlap with the Si:As and lead to non-trivial dynamics.

The average magnetization of the readouts (Si:As) is given by

[TABLE]

The average spin flip probability is then defined as

[TABLE]

where $\ket{\uparrow}\bra{\uparrow}_{i}$ is the projector onto the $i$-th spin being in spin up.

To make a differential measurement, we compare the average spin flip probability of the readout dopants when the control dopants are in the ground state with the average spin flip probability when the control dopants are collectively excited to the 2p+- state. As a result, we then expect

[TABLE]

We can deduce that the non-trivial dynamics of the spin flip probability must have been due to interactions between species as the species themselves are in an eigenstate of their respective Heisenberg Hamiltonian.

In the following, we will now predict the outcome of the experiment using an exact diagonalization technique.

IV.2 Dynamics of the arsenic spin flip probability within MACE

In order to simulate the experimental outcome, i.e. the magnetization dynamics, we use the Moving Average Cluster Expansion(MACE) technique which has been shown to capture the magnetization dynamics in disordered long-range interacting quantum spin systems realized by cold dipolar molecules Hazzard et al. (2014) and Rydberg atoms Piñeiro Orioli et al. (2018).

MACE assumes that in a system where spins are randomly placed in space, contributions to the local magnetization dynamics of a particular spin are only made by the spins with which it has the largest interaction exchange energies. Therefore, it is enough to exactly diagonalize this cluster and then average over all such clusters. Convergence is then sought by increasing the cluster size. In our system, convergence was found for cluster sizes 5 and above due to the short-range nature of the interactions. This is in contrast to the large cluster sizes needed in algebraically decaying interactions Schuckert et al. (2018). Plotted on Fig. 8 are results from cluster size 8 with jackknife error estimates in grey Ambegaokar and Troyer (2010).

Neglecting the effect of the magnetic field on the strength of the exchange interaction, we calculated that there is no visible change in the magnetisation of the arsenic for fields up to 10 Tesla. We have also checked that there are no significant differences between defining clusters with the largest exchange interactions or the nearest neighbours of the spin of interest.

The difference in the dynamics of the Si:As spin-flip probability between Si:P being in the orbital ground state or the excited state can be seen in Fig. 8 to be 13% within the 200ps experimentally detected decoherence time of the Si:P 2p+- excited state, for doping densities of both As and P of $7\times 10^{9}$ dopants per cm2, which is accessible with current doping techniques.

This shows that already with only global measurements at hand, non-trivial quantum many-body dynamics could be probed within this dual-doped system. Comparing measurements where Si:P is in the ground and excited state would furthermore show that substantial control over the dynamics of the Si:As spins can be implemented within this scheme, therefore showing that one of the prime requirements of the SFG entangling gate are within reach of current experimental platforms.

Conclusion

In this paper we have established the optimum densities of optically controlled entangling gates achievable in randomly doped silicon. In order to determine which spacing between dopant electrons can realistically create entangling gates, we created spatial maps of the Heisenberg interaction exchange energy between same-species and dual species dopants. We focused on Si:P and Si:As because they are, to date, the most well-understood dopants for diffusion onto a silicon surface in ultra-high vacuum.

Respecting these conditions we obtained matching results in both the nearest neighbour treatment of a Poisson point process and a Monte-Carlo simulation. Densities of entangling gates in 2D were in the low $10^{9}$ dopants per cm2 and in 3D were in the low $10^{13}$ dopants per cm3, corresponding to relatively low percentages of active dopants. However, by considering same-species ‘Heisenberg’ gates, similar densities corresponded to percentages up to 27% of dopants involved in ‘successful’ gate structures. We also showed that in 3D there arise situations in which gates with three and four qubits reach significant densities. Finally that by dividing the donors into two parallel 2D layers the densities of gates increase still further.

Finally, we proposed a proof of principle experiment aimed at demonstrating the onset of two-body (entangling) interactions caused by the optical excitation of one of the species in a 2D randomly doped sample. The magnetisation dynamics of the Si:As donor electron spins were calculated using the Moving Average Cluster Expansion technique, depending on the orbital state of the Si:P donor electrons. For low but equal densities of Si:P and Si:As, we found a large difference in the time evolution of the spin flip probability of more than 10% between the cases where phosphorous was in the ground or excited state.

This work can be extended in various ways. Single donor placement techniques such as hydrogen lithography could be used to create two or three dopant structures in silicon to verify the maps of the exchange interaction energy calculated here. The theory we have developed can be applied to different dopant species such as selenium or acceptors such as boron, could be modified to accommodate hard-sphere configurations and which would be relevant for Rydberg atoms which deviate from pure Poissonian statistics in the blockaded regime Saffman et al. (2010).

Acknowledgements.

We thankfully acknowledge discussions with K. Stockbridge, B. N. Murdin, N. J. Curson, G. Matmon and R. Crane. A.S. thanks the London Centre for Nanotechnology for hospitality. We gratefully acknowledge the financial support from the UK Engineering and Physical Sciences Research Council [COMPASSS/ADDRFSS, Grant No. EP/M009564/1]. A.S. acknowledges financial support from the International Max Planck Research School for Quantum Science and Technology (IMPRS-QST) funded by the Max Planck Gesellschaft(MPG).

Appendix A Exchange calculation with multivalley wavefunctions

In order to calculate the exchange interaction between two donors we use the Heitler-London approximation Koiller et al. (2001)

[TABLE]

where $\bm{R}=\bm{R}_{2}-\bm{R}_{1}$ is the separation vector between the two donors and $\epsilon_{r}=11.4$ is the dielectric constant of silicon. The wavefuction of each donor can be either the $1s\text{A}_{1}$ ground state, the $2p_{0}$ or $2p_{\pm}$ excited states. In the multivalley effective mass theory all of these wavefunctions can be expanded as

[TABLE]

where $j=1,2,3$ indicates the $1s\text{A}_{1}$ ground state, the $2p_{0}$ and $2p_{\pm}$ excited states, respectively; $\mu=\pm x,\pm y,\pm z$ indicates the valleys of silicon’s conduction band minima, $\phi_{\mu}(\bm{r})=e^{i\bm{k}_{\mu}.\bm{r}}u_{\mu}(\bm{r})$ are the Bloch functions at the minima, $|\bm{k}_{\mu}|=0.84\times 2\pi/a_{0}$ where $a_{0}=0.543$nm is the lattice constant. The $F_{\mu}(\bm{r})$ are the envelope functions. For the $1s\text{A}_{1}$ state the envelope function of the $+z$ valley is Kohn and Luttinger (1955)

[TABLE]

where $a_{1}=2.42$nm and $b_{1}=1.39$ nm. The factor $\alpha$ accounts for the contraction of the ground state Thomas et al. (1981) due to the central cell correction (CCC) and is given by $\alpha=\sqrt{E_{\text{SV}}/E_{B}}$ where $E_{\text{SV}}=29.7$ meV is the binding energy obtained from a single valley theory without the CCC, as done in Ref. Kohn and Luttinger (1955), and $E_{B}$ is the experimental binding energy which is $45.58$meV for Si:P and $53.77$meV for Si:As Ramdas and Rodriguez (1981). The other envelope functions are obtained by using $F_{\mu}=F_{-\mu}$ and cyclic permutations of $x,y,z$.

The excited state energies and wavefunctions of Si:P are identical to those of Si:As Ramdas and Rodriguez (1981); Li et al. (2018), and is dependent on the polarization of the light field Chick et al. (2017). For polarization with the unit vector $\bm{\epsilon}=[\epsilon_{x},\epsilon_{y},\epsilon_{z}]$ the excited states are also given by Eq. (20) but the envelope functions are now

[TABLE]

where $a_{2}=3.68$nm, $b_{2}=2.23$nm, $a_{3}=5.45$nm, and $b_{3}=3.35$nm. The other envelope functions can again be derived using $F_{\mu}=F_{-\mu}$ and cyclic permutations of $x,y,z$.

We can further expand the Bloch function in Eq. (20) in terms of plane waves as $\phi_{\mu}(\bm{r})=\sum_{\bm{G}}c_{\bm{G}}e^{i\bm{G}.\bm{r}}$ where $\bm{G}$ is the reciprocal lattice vector. Substituting this into Eq. (20) and then Eq. (A), neglecting the fast oscillating terms in the resulting integrand, and using the equality $\sum_{\bm{G}}|c_{\bm{G}}|^{2}=1$ we arrive at

[TABLE]

where

[TABLE]

These highly oscillating integrals are evaluated numerically with the Vegas Package Hahn (2005).

Appendix B Monte Carlo Simulation

The Monte-Carlo simulation was written in php with a MySQLi database. The front-end is written in html/css so that simulations can be run through a browser and made available on the internet. It scales optimally with the number of dopants, i.e. O(nc + nr), where nc is the number of control dopants and nr is the number of readout dopants, such that it allows for computations on large densities and such that the wait time is minimal. Optimal scaling stems principally from the avoidance of full Pythagoras computations. This is achieved with four different techniques. As regards the first step in the whole algorithm which is to analyse the controls in order to identify all viable controls, the avoidance of Pythagoras is achieved by first partitioning the space according to the expected mean nearest neighbour distance, second by using Pythagoras only for points lying between the enclosing squares (Fig. 9), third by skipping to the next point as soon as conditions are found to be breached and fourth by reciprocating information between adjacent points.

As regards the second step which is to work within the viable control distribution and the full readout distribution, avoidance of Pythagoras is achieved by using the search sequence represented diagrammatically in Fig. 10, referred to as a Haystack search.

Only the readouts close enough to a viable control and which might therefore constitute or interfere with an entangling gate are analysed. Consequently slightly more than half the readout dopants are never involved in the algorithm.