The compositional homogeneity of the metal particle during vapor–liquid–solid growth of nanowires

The vapor–liquid–solid (VLS) mechanism is probably the most versatile method to fabricate semiconductor nanowires and several investigations assume a compositionally homogeneous catalyst particle. In this investigation we address the compositional homogeneity of the catalyst particle during growth of nanowires. Using diffusion calculations, we show that the particle is indeed homogeneous during VLS growth, but can have a strong concentration gradient during vapor–solid–solid growth, that is, growth with a solid particle. We also show that the response to a concentration change is extremely fast, meaning that if the concentration at the surface of the particle changes, the entire particle reaches this new concentration effectively instantaneously.

Scientific RepoRtS | (2020) 10:11041 | https://doi.org/10.1038/s41598-020-67618-x www.nature.com/scientificreports/ thus sets the rate that atoms from the liquid incorporate into the solid. The measure is a dimensionless number, χ , that we define as where c max and c min are the maximum and minimum concentrations of nanowire-constituting atoms in the particle, respectively. From this definition it is clear that 0 ≤ χ ≤ 1 , where χ = 0 indicates a compositionally homogeneous particle and χ = 1 accounts for the maximum attainable concentration gradient through the particle. In order to find c max and c min we aim to solve the steady state diffusion equation, with the geometry and the boundary conditions defined in the caption of Fig. 1. In Eq. (2), ∇ 2 is the Laplacian and D and c are the diffusivity and the concentration of the nanowire species in the liquid particle, respectively. Thus, we consider the following boundary value problem with Ŵ = G/(D�) , where G is the axial growth rate of the nanowire and is the molecular volume of the nanowire species in the solid phase. Since c is, per definition, a harmonic function 11 , its maximum and minimum values are attained on the boundary of the domain. Using the symmetry of the domain it is easy to see that c max = c S and c min = c(0) , the latter being the concentration at the origin (the center of the growth interface), which we from now on denote c 0 . The homogeneity coefficient defined in Eq. (1) becomes Here we note that both of the concentrations in Eq. (6) are excess concentrations, meaning that the true concentrations are c S + c eq at the surface and c 0 + c eq is the minimum interface concentration, where c eq is the equilibrium solubility of the species in the metal particle. Most importantly, we also assume that the precursor material supply is efficient enough to ensure a constant surface concentration, c S , during some finite time interval, which is longer than the average time it takes for the species to diffuse through the particle.
The problem at hand is to determine an expression for χ in terms of R , Ŵ , and c S . However, with the geometry and boundary conditions according to Fig. 1, Eqs. (3)(4)(5) have no analytic solution using the standard technique of separation of variables. So, in order to determine an expression for c 0 , we first solve Eqs. (3)(4)(5) for simpler geometries, with the hope of finding a suitable generalization. These analytic calculations are outlined in the next section. Here we only mention the general result, which is that c 0 = c S − aŴR , where a is a constant with a value

Figure 1.
Schematic illustration of the geometry of the hemispherical catalyst particle with the coordinate system and boundary conditions indicated. The notation for the different parts of the domain are: B + (R) is the interior of the hemisphere, S + (R) is the curved surface of the hemisphere where c = c S , and U(R) is the growth interface where ∂c/∂z = Ŵ . In the center of the growth interface (the origin), c = c 0 . The subscript " + " refers to the upper hemispace.
Scientific RepoRtS | (2020) 10:11041 | https://doi.org/10.1038/s41598-020-67618-x www.nature.com/scientificreports/ that depends on the geometry of the catalyst particle. Inserting this into Eq. (6), we can express the homogeneity coefficient as a proportionality, In Table 1, we have listed the values of a for the different dimensionalities and geometries considered. For the hemisphere, we show that a = 1/2 in a separate section. The readers without interest in the details of our calculations can skip directly to the "Discussion of catalyst particle homogeneity" section without missing anything of essence to our conclusions.

Analytic calculations
Here we show the analytic solutions to Eqs. (3)(4)(5) for a few simple geometries and we start with the trivial one dimensional case: The solution to this boundary value problem is c(z) = c S − Ŵ(R − z) and the homogeneity measure in Eq. (7) is given by χ = ŴR/c S , that is, a = 1.
Next we solve Eqs. (3)(4)(5) in two dimensions, that is, for a semicircle instead of a hemisphere as in Fig. 1. With boundary conditions corresponding to those in Fig. 1, but with reduced dimension, we get the solution, where we have used planar polar coordinates so that r = √ x 2 + z 2 and θ = arctan (z/x) . The smallest interface concentration is also in this case found in the origin, c(0, θ ) = c S − 2ŴR/π , resulting in the homogeneity measure χ = 2ŴR/(π c S ) , that is, a = 2/π . While these two cases can serve as approximations to the diffusion profile with hemispherical boundary conditions ( Fig. 1), they are also interesting in their own right. The one-dimensional case corresponds to VLS film growth 12 and the two-dimensional case could be relevant for VLS growth of flat, fin-or sail-like structures 13 .
For the sake of completeness we write down the solution to the two dimensional boundary value problem for a rectangular boundary with sides of height R and a top edge of width 2R . At the sides and the top edge c = c S and at the bottom edge ∂c/∂z = Ŵ . With these boundary values, the solution to Eq. (2) is where n = (n + 1/2)π . The smallest interface concentration is given by Finally, we generalize the rectangular boundary and solve Eq. (2) in three dimensions with cylindrical boundary conditions. That is, we assume that the catalyst particle is shaped like a cylinder of radius R and height R . The concentration at the top circular area and at the side surface is c S and at the bottom circular area (the nanowire growth interface) we have the flux boundary condition ∂c/∂z = Ŵ , similar to the previous cases. With these boundary conditions, the solution to Eq. (2) is given by where J 0 and J 1 are the zeroth and first order Bessel functions of the first kind, respectively. The parameter α 0n is the nth zero of J 0 . The minimum concentration at the growth interface can be calculated as ≈ c S − 0.524ŴR , resulting in χ ≈ 0.524ŴR/c S ( a ≈ 0.524).

Scaling analysis
Since we have seen that the homogeneity index can be written in the form of Eq. (7) for all the investigated cases, we use a scaling approach to show that this proportionality is general and therefore valid for the more realistic geometry in Fig. 1. The problem at hand, Eqs. (3)(4)(5), contains a characteristic length scale, R , and a characteristic concentration, c S . Based on this we introduce the dimensionless variables x = R −1 x and u(x) = c −1 S c(Rx) . With these variable changes, Eqs. (3-5) transform into where we have dropped the overbars and derivatives are taken with respect to the dimensionless variables. The parameter γ is given by γ = RŴ/c S . In terms of the new variables, the homogeneity coefficient becomes We will now express χ in terms of γ , which is the only parameter left in the problem. Let u γ denote the solution to Eqs. (12)(13)(14) for a given scaled flux γ . Specifically, if γ = 0 then u 0 (x) = 1 in B + (1) . Suppose that we have found u 1 (x) , then we can construct u γ as which is clearly harmonic and satisfies the boundary conditions. Since it is a solution, it is the unique solution 14 . (15) we arrive at which proves Eq. (7), since γ = ŴR/c S , and the proportionality constant a can be identified as a = 1 − u 1 (0).
It is easy to verify that the proportionality in Eq. (17) holds for any reasonable domain if the part of the boundary where the flux γ is defined can be described using one characteristic length scale. However, a will be different for differently shaped domains and it can only be analytically calculated in certain special cases, as we have seen. In Fig. 2 we show numeric calculations of the proportionality constant, a , which serves as a shape factor, for different values of the contact angle, θ , describing the shape of the liquid metal particle. The value of a remains finite for the full range of θ , with a approaching zero in the small θ limit and a ≈ 0.63 as θ approaches 180°. We also see that a = 0.5 at θ = 90 • and in the next section we will show that this is an exact result. Figure 2. The proportionality constant a as a function of the contact angle θ as defined in the inset, calculated by finite element modelling and explicit calculation of the homogeneity index according to Eq. (6). Note that in the limit as θ → 180 • , the radius of the largest modelled particle is almost 100 times the wire radius R , but a nonetheless remains limited to 0.63.

Homogeneity coefficient for the hemisphere
We now turn to the case of the hemisphere. In order to compute a in this case, the new dependent variable w(x) = u 1 (x) − z is introduced. Notice that the function w is harmonic and is a solution to the boundary value problem The boundary condition on U(1) is homogeneous, which allows us to symmetrize the problem about the plane z = 0 . In other words, we seek a solution w to the following Dirichlet problem on the entire unit sphere, It is well-known that such a problem has a unique solution w , which is infinitely differentiable on B(1) and extends continuously onto the boundary S(1) . Furthermore, the reflected function w x, y, −z is also a solution to Eqs. (21)(22). Since the solution is unique it follows that w is an even function with respect to reflection about the plane z = 0 and therefore ∂w/∂z = 0 for z = 0 . If we define w to be the restriction of w to B + (1) , then w is the desired solution to Eqs. (18)(19)(20).
It is clear that u 1 (0) = w(0) =w(0) and the latter can be computed using the mean value theorem for harmonic functions 11 : where dσ is the surface-area measure on S(1) and |S(1)| = 4π (the surface-area of the unit sphere). We know that w = 1 − |z| on S(1) and inserting this in Eq. (23) and using, for instance, spherical coordinates it is a trivial task to compute w(0) = 1/2 . From this follows that a = 1 −w(0) = 1/2 for the hemisphere in Fig. 1.

Discussion of catalyst particle homogeneity
In this section, we will discuss the homogeneity coefficient, χ , and relate it to nanowire growth experiments, both VLS and vapor-solid-solid (VSS) growth. We will also relate χ to other dimensionless numbers.
Since the previously outlined analytically solvable cases can serve as approximations to the more realistic problem, the values of a are collected in Table 1, for easy comparison. Here we see that the one-dimensional approximation, Eq. (8), gives an error of a factor of two, which can still be acceptable, given the extreme simplicity of this approximation. The cylinder approximation overestimates χ by only 5% and the semicircle one with 27%. In the next section we will indeed use the one-dimensional approximation when calculating the time it takes to refill the particle.
Before we can discuss the homogeneity of the catalyst particle we need to estimate the parameters in Eq. (7). Typical growth temperatures for GaAs nanowires are 400-700 °C and the diffusivity of metal atoms in metal solvents are all in the range 10 -9 -10 -8 m 2 /s, depending on temperature and materials combination 15 . So, for our order of magnitude estimation it will suffice to set D Ga ≈ 5 × 10 -9 m 2 /s for the diffusivity of Ga in a liquid Au-Ga alloy. This agrees well with the DFT calculation of the Ga diffusivity in liquid Au, 3 × 10 -9 m 2 /s, presented in Ref. 16 . We also set the diffusivity of As in the Au-Ga liquid to D As ≈ 5 × 10 -9 m 2 /s. This is consistent with the approximation used by Roy and Chhabra 15 , where D AB is the diffusivity of solute A in solvent B, d A and d B are the respective atomic diameters, and D BB is the self-diffusivity of B. Combining Eq. (24) for the two cases: diffusion of Ga in Au-Ga and diffusion of As in Au-Ga, we can eliminate the atomic diameter and the self-diffusivity of Au. Then, since the van der Waals radii of Ga and As are almost equal ( r Ga = 1.87 Å and r As = 1.85 Å) 17 , the approximation D Ga ≈ D As is justified. Our estimation of D As also agrees with the diffusivities of about 4 × 10 -9 m 2 /s, extracted from liquid phase epitaxy experiments 18,19 but also measured by Dawson (see Ref. 18 ). Next we approximate the surface concentration as c S = x S ρ , where x S is the atomic fraction of the nanowire species, and ρ is the atomic density of liquid Au, ρ = 5.3 × 10 28 m -3 . The volume of a Ga-As pair in the solid is given by = 4.52 × 10 -29 m 3 . Now we can calculate χ for the growth of Au-catalyzed GaAs nanowires growing by the VLS mechanism. In order to find an upper limit for χ , we insert the highest axial growth rate that we are aware of, G ≈ 1 μm/s, observed for aerotaxy growth 10 . We use a large nanowire radius, R = 100 nm and we let the growth be As limited, with an As concentration at the surface of 1 at%, x S = 0.01. For this set of parameters we get χ = 4 × 10 -4 , that is, the catalyst particle is very close to being homogeneous also for this extremely high growth rate. This means that it is always safe to assume that the particle is compositionally homogeneous for VLS growth and that the axial Scientific RepoRtS | (2020) 10:11041 | https://doi.org/10.1038/s41598-020-67618-x www.nature.com/scientificreports/ growth rate is not limited by diffusion through the liquid catalyst particle. In Fig. 3 we show a numerical calculation of the concentrations in the catalyst particle for this set of parameters. Here it is clearly seen that the variation in concentration is extremely small and that the minimum concentration is at the center of the growth interface.
For the same set of parameters, we can in fact estimate what the growth rate would be if diffusion through the liquid were rate limiting. In this case we set χ = 1 and solve for the growth rate, resulting in the enormous growth rate G = 240 μm/s. It seems unlikely that crystalline nanowires can form with this high growth rate. On the other hand, diffusion through the particle can be rate limiting for VSS growth, where the diffusivity is orders of magnitude smaller. Persson et al. 20 investigated gold catalyzed growth of GaAs nanowires using chemical beam epitaxy, and concluded that the growth proceeded by the VSS mechanism and was limited by Ga diffusion through the solid catalyst particle. Indeed, using the data from this investigation, we estimate a homogeneity index of χ ≈ 1 , indicating a strong concentration gradient in the particle, consistent with diffusion controlled growth. On the other hand, Koryakin et al. 21 have investigated gold catalyzed growth of InAs nanowires at very low temperatures and also in this case, VSS growth was concluded. However, their relatively high diffusivity and their low growth rate led to a very low homogeneity index, χ ≈ 10 −5 . That is, the particles seem homogeneous and diffusion through the bulk of the particle cannot be rate limiting. This is not in conflict with the conclusions of Koryakin et al., who found that the growth is limited by diffusion of As along the nanowire-catalyst interface, which limits the nucleation rate of new layers 21 .
So far we have used a steady state model, meaning that the flux into the particle equals the flux out of the particle, which is equivalent to the growth rate. In certain experimental situations this steady state is broken so that the supply of material is smaller than the growth rate. This happens if the total amount of excess material in the seed particle is not large enough to complete one layer, that is if c S < 2 √ 3/ a 2 L R for a cylindrical nanowire growing in a {111}-direction with a hemispherical particle, where a L is the lattice constant. This corresponds to x S < 2 √ 3/ a 2 L Rρ , which evaluates to approximately 0.002 for R = 100 nm and 0.02 for R = 10 nm. Thus, for sufficiently thin wires, the excess amount of material in the catalyst particle will not suffice to complete a layer.
This leads to a situation where the layer nucleates and grows using both the (small) initial supply in the particle and the flux from the vapor phase. When this small initial supply is consumed, the concentration in the particle reaches a certain critical concentration related to a minimum in the Helmholtz free energy 22 , which we here estimate as c eq ( c 0 = 0 ), which in any case should be the lower limit of this concentration. After this the growth rate of the remaining part of the layer is limited by the flux from the vapor phase and will thus proceed slower than initially. For VLS growth, this effect is known as the "stopping effect" 23 . As we will show in the next section, the timescale for diffusion through the liquid is fast enough for steady state diffusion to set in almost immediately after the concentration is changed. Here we also mention that even in the case when the amount of material in the particle is sufficiently large so that there is no "stopping effect", there can still be some temporal variation of the concentration and thus of the supersaturation, which gives rise to anti-correlation of nucleation events, so called nucleation anti-bunching 24 .
Finally, we discuss χ in terms of other dimensionless numbers. The Damköhler number is defined as the reaction rate constant divided by the diffusion rate constant, Da = k r /k d 25 . If the nanowire growth rate can be described as a pseudo first-order process, the reaction rate constant can be written as k r = G/(�c 0 ) . The diffusion rate constant is given by k d = D/(aR) . This leads to Using Eq. (6) we substitute c 0 = c S (1 − χ) in Eq. (25) to arrive at Da = χ/(1 − χ) , or As the Damköhler number measures the growth rate in comparison to the diffusion rate, this can be a convenient route to estimate the catalyst particle homogeneity. We immediately see that if Da ≪ 1 , which it is for www.nature.com/scientificreports/ VLS growth, then χ ≈ Da. On the other hand, for materials systems where Da ≫ 1 , χ ≈ 1 and the growth is limited by diffusion through the particle. Another relevant dimensionless number is the mass transfer Biot number, Bi m . It is defined as the mass transfer rate at the interface divided by the mass transfer rate in the bulk 26 . Since the interface mass transfer in this system is identical to the growth rate, we have that Bi m = Da.

Diffusion time
In this section we calculate the time it takes to diffuse through the particle and refill it again to some small supersaturation after a stopping event, so that growth can proceed. To calculate this time, we use the one-dimensional approximation and we change the coordinate system as compared to Eq. (8), so that the particle surface is located at z = 0 and that the interface is located at z = R . The concentration is still the excess concentration (so that c = 0 means that the true concentration is c eq ). Since we will calculate the refill time, we solve the time dependent diffusion equation with a zero flux boundary condition at the interface, emulating a situation where growth has stopped, according to where ζ = z/R and r = R/ 2 √ Dt . In Fig. 4, we plot (c S − c 0 )/c S , where c 0 = c(R, r) , the concentration at the growth interface, as a function of r , our rescaled time variable. We see that as r decreases, or time increases, c 0 approaches c S . The inset shows the concentration profiles, that is (c S − c)/c S as a function of ζ (where ζ = 0 at the surface and ζ = 1 at the growth interface), for different values of r . Using Fig. 4, we can calculate the time it would take for the excess interface concentration to increase from 0 to just below c S , that is, the time it takes to make the particle homogeneous. If we for instance choose c 0 = 0.99c S as a homogeneity criterion (corresponding to χ = 0.01 ), then for D = 5 × 10 -9 m 2 /s and R = 100 nm, we get t = 4 µs. If we instead require c 0 = 0.999c S (or χ = 0.001 ), we get t = 6 µs for the same choice of parameters. It is interesting to note that this time is independent of the value of c S , depending only on the required ratio, c 0 /c S , and for relevant ratios for the homogeneity condition, it depends only weakly on this ratio.
Based on these time estimations we conclude that diffusion through the particle indeed occurs on a much faster time scale than it takes to grow a layer (0.5 ms/layer for Aerotaxy, typically much longer for MOVPE). This implies that after any change in surface concentration, the interface concentration approaches the surface concentration almost instantaneously, which in turn implies that the particle is always compositionally homogeneous, even if its concentration varies with time. Since we have used the one-dimensional approximation we expect the calculated times to be overestimated by a factor of two, which would make the homogeneity argument even stronger, and is in any case accurate enough for an order of magnitude estimation.

conclusions
We have investigated the compositional homogeneity of the metal catalyst particle during VLS (vapor-liquid-solid) growth of nanowires. We have introduced a homogeneity measure in the form of a dimensionless number and using steady state diffusion calculations we show that the catalyst particle is homogeneous during VLS growth but that there can be a large concentration gradient through the particle in the case of VSS growth, that is, growth with a solid particle. We have also performed time dependent calculations, which show that the response to a concentration change can be considered instantaneous on the time scale relevant for VLS growth. To conclude, the catalyst particle is homogeneous even if the concentration can vary with time, depending on the growth conditions and the size of the particle. The growth rate is not limited by diffusion through the particle for VLS growth of nanowires.

Methods
The analytical, steady state solutions to the diffusion equation, Eq. (2), were obtained using the Fourier method, that is, separation of variables. The time dependent solution to the zero flux boundary problem, Eq. (27), was obtained by summing and subtracting translated and reflected as well as only translated error function solutions, which also satisfy the diffusion equation, so that the derivative at the specified location becomes zero and that the surface concentration still has the desired value. This technique has been described by Crank 27 .
[erf(−ζ r + 2(2n + 1)r) − erf(ζ r + 2(2n + 1)r)]+ The numerical solution to the hemispherical problem was performed with finite element modelling (FEM), using COMSOL Multiphysics software 28 . Equations (3)(4)(5) were here solved numerically using a rotated twodimensional version of the geometry defined in Fig. 1, and a typical result for realistic parameter values is shown in Fig. 3 (solved in 3D for this illustrative plot). The parameter a for varying shapes of the seed particle (spherical segments with contact angle θ ) plotted in Fig. 2 was also calculated using FEM in a rotated 2D geometry, with the number of mesh elements at the particle-wire interface constrained to 100 regardless of particle size.  √ Dt . Note that c 0 is essentially zero for very short times and tends to c S in the long time limit. The inset shows the concentration profile for different values of r . Note that ζ = 0 at the surface of the particle and ζ = 1 at the growth interface.