Surface‑Diffusion‑Driven Morphology Evolution of Pit‑Patterned Si(001) Substrates: A Phase‑Field Study

Background

Lattice‑mismatched heteroepitaxy of semiconductors such as Ge/Si or InGaAs/GaAs frequently leads to the formation of three‑dimensional islands via the Stranski–Krastanow (SK) growth mode. Although self‑assembly can produce quantum dots, the random nucleation inherent to this process often yields a broad size and shape distribution that limits device performance.

Over the past decades, a variety of techniques has been developed to steer heteroepitaxial growth toward ordered structures. Among them, pit‑patterned substrates have emerged as one of the most versatile approaches for achieving both high lateral ordering and precise size control of islands.

Pit patterns are fabricated using top‑down methods such as nanoimprint lithography, electron‑beam lithography combined with reactive‑ion or wet chemical etching, and nanoindentation. These techniques allow the design of ordered pit arrays with nanometer precision, which, under optimized growth conditions, give rise to nearly perfect lateral ordering of islands.

Because the pit shape influences the total surface energy and, consequently, island nucleation, accurate control of pit morphology is essential. At elevated temperatures, capillarity drives the evolution of the pits toward energetically favorable configurations, sometimes leading to complete healing. Therefore, annealing or additional material deposition is routinely employed to stabilize metastable pit shapes. Even after stabilization, further evolution can occur during the heteroepitaxial growth phase.

In this study we describe the evolution of pit‑patterned substrates driven by surface‑energy reduction via surface diffusion. We adopt a phase‑field formulation that captures length and time scales relevant to experiments and has previously been applied to diffusion‑limited kinetics in heteroepitaxial systems. The model is capable of reproducing equilibrium shapes with realistic anisotropic surface energies and has been validated against experimental data.

We focus on pit‑patterned Si(001) surfaces, a system that has been extensively investigated in the literature. The paper is organized as follows: the phase‑field model is introduced in the next section; the smoothing of Si(001) pits and the kinetic pathway toward equilibrium are discussed in Section 3; the application to Ge overgrowth, which induces a pit‑rotation, is presented in Section 4; and the paper concludes with a discussion of the implications of our findings.

Methods

Phase‑Field Model

The phase‑field model introduces a continuous order parameter φ that varies smoothly from φ = 1 (solid) to φ = 0 (vacuum). The evolution of φ is governed by a degenerate Cahn–Hilliard equation that captures surface diffusion while respecting the strong anisotropy of the surface energy. The total free energy functional reads

F = ∫Ω γ(n̂) \,[ ε/2 |∇φ|² + 1/ε B(φ) ] d𝑟 + ∫Ω β/(2ε) [ -ε∇²φ + 1/ε B′(φ) ]² d𝑟,

where Ω is the simulation domain, ε the interface thickness, B(φ)=18φ²(1−φ)² the double‑well potential, and β the Willmore regularization parameter that rounds sharp corners in the strong‑anisotropy regime. The chemical potential μ=δF/δφ drives the evolution through

∂φ/∂t = D ∇·[ M(φ) ∇μ ],

with D the diffusion coefficient and M(φ)=36/ε φ²(1−φ)² the surface mobility. The expression for μ includes the anisotropic surface‑energy density γ(n̂) and the Willmore regularization term, ensuring a physically realistic description of facet evolution.

Anisotropic Surface Energy

We model the surface‑energy density as

γ(n̂) = γ₀ [ 1 − Σᵢ αᵢ (n̂·mᵢ)^{wᵢ} Θ(n̂·mᵢ) ],

where the vectors mᵢ correspond to the crystallographic directions of minimum energy (〈001〉, 〈113〉, 〈110〉, 〈111〉 for silicon). The depth of each minimum αᵢ is derived from the relative surface‑energy values γᵢ reported in the literature, while the width parameter wᵢ controls the angular extent of the facet. For silicon we set γ₀ = 1, α_{〈001〉}=0.15, and w values of 50 (except w_{〈113〉}=100). The Willmore regularization parameter β is chosen as 0.005 to resolve the rounded corners while maintaining computational feasibility.

Initial Morphology and Simulation Setup

We construct an initial pit geometry on a (001) surface defined by a circular depression of radius L and depth H, smoothly connected to the surrounding plane using a Gaussian tail with width σ=L/2. The signed distance d(r) from the interface is inserted into the hyperbolic‑tangent profile φ(r)=½[1−tanh(3 d(r)/ε)] to generate the initial condition. A typical domain size of L = 1 (in arbitrary units) and interface thickness ε = 0.2 are used. The finite‑element toolbox AMDiS is employed with a semi‑implicit time integration scheme and mesh refinement at the interface. Periodic boundary conditions are applied in the lateral directions, and no‑flux conditions are used along the vertical axis. The diffusion coefficient D is set to unity, so the simulation time is reported in arbitrary units.

Results and Discussion

Smoothing of Si(001) Pits

Figure 2 shows the early‑time evolution of pits with aspect ratios R=0.25 and R=0.5. The initially featureless profiles spontaneously develop facets corresponding to the minima of the silicon surface‑energy density. For the shallow pit (R=0.25), a square 〈001〉 base is surrounded by 〈113〉 edges and narrow 〈110〉 triangles. For the deeper pit (R=0.5), additional 〈111〉 facets appear near the bottom, and the 〈110〉 facets widen. These predictions agree quantitatively with experimental observations of pit‑shaped Si(001) surfaces.

Figure 3 traces the long‑time evolution of the deeper pit (R=0.5). The sequence of perspective and top views illustrates the progressive disappearance of steep 〈111〉 facets, the coarsening of adjacent 〈113〉 facets, and the eventual vanishing of 〈110〉 facets, which together yield a flattened surface. The total surface energy decreases monotonically, as shown in Figure 4, and the kinetics follow a common pathway for all R≤0.5. The characteristic time scale for a pit of tens of nanometers at 1100–1200 °C is on the order of hours, consistent with literature values for silicon surface diffusion.

Mimicking the Shape Change due to Ge Overgrowth

To assess the impact of depositing a material with a different equilibrium facet set, we simulated the overgrowth of Ge on a pit patterned Si(001) substrate with R≈0.1. Starting from the morphology at t≈5.0 (Fig. 3), we replaced the silicon surface‑energy density with one that includes a deep minimum along the 〈105〉 direction, characteristic of Ge/Si(001) systems. The width parameter w_{〈105〉} was increased to 500 to capture the shallow 〈105〉 facets.

Figure 5 demonstrates that the surface rapidly reorganizes: 〈105〉 facets nucleate between the remaining 〈113〉 facets, expanding at the expense of the latter. The final pit is bounded exclusively by 〈105〉 facets, and the top view shows a 45° rotation of the pit outline relative to the initial orientation. This rotation has been reported experimentally during Ge deposition and is now shown to arise purely from surface‑energy minimization.

Conclusions

We have developed a continuum phase‑field framework that captures the surface‑diffusion‑driven evolution of pit‑patterned Si(001) substrates. By explicitly incorporating the strong anisotropy of the silicon surface energy, the model reproduces the formation of metastable facets, the eventual flattening of the pit, and the kinetic pathways that lead to these states. Extending the surface‑energy description to include 〈105〉 facets demonstrates that the rotation of pit outlines observed during Ge overgrowth can be explained entirely by surface‑energy reduction. The model is general and can be readily adapted to other substrates by re‑parameterizing the anisotropy, offering a predictive tool for designing pit geometries that guide heteroepitaxial island placement with nanometer precision.