Solid-phase silicon homoepitaxy via shear-induced amorphization and recrystallization

The development of epitaxy techniques for localized growth of crystalline silicon nanofilms and nanostructures has been crucial to recent advances in electronics and photonics. A precise definition of the crystal growth location, however, requires elaborate pre-epitaxy processes for substrate patterning. Our molecular dynamics simulations reveal that homoepitaxial silicon nanofilms can be directly deposited by a crystalline silicon tip rubbing against the substrate, thus enabling geometrically controlled crystal growth with no need for substrate pre-patterning. We name this solid-phase epitaxial growth triboepitaxy as it solely relies on shear-induced amorphization and recrystallization that occur even at low temperature at the sliding interface between two silicon crystals. The interplay between the two concomitant, shear-induced processes is responsible for the formation of an amorphous sliding interface with constant nanometric thickness. If the two elastically anisotropic crystals slide along different crystallographic orientations, the amorphous layer can move unidirectionally perpendicular to the sliding plane, causing the crystal with lowest elastic energy per atom to grow at the expenses of the other crystal. As triboepitaxial growth is governed by the shear elastic response of the two crystals along the sliding direction, it can be implemented as a mechanical scanning-probe lithography method in which epitaxial growth is controlled by tuning the crystallographic misorientation between tip and substrate, the tip's size or the normal force. These results suggest a radically new way to conceive nanofabrication techniques that are based on tribologically induced materials transformations.

The development of nanometer-accurate techniques for localized epitaxial growth of singlecrystalline films and nanostructures is essential for modern electronic, optoelectronic and photonic applications. 1,2 To precisely control the location, shape and size of the grown crystalline structures, complex processes, based for instance on nanolithography, 3 are necessary to predefine the growth patterns on the substrate surface. 1 A direct deposition of crystalline patterns by means of mechanical scanning-probe lithography techniques 4 would help circumvent this issue as the deposition process would be localized at the contact between writing tip and substrate and, in addition, its accuracy would not be diffraction-limited.
However, the implementation of mechanical or tribological scanning-probe nanolithography in silicon generally relies on nanoscratching to remove surface atomic layers, 5,6 or to induce localized surface amorphization 7 or oxidation. 8 Despite producing nanoscale features on silicon substrates, these techniques have not been used for localized deposition of crystalline silicon films. Yet our previous studies on shear-induced amorphization of diamond 9 and silicon 10 suggest that the direct deposition of amorphous silicon (a-Si) films from a silicon tip onto a silicon substrate might be possible. Such a-Si films could be recrystallized in a following step by thermally induced solid-phase epitaxy. 11 In order to combine these two steps into a single mechanical process, we propose that the same plastic shear deformation that is responsible for amorphization can also enable recrystallization, since, similarly to temperature, it provides the atomic mobility necessary to overcome the energy barriers for recrystallization of the metastable amorphous phase. Indeed, as we reported recently, 10 a nanoscale amorphous layer forms at the sliding interface between two silicon crystals. A competition between shear-induced amorphization 9 and recrystallization 10 at the amorphous-crystal interfaces delimiting the a-Si region results in a constant thickness of the amorphous layer, whose position changes stochastically with no unidirectional drift if the two crystals have identical orientation. In this article, we use molecular dynamics (MD) simulations to explore conditions that render the movement of the amorphous layer unidirectional, thus leading to the growth of one crystal at the expense of the other one.
We find that this tribologically induced homoepitaxial process is possible if the two sliding crystals have different crystallographic orientations and we name it "triboepitaxy", a combination of nanotribology and solid-phase epitaxy. We show that, for a variety of surface orientations, the growth direction is determined by the elastic energy per atom : the crystal with lowest grows. Hence, triboepitaxial growth can be controlled by exploiting elastic anisotropy, i.e. by tuning the relative orientation of the two crystals with respect to the sliding direction. Finally, we propose that triboepitaxy could enable a novel scanning-probe nanolithography technique for the deposition of crystalline silicon nanostructures with no need for elevated temperatures or contaminating species to catalyse recrystallization. 3 In this case, triboepitaxial growth can also be controlled by exploiting elastic finite-size effects, e.g. by tuning shape and size of the writing tip.
We perform classical reactive MD simulations 12,13 to study the evolution of a shear-induced a-Si phase between two diamond-cubic silicon crystals. 9,10 Different interatomic potentials, [13][14][15] including a recent potential based on machine learning of density-functional-theory results, 15,16 yield comparable results (Supporting Information S1). Figure 1a displays the interface between two silicon (110) surfaces before sliding along the [ direction (periodic boundary conditions are applied in the sliding plane; analogous simulations for surfaces are described in Supporting Information S6.1). After imposing an external normal pressure and thermalizing the system to room temperature by means of specially tailored barostat and thermostat, 9 a top layer of atoms is rigidly moved at a constant velocity along the sliding direction, while a bottom layer is kept fixed.  Upon sliding, an a-Si layer forms (snapshot in Figure 1b) at the interface between the two crystals. It entirely accommodates the plastic shear deformation of the system 10 and its thickness increases with sliding distance and with applied normal pressure (Figure 1c for ). As previously reported, 10 the crystalline-to-amorphous transition is mechanically induced, local temperatures remain below K and a melting transition can be excluded. Interestingly, as the sliding distance increases, saturates at a constant value . Since mechanical amorphization alone would lead to a strict monotonous increase of , 9,10 saturation indicates the presence of a competing process, namely shear-induced recrystallization. 10  A further evidence of the recrystallization process is provided by the blue line in Figure 1c.
Here, we perform a sliding simulation at high pressure ( GPa) and suddenly decrease the normal pressure to 5 GPa at a sliding distance 200 nm. This causes a rapid decrease of the a-Si thickness from nm to nm, i.e. a 2.1-nm growth of the crystals. In accordance with the dependency , is exactly the a-Si thickness the system would have reached had the simulation been performed at from the beginning.
While no net crystal growth occurs for symmetric tribopartners (i.e. same crystal orientation and sliding direction), we now show that shear-induced growth of one crystal at the expense of the other can be achieved by breaking the symmetry of the sliding system. Figure 1e displays a Si(110) sliding system in which the lower crystal is rotated by with respect to the upper one, so that sliding proceeds along different crystallographic directions for the two crystals. Now, the lower crystal grows rapidly, as shown by the upwards migration of and in Figure 1f, while remains roughly constant at (Movie S1). Growth is conveniently characterized by the rate , which is and approximately constant in Figure   1f.
To elucidate the mechanism underlying the observed triboepitaxial growth, we first focus on the shear stress . Although the a-Si region migrates, the time-averaged shear stress remains constant. Since also and remain roughly constant, we conclude that the properties of the crystal-amorphous-crystal transition region remain unaltered throughout the simulation. Figure 2a shows a saw-tooth-shaped characteristic of stick-slip instabilities.
During stick phases, the system deforms elastically. When exceeds the interfacial shear strength, the crystal-amorphous-crystal transition region undergoes plastic slip and its atoms gain mobility. In Figure 2a, we follow the typical trajectory of an atom directly involved in triboepitaxy during several stick-slip phases (Movie S2). Initially, the atom is part of the upper crystal. A series of plastic slip events cause the atom to be dragged into the a-Si region and to move to the lower amorphous-crystal interface where it eventually recrystallizes. To understand what drives the migration of the a-Si interface, we note that triboepitaxy resembles grain-boundary migration in metals under shear [17][18][19] and that thermally induced migration of grain boundaries under anisotropic elastic strain is driven by the elastic energy density difference between grains. 20,21 We show in the following that also determines the growth direction in shear-induced migration. An arbitrary steady-state configuration from the MD trajectory in Figure 1f is selected and relaxed to and  In conclusion, the triboepitaxy concept presented in this work lends itself to immediate experimental validation. Since triboepitaxial growth is simply governed by the elastic properties of the two sliding crystals, predictions on the crystal growth direction are straightforward and the process can be controlled by easily accessible parameters and implemented using existing technologies. 4 A first step in this direction could be made by using atomic force microscopy probes inside a transmission electron microscopy. 24 Adhesion experiments performed under these conditions suggest that surface-passivating species that prevent the formation of covalent bonds across the tribological interface are easily removed by sliding. 24 To further investigate the generality of the proposed mechanisms, simulations are underway to investigate triboepitaxy at sliding contacts between two differently oriented surfaces (Supporting Information S6.4). In such cases, due to the different elastic responses of the two crystals to normal pressure, the applied normal force can become an additional parameter to control crystal growth direction and rate. Since one of the main ideas behind triboepitaxy is that shear-induced plasticity can replace temperature to help atoms overcome energy barriers for recrystallization, the method could be tested on other materials showing shear-induced amorphization (e.g. diamond 9 ) or materials on which solid-state epitaxy was applied successfully (e.g. germanium 11 ). More generally, this study reveals a new way to apply tribological concepts in the context of nanolithography and nanofabrication. While mechanical scanning-probe techniques are usually based on nanomachining, 4 i.e. removal of material from a surface, we propose that tribologically induced phase transitions can be exploited for direct deposition of nanostructures or selective microstructural modification of a surface.

Computational Methods
Interatomic potentials. All atomistic simulations described in the article are performed with a screened version of the bond-order potential by Kumagai and coworkers. 12,13 Simulations employing different interatomic potentials that can describe shear-induced amorphization/crystallization processes and elastic anisotropy deliver analogous results (see Supporting Information S1 for a comparison to the Stillinger-Weber potential 14 and a more a recent potential based on machine learning of density-functional-theory results (Gaussian Approximation Potential, GAP) 15 ). interatomic potentials that can suitably describe these processes deliver results that are analogous to those described in the article. Figure S1 summarizes the results of molecular dynamics (MD) simulations of triboepitaxial growth that we obtained using the Stillinger-Weber potential 4 and a more sophisticated machine-learning-based potential (Gaussian Approximation Potential, GAP) 5 . . The simulation temperature is K and the normal pressure is

Molecular dynamics sliding simulations.
. The system size for the Stillinger-Weber simulation is identical to the one used for the screened Kumagai simulation (see Table S1), but rescaled to the optimized Stillinger-Weber lattice constant Å. For the GAP potential the system size is reduced to Å 3 ( atoms, adatoms, lattice constant Å) due to the higher computational cost. For the same reason, a higher sliding velocity is used in this case. The snapshots show the atomic configuration at different , the colors of the atoms refer to the crystal the atoms belong to at .

S2 Overview of the systems used for the simulations of periodic interfaces
Crystal 1 Crystal 2 System size Strain crystal 1 Strain crystal 2 Atoms *The results for the -system are obtained using the setup with atoms, except for the those shown in Fig. 2d, where the symmetric setup ( ) with atoms is used in order to ensure that the system is perfectly inversed at Table S1. Overview of the systems used for the simulations of periodic interfaces. The individual crystals are characterized using the notation , where the integer numbers between round brackets are the Miller indices of the sliding surface while those between square brackets indicate the sliding direction. refers to the length of the cell vector along the shear direction while refers to the direction in the interface plane normal to the shear direction. In each system, adatoms are initially added to facilitate the phase transitions. Negative values indicate a compressive strain.

S3 Role of strain due to lattice mismatch in periodic simulations
In general, in order to fit two crystals with different sliding direction in the same periodic simulation cell, the crystals have to be strained along the directions of the sliding plane to align their periodicities. This spurious strain changes the elastic response of the crystals and thus affects . To mitigate this effect, we choose multiples of the unit cells such that the necessary strain for the alignment is minimal compatibly with an affordable computational cost. Moreover, we distribute the strain equally in both crystals (Table S1). Although the magnitude of the mismatch strain in our sliding systems is relatively low, we verify that the interface migration is not an artifact of the mismatch strain. We consider the system and perform two additional MD simulations: one in which the mismatch strain is fully applied to the upper crystal and another one for the inverse case. Figure S2 shows the position of the lower amorphous-crystal interface as a function of the sliding distance.
Although imposing the spurious mismatch strain slightly increases the elastic energy of the strained crystal, this has no significant influence on triboepitaxial growth and in all three cases the growth direction and the growth rate are the same.
This test is especially interesting for the case in which the whole strain in accommodated in the -oriented crystal. Here, without shear, the -oriented crystal is energetically more stable than the -oriented crystal. If this interface migration was driven by temperature, without applied shear, we would expect a growth of the -oriented crystal. However, with increasing , the sign of changes and the -oriented becomes the most stable crystal and grows. To circumvent any artifact, the mismatch strain is equally distributed in both crystals.
An additional test shows that the sign of from quasistatic simulations of the tribosystem including strain (Fig. 2c) is identical to the sign of as obtained by simple shear deformations of the separate crystals (without any strain induced by lattice mismatch). This confirms that the energetic ordering of the crystals under shear is not affected by mismatch strain. Figure S2. Role of strain due lattice mismatch in periodic simulations. Position of the lower amorphous-crystal interface as a function of for the -sliding system for different distributions of the strain necessary to align the periodicities of the upper and the lower crystal. The strain is either equally distributed between both crystals ("Equal strain") or accommodated entirely in one of the crystals, while the other crystal remains in its optimal configuration ("Strain in / "). The normal load is , K and m/s.

S4 Dependence of triboepitaxy on the sliding velocity
Irrespective of the sliding velocity, the magnitude of the crystal growth depends solely on the sliding distance for given temperature and normal pressure. Figure S3a

S5 as a function of the sliding distance
In Fig. 2c and 2d, is evaluated by selecting a snapshot from the MD trajectory after the steady-state is reached. In order to show that the choice of the particular steady-state snapshot does not affect the sign of (and even its magnitude), we perform the quasistatic shear simulations of Fig. 2c for equidistant snapshots extracted from the MD trajectory between and nm. For each simulation, we determine three energy levels at as shown in Fig. S4a: (i)   (solid line) as obtained by quasistatic shear simulations. The initial configuration for this calculation is the atomic configuration at nm (Fig. 1f). The gray dashed line is a guide to the eye highlighting that at both crystals are energetically degenerate. The dashed blue line refers to the energy level of the lower crystal under shear, the dashed red line to the one of the upper crystal. The inset is a zoom-out, where additionally the maximal energy density of the amorphous zone is marked with a dashed purple line. Panel (b) shows the evolution of these energy levels as a function of the sliding distance . All energies are determined at the average shear stress of the dynamic simulation.

S6.1 Shear-induced amorphization and recrystallization for -oriented crystals
To show that concurring amorphization and recrystallization processes also occur for other surface orientations and give rise to a steady-state with constant , we reproduce the results described in Fig. 1 for a sliding contact between two -oriented surfaces (Fig. S5). Despite different amorphization/crystallization kinetics, the results are analogous to those we obtained for two -oriented crystals. In particular:  Upon sliding with different external loads, the symmetric system (i.e. same sliding direction for both crystals) undergoes mechanical amorphization and the thickness of the interfacial disordered layer increases with increasing ;  grows with increasing sliding distance until it saturates at a constant value ;  is a characteristic feature of the system under shear for given normal load and temperature;  Continuously ongoing crystallization and amorphization events lead to a stochastic movement of the amorphous-crystal interfaces for the symmetric system;

S6.2 Triboepitaxy on other surface orientations
To assess the generality of the triboepitaxial growth process we repeat the periodic simulations of Fig. 1e,f for additional surface orientations. For -oriented interfaces amorphization is strongly suppressed under the considered pressures due to the high stability of the shuffle plane 1 (similar to the diamond case 6 ) and we do not observe interface migration.
For , and surface orientations, tribo-interfaces consisting of two crystals with the same orientation but sliding along different crystallographic directions are constructed as shown in Fig. S6a. For all cases, we observe the formation of an interfacial amorphous layer with a constant height that migrates upwards with an approximately constant rate (Fig. S6b). Figure S6c shows the energy density difference between the upper and the lower crystal as a function of the shear stress . We find that for any the upper crystal is energetically less stable than the lower crystal and that the growth direction correlates with the sign of . Apart for one case (red line), we also observe a correlation between and . We note, however, that different amorphization/crystallization kinetics on different crystallographic orientations 7,8 are likely to affect the growth rate and, therefore, the correlation between growth rate and elastic energy difference.

S6.3 Suitable orientations for triboepitaxial growth in an experiment
The most effective way to grow crystalline structures on a silicon substrate using the experiment suggested in Fig. 3d is to use tip/substrate systems in which both the finite size effects of the tip and the bulk elastic anisotropy favor the growth of the substrate crystal. One such example is the system presented in the article: both surfaces have a orientation; the oscillatory motion for the tip is along the direction and the equivalent direction, while the oscillatory motion for the substrate occurs along the and the equivalent direction.
It is important to point out that triboepitaxy is not restricted to these specific crystallographic orientations. Here, we discuss other suitable orientations for tip/substrate systems, where both the tip and the substrate have the same crystallographic surface orientation.
The case with perhaps the highest relevance for electronic applications is the triboepitaxial growth on surfaces. One possibility to grow such surfaces is to use a -oriented tip with oscillations along the and directions for the tip and along the and directions of the substrate. The elastic anisotropy alone leads to a rather slow growth rate of about if compared to other surface orientations (see Fig. S6b). Nonetheless, with the typical experimental velocities nm/s 9,10 significant growth velocities nm/s can be achieved, which could be further increased by finite size effects. Alternatively, we find in preliminary simulations of silicon contacts with different surface orientations that the periodic -system shows a higher rate and could also be used to grow surfaces.
The so far discussed low-index surface orientations are symmetric with respect to an inversion of the sliding direction. This is not necessarily the case for higher index surfaces. Therefore, to be able to efficiently perform triboepitaxial growth on such surfaces using the experiment in Fig. 3d, one should make sure that the growth occurs in both sliding directions during oscillatory motion. Figure S7 shows how the energy per atom as a function of changes upon inversion of the deformation direction for the higher index surface orientations considered previously (Fig. S6).

S6.4 Triboepitaxy for interfaces between crystals with different surface orientation:
Pressure as an additional parameter to control the growth direction All simulations described in the article consider interfaces consisting of surfaces with identical crystallographic orientations. Preliminary simulations show that triboepitaxy can also be observed by pairing surfaces with different crystallographic orientations. In this context, the combination of crystals with very different elastic responses to normal load is particularly interesting, as this directly influences .
Based on simple shear deformations of crystals under compressive strain, we identify two systems that show an inversion of upon an increase of : A sliding system and a -sliding system. In both cases, for the crystal grows (dark green curves in Fig. S8a,b) -as one would expect from the sign of (Fig. S8c). With increasing normal load, and the migration rate decrease until both quantities eventually become negative and the crystal growth direction is inverted (red curves in Fig. S8a,b). This reveals that due to different responses of the two crystals to the compressive load, the latter can be used to control the growth direction. Figure S8c shows that and the migration rate are correlated and that in the system the growth direction is fully determined by the sign of . However, for the -case, the curve does not pass through the origin.

S7 Robustness of triboepitaxy with respect to temperature and pressure variations
To understand how triboepitaxy depends on temperature and normal load, we perform two additional series of simulations using the tribological systems shown in Fig. S9a.
Firstly, we reduce the target temperature of the thermostats from K (Fig. S9b) to K (Fig.   S9c) while maintaining the same normal load . In the latter case, the temperature at the interface is typically around -K. Since the triboepitaxial process is mechanically driven, it should take place even at very low temperature as long as the shear rate in the interface amorphous region is finite (the average shear rate is given by ). For all interfaces, we observe that the decrease in temperature does not affect the growth direction, while it causes a decrease in the growth rate.
Secondly, we increase the normal load from to GPa (Fig. S9d), while keeping the temperature constant at K. We find that the growth rates are not significantly affected by the pressure variation. Since both crystals in each interface have the same surface orientation, the elastic response to normal load is approximately identical. This can change drastically when crystals with different surface orientations are paired (a discussion can be found in Section S6.4).
The important aspect is that in all examined cases the growth direction is determined by the sign of as calculated by quasistatic shear simulations in spite of substantial variations of temperature and normal pressure. In conclusion, the reported phenomena are not restricted to a narrow window of specific operation conditions and we expect them to be robust with respect to temperature and normal load variations. Figure S9. Temperature and pressure dependence of triboepitaxy. Vertical position of the lower amorphous-crystal interface for the systems depicted in (a) as a function of the sliding distance . Panel (b) is a reproduction of Fig. S6b. Panel (c) shows the results when the target temperature of the thermostat is reduced to K. The inset depicts the same curves for a longer sliding distance to show that triboepitaxy occurs for all cases after sufficiently long sliding distances. In panel (d) the normal load is increased to GPa with K. The sliding velocity is m/s for all simulations. The color code of (a) applies to (b)-(d).

S8 Analytical model for the effective shear modulus of the tip
The analytical expression for the effective shear modulus of the tip as a function of its length consists a of bulk contribution and a surface contribution . The surface thickness (Fig. 3c) is the thickness of the tip's surface region in which the elastic response differs from the bulk elastic response, i.e. the response of the infinite tip ( ). In Ref. 11, where an analogous model is proposed for quartz plates, is set to Å to separate the surface and bulk contributions to the effective Young's modulus of the system. To estimate an upper bound for , we refer to the typical decay length of the stress/strain field ahead of an atomically sharp crack in a silicon crystal of about Å 12 .
Applying the constraint Å < < Å, and considering that can only be varied by multiples of the distance between equivalent surface planes ( Å) to allow for the calculation of , the best fit to the values calculated by means of quasistatic shear deformation is provided by Å (Fig. S10). Due to the simplicity of the analytic model, which just serves as a qualitative model for finite-size elastic effects, the resulting value is a qualitative estimate, albeit physically reasonable. This, however, does not affect the conclusions drawn in the article, which are uniquely based on the simulations results. Figure S10. Effective shear modulus of the tip (Fig. 3c) as a function of its length . Black squares show from quasistatic shear deformation. The curves show the interpolation between the shear modulus of the bulk region GPa and the shear modulus of the surface region for different values of the surface thickness . The solid curve shows the best fit to the data points for Å.