3D Ordering at the Liquid–Solid Polar Interface of Nanowires

The nature of the liquid–solid interface determines the characteristics of a variety of physical phenomena, including catalysis, electrochemistry, lubrication, and crystal growth. Most of the established models for crystal growth are based on macroscopic thermodynamics, neglecting the atomistic nature of the liquid–solid interface. Here, experimental observations and molecular dynamics simulations are employed to identify the 3D nature of an atomic‐scale ordering of liquid Ga in contact with solid GaAs in a nanowire growth configuration. An interplay between the liquid ordering and the formation of a new bilayer is revealed, which, contrary to the established theories, suggests that the preference for a certain polarity and polytypism is influenced by the atomic structure of the interface. The conclusions of this work open new avenues for the understanding of crystal growth, as well as other processes and systems involving a liquid–solid interface.


DOI: 10.1002/adma.202001030
Understanding the nature of interfaces is key to determining and engineering their properties. However, while liquid-solid interfaces are fundamental to many applications and interactions, relatively little is known about them at the atomic level. Here, we use vapor-liquid-solid (VLS) grown nanowires (NWs) as a model system to identify and interpret the atomistic nature of a liquid-solid interface in three dimensions. VLS is a widely used method of crystal growth relying on the precipitation of a solid from a supersaturated nanodroplet. [1] Recent in situ investigations have revealed that precipitation/growth proceeds in a layer-by-layer fashion, similarly to thin film epitaxy. [2,3] Growth of NWs can occur in different growth directions, crystal structures, and polarities, [4][5][6] currently understood by relating the prevalence of these properties to the wetting characteristics of the liquid on the tip of the NWs, [5,7] where the role of the liquid droplet contact angle is widely discussed but remains controversial. [4,8] So far, the liquidsolid boundary is considered as a clear-cut, binary interface, with no regard for its structure. Modeling has mostly reasoned in terms of macroscopic parameters, including contact angle, the surface energies at the solid-vapor and liquid-solid interfaces and the chemical potentials. [9,10] The atomistic nature of the participating parts has rarely been considered.
The idea that a solid surface can induce local order to an adjacent liquid phase has previously been demonstrated in simulations of interfaces between solids and aqueous solutions, as well as liquid metals in contact with their solid compounds. [11][12][13] In addition, there have been some experimental attempts to study such ordering. [14][15][16] In particular, grazing incident X-ray diffraction has been used to study the ordering of liquid AuSi alloy, [17,18] and in situ X-ray diffraction has been employed to study the liquid ordering at the interface between P-terminated crystalline InP substrate and liquid InAu phase. [19] While this second study was targeted at understanding VLS-NW growth, due to the geometrical difficulties for probing NWs using X-rays, it was limited to investigating an interface in which a bulk liquid completely covers a bulk substrate, as opposed to the measurement of a nanoscale droplet in contact with the crystalline NW. It further applied an assumption of a fixed spacing between the ordered liquid layers. In addition, the composition of the liquid phase was tuned to match that observed in the gold-catalyzed growth of InP NWs; therefore, it is not clear if the results can apply to the self-catalyzed growth of III-V NWs. Using electron microcopy, Kaplan and co-workers studied the longitudinal ordering of liquid Al at the interface with Al 2 O 3 , which can extend up to five layers depending on the facets. [14] However, thus far, the effect of such ordering has been ignored in the VLS growth.
In this work, we combine scanning transmission electron microscopy (STEM), electron energy loss spectroscopy (EELS), atomic simulations using machine learning potentials and molecular dynamics (MD) to demonstrate, for the first time, the polaritydependent 3D nature of the liquid ordering at the solid interface.
The study is performed in the frame of the VLS growth of NWs. We elucidate the interplay between the ordering and the formation of a new bilayer, the preference for a certain polarity and the effect of ordering on defect formation. While our study concentrates on the Ga-assisted growth of GaAs, the results can extend to other relevant cases such as Au-assisted NW growth, bulk crystal growth from the melt, and even liquid-solid interaction in the context of biological environments.
We start by providing experimental evidence of liquid ordering at the interface with the solid in a VLS-grown NW using aberration-corrected STEM. Our study focuses on the Ga(l)-GaAs(s) interface of A and B polar GaAs NWs obtained by the Ga-assisted method. [20,21] A schematic drawing of the process is provided in Figure 1A. A Ga droplet preferentially gathers As 4 molecules, which dissolve in the droplet. The chemical bonding between As and Ga results in the growth of GaAs at the bottom of the droplet. [20] Along the [111] direction, GaAs is organized as indivisible bi-layers of As-Ga pairs, known as dumbbells. Due to their different electron affinity and the relative position with respect to the (111) surface, each bilayer exhibits an intrinsic electric field, called polarity. Positive (A) and negative (B) polarities correspond to termination by Ga and As, respectively.
Representative aberration-corrected high-angle annular darkfield (HAADF) STEM images of the interface for two NWs observed along the [11 0] zone axis of zinc-blende (ZB) and [112 0] of wurtzite (WZ), with A and B-polarities, are shown in Figure 1B,C. We mark the As and Ga atoms in blue and red, respectively. Polarity determination for these NWs is offered in Figure S1 of the Supporting Information. The A-polar NW exhibits a pure ZB structure, whereas the B-polar NW exhibits a mixed-phase structure, finishing with WZ. We attribute this to the change in conditions upon growth termination and the known tendency toward polytypic growth along this polarity. [5] Thanks to their intensity dependence on the atomic number as ≈Z 1.7-2 , HAADF images provide an easily interpretable and chemically sensitive representation of the structure. As expected, the solid shows strong peaks corresponding to the atomic columns, whereas, further from the interface, the liquid shows a uniform intensity because of its disordered nature. However, close to the interface, the liquid shows intensity fluctuations corresponding to an ordering. For the A-polar interface, this is limited to one additional layer, while several layers are observed for the B-polar interface. This is confirmed in the integrated HAADF intensity profiles shown in the right of panels of Figure 1B,C and we attribute it to a longer-range ordering on top of the B-polar interface. The first layer seems to be more clearly structured, with further layers becoming gradually more amorphous. This is in contrast to the findings by Krogstrup et al., [19] in which the first three layers show similarly strong ordering, with an abrupt reduction in ordering for the fourth layer. One should note that, while high-resolution transmission electron microscopy has been employed before to visualize the ordering of a liquid in contact with a crystalline solid, [14] aberration-corrected HAADF-STEM offers a more directly interpretable image contrast, free of both thickness/focus related contrast inversions and delocalization effects. This derives from the incoherent nature of this imaging mode.
We further analyze the liquid-solid interface using EELS hyperspectral mapping to probe the bulk plasmon response around the interface. Being related to the valence/free electron density and characteristics of a material, this technique can discriminate between different chemical/structural phases. Compared to analysis by core-loss EELS or energy dispersive X-ray spectroscopy, the required electron beam dose is orders of magnitude lower, allowing us to apply the technique without destroying the fragile structure of the interface.
Low-loss EELS spectrum images of the interface were recorded using an atomically sized probe. As an example, Figure 2A shows a HAADF-STEM image acquired simultaneously with the EELS signal from a spectrum image of an A-polar GaAs NW in contact with a Ga droplet. Figure 2B depicts the corresponding spectral evolution integrated across a region of the interface, represented by the blue area shown in Figure 2A.
The bulk plasmon excitation shifts smoothly from a broad peak in the GaAs to a sharper peak in the liquid Ga, over a spatial distance of ≈10 nm. Intensity oscillations in the GaAs corresponding to the atomic planes are clearly observable due to signal convolution with an elastic scattering contribution. Reference spectra extracted away from the interface at positions deep in the Ga and GaAs phases (marked P1 and P3 in Figure 2A), and a spectrum extracted from the interface (at position P2 in Figure 2A) are presented in Figure 2C,D, respectively. The spectrum at P2 might be assumed to be a linear combination of the contributions of the liquid and solid phases. However, this cannot correctly describe the P2 spectrum: a non-negligible residual is found, as presented in Figure S2 of the Supporting Information. We interpret this residue as the interface contribution. To quantitatively study this interface contribution, we use independent component analysis (ICA), which identifies a finite number of spectral components that can reconstruct the spectrum image. [22][23][24] Although care must be taken when interpreting the components obtained from ICA, it has been successfully applied to unmix independent components in low loss EELS data. [22] Once scaled in intensity, the first two components in our ICA analysis, shown in Figure 2C and labeled S and L, match spectra P1 and P3 well and can be safely interpreted as representing the solid GaAs and liquid Ga components. A third ICA component is also found, which is interpreted as the "ordered liquid" (OL) contribution, dominating the P2 interface signal, see Figure 2D. A reconstruction of the entire dataset using this three-component model shows a negligible residual, demonstrated in Figure S2 of the Supporting Information, confirming that these three components are sufficient to describe the data. This is further emphasized from an analysis of the scree plot in Figure S2 of the Supporting Information, a common tool in multivariate statistical analysis used to identify the number of significant sources of information in a dataset, which shows that the first three components represent the majority of variations in the spectrum image. [24] Considering the confined width of the OL region between GaAs and Ga, combined with the delocalized nature of the plasmon oscillations, the low-loss EELS response of the ordered liquid region cannot be isolated by a direct measurement. While the effect of an ordered liquid region on the overall low-loss EEL spectra of liquid-solid interfaces has been shown previously, [25] our use of ICA isolates its response for the first time. A similar analysis of B-polar WZ interfaces is demonstrated in the Supporting Information, leading to an ordered liquid spectral contribution similar to the A-polar component in the range around the plasmon peak. From these results, it is possible to determine a "bulk plasmon" energy for the OL of 14.6 eV (compared to 13.7 and 15.9 eV for Ga and GaAs), with a full width at half maximum (FWHM) of 2.35 eV (compared with 1.95 and 5.3 eV for Ga and GaAs, respectively). Thus, the peak position and FWHM of the plasmon for the OL sit between those of the liquid Ga and solid GaAs. Because of the single, Lorentzian-shaped bulk plasmon peak for each component, the response is consistent with that given by the Drude model, in which the peak position and its FWHM can be linked to the free electron density, n, and its damping constant, Γ. [26] The results indicate that both n and Γ decrease from the solid GaAs through the ordered liquid to the liquid Ga. The extent of damping depends on the lattice and band structure. For instance, damping can be caused by the transfer of energy from the plasmon resonance to single-electron transitions (i.e., creation of e − -h + pairs). Conversely, the closer the electronic nature of a material approximates a free electron gas, the smaller the damping constant is. It is therefore logical that the solid semiconducting GaAs NW has a larger Γ than the metallic liquid Ga, since the lattice and band structure of the former will increase the probability of energy transfer to single electron transitions. A similar difference is seen between semiconducting Si and metals such as Al.
From this analysis, one can conclude that the ordered liquid has a distinct electronic nature, which can intuitively be correlated with its semistructured nature. Considering this and given the fact that the ordered liquid cannot exist as a bulk phase outside of the interface, we propose that it can be considered as an interfacial complexion. [27][28][29] The existence of complexion could have profound consequences for crystal growth. To understand this, the 3D nature of the ordering and the influence of the polarity on its structure are key aspects, which we shall address next.
We now investigate the nature of the ordering using MD simulations. For this purpose, we use the 2D projection of the structure in STEM imaging to validate the MD results.
Given the size and time scales needed for these simulations, brute-force ab initio MD is prohibitively expensive. To circumvent this, we have trained a neural network potential (NNP) [30,31] based on a relatively small number of reference density functional theory (DFT) calculations. [32,33] Then, we have used this machine-learning model to drive the dynamics. [31,34] MD details in addition to the calculations and validation of the NN based on a multiple-timestep integration [35,36] are discussed in the Supporting Information. Figure 3A,B compare the projected linear density obtained across the simulation cell with intensity line profiles derived from experimental images for A-and B-polar cases, respectively. For any given polarity, simulations do not show significant differences in the liquid ordering with respect to the ZB or WZ solid phases. The experimental curves are obtained by projecting the intensity profiles from Figure 1B,C along the axis normal to the surface. In the solid, we obtain regularly ordered peaks at the positions of the dumbbells. The liquid also exhibits some peaks in the density profile, characteristic of atomic-level ordering. The range of the ordering is different for the two polarities. As noted earlier, in STEM images, while the B-polar order is observed for four layers, it does not extend beyond the first layer for the A-polar case. Similar to the experimental observations, the simulations show that the ordering gradually diminishes when the distance from the crystalline phase is increased, which is in contrast to the conclusions given by Algra et al. [19] for InP in contact with liquid InAu. This tendency correlates well with the simulations, which show that the in-plane order of the A-polar case is more diffuse than that of B-polar, as demonstrated in Figure S6 of the Supporting Information. Note that an exact correlation of the projected linear density obtained from the simulations with the HAADF-STEM integrated intensities is not expected: although the latter scales approximately with the former, a full quantum mechanical image simulation using the atomistic model derived from MD would be necessary for a quantitative comparison. The key observation, however, is that there is a quantitative match between the MD predictions of the spacing between ordered layers with the experimental observations. This is demonstrated in Table S3 of the Supporting Information, in which the experimental observations and MD predictions of the spacing between ordered layers are compared for a B-polar interface. We therefore conclude that MD based on the NNP correctly predicts the nature of the liquid ordering.
The 3D representations from the simulations indicate that the order in the liquid extends fully within the plane parallel to the interface, as illustrated in Figure S6 of the Supporting Information, consistent with the fact that this order is visible in HAADF images. One interesting observation is that contrary to what is assumed by Algra et al., [19] the spacings between consecutive ordered liquid layers are not identical. As clear from Table S3 of the Supporting Information, these layers gradually become closer to each other as the distance from the crystalline solid is increased. We further find that the order at the liquidsolid interface is dominated by the tendency of Ga atoms to form Ga-Ga pairs. Ga atoms form pairs with the solid when the surface is Ga-terminated, while it forms Ga-Ga pairs in the liquid when the solid is As-terminated. This tendency for dimerization is well-known in liquid Ga, which is often referred to as a molecular metal. [37] We should note that besides this remarkable ordering at the liquid-solid interface, the liquid remains a liquid as the position of the atoms is dynamic and thus varies with time. In Figure SI9 of the Supporting Information, we provide a comparison of the Ga-Ga radial pair distribution functions of the solid, liquid, and interface phases as an illustration.
The simulations shown previously were performed at 300 K, corresponding to the temperature used to acquire the STEM images. To make the link with the growth process, we performed the same simulations at the growth temperature (900 K). Figure 4 depicts the ordering of the interface at 900 K for both polarities where the spatial variation in the Ga atomic density is indicated by red isosurfaces. Figure 4A depicts a 3D view of the ensemble for a B-polar WZ solid in contact with liquid Ga. Figure 4B provides top and side views for the A-and B-polar solids. In both cases, we find ordering within the plane, however with significant differences in the arrangement of Ga atoms as a function of the polarity.
On the B-polar surface, Ga is adsorbed right on top of the terminal As. This is consistent with the large electronegativity difference between As and Ga that is likely to induce strong electrostatic interactions. In-plane ordering is also present at the A-polar surface. Here, instead, Ga atoms in the liquid pair with the terminal Ga atoms in the solid, consistent with the dimerization tendency of Ga. As expected, at the higher temperature, the range of the order has decreased relative to observations at 300 K, with simulations indicating the presence of at most one ordered layer on top of the solid. Most of the qualitative features of the liquid ordering, however, are preserved between the two temperatures, suggesting that experiments performed at 300 K can provide insights into the relevant mechanisms for growth at higher temperatures.
We move now to the microscopic picture of GaAs growth. Any new layer of GaAs forms new Ga-As pairs. Ga and As atoms are shown in red and blue in Figure 4, while the positions leading to ZB and WZ configurations are indicated by squares and hexagons, respectively. The inclined dashed lines provide a guide to define the atomic positions that lead to the relevant crystalline structures in Figure 4B,C. In the A-polar case, As must displace and occupy the position of ordered Ga in the liquid, with which it may form a dumbbell. The Ga atoms have two choices, as indicated in Figure 4 by the square and hexagon. Depending on the position selected, ZB (ABCABC stacking) or a twin (ABA stacking) is formed. These two positions are not equivalent in terms of their first and second nearest neighbors configuration. In the ABC stacking, Ga is found at the middle of a projected hexagon, while in the ABAB it is at a vertex. Further calculations beyond the scope of this work should corroborate if there is a position that is more energetically favorable. Still, to form a new bilayer, As atoms should first displace the ordered Ga layer on top of the NW, and this should increase the formation barrier. This is consistent with the difficulty in synthesizing A-polar GaAs NWs, [5,21] suggesting that the slower process may help the growing layer achieve the more thermodynamically stable ABC stacking.
The associated low probability of twinning also reduces the probability of WZ formation, as both are closely connected. [38] Experimentally, the A-polar GaAs NWs exhibit mostly a ZB structure with an absence of transverse twins. [5,21] During the formation of a new bilayer for the B-polar surface, the Ga atoms adsorbed on top of the As-terminated surface can stay at their position. The incoming As atoms occupy empty positions in the second row. The fact that the growth can proceed without displacing Ga atoms is consistent with the observation of a more facile growth for this polarity. The higher growth rate can also partly explain the higher propensity to introduce stacking defects and polytypism in B-polar GaAs NWs. [5] These results could explain why in polar semiconductors, growth in a certain polarity is preferred and how polarity determines the tendency for polytypism.
In conclusion, we have demonstrated that interfaces between a solid GaAs NW and the liquid Ga droplet exhibit a 3D order that depends on the polarity of the solid and the temperature. This order can be explained by the electrostatic interactions within the liquid and the solid termination. Detailed simulations of the interface, enabled by the combination of first-principles calculations and machine-learning, provide an atomistic picture of the growth process from the liquid, the likelihood for growth in a certain polarity and defect formation. This work lays the foundations for further atomistic studies of crystal growth processes. While our results have been obtained for the Ga(l)-GaAs(s) system, we believe they are of a general character. Similar ordering features likely play an important role in the crystal growth of other materials and in more general solidification processes, raising new questions for groups observing NW growth by in situ electron microscopy. This work also opens analogous perspectives in other fields such as catalysis and (photo)electrochemistry.

Experimental Section
NW Growth: GaAs NWs were grown on GaAs (100) substrates as described by Zamani et al. [21] HF etching was used to remove the substrate native oxide. Silicon oxide layers were obtained by spin-coating a solution of hydrogen silsesquioxane and methyl isobutyl ketone. The spin-coated layers were annealed at 300 °C for 10 min. The oxide thicknesses used for the growth was around 3.3 nm. Two additional annealing steps at 150 and 300 °C, lasting 2 h each were performed inside the MBE chamber. The growth conditions for the NWs were: a V/III ratio of 1.9, equivalent 2D Ga-limited GaAs growth rate of 0.24 Å s −1 , at a substrate temperature of 640 °C.
HAADF-STEM Imaging: The STEM samples were prepared by scraping as-grown samples onto Cu TEM grids with amorphous holey carbon films. The STEM imaging was performed using a double aberration-corrected FEI Titan Themis 60-300, operated at 300 kV high tension (HT). To minimize both electron-beam induced damage to the interface and scan distortions from, e.g., sample drift, the HAADF-STEM images were acquired as 24-frame image series, which were then averaged while correcting for rigid and nonrigid displacements using the SmartAlign plug-in for DigitalMicrograph. [39] Consecutive frames were acquired at a 90° rotations using a low beam current (<40 pA). The electron probe convergence semiangle and HAADF detector inner collection semiangles were 20 and 50.5 mrad. Fischione photomultiplier tube was used as the HAADF detector. Pixel dwell time of 0.5 µs and pixel step size ≈6 pm were employed. The images were captured with 2k × 2k size.
EELS Hyperspectral Imaging: The data were acquired with the same microscope, HT and probe convergence semiangle, using a Gatan GIF Quantum ERS spectrometer with Ultrascan 2k × 2k detector. In order to avoid damaging the interface, the probe current was set to 75 pA and a scan step size of ≈0.9 Å was applied. The energy dispersion was set to 0.05 eV per ch, while the zero-loss peak FWHM was 1.05 eV. A spectrometer entrance aperture of 2.5 mm, corresponding to a collection semiangle of 47 mrad was employed. For the measurements on the A-polar NW discussed in Figure 2 and in the Supporting Information, per spectrum pixel time was 0.25 ms, while this value was 0.5 ms for the measurement on B-polar NW reported in the Supporting Information. For all the EELS data, detector was binned by 130× on nondispersive axis and it was operating in high speed mode. "High quality dark correction" was applied. Processing the data started with removing spectral spikes from cosmic rays, centering the zero-loss peak to compensate for HT fluctuations and removing it, removing plural scattering and the plasmon peak replica. Then, second-order ICA was performed using the DigitalMicrograph plugin described by Lucas et al., [24] details of which are presented in the Supporting Information.
Molecular Dynamics: To investigate the ordering at the interface at the two (111) surfaces, an orthorhombic supercell composed of a central solid GaAs section (144 atoms, corresponding to 6 layers of 24 atoms), in contact with liquid Ga (192 atoms) on both of its surfaces, totaling 336 atoms, was used.
The initial lattice parameter for the solid part was set to the one obtained from DFT calculations, whereas the initial density of liquid Ga was set to that obtained with independent simulations in a smaller (96 atoms) box at the objective temperature. All the simulations were run using i-PI [36] in combination with large-scale atomic/molecular massively parallel simulator, [34] and n2p2 [31] to evaluate the NNP. First, the system was equilibrated in the NσT ensemble, allowing the cell degrees of freedom to change independently. After equilibration, production simulations were run in the NVT ensemble, using the average lattice parameters, at the temperatures indicated in the text.
The temperatures were controlled using a combination of a generalized Langevin [40] and stochastic velocity rescaling [41] thermostats. Pressures, where applicable, were constrained using an anisotropic barostat. [42] Simulations were run with a timestep of 4 fs (at 300 K) and 2 fs (at 900 K), for a total of 10 ns.
Neural Network and Ab Initio Calculations: A neural network potential of the Behler-Parrinello kind [30] was trained to reproduce the results obtained with DFT calculations. A total of 970 structures were used to fit the potential, with a 90/10 split between training and test. At the end of the training procedure, the average error in the test set was 2.4 meV per atom for energies and 120 meV Å −1 for forces.
The DFT calculations were carried out with the plane-wave code Quantum-Espresso, [32] using the Perdew-Burke-Ernzerhof exchangecorrelation functional, [43] and ultrasoft pseudopotentials [33] from the SSSP library (version 0.7). [44] The wavefunction energy cutoff was set to 50 Ry for all calculations, whereas the kinetic energy cutoff was set to 400 Ry.
Energy and forces for the structures used as training points were computed with a converged (<1 meV per atom absolute convergence) 3 × 3 × 1 Monkhorst-Pack k-point grid. [45]

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.