Universality and dependence on initial conditions in the class of the nonlinear molecular beam epitaxy equation
We report extensive numerical simulations of growth models belonging to the nonlinear molecular beam epitaxy (nMBE) class, on flat (fixed-size) and expanding substrates (ES). In both and , we find that growth regime height distributions (HDs), and spatial and temporal covariances are universal, but are dependent on the initial conditions, while the critical exponents are the same for flat and ES systems. Thus, the nMBE class does split into subclasses, as does the Kardar-Parisi-Zhang (KPZ) class. Applying the “KPZ ansatz” to nMBE models, we estimate the cumulants of the HDs. Spatial covariance for the flat subclass is hallmarked by a minimum, which is not present in the ES one. Temporal correlations are shown to decay following well-known conjectures.
Scaling invariance and universality, two pillars of the theory of phase transitions and critical phenomena, have been also very important in the study of nonequilibrium systems Ódor (2008). One of the most prominent examples is the dynamics of growing interfaces, whose width increases in time as (while the correlation length parallel to the substrate scales as ), and with the system size as (when ). A set of exponents , and defines an universality class and, interestingly, only a few classes exist, which are determined by some fundamental symmetries Barabasi and Stanley (1995). For example, interfaces evolving under tension and growing in the direction of its local normal are expected to belong to the Kardar-Parisi-Zhang (KPZ) class, being described at a coarse-grained level by the KPZ equation Kardar et al. (1986)
The Edwards-Wilkinson (EW) Edwards and Wilkinson (1982) equation (class) is given by . On the other hand, when the growth is dominated by the surface diffusion of adatoms, as is the case in molecular beam epitaxy (MBE), it is expected to fall into the nonlinear MBE (nMBE) class, associated with the equation by Villain Villain (1991) and Lai and Das Sarma Lai and Sarma (1991)
or in its linear counterpart (with ). In all these growth equations, is the height at substrate position and time ; , and , with , are constants and is a white-noise, with and variance Barabasi and Stanley (1995).
Recent theoretical Sasamoto and Spohn (2010); *Amir; *Calabrese; *Imamura, experimental Takeuchi and Sano (2010); *TakeuchiSP and numerical [Forrecentsurveyofliteraturesee; e.g.; ]healytake2015 works on KPZ systems have changed our view of KPZ universality by demonstrating that this class splits into subclasses depending on initial conditions (ICs), or surface geometry. More specifically, while the scaling exponents (, and ) are the same for KPZ growth starting from a flat substrate (flat IC/geometry) or from a seed - so that the active growing zone expands in time - (usually called curved or droplet geometry), the (1-point) height distributions (HDs) and (2-point) spatial and temporal correlators are different, but universal in each IC/geometry. In , the height fluctuations are given by Tracy-Widom distributions Tracy and Widom (1994), and spatial covariances are associated with Airy processes Prähofer and Spohn (2002); *Sasa2005; *Borodin.etal-CPAM2008. In higher dimensions, universality and IC dependence of KPZ HDs have been demonstrated numerically Halpin-Healy (2012); *Oliveira13; Alves et al. (2014); Halpin-Healy and Takeuchi (2015) and confirmed experimentally for the () flat subclass Almeida et al. (2014); *healy2014; *Almeida15.
Despite the importance of the nMBE class - since MBE is the main technique for thin film deposition - basically nothing is known about universality of (growth regime) HDs and IC dependence in these systems. In order to decrease this abyss between KPZ and nMBE classes, in this work we present a detailed numerical analysis of nMBE models studied on flat substrates of fixed-size (flat IC) and enlarging sizes (ES IC, which mimic the curved geometry Carrasco et al. (2014)). Results from large scale simulations, in and , demonstrate that universal and IC-dependent HDs and correlators also exist in nMBE growth. Beyond the obvious application of the flat subclass (for MBE growth on flat substrates), we note that the ES one might be appealing for deposition on textured substrates. A prominent example, which is very important for several applications Savin et al. (2015), is etched Si(100) surfaces, where inverted pyramid holes can be formed Wang et al. (2015) and, depending on the growth conditions and Si-adsorbate affinity, expanding surfaces might be observed while growth proceeds inside the holes.
To distill the universality of HDs, let us consider the so-called “KPZ ansatz” Prähofer and Spohn (2000)
where (the asymptotic growth velocity), (setting the amplitude of ) and (a stochastic correction) are non-universal (system-dependent) parameters, while is a random variable yielding the height fluctuations (which are universal in the KPZ class). A simple analysis of Eq. 2, considering periodic boundary conditions (PBC), shows that the mean height is always given by . Comparing this with Eq. 3, one sees that is equal to the deposition flux (), while the mean of the nMBE HDs is null (i. e., ), as well as are corrections in . This implies that the shift observed in the mean of KPZ HDs does not exist in the nMBE ones, since . The exponents and , and so , are exactly known from two-loop renormalization, where is a correction to the one-loop result Janssen (1997). Following a dimensional analysis of the nMBE equation, as done in Ref. Amar and Family (1992); *AmarFamily2 for 1-loop exponents (), we find here the scaling of the variance of HDs (for 2-loops) as
where is the correlation length and sets the roughness amplitude at the steady state regime (where ). In the growth regime, , so that . Comparing this with Eq. 3, one may identify and .
The standard discrete model in the nMBE class is the conserved restricted solid-on-solid (CRSOS) model Kim et al. (1994), where a (randomly deposited) particle aggregates in a site (i. e., ) if the restriction is satisfied for all nearest-neighbors (NN) . Otherwise, it is deposited at the nearest site of satisfying the restriction Kim et al. (1994). Theoretical calculations Park et al. (2001); *ParkKimPark2 for this model with (hereafter called CRSOS1), in , have demonstrated that it is described by the nMBE equation, in the hydrodynamic limit, with parameters , , and . Therefore, and for this model 111The one-loop () values are and , differing only slightly from the two-loop ones., which will be used as a benchmark in our analyses. Another classical nMBE model is the one from Das Sarma and Tamborenea (DT) Sarma and Tamborenea (1991), where the freshly (randomly) deposited particle, in a site , can move to its NN sites in order to increase the number of lateral neighbors. While the scaling of the original DT model is featured by strong corrections, a version with noise reduction, where an aggregation occurs at a given site only after deposition is attempted at that site, displays scaling exponents in good agreement with the nMBE class in Punyindu and Sarma (1998). Data for are presented in the following 222We have verified that simulations for yield the same asymptotic results.. Extensive simulations of the CRSOS model on substrates of fixed lateral sizes up to () and () were carried out for , (CRSOS2) and (CRSOS4). The DT model is investigated in for the same sizes. Furthermore, these models are also studied on enlarging substrates, using the method introduced by us in Ref. Carrasco et al. (2014). In this case, the growth starts on (flat) substrates of lateral size , which expand (in each dimension) at a constant rate by randomly duplicating columns. Here, one sets in and and in . In all models, PBCs are considered, and the deposition flux is defined as one particle per site per time unit, so that .
Effective growth exponents for CRSOS models with flat and ES ICs are compared in Fig. 1a. The convergence to the same asymptotic value demonstrates that the substrate enlargement preserves the roughness scaling properties, as expected Carrasco et al. (2014); Masoudi et al. (2012); *Escudero2. In general, one observes a very slow convergence of , which is still a bit smaller than the two-loop exponent even after very long times. This suggests the existence of additional corrections in the ansatz, so that , where , but and/or the covariance . Hence, . Indeed, by plotting versus time for the CRSOS1 model (see Fig. 1b), instead of a constant one finds a slightly decreasing behavior consistent with , so that if or otherwise. In any case, the extrapolation of to give us the variance of the HDs (for the CRSOS1 model). Higher order cumulants are determined in the same way, from , as shown in Fig. 1b, for and . The asymptotic cumulants for both ICs are summarized in Tab. 1. While in both cases, mild and considerable differences exist in and , respectively, demonstrating that the HDs are IC-dependent. In our analysis, we are assuming that is the same for flat and ES ICs Carrasco et al. (2014).
Since the parameter is known only for the CRSOS1 model, to confirm the universality of the HDs, we investigate the (adimensional) cumulant ratios: skewness and kurtosis . In the flat case, corrections and are found in and , respectively (see Fig. 1c). For ES, the exponents seem consistent with twice the ones for flat IC, but the extrapolated values are almost the same if we assume identical corrections. The asymptotic values of and for all investigated models, in the same dimension and IC, agree quite well, as shown in Tab. 2, confirming the universality of the HDs, as well as their IC dependence. Interestingly, is always very close to zero. Moreover, for flat ICs, is almost the same for and , so that these HDs have quite similar shapes, while in the ES case a decreasing is observed. This contrasts with the KPZ HDs, whose and are increasing functions of Alves et al. (2014), and it is possibly related to the fact that the nonlinearity in nMBE growth becomes irrelevant at its upper critical dimension Barabasi and Stanley (1995), where and are expected to vanish. We recall that the corresponding values of and for KPZ HDs (with flat and ES ICs in and ) fall into the ranges and Prähofer and Spohn (2000); Halpin-Healy (2012); Oliveira et al. (2013), being considerably larger than the ones in Tab. 2. Larger ratios ( and in and in ) have also been reported for the steady state HDs of the CRSOS model Reis (2004); *Tiagorug07, while a much smaller skewness () was recently found in a (one-loop) renormalization analysis of the nMBE equation in this regime Singha and Nandy (2016).
Although, without knowing , we cannot determine for all models, the product can be estimated, as done in Fig. 1d. Then, assuming the universality of the ’s in Tab. 1, one readily obtains (CRSOS4) and (DT, with ) in one dimension. The reliability of such estimates is confirmed by the nice data collapse shown in Fig. 2a, where the HDs , with , for different models are compared. We remark that these collapses confirm that is the same for fixed-size and enlarging substrates. Additional evidence of this is provided by the universality of the “cross-subclass” Halpin-Healy (2013) variance ratios , as shown in Tab. 2. To compare the HDs, we use the variable , which turns out to be simply , so that flat and ES ’s have variances 1 and , respectively. Again, a very good collapse is found (see Fig. 2b), which confirms that HDs are also universal and IC-dependent.
Now, we turn to the analysis of the spatial covariance
where , is a scaling function and is the same as defined above, in . Figures 3a and 3b show the rescaled for all investigated models in and , respectively. Interestingly, the curves for flat ICs cross the zero and have a minimum in the negative region, indicating the existence of a characteristic length in the interfaces, which is not present when the substrate expands. Since the ’s are not known for the CRSOS4 and DT models (in ), we determine them by making the minima of their curves (in flat case) to coincide with the one for the CRSOS1 model. This yields (CRSOS4) and (DT). The constant is obtained in the same way for curves, but shifting all minima to 1. Moreover, in this dimension one uses (obtained from simulations), instead of in the rescaling, so that and . The good collapse of rescaled curves confirms that and are universal, but . Hence, different processes exist for generating the flat and ES nMBE interfaces. It is worthy noting that the scaling functions [’s] for and are very similar (for a given IC), when appropriately rescaled. For instance, in the ES subclass, one finds approximately , for large , in both dimensions.
As an aside, from estimates of ’s and ’s in , one finds (CRSOS4) and (DT, ). Moreover, disregarding the (small) two-loop correction in , one obtains (CRSOS4) and (DT, ).
We also investigate the temporal covariance
Once more, the nice data collapse displayed in Fig. 4 demonstrates that universal IC-dependent scaling functions exist in the nMBE class. In , we have used , rather than in rescaling. Only data for ES ICs do not collapse well [see the inset of Fig. 4b], due to strong finite-time corrections , but when extrapolating the (rescaled) curves to , a very good agreement is obtained, as the main plot of Fig. 4b shows. A similar procedure has been employed to analyze the universality of in the KPZ class Carrasco et al. (2014). Substantially, in both dimensions, we find a power law decay , with exponents (flat) and (ES), in striking agreement with conjectures by Kallabis and Krug Kallabis and Krug (1999) and Singha Singha (2005), respectively.
In summary, we have demonstrated that 1-point height fluctuations in the nMBE class evolve, in the growth regime, according to the “KPZ ansatz” (Eq. 3) with universal and IC-dependent HDs. Moreover, 2-point spatial and temporal correlators are also IC-dependent. Therefore, the nMBE class splits into subclasses sharing the same critical exponents, similarly to KPZ systems. The absence of such splitting in HDs of linear classes, which are Gaussian for flat and ES ICs 333We have confirmed this by numerically integrating Eqs. 1 and 2 with in , suggests that this is a feature of nonlinear interfaces, possibly due to the lack of an up-down reflection symmetry in them. We claim that our findings will be very useful to confirm the universality class of growing systems, along the same lines of Refs. Takeuchi and Sano (2010); *TakeuchiSP; Almeida et al. (2014); *healy2014; *Almeida15; Brandt et al. (2015), especially because effective local roughness exponents close to the nMBE value () have been found in grained/mounded films Barabasi and Stanley (1995); Oliveira and Reis (2007b); *TiagoGraos2, but they can be a simple consequence of a geometric effect Oliveira and Reis (2007b); *TiagoGraos2. From a theoretical side, our results will certainly motivate and guide analytical works toward exact solutions of the nMBE equation and related discrete models.
Acknowledgements.We acknowledge support from CNPq, CAPES and FAPEMIG (Brazilian agencies), and thank S. O. Ferreira for helpful discussions. T.J.O. appreciates the kind hospitality of the group of Prof. James Evans at Iowa State University, where part of this work was done.
- Ódor (2008) G. Ódor, Universality in nonequilibrium lattice systems: theoretical foundations (World Scientific, Singapore, 2008).
- Barabasi and Stanley (1995) A.-L. Barabasi and H. E. Stanley, Fractal Concepts in Surface Growth (Cambridge University Press, Cambridge, England, 1995).
- Kardar et al. (1986) M. Kardar, G. Parisi, and Y.-C. Zhang, Phys. Rev. Lett. 56, 889 (1986).
- Edwards and Wilkinson (1982) S. F. Edwards and D. R. Wilkinson, Proc. R. Soc. London, Ser. A 381, 17 (1982).
- Villain (1991) J. Villain, J. Phys. I 1, 19 (1991).
- Lai and Sarma (1991) Z.-W. Lai and S. D. Sarma, Phys. Rev. Lett. 66, 2348 (1991).
- Sasamoto and Spohn (2010) T. Sasamoto and H. Spohn, Phys. Rev. Lett. 104, 230602 (2010).
- Amir et al. (2011) G. Amir, I. Corwin, and J. Quastel, Commun. Pure Appl. Math. 64, 466 (2011).
- Calabrese and Le Doussal (2011) P. Calabrese and P. Le Doussal, Phys. Rev. Lett. 106, 250603 (2011).
- Imamura and Sasamoto (2012) T. Imamura and T. Sasamoto, Phys. Rev. Lett. 108, 190603 (2012).
- Takeuchi and Sano (2010) K. A. Takeuchi and M. Sano, Phys. Rev. Lett. 104, 230601 (2010).
- Takeuchi et al. (2011) K. A. Takeuchi, M. Sano, T. Sasamoto, and H. Spohn, Sci. Rep. 1, 34 (2011).
- Halpin-Healy and Takeuchi (2015) T. Halpin-Healy and K. A. Takeuchi, J. Stat. Phys. 160, 794 (2015).
- Tracy and Widom (1994) C. Tracy and H. Widom, Commun. Math. Phys. 159, 151 (1994).
- Prähofer and Spohn (2002) M. Prähofer and H. Spohn, J. Stat. Phys. 108, 1071 (2002).
- Sasamoto (2005) T. Sasamoto, J. Phys. A: Math. Theor. 38, L549 (2005).
- Borodin et al. (2008) A. Borodin, P. L. Ferrari, and T. Sasamoto, Commun. Pure Appl. Math. 61, 1603 (2008).
- Halpin-Healy (2012) T. Halpin-Healy, Phys. Rev. Lett. 109, 170602 (2012).
- Oliveira et al. (2013) T. J. Oliveira, S. G. Alves, and S. C. Ferreira, Phys. Rev. E 87, 040102 (2013).
- Alves et al. (2014) S. G. Alves, T. J. Oliveira, and S. C. Ferreira, Phys. Rev. E 90, 020103(R) (2014).
- Almeida et al. (2014) R. A. L. Almeida, S. O. Ferreira, T. J. Oliveira, and F. D. A. A. Reis, Phys. Rev. B 89, 045309 (2014).
- Halpin-Healy and Palasantzas (2014) T. Halpin-Healy and G. Palasantzas, Europhys. Lett. 105, 50001 (2014).
- Almeida et al. (2015) R. A. L. Almeida, S. O. Ferreira, I. R. B. Ribeiro, and T. J. Oliveira, Europhys. Lett. 109, 46003 (2015).
- Carrasco et al. (2014) I. S. S. Carrasco, K. A. Takeuchi, S. C. Ferreira, and T. J. Oliveira, New J. Phys. 16, 123057 (2014).
- Savin et al. (2015) H. Savin, P. Repo, G. von Gastrow, P. Ortega, E. Calle, M. Garín, and R. Alcubilla, Nature Nanotech. 10, 624 (2015).
- Wang et al. (2015) Y. Wang, L. Yang, Y. Liu, Z. Mei, W. Chen, J. Li, H. Liang, A. Kuznetsov, and D. Xiaolong, Sci. Rep. 5, 10843 (2015).
- Prähofer and Spohn (2000) M. Prähofer and H. Spohn, Phys. Rev. Lett. 84, 4882 (2000).
- Janssen (1997) H. K. Janssen, Phys. Rev. Lett. 78, 1082 (1997).
- Amar and Family (1992) J. G. Amar and F. Family, Phys. Rev. A 45, R3373 (1992).
- Arnar and Family (1992) J. G. Arnar and F. Family, Phys. Rev. A 45, 5378 (1992).
- Kim et al. (1994) Y. Kim, D. K. Park, and J. M. Kim, J. Phys. A 27, L533 (1994).
- Park et al. (2001) S.-C. Park, D. Kim, and J.-M. Park, Phys. Rev. E 65, 015102(R) (2001).
- Park et al. (2002) S.-C. Park, J.-M. Park, and D. Kim, Phys. Rev. E 65, 036108 (2002).
- (34) The one-loop () values are and , differing only slightly from the two-loop ones.
- Sarma and Tamborenea (1991) S. D. Sarma and P. Tamborenea, Phys. Rev. Lett. 66, 325 (1991).
- Punyindu and Sarma (1998) P. Punyindu and S. D. Sarma, Phys. Rev. E 57, R4863 (1998).
- (37) We have verified that simulations for yield the same asymptotic results.
- Masoudi et al. (2012) A. A. Masoudi, S. Hosseinabadi, J. Davoudi, M. Khorrami, and M. Kohandel, J. Stast. Mech. 2012, L02001 (2012).
- Escudero (2009) C. Escudero, J. Stat. Mech. 2009, P07020 (2009).
- Reis (2004) F. D. A. A. Reis, Phys. Rev. E 70, 031607 (2004).
- Oliveira and Reis (2007a) T. J. Oliveira and F. D. A. A. Reis, Phys. Rev. E 76, 061601 (2007a).
- Singha and Nandy (2016) T. Singha and M. K. Nandy, J. Stat. Mech. 2016, 023205 (2016).
- Halpin-Healy (2013) T. Halpin-Healy, Phys. Rev. E 88, 042118 (2013).
- Kallabis and Krug (1999) H. Kallabis and J. Krug, Europhys. Lett. 45, 20 (1999).
- Singha (2005) S. B. Singha, J. Stat. Mech. 2005, P08006 (2005).
- (46) We have confirmed this by numerically integrating Eqs. 1 and 2 with in .
- Brandt et al. (2015) I. S. Brandt, V. C. Zoldan, V. Stenger, C. C. Plá Cid, A. A. Pasa, T. J. Oliveira, and F. D. A. A. Reis, J. Appl. Phys. 118, 145303 (2015).
- Oliveira and Reis (2007b) T. J. Oliveira and F. D. A. A. Reis, J. Appl. Phys. 101, 063507 (2007b).
- Oliveira and Reis (2011) T. J. Oliveira and F. D. A. A. Reis, Phys. Rev. E 83, 041608 (2011).