Analysis of Erbium and Vanadium Diffusion in Porous Silicon Carbide

Abstract

Experimental data on diffusion of erbium and vanadium in porous and nonporous silicon carbide at 1700 and 2200°C have been used for modelling diffusion in porous SiC. It is shown that the consideration of pore structure modification under annealing via vacancy redistribution allows for satisfactory description of dopant diffusion. As expected, important contribution to the diffusion in the porous medium is found to be made by the walls of the pores: in SiC, the vacancy surface diffusion coefficient on the walls appears to exceed that in the bulk of the material by an order of magnitude. When thermal treatment transforms pore channels into closed voids, pathways for accelerated diffusion cease to exist and diffusion rates in porous and nonporous SiC become similar.

1. Introduction

Silicon carbide (SiC) is a semiconductor that possesses properties highly suitable for high-power, high-frequency, and high-temperature microelectronic applications. For a long time, selective doping in SiC technology has been implemented mostly with ion implantation. Diffusion doping was considered to be impractical due to low diffusion coefficients of impurities at temperatures below 1800°C, and the interest in diffusion in SiC was generally related to applications requiring “drive-in” annealing after implantation (see, e.g., [1–3]). Papers devoted to “classic” thermal diffusion in SiC were mostly limited to early years [4, 5]. However, recent achievements in fabrication of SiC-based microelectronic devices using diffusion doping [1, 2, 6] stirred renewed interest in diffusion in SiC.

Diffusion in SiC is believed to proceed mostly via kick-out or vacancy-related mechanisms [1–5]. Vacancies define many physical properties of SiC [7, 8], and one of the current subjects of research is their contribution to the formation and properties of porous SiC, a material with many prospective applications [9, 10]. The very formation of porous structure in SiC is believed to proceed via vacancy migration and clustering [11], and there are all grounds to believe that experimentally observed modification of the structure under thermal annealing [12, 13] is also due to vacancy redistribution. In this respect, diffusion in porous SiC represents a subject of special interest, as pores serve as the sources of vacancies, which should contribute to diffusion. In this paper, we consider diffusion in SiC with porous structure formed by the vacancy mechanism, and model a dynamic process, where diffusion of dopants proceeds against the background of fluxes of vacancies responsible for pore ripening during high-temperature annealing.

2. Experimental Details

To form the porous structure, 6H–SiC wafers were subjected to surface anodization in HF. The anodization mode was chosen such as to ensure formation of pores via vacancy-related mechanism, with no chemical etching of the surface involved [11]. As a result of anodization, porous layers were formed with thickness significantly greater than the expected diffusion depth. The dopants were erbium and vanadium. Erbium doping of SiC is attractive for infrared applications relying on radiative intrashell transitions in Er ions placed in semiconductor lattice [14], while vanadium is a dopant of choice for fabricating semi-insulating SiC [15]. In each experiment, the dopant was simultaneously introduced into SiC layer with a “nanoporous” structure, and into a reference sample of nonporous SiC. The diffusion was carried out in the so-called sandwich-type cell. Conditions in the cell prevent the crystal surface from decomposition in the course of high-temperature annealing, thus making possible quasi-equilibrium diffusion doping of nonporous SiC at temperatures 𝑇≤2600°C [5]. The diffusion was carried out at 𝑇=1700 and 2200°C for a period of time of 30 and 120 min. The dopant profiles were investigated using secondary ion mass spectroscopy (SIMS).

3. Experimental Results

To follow changes in the porous structure induced by the high temperature, a control set of samples was first subjected to thermal anneal under the same conditions, which were later used for diffusion. Figure 1 shows a typical example of the modification of porous structure in 6H-SiC as a result of anneal at 𝑇=1700°C for 120 and 240 min. As can be seen, initial porous structure, which appeared in a form of branched thin pores with the pore stem diameter of 30 to 40 nm (a), after 120 min of anneal evolved into a set of “cylinders” with the diameter of 100 to 200 nm (b). Further annealing transformed the “cylinders” into closed voids, with some voids showing tendency to faceting (c). The pattern of porous structure after anneal at 2200°C for 120 min was similar to that shown in Figure 1(c), as with temperature increasing the modification of porous structure proceeds faster [13]. According to Cheremskoy et al. [16], the closing and faceting of the voids confirm that pore modification proceeds according to ripening mechanism, and it is due to vacancy redistribution.

(a)

(b)

(c)

Figure 1 

Transmission (a) and scanning ((b) and (c)) electron microscopy images of the cleft edge of porous layer in SiC: straight after anodization (a), after 1700°C anneal for 120 min (b), and after 1700°C anneal for 240 min (c).

SIMS measurements with depth profiling, performed both after the above anneals and dopant diffusion, did not reveal any serious deviation from SiC stoichiometry in the processed samples. The examples of dopant diffusion profiles are presented in Figure 2. As was shown elsewhere [17], only a small part of such profiles related to the area close to the surface could be fitted with complementary error function, which testified to the concentration-dependent diffusion. Figure 2(a) presents erbium profiles in porous and nonporous SiC after diffusion at 2200°C for 30 and 120 min. As is seen, the erbium surface concentration in the samples was close to 3 × 10¹⁸ cm³. For samples diffused for 30 min, we observed a difference in the diffusion depth in porous and nonporous specimens. For samples diffused for 120 min, this difference became somewhat larger, but considering the difference in the surface concentration, not as large as could be expected if the “classic” inverse square root diffusion law was applicable.

(a)

(b)

Figure 2 

Experimental profiles showing erbium (a) and vanadium (b) concentration versus depth after diffusion at 2200°C for 30 min (curves 1 and 2) and 120 min (curves 1′ and 2′), in: 1: nonporous SiC, 2: porous SiC.

A similar situation was observed with vanadium diffusion at 2200°C. For all the diffused samples, the vanadium surface concentration was ~2 × 10¹⁷ cm⁻³ (Figure 2(b)). The short-time (30 min) diffusion yielded diffusion depth in porous SiC larger than that in nonporous material. Increasing the anneal time to 120 min, the difference in diffusion depth in porous and nonporous specimens remained almost the same.

4. Modelling

For modelling, let us consider diffusion of impurity atoms into a porous layer. The diffusion starts at the moment of time 𝑡=0 and proceeds through an interface with 𝑥=0. It can be described with the second Fick law written in the following form: 𝜕𝐶𝑥,𝑦,𝑧,𝑡𝐷𝜕𝑡=div_𝐶𝐶𝐷⋅grad𝑥,𝑦,𝑧,𝑡+div^𝐶𝑆𝜇,𝑉𝑘𝑇⋅grad𝑥,𝑦,𝑧,𝑡 (1) with the following boundary and initial conditions: 𝐶0,𝑦,𝑧,0=𝐶₀,𝐶𝑥>0,𝑦,𝑧,0=0,𝜕𝐶𝑥,𝑦,𝑧,𝑡𝜕𝑥_{𝑥=𝐿𝑥}=𝜕𝐶𝑥,𝑦,𝑧,𝑡𝜕𝑦_𝑦=0=𝜕𝐶𝑥,𝑦,𝑧,𝑡𝜕𝑦_{𝑦=𝐿𝑦}=𝜕𝐶𝑥,𝑦,𝑧,𝑡𝜕𝑧_𝑧=0=𝜕𝐶𝑥,𝑦,𝑧,𝑡𝜕𝑧_{𝑧=𝐿𝑧}𝐶=0,0,𝑦,𝑧,𝑡=𝐶₀, (2) where 𝑘=1.38⋅10⁻²³JK⁻¹ is the Boltzmann constant, 𝑉 is molar volume; 𝜇(𝑥,𝑦,𝑧,𝑡)=𝑅⋅𝑇⋅ln(𝑉₂/𝑉₁) is the chemical potential, 𝑉₁ and 𝑉₂ are the initial and final pore volumes, respectively, 𝑅=8.31Jmole⁻¹K⁻¹ is the molar gas constant, 𝐿_𝑥, 𝐿_𝑦, and 𝐿_𝑧 are the dimensions of doped sample, 𝐶(𝑥,𝑦,𝑧,𝑡) is the impurity concentration, and 𝐷_𝐶 and 𝐷_𝐶𝑆 are the diffusion coefficients in the bulk and on the surface (pore walls) of the material, respectively (see, e.g., [18]). As follows from the experimental data, we should consider the dependence of 𝐷_𝐶 and 𝐷_𝐶𝑆 on the vacancy and dopant concentration. These dependences can be described using the functions suggested by Gotra [19], Zorin et al. [20], and Ryssel and Ruge [21], respectively, 𝐷_𝐶𝑉=𝐷_𝐶𝑉𝐿(𝑇)1+𝜉_𝑉𝐶^𝜂𝑥,𝑦,𝑧,𝑡𝑃^𝜂×𝑥,𝑦,𝑧,𝑇1+𝜍_𝑉𝑉𝑥,𝑦,𝑧,𝑡𝑉_∞,𝐷_𝐶𝑆=𝐷_𝐶𝑆𝐿(𝑇)1+𝜉_𝑆𝐶^𝜂𝑥,𝑦,𝑧,𝑡𝑃^𝜂×𝑥,𝑦,𝑧,𝑇1+𝜍_𝑆𝑉𝑥,𝑦,𝑧,𝑡𝑉_∞. (3) Here, 𝐷_𝐶𝐿(𝑇) and 𝐷_𝐶𝑆𝐿(𝑇) (where index 𝑉 relates to diffusion in bulk and 𝑆 relates to surface diffusion) are linear (i.e., independent on dopant concentration, the linearity being indicated by index 𝐿) coefficients of the approximations. These (as well as solubility limits) have been extracted from the experimental data presented earlier [17]. The second multipliers in (3) describe concentration dependence of diffusion coefficients [19]. 𝑃(𝑥,𝑦,𝑧,𝑇) is the dopant solubility limit. Value of the parameter 𝜂 (which is equal to an integer value in the range of 1 to 3 [19]) depends on the interaction between dopant atoms and defects [19, 22]. The third multipliers in (3) describe vacancy dependences of the diffusion coefficients [20, 21]. 𝑉_∞ is the equilibrium vacancy distribution. Four fitting parameters, 𝜉_𝑉, 𝜉_𝑆, 𝜁_𝑉, 𝜁_𝑆 can be considered as scale multipliers of the limit of solubility of dopant and equilibrium vacancy concentration. The vacancy concentration can be described as follows: 𝜕𝑉𝑥,𝑦,𝑧,𝑡𝐷𝜕𝑡=div_𝑉(𝑇)𝑉𝐷⋅grad𝑥,𝑦,𝑧,𝑡+div^𝑉𝑆(𝑇)𝜇𝑉𝑘𝑇⋅grad𝑥,𝑦,𝑧,𝑡 (4) with the initial and boundary conditions: 𝑉𝑥,𝑦,𝑧,𝑡_𝑆1+^→_𝑉𝑛𝑡=𝑉_∞1+2𝛾𝜔𝑘𝑇𝑥²₁+𝑦²₁+𝑧²₁,𝑉𝑥,𝑦,𝑧,𝑡_𝑆2=𝑉_∞,𝑉𝑥,𝑦,𝑧,0=𝑓_𝑉𝑥,𝑦,𝑧=𝑉_∞2𝛾𝜔𝑘𝑇𝑥²₁+𝑦²₁+𝑧²₁×𝑒^{−𝑥/𝑥1}𝑒^{−𝑦/𝑦1}𝑒^{−𝑧/𝑧1}−1, (5) where 𝑆₁ designates initial pore surface, 𝑆₂ stands for the surface of the outward boundary of the porous area, 𝑉_𝑛 is the absolute value of the velocity of the movement of the surface of the growing or shrinking pore [18], 𝐷_𝑉 and 𝐷_𝑉𝑆 are volumetric and surface diffusion coefficients of vacancies; 𝛾 is the specific surface energy [16], and 𝜔=𝑎³, where 𝑎 is interatomic distance. The appropriateness of using (4) for describing the evolution of a pore under annealing was considered in detail by Kitayama et al. [18]. The initial and boundary conditions for this equation represent three-dimensional generalization of those given by Cheremskoy et al. [16], with the consideration of pore size changes.

A typical vacancy distribution in the vicinity of a pore, which was obtained as a result of modelling using (4) is shown in Figure 3. Here, “x” is the coordinate across the pore: the maximum value of vacancy concentration corresponds to the center of the pore. The initial pore was regarded as a cylinder with the radius of 40 nm. The parameters of modelling were chosen so that the simulated annealing would lead to the growth of pore “cylinder” with the diameter of ~200 nm, which corresponds to the experimental data presented in Figure 1(b). The figure shows how annealing leads to the increase in pore diameter, and confirms qualitative adequacy of our model for consideration of pore modification under heat treatment on the microscopic scale.

Figure 3 

Modelled vacancy distribution in the vicinity of a pore in porous SiC before (curve 1) and after (curve 2) anneal at 1700°C for 120 min.

To consider the diffusion, let us transform the diffusion coefficients into the following form: 𝐷_𝐿(𝑥,𝑦,𝑧,𝑇)=𝐷_0𝐿[1+𝜀_𝐿𝑔_𝐿(𝑥,𝑦,𝑧,𝑇)] and 𝐷_𝑉(𝑥,𝑦,𝑧,𝑇)=𝐷_0𝑉[1+𝜀_𝑉𝑔_𝑉(𝑇)], where 𝐷_0𝐿 and 𝐷_0𝑉 are average values of diffusion coefficients, 0≤𝜀_𝐿<1, 0≤𝜀_𝑉<1, |𝑔_𝐿(𝑥,𝑦,𝑧,𝑇)|≤1, |𝑔_𝑉(𝑇)|≤1, (see, e.g., [22–24]). We determined solutions of (1) and (4) as the following power series: 𝐶𝑥,𝑦,𝑧,𝑡=^∞_𝑖=0𝜀^𝑖_𝐿^∞_𝑗=0𝜉^𝑗𝐶_𝑖𝑗,𝑉=𝑥,𝑦,𝑧,𝑡𝑥,𝑦,𝑧,𝑡^∞_𝑖=0𝜀^𝑖_𝑉𝑉_𝑖𝑥,𝑦,𝑧,𝑡. (6) Substitution of the series (6) into (1), (4), and appropriate boundary and initial conditions gives us possibility to obtain equations for zero-order approximation of dopant and vacancy concentration 𝐶₀₀(𝑥,𝑦,𝑧,𝑡) and 𝑉₀(𝑥,𝑦,𝑧,𝑡), their corrections 𝐶_𝑖𝑗(𝑥,𝑦,𝑧,𝑡) and 𝑉_𝑖(𝑥,𝑦,𝑧,𝑡)(𝑖≥1,𝑗≥1), and the boundary and initial conditions for all functions 𝐶_𝑖𝑗(𝑥,𝑦,𝑧,𝑡) and 𝑉_𝑖(𝑥,𝑦,𝑧,𝑡)(𝑖≥0,𝑗≥0): 𝜕𝐶₀₀𝑥,𝑦,𝑧,𝑡𝜕𝑡=𝐷_0𝐶Δ𝐶₀₀𝐷𝑥,𝑦,𝑧,𝑡+div^𝐶𝑆𝜇,𝑉𝑘𝑇⋅grad𝑥,𝑦,𝑧,𝑡𝜕𝐶_𝑖0𝑥,𝑦,𝑧,𝑡𝜕𝑡=𝐷_0𝐶Δ𝐶_𝑖0𝑥,𝑦,𝑧,𝑡+𝐷_0𝐶div𝑔_𝐿𝐶𝑥,𝑦,𝑧,𝑇×grad₀₀𝑥,𝑦,𝑧,𝑡,𝑖≥1,𝜕𝐶₀₁𝑥,𝑦,𝑧,𝑡𝜕𝑡=𝐷_0𝐶Δ𝐶₀₁𝑥,𝑦,𝑧,𝑡+𝐷_0𝐶𝐶div^𝛾00𝑥,𝑦,𝑧,𝑡𝑃^𝛾𝐶𝑥,𝑦,𝑧,𝑇grad₀₀𝑥,𝑦,𝑧,𝑡,𝜕𝐶₀₂𝑥,𝑦,𝑧,𝑡𝜕𝑡=𝐷_0𝐶div𝐶₀₁𝐶𝑥,𝑦,𝑧,𝑡^𝛾−100𝑥,𝑦,𝑧,𝑡𝑃^𝛾𝐶𝑥,𝑦,𝑧,𝑇×grad₀₀𝑥,𝑦,𝑧,𝑡+𝐷_0𝐶Δ𝐶₀₂𝑥,𝑦,𝑧,𝑡+𝐷_0𝐶𝐶div^𝛾00𝑥,𝑦,𝑧,𝑡𝑃^𝛾𝐶𝑥,𝑦,𝑧,𝑇×grad₀₁𝑥,𝑦,𝑧,𝑡,𝜕𝐶₁₁𝑥,𝑦,𝑧,𝑡𝜕𝑡=𝐷_0𝐶Δ𝐶₁₁𝑥,𝑦,𝑧,𝑡+𝐷_0𝐶div𝑔_𝐿𝐶𝑥,𝑦,𝑧,𝑇×grad₀₁𝑥,𝑦,𝑧,𝑡+𝐷_0𝐶𝐶div^𝛾00𝑥,𝑦,𝑧,𝑡𝑃^𝛾𝐶𝑥,𝑦,𝑧,𝑇×grad₁₀𝑥,𝑦,𝑧,𝑡,𝜕𝑉₀𝑥,𝑦,𝑧,𝑡𝜕𝑡=𝐷_0𝑉Δ𝑉₀𝐷𝑥,𝑦,𝑧,𝑡+div^𝑉𝑆𝜇,𝑉𝑘𝑇⋅grad𝑥,𝑦,𝑧,𝑡𝜕𝑉_𝑖𝑥,𝑦,𝑧,𝑡𝜕𝑡=𝐷_0𝑉Δ𝑉_𝑖0𝑥,𝑦,𝑧,𝑡+𝐷_0𝑉div𝑔_𝑉𝑉𝑥,𝑦,𝑧,𝑇×grad_𝑖−1𝑥,𝑦,𝑧,𝑡,𝑖≥1,𝜕𝐶_𝑖𝑗𝑥,𝑦,𝑧,𝑡𝜕𝑥_𝑥=0=𝜕𝐶_𝑖𝑗𝑥,𝑦,𝑧,𝑡𝜕𝑥_{𝑥=𝐿𝑥}=𝜕𝐶_𝑖𝑗𝑥,𝑦,𝑧,𝑡𝜕𝑦_𝑦=0=𝜕𝐶_𝑖𝑗𝑥,𝑦,𝑧,𝑡𝜕𝑦_{𝑦=𝐿𝑦}=𝜕𝐶_𝑖𝑗𝑥,𝑦,𝑧,𝑡𝜕𝑧_𝑧=0=𝜕𝐶_𝑖𝑗𝑥,𝑦,𝑧,𝑡𝜕𝑧_{𝑧=𝐿𝑧}𝐶=0,₀₀0,𝑦,𝑧,0=𝐶₀,𝐶_𝑖𝑗,𝐶0,𝑦,𝑧,0=0𝑖≥1,𝑗≥1_𝑖𝑗,𝐶𝑥>0,𝑦,𝑧,0=0𝑖≥0,𝑗≥0₀₀0,𝑦,𝑧,𝑡=𝐶₀,𝐶_𝑖𝑗,𝑉0,𝑦,𝑧,𝑡=0𝑖≥1,𝑗≥1₀𝑥,𝑦,𝑧,𝑡_𝑆2=𝑉_∞,𝑉_𝑖𝑥,𝑦,𝑧,𝑡_𝑆2(),𝑉=0𝑖≥0₀₀𝑥,𝑦,𝑧,𝑡_𝑆1+^→_𝑉𝑛𝑡=𝑉_∞1+2𝛾𝜔𝑘𝑇𝑥²₁+𝑦²₁+𝑧²₁,𝜕𝑉_𝑖𝑥,𝑦,𝑧,𝑡𝜕𝑥_𝑥=0=𝜕𝑉_𝑖𝑥,𝑦,𝑧,𝑡𝜕𝑥_{𝑥=𝐿𝑥}=𝜕𝑉_𝑖𝑥,𝑦,𝑧,𝑡𝜕𝑦_𝑦=0=𝜕𝑉_𝑖𝑥,𝑦,𝑧,𝑡𝜕𝑦_{𝑦=𝐿𝑦}=𝜕𝑉_𝑖𝑥,𝑦,𝑧,𝑡𝜕𝑧_𝑧=0=𝜕𝑉_𝑖𝑥,𝑦,𝑧,𝑡𝜕𝑧_{𝑧=𝐿𝑧}𝑉=0,_𝑖𝑗𝑥,𝑦,𝑧,𝑡_𝑆1+^→_𝑉𝑛𝑡,𝑉=0𝑖≥0,𝑗≥0_𝑖𝑗,𝑉𝑥,𝑦,𝑧,0=0𝑖≥0,𝑗≥0₀₀𝑥,𝑦,𝑧,0=𝑉_∞2𝛾𝜔𝑘𝑇𝑥²₁+𝑦²₁+𝑧²₁×𝑒^{−𝑥/𝑥1}𝑒^{−𝑦/𝑦1}𝑒^{−𝑧/𝑧1}−1. (7)(8)(9) Solutions of (7) with conditions (8) are presented in the appendix. Though the values of the vacancy and dopant concentration can be determined with more accuracy, second-order approximation usually gives satisfactory qualitative and some quantitative results (see, e.g., [22–24]).

5. Discussion

Let us analyze the experimental diffusion profiles using the relations obtained in the previous section. We will restrict our consideration to one dimension only, which is inward the porous layer (direction x). Considering the experimental results, we may assume that the results of 30 min-long anneals were somehow affected by transient processes, yet 120 min-long diffusion proceeded under “steady-state” conditions. This assumption is supported by the fact that the propagation of the diffusion fronts as presented in Figure 2 did not obey the inverse square root law. We fitted experimental profiles shown in Figure 2 by minimization of mean-squared error between experimental and theoretical results. The solid lines correspond to erbium profiles in porous and nonporous SiC after diffusion at 1700°C according to the procedure developed above for 120 min, that is, long diffusion anneals. The parameters obtained are listed in Table 1 and the fitting curves are shown in Figure 4. Values of fitting parameters are presented in the Table 2. The parameters 𝜉_𝑉, 𝜉_𝑆, 𝜁_𝑉, and 𝜁_𝑆 give us corrections to experimental values of dopant solubility limit and equilibrium vacancy distribution. Their usefulness is determined by the fact that one can obtain quite a large difference between experimental data taken from different references.

Table 1 

Diffusion parameters obtained from the fitting.

Table 2 

Fitting parameters for Figure 4.

(a)

(b)

(c)

Figure 4 

Distribution in non-porous SiC (curves 1) and porous SiC (curves 2) of: erbium after diffusion at 1700°C for 120 min (a), vanadium after diffusion at 2200°C for 120 min (b), and erbium after diffusion at 2200°C for 120 min (c). Points show experimental data, solid curves show results of fitting.

It follows from the results obtained that the vacancy diffusion coefficient on the surface of the pore walls 𝐷_𝑉𝑆 exceeds that in the bulk of the material 𝐷_𝑉 approximately by an order of magnitude. Erbium appears to diffuse in SiC much faster than vanadium, seemingly due to stronger concentration dependence of its diffusion coefficient, that is, concentration dependence of diffusion coefficient leads to higher diffusion coefficient for erbium in comparison with that of vanadium. Summarizing the whole set of diffusion data obtained (including our earlier data [17]), it appears that when a short-term diffusion anneal is used, diffusion in porous SiC proceeds faster than that in nonporous SiC simply because porous structure does not have enough time to evolve into closed “cylinders.” For long-term anneals, pathways for accelerated diffusion (pore walls) disappear as pore channels get transformed into closed voids, and diffusion coefficients in porous and nonporous SiC become close. This effect becomes more pronounced with increasing diffusion temperature. At low temperatures (𝑇<1000°C), the diffusion in porous SiC always proceeds much faster than that in nonporous SiC due to weak pore structure modification.

6. Conclusion

Modelling based on experimental data on the diffusion of erbium and vanadium in porous silicon carbide at 1700 and 2200°C shows that this process is greatly affected by the diffusion of vacancies, which is a typical process in the material with vacancy-related porosity under thermal treatment. Consideration of the vacancy redistribution under annealing allows for satisfactory describing the dopant diffusion in porous SiC and yields basic diffusion parameters. The vacancy diffusion coefficient on the surface of the pore walls appears to exceed that in the bulk of the material by an order of magnitude. The values of “bulk” diffusion coefficients obtained in this work agree with the data obtained for typical dopants in SiC (see, e.g., [25]). Generally, it appears that the consideration of vacancy redistribution during high-temperature annealing allows for correct description of diffusion in a porous semiconductor.

Appendix

Equations for the functions 𝐶_𝑖𝑗(𝑥,𝑦,𝑧,𝑡) and 𝑉_𝑖(𝑥,𝑦,𝑧,𝑡)(𝑖≥0,𝑗≥0) have been solved by the standard Fourier approach (see, e.g., the work by Carslaw and Jaeger [26]). Solutions of (7) with conditions (8) can be written in the form: 𝐶₀₀𝑥,𝑦,𝑧,𝑡=−2𝜋𝐷_0𝐶𝐿²_𝑥𝐿_𝑦𝐿_𝑧^∞_𝑛=1𝑛𝐹_𝑛𝐶𝑐_𝑛(𝑥)𝑐_𝑛𝑦𝑐_𝑛(𝑧)𝑒_𝑛𝐶(𝑡)×^𝑡₀𝑒_𝑛𝐶(−𝜏)^𝐿𝑥₀𝑠_𝑛(𝑢)^𝐿𝑦₀𝑐_𝑛(𝑣)×^𝐿𝑧₀𝑐_𝑛(𝑤)𝐷_𝐶𝑆(𝑢,𝑣,𝑤,𝑇)()𝑉𝑘𝑇𝜕𝜇𝑢,𝑣,𝑤,𝜏−𝜕𝑢𝑑𝑤𝑑𝑣𝑑𝑢𝑑𝜏2𝜋𝐷_0𝐶𝐿_𝑥𝐿²_𝑦𝐿_𝑧^∞_𝑛=1𝐹_𝑛𝐶𝑐_𝑛(𝑥)𝑐_𝑛𝑦𝑐_𝑛(𝑧)×^𝑡₀𝑒_𝑛𝐶(−𝜏)^𝐿𝑥₀𝑐_𝑛(𝑢)^𝐿𝑦₀𝑠_𝑛(𝑣)×^𝐿𝑧₀𝑐_𝑛(𝑤)𝐷_𝐶𝑆(𝑢,𝑣,𝑤,𝑇)()𝑉𝑘𝑇𝜕𝜇𝑢,𝑣,𝑤,𝜏𝜕𝑣𝑑𝑤𝑑𝑣𝑑𝑢𝑑𝜏×𝑛𝑒_𝑛𝐶(𝑡)−2𝜋𝐷_0𝐶𝐿_𝑥𝐿_𝑦𝐿²_𝑧^∞_𝑛=1𝑛𝐹_𝑛𝐶𝑐_𝑛(𝑥)×𝑐_𝑛𝑦𝑐_𝑛(𝑧)𝑒_𝑛𝐶(𝑡)^𝑡₀𝑒_𝑛𝐶(−𝜏)^𝐿𝑥₀𝑐_𝑛(𝑢)^𝐿𝑦₀𝑐_𝑛(𝑣)×^𝐿𝑧₀𝑠_𝑛(𝑤)𝐷_𝐶𝑆(𝑢,𝑣,𝑤,𝑇)()𝑉𝑘𝑇𝜕𝜇𝑢,𝑣,𝑤,𝜏𝜕𝑤𝑑𝑤𝑑𝑣𝑑𝑢𝑑𝜏, (A.1) where 𝑒_𝑛𝐶(𝑡)=exp[−𝜋²𝑛²𝐷_0𝐿𝑡(𝐿_𝑥⁻²+𝐿_𝑦⁻²+𝐿_𝑧⁻²)],𝐹_𝑛𝐶=∫^𝐿𝑥₀𝑐_𝑛(𝑢)∫^𝐿𝑦₀𝑐_𝑛(𝑣)∫^𝐿𝑧₀𝑓_𝐶(𝑢,𝑣,𝑤)𝑐_𝑛(𝑤)𝑑𝑤𝑑𝑣𝑑𝑢,𝑐_𝑛(𝑥)=cos(𝜋𝑛𝑥/𝐿), 𝐶_𝑖0𝑥,𝑦,𝑧,𝑡=−2𝜋𝐷_0𝐶𝐿²_𝑥𝐿_𝑦𝐿_𝑧^∞_𝑛=1𝑛𝐹_𝑛𝐶𝑐_𝑛(𝑥)𝑐_𝑛𝑦𝑐_𝑛(𝑧)𝑒_𝑛𝐶(𝑡)×^𝑡₀𝑒_𝑛𝐶(−𝜏)^𝐿𝑥₀𝑠_𝑛(𝑢)^𝐿𝑦₀𝑐_𝑛(𝑣)×^𝐿𝑧₀𝑐_𝑛(𝑤)𝑔_𝐿()×𝑢,𝑣,𝑤,𝑇𝜕𝐶_𝑖−10(𝑢,𝑣,𝑤,𝜏)−𝜕𝑢𝑑𝑤𝑑𝑣𝑑𝑢𝑑𝜏2𝜋𝐷_0𝐶𝐿_𝑥𝐿²_𝑦𝐿_𝑧^∞_𝑛=1𝐹_𝑛𝐶𝑐_𝑛(𝑥)𝑐_𝑛𝑦×𝑐_𝑛(𝑧)^𝑡₀𝑒_𝑛𝐶(−𝜏)^𝐿𝑥₀𝑐_𝑛(𝑢)^𝐿𝑦₀𝑠_𝑛(𝑣)×^𝐿𝑧₀𝑐_𝑛(𝑤)𝑔_𝐿()×𝑢,𝑣,𝑤,𝑇𝜕𝐶_𝑖−10(𝑢,𝑣,𝑤,𝜏)𝜕𝑣𝑑𝑤𝑑𝑣𝑑𝑢𝑑𝜏×𝑛𝑒_𝑛𝐶(𝑡)−2𝜋𝐷_0𝐶𝐿_𝑥𝐿_𝑦𝐿²_𝑧^∞_𝑛=1𝑛𝐹_𝑛𝐶𝑐_𝑛(𝑥)𝑐_𝑛𝑦𝑐_𝑛(𝑧)×𝑒_𝑛𝐶(𝑡)^𝑡₀𝑒_𝑛𝐶(−𝜏)^𝐿𝑥₀𝑐_𝑛(𝑢)^𝐿𝑦₀𝑐_𝑛(𝑣)^𝐿𝑧₀𝑠_𝑛(𝑤)×𝑔_𝐿(𝑢,𝑣,𝑤,𝑇)𝜕𝐶_𝑖−10(𝑢,𝑣,𝑤,𝜏)𝜕𝑤𝑑𝑤𝑑𝑣𝑑𝑢𝑑𝜏, (A.2) where 𝑖≥1,𝑠_𝑛(𝑥)=sin(𝜋𝑛𝑥/𝐿), 𝐶₀₁𝑥,𝑦,𝑧,𝑡=−2𝜋𝐷_0𝐶𝐿²_𝑥𝐿_𝑦𝐿_𝑧^∞_𝑛=1𝑛𝐹_𝑛𝐶𝑐_𝑛(𝑥)𝑐_𝑛𝑦𝑐_𝑛(𝑧)𝑒_𝑛𝐶(𝑡)×^𝑡₀𝑒_𝑛𝐶(−𝜏)^𝐿𝑥₀𝑠_𝑛(𝑢)^𝐿𝑦₀𝑐_𝑛(𝑣)^𝐿𝑧₀𝑐_𝑛(𝑤)×𝐶^𝛾₀₀(𝑢,𝑣,𝑤,𝜏)𝑃^𝛾(𝑢,𝑣,𝑤,𝑇)𝜕𝐶₀₀(𝑢,𝑣,𝑤,𝜏)−𝜕𝑣𝑑𝑤𝑑𝑣𝑑𝑢𝑑𝜏2𝜋𝐷_0𝐶𝐿_𝑥𝐿²_𝑦𝐿_𝑧^∞_𝑛=1𝐹_𝑛𝐶𝑐_𝑛(𝑥)𝑐_𝑛𝑦×𝑛^𝑡₀𝑒_𝑛𝐶(−𝜏)^𝐿𝑥₀𝑐_𝑛(𝑢)^𝐿𝑦₀𝑠_𝑛(𝑢)^𝐿𝑧₀𝑐_𝑛(𝑤)×𝐶^𝛾₀₀(𝑢,𝑣,𝑤,𝜏)𝑃^𝛾(𝑢,𝑣,𝑤,𝑇)𝜕𝐶₀₀(𝑢,𝑣,𝑤,𝜏)𝜕𝑣𝑑𝑤𝑑𝑣𝑑𝑢𝑑𝜏×𝑐_𝑛(𝑧)𝑒_𝑛𝐶(𝑡)−2𝜋𝐷_0𝐶𝐿_𝑥𝐿_𝑦𝐿²_𝑧^∞_𝑛=1𝑛𝑐_𝑛(𝑥)𝑐_𝑛𝑦𝑐_𝑛(𝑧)𝑒_𝑛𝐶(𝑡)×^𝑡₀𝑒_𝑛𝐶(−𝜏)^𝐿𝑥₀𝑐_𝑛(𝑢)^𝐿𝑦₀𝑐_𝑛(𝑢)^𝐿𝑧₀𝑠_𝑛(𝑤)×𝐹_𝑛𝐶𝐶^𝛾₀₀(𝑢,𝑣,𝑤,𝜏)𝑃^𝛾(𝑢,𝑣,𝑤,𝑇)𝜕𝐶₀₀(𝑢,𝑣,𝑤,𝜏)𝐶𝜕𝑣𝑑𝑤𝑑𝑣𝑑𝑢𝑑𝜏,₀₂𝑥,𝑦,𝑧,𝑡=−2𝜋𝐷_0𝐶𝐿²_𝑥𝐿_𝑦𝐿_𝑧^∞_𝑛=1𝑛𝐹_𝑛𝐶𝑐_𝑛(𝑥)𝑐_𝑛𝑦𝑐_𝑛(𝑧)𝑒_𝑛𝐶(𝑡)×^𝑡₀𝑒_𝑛𝐶(−𝜏)^𝐿𝑥₀𝑠_𝑛(𝑢)^𝐿𝑦₀𝑐_𝑛(𝑣)×^𝐿𝑧₀𝐶^𝛾−1₀₀(𝑢,𝑣,𝑤,𝜏)𝑃^𝛾()×𝐶𝑢,𝑣,𝑤,𝑇₀₁(𝑢,𝑣,𝑤,𝜏)𝜕𝐶₀₀(𝑢,𝑣,𝑤,𝜏)𝜕𝑢+𝐶₀₀(𝑢,𝑣,𝑤,𝜏)𝜕𝐶₀₁(𝑢,𝑣,𝑤,𝜏)𝜕𝑢×𝑐_𝑛(𝑤)𝑑𝑤𝑑𝑣𝑑𝑢𝑑𝜏−2𝜋𝐷_0𝐶𝐿_𝑥𝐿²_𝑦𝐿_𝑧^∞_𝑛=1𝑛𝐹_𝑛𝐶𝑐_𝑛(𝑥)×𝑐_𝑛𝑦𝑐_𝑛(𝑧)𝑒_𝑛𝐶(𝑡)^𝑡₀𝑒_𝑛𝐶(−𝜏)^𝐿𝑥₀𝑐_𝑛(𝑢)×^𝐿𝑦₀𝑠_𝑛(𝑣)^𝐿𝑧₀𝑐_𝑛(𝑤)𝐶^𝛾−1₀₀(𝑢,𝑣,𝑤,𝜏)𝑃^𝛾()×𝐶𝑢,𝑣,𝑤,𝑇₀₁(𝑢,𝑣,𝑤,𝜏)𝜕𝐶₀₀(𝑢,𝑣,𝑤,𝜏)𝜕𝑣+𝐶₀₀(𝑢,𝑣,𝑤,𝜏)𝜕𝐶₀₁(𝑢,𝑣,𝑤,𝜏)−𝜕𝑣𝑑𝑤𝑑𝑣𝑑𝑢𝑑𝜏2𝜋𝐷_0𝐶𝐿_𝑥𝐿_𝑦𝐿²_𝑧^∞_𝑛=1𝑛𝐹_𝑛𝐶𝑐_𝑛(𝑥)𝑐_𝑛𝑦×𝑐_𝑛(𝑧)𝑒_𝑛𝐶(𝑡)^𝑡₀𝑒_𝑛𝐶(−𝜏)^𝐿𝑥₀𝑐_𝑛(𝑢)×^𝐿𝑦₀^𝐿𝑧₀𝐶^𝛾−1₀₀(𝑢,𝑣,𝑤,𝜏)𝑃^𝛾()×𝐶𝑢,𝑣,𝑤,𝑇₀₁(𝑢,𝑣,𝑤,𝜏)𝜕𝐶₀₀(𝑢,𝑣,𝑤,𝜏)𝜕𝑤+𝐶₀₀(𝑢,𝑣,𝑤,𝜏)𝜕𝐶₀₁(𝑢,𝑣,𝑤,𝜏)𝜕𝑤×𝑠_𝑛(𝑤)𝑑𝑤𝑐_𝑛(𝑣)𝐶𝑑𝑣𝑑𝑢𝑑𝜏,₁₁𝑥,𝑦,𝑧,𝑡=−2𝜋𝐷_0𝐶𝐿²_𝑥𝐿_𝑦𝐿_𝑧^∞_𝑛=1𝑛𝐹_𝑛𝐶𝑐_𝑛(𝑥)𝑐_𝑛𝑦𝑐_𝑛(𝑧)𝑒_𝑛𝐶(𝑡)×^𝑡₀𝑒_𝑛𝐶(−𝜏)^𝐿𝑥₀𝑠_𝑛(𝑢)^𝐿𝑦₀𝑐_𝑛(𝑣)×^𝐿𝑧₀𝑔_𝐿(𝑢,𝑣,𝑤,𝑇)𝜕𝐶₀₁(𝑢,𝑣,𝑤,𝜏)+𝐶𝜕𝑢_𝛾00(𝑢,𝑣,𝑤,𝜏)𝑃^𝛾(𝑢,𝑣,𝑤,𝑇)𝜕𝐶₁₀(𝑢,𝑣,𝑤,𝜏)𝜕𝑢+𝐶₁₀()𝐶𝑢,𝑣,𝑤,𝜏^𝛾−100(𝑢,𝑣,𝑤,𝜏)𝑃^𝛾()×𝑢,𝑣,𝑤,𝑇𝜕𝐶₀₀(𝑢,𝑣,𝑤,𝜏)𝜕𝑢+𝑔_𝐿()×𝐶𝑢,𝑣,𝑤,𝑇_𝛾00(𝑢,𝑣,𝑤,𝜏)𝑃^𝛾(𝑢,𝑣,𝑤,𝑇)𝜕𝐶₀₀(𝑢,𝑣,𝑤,𝜏)𝜕𝑢×𝑐_𝑛(𝑤)𝑑𝑤𝑑𝑣𝑑𝑢𝑑𝜏−2𝜋𝐷_0𝐶𝐿_𝑥𝐿²_𝑦𝐿_𝑧^∞_𝑛=1𝑛𝐹_𝑛𝐶𝑐_𝑛(𝑥)𝑐_𝑛×𝑦𝑐_𝑛(𝑧)𝑒_𝑛𝐶(𝑡)^𝑡₀𝑒_𝑛𝐶(−𝜏)^𝐿𝑥₀𝑐_𝑛(𝑢)×^𝐿𝑦₀𝑠_𝑛(𝑣)^𝐿𝑧₀𝑐_𝑛(𝑤)×𝑔_𝐿(𝑢,𝑣,𝑤,𝑇)𝜕𝐶₀₁(𝑢,𝑣,𝑤,𝜏)+𝐶𝜕𝑣_𝛾00(𝑢,𝑣,𝑤,𝜏)𝑃^𝛾(𝑢,𝑣,𝑤,𝑇)𝜕𝐶₁₀(𝑢,𝑣,𝑤,𝜏)𝜕𝑣+𝐶₁₀()𝐶𝑢,𝑣,𝑤,𝜏^𝛾−100(𝑢,𝑣,𝑤,𝜏)𝑃^𝛾()×𝑢,𝑣,𝑤,𝑇𝜕𝐶₀₀(𝑢,𝑣,𝑤,𝜏)𝜕𝑣+𝑔_𝐿()×𝑢,𝑣,𝑤,𝑇𝜕𝐶₀₀(𝑢,𝑣,𝑤,𝜏)𝐶𝜕𝑣^𝛾00(𝑢,𝑣,𝑤,𝜏)𝑃^𝛾()−𝑢,𝑣,𝑤,𝑇𝑑𝑤𝑑𝑣𝑑𝑢𝑑𝜏2𝜋𝐷_0𝐶𝐿_𝑥𝐿_𝑦𝐿²_𝑧^∞_𝑛=1𝑛𝑐_𝑛(𝑥)𝑐_𝑛𝑦𝑐_𝑛(𝑧)𝑒_𝑛𝐶(𝑡)×^𝑡₀𝑒_𝑛𝐶(−𝜏)𝐹_𝑛𝐶^𝐿𝑥₀𝑠_𝑛(𝑢)×^𝐿𝑦₀𝑐_𝑛(𝑣)^𝐿𝑧₀𝑐_𝑛(𝑤)×𝑔_𝐿()×𝑢,𝑣,𝑤,𝑇𝜕𝐶₀₁(𝑢,𝑣,𝑤,𝜏)+𝜕𝑤𝜕𝐶₁₀(𝑢,𝑣,𝑤,𝜏)×𝐶𝜕𝑤_𝛾00(𝑢,𝑣,𝑤,𝜏)𝑃^𝛾(𝑢,𝑣,𝑤,𝑇)+𝐶₁₀()×𝐶𝑢,𝑣,𝑤,𝜏_𝛾−100(𝑢,𝑣,𝑤,𝜏)𝑃^𝛾(𝑢,𝑣,𝑤,𝑇)𝜕𝐶₀₀(𝑢,𝑣,𝑤,𝜏)+𝜕𝑤𝜕𝐶₀₀(𝑢,𝑣,𝑤,𝜏)𝑔𝜕𝑤^𝐿()×𝐶𝑢,𝑣,𝑤,𝑇_𝛾00(𝑢,𝑣,𝑤,𝜏)𝑃^𝛾(𝑉𝑢,𝑣,𝑤,𝑇)𝑑𝑤𝑑𝑣𝑑𝑢𝑑𝜏,₀𝑥,𝑦,𝑧,𝑡=−2𝜋𝐷_0𝑉𝐿²_𝑥𝐿_𝑦𝐿_𝑧^∞_𝑛=1𝑛𝐹_𝑛𝑉𝑐_𝑛(𝑥)𝑐_𝑛𝑦𝑐_𝑛(𝑧)×𝑒_𝑛𝑉(𝑡)^𝑡₀𝑒_𝑛𝑉(−𝜏)^𝐿𝑥₀𝑠_𝑛(𝑢)^𝐿𝑦₀𝑐_𝑛(𝑣)×^𝐿𝑧₀𝑐_𝑛(𝑤)𝐷_𝑉𝑆(𝑢,𝑣,𝑤,𝑇)×()𝑉𝑘𝑇𝜕𝜇𝑢,𝑣,𝑤,𝜏𝜕𝑢𝑑𝑤𝑑𝑣𝑑𝑢𝑑𝜏−2𝜋𝐷_0𝑉𝐿_𝑥𝐿²_𝑦𝐿_𝑧×^∞_𝑛=1𝐹_𝑛𝑉𝑐_𝑛(𝑥)𝑐_𝑛𝑦𝑐_𝑛(𝑧)^𝑡₀𝑒_𝑛𝐶()×−𝜏_𝐿𝑥0𝑐_𝑛(𝑢)^𝐿𝑦₀𝑠_𝑛(𝑣)^𝐿𝑧₀𝑐_𝑛(𝑤)×𝐷_𝑉𝑆(𝑢,𝑣,𝑤,𝑇)()𝑉𝑘𝑇𝜕𝜇𝑢,𝑣,𝑤,𝜏𝜕𝑣𝑑𝑤𝑑𝑣𝑑𝑢𝑑𝜏×𝑛𝑒_𝑛𝑉(𝑡)−2𝜋𝐷_0𝑉𝐿_𝑥𝐿_𝑦𝐿²_𝑧×^∞_𝑛=1𝑛𝐹_𝑛𝑉𝑐_𝑛(𝑥)𝑐_𝑛𝑦𝑐_𝑛(𝑧)𝑒_𝑛𝑉(𝑡)×^𝑡₀𝑒_𝑛𝑉(−𝜏)^𝐿𝑥₀𝑐_𝑛(𝑢)^𝐿𝑦₀𝑐_𝑛(𝑣)×^𝐿𝑧₀𝑠_𝑛(𝑤)𝐷_𝑉𝑆(𝑢,𝑣,𝑤,𝑇)×()𝑉𝑘𝑇𝜕𝜇𝑢,𝑣,𝑤,𝜏𝑉𝜕𝑤𝑑𝑤𝑑𝑣𝑑𝑢𝑑𝜏,_𝑖𝑥,𝑦,𝑧,𝑡=−2𝜋𝐷_0𝑉𝐿²_𝑥𝐿_𝑦𝐿_𝑧^∞_𝑛=1𝑛𝐹_𝑛𝑉𝑐_𝑛(𝑥)𝑐_𝑛𝑦𝑐_𝑛(𝑧)×𝑒_𝑛𝑉(𝑡)^𝑡₀𝑒_𝑛𝑉(−𝜏)^𝐿𝑥₀𝑠_𝑛(𝑢)^𝐿𝑦₀𝑐_𝑛(𝑣)×^𝐿𝑧₀𝑐_𝑛(𝑤)𝑔_𝑉(𝑢,𝑣,𝑤,𝑇)𝜕𝑉_𝑖−1(𝑢,𝑣,𝑤,𝜏)−𝜕𝑢𝑑𝑤𝑑𝑣𝑑𝑢𝑑𝜏2𝜋𝐷_0𝑉𝐿_𝑥𝐿²_𝑦𝐿_𝑧^∞_𝑛=1𝐹_𝑛𝑉𝑐_𝑛(𝑥)×𝑐_𝑛𝑦𝑐_𝑛(𝑧)^𝑡₀𝑒_𝑛𝑉(−𝜏)^𝐿𝑥₀𝑐_𝑛(𝑢)×^𝐿𝑦₀𝑠_𝑛(𝑣)^𝐿𝑧₀𝑐_𝑛(𝑤)𝑔_𝑉()×𝑢,𝑣,𝑤,𝑇𝜕𝑉_𝑖−1(𝑢,𝑣,𝑤,𝜏)𝜕𝑣𝑑𝑤𝑑𝑣𝑑𝑢𝑑𝜏×𝑛𝑒_𝑛𝑉(𝑡)−2𝜋𝐷_0𝑉𝐿_𝑥𝐿_𝑦𝐿²_𝑧^∞_𝑛=1𝑛𝐹_𝑛𝑉𝑐_𝑛(𝑥)𝑐_𝑛𝑦𝑐_𝑛(𝑧)×𝑒_𝑛𝑉(𝑡)^𝑡₀𝑒_𝑛𝑉(−𝜏)^𝐿𝑥₀𝑐_𝑛(𝑢)×^𝐿𝑦₀𝑐_𝑛(𝑣)^𝐿𝑧₀𝑠_𝑛(𝑤)𝑔_𝑉()×𝑢,𝑣,𝑤,𝑇𝜕𝑉_𝑖−1(𝑢,𝑣,𝑤,𝜏)𝜕𝑤𝑑𝑤𝑑𝑣𝑑𝑢𝑑𝜏. (A.3)

Acknowledgments

This work was partly supported by the grant of the President of Russia (Project no. MK-548.2010.2).