Using non-homogeneous point process statistics to find multi-species event clusters in an implanted semiconductor

The Poisson distribution of event-to- i th -nearest-event radial distances is well known for homogeneous processes that do not depend on location or time. Here we investigate the case of a non-homogeneous point process where the event probability (and hence the neighbour configuration) depends on location within the event space. The particular non-homogeneous scenario of interest to us is ion implantation into a semiconductor for the purposes of studying interactions between the implanted impurities. We calculate the probability of a simple cluster based on nearest neighbour distances, and specialise to a particular two-species cluster of interest for qubit gates. We show that if the two species are implanted at different depths there is a maximum in the cluster probability and an optimum density profile.


Introduction
Individual interacting impurity atoms can be important for donor qubit gates, such as that proposed by Stoneham et al [1], while an important class of theoretical physics problems is produced by the Hubbard model, which relies on hopping and magnetic interactions between neighbours in chains [2]. In the case of donor impurities in a semiconductor, deterministic placement using scanning probe tips has improved greatly in recent years, but is currently limited to a small number of species of impurity (principally phosphorus and arsenic [3] in silicon [4,5] and germanium [6], and Mn in GaAs [7]). Ion implantation methods can also be used to create impurity layers in semiconductors with merits including flexibility with regards to the numerous available implantable species and far faster device fabrication times which are less costly and more easily scalable. These merits clearly come at the cost of much less precision. Given the stochasticity of the donor placement it is important to look at the effects of the implant distribution on the neighbour-neighbour distances, and hence the probability of observable interactions. Contemporary work in this area [8] has focused on analytically understanding the interactions between donors, the dependence these interactions have on donor spacing and using the results of homogeneous Poisson point process statistics, optimising for these interactions. Here we generalise the statistics to include inhomogeneity of the impurities and optimise the event density profiles for a comparable event cluster definition. We show that in the case of a Gaussian distribution of events (which is a good approximation of the distribution of impurities after ion implantation) an analytic solution for the nonhomogeneous nearest neighbour distribution exists. We also show that the numerical optimisation involved can be accelerated by introducing an appropriate heuristic.
Many physical problems involving stochastic probability have been studied which make use of point process statistics. They have been used to model distributions of events ranging from plants in a field [9], the locations of cellular network base stations [10] to the distribution of astronomical bodies [11]. The ability to model a distribution of points in an event space and be able to quantify irregularities such as clustering of points helps provide an insight into correlations between events and the consequences resulting from such a distribution. Using the well understood construct that is the Poisson point process [12,13], the nearest neighbour distribution of an impurity species has been used to model optical properties of donors in silicon [14] due to nearest neighbour interactions and also to calculate the probability of finding large clusters of donors [15] in homogeneously doped bulk semiconductors.
It is useful to discuss clusters comprising two different species so that we may expect to detect the effect that excitation of one of the species has on its interaction with the other. Species specific control/detection might make use of optical or electronic resonances. It seems reasonable to expect that if density profiles of species A and B are implanted at a different depth then we can control the distributions for the  A B separation, and by controlling the peak densities of each profile we might control the  A A and  B Bseparations. It is clear that if the homospecies separations (i.e.  A A or  B B) are small compared with the width of the implant profile then the likelihood of an A−B interaction will decrease relative to homospecies interactions and signal will be lost. This suggests low density sheets are ideal for inter-species interactions. Conversely, if the densities are very low then the  A B separation will be controlled by the density rather than the separation of the sheets. There is clearly an optimum to be found. The statistics developed in this work, though applied to the problem of interacting donor impurities for quantum technology, has been presented so that readers interested in structure between multiple species of events described by non-homogeneous density distributions may easily apply the ideas contained to their work.
In this work we analyse the distribution of nearest neighbour (NN) distances existing in a two-species point process where the density of those species are non-homogeneous in depth. When the event density varies as a Gaussian profile in one spatial dimension a solution to the distribution of ithnearest neighbours is shown to be analytically solvable. We have extended the calculation to the probability of finding particular multi-species cluster configurations and show that optimising the density profiles in favour of such clusters is a numerically feasible task which may be improved by an approximation. The optimised probabilities found with the aid of such an approximation were checked against those using the analytical solution which was also supported by simulation of the deemed optimal doping parameters.

Non-homogeneous Poisson point process
In this section we give some definitions of symbols useful later, and show the relationship between nonhomogeneous and homogeneous Poisson distributions in order to indicate the method for investigating neighbour-distances.
If the distribution of events is non-homogeneous, the density n x ( ) and the expected number of events d d = N n x x ( ) in the infinitesimal volume dx varies with the location, x. If δN is small then the probability of an event in dx is equal to δN . We may now find the probability that there are no events within a larger volume V by dividing it up into elemental volumes. The expected number of events in V is Assuming n x ( ) is well-behaved, we may choose the size of the ithelement dx i ( )so that the product d n x x i i ( ) is a constant. The probability of an event within dx i is then the same for every element, and it follows (see appendix A) that the probability of m events enclosed in the larger volume V is given by the probability mass function for a Poisson distribution It is tempting to try to use the probability d n x x ( ) (of an event in the elemental volume dx around x) to define a probability density function (PDF). This is best avoided because the probability of an event in a larger volume is not given by the integral over that volume ò }from comparison of equations (1) and (2) (unless V and N are very small, obviously). This is due to the fact that the possibility of two or more events is nonnegligible for a large volume. Later, we shall be concerned with both questions of counting types of events in some volume where an integration like equation (1) is needed, and of the probability of occupation by specific numbers of events (0, 1, 2 or more events etc) in some volume where equation (2) is needed.
Previously, clusters of impurities with homogeneous density have been discussed in terms of the distribution of neighbour-neighbour distance [15]. In order to put our discussion into this context we give the nonhomogeneous case, which follows immediately from equation (2). We define the probability d )has its ithnearest event of species A at a radial distance between d  + r r r. We shall refer to  p r A x i ( ) as the nearest neighbour probability density function (NNPDF). In previous literature the NNPDF is the precursor to what is referred to as the 'void nearest neighbour distribution function' [12] which is simply the cumulative distribution of the NNPDF as defined here. The term 'void' is used since there is no event specified at the point x whose neighbour is being found. To calculate the NNPDF, consider the sphere V r x ; sphere ( )centred on x of radius r, and the infinitesimal shell of thickness dr around it. The probability of finding the first nearest event within the shell is equal to the product of the probabilities of having no events within the sphere and one within the shell (since these two conditions are independent): where the probability of an event in the shell is . Hence the NNPDF for the ithnearest neighbour is The distribution around a void can be extended to the distribution for the neighbours around an event by taking into account the density of events in the infinitesimal volume dx at x. We can recover the homogeneous results in 3D bulk doped or perfect 2D delta-doped layers in semiconductors, etc. For example, in 3D we replace . Hence all terms in equations (2) and (3) are independent of x and we obtain the familiar 3D homogeneous neighbour-neighbour distributions,

Poisson point process with expectation varying in 1D
In this work we are particularly concerned with impurities that have been implanted from the surface. We therefore specialize to the case with inhomogeneity of event density in only one dimension, specifically when the density n A is a function of z only. This is the case with broad-area ion implants which produce a finite spread of penetration perpendicular to the surface and are homogeneous in the plane. The expected number of events per unit area in To find the expected number of events in the sphere of radius r, one integrates over the thin discs perpendicular to z (which have constant density). For a sphere V r z where the limits of the integral are such that the area in the square brackets is positive, i.e. from z−r to z+r. Differentiating under the integral sign, for use in equation (3).

Neighbour separations with Gaussian density profiles
In the case of Gaussian density profiles which are reasonable approximations for typical mono-energetic implants, is the r.m.s. thickness of the density profile. When substituted into (3) for i=1 this gives the analytical solution for the probability of a first nearest event within d  + r r rmeasured from a starting point at depth z and hence the NNPDF: Having found the NNPDF we can further give the number of A atoms (per unit area) between since the density of Aʼs at z and the probability density for the surrounding Aʼs are independent. This is easily generalized to a multiple species situation: the density of Aʼs (per unit area) at z with an ithnearest species B neighbour at r is =   n r z p r n z , . 8 whereby the two Gaussian density profiles for n z A B , ( ) have a unit separation in depth, and each has an r.m.s width of 1 2 i.e.
unit of length, and unit height i.e.
. Referring to figure 1 this is equivalent to integrating over horizontal slices or projecting onto the r-axis. The probability of, say, the first nearest neighbour B having a separation of 1 unit from an A can be optimised by varying the separation and density of the layers. In the next section we consider the effect of adding constraints on the next nearest neighbours.

Density of specific cluster configurations
The total number of useful clusters, N good , is a question of counting events and may be found from an integration like equation (1). The number of useful clusters in the elemental volume dx around x is given by the number of Aʼs in the elemental volume multiplied by the probability that each is part of a useful cluster, Hence the total is . The centre of mass of the NNDF for i=1 has a z position close to the peak of the A density and an r position defined by the depth separation between the A and B layers. As expected, the 20th nearest neighbours are further away than the 1st nearest. Vertical slices may be used to find the expected distance to a B for Aʼs of specific depth, and integrating horizontal slices (projecting onto the r axis) produces the total probability density irrespective of the depth of A.
We have used the shorthand notation P Good Cluster Ax { } for the conditional probability of a useful cluster configuration around x when it is given that a central A exists at point x. This probability depends on the probability of specified numbers of events in regions of space around the A, and must be found from computations like equation (2). In the case of a 1D non-homogeneous density variation along z, the number of useful clusters (per unit area) is

Cluster probability for donor qubit gates
Here we investigate the calculation of P Good Cluster Az { } for the example case of a simple qubit gate made with donors in silicon.
We are interested in multi-species clusters that have specifications on the separations. We imagine a pair of qubits that carry quantum information in their spin. The gate operation is performed by controlling the entanglement between two impurity electrons. By changing the state of the impurity of species A which we name the control, it will get entangled with a nearby impurity electron of species B which we name the target. This controlled change of state might be by optical excitation or use of electric fields etc. To facilitate the controlled interaction they must be separated by an appropriate radial distance, equivalent to defining a radial interval rather than an infinitesimal shell. We imagine such changes in interaction range occur in other fields such as in ecological networks of consumer species and resource species where there is e.g. a seasonal change in the interaction strength. We now add a specification on the next-nearest neighbours because the control and target impurities should be sufficiently isolated from the environment, i.e. other impurities, that they do not decohere. In ecology this might be analogous to the effects of competition.
With the results of the previous section we can calculate the probability density for an A donor having its first nearest (target) B at the optimal distance using equations (8) and (6) (or (3) for non Gaussian event density profiles). We now examine the combination of this condition with specification that the second nearest B and first nearest A are out of range. We allow for some tolerance on the useful target B distance.
In the simplest specification of our useful cluster for which the probability is P Good Cluster Ax { } , we define a useful cluster as one in which the A control atom has: 1. its nearest A outside the range >  r r ii. region '3' on figure 2 includes exactly one B. This region will be referred to as V Ax target , the volume around the x ( ) (other than those at x and ¢ x , which we consider to be elemental volumes so small that it does not affect the required expectation). Since a useful cluster having its B at ¢ x is mutually exclusive with a useful cluster having its B at  x we can add these probabilities, i.e. integrate over the allowed range of ¢ x . In the case of the 1D non-homogeneous problem the probability of a B in an elemental ring at cylindrical coordinates ¢ z and ¢ r c (from the vertical axis containing the central A control atom) is p d d ¢ ¢ ¢ ¢ r n z r z ( ( ) )is a slice through the intersection of two spheres, which is surprisingly complicated but may be written analytically. Even so, the integral in equation (11) is a nested triple integral with complicated bounds. In cases where many calculations of P Good cluster z { } are required, such as in our problem of optimising the species density profiles, it is helpful to produce a heuristic method that accelerates the numerical calculation of this probability.
Heuristic method to approximate the best case cluster probability So long as it is given that there is only B 1 within the region V Az target , then the probability of finding B 1 between ave ( )(as usual the subscript indicates it is given that there is an A control atom at z). We may now use an approximate version of equation (11): There are a number of reasonable but different choices for calculating the most important location of the target ¢ Z Az and ¢ R Az for use in equation (13). We tried finding the expectation radius using (   where ¢ C z ( ) is the circumference of the small circle of the sphere ¢ R Az through ¢ z . We also tried finding the expectation depth and then the expectation spherical radius given this depth. Finally, we tried looking for the most likely depth (where the maximum n B (z) occurs) and most likely radius (at the outer edge of V Az target ). The latter is the easiest to find and requires no integration yet dramatically underestimates P Good cluster Az { } . We found that calculating ¢ R Az first gave the best agreement with (11) for the good cluster configuration we were interested in.
An approximate solution to this probability which is less computationally intensive accelerates the process of numerically optimising that probability. The closer the approximate solution is to the optimum found using the vigorous method, the more efficiently one can converge to an optimum Gaussian doping profile.

Results of optimising cluster probability
We provide a numerical example of the cluster optimization in the case of the silicon donor qubit gate using separation tolerances estimated from consideration of exchange interactions [8]. The separation range for the control to target distance is from =  r 15nm , and the separation of the two layers μ. Here we make the assumption that both profiles can be implanted at different depths with the same width. In practise, independent control of layer depth and width is not achievable with ion implantation and the depth profile for a particular implant species and target depends principally on the implant energy. The profile of impurities will also change during necessary post-processing such as diffusion during annealing and whether each species is annealed as implanted or if the system is annealed only once after all implants. As these specifics are dependent on the species and target of interest we leave possible simplifications to the reader. It has been shown previously that Bi implanted at a high energy into silicon (therefore creating many lattice defects) can be electrically activated to a high enough quality that the donor electron spin states are measurable [16]. This demonstrates the feasibility of ion implantation as the means creating the qubits of interest here. The following optimisation can be used to determine a best case implant profile if the final donor profiles are expected to be Gaussian. If the final active donor profile can be measured, it may be used directly with equations (10) and (11) to determine the final viable cluster yield.
We first maximised the total density of good clusters, N D good 2 for a given combination of width and separation using the heuristic as detailed in the previous section. This is shown in figure 3(a) For sufficiently large layer widths layer separation clearly has no effect on either N D tend to the optimum homogeneous bulk densities [8]. We see that for layers spaced far apart the density of viable clusters tends to zero as the layer width is reduced, as expected as it becomes increasingly difficult to obtain an  A B distance within the allowed range. For very narrow layers, there is an obvious optimum layer separation of 15 nm which can be deduced from geometrical configurations relating to the definition of our 'good configuration'. This distance is the same as the minimum separation of control from target. For intermediate layer widths of around 10 nm the lines cross, and the optimum is now obtained for layers of zero separation, i.e. the target and control layers should be at the same depths.
The areal density N D good 2 is not necessarily the most useful optimization objective function. Even with a small density, the total number of good clusters may be increased simply by increasing the sample size. The ratio of signal to background in an experiment might be improved if instead we maximise the fraction of donors that are involved in a good cluster. For example, we imagine an experiment where we detect the effects of the interaction by measuring the effect on the spin of the Aʼs after exciting the Aʼs (which produces an effect only for those Aʼs that are part of a good cluster). In our simple cluster of interest there is only one A per cluster, so the number of Aʼs involved in a good cluster is just equal to the number of good clusters, and the signal to background will be optimized by maximising the fraction N n D A D good describing the ratio of signal to background shown in figure 3(b). 3(b) shows that if one were able to fabricate atomically flat layers separated in depth then this fraction is optimised when such layers are separated by 15 nm. An interesting situation arises if there is a lower limit on the possible width of the density profiles. This is the case when ion implanting species into a lattice. Depending on the implantation specifics, it is difficult to make very thin layers due to ion straggle. Figure 3 shows that if the layers cannot be fabricated with widths less than 10 nm, the optimum configuration is to have the two species coplanar (m = 0) with widths of »15 nm (not to be as thin as possible) achieving a reasonable 10%. The fractions in both of these cases are an improvement over the optimised, bulk dope case achieving a good cluster fraction of 9%. The results of 3 doping configurations were simulated using a brute force (Monte Carlo) approach (filled squares with error) confirming the densities as calculated by (10).
The quality of the optimization is shown in figure 4, in which we compare it with optimization using the full solution of equations 11-12, and also using a brute force (Monte Carlo) method. Agreement is excellent for cases examined here where the two species profiles overlap in depth. We found that the heuristic method does not agree with the full solution when the layers are thinner than the separation between layers as shown in figure 5. Under these circumstances the optimum B density can differ by as much as 2 orders of magnitude.  . Focussing on the co-planar (μ=0) configuration we show how close our heuristic method agrees with the full result for the density of viable clusters in that configuration. Both methods are shown to be in agreement with a brute force approach (filled squares with error). Here the error in the mean fraction of A impurities that gate viable clusters is more visible. This error was minimised through repeated simulation. Like the full solution, optimising using the brute force approach is considerably more impractical than the solution found using the suitable heuristic.

Conclusion
This work has demonstrated how the ideas used in non-homogeneous Poison point process statistics may be used not only to describe the distribution of ithnearest neighbour distances between events that exist nonhomogeneously in space but also to calculate the probabilities of events forming more complex structures existing in such a process. When the event species density inhomogeneity is described by a Gaussian probability density function of depth and with complete spatial randomness in the xy-plane an analytic solution to the NNPDF is yielded. This density surface gives a good representation of event-event separations existing in nonhomogeneous processes. for the purposes of research into solid state qubit gates, we calculate the frequency of appearance of suitable structures for qubit interactions in random doped samples. The introduction of a heuristic for calculating the probabilities of such structures improves the ability to optimise the dopant profiles for these interactions. Such an optimisation will influence the fabrication requirements of test-bed solid state quantum devices whereby distance dependant interactions with nearest neighbours are of paramount importance.
The sizes of the i th elements dx i ( )that make up the finite domain V in which there are an expected N(V ) events can be varied such that the expected number of events within those elements is constant.
where M is the number of elemental volumes. As N is finite, we can make M large enough and dN small enough that we may neglect the possibility of two or more events in any element; hence the probability of no events is d -N 1 for each element. Since the incidence of events in each element is independent from the others, the probability of no events in V is ). This is the same as the usual Poisson probability of zero events when the expectation is N, in spite of the fact that n x ( ) is non-homogeneously distributed within V. We may further find the probability of one event in the volume V by considering one additional elemental volume added to it; ). Continuing similar arguments for higher numbers of events, m, we again recover the standard Poisson distribution (2) with expectation N (even though n x ( ) is non-homogeneous).