PHYSICAL REVIEW E 102, 042802 (2020) Scaling of roughness and porosity in thin film deposition with mixed transport mechanisms and adsorption barriers Fábio D. A. Aarão Reis ,1,* Daniel O. Mallio ,1,† Jose Luis Galindo ,2 and Rafael Huertas 1 2,‡ Instituto de Física, Universidade Federal Fluminense, Avenida Litorânea s/n, 24210-340 Niterói Rio de Janeiro, Brazil 2 Departament of Optics, TexColImag Group, University of Granada, Granada 18071, Spain (Received 6 July 2020; accepted 13 October 2020; published 28 October 2020) Thin film deposition with particle transport mixing collimated and diffusive components and with barriers for adsorption are studied using numerical simulations and scaling approaches. Biased random walks on lattices are used to model the particle flux and the analogy with advective-diffusive transport is used to define a Peclet number P that represents the effect of the bias towards the substrate. An aggregation probability that relates the rates of adsorption and of the dominant transport mechanism plays the role of a Damkohler number D, where D  1 is set to describe moderate to low adsorption rates. Very porous deposits with sparse branches are obtained with P ≪ 1, whereas low porosity deposits with large height fluctuations at short scales are obtained with P ≫ 1. For P  1 in which the field bias is intense, an initial random deposition is followed by KardarParisi-Zhang (KPZ) roughening. As the transport is displaced from those limiting conditions, the interplay of the transport and adsorption mechanisms establishes a condition to produce films with the smoothest surfaces for a constant deposited mass: with low adsorption barriers, a balance of random and collimated flux is required, whereas for high barriers the smoothest surfaces are obtained with P ∼ D1/2 . For intense bias, the roughness is shown to be a power law of P/D, whose exponent depends on the growth exponent β of the KPZ class, and the porosity has a superuniversal scaling as (P/D)−1/3 . We also study a generalized ballistic deposition model with slippery particle aggregation that shows the universality of these relations in growth with dominant collimated flux, particle adsorption at any point of the deposit, and negligible adsorbate diffusion, in contrast with the models where aggregation is restricted to the outer surface. DOI: 10.1103/PhysRevE.102.042802 I. INTRODUCTION Thin porous films attract interest in several scientific areas due to their variety of applications, such as production of battery electrodes [1], solar cells [2], membranes for ion or molecule separation [3], and superhydrophobic coatings [4]. This motivates the study of several models for their deposition by vapor or solution techniques. Those models must account for the different transport mechanisms involved in the flux of materials to the deposits and for the different relaxations mechanisms taking place during or after the growth. Two widely studied lattice models that represent limiting cases of particle flux are ballistic deposition (BD) [5,6] and diffusion-limited deposition (DLD) [7]. BD considers a collimated flux and aggregation of particles at the first contact with the deposit; it was originally proposed for sedimentary rock formation [5]. DLD considers a purely diffusive particle flux and aggregation at the first contact with the deposit growing on a planar substrate [7,8]; it is an extension of the diffusionlimited aggregation model [9] in which the deposit grows from a point seed. The interplay of collimated and diffusive * fdaar@protonmail.com danielmallio@if.uff.br ‡ rhuertas@ugr.es † 2470-0045/2020/102(4)/042802(12) transport to the deposits was subsequently described by biased random walks in models with mechanisms similar to BD and DLD [10–21]. Extensions of BD and DLD were also proposed to account for post-deposition particle relaxation [13,22–28] or for finite rates of aggregation in the neighborhood of the deposit [7,29–34]. Possible applications include, but are not limited to, film deposition by electrochemical, electrospray, and nanocolloid drying techniques [13,26,27,34,35], in which transport is affected by fluctuating forces of the medium and by external fields. Recently, Ref. [18] introduced a deposition model with the biased random walk approach and with a constraint between particle adsorption and its motion near the deposit. It is termed potential adherence rule model (PARM) and assumes that a moving particle aggregates only if it attempts to hop onto a previously aggregated particle. This condition may be interpreted as an effect of an energy barrier which has to be overcome for the adsorption to occur and implies a finite adsorption rate, in contrast with the infinitely large rates of BD and DLD. An alternative interpretation of PARM is that the attachment of a reaction product occurs if the energy barrier for a surface reaction is overcome. Surface and bulk properties of two-dimensional deposits grown with PARM were formerly studied for broad ranges of the Peclet number (P) [18,20], which is defined in an analogy with advectivediffusive transport to relate the rates of collimated and random particle motion. This interpretation also shows that the model 042802-1 ©2020 American Physical Society REIS, MALLIO, GALINDO, AND HUERTAS PHYSICAL REVIEW E 102, 042802 (2020) TABLE I. Summary of the main features of the deposition models. Model Particle transport to the deposit Possible adsorption positions BD [5] Collimated DLD [7] RD [6] PARM [18] PARMA (this paper) Diffusive Collimated Mixed Mixed Top of column or nearest neighbor (NN) of top Any Top of column Any Any GBD [33] BD-RD [29] Collimated Collimated Any Top of column or NN of top can be applied to deposition of particles suspended in flowing solutions. The aim of this paper is to study simple thin film deposition models with barriers for the adsorption of incoming particles and the mix of collimated and diffusive flux. Our first step is to perform a numerical study of PARM in three dimensions. Next we introduce an extension of that model, which is termed PARMA (PARM with additional adsorption barrier), that represents aggregation rates smaller that those of PARM. In PARMA, the interplay of the adsorption and of the dominant transport mechanism can be described by a Damkohler number D, which is constrained to D  1 to represent moderate or slow adsorption. The relevance of the competition of diffusion and reactions to determine the morphology of films was actually stressed in recent works on electrochemical [36] and electrostatic spray [1] deposition. Here we also study a generalized BD model (GBD) introduced in Ref. [30] and further developed in Ref. [33], which considers collimated particle flux but finite probability of lateral aggregation at each contact with the deposit. Among the results, we show that films grown with a balance of random and biased transport attain a minimal surface roughness for a given deposited mass. The numerical study of the cases of dominant advection shows that the films have Kardar-Parisi-Zhang (KPZ) [37] roughening and shows power law relations for the roughness and for the (low) porosity as a function of P and D. Theoretical justifications for those results are then provided by scaling approaches, which show the conditions for production of the smoothest film surfaces and determine the exponents of the relations for the roughness and for the porosity with large P/D. This includes the observation of a superuniversal scaling of the porosity in terms of P/D in the general context of PARMA. Equivalent relations are extended to GBD. The rest of this paper is organized as follows. In Sec. II, we present the models studied here, define the dimensionless numbers P and D, the main quantitites to be measured, and present simulation details. In Sec. III, we present the numerical results. In Sec. IV, we present theoretical arguments to explain the condition of minimal roughness in PARM and PARMA and to explain the scaling of roughness and porosity in cases of dominant collimated transport, including a discussion of GBD. In Sec. V, we summarize our results and present our conclusions. Adsorption condition First contact First contact First contact Hop to occupied site Hop to occupied site and additional barrier Finite probability Finite probability II. MODELS AND METHODS The deposition models in dimension d = 3 are defined in the region z  0 of a simple cubic lattice. The edge of a site has length a, and each particle occupies a single lattice site. The initially flat substrate is located at z = 0 and all points with z > 0 are empty. The lateral size of the lattice is L and periodic boundary conditions are considered in the lateral directions x and y. The set of particles with the same (x, y) position is termed a column of the deposit, and the height variable h(x, y) is the position z of the topmost particle in that column. Deposition in d = 2 is defined in the region z  0 of the xz plane; other definitions are analogous to those in d = 3. The models are of limited mobility, i.e., when a new particle is deposited, all previously deposited particles are already at their final aggregation positions. In Table I, we summarize the main features of the models discussed in this paper. A. PARM PARM is defined as in Refs. [18,20,38]. Each new particle is released at a randomly chosen column in a height z much larger than all values of h(x, y) [h(x) in d = 2]. It executes a sequence of hops to NN sites with probabilities depending on the hop directions. In each hop attempt, the diffusive and the collimated components are chosen with probabilities κ and 1 − κ, respectively. If diffusion is chosen, the hop direction is then randomly chosen among the 2d NNs; this gives a probability κ/(2d ) of attempting to hop to each NN. If the collimated component is chosen, the particle attempts to hop to the site immediately below (direction −z). In any case, if the chosen site is empty, the particle moves to that site; otherwise, the particle permanently aggregates at its current position. Figure 1(a) illustrates these rules in d = 2 to facilitate visualization. The collimated and diffusive components may be interpreted as effects of a constant field and a fluctuating field, respectively. If growth occurs in a solution that flows in the z direction, the interpretation of a competition of advective and diffusive transport is also possible as suggested in the original works on PARM [18,20]. The aggregation in PARM is interpreted here as an effect of an adsorption barrier, which is overcome when the particle motion at the current site leads 042802-2 SCALING OF ROUGHNESS AND POROSITY IN THIN … (a) PHYSICAL REVIEW E 102, 042802 (2020) in the limiting cases κ = 0 and κ = 1. However, in the totally collimated case κ = 0, no lateral hop is allowed, so the adsorption always takes place at the top of the column of incidence (the probability q can only retard the adsorption after the incident particle reaches the top of that column). Thus, PARMA is also equivalent to RD in this condition. Here we described the additional barrier for aggregation in PARMA by a probability. However, for a specific application, q may be related to an energy barrier and the temperature, e.g., by an Ahrrenius law. (b) PARM GBD FIG. 1. (a) Illustration of the deposition of two particles in PARM. Each incident particle (red square) is released far above the deposit, executes a biased random motion among NN sites (solid line), and is adsorbed at the position (dashed square) from which it attempts to hop onto a previously deposited particle (blue squares). The last hop attempt is indicated by an arrow and dots (green) indicate the positions of the trajectory in which the particle has an occupied NN but does not attempt to hop onto that NN. (b) Illustration of the deposition of two particles in GBD. Each incident particle is released above the deposit, has a vertical trajectory, and aggregates laterally with probability p; each aggregation attempt is represented by a green line connecting with a NN deposited particle. The particle at the left aggregates at the fourth attempt, and the particle at the right aggregates when it reaches the top of the column. to a direct collision with a previously aggregated particle. Within this interpretation, we are considering that the energy fluctuation of the incident particle is sufficient to overcome that barrier. If the flux is only diffusive (κ = 1), PARM is similar to DLD, but not equivalent. In PARM, the incident particle may aggregate at the first contact with the deposit with a finite probability because the adsorption depends on choice of hop direction, but in DLD there is an infinitely large adsorption rate at the neighborhood of the deposit. In the case of purely collimated flux (κ = 0), no lateral hop is possible during the particle trajectory, so the incoming particle always aggregates at the top of the column of incidence. This produces compact deposits with uncorrelated local heights so that PARM becomes equivalent to the random deposition (RD) model [6]. C. Peclet and Damkohler numbers In processes with advective and diffusive particle transports, the Peclet number P is defined as the ratio between the respective rates of those processes. The concept may be extended to processes with electrophoresis or thermophoresis towards a hard wall as discussed in Ref. [39]. Following this reasoning, Ref. [10] defined the Peclet number of the biased random walk model considered here as the ratio of the characteristic rates of collimated and diffusive transport within the distance a. First, it was shown that the diffusion coefficient is D = κ (2 − κ )a2 /(2 dτ ) and that the velocity towards the substrate is v = (1 − κ )a/τ , where τ is the time of a single hop. These relations give P= (1) The purely advective case κ = 0 corresponds to P → ∞, and the purely diffusive case κ = 1 corresponds to P = 0. The Damkohler number D describes the competition of transport and reaction (adsorption) in PARMA [40]. For a particle in the neighborhood of the deposit, the probability of adsorption to a given NN is q/(2d ) (2d is the number of NNs). The residence time of the incident particle in a site is τ , which gives a reaction rate q/(2 dτ ) (probability per unit time). The relevant transport rate depends on the dominant mechanism, which is controlled by the Peclet number. For P > 1, the flux is dominated by the bias, so the Damkohler number is defined as the ratio between the reaction rate and the rate of vertical motion in the length scale a, B. PARMA D= PARMA represents larger energy barriers for adsorption, which may not be overcome with the average kinetic energy of the hopping particle. The rules for the particle flux remain the same as in PARM, i.e., with a mixture of collimated and diffusive components. However, when the mobile particle attempts to hop to an occupied site, it aggregates at the current position with probability q < 1, otherwise it continues to move from the current position with the same transport rules. For instance, consider the hopping particles in Fig. 1(a). In PARMA, the last indicated hop attempts lead to aggregation at the marked position (dashed square) with probability q, but with probability 1 − q the particle continues to hop starting from the marked position. The same probability q is applicable for adsorption in further attempts to hop to an occupied site. PARM corresponds to the case q = 1. PARMA has differences from BD and DLD because adsorption may not occur at the first contact with the deposit 2d (1 − κ ) v/a . = D/a2 κ (2 − κ ) q/(2dτ ) q = . v/a 2d (1 − κ ) (2) The condition P > 1 in d = 3 implies κ < 0.838, so Eq. (2) always gives D ∼ q, i.e., other numerical factors relating D and q are on the order of 1. For P < 1, in a regime dominated by diffusion, D is defined as the ratio between the reaction rate and the diffusion rate in the length scale a, D= q/(2 dτ ) q . = 2 D/a κ (2 − κ ) (3) For d = 3, 0.838  κ  1 in this regime, so Eq. (3) also gives D ∼ q for any value of κ. For these reasons, the order of magnitude of the Damkohler number is always that of the adsorption probability q. In PARM, we have D ∼ 1, which means that adsorption and transport are balanced. In PARMA with q ≪ 1, adsorption is slow compared to transport. These models exclude cases of very fast adsorption or of a very fast sequence of reaction and 042802-3 REIS, MALLIO, GALINDO, AND HUERTAS PHYSICAL REVIEW E 102, 042802 (2020) product adsorption; the excluded cases have transport-limited conditions represented, e.g., by the standard adherence rule model studied in Refs. [18,20,38], whose limiting behaviors are exactly equal to the DLD and BD models. D. GBD GBD is defined as in Refs. [30,33,41]. Particles are sequentially released at randomly chosen columns and heights larger than all h(x, y)’s. They execute fully biased walks towards the deposit, i.e., walks with hops only in the −z direction. After each hop, if the particle has n NN occupied sites, it can aggregate at the current position with probability np, i.e., with a probability p for each NN; with probability 1 − np, the particle hops to the site immediately below it. If the particle reaches the top of the column of incidence [position z = h(x, y) + a], it permanently aggregates there because lateral hops are not allowed. Thus, pore formation is possible only for p > 0. Figure 1(b) illustrates the GBD rules in d = 2 to facilitate the visualization. For p = 1, the model is equivalent to BD in which aggregation occurs at the first contact with a NN occupied site. For p = 0, it is equivalent to RD because lateral aggregation is not possible. A probability of aggregation 0 < p < 1 may also be interpreted as an effect of an adsorption barrier, but GBD differs from PARM and PARMA with diffusive transport (κ > 0). Note that GBD is also not equivalent to the competitive model in which the rule for particle aggregation is chosen between that of BD (probability p), i.e., aggregation at the first contact and that of RD (probability 1 − p), i.e., aggregation at the top of the column of incidence [29,42]. GBD is a type of slippery aggregation model [41] in which the incident particle can be adsorbed at several points whereas it slips in contact with deposited NNs. E. Quantities of interest The growth in PARM and PARMA occurs in timescales that depend on the type of flux. However, since we are mainly interested in the structural properties of the deposits, we analyze their evolution in terms of an effective number of deposited monolayers, which is denoted as T . This number is the ratio between the number of deposited particles N and the number of columns (L/a)d−1 , T = N . (L/a)d−1 (4) T may also be interpreted as the thickness of a compact film with the same mass. The thickness of a porous deposit H is always larger than that of a compact film with the same mass T . The dimensionless thickness H is defined as the average of the heights of the outer surface,   1  H≡ h(x, y) . (5) (L/a)2 x,y The main quantity that characterizes the fluctuations of the outer surface is the roughness, whose dimensionless value is defined as W ≡  1  [h(x, y) − H]2 (L/a)2 x,y 1/2  . (6) These definitions apply to deposits in d = 3; in d = 2, sums are over L/a columns in the x direction. The main quantity that characterizes the inner part of a porous deposit is the porosity, which is defined as the volume fraction of the pore space. The deposit is limited by the substrate, the lateral sides, and the outer surface, so its volume is H (L/a)2 and the total porosity T is complementary to the volume fraction of the solid, T ≡ 1 − T (L/a)2 T =1− . 2 H H (L/a) (7) As the deposits grow, we observe a slightly nonlinear evolution of H with T , which is related to small variations in the porosity of different layers of the film [43]. In a certain interval of deposition time, corresponding to the number of monolayers in the interval T1  T  T2 , the local porosity is given by (T1 , T2 ) = 1 − T2 − T1 . H (T2 ) − H (T1 ) (8) In the three-dimensional films studied here, the local porosity has an approximately constant value if T1 and T2 are both close to the maximal thickness. In the study of transport of fluids or solutes in porous deposits, it is important to distinguish the total porosity and the effective porosity, the latter comprising only the pores that participate in the transport between two specified frontiers. This is not of direct interest for the present paper, so we only analyze the total porosity. Moreover, we are particularly interested in deposits produced in regimes of dominant advective transport, which have low porosity and mostly isolated pores. F. Simulation details Simulations of PARM in d = 3 were performed in L = 128a, L = 256a, and L = 512a, in most cases up to T = 103 . The simulation in three different sizes is necessary to avoid finite size effects, which are particularly strong for small P. Our numerical analysis will be based on data from L = 256a, which we ensure that are not contaminated by such effects. Deposition was simulated from P ∼ 1 (balanced diffusion and advection) to P ∼ 103 (dominant advection). For each P, average quantities were calculated from 102 different deposits. PARMA was simulated with adsorption probability q = 0.25 in d = 3 and sizes L = 128a, L = 256a, and L = 512a. Average quantities were calculated from 102 different deposits for each P in the range of P ∼ 1 to P ∼ 104 . Several results for PARM in d = 2 were formerly presented in Refs. [18,20]. Here we use part of those data to investigate the dependence of the roughness and of the porosity on P. The simulations were performed in L = 400a and sufficiently large heights to allow the deposition of T = 2.5 × 103 monolayers. For each P, average quantities were calculated from 25 different deposits. 042802-4 SCALING OF ROUGHNESS AND POROSITY IN THIN … PHYSICAL REVIEW E 102, 042802 (2020) GBD was simulated in d = 3 with L = 1024a, which is sufficiently large to avoid finite size effects. For 0.005  p  0.1, average quantities were obtained from 103 different deposits up to T = 104 . III. NUMERICAL RESULTS A. General features of PARM in d = 3 Figure 2 shows top portions of deposits grown with six different values of P. Figures 3(a)–3(f) show cross sections of the same deposits. FIG. 3. Top portions of deposits with T = 200 layers grown with PARM with the indicated values of P. FIG. 2. Top portions of deposits grown with PARM with the indicated values of P and constant T = 200. For P < 1, the deposits have many branches and high porosity; this is the reason why several isolated particles or clusters are observed in the panels of P = 0.06 and 0.301 in Figs. 2 and 3. The pore system is clearly connected along the vertical direction. As P increases, the porosity decreases and a large fraction of the pores are isolated. For large P, the deposits are much more compact; their pores are elongated in the vertical direction, which can be observed by close inspection of the images of the films obtained with P = 298 in Figs. 2 and 3. The global height fluctuations in the films grown in limiting cases, i.e., with P ≪ 1 and P ≫ 1, are apparently larger than the fluctuations in the films grown with intermediate values of P, i.e., P ∼ 1. The decrease in the porosity with increasing P, for a constant number of deposited layers, corresponds to an increase in the film density; see Eq. (7). This was experimentally observed, for instance, in recent electrospray deposition of nanoparticulated TiO2 films [44]. 042802-5 REIS, MALLIO, GALINDO, AND HUERTAS PHYSICAL REVIEW E 102, 042802 (2020) FIG. 5. (a) Surface roughness at T = 103 and (b) average porosity between T1 = 800 and T2 = 1000 as a function of P for films grown with PARM in d = 3. Dashed lines are linear fits of the data for P > 102 . The uncertainties in both quantities are smaller than the size of the symbols. FIG. 4. Surface roughness of PARM in d = 3 as a function of the number of deposited monolayers for: (a) low and intermediate values of P and (b) intermediate and large P’s. The uncertainties are smaller than the size of the symbols. B. Roughness and porosity of PARM in d = 3 Figure 4(a) shows the roughness evolution in cases with balanced collimated and diffusive flux (P ∼ 1) in lattice size L = 256a. For the smallest P, the ranges of T in the plots are small because finite size effects appear for larger thicknesses. For a constant and large number of monolayers T , i.e., constant number of deposited particles, it is clear that W has a nonmonotonic variation with P with minimal values obtained for P ∼ 6. Figure 4(b) shows the roughness evolution in cases of highly biased transport (P ≫ 1) in the same lattice size. For a quantitative analysis of the roughness evolution, we recall that kinetic roughening theory predicts β W ∼T , (9) with β < 0.5 for systems with stable and correlated growth and β = 0.5 for RD [6,45]. Effective values of the exponent β may be extracted from the slopes of log10 W × log10 T plots like those in Figs. 4(a) and 4(b). The growth with P = 0 is similar to that of DLD, which is unstable, i.e., W increases faster than in Eq. (9). However, the partially collimated flux for any P ≪ 1 is expected to stabilize the growth [12]. This may lead to long crossovers and explain the variations in the slopes of the plots for different P in Fig. 4(a). The effective slope slowly increases with P in this regime, which parallels the recent observation of larger effective growth exponents β for larger flow rates observed in radial imbibition of a porous medium [46]. Figure 4(b) shows the roughness evolution for intermediate and large P. As P increases, we observe a longer regime with effectively RD, so the initial slopes of the plots are close to 0.5. For large T , the growth crosses over to a correlated regime. The plots in Fig. 4(b) have slopes close to the value β = 0.242 of KPZ scaling [47,48] for all P  10’s. KPZ scaling is theoretically predicted for the growth of porous deposits with stable surface roughening [12,13,49]. However, in models where the height fluctuations are large, deviations are frequently observed in numerical simulations; this is the case, for instance, of BD in d = 3 [50]. For this reason, it is somewhat surprising that the KPZ exponents are obtained here for P ≫ 1 with good accuracy at relatively short deposition times. For P  1, theoretical arguments developed for related systems [12,13] may be extended to PARM and suggest that KPZ is the asymptotic class of roughening for all P’s. The deviations observed here, however, show that much longer simulations and more powerful methods of analysis (e.g., as in Ref. [51]) would be necessary for a numerical confirmation of this feature. Our next step is to analyze the effect of P on the roughness and on the porosity of films with a constant number of deposited layers. For T ∼ 103 , Fig. 4(b) shows approximately parallel curves for different values of P. It suggests a power law roughness-flux relation for constant T as W ∼ Pγ , (10) with a positive exponent γ . Figure 5(a) shows the roughness measured at T = 103 as a function of P, which confirms the power law increase and gives the estimate γ = 0.160 ± 0.001 for the largest values of P. Note that the corresponding 042802-6 SCALING OF ROUGHNESS AND POROSITY IN THIN … PHYSICAL REVIEW E 102, 042802 (2020) FIG. 6. (a) Surface roughness at T = 2.5 × 103 and (b) total porosity as a function of P for films grown with PARM in d = 2. Dashed lines are linear fits of the data for P > 102 . The uncertainty in the porosity is smaller than the size of the symbols. deposits have the same number of particles, but their average thicknesses H are different. Figure 5(b) shows the local porosity measured between T = 800 and T = 103 for several values of P. It monotonically decreases with P, which is explained by the decrease in the frequency of lateral aggregations. For larger P, Fig. 5(b) suggests a porosity-flux relation as  ∼ P −λ , (11) with a positive exponent λ. The linear fit in Fig. 5(b) gives λ = 0.337 ± 0.008. C. Roughness and porosity of PARM in d = 2 The evolution of the roughness and of the porosity with the number of deposited particles N was discussed in Refs. [18,20]. For fixed N, the nonmonotonic variation of W with P was observed, similar to that in d = 3 [Fig. 4(a)]; the minimum value of W was obtained for P ∼ 2. Here we extend the analysis of the data shown in Ref. [18] to determine the effect of P on the surface and bulk structure in conditions of dominant collimated flux. In Fig. 6(a), we show the roughness of the thickest deposits with T = 2.5 × 103 monolayers as a function of P in the range of 101 < P < 103 . The linear fit of those data shows consistency with Eq. (10) with γ = 0.111 ± 0.004. In Fig. 6(b), we show the total porosity T in the same range of P. It it consistent with Eq. (11) with λ = 0.343 ± 0.007. D. Roughness and porosity of PARMA in d = 3 Figure 7(a) shows the roughness evolution with adsorption probability q = 0.25 for 1  P  10. The plot is restricted to T < 300 because finite size effects are observed at larger thicknesses for the smallest P. For constant T , the minimal roughness is obtained for P ∼ 3, which is smaller than the value of P that provides the minimal roughness in PARM (q = 1) [Fig. 4(a)]. Figure 7(b) shows the roughness evolution for q = 0.25 and P ≫ 1, i.e., strongly collimated flux. Comparison with Fig. 4(b) shows that the initial uncorrelated regime (RD) lasts longer than that of PARM. For instance, for P = 598, this regime ends at T < 102 in PARM, but it clearly ends at T > 102 in PARMA. For large thicknesses, there is also evidence of KPZ scaling; in the largest P, this observation is FIG. 7. Surface roughness of PARMA in d = 3 for q = 0.25 and L = 256a as a function of the number of deposited monolayers: (a) low and intermediate values of P and (b) intermediate and large P’s. The uncertainties are smaller than the size of the symbols. possible because the simulations of PARMA were extended to T = 5 × 103 . For fixed T , the roughness monotonically increases with P in this regime. Figure 8(a) shows the roughness measured at T = 5 × 103 as a function of P. The linear fit of the data for P > 3 × 102 shows agreement with Eq. (10) and gives γ = 0.190 ± 0.001. Figure 5(b) shows the local porosity of the same deposits measured between T = 4 × 103 and T = 5 × 103 . The linear fit for large P agrees with Eq. (11) with λ = 0.325 ± 0.009. E. GBD in three dimensions Here we focus on the cases of small p, in which adsorption with pore formation is not frequent. This regime parallels that of large P in the models with noncollimated flux. Figure 9(a) shows the evolution of the roughness in these conditions. The initial RD regime is also followed by a crossover and an asymptotic KPZ regime, which is confirmed by the effective slopes of the plots near β ≈ 0.242 [47,48]. For fixed T , as 042802-7 REIS, MALLIO, GALINDO, AND HUERTAS PHYSICAL REVIEW E 102, 042802 (2020) An interesting feature here is that similar values of the exponent λ of the porosity-flux relation were obtained in all models in d = 2 and d = 3. However, numerical estimates of the exponent γ of the roughness-flux relation have larger relative differences. IV. SCALING APPROACHES A. Smoothest surfaces in PARM FIG. 8. (a) Surface roughness at T = 5 × 103 and (b) local porosity between T = 4 × 103 and 5 × 103 as a function of P for films grown with PARMA in d = 3 with q = 0.25. Dashed lines are linear fits of the data for P > 102 . The uncertainties are smaller than the size of the symbols. p decreases, W increases because the correlations in the column heights are weaker. Simulations for larger values of p [not shown in Fig. 9(a)] show that the roughness has small variations for a constant thickness and p > 0.1. Figure 9(b) shows the roughness at T = 104 as a function of p. The linear fit suggests a power law relation between those quantities for constant T as W ∼ p−γ1 , (12) with γ1 = 0.135 ± 0.005. Figure 9(c) shows the total porosity T obtained in the deposition of T = 104 layers as a function of p. The linear fit also suggests a power law relation as T ∼ pλ1 , (13) with λ1 = 0.302 ± 0.003. FIG. 9. (a) Roughness as a function of the number of deposited layers in GBD in d = 3 for L = 1024a and several values of p. The dashed line has slope 0.242. (b) Roughness and (c) total porosity measured at T = 104 as a function of p. The dashed lines are linear fits of the data in these plots. When films with a constant number of deposited layers T are compared, there is a value of P in which the roughness is minimal. For applications of this type of model, it is important to understand the conditions for formation of films with the smoothest possible surfaces. We begin analyzing the case of dominant diffusion P ≪ 1 in which transport and aggregation follow rules similar to DLD. The growth in DLD is unstable because the highest branches of the deposit capture most of the incoming particles, shadowing the lowest points; the height fluctuations become very large [see, e.g., the case P ≪ 1 in Fig. 3(a)]. The aggregation in PARM is slightly different because it does not occur at the first contact with the deposit. For instance, in the trajectory at the right side of Fig. 1(a), the particle stood at three sites in which it had an aggregated NN, but it only aggregated at the fourth of those sites because it attempted to hop to an occupied NN. For each occupied NN that the mobile particle meets, there is a probability of aggregation 1/(2d ). Thus, PARM with κ = 1 (P = 0) is equivalent to DLD with a probability 1/(2d ) of aggregation to a NN occupied site. For nonzero P, a small bias can suppress the growth instability, i.e., it can lead to a slower increase in the roughness at long times as in Eq. (9). This parallels other systems in which instabilities are suppressed by biased fluxes [12,13]. For this reason, in the diffusion dominated regime P ≪ 1 and for constant T, W decreases as P increases; Fig. 4(b). Now consider the case of fully biased motion P → ∞ in which the particle flux is collimated. The aggregation always occurs at the top of the column of incidence, which corresponds to RD [6]. For finite P ≫ 1, the lateral aggregation is the only mechanism in PARM that leads to correlations between neighboring columns, and this lateral aggregation is associated with the diffusional transport. The correlations along the surface lead to roughness scaling as in Eq. (9) with β < 1/2. For large T , it explains the decrease in W with the decrease in P in this advection dominated regime. The arguments in the above paragraphs mean that, from one side, W decreases because the intensity of the bias increases, and from the other side, W decreases because diffusion becomes more relevant. This implies that the roughness reaches a minimal value when the diffusive and the collimated components of the flux are balanced, i.e., for P ∼ 1. (14) Those mechanisms are independent of the spatial dimension, which is consistent with the numerical results in d = 2 [18] and d = 3 (Sec. III B). 042802-8 SCALING OF ROUGHNESS AND POROSITY IN THIN … PHYSICAL REVIEW E 102, 042802 (2020) B. Smoothest surfaces in PARMA Here we extend the previous arguments to PARMA in which we also have to account for the (possibly low) adsorption probability q. Recall that D ∼ q in all cases. For P ≪ 1, the shadowing effect of the highest branches of the deposit is reduced when q < 1, which helps to control the instability. This is equivalent to an effectively larger P. When the incident particle is near the deposit, its vertical spread is proportional to the square root of the number of hops. This number is proportional to the number ∼1/q of contacts with the deposit before aggregation. Thus, the vertical spread is ∼1/q1/2 and the effective P is expected to increase by this factor. Now consider the case P ≫ 1. The probability q < 1 leads to longer RD regimes as observed in Sec. III D, which is also equivalent to increasing P. In the trajectory of the incident particle near the deposit, the number of contacts before aggregation is ∼1/q, so the spread in the particle position is ∼1/q1/2 and we expect that this factor increases the effective value of P. For these reasons, the condition of minimal roughness in PARMA is P/q1/2 ∼ 1. This is consistent with our numerical results: the value of P of minimal roughness decreases by a factor of ∼2 [from 5.63 in Fig. 4(a) to 2.86 in Fig. 7(a)] as q decreases by a factor of 4 (from 1 to 0.25). The condition can be written in terms of P and D as P ∼ D1/2 . (15) Note that numerical factors on the order of 1 were omitted in this relation and that such factors may be model dependent. It is important to recall that the instability control of DLD by a finite probability of aggregation (or by a finite reaction rate) is consistent with theoretical predictions [52] and was experimentally observed in a recent work on electrodeposition of Ni and NiW films [53]. Eq. (15) actually suggests that a very small D (very slow aggregation rate) may lead to small roughness even if P is small (diffusion-dominated flux). The experimental realization of the regime with highly collimated flux of PARM and PARMA in which the roughness increases with P is more difficult because these models are applicable only to systems in which adsorbate diffusion (i.e., relaxation after adsorption) is negligible. C. Roughness and porosity in PARM with large Peclet numbers Here we determine the exact values of the exponents γ and λ that relate roughness and porosity with large Peclet numbers in PARM. Subsequently, the approach is extended to PARMA. The reasoning partly rephrases the arguments of Refs. [32,42,54]. In the cases of κ ≪ 1, Eq. (1) gives P≈ d ≫ 1. κ (16) Fully biased flux (κ = 0) leads to RD in which the the height difference between neighboring columns randomly increases. Large fluctuations of heights of neighboring columns also appear in the deposits grown with 0 < κ ≪ 1 as shown in Fig. 3(f). This geometry is schematically illustrated in Fig. 10 in which columns A and B have a large height difference δh. δh A B FIG. 10. Scheme of two neighboring columns of the deposit (light blue) with a large height difference. An incident particle (red) moves along column B and may aggregate at any of the points marked with a cross. The top particles of both columns are highlighted (dark blue). It corresponds to a fluctuation, δN = δh , a (17) in the number of particles accumulated in those columns. Let δT be the number of deposited monolayers in the time that the uncorrelated height difference δh developed. It is related to the fluctuation in the number of particles as δN ≈ (δT )1/2 (18) (recall that T and N are dimensionless quantities). When an incident particle is close to the deposit, its trajectory is almost vertical. If it is moving in column B of Fig. 10, it may aggregate near one of the NNs of column A before reaching the top of column B. The δN possible aggregation points are indicated by crosses in Fig. 10. In each point, the probability of aggregation is κ/(2d ), which is the probability of a hop attempt onto a NN of column A. Thus, the probability that the incident particle aggregates to any of those NNs is κ P1 = δN . (19) 2d This reasoning is applicable to a single incident particle, but the height difference δh was developed whereas δT monolayers were deposited. In this time interval, the average number of particles deposited in any column is δT and the total probability of lateral aggregation of the particles that reach column B is   κ δh 3 Pag ≈ δT P1 ≈ , (20) 2d a where we used Eqs. (17)–(19). When Pag ∼ 1, we expect that, at least, one lateral aggregation occurs in column B. This is sufficient to create a correlation in the heights of columns A and B. Using Eq. (20), it occurs when the height fluctuation reaches a crossover value δhc given by  −1/3 κ δhc ∼ a . (21) 2d 042802-9 REIS, MALLIO, GALINDO, AND HUERTAS PHYSICAL REVIEW E 102, 042802 (2020) The number of deposited monolayers to reach the crossover is  −2/3 κ δTc ∼ , (22) 2d where we used Eqs. (17) and (18). Lateral aggregation is the mechanism that generates the excess velocity characteristic of KPZ roughening [49]. Consequently, a crossover from RD to KPZ is expected after the deposition of δTc monolayers. Since δTc decreases with κ, it increases with P [Eq. (16)], which explains the longer regions with RD for larger P observed in Figs. 4(b) and 7(b). The elongated pores in Fig. 3(f) are the consequence of the large height differences that develop before an event of lateral aggregation occurs. The number of pore sites produced in this event is on the order of δNc ∼ δhc /a [Eq. (17)] because the aggregation can occur at any point in the particle trajectory. In a given column, these pores were formed when δTc particles were deposited; thus, Eqs. (21) and (22) lead to a porosity,  1/3 δhc /a κ δhc /a ≈ ∼ ∼ P −1/3 , (23) ≈ δTc + δhc /a δTc 2d where the last step used Eq. (16). Equation (23) shows that λ = 1/3 in Eq. (11). This is a superuniversal exponent because it has the same value in any spatial dimension. Our numerical estimates are in good agreement with this prediction. The roughness of a KPZ interface varies as in Eq. (9). For large P, the numerical factors omitted in Eq. (9) have to be consistent with the initial RD in which height fluctuations on the order of δhc are attained after the deposition of δTc monolayers. This is possible with   T β WPARM ∼ δhc . (24) δTc Using Eqs. (16), (21), and (22), we obtain WPARM ∼ P 1/3−2β/3 . (25) This result implies γ = 1/3 − 2β/3 in Eq. (10). In d = 2, β = 1/3 [37], which gives γ = 1/9 ≈ 0.111; this is in good agreement with the numerical estimated obtained in Sec. III C. In d = 3, β ≈ 0.242 [47,48], which gives γ ≈ 0.172; the estimate in PARM (Sec. III B) is slightly smaller than this value. D. Roughness and porosity in PARMA with large Peclet numbers FIG. 11. Local porosity of PARM and PARMA as a function of P/D. The dashed line has a slope −1/3. data of Figs. 5(b) and 8(b) and setting D = q. This confirms the scaling in Eq. (26). The effect of a small aggregation rate q on the roughness is described by the same reasoning. Consequently, Eq. (25) is generalized by replacing P by the ratio P/D as  1/3−2β/3 P . (27) WPARMA ∼ D E. GBD with small aggregation probability The above scaling arguments can be extended to GBD. The probability of aggregation κ/(2d ) of PARM, which follows from the lateral hops, is replaced by p in GBD. Thus, the porosity of the deposit is obtained with this substitution in Eq. (23),  ∼ p1/3 . This means λ1 = 1/3 in Eq. (13), which is close to the numerical estimate obtained in Sec. III E. The roughness at constant T is obtained by replacing P in Eq. (25) by p−1 [see Eq. (16)], WGBD ∼ p2β/3−1/3 . Consider the case of PARMA with large P, i.e., small κ. The probability of lateral aggregation, which is small in PARM, is reduced by the factor q in PARMA. Thus, we can replace κ by qκ in the scaling relations (19)–(23). Since the interplay between diffusion and aggregation can be described by D ∼ q, we have qκ ∼ D/P. The general scaling relation expected for the porosity is  −1/3 P ∼ (26) D for large ratios P/D. Figure 11 shows the local porosity of PARM and PARMA as a function of P/D, using the same (28) (29) This leads to γ1 ≈ 0.172 in Eq. (12) for d = 3. The numerical estimate obtained in Sec. III E is more than 20% smaller. However, such deviations are usually expected in numerical estimates of small scaling exponents. Equations (28) and (29) are different from the relations obtained for the porosity and roughness in the competitive models involving BD and RD [29,42] and in related models, such as bidisperse ballistic aggregation [55] or ballistic deposition with bond breaking [56]. For instance, those models show porosity scaling as ∼p1/2 [29,42,55] or as ∼p [56], where p is the probability of a ballistic component, which may lead to lateral aggregation. In those models, the lateral 042802-10 SCALING OF ROUGHNESS AND POROSITY IN THIN … PHYSICAL REVIEW E 102, 042802 (2020) aggregation can occur only at the side of a top particle (of a NN column), i.e., it depends only on the configuration of the outer surface of the deposit. Instead, in GBD, PARM, and PARMA, the aggregation may occur at any point of a wall separating a high and a low column as in Fig. 10. Thus, GBD is the simplest realization of a deposition model with highly collimated transport and with particle adsorption allowed at any contact with the deposit. V. SUMMARY AND CONCLUSION We studied thin film deposition models with particle flux mixing diffusive and collimated components and with aggregation depending on the particle motion near the deposit. In the first model with these features, PARM [18], aggregation of the incident particle occurs only when it attempts to hop to a site occupied by a previously deposited particle. The extended model PARMA considers an additional barrier for the adsorption. The properties of the porous deposits depend on P, which characterizes the transport at a length scale a equal to the particle size, and on D, which characterizes adsorption and transport rates at that scale. PARM has D ∼ 1, and PARMA has D  1, so all models consider adsorption limited conditions or a balance of adsorption and transport rates; this contrasts with models of instantaneous sticking of incident particles, such as BD and DLD. For P ≪ 1, the films are highly porous due to the preferential capture of incoming particles by the highest branches of the deposits. As P increases, this instability is suppressed; for a constant number of deposited monolayers, the film roughness decreases. For P ≫ 1, there is an initial uncorrelated growth, i.e., random deposition, but KPZ scaling is observed with small corrections at long times. For a constant number of deposited layers, the film roughness has a minimal value when growth occurs with P ∼ D1/2 . In PARM, the condition P ∼ 1 reflects a balance of the diffusive and collimated components in the particle flux. In PARMA, the effect of D accounts for the fact that an increase in the adsorption barrier (decrease of D) works against the lateral aggregations that correlates the local heights and smoothens the surface. Using a scaling approach, we also explained the power law relations between the roughness and P and between the porosity and P in PARM with dominant collimated flux. The exponents for the roughness are related to the growth exponent [1] C. Zhu, Y. Fu, and Y. Yu, Adv. Mater. 31, 1803408 (2019). [2] L. M. Peter, Electrochem. Commun. 50, 88 (2015). [3] S. Mondal, I. M. Griffiths, and G. Z. Ramon, J. Membr. Sci. 588, 117166 (2019). [4] J. Tam, G. Palumbo, and U. Erb, Materials 9, 151 (2016). [5] M. J. Vold, J. Phys. Chem. 63, 1608 (1959). [6] A. Barabási and H. E. Stanley, Fractal Concepts in Surface Growth (Cambridge University Press, New York, 1995). [7] P. Meakin, Phys. Rev. A 27, 2616 (1983). [8] P. Meakin, Phys. Rev. B 30, 4207 (1984). [9] T. A. Witten and L. M. Sander, Phys. Rev. Lett. 47, 1400 (1981). β of the KPZ class and are in reasonable agreement with simulation results in d = 2 and d = 3. The porosity shows a superuniversal scaling as P −1/3 , i.e., independent of the spatial dimension. In PARMA, the same scaling relations are obtained by replacing P by P/D. We also studied a generalized ballistic deposition model in which the incident particle trajectory is straight, but the probability of aggregation to lateral neighbors along that trajectory is p < 1. For small p, the roughness and the porosity have similar relations as those of PARM with P replaced by p−1 , which is confirmed numerically. The difference between these models and previous competitive models involving ballistic and random deposition is related to the possibility of adsorption at any contact point of the particle trajectory. Thus, that generalized ballistic deposition model provides the simplest example of a class of models with the so-called slippery particle motion [30]. The results obtained here may be relevant for applications in which the transport properties in a gas or in a solution can be tuned to control the film properties. To predict the conditions for producing films with smoother surfaces, the relation involving the Peclet and the Damkohler numbers may be particularly useful. Moreover, the superuniversal scaling of the porosity may be used to interpret properties of low porosity films. The scaling approaches developed here may also be useful in the study of other models with mixed transport properties and with simple subsurface aggregation conditions, such as those of GBD. Extensions of this type of model should also include relaxation of deposited particles, which may lead to significant changes in the deposit patterns as illustrated in recent works [26–28]. ACKNOWLEDGMENTS F.D.A.A.R. acknowledges support from the Brazilian agencies FAPERJ (Grant No. E-26/202.881/2018), CNPq (Grant No. 305391/2018-6), and CAPES (Grant No. 88887.310427/2018-00-PrInt). D.O.M. acknowledges support from CAPES (Grant No. 88882.332193/2018-01). R.H acknowledges support to Project No. FIS2017-89258-P from Ministerio de Economía, Industria y Competitividad, Agencia Estatal de Investigación, Spain, and from the European Union FEDER (European Regional Development Funds). [10] M. Tassopoulos, J. A. O’Brien, and D. E. Rosner, AIChE J. 35, 967 (1989). [11] A. Sánchez, M. J. Bernal, and J. M. Riveiro, Phys. Rev. E 50, 2427(R) (1994). [12] M. Castro, R. Cuerno, A. Sánchez, and F. Domínguez-Adame, Phys. Rev. E 57, 2491(R) (1998). [13] M. Castro, R. Cuerno, A. Sanchez, and F. Dominguez-Adame, Phys. Rev. E 62, 161 (2000). [14] S. C. Ferreira, S. G. Alves, A. Faissal Brito, and J. G. Moreira, Phys. Rev. E 71, 051402 (2005). [15] D. Rodríguez-Pérez, J. L. Castillo, and J. C. Antoranz, Phys. Rev. E 72, 021403 (2005). 042802-11 REIS, MALLIO, GALINDO, AND HUERTAS PHYSICAL REVIEW E 102, 042802 (2020) [16] S. G. Alves and S. C. Ferreira, Jr., Phys. Rev. E 73, 051401 (2006). [17] D. Rodríguez-Pérez, J. L. Castillo, and J. C. Antoranz, Phys. Rev. E 76, 011407 (2007). [18] J. L. Galindo and R. Huertas, J. Appl. Phys. 114, 064905 (2013). [19] T. J. Oliveira and F. D. A. Aarão Reis, J. Stat. Mech. (2014) P09006. [20] J. L. Galindo, R. Huertas, A. Carrasco-Sanz, A. Lapresta, J. Galindo, and E. Vasco, J. Appl. Phys. 120, 034902 (2016). [21] M. J. Kartha, Phys. Lett. A 381, 556 (2017). [22] Y. P. Pellegrini and R. Jullien, Phys. Rev. Lett. 64, 1745 (1990). [23] S. Sadhukhan, T. Dutta, and S. Tarafdar, J. Stat. Mech. (2007) P06006. [24] F. D. A. A. Reis, D. di Caprio, and A. Taleb, Phys. Rev. E 96, 022805 (2017). [25] Z. Xun, Z. Zhang, Y. Chen, L. Wu, and G. Tang, Physica A 471, 569 (2017). [26] D. di Caprio, A. Taleb, and F. D. A. Aarão Reis, J. Phys. Chem. C 122, 21418 (2018). [27] F. Hao, A. Verma, and P. P. Mukherjee, ACS Appl. Mater. Interfaces 10, 26320 (2018). [28] B. S. Vishnugopi, F. Hao, A. Verma, and P. P. Mukherjee, Phys. Chem. Chem. Phys. 22, 11286 (2020). [29] C. M. Horowitz and E. V. Albano, J. Phys. A 34, 357 (2001). [30] A. Robledo, C. N. Grabill, S. M. Kuebler, A. Dutta, H. Heinrich, and A. Bhattacharya, Phys. Rev. E 83, 051604 (2011). [31] K. Banerjee, J. Shamanna, and S. Ray, Phys. Rev. E 90, 022111 (2014). [32] F. D. A. Aarão Reis, Phys. Rev. E 91, 062401 (2015). [33] B. Mal, S. Ray, and J. Shamanna, Phys. Rev. E 93, 022121 (2016). [34] Y. Kameya, J. Nanopart. Res. 19, 214 (2017). [35] F. Higuera, J. Aerosol Sci. 118, 45 (2018). [36] T. Yang and Y. Han, Cryst. Growth Des. 16, 2850 (2016). [37] M. Kardar, G. Parisi, and Y. C. Zhang, Phys. Rev. Lett. 56, 889 (1986). [38] J. L. Galindo, F. D. A. Aarão Reis, R. Huertas, A. CarrascoSanz, J. Galindo, V. Medina, and E. Vasco, Thin Solid Films 690, 137448 (2019). [39] D. E. Rosner, Transport Processes in Chemically Reacting Flow Systems (Dover, Mineola, NY, 2000). [40] H. S. Fogler, Elements of Chemical Reaction Engineering, 4th ed., Prentice Hall PTR International Series in the Physical and Chemical Engineering Sciences (Prentice Hall, Upper Saddle River, NJ, 2006). [41] P. Nath and D. Jana, Phys. Sci. Rev. 4, 20170109 (2019). [42] F. D. A. Aarão Reis, Phys. Rev. E 73, 021605 (2006). [43] J. Krug, J. Phys. A 23, L987 (1990). [44] J. Tang and A. Gomez, Aerosol Sci. Technol. 51, 755 (2017). [45] J. Krug, Adv. Phys. 46, 139 (1997). [46] Y.-J. Chen, S. Watanabe, and K. Yoshikawa, J. Phys. Chem. C 119, 12508 (2015). [47] J. Kelling and G. Ódor, Phys. Rev. E 84, 061150 (2011). [48] J. Kelling, Géza Ódor, and S. Gemming, J. Phys. A: Math. Theor. 51, 035003 (2018). [49] W. E. Hagston and H. Ketterl, Phys. Rev. E 59, 2699 (1999). [50] S. G. Alves, T. J. Oliveira, and S. C. Ferreira, Phys. Rev. E 90, 052405 (2014). [51] E. E. M. Luis, T. A. de Assis, S. C. Ferreira, and R. F. S. Andrade, Phys. Rev. E 99, 022801 (2019). [52] R. Cuerno and M. Castro, Phys. Rev. Lett. 87, 236103 (2001). [53] P. A. Orrillo, S. N. Santalla, R. Cuerno, L. Vázquez, S. B. Ribotta, L. M. Gassa, F. J. Mompean, R. C. Salvarezza, and M. E. Vela, Sci. Rep. 7, 17997 (2017). [54] C. M. Horowitz and E. V. Albano, Phys. Rev. E 73, 031111 (2006). [55] F. A. Silveira and F. D. A. Aarão Reis, Phys. Rev. E 75, 061608 (2007). [56] J. S. O. Filho, T. J. Oliveira, and J. A. Redinz, Physica A 392, 2479 (2013). 042802-12