Abstract
In the midst of the epitaxial circuitry revolution in silicon technology, we look ahead to the next paradigm shift: effective use of the third dimension  in particular, its combination with epitaxial technology. We perform ab initio calculations of atomically thin epitaxial bilayers in silicon, investigating the fundamental electronic properties of monolayer pairs. Quantitative band splittings and the electronic density are presented, along with effects of the layers’ relative alignment and comments on disordered systems, and for the first time, the effective electronic widths of such device components are calculated.
Keywords:
Ab initio; Density functional theory; Bilayers; Phosphorus in siliconBackground
We are currently living through a transition in electronic circuitry from the classical to the quantum domain. With Moore’s Law on the way out, thanks to the recent unveiling of ohmic 2 nm epitaxial nanowires [1] and epitaxially gated singleatom quantum transistors [2], the challenge for scientists becomes finding new ways to increase the density and speed of devices as we can no longer rely on being able to shrink their components.
Farsighted speculation has already been abundant for many years regarding efficient use of the third dimension in device architecture [36], breaking the twodimensional paradigm of current electronics manufacturing techniques. Recent germaniumbased works [7,8] illustrated fundamental physics required for full 3D device implementation and heralded the creation of multiple stacked δlayers of dopants within a semiconductor. Each of these layers could potentially display atomically abrupt doped regions for inplane device function and control. Multiple layers of this nature have indeed been created in Ge [9]. The P in Ge atomic layer deposition technique parallels phosphorus in silicon 1/4 monolayer (ML) doping (Si: δP), created using scanning tunnelling microscope lithography, with a few minor technological improvements (annealling temperatures, amongst others) [8]. In contrast, one major advantage of improvements to silicon technology is that uptake may be far easier, given the ubiquity of silicon architecture in the present everyday life. We may therefore expect to see, in the near future, Si: δP systems of similar construction.
The time is thus ripe to attend to possible threedimensional architectures built from phosphorus in silicon. Although Si:P singledonor physics is well understood, and several studies have been completed on singlestructure epitaxial Si: δP circuit components (such as infinite single monolayers [1017], single thicker layers [18,19], epitaxial dots [20], and nanowires [1,21]), a true extension studying interactions between device building blocks in the third dimension is currently missing.
The description of experimental devices is a thorny problem involving the tradeoff between describing quantum systems with enough rigour and yet taking sufficient account of the disorder inherent to manufacturing processes. A first approach might therefore be to study the simplest case of interacting device components, namely two Pdoped single monolayers (bilayers) [22,23]. Given the computational limitations of ab initio modelling it is currently not possible to treat the disordered multilayer system in full. Two approaches suggest themselves. In [23] the approach was to simplify the description of the deltalayer in order to study disorder through a mixed atom pseudopotentials approach. Here, we instead develop a rigorous model of an idealised, perfectly ordered multilayer system in order to make connections to an understanding gleaned from both the mixedatom approach and from other idealised models. The two approaches are complementary: alone, neither achieves a complete description, but together, they offer good comparisons from which one may draw the firmest conclusions available regarding experimental devices. The second approach, dwelt upon in this work, also offers descriptions of systems that should become available with improvements to the manufacturing processes mentioned above. As such, this is the focus of our discussion.
Whilst singlemonolayer studies converge properties by increasingly isolating the layers [11,14,16], at closer separations, it is impossible to divorce specific interactions between two layers from those between all of their (infinite) periodic replications. Further, effects arising due to atomicscale mismatches in each layer’s doping locations cannot be seen when the neighbouring layer is a perfect replica. Building upon the methodology established whilst investigating single δ layers [16], expanded upon when considering thicker layers comprised of multiple adjacent δ layers [19], and further extended to consider δdoped nanowires [21], here, we model Si: δP bilayers, varying both their vertical separation (Figure 1a) and their relative inplane alignment (Figure 1b).
Figure 1. Model schematics. (a) TypeA bilayer system: tetragonal cell (lines), donors (P _{1}, P _{2}), periodic images (translucent circles), and effective donor layers (translucent sheets). Varying separation within bilayers (arrows). (b) Secondlayer dopant (inplane) positions: P _{1} projection (black circle), coplanar Si atoms (circles), typeA, B, and C positions, other monolayers’ atoms’ projections (dashed circles), and periodic boundary (square).
Methods
δ layers of P are created on Si (001) terraces before being epitaxially coated with further Si [2427]. It is easy to envision this coating process being monitored and halted at a desired buffer thickness, before a new δ layer of P is created (and/or patterned). Single δ layer findings [16] suggest that layers interact when less than 80 monolayers (approximately 10.9 nm) of silicon separate them, and that at 80 ML, their properties converge with respect to silicon cladding depth. In that model, periodic replications of the layers were identical by construction, with no possibility of any deviation. Here, we explicitly allow for such differences by including a second layer in the model.
c(2×2) cells including two δlayers at N ML separation and 80 ML of Si cladding were built (N ∈ {4,8,16,40,60,80}). Doping into a new layer can be accomplished at several locations [19]. For Nmod(4) = 0 systems, this can occur in three ways (Figure 1b): directly above the original dopant (type A), at either position nearest A in the plane (type B), or at maximal inplane separation (type C). Note that B is not the nearest neighbour of either A or C in the silicon lattice (see Figure 1b).
We note that Nmod(4) ∈ {1,2,3} systems exhibit new position types, requiring further modelling. Although such investigation would greatly inform the ongoing discussion of disorder in δdoped systems, due to computational resource constraints, they are not considered here.
Models were replicated as A_{N}, B_{N}, C_{N}, and undoped (for bulk properties comparison without bandfolding complication) structures. Electronic relaxation was undertaken, with opposite donor spins initialised for each layer and various properties calculated. The general method of [16] using SIESTA[28], and energy convergence of 10^{6} eV, was used with two exceptions: an optimised double ζ with polarisation (DZP) basis [19] (rather than the default) was employed for all calculations, and the C_{80} model was only converged to 2 × 10^{4} in density (and 10^{6} eV in energy) due to intractability. Band structures had at least 25 points between highsymmetry locations.
The choice of a DZP basis over a single ζ with polarisation (SZP) basis was discussed in [16], where it was found for single δ layers to give valley splittings in far better agreement with those calculated via planewave methods. In the recent study by Carter et al. [23], less resourceintensive methods were employed to approximate the disorderedbilayer system, however, here we employ the DZP basis to model the completely ordered system.
Results and discussion
Benchmarking of N = 80 model
Although we used the general method of [16], as we used the optimised basis of [19], we benchmark our A_{80} model with their 80 ML single δlayer (δ_{1}) calculation rather than those of [16]. (Lee et al. [18] also used the same general method.) Our supercell being precisely twice theirs, apart from having spin freedom between layers, results should be near identical.
Figure 2 is the A_{80} band structure. Agreement is very good; band shapes are similar, and the structure is nearly identical. A closer look reveals that A_{80} has two bands to the δ_{1}’s one, as we should expect – A_{80} has two dopant layers to δ_{1}’s one. Due to 80 ML of Si insulation, the layers behave independently, resulting in degenerate eigenspectra. Comparison of band minima shows quantitative agreement within 20 meV; the discrepancy is likely a combination of numerical differences in the calculations (generally accurate to approximately 5 meV), the additional spin degree of freedom (which may allow less repulsion between the layers), and band folding from the extension of the bilayer supercell in z.
Band structures and splittings
Band structures for other models were calculated in the same fashion. Comparisons of band minima are shown in Table 1. Within types, the band minima change drastically as N shrinks and the δ sheets come closer together. The natural progression of this is to the δ_{2} results [19], where two layers are directly adjacent (although the location of the dopant in the second layer will be different, as mentioned above, due to the nature of the silicon lattice).
Table 1. Bilayer models’ band minima energies, Fermi levels, and differences between band minima
In the largeseparation limit (N ≥ 40), the values across types (same N) are quite similar. The full band structures (60, 80 not shown here) are effectively identical from the valence band maximum (VBM) to well above the Fermi level. We focus upon the occupied spectra from VBM to E_{F}: as N decreases, differences due to small changes in donor position become apparent. In particular, we find (see Figure 3) that the C_{4} model exhibits drastically wider splitting between the first two bands than A_{4}, which in turn is significantly wider than B_{4}. N ≥ 40 models show occupation of four bands; a fifth (with minimum away from Γ) dips below E_{F} for N = 16 and 8. (For N = 4, the minimum shifts to be at Γ.) The tetragonal symmetry means that this fifth band is fourfold degenerate, so these models have four further, for a total of eight, channels open for conduction, until they merge by N = 4. These fifth bands, however, do not penetrate very far below the Fermi level and are henceforth ignored.
Figure 3. Band structure of N ≤ 40 models, from M to Γ to X . The valence band maxima have been set to zero energy.
As has been noted before [14,16], the specific ordering of donors and symmetries inherent in (or broken by) their placement have great effect upon band energies. Whilst for single layers, valley splitting was paramount [15,16]; here, we introduce the additional possibilities of Coulombic interaction with faraway dopants and quantum interactions with nearneighbour dopants.
Upon closer inspection, holding too closely to singlelayer valley splitting proves to be a somewhat naïve way of discussing of the band structures of Figure 3. When all models are compared from N = 80 down, it is easily seen that bands come in pairs in the bilayer models, and therefore, at N = 80, the equivalent of singlelayer valley splitting is the gap between bands one and three (type 2 in Table 1). Due to their large spatial separation, electrons inhabiting bands one and two will overlap only to a negligible extent and, hence, share the same energy here. (This type 1 separation corresponds to interlayer effects  see ‘Consideration of disorder’ section for further discussion.)
As N →4, however, the layers approach and interact; for the Ctype model, bands two and three quite clearly cross each other, and it is possible that some mixing of states occurs  which might well be utilised for information transfer between circuit components in a threedimensional device design; consider two wires crossing at close distance (N < 16) in order to share a state between them.
In fact, the differences columns of Table 1 show that the valley splitting is not particularly perturbed until the layers are quite close to each other (A_{4}, B_{8}, and C_{4}), whilst bands which are effectively degenerate at N = 80 are not for N ≤ 16. The layers are interacting, affecting the multielectronic wavefunction under these closeapproach conditions. At N = 4, it is currently impossible to say which contributes more to the band structure.
Within the approximate treatment in [23] it was concluded that the valley splitting in the interacting deltalayers is the same as that for the individual deltalayer. Here we find that in the DZP approach the valley splitting of 119 meV for the interacting deltalayers is about 30% larger than for the individual deltalayer [19]. Of course, Carter et al. themselves acknowledge that their reduced basis functions are not complete enough to represent the ideal system; the SZP results on disordered systems could not have predicted such a difference. We therefore suggest that their estimate of splitting of 63 meV be revised upwards somewhat; the 30% difference seen between ideal single and double layers may be thought of as an upper bound, since the influence of disorder may well counter that of introducing the second layer.
Density of states and conduction
Figure 4 shows the electronic densities of states (DOS) of the A_{N} models. As evidenced by the changes in the band minima, lower N leads to occupation further into the band gap. In all cases, the occupation is maintained across E_{F}, indicating that the structures are conductive. The DOS of highN models are in good agreement with each other, confirming that these layers are well separated, whilst those of smaller N show shifts of density peaks relative to each other and to A_{80}.
Figure 4. Densities of states of A_{N} models. A_{4} (solid black), A_{8} (dashed black), A_{16} (dotted black), A_{40} (solid red), A_{60} (dashed red), A_{80} (dotted red), and bulk (shaded grey). Twentyfive meV Gaussian smearing applied for visualisation purposes.
Less affected by donor placement than the band structure, the DOS shows negligible difference between types by N = 16 (Figure 5). Changes between the DOS of N = 1680 models (not shown) therefore arise solely from the interlayer distance. When one considers the interdonor separation length, consisting of N layers’ separation and a component describing the inplane separation due to model type, this separation length is far more sensitive to variations of type when the interlayer separation is short. At N = 4, there is already a significant scale difference between the two vector components’ magnitudes which is only exacerbated by increasing N.
Figure 5. Densities of states of (a) N = 4, (b) N = 8, and (c) N = 16 models. Types A (black solid lines), type B (blue dashed lines), type C (red dotted lines), and bulk (grey shaded backgrounds). Energy zero is set to the VBM, Gaussian smearing of 25 meV applied for visualisation purposes.
The perpendicular electronic crosssection
Electronic crosssections are inferred from the local densities of states (LDOS; integrated from VBM to E_{F}) and may be useful in planning classical devices. A_{N} models are shown in Figure 6a, where isolation of wellseparated and interaction between closely spaced layers are obvious. Significant density overlap begins between N = 8 and 16.
Figure 6. Local density of states: side view. (a) Charge density (by LDOS) of A_{N} models, lineaveraged along the [110] direction; (b) contour plot of C_{N} models’ Ψ_{gap}, maximum along [110] taken for each point. All data normalised to [0,1].
Figure 6b depicts the worstcase overlap of the gapstates’ wavefunction (modulus). By N = 40, we see (for quantum information purposes) nonnegligible overlap (>2%) between the layers. Conversely, N ≥ 80 models show that Ψ_{gap} falls off to less than e^{5}. By N = 8, Ψ_{gap} ≥ e^{2} between the layers. This information will be crucial in assessing future quantum device designs.
Interestingly, the falloff from the center of the N = 4 model is decidedly similar to the falloff of the wellseparated layers of the N = 80 model, as Figure 7 illustrates. The bilayer density is slightly higher in the central nanometre and almost negligibly lower in the tail regions. Unlike the δ_{2} model [19], which featured doping in two adjacent layers of the Si crystal, the charge density is not pulled inwards much more than a simple combination of two single layers would suggest.
Figure 7. Single layer versus bilayer density profiles. Average of A_{80} layer profiles about their centers (dotted black), A_{80} average profile shifted to center on bilayer positions (solid black), summed shifted profiles (dashed blue), and planeaveraged A_{4} profile (solid red). Data were planeaveraged, collapsed to [001], and normalised such that charge density integrated to one.
Inplane density maps
Inplane density maps will be of interest when considering transport and also when considering disorder. Figure 8a shows the inplane charge density for all models. Inplane alignment does indeed have a great effect upon the charge density; A_{N} models exhibit large lowdensity central regions (away from the donors) whilst B_{N} have highdensity pathways in one direction, and C_{N} show the greatest extent of highdensity regions.
Figure 8. Local density of states: topdown view. (a) Charge density (all models), lineaveraged along [001] and normalised such that their values’ ranges are each [0,1]. (b) Charge densities of N ∈ {4,80} models, normalised to Ψ^{2} = 1. Differences also shown, on two scales.
To focus on bilayerspecific effects, N = 4 and 80 models were rescaled, and their differences are shown in Figure 8b. The electronic density reorganises as the layers approach, in a typedependent manner. The magnitude of the rearrangement is ≤ 20% of the singlelayer density.
Consideration of disorder
As mentioned earlier, though the main focus of this work is perfectly ordered systems, recent attention has been given to disorder. Here, we consider how these ordered results can contribute to that discussion. As it is useful to recall which calculations have been previously performed in the literature, Table 2 summarises the state of the field and introduces terminology to distinguish between the various models.
Table 2. Listing of ab initio works in this field covering systems with 1/4 ML phosphorus density
Interacting δ layers have recently been studied from the point of view that current experimental systems involve some inherent level of disorder [23]. Whilst it is recognised that a complete DZP model of interacting quasidisordered bilayers is currently intractable (let alone incorporating disorder on any realistic scale), they offered the rational approach of contrasting a DZP model of a single quasidisordered δ layer against an SZP model and then extending the SZP model to cover a quasidisordered bilayer. The reasonable assumption there was that the differences between SZP and DZP models should be similar in both cases, and the valley splittings of the (missing) DZP model of a quasidisordered bilayer could thus be inferred. They also considered the approach of using a DZP basis and mixed pseudopotential to describe the disorder; this approach is vastly cheaper computationally and purports to inform us about the splittings due to the presence of the second layer. It is supported by SZP mixed and explicit pseudopotential results in which these interlayer splittings are preserved.
The approach taken in this paper, of calculating the properties of an explicitly ordered bilayer system using a DZP basis, complements that previous work. We can equivalently make comparisons between the ordered singlelayer systems of [19] (δDZPord) and ordered doublelayer systems as calculated with DZP bases here (δδDZPord), and between the δDZPord systems of [19] and the (DZP) quasidisordered singlelayer system (δDZPdis) presented in [23], in order to draw inferences about the (intractable, missing) δδDZPdis model, without at any stage compromising the accuracy of the results by using a lesscomplete basis set. (We shall now proceed to drop the ‘DZP’ from the labels, since it is ubiquitous here.)
One important point in the consideration of disorder from these ideal models is that, at the lowest separation distances, the crystalline order and alignment of the layers is greatly influencing their band structure. In a disordered system, the alignment effects would largely be negated, or averaged out, since one would expect to encounter all possible arrangements. We therefore limit ourselves to discussing averages of splittings.
The δord layers show valley splittings (VS) of 92 meV, as compared to the 120(±10%) meV of the δδord bilayer systems presented here (apart from separations of less than 8 monolayers). The δdis system showed a valley splitting of 63 meV, indicating that we might expect a reduction of valley splitting of up to 32% due to the (partial) inclusion of disorder. We can then infer that the valley splitting in the δδdis systems should be around 81 meV, unless their separations are small (see Table 3).
Table 3. Model properties and prediction of disordered splittings
We can estimate the interlayer splitting by taking the differences between bands 1 and 2 and bands 3 and 4 (except at low separation). Averaged values for these are also presented in Table 3. Unfortunately, beyond the SZP models, we have no further information as to the likely behaviour of the δδdis model at the DZP level in this regard, as there can be no interlayer splitting in the isolate singlelayer models to compare against.
It is clear from Table 3 that the estimated values for the valley splitting differ from those predicted by the SZP approach (63 meV for all but ‘extremely close separations’). We are in agreement with the finding that narrow separations affect the value greatly. Even allowing for the possibility of overestimation of the valley splitting here (the δord value was 92 meV) only adjusts the estimated δδord value by 8 meV, not the 20 required to match the values obtained using the SZP approach.
Obviously, the extension to a full DZP model has brought to light behaviours at small separation not evident from the SZP approach, and further work is required to elucidate these as computational resources improve.
Conclusions
We have modelled Si: δP bilayers, varying their separation and inplane alignment. Whilst layers behave independently at large separations (above 40 ML), they interact when brought close together: band structures are affected considerably; variation in the energy splitting between the first two occupied bands for N = 4 is considerable, and this variation must be taken into account in any future models of disorder in such closely spaced layers; inplane charge densities shift by ≤20%. Outofplane charge densites overlap to varying extent; wavefunction moduli are more sensitive. For 8 ≤ N ≤ 16, four new conduction channels open, making eight in total. Consequences for device design will depend heavily on the desired purpose; detailed information has been presented for several possible issues to facilitate successful design and operation of future threedimensional devices, be they classical or quantum in nature. Finally, despite single ζ with polarisation results indicating that valley splittings are the same in single and double δlayered systems, our results indicate otherwise at double ζ with polarisation level (previously shown to be adequately complete), with implications for the ongoing discussion of disordered systems of this type.
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
DWD, MCP, and LCLH planned the study. DWD, MCP and AB performed the calculations. All authors analysed the results and wrote the manuscript. All authors read and approved the final manuscript.
Acknowledgements
The authors acknowledge funding by the ARC Discovery grant DP0986635. This research was undertaken on the NCI National Facility, Canberra, Australia, supported by the Australian Commonwealth Government.
References

Weber B, Mahapatra S, Ryu H, Lee S, Fuhrer A, Reusch TCG, Thompson DL, Lee WCT, Klimeck G, Hollenberg LCL, Simmons MY: Ohm’s law survives to the atomic scale.
Science 2012, 335:6467. PubMed Abstract  Publisher Full Text

Fuechsle M, Miwa JA, Mahapatra S, Ryu H, Lee S, Warschkow O, Hollenberg LCL, Klimeck G, Simmons MY: A singleatom transistor.
Nat Nanotechnol 2012, 7:242246. PubMed Abstract  Publisher Full Text

Eisele I: Deltatype doping profiles in silicon.
Appl Surf Sci 1989, 36:3951. Publisher Full Text

Kiunke W, Hammerl E, Eisele I: Electrical transport between delta layers in silicon.
J Appl Phys 1992, 72(8):3602. Publisher Full Text

RodriguezVargas I, GaggeroSager LM: Subband and transport calculations in double ntype δdoped quantum wells in Si.
J Appl Phys 2006, 99:033702. Publisher Full Text

Cartwright J: Intel enters the third dimension.
Nature News 2011.
doi:10.1038/news.2011.274

Scappucci G, Capellini G, Klesse WM, Simmons MY: Dualtemperature encapsulation of phosphorus in germanium δlayers toward ultrashallow junctions.
J Cryst Growth 2011, 316:8184. Publisher Full Text

Scappucci G, Cappellini G, Johnston B, Klesse WM, Miwa JA, Simmons MY: A complete fabrication route for atomicscale, donorbased devices in singlecrystal germanium.
Nano Lett 2011, 11:22722279. PubMed Abstract  Publisher Full Text

Scappucci G, Capellini G, Klesse WM, Simmons MY: Phosphorus atomic layer doping of germanium by the stacking of multiple δ layers.
Nanotechnology 2011, 22:375203. PubMed Abstract  Publisher Full Text

Qian G, Chang YC, Tucker JR: Theoretical study of phosphorous δdoped silicon for quantum computing.

Carter DJ, Warschkow O, Marks NA, McKenzie DR: Electronic structure models of phosphorus δdoped silicon.

Ryu H, Lee S, Klimeck G: A study of temperaturedependent properties of ntype δdoped si bandstructures in equilibrium. In Proceedings of the 13th International Workshop on Computational Electronics. Beijing: Tsinghua University; 2009:14.
arXiv:1003.4926v1 [condmat.mtrlsci]

Ryu H, Lee S, Weber B, Mahapatra S, Simmons MY, Hollenberg LCL, Klimeck G: Quantum transport in ultrascaled phosphorusdoped silicon nanowires.

Carter DJ, Marks NA, Warschkow O, McKenzie DR: Phosphorus δdoped silicon: mixedatom pseudopotentials and dopant disorder effects.
Nanotechnology 2011, 22:065701. PubMed Abstract  Publisher Full Text

Drumm DW, Hollenberg LCL, Simmons MY, Friesen M: Effective mass theory of monolayer δ doping in the highdensity limit.
Phys Rev B 2012, 85(15):155419.
arXiv:1201.3750v1 [condmat.mtrlsci]

Drumm DW, Smith JS, Budi A, Per MC, Russo SP, Hollenberg LCL: Ab initio electronic properties of monolayer phosphorus nanowires in silicon.
Phys Rev Lett 2013, 110:126802. PubMed Abstract  Publisher Full Text

Smith JS, Cole JH, Russo SP: Electronic properties of δdoped Si:P and Ge:P layers in the highdensity limit using a ThomasFermi method.

Lee S, Ryu H, Campbell H, Hollenberg LCL, Simmons MY, Klimeck G: Electronic structure of realistically extended atomistically resolved disordered Si:P δdoped layers.

Budi A, Drumm DW, Per MC, Tregonning A, Russo SP, Hollenberg LCL: Electronic properties of multiple adjacent δdoped Si:P layers: the approach to monolayer confinement.

Fuechsle M, Mahapatra S, Zwanenburg FA, Friesen M, Eriksson MA, Simmons MY: Spectroscopy of fewelectron singlecrystal silicon quantum dots.
Nat Nanotechnol 2010, 5:502505. PubMed Abstract  Publisher Full Text

Drumm DW, Budi A, Per MC, Russo SP, Hollenberg LCL: Ab initio calculation of valley splitting in monolayer δdoped phosphorus in silicon.
Nanoscale Research Letters 2013, 8:arXiv:1201.3751v1 [condmat.mtrlsci].

Drumm DW: Physics of lowdimensional nano structures. PhD thesis, The University of Melbourne, 2012

Carter DJ, Warschkow O, Marks NA, Mackenzie DR: Electronic structure of two interacting phosphorus δdoped layers in silicon.

Tucker JR, Shen TC: Prospects for atomically ordered device structures based on STM lithography.

O’Brien JL, Schofield SR, Simmons MY, Clark RG, Dzurak AS, Curson NJ, Kane BE, McAlpine NS, Hawley ME, Brown GW: Towards the fabrication of phosphorus qubits for a silicon quantum computer.

Shen TC, Ji JY, Zudov MA, Du RR, Kline JS, Tucker JR: Ultradense phosphorous delta layers grown into silicon from PH_{3} molecular precursors.
Appl Phys Lett 2002, 80(9):15801582. Publisher Full Text

Fuechsle M, Ruess FJ, Reusch TCG, Mitic M, Simmons MY: Surface gate and contact alignment for buried, atomically precise scanning tunneling microscopyppatterned devices.
J Vac Sci Tech B 2007, 25(6):25622567. Publisher Full Text

Artacho E, Anglada E, Diéguez O, Gale JD, Garciá A, Junquera J, Martin P, Ordejón RM, Pruneda JM, SánchezPortal D, Soler JM: The SIESTA method; developments and applicability.
J Phys Condens Matter 2008, 20:064208. PubMed Abstract  Publisher Full Text