Modeling the temporal evolution and stability of thin evaporating films for wafer surface processing Open Access

Thin liquid films play a crucial role in many applications. As soon as a surface comes in contact with humid air, a thin water layer will adsorb on it. These layers can be undesirable, e.g., in vacuum chambers, where significant efforts are made to remove the water.¹ On the other hand, thin liquid films are needed in applications such as coating, particle deposition, wafer bonding, and the cooling of electronic devices.^2–5 

The models used to describe thin films are often based on the assumption that the pure liquid is in equilibrium with its own pure vapor.^6–8 In reality, this is mostly not the case. There are some models that assume the thin film to be in contact with a mixture of its own vapor and air,^9,10 and only a few of them include the disjoining pressure.¹¹ However, the disjoining pressure is needed to study thin films down to their equilibrium thicknesses (typically a few monolayers) because it describes the interaction between the solid-liquid and the liquid-vapor interfaces through the ultrathin liquid phase.

The goal of the present paper is to combine literature-known models for the evaporation of thin water films,¹¹ the disjoining pressure,¹² and the spin drying of wafers¹³ to predict the temporal evolution of water layers. As an example, such a model is useful to study water on wafers that are rinsed and spin dried after surface treatments and are going to be bonded, where the amount of water plays a crucial role for the resulting bond strength. In principle, the approach shown in this paper to simulate the whole technological process from rinsing to reaching the equilibrium layer thickness is versatile and can also be applied to models describing other pre-treatments than spin drying of the substrate. The disjoining pressure for water on the substrate, here LiTaO₃ (LTO), is determined from Gibbs free energy calculated by ab initio density functional theory (DFT).

This paper is structured as follows: In Sec. II, the models for evaporation and spin drying are introduced. Furthermore, we show how the disjoining pressure is calculated. The simulation of adsorption of water on LTO as well as the Gibbs free energy by means of ab initio DFT is covered in Sec. III. Section IV summarizes the numerical details of the implementation of the models. The results are shown in Sec. V. We derive an analytical approximation for the time needed to reach equilibrium and compare it with the numerical solution. The spin drying model is used as an input parameter to obtain the initial layer thickness for the evaporation model. As an example, we show how the combined models can be used to minimize the time needed to reach the equilibrium water layer thickness on LTO wafers that are spin dried after rinsing. Finally, Sec. VI summarizes our work.

We start with a brief summary of the equations used to describe evaporation of a thin water film, i.e., the diffusion equation and the mass balance at the water–air interface. Afterward, we will derive an upper bound for the water layer thickness after spin drying. This can be used as an input parameter for the evaporation model. To capture the influence of the solid substrate, the disjoining pressure is incorporated in the boundary condition of the diffusion equation. Its calculation based on ab initio DFT data for the adsorption of water on a LTO surface will be outlined. Finally, we show how to calculate the equilibrium thickness of the water layer.

To study the evaporation of a thin water layer, we use the one-dimensional model proposed by Ajaev et al.¹¹ For the sake of brevity, we will only summarize the set of equations and the corresponding initial and boundary conditions. A sketch of the geometry is shown in Fig. 1: a thin water layer with thickness h on top of a substrate. The solid–liquid interface is located at z = 0. The water is in contact with air, of which the relative humidity H is fixed at a point far away from the layer, i.e., at z = L. The one-dimensional diffusion equation describes the vapor above the water layer,

where c is the vapor concentration, z is the vertical coordinate, t is the time, and D is the diffusivity of water vapor in air. The thickness h of the water layer is determined by the conservation of mass. Utilizing Fick’s first law leads to

where ρ is the density of liquid water. As h ≪ L, we can solve Eq. (1) in the domain 0 ≤ z ≤ L and evaluate the right-hand side of Eq. (2) at z = 0. The advantage of this fixed diffusion domain is that the simulation can be performed without a moving boundary.¹¹ The boundary conditions for the concentration (already formulated for the fixed domain) are

where c_sat is the saturation concentration, V_m is the molar volume of water, R is the universal gas constant, and T is the temperature. The inclusion of the disjoining pressure Π = Π(h(t)) into the evaporation model takes into account the influence of very thin layers (in the order of several monolayers). The calculation of Π is described in detail in Sec. II C. Note that Π → 0 for h → ∞, i.e., both boundary conditions are constant for thick layers. The initial conditions are given as

where h₀ is the initial layer thickness. This initial thickness depends on the pretreatment of the substrate and the model to determine h₀ has to be chosen accordingly. As an example, we will use the spin drying model described in Sec. II B.

The spin drying model proposed by Emslie et al.¹³ is based on the equilibrium of viscous and centrifugal forces. They showed that for a Newtonian fluid layer with initial uniform thickness h_init on top of an infinite rotating disk, the thickness of the layer after spinning time t_spin is given by

The parameter κ = ρω²/(3η) contains the fluid density ρ, the angular velocity ω, and the dynamic viscosity η. The upper bound for the layer thickness after spin drying (in the limit of $hinit2≫4κtspin$⁠) is

Note that this upper bound is independent of h_init. We will use h₀ = h_spin as initial condition for the evaporation model. With this, we get an upper bound for the evaporation time dependent on the spinning parameters.

The disjoining pressure Π = Π(h) is used to capture the influence of the interaction of the solid-liquid and liquid-vapor interfaces through the very thin liquid layer and is defined as the change in Gibbs free energy G with liquid layer thickness h and per unit area at constant cross-sectional area A, temperature T, and pressure p (see, e.g., Ref. 14),

Following Churaev and Sobolev¹² and Starov,¹⁵ Π consists of three components,

1. a molecular component given by

   $Πm=−AH6πh3,$

   (8)

   where A_H is the Hamaker constant;
2. an electrostatic component given by

   $Πe=−ε(Δϕ)28πh2,$

   (9)

   where ɛ is the permittivity of water and Δϕ is the potential difference between the substrate and the water–air interface; and
3. a structural component given by

   $Πs=Kexp−hλ.$

   (10)

   As stated by Starov,¹⁶ there is no way to determine Π_s based on a firm theoretical background. Therefore, only this semiempirical equation can be used to study this component of the disjoining pressure. In contrast to λ, which has a clear physical meaning as decay length or correlation length (typically in the range of 1.2 up to 3.1 Å for water at room temperature¹⁷), the parameter K can be extracted only from fitting to experimental data.

In the present paper, we study pure water on a LTO surface because, in general, deionized water is used for wafer rinsing. Hence, the electrostatic component has to be neglected because there is only a very low concentration of ions in the liquid, which could be responsible for the formation of the potential difference Δϕ. The disjoining pressure used in the following can be written as

To determine the parameters A_H, K, and λ, we equate Eqs. (7) and (11) and solve the resulting differential equation for G. This expression for G is a function of the layer thickness h and will be fitted to ab initio DFT data of LTO surfaces with different amounts of water adsorbed. Details about the ab initio DFT calculations as well as the fitting process are elaborated in Sec. III.

Equilibrium is reached when the boundary conditions (3a) and (3b) are equal. This leads to

Inserting Eq. (11) for the disjoining pressure and rearranging give the final equation, which is numerically solved for h_eq,

The adsorption of water on the LTO surface is studied using spin-polarized ab initio DFT as implemented in the Quantum Espresso (QE) package.¹⁸ The exchange and correlation interactions are modeled by the Perdew–Burke–Ernzerhof (PBE) functional within the generalized gradient approximation.¹⁹ Ultrasoft pseudopotentials²⁰ with a plane-wave basis set are expanded up to a kinetic energy cutoff of 30 Ry for wave functions and 240 Ry for charge density, respectively. To account for dispersion interactions, a nonlocal van der Waals density functional (vdW-DF)²¹ as implemented in QE is used. Atomic structures are optimized using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) quasi-Newton algorithm until the total energy changes less than 0.126 meV and all components of all forces are smaller than 25.7 meV/Å.

LTO is a trigonal crystal belonging to the point group 3m and the space group R3c. The parameters of the LTO unit cell as computed by ab initio DFT are a = b = 5.093 Å and c = 13.705 Å. These parameters agree well with the experimental values of a = b = 5.154 Å and c = 13.783 Å.²² The optimized LTO crystal was then used to create the LTO (0411) surface, as shown in Fig. 2(a). The (0411) orientation of LTO corresponds to a 41.69° X–Y-Cut, which is a preferred wafer cut for several applications.²³ 

To verify that ab initio DFT correctly describes liquid water, we have calculated the formation energy. The liquid water model is cubic and consists of 96 H₂O molecules, corresponding to a density of 0.996 g/cm³. Ab initio MD simulations are performed at 500 K to generate the initial structure of liquid water. The final structure is optimized using the ab initio DFT and BFGS methods. The formation energy of liquid water, E_f, is calculated as

where E_96H2O(l) is the total energy of 96 H₂O molecules in a cubic model and E_H2O(g) is the total energy of a single H₂O molecule. The value of E_f as calculated by ab initio DFT is −0.53 eV. For comparison, the experimental enthalpy of the reaction H₂O(g) → H₂O(l) at 273 K is −0.47 eV.²⁴ 

To model the adsorption of water, some molecules are randomly distributed above the LTO surface. Ab initio DFT calculations are then performed to optimize the adsorption structure and compute the adsorption energy. An example of the adsorption structure of 16 H₂O molecules on the LTO surface is shown in Fig. 2(b). The water adsorption energy E_ads is defined as

where E_surf+nH2O is the total energy of the surface with n adsorbed H₂O molecules on it, E_surf is the total energy of the clean surface, and E_H2O(g) is the total energy of a single water molecule.

Once the adsorption energy E_ads and formation energy E_f are calculated, the Gibbs free energy of adsorption can be computed as

where the chemical potentials of liquid and gaseous water, $μH2O(l)o$ and $μH2O(g)o$⁠, were fitted to values taken from Barin.²⁴ θ denotes the water coverage in monolayers and p_r = p_H2O/p^o is the water partial pressure divided by the total ambient pressure. The derivation of Eq. (16) is explained in the  Appendix A.

To compute the disjoining pressure parameters, we follow the procedure described in Sec. II C. Because Π is defined as derivative of G with respect to the liquid layer thickness, we can identify G with ΔG_ads and Eq. (7) becomes²⁵ 

Here, we used h = h_mlθ, where h_ml ≈ 3.11 Å is the height of one monolayer, and multiplied ΔG_ads with the amount of substance of one monolayer of water to keep the unit of Π. Equating Eqs. (11) and (17) and solving the differential equation lead to

with C being the integration constant. This expression is used to fit the data points calculated with Eq. (16). As an example, Fig. 3 shows the Gibbs free energy for adsorption of water on LTO calculated by ab initio DFT and the corresponding fit for T = 23 °C, p^o = 1 atm, and H = 50%. The resulting fit parameters are summarized in Table I.

We end this section with the calculation of the Hamaker constant for comparison with the fit result. For the interaction of the air–water surface with the water–LTO surface,²⁶ 

with Boltzmann’s constant k_B, the reduced Planck constant ℏ, and the electronic adsorption frequency ν—which has a value of 3 × 10¹⁵ s⁻¹ for water.²⁷ The refractive indices n_i as well as the dielectric constants ɛ_i for air (i = 1), LTO (i = 2), and water (i = 3) are summarized in Table II. For T = 23 °C, we get A_H = −6.68 × 10⁻²⁰ J, which is in the same order of magnitude as the fitted value.

To solve the model consisting of Eqs. (1) and (2) , we use SciPy’s odeint function, which is based on the LSODA solver provided by the FORTRAN library ODEPACK.²⁹ The solver can handle both nonstiff and stiff problems and switches automatically between the suitable methods. For nonstiff problems, Adams method is used, while backward differentiation formula (BDF) method is used in the case of stiff problems. The automatic switching between these methods is controlled by comparing the time step size needed to perform the next step within the requested accuracy (here, we used the default values given by SciPy’s odeint function) as well as the cost per step of each method. In general, the step size is selected based on asymptotic local error analysis. The local error is estimated from the difference between the predicted and corrected values of the solution.³⁰ In addition, the algorithm checks if a change of the order of the method would be beneficial, i.e., leads to greater step sizes. More details about the switching algorithm and the determination of the step size are given in the original publication by Petzold.³¹ The system is discretized in space by second order finite differences. If not stated otherwise, we use L = 1 cm as diffusion domain length and Δz = 10 µm as grid step size. The needed properties of water are calculated based on the formulas given by Kretschmar and Wagner.³² Fitting of the Gibbs free energy ab initio DFT data is done with the lmfit package.³³ 

After briefly showing the outcome of the spin drying model, we discuss in detail the evaporation model. This results in an example showing the optimization of the total time needed to reach equilibrium by evaporation including the spin drying process.

To solve the evaporation model numerically, an input regarding the initial water layer thickness is needed. Based on Eq. (6), the water layer thickness after spin drying for different spinning parameters at T = 23 °C and p^o = 1 atm is shown in Fig. 4. As expected, the thickness decreases with increasing spinning time and angular velocity, respectively. Later, we will use this thickness as initial value in Eq. (4a). The spin drying model is valid as long as the water layer thickness h_spin after spin drying is thick enough so that the disjoining pressure vanishes. As seen in Fig. 4, this is the case because h_spin is in the range of 0.2–4 µm for the given spinning parameters.

Before solving the evaporation model for arbitrary initial thicknesses, an analytical approximation for thick initial layers is presented. In this case, the disjoining pressure vanishes. Then, the diffusion equation has constant boundary conditions,

and the initial concentration profile is given as

In this case, an analytical solution of the diffusion equation can be found (cf., e.g., Ref. 34),

and the layer thickness calculated with Eq. (2) is

As we assumed a thick initial layer, we can expect high values of t_eq, which is the time needed to reach the equilibrium layer thickness h_eq. Assuming t_eq ≫ L²/D and h₀ ≫ h_eq leads to

Note that t_eq only depends on the fluid but not on the substrate in the limit of thick films. In Sec. V B 2, we will solve the evaporation model for thin films under the influence of disjoining pressure and compare the results with approximation (24).

Now, we want to study the transition to ultrathin layers and hence have to include the disjoining pressure. The evaporation process is assumed to be isothermal and the ambient pressure is p^o = 1 atm = const. For T = 23 °C and H = 50%, the parameters of the disjoining pressure (11) can be found in Table I. These parameters are constant over the whole simulation, because T, p^o, and H are constant. The equilibrium height given by solving Eq. (13) numerically is h_eq ≈ 2.0 monolayers = 6.2 Å.

Figure 5(a) shows the result of the evaporation model for different initial heights h₀. Because h approaches h_eq asymptotically, we have to define a condition for reaching the equilibrium state. Therefore, we assume equilibrium if h and h_eq differ no more than 1%. The times t_eq needed to reach the equilibrium are marked by dashed vertical lines. In Fig. 5(b), the corresponding vapor concentrations at the water–air interface are shown. These concentrations are first equal to the saturation concentration c_sat and then approach Hc_sat as soon as the influence of disjoining pressure comes into play.

Figure 6 shows the equilibrium time obtained by solving the evaporation model numerically (solid lines) for different lengths L of the diffusion domain for p^o = 1 atm, T = 23 °C, and H = 50%. The dashed lines mark the corresponding approximation for thick initial layers, Eq. (24). For high values of h₀, the numerical solution always approaches the analytical approximation. As we used t_eq ≫ L²/D for derivation of Eq. (24), higher values of h₀ are needed for larger L to reach the approximation. For very small initial thicknesses, the solid lines overlap, i.e., the results of the evaporation model become independent of L.

We want to close this section with an example: the optimization of reaching equilibrium thickness including the spin drying process regarding total time. For this, we use the results of the spin drying model (cf. Fig. 4) as an input for the evaporation model, i.e., as initial thickness h₀. Note that the spin drying model neglects evaporation during spinning. Figure 7 shows the total time, i.e., spin drying time + evaporation time, needed to reach the equilibrium layer thickness as a function of the spin drying parameters ω and t_spin for p^o = 1 atm, T = 23 °C, and H = 50%. For each value of ω (in steps of 100 rpm), we marked the value of t_spin (in steps of 1 s) for which the total time reaches its minimum with a white cross. As the corresponding initial thicknesses are at least 0.5 µm, we can use the analytical approximation (24) to calculate t_eq (cf. Figs. 4 and 6, here L = 1 cm). Substituting Eq. (6) for h₀ in Eq. (24) with κ = ρω²/(3η) gives the total time

which has a minimum at

marked with a solid white line in Fig. 7. Note that Eq. (26) depends only on the fluid and spinning parameters but not on the substrate properties. This changes if approximation (24) does not hold, e.g., for increasing L. In this case, the evaporation model has to be solved numerically and the influence of the substrate enters via the disjoining pressure. Note that spin drying is only an example of possible pretreatments of the substrate. By substitution of the spin drying model, other technological processes can be described. In this way, different optimizations can be performed to help reduce the process time and hence save costs.

An evaporation model for thin fluid films on solid substrates based on diffusion and mass balance was studied. Inclusion of disjoining pressure allows the description of ultrathin layers down to the equilibrium thickness, which is typically a few monolayers. For thick initial layers, where the disjoining pressure vanishes, an analytical approximation was derived. In this case, evaporation only depends on the fluid and ambient parameters (because all information about the substrate is described with disjoining pressure only).

As an illustrative example, we combined the evaporation model with a spin drying model to study water on LiTaO₃ (LTO). In this way, we simulated the whole temporal evolution of a water layer from wafer rinsing and spinning until reaching the equilibrium layer thickness. Ab initio DFT simulations were performed to compute the adsorption energy of water on LTO. These results were used to calculate the disjoining pressure. For thick initial layers, we found an expression for the spinning time, which leads to a minimal total time for reaching equilibrium. The approach described in this paper can be applied to many other technological processes by substitution of the spin drying model or the solid material to account for different pretreatments and substrates, respectively. In principle, liquids other than water can also be studied.

The authors have no conflicts to disclose.

Max Huber: Conceptualization (supporting); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (equal); Project administration (supporting); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Xiao Hu: Conceptualization (supporting); Formal analysis (supporting); Investigation (supporting); Methodology (equal); Writing – original draft (supporting); Writing – review & editing (equal). Andreas Zienert: Conceptualization (lead); Methodology (supporting); Project administration (supporting); Writing – review & editing (equal). Jörg Schuster: Conceptualization (lead); Project administration (lead); Supervision (equal); Writing – review & editing (equal). Stefan E. Schulz: Conceptualization (lead); Project administration (lead); Supervision (equal); Writing – review & editing (equal).

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

The Gibbs free energy of adsorption of water on a surface is calculated as

where μ_surf+H2O, μ_surf, and μ_H2O(g) are the chemical potentials of the surface with adsorbed water, the clean surface, and gaseous water, respectively. If we assume the gaseous water as an ideal gas and the adsorption process as isothermal, then

Here, $μH2O(g)o(T)$ is the standard chemical potential of gaseous water, which can be directly obtained from a thermodynamic database. p_H2O and p^o are the water partial pressure and the standard pressure, respectively. Additionally,

and

The calculation of the standard chemical potentials $μsurf+H2Oo(T)$ and $μsurfo(T)$ will be described in  Appendix B. θ denotes the water coverage. Inserting Eqs. (A2)–(A4) into (A1), we obtain

where

is the standard Gibbs free energy of adsorption of water on the surface. At equilibrium, ΔG_ads = 0 and, therefore,

with equilibrium constant K_eq. Note that K_eq corresponds to Henry’s adsorption model.

At a given temperature T, the standard chemical potential of species i can be written as

where Δμ_i(T) contains the entropy and enthalpy corrections at T. μ_i(0) is the chemical potential at T = 0 K and is equal to the ab initio DFT total energy E_i. Therefore,

and

To determine the corrections Δμ_surf+H2O(T) and Δμ_surf(T), we assume that they are the same as in the case of liquid water, i.e.,

With this, we get

which can be rewritten using Eq. (B1),

The standard chemical potential $μH2O(l)o(T)$ of liquid water can be obtained from a database. The standard Gibbs free energy of adsorption (A6) can now be written as

We define the formation energy of liquid water as

and the adsorption energy of water on the surface as

With these definitions,

and we end up with the final expression for the Gibbs free energy of adsorption (A5),