Countable families of solutions of a limit stationary semilinear fourth-order Cahn--Hilliard equation I. Mountain pass and Lusternik--Schnirel'man patterns in R^N

Solutions of the stationary semilinear Cahn--Hilliard equation -\Delta^2 u - u -\Delta(|u|^{p-1}u)=0 in R^N, with p>1, which are exponentially decaying at infinity, are studied. Using the Mounting Pass Lemma allows us the determination of two different solutions. On the other hand, the application of Lusternik--Schnirel'man (L--S) Category Theory shows the existence of, at least, a countable family of solutions. However, through numerical methods it is shown that the whole set of solutions, even in 1D, is much wider. This suggests that, actually, there exists, at least, a countable set of countable families of solutions, in which only the first one can be obtained by the L--S min-max approach.


Introduction and motivation for main problems
1.1. Models and preliminaries. This paper studies some multiplicity properties of the steady-states of the following fourth-order parabolic equation of the Cahn-Hilliard (C-H) type: The problem is completed with bounded smooth initial data, (1.2) u(x, 0) = u 0 (x) in R N .
Assuming that data u 0 (x) are sufficiently fast exponentially decaying at infinity, the same behaviour holds for the unique classic solution of (1.1), at least locally in time, since, for p > 1, u(x, t) may blow-up in finite time; see key references and results in [6] and [14]. The classic Cahn-Hilliard equation describes the dynamics of pattern formation in phase transition in alloys, glasses, and polymer solutions. When a binary solution is cooled sufficiently, phase separation may occur and then proceed in two ways: either nucleation, in which nuclei of the second phase appear randomly and grow, or, in the so-called spinodal decomposition, the whole solution appears to nucleate at once and then periodic or semi-periodic structures appear. Pattern formation resulting from phase transition has been observed in alloys, glasses, and polymer solutions. Parallelly, these types of equations possesses a great interest in biology and, after certain transformations (see below) as prototype for nerve conduction in the form of a FitzHugh-Nagumo system (cf. [17,25]).
Cahn-Hilliard equations-type such as (1.1) have been studied by many authors during the last decades. Among some of the works where several aspects related to these equations are [2,33,34,35,36]. Additionally, one can check the surveys in [14,20], where necessary aspects of global existence and blow-up of solutions for (1.1) are discussed in sufficient detail.
From the mathematical point of view, the non-stationary equation (1.1) involves a fourth order elliptic operator and it contains a negative viscosity term. In a more general (and applied) setting, the unknown function is a scalar u = u(x, t), x ∈ R N , t ∈ R + , and the equation reads This equation has been extensively studied in the past years but many questions still remain unanswered, especially in relation to multiplicity problems. Here, as customary, This yields the following C 1 -functional associated with (1.5): Solutions of the equation (1.5) are then obtained as critical points of the functional (1.6). The non-local operator (−∆) −1 is a positive linear integral operator from L q (R N ) to L 2 (R N ), within the Sobolev's range, i.e., 1 ≤ q ≤ 2N N −2 . Moreover, the operator (−∆) − 1 2 can be correctly defined as the square root of the operator (−∆) −1 and it will also be referred to as a non-local linear operator.
Also, since the problem is set in R N we are defining this operator in a class of exponentially decaying functions; see below the details and conditions to have the weak expression of the problem with these exponentially decay solutions.
A similar fourth-order problem was studied in [6], where further references can be found. However, in [6] the stationary equation (1.4) was considered in a bounded smooth domain Ω ∈ R N , with homogeneous Navier-type boundary conditions (1.7) u = ∆u = 0 on ∂Ω.
In particular, it was shown that this problem (1.4), (1.7) admits a countable family of solutions (critical points/values). Moreover, in [19, § 6] for a different variational fourth-order problem in R, it was shown that a wider countable set of countable families of solutions can be expected, where only the first infinite family is the Lusternik-Schnirel'man one. Hence, due to the analysis performed in [18,19] we need to study, in a general multidimensional geometry, existence and multiplicity for the elliptic equation (1.4) in a class of functions properly decaying at infinity (in fact exponentially), (1.8) lim |x|→∞ u(x) = 0.
Consequently, to get the results obtained here one can proceed following two different approaches. Firstly, bearing in mind solutions satisfying (1.8), with a fast decay at infinity, one can choose a sufficiently large radius R > 0 for the ball B R and consider the variational problem for the functional F(u) denoted by (1.6) in W 1,2 0 and assuming Dirichlet boundary conditions on S R = ∂B R . Thus, both spaces L 2 (B R ) and L p+1 (B R ) are compactly imbedded into W 1,2 0 (B R ) in the subcritical Sobolev's range (1.9) 1 < p ≤ p * = 1 + N +2 N −2 ≡ 2N N −2 = p S , N ≥ 3 (p * = +∞ for N = 1, 2). In other words, Note that, here, for the fourth-order elliptic operator in (1.4), p S = 2N N −2 , with N ≥ 3, because this operator has the representation so that, the necessary embedding features are governed by a standard second-order one ∆u + |u| p−1 u. The next step is passing to the limit as R → +∞, by using some uniform (in R) bounds on such families of solutions in B R . Since the category of the functional subset increases without bound in such a limit, eventually, we then expect to arrive at, at least, a countable family of various critical values/points. However, here we perform a variational study directly in R N , which was done previously for many fourth-order ODEs and elliptic equations; see [23,24,38] as key examples (though those equations, mainly, contain coercive operators, with "non-oscillatory" behaviour at infinity).
Namely, by a linearized analysis we first check that equation (1.4) provides us with a sufficient "amount" of exponentially decaying solutions at infinity. Obviously, in any bounded class of such functions, and in a natural functional viewing, since, loosely speaking, nothing happens at infinity (effectively, the solutions vanish there), the variational problem can be treated as the one in a bounded domain. So that the embedding in (1.10) comes in charge in some sense. Remark 1.1. Throughout this paper we shall consider radial solutions. Hence, we take into consideration a well known result about continuous Sobolev's embedding (see for instance [1]), which are compact replacing E by the radial subspace W 2,2 rad (R N ) and if in addition 2 ≤ N and p < p * (see [32]) where Hence, there is a constant S N > 0 such that . Note that for the second order case one has the continuous embedding (1.11) in the subcritical Sobolev's range (1.9).
Let us briefly summarize what we obtain here. First, we perform an analysis based on the application of the Mountain Pass Theorem in order to ascertain the existence of one solution and, additionally, the existence of more than one for the equation (1.4). To ascertain the existence of a second solution we use the auxiliary equation where u * represents a solution of the problem (1.5). Note that, since u * is a solution of (1.5) if u is a solution of (1.13) then, u + u * will be also a solution of (1.5). Thus, by construction and proving the existence of a solution for the equation (1.13), we find the existence of a second solution for (1.5) applying again a Mountain Pass argument. Indeed, performing this analysis and thanks to the equivalence between (1.4) and (1.5), we finally obtain the existence of at least a second solution for our main equation (1.4).
Thus, we state (proved in Section 3) the following.
Theorem 1.1. Suppose N ≥ 2 and radial solutions in W 1,2 rad (R N ) with exponential decay. Then, the non-local equation (1.5) possesses at least two solutions in W 1,2 rad (R N ). As far as we know, to get the existence of more than one solution, this seems to be the best available approach, since, for this type of higher-order PDEs we have a big lack of classical methodology and PDE theory.
In addition, secondly, we apply a L-S-fibering approach to get a countable family of solutions (critical points), though without any detailed information about how they look like.
Therefore we state, and prove in Section 4, the following.
Theorem 1.2. Suppose N ≥ 2 and radial solutions in W 1,2 rad (R N ) with exponential decay. Then, there is a countable family of solutions for the non-local equation (1.5) of the L-S type.
Finally, we apply advanced numerical methods to describe general "geometric" structure of various solutions assuming symmetric for even profiles and anti-symmetric conditions for odd profiles (see below). In particular, we introduce some chaotic patterns for equation (1.4) for p = 3 and p = 2, showing some profiles that become very chaotic away from the point of symmetry.
These numerical experiments suggest that the whole set of solutions, even in 1D, is much wider. However, these numerical experiments, together with some analytical approaches and estimates, will be extended and analysed with more detail in [4]. Indeed, it will be proved there that there exists, at least, a countable set of countable families of critical points, in which only the first one can be obtained by the L-S min-max approach.
Also, we observe that performing numerical experiments, shooting smoothly from x = 0 with u (0) = u (0) = u (0) = 0 and varying u(0), those chaotic patterns become more periodic when u(0) increases. Indeed, this fact appears to be sooner for p = 3 than for p = 2.
Moreover, it should be mentioned that this kind of transition behaviour is seen in similar phase solidification fourth-order equations, such as the Kuromoto-Sivashinsky and Swift-Hohenberg equations [37], [13], as a critical order parameter increases.
1.3. Previous related results. Recall again that, in [6], existence and multiplicity results were obtained for the steady-states of the unstable C-H equation of the form (1.14) −∆ 2 u + γu − ∆(|u| p−1 u) = 0 in Ω (p > 1), for a real parameter γ, where the multiplicity essentially depended on this parameter, which affected the category of the functional subset associated with the principle linear operator. The analysis was based on variational methods such as the fibering method, potential operator theory and Lusternik-Schnirel'man category-genus theory, and others, such as homotopy approaches or perturbation theory existence. Specifically, it was obtained that, depending on the parameter γ, there exists a different number of stationary solutions, i.e., • If the parameter γ ≤ Kλ 1 , with K > 0 a positive constant and λ 1 > 0 the first eigenvalue of the bi-harmonic operator, i.e., ∆ 2 ϕ 1 = λ 1 ϕ 1 with Navier boundary conditions (1.7), then there exists at least one solution for the equation (1.14), and; • When the parameter is greater than the first eigenvalue of the bi-harmonic equation λ 1 , multiplied by the positive constant K, then there will not be any solution at all, if one assumes only positive solutions. However, for oscillatory solutions of changing sign the number of possible solutions increases with the value of the parameter γ. In fact, when the parameter γ goes to infinity, one has an arbitrarily large number of distinct solutions.
Remark 1.2. Note that the previous distinction in terms of the number of solutions is related to pattern formation of problems such as the Kuramoto-Sivashinsky equation; see [22,37] for further references and details.
Furthermore, the so-called Mountain Pass Theorem [9] has been previously applied to find a second solution. Indeed, in [41] the authors analyzed the existence of a second solution for the fourth-order elliptic problem where δ ai represents the delta function supported at a i ∈ R N and Other methods utilized to solve similar problems to (1.4) such as where µ 1 and µ 2 are the squares of the roots of the characteristic polynomial and apply the Strong Maximum Principle having a positive solution if we assume that u is a homoclinic orbit (see [38,40] for any further details with N = 1). Moreover, in this case one can ascertain certain qualitative properties for that homoclinic orbit, such as the existence of a precisely critical point for u, a priori estimates and the symmetry of u with respect to that critical point. However, if 0 < β < √ 2, one cannot obtain such a decoupling. Also, similar fourth-order equations to (1.16) have been analyzed in [23,24] in the dimensional space R 4 and via Hamiltonian methods obtaining the connection between the critical points of the Hamiltonian.
On the other hand, note that problem (1.4) can be written as the following elliptic system: where v := (−∆) −1 u, which gives a different perspective to the problem in hand. The solutions of the system (1.18) represent the steady-states of a reaction-diffusion system of the form with a great interest in biology and as prototype for nerve conduction in the form of a FitzHugh-Nagumo system (cf. [17,25]). Mathematically, we note that system (1.18) is weakly coupled. However, this coupling is non-cooperative, in other words, the coupled terms have different signs. This means that the Maximum Principle is not valid for the system (1.18). For a cooperative system where the Maximum Principle can be applied to a very similar system to (1.18) see [3].
Basically, so far, one can see that, due to the lack of comparison methods, Maximum Principle, and several others classical methods in the analysis of higher-order PDEs. Therefore, in general the methodology is very limited and restricted to very specific examples.
1.4. Further extensions. Our results can be applied to other C-H models. For example note that, for the "true" fourth-order semilinear operator, the critical Sobolev's exponent is different: Obviously, a critical Sobolev's range as in (1.20) occurs for a different sixth-order C-H equation However, if the unstable nonlinear diffusion operator is of fourth order, as in we again arrive at the "second-order" Sobolev range as in (1.9). Equations (1.21) and (1.22) can be studied in similar lines, but some aspects become more technical, though not affecting the principal conclusions and results.

2.
Preliminary results: exponentially decaying patterns in R N 2.1. Exponentially decaying patterns in R N . The preliminary conclusions presented here formally allow us to consider our equations (1.4) in the whole R N , unlike as in [6] where the problem was assumed to be in a bounded domain Ω ⊂ R N .
Indeed, for the functional (1.6) we deal with the integrals over R N and, actually, with the functional setting over a certain weighted Sobolev space 1 , instead of W 2,2 0 (Ω) as assumed in [6]. Such a functional setting of the problem in R N is key in what follows. In fact, a proper functional setting assumes certain admissible asymptotic decay of solutions at infinity, which, for (1.4), is governed by the corresponding linearized operator.
Thus, considering (1.4) in the radial geometry, with u = u(r) and r = |x| ≥ 0, we then obtain Next, as usual, calculating the admissible decaying asymptotics from (2.1), using a two scale WKBJ-type asymptotics, as a first approximation (sufficient for our purposes), we use an exponentially pattern of the form u(r) ≈ r δ e ar (as r → ∞) in (2.1) leading easily to the following characteristic equation: . To be precise, note that the first equation in (2.2) comes from the homogeneity of the leading terms in (2.1). The second equality in (2.2) comes from a similar argument after evaluating the next leading terms on the left-hand side in (2.1). This yields a two-dimensional exponential bundle: where C 1,2 ∈ R are two arbitrary parameters of this linearized bundle. Through this asymptotic analysis we shall be able to show some patterns after performing a shooting problem in Section 5.
Let us now perform a preliminary analysis. In this radial geometry, any regular bounded solution of (1.4) must satisfy two boundary conditions at the origin (making the bi-Laplacian non-singular) Hence, using a standard shooting strategy from r = +∞, algebraically, at least two parameters are needed to satisfy both (2.4). Looking again at (2.3), where there exist two parameters C 1,2 ∈ R, we observe that matching with two symmetry boundary conditions (2.4) yields a well-posed and well-balanced algebraic "2D-2D shooting problem".
To justify existence of such solutions (critical points), we then need carefully apply different variational techniques. Of course, for the elliptic problem in R N , one needs a more technical and delicate "asymptotic separation" analysis. Namely, one needs to resolve the following "separation" procedure (as is well-known, for the bi-Laplacian, a standard separation of variables is not available): with R(r) obtained above, it then leads to a complicated equations on Y (σ). In general, such an asymptotic separation procedure is expected to determine a sufficiently wide infinite-dimensional asymptotic bundle of solutions with an exponential decay at infinity. However, we do not need such a full and rather technical analysis. We must admit that we still do not know that whether a countable family of L-S critical points are radially symmetric solutions or not. If the former is true, then the above radial analysis is sufficient. In general, using our previous experience, we expect that min-max critical points are not all radially symmetric, but cannot prove that.
Finally, we note that, obviously, this is not the case in 1D. Indeed, in addition to the symmetric (even) profiles satisfying (2.4), there exist others dipole-like/anti-symmetric (odd) solutions such that We can also check existence of such solutions numerically; see Section 5.
2.2. Spectral theory in R N . Now, we consider the linear spectral problem for the corresponding linearized non-local equation (1.5), i.e., which was analyzed, however, in a bounded domain, in full detailed in [17]. In problem (2.8) λ β represents the βth-eigenvalue associated with the eigenfunction ψ β . In the following we show and prove several properties of the spectrum of the linear eigenvalue problem (2.8). Subsequently, we will apply them, in particular, in order to get estimations of the category; see section 4.
As usual, we begin checking that the linearized equation (2.8) admits solutions with a proper exponential decay at infinity. Indeed, writing it down in the equivalent form we have that a proper exponential decay at infinity holds (here, as before, r = |x| 1, so we are again restricted to the radial setting of this linear spectral problem).
Moreover, due to the positivity of the operator L = ∆ 2 + Id on the left-hand side of (2.9), λ β = 0 is not an eigenvalue, so that λ β > 0, for any β. Secondly, owing to [11,Theorem V I.8] and the compactness of the inverse integral operator the spectrum might contain either infinitely many isolated real eigenvalues or a finite number of isolated eigenvalues, formed by a monotone sequence of eigenvalues In other words, for a sufficiently large α, the resolvent of the operator L + αI is compact and, hence, the spectrum is discrete and, of course, real, by symmetry (self-adjoint). Note as well that when the operator is self-adjoint the method shown in [21] can be used to prove that there are infinitely many eigenvalues. Thus, we have the following result: Proposition 2.1. The operator L admits a discrete set of eigenvalues that tend to +∞ and there exists at least a solution ψ ∈ W 1,2 rad (R N ). Remarks.
• The strict positivity of eigenvalues, for eigenfunctions with an exponential decay at infinity follows from the equality obtained by multiplying 2.9 by ψ β and integrating in L 2 . • Moreover, the corresponding associated family of eigenfunctions {ψ β } is a complete orthogonal set in L 2 . • According to our analysis above, to get exponentially decaying solutions, the real eigenvalues in (2.8) must satisfy (2.11) λ β > 0, for any β.
Note that problem (2.9) admits a positive first eigenvalue λ 1 characterized by the Raileigh quotient Moreover, we also find the following simple observation: If λ β > 0 is an eigenvalue of (2.9) for any β, then Proof. Indeed, from 2.10, integrating by parts and applying the Hölder's inequality yields , which proves the proposition. Furthermore, the natural space for the eigenfunctions of problem (2.9) is W 2,2 rad (R N ), i.e., the closure of C ∞ 0 -functions with respect to the norm is a reflexive Banach space. Under those assumptions we have the following variational expression of the problem (2.9): is an eigenfunction of the problem (2.9) associated with the eigenvalue λ β . Moreover, a weak formulation of the problem requires the following result: Then, there is a sequence {R n } ⊂ R + , with R n → ∞ as n → +∞, such that lim n→+∞ ∂B Rn u ∂u ∂n dS = 0, lim n→+∞ ∂B Rn ∇u ∂∇u ∂n dS = 0 and lim n→+∞ ∂B Rn u ∂∆u ∂n dS = 0, where B Rn is the ball of radius R n > 0, centered at the origin, and n denotes the unitary outward normal vector on ∂B Rn .
Note that, since, for λ β > 0 for any β, we always deal with solutions with exponential decay at infinity (see above), the previous assumptions are, hence, achievable.

Mountain Pass Theorem and existence of at least two solutions
In this section, we apply the celebrated Mountain Pass Theorem (cf. [9,39] for details of this highly cited theorem) to ascertain the existence of a solution and, consequently, of at least a second solution, for the problem (1.4) Recall that, in R N , with any N ≥ 2, we are restricted to a class of radially symmetric solutions, for which we know their exponential decay at infinity. For N = 1, we can deal with both even and odd patterns. However, in general, restrictions to both symmetries or not, makes no difference in the variational analysis. Nevertheless, in order to have the compact Sobolev's embedding (1.11) in the subcritical range (1.9) our results in this section are restricted to N ≥ 2.
) a radial function. However, by elliptic regularity, those weak solutions u are also strong solutions of the equation (1.5); see [10].
Then, we look for the existence of critical points for the functional F(u) (1.6). As previously proved in [6] this functional is weakly lower semicontinuous and its Fréchet derivative has the expression such that the directional derivative (Gateaux's derivative) of the functional (1.6) is the following d dt F(u + tϕ) |t=0 = D u F(u), ϕ = D u F(u)ϕ. Moreover, we know that the functional is C 1 and the critical points of the functional (1.6) denoted by C := {u ∈ W 1,2 rad (R N ) : D u F(u)ϕ = 0, for any ϕ ∈ W 1,2 rad (R N )} are weak solutions in H 1 rad (R N ) for the equation (3.10), i.e., D u F(u)ϕ = 0.
Thus, u ∈ C if and only if where F is called the gradient of F at u. Again, by classic elliptic regularity (Schauder's theory; see [10] for further details) we will then always obtain classical solutions for such equations.
The main ingredient we are applying in getting the existence of a solution is the celebrated Mountain Pass Theorem due to Ambrosetti-Rabinowitz [9]. Before stating the Mountain Pass Theorem, we show a very important and necessary condition in order to satisfy the theorem. This is the so-called Palais-Smale condition.

Palais-Smale condition (PS)
. Let E be a Banach space and {u n } ⊂ E a sequence such that (3.2) F(u n ) is bounded and F (u n ) → 0, as n → ∞, then {u n } is pre-compact, i.e., {u n } has a convergence subsequence. In particular, F(u n ) → c and F (u n ) → 0 ⇒ {u n } has a convergence subsequence.
The (PS) condition is a convenient way of imposing some kind of compactness into the functional F. Indeed, this (PS) condition implies that the set of critical points at the level value c is compact for any c ∈ R.
Remark 3.1. Since the functional F is C 1 it is easily proved that if there exists a minimizing sequence {u n } of solutions of the equation (1.5) weakly convergent in H 1 rad (R N ) to certain u 0 ∈ H 1 rad (R N ) and such that F (u n ) → 0 then we can assure that u 0 is a critical point, i.e., F (u 0 ) = 0.
Then, the functional F possesses a critical value c ≥ α > 0 characterized as Remark 3.2. As Ambrosetti-Rabinowitz mentioned, on a heuristic level, the theorem says that, if a pair of points in the graph of F are separated by a mountain range there must be a mountain pass containing a critical point between them. Also, although the statement of the theorem does not imply it, normally in the applications the origin u = 0 is a local minimum for the functional F. As we will see below that is our case. Most of the critical points will be maxima or minima. However, we cannot assure directly that those critical points are global maxima or minima, we shall need to work a bit harder to obtain that. Thus, we first show that for our problem the trivial solution u = 0 is a local minimum. Proof. Take a function g ∈ H 1 rad (R N ) normalized in the following way Then, taking a real number t sufficiently close to 0, using the expression of the functional F and applying the compact Sobolev's embedding (1.11), (1.9) we find that for a constant K > 0 that depends on the Sobolev's constant S N . Then, we arrive at with t sufficiently close to zero, i.e., t p−1 < p+1 2K .
As a first step we prove that the (PS) condition is satisfied by the functional F (1.6). Proof. Let {u n } be a sequence such that F(u n ) → c and F (u n ) → 0. Then, since F (u n ) → 0, for any ε > 0, there exists a subsequence of {u n } (denoted again by {u n }) such that . Furthermore, due to (3.2) we have that F(u n ) is bounded (or F(u n ) → c), i.e., |F(u n )| ≤ c, (3.5) or, 1 Hence, for a positive constant µ ∈ R (to be chosen below) and using (3.4), (3.5) µF(u n ) − D u F(u n ), u n = c + o(1) u n W 1,2 rad (R N ) = c + ε u n W 1,2 rad (R N ) . Indeed, we actually have that . Subsequently, for µ > 2 we arrive at the inequality

Now, analysing the function
and since G(0) = 0 one can easily realize that G(s) ≥ 0, for any s ≥ 0 if µ < p + 1.
Consequently, for any appropriate ε > 0 sufficiently small, and 2 < µ < p + 1 we get the boundedness of the norms in W 1,2 rad (R N ) for the elements of the subsequence {u n }. Note that p > 1.
Subsequently, due to the Hardy-Littlewood-Sobolev's inequality we find that with K a positive constant.
Thus, we can assure that there exists ρ, α > 0 such that Additionally, assuming the compact Sobolev's embedding (1.11) and satisfying the decay condition at infinity (1.8), we find that as t → ∞ and for a positive constant K proving the second condition of the lemma. Note that the positive constant K comes after using the Hardy-Littlewood-Sobolev inequality (3.6) with p * corresponding to the Sobolev range (1.9).
These results provide us with the existence of at least a solution for the non-local equation (1.5) and, hence, the existence of stationary solutions for the Cahn-Hilliard equation (1.4).

3.3.
Existence of at least a second solution. Now, we follow a similar argument line as performed above to get the existence of another solution for the non-local equation (1.5). Therefore, to ascertain the existence of at least a second solution, we will work on seeking for solutions of the equation such that u * is solution of the equation (1.5) and, hence, to the stationary Cahn-Hilliard equation (1.4). Note that, since u * is a solution of (1.5) if u = 0 is a solution of (3.10) then, u + u * will be also a solution of (1.5) different from u * As seen above the critical points of the functional F(u) (1.6) correspond to weak solutions of the equation (1.5). However, by the elliptic regularity, u * is also a strong solution of the equation (1.5). Now, for (3.10), we look for the existence of critical points for the functional denoted by Furthermore, note that this functional is weakly lower semicontinuous and its Fréchet derivative has the expression such that the directional derivative of the functional (3.11) is To arrive at the expression for those derivatives and equalities above for the functional J (u) (3.11) follow previous arguments shown for the functional F(u) (1.6) in [6]. For several additional properties we stress the work made in Sato-Watanabe [41]. In particular, Lemmas 5.1 and 5.2, where this argument to ascertain a second solution was used for a problem of the type (1.15).
Again, by classic elliptic regularity the solutions of the equation (3.10) are also classical solutions. Then, as previously mentioned, we apply again Mountain Pass Theorem to ascertain the existence of a second solution for the functional J denoted by (3.11).
Furthermore, note that, again, since the functional J is C 1 , if there exists a minimizing sequence {u n } of solutions of the equation (3.10) weakly convergent in H 1 rad (R N ) to certain u 0 ∈ H 1 rad (R N ), such that J (u n ) → 0 then we can assure that u 0 is a critical point J (u 0 ) = 0.
Also, since the nonlinearity of the equation (1.5) satisfies this implies that (3.10) possesses again the trivial solution u = 0. Indeed, it is not difficult to see that u = 0 is a local minimum just applying similarly Lemma 3.1 since J (0) = 0. Before proving the Mountain Pass Theorem for the perturbed problem (3.10) with the associated functional (3.11) we prove some properties for several auxiliary functions that will be used later. T (s) = (a + s) p s − a p s − 2+µ p+1 (a + s) p+1 + 2+µ p+1 a p+1 + (2 + µ)a p s + µ 2 pa p−1 s 2 , satisfies that T (s) ≥ 0, for any s ≥ 0 and µ = min{1, p − 1}.
Subsequently, we prove the (PS) condition for functional J (3.11).
Proof. Let {u n } be a sequence that satisfies the PS condition (3.2). Since J (u n ) → 0, for any ε > 0, there exists a subsequence of {u n } (denoted again by {u n }) such that . Additionally, due to the boundedness of the functional J (u n ) we have that Subsequently, for an appropriate positive constant µ ∈ R Now, we analyse the terms Indeed, thanks to Lemma 3.4 (i) we find that 20) and, hence, substituting (3.20) into (3.19) yields to Next, we apply the spectral theory of the non-local weighted eigenvalue problem with a weight of the form a(x) = p|u * | p−1 , and where λ * 1 is the first eigenvalue of the problem (3.21) associated with the eigenfunction ψ 1 . Then, it is easy to prove (see [17] for any further details about the spectral theory of certain classes of non-local operators of this form in bounded domains) that Note that, to extend the spectral theory developed in [17] for bounded domains, one needs to follow the analysis carried out above in Section 2 for these non-local problems, assuming exponential decay for the solutions.
Then, using the weighted eigenvalue problem (3.21) we find that Note that µ > 0 and due to Proposition 2.2 we can show that λ * 1 > 2 so that , and, hence, thanks to the Sobolev's embedding (1.11) (see Talenti [42] for best possible constants) and rearranging terms we arrive at , for a positive constant K. Therefore, we finally get the boundedness of the norms in W 1,2 rad (R N ) for the elements of the subsequence {u n }.
To conclude the prove we show the strong convergence of the subsequence {u n }. Indeed, applying the Hardy-Littlewood-Sobolev's inequality (3.6), inequality (3.17) and the compact Sobolev's embedding (1.11), for N ≥ 2 and 1 < p < p * for (1.9), we find that , for a positive constant K and, hence, proving the (PS) condition.
To conclude the proof of the conditions of the Mountain Pass Theorem 3.1, we prove the following result.
Additionally, assuming the compact Sobolev's embedding (1.11), as well as the Hardy-Littlewood-Sobolev's inequality (3.6) and satisfying the decay condition at infinity (1.8), we can see that, applying Hölder's inequality, Remark 3.5. Here we only obtain the existence of two solutions. Moreover, since those solutions could be oscillatory of changing sign we cannot actually compare them, nor knowing which ones they are in the subsequent sequence of solutions obtained via Lusternik-Schnirel'man's theory.

Towards to a first countable family of L-S critical points
In order to estimate the number of critical points of a functional, we shall need to apply Lusternik-Schnirel'man's (L-S) classic theory of calculus of variations. Thus, the number of critical points of the functional (1.6) will also depend on the category of a functional subset (see below some details).
There are very important applications of minimax methods to establish the existence of multiple critical points of functionals, which are invariant under a group of symmetries G, in the sense that F(hu) = F(u), ∀h ∈ G, and u ∈ W 1,2 (R N ) (a Banach space for our particular case). Note that, normally, the group of symmetries is G := {id, −id}.
As a simple case, suppose the functional F to be even, then F(u) = F(−u) for any u in the appropriate Banach space. One of the most famous methods to get these multiplicity results is due to Lusternik-Schnirel'man for symmetric functionals (see [39,Chapters 8,9,10] for further details and [6] for a discussion of this methodology for a similar problem to the one under consideration here). Basically, this topological theory for potential compact operators is a natural extension of the standard minimax principles which characterize the eigenvalues of linear compact self-adjoint operators. Applying the Calculus of Variations to the eigenvalue problem in an appropriate functional setting, one can see that the critical values of the functional involved are precisely the eigenvalues of the problem. Indeed, performing this characterization in the unit sphere S N −1 one has that the eigenvalues are the critical points of the functional associated to a linear operator L in the unit ball In our particular case, this functional subset is the following: in the spirit of the eigenvalues of linear operators in the unit ball (4.1). According to the L-S approach (see [10,28], etc.), in order to obtain the critical points of a functional on the corresponding functional subset, R 0 , one needs to estimate the category ρ of that functional subset. Thus, the category will provide us with the (minimal) number of critical points that belong to the subset R 0 . Namely, similar to [6,18], we, formally, may use a standard result: Lemma 4.1. The category of the manifold R 0 , denoted by ρ(R 0 ), is given by the number of eigenvalues (with multiplicities) of the corresponding linear eigenvalue problem satisfying: Proof. Let λ β be the β-eigenvalue of the linear bi-harmonic problem (4.4) such that taking into consideration the multiplicity of the β-eigenvalue, under the natural "normalizing" constraint k≥1 a k = 1.
Here, (4.5) represents the associated eigenfunctions to the eigenvalue λ β and Therefore, R 0 contains a sphere of an arbitrary bounded dimension. Hence, its category is then infinite.
One can say that the functional (1.6) is between at least two values, a maximum and a minimum one, Therefore, having at least two positive critical points for such a functional and since the L-S characterisation provides us with a lower bound for solutions, but not exactly how many, we should not ruled out the situation in which there are infinitely many critical points. Note that ρ(R 0 ) measures, at least, a lower bound of the total of number of L-S critical points. Moreover, thanks to the spectral theory shown above in Section 2 we have the sufficient spectral information about the eigenvalue problem (4.4) to obtain a sharp estimate of the category (4.3).
with β standing for the category of R 0 .
Hence, if A 1 contains a subset A 0 ∈ A β such that sup u∈A0 F(u) ≤ sup u∈A1 F(v) < c β+1 + ε, and which completes the proof.
Roughly speaking, since the dimension of the sets A belonging to the classes of sets A β increases with β such that A 1 ⊃ A 2 ⊃ . . . ⊃ A β , this guarantees that the critical points delivering critical values (4.8) are all different. Hence, to get those critical values we need to estimate the category ρ of that set R 0 .
Additionally, for this functional we show the following particular result (see [39,Chap. 9] for any further details) which provides us with a countable family of critical points for the functional F (1.6) following the spirit of the Mountain Pass Theorem. Remark 4.1. Thus, we find a countable family of critical points of the functional F defined by (4.8) such that c β = F(u β ) with u β a weak solution of the problem (1.5). However, we cannot assure how many exactly since so far we only acknowledge the existence of two solutions obtained through the Mountain Pass arguments performed in the previous section. Hence, suppose u β (in fact a critical point of (1.6)) is the function on which the subsequent infimum is achieved having at least β critical points, in other words γ(R 0 ) ≥ β. Now, let us take a two hump structure (as done in [18] with sufficiently large |a|, If necessary, we also perform a slight modification ofû(x) to have exponentially decay solutions at infinity. Thus, sinceû belongs to R 0 we have Observe that Thus, for anyû = M u β with u β ∈ R 0 , such that we have that and, hence, meaning that, in the present case, a two-hump structure cannot be a L-S β-solution. In particular for β = 1 or β = 2 we will have just one or two solutions. Actually the ones obtained in Section 3, which are the only ones we know there existence explicitly.

Numerical analysis in 1D and in the radial geometry
Considering radial geometry as discussed in section 2.1, (1.4) takes the form This is considered together with the far-field asymptotic behaviour (2.3), written in the form where we conveniently use the constants k 1,2 in place of C 1,2 , together with symmetry (2.4) or antisymmetry (2.7) conditions at the origin. Using the far-field behaviour (5.2), a 2-D shooting problem may be formulated where the parameters k 1,2 are determined by satisfying the symmetry condition (2.4) at the origin for the even profiles and antisymmetry condition (2.7) for the odd profiles. Matlab's ode15s solver is used with tight error tolerances (RelTol=AbsTol=10 −13 ).      We complete this numerical introduction with illustration of a different class of solutions to (5.1). We may perform numerical experiments, shooting smoothly from x = 0 with u (0) = u (0) = u (0) = 0 and varying u(0). Figures 6 and 7 are selected profiles in one-dimension (N = 1) in the illustrative parameter cases p = 3 and p = 2 respectively. For sufficiently small u(0) seemingly chaotic patterns are obtained, which emerge into a more periodic structure as u(0) increases. This emergence appears sooner for p = 3 than p = 2 when increasing the size of u(0). Such transition behaviour is seen in similar phase solidification fourth-order equations, such as the Kuromoto-Sivashinsky and Swift-Hohenberg equations [37], [13], as a critical order parameter increases.