**1. Introduction**

The conventional theory of heat conduction starts from the Fourier law and, coupled with the energy conservation law, implies the conventional parabolic heat conduction equation. In bodies exhibiting complex features, the usual Fourier law and heat conduction equation are not precise enough, and nonstandard theories, in which these equations are substituted by extended (nonlocal) equations, are elaborated.

Time nonlocal equations give an account of memory. Gurtin and Pipkin [1] suggested the general relation between the heat flux vector **q** and the temperature gradient:

$$\mathbf{q} = -k \int\_0^\infty K(\tau) \,\mathrm{grad}\, T(t-\tau) \,\mathrm{d}\tau. \tag{1}$$

In Equation (1), *k* can be treated as the generalized thermal conductivity, *T* denotes temperature, *t* stands for time, and *<sup>K</sup>*(*τ*) is the weight function.

Nigmatullin [2,3] proposed a similar version of Equation (1):

$$\mathbf{q} = -k \int\_0^t K(t - \tau) \,\mathrm{grad}\, T(\tau) \,\mathrm{d}\tau,\tag{2}$$

and obtained for temperature the equation with memory:

$$\frac{\partial T}{\partial t} = a \int\_0^t K(t - \tau) \,\Delta T(\tau) \,\mathrm{d}\tau. \tag{3}$$

In Equation (3), *a* can be treated as the generalized thermal diffusivity coefficient.

Different choices of the weight function *K*(*t* − *τ*) ("full sclerosis", "short-tail" memory, "middle-tail" memory, "long-tail" memory, and "full memory") were analized in [4–8].

Though different time nonlocal generalizations of the Fourier law have a long history, they still attract the attention of researchers (see, for example, Reference [9], and the extensive discussion in [10]).

Fractional calculus (the theory of integrals and derivatives of non-integer order) has a large number of uses in engineering, geophysics, physics, geology, chemistry, rheology, biology, finance, and medicine (see [11–21], among many others).

Equation (2), written with the "long-tail" power kernel [4–8], can be represented in terms of fractional calculus:

$$\mathbf{q}(t) = -kD\_{RL}^{1-\underline{u}} \text{grad } T(t), \qquad 0 < \underline{u} \le 1,\tag{4}$$

$$\mathbf{q}(t) = -kI^{\mathfrak{n}-1} \text{grad } T(t), \qquad 1 < \mathfrak{a} \le 2. \tag{5}$$

In Equations (4) and (5), *Iα f*(*t*) and *<sup>D</sup><sup>α</sup>RL f*(*t*) are the Riemann–Liouville fractional integral and derivative of the order *α* [22–24]:

$$d^a f(t) = \frac{1}{\Gamma(a)} \int\_0^t (t - \tau)^{a-1} f(\tau) \, d\tau, \qquad a > 0,\tag{6}$$

$$D\_{RL}^{a}f(t) = \frac{\mathbf{d}^{\rm m}}{\mathbf{d}t^{m}} \left[ \frac{1}{\Gamma(m-a)} \int\_{0}^{t} (t-\tau)^{m-a-1} f(\tau) \,d\tau \right], \qquad m-1 < a < m. \tag{7}$$

In Equations (6) and (7), <sup>Γ</sup>(*α*) denotes the Gamma function.

It should be emphasized that in fractional calculus, there is no sharp boundary between integrals and derivatives. This is why some authors [11,22] do not use a separate notation for the fractional integral *Iα f*(*t*), which is designated as *D*−*<sup>α</sup> RL f*(*t*) with *α* > 0. Making use of this notation, Equations (4) and (5) can be depicted as one time-nonlocal relation:

$$\mathbf{q}(t) = -kD\_{RL}^{1-a}\mathbf{grad}\,T(t), \qquad 0 < a \le 2\,. \tag{8}$$

Coupled with the law of conservation of energy, the constitutive relation (8) implies the time fractional heat conduction equation:

$$\frac{\partial^a T}{\partial t^a} = a \,\Delta T, \qquad 0 < a \le 2,\tag{9}$$

with the Caputo fractional derivative:

$$\frac{\mathrm{d}^{\mathrm{af}}f(t)}{\mathrm{d}t^{\alpha}} = \frac{1}{\Gamma(m-a)} \int\_{0}^{t} (t-\tau)^{m-a-1} \frac{\mathrm{d}^{m}f(\tau)}{\mathrm{d}\tau^{m}} d\tau, \qquad m-1 < a < m. \tag{10}$$

The interested reader can find the details of obtaining the time fractional heat conduction Equation (9) from the constitutive relation (8) in [25].

Equations (4) and (9) for 0 < *α* < 1 describe the so-called subdiffusion (slow diffusion), which is characterized by the mean-squared displacement *x*<sup>2</sup> ∼ *t<sup>α</sup>*.

Recall the rules of the Laplace transform of fractional integrals and derivatives [22–24]:

$$
\mathcal{L}\left\{l^a f(t)\right\} = \frac{1}{s^a} f^\*(s),
\tag{11}
$$

$$\mathcal{L}\left\{\mathcal{D}\_{RL}^{a}f(\mathbf{t})\right\} = s^{a}f^{\*}(\mathbf{s}) - \sum\_{k=0}^{m-1} D^{k}I^{m-a}f(0^{+})s^{m-1-k}, \qquad m-1 < a < m,\tag{12}$$

*Symmetry* **2019**, *11*, 800

$$\mathcal{L}\left\{\frac{\mathbf{d}^a f(t)}{\mathbf{d}t^a}\right\} = s^a f^\*(s) - \sum\_{k=0}^{m-1} f^{(k)}(0^+) s^{a-1-k}, \qquad m-1 < a < m. \tag{13}$$

In Equations (11)–(13), the asterisk marks the transform, and *s* stands for the Laplace transform variable.

A variety of boundary conditions for the time fractional heat conduction equation were analyzed in [8,20,26,27]. If the surfaces of two bodies are in perfect thermal contact, the temperatures at their joint boundary *S* and the heat fluxes through this boundary are the same for both bodies, and we arrive at the following boundary conditions:

$$\left.T\_1\right|\_S = T\_2\Big|\_{S'}\tag{14}$$

$$k\_1 D\_{RL}^{1-a} \frac{\partial T\_1}{\partial n} \bigg|\_{S} = k\_2 D\_{RL}^{1-\beta} \frac{\partial T\_2}{\partial n} \bigg|\_{S} \qquad 0 < a \le 2, \qquad 0 < \beta \le 2,\tag{15}$$

where subscripts 1 and 2 mark bodies 1 and 2, respectively, and *n* is the common normal at the joint surface. In Equation (15), it is assumed that in the general case, the heat conduction in one of the bodies is described by the equation with the Caputo derivative of the order *α*, whereas in the second body, the process is governed by the equation with the Caputo derivative of the order *β*.

In previous publications [27–29], the problems of fractional heat conduction in composed bodies were investigated for one spatial coordinate. In the present paper, heat conduction with the Caputo fractional derivative is considered in two joint half-planes under the conditions of perfect thermal contact in the case of two spatial coordinates. The fundamental solutions to the Cauchy and source problems are studied. The Fourier and Laplace transforms are employed. The Fourier transforms are inverted analytically, and the Laplace transform is inverted numerically, employing the Gaver–Stehfest method [30–33]. The interested reader is also referred to the recent paper [34] on numerical inversion of the Laplace transform for solving fractional deifferenetial equations, where additional references can be found.

### **2. The Fundamental Solution to the Cauchy Problem**

Consider the time fractional heat conduction equations in two half-planes:

$$\begin{split} \frac{\partial^{a}T\_{1}(\mathbf{x},y,t)}{\partial t^{a}} &= a\_{1} \left[ \frac{\partial^{2}T\_{1}(\mathbf{x},y,t)}{\partial \mathbf{x}^{2}} + \frac{\partial^{2}T\_{1}(\mathbf{x},y,t)}{\partial y^{2}} \right], \\ 0 &< \mathbf{x} < \infty, \quad -\infty < y < \infty, \quad 0 < t < \infty, \quad 0 < a \le 2, \end{split} \tag{16}$$

$$\frac{\partial^{\beta}T\_{2}(\mathbf{x},\mathbf{y},t)}{\partial t^{\beta}} = a\_{2} \left[ \frac{\partial^{2}T\_{2}(\mathbf{x},\mathbf{y},t)}{\partial \mathbf{x}^{2}} + \frac{\partial^{2}T\_{2}(\mathbf{x},\mathbf{y},t)}{\partial \mathbf{y}^{2}} \right],\tag{17}$$
 
$$-\infty < \mathbf{x} < 0, \quad -\infty < \mathbf{y} < \infty, \quad 0 < t < \infty, \quad 0 < \beta \le 2.$$

We assume the initial conditions:

$$t = 0: \qquad T\_1(\mathbf{x}, y, t) = p\wp \delta\left(\mathbf{x} - l\right)\delta\left(y\right), \quad 0 < \mathbf{x} < \infty, \quad -\infty < y < \infty, \quad 0 < a \le 2,\tag{18}$$

$$\mathbf{t} = \mathbf{0}: \qquad \frac{\partial T\_1(\mathbf{x}, y, t)}{\partial t} = \mathbf{0}, \quad \mathbf{0} < \mathbf{x} < \infty, \qquad -\infty < y < \infty, \quad 1 < a \le 2,\tag{19}$$

$$\mathbf{x}, t = 0: \qquad T\_2(\mathbf{x}, y, t) = 0, \qquad -\infty < \mathbf{x} < 0, \quad -\infty < y < \infty, \quad 0 < \beta \le 2,\tag{20}$$

$$\mathbf{t} = \mathbf{0}: \qquad \frac{\partial T\_2(\mathbf{x}, y, t)}{\partial t} = \mathbf{0}, \quad -\infty < \mathbf{x} < \mathbf{0}, \quad -\infty < y < \infty, \quad 1 < \beta \le 2,\tag{21}$$

and the conditions of perfect thermal contact:

$$\mathbf{x} = \mathbf{0}: \qquad T\_1(\mathbf{x}, y, t) = T\_2(\mathbf{x}, y, t), \quad -\infty < y < \infty, \quad 0 < t < \infty, \quad 0 < \mathbf{a} \le 2, \quad 0 < \boldsymbol{\beta} \le 2,\tag{22}$$

$$\mathbf{x} = \mathbf{0}: \qquad k\_1 D\_{RL}^{1 - \alpha} \frac{\partial T\_1(\mathbf{x}, y, t)}{\partial \mathbf{x}} = k\_2 D\_{RL}^{1 - \beta} \frac{\partial T\_2(\mathbf{x}, y, t)}{\partial \mathbf{x}},$$

$$-\infty < y < \infty, \quad 0 < t < \infty, \quad 0 < \mathbf{a} \le 2, \quad 0 < \beta \le 2. \tag{23}$$

The temperatures should also satisfy the zero conditions at infinity:

$$\lim\_{\mathbf{x}\to\infty}T\_1\left(\mathbf{x},\mathbf{y},t\right)=0,\qquad \lim\_{\mathbf{y}\to\pm\infty}T\_1\left(\mathbf{x},\mathbf{y},t\right)=0,\tag{24}$$

$$\lim\_{\chi \to -\infty} T\_2\left(x, y, t\right) = 0, \qquad \lim\_{y \to \pm \infty} T\_2\left(x, y, t\right) = 0,\tag{25}$$

$$\lim\_{t \to \infty} T\_1\left(\mathbf{x}, y, t\right) = 0, \qquad \lim\_{t \to \infty} T\_2\left(\mathbf{x}, y, t\right) = 0. \tag{26}$$

We introduce the unknown function:

$$\log \left( y, t \right) = k\_1 D\_{RL}^{1 - a} \frac{\partial T\_1(x, y, t)}{\partial x} \Big|\_{x = 0} = k\_2 D\_{RL}^{1 - \beta} \frac{\partial T\_2(x, y, t)}{\partial x} \Big|\_{x = 0}. \tag{27}$$

In the subsequent text, we employ the following cos-Fourier transforms with respect to the coordinate *x*: for *x* > 0:

$$\mathcal{F}\_{\mathfrak{C}}\left\{f(\mathbf{x})\right\} = \tilde{f}(\mathfrak{f}) = \int\_{0}^{\infty} f(\mathbf{x}) \cos \left(\mathbf{x}\_{\mathfrak{V}}^{\mathfrak{r}}\right) d\mathbf{x}\_{\mathfrak{r}} \tag{28}$$

$$\mathcal{F}\_{\varepsilon}^{-1}\left\{\widetilde{f}(\xi)\right\} = f(\mathbf{x}) = \frac{2}{\pi} \int\_{0}^{\infty} \widetilde{f}(\xi) \cos \left(\mathbf{x}\xi\right) \,\mathrm{d}\xi,\tag{29}$$

$$\mathcal{F}\_c \left\{ \frac{\mathrm{d}^2 f(\mathbf{x})}{\mathrm{d} \mathbf{x}^2} \right\} = -\xi^2 \tilde{f}(\xi) - \left. \frac{\mathrm{d} f(\mathbf{x})}{\mathrm{d} \mathbf{x}} \right|\_{\mathbf{x} = 0^+} \, ^{ \prime} \tag{30}$$

and for *x* < 0:

$$\mathcal{F}\_{\varepsilon}\left\{f(\mathbf{x})\right\} = \tilde{f}(\xi) = \int\_{-\infty}^{0} f(\mathbf{x}) \cos\left(\mathbf{x}\xi\right) d\mathbf{x},\tag{31}$$

$$\mathcal{F}\_{\mathfrak{c}}^{-1}\left\{\widetilde{f}(\mathfrak{f})\right\} = f(\mathfrak{x}) = \frac{2}{\pi} \int\_{-\infty}^{0} \widetilde{f}(\mathfrak{f}) \, \cos \left(\mathfrak{x} \mathfrak{f}\right) \, \mathrm{d}\mathfrak{f},\tag{32}$$

$$\mathcal{F}\_{\mathfrak{C}}\left\{\frac{\mathrm{d}^{2}f(\mathbf{x})}{\mathrm{d}\mathbf{x}^{2}}\right\}=-\mathfrak{J}^{2}\widetilde{f}(\mathfrak{J})+\left.\frac{\mathrm{d}f(\mathbf{x})}{\mathrm{d}\mathbf{x}}\right|\_{\mathbf{x}=\mathbf{0}^{-}}\,.\tag{33}$$

Applying to the initial boundary value problem (16)–(27) the Laplace transform with respect to time *t*, the exponential Fourier transform with respect to the coordinate *y*, and the foregoing cos-Fourier transforms with respect to the coordinate *x*, in the transform domain we ge<sup>t</sup> is:

$$\stackrel{\circ}{T}\_{1}^{\circ \circ} \left( \stackrel{x}{\sharp}, \eta, s \right) = \frac{p\_0}{\sqrt{2\pi}} \cos \left( l \xi \right) \frac{s^{a-1}}{s^a + a\_1 \left( \xi^2 + \eta^2 \right)} - \frac{a\_1}{k\_1} \, \hat{\phi}^\* \left( \eta, s \right) \frac{s^{a-1}}{s^a + a\_1 \left( \xi^2 + \eta^2 \right)} \, \tag{34}$$

$$
\widetilde{\hat{T}}\_2^\*\left(\xi,\eta,s\right) = \frac{a\_2}{k\_2} \widetilde{\varphi}^\*\left(\eta,s\right) \frac{s^{\mathfrak{G}-1}}{s^{\mathfrak{G}} + a\_2\left(\xi^2 + \eta^2\right)}.\tag{35}
$$

All the Fourier transforms are marked by the tilde, *ξ* denotes the cos-Fourier transform variable, and *η* stands for the exponential Fourier transform variable.

Inversion of the cos-Fourier transforms with taking into account the integral below [35]:

$$\int\_0^\infty \frac{1}{\xi^2 + c^2} \cos \left( \mathbf{x}\_\xi^x \right) \mathbf{d}\_\xi^x = \frac{\pi}{2c} \exp \left( -c|\mathbf{x}| \right), \quad c > 0,\tag{36}$$

gives:

$$\begin{split} \bar{T}\_{1}^{\*}\left(\mathbf{x},\eta,s\right) &= \frac{p\_{0}}{\sqrt{2\pi}} \frac{s^{a-1}}{2a\_{1}} \frac{1}{\sqrt{\eta^{2}+s^{a}/a\_{1}}} \Bigg\{ \exp\left[-\sqrt{\eta^{2}+\frac{s^{a}}{a\_{1}}}\left(\mathbf{x}+l\right)\right] + \exp\left[-\sqrt{\eta^{2}+\frac{s^{a}}{a\_{1}}}\left|\mathbf{x}-l\right|\right] \Bigg\} \\ &- \frac{s^{a-1}}{k\_{1}} \bar{\boldsymbol{\eta}}^{\*}\left(\boldsymbol{\eta},s\right) \frac{1}{\sqrt{\eta^{2}+s^{a}/a\_{1}}} \exp\left(-\sqrt{\eta^{2}+\frac{s^{a}}{a\_{1}}}\mathbf{x}\right), \quad \mathbf{x}>0, \end{split} \tag{37}$$

$$\tilde{T}\_2^\*\left(\mathbf{x}, \eta, s\right) = \frac{s^{\mathfrak{G}-1}}{k\_2} \,\hat{\varphi}^\*\left(\eta, s\right) \frac{1}{\sqrt{\eta^2 + s^{\mathfrak{G}}/a\_2}} \exp\left(-\sqrt{\eta^2 + \frac{s^{\mathfrak{G}}}{a\_2}} |\mathbf{x}|\right), \quad \mathbf{x} < 0. \tag{38}$$

The perfect thermal contact boundary condition (19) written in the transform domain:

$$
\hat{T}\_1^\*(0, \eta, s) = \hat{T}\_2^\*(0, \eta, s) \tag{39}
$$

makes it possible to determine the function: *ϕ*∗ (*η*,*<sup>s</sup>*):

$$\widetilde{\varphi}^\*(\eta, s) = \frac{p\_0 k\_1 k\_2}{\sqrt{2\pi} a\_1} \frac{\sqrt{\eta^2 + s^\beta / a\_2}}{s^{\beta - a} k\_1 \sqrt{\eta^2 + s^a / a\_1} + k\_2 \sqrt{\eta^2 + s^\beta / a\_2}} \exp\left(-\sqrt{\eta^2 + \frac{s^a}{a\_1}} l\right), \tag{40}$$

and obtain the expressions for temperatures:

$$\begin{split} \widetilde{T}\_{1}^{\*}\left(\mathbf{x},\eta,s\right) &= \frac{p\_{0}}{\sqrt{2\pi}a\_{1}} \frac{s^{s-1}}{\sqrt{\eta^{2}+s^{s}/a\_{1}}} \bigg\{ \frac{1}{2} \exp\left[-\sqrt{\eta^{2}+\frac{s^{s}}{a\_{1}}}\left(\mathbf{x}+l\right)\right] \\ &+ \frac{1}{2} \exp\left[-\sqrt{\eta^{2}+\frac{s^{s}}{a\_{1}}}\left|\mathbf{x}-l\right|\right] \\ &- \frac{k\_{2}\sqrt{\eta^{2}+s^{s}/a\_{2}}}{s^{\theta-a}k\_{1}\sqrt{\eta^{2}+s^{s}/a\_{1}}+k\_{2}\sqrt{\eta^{2}+s^{\theta}/a\_{2}}} \exp\left[-\sqrt{\eta^{2}+\frac{s^{s}}{a\_{1}}}\left(\mathbf{x}+l\right)\right] \bigg\}, \quad \mathbf{x}>0, \\\\ \widetilde{T}\_{2}^{\*}\left(\mathbf{x},\eta,s\right) &= \frac{p\_{0}}{\sqrt{2\pi}a\_{1}} \frac{k\_{1}s^{\theta-1}}{s^{\theta-a}k\_{1}\sqrt{\eta^{2}+s^{s}/a\_{1}}+k\_{2}\sqrt{\eta^{2}+s^{\theta}/a\_{2}}} \\ &\times \exp\left(-\sqrt{\eta^{2}+\frac{s^{\theta}}{a\_{2}}}\left|\mathbf{x}\right| - \sqrt{\eta^{2}+\frac{s^{s}}{a\_{1}}}l\right), \quad \mathbf{x}<0. \end{split} \tag{42}$$

The inverse transforms give:

$$\begin{split} T\_{1}\left(\mathbf{x},y,t\right) &= \frac{p\_{0}}{2\pi a\_{1}} \mathcal{L}^{-1} \left\langle \int\_{-\infty}^{\infty} \frac{s^{a-1}}{\sqrt{\eta^{2} + s^{a}/a\_{1}}} \left\{ \frac{1}{2} \exp\left[-\sqrt{\eta^{2} + \frac{s^{a}}{a\_{1}}} \left(\mathbf{x} + l\right)\right] \right. \\ &+ \frac{1}{2} \exp\left[-\sqrt{\eta^{2} + \frac{s^{a}}{a\_{1}}} \left|\mathbf{x} - l\right|\right] \\ &- \frac{k\_{2}\sqrt{\eta^{2} + s^{\beta}/a\_{2}}}{s^{\beta-a}k\_{1}\sqrt{\eta^{2} + s^{a}/a\_{1}} + k\_{2}\sqrt{\eta^{2} + s^{\beta}/a\_{2}}} \exp\left[-\sqrt{\eta^{2} + \frac{s^{\kappa}}{a\_{1}}} \left(\mathbf{x} + l\right)\right] \right\} \cos\left(y\eta\right) \,\mathrm{d}\eta \right\rangle, \; \mathbf{x} > 0, \end{split} \tag{43}$$

$$\begin{split} T\_2\left(\mathbf{x}, y, t\right) &= \frac{p\_0 k\_1}{2\pi a\_1} \mathcal{L}^{-1} \left\{ \int\_{-\infty}^{\infty} \frac{s^{\beta - 1}}{s^{\beta - a} k\_1 \sqrt{\eta^2 + s^a/a\_1} + k\_2 \sqrt{\eta^2 + s^{\beta}/a\_2}} \\ &\times \exp\left( -\sqrt{\eta^2 + \frac{s^{\beta}}{a\_2}} |\mathbf{x}| - \sqrt{\eta^2 + \frac{s^{\pi}}{a\_1}} l \right) \cos\left( y\eta \right) \mathbf{d}\eta \right\}, \quad \mathbf{x} < 0. \end{split} \tag{44}$$

Temperatures in joint half-planes are presented in Figures 1–6 for varied combinations of parameters. Calculations have been carried out with the following nondimensional quantities:

$$\bar{\mathbf{x}} = \frac{\mathbf{x}}{l}, \quad \bar{y} = \frac{y}{l}, \quad \kappa = \frac{\sqrt{a\_1}}{l} \frac{t\_0^{a/2}}{l}, \quad \bar{T} = \frac{a\_1 t\_0^a}{p\_0} T,\tag{45}$$

where *t*0 is the characteristic time of the process. In all the figures, the values *α* = *β* = 0.5 and *κ* = 1 have been taken.

**Figure 1.** Behavior of the solution (43) and (44) versus time; *x*¯ = 1, *y* = 0, *k*1 = *k*2.

**Figure 2.** Behavior of the solution (43) and (44) versus time; *x*¯ = 1, *y* = 0, *a*1 = *a*2.

**Figure 3.** Dependence of the solution (43) and (44) on the coordinate *x*; *y*¯ = 0, ¯*t* = 1, *k*1 = *k*2.

**Figure 4.** Dependence of the solution (43) and (44) on the coordinate *x*; *y*¯ = 0, ¯*t* = 1, *a*1 = *a*2.

**Figure 5.** Dependence of the solution (43) and (44) on the coordinate *y*; *x*¯ = 1, ¯*t* = 1, *k*1 = *k*2.

**Figure 6.** Dependence of the solution (43) and (44) on the coordinate *y*; *x*¯ = 1, ¯*t* = 1, *a*1 = *a*2.

### **3. The Fundamental Solution to the Source Problem**

Next, we consider the time fractional heat conduction equation with the source term in the domain

*x* > 0:

$$\begin{split} \frac{\partial^{a}T\_{1}(\mathbf{x},\mathbf{y},t)}{\partial t^{a}} &= a\_{1} \left[ \frac{\partial^{2}T\_{1}(\mathbf{x},\mathbf{y},t)}{\partial \mathbf{x}^{2}} + \frac{\partial^{2}T\_{1}(\mathbf{x},\mathbf{y},t)}{\partial \mathbf{y}^{2}} \right] + \mathbf{g}\_{0} \delta\left(\mathbf{x}-l\right) \delta\left(\mathbf{y}\right) \delta\left(t\right), \\ &0 < \mathbf{x} < \infty, \quad -\infty < \mathbf{y} < \infty, \quad 0 < t < \infty, \quad 0 < a \le 2, \end{split} \tag{46}$$

and the corresponding equation in the region *x* < 0:

$$\begin{split} \frac{\partial^{\beta}T\_{2}(\mathbf{x},\mathbf{y},t)}{\partial t^{\beta}} &= a\_{2} \left[ \frac{\partial^{2}T\_{2}(\mathbf{x},\mathbf{y},t)}{\partial \mathbf{x}^{2}} + \frac{\partial^{2}T\_{2}(\mathbf{x},\mathbf{y},t)}{\partial \mathbf{y}^{2}} \right], \\ -\infty &< \mathbf{x} < 0, \quad -\infty < \mathbf{y} < \infty, \quad 0 < t < \infty, \quad 0 < \beta \le 2. \end{split} \tag{47}$$

We assume the zero initial conditions:

$$t = 0: \qquad T\_1(\mathbf{x}, \mathbf{y}, t) = 0, \qquad 0 < \mathbf{x} < \infty, \quad -\infty < \mathbf{y} < \infty, \qquad 0 < \mathbf{z} \le 2,\tag{48}$$

$$t = 0: \qquad \frac{\partial T\_1(\mathbf{x}, y, t)}{\partial t} = 0, \quad 0 < \mathbf{x} < \infty, \qquad -\infty < y < \infty, \quad 1 < a \le 2,\tag{49}$$

$$t = 0: \qquad T\_2(\mathbf{x}, y, t) = 0, \qquad -\infty < \mathbf{x} < 0, \quad -\infty < y < \infty, \quad 0 < \beta \le 2,\tag{50}$$

$$t = 0: \qquad \frac{\partial T\_2(\mathbf{x}, y, t)}{\partial t} = 0, \quad -\infty < \mathbf{x} < 0, \quad -\infty < y < \infty, \quad 1 < \beta \le 2,\tag{51}$$

and the conditions of perfect thermal contact at the joint boundary:

$$\mathbf{x} = \mathbf{0}: \qquad T\_1(\mathbf{x}, y, t) = T\_2(\mathbf{x}, y, t), \quad -\infty < y < \infty, \quad 0 < t < \infty, \quad 0 < a \le 2, \quad 0 < \beta \le 2,\tag{52}$$

$$\mathbf{x} = \mathbf{0}: \qquad k\_1 D\_{RL}^{1 - a} \frac{\partial T\_1(\mathbf{x}, y, t)}{\partial \mathbf{x}} = k\_2 D\_{RL}^{1 - \beta} \frac{\partial T\_2(\mathbf{x}, y, t)}{\partial \mathbf{x}},$$

$$-\infty < y < \infty, \quad 0 < t < \infty, \quad 0 < a \le 2, \quad 0 < \beta \le 2.\tag{53}$$

Using the integral transform technique, we finally obtain:

$$\begin{split} T\_1(\mathbf{x}, y, t) &= \frac{\mathcal{G}\_0}{2\pi a\_1} \mathcal{L}^{-1} \left\langle \int\_{-\infty}^{\infty} \frac{1}{\sqrt{\eta^2 + s^2/a\_1}} \left\{ \frac{1}{2} \exp\left[ -\sqrt{\eta^2 + \frac{s^2}{a\_1}} \left( \mathbf{x} + l \right) \right] \right\} \right. \\ &+ \frac{1}{2} \exp\left[ -\sqrt{\eta^2 + \frac{s^2}{a\_1}} \left| \mathbf{x} - l \right| \right] \\ &- \frac{k\_2 \sqrt{\eta^2 + s^2/a\_2}}{s^{\beta - a} k\_1 \sqrt{\eta^2 + s^2/a\_1} + k\_2 \sqrt{\eta^2 + s^2} / a\_2} \exp\left[ -\sqrt{\eta^2 + \frac{s^2}{a\_1}} \left( \mathbf{x} + l \right) \right] \right\} \cos \left( y\eta \right) \, \mathrm{d}\eta \Big], \, \mathbf{x} > 0, \\\\ T\_2 \left( \mathbf{x}, y, t \right) &= \frac{\mathcal{G} \mathbb{B}^{\mathbf{k}\_1}}{2 \pi a\_1} \mathcal{L}^{-1} \left\{ \int\_{-\infty}^{\infty} \frac{s^{\beta - a}}{s^{\beta - a} k\_1 \sqrt{\eta^2 + s^2/a\_1} + k\_2 \sqrt{\eta^2 + s^2} / a\_2} \\ &\times \exp\left( -\sqrt{\eta^2 + \frac{s^{\beta}}{a\_2}} \left| \mathbf{x} \right| - \sqrt{\eta^2 + \frac{s^{\beta}}{a\_1}} \, l \right) \cos \left( y\eta \right) \, \mathrm{d}\eta \right\}, \quad \mathbf{x} < 0. \end{split} \tag{55}$$

The numerical results according to Equations (54) and (55) with *α* = *β* = 0.5, *κ* = 1 are presented in Figures 7–12. The nondimensional temperature is introduced as:

$$
\bar{T} = \frac{a\_1 t\_0}{\mathcal{g}\_0},
\tag{56}
$$

while other nondimensional quantities are the same as in Equation (45).

**Figure 7.** Behavior of the solution (54) and (55) versus time; *x*¯ = 1, *y* = 0, *k*1 = *k*2.

**Figure 8.** Behavior of the solution (54) and (55) versus time; *x*¯ = 1, *y* = 0, *a*1 = *a*2.

**Figure 9.** Dependence of the solution (54) and (55) on the coordinate *x*; *y*¯ = 0, ¯*t* = 1, *k*1 = *k*2.

**Figure 10.** Dependence of the solution (54) and (55) on the coordinate *x*; *y*¯ = 0, ¯*t* = 1, *a*1 = *a*2.

**Figure 11.** Dependence of the solution (54) and (55) on the coordinate *y*; *x*¯ = 1, ¯*t* = 1, *k*1= *k*2.

**Figure 12.** Dependence of the solution (54) and (55) on the coordinate *y*; *x*¯ = 1, ¯*t* = 1, *a*1 = *a*2.
