**3. HAM Solutions Methodology**

The homotopy analysis method (HAM) was applied to solve Equations (15)–(18). Shijun Liao developed this technique in 1992. It is often valid, regardless of whether there are a limited number of parameters or otherwise. It can be used to solve both weakly and strongly nonlinear problems. It offers a wide range of options for selecting the base functions of solutions, as well as discretion in choosing the linear operators. However, it provides a convenient method for ensuring the convergence of series solutions. Therefore, this method differs from other techniques, with examples like Adomain decomposition and the delta expansion methods. In the introduction section, some studies on the approach were presented.

Taking the initial guesses of the *f*(*ζ*),*θ*(*ζ*), *φ*(*ζ*), and *χ*(*ζ*) with the auxiliary linear operators respectively as

$$f\_0(\zeta) = 1 - e^{-\frac{\tau}{\zeta}},\ \theta\_0(\zeta) = \left(\frac{B\_i}{1 + B\_i}\right) e^{-\frac{\zeta}{\zeta}},\ \phi\_0(\zeta) = -\left(\frac{Nt}{Nb}\right) e^{-\frac{\zeta}{\zeta}},\ \chi\_0(\zeta) = e^{-\frac{\tau}{\zeta}}.\tag{24}$$

and

$$\mathcal{L}\_f = f^{\prime\prime\prime} - f^{\prime}, \quad \mathcal{L}\_\theta = \theta^{\prime\prime} - \theta, \quad \mathcal{L}\_\phi = \phi^{\prime\prime} - \phi, \quad \mathcal{L}\_\chi = \chi^{\prime\prime} - \chi. \tag{25}$$

the properties are satisfied as given below

$$\begin{aligned} \mathcal{L}\_f(\Lambda\_1 + \Lambda\_2 e^{\overline{\zeta}} + \Lambda\_3 e^{-\overline{\zeta}}) &= 0, & \mathcal{L}\_\theta(\Lambda\_4 e^{\overline{\zeta}} + \Lambda\_5 e^{-\overline{\zeta}}) &= 0, \\ \mathcal{L}\_\phi(\Lambda\_6 e^{\overline{\zeta}} + \Lambda\_7 e^{-\overline{\zeta}}) &= 0, & \mathcal{L}\_\chi(\Lambda\_8 e^{\overline{\zeta}} + \Lambda\_9 e^{-\overline{\zeta}}) &= 0 \end{aligned} \tag{26}$$

with arbitrary constants Λ*i*, *i* ∈ [1, 9]. The Zeroth order form of the problem is given by

$$(1-p)\mathcal{L}\_f[f(\mathbb{Q};p)-f\_0(\mathbb{Q})] = ph\_f \mathbf{N}\_f[f(\mathbb{Q},p), \theta(\mathbb{Q},p), \phi(\mathbb{Q},p), \chi(\mathbb{Q},p)],\tag{27}$$

$$|(1-p)\mathcal{L}\_{\theta}[\theta(\zeta;p)-\theta\_{0}(\zeta)]| = ph\mathsf{N}\_{\theta}[\theta(\zeta,p),f(\zeta,p),\phi(\zeta,p)],\tag{28}$$

$$(1 - p)\mathcal{L}\_{\Phi}[\phi(\zeta, p) - \phi\_0(\zeta)] = p h\_{\Phi} \mathbf{N}\_{\Phi}[\phi(\zeta, p), \theta(\zeta, p), f(\zeta, p)],\tag{29}$$

$$\mathbb{P}(1-p)\mathcal{L}\_{\mathbb{X}}[\chi(\mathbb{zeta},p)-\chi\_{0}(\mathbb{zeta})] = ph\_{\mathbb{X}}\mathbf{N}\_{\mathbb{X}}[\chi(\mathbb{zeta},pp),\Phi(\mathbb{zeta},p),f(\mathbb{zeta},p)],\tag{30}$$

with *p* ∈ [0, 1] as the embedded parameter, and nonlinear operators **N***<sup>f</sup>* , **N***θ*, **N***φ*, and **N***<sup>χ</sup>* obtained by using Equations (15)–(18).

The problems' equivalent *m* order of the deformation are

$$L\_f[f\_m(\zeta, p) - \eta\_m f\_{m-1}(\zeta)] = h\_f \mathcal{R}\_{f,m}(\zeta),\tag{31}$$

$$
\mathcal{L}\_{\theta}[\theta\_m(\zeta, p) - \eta\_m \theta\_{m-1}(\zeta)] = h\_{\theta} \mathcal{R}\_{\theta, m}(\zeta), \tag{32}
$$

$$
\mathcal{L}\_{\Phi}[\phi\_{\mathfrak{m}}(\zeta, p) - \eta\_{\mathfrak{m}}\phi\_{\mathfrak{m}-1}(\zeta)] = h\_{\Phi}\mathcal{R}\_{\Phi, \mathfrak{m}}(\zeta), \tag{33}
$$

$$
\mathcal{L}\_{\chi}|\chi\_{m}(\zeta,p) - \eta\_{m}\chi\_{m-1}(\zeta)| = h\_{\chi}\mathcal{R}\_{\chi\times m}(\zeta),\tag{34}
$$

$$\begin{cases} f\_m = S, f\_m' = 0, \theta\_m' - B\_1 \theta\_m = 0, Nb\theta\_m' + Nt\theta\_m' = 0, \chi\_m = 0, at\,\zeta = 0\\ f\_m' = 0, \theta\_m = 0, \phi\_m = 0, \chi\_m = 0 as\zeta \to \infty. \end{cases} \tag{35}$$

$$\eta\_{\mathcal{W}} = \begin{cases} 0, & \text{if } m \le 1 \\ 1, & \text{if } m > 1, \end{cases} \tag{36}$$

where R*<sup>m</sup> <sup>f</sup>* (*ζ*),R*<sup>m</sup> <sup>θ</sup>* (*ζ*),R*<sup>m</sup> <sup>φ</sup>* (*ζ*),R*<sup>m</sup> <sup>χ</sup>* (*ζ*) can be obtained using Equations (15)–(18). The general solutions are given by

$$f\_m(\zeta) = f\_m^\kappa(\zeta) + \Lambda\_1 + \Lambda\_2 e^{\overline{\zeta}} + \Lambda\_3 e^{-\overline{\zeta}},\tag{37}$$

$$
\theta\_m(\zeta) = \theta\_m^s(\zeta) + \Lambda\_4 e^{\zeta} + \Lambda\_5 e^{-\zeta},
\tag{38}
$$

$$
\phi\_m(\zeta) = \phi\_m^s(\zeta) + \Lambda \varsigma e^{\overline{\zeta}} + \Lambda \gamma e^{-\overline{\zeta}},
\tag{39}
$$

$$
\chi\_m(\zeta) = \chi\_m^s(\zeta) + \Lambda\_8 e^{\overline{\zeta}} + \Lambda\_9 e^{-\overline{\zeta}},
\tag{40}
$$

where (*f <sup>s</sup> <sup>m</sup>*(*ζ*), *θ<sup>s</sup> <sup>m</sup>*(*ζ*), *φ<sup>s</sup> <sup>m</sup>*(*ζ*), *χ<sup>s</sup> <sup>m</sup>*(*ζ*)) are special solutions.
