**3. Stability Analysis**

The system of Equations (8)–(10) along with the boundary conditions (11), is capable of generating more than one solution and ultimately permits the requirement analysis of the flow to identify the reliable and feasible solution. Going through the outstanding work done by Merkin [53] and Merrill et al. [54] in the stability analysis, the application of an unstable form of the boundary layer problem is analyzed by using the time variable and a dimensionless time variable denoted by τ [55]. Next, we consider the following new similarity variables:

 $\mu = \frac{cx}{1 + at}$  $\frac{\partial f}{\partial \eta}(\eta, \tau)$ ,  $\upsilon = \frac{cy}{1 + at}$  $\frac{\partial g}{\partial \eta}(\eta, \tau)$ ,  $w = -\sqrt{\frac{\upsilon c}{1 + at}}[f(\eta, \tau) + g(\eta, \tau)]$ ,  $\theta(\eta) = \frac{(T - T\_{\text{ov}})}{(T\_w - T\_{\text{ov}})}$ ,  $\eta = \sqrt{\frac{c}{\nu(1 + at)}}z$ ,  $\tau = \frac{ct}{1 + at}$ .

Using the similarity variables of Equation (15) into Equations (8)–(10), the altered differential equations and the boundary conditions are as follows:

$$\begin{cases} \frac{\mu\_{\rm hf}/\mu\_f}{\rho\_{\rm hf}/\rho\_f} \frac{\partial^3 f}{\partial \eta^3} + \left( A \frac{1}{2} \eta + f + g \right) \frac{\partial^2 f}{\partial \eta^2} + A \frac{\partial f}{\partial \eta} - \left( \frac{\partial f}{\partial \eta} \right)^2 - A \left( \lambda + \chi \right) + \left( \lambda + \chi \right)^2\\ -M \left( \frac{\partial f}{\partial \eta} - \lambda - \chi \right) - \left( 1 - A \tau \right) \frac{\partial^2 f}{\partial \eta \partial \tau} = 0, \end{cases} \tag{16}$$

$$\begin{cases} \frac{\mu\_{\rm mf}/\mu\_f}{\rho\_{\rm mf}/\rho\_f} \frac{\partial^3 \chi}{\partial \eta^3} + \left( A\_2^1 \eta + f + g \right) \frac{\partial^2 \chi}{\partial \eta^2} + A \frac{\partial \chi}{\partial \eta} - \left( \frac{\partial \chi}{\partial \eta} \right)^2 - A \left( \lambda - \gamma \right) + \left( \lambda - \gamma \right)^2\\ -M \left( \frac{\partial \chi}{\partial \eta} - \lambda + \gamma \right) - \left( 1 - A \tau \right) \frac{\partial^2 \chi}{\partial \eta^2 \tau} = 0, \end{cases} \tag{17}$$

$$\frac{1}{\Pr} \frac{k\_{\rm hnf}/k\_f}{\rho \mathbb{C}\_{\rm p\rm hnf}/\rho \mathbb{C}\_{pf}} \frac{\partial^2 \theta}{\partial \eta^2} + \left( A \Big(\frac{1}{2}\eta - \tau\Big) + f + g + 1 \right) \frac{\partial \theta}{\partial \eta} - (1 - A\tau) \frac{\partial \theta}{\partial \tau} = 0,\tag{18}$$

$$\begin{array}{l} f(0,\tau) = \mathcal{g}(0,\tau) = 0, \frac{\partial f}{\partial \eta}(0,\tau) = \frac{\partial \mathcal{g}}{\partial \eta}(0,\tau) = \varepsilon, \; \mathcal{O}(0,\tau) = 1, \\\ \frac{\partial f}{\partial \eta}(\eta,\tau) \to \lambda + \gamma, \frac{\partial \mathcal{g}}{\partial \eta}(\eta,\tau) \to \lambda - \gamma, \mathcal{O}(\eta,\tau) \to 0, \; \text{as } \eta \to \infty. \end{array} \tag{19}$$

Weidman et al. [55] highlighted that the stability of solutions is introduced by determining the system's decay or initial growth. This can be achieved via considering the following perturbation expressions of the primary flow *f* = *f*0(η), *g* = *g*0(η) and θ = <sup>θ</sup>0(η) with the resulting equation:

$$\begin{cases} f(\eta,\tau) = f\_0(\eta) + e^{-\omega\tau} F(\eta), \mathcal{g}(\eta,\tau) = \mathcal{g}\_0(\eta) + e^{-\omega\tau} G(\eta), \\ \theta(\eta,\tau) = \theta\_0(\eta) + e^{-\omega\tau} H(\eta), \end{cases} \tag{20}$$

where ω is an unknown parameter of the eigenvalue, *<sup>F</sup>*(η), *<sup>G</sup>*(η) and *<sup>H</sup>*(η) are comparatively slight to *f*0(η), *g*0(η) and <sup>θ</sup>0(η). Substituting (20) into Equations (16)–(18) we attained the following system of equations:

$$\begin{cases} \frac{\mu\_{\rm Imf} \wedge \mu\_f}{\rho\_{\rm Imf} \wedge \rho\_f} \frac{\partial^3 F}{\partial \eta^3} + \left( A\_2^1 \eta + f\_0 + \mathbf{g}\_0 \right) \frac{\partial^2 F}{\partial \eta^2} + \left( A - M + (1 - A \tau) \omega - 2 \frac{\partial f\_0}{\partial \eta} \right) \frac{\partial F}{\partial \eta} \\ + (F + \mathbf{G}) \frac{\partial^2 f\_0}{\partial \eta^2} = 0, \end{cases} \tag{21}$$

$$\begin{cases} \frac{\mu\_{\rm hf}/\mu\_f}{\partial \mu\_{\rm hf}/\rho\_f} \frac{\partial^3 \zeta\_0}{\partial \eta^2} + \left( A\_2^1 \eta + f\_0 + \mathcal{g}\_0 \right) \frac{\partial^2 \zeta\_0}{\partial \eta^2} + \left( A - M + (1 - A\tau)\omega - 2 \frac{\partial \mathcal{g}\_0}{\partial \eta} \right) \frac{\partial \zeta\_0}{\partial \eta} \\ + (F + G) \frac{\partial^2 \mathcal{g}\_0}{\partial \eta^2} = 0, \end{cases} \tag{22}$$

$$\frac{1}{\Pr} \frac{k\_{\rm lmf}/k\_f}{(\rho \mathbb{C}\_p)\_{\rm lmf}/(\rho \mathbb{C}\_p)\_f} \frac{\partial^2 H}{\partial \eta^2} + (F + G) \frac{\partial \theta\_0}{\partial \eta} + \left( A \frac{1}{2} \eta + f\_0 + \mathfrak{g}\_0 \right) \frac{\partial H}{\partial \eta} + aH = 0,\tag{23}$$

subject to the boundary conditions:

 $F(0) = G(0) = 0$ ,  $\frac{\partial F}{\partial \eta}(0) = \frac{\partial G}{\partial \eta}(0) = 0$ ,  $H(0) = 0$ ,  $\frac{\partial F}{\partial \eta}(\infty) \to 0$ ,  $\frac{\partial G}{\partial \eta}(\infty) \to 0$ ,  $H(\infty) \to 0$ .

The stability of the steady-state flow and heat transfer solutions *f*0(η), *g*0(η) and <sup>θ</sup>0(η) was implemented by setting τ = 0 with *F* = *f*0(η), *G* = *g*0(η) and *H* = <sup>θ</sup>0(η) in (21)–(24). Finally, the initial development or the solution decay of the solution (20) is identified. The value of ω is obtained by solving the following eigenvalue problem:

$$\frac{\mu\_{\text{Im}f}/\mu\_f}{\rho\_{\text{Im}f}/\rho\_f} \mathbf{F}\_0^{\prime\prime} + \left(A\frac{1}{2}\eta + f\_0 + \mathfrak{g}\_0\right) \mathbf{F}\_0^{\prime\prime} + \left(A - M + \omega - 2f\_0^{\prime}\right) \mathbf{F}^{\prime} + (\mathcal{F} + \mathcal{G})f\_0^{\prime\prime} = 0,\tag{25}$$

$$\frac{\mu\_{\text{Im}f}/\mu\_f}{\rho\_{\text{Im}f}/\rho\_f} \mathcal{G}\_0''' + \left(A\frac{1}{2}\eta + f\_0 + g\_0\right) \mathcal{G}\_0'' + \left(A - M + \omega - 2g\_0'\right) \mathcal{G}' + (F + G)\mathcal{g}\_0'' = 0,\tag{26}$$

$$\frac{1}{\Pr} \frac{k\_{\rm Inf}/k\_f}{(\rho \mathbf{C}\_p)\_{\rm Inf}/(\rho \mathbf{C}\_p)\_f} H' + (\mathbf{F} + \mathbf{G}) \theta \mathbf{o}' + \left( A \frac{1}{2} \eta + f \mathbf{o} + \mathbf{g} \right) H' + aH = 0,\tag{27}$$

along with the conditions:

$$\begin{array}{l} F\_0(0) = G\_0(0) = 0, \; F\_0'(0) = G\_0'(0) = 0, \; H\_0(0) = 0, \\\ \; F\_0'(\infty) \to 0, \; G\_0'(\infty) \to 0, \; H\_0(\infty) \to 0. \end{array} \tag{28}$$

The range of possible eigenvalues can be determined by relaxing a boundary condition on *F* 0(η) according to the previous research by Harris et al. [56]. In this study, the condition *F* 0(∞) → 0 is relaxed, and the linear eigenvalue problem (25)–(28) is solved together with the new boundary condition *F* 0 (0) = 1 for a fixed value of ω1. The flow is considered unstable if the smallest eigenvalue ω1. is negative, which indicates an initial growth of disturbances occurred, while a positive value of the smallest eigenvalue signifies that the flow is physically achievable and stable.
