Continuum models of focused electron beam induced processing

Review

Introduction to continuum models of focused electron beam induced processing (FEBIP)

Continuum FEBIP models enable the simulation of process rates that govern focused electron beam induced etching (FEBIE), deposition (FEBID) [1-16] and surface functionalization [17] techniques. They are typically used to simulate growth rates and nanostructure geometries as a function of experimental parameters, and to help elucidate the underlying growth mechanisms. Continuum FEBIP models are comprised of differential equations for the rates of change of concentrations of all surface-adsorbed species thought to be involved in the deposition or etch kinetics. The rate equations are functions of time and space, and require specification of the molecular properties of each adsorbate and the electron flux profile(s) at the solid–vacuum interface. Simple continuum FEBIP models can be solved analytically, yielding governing laws delineating the so-called “reaction-rate” and “mass transport” limited process regimes, and resolution scaling laws. Numerical models can account for adsorbate diffusion and enable modeling of processes such as simultaneous FEBIE and FEBID performed using a mixture of precursor gases. Here we provide software that can be used to simulate a wide range of processes reported in the FEBIP literature, and review the underlying continuum FEBIP models (recent general reviews of FEBIP can be found in [4,10,18-21]). We begin with a discussion of the reaction rate limited regime and the most common continuum model input parameters: initial adsorbate coverage, electron flux profile and the gas flux distribution produced by a capillary-style gas injection system. We then cover simple continuum models that are valid in the reaction rate limited regime (where net adsorbate transport via surface diffusion is negligible) and can be used to model FEBIP performed using continuous and pulsed electron beams, physisorbed and chemisorbed precursor molecules, gas mixtures, and multiple reaction products. Finally, we cover a number of models that account for surface diffusion and can be used to model FEBIP in both the reaction and mass transport limited regimes. Throughout, we emphasize the underlying assumptions and limitations inherent to each model.

Before beginning our discussion, we note that the terms “rate”, “concentration” and “flux” are always used to describe quantities with units of reciprocal time [s⁻¹], reciprocal area [Å⁻²] and their product [Å⁻²s⁻¹], respectively. For example, the concentration of adsorbate species ‘a’, N_a, their desorption rate k_a and desorption flux N_ak_a have units of [adsorbates/Ų], [molecules/s] and [molecules/Ų/s], respectively. Frequently used symbols in this review are defined in Table 1. The term “growth” is applied to both FEBIE and FEBID (positive and negative growth rates refer to the growth of deposits and etch pits, respectively). We limit our discussion to FEBIP performed using a stationary, continuous or pulsed, radially symmetric electron beam (i.e., models implemented in cylindrical coordinates). Examples of models of scanned beams can be found in [22-24].

Reaction rate limited growth regime

In order to reach clear conclusions it is often desirable to perform simulations and experiments under simplified conditions where one or more processes are negligible. In the case of FEBIP models, a useful simplification occurs when the net transport of adsorbates via surface diffusion is negligible. This condition can be met by performing FEBIP in the so-called reaction rate limited growth regime (also called the electron-limited growth regime) where the extent of adsorbate depletion caused by electrons is negligible [4,10,19-21]. The strict definition of reaction rate limited growth is that the electron-induced dissociation rate is much smaller than the sum of the adsorption rate and the thermal desorption rate, i.e., using the symbols defined below and in Table 1, for adsorbate species ‘a’, . Conversely, the mass transport limited growth regime (also called the adsorbate-limited growth regime) is defined as . In the mass transport limited regime, the magnitude of the electron flux is sufficiently high to cause significant depletion of precursor adsorbates, and adsorbate diffusion into the area irradiated by electrons makes an important contribution to the growth rate.

In practice, the reaction rate limited regime can be identified simply by measuring or simulating the steady state growth rate as a function of electron flux, as illustrated in Figure 1 (in the figures, numerical values are excluded from axis labels when the plots are used to illustrate general trends that occur under a wide range of FEBIP conditions). Here, the linear portion of the curve corresponds to the reaction rate limited growth regime. In the mass transport limited growth regime, the deposit/pit shapes can provide information on the role of diffusion in the FEBIP process. Inclusion of diffusion in continuum FEBIP models adds a layer of complexity to the modelling and to interpretation of the model outputs. It is necessary only if diffusion is a major contribution to the growth rate. We will discuss models that incorporate diffusion in a self-contained section below, but will limit the majority of our discussion to diffusion-less models.

Initial adsorbate coverage

FEBIP models require specification of the initial concentration of surface-adsorbed precursor molecules , i.e., the steady-state concentration of adsorbates in the absence of electron irradiation. can be found by solving for the difference between the flux of molecules adsorbing from and returning to the gas phase. In the gas phase, the molecule flux is given by:

where P_a is the pressure of the precursor gas for adsorbate ‘a’, m_a is the mass of a gas molecule, k_B is Boltzmann’s constant and T_g is the gas temperature. The simplest case of gas-molecule adsorption onto a substrate surface is that of physisorption, described by a single potential well at the surface as shown in Figure 2. The flux Λ_a of precursor molecules physisorbing to vacant surface sites is given by:

where s_a is the sticking coefficient (in the limit of zero surface coverage), and Θ is the adsorbate coverage, i.e., the fraction of surface sites occupied by physisorbed gas molecules:

where A_a is the area of a surface site occupied by adsorbate ‘a’, and 1/A_a is the maximum possible concentration of species ‘a’.

We note that Equation 2 describes non-activated Langmuir adsorption of a single molecular species ‘a’. The Langmuir model limits the surface coverage to one monolayer (hence the term (1 − Θ) in Equation 2), and can be modified to account for other adsorption behavior such as multilayer adsorption and thermally activated chemisorption. We also note that most FEBIP models assume that s_a is independent of temperature.

The thermal desorption rate k_a of the physisorbed species ‘a’ is given by:

where τ_a is the adsorption time (i.e., adsorbate residence time), κ_a is the desorption attempt frequency, E_a the desorption energy, (i.e., the depth of the potential well shown in Figure 2), and T is the temperature of the substrate surface. The thermal desorption flux is given by N_ak_a, and the adsorbate concentration is found by solving

for N_a(t):

In Equation 6, t = 0 represents the time at which the gas flux F_a is turned on in the model, i.e., the time at which the gas pressure is changed from 0 to P_a. A typical time-evolution of N_a in the absence of electron irradiation is shown in Figure 3. As t → ∞, the surface coverage reaches a steady-state equilibrium value, , which is the initial value that is input into FEBIP models (in which t = 0 represents the time at which the electron flux is turned on):

Electron flux profile

Spatially-resolved FEBIP models require specification of the electron flux profile f(r). Focused electron beams are usually approximated by a Gaussian function:

where Ω is the standard deviation (i.e., full width at half maximum, FWHM = ) of the Gaussian beam. Defocused beams typically have a tophat shape that can be approximated by:

where f₀ is the (maximum) flux at the beam axis (r = 0), β defines the abruptness of the edge of the tophat, and ξ is the beam radius. Tophat beams are useful for quantitative analyses and modeling of experimental data because they are easy to measure and control with a high degree of accuracy [2,3,6,25,26].

A Gaussian and two tophat electron flux profiles are shown in Figure 4. These profiles can be used to approximate those encountered in FEBIP. Actual flux profiles have contributions from primary, backscattered and secondary electrons, each of which has a unique spatial profile and a unique energy distribution [19].

Gas flow from a capillary-style gas injection system (GIS)

FEBIP precursor gases are injected into a specimen chamber using one of two methods. In the first method, the entire vacuum chamber, or a sub-chamber [6] is filled with a precursor gas, as is done in environmental electron microscopy [27-30]. Such vacuum systems can be configured so that the gas pressure is uniform across the substrate surface, and can be measured accurately by conventional pressure gauges. However, in the vast majority of FEBIP setups, a gas injection capillary is used to inject the precursor gas into a chamber that is pumped continuously by a high-vacuum pumping system. The low conduction of the capillary makes it the element that limits the flow rate and shapes the flux profile of such gas injection systems. The capillaries are useful because they simultaneously assure a high local molecule flux at the substrate surface, and a low vacuum chamber background pressure, as required for robust operation of the vacuum system. However, a disadvantage of a capillary-style GIS is that the gas pressure varies dramatically throughout the chamber, and the gas pressure at the substrate surface cannot be measured accurately using conventional pressure gauges. Hence, in this section, we describe a Monte-Carlo simulator developed at Empa (the “GIS simulator”) for calculating gas pressure distributions generated by a capillary-style GIS. The code can be used to calculate the gas molecule flux at the beam impact point on the substrate surface (which is an input parameter to all FEBIP models).

Gas molecule impingement rates on the substrate have lateral distributions that depend on the capillary nozzle geometry, the angle and distance between capillary and substrate, and the molecule flow regime. The GIS simulator code enables the selection of some pre-defined nozzle geometries, capillary inner and outer diameters, angles and distances to the substrate, and either molecular or transient flow regime. By default it maps the impingement distribution on planar substrates, the molecule flux along the tube, and the angular nozzle effusion distribution. The code was validated by capillary injection experiments for both flow regimes [31,32].

Due to the Monte Carlo code implementation, the GIS simulator does not account for gas pressure gradients inside the capillary and above the substrate. The code can be used in an executable version or an editable version (http://www.empa.ch/febipcode) to include specific shadow effects by non-planar substrate geometries (such as deposits grown by FEBID), surface gas phase collisions, or other capillary geometries.

Flow regimes

To make a correct choice of the flow regime for the simulation one needs to determine the mean free path λ between molecule–molecule collisions and the Knudsen number (the ratio of λ to the capillary diameter d):

The Knudsen number specifies the ratio of wall collisions to molecule–molecule collisions. A value of 10 < Kn < ∞ signifies rare flow conditions under which the molecule gas flow distribution is determined predominantly by tube wall collisions (molecular flow regime). In the range 0.1 < Kn < 10 collisions between molecules become more important for shaping the flow (transient flow regime). When Kn < 0.1 the flow becomes viscous and can be treated by continuum models; this regime is not covered by the GIS simulator.

The calculation of the mean free path along a capillary is not trivial as the pressure P along the capillary is not constant. There will be a pressure gradient along capillaries directing net flow from the precursor reservoir (P≈ vapour pressure) to the vacuum chamber (P≤ 10⁻⁴ mbar) [33]. The commonly known relation , where δ is the molecule diameter and other variables are as defined in Equation 1, can only be used if the pressure close to the capillary exit is known.

Mean free paths for a given GIS can be obtained by monitoring the precursor consumption rate Q (molecules per unit time) experimentally and calculating the mean free path using

where R is the ideal gas constant and M is the molar mass [33]. The Knudsen number can then be calculated by inserting Equation 11 into Equation 10. The consumption rate of the precursor can be monitored by measuring its mass or volume change due to evaporation during FEBIP. Alternatively, mass flow controllers can be used to supply a defined flow rate Q (sccm). Figure 5 shows the dependence of the Knudsen number on the pipe diameter for two flow rates of H₂ and H₂O. At a fixed flow rate and capillary diameter, the Knudsen numbers and mean free paths scale with . Values for a few molecules are given in Table 2 together with their vapour pressure and monolayer adsorbate concentration n₀ = 1/A_a = 1.154δ⁻².

To estimate an upper limit for the impingement rate, consider a substrate placed directly in front of a capillary with an exit area of 1 mm² and a flow rate of 1 sccm. This results in a molecule flux F_a ≈ 4.5 × 10¹⁹ molecules cm⁻²s⁻¹ leaving the capillary, or F_a ≈ 5 × 10⁵ monolayers per second (taking 10¹⁴ cm⁻² as monolayer coverage). At a minimum, impingement rates of FEBIP molecules should be greater than the impingement rates of residual gases (at least one monolayer per second at a background pressure of 10⁻⁶ mbar). Of course, the molecule flux will rapidly decrease with distance between the capillary and the substrate.

Wall uptake

The GIS simulator incorporates an uptake factor for capillary wall collisions. Setting this value to zero means that all molecules colliding with the wall will immediately desorb; setting the value to 1 means that all molecules colliding with the wall will adsorb permanently (i.e., only molecules without wall collisions will exit the capillary). In most cases the uptake factor can be set to zero, as molecule condensation on capillary walls is normally avoided through the pressure gradient and keeping the reservoir at a temperature that is lower than or equal to that of the capillary walls. However, some organometallic molecules can spontaneously decompose in contact with the wall material and thus provide a surface for continuous autocatalytic decomposition through successive molecule–wall collisions. (Metal carbonyls are known for autocatalytic decomposition but data quantifying the uptake coefficients is sparse. A change in colour of the capillary points to such a mechanism). Analogous to the atomic layer deposition process [45], a more likely scenario is that at room temperature the molecules decompose partially via chemisorption upon collisions with unoccupied adsorption sites on the wall. This process is self-limiting after a short transient time during which all wall sites are occupied and the wall is passivated for further molecule decomposition.

Nozzle geometries

Two nozzle geometries are incorporated in the GIS simulator [31], and users can implement other geometries. Here we release the source code of version 1.5 containing the geometries shown in Figure 6 for straight capillaries and for bevelled capillaries with access holes. This C++ code is the base of the executable GIS simulator tool downloadable at http://www.empa.ch/febipcode. We would like to note that the code is not professionally commented nor written but has been checked against analytical solutions (for simple geometries) and against experiments. The current version of the code does not consider desorption of molecules from the substrate.

Simulation method

The physics behind the GIS simulator is described in detail in [31,33,46] and we give only a brief summary here. The code is based on the so called test-particle Monte Carlo method, which works in the molecular and transient flow regimes. The molecule flux distribution is obtained as a function of the nozzle geometry and the arrangement of nozzle and substrate. Molecule trajectories are computed consecutively for 10⁶ to 10⁷ molecules. This approach is strictly only valid for molecular flow conditions where molecule trajectories are independent of each other, and collisions occur only with the inner tube wall. However, it was shown experimentally [31] and by comparison with direct simulation Monte Carlo [33] that the test-particle Monte Carlo method with intermolecular collisions gives excellent results in the following way (points 3 and 4 are simplifications, which can be changed by the user if required):

The capillary length must be entered into the simulation. For long capillaries this is an enormous computational effort. A reduced capillary length L_r can be entered instead, still giving accurate results. For molecular flow conditions, the spatial distribution of impinging molecules on the substrate is not significantly altered for lengths L_r > 15d. For transient flow conditions, the reduced length depends on λ with L_r = 3λ to 5λ [31,33].

Simulation examples with the GIS simulator

It is most useful to know how the nozzle geometry and the arrangement of the capillary with respect to the substrate influences the precursor flux at the substrate. The dependence of precursor flow on nozzle geometry (straight nozzle, bevelled nozzle with FEB access hole, closed capillary with FEB access holes), capillary angle, and shadowing by 3D objects have been studied by Friedli et al. [31,33]. More recently, it was shown that certain FEBIP scan strategies also lead to shadowing effects which can cause disruptions in surface flatness [23]. A specific simulation of a gas-flow distribution on a cantilever-based mass sensor enabled the estimation of the residence time of Me₃PtCpMe on SiO₂ [47]. With the released code the reader can include new nozzle geometries or substrate morphologies, for example a conical nozzle geometry as shown in Figure 7. With respect to the straight tube geometry the conical geometries are less practical for FEBIP in the molecular flow regime: Although cones “focus” the exiting molecules, they laterally reduce the maximum impinging flux area, physically cover the maximum flux area, and reduce the maximum impinging flux value slightly as shown by the J/J_tot contours. This behavior can be understood from the cosine-law redistribution of molecule collisions with the inner tube wall dominating the exiting precursor distribution in the molecular flow regime in contrast to the viscous flow experienced daily in water taps. Nonetheless, keeping the inner tube geometry straight and making the outer geometry conical could solve some space restrictions in the scanning electron microscope or allow to bring the nozzle closer to the substrate to increase the molecule impingement rate.

Single gas species FEBIP model: etching or deposition

In the following sections, we outline continuum FEBIP models that can be utilized to simulate a wide range of experiments reported in the FEBIP literature. We start with the simple case of a precursor gas comprised of a single molecular species “a” (i.e., an etch or deposition precursor) that physisorbs to the substrate surface. In this case, the continuum FEBIP model is based on an expression for the rate of change of the adsorbate concentration N_a at each point r on the substrate surface. The expression for ∂N_a/∂t is obtained by extending Equation 5 to account for the dissociation of adsorbates by electrons [16,19]:

This equation is a sum of fluxes representing gas molecule arrival at the surface through adsorption (Λ_a), adsorbate removal through thermal desorption (−N_ak_a), and adsorbate conversion to fragment species through electron induced dissociation (−∂N_α/∂t). In Equation 12, t = 0 represents the time at which electron irradiation is activated in the model, and the initial adsorbate concentration is given by Equation 7.

Most models of FEBIP assume that the adsorbate dissociation flux, ∂N_α/∂t, is proportional to the product of the electron flux f and N_a:

where n_α is the number of fragments α generated per dissociation event, and σ_α is the effective cross-section [3,18] for the generation of fragments that volatilize the substrate (in FEBIE) or deposit onto the substrate (in FEBID).

The steady state adsorbate concentration can be found by solving Equation 12 in the limit t→∞:

Substituting into Equation 13 gives the steady state growth rate . The vertical growth rate ∂h/∂t, which is easily measured experimentally, is proportional to :

where h is the deposit height (or etch pit depth), ι is ±1 for deposition and etching, respectively, and V_D is the volume of a single molecule removed from or added to the substrate in the etch or deposition reaction.

The electron flux f can either be fixed at f₀, or specified by an r-dependent function such as Equation 8 or Equation 9. Setting f to f₀ is justified when simulating the growth rate or the deposit height (etch pit depth) at the beam axis (r = 0). We emphasize that the above equations are valid only in the reaction rate limited growth regime where net transport of adsorbates through diffusion is negligible. We also note that electron-stimulated desorption (ESD) [48-51] is assumed to be negligible. This assumption is typically justified because, for most adsorbates, the ESD cross-section σ_E lies in the range from 10⁻⁷ to 10⁻² Ų [50] (i.e., in general, σ_E << σ_α). However, if necessary, ESD can be incorporated in the above FEBIP model simply by adding the term (−σ_EfN_a) to Equation 12 [52].

Equation 15 can be used to calculate FEBIP growth rates as a function of experimental parameters such as the precursor gas pressure, electron flux and the substrate temperature. For example, Figure 8 shows a set of spatially resolved steady-state vertical growth rates calculated using a Gaussian electron flux profile and substrate temperatures of 300, 320, 340, 360 and 380 K. The growth rate decreases exponentially with temperature due to an increase in the thermal desorption rate k_a (given by Equation 4).

Dimensionless FEBIP models

We now introduce a number of dimensionless parameters, originally employed by Utke et al. [15], which are useful for describing adsorbate kinetics in FEBIP [1] and giving concise scaling laws for the lateral resolution of the FEBIP process.

Irradiative depletion is a dimensionless parameter that quantifies the adsorbate concentration at the beam centre relative to the non-irradiated area. It is proportional to f₀ and can be expressed as [53]:

The effective residence time in the absence of electron irradiation (or outside the irradiated area) τ_out,a is given by:

where τ_a is given by Equation 4, and n₀ is the maximum possible adsorbate concentration (corresponding to one monolayer) in the absence of electron irradiation (i.e., n₀ = 1/A_a).

The steady-state concentration of adsorbates under an electron beam, given by Equation 14, and the steady-state vertical growth rate defined by Equation 15 can now be reformulated as:

In the above, is a dimensionless (normalized) version of the Gaussian or tophat electron electron flux profile given by Equation 8 or Equation 9, respectively, and is given by:

The lateral resolution can be quantified by the dimensionless parameter , expressed as the ratio of the diameter of the deposit (FWHM_D) and the electron beam (FWHM):

In the case of negligible surface diffusion, can be derived analytically yielding the -vs- scaling law of FEBIP resolution for stationary Gaussian [15] and tophat electron beams:

The resolution scaling behavior defined by the above equations is shown in Figure 9. In the limit of very low irradiative depletion the resolution parameter is equal to 1, giving the highest possible FEBIP resolution (i.e., the smallest possible deposit size in the case of FEBID). In other words, the reaction limited (electron-limited) regime gives better lateral resolution than the mass transport limited (adsorbate-limited) regime. Deposits generated by tophat beams exhibit a less pronounced dependency of lateral resolution on irradiative depletion than those generated by Gaussian beams, and the dependence of on is a function of the blurring parameter β (see Equation 9). However, tophat profiles typically have a much larger beam diameter 2ξ compared to the FWHM of a Gaussian.

Pulsed exposure: To illustrate the use of dimensionless parameters in FEBIP models, we consider the case of FEBIP performed using a pulsed, Gaussian electron beam. Specifically, we consider the adsorbate kinetics during an electron pulse, under the assumption that that the surface is replenished entirely between consecutive electron pulses. When the electron beam is turned on, the adsorbate concentration N_a initially decays with time due to electron induced dissociation. The decay rate depends on the effective residence time τ_a,out and the irradiative depletion parameter . The resulting time-evolution of N_a is given by the solution to Equation 12, which can be expressed as:

The corresponding vertical growth rate is proportional to N_a(t) and the electron flux, and is given by:

The lateral size of the growing structure or etch pit will evolve in time, yielding the -vs-t law of pulsed FEBIP resolution, expressed as the lateral resolution , (neglecting surface diffusion) versus the electron beam exposure time t [1]:

In the limit t → 0, . In the steady state (t → ∞) this converges to the scaling law given by Equation 22.

We note that Equation 8 and Equation 9 for the beam profile do not account for lateral contributions of backscattered and forward scattered electrons as well as their generated secondary electrons to the final shape of the deposit or etch pit. The above model also does not account for the change of those contributions on the developing three-dimensional pillar and etch-pit geometry. These simplifications limit the applicability of the resolution scaling laws given by Equation 22, Equation 23, and Equation 26 as well as Equation 79 and Equation 80 to aspect ratios of roughly 1 to 4 for deposits and etch pits. A more detailed discussion on the applicability range can be found in [10]. For high-aspect ratio structures either a full Monte Carlo electron trajectory approach can be chosen [54] or the continuum equations outlined in our present work need to be solved on a curvilinear reference surface, see [13] and [55].

Multiple gas species FEBIP model: Simultaneous etching and deposition

Deposition and etching arising from electron-induced dissociation of multiple adsorbed species is important in a number of experimental processes. Deliberate use of a deposition and etch precursor gas mixture, as shown in Figure 10, is often advantageous, such as in the deposition of high purity materials using a FEBID precursor and an oxygen-containing background gas [56-63], or a mixture of two FEBID precursors [4,8,64-66]. Unintentional deposition of carbonaceous films through electron-activated cross linking of hydrocarbon contaminants is a common problem in many etch processes [6,9,14].

Simultaneous FEBID and FEBIE can be simulated using models such as those described in [11], with adjustment to account for system-specific details such as the surface site competition behavior of the etch and deposit precursor molecules, adsorbate–adsorbate interactions, and whether the etch precursor can volatilize all or only some fraction of the deposited material. A specific example of this is H₂O-mediated FEBIE of deposits formed from organometallic precursors, where only the deposited carbon (and not the metal) is etched.

Here, we illustrate how simultaneous etching and deposition can be incorporated into a continuum FEBIP model for the case of a mixture of H₂O and C₁₀H₂₂ as the precursors for etching and deposition of carbon, respectively [11]. In this model, we assume that both etch and deposition precursor adsorbates (denoted by the subscripts “e” and “d”, respectively) physisorb to the substrate surface (Figure 2) in a competitive process, with the total surface coverage Θ being limited to 1 ML. The fluxes and coverages of the two species are given by the following set of Equation 27–Equation 31:

Here, Equation 27 and Equation 28 are directly analogous to Equation 1, Equation 29 and Equation 30 are analogous to Equation 2, and Equation 31 is equivalent to Equation 3, modified to account for the fact that the adsorbates “e” and “d” physisorb to the surface and compete for the same surface sites. Growth rates can be found by solving differential equations for the rate of change of the concentration of the species “e”, “d” and “D”, where the latter denotes molecules that are deposited as a result of FEBID:

Equation 32 and Equation 33 are analogous to Equation 12, modified to account for the fact that FEBIE can volatilize the adsorbates “d” (i.e., the flux (∂N_ε/∂t)σ_rN_d represents FEBIE of the deposition precursor adsorbates). Equation 34 accounts for the simultaneous deposition of carbon through FEBID (∂N_δ/∂t) and volatilization of the deposit via EBIE [(∂N_ε/∂t)(1 − σ_rN_d)σ_rDN_D ], where (∂N_ε/∂t)σ_rN_d discounts those etch precursor fragments ε that are consumed in volatilizing the deposition precursor adsorbates “d”, and N_D must be capped at 1/A_d so that only molecules in the top monolayer of the deposit are available for volatilization by the fragments ε. The reaction cross-section σ_r (and σ_rD) accounts for the effectiveness of collisions between the etch precursor fragments ε and the adsorbates “d” (or deposited molecules “D”) in contributing to etching [11]. Equation 35 is the vertical growth rate, and is directly analogous to Equation 15.

The fluxes ∂N_ε/∂t and ∂N_δ/∂t represent electron-induced dissociation of the adsorbates “e” and “d”, respectively, and are given by:

where the symbols have the same meanings as in Equation 13. The thermal desorption rates k_e and k_d have the same functional form as Equation 4:

Initial concentrations of the species “e” and “d” (i.e., N_e0 and N_d0, respectively) can be found by solving Equation 32 and Equation 33 for the case of zero electron flux. The resulting time-evolution of surface coverage is more complex than that of the single gas species model shown in Figure 3, due to the nature of the surface site competition effect implemented in Equation 29–Equation 31. This is illustrated in Figure 11a for the case where P_e >> P_d and E_e << E_d. Initially, N_e increases much more rapidly than N_d because P_e >> P_d (and hence F_e >> F_d). However, at times greater than ca. 10⁻⁴ s the increase in N_d causes a corresponding decrease in N_e because the two species compete for the same surface sites and E_d >> E_e (and hence τ_d<< τ_e).

The utility of the simultaneous FEBIE/FEBID model defined by Equation 32, Equation 33 and Equation 34 is illustrated in Figure 11b, which shows a number of deposits simulated using a Gaussian beam (given by Equation 8), with a maximum flux (f₀) of 1 × 10³, 2 × 10³ and 5 × 10³ Å⁻²s⁻¹. At the lowest electron flux (f₀ = 1 × 10³ Å⁻²s⁻¹), the deposit shape is approximately Gaussian because the extent of adsorbate depletion by the beam is low. At the intermediate flux (f₀ = 2 × 10³ Å⁻²s⁻¹), the central region of the deposit contains an indent caused by preferential depletion of the FEBID precursor adsorbates “d”. At the highest electron flux (f₀ = 5 × 10³ Å⁻²s⁻¹), FEBIE dominates in the central region of the deposit due to complete depletion of “d”. The preferential depletion of species “d” is a direct consequence of the slower replenishment rate caused by the fact that P_e >> P_d. The difference between the replenishment rates of “e” and “d” is illustrated in Figure 11a (at times shorter than ca. 10⁻⁵ s).

We emphasize that the changes in deposit geometry seen in Figure 11b are caused by a competition between simultaneous FEBIE and FEBID processes, and are not caused by adsorbate diffusion (which is ignored by Equation 27–Equation 37). This competition and the resulting phenomenon of electron flux controlled switching between etching and deposition have been discussed in detail in [11,14]. The additional effects of adsorbate diffusion can be investigated using a version of the above model that incorporates diffusion (discussed below).

Multiple adsorption states: Thermally activated chemisorption

The FEBID growth rate typically decreases with increasing substrate temperature (as illustrated by Figure 8) due to thermal depopulation of the physisorbed state shown in Figure 2. However, recently it was shown that activated chemisorption can give rise to a more complex dependence of the growth rate on T due to thermal population and depopulation of chemisorbed states [5,67]. The simplest case of a system that can undergo physisorption and activated chemisorption can be described by the potential energy diagram shown in Figure 12 [68], comprised of a physisorbed state (denoted “p”) and a chemisorbed state (denoted “c”) separated by an activation barrier for conversion between these states. The respective physisorption and chemisorption potential wells have depths E_p and E_c, and the activation barrier has height E_conv.

Chemisorption is typically a dissociative process in which the activation barrier represents the energy needed to fragment the precursor molecule. For example, O₂ chemisorbs onto numerous surfaces through the thermally activated reaction O₂ → 2O, where O₂ is in the gas phase or in the physisorbed state, and the two O atoms are in the chemisorbed state. Similarly, XeF₂ can fluorinate many surfaces through decomposition pathways that generate chemisorbed F [17,69-71], and the FEBID precursor tetraethyl orthosilicate (TEOS, Si(OC₂H₅)4) chemisorbs onto SiO₂ as a mixture of (SiO)₂Si(OC₂H₅)₂ and (SiO)Si(OC₂H₅)₃ [72-74].

Chemisorption of gas phase precursor molecules can take place if their thermal energy is sufficiently high to surmount the barrier E_conv − E_p (see Figure 12):

where s_c is the sticking coefficient for activated sticking of gas molecules into the chemisorbed state (in the limit of zero surface coverage), s₀ is the preexponential factor (i.e., the sticking coefficient for the limiting case of E_conv − E_p = 0), and T_g is the temperature of the precursor gas. Similarly, molecules in the physisorbed state can surmount E_conv by gaining thermal energy from the substrate and populate the chemisorbed state at a rate k_conv:

This pathway is more relevant to FEBID than chemisorption from the gas phase (i.e., Equation 40) because FEBID is typically implemented as a cold-wall deposition technique, whereby T_g is dominated by the temperature of the capillary used to deliver the precursor gas, and the substrate temperature T is used to control growth kinetics.

The thermal desorption rates of physisorbed (k_p) and chemisorbed species (k_c) are given by:

Activated chemisorption can be incorporated into the continuum FEBID model by accounting for the above transitions between the gas phase, physisorbed and chemisorbed states by replacing Equation 12 with a pair of equations for the rate of change of concentration of physisorbed molecules (∂N_p/∂t) and chemisorbed molecules (∂N_c/∂t) [5]:

and by replacing Equation 13 with:

where n_p and n_c are the number of fragments generated per dissociation event of the physisorbed and chemisorbed molecules, respectively, assuming that both contribute to FEBID. Hence, the vertical growth rate is given by:

which assumes that the volumes of molecules deposited by the dissociation of physisorbed and chemisorbed adsorbates are both equal to V_D.

The fluxes Λ_p, Λ_c and Λ_conv correspond to physisorption from the gas phase, chemisorption from the gas phase, and transitions from the physisorbed to chemisorbed states, respectively:

where F_a is the flux of precursor molecules “a”. The term (1 − s_c) excludes gas molecules that are trapped in the physisorption potential well, and Θ_p and Θ_c are the coverages of physisorbed and chemisorbed adsorbates, each limited to 1 ML by:

We note that energetic electrons can cause partial decomposition of precursor molecules and hence induce transitions from the physisorbed to the chemisorbed state. This effect is neglected here, but can be incorporated into FEBIP models as in [17] for the case of fluorination caused by decomposition of XeF₂.

The initial concentration of physisorbed and chemisorbed species is found by solving Equation 44 and Equation 45 for the case of zero electron flux. An example of the pressure dependence of and is shown in Figure 13a for a particular set of the energies E_p, E_conv and E_c (corresponding to the chemisorption of TEOS onto SiO₂ [5]), and substrate temperatures of 200, 300 and 350 K. The residence times of physisorbed and chemisorbed adsorbates are given by and , which are governed by the depths of the chemisorption and physisorption potential wells shown in Figure 12. At a substrate temperature of 350 K, chemisorbed states are rapidly populated through the physisorbed state (Equation 41), but depopulate slowly as E_conv << E_c (i.e., k_conv >> k_c). Conversely, physisorbed states are rapidly depopulated because E_p is smaller than both E_conv and E_c. Consequently, saturates at the monolayer limit of 1/A_c at a much lower pressure than saturates at the monolayer limit 1/A_p. At the lower substrate temperatures of 300 and 200 K, k_conv, k_p and k_c are reduced through Equation 41, Equation 42 and Equation 43, yielding the adsorption isotherms shown in Figure 13a.

It is important to note that the rate of chemisorption is always limited by the activation barrier shown in Figure 12. Hence, the rate at which physisorbed states are populated (and replenished during FEBID) is always greater than the rate at which chemisorbed states are populated. This is illustrated by the plots of N_p(t) and N_c(t) shown in Figure 13b where, starting with a depopulated surface, the physisorbed species take ca. 10⁻⁶ s to populate the surface and reach a steady state, whereas it takes the chemisorbed species ca. 1 s to reach a steady state.

Equation 47 can be used to calculate the temperature-dependence of FEBID, as shown in Figure 14 for the case of TEOS [5]. The contribution to growth rate made by physisorbed adsorbates (i.e., fσ_pN_p in Equation 46) decreases with increasing T due to an increase in the thermal desorption rate k_p (and hence a decrease in the adsorbate coverage N_p). The chemisorption component (i.e., fσ_cN_c in Equation 46) is negligible at room temperature because E_conv >> k_BT_g and E_c − E_p >> k_BT. It is characterized by a peak produced by thermal conversion of adsorbates from the physisorbed to the chemisorbed state (i.e., the low-temperature tail of the peak) and desorption from the chemisorbed state (i.e., the high-temperature tail), respectively. The general temperature dependence seen in Figure 14 (i.e., a decrease in the physisorption component followed by a peak caused by chemisorption) exists only if k_conv << k_p << k_c [5]. A practical benefit of performing FEBID at elevated temperatures (whereby, in this case, T ≥ 350 K) is that the deposits have much higher purity than those fabricated at room temperature [5].

Role of multiple reaction products in electron beam induced etching

Etch processes such as XeF₂-mediated FEBIE of SiO₂ [12] and NF₃ FEBIE of Si [2] proceed through multiple chemical pathways that involve a number of reaction products. Here, we use the case of NF₃ FEBIE of Si [2] to illustrate how such processes can be incorporated in continuum FEBIE models.

FEBIE of Si using NF₃ as the precursor gas involves the reaction products SiF, SiF₂, SiF₃, and SiF₄, each of which can be dissociated by electrons, or can desorb from the substrate surface. For example, SiF₂ can either gain F to form SiF₃, dissociate to form SiF and F, or it can desorb (and hence give rise to etching by removing a Si atom):

Multi-step FEBIE reactions can be modeled using a set of differential rate equations that account for each molecular species at the substrate surface, which in this case are the NF₃ precursor adsorbates (denoted by “a”), F radicals (α), and the reaction products SiF, SiF₂, SiF₃, and SiF₄:

The symbols in Equation 53 and Equation 54 represent the same quantities as the corresponding symbols in Equation 12 and Equation 13 (see Table 1). In Equation 55, N_n (n = 1–4) represents the concentration of SiF_n molecules, σ_r is the cross-section for the electron-induced scission of the Si–F bond, and k_n is the desorption rate of SiF_n:

where E_n is the binding energy of SiF_n. N_η is the concentration of surface sites at which F can bond to a Si atom, and Θ_η is the coverage of surface sites occupied by F (see below in Figure 15):

The term (1/n)A_Si limits the concentration of fluorinated Si atoms to one monolayer, m is an integer with lower and upper limits of 0 and 3 because an unfluorinated Si atom (designated by m = 0) can react with F to form SiF, whereas SiF₄ (designated by m = 4) cannot gain an additional fluorine atom. The integer n is bound by 1 and 4 because the total coverage of sites occupied by F must account for SiF, SiF₂, SiF₃ and SiF₄ species. N₀ (the concentration of unfluorinated Si atoms at the surface) is given by:

where A_Si is the area of a single Si surface site, and N_Si is the concentration of fluorinated Si sites:

The vertical etch rate is governed by the desorption rate of SiF_n molecules:

Models such as the above (Equation 53–Equation 61) are typically based on a number of simplifying assumptions in addition to those made in simple FEBID models. This particular model assumes that all available F radicals generated by electron-induced dissociation of NF₃ are converted to SiF_n by reacting with SiF_n−1, and that any surplus fluorine atoms desorb as shown in Figure 15 (the desorption is assumed instant, hence the absence of a desorption flux in Equation 54). The probability of a reaction between F and each species SiF_n−1 is assumed equal and the total F coverage is limited to 1 ML, hence the term (N_n−1/N_η)(1 − Θ_η) in Equation 55 that governs the partitioning of the available F radicals illustrated by the flow chart shown in Figure 15. It is also assumed that σ_r is independent of n. Hence, the flux σ_rfN_n+1 represents the creation of SiF_n through electron-induced dissociation of SiF_n+1, and −N_nσ_rf the consumption of SiF_n through electron-induced dissociation of SiF_n. Finally, it is assumed that electrons are not consumed in any of the reactions, and that electron-stimulated desorption of all species is negligible.

The most significant consequence of the model defined by Equation 53–Equation 61 is that etching is inhibited at high electron fluxes, as seen in Figure 16a, due to electron-induced dissociation of SiF_n. Specifically, this effect dominates when the dissociation rate of species n is much greater greater than the corresponding thermal desorption rate (i.e., σ_rf >> k_n). The simulations presented here were performed using model input parameters appropriate for cryogenic NF₃-mediated FEBIP of Si performed at a substrate temperature of 100 K [2]. The model can be used to calculate a number of quantities such as the concentrations of NF₃, SiF, SiF₂, SiF₃ and SiF₄ molecules shown in Figure 16b, and the concentrations of unfluorinated (N₀) and fluorinated (N_Si) Si sites shown in Figure 16c.

Special cases

The models outlined above can be used (either directly or in modified form) to simulate most processes reported in the FEBIP literature. Here, we summarize a few special cases for which the models require modification. We include these for completeness, and refer the reader to the cited papers for detailed descriptions of the changes that must be made to the above FEBIP models.

Adsorbate depletion in high aspect ratio pits [7]: The replenishment rate of precursor molecules consumed in FEBIE can be limited by the gas flow conductance of a growing etch pit. The replenishment rate affects the adsorbate concentration N which, in turn, affects the etch rate through equations such as Equation 15. Etch-pit conductance can be the most significant process limiting the FEBIE rate when fabricating high aspect ratio pits. As an etch pit grows during FEBIE, the conductance decreases, causing the etch rate to decrease with time as discussed in detail in [7].

Dynamic surface site activation during FEBIE [3]: Surface site activation caused by electron beam restructuring of the substrate (i.e., beam damage) can give rise to FEBIE of materials that can not be etched in their unmodified state. It can also alter etch kinetics, and cause the etch rate to increase with electron beam irradiation time.

Electron beam induced surface functionalization and spontaneous decomposition of precursor molecules at the substrate surface: These two processes have been modeled for the case of XeF₂ [17] which can fragment through a dissociative chemisorption pathway, leading to fluorination of many surfaces [69-71,75]. The model in [17] is a variant of the above model of thermally activated chemisorption defined by Equation 40–Equation 52.

Electron beam dwell time as a control parameter of the composition of materials deposited using a mixture of two precursor gases [8]: This situation, depicted in Figure 10, can be modeled by solving rate equations for two FEBID precursor adsorbates, denoted by “A” and “B”, that compete for adsorption sites on the surface. The adsorbates are dissociated by electrons, producing non-volatile reactions products D_A and D_B (i.e., the deposit), and volatile reaction products V_A and V_B:

The evolution of the surface densities n_A and n_B of the two adsorbates is described by the following set of inhomogeneous first-order differential equations [8]:

The generic solution of this set as a function of electron-beam exposure (dwell time) is given by:

with the constants n_dA,B, k_d, Δn_A,B, κ defined in [8]. The dissociation yields Y_A,B are then obtained by integrating over the electron beam dwell time:

A graphical representation of Equation 66 and Equation 67 is shown in Figure 17. As the molecule fluxes, residence times, and dissociation cross-sections of molecules “A” and “B” are very likely different from each other, it can be seen that the composition (given by the magnitude Z on the right hand axis) of the deposits can be tuned by changing the electron beam dwell time per pixel. Equation 64 and Equation 65 hold for micrometer-sized deposits. For small-scale deposits surface diffusion can become important and different yields and compositions are obtained as function of exposure time, see Gabureac et al. [76,77]. Including surface diffusion into the multi-adsorbate model is straightforward, however, is at the expense of having a closed analytical solution.

Diffusion

All of the above models are, strictly speaking, valid only in the reaction rate limited growth regime. They can be used to identify conditions under which depletion becomes significant, causing the system to transition into the mass transport limited regime (e.g., under intermediate and high electron flux conditions in Figure 11b, and the high electron flux portions of the curves in Figure 16). However, when the extent of depletion is significant (or, more precisely, when the net transport of adsorbates through diffusion is significant), the models must be modified to account for the diffusion of all mobile species at the surface. In the case of FEBID or FEBIE performed using a single species of physisorbed precursor molecules “a”, this is achieved simply by adding a diffusion term to Equation 12 [19]:

D_a is the diffusion coefficient, given by:

where is the energy barrier for diffusion and is the preexponential factor (i.e., the diffusivity in the limit ).

The above approach for incorporating diffusion can also be applied to more complex FEBIP models, as we illustrate below. We note, however, that all of the continuum models discussed here assume a flat substrate surface. That is, the models track the time-evolution of deposits and etch pits made by FEBIP, but do not account for the effects of the resulting changes in surface geometry on adsorbate kinetics (e.g., adsorption to and diffusion along the evolving pillar sidewalls or etch pit sidewalls).

Gas mixtures

Equation 32–Equation 35 are used to model a gas mixture comprised of an etch precursor “e” and a deposition precursor “d” that simultaneously etch and deposit a material such as carbon. To account for the diffusion of the adsorbates “e” and “d” at the substrate surface, Equation 32 and Equation 33 are replaced by [11]:

where D_e and D_d are the respective diffusion coefficients which have the same functional form as Equation 69:

Equation 34 and Equation 35 remain unchanged since localized FEBIP requires the substrate temperature to be sufficiently low for diffusion of the deposited species “D” (e.g., carbon) to be negligible.

Figure 18 shows deposit geometries simulated using Equation 70 and Equation 71, a Gaussian electron-beam profile and substrate temperatures of 285, 290 and 295 K. The dip in the center of each deposit is caused by the fact that etching dominates near the beam axis where the deposition precursor is preferentially depleted, while deposition dominates in the adjacent regions of low electron flux. The deposit geometry changes with temperature because an increase in T causes: (i) a decrease in N_e and N_d through thermal desorption (Equation 38 and Equation 39), and (ii) an increase in N_e and N_d through diffusion (Equation 72 and Equation 73). The dramatic change in geometry with increasing T shown in Figure 18 is largely due to reduced depletion of adsorbates “d” near the beam axis due to an increase in the rate of diffusion (Equation 73).

Thermally activated chemisorption

The model of activated chemisorption described by the potential energy diagram shown in Figure 12 can be extended to account for diffusion by adding diffusion terms to Equation 44 and Equation 45:

where D_p and D_c are the diffusion coefficients for physisorbed and chemisorbed adsorbates, respectively.

Figure 19 shows FEBID rates simulated using Equation 74 and Equation 75, a Gaussian electron-beam profile and substrate temperatures of 300 K and 500 K. At 300 and 500 K, FEBID is dominated by electron-induced dissociation of physisorbed and chemisorbed gas molecules, respectively (see Figure 14), because of preferential population of the two respective states shown in Figure 12.

Dimensionless FEBIP models

In this section, we incorporate surface diffusion into the dimensionless FEBIP model, introduced above through Equation 16–Equation 26.

The magnitude of adsorbate surface diffusion with respect to the FWHM size of the beam can be correlated by introducing the dimensionless surface diffusion replenishment parameter:

Substituting the irradiative depletion from Equation 16 and the surface diffusion replenishment to the adsorbate rate equation, it becomes solely dependent on those two parameters and can be solved numerically [79]. Very small values of the surface diffusion replenishment parameter (e.g., = 0.052) have a negligible effect on the lateral growth rate, as shown by row 1 of Figure 20 and Figure 21. However, large values (e.g., = 0.52) can alter the deposit geometries significantly, particularly in the case of a tophat electron beam (see rows 2 and 3 of Figure 20 and Figure 21). If the contribution of surface diffusion is very high ( >> 0.1), adsorbate replenishment through this process dominates the growth kinetics. For example, when = 100, the deposit height is increased by a factor of 10 relative to the situation where adsorbate replenishment occurs only through adsorption from the gas phase.

The resolution parameter for a stationary Gaussian beam can be obtained from the deposit shape, using a numerical solution of Equation 68, and the expression for the growth rate given by Equation 15. The numerical result is very well approximated by the -vs- scaling law of FEBIP resolution formulated in [1]:

The above formula shows that even a relatively small surface diffusion contribution can lead to a decreasing in the deposit and etch pit size, and thus an improved lateral resolution. Once surface diffusion overcomes irradiative depletion, the deposit or etch pit diameter will approach that of the beam ( = 1), as discussed in [1].

Pulsed exposure: Equation 68 can be solved numerically, and the corresponding general scaling law of FEBIP resolution, including irradiative depletion, surface diffusion, and exposure dwell time, can be approximated by [1]:

The dependencies of deposit geometry on the electron beam exposure time t, irradiative depletion and surface diffusion path are shown in Figure 20 and Figure 21 for Gaussian and tophat beams, respectively. At short exposure times, depletion is negligible and the deposit shapes reflect the electron flux profiles. At long exposure times, the deposit geometry is determined by the extent of adsorbate depletion, and the relative contributions to adsorbate replenishment through diffusion and adsorption from the gas phase [1].

Conclusion

In summary, we have reviewed continuum modeling techniques developed by the authors that can be used to simulate most processes reported in the FEBIP literature. Accompanying this article, we have released a software implementing most of these techniques at http://www.empa.ch/febipcode.

The software release consists of the following:

FEBIP models implemented in Mathematica notebooks:

1. Single physisorbed gas species initial coverage calculator, used to generate Figure 3.

2. Single physisorbed gas species FEBIP model based on Equation 12. This model was used to generate Figure 8.

3. Gas mixture FEBIP model based on Equation 32–Equation 39 and used to generate Figure 11a. For this simulation, initial coverages can be calculated by setting the electron flux to zero.

4. Gas mixture FEBIP model that incorporates diffusion, based on Equation 70–Equation 73 and used to generate Figure 18 and Figure 11b (in the latter, the diffusion coefficients were made negligible).

5. Thermally activated chemisorption initial coverage calculator used to generate Figure 13.

6. Thermally activated chemisorption, based on Equation 44–Equation 52 and used to generate Figure 14.

7. Thermally activated chemisorption with diffusion, based on Equation 74–Equation 77 and used to generate Figure 19.

8. EBIE with multiple reaction products, based on Equation 53–Equation 61 and used to generate Figure 16. Initial coverages can be calculated by setting the electron flux to zero.

Gaussian (Equation 8) or tophat (Equation 9) focused electron beams can be used by selecting the appropriate equations in each notebook.