FEBID Continuum Modelling

• FEBID continuum modelling is a mesoscale framework using coupled rate equations and diffusion models to simulate precursor adsorption, dissociation, and nanostructure growth.
• It integrates detailed electron transport—including bulk and surface energy losses—and secondary electron generation to accurately drive precursor dissociation.
• Extended kinetics models incorporate ligand co-deposition, predicting growth rate peaks and compositional variations based on dwell time and electron flux.

Searching arXiv for the specified FEBID continuum modelling papers to ground the article. Focused electron-beam-induced deposition (FEBID) continuum modelling denotes a class of mesoscale and continuum descriptions in which nanostructure growth under a focused electron beam is represented through coupled surface-kinetic and transport fields rather than by atomistic trajectory-by-trajectory chemistry alone. In the formulation emphasized for FEBID, the evolving precursor surface coverage $\theta(\mathbf{r},t)$ is governed by adsorption, desorption, diffusion, and electron-induced dissociation, while the deposition rate is linked to the dissociation term through fragment yield and sticking probability (Salvat-Pujol et al., 2015). Recent extensions broaden this framework from precursor-only kinetics to coupled precursor–ligand kinetics, enabling prediction not only of growth rate and shape but also of metallic composition as a function of dwell time, electron flux, and ligand residence/dissociation behavior (Jurczyk et al., 29 Sep 2025).

1. Continuum description of FEBID growth

FEBID grows nanostructures by scanning a focused keV electron beam over a substrate bearing adsorbed organometallic precursor molecules (Salvat-Pujol et al., 2015). In continuum models, deposit evolution is described by coupled rate-equation and diffusion–reaction fields for precursor surface coverage and deposition rate (Salvat-Pujol et al., 2015). Within that picture, the electron transport problem enters as a source term: near-surface primary electrons, backscattered electrons, and especially low-energy secondary electrons drive precursor dissociation (Salvat-Pujol et al., 2015).

A standard continuum form for the precursor coverage is

$$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$

where $F$ is the local adsorption flux, $k_{\mathrm{des}}$ the thermal desorption rate, $D$ the surface diffusion coefficient, and $k_{\mathrm{diss}}(\mathbf{r})$ the electron-induced dissociation rate (Salvat-Pujol et al., 2015). The dissociation rate is driven by the local electron energy flux,

$$k_{\mathrm{diss}}(\mathbf{r})=\int_0^{E_{\max}} dE\;\sigma_{\mathrm{diss}}(E)\,\Phi_e(E,\mathbf{r}),$$

with $\sigma_{\mathrm{diss}}(E)$ the dissociation cross section and $\Phi_e(E,\mathbf{r})$ the surface electron flux spectrum including primary, backscattered, and secondary-electron contributions (Salvat-Pujol et al., 2015). The deposition rate is then linked to dissociation via

$$R_{\mathrm{dep}}(\mathbf{r})=Y_{\mathrm{frag}}\,p_{\mathrm{stick}}\;k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$

optionally modified by re-etching or desorption of volatile ligands (Salvat-Pujol et al., 2015).

A closely related single-species continuum formulation writes the precursor kinetics in Langmuir form with one-monolayer maximum coverage,

$$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$0

or, when diffusion is neglected,

$$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$1

with $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$2, $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$3, and $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$4 (Jurczyk et al., 29 Sep 2025). This formalism underlies the classical continuum model used to predict nanoprint shape and growth rate (Jurczyk et al., 29 Sep 2025).

The two papers describe complementary layers of the same continuum hierarchy. One focuses on the electron-transport source term $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$5 and its near-surface corrections (Salvat-Pujol et al., 2015); the other extends the surface-kinetic state variables from intact precursor alone to precursor plus detached ligand, thereby making composition itself a model output (Jurczyk et al., 29 Sep 2025).

2. Electron-transport source terms and the role of surface excitations

A quantitatively correct source term for FEBID requires that electron energy losses and angular and energy distributions be accurately represented in the near-surface region, where surface excitations of the solid’s collective modes strongly affect reflected electron energy-loss spectra (REELS) and secondary-electron emission (Salvat-Pujol et al., 2015). Surface excitations are collective charge oscillations induced by charged projectiles near an interface; they occur on both sides of a planar boundary, in the solid and in the vacuum, because the surface charge density responds to the projectile’s field regardless of which side the projectile occupies (Salvat-Pujol et al., 2015).

For primary energies up to a few keV, surface excitations contribute a sizable fraction of REELS intensity and control secondary-electron generation irrespective of the primary energy (Salvat-Pujol et al., 2015). In REELS at 1 keV, including surface excitations adds low-loss features and intensity: about 20% of the first few tens of eV of loss for Si and 15% for Cu (Salvat-Pujol et al., 2015). Vacuum-side surface losses often make up more than half of the total surface-excitation intensity in REELS (Salvat-Pujol et al., 2015). The probability for surface loss per crossing scales approximately with the surface dwell time, $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$6, remaining significant up to a few keV and enhanced near normal incidence; an in–out asymmetry between incoming and outgoing electrons is most pronounced near normal incidence and around 1 keV (Salvat-Pujol et al., 2015).

These near-interface effects have direct consequences for continuum FEBID source terms. If surface losses are omitted, the near-surface electron flux spectrum $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$7 is distorted in exactly the energy range most relevant to electron-driven dissociation, because secondary electrons with energies $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$8 eV dominate precursor dissociation and surface plasmon decay feeds slow secondary electrons in that range (Salvat-Pujol et al., 2015). This suggests that continuum FEBID models are sensitive not merely to total electron dose but to a detailed, interface-corrected spectral decomposition of the local electron flux.

Surface excitations also affect angular transport indirectly. They increase energy loss, modify effective scattering geometries near the interface, tend to increase backscattering, and can fill elastic minima through small deflections accompanying inelastic events (Salvat-Pujol et al., 2015). A plausible implication is that spatial broadening and proximity effects in FEBID cannot be inferred reliably from bulk transport alone when the growth zone remains within a few nanometers of the vacuum–solid boundary.

3. Dielectric-response formalism for bulk and surface losses

The electron-transport component is described semiclassically through the dielectric formalism using the complex dielectric function $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$9, where $F$0 is the momentum transfer and $F$1 the energy transfer (Salvat-Pujol et al., 2015). For a projectile of charge $F$2 and speed $F$3, the induced field yields the stopping power and the differential inelastic inverse mean free path (DIIMFP) (Salvat-Pujol et al., 2015).

For bulk losses, the DIIMFP per unit path length is

$F$4

with kinematic limits

$F$5

The bulk inelastic inverse mean free path is

$F$6

and the bulk stopping power is

$F$7

The dielectric function can be obtained by extending optical data $F$8 to finite $F$9, for example through Lindhard–Mermin or superpositions of Drude–Lindhard oscillators (Salvat-Pujol et al., 2015).

Near an interface, boundary conditions require a decomposition into bulk and surface contributions. Using the image-charge or extended-pseudomedia method, a commonly used semiclassical form for the surface DIIMFP is

$k_{\mathrm{des}}$0

where the factor $k_{\mathrm{des}}$1 captures the evanescent character of surface modes away from the interface, the loss function $k_{\mathrm{des}}$2 describes surface-mode coupling, and the geometry factor $k_{\mathrm{des}}$3 accounts for incidence or emergence and in–out asymmetry (Salvat-Pujol et al., 2015). In the simplest crossing approximation,

$k_{\mathrm{des}}$4

so the per-crossing probability scales with the surface dwell time (Salvat-Pujol et al., 2015).

Integrating across a near-surface zone of thickness $k_{\mathrm{des}}$5 gives the per-crossing surface loss probability spectrum,

$k_{\mathrm{des}}$6

and the effective inelastic mean free path near the surface becomes

$k_{\mathrm{des}}$7

with

$k_{\mathrm{des}}$8

In Monte Carlo sampling, elastic steps use differential elastic cross sections, while inelastic steps use bulk DIIMFP in the bulk and bulk plus surface DIIMFP within the near-surface zone (Salvat-Pujol et al., 2015).

Within FEBID continuum modelling, this formalism provides the transport kernels from which $k_{\mathrm{des}}$9 is constructed. The central methodological point is that the surface contribution is not a perturbative embellishment but a required term of the effective transport operator in the interfacial region.

4. Secondary-electron emission and dissociation-relevant spectra

Secondary-electron generation couples the energy-loss spectrum of primary and backscattered electrons to the low-energy electrons that actually drive much of FEBID chemistry (Salvat-Pujol et al., 2015). A practical link between inelastic losses and secondary-electron yield is

$D$0

where $D$1 and $D$2 are the average numbers of secondary electrons generated per loss event of energy $D$3 in the bulk and via surface modes, respectively (Salvat-Pujol et al., 2015).

A simple model that reproduces coincidence data launches one secondary electron per energy-loss event with initial energy $D$4,

$D$5

and applies this to all carriers—primary, backscattered, and secondary electrons—with DIIMFPs valid down to approximately 1 eV for transport (Salvat-Pujol et al., 2015). The emitted secondary-electron energy distribution including escape is

$D$6

with escape probability approximated by

$D$7

where $D$8 accounts for transmission through the surface barrier and $D$9 is the effective attenuation length including bulk plus surface losses for the secondary electrons themselves (Salvat-Pujol et al., 2015).

The importance of this formulation for FEBID is empirical as well as theoretical. Realistic secondary-electron spectra require surface losses for primaries, backscattered electrons, and the secondary electrons; otherwise the onset and structure of observed peaks are not reproduced (Salvat-Pujol et al., 2015). Including surface losses yields agreement with measurements on Si, Al, and Ag at 100 eV primary energy, whereas omitting them causes simulated coincidence spectra to fail in reproducing observed onset and structure (Salvat-Pujol et al., 2015). Surface plasmon decay boosts the secondary-electron yield and enhances proximity effects, and an increase of both the secondary-electron yield and the backscatter fraction is expected when surface losses are included (Salvat-Pujol et al., 2015).

Because $k_{\mathrm{diss}}(\mathbf{r})$0 weights $k_{\mathrm{diss}}(\mathbf{r})$1 by the dissociation cross section, and because dissociation often concentrates in the low-energy regime where secondary electrons dominate, the treatment of secondary-electron generation and escape is a constitutive part of continuum FEBID modelling rather than a post-processing step. This suggests that any continuum model aimed at quantitative deposition-rate or composition prediction inherits the quality limits of its low-energy transport treatment.

5. Extended continuum kinetics for ligand co-deposition and composition

A more recent extension of FEBID continuum modelling incorporates detached ligand kinetics explicitly, replacing the traditional assumption that volatile ligand fragments desorb instantaneously (Jurczyk et al., 29 Sep 2025). In the traditional model, gaseous organometallic precursor molecules $k_{\mathrm{diss}}(\mathbf{r})$2 adsorb, diffuse, desorb thermally, and dissociate under electrons; the volatile ligand fragments are assumed to desorb instantly, while the nonvolatile metal-containing fragment remains and contributes to growth (Jurczyk et al., 29 Sep 2025). The extended model introduces a second adsorbed species, the detached ligand $k_{\mathrm{diss}}(\mathbf{r})$3, together with a second electron-induced reaction leading to co-deposited residue (Jurczyk et al., 29 Sep 2025).

The two electron-induced reactions are

$k_{\mathrm{diss}}(\mathbf{r})$4

and

$k_{\mathrm{diss}}(\mathbf{r})$5

where $k_{\mathrm{diss}}(\mathbf{r})$6 denotes the deposited metal fragment or metal-containing fragment, $k_{\mathrm{diss}}(\mathbf{r})$7 the nonvolatile residue from ligand dissociation, and $k_{\mathrm{diss}}(\mathbf{r})$8 a volatile fragment assumed to desorb instantly (Jurczyk et al., 29 Sep 2025).

Neglecting diffusion to preserve analytical tractability, the coupled coverages satisfy

$k_{\mathrm{diss}}(\mathbf{r})$9

$$k_{\mathrm{diss}}(\mathbf{r})=\int_0^{E_{\max}} dE\;\sigma_{\mathrm{diss}}(E)\,\Phi_e(E,\mathbf{r}),$$0

where $$k_{\mathrm{diss}}(\mathbf{r})=\int_0^{E_{\max}} dE\;\sigma_{\mathrm{diss}}(E)\,\Phi_e(E,\mathbf{r}),$$1 accounts for relative footprint area, the blocking term $$k_{\mathrm{diss}}(\mathbf{r})=\int_0^{E_{\max}} dE\;\sigma_{\mathrm{diss}}(E)\,\Phi_e(E,\mathbf{r}),$$2 expresses site competition, and $$k_{\mathrm{diss}}(\mathbf{r})=\int_0^{E_{\max}} dE\;\sigma_{\mathrm{diss}}(E)\,\Phi_e(E,\mathbf{r}),$$3, $$k_{\mathrm{diss}}(\mathbf{r})=\int_0^{E_{\max}} dE\;\sigma_{\mathrm{diss}}(E)\,\Phi_e(E,\mathbf{r}),$$4 define ligand dissociation and desorption rates (Jurczyk et al., 29 Sep 2025). The initial conditions for pulsed exposure are

$$k_{\mathrm{diss}}(\mathbf{r})=\int_0^{E_{\max}} dE\;\sigma_{\mathrm{diss}}(E)\,\Phi_e(E,\mathbf{r}),$$5

The analytical solutions are sums of exponentials with characteristic rates $$k_{\mathrm{diss}}(\mathbf{r})=\int_0^{E_{\max}} dE\;\sigma_{\mathrm{diss}}(E)\,\Phi_e(E,\mathbf{r}),$$6 and steady-state coverages

$$k_{\mathrm{diss}}(\mathbf{r})=\int_0^{E_{\max}} dE\;\sigma_{\mathrm{diss}}(E)\,\Phi_e(E,\mathbf{r}),$$7

$$k_{\mathrm{diss}}(\mathbf{r})=\int_0^{E_{\max}} dE\;\sigma_{\mathrm{diss}}(E)\,\Phi_e(E,\mathbf{r}),$$8

with $$k_{\mathrm{diss}}(\mathbf{r})=\int_0^{E_{\max}} dE\;\sigma_{\mathrm{diss}}(E)\,\Phi_e(E,\mathbf{r}),$$9 (Jurczyk et al., 29 Sep 2025). The stationary coverage ratio is

$\sigma_{\mathrm{diss}}(E)$0

Time-averaged yields per electron over a dwell time $\sigma_{\mathrm{diss}}(E)$1 are defined as

$\sigma_{\mathrm{diss}}(E)$2

and the growth rate is

$\sigma_{\mathrm{diss}}(E)$3

with $\sigma_{\mathrm{diss}}(E)$4 and $\sigma_{\mathrm{diss}}(E)$5 the atomistic volumes per incorporated species (Jurczyk et al., 29 Sep 2025). A central prediction is that the traditional single-species model yields a growth rate that monotonically decreases with dwell time through precursor depletion, whereas ligand co-deposition can generate a peak in total growth rate as ligand-derived residue rises while precursor-derived deposition falls (Jurczyk et al., 29 Sep 2025).

Composition enters through the atomic fractions

$\sigma_{\mathrm{diss}}(E)$6

or, more generally,

$\sigma_{\mathrm{diss}}(E)$7

where $\sigma_{\mathrm{diss}}(E)$8 denotes the metal share per precursor-dissociation event and $\sigma_{\mathrm{diss}}(E)$9 the number of atoms contributing to the residue per ligand-dissociation event (Jurczyk et al., 29 Sep 2025). In steady state,

$\Phi_e(E,\mathbf{r})$0

so the steady-state composition depends solely on the detached ligand kinetics, not on the pristine precursor kinetics (Jurczyk et al., 29 Sep 2025). This is one of the clearest conceptual consequences of the extension: precursor supply controls how much material grows, while ligand desorption and ligand dissociation control what fraction of that material is metallic.

6. Validation, parameterization, and practical implementation

Both the transport-focused and ligand-kinetics-focused formulations are presented as predictive modelling frameworks that require parameterization and validation against electron-spectroscopic and growth data. For the electron-transport part, recommended dielectric parameterizations use optical $\Phi_e(E,\mathbf{r})$1 data extended to finite $\Phi_e(E,\mathbf{r})$2 through Lindhard–Mermin or Drude–Lindhard oscillator fits, with elastic differential cross sections from partial-wave calculations for the substrate and deposit materials (Salvat-Pujol et al., 2015). Bulk inelastic mean free paths at keV energies are typically 10–100 nm depending on material, while the surface contribution is added within about 1–2 nm from the interface; a near-surface zone of $\Phi_e(E,\mathbf{r})$3 Å is used in simulations (Salvat-Pujol et al., 2015). Common substrates and deposits include Si, SiO$\Phi_e(E,\mathbf{r})$4, metals such as Cu, and FEBID deposits from W-, Co-, and Pt-containing precursors on SiO$\Phi_e(E,\mathbf{r})$5 (Salvat-Pujol et al., 2015).

A recommended algorithmic recipe for computing $\Phi_e(E,\mathbf{r})$6 begins by choosing the dielectric function and elastic cross sections, then performing Monte Carlo precomputation over the relevant substrate-plus-deposit geometry, using bulk DIIMFP in the bulk and adding surface DIIMFP within $\Phi_e(E,\mathbf{r})$7 Å of interfaces with in–out asymmetry and both solid- and vacuum-side kernels (Salvat-Pujol et al., 2015). Secondary electrons are generated at each loss event with $\Phi_e(E,\mathbf{r})$8, transported with elastic scattering and DIIMFP down to about 1 eV including surface losses, and the surface electron flux spectrum is tallied on a grid for coupling into $\Phi_e(E,\mathbf{r})$9 (Salvat-Pujol et al., 2015). An analytical or semi-analytical alternative replaces full Monte Carlo by continuous slowing down with an effective mean free path $$R_{\mathrm{dep}}(\mathbf{r})=Y_{\mathrm{frag}}\,p_{\mathrm{stick}}\;k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$0 and angular spread from elastic multiple scattering, calibrated against REELS and secondary-electron yield data (Salvat-Pujol et al., 2015).

Validation of the transport model proceeds by comparing simulated REELS with and without surface excitations to measurements. At 1 keV, including surface losses reproduces low-loss features and adds about 15–20% intensity in the first tens of eV for Si and Cu, while validation of secondary-electron spectra against $$R_{\mathrm{dep}}(\mathbf{r})=Y_{\mathrm{frag}}\,p_{\mathrm{stick}}\;k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$1 coincidence data shows that only models including surface losses for primaries, backscattered electrons, and secondary electrons reproduce experimental peak onsets and structures (Salvat-Pujol et al., 2015).

For the ligand co-deposition model, parameterization uses site densities, molecular fluxes, dissociation cross sections, residence times, and atomistic volumes. An illustrative parameter set uses $$R_{\mathrm{dep}}(\mathbf{r})=Y_{\mathrm{frag}}\,p_{\mathrm{stick}}\;k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$2 m$$R_{\mathrm{dep}}(\mathbf{r})=Y_{\mathrm{frag}}\,p_{\mathrm{stick}}\;k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$3s$$R_{\mathrm{dep}}(\mathbf{r})=Y_{\mathrm{frag}}\,p_{\mathrm{stick}}\;k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$4, $$R_{\mathrm{dep}}(\mathbf{r})=Y_{\mathrm{frag}}\,p_{\mathrm{stick}}\;k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$5 m$$R_{\mathrm{dep}}(\mathbf{r})=Y_{\mathrm{frag}}\,p_{\mathrm{stick}}\;k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$6, $$R_{\mathrm{dep}}(\mathbf{r})=Y_{\mathrm{frag}}\,p_{\mathrm{stick}}\;k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$7 s$$R_{\mathrm{dep}}(\mathbf{r})=Y_{\mathrm{frag}}\,p_{\mathrm{stick}}\;k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$8, $$R_{\mathrm{dep}}(\mathbf{r})=Y_{\mathrm{frag}}\,p_{\mathrm{stick}}\;k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$9 m$$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$00s$$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$01, $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$02 m$$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$03, and $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$04 s, together with ligand parameter sets that produce steady-state metal fractions of approximately $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$05 and $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$06 (Jurczyk et al., 29 Sep 2025). For Cr(CO)$$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$07 at 3 kV and 21 pA, with dwell times from 50 ns to 0.1 s and total exposure 0.1 s, the model reproduces the experimentally observed peak in total growth rate around $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$08–$$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$09 $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$10s, with the metal contribution monotonically decreasing and the ligand-residue contribution rising and then decreasing (Jurczyk et al., 29 Sep 2025). For silver carboxylates under a CASINO-derived flux profile at 20 kV and 0.5–0.6 nA with FWHM approximately 220 nm, the model captures the experimentally observed increase of metal content from center to halo and its sensitivity to the product $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$11 (Jurczyk et al., 29 Sep 2025).

The following table summarizes the roles of the two cited works within FEBID continuum modelling.

Aspect	Transport-centred formulation	Ligand-kinetics extension
Primary state variables	$$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$12, $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$13	$$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$14, $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$15
Main physical emphasis	Bulk and surface energy losses, REELS, SE emission	Detached ligand dissociation/desorption and composition
Core prediction target	Correct near-surface source term for $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$16	Growth rate and metallic fraction versus $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$17, $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$18, $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$19
Key validation	REELS and coincidence SE spectra	Growth-rate peak for Cr(CO)$$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$20 and Ag center/halo composition

Taken together, these formulations indicate a layered modelling program: transport determines the dissociation-driving electron spectrum, while surface kinetics determines how that spectrum is converted into deposition rate, morphology, and composition.

7. Limitations, misconceptions, and open directions

Several pitfalls are identified explicitly. On the transport side, problematic simplifications include using bulk-only DIIMFP everywhere, ignoring the extent of the near-surface zone, extending elastic and inelastic cross sections below their formal validity without calibration, mis-parameterizing $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$21 by relying on optical-only $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$22 without finite-$$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$23 extension, and neglecting multilayer interfaces such as thin adsorbate films or oxide layers whose composite surface modes shift the loss function (Salvat-Pujol et al., 2015). Omitting vacuum-side losses or in–out asymmetry underestimates low-loss intensity and secondary-electron yield, while neglecting the dwell-time scaling with angle distorts spatial distributions near edges and slopes (Salvat-Pujol et al., 2015).

On the kinetic side, the ligand co-deposition model assumes Langmuir monolayer adsorption, no surface diffusion in the analytical treatment, a two-reaction pathway, no autocatalysis, instantaneous desorption of the volatile product, and effective constant parameters $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$24 and $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$25 (Jurczyk et al., 29 Sep 2025). It therefore does not encompass multilayer adsorption relevant for Co$$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$26(CO)$$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$27, autocatalysis seen in Co and Fe systems, or full fragmentation networks leading to carbide or oxide formation (Jurczyk et al., 29 Sep 2025). Electron transport is reduced to a flux $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$28, so detailed secondary-electron spectra and multi-electron effects are not resolved within that kinetic extension itself (Jurczyk et al., 29 Sep 2025).

A recurrent misconception is that FEBID continuum modelling is equivalent to fitting a local dose–growth curve. The cited work on transport indicates instead that the source term must resolve near-surface spectral and angular structure, including both solid- and vacuum-side surface excitations (Salvat-Pujol et al., 2015). A second misconception is that composition follows directly from precursor depletion. The ligand co-deposition framework shows that steady-state metallic composition is governed by detached ligand kinetics, specifically $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$29 and $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$30, rather than by the precursor kinetics alone (Jurczyk et al., 29 Sep 2025).

Open problems are also identified. For transport, these include improved dielectric models with nonlocal and band-structure effects at low energies, surface roughness and three-dimensional morphology effects on surface modes, multilayer substrate–deposit interfaces, improved secondary-electron generation kernels beyond the one-secondary-electron-per-loss approximation, and more accurate low-energy elastic and inelastic cross sections below about 50 eV (Salvat-Pujol et al., 2015). For kinetics, proposed extensions include adding diffusion terms and solving 2D or 3D PDEs numerically, incorporating multilayer adsorption isotherms and replenishment under gas-jet geometries, constructing reaction networks for multi-path fragmentation, coupling to electron-transport solvers for spatially and energetically resolved $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$31, and measuring $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$32 and $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$33 with stronger temperature dependence characterization and deterministic fitting frameworks (Jurczyk et al., 29 Sep 2025).

A plausible implication of reading the two frameworks together is that the next stage of FEBID continuum modelling lies in their tighter integration: dielectric-response-based transport with surface excitations provides the physically correct $$\frac{\partial \theta(\mathbf{r},t)}{\partial t} = F(\mathbf{r},t) - k_{\mathrm{des}}\,\theta(\mathbf{r},t) + D\,\nabla^2\theta(\mathbf{r},t) - k_{\mathrm{diss}}(\mathbf{r})\,\theta(\mathbf{r},t),$$34, while multi-species surface kinetics translates that flux into growth rate, morphology, and composition. In that integrated view, the decisive modelling domain is the few-nanometer near-surface zone where surface excitations govern low-energy electron transport and detached-ligand chemistry governs steady-state purity (Salvat-Pujol et al., 2015).