Phase-field modeling and n-point polytope characterization of nanostructured protuberances formed during vapor-deposition of phase-separating alloy films

Distinct surface features such as grooves and protrusions formed in monolithic and alloyed films can influence their optical, electronic, and mechanical properties. The presence of these nanometer-sized defects in films, which are a key component of electronic devices and integrated circuits, often causes catastrophic failure. Therefore, it is important to obtain a fundamental understanding of the mechanisms by which surface irregularities form and evolve in alloy films.

Over the years, numerous studies on the growth and control of surface features^1–4 have been conducted. In monolithic films, the formation and growth of hillocks have been attributed to residual stresses arising due to a difference in the coefficient of thermal expansion between the film and the substrate.^1,5,6 Beyond residual stresses, the formation of irregular features on the film surface can be influenced by several other factors such as the lattice misfit strain,⁷ the energy of the substrate–film interface,^1,8 and the film thickness.⁹ 

More recent studies explored the morphological evolution of surface irregularities in binary metallic alloy films during the deposition and heat treatment,^10–12 implying that factors beyond the substrate and residual stresses may be responsible for their occurrence. Xue et al.¹³ attributed the formation of surface irregularities in Cu–W co-sputtered films to surface diffusion of atoms. However, their interpretation does not account for the impact of phase-separated Cu- and W-rich domains on surface irregularities. Powers et al.¹⁴ proposed that dissimilar mobilities between the alloying elements can influence hillock formation in immiscible phase-separating films of Cu–Ta and Cu–Mo–Ag. This finding was corroborated through SEM characterization of Cu–Ta vapor co-deposited films, which indicated a greater tendency of hillock formation with an increase in the deposition temperature. While it can be argued that adatom diffusion can exacerbate the size surface protuberances due to differing interdiffusion coefficients of the co-deposited elements, the underlying driving force that is responsible for the nucleation of such surface instabilities requires deeper analyses. Moreover, a comprehensive study analyzing the role of thermodynamic and kinetic factors on the morphological evolution of surface protuberances has not been performed before. Previous experiments and modeling have revealed the occurrence of nanoscale concentration modulations in immiscible phase-separating alloy films of Cu–Ta and Cu–Mo.^15–18 The co-deposited film nanostructures comprised of vertical concentration modulations (VCMs), lateral concentration modulations (LCMs), rod-shaped domains, or a mixture of variants, depending on the alloy composition and process parameters such as deposition rate and temperature. Raghavan et al.¹⁹ extended the theory of Srolovitz and co-workers^20,21 by demonstrating that the nanostructural transitions were governed by interfacial energy minimization as the deposition rate was increased. Therefore, the role of interfacial energies, and in particular, the influence of surface contact angles and energy on the formation of surface protuberances in co-deposited films warrants further investigation. The energy of the interface between phase-separated domains of Cu–Ta, Ag–W, and Ag–Cr films is comparable to or greater than the surface energy of the alloyed components due to the incoherent nature of the interphase boundaries.^22–24 This allows for dihedral angles (or contact angles), $θ$ at the film-vapor interface to be less than $180°$ under equilibrium conditions. Based on the Young equation described in Sec. II, the equilibrium contact angles for immiscible phase-separating films of Cu–Ta, Ag–Cr, and Ag–W are calculated to equal $112°$⁠, $92°$⁠, and $76°$⁠, respectively.^22–24 Therefore, it is necessary to investigate whether or not the initiation and/or growth of surface protuberances is influenced by the energy of the phase-separating domain boundaries.

While a fundamental understanding of the multiphysics leading to protuberances in vapor-deposited phase-separating films is essential, rigorous quantification of the nanostructured modulations and features, including defects, using statistical correlation tools is required for a comprehensive understanding. Prior studies focused on statistical analyses to assert relationships between the hillock size and spatial distribution in stress-relieved monolithic films.^9,25,26 Other studies have used atomistic and finite element modeling techniques to better understand the mechanisms leading to hillock formation in monolithic films.^7,27 However, as these previously reported techniques for characterizing the evolution of nanostructured features are limited to a selected list of statistical measures (e.g., size and shape distribution, average hillock distance, and orientation relationships), they are unable to expose the underlying symmetry of emerging nanoscaled features and their associated length scales during the morphological evolution of vapor-deposited films. Moreover, the availability of literature on applying numerical or statistical methods for understanding the evolution of surface protuberances in immiscible alloy films is sparse. Stewart and Dingreville²⁸ recently developed an advection-based model of vapor-deposited phase-separating films. However, this study was primarily limited to 2D, and any influence of the variability in surface contact angles, which can significantly alter the morphological evolution of protuberances, was not addressed. Limiting the simulations to 2D intrinsically restricts the characterization of hierarchical symmetries of the amplifying protrusions, which can infer a host of insights into the growth mechanisms. This limitation highlights a pressing need to couple 3D phase-field simulations of nanostructure evolution with advanced characterization techniques that leverage statistical probability functions capable of identifying symmetrical patterns while simultaneously quantifying the corresponding length scales as the protuberances evolve on the surface of vapor-deposited films. Here, we propose a 3D phase-field approach that assesses the role of surface and interfacial energies in the formation of protuberances on the film surface. Our calculations account for the variability in contact angles, seed morphology, and adatom kinetics that are found to significantly influence the morphological evolution of surface protrusions. We then characterize the evolving surface instabilities via a novel set of statistical morphological descriptors, i.e., the n-point polytope functions $Pn$⁠. The $Pn$ function is associated with the probability of finding n-point distributions specified by a regular n-polytope (polygons in 2D) in the phase-separating domains.²⁹ It has previously been shown that $Pn$ functions can effectively decompose the structural features of interest into a set of polytope basis, allowing for easy detection of any underlying symmetries or emerging features and their associated length scales during nanostructural evolution.^29,30 In the following study, we isolate unique surface patterns based on the relative prevalence of these dominant length scales. We also simulate film deposition in 3D at distinct surface contact angles that enable the formation of surface grooves as the pure domains phase separate by starting from a random seed comprising of small concentration fluctuations (Fig. 1). Finally, we apply statistical $Pn$ analysis to these computationally generated films to identify and compare the nanostrutured protuberances emerging from these films with the micrographs of sputtered Cu–W and Cu–Mo–Ag films.¹⁴ 

In this section, we discuss a ternary Cahn–Hilliard model to explore the influence of surface contact angles on the nanostructural evolution of binary immiscible alloys films of A-B type, in which the phase separation results in the formation of pure A and B domains, while the vapor is considered as the third phase.^31–33 The evolution of the film nanostructure is driven by a phenomenological minimization of the total free energy. The total free energy functional can be written as

where $Nv$ is the number of molecules per unit volume (assumed independent of composition and position) and $κi$ (⁠$i=A,B$⁠, and $ν$⁠) are the gradient energy coefficients associated with the concentration field. The order parameters, $ϕA$ and $ϕB$ denote the scaled concentrations of the alloying elements, A and B in the film, while $ϕv$ corresponds to the vapor phase, which is in contact with the film’s surface. In the interest of mass conservation, the sum of the order parameters is constrained to be unity by imposing

The chemical free energy per molecule, $f(ϕA,ϕB,ϕv)$⁠, corresponds to a regular solution such that

where $χij$ (⁠$i,j=A,B,v$⁠; $i≠j$⁠) is the pairwise interaction energy between the components, $kB$ is the Boltzmann constant, and $T$ is the absolute temperature. An expanded form of Eq. (3) can be written as

We can eliminate one of the field variables by replacing $ϕv$ with $1−ϕA−ϕB$ in Eq. (4) and rewrite as

The chemical potential of A-rich and B-rich phases are obtained by deriving the variational derivatives of the free energy functional in Eq. (3) with respect to $ϕA$ and $ϕB$ such that

and

where $κAA=κA+κv$⁠, $κBB=κB+κv$⁠, and $κAB=κBA=κv$⁠.

The kinetics of phase separation is obtained via a continuity equation given by

where $Ji′$ is the total flux of each component in the system. We adopt a formulation that incorporates the net vacancy flux, $J°$⁠, during the diffusion process as reported earlier by Kramer et al.^34,35 and Bhattacharyya et al.,^36,37 to obtain the net flux of each component, $i$⁠, given by

$Ji$ is given by

where $μi$ and $Mi$ are the chemical potential per site and Onsager coefficient of the $i$th component, respectively. The vacancy flux, $J°$⁠, is given by

By substituting Eqs. (11) and (9) in Eq. (8) and using the Gibbs–Duhem relationship,^35,37 we arrive at the following expressions for the temporal evolution of A and B-rich phases:

and

where $MAA$ and $MBB$ are the atomic mobilities of A and B atoms in non-A-rich and non-B-rich phases, respectively, while $MAB$ and $MBA$ are mobilities of A atoms in B-rich phase and B atoms in A-rich phase, respectively.^35,36 These are coupled to the diffusion coefficients of the alloying components, $Di$⁠, by extending the Nernst–Einstein relation to a ternary system³⁸ 

and

We solve Eqs. (12) and (13) by first non-dimensionalizing the parameters using the relation $l∗=(κi/2kBT)1/2Δx$ and $t∗=(kBT/Mii∗l∗2)$⁠, where $l∗$ and $t∗$ are the characteristic length and time, respectively, and $Mii∗$ is the dimensional value of mobility for phases $i=A,B$⁠. The dimensionless form of the equations is solved via an explicit finite difference scheme for temporal and spatial derivatives.

In vapor-deposited films where the nanostructural evolution is characterized by the formation of phase-separating domains, interfacial energies play a dominant role in the morphological self-assembly. The role of excess energies at the interphase and surface boundaries can be quantified by measuring contact angle, $θ$⁠, at the steady state, which is governed by a balance of forces at the triple junction of A, B, and vapor phases. This is encapsulated within Young’s equation as

where $θ$ is the surface contact angle, $σAB$ is the interfacial energy between the phase-separated pure domains of A and B, while $σfv$ equals the energy of film surface which is in contact with the vapor phase. Here, $f$ denotes either the A-rich or the B-rich phase, depending on which of the two constitutes the surface that is in contact with the vapor. We would like to clarify that the interfacial energies of either domain that are in contact with the vapor are assumed to be equal, which in turn leads to the formation of symmetrical surface grooves.

The excess energy at the interface between two phases, defined by components, A and B, is given by the expression³¹ 

where $Nv$ denotes the number of molecules per unit volume, and $λ$ is the interaction distance, which is related to the intermolecular distance and is assumed to be a constant for the alloy. $ϕ$ is the composition field that varies in value from $ϕA$ within the A-rich phase to $ϕB$ within the B-rich phase. This equation is solved by approximating the interfacial energy based on the bulk free energy at its maximum value,³¹ given by

where $ϕe$ is the equilibrium value attained by $ϕ$ when the system reaches a steady state. In a binary system, a single order parameter, $ϕ$⁠, which varies from A-rich and B-rich phases, is sufficient to define its compositional layout. However, in a three-component system with an additional vapor phase, the excess energy at the interface would depend on not one but two order parameters that operate independently from one another. This renders the calculation of interfacial energy from Eq. (17), and therein, the contact angle, non-trivial. However, a relationship between the excess energy associated with the interphase or surface boundaries and the interaction parameters corresponding to the alloy is known as per Eq. (18). We, therefore, tune the excess energies by applying different values for $χAv$ and $χBv$ while maintaining $χAB$ a constant. The ratio of the excess energies at the surface and interphase boundaries would then yield distinct dihedral angles depending on the force balance given by Eq. (16). We verify the exact value of the dihedral angles for every case by deducing the angles obtained at steady-state from simulations of bulk phase separations of a three-phase system that comprises A-rich, B-rich, and vapor phases, as discussed in the  Appendix.

The Cahn–Hilliard equations (12) and (13) are solved in a three-dimensional computational domain whose dimension varies along the deposition axis prior to attaining a maximum size of $150×150×150$ grid points. The assigned energy scale, which depends on the interfacial energy of the phase-separating domains, yields a characteristic length scale, $l∗$⁠, equal to $1.1$ nm. Although the model developed in this work can be generically applied to analyze nanostructural evolution in any immiscible alloy, the chosen values of interfacial energies correspond to those of Cu–Mo and Cu–Ta available in the literature.^15,16 Other simulation parameters are non-dimensionalized as per the method outlined in our previous work.¹⁹ The set of non-dimensionalized deposition parameters used to perform the simulations reported in this article are listed in Table I.

The initial condition for simulations reported below comprises of seed that either has a perturbed topology with VCM or is topologically flat and composed of Langevin-type compositional noise whose mean value preserves the overall film composition. The top half of the simulation domain corresponds to a vapor phase. To mimic the process of deposition, a layer of grid points with a topology identical to the film surface, given by $ϕv=0.5$⁠, and comprises of Langevin noise is overlaid after periodic intervals. The duration of the interval assumed to be constant corresponds to the simulated deposition rate. Apparently, our method warrants tracking the topology of the undulated film surface prior to every deposition event. To maintain the relative concentration ratios conserved, a layer of thickness $Δx$ of the vapor phase is also added to the top of the simulation domain following every deposition event, as shown in Fig. 2. The free energy functional (1) does not consider the influence of misfit strains that may arise due to differences in the crystal structure of phase-separating domains. Therefore, the above formulation of phase-field equations is best suited for simulating nanostructural evolution in those films where lattice strains are negligible.

For the characterization of film nanostructures, we employ a recently developed set of statistical morphological descriptors, i.e., the $n$-point polytope functions $Pn$^29,30 to quantify the nanostructural evolution characteristics of the vapor-deposited immiscible alloy films. Without loss of generality, we consider the $Pn$ functions defined for a heterogeneous material (represented as a random field) in which the different structural features are segmented and grouped into different “phases.” The definition of $Pn$ is then given as follows:

$Pn(r)≡$ Probability that all of the $n$ vertices of a randomly selected regular $n$-point polytope with edge length $r$ fall into the phase of interest.

Here, we focus on the $n$-polygon functions (see Fig. 3), involving $n$-point regular polygons, for which the vertex (edge) number $n$ can take any positive integer values. We note that the $n$-point polygonal functions can be computed from both single 2D slice and 3D realizations by averaging several 2D slices along the distinct planes. It is apparent from the above definition that $Pn$ functions belong to the subset of corresponding standard $n$-point probability functions $Sn$⁠, which gives the probability of occurrence of specific $n$-point configurations in the phase of interest.^39–42 For $n=2$⁠, the 2-point polytope function $P2$ is identical to the standard 2-point correlation function $S2$ in the case of statistically homogeneous and isotropic systems, which has been extensively employed to quantify complex material systems and microstructural evolution.^43–46 

Figure 3 schematically illustrates the stochastic events that contribute to the $Pn$ functions. For $r=0$⁠, the polygon reduces to a single point, and $Pn(r=0)$ corresponds to the probability of a randomly selected point falling in a the phase of interest, $η$⁠, which equals its volume fraction. For non-zero $r$ values, $Pn(r)$ provides $n$-point spatial correlations in the phase of interest. For extremely large $r$ values (e.g., $r→∞$⁠), the probabilities of finding the vertices within phases of interest are almost independent of one another; thus, we have $Pn(r→∞)≈ηn$⁠.

The probability-based definition of the $Pn$ functions allows us to easily compute these functions from microstructural data, including both 2D images and 3D digital representations of the material nanostructures. To compute the value of $Pn(r)$ with $r=a$⁠, where $a$ is the edge length of the corresponding regular n-polytope (e.g., a $n$-polygon), the following procedure is used: the $n$-polygons of edge length $a$ are generated $M$ number of times such that their center and orientation within the nanostructure are randomly selected. Every time a polygon is embedded on the micrograph, we determine if the corresponding vertices fall into the phase of interest (i.e., a “successful” event), and the total number of successful events $Ms$ out of a total of $M$ trials is subsequently recorded. The rate of success $Pn(r=a)=Ms/M$ is computed, which equals the probability that a randomly selected $n$-polygon having all its vertices falling into the phase of interest. The above procedure is repeated for different $r$ values to deduce an entire $Pn(r)$ function.

We note that by definition, the $Pn$ functions can easily capture the underlying (hidden) symmetry of the pattern of interest. In our following analysis of simulated films, we apply the $Pn$ function to the upper-most concentration modulations of the protruding nanostructures as visualized in the plan-view by segmenting the top-most A-rich (blue) and B-rich (red) layered phase domains and applying the n-polytope functions to each of the phase domains. We use the characteristic lengthscale of $l∗=1.1$ nm to translate the film dimensions in our plots of the $Pn$ function. In the following article, we mainly report the results of $P2$ (⁠$S2$⁠) and $P4$ analysis. The former is a commonly used generic descriptor, while the latter captures the fourfold symmetry of the deposited films, which is inherited from the symmetries that are present in the initial seed morphology.

We first discuss the characteristic protuberances formed at measured $θ≈32°,41°$⁠, and $51°$ from a seed layer with a perturbed topology with alternating layers of phase-separated domains. Our choice of the seed morphology is based on commonly observed vertical concentration modulations in previous TEM characterization of sputtered films.¹⁶ Next, we leverage the $Pn$ function to analyze the morphological self-assembly and the evolution of pattern symmetries. Finally, we simulate and characterize the films starting from seed with a random topology and comprising of Langevin noise, and compare with experimentally observed hillock formations in sputtered films of Cu–Ta and Cu–Mo–Ag.¹⁴ 

To better understand the underlying mechanisms by which protuberances evolve, we simulate the nanostructural evolution during the vapor deposition of binary immiscible films using the phase-field model described in Sec. II. A schematic diagram in Fig. 2 shows the hemispherical seed morphology composed of equidistant perturbations of radius $17Δx$⁠, such that every protrusion has four nearest neighbors. The alternating layers of pure A and B phases maintain an overall film composition of 50 at. % B. Deposition is simulated by adding a noisy layer of thickness equal to $Δx$ over the film’s surface at regular time intervals, which is numerically equal to the inverse of the deposition rate, $R$ [Fig. 2(b)]. $R$ is non-dimensionalized and equals $0.1$ for all simulations reported in this study. Since the formation of nanostructured domains is known to be guided by an interplay of deposition rate and atomic mobility, which is related to temperature, it is likely that the contact angles specified by Eq. (16) are never attained during the course of deposition simulations. However, the following study analyzes how the morphological evolution of protuberances is altered at distinct ratios of the surface energy of the film with respect to the vapor phase and the interfacial energy of the A-B interface.

When the contact angle, $θ$⁠, is assigned to be equal to $32°$⁠, we observe the evolution of distinct nanostructured patterns as the deposition progresses. At $t=12$⁠, the seed topology is largely retained, barring a minor coarsening of the phase-separated domains below, as seen in Fig. 4(a). At $t=612$⁠, the surface flattens with the top-most layer of adjacent protrusions fusing with one another, as shown in Fig. 4(b). However, the average width of the protrusions remains unchanged. As the deposition continues at $t=1212$⁠, we observe the perturbed seed regions (first-tier protrusions) transition into new second-tier protrusions or protuberances. These protuberances have a checkerboard morphology with their edges fusing with one another, leaving pits through which remnants of the original protrusions are visible in the plan-view [Fig. 4(c)]. At $t=1812$⁠, we see another transition leading to the formation of a third-tier of protuberances over the original first-tier protrusions. These protuberances appear to form isolated islands without any interlinking of domains at the corners. Finally, at $t=2412$ [Fig. 4(e)], the isolated islands give way to fused protuberances with domains interlinked at the corners, manifesting as an open honeycomb nanostructure. No new patterns emerge after this point as the deposition continues, which indicates that the surface protuberances of the film nanostructure have stabilized.

Next, we simulate a film with a lesser skew between the surface and interfacial energies, which results in a contact angle of $41°$⁠. In Fig. 5, we can observe that after $t=12$⁠, where the seed geometry is retained, the film surface comprises of shallower hills and pits as compared to Fig. 4. At $t=612$⁠, the first-tier protrusions are readily substituted by a continuous domain of the A-rich phase that stretches over the surface of the film. However, the surface topology remains undulated with tiny islands of the B-rich phase emerging in the pits. At $t=1212$⁠, the continuous domains give way to a honeycomb-like nanostructured pattern, with narrow pits forming over the region corresponding to former valleys of the first-tier protrusions [Fig. 5(c)]. As the deposition continues, the honeycomb-like film nanostructure is found to be stabilized at $t=1812$ and $2412$⁠. When the contact angle is further increased to $51°$⁠, the nanostructural evolution is altered as depicted in Fig. 6. After $t=612$⁠, a continuous layer of A-rich domain forms that completely envelopes the original first-tier protrusions formed earlier at $t=12$ [Figs. 6(a) and 6(b)]. Surface undulations are smaller compared to Fig. 5, which corresponds to a smaller contact angle of $41°$⁠. At $t=1212$⁠, the honeycomb nanostructure emerges, albeit with narrower pits and widely spaced hills. Between $t=1812$ and 2412, this honeycomb nanostructure, which is coarser than the one shown in Fig. 5(e), stabilizes.

We now employ the polytope functions, $P2$ and $P4$⁠, to quantify the simulated nanostructural evolution in vapor-deposited phase-separating films. In particular, for a fixed contact angle $θ$⁠, we select three representative snapshots of the top surface of the nanostructured film as the deposition continues and segment the phase of interest. Without any loss of generality, we mainly focus on the blue (A-rich) phase in our analysis. Given that the film composition is equiatomic, the blue and the red (B-rich) phases are complementary to each other. Therefore, analyzing one of the phases is sufficient to statistically quantify the nanostructured film pattern.

Figures 7(a)–7(c) show the computed functions from the simulated thin film nanostructures for all the three contact angles $θ≈32°$⁠, $41°$⁠, and $51°$⁠, respectively. Varying the contact angle alters the evolution pathway, specifically, for $θ≈32°$⁠, where the morphology of the blue phase evolves from an interconnected network to intermediate compact clusters, which later stabilize into a connected complex topology. For the other two angles, regular nanostructures comprising pits and hills develop as the deposition proceeds, although the spacings increase at a larger contact angle. Nonetheless, the $Pn$ functions reveal universal hidden correlations between these distinct nanostructural evolution pathways. In particular, the peaks in $P2$ and $P4$⁠, respectively, indicate characteristic length scales associated with strong two-point and 4-point correlations. Apparently, these length scales all coincide with one another, regardless of the underlying geometry and topology of the film nanostructure. Also, the strong peaks in $P4$ at $r≈30$ and 60 nm indicate a predominant fourfold symmetry, which is inherited from the seed morphology and preserved during the evolution of the film nanostructure.

In order to analyze the impact of seed morphology on the formation of surface protrusions in vapor-deposited films, we simulate the nanostructural evolution in films by starting from random seed morphology, where a noisy seed layer of thickness $10Δx$ acts as a template upon which new layers are deposited at time intervals that correspond to the rate, $R$⁠. Phase-field simulations of the nanostructural evolution during the physical vapor deposition were carried out at distinct $θ$ values approximately equaling $32°$⁠, $41°$⁠, and $51°$⁠. Figure 8(a) shows a film nanostructure deposited at $θ≈32°$ starting from a random seed, wherein the side and top views depict deep grooves and periodically spaced protrusions. Increasing the contact angle reduces the pit depth while marginally increasing the protrusion spacing, as evident from the formation of shallower grooves in Figs. 8(b) and 8(c). This observation is duly supported by a decrease in film densities, $ρ$⁠, as a function of film thickness, shown by the plot in Fig. 8(d), which is obtained by post-processing the computational film nanostructure simulated at the three contact angles. Here, $ρ$⁠, which corresponds to the normalized density of the film, is computed for every 2D slice such that the normal vectors to individual slices are aligned along the deposition axis. Thus, the value of $ρ$ varies from a maximum value of 1.0 in basal regions of the film where voids or protuberances are absent to a lower limit below 0.1 that corresponds to the film surface. When $θ≈32°$⁠, $ρ$ decreases steeply from a maximum value of $1.0$ at $h=30Δx$⁠. For the film nanostructure corresponding to $θ≈41°$⁠, $ρ$ decreases at a relatively gradual rate, whereas for $θ≈51°$⁠, the drop from $ρ=1.0$ is most gentle with respect to the other two cases until $h=80Δx$⁠, before it reaches 0.64. It has to be noted that for all the three cases that pertain to distinct contact angles, $ρ$ plummets beyond certain $h$ that corresponds to the root of the protuberances.

To characterize the nanostructured protuberances reported in Fig. 8, we quantify the topography of the deposited films using the 2-point function $P2$⁠, as shown by the plot in Fig. 9(a). Since the nanostructures in this case, which emerged from the random seeds, do not possess long-range correlations, the associated $P2$ functions exhibit the well-known damped oscillations.^39,40 Apparently, the length scale associated with the first local minimum and maximum corresponds to the average width of the domains and the inter-domain spacing, respectively. It can be observed that although changing the contact angle does not significantly affect the domain size and spacing from a plan-view perspective, the oscillations show greater damping at larger $θ$s. A better understanding of these nanostructured protuberance can be attained by quantifying the disordered layered morphology resulting from the random seeds from a side-view perspective. Since the layered morphology exhibits directional anisotropy due to the preferred alignment of interfaces perpendicular to the deposition axis, we choose the directional $P2$ functions along the orthogonal directions, respectively. In Fig. 9(b), the directional $P2$ functions plotted in the vertical direction that aligns with the deposition axis clearly indicate a preferred alignment of interfaces into VCM, described previously by Ankit et al. and others.^15,16,19 Moreover, the wavelength associated with the oscillations, i.e., the average distance between successive peaks, increases as the contact angle $θ$ increases, indicating effective coarsening of the layers. Along horizontal directions that are orthogonal to the deposition axis, the $P2$ functions are long-ranged and exhibit weaker oscillations.

In order to compare the nanostructural characteristics of the simulated films with those observed in vapor-deposited experiments, we leverage the plan-view SEM micrographs of Cu–Ta and Cu–Mo–Ag sputtered films reported earlier by Powers et al.¹⁴ The contour lines plotted on the SEM micrographs in Figs. 10(a)–10(c) demarcate hillock roots. Similarly, we algorithmically map and segment the plan-view images of simulated films in Fig. 8 to a depth of $15Δz$ to reveal protuberance spacing, as shown in Fig. 10. On comparing the simulated nanostructures [Figs. 10(d)–10(f)] with the previously reported experimental micrographs [Figs. 10(a)–10(c)], we observe several similarities. In particular, Figs. 10(b) and 10(c) show hillocks that appear as disconnected islands, similar to the ones in Figs. 10(d)–10(f). As per the experimental findings,¹⁴ the hillocks become more distinct as the deposition temperature is increased. However, it is to be noted that the temperature is assumed to be the same for all the phase-field simulations reported in this article, while altering the contact angles leads to an observation similar to the experiments. It is also worth clarifying that due to a predominantly non-equilibrium nature of the deposition process, the contact angles that actually develop at the film surface may differ with respect to the prescribed one whose value is derived based on Eq. (16). Our assertion is duly supported by the sequence of nanostructural evolution shown in Figs. 4–6, where the seed undulations first flatten out before secondary protrusions form that are distinct in shape. While it is principally difficult to determine the exact values of $θ$ that pertain to non-equilibrium conditions during the deposition process, based on our simulations, we can certainly conclude that prescribing a smaller contact angle enhances the propensity of protuberance formation, as observed in Cu–Ta and Ag–W films. From a materials design perspective, this naturally boils down to selecting an appropriate immiscible alloy such that surface contact angles are in accordance with the desired surface roughness.

It is well known that the interface between phase-separated domains has orientation dependencies. The energetically stable domain orientations tend to develop during post-deposition annealing of films where there is sufficient time for diffusion occur. However, in non-equilibrium processes, such as sputtering, that have been reported in this work, energetically favored orientations are typically never attained primarily due to a persistent flux of atoms at the film surface, which competes with atomic diffusion that is responsible for phase separation. This gives rise to surface undulations which may amplify over time, especially in those films where the ratio of energies pertaining to domain boundaries and surface lends to a smaller contact angle. Otherwise, the incipient surface modulations tend to flatten out as the deposition proceeds.

Finally, it is worth clarifying that the current phase-field simulations do not account for residual stresses that can influence protuberance growth in vapor-deposited films.^14,17 Additionally, the influence of lattice mismatch, grain boundaries, and the deposition rates on protuberance formation are worth exploring. Therefore, in the future, we will incorporate the multiphysics associated with grain boundaries, lattice mismatch, and residual stresses to comprehend the evolution of nanostructured domains in immiscible alloy films.

In this work, we have formulated a phase-field model for simulating the evolution of periodic and irregular nanostructured protuberances in phase-separating films during the physical vapor deposition of binary immiscible alloys. When the deposition was simulated starting from predefined seed morphology, phase-separated domains self-assembled into VCMs leading up to the formation of protuberances that organized into a characteristic honeycomb nanostructure. When the ratio of interfacial to surface energies of the film is adjusted to prescribe a contact angle of $32°$⁠, the nanostructural evolution of the protuberances was found to occur in three stages; the seed modulations first amplified and then flattened out to form a checkerboard pattern, which ultimately stabilized as a honeycomb nanostructure. For larger contact angles, distinct stages of growth are limited to two; the seed modulations first amplify before self-assembling into interconnected protuberances in a honeycomb-like array. The contact angle is also found to alter the packing density of these nanostructured protuberances. At $θ≈32°$⁠, the spacing between the protuberances is large, as they barely interlink at the corners to define an open honeycomb nanostructure. At $θ≈41°$⁠, the packing of protuberances is denser due to an enhanced interlinking of corners. The packing is most dense for $θ≈51°$⁠, where the intermittent holes within the honeycomb pattern appear to be much narrower when compared to the corresponding nanostructure simulated at a contact angle of $41°$⁠. $Pn$ analysis of the evolving protuberances at different contact angles reveals that the emerging nanostructured domains of A-rich and B-rich phases are complementary as inferred from the top-view. The $P2$ function reveals twofold symmetry and indicates characteristic length scales peaking at $r≈54$ nm irrespective of the $θ$ revealing a hidden correlation omnipresent at all angles. However, the width of peaks is greater as $θ$ increases, indicating a greater spread of the length scales at larger contact angles. The $P4$ function reveals fourfold symmetry with strong peaks at $r=29$ and 57 nm, which is inherited from the seed morphology.

Vapor co-deposition with a noisy seed that mimics an actual vapor deposition process reveals widely spaced protuberances with deep grooves when the prescribed angle, $θ≈32°$⁠. The groove depths decrease when $θ$ is increased to $41°$ and $51°$⁠, respectively, albeit with a marginal increase in the protuberance spacing. $Pn$ analysis of these computational film nanostructures reveals a decrease in the domain width and inter-domain spacings with an increase in $θ$⁠. When we compare our simulations with SEM micrographs of immiscible Cu–Ta and Cu–Mo–Ag films,¹⁴ we are able to identify several similarities in size and spacing of the nanostructured features.

Our work reported in this article strongly corroborates the hypothesis that surface instigation of protuberances is closely linked to the ratios between domain boundary and surface energies. Smaller contact angles between the film and the vapor phases promote amplification of surface modulations leading to the formation of steep protuberances. In comparison, larger contact angles lead to shallower and more widely spaced protuberances. We anticipate that the fundamental insights gained from this study will provide guidance for depositing films with preferred surface morphologies for optical applications or prevent surface irregularities in cases where they are undesirable, for example, in semiconductor and solar cell interconnects.

Financial support from the NSF MEP Program via Award No. 1763128 is gratefully acknowledged.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

It is well known that the energetics of phase separation in a regular solution is governed by the interaction parameters, $χ$⁠.³¹ Therefore, based on Young’s equation [Eq. (16)], the dihedral angle is also dependent on the relative values of the interaction parameters, $χAB$⁠, $χAv$⁠, and $χBv$⁠. To validate our phase-field model so that it accurately captures the physics of force balance at phase triple junctions, we analyze the surface grooving phenomenon in steady-state films comprising two equal-sized domains of pure A and pure B components in contact with the vapor. Simulations are performed in a large numerical 2D domain of size $2000×400$ grid points, divided equally into domains of pure A, B, and vapor components, as shown in Fig. 11. The interaction parameters, $χAB$⁠, $χAv$⁠, and $χBv$⁠, are varied to change the excess energies at the interface between A-rich and B-rich phases, and the surface between film and vapor. No flux boundary conditions are assigned along the two principal axes, and the films are allowed to anneal until a groove appears at the triple junction. The value of $χAB$ is maintained constant at a non-dimensionalized value of $3.5$⁠, while $χAv$ and $χBv$ are assigned values of $2.7$ and $2.75$–$2.8$⁠. For calibration, we also perform a simulation corresponding to $χAB=χAv=χBv=3.5$⁠. However, to retain absolute symmetry, we also require that $κAA=κBB=2κAB$⁠. We assign gradient parameters such that $κAA$ and $κBB$ are less than $2κAB$ across all simulated depositions, so the resulting dihedral angle is expected to be close but slightly below $120°$ when the equality condition for $χ$ is enforced. The dihedral angles that are obtained from curve-fitting analyses along with their representative thermodynamic parameters are listed in Table II.