On approach to increase density of field-effect heterotransistors in the framework of three-stage differential amplifier. Influence of mismatch-induced stress and porosity of materials on technological process

In the present time several actual problems of the solid state electronics (such as increasing of performance, reliability and density of elements of integrated circuits: diodes, field-effect and bipolar transistors) are intensively solving [1– 5]. To increase the performance of these devices it is attracted an interest determination of materials with higher values of charge carriers mobility [6– 9]. One way to decrease dimensions of elements of integrated circuits is manufacturing them in thin film heterostructures [2– 4, 11]. In this case it is possible to use inhomogeneity of heterostructure and necessary optimization of doping of electronic materials [11] and development of epitaxial technology to improve these materials (including analysis of mismatch induced stress) [12– 14]. An alternative approaches to increase dimensions of integrated circuits are using of laser and microwave types of annealing [15, 16].

In the framework of the paper we introduce an approach to optimize manufacture of field-effect heterotransistors. The approach gives a possibility to decrease their dimensions with increasing their density in the framework of a three-stage differential amplifier [17]. We also consider possibility to decrease mismatch-induced stress to decrease quantity of defects, generated due to the stress. In this paper we consider a heterostructure, which consist of a substrate and an epitaxial layer (see Figure 1).

We also consider a buffer layer between the substrate and the epitaxial layer. The epitaxial layer includes into itself several sections, which were manufactured by using another materials. These sections have been doped by diffusion or ion implantation to manufacture the required types of conductivity (p or n). These areas became sources, drains and gates (see Figure 1). After this doping it is required annealing of dopant and/or radiation defects. Main aim of the present paper is analysis of redistribution of dopant and radiation defects to determine conditions, which correspond to decreasing of elements of the considered amplifier and at the same time to increase their density. At the same time we consider a possibility to decrease mismatch-induced stress.

1. Method of solution

To solve our aim we determine and analyzed spatio-temporal distribution of concentration of dopant in the considered heterostructure. We determine the distribution by solving the second Fick’s law in the following form [1, 18– 21].

\[\begin{align} \label{GrindEQ__1_} \frac{\partial C(x,y,z,t)}{\partial t} &= \frac{\partial}{\partial x} \left[ D \frac{\partial C(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ D \frac{\partial C(x,y,z,t)}{\partial y} \right] + \frac{\partial}{\partial z} \left[ D \frac{\partial C(x,y,z,t)}{\partial z} \right] \nonumber\\ &\quad + \Omega \frac{\partial}{\partial x} \left[ \frac{D_S}{kT} \nabla_S \mu_1(x,y,z,t) \int_0^{L_z} C(x,y,W,t) \, dW \right] \nonumber\\ &\quad + \Omega \frac{\partial}{\partial y} \left[ \frac{D_S}{kT} \nabla_S \mu_1(x,y,z,t) \int_0^{L_z} C(x,y,W,t) \, dW \right] \nonumber\\ &\quad + \frac{\partial}{\partial x} \left[ \frac{D_{CS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ \frac{D_{CS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial y} \right] + \frac{\partial}{\partial z} \left[ \frac{D_{CS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial z} \right]\nonumber\\ &\quad +\frac{\partial }{\partial \, \, x} \left[\frac{D_{C\, S} }{\bar{V}\, k\, T} \frac{\partial \, \, \mu _{2} \left(x,y,z,t\right)}{\partial \, \, x} \right]+\frac{\partial }{\partial \, \, y} \left[\frac{D_{C\, S} }{\bar{V}\, k\, T} \frac{\partial \, \, \mu _{2} \left(x,y,z,t\right)}{\partial \, \, y} \right]+\frac{\partial }{\partial \, \, z} \left[\frac{D_{C\, S} }{\bar{V}\, k\, T} \frac{\partial \, \, \mu _{2} \left(x,y,z,t\right)}{\partial \, \, z} \right], \end{align}\tag{1}\] with boundary and initial conditions \[\begin{aligned} \left. \frac{\partial C(x,y,z,t)}{\partial x} \right|_{x=0} &=& 0, \\ \left. \frac{\partial C(x,y,z,t)}{\partial x} \right|_{x=L_x} &=& 0,\\ \left. \frac{\partial C(x,y,z,t)}{\partial y} \right|_{y=0} &=& 0, \\ C(x,y,z,0) &=& f_C(x,y,z), \\ \left. \frac{\partial C(x,y,z,t)}{\partial y} \right|_{y=L_y} &=& 0, \\ \left. \frac{\partial C(x,y,z,t)}{\partial z} \right|_{z=0} 7&=& 0, \\ \left. \frac{\partial C(x,y,z,t)}{\partial z} \right|_{z=L_z} &=& 0. \end{aligned}\]

Here C(x,y,z,t) is the spatio-temporal distribution of concentration of dopant; \(\Omega\) is the atomic volume of dopant; \(\nabla\)\({}_{s}\) is the symbol of surficial gradient; \(\int _{0}^{L_{z} }C\, \left(x,y,z,t\right)\, d\, z\) is the surficial concentration of dopant on interface between layers of heterostructure (in this situation we assume, that Z-axis is perpendicular to interface between layers of heterostructure); \(\mu\)\({}_{1}\)(x,y,z,t) and \(\mu\)\({}_{2}\)(x,y,z,t) are the chemical potential due to the presence of mismatch-induced stress and porosity of material; D and D\({}_{S}\) are the coefficients of volumetric and surficial diffusions. Values of dopant diffusions coefficients depends on properties of materials of heterostructure, speed of heating and cooling of materials during annealing and spatio-temporal distribution of concentration of dopant. Dependences of dopant diffusions coefficients on parameters could be approximated by the following relations [22– 24]. \[D_{C} =D_{L} \left(x,y,z,T\right)\, \left[1+\xi \frac{C^{\gamma } \left(x,y,z,t\right)}{P^{\gamma } \left(x,y,z,T\right)} \right]\, \left[1+\varsigma _{1} \frac{V\, \left(x,y,z,t\right)}{V^{*} } +\varsigma _{2} \frac{V^{2} \left(x,y,z,t\right)}{\left(V^{*} \right)^{2} } \right],\nonumber\ \] \[D_{S} =D_{S\, L} \left(x,y,z,T\right)\, \left[1+\xi _{S} \frac{C^{\gamma } \left(x,y,z,t\right)}{P^{\gamma } \left(x,y,z,T\right)} \right]\, \left[1+\varsigma _{1} \frac{V\, \left(x,y,z,t\right)}{V^{*} } +\varsigma _{2} \frac{V^{2} \left(x,y,z,t\right)}{\left(V^{*} \right)^{2} } \right]. \tag{2}\]

Here D\({}_{L}\)\({}_{\ }\)(x,y,z,T) and D\({}_{LS}\)\({}_{\ }\)(x,y,z,T) are the spatial (due to accounting all layers of heterostruicture) and temperature (due to Arrhenius law) dependences of dopant diffusion coefficients; T is the temperature of annealing; P\({}_{\ }\)(x,y,z,T) is the limit of solubility of dopant; parameter \(\gamma\) depends on properties of materials and could be integer in the following interval \(\gamma\)\({}_{\ }\)\(\in\)[1,3] [22]; V\({}_{\ }\)(x,y,z,t) is the spatio-temporal distribution of concentration of radiation vacancies; V\({}^{*}\) is the equilibrium distribution of vacancies. Concentrational dependence of dopant diffusion coefficient has been described in details in [22]. Spatio-temporal distributions of concentration of point radiation defects have been determined by solving the following system of equations [18– 21, 23, 24].

\[\begin{align} \label{GrindEQ__3_} \frac{\partial I(x,y,z,t)}{\partial t} &= \frac{\partial}{\partial x} \left[ D_I(x,y,z,T) \frac{\partial I(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ D_I(x,y,z,T) \frac{\partial I(x,y,z,t)}{\partial y} \right] \nonumber\\ &\quad + \frac{\partial}{\partial z} \left[ D_I(x,y,z,T) \frac{\partial I(x,y,z,t)}{\partial z} \right] – k_{I,I}(x,y,z,T) I^2(x,y,z,t) \nonumber\\ &\quad – k_{I,V}(x,y,z,T) I(x,y,z,t) V(x,y,z,t) \nonumber\\ &\quad + \Omega \frac{\partial}{\partial x} \left[ \frac{D_{IS}}{kT} \nabla_S \mu(x,y,z,t) \int_0^{L_z} I(x,y,W,t) \, dW \right] \nonumber\\ &\quad + \Omega \frac{\partial}{\partial y} \left[ \frac{D_{IS}}{kT} \nabla_S \mu(x,y,z,t) \int_0^{L_z} I(x,y,W,t) \, dW \right] \nonumber\\ &\quad + \frac{\partial}{\partial x} \left[ \frac{D_{IS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ \frac{D_{IS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial y} \right] + \frac{\partial}{\partial z} \left[ \frac{D_{IS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial z} \right] \end{align}\tag{3}\]

\[\begin{aligned} \frac{\partial V(x,y,z,t)}{\partial t} &= \frac{\partial}{\partial x} \left[ D_V(x,y,z,T) \frac{\partial V(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ D_V(x,y,z,T) \frac{\partial V(x,y,z,t)}{\partial y} \right] \\ &\quad + \frac{\partial}{\partial z} \left[ D_V(x,y,z,T) \frac{\partial V(x,y,z,t)}{\partial z} \right] – k_{V,V}(x,y,z,T) V^2(x,y,z,t) \\ &\quad – k_{I,V}(x,y,z,T) I(x,y,z,t)V(x,y,z,t) \\ &\quad + \Omega \frac{\partial}{\partial x} \left[ \frac{D_{VS}}{kT} \nabla_S \mu(x,y,z,t) \int_0^{L_z} V(x,y,W,t) \, dW \right] \\ &\quad + \Omega \frac{\partial}{\partial y} \left[ \frac{D_{VS}}{kT} \nabla_S \mu(x,y,z,t) \int_0^{L_z} V(x,y,W,t) \, dW \right] \\ &\quad + \frac{\partial}{\partial x} \left[ \frac{D_{VS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ \frac{D_{VS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial y} \right] + \frac{\partial}{\partial z} \left[ \frac{D_{VS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial z} \right] \end{aligned}\] with boundary and initial conditions \[\left. \frac{\partial \, \, I\, \left(x,y,z,t\right)}{\partial \, x} \right|_{x=0} =0,\nonumber\ \] \[\left. \frac{\partial \, \, I\, \left(x,y,z,t\right)}{\partial \, x} \right|_{x=L_{x} } =0,\nonumber\ \] \[ \left. \frac{\partial \, \, I\, \left(x,y,z,t\right)}{\partial \, y} \right|_{y=0} =0,\nonumber\ \] \[\left. \frac{\partial \, \, I\, \left(x,y,z,t\right)}{\partial \, y} \right|_{y=L_{y} } =0,\nonumber\ \] \[\left. \frac{\partial \, \, I\, \left(x,y,z,t\right)}{\partial \, z} \right|_{z=0} =0,\nonumber\ \] \[ \left. \frac{\partial \, \, I\, \left(x,y,z,t\right)}{\partial \, z} \right|_{z=L_{z} } =0,\nonumber\ \] \[\left. \frac{\partial \, \, V\, \left(x,y,z,t\right)}{\partial \, x} \right|_{x=0} =0,\nonumber\ \] \[\left. \frac{\partial \, \, V\, \left(x,y,z,t\right)}{\partial \, x} \right|_{x=L_{x} } =0,\nonumber\ \] \[\left. \frac{\partial \, \, V\, \left(x,y,z,t\right)}{\partial \, y} \right|_{y=0} =0,\nonumber\ \] \[\left. \frac{\partial \, \, V\, \left(x,y,z,t\right)}{\partial \, y} \right|_{y=L_{y} } = 0,\nonumber\ \] \[\left. \frac{\partial \, \, V\, \left(x,y,z,t\right)}{\partial \, z} \right|_{z=0} =0,\nonumber\ \] \[\left. \frac{\partial \, \, V\, \left(x,y,z,t\right)}{\partial \, z} \right|_{z=L_{z} } =0,\nonumber\ \] \[ {I }( {x}, {y}, {z},0) {=} {f{}_{I}}( {x}, {y}, {z}), {V }( {x}, {y}, {z},0)\nonumber\ \] \[ =f_{V}( {x}, {y}, {z}), V\, \left(x_{1} +V_{n} t,y_{1} +V_{n} t,z_{1} +V_{n} t,t\right) =V_{\infty } \left(1+\frac{2\, \ell \, \omega }{k\, T\, \sqrt{x_{1}^{2} +y_{1}^{2} +z_{1}^{2} } } \right). \tag{4}\]

Here I (x,y,z,t) is the spatio-temporal distribution of concentration of radiation interstitials; I\({}^{*}\) is the equilibrium distribution of interstitials; D\({}_{I}\)(x,y,z,T), D\({}_{V}\)(x,y,z,T), D\({}_{IS}\)(x,y, z,T), D\({}_{VS}\)(x,y,z,T) are the coefficients of volumetric and surficial diffusions of interstitials and vacancies, respectively; terms V\({}^{2}\)(x,y,z,t) and I\({}^{2}\)(x,y,z,t) correspond to generation of divacancies and diinterstitials, respectively (see, for example, [24] and appropriate references in this book); k\({}_{I}\)\({}_{,}\)\({}_{V}\)(x,y,z,T), k\({}_{I}\)\({}_{,}\)\({}_{I}\)(x,y,z,T) and k\({}_{V}\)\({}_{,}\)\({}_{V}\)(x,y,z,T) are the parameters of recombination of point radiation defects and generation of their complexes; k is the Boltzmann constant; \(\omega\) = a\({}^{3}\), a is the interatomic distance; \(\ell\) is the specific surface energy. To account porosity of buffer layers we assume, that porous are approximately cylindrical with average values \(r=\sqrt{x_{1}^{2} +y_{1}^{2} }\) and z\({}_{1}\) before annealing [21]. With time small pores decomposing on vacancies. The vacancies absorbing by larger pores [25]. With time large pores became larger due to absorbing the vacancies and became more spherical [25]. Distribution of concentration of vacancies in heterostructure, existing due to porosity, could be determined by summing on all pores, i.e. \[V\left(x,y,z,t\right)=\sum _{i=0}^{l}\sum _{j=0}^{m}\sum _{k=0}^{n}V_{p} \left(x+i\, \alpha ,y+j\, \beta ,\, z+k\, \chi ,t\right) , R=\sqrt{x^{2} +y^{2} +z^{2} } .\] Here \(\alpha\), \(\beta\) and \(\chi\) are the average distances between centers of pores in directions x, y and z; l, m and n are the quantity of pores inappropriate directions.

Spatio-temporal distributions of divacancies \(\Phi\)\({}_{V}\) (x,y,z,t) and diinterstitials \(\Phi\)\({}_{I}\) (x,y,z, t) could be determined by solving the following system of equations [23, 24].

\[\begin{align} \label{GrindEQ__5_} \frac{\partial \Phi_I(x,y,z,t)}{\partial t} &= \frac{\partial}{\partial x} \left[ D_{\Phi_I}(x,y,z,T) \frac{\partial \Phi_I(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ D_{\Phi_I}(x,y,z,T) \frac{\partial \Phi_I(x,y,z,t)}{\partial y} \right] \nonumber\\ &\quad + \frac{\partial}{\partial z} \left[ D_{\Phi_I}(x,y,z,T) \frac{\partial \Phi_I(x,y,z,t)}{\partial z} \right] \nonumber\\ &\quad + \Omega \frac{\partial}{\partial x} \left[ \frac{D_{\Phi_IS}}{kT} \nabla_S \mu_1(x,y,z,t) \int_0^{L_z} \Phi_I(x,y,W,t) \, dW \right] \nonumber\\ &\quad + \Omega \frac{\partial}{\partial y} \left[ \frac{D_{\Phi_IS}}{kT} \nabla_S \mu_1(x,y,z,t) \int_0^{L_z} \Phi_I(x,y,W,t) \, dW \right] \nonumber\\ &\quad + k_{I,I}(x,y,z,T) I^2(x,y,z,t) + k_I(x,y,z,T) I(x,y,z,t) \nonumber\\ &\quad + \frac{\partial}{\partial x} \left[ \frac{D_{\Phi_IS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ \frac{D_{\Phi_IS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial y} \right] + \frac{\partial}{\partial z} \left[ \frac{D_{\Phi_IS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial z} \right] \end{align} \tag{5}\]

\[\begin{aligned} \frac{\partial \Phi_V(x,y,z,t)}{\partial t} &= \frac{\partial}{\partial x} \left[ D_{\Phi_V}(x,y,z,T) \frac{\partial \Phi_V(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ D_{\Phi_V}(x,y,z,T) \frac{\partial \Phi_V(x,y,z,t)}{\partial y} \right] \\ &\quad + \frac{\partial}{\partial z} \left[ D_{\Phi_V}(x,y,z,T) \frac{\partial \Phi_V(x,y,z,t)}{\partial z} \right] \\ &\quad + \Omega \frac{\partial}{\partial x} \left[ \frac{D_{\Phi_V S}}{kT} \nabla_S \mu_1(x,y,z,t) \int_0^{L_z} \Phi_V(x,y,W,t) \, dW \right] \\ &\quad + \Omega \frac{\partial}{\partial y} \left[ \frac{D_{\Phi_V S}}{kT} \nabla_S \mu_1(x,y,z,t) \int_0^{L_z} \Phi_V(x,y,W,t) \, dW \right] \\ &\quad + k_{V,V}(x,y,z,T) V^2(x,y,z,t) + k_V(x,y,z,T) V(x,y,z,t) \\ &\quad + \frac{\partial}{\partial x} \left[ \frac{D_{\Phi_V S}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ \frac{D_{\Phi_V S}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial y} \right] + \frac{\partial}{\partial z} \left[ \frac{D_{\Phi_V S}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial z} \right] \end{aligned}\] with boundary and initial conditions \[\begin{aligned} \label{GrindEQ__6_} \left. \frac{\partial \Phi _{I} \left(x,y,z,t\right)}{\partial x} \right|_{x=0} =&0,\nonumber\\ \left. \frac{\partial I \left(x,y,z,t\right)}{\partial x} \right|_{x=L_{x} } =&0,\nonumber\\ \left. \frac{\partial I \left(x,y,z,t\right)}{\partial y} \right|_{y=0} =&0,\nonumber\\ \left. \frac{\partial I \left(x,y,z,t\right)}{\partial y} \right|_{y=L_{y} } =&0,\nonumber\\ \left. \frac{\partial \Phi _{I} \left(x,y,z,t\right)}{\partial z} \right|_{z=0} =&0,\nonumber\\ \left. \frac{\partial I \left(x,y,z,t\right)}{\partial z} \right|_{z=L_{z} } =&0,\nonumber\\ \left. \frac{\partial \Phi _{V} \left(x,y,z,t\right)}{\partial x} \right|_{x=0} =&0,\nonumber\\ \left. \frac{\partial V \left(x,y,z,t\right)}{\partial x} \right|_{x=L_{x} } =&0,\nonumber\\ \left. \frac{\partial V \left(x,y,z,t\right)}{\partial y} \right|_{y=0} =&0, \nonumber\\ \left. \frac{\partial V \left(x,y,z,t\right)}{\partial y} \right|_{y=L_{y} } =&0,\nonumber\\ \left. \frac{\partial V \left(x,y,z,t\right)}{\partial z} \right|_{z=0} =&0,\nonumber\\ \left. \frac{\partial \Phi _{V} \left(x,y,z,t\right)}{\partial z} \right|_{z=L_{z} } =&0,\nonumber\\ \Phi_{I}({x},{y},{z},0)=&f_{I}( {x}, {y}, {z}),\nonumber\\ {\Phi_{V} }( {x}, {y}, {z},0)=&f_{V}( {x}, {y}, {z}). \end{aligned} \tag{6}\]

Here D\({}_{I}\)(x,y,z,T), D\({}_{V}\)(x,y,z,T), D\({}_{IS\ }\)(x,y,z,T) and D\({}_{VS}\)(x,y,z,T) are the coefficients of volumetric and surficial diffusions of complexes of radiation defects; k\({}_{I}\)(x,y,z,T) and k\({}_{V}\)(x,y,z,T) are the parameters of decay of complexes of radiation defects.

Chemical potential \(\mu\)\({}_{1}\) in Eq.(1) could be determine by the following relation [18];

\[\label{GrindEQ__7_} \mu_1=E(z)\sigma_{ij} [uij(x,y,z,t)+uji(x,y,z,t)]/2, \tag{7}\] where E(z) is the Young modulus, \(\sigma\)\({}_{ij}\) is the stress tensor; \(u_{ij} =\frac{1}{2} \left(\frac{\partial \, u_{i} }{\partial \, x_{j} } +\frac{\partial \, u_{j} }{\partial \, x_{i} } \right)\) is the deformation tensor; u\({}_{i}\), u\({}_{j}\) are the components u\({}_{x}\)(x,y,z,t), u\({}_{y}\)(x,y,z,t) and u\({}_{z}\)(x,y,z,t) of the displacement vector \(\vec{u}\, \left(x,y,z,t\right)\); x\({}_{i}\), x\({}_{j}\) are the coordinate x, y, z. The Eq. (3) could be transform to the following form \[\begin{aligned} \mu \, \, \left(x,y,z,t\right)=&\left[\frac{\partial \, u_{i} \left(x,y,z,t\right)}{\partial \, x_{j} } +\frac{\partial \, u_{j} \left(x,y,z,t\right)}{\partial \, x_{i} } \right]\, \left\{\frac{1}{2} \left[\frac{\partial \, u_{i} \left(x,y,z,t\right)}{\partial \, x_{j} } +\frac{\partial \, u_{j} \left(x,y,z,t\right)}{\partial \, x_{i} } \right]\right. \\ &\left. -\varepsilon _{0} \delta _{ij} +\frac{\sigma \, \left(z\right)\, \delta _{ij} }{1-2\, \sigma \, \left(z\right)} \left[\frac{\partial \, u_{k} \left(x,y,z,t\right)}{\partial \, x_{k} } -3\, \varepsilon _{0} \right]-K\, \left(z\right)\, \beta \, \left(z\right)\, \left[T\, \left(x,y,z,t\right)-T_{0} \right]\, \, \delta _{ij} \right\}\frac{\Omega }{2} E\, \left(z\right), \end{aligned}\] where \(\sigma\) is Poisson coefficient; \(\epsilon\)\({}_{0}\) = (a\({}_{s}\)–a\({}_{EL}\))/a\({}_{EL}\) is the mismatch parameter; a\({}_{s}\), a\({}_{EL}\) are lattice distances of the substrate and the epitaxial layer; K is the modulus of uniform compression; \(\beta\) is the coefficient of thermal expansion; T\({}_{r}\) is the equilibrium temperature, which coincide (for our case) with room temperature. Components of displacement vector could be obtained by solution of the following equations [19]; \[\rho \, \, \left(z\right)\frac{\partial ^{2} \, u_{x} \left(x,y,z,t\right)}{\partial \, t^{2} } =\frac{\partial \, \sigma _{xx} \left(x,y,z,t\right)}{\partial \, x} +\frac{\partial \, \sigma _{xy} \left(x,y,z,t\right)}{\partial \, y} +\frac{\partial \, \sigma _{xz} \left(x,y,z,t\right)}{\partial \, z}\] \[\rho \, \, \left(z\right)\frac{\partial ^{2} \, u_{y} \left(x,y,z,t\right)}{\partial \, t^{2} } =\frac{\partial \, \sigma _{yx} \left(x,y,z,t\right)}{\partial \, x} +\frac{\partial \, \sigma _{yy} \left(x,y,z,t\right)}{\partial \, y} +\frac{\partial \, \sigma _{yz} \left(x,y,z,t\right)}{\partial \, z}\] \[\rho \, \, \left(z\right)\frac{\partial ^{2} \, u_{z} \left(x,y,z,t\right)}{\partial \, t^{2} } =\frac{\partial \, \sigma _{zx} \left(x,y,z,t\right)}{\partial \, x} +\frac{\partial \, \sigma _{zy} \left(x,y,z,t\right)}{\partial \, y} +\frac{\partial \, \sigma _{zz} \left(x,y,z,t\right)}{\partial \, z} ,\] where \(\sigma _{ij} =\dfrac{E\, \left(z\right)}{2\, \left[1+\sigma \, \left(z\right)\right]} \left[\dfrac{\partial \, u_{i} \left(x,y,z,t\right)}{\partial \, x_{j} } +\dfrac{\partial \, u_{j} \left(x,y,z,t\right)}{\partial \, x_{i} } -\dfrac{\delta _{ij} }{3} \dfrac{\partial \, u_{k} \left(x,y,z,t\right)}{\partial \, x_{k} } \right]+K\, \left(z\right)\, \delta _{ij} \times\) \(\times \dfrac{\partial \, u_{k} \left(x,y,z,t\right)}{\partial \, x_{k} } -\beta \, \left(z\right)\, K\, \left(z\right)\, \left[T\left(x,y,z,t\right)-T_{r} \right]\), \(\rho\) (z) is the density of materials of heterostructure, \(\delta\)\({}_{ij}\) Is the Kronecker symbol. With account the relation for \(\sigma\)\({}_{ij}\) last system of equation could be written as \[\begin{aligned} \rho(z) \frac{\partial^2 u_x(x,y,z,t)}{\partial t^2} &= \left\{ K(z) + \frac{5E(z)}{6[1+\sigma(z)]} \right\} \frac{\partial^2 u_x(x,y,z,t)}{\partial x^2} + \left\{ K(z) – \frac{E(z)}{3[1+\sigma(z)]} \right\} \frac{\partial^2 u_y(x,y,z,t)}{\partial x \partial y} \\ &\quad + \frac{E(z)}{2[1+\sigma(z)]} \left( \frac{\partial^2 u_y(x,y,z,t)}{\partial y^2} + \frac{\partial^2 u_z(x,y,z,t)}{\partial z^2} \right) \\ &\quad + \left\{ K(z) + \frac{E(z)}{3[1+\sigma(z)]} \right\} \frac{\partial^2 u_x(x,y,z,t)}{\partial x \partial z} – K(z) \beta(z) \frac{\partial T(x,y,z,t)}{\partial x},\end{aligned}\] \[\begin{align} \label{GrindEQ__8_} \rho(z) \frac{\partial^2 u_y(x,y,z,t)}{\partial t^2} &= \frac{E(z)}{2[1+\sigma(z)]} \left( \frac{\partial^2 u_y(x,y,z,t)}{\partial x^2} + \frac{\partial^2 u_x(x,y,z,t)}{\partial x \partial y} \right) – \frac{\partial T(x,y,z,t)}{\partial y} \nonumber\\ &\quad + \left\{ K(z) \beta(z) + \frac{\partial}{\partial z} \left( \frac{E(z)}{1+\sigma(z)} \left[ \frac{\partial u_y(x,y,z,t)}{\partial z} + \frac{\partial u_z(x,y,z,t)}{\partial y} \right] \right) \right\} \nonumber\\ &\quad + \left\{ \frac{5E(z)}{12[1+\sigma(z)]} + K(z) \right\} \frac{\partial^2 u_y(x,y,z,t)}{\partial y^2} + \left\{ K(z) – \frac{E(z)}{6[1+\sigma(z)]} \right\} \frac{\partial^2 u_y(x,y,z,t)}{\partial y \partial z}.\end{align} \tag{8}\]

\[\begin{aligned} \rho(z) \frac{\partial^2 u_z(x,y,z,t)}{\partial t^2} &= \frac{E(z)}{2[1+\sigma(z)]} \left[ \frac{\partial^2 u_z(x,y,z,t)}{\partial x^2} + \frac{\partial^2 u_z(x,y,z,t)}{\partial y^2} + \frac{\partial^2 u_x(x,y,z,t)}{\partial x \partial z} + \frac{\partial^2 u_y(x,y,z,t)}{\partial y \partial z} \right] \\ &\quad + \frac{\partial}{\partial z} \left\{ K(z) \left[ \frac{\partial u_x(x,y,z,t)}{\partial x} + \frac{\partial u_y(x,y,z,t)}{\partial y} + \frac{\partial u_z(x,y,z,t)}{\partial z} \right] \right\} \\ &\quad + \frac{1}{6} \frac{\partial}{\partial z} \left\{ \frac{E(z)}{1+\sigma(z)} \left[ 6\frac{\partial u_z(x,y,z,t)}{\partial z} – \frac{\partial u_x(x,y,z,t)}{\partial x} – \frac{\partial u_y(x,y,z,t)}{\partial y} – \frac{\partial u_z(x,y,z,t)}{\partial z} \right] \right\} \\ &\quad – K(z) \beta(z) \frac{\partial T(x,y,z,t)}{\partial z}. \end{aligned}\]

Conditions for the system of Eq. (8) could be written in the form \[\frac{\partial \, \, \vec{u}\left(0,y,z,t\right)}{\partial \, x} =0; \frac{\partial \, \vec{u}\, \left(L_{x} ,y,z,t\right)}{\partial \, x} =0; \frac{\partial \, \, \vec{u}\, \left(x,0,z,t\right)}{\partial \, y} =0; \frac{\partial \, \, \vec{u}\, \left(x,L_{y} ,z,t\right)}{\partial \, y} =0;\] \[\frac{\partial \, \, \vec{u}\, \left(x,y,0,t\right)}{\partial \, z} =0; \frac{\partial \, \, \vec{u}\, \left(x,y,L_{z} ,t\right)}{\partial \, z} =0; \vec{u}\, \left(x,y,z,0\right)=\vec{u}_{0} ; \vec{u}\, \left(x,y,z,\infty \right)=\vec{u}_{0} .\]

We determine spatio-temporal distributions of concentrations of dopant and radiati-on defects by solving the Eqs.(1), (3) and (5) framework standard method of averaging of function corrections [26]. Previously we transform the Eqs.(1), (3) and (5) to the following form with account initial distributions of the considered concentrations

\[\begin{align} \frac{\partial C(x,y,z,t)}{\partial t} &= \frac{\partial}{\partial x} \left[ D \frac{\partial C(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ D \frac{\partial C(x,y,z,t)}{\partial y} \right] + \frac{\partial}{\partial z} \left[ D \frac{\partial C(x,y,z,t)}{\partial z} \right] \nonumber\\ &\quad + f_C(x,y,z)\delta(t) + \frac{\partial}{\partial x} \left[ \frac{D_{CS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ \frac{D_{CS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial y} \right] \nonumber\\ &\quad + \frac{\partial}{\partial z} \left[ \frac{D_{CS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial z} \right] + \Omega \frac{\partial}{\partial x} \left[ \frac{D_S}{k T} \nabla_S \mu(x,y,z,t) \int_{0}^{L_{\xi}} C(x,y,W,t) \, dW \right] \nonumber\\ &\quad + \Omega \frac{\partial}{\partial y} \left[ \frac{D_S}{k T} \nabla_S \mu(x,y,z,t) \int_{0}^{L_{\xi}} C(x,y,W,t) \, dW \right] \end{align} \tag{1a}\]

\[\begin{align} \frac{\partial I(x,y,z,t)}{\partial t} &= \frac{\partial}{\partial x} \left[ D_I(x,y,z,T) \frac{\partial I(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ D_I(x,y,z,T) \frac{\partial I(x,y,z,t)}{\partial y} \right] \\ &\quad + \frac{\partial}{\partial z} \left[ D_I(x,y,z,T) \frac{\partial I(x,y,z,t)}{\partial z} \right] + \Omega \frac{\partial}{\partial x} \left[ \frac{D_{IS}}{kT} \nabla_S \mu_1(x,y,z,t) \int_{0}^{L_z} I(x,y,W,t) \, dW \right] \\ &\quad + \Omega \frac{\partial}{\partial y} \left[ \frac{D_{IS}}{kT} \nabla_S \mu_1(x,y,z,t) \int_{0}^{L_z} I(x,y,W,t) \, dW \right] \\ &\quad – k_{I,I}(x,y,z,T) I^2(x,y,z,t) – k_{I,V}(x,y,z,T) I(x,y,z,t) V(x,y,z,t) + f_I(x,y,z)\delta(t) \end{align}\tag{3a}\] \[\begin{align} \frac{\partial V(x,y,z,t)}{\partial t} &= \frac{\partial}{\partial x} \left[ D_V(x,y,z,T) \frac{\partial V(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ D_V(x,y,z,T) \frac{\partial V(x,y,z,t)}{\partial y} \right]\\&\quad + \frac{\partial}{\partial z} \left[ D_V(x,y,z,T) \frac{\partial V(x,y,z,t)}{\partial z} \right] + \Omega \frac{\partial}{\partial x} \left[ \frac{D_{VS}}{kT} \nabla_S \mu_1(x,y,z,t) \int_0^{L_z} V(x,y,W,t) \, dW \right]\\&\quad + \Omega \frac{\partial}{\partial y} \left[ \frac{D_{VS}}{kT} \nabla_S \mu_1(x,y,z,t) \int_0^{L_z} V(x,y,W,t) \, dW \right] – k_{V,V}(x,y,z,T) V^2(x,y,z,t)\\ &\quad – k_{I,V}(x,y,z,T) I(x,y,z,t)V(x,y,z,t) + f_V(x,y,z)\delta(t) + \frac{\partial \Phi_I(x,y,z,t)}{\partial t} \\&= \frac{\partial}{\partial x} \left[ D_{\Phi_I}(x,y,z,T) \frac{\partial \Phi_I(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ D_{\Phi_I}(x,y,z,T) \frac{\partial \Phi_I(x,y,z,t)}{\partial y} \right] \\ &\quad + \frac{\partial}{\partial z} \left[ D_{\Phi_I}(x,y,z,T) \frac{\partial \Phi_I(x,y,z,t)}{\partial z} \right] + \Omega \frac{\partial}{\partial x} \left[ \frac{D_{\Phi_I,S}}{kT} \nabla_S \mu_1(x,y,z,t) \int_0^{L_z} \Phi_I(x,y,W,t) \, dW \right] \\ &\quad + \Omega \frac{\partial}{\partial y} \left[ \frac{D_{\Phi_I,S}}{kT} \nabla_S \mu_1(x,y,z,t) \int_0^{L_z} \Phi_I(x,y,W,t) \, dW \right] \\ &\quad + \frac{\partial}{\partial x} \left[ \frac{D_{\Phi,S}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ \frac{D_{\Phi,S}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial y} \right] + \frac{\partial}{\partial z} \left[ \frac{D_{\Phi,S}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial z} \right] \\ &\quad + k_{I,I}(x,y,z,T) I^2(x,y,z,t) + f_{\Phi}(x,y,z)\delta(t) \end{align} \tag{5a}\]

\[\begin{aligned} \frac{\partial \Phi_V(x,y,z,t)}{\partial t} &= \frac{\partial}{\partial x} \left[ D_{\Phi_V}(x,y,z,T) \frac{\partial \Phi_V(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ D_{\Phi_V}(x,y,z,T) \frac{\partial \Phi_V(x,y,z,t)}{\partial y} \right] \\ &\quad + \frac{\partial}{\partial z} \left[ D_{\Phi_V}(x,y,z,T) \frac{\partial \Phi_V(x,y,z,t)}{\partial z} \right] + \Omega \frac{\partial}{\partial x} \left[ \frac{D_{\Phi_V S}}{kT} \nabla_S \mu_1(x,y,z,t) \int_0^{L_z} \Phi_V(x,y,W,t) \, dW \right] \\ &\quad + \Omega \frac{\partial}{\partial y} \left[ \frac{D_{\Phi_V S}}{kT} \nabla_S \mu_1(x,y,z,t) \int_0^{L_z} \Phi_I(x,y,W,t) \, dW \right] + k_I(x,y,z,T) I(x,y,z,t) \\ &\quad + \frac{\partial}{\partial x} \left[ \frac{D_{\Phi_V S}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ \frac{D_{\Phi_V S}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial y} \right] + \frac{\partial}{\partial z} \left[ \frac{D_{\Phi_V S}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial z} \right] \\ &\quad + k_{V,V}(x,y,z,T) V^2(x,y,z,t) + f_{\Phi_V}(x,y,z) \delta(t). \end{aligned}\]

Farther we replace concentrations of dopant and radiation defects in right sides of Eqs. (1a), (3a) and (5a) on their not yet known average values \(\alpha\)\({}_{1}\). In this situation we obtain equations for the first-order approximations of the required concentrations in the following form

\[\begin{align} \frac{\partial C_1(x,y,z,t)}{\partial t} &= \alpha_{1C} \Omega \frac{\partial}{\partial x} \left[ z \frac{D_S}{kT} \nabla_S \mu_1(x,y,z,t) \right] + \alpha_{1C} \Omega \frac{\partial}{\partial y} \left[ z \frac{D_S}{kT} \nabla_S \mu_1(x,y,z,t) \right] \\ &\quad + f_C(x,y,z)\delta(t) + \frac{\partial}{\partial x} \left[ \frac{D_{CS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ \frac{D_{CS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial y} \right] \\ &\quad + \frac{\partial}{\partial z} \left[ \frac{D_{CS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial z} \right] \end{align}\tag{1b}\] \[\begin{align} \frac{\partial I_1(x,y,z,t)}{\partial t} &= \alpha_{1I} z \Omega \frac{\partial}{\partial x} \left[ \frac{D_{IS}}{kT} \nabla_S \mu(x,y,z,t) \right] + \alpha_{1I} \Omega \frac{\partial}{\partial y} \left[ z \frac{D_{IS}}{kT} \nabla_S \mu(x,y,z,t) \right] \\ &\quad + \frac{\partial}{\partial x} \left[ \frac{D_{IS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ \frac{D_{IS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial y} \right] + \frac{\partial}{\partial z} \left[ \frac{D_{IS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial z} \right] \\ &\quad + f_I(x,y,z)\delta(t) – \alpha_{1I}^2 k_{I,I}(x,y,z,T) – \alpha_{1I} \alpha_{1V} k_{I,V}(x,y,z,T) \end{align}\tag{3b}\] \[\begin{align} \frac{\partial V_1(x,y,z,t)}{\partial t} &= \alpha_{1V} z \Omega \frac{\partial}{\partial x} \left[ \frac{D_{VS}}{kT} \nabla_S \mu_1(x,y,z,t) \right] + \alpha_{1V} \Omega \frac{\partial}{\partial y} \left[ z \frac{D_{VS}}{kT} \nabla_S \mu_1(x,y,z,t) \right] \\ &\quad + \frac{\partial}{\partial x} \left[ \frac{D_{VS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ \frac{D_{VS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial y} \right] + \frac{\partial}{\partial z} \left[ \frac{D_{VS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial z} \right] \\ &\quad + f_V(x,y,z)\delta(t) – \alpha_{1V}^2 k_{V,V}(x,y,z,T) – \alpha_{1V} \alpha_{1I} k_{I,V}(x,y,z,T) \end{align}\] \[\begin{align} \frac{\partial \Phi_1(x,y,z,t)}{\partial t} &= \alpha_{\Phi_1} z \Omega \frac{\partial}{\partial x} \left[ \frac{D_{\Phi_I S}}{kT} \nabla_S \mu_1(x,y,z,t) \right] + \alpha_{\Phi_1} \Omega \frac{\partial}{\partial y} \left[ z \frac{D_{\Phi_I S}}{kT} \nabla_S \mu_1(x,y,z,t) \right] \\ &\quad + \frac{\partial}{\partial x} \left[ \frac{D_{\Phi_I S}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ \frac{D_{\Phi_I S}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial y} \right] + \frac{\partial}{\partial z} \left[ \frac{D_{\Phi_I S}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial z} \right] \\ &\quad + f_{\Phi_I}(x,y,z)\delta(t) + k_I(x,y,z,T)I(x,y,z,t) + k_{I,I}(x,y,z,T)I^2(x,y,z,t) \end{align} \tag{5b}\]

\[\begin{aligned} \frac{\partial \, \, \Phi _{1V} \left(x,y,z,t\right)}{\partial \, \, t} =&\alpha _{1\Phi _{V} } z\, \Omega \frac{\partial }{\partial \, x} \left[\frac{D_{\Phi _{V} S} }{k\, T} \nabla _{S} \mu _{1} \left(x,y,z,t\right)\right]+\alpha _{1\Phi _{V} } z\, \Omega \frac{\partial }{\partial \, y} \left[\frac{D_{\Phi _{V} S} }{k\, T} \nabla _{S} \mu _{1} \left(x,y,z,t\right)\right]+\\ &+\frac{\partial }{\partial \, x} \left[\frac{D_{\Phi _{V} S} }{\bar{V}\, k\, T} \frac{\partial \, \mu _{2} \left(x,y,z,t\right)}{\partial \, x} \right]+\frac{\partial }{\partial \, y} \left[\frac{D_{\Phi _{V} S} }{\bar{V}\, k\, T} \frac{\partial \, \mu _{2} \left(x,y,z,t\right)}{\partial \, y} \right]+\frac{\partial }{\partial \, z} \left[\frac{D_{\Phi _{V} S} }{\bar{V}\, k\, T} \frac{\partial \, \mu _{2} \left(x,y,z,t\right)}{\partial \, z} \right]+\\ &+f_{\Phi _{V} } \, \left(x,y,z\right)\, \delta \left(t\right)+k_{V} \left(x,y,z,T\right)\, V\, \left(x,y,z,t\right)+k_{V,V} \left(x,y,z,T\right)\, V^{2} \left(x,y,z,t\right). \end{aligned}\]

Integration of the left and right sides of the Eqs. (1b), (3b) and (5b) on time gives us possibility to obtain relations for above approximation in the final form

\[\begin{align} C_{1} \left(x,y,z,t\right)=&\alpha _{1C} \Omega \frac{\partial }{\partial \, x} \int _{0}^{t}D_{S\, L} \left(x,y,z,T\right)\, \frac{z}{k\, T} \left[1+\varsigma _{1} \frac{V\, \left(x,y,z,\tau \right)}{V^{*} } +\varsigma _{2} \frac{V^{2} \left(x,y,z,\tau \right)}{\left(V^{*} \right)^{2} } \right]\, \, \\ &\left. \times \, \nabla _{S} \mu _{1} \left(x,y,z,\tau \right)\, \left[1+\frac{\xi _{S} \alpha _{1C}^{\gamma } }{P^{\gamma } \left(x,y,z,T\right)} \right]\, d\, \tau \right\}+\alpha _{1C} \frac{\partial }{\partial \, y} \int _{0}^{t}D_{S\, L} \left(x,y,z,T\right)\, \left[1+\frac{\xi _{S} \alpha _{1C}^{\gamma } }{P^{\gamma } \left(x,y,z,T\right)} \right]\, \, \\& +\times \, \, \Omega \, \nabla _{S} \mu _{1} \left(x,y,z,\tau \right)\, \frac{z}{k\, T} \left[1+\varsigma _{1} \frac{V\, \left(x,y,z,\tau \right)}{V^{*} } +\varsigma _{2} \frac{V^{2} \left(x,y,z,\tau \right)}{\left(V^{*} \right)^{2} } \right]\, \, d\, \tau +\frac{\partial }{\partial \, \, x} \int _{0}^{t}\frac{D_{C\, S} }{\bar{V}\, k\, T} \times \\& \times \frac{\partial \, \, \mu _{2} \left(x,y,z,\tau \right)}{\partial \, \, x} \, d\, \tau +\frac{\partial }{\partial \, \, y} \int _{0}^{t}\frac{D_{C\, S} }{\bar{V}\, k\, T} \frac{\partial \, \, \mu _{2} \left(x,y,z,\tau \right)}{\partial \, \, y} \, d\, \tau +\frac{\partial }{\partial \, \, z} \int _{0}^{t}\frac{D_{C\, S} }{\bar{V}\, k\, T} \frac{\partial \, \, \mu _{2} \left(x,y,z,\tau \right)}{\partial \, \, z} \, d\, \tau +\\ &+f_{C} \, \left(x,y,z\right) \end{align}\tag{1c}\] \[\begin{align}I_{1} \left(x,y,z,t\right)=&\alpha _{1I} z\, \Omega \frac{\partial }{\partial \, x} \int _{0}^{t}\frac{D_{IS} }{k\, T} \nabla _{S} \mu _{1} \left(x,y,z,\tau \right)\, d\, \tau +\alpha _{1I} z\, \Omega \frac{\partial }{\partial \, y} \int _{0}^{t}\frac{D_{IS} }{k\, T} \nabla _{S} \mu _{1} \left(x,y,z,\tau \right)\, \, d\, \tau \\ &+\frac{\partial }{\partial \, x} \int _{0}^{t}\frac{D_{I\, S} }{\bar{V}\, k\, T} \frac{\partial \, \mu _{2} \left(x,y,z,t\right)}{\partial \, x} d\, \tau +\frac{\partial }{\partial \, y} \int _{0}^{t}\frac{D_{I\, S} }{\bar{V}\, k\, T} \frac{\partial \, \mu _{2} \left(x,y,z,t\right)}{\partial \, y} d\, \tau +\frac{\partial }{\partial \, z} \int _{0}^{t}\frac{D_{I\, S} }{\bar{V}\, k\, T} \frac{\partial \, \mu _{2} \left(x,y,z,t\right)}{\partial \, z} d\, \tau \\ &+f_{I} \, \left(x,y,z\right)-\alpha _{1I}^{2} \int _{0}^{t}k_{I,I} \left(x,y,z,T\right)\, d\, \tau -\alpha _{1I} \alpha _{1V} \int _{0}^{t}k_{I,V} \\ &\left(x,y,z,T\right)\, \, d\, \tau \end{align}\tag{3c}\] \[\begin{aligned} V_{1} \left(x,y,z,t\right)=&\alpha _{1V} z\, \Omega \frac{\partial }{\partial \, x} \int _{0}^{t}\frac{D_{IS} }{k\, T} \nabla _{S} \mu _{1} \left(x,y,z,\tau \right)\, d\, \tau +\alpha _{1V} z\, \Omega \frac{\partial }{\partial \, y} \int _{0}^{t}\frac{D_{IS} }{k\, T} \nabla _{S} \mu _{1} \left(x,y,z,\tau \right)\, \, d\, \tau \\ &+\frac{\partial }{\partial \, x} \int _{0}^{t}\frac{D_{V\, S} }{\bar{V}\, k\, T} \frac{\partial \, \mu _{2} \left(x,y,z,t\right)}{\partial \, x} d\, \tau +\frac{\partial }{\partial \, y} \int _{0}^{t}\frac{D_{V\, S} }{\bar{V}\, k\, T} \frac{\partial \, \mu _{2} \left(x,y,z,t\right)}{\partial \, y} d\, \tau +\frac{\partial }{\partial \, z} \int _{0}^{t}\frac{D_{V\, S} }{\bar{V}\, k\, T} \frac{\partial \, \mu _{2} \left(x,y,z,t\right)}{\partial \, z} d\, \tau \\ &+f_{V} \, \left(x,y,z\right)-\alpha _{1V}^{2} \int _{0}^{t}k_{V,V} \left(x,y,z,T\right)\, d\, \tau -\alpha _{1I} \alpha _{1V} \int _{0}^{t}k_{I,V} \left(x,y,z,T\right)\, \, d\, \tau \\ &\Phi _{1I} \left(x,y,z,t\right)=\alpha _{1\Phi _{I} } z\, \Omega \frac{\partial }{\partial \, x} \int _{0}^{t}\frac{D_{\Phi _{I} S} }{k\, T} \nabla _{S} \mu _{1} \left(x,y,z,\tau \right)\, \, d\, \tau +\, \Omega \frac{\partial }{\partial \, x} \int _{0}^{t}\frac{D_{\Phi _{I} S} }{k\, T} \nabla _{S} \mu _{1} \left(x,y,z,\tau \right)\, \, d\, \tau \times \\ &\times \, \alpha _{1\Phi _{I} } z+f_{\Phi _{I} } \, \left(x,y,z\right)+\frac{\partial }{\partial \, x} \int _{0}^{t}\frac{D_{\Phi _{I} S} }{\bar{V}\, k\, T} \frac{\partial \, \mu _{2} \left(x,y,z,\tau \right)}{\partial \, x} d\, \tau +\frac{\partial }{\partial \, y} \int _{0}^{t}\frac{D_{\Phi _{I} S} }{\bar{V}\, k\, T} \frac{\partial \, \mu _{2} \left(x,y,z,\tau \right)}{\partial \, y} d\, \tau +\\ &+\frac{\partial }{\partial \, z} \int _{0}^{t}\frac{D_{\Phi _{I} S} }{\bar{V}\, k\, T} \frac{\partial \, \mu _{2} \left(x,y,z,\tau \right)}{\partial \, z} d\, \tau +\int _{0}^{t}k_{I} \left(x,y,z,T\right)\, I\, \left(x,y,z,\tau \right)\, d\, \tau \\ &+\int _{0}^{t}k_{I,I} \left(x,y,z,T\right)\, I^{2} \left(x,y,z,\tau \right)\, d\, \tau \end{aligned} \tag{5c}\]

\[\begin{aligned} \Phi _{1V} \left(x,y,z,t\right)=&\alpha _{1\Phi _{V} } z\, \Omega \frac{\partial }{\partial \, x} \int _{0}^{t}\frac{D_{\Phi _{V} S} }{k\, T} \nabla _{S} \mu _{1} \left(x,y,z,\tau \right)\, \, d\, \tau +\, \Omega \frac{\partial }{\partial \, x} \int _{0}^{t}\frac{D_{\Phi _{V} S} }{k\, T} \nabla _{S} \mu _{1} \left(x,y,z,\tau \right)\, \, d\, \tau \times \\ &\times \, \alpha _{1\Phi _{V} } z+f_{\Phi _{V} } \, \left(x,y,z\right)+\frac{\partial }{\partial \, x} \int _{0}^{t}\frac{D_{\Phi _{V} S} }{\bar{V}\, k\, T} \frac{\partial \, \mu _{2} \left(x,y,z,\tau \right)}{\partial \, x} d\, \tau +\frac{\partial }{\partial \, y} \int _{0}^{t}\frac{D_{\Phi _{V} S} }{\bar{V}\, k\, T} \frac{\partial \, \mu _{2} \left(x,y,z,\tau \right)}{\partial \, y} d\, \tau \\ &+\frac{\partial }{\partial \, z} \int _{0}^{t}\frac{D_{\Phi _{V} S} }{\bar{V}\, k\, T} \frac{\partial \, \mu _{2} \left(x,y,z,\tau \right)}{\partial \, z} d\, \tau +\int _{0}^{t}k_{V} \left(x,y,z,T\right)\, V\, \left(x,y,z,\tau \right)\, d\, \tau \\ &+\int _{0}^{t}k_{V,V} \left(x,y,z,T\right)\, V^{2} \left(x,y,z,\tau \right)\, d\, \tau . \end{aligned}\]

We determine average values of the first-order approximations of concentrations of dopant and radiation defects by the following standard relation [26]; \[\label{GrindEQ__9_} \alpha _{1\rho } =\frac{1}{\Theta L_{x} L_{y} L_{z} } \int _{0}^{\Theta }\int _{0}^{L_{x} }\int _{0}^{L_{y} }\int _{0}^{L_{z} }\rho _{1} \left(x,y,z,t\right)\, d\, z\, d\, y\, d\, x\, d\, t . \tag{9}\]

Substitution of the relations (1c), (3c) and (5c) into relation (9) gives us possibility to obtain required average values in the following form \[\begin{aligned} \alpha _{1C} =&\frac{1}{L_{x} L_{y} L_{z} } \int _{0}^{L_{x} }\int _{0}^{L_{y} }\int _{0}^{L_{z} }f_{C} \left(x,y,z\right)\, \, d\, z\, d\, y\, d\, x ,\\ \alpha _{1I} =&\sqrt{\frac{\left(a_{3} +A\right)^{2} }{4\, a_{4}^{2} } -4\, \, \left(B+\frac{\Theta \, a_{3} B+\Theta ^{2} L_{x} L_{y} L_{z} a_{1} }{a_{4} } \right)} \frac{a_{3} +A}{4\, a_{4} } ,\\ \alpha _{1V} =&\frac{1}{S_{IV\, 00} } \left[\frac{\Theta }{\alpha _{1I} } \int _{0}^{L_{x} }\int _{0}^{L_{y} }\int _{0}^{L_{z} }f_{I} \, \left(x,y,z\right)\, \, d\, z\, d\, y\, d\, x -\alpha _{1I} S_{II\, 00} -\Theta L_{x} L_{y} L_{z} \right], \end{aligned}\] where \[\begin{aligned} S_{\rho \rho \, \, ij} =&\int _{0}^{\Theta }\left(\Theta -t\right)\int _{0}^{L_{x} }\int _{0}^{L_{y} }\int _{0}^{L_{z} }k_{\rho ,\rho } \left(x,y,z,T\right)\, I_{1}^{i} \left(x,y,z,t\right)\, V_{1}^{j} \left(x,y,z,t\right)\, d\, z\, d\, y\, d\, x\, d\, t ,\\ a_{4} =&S_{II\, 00} \times \, \, \left(S_{IV\, 00}^{2} -S_{II\, 00} S_{VV\, 00} \right),\\ a_{3} =&S_{IV00} S_{II00} +S_{IV00}^{2} -S_{II00} S_{VV00} ,\\ a_{2} =&\int _{0}^{L_{x} }\int _{0}^{L_{y} }\int _{0}^{L_{z} }f_{V} \, \left(x,y,z\right)\, d\, z\, d\, y\, d\, x \times S_{IV00} S_{IV00}^{2} +S_{IV00} \Theta L_{x}^{2} L_{y}^{2} L_{z}^{2}\\ &+2\, S_{VV00} S_{II00} \int _{0}^{L_{x} }\int _{0}^{L_{y} }\int _{0}^{L_{z} }f_{I} \, \left(x,y,z\right)\, d\, z\, d\, y\, d\, x -\Theta \, L_{x}^{2} L_{y}^{2} L_{z}^{2} S_{VV00}\\ &-S_{IV00}^{2} \int _{0}^{L_{x} }\int _{0}^{L_{y} }\int _{0}^{L_{z} }f_{I} \, \left(x,y,z\right)\, d\, z\, d\, y\, d\, x , \\ a_{1} =&S_{IV00} \int _{0}^{L_{x} }\int _{0}^{L_{y} }\int _{0}^{L_{z} }f_{I} \, \left(x,y,z\right)\, d\, z\, d\, y\, d\, x ,\\ a_{0} =&S_{VV00} \times \left[\int _{0}^{L_{x} }\int _{0}^{L_{y} }\int _{0}^{L_{z} }f_{I} \, \left(x,y,z\right)\, d\, z\, d\, y\, d\, x \right]^{2} , \\ A=&\sqrt{8\, y+\Theta ^{2} \frac{a_{3}^{2} }{a_{4}^{2} } -4\, \Theta \, \frac{a_{2} }{a_{4} } } ,\\ B=&\frac{\Theta \, a_{2} }{6\, a_{4} } +\sqrt[{3}]{\sqrt{q^{2} +p^{3} } -q} -\sqrt[{3}]{\sqrt{q^{2} +p^{3} } +q} ,\\ q=&\frac{\Theta ^{3} a_{2} }{\, 24\, a_{4}^{2} } \left(4a_{0} -\Theta \, L_{x} L_{y} L_{z} \frac{a_{1} a_{3} }{a_{4} } \right)-\Theta ^{2} \frac{a_{0} }{8a_{4}^{2} } \left(4\, \Theta \, a_{2} -\Theta ^{2} \frac{a_{3}^{2} }{a_{4} } \right)-\frac{\Theta ^{3} a_{2}^{3} }{54\, a_{4}^{3} } -L_{x}^{2} L_{y}^{2} L_{z}^{2} \frac{\Theta ^{4} a_{1}^{2} }{8\, a_{4}^{2} } ,\\ p=&\Theta ^{2} \frac{4a_{0} a_{4} -\Theta \, L_{x} L_{y} L_{z} a_{1} a_{3} }{12\, a_{4}^{2} } -\frac{\Theta \, a_{2} }{18\, a_{4} } ,\\ \alpha _{1\Phi _{I} } =&\frac{R_{I1} }{\Theta L_{x} L_{y} L_{z} } +\frac{S_{II\, 20} }{\Theta L_{x} L_{y} L_{z} } +\frac{1}{L_{x} L_{y} L_{z} } \int _{0}^{L_{x} }\int _{0}^{L_{y} }\int _{0}^{L_{z} }f_{\Phi _{I} } \, \left(x,y,z\right)\, \, d\, z\, d\, y\, d\, x , \\ \alpha _{1\Phi _{V} } =&\frac{R_{V1} }{\Theta L_{x} L_{y} L_{z} } +\frac{S_{VV\, 20} }{\Theta L_{x} L_{y} L_{z} } +\frac{1}{L_{x} L_{y} L_{z} } \int _{0}^{L_{x} }\int _{0}^{L_{y} }\int _{0}^{L_{z} }f_{\Phi _{V} } \, \left(x,y,z\right)\, \, d\, z\, d\, y\, d\, x , \end{aligned}\] where \(R_{\rho \, \, i} =\int _{0}^{\Theta }\left(\Theta -t\right)\int _{0}^{L_{x} }\int _{0}^{L_{y} }\int _{0}^{L_{z} }k_{I} \left(x,y,z,T\right)\, I_{1}^{i} \left(x,y,z,t\right)\, \, d\, z\, d\, y\, d\, x\, d\, t\).

We determine approximations of the second and higher orders of concentrations of dopant and radiation defects framework standard iterative procedure of method of averaging of function corrections [26]. Framework this procedure to determine approximations of the n-th order of concentrations of dopant and radiation defects we replace the required concentrations in the Eqs. (1c), (3c), (5c) on the following sum \(\alpha\)\({}_{n}\)+\(\rho\) \({}_{n}\)\({}_{-1\ }\)(x,y,z,t). The replacement leads to the following transformation of the appropriate equations

\[\begin{align} \frac{\partial C_2(x,y,z,t)}{\partial t} &= \frac{\partial}{\partial x} \left( \left\{ 1 + \zeta \frac{[\alpha_{2C} + C_1(x,y,z,t)]^{\gamma}}{P^{\gamma}(x,y,z,T)} \right\} \left[ 1 + \zeta_1 \frac{V(x,y,z,t)}{V^*} + \zeta_2 \frac{V^2(x,y,z,t)}{(V^*)^2} \right] D_L(x,y,z,T) \frac{\partial C_1(x,y,z,t)}{\partial x} \right) \\ &\quad + \frac{\partial}{\partial y} \left( \left\{ 1 + \zeta \frac{[\alpha_{2C} + C_1(x,y,z,t)]^{\gamma}}{P^{\gamma}(x,y,z,T)} \right\} \left[ 1 + \zeta_1 \frac{V(x,y,z,t)}{V^*} + \zeta_2 \frac{V^2(x,y,z,t)}{(V^*)^2} \right] D_L(x,y,z,T) \frac{\partial C_1(x,y,z,t)}{\partial y} \right) \\ &\quad + \frac{\partial}{\partial z} \left( \left\{ 1 + \zeta \frac{[\alpha_{2C} + C_1(x,y,z,t)]^{\gamma}}{P^{\gamma}(x,y,z,T)} \right\} \left[ 1 + \zeta_1 \frac{V(x,y,z,t)}{V^*} + \zeta_2 \frac{V^2(x,y,z,t)}{(V^*)^2} \right] D_L(x,y,z,T) \frac{\partial C_1(x,y,z,t)}{\partial z} \right) \\ &\quad + f_C(x,y,z)\delta(t) + \frac{\partial}{\partial x} \left[ \frac{D_{CS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ \frac{D_{CS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial y} \right] + \frac{\partial}{\partial z} \left[ \frac{D_{CS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial z} \right] \\ &\quad + \Omega \frac{\partial}{\partial x} \left\{ \frac{D_S}{k T} \nabla_S \mu_1(x,y,z,t) \int_0^{L_z} [\alpha_{2C} + C(x,y,W,t)] \, dW \right\} \\ &\quad + \Omega \frac{\partial}{\partial y} \left\{ \frac{D_S}{k T} \nabla_S \mu_1(x,y,z,t) \int_0^{L_z} [\alpha_{2C} + C(x,y,W,t)] \, dW \right\} \end{align}\tag{1d}\] \[\begin{align} \frac{\partial I_2(x,y,z,t)}{\partial t} &= \frac{\partial}{\partial x} \left[ D_I(x,y,z,T) \frac{\partial I_1(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ D_I(x,y,z,T) \frac{\partial I_1(x,y,z,t)}{\partial y} \right] \\ &\quad + \frac{\partial}{\partial z} \left[ D_I(x,y,z,T) \frac{\partial I_1(x,y,z,t)}{\partial z} \right] – k_{I,I}(x,y,z,T) [\alpha_{1I} + I_1(x,y,z,t)]^2 – k_{I,V}(x,y,z,T) \\ &\quad \times [\alpha_{1I} + I_1(x,y,z,t)][\alpha_{1V} + V_1(x,y,z,t)] \\ &\quad + \Omega \frac{\partial}{\partial x} \left\{ \nabla_S \mu(x,y,z,t) \int_0^{L_z} [\alpha_{2I} + I_1(x,y,W,t)] \, dW \right\} \\ &\quad + \Omega \frac{\partial}{\partial y} \left\{ \frac{D_{IS}}{kT} \nabla_S \mu(x,y,z,t) \int_0^{L_z} [\alpha_{2I} + I_1(x,y,W,t)] \, dW \right\} \\ &\quad + \frac{\partial}{\partial x} \int_0^t \frac{D_{IS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial x} \, d\tau + \frac{\partial}{\partial y} \int_0^t \frac{D_{IS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial y} \, d\tau \\ &\quad + \frac{\partial}{\partial z} \int_0^t \frac{D_{IS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial z} \, d\tau \end{align}\tag{3d}\] \[\begin{align} \frac{\partial V_2(x,y,z,t)}{\partial t} &= \frac{\partial}{\partial x} \left[ D_V(x,y,z,T) \frac{\partial V_1(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ D_V(x,y,z,T) \frac{\partial V_1(x,y,z,t)}{\partial y} \right] \\ &\quad + \frac{\partial}{\partial z} \left[ D_V(x,y,z,T) \frac{\partial V_1(x,y,z,t)}{\partial z} \right] – k_{V,V}(x,y,z,T) [\alpha_{1V} + V_1(x,y,z,t)]^2 \\ &\quad – k_{I,V}(x,y,z,T) [\alpha_{1I} + I_1(x,y,z,t)][\alpha_{1V} + V_1(x,y,z,t)] \\ &\quad + \Omega \frac{\partial}{\partial x} \left\{ \nabla_S \mu(x,y,z,t) \int_0^{L_z} [\alpha_{2V} + V_1(x,y,W,t)] \, dW \right\} \\ &\quad + \Omega \frac{\partial}{\partial y} \left\{ \frac{D_{VS}}{kT} \nabla_S \mu(x,y,z,t) \int_0^{L_z} [\alpha_{2V} + V_1(x,y,W,t)] \, dW \right\} \\ &\quad + \frac{\partial}{\partial x} \int_0^t \frac{D_{VS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial x} \, d\tau + \frac{\partial}{\partial y} \int_0^t \frac{D_{VS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial y} \, d\tau \\ &\quad + \frac{\partial}{\partial z} \int_0^t \frac{D_{VS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial z} \, d\tau \end{align}\] \[\begin{align} \frac{\partial \Phi_{2I}(x,y,z,t)}{\partial t} &= \frac{\partial}{\partial x} \left[ D_{\Phi_I}(x,y,z,T) \frac{\partial \Phi_{1I}(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ D_{\Phi_I}(x,y,z,T) \frac{\partial \Phi_{1I}(x,y,z,t)}{\partial y} \right] \\ &\quad + \Omega \frac{\partial}{\partial x} \left\{ \frac{D_{\Phi_IS}}{kT} \nabla_S \mu(x,y,z,t) \int_0^{L_z} [\alpha_{2\Phi_I} + \Phi_{1I}(x,y,W,t)] \, dW \right\} + k_{I,I}(x,y,z,T) I^2(x,y,z,t) \\ &\quad + \Omega \frac{\partial}{\partial y} \left\{ \frac{D_{\Phi_IS}}{kT} \nabla_S \mu(x,y,z,t) \int_0^{L_z} [\alpha_{2\Phi_I} + \Phi_{1I}(x,y,W,t)] \, dW \right\} + k_I(x,y,z,T) I(x,y,z,t) \\ &\quad + \frac{\partial}{\partial x} \left[ \frac{D_{\Phi_IS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ \frac{D_{\Phi_IS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial y} \right] + \frac{\partial}{\partial z} \left[ \frac{D_{\Phi_IS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial z} \right] \\ &\quad + \frac{\partial}{\partial z} \left[ D_{\Phi_I}(x,y,z,T) \frac{\partial \Phi_{1I}(x,y,z,t)}{\partial z} \right] + f_{\Phi_I}(x,y,z) \delta(t) \end{align} \tag{5d}\]

\[\begin{aligned} \frac{\partial \Phi_{2V}(x,y,z,t)}{\partial t} &= \frac{\partial}{\partial x} \left[ D_{\Phi_V}(x,y,z,T) \frac{\partial \Phi_{1V}(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ D_{\Phi_V}(x,y,z,T) \frac{\partial \Phi_{1V}(x,y,z,t)}{\partial y} \right] \\ &\quad + \Omega \frac{\partial}{\partial x} \left\{ \frac{D_{\Phi_V S}}{kT} \nabla_S \mu(x,y,z,t) \int_0^{L_z} [\alpha_{2\Phi_V} + \Phi_{1V}(x,y,W,t)] \, dW \right\} + k_{V,V}(x,y,z,T)V^2(x,y,z,t) \\ &\quad + \Omega \frac{\partial}{\partial y} \left\{ \frac{D_{\Phi_V S}}{kT} \nabla_S \mu(x,y,z,t) \int_0^{L_z} [\alpha_{2\Phi_V} + \Phi_{1V}(x,y,W,t)] \, dW \right\} + k_V(x,y,z,T)V(x,y,z,t) \\ &\quad + \frac{\partial}{\partial x} \left[ \frac{D_{\Phi_V S}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ \frac{D_{\Phi_V S}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial y} \right] + \frac{\partial}{\partial z} \left[ \frac{D_{\Phi_V S}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial z} \right] \\ &\quad + \frac{\partial}{\partial z} \left[ D_{\Phi_V}(x,y,z,T) \frac{\partial \Phi_{1V}(x,y,z,t)}{\partial z} \right] + f_{\Phi_V}(x,y,z)\delta(t) \end{aligned}\]

Integration of the left and the right sides of Eqs. (1d), (3d) and (5d) gives us possibility to obtain relations for the required concentrations in the final form

\[\begin{align} C_2(x,y,z,t) &= \frac{\partial}{\partial x} \int_0^t \left\{ \left( 1 + \zeta \frac{[\alpha_{2C} + C_1(x,y,z,\tau)]^{\gamma}}{P^{\gamma}(x,y,z,T)} \right) \left( 1 + \zeta_1 \frac{V(x,y,z,\tau)}{V^*} + \zeta_2 \frac{V^2(x,y,z,\tau)}{(V^*)^2} \right) \right\} \\ &\quad \times D_L(x,y,z,T) \frac{\partial C_1(x,y,z,\tau)}{\partial x} \, d\tau \\ &\quad + \frac{\partial}{\partial y} \int_0^t D_L(x,y,z,T) \left( 1 + \zeta_1 \frac{V(x,y,z,\tau)}{V^*} + \zeta_2 \frac{V^2(x,y,z,\tau)}{(V^*)^2} \right) \frac{\partial C_1(x,y,z,\tau)}{\partial y} \, d\tau \\ &\quad + \frac{\partial}{\partial z} \int_0^t D_L(x,y,z,T) \left( 1 + \zeta_1 \frac{V(x,y,z,\tau)}{V^*} + \zeta_2 \frac{V^2(x,y,z,\tau)}{(V^*)^2} \right) \frac{\partial C_1(x,y,z,\tau)}{\partial z} \\ &\quad \times \left( 1 + \zeta \frac{[\alpha_{2C} + C_1(x,y,z,\tau)]^{\gamma}}{P^{\gamma}(x,y,z,T)} \right) d\tau + f_C(x,y,z) \\ &\quad + \Omega \frac{\partial}{\partial x} \int_0^t \frac{D_S}{kT} \nabla_S \mu(x,y,z,\tau) \int_0^{L_z} [\alpha_{2C} + C_1(x,y,W,t)] \, dW \, d\tau \\ &\quad + \Omega \frac{\partial}{\partial y} \int_0^t \frac{D_S}{kT} \nabla_S \mu(x,y,z,\tau) \int_0^{L_z} [\alpha_{2C} + C_1(x,y,W,t)] \, dW \, d\tau \\ &\quad + \frac{\partial}{\partial x} \left[ \frac{D_{CS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ \frac{D_{CS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial y} \right] + \frac{\partial}{\partial z} \left[ \frac{D_{CS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial z} \right] \end{align}\tag{1e}\] \[\begin{align} I_2(x,y,z,t) &= \frac{\partial}{\partial x} \int_0^t D_I(x,y,z,T) \frac{\partial I_1(x,y,z,\tau)}{\partial x} \, d\tau + \frac{\partial}{\partial y} \int_0^t D_I(x,y,z,T) \frac{\partial I_1(x,y,z,\tau)}{\partial y} \, d\tau \\ &\quad + \frac{\partial}{\partial z} \int_0^t D_I(x,y,z,T) \frac{\partial I_1(x,y,z,\tau)}{\partial z} \, d\tau – \int_0^t k_{I,I}(x,y,z,T) [\alpha_{2I} + I_1(x,y,z,\tau)]^2 \, d\tau \\ &\quad – \int_0^t k_{I,V}(x,y,z,T) [\alpha_{2I} + I_1(x,y,z,\tau)][\alpha_{2V} + V_1(x,y,z,\tau)] \, d\tau \\ &\quad + \frac{\partial}{\partial x} \int_0^t \nabla_S \mu(x,y,z,\tau) \times \Omega \frac{D_{IS}}{kT} \int_0^{L_z} [\alpha_{2I} + I_1(x,y,W,\tau)] \, dW \, d\tau \\ &\quad + \frac{\partial}{\partial y} \int_0^t \nabla_S \mu(x,y,z,\tau) \times \Omega \frac{D_{IS}}{kT} \int_0^{L_z} [\alpha_{2I} + I_1(x,y,W,\tau)] \, dW \, d\tau \\ &\quad + \frac{\partial}{\partial x} \left[ \frac{D_{IS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ \frac{D_{IS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial y} \right] + \frac{\partial}{\partial z} \left[ \frac{D_{IS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial z} \right] \\ &\quad + f_I(x,y,z) \end{align}\tag{3e}\] \[\begin{align} V_2(x,y,z,t) &= \frac{\partial}{\partial x} \int_0^t D_V(x,y,z,T) \frac{\partial V_1(x,y,z,\tau)}{\partial x} \, d\tau + \frac{\partial}{\partial y} \int_0^t D_V(x,y,z,T) \frac{\partial V_1(x,y,z,\tau)}{\partial y} \, d\tau \\ &\quad + \frac{\partial}{\partial z} \int_0^t D_V(x,y,z,T) \frac{\partial V_1(x,y,z,\tau)}{\partial z} \, d\tau \\ &\quad – \int_0^t k_{V,V}(x,y,z,T) [\alpha_{2V} + V_1(x,y,z,\tau)]^2 \, d\tau \\ &\quad – \int_0^t k_{I,V}(x,y,z,T) [\alpha_{2I} + I_1(x,y,z,\tau)][\alpha_{2V} + V_1(x,y,z,\tau)] \, d\tau \\ &\quad + \Omega \frac{\partial}{\partial x} \int_0^t \nabla_S \mu(x,y,z,\tau) \int_0^{L_z} [\alpha_{2V} + V_1(x,y,W,\tau)] \, dW \, d\tau \\ &\quad + \Omega \frac{\partial}{\partial y} \int_0^t \nabla_S \mu(x,y,z,\tau) \int_0^{L_z} [\alpha_{2V} + V_1(x,y,W,\tau)] \, dW \, d\tau \\ &\quad + \frac{\partial}{\partial x} \left[ \frac{D_{VS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial x} \right] + \frac{\partial}{\partial y} \left[ \frac{D_{VS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial y} \right] + \frac{\partial}{\partial z} \left[ \frac{D_{VS}}{V k T} \frac{\partial \mu_2(x,y,z,t)}{\partial z} \right] \\ &\quad + f_V(x,y,z) \end{align}\] \[\begin{align} \Phi_{2I}(x,y,z,t) &= \frac{\partial}{\partial x} \int_0^t D_{\Phi_I}(x,y,z,T) \frac{\partial \Phi_{1I}(x,y,z,\tau)}{\partial x} \, d\tau + \frac{\partial}{\partial y} \int_0^t D_{\Phi_I}(x,y,z,T) \frac{\partial \Phi_{1I}(x,y,z,\tau)}{\partial y} \, d\tau \\ &\quad + \frac{\partial}{\partial z} \int_0^t D_{\Phi_I}(x,y,z,T) \frac{\partial \Phi_{1I}(x,y,z,\tau)}{\partial z} \, d\tau \\ &\quad + \Omega \frac{\partial}{\partial x} \int_0^t \nabla_S \mu(x,y,z,\tau) \frac{D_{\Phi_I S}}{kT} \int_0^{L_z} [\alpha_{2\Phi_I} + \Phi_{1I}(x,y,W,\tau)] \, dW \, d\tau \\ &\quad + \Omega \frac{\partial}{\partial y} \int_0^t \nabla_S \mu(x,y,z,\tau) \frac{D_{\Phi_I S}}{kT} \int_0^{L_z} [\alpha_{2\Phi_I} + \Phi_{1I}(x,y,W,\tau)] \, dW \, d\tau \\ &\quad + \int_0^t k_{I,I}(x,y,z,T) I^2(x,y,z,\tau) \, d\tau + \int_0^t k_I(x,y,z,T) I(x,y,z,\tau) \, d\tau \\ &\quad + \frac{\partial}{\partial x} \int_0^t \frac{D_{\Phi_I S}}{V k T} \frac{\partial \mu_2(x,y,z,\tau)}{\partial x} \, d\tau + \frac{\partial}{\partial y} \int_0^t \frac{D_{\Phi_I S}}{V k T} \frac{\partial \mu_2(x,y,z,\tau)}{\partial y} \, d\tau \\ &\quad + \frac{\partial}{\partial z} \int_0^t \frac{D_{\Phi_I S}}{V k T} \frac{\partial \mu_2(x,y,z,\tau)}{\partial z} \, d\tau + f_{\Phi_I}(x,y,z) \end{align} \tag{5e}\]

\[\begin{aligned} \Phi_{2V}(x,y,z,t) &= \frac{\partial}{\partial x} \int_0^t D_{\Phi_V}(x,y,z,T) \frac{\partial \Phi_{1V}(x,y,z,\tau)}{\partial x} \, d\tau + \frac{\partial}{\partial y} \int_0^t D_{\Phi_V}(x,y,z,T) \frac{\partial \Phi_{1V}(x,y,z,\tau)}{\partial y} \, d\tau \\ &\quad + \frac{\partial}{\partial z} \int_0^t D_{\Phi_V}(x,y,z,T) \frac{\partial \Phi_{1V}(x,y,z,\tau)}{\partial z} \, d\tau \\ &\quad + \Omega \frac{\partial}{\partial x} \int_0^t \nabla_S \mu(x,y,z,\tau) \frac{D_{\Phi_V S}}{kT} \int_0^{L_z} [\alpha_{2\Phi_V} + \Phi_{1V}(x,y,W,\tau)] \, dW \, d\tau \\ &\quad + \Omega \frac{\partial}{\partial y} \int_0^t \nabla_S \mu(x,y,z,\tau) \frac{D_{\Phi_V S}}{kT} \int_0^{L_z} [\alpha_{2\Phi_V} + \Phi_{1V}(x,y,W,\tau)] \, dW \, d\tau \\ &\quad + \int_0^t k_{V,V}(x,y,z,T) V^2(x,y,z,\tau) \, d\tau + \int_0^t k_V(x,y,z,T) V(x,y,z,\tau) \, d\tau \\ &\quad + \frac{\partial}{\partial x} \int_0^t \frac{D_{\Phi_V S}}{V k T} \frac{\partial \mu_2(x,y,z,\tau)}{\partial x} \, d\tau + \frac{\partial}{\partial y} \int_0^t \frac{D_{\Phi_V S}}{V k T} \frac{\partial \mu_2(x,y,z,\tau)}{\partial y} \, d\tau \\ &\quad + \frac{\partial}{\partial z} \int_0^t \frac{D_{\Phi_V S}}{V k T} \frac{\partial \mu_2(x,y,z,\tau)}{\partial z} \, d\tau + f_{\Phi_V}(x,y,z) \end{aligned}\]

Average values of the second-order approximations of required approximations by using the following standard relation [26]; \[\label{GrindEQ__10_} \alpha _{2\rho } =\frac{1}{\Theta L_{x} L_{y} L_{z} } \int _{0}^{\Theta }\int _{0}^{L_{x} }\int _{0}^{L_{y} }\int _{0}^{L_{z} }\left[\rho _{2} \left(x,y,z,t\right)-\rho _{1} \left(x,y,z,t\right)\right]\, \, d\, z\, d\, y\, d\, x\, d\, t . \tag{10}\]

Substitution of the relations (1e), (3e), (5e) into relation (10) gives us possibility to obtain relations for required average values \(\alpha\)\({}_{\ }\)\({}_{2}\)

\[\begin{aligned} \alpha2C=&0,\\ \alpha2I\ =&0, \\ \alpha2V\ =&0,\\ \alpha _{2V} =&\sqrt{\frac{\left(b_{3} +E\right)^{2} }{4\, b_{4}^{2} } -4\, \, \left(F+\frac{\Theta \, a_{3} F+\Theta ^{2} L_{x} L_{y} L_{z} b_{1} }{b_{4} } \right)} -\frac{b_{3} +E}{4\, b_{4} } ,\\ \alpha _{2I} =&\frac{C_{V} -\alpha _{2V}^{2} S_{VV00} -\alpha _{2V} \left(2\, S_{VV01} +S_{I\, V10} +\Theta L_{x} L_{y} L_{z} \right)\, \, -S_{V\, V02} -S_{I\, V11} }{S_{IV01} +\alpha _{2V} S_{IV00} } , \end{aligned}\] where \[\begin{aligned} b_{4} =&\frac{1}{\Theta L_{x} L_{y} L_{z} } S_{IV00}^{2} S_{V\, V00} -\frac{1}{\Theta L_{x} L_{y} L_{z} } S_{VV00}^{2} S_{II00},\\ b_{3} =&-\frac{S_{II00} S_{VV00} }{\Theta L_{x} L_{y} L_{z} } \left(2S_{VV01} +S_{IV10} +\Theta L_{x} L_{y} L_{z} \right)+\frac{S_{IV00} S_{VV00} }{\Theta L_{x} L_{y} L_{z} } \left(S_{IV01} +2S_{II10} +S_{IV01} +\Theta L_{x} L_{y} L_{z} \right)\\ &+\frac{S_{IV00}^{2} }{\Theta L_{x} L_{y} L_{z} } \left(2S_{VV01} +S_{IV10} +\Theta L_{x} L_{y} L_{z} \right)-\frac{S_{IV00}^{2} S_{IV10} }{\Theta ^{3} L_{x}^{3} L_{y}^{3} L_{z}^{3} } ,\\ b_{2} =&\frac{S_{II00} S_{VV00} }{\Theta L_{x} L_{y} L_{z} } \left(S_{VV02} +S_{IV11} +C_{V} \right)-\left(S_{IV10} -2S_{VV01} +\Theta L_{x} L_{y} \times L_{z} \right)^{2}\\ &+\frac{S_{IV01} S_{VV00} }{\Theta L_{x} L_{y} L_{z} } \left(\Theta L_{x} L_{y} L_{z} +2S_{II10} +S_{IV01} \right)\\ &+\frac{S_{IV00} }{\Theta L_{x} L_{y} L_{z} } \left(S_{IV01} +2S_{II10} +2S_{IV01} +\Theta L_{x} L_{y} \times L_{z} \right)\, \left(2S_{VV01} +\Theta L_{x} L_{y} L_{z} +S_{IV10} \right)\\ &-\frac{S_{IV00}^{2} }{\Theta L_{x} L_{y} L_{z} } \left(C_{V} -S_{VV02} -S_{IV11} \right)+\frac{C_{I} S_{IV00}^{2} }{\Theta ^{2} L_{x}^{2} L_{y}^{2} L_{z}^{2} } -\frac{2\, S_{IV10} }{\Theta L_{x} L_{y} L_{z} } \times S_{IV00} S_{IV01} ,\\ b_{1} =&S_{II00} \frac{S_{IV11} +S_{VV02} +C_{V} }{\Theta L_{x} L_{y} L_{z} } \left(2S_{VV01} +S_{IV10} +\Theta L_{x} L_{y} L_{z} \right)\\ &+\frac{S_{IV01} }{\Theta L_{x} L_{y} L_{z} } \left(\Theta L_{x} L_{y} \times L_{z} +2S_{II10} +S_{IV01} \right)\, \left(2S_{VV01} +S_{IV10} +\Theta LL_{y} L_{z} \right)-\frac{S_{IV10} S_{IV01}^{2} }{\Theta L_{x} L_{y} L_{z} } \\ &-\frac{S_{IV00} }{\Theta L_{x} L_{y} L_{z} } \left(3S_{IV01} +2S_{II10} +\Theta L_{x} L_{y} L_{z} \right)\, \left(C_{V} -S_{VV02} -S_{IV11} \right)+2\, C_{I} S_{IV00} S_{IV01} ,\\ b_{0} =&\frac{S_{II00} }{\Theta L_{x} L_{y} L_{z} } \left(S_{IV00} +S_{VV02} \right)^{2} -\frac{S_{IV01} }{L_{x} L_{y} L_{z} } \times \frac{1}{\Theta } \left(\Theta L_{x} L_{y} L_{z} +2\, S_{II10} +S_{IV01} \right)\, \left(C_{V} -S_{VV02} -S_{IV11} \right)\\ &+2C_{I} S_{IV01}^{2} -S_{IV01} \frac{C_{V} -S_{VV02} -S_{IV11} }{\Theta L_{x} L_{y} L_{z} } \times \frac{1}{\Theta } \left(\Theta L_{x} L_{y} L_{z} +2\, S_{II10} +S_{IV01} \right)\, \left(C_{V} -S_{VV02} -S_{IV11} \right)\\ &+2C_{I} S_{IV01}^{2} -S_{IV01} \frac{C_{V} -S_{VV02} -S_{IV11} }{\Theta L_{x} L_{y} L_{z} } \times \, S_{IV01} \left(\Theta L_{x} L_{y} L_{z} +2\, S_{II10} +S_{IV01} \right),\\ C_{I} =&\frac{\alpha _{1I} \alpha _{1V} }{\Theta L_{x} L_{y} L_{z} } S_{IV00} +\frac{\alpha _{1I}^{2} S_{II00} }{\Theta L_{x} L_{y} L_{z} } -\frac{S_{II20} S_{II20} }{\Theta L_{x} L_{y} L_{z} } -\frac{S_{IV11} }{\Theta L_{x} L_{y} L_{z} } ,\\ C_{V} =&\alpha _{1I} \alpha _{1V} S_{IV00} +\alpha _{1V}^{2} S_{VV00} -S_{VV02} -S_{IV11} ,\\ E=&\sqrt{8\, y+\Theta ^{2} \frac{a_{3}^{2} }{a_{4}^{2} } -4\, \Theta \, \frac{a_{2} }{a_{4} } } ,\\ F=&\frac{\Theta \, a_{2} }{6\, a_{4} } +\sqrt[{3}]{\sqrt{r^{2} +s^{3} } -r} -\sqrt[{3}]{\sqrt{r^{2} +s^{3} } +r} ,\\ r=&\frac{\Theta ^{3} b_{2} }{\, 24\, b_{4}^{2} } \left(4b_{0} -\Theta \, L_{x} L_{y} L_{z} \frac{b_{1} b_{3} }{b_{4} } \right)-\frac{\Theta ^{3} b_{2}^{3} }{54\, b_{4}^{3} } -b_{0} \frac{\Theta ^{2} }{8b_{4}^{2} } \times \left(4\, \Theta \, b_{2} -\Theta ^{2} \frac{b_{3}^{2} }{b_{4} } \right)-L_{x}^{2} L_{y}^{2} L_{z}^{2} \frac{\Theta ^{4} b_{1}^{2} }{8\, b_{4}^{2} } ,\\ s=&\Theta ^{2} \frac{4b_{0} b_{4} -\Theta \, L_{x} L_{y} L_{z} b_{1} b_{3} }{12\, b_{4}^{2} } -\frac{\Theta \, b_{2} }{18\, b_{4} } . \end{aligned}\]

Farther we determine solutions of Eq. (8), i.e. components of displacement vector. To determine the first-order approximations of the considered components framework method of averaging of function corrections we replace the required functions in the right sides of the equations by their not yet known average values \(\alpha\)\({}_{i}\). The substitution leads to the following result \[\begin{aligned} \rho \, \left(z\right)\frac{\partial ^{2} \, u_{1x} \left(x,y,z,t\right)}{\partial \, t^{2} } =&-K\, \left(z\right)\, \beta \, \left(z\right)\frac{\partial \, T\, \left(x,y,z,t\right)}{\partial \, x} , \rho ,\\ \left(z\right)\frac{\partial ^{2} \, u_{1y} \left(x,y,z,t\right)}{\partial \, t^{2} } =&\rho \, \left(z\right)\frac{\partial ^{2} \, u_{1y} \left(x,y,z,t\right)}\partial ,\\ t^{2} =&-K\, \left(z\right)\, \beta \, \left(z\right)\frac{\partial \, T\, \left(x,y,z,t\right)}{\partial \, y} , \\ \rho \left(z\right)\frac{\partial ^{2} \, u_{1z} \left(x,y,z,t\right)}{\partial \, t^{2} } =&-K\, \left(z\right)\, \beta \, \left(z\right)\frac{\partial \, T\, \left(x,y,z,t\right)}{\partial \, z} . \end{aligned}\]

Integration of the left and the right sides of the above relations on time t leads to the following result \[\begin{aligned} u_{1x} \left(x,y,z,t\right)=&u_{0x} +K\, \left(z\right)\, \frac{\beta \, \left(z\right)}{\rho \, \left(z\right)} \frac{\partial }{\partial \, x} \int _{0}^{t}\int _{0}^{\vartheta }\, T\, \left(x,y,z,\tau \right)\, \, d\, \tau \, d\, \vartheta \\ & -K\, \left(z\right)\, \frac{\beta \, \left(z\right)}{\rho \, \left(z\right)} \frac{\partial }{\partial \, x} \int _{0}^{\infty }\int _{0}^{\vartheta }\, T\, \left(x,y,z,\tau \right)\, \, d\, \tau \, d\, \vartheta ,\\ u_{1y} \left(x,y,z,t\right)=&u_{0y} +K\, \left(z\right)\, \frac{\beta \, \left(z\right)}{\rho \, \left(z\right)} \frac{\partial }{\partial \, y} \int _{0}^{t}\int _{0}^{\vartheta }\, T\, \left(x,y,z,\tau \right)\, \, d\, \tau \, d\, \vartheta \\ &-K\, \left(z\right)\, \frac{\beta \, \left(z\right)}{\rho \, \left(z\right)} \frac{\partial }{\partial \, y} \int _{0}^{\infty }\int _{0}^{\vartheta }\, T\, \left(x,y,z,\tau \right)\, \, d\, \tau \, d\, \vartheta ,\\ u_{1z} \left(x,y,z,t\right)=&u_{0z} +K\, \left(z\right)\, \frac{\beta \, \left(z\right)}{\rho \, \left(z\right)} \frac{\partial }{\partial \, z} \int _{0}^{t}\int _{0}^{\vartheta }\, T\, \left(x,y,z,\tau \right)\, \, d\, \tau \, d\, \vartheta \\ &-K\, \left(z\right)\, \frac{\beta \, \left(z\right)}{\rho \, \left(z\right)} \frac{\partial }{\partial \, z} \int _{0}^{\infty }\int _{0}^{\vartheta }\, T\, \left(x,y,z,\tau \right)\, \, d\, \tau \, d\, \vartheta . \end{aligned}\]

Approximations of the second and higher orders of components of displacement vector could be determined by using standard replacement of the required components on the following sums \(\alpha\)\({}_{i}\)+u\({}_{i}\)(x,y,z,t) [26]. The replacement leads to the following result \[\begin{aligned} \rho \left(z\right)\frac{\partial ^{2} u_{2x} \left(x,y,z,t\right)}{\partial t^{2} } =&\left\{K \left(z\right)+\frac{5E \left(z\right)}{6 \left[1+\sigma \left(z\right)\right]} \right\}\frac{\partial ^{2} u_{1x} \left(x,y,z,t\right)}{\partial x^{2} } +\left\{K \left(z\right)-\frac{E \left(z\right)}{3 \left[1+\sigma \left(z\right)\right]} \right\}\frac{\partial ^{2} u_{1y} \left(x,y,z,t\right)}{\partial x \partial y} \\ &\times \frac{\partial ^{2} u_{1y} \left(x,y,z,t\right)}{\partial x \partial y} +\frac{E \left(z\right)}{2 \left[1+\sigma \left(z\right)\right]} \left[\frac{\partial ^{2} u_{1y} \left(x,y,z,t\right)}{\partial y^{2} } +\frac{\partial ^{2} u_{1z} \left(x,y,z,t\right)}{\partial z^{2} } \right]\\ &-\frac{\partial T \left(x,y,z,t\right)}{\partial x} \times K\left(z\right) \beta \left(z\right)+\left\{K \left(z\right)+\frac{E \left(z\right)}{3 \left[ 1+\sigma \left(z\right)\right]} \right\}\frac{\partial ^{2} u_{1z} \left(x,y,z,t\right)}{\partial x \partial z} \\ \rho \left(z\right)\frac{\partial ^{2} u_{2y} \left(x,y,z,t\right)}{\partial t^{2} } =&\frac{E \left(z\right)}{2\left[1+\sigma \left(z\right)\right]} \left[\frac{\partial ^{2} u_{1y} \left(x,y,z,t\right)}{\partial x^{2} } +\frac{\partial ^{2} u_{1x} \left(x,y,z,t\right)}{\partial x \partial y} \right]- \frac{\partial T \left(x,y,z,t\right)}{\partial y} \\ &\times K \left(z\right) \beta \left(z\right)+\frac{\partial }{\partial z} \left\{\frac{E \left(z\right)}{2 \left[1+\sigma \left(z\right)\right]} \left[\frac{\partial u_{1y} \left(x,y,z,t\right)}{\partial z} +\frac{\partial u_{1z} \left(x,y,z,t\right)}{\partial y} \right]\right\}\\ &+\frac{\partial ^{2} u_{1y} \left(x,y,z,t\right)}{\partial y^{2} } \times \left\{\frac{5E \left(z\right)}{12 \left[1+\sigma \left(z\right)\right]} +K \left(z\right)\right\}+\left\{K \left(z\right)-\frac{E \left(z\right)}{6 \left[1+\sigma \left(z\right)\right]} \right\}\frac{\partial ^{2} u_{1y} \left(x,y,z,t\right)}{\partial y \partial z} \\ &+K \left(z\right)\frac{\partial ^{2} u_{1y} \left(x,y,z,t\right)}{\partial x \partial y} \\ \rho \left(z\right)\frac{\partial ^{2} u_{2z} \left(x,y,z,t\right)}{\partial t^{2} } =&\frac{E \left(z\right)}{2 \left[1+\sigma \left(z\right)\right]} \left[\frac{\partial ^{2} u_{1z} \left(x,y,z,t\right)}{\partial x^{2} } +\frac{\partial ^{2} u_{1z} \left(x,y,z,t\right)}{\partial y^{2} } +\frac{\partial ^{2} u_{1x} \left(x,y,z,t\right)}{\partial x \partial z} + \frac{\partial ^{2} u_{1y} \left(x,y,z,t\right)}{\partial y \partial z} \right]\\ &+\frac{\partial }{\partial z} \left\{K \left(z\right) \left[\frac{\partial u_{1x} \left(x,y,z,t\right)}{\partial x} +\frac{\partial u_{1y} \left(x,y,z,t\right)}{\partial y} +\frac{\partial u_{1x} \left(x,y,z,t\right)}{\partial z} \right]\right.\\ &+\frac{E \left(z\right)}{6 \left[1+\sigma \left(z\right)\right]} \frac{\partial }{\partial z} \left[6\frac{\partial u_{1z} \left(x,y,z,t\right)}{\partial z} -\frac{\partial u_{1x} \left(x,y,z,t\right)}{\partial x} -\frac{\partial u_{1y} \left(x,y,z,t\right)}{\partial y} -\frac{\partial u_{1z} \left(x,y,z,t\right)}{\partial z} \right]\\ &\left. \left. -\frac{\partial u_{1x} \left(x,y,z,t\right)}{\partial x} -\frac{\partial u_{1y} \left(x,y,z,t\right)}{\partial y} -\frac{\partial u_{1z} \left(x,y,z,t\right)}{\partial z} \right]\right\}\frac{E \left(z\right)}{1+\sigma \left(z\right)}\\ & -K \left(z\right) \beta \left(z\right)\frac{\partial T \left(x,y,z,t\right)}{\partial z} . \end{aligned}\]

Integration of the left and right sides of the above relations on time t leads to the following result

\[\begin{aligned} u_{2x}(x,y,z,t) &= \frac{1}{\rho(z)} \left\{ \left( K(z) + \frac{5E(z)}{6[1+\sigma(z)]} \right) \frac{\partial^2}{\partial x^2} \int_0^t \int_0^\vartheta u_{1x}(x,y,z,\tau) \, d\tau d\vartheta \right. \\ &\quad + \frac{1}{\rho(z)} \left( K(z) – \frac{E(z)}{3[1+\sigma(z)]} \right) \frac{\partial^2}{\partial x \partial y} \int_0^t \int_0^\vartheta u_{1y}(x,y,z,\tau) \, d\tau d\vartheta \\ &\quad + \frac{E(z)}{2\rho(z)[1+\sigma(z)]} \frac{\partial^2}{\partial y^2} \int_0^t \int_0^\vartheta u_{1y}(x,y,z,\tau) \, d\tau d\vartheta + \frac{1}{1+\sigma(z)} \frac{\partial^2}{\partial z^2} \int_0^t \int_0^\vartheta u_{1z}(x,y,z,\tau) \, d\tau d\vartheta \\ &\quad + \frac{1}{\rho(z)} \frac{\partial^2}{\partial x \partial z} \int_0^t \int_0^\vartheta u_{1z}(x,y,z,\tau) \, d\tau d\vartheta \\ &\quad + \left( K(z) + \frac{E(z)}{3[1+\sigma(z)]} \right) \beta(z) \frac{1}{\rho(z)} \frac{\partial}{\partial x} \int_0^t \int_0^\vartheta T(x,y,z,\tau) \, d\tau d\vartheta \\ &\quad – \frac{\partial^2}{\partial x^2} \int_0^\infty \int_0^\vartheta u_{1x}(x,y,z,\tau) \, d\tau d\vartheta \left( \frac{K(z)+\frac{5E(z)}{6[1+\sigma(z)]}}{\rho(z)} \right) \\ &\quad – \frac{\partial^2}{\partial x \partial y} \int_0^\infty \int_0^\vartheta u_{1y}(x,y,z,\tau) \, d\tau d\vartheta \left( K(z) – \frac{E(z)}{3[1+\sigma(z)]} \right) \\ &\quad – \frac{1}{2\rho(z)} \frac{E(z)}{1+\sigma(z)} \left[ \frac{\partial^2}{\partial y^2} \int_0^\infty \int_0^\vartheta u_{1y}(x,y,z,\tau) \, d\tau d\vartheta + \frac{\partial^2}{\partial z^2} \int_0^\infty \int_0^\vartheta u_{1z}(x,y,z,\tau) \, d\tau d\vartheta \right] \\ &\quad – \frac{1}{\rho(z)} \left( K(z) + \frac{E(z)}{3[1+\sigma(z)]} \right) \frac{\partial^2}{\partial x \partial z} \int_0^\infty \int_0^\vartheta u_{1z}(x,y,z,\tau) \, d\tau d\vartheta \\ &\quad \left. + u_{0x} + K(z) \frac{\beta(z)}{\rho(z)} \frac{\partial}{\partial x} \int_0^\vartheta \int_0^\vartheta T(x,y,z,\tau) \, d\tau d\vartheta \right\} \end{aligned}\]

\[\begin{aligned} u_{2y}(x,y,z,t) &= \frac{E(z)}{2 \rho(z) [1+\sigma(z)]} \left[ \frac{\partial^2}{\partial x^2} \int_0^t \int_0^\vartheta u_{1x}(x,y,z,\tau) \, d\tau d\vartheta + \frac{\partial^2}{\partial x \partial y} \int_0^t \int_0^\vartheta u_{1x}(x,y,z,\tau) \, d\tau d\vartheta \right] \\ &\quad + \left( \frac{1}{1+\sigma(z)} + \frac{K(z)}{\rho(z)} \right) \frac{\partial^2}{\partial y^2} \int_0^t \int_0^\vartheta u_{1y}(x,y,z,\tau) \, d\tau d\vartheta \\ &\quad + \frac{1}{\rho(z)} \left\{ \left( \frac{5E(z)}{12[1+\sigma(z)]} + K(z) \right) \frac{\partial^2}{\partial y^2} \int_0^t \int_0^\vartheta u_{1y}(x,y,z,\tau) \, d\tau d\vartheta \right\} \\ &\quad + \frac{1}{2\rho(z)} \frac{\partial}{\partial z} \left( \frac{E(z)}{1+\sigma(z)} \left[ \frac{\partial}{\partial z} \int_0^t \int_0^\vartheta u_{1y}(x,y,z,\tau) \, d\tau d\vartheta \right] \right) \\ &\quad + \frac{1}{2\rho(z)} \frac{\partial}{\partial z} \left( \frac{\partial}{\partial y} \int_0^t \int_0^\vartheta u_{1z}(x,y,z,\tau) \, d\tau d\vartheta \right) – K(z) \frac{\beta(z)}{\rho(z)} \int_0^t \int_0^\vartheta T(x,y,z,\tau) \, d\tau d\vartheta \\ &\quad – \left( \frac{E(z)}{6[1+\sigma(z)]} \right) \left[ \frac{\partial^2}{\partial x^2} \int_0^\infty \int_0^\vartheta u_{1x}(x,y,z,\tau) \, d\tau d\vartheta + \frac{\partial^2}{\partial x \partial y} \int_0^\infty \int_0^\vartheta u_{1y}(x,y,z,\tau) \, d\tau d\vartheta \right] \\ &\quad – K(z) \frac{1}{\rho(z)} \left[ \frac{\partial^2}{\partial y \partial z} \int_0^\infty \int_0^\vartheta u_{1y}(x,y,z,\tau) \, d\tau d\vartheta \right] + \frac{5E(z)}{12[1+\sigma(z)]} \frac{\partial}{\partial y} \int_0^\infty \int_0^\vartheta u_{1z}(x,y,z,\tau) \, d\tau d\vartheta \\ &\quad + u_{0y} \end{aligned}\]

\[\begin{aligned} u_z(x,y,z,t) &= \frac{E(z)}{2[1+\sigma(z)]} \left[ \frac{\partial^2}{\partial x^2} \int_0^\infty \int_0^\vartheta u_{1z}(x,y,z,\tau) \, d\tau d\vartheta + \frac{\partial^2}{\partial y^2} \int_0^\infty \int_0^\vartheta u_{1z}(x,y,z,\tau) \, d\tau d\vartheta \right] \\ &\quad + \frac{\partial^2}{\partial x \partial z} \int_0^\infty \int_0^\vartheta u_{1x}(x,y,z,\tau) \, d\tau d\vartheta + \frac{\partial^2}{\partial y \partial z} \int_0^\infty \int_0^\vartheta u_{1y}(x,y,z,\tau) \, d\tau d\vartheta \\ &\quad + \frac{1}{\rho(z)} \frac{\partial}{\partial z} \left\{ K(z) \left[ \frac{\partial}{\partial x} \int_0^\infty \int_0^\vartheta u_{1x}(x,y,z,\tau) \, d\tau d\vartheta + \frac{\partial}{\partial y} \int_0^\infty \int_0^\vartheta u_{1x}(x,y,z,\tau) \, d\tau d\vartheta \right] \right\} \\ &\quad + \frac{1}{\rho(z)} \frac{\partial}{\partial z} \left\{ \frac{1}{6} \frac{E(z)}{1+\sigma(z)} \left[ 6 \frac{\partial}{\partial z} \int_0^\infty \int_0^\vartheta u_{1z}(x,y,z,\tau) \, d\tau d\vartheta \right] \right\} \\ &\quad – \frac{\partial}{\partial x} \int_0^\infty \int_0^\vartheta u_{1x}(x,y,z,\tau) \, d\tau d\vartheta – \frac{\partial}{\partial y} \int_0^\infty \int_0^\vartheta u_{1y}(x,y,z,\tau) \, d\tau d\vartheta \\ &\quad – \frac{\partial}{\partial z} \int_0^\infty \int_0^\vartheta u_{1z}(x,y,z,\tau) \, d\tau d\vartheta – K(z) \frac{\beta(z)}{\rho(z)} \frac{\partial}{\partial z} \int_0^\infty \int_0^\vartheta T(x,y,z,\tau) \, d\tau d\vartheta + u_{0z}. \end{aligned}\]

In the framework of this paper we determine concentration of dopant, concentrations of radiation defects and components of displacement vector by using the second-order approximation in the framework of the method of averaging of function corrections. This approximation is usually enough good approximation to make qualitative analysis and to obtain some quantitative results. All obtained results have been checked by comparison with results of numerical simulations.

2. Discussion

In this section we analyzed dynamics of redistributions of dopant and radiation defects during annealing and under influence of mismatch-induced stress and modification of porosity. Typical distributions of concentrations of dopant in heterostructures are presented on Figures 2 and 3 for diffusion and ion types of doping, respectively. These distributions have been calculated for the case, when value of dopant diffusion coefficient in doped area is larger, than in nearest areas. The figures show, that inhomogeneity of heterostructure gives us possibility to increase compactness of concentrations of dopants and at the same time to increase homogeneity of dopant distribution in doped part of epitaxial layer. However framework this approach of manufacturing of bipolar transistor it is necessary to optimize annealing of dopant and/or radiation defects. Reason of this optimization is following. If annealing time is small, the dopant did not achieve any interfaces between materials of heterostructure. In this situation one cannot find any modifications of distribution of concentration of dopant. If annealing time is large, distribution of concentration of dopant is too homogenous. We optimize annealing time framework recently introduces approach [27– 35]. Framework this criterion we approximate real distribution of concentration of dopant by step-wise function (see Figures 4 and 5). Farther we determine optimal values of annealing time by minimization of the following mean-squared error \[\label{GrindEQ__15_} U=\frac{1}{L_{x} L_{y} L_{z} } \int _{0}^{L_{x} }\int _{0}^{L_{y} }\int _{0}^{L_{z} }\left[C\, \left(x,y,z,\Theta \right)-\psi \, \left(x,y,z\right)\right]\, \, d\, z\, d\, y\, d\, x , \tag{11}\] where \(\psi\) (x,y,z) is the approximation function. Dependences of optimal values of annealing time on parameters are presented on Figures 6 and 7 for diffusion and ion types of doping, respectively. It should be noted, that it is necessary to anneal radiation defects after ion implantation. One could find spreading of concentration of distribution of dopant during this annealing. In the ideal case distribution of dopant achieves appropriate interfaces between materials of heterostructure during annealing of radiation defects. If dopant did not achieves any interfaces during annealing of radiation defects, it is practicably to additionally anneal the dopant. In this situation optimal value of additional annealing time of implanted dopant is smaller, than annealing time of infused dopant.

Farther we analyzed influence of relaxation of mechanical stress on distribution of dopant in doped areas of heterostructure. Under following condition \(\epsilon\)\({}_{0}\)\(\mathrm{<}\) 0 one can find compression of distribution of concentration of dopant near interface between materials of heterostructure. Contrary (at \(\epsilon\)\({}_{0}\)\(\mathrm{>}\)0) one can find spreading of distribution of concentration of dopant in this area. This changing of distribution of concentration of dopant could be at least partially compensated by using laser annealing [35]. This type of annealing gives us possibility to accelerate diffusion of dopant and another processes in annealed area due to inhomogenous distribution of temperature and Arrhenius law. Accounting relaxation of mismatch-induced stress in heterostructure could leads to changing of optimal values of annealing time. At the same time modification of porosity gives us possibility to decrease value of mechanical stress. On the one hand mismatch-induced stress could be used to increase density of elements of integrated circuits. On the other hand could leads to generation dislocations of the discrepancy. Figures 8 and 9 show distributions of concentration of vacancies in porous materials and component of displacement vector, which is perpendicular to interface between layers of heterostructure.

3. Conclusion

In this paper we model redistribution of infused and implanted dopants with account relaxation mismatch-induced stress during manufacturing field-effect heterotransistors in the framework of the three-stage differential amplifier. We formulate recommendations for optimization of annealing to decrease dimensions of transistors and to increase their density. We formulate recommendations to decrease mismatch-induced stress. We introduced an analytical approach to model diffusion and ion types of doping with account concurrent changing of parameters in space and time. At the same time the approach gives us possibility to take into account nonlinearity of considered processes.