A mathematical analysis of an extended MHD Darcy–Forchheimer type fluid

The presented analysis has the aim of introducing general properties of solutions to an Extended Darcy–Forchheimer flow. The Extended Darcy–Forchheimer set of equations are introduced based on mathematical principles. Firstly, the diffusion is formulated with a non-homogeneous operator, and is supported by the addition of a non-linear advection together with a non-uniform reaction term. The involved analysis is given in generalized Hilbert–Sobolev spaces to account for regularity, existence and uniqueness of solutions supported by the semi-group theory. Afterwards, oscillating patterns of Travelling wave solutions are analyzed inspired by a set of Lemmas focused on solutions instability. Based on this, the Geometric Perturbation Theory provides linearized flows for which the eigenvalues are provided in an homotopy representation, and hence, any exponential bundles of solutions by direct linear combination. In addition, a numerical exploration is developed to find exact Travelling waves profiles and to study zones where solutions are positive. It is shown that, in general, solutions are oscillating in the proximity of the null critical state. In addition, an inner region (inner as a contrast to an outer region where solutions oscillate) of positive solutions is shown to hold locally in time.

(1) V = v y, t , 0 , divV = 0 → ∂v ∂x = 0. (2) v(y, t) =v 0 (y). Based on the Darcy-Forchheimer model driven by Eq. (3), it is the intention to provide an extended model introduced under certain theoretical approaches. Previously, it is the intention to provide a historical context in what regards with the different analytical paths followed along this article.

OPEN
To start with, consider the seminal work firstly introduced by Kolmogorov, Petrovskii and Piskunov 17 , in combustion theory, and by Fisher 18 , in biology to predict genes interaction. These works introduced key ideas on how to face non-linear reaction-diffusion equations. As a basic principle, solutions were formulated under the scope of an incipient theory known as Travelling waves dynamics. The main question was focused on the conversion of a Partial Differential problem into an Ordinary Differential one, but it introduced a philosophical concern on the way solutions behave. Indeed, it was observed that depending on the selected Travelling wave speed, solutions may be oscillating or monotonic. Then, the question was mainly focus on selecting an appropriate travelling speed replicating the correct profile behaviour. Since the first formulation, the Fisher-KPP model has been applied to different ubiquitous applications: from ecology or biology (refer to [19][20][21] to non-linear fluids 22 . The proposed problem along this article can be named as an Extended Darcy-Forchheimer, in the sense that further theoretical terms are introduced on the classical formulation in Eq. (3). Firstly, an heterogeneous diffusion is introduced to account for potential heterogeneities in the proximity of the null critical state. Such heterogeneous diffusion is formulated with a high order operator. As an example on a process deriving into a fourth order operator, the reader can refer to 23 ) where the authors proposed an energy approach to understand the patterns formation in materials. In addition, high order operators have been of wide use in physics to control unstable patterns in the proximity of critical points given by bistable systems (refer to [24][25][26]. Mathematically, such operators are an intense area of research. As a representative example, the Giorgi's conjecture has been proved for an Allen-Cahn equation formulated with a high order operator (see 27 ). The fact of introducing heterogeneous diffusion shall be contemplated in accordance with the problem particularities to predict. Along this analysis, a fourth order operator has been admitted, nonetheless other possible diffusion types can be assumed depending of the reality to model (for instance, refer to the p-Laplacian approach in 28 together with the p-Laplacian non-Lipschitz 29 ). Finally and in relation with the advection term introduced, some interesting results to show well-posedness of solutions have been proved by Montaru in 30 .
Based on the last discussion, the proposed problem (P) is formed of three main terms: a heterogeneous diffusion in the form of a high order operator, a non-linear advection and a non-uniform reaction term of Darcy-Forchheimer type.
As discussed, the introduction of a high order operator aims the modelling of heterogeneous diffusion in a given media. Each type of fundamental non-linear diffusion has certain properties of potential interest in modelling. For instance, the Porous Medium Equation and the p-Laplacian driven diffusion have the property known as finite propagation along compact supports in finite mass functions. In this case, the high order diffusion induces a set of oscillating solutions whose properties permit to model diffusion in non-homogeneous media. Furthermore, the introduction of a non-linear advection has the objective of modelling specific flow movement induced by external means and encountered through the coefficients c, q. The non-uniform reaction depends on the spatial variable y through the term |y| µ and is postulated to model the environmental medium variation. Indeed, along the integration domain y, the media can change leading to a temporal change in the velocity (fluid local acceleration). Note that many functions can be provided instead of |y| µ , nonetheless the proposed one permits to account for radial topological variations typical of annular geometries.
The search of precise profiles of solutions is given in the travelling wave domain. This will permit to account for a system of ordinary differential equations, with certain further considerations of an appropriate travelling wave speed in combination with the other parameters introduced q, c, µ . Note that the travelling waves approach to non-linear reaction and diffusion models have been considered in different scenarios. The reader can refer to [31][32][33][34][35] . The mathematical analysis along this paper starts with the formulation of the problem as an abstract evolution supported by a strongly continuous semigroup. This permits to account for the analysis of regularity, existence and uniqueness of solutions. Afterwards, instability of solution profiles are studied in the Travelling wave domain. To this end, a set of Lemmas are introduced supported by the definition of generalized norms in a Hilbert-Sobolev space. Afterwards, a numerical exploration provides different types of solution profiles for different values in the involved model parameters c, q, µ . To make the numerical exploration possible, the solver bvp4c in MATLAB is used. Such solver is based on an implicit Runge-Kutta algorithm with interpolant extensions 36 together with a collocation method at the borders.
In addition, a mollifying norm is introduced. This is a relevant definition as the exponential mollifier permits to study possible patterns of extreme oscillating solutions. Definition 2 Given j ∈ R + and given 0 ≤ t < ∞ , the mollifying space H j is defined as per the norm: where the exponential kernel satisfies the A p -condition for (p = 1)37.

Regularity and bound of solutions. Firstly, consider the fundamental problem:
where L is a general operator in R N (including with no loss of results, N = 1): Then, the following bound properties hold:

Lemma 2.1
For v 0 ∈ L 2 (R) , the following bound holds: Admit v 0 ∈ L 2 (R) ∩ H j (R) , then: In addition, the following holds: Proof Given the fundamental equation (8), an abstract evolution is represented as: Performing the Fourier transform ( f − domain): Then: Now, admit the mollifying-weighted Sobolev space H j with the norm (Eq. 7): Scientific Reports | (2022) 12:5228 | https://doi.org/10.1038/s41598-022-08623-0 www.nature.com/scientificreports/ Note that as per initial assumption v 0 ∈ L 2 (R) , then: After standard arrangements: To conclude: for any given t ≥  38, p. 79 ). Then, it is concluded on the regularity of the derivatives up to the third order while the fourth order derivative shall be considered as a control parameter. The mollifying norm can be used as a bounding condition provided the fourth derivative is finite, leading to a finite ϒ.
We have proved that the space of mollifiers H j is bounded by the Lebesgue L 2 . In addition, the Sobolev space of oscillating solutions H 4 α has been shown bounded by the mollifying space H j upon a scaling ϒ .
The high order linear spatial operator − 2 (or in R , D 4 x ) can be regarded as the infinitesimal generator of a one parametric (being the parameter τ with 0 < τ ≤ t ) strongly continuous semi-group. Based on this, the following abstract form holds (note that the formulation is provided in R N for N = 1, 2, 3...): A solution to the fundamental v t = − 2 v with Dirac pulse as initial data v(y, 0) = δ(y) is given by: A kernel N(y, t) is determined making the inverse transformation: The last integral is convergent for f ∈ R . As a consequence, a kernel N(y, t) exists. Based on such kernel the following mapping is defined in the space H 4 α : such that for 0 < τ ≤ t , it holds: The coming theorem provides regularity conditions for the mapping � v 0 ,τ (v).

Theorem 2.1
The single parametric operator � v 0 ,τ (being τ the free parameter) is bounded in the space H 4 α (R) with the norm (Eq. 5).
Proof As a previous first step, the following inequality B 0 �v 0 � α ≤ �v� α (being B 0 > 0 a suitable constant) is shown:

Uniqueness.
To show uniqueness, the following theorem holds: Proof To show the proposed Theorem, consider: N and ∇N are bounded. Indeed, both are convergent functions in R (see Eq. (23)), then: for any value in the involved parameters s, ν.
Each of the involved integrals in Eq. (29) shall be assessed. Previously, the following supporting function is defined: For two fixed values of ǫ and s = T , the previous supporting function is bounded: In the proximity of null critical point any solution converges (this is further shown in the travelling waves analysis to come in the following Section), then it is possible to define M = max{|A + Bv 1 |, |A + Bv 2 |}: where the term �|y| µ � α < ∞ in the norm (Eq. 5). For any ball centered in 0 < τ ≤ t with radium given by any multiple of τ − s , uniqueness holds in the

Travelling waves
Travelling waves patterns of solutions are provided by the change v(y, t) = �(θ), θ = y − t ∈ R . Note that refers to the travelling wave velocity and � ∈ L 2 loc (R) ∩ L ∞ (R) . As well it may comply with (4) is formulated in the travelling wave domain as per the following expression: To control the non-uniform reaction, a bound truncation is defined to the term |y| µ as per: After the truncation defined, the following problem is considered (named P ′ ): Note that the truncated term |y| µ ǫ is bounded and it permits to define internal balls for |y| ≤ ǫ . The problem P ′ can be mentioned as the controlled problem, as the reaction evolves according to the controlling parameter ǫ . Consequently and based on the proposed problem P ′ , the following Lemma, to characterize the travelling wave motion, holds: The travelling wave stops for c =  www.nature.com/scientificreports/ As a consequence: As pointed, the integral is determined between −∞ and +∞ such that asymptotically, the following holds: Therefore: By standard procedure, it is easy to check that: ′ = − 1 2 . Consider now the advection term integration: Now, the integral on the left reads: so that: Eventually, the following integral is assessed: Compiling the different integrals assessed: Note that problem (4) is recovered for ǫ → ∞ . As ǫ > 0 , the Lemma postulations are easily checked by making: Idem for < and = in the previous expression.

Characterization of travelling waves oscillating behaviour. The pursuit intention consists on show-
ing and characterizing the instability behaviour of travelling waves profiles in the proximity of the null critical state = 0 . The proposed approach is inspired by a set of Lemmas and Theorems employed to characterize unstable patterns in Kuramoto-Sivashinsky equation 39 , in a Cahn-Hilliard equation 31 and in a high and odd order operator 32 . The proposal followed along our analysis varies compared to that in the cited references to account for the problem P ′ particularities.

Lemma 3.2 Given the Hilbert-Sobolev spaces H 4 , H 4
α together with the Lebesgue L 2 . The following bounds hold: Proof Admit the expression (7), then: as a consequence D 1 = 1 . The next inequality formulated is shown based on the norm (Eq. 5). Indeed: a.e. in R . It suffices to consider D 2 = 1 .
To show convergence principles in the oscillating travelling wave profile close the null equilibrium = 0 , the following perturbation function is introduced: κ(y, t) = v(y, t) − �(y, t) . The problem P (Eq. 4) is re-formulated in virtue of the new function κ(y, t) and the Travelling wave profile �(y, t) close the condition κ = 0.
It shall be considered that the involved terms during the perturbation analysis is given only in the derivatives, independently of the order, as any oscillating solution changes the sign of derivatives during evolution. Convergence is, hence, shown in the case of encompassed motion between the actual solution v and the travelling wave .
Admit the problem: Note that oscillating profiles are given due to the high order operator. To ensure this occurs, the term �(κ) is shown to be bounded in the following Lemma:

Proof
Then, the following holds: It shall be noted that the previous Lemma can be proved following a similar procesure for � : H α → L 2 , provided there exist ϕ > 1 , r > 0 and D 4 > 0 such that �(κ) L 2 < D 4 κ ϕ α , for 0 < κ ϕ < r . In this last case Note that � κ 0 ,τ has been proved as a bounded operator in Theorem 2.1. Furthermore, the mapping operator � κ 0 ,τ is the infinitesimal generator of a strongly continuous semi-group that admits an exponential representation to the linear operator L 0 . Based on this, the following Lemma enunciates: Lemma 3.4 The heterogeneous diffusion given by linear L 0 generates a strongly continuous semigroup abstractly represented by e tL 0 that satisfies: www.nature.com/scientificreports/ Proof First, and given the fundamental problem κ t = −κ (4) , the following exponential solution holds in the Fourier domain (f): Then: Then: A finite E 1 holds upon integration with t in (0, 1]. Proceeding in a similar way for E 2 : In addition, the following holds: A finite E 2 holds after integration with regards to t in (0, 1]. The next intention is to characterize the spectrum of L (as defined in expression (9) and for N = 1 ) close the null critical point. To this end, the following Lemma enunciates: Lemma 3.5 The local (local means close the equilibrium null state) spectrum of L (see expressions (8) and (9)) in H 4 has at least one eigenvalue (s) such that Re(s) > 0.
Proof The spectrum of L in the proximity of the null critical state may be shown making use of Evans functions (see 40 for a detail review of such theory and the relation with the characteristic polynomial of a matrix). In the presented analysis, the Eq. (38) is firstly converted into a matrix representation and an homotopy graph is provided close the null critical point to study the influence of the travelling wave speed.
The characteristic polynomial for the matrix is: Admit the approach � 1 < ǫ → 0 , then it is easy to check that R(s) has at least one root such that Re(s) > 0 . Indeed, given an arbitrary ǫ in the expression |y| µ ǫ , for 0 < s << 1 , R(s) < 0 while for s >> 1 , R(s) > 0 . The roots behaviour to R(s) can be guessed admitting |y| Geometric perturbation theory. Firstly, consider the following fundamental manifold defined in accordance with expression (38):   The manifold M ǫ shall be shown as a normal hyperbolic. To this end, the Fenichel invariant manifold Theorem 41 is considered as established in 42,43 . The eigenvalues of M ǫ in the proximity of the critical point, and in the transversal direction to the tangent space shall have non-zero real part. After computation, the eigenvalues for M ǫ are given as: The eigenvector to the zero eigenvalue is [ǫ, 0, 0, 0] that is tangent to M ǫ . Therefore, all the eigenvalues transversal to M ǫ have non-zero real part leading to state that M ǫ is a hyperbolic manifold. Now, the manifold M ǫ shall be proved to be locally invariant under the flows defined (particularly considering the flow for φ ′ 4 ). To this end and according to reference 43 it shall be shown that for all r > 0 , for all open interval J with a constant a ∈ J and for any given value of i ∈ N , there exists a δ such that for ǫ ∈ (0, δ) , the manifold M ǫ is kept invariant. To this end, admit i ≥ 1 and the following set of functions: A value for r > 0 can be selected in the consideration that M ǫ ∩ B r (0) is a sufficiently large interval so as to permit the evolution of the travelling wave profile along its domain. The assessment of δ is based on determining the distance between the flows for M ǫ and M 0 . To this end, consider that the involved functions in each flow are measurable (at least a.e.) in B r (0) with the norm (Eq. 5): In the null critical state, it can be considered the asymptotic approach φ 2 = 0 such that for a small φ 1 : For φ 1 < |ǫ| << 1 , admit the definition of δ(t) = cq�φ q−1 1 − ǫ q−1 � α + |θ + t| µ �A + Bφ 1 � α , for each t in B r (0) . Under these conditions, δ is finite and consequently M ǫ is close M 0 to preserve the normal hyperbolic condition.
Note that the linearized flows in the invariant manifold M ǫ permit to determine asymptotic travelling wave profiles in B r (0) by direct computation of a linear system of ODEs and given specific values in the involved parameters c, q and (see Figs. 5 and 6).

Construction of travelling wave profiles.
The existence and construction of travelling wave profiles can be made considering the common transversality between the manifolds M ǫ and M µ . For each manifold, a profile is obtained in the asymptotic condition and close the null critical point.
(72) 0, + cqǫ q−1 1 3 , + cqǫ q−1 1 3 e i2̺/3 , + cqǫ q−1 1 3 e i4̺/3 . www.nature.com/scientificreports/ As represented in Figs. 5 and 6 , it is possible to conclude on the existence of complex conjugate eigenvalues leading to oscillating bundles of solution. This suggests that the oscillating property is inherent to our problem and exists for any searched solutions. One of the main consequences of this last statement is the impossibility to find an appropriate travelling wave speed for which oscillations vanish in accordance with the Fisher-KPP philosophy introduced in 17 . Alternatively, the exploration of a specific speed for which the profile is positive can be translated into searching an appropriate speed for which the first travelling wave minimum is positive. Probably, this may happen only locally in time given the non-uniform reaction term.
The obtaining process of travelling waves profiles is performed with the solver bvp4c in the Sofware MATLAB. This solver is based on an implicit Runge-Kutta with interpolant extensions. In addition, the solver employs a collocation method to build the travelling profiles for particular conditions at the borders 36 . In our analysis, we consider a step like function as initial data � 0 (θ) = H(−θ) . This function enables us to study the evolution of a positive mass (at θ < 0 ) and a null one (for θ > 0 ). Furthermore, a heteroclinic orbit holds between the two states as � 0 (θ < 0) = 1 and � 0 (θ > 0) = 0 . To vanish the influence of the pseudo-boundary, required by the collocation procedure, the integration domain is admitted as large, i.e. θ ∈ (−1000, 1000).
Our main target is to explore the existence of a value in the scaled variable ̺ = t for which the first travelling minimum is positive. For the sake of simplicity and without loss of generality, the numerical scheme has been obtained for = 1 , introducing a wide interval of t values until the first minimum in the profile is positive. The numerical exploration suggests that in virtue of the non-uniform reaction, it is not possible to find a unique travelling wave speed (valid for any time) for which the first minimum is positive, and hence, ensuring a stationary inner region of positive solutions. Consequently, any positive region is explored only locally in time.
The numerical assessments have been pursuit for specific values in the involved parameters (see Figs. 7 and 8 ). In particular, c = 1 , q = 2 , A = B = 1 and µ = 2 . Note that for the selected values the convection coefficient complies with expression (39), then the profile moves from the left to the right.  www.nature.com/scientificreports/ Considering the first minimum positivity in Fig. 9, it is possible to obtain a spatial region of purely monotone solutions and for which a maximum principle holds locally in time. This spatial inner region is assessed based on the localization of θ = 1.88219 , the local time t = 0.235 and the travelling wave change θ = y − t , such that:

Conclusions
The proposed analysis has been focused on the mathematical treatment of an Extended Darcy-Forchheimer flow. Such extended problem has been introduced based on mathematical principles, having been such mathematical scope the main target of the presented study. To this end, we started by the definition of generalized functional spaces to characterize the oscillating behaviour of the high order operator. To this end, the H 4 α Hilbert-Sobolev space has permitted to account for the properties of solutions up to its fourth order derivative. In addition, the mollifier space H j has permitted to account for global mollification properties in the exponential mollifyingkernel. Afterwards, the problem has been formulated making use of an abstract evolution, such that existence and uniqueness of solutions are shown upon and based on the generalized Hilbert-Sobolev spaces defined. Solutions have been shown to be inherently unstable in virtue of the high order operator. To this end, a certain series of Lemmas (from Lemma 3.2 to 3.5) have been shown to apply in our case. Finally, a numerical exploration has been promoted in the search of travelling waves profiles and regions where positivity holds. In this sense and as a main finding, it has been observed that solutions are oscillating in the proximity of the null critical point (and hence being negative since the first minimum) except for a local time that has been precised assessed for specific values in the problem parameters ( c = 1 , q = 2 , A = B = 1 and µ = 2 ). It shall be noted that the numerical explorations have been provided for certain values in the involved parameters p, c, q, A and B but conclusions

Data availability
The author states that all data used is accessible.