A heterostructure is, in its broadest sense, any stack of solid layers in which at least one material parameter chemical composition, crystal structure, doping, or orientation varies along the growth direction. In modern device engineering, the term is most often reserved for epitaxial stacks in which each layer is a single crystal that is crystallographically registered with its neighbour. The coherent registration of two layers with different equilibrium lattice constants, a_s (substrate) and a_EL (epitaxial layer), is parameterised by the mismatch parameter:

ε₀ = (a_s − a_EL) / a_EL , (1)

which quantifies the in-plane strain that must be accommodated by elastic deformation of one or both layers. When |ε₀| is sufficiently small and the epitaxial layer is sufficiently thin, the strain is stored elastically, and the interface remains coherent. Beyond a critical thickness, first analysed by Matthews and Blakeslee [14] and later refined by People and Bean [4], the stored elastic energy exceeds the line energy of a misfit dislocation, and plastic relaxation sets in [28]. The dislocations nucleated at the interface thread through the epitaxial layer, where they degrade carrier mobility, increase leakage currents, and reduce the radiative efficiency of optoelectronic devices. Recent synchrotron X-ray diffraction and Raman mapping studies on GaN-on-Si epitaxial layers have directly correlated residual stress with threading-dislocation density and have quantified the contribution of superlattice buffers to stress relaxation [29], [30].

Beyond the purely geometrical contribution associated with ε₀, the stress state of an epitaxial stack is also influenced by the difference in thermal expansion coefficients between the substrate and the layer. VPE is carried out at elevated temperatures, typically between 600 °C and 1100 °C depending on the material system. On cooling from the growth temperature to room temperature, the differential thermal contraction introduces an additional thermo-elastic stress that adds to the lattice-mismatch contribution [15], [31]. The total stress at the operating temperature of the device is therefore the superposition of a growth-temperature mismatch component and a thermo-elastic component accumulated during cooling. The two contributions may have opposite signs, so that a layer which is compressively strained at the growth temperature may end up in tension at room temperature, or vice versa, a fact that complicates the design of strain-engineered heterostructures and is exploited, for example, in step-graded AlGaN buffers for GaN-on-SiC technology [7], [15].

The interaction of ionising radiation with semiconductor crystals has been studied intensively since the 1950s and is the subject of several authoritative monographs [32], [33]. When energetic particles electrons, ions, or photons impinge on a crystalline solid, they transfer energy to the lattice through a cascade of elastic and inelastic collisions. If the energy transferred to a lattice atom exceeds the displacement threshold, the atom leaves its site and a Frenkel pair is created, consisting of a vacancy and an interstitial. The primary point defects subsequently migrate, recombine, and aggregate into more complex defect structures: divacancies, diinterstitials, and larger clusters. The kinetics of these processes are governed by the local temperature, the local stress state, and the concentrations of the interacting species. Recent reviews have emphasised the multiscale nature of the problem, linking ab-initio defect energetics, kinetic Monte Carlo simulations, and continuum rate-equation models into a unified description of point-defect kinetics [34–36].

Radiation processing of semiconductors is normally viewed as a source of damage to be avoided or removed by annealing. The present work inverts this perspective: the damaged layer is deliberately introduced into the substrate, and its subsequent evolution during high-temperature VPE is exploited as a stress-relief mechanism. The conceptual basis for this inversion is the elastodiffusive coupling that links the local concentration of point defects to the chemical potential and, through it, to the stress tensor. By engineering the spatial and temporal distribution of the defect population, one may, in principle, redistribute the elastic energy stored in the heterostructure and reduce the peak stress at the critical interface. Ion-implantation-driven defect engineering has already been demonstrated in functional oxide thin films, where it has been used to manipulate carrier recombination dynamics and to enhance charge separation [37], and a closely related principle has been applied to compensate film stress in silicon mirror substrates [38]. The framework developed here extends these ideas to the problem of epitaxial stress reduction.

Closed-form analytical solutions to nonlinear, coupled, time-dependent boundary-value problems of the kind encountered here are rarely available. Numerical simulation is the method of choice when high accuracy is required for a specific geometry, but it obscures the parametric dependences that govern the physics of the problem and provides little insight into the sensitivity of the solution to material parameters. Approximate analytical methods occupy an intermediate position: they sacrifice some accuracy in exchange for explicit functional dependences that can be inspected, differentiated, and used for design optimisation. Among recent alternatives, modern metaheuristic and machine-learning-based approaches have proven useful for the inverse problem of parameter identification in similar nonlinear transport models [21–27]; such approaches are complementary to the forward analytical treatment pursued here.

The MAFC belongs to this intermediate class. It was introduced by Sokolov [17] in the 1950s as a technique for approximating the solutions of nonlinear differential and integral equations [17], [19]. The central idea of the method is to replace, at each iteration order, the unknown function in the nonlinear terms of the equation by an unknown constant, the average of the function over an appropriate domain. The original nonlinear problem is thereby reduced, at each order, to a sequence of linear problems that can be solved in closed form. The unknown constants are then determined self-consistently by an orthogonality condition, which guarantees that the residual of the original equation is, on average, minimised over the chosen domain. Convergence of the iteration is typically rapid for weakly nonlinear problems, and the second-order approximation is usually sufficient for both qualitative analysis and quantitative estimates [20], [39].

The analysis presented in the following sections is restricted to a single heterostructure–reactor configuration: a cylindrical substrate of radius R, hosted in a horizontal VPE reactor with a sloping holder, on which a single epitaxial layer is grown from a gas mixture of reagents diluted in an inert carrier. The radiation processing is assumed to be carried out ex situ, before the substrate is loaded into the reactor, and to produce a known

initial distribution of point defects. The growth temperature is assumed to be uniform in the gas phase and constant in time. These assumptions are idealisations of the actual process, but they capture the essential physics of the problem and allow the analytical machinery of the MAFC to be exercised in full. A schematic of the modelled reactor is shown in Fig. 1.

Longitudinal section showing the gas-flow direction, the heated reactor wall, the substrate, and the sloping holder approximated by straight lines with slope angle φ₁, and b. Cross-section of the holder: the slope angle φ₁ and the reactor radius R define the four-fold angular symmetry exploited throughout the analysis.

The reactor is modelled as a horizontal cylinder of radius R whose axis coincides with the z-axis of a cylindrical coordinate system (r, φ, z). The substrate is held by a sloping holder that subtends an opening half-angle φ₁ on each side of the reactor axis (Fig. 1). The choice of a sloping holder, rather than a flat one, is motivated by the fact that the former produces a more uniform boundary layer over the substrate and a more controllable mass transport of the reagent species to the growth surface, in agreement with recent CFD optimisation studies of horizontal HVPE reactors [11–13]. The symmetry of the configuration with respect to the planes φ = ±φ₁, φ = π ± φ₁ is exploited throughout the analysis to reduce the computational domain. The substrate is treated as an elastic isotropic solid with position-dependent Young modulus E(z), Poisson coefficient ν, density ρ(z), and thermal expansion coefficient β; the epitaxial layer is treated as a thin film whose elastic and transport properties coincide with those of the substrate in the limiting case of lattice matching and are otherwise characterised by an additional mismatch parameter ε₀.

Let I(r, φ, z, t) and V(r, φ, z, t) denote the spatio-temporal concentrations of radiation-induced interstitials and vacancies, with equilibrium concentrations V* (assumed, for simplicity, equal for the two species). The evolution of these concentrations is governed by the coupled diffusion–recombination equations [1], [32], [33].

2 (2)

= ∇ · (DI∇I) − k{I,V}IV − 2 k{ } + GI − ∇ · (D{IS}∇ I),

I,I I S

2 (3)

= ∇ · (DV∇V) − k{I,V}IV − 2 k{ } + GV − ∇ · (D{VS}∇ V),

V,V V S

where D_I(r, φ, z, T) and D_V(r, φ, z, T) are the volumetric diffusion coefficients of interstitials and vacancies, D_{IS} and D_{VS} the corresponding surface diffusion coefficients, k_{I, V}, k_{I, I}, and k_{V, V} the recombination and pair-generation rate parameters, and G_I, G_V the volumetric generation rates of the two species (non-zero only during the radiation-processing stage and zero during the VPE growth stage considered here). The surface-diffusion terms, written with the surface gradient operator ∇_S, account for the enhanced mobility of defects along the free surfaces and interfaces of the heterostructure. The boundary and initial conditions for the systems (2)–(3) read:

I(r, φ, z, 0) = f_I(r, φ, z), V(r, φ, z, 0) = f_V(r, φ, z), (4)

I(0, φ, z, t) < ∞, V(0, φ, z, t) < ∞, (5)

∂I/ ∂r|_{r = R} = 0, ∂V/ ∂r|_{r = R} = 0, (6)

I(r, −φ₁, z, t) = I(r, φ₁, z, t) = I(r, π − φ₁, z, t) = I(r, π + φ₁, z, t), (7)

V(r, −φ₁, z, t) = V(r, φ₁, z, t) = V(r, π − φ₁, z, t) = V(r, π + φ₁, z, t). (8)

Condition (4) specifies the initial defect distributions inherited from the radiation-processing stage; Condition (5) enforces regularity at the reactor axis; Condition (6) imposes zero radial flux at the reactor wall; and Conditions (7)–(8) encode the four-fold angular symmetry imposed by the geometry of the holder. The formation of more complex defects, divacancies Φ_V and diinterstitials Φ_I, is described by an analogous system of equations coupled to Conditions (2)–(3):

∂Φ_I/ ∂t = ∇ · (D_{ΦI} ∇Φ_I) + k_{I, I} I² − k_I Φ_I − ∇ · (D_{ΦIS} ∇_S Φ_I), (9)

∂Φ_V/ ∂t = ∇ · (D_{ΦV} ∇Φ_V) + k_{V, V} V² − k_V Φ_V − ∇ · (D_{ΦVS} ∇_S Φ_V), (10) where D_{ΦI} and D_{ΦV} are the volumetric diffusion coefficients of the complexes, D_{ΦIS} and D_{ΦVS} the corresponding surface-diffusion coefficients, and k_I, k_V the decay-rate parameters of the complexes. The boundary and initial conditions for the complex concentrations mirror those of the point defects, with the additional specification Φ_I(r, φ, z, 0) = Φ_V(r, φ, z, 0) = 0 (the complexes are absent immediately after the radiation-processing step and are generated entirely during the high-temperature stage). The structure of Eqs. (2)–(10) is consistent with the rate-equation models of point-defect kinetics reviewed in [34], [35], in which the recombination of Frenkel pairs is treated as a short-range process, and the formation of complexes is described by second-order reaction kinetics.

The coupling between the defect concentrations and the elastic state of the heterostructure is mediated by the chemical potential μ₁. Following the framework developed in [33], [40], the chemical potential may be written as

μ₁ = E(z) · Ω · σ_{ij} · [u_{ij}(r, φ, z, t) + u_{ji}(r, φ, z, t)] / 2 , (11)

where E(z) is the position-dependent Young modulus, Ω the atomic volume, σ_{ij} the stress tensor, and u_{ij} the deformation tensor. After substitution of the constitutive relation for a linear elastic isotropic solid, Eq. (11) may be rewritten in the explicit form

μ₁ = E(z) · Ω · [ ε₀ + (∂u/ ∂z) − β (T − T_r) ] / (1 − ν) , (12)

where ν is the Poisson coefficient, ε₀ = (a_s − a_EL)/a_EL is the mismatch parameter introduced in Eq. (1), β is the coefficient of linear thermal expansion, T_r is the reference (room) temperature, and u(r, φ, z, t) denotes the relevant component of the displacement vector. The dependence of μ₁ on the local defect concentrations is implicit in Eq. (12) through the strain term ∂u/∂z, which is itself determined by the elastodiffusive equilibrium. The components of the displacement vector, u_r, u_φ, u_z, satisfy the equations of dynamic elasticity [40], [41].

ρ(z) · ∂²u_i/ ∂t² = ∂σ_{ij}/ ∂x_j + ρ(z) · f_i , (13)

where ρ(z) is the position-dependent material density, δ_{ij} the Kronecker symbol, and f_i the body-force density. With the constitutive relation for σ_{ij} substituted into Eq. (13), the system takes the form

ρ(z) · ∂²u_i/ ∂t² = (λ + μ) · ∂/ ∂x_i (∂u_j/. ∂x_j) + μ · ∇²u_i + ρ(z) · f_i , (14)

in which λ and μ are the Lamé coefficients of the heterostructure. The boundary and initial conditions for the Systems (13)–(14) impose zero normal stress on the free surfaces, continuity of displacement and traction at the internal interfaces, and zero initial displacement and velocity. These conditions complete the elastodiffusive block of the model.

The growth of the epitaxial layer from the gas phase is accompanied by heat transfer between the heated walls of the reactor, the flowing gas mixture, and the substrate. The temperature field T(r, φ, z, t) is described by the second Fourier law [31], [42].

ρ · c · (∂T/ ∂t + v · ∇T) = ∇ · (λ_T ∇T) + p(r, φ, z, t) , (15)

where v is the gas-flow velocity, c is the heat capacity of the gas mixture, λ_T is the thermal conductivity, and p(r, φ, z, t) is the volumetric power density deposited into the substrate–holder system. The thermal conductivity of the gas can be expressed in terms of molecular parameters as λ_T = (1/3) · ρ · c_v · ⟨v_m⟩ · ℓ, where ⟨v_m⟩ is the mean molecular speed, ℓ the mean free path, and c_v the specific heat at constant volume [31]. The gas-flow velocity v and the concentration C of the reagent mixture are governed by the Navier–Stokes and the second Fick equations, respectively:

ρ · (∂v/ ∂t + (v · ∇) v) = −∇P + ρ · ν · ∇²v + ρ · g , (16)

∂C/ ∂t + (v · ∇) C = ∇ · (D ∇C) , (17)

in which P is the pressure, ν the kinematic viscosity of the gas, D the diffusion coefficient of the reagent in the carrier, and g the gravitational acceleration. The flow is assumed to be laminar, and the substrate radius R is taken to be much larger than the thicknesses of the diffusion and near-boundary layers, so that the curvature of the boundary layer can be neglected in leading order. This regime is consistent with the laminar-flow assumption employed in the CFD studies of horizontal HVPE reactors reported in [11], [12]. In the limiting-flow regime, when every molecule of the depositing species that reaches the substrate is incorporated into the growing layer, the boundary and initial conditions read:

C(r, φ, −L, t) = C₀, C(r, φ, z, 0) = C₀ · δ(z + L), (18)

T(r, φ, z, 0) = T_r, T(r, φ, −L, t) = T_wall, ∂T/ ∂r|_{r = R} = 0, (19)

v_r(r, φ, −L, t) = 0, v_z(r, φ, −L, t) = V₀, v_z(r, φ, L, t) = V₀, (20)

where C₀ is the inlet concentration of the reagent, T_wall the temperature of the reactor wall, V₀ the axial velocity imposed at the inlet, L the half-length of the reactor, δ the Dirac delta, and the angular symmetry conditions analogous to Eqs. (7)–(8) are imposed on C, T, and on each component of v. The full system of Eqs. (2)–(20), together with the auxiliary relations defining the material parameters, constitutes the complete mathematical model of the heterostructure–reactor system.

The MAFC solutions derived above describe the full spatio-temporal evolution of the defect concentrations during the high-temperature VPE stage. The principal qualitative features of the solutions are summarised in Fig. 2 and may be described as follows. First, the concentrations of interstitials and vacancies decrease monotonically with time as the defects recombine and diffuse out of the damaged region; the rate of decrease is controlled by the recombination parameters k_{I, V}, k_{I, I}, k_{V, V} and by the diffusion coefficients D_I, D_V at the growth temperature, in agreement with the short-range recombination picture of Frenkel-pair kinetics discussed in [35], [36]. Second, the concentrations of the complexes Φ_I and Φ_V initially increase as the point defects aggregate, reach a maximum at a time of order (k_{I, I} I²)^{-1}, and subsequently decrease as the complexes themselves decay. Third, the spatial profiles of I and V broaden with time as diffusion redistributes the defects from their initial localised distribution into the bulk of the substrate. Comparison of the stressed and unstressed profiles (Fig. 2a) reveals that the elastodiffusive coupling has a measurable influence on the vacancy distribution: in the stressed heterostructure, the vacancy concentration near the interface is reduced relative to the unstressed case, reflecting the stress-driven contribution to the chemical potential in Eq. (12). The reduction is modest in absolute terms of the order of a few per cent for typical parameter values but it is systematic and grows with the magnitude of the mismatch parameter ε₀. The implication is that the stress field itself participates in the redistribution of the point defects and, conversely, that the redistribution of the defects feeds back into the stress field through the chemical-potential term. This two-way coupling is precisely the mechanism that the present analytical framework is designed to capture and that purely numerical, fixed-stress simulations tend to overlook.

The displacement-vector component u_z along the growth direction is shown, in normalised form, in Fig. 2b for the unirradiated and irradiated cases. In the unirradiated heterostructure (curve 1), u_z varies monotonically across the interface, reflecting the abrupt change in the lattice constant and the absence of any relaxation mechanism other than elastic deformation. In the irradiated heterostructure (curve 2), the profile of u_z is markedly smoother: the gradient ∂u_z/∂z at the interface is reduced, and the displacement component varies more gradually across the transition region. The smoothing of the displacement profile is the direct geometric signature of the reduction of the mismatch-induced stress, since the stress is proportional to the strain, which is in turn the gradient of the displacement. This trend is consistent with the experimental observation that Ar-ion implantation pretreatment of sapphire substrates reduces the total dislocation density in subsequently grown GaN layers by approximately 65 % [16], and with the broader literature on ion-implantation-driven defect engineering of thin films [37], [38].

The reduction of the interfacial gradient ∂u_z/∂z, and hence of the peak mismatch-induced stress, is the central quantitative result of the present analysis. For typical values of the model parameters substrate thickness of order 0.5 mm, epitaxial-layer thickness of order 1–10 μm, growth temperature of order 1000 K, and initial defect concentrations of order 10¹⁸ cm⁻³ the reduction of the peak stress is estimated to be in the range of 10–30 %, depending on the precise values of the radiation dose, the diffusion coefficients, and the recombination parameters. The dose dependence of the stress reduction is illustrated in Fig. 3: the reduction grows monotonically with the dose at low doses, saturates in the intermediate regime, and eventually declines at very high doses as the residual damage accumulated in the substrate begins to dominate. The existence of an optimal dose is a non-trivial prediction of the model and constitutes one of its main engineering implications.

The MAFC solutions of the gas-dynamic and heat-transfer equations describe the spatio-temporal evolution of the velocity, temperature, and reagent-concentration fields in the reactor; representative radial profiles of the velocity components and a typical axial temperature profile are shown in Fig. 4. The axial velocity component v_z approaches the inlet value V₀ monotonically along the reactor axis, with a boundary-layer correction near the sloping walls of the holder that depends on the slope angle φ₁ and on the kinematic viscosity ν. The radial component v_r is induced by the slope of the holder and is responsible for the redistribution of the reagent across the wafer; its magnitude scales as V₀ · tan(φ₁) near the holder surface. The azimuthal component v_φ, which would be absent in a strictly two-dimensional treatment, is generated by the rotation of the substrate and contributes to the azimuthal uniformity of the deposition. The qualitative shape of these profiles is consistent with the CFD results reported for horizontal HVPE reactors in [11–13].

The temperature field T(r, φ, z, t) approaches a quasi-steady profile on a time scale of order L²/κ_T, where κ_T = λ_T / (ρ c) is the thermal diffusivity of the gas. For typical VPE conditions, this time scale is much shorter than the growth time, so that the temperature distribution can be treated as quasi-steady during most of the growth. The reagent-concentration field C(r, φ, z, t) exhibits a boundary-layer structure near the substrate, with a depletion region whose thickness scales as (D · L / V₀)^{1/2}. The uniformity of C over the substrate surface, which directly determines the uniformity of the deposition rate, is controlled by the slope angle φ₁ and by the rotation rate ω of the substrate, again in qualitative agreement with the parametric studies reported in [11], [12].

Interstitial concentration I Monotonic decay, accelerated by stress-driven diffusion Vacancy concentration V Reduced near interface in stressed layer Displacement component u_z Smoother profile across interface Mismatch-induced stress σ Reduced peak magnitude at interface Interface sharpness Morphological smoothing of the transition region

concentration ≈ 10¹⁸ cm⁻³).

The reduction of the mismatch-induced stress by preliminary radiation processing of the substrate can be understood as the result of three coupled mechanisms that operate simultaneously during the high-temperature VPE stage. The first mechanism is the stress-driven redistribution of the point defects themselves. According to Eqs. (11)–(12), the local chemical potential of a point defect carries an explicit contribution from the stress tensor. In a region of high compressive stress, such as the immediate neighbourhood of a coherent interface with positive mismatch, the chemical potential of an interstitial is raised, while that of a vacancy is lowered. The resulting gradient of the chemical potential drives interstitials away from, and vacancies towards, the high-stress region, thereby redistributing the stored elastic energy over

The second mechanism is the accelerated recombination of the radiation-induced Frenkel pairs at the elevated growth temperature. The high thermal energy available during VPE lowers the activation barriers for vacancy–interstitial recombination, and the local recombination rate k_{I, V} I V increases accordingly, in line with the experimental and atomistic-modelling evidence reviewed in [34–36]. As the point-defect concentrations decay, the local free energy of the damaged region decreases, and the elastic energy stored in the heterostructure is redistributed through the elastodiffusive coupling. The third mechanism is the morphological smoothing of the interface, which results from the enhanced atomic mobility in the defect-rich near-surface region of the substrate. A smoother interface implies a smaller local curvature and, in turn, a smaller peak stress concentration at the boundary between the two layers.

It is instructive to compare the present approach with the more established stress-relief strategies reviewed in the Background section. Graded buffer layers, such as the step-graded AlGaN buffers analysed in [15], reduce the mismatch stress by distributing the lattice-constant mismatch over a thick intermediate layer; they are highly effective but require a substantial increase in the total thickness of the stack and in the growth time. Compliant and porous substrates reduce the effective stiffness of the substrate and thereby the magnitude of the stress transmitted to the epitaxial layer [44]; they require specialised substrate preparation and are sensitive to mechanical damage. Post-growth thermal annealing promotes stress relaxation by dislocation glide; it is widely used but does not address the stress generation mechanism at the growth stage.

The radiation-processing approach analysed in the present work is complementary to these strategies. It does not require any modification of the growth equipment or of the substrate geometry, and it can be combined with the use of graded buffers or post-growth annealing to achieve a cumulative stress reduction. Its principal cost is the additional processing step, the irradiation of the substrate, and the requirement of a controlled radiation dose. The optimisation of the dose, the energy, and the species of the irradiating particles, in conjunction with the growth parameters, is a natural direction for future experimental work. Modern metaheuristic and machine-learning-based optimisation techniques, including the SDMB sampling algorithm for imbalanced data [22] and the farmer-ants optimisation algorithm for discrete problems [23], provide promising tools for the inverse problem of identifying the optimal irradiation protocol from limited experimental data; related approaches have already been used for image-segmentation [24], categorical classification [25], IoT-driven system monitoring [26], and bioprocess optimization [27], and could be directly transferred to the present context.