## 2 Analytical Solution  
To solve the PDE (2), we use the well-known reduction technique, where we define a new variable \(\eta=\frac{\pi}{i\theta}~{}\in\mathbb{R}\), which is a combination of the spatial and temporal variable. Then, we try to find the solution \(u(x,~{}t)\) with the self-similar Ansatz in the form of \(u(x,~{}t)=t^{-\alpha}f\big{(}x/t^{\beta}\big{)}\), where \(\alpha\) and \(\beta\) are arbitrary real constants, and \(f(\eta)\) is the shape function with existing first and second continuous derivatives with respect to \(\eta\). Evaluating the first temporal and second spatial derivative of this Ansatz and substituting these back to the PDE (2) yields  
\[-\alpha t^{-\alpha-1}f(\eta)-\beta\eta t^{-\alpha-1}f^{\prime}(\eta)=D\Big{(} m\eta^{m-1}t^{-\beta}t^{-\alpha-\beta}f^{\prime}(\eta)+\eta^{m}t^{-\alpha-2 \beta}f^{\prime}\Big{)} \tag{3}\]  
where prime means derivation in respect to the reduced variable \(\eta\). To complete the reduction mechanism, the explicit presence of the time variable has to be eliminated. From this requirement, we arrive at the following constraints for the parameters, \(\alpha\) = arbitrary real number, \(\beta=1/2\), while \(m\) remains an arbitrary real parameter. Now, we have the ordinary differential equation (ODE) for \(f(\eta)\)  
\[f(\eta) = \mathrm{e}\,\frac{\mathrm{e}^{\frac{\eta^{-\alpha+2}}{4D(m-2)}} \left[\frac{c_{1}}{\sqrt{\eta}}\cdot M_{\frac{4\alpha-1}{2m-4}},\,\frac{m-1}{ 2m-4}\left(\frac{\eta^{2-m}}{2D(m-2)}\right)\ +\ \frac{c_{2}}{\sqrt{\eta}}\cdot W_{\frac{4\alpha-1}{2m-4}},\,\frac{m-1}{2m-1} \left(\frac{\eta^{2-m}}{2D(m-2)}\right)\right]\] \[\text{where }M\text{ and }W\text{ are the Whittaker functions \@@cite[cite]{[\@@bibref{}{Buffaker}{}{}]},\@@cite[cite]{[\@@bibref{}{Buffaker}{}{}]}. There is a connection to the former results, which can be seen if one checks the following formulas which express the Whittaker functions [36] in terms of the Kummer functions \(M\) and \(U\):  
\[M_{\kappa,\mu}(z)=\mathrm{e}^{-\frac{z}{2}z^{\mu+\frac{1}{2}}M\bigg{(}\mu- \kappa+\frac{1}{2},1+2\mu;z\bigg{)}\,\ W_{\kappa,\mu}(z)=\mathrm{e}^{-\frac{z}{2}z^{\mu+\frac{1}{2}}U\bigg{(}\mu- \kappa+\frac{1}{2},1+2\mu;z\bigg{)} \tag{5}\]  
Due to the exponential factor in Equation (5), the Whittaker functions have a quicker decay than the Kummer functions. The shape functions for some parameter values are presented in Figure 1. The time development of the function \(u\) for two given parameter sets are shown in Figures 2 and 3. It can be shown with a careful parameter analysis that for negative values of \(\alpha\), the solutions have an exponential growth at large time and spatial coordinates which may be considered non-physical; thus, we exclude them from further numerical investigations. Diffusion processes where the concentration or the number of particles explodes violate energy and matter conservation laws. On the other hand, for some values of the parameters and variables, the values of \(u\) can be complex, which may be considered non-physical. We also note that because our analytical solution is valid on the whole real axis, the boundary conditions need to be specified only when the analytical solution is going to be reproduced by numerical methods. Figure 1: The solution of Equation (3) with \(D=1\); \(\alpha=1/2\); \(c_{1}=1\); \(c_{2}=0\): the black, red, blue, and green lines are for (\(m=0\); \(1/2\); \(1\); \(5/2\)). In our numerical experiments, we always choose \(D=1,\;\;c_{1}=0\) and set the value of \(c_{2}\) to be a normalization constant. It means that we are going to numerically reproduce the following reference solution:  
\[u(x,t)=t^{-\alpha}f\bigg{(}\frac{x}{t^{\beta}}\bigg{)}\;=\;c_{2}\sqrt{\frac{t^{ \frac{1}{2}-2\alpha}}{x}}\cdot\mathrm{e}^{\frac{(x/\sqrt{t})^{2}-m}{4(m-2)} \cdot}W_{\frac{4m-1}{2m-4},\;\frac{m-1}{2m-4}}\left(\frac{\left(x/\sqrt{t} \right)^{2-m}}{2(m-2)}\right) \tag{6}\]  
In Figures 2 and 3, we present 3D plots of this \(u\) function for those two cases, which will be reproduced numerically in the latter sections.