\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{cite}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2015 (2015), No. 90, pp. 1--29.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2015 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2015/90\hfil Branching analysis of a countable family]
{Branching analysis of a countable family of global
similarity solutions of a fourth-order thin film equation}

\author[P. Alvarez-Caudevilla, V. A. Galaktionov \hfil EJDE-2015/90\hfilneg]
{Pablo Alvarez-Caudevilla, Victor A. Galaktionov}

\address{Pablo Alvarez-Caudevilla \newline
Universidad Carlos III de Madrid,
Av. Universidad 30, 28911-Legan\'es, Spain}
\email{pacaudev@math.uc3m.es Phone +34-916249099}

\address{Victor A. Galaktionov \newline
Department of Mathematical Sciences, University of Bath,
 Bath BA2 7AY, UK}
\email{vag@maths.bath.ac.uk Phone +44-1225826988}

\thanks{Submitted April 11, 2014. Published April 10, 2015.}
\subjclass[2000]{35K55, 35B32, 35G20, 35K41, 35K65}
\keywords{Thin film  equation;  Cauchy problem; source-type
similarity solutions;  \hfill\break\indent 
finite interfaces; oscillatory
sign-changing behaviour; Hermitian spectral theory; branching}

\begin{abstract}
 The main goal in this article is to justify that source-type and
 other global-in-time similarity solutions of the Cauchy problem
 for the fourth-order thin film equation
 $$ 
 u_t=-\nabla \cdot (|u|^n \nabla \Delta u) \quad
 \text{in $\mathbb{R}^N \times \mathbb{R}_ + $ where $n>0$,
 $N \ge 1$},
 $$
 can be obtained by a continuous deformation (a homotopy path) as
 $n \to 0^+$. This is done by reducing to similarity solutions 
(given by eigenfunctions  of a rescaled  linear operator $\mathbf{B}$) 
of the classic \emph{bi-harmonic  equation}
 $$ 
 u_t = - \Delta^2 u \quad\text{in $\mathbb{R}^N \times \mathbb{R}_ +$, where
 $\mathbf{B}=-\Delta^2 +\frac 14 y \cdot \nabla+ \frac N4 I$}.
 $$
 This approach leads to a countable family of various global similarity patterns
 of the thin film equation, and  describes their oscillatory sign-changing behaviour
 by using the known asymptotic properties  of the fundamental solution of
 bi-harmonic equation.
 The branching from  $n=0^+$ for thin film equation  requires Hermitian spectral
 theory for a  pair $\{\mathbf{B}, \mathbf{B}^*\}$ of non-self adjoint operators
 and leads to a number of difficult mathematical problems. These include, as a key
 part, the problem of multiplicity of solutions, which is under  particular scrutiny.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{remark}[theorem]{Remark}
\allowdisplaybreaks

\section{Introduction: TFEs, connections with classic PDE theory, layout}
 \label{S1}

\subsection{Main models, applications, and preliminaries}

    We study the global-in-time behaviour of  solutions of
the fourth-order quasilinear evolution equation of parabolic type,
called the \emph{thin film equation} (TFE-4),
\begin{equation}
\label{i1}
    u_{t} = -\nabla \cdot(|u|^{n} \nabla \Delta u)
 \quad \text{in } \mathbb{R}^N \times \mathbb{R}_ +\,,
\end{equation}
where $\nabla=\operatorname{grad}_x$, $\Delta=\nabla \cdot \nabla$ stands for the Laplace
operator in $\mathbb{R}^N$, and $n>0$ is a real parameter.
The TFE-4 \eqref{i1} is written for solutions of changing
sign, which can occur in the \emph{Cauchy problem} (CP) and
also in some \emph{free-boundary problems} (FBPs); see proper settings shortly.

Fourth- and sixth-order TFEs with a similar form to \eqref{i1} such as
(TFE-6)
\begin{equation}
\label{i1166}
    u_{t} = \nabla \cdot (|u|^{n} \nabla \Delta^2 u)\,,
\end{equation}
 as well as more complicated doubly nonlinear degenerate parabolic
 models  (see \cite{AEG1,Ki01} for other typical examples),
 have various applications in thin film,
lubrication theory, and in several other hydrodynamic-type
problems.

These equations model the
dynamics of a thin film of viscous fluid, such as the spreading of a
liquid film along a surface, where $u$ stands for the height of the
film (then clearly $u \ge 0$ that naturally leads to a FBP
setting).

Specifically, when $n=3$ we are dealing with a problem
in the context of lubrication theory for thin viscous films that
are driven by surface tension and when $n=1$ with Hele--Shaw
flows. We refer e.g. to \cite{EGK1, EGK3, Gia08, Grun04} for
recent  surveys and for extended lists of references
concerning physical derivations of various models, key
mathematical results and further applications. Moreover, since the 1980s
such equations also play quite a special role in nonlinear PDE
theory.

It is worth mentioning that \emph{nonnegative
solutions} with compact support  of various FBPs are mostly
physically relevant, and that the pioneering mathematical
approaches by Bernis and Friedman in 1990 \cite{BF1} were
developed mainly for such solutions in the one dimensional case.

Furthermore, some very important extensions to the N-dimensional case were
achieved by Beretta, Bertsch and Dal Passo in \cite{BBP}.


 However, \emph{solutions of changing sign} have already been under
 scrutiny for a few years
(see \cite{BW06, EGK2, EGK4}), which in particular, can have some
biological motivations as stated in personnal communications with
Professor J. R. King, to say nothing of general PDE theory.

It turned out that  these classes of the so-called
``oscillatory solutions of changing sign" of \eqref{i1} were
rather difficult to tackle rigorously by standard and classical
methods. Principally, due to the additional fact that any kind of detailed analysis
for higher-order equations is much more difficult than for
second-order counterparts (such as the notorious classic \emph{porous medium
equation} $u_{t} = \Delta (|u|^{n-1}u)$) in view of the lack of maximum
principle, comparison and order preserving semigroups and
potential properties of the operators involved.

Thus, practically
all the existing methods for monotone or variational operators
are not applicable to the TFE-4 \eqref{i1}.
Moreover, for TFEs such as \eqref{i1} even their self-similar radial (i.e. ODE)
representatives can lead to several surprises in trying to
describe sign-changing features close to interfaces; see
\cite{EGK2} for a collection of such hard properties.

On the other hand, it also turned
out that, for better understanding of such singular
oscillatory properties of solutions of the CP for \eqref{i1}, it
is  fruitful to consider the (homotopic) limit $n \to 0^+$
(see \cite{TFE4PV} for an extensive analysis of this homotopic approach)
 owing to Hermitian spectral theory developed in
\cite{EGKP} for a pair $\{\mathbf{B},\mathbf{B}^*\}$ of linear rescaled
operators  for $n=0$, i.e. for the \emph{bi-harmonic equation}
\begin{equation}
\label{s1}
 u_{t} = -\Delta^2 u\quad \text{in } \mathbb{R}^N \times \mathbb{R}_ +\,,
    \quad \mathbf{B}=-\Delta^2 +\frac 14 y \cdot \nabla+ \frac N4\, I, \quad
    \mathbf{B}^*=-\Delta^2 - \frac 14  y \cdot \nabla,
\end{equation}
which will be key for our further analysis and whose solutions are
$C^\infty$, have infinite speed of propagation
and oscillates infinitely near the interfaces.

\subsection{Main results}


In this article, our goal is, using this continuity/homotopy deformation
approach as ``$n \to 0^+$", to reduce the nonlinear eigenvalue problem
 (described in detail later on)
\begin{equation} \label{eigp}
   - \nabla \cdot (|f|^n \nabla \Delta f) + \beta y \cdot \nabla f +
 \alpha f=0 \quad\text{in }\mathbb{R}^N,
 \end{equation}
 to the linear eigenvalue problem
with the non-self-adjoint operator $\mathbf{B}$,
after writing the solutions of equation \eqref{i1} by
\begin{equation*}
 u(x,t):= t^{-\alpha} f(y), \quad \text{with }
    y:=\frac{x}{t^\beta},\quad \beta= \frac{1-n \alpha}4, \quad \text{for }  t>0,
\end{equation*}
with $f$ as the similarity profiles satisfying the nonlinear eigenvalue
 problem \eqref{eigp}.

Thus, we shall focus our
analysis on the Cauchy problem \eqref{i1} for exponents $n>0$,
which are assumed to be sufficiently small, studying large time behaviour
of the solutions of \eqref{i1},
i.e. \emph{source-type solutions} or \emph{global asymptotic behaviour}
(as $t \to \infty$)
with conservation of mass, in the case on a non-zero mass, i.e.,
\begin{equation*}
    \int u(x,t)\, \mathrm{d}x \neq 0.
\end{equation*}

Thus, we perform a systematic analysis
 of the behaviour of the similarity solutions
through a so-called \emph{homotopic approach} (branching from
$n=0$)
via branching theory, by using the Lyapunov--Schmidt methods,
obtaining relevant results and properties for the solutions of the
self-similar equation associated with \eqref{i1} and, hence, for
the proper solutions of \eqref{i1}.

Loosely speaking,
this approach is characterized as follows: good proper (similarity
or not) solutions of the Cauchy problem for the TFE \eqref{i1} are
those that can be continuously deformed (via a homotopic path) as
$n \to 0^+$ to the corresponding solutions of the {bi-harmonic
equation} \eqref{s1},
which will play a crucial role in the subsequent analysis.
Note that in \cite{PVcauchy} this analysis has been carried out
directly on the parabolic TFE-4 \eqref{i1}.

Then, under the previous conditions we are able to show the following result.

\begin{theorem} \label{thm1.1}
For sufficiently small $n>0$
\begin{equation}  \label{main}
 \text{there exists a countable set of
 solutions $\{f_k, \alpha_k\}_{|\sigma|=k \ge 0}$,}
\end{equation}
which depend directly on the dimension of
the eigenfunctions of the linear operator $\mathbf{B}$, where $\sigma$ is a multiindex
in $\mathbb{R}^N$ to numerate the pairs.
\end{theorem}

\begin{remark} \label{rmk1.1} \rm
Note that this continuity with respect to the parameter $n$
was already observed by Bernis, Hulshof \& Quir\'os in \cite{BHQ}.
For the one dimensional case and
assuming non-negative solutions for \eqref{i1} they ascertained that
the limit when $n \to 0^+$ cannot be the
CP due to the oscillatory properties of \eqref{s1}. However, since we are
supposing changing sign solutions for the the equation \eqref{i1}
our limiting problem might be the
CP \eqref{s1}.
\end{remark}

Furthermore, our homotopic-like approach
is based upon the spectral properties known for the linear
counterpart \eqref{s1} of the TFE \eqref{i1}. Moreover, owing to
the oscillatory character of the solutions of the bi-harmonic
equation \eqref{s1}, being
a ``limit case" of the TFE \eqref{i1}, close to the interfaces
this homotopy study exhibits a typical difficulty concerning the
desired structure of the transversal zeros of solutions, at least
for small $n>0$. The proof of such a transversality zero property is a
difficult open problem, though qualitatively, this was rather well
understood, \cite{EGK1}.


\begin{remark} \label{rmk1.2} \rm
It is worth mentioning that, unlike the FBPs, studied in hundreds
of papers since the 1980s (see \cite{Gia08} and \cite{EGK2} for key
references and  alternative versions of uniqueness approaches),
thin film theory for the Cauchy problem for \eqref{i1} or
\eqref{i2} has recently led to a number of difficult open problems and
is not still fully developed; see the above references as a guide
to main difficulties and ideas.

In fact, the concept of  proper solutions is still rather obscure for the
Cauchy problem, since any classic or standard notions of
weak-mild-generalized-\dots solutions fail in the CP setting.
Therefore, the complete definition  of the solutions of the CP for
\eqref{i1} is still unclear, although the results obtained in \cite{PVcauchy} provide
 us with a great improvement in that respect.
\end{remark}

\subsection{Previous related results, further extensions and layout}

Observe that this approach certainly has great similarities to a previous
one developed in the last two decades of the twentieth century for second-order
operators. It is well understood that for any $n>1$, (PME-2)
\begin{equation}
 \label{PME2}
    u_{t} = \Delta (|u|^{n-1}u)
\end{equation}
has a family of exact self-similar compactly
supported source-type solutions (the ZKB ones from 1950s), which
describes the large time behaviour of compactly supported solutions
with conservation of mass.


In addition, \eqref{PME2} also admits a countable
family of other similarity solutions; see \cite{GHCo} for key
references and most recent results.

The PME-2 \eqref{PME2} can be interpreted as a nonlinear degenerate
version of the classic \emph{heat equation} for $n=0$,
  $$
   u_t=\Delta u \quad \text{in }\mathbb{R}^N \times \mathbb{R}_+.
   $$
 Note that passing to the limit $n \to 0^+$ in \eqref{PME2}
for nonnegative solutions   was considered a difficult
 mathematical problem in the 1970s-80s, which exhibited typical
 (but clearly simpler than in the TFE case) features
 of a ``homotopy" transformation of PDEs.
 This study was initiated by Kalashnikov in 1978
 \cite{Kal78} for the one-dimensional case.


 Further detailed results in $\mathbb{R}^N$ were obtained in
 \cite{Ben81}; see also \cite{Cock99}. More recent
 estimates were obtained in \cite{Pan08, Pan07} for the 1D PME-2 \eqref{PME2}
 establishing the rate of convergence of solutions
 as $n \to 0^\pm$, such as $O(n)$ as $n \to 0^-$ (i.e. from $n<0$,
the fast diffusion range, where solutions are smoother)
  in $L^1(\mathbb{R})$ \cite{Pan08}, and $O(n^2)$ as $n \to 0^+$ in
 $L^2(\mathbb{R} \times (0,T))$ \cite{Pan07}.

 However, most of such convergence
 results are obtained for \emph{nonnegative} solutions of
 \eqref{PME2}.
 For solutions of changing sign, even for this second-order equation \eqref{PME2},
 there are some open problems; see \cite{GHCo} for
 references and further details.

Thus,  in the twenty-first century, higher-order TFEs such as
\eqref{i1}, though looking
  like a natural counterpart/extension of the PME-2
 \eqref{PME2}, corresponding mathematical
TFE theory is more complicated with several remaining problems still open.


Finally,  to summarize let us also mention that
 higher-order semilinear and quasilinear
parabolic equations occur in applications to thin film theory,
nonlinear diffusion, lubrication theory, flame and wave
propagation (the Kuramoto--Sivashinsky equation and the extended
Fisher--Kolmogorov equation), phase transition at critical
Lifshitz points and bi-stable systems (see Peletier-Troy
\cite{PT} for further details, models and results).

Furthermore, the analysis carried out in this paper could be extended to
more complicated models such as
the \emph{unstable fourth-order thin film equation} (the
unstable TFE-4):
\begin{equation}
\label{i2}
    u_{t} = -\nabla \cdot(|u|^{n} \nabla \Delta u)-\Delta(|u|^{p-1}u)\,,
\end{equation}
with the unstable homogeneous second-order diffusion term, where
$p>1$ is a fixed exponent; see \cite{EGK1} for physical
motivations, references,  and other basics. Here, \eqref{i2}
represents a fourth-order nonlinear parabolic equation with the
backward (unstable) diffusion term in the second-order operator. Note, that
the fourth-order term reflects surface tension effects and the
second-order term can reflect gravity, van der Waals interactions,
thermocapillarity effects, or geometry of the solid substrate.

Although, the analysis of equations such as \eqref{i2} will be the ultimate goal,
they are not within the scope of this paper.
Thus, the layout of the paper is as follows:
\begin{itemize}
\item[(I)] Problem setting and spectral theory for the linear equation \eqref{s1}.
Sections \ref{S2} and \ref{S3}

\item[(II)] Proof of the main results of the paper. Study of a countable
 family of global self-similar solutions of
\eqref{i1} via their branching from eigenspaces at $n=0^+$,
Section \ref{S4}.
\end{itemize}


\section{Problem setting and self-similar solutions}
\label{S2}

\subsection{The FBP and CP}

For  both the FBP and the CP, the solutions are assumed to satisfy
standard free-boundary conditions or boundary conditions at infinity:
\begin{equation}\label{i3}
\begin{gathered}
u=0, \quad \text{zero-height},\\
    \nabla u=0, \quad \text{zero contact angle},\\
    -\mathbf{n} \cdot \nabla (|u|^{n}  \Delta u)=0, \quad
    \text{conservation of mass (zero-flux)},
\end{gathered}
\end{equation}
at the singularity surface (interface) $\Gamma_0[u]$, which is the
lateral boundary of
\begin{equation*}
    \operatorname{supp}u \subset \mathbb{R}^N \times \mathbb{R}_ +,\quad N \geq 1\,,
\end{equation*}
where $\mathbf{n}$ stands for the unit outward normal to
$\Gamma_0[u]$.  Note that, for sufficiently smooth interfaces,
the condition on the flux can be read as
\begin{equation*}
    \lim_{\operatorname{dist}(x,\Gamma_0[u])\downarrow 0}
    -\mathbf{n} \cdot \nabla (|u|^{n}  \Delta u)=0.
\end{equation*}
For the FBP, dealing with \emph{nonnegative} solutions, this
setting is assumed to define a unique solution. However, this
uniqueness result is known in 1D only; see \cite{Gia08},
where the interface equation was included into the problem
setting. We also refer to \cite[\S~6.2]{EGK2}, where a ``local"
uniqueness is explained via \emph{von Mises} transformation, which
fixes the interface point. For more difficult, non-radial
geometries in $\mathbb{R}^N$, there is no hope of getting any uniqueness for
the FBP, in view of possible very complicated shapes of supports
leading to various ``self-focusing" singularities of interfaces at
some points, which can dramatically change the required regularity
of solutions.

For the CP, the assumption of non-negativity is got rid of, and
solutions become oscillatory close to interfaces. It is then key,
for the CP, that the solutions are expected to be
``smoother" at the interface than those for the FBP, i.e. \eqref{i3}
are not sufficient to define their regularity. These \emph{maximal
regularity} issues for the CP, leading to oscillatory solutions,
are under scrutiny in \cite{EGK2}. However, since as far as we know there is no
knowledge of how the solutions for these problems should be,
little more can be said about it.

Moreover, we denote by
\begin{equation*}
    M(t):=\int u(x,t) \, \mathrm{d}x,
\end{equation*}
the mass of the solution, where integration is performed over
smooth support ($\mathbb{R}^N$ is allowed for the CP only). Then,
differentiating $M(t)$ with respect to $t$ and applying the
divergence theorem (under natural regularity assumptions on
solutions and free boundary) we have that
\begin{equation*}
J(t):=  \frac{\mathrm{d}M}{\mathrm{d}t}
= -  \int_{\Gamma_0\cap\{t\}}\mathbf{n} \cdot \nabla
     (|u|^{n}  \Delta u )\, .
\end{equation*}
The mass is conserved if
   $ J(t) \equiv 0$, which is assured by the flux condition in
   \eqref{i3}.


The problem is completed with bounded, smooth, integrable,
compactly supported initial data
\begin{equation}\label{i4}
    u(x,0)=u_0(x) \quad \text{in } \Gamma_0[u] \cap \{t=0\}.
\end{equation}

In the CP for \eqref{i1} in $\mathbb{R}^N \times \mathbb{R}_ +$, one needs to pose
bounded compactly supported initial data \eqref{i4} prescribed in
$\mathbb{R}^N$. Then,  under the same zero flux condition at finite
interfaces (to be established separately), the mass is preserved.

\subsection{Global similarity solutions: towards a nonlinear eigenvalue problem}

 We now begin to specify the self-similar
solutions of the equation \eqref{i1}, which
are admitted due to its natural scaling-invariant  nature.
In the case of the mass being conserved, we have global in time
source-type solutions.


Using the following scaling in \eqref{i1}
\begin{gather*}
    x:= \mu \bar x,\quad t:= \lambda \bar t,\quad u:= \nu \bar u, \quad
    \text{with}, \\
 \frac{\partial u}{\partial t}= \frac{\nu}{\lambda} \frac{\partial \bar u}{\partial t},\quad
 \frac{\partial u}{\partial x_{i}}= \frac{\nu}{\mu} \frac{\partial \bar u}{\partial x_{i}},\quad
 \frac{\partial^2 u}{\partial x_{i}^2}= \frac{\nu}{\mu^2} \frac{\partial^2 \bar u}{\partial
    x_{i}^2},
\end{gather*}
 and substituting those expressions in \eqref{i1} yields
\begin{equation*}
    \frac{\nu}{\lambda} \frac{\partial \bar u}{\partial t}=-
    \frac{\nu^{n+1}}{\mu^4} \nabla \cdot (|\bar u|^{n} \nabla \Delta \bar u)\,.
\end{equation*}
To keep this equation invariant, the following  must be fulfilled:
\begin{equation}
\label{sf2}
    \frac{\nu}{\lambda}=\frac{\nu^{n+1}}{\mu^4},
\end{equation}
so that
\begin{equation*}
    \mu := \lambda^\beta \Rightarrow \nu := \lambda^{\frac{4\beta-1}{n}} \quad
    \text{and}\quad
    u(x,t):= \lambda^{\frac{4\beta-1}{n}} \bar u(\bar x,\bar t) = \lambda^{\frac{4\beta-1}{n}}
    \bar u(\frac{x}{\mu},\frac{t}{\lambda}).
\end{equation*}
Consequently,
\begin{equation*}
    u(x,t):= t^{\frac{4\beta-1}{n}} f(\frac{x}{t^\beta}),
\end{equation*}
where $t=\lambda$ and $f(x/t^\beta)=\bar u(x/t^\beta,1)$.
Owing to \eqref{sf2}, we obtain
\begin{equation*}
   n\alpha +4\beta=1,
\end{equation*}
which links the parameters $\alpha$ and $\beta$. Hence, substituting
\begin{equation}
\label{sf3}
    u(x,t):= t^{-\alpha} f(y), \quad \text{with}\quad
    y:=\frac{x}{t^\beta},\quad \beta= \frac{1-n \alpha}4,
\end{equation}
into \eqref{i1} and rearranging terms, we find that the function
$f$ solves a quasilinear elliptic equation of the form
\begin{equation}
\label{sf4}
    \nabla \cdot(|f|^{n} \nabla \Delta f)=\alpha f+\beta y \nabla \cdot f\,.
\end{equation}
 Finally,
 thanks to the above relation between $\alpha$ and $\beta$, we find a
\emph{nonlinear eigenvalue problem} of the form
\begin{equation} \label{sf5}
 -\nabla \cdot(|f|^{n} \nabla \Delta f) +\frac{1-\alpha n}{4} y \nabla \cdot f +\alpha
    f=0,     \quad f \in C_0(\mathbb{R}^N)\, ,
\end{equation}
 where we add to the equation \eqref{sf4} a natural assumption that
 $f$ must be compactly supported (and, of course, sufficiently
 smooth at the interface, which is an accompanying question to be
 discussed as well).

Thus, for such degenerate elliptic equations,
  the functional setting in \eqref{sf5} assumes that we are
 looking for  (weak) \emph{compactly supported} solutions $f(y)$ as
 certain ``nonlinear eigenfunctions" that hopefully occur for special values of nonlinear eigenvalues
  $\{\alpha_k\}_{k \ge 0}$. Our goal is to justify
  that, labelling  the eigenfunctions via a multiindex $\sigma$,
 \begin{equation}  \label{sf51}
 \parbox{10cm}{Equation \eqref{sf5} possesses  a countable set of
 eigenfunction/value pairs $\{f_k,\alpha_k\}_{|\sigma|=k \ge 0}$.}
  \end{equation}

 Concerning  the well-known properties of finite propagation for TFEs,
 we refer to papers  \cite{EGK1}--\cite{EGK4}, where a large amount of earlier
 references are available; see also \cite{GMPSobI} for more recent
 results and references in this elliptic area.

 However, one should observe that there are still
a few entirely rigorous results, especially those that are
attributed to the Cauchy problem for TFEs; for example \cite{PVcauchy}.

In the linear case $n=0$,
 the condition $f \in C_0(\mathbb{R}^N)$ is naturally replaced by the requirement that the
 eigenfunctions $\psi_\beta(y)$ exhibit typical exponential decay at
 infinity, a property that is reinforced by introducing  appropriate weighted
$L^2$-spaces.
Actually,  using the homotopy limit $n \to 0^+$, we will be obliged for
 small $n>0$,  instead of $C_0$-setting in
\eqref{sf5}, to use the following weighted $L^2$-space:
  \begin{equation}   \label{WW11}
  f \in L^2_\rho(\mathbb{R}^N), \quad \text{where }  \rho(y)={\mathrm e}^{a |y|^{4/3}},
  \; a>0 \text{ small }.
  \end{equation}

Note that, in the case of the Cauchy problem with conservation of
mass making use of the self-similar solutions \eqref{sf3}, we have
that
\[   M(t):=\int_{\mathbb{R}^N} u(x,t) \, \mathrm{d}x=t^{-\alpha} \int_{\mathbb{R}^N}
     f\big(\frac{x}{t^\beta}\big)
    \, \mathrm{d}x = t^{-\alpha+\beta N} \int_{\mathbb{R}^N} f(y)
    \, \mathrm{d}y,
\]
 where the actual integration is performed over the support
$\operatorname{supp} f$ of the nonlinear eigenfunction.
 Then, as is well known,
if $\int f \not = 0$,
  the exponents are calculated giving the first explicit nonlinear eigenvalue:
\begin{equation}\label{alb1}
-\alpha + \beta N=0 \Rightarrow
    \alpha_0(n)=\frac{N}{4+Nn} \quad \text{and}  \quad \beta_0(n)=\frac{1}{4+Nn}.
\end{equation}

\section{Hermitian spectral theory of the linear rescaled operators}
\label{S3}

 In this section, we establish the spectrum $\sigma(\mathbf{B})$
of the linear operator $\mathbf{B}$ obtained from the rescaling of the
linear counterpart of \eqref{i1}, i.e. the bi-harmonic equation
\eqref{s1}, which will be  essentially used in what follows.

\subsection{How the operator $\mathbf{B}$ appears: a linear eigenvalue problem}

Let $u(x,t)$ be the unique solution of the CP for the linear
parabolic bi-harmonic equation \eqref{s1} with the initial data
(the space as in \eqref{WW11} to be more properly introduced
shortly)
\begin{equation*}
    u_0 \in L_{\rho}^2(\mathbb{R}^N),
\end{equation*}
given by the convolution Poisson-type integral
\begin{equation}\label{s2}
    u(x,t)=b(t)* u_0 \equiv t^{-N/4} \int_{\mathbb{R}^N} F((x-z)t^{-1/4})
     u_0(z)\, \mathrm{d}z.
\end{equation}
Here, by scaling invariance of the problem, in a similar way as
was done in the previous section for \eqref{i1}, the unique
fundamental solution of the operator $\frac{\partial}{\partial t} + \Delta^2$ has
the self-similar structure
\begin{equation}
\label{s3}
   b(x,t)=t^{-N/4} F(y), \quad y:=\frac{x}{t^{1/4}} \quad  (x\in \mathbb{R}^N).
\end{equation}
Substituting $b(x,t)$ into \eqref{s1} yields that the rescaled
fundamental kernel $F$ in \eqref{s3} solves the linear elliptic
problem
\begin{equation}\label{s4}
  \mathbf{B}F \equiv -\Delta_y^2 F + \frac{1}{4} y \cdot \nabla_y F  +\frac{N}{4} F=0
    \quad \text{in} \quad \mathbb{R}^N,\quad \int_{\mathbb{R}^N} F(y) \, \mathrm{d}y=1.
\end{equation}
$\mathbf{B}$ is a non-symmetric linear operator, which is bounded
from $H_{\rho}^4(\mathbb{R}^N)$ to $L_{\rho}^2(\mathbb{R}^N)$ with the exponential
weight as in \eqref{WW11}.
Here,  $a\in (0,2d)$ is any positive constant, depending on the
parameter $d>0$, which characterises  the exponential decay of the
kernel $F(y)$:
 \begin{equation}
  \label{F11}
 |F(y)| \le D {\mathrm e}^{-d |y|^{4/3}} \quad \text{in} \quad
 \mathbb{R}^N \quad \big(D>0, \quad d=3 \cdot 2^{-11/3}\big).
   \end{equation}
 By  $F$ we denote the oscillatory rescaled kernel as
the only solution of \eqref{s4}, which has exponential decay,
oscillates as $|y|\to \infty$, and satisfies the  standard
pointwise estimate \eqref{F11}.

Thus, we need to solve the corresponding \emph{linear eigenvalue
problem}:
 \begin{equation}  \label{LP1}
  \mathbf{B} \psi = \lambda \psi \quad\text{in }\mathbb{R}^N, \quad \psi \in L^2_\rho(\mathbb{R}^N).
   \end{equation}
One can see that the nonlinear
   one \eqref{sf5} formally reduces to \eqref{LP1} at $n=0$ with
   the following shifting of the corresponding eigenvalues:
    $$
    \lambda=-\alpha + \frac N4.
   $$
 In fact, this is the main reason to calling \eqref{sf5} a
 nonlinear eigenvalue problem and, crucially, the discreteness of
 the real spectrum of the linear one \eqref{LP1} will be shown to
 be inherited by the nonlinear problem.

\subsection{Functional setting and semigroup expansion}

Thus,  we solve \eqref{LP1} and calculate the spectrum of $\sigma(\mathbf{B})$ in the
weighted space $L_{\rho}^2(\mathbb{R}^N)$. We then need the following
Hilbert space:
\begin{equation*}
    H_{\rho}^4(\mathbb{R}^N) \subset L_{\rho}^2(\mathbb{R}^N) \subset L^2(\mathbb{R}^N).
\end{equation*}
The Hilbert space $H_{\rho}^4(\mathbb{R}^N)$ has the following inner product:
\begin{equation*}
\langle v,w \rangle_{\rho} := \int_{\mathbb{R}^N} \rho(y) \sum_{k=0}^{4}
     D^{k} v(y) \overline{{D^{k} w(y)}}
     \, \mathrm{d}y,
\end{equation*}
where $D^{k} v$ stands for the vector
    $\{D^{\beta} v\,,\,|\beta|=k\}$,
and the norm
\begin{equation*}
   \| v \|_{\rho}^2 := \int_{\mathbb{R}^N} \rho(y) \sum_{k=0}^{4}
     |D^{k} v(y)|^2 \, \mathrm{d}y.
\end{equation*}
 Next, introducing the rescaled variables
\begin{equation}\label{s6}
    u(x,t)=t^{-N/4} w(y,\tau), \quad
 y:=\frac{x}{t^{1/4}}, \quad \tau= \ln t :
    \mathbb{R}_ + \to \mathbb{R},
\end{equation}
we find that the rescaled solution $w(y,\tau)$ satisfies the
evolution equation
\begin{equation}\label{s7}
    w_{\tau} = \mathbf{B}w\,,
\end{equation}
since, substituting the representation of $u(x,t)$ \eqref{s6} into
\eqref{s1} yields
\begin{equation*}
   -\Delta_y^2 w + \frac{1}{4}  y \cdot \nabla_y w  +\frac{N}{4}  w
= t \frac{\partial w}{\partial t} \frac{\partial \tau}{\partial t}.
\end{equation*}
Thus, to keep this invariant, the following should be satisfied:
\[
t  \frac{\partial \tau}{\partial t}=1  \Rightarrow \tau = \ln t,
\]
i.e.,  as defined in \eqref{s6}.
Hence, $w(y,\tau)$ is the solution of the Cauchy problem for the
equation \eqref{s7} and with the following initial condition at
$\tau=0$, i.e. at $t=1$:
\begin{equation}\label{s8}
    w_0(y) = u(y,1)\equiv b(1) *  u_0 = F *  u_0\, .
\end{equation}
Then, the linear operator $\frac{\partial}{\partial \tau} - \mathbf{B}$ is also a
rescaled version of the standard parabolic one
$\frac{\partial}{\partial t} + \Delta^2$. Therefore, the corresponding semigroup
${\mathrm e}^{\mathbf{B} \tau}$ admits an explicit integral representation.
This helps to establish some properties of the operator $\mathbf{B}$ and
describes other evolution features of the linear flow.

Indeed, from \eqref{s2} we find the following explicit representation of the
semigroup:
\[
w(y,\tau)=\int_{\mathbb{R}^N} F \big(y-z{\mathrm e}^{-\frac{\tau}{4}}\big)\, u_0(z) \,
    \mathrm{d}z \equiv {\mathrm e}^{\mathbf{B} \tau} w_0,
\]
where $x=t^{1/4}y$,  $\tau=\ln t$.
Subsequently, consider Taylor's power series of the analytic
kernel\footnote{We hope that returning here to the multiindex $\beta$
instead of $\sigma$ in \eqref{sf51} will not lead to a confusion with the
exponent $\beta$ in self-similar scaling \eqref{sf3}.}
\begin{equation}
\label{s10}
F\big(y-z {\mathrm e}^{-\frac{\tau}{4}}\big)=\sum_{(\beta)} {\mathrm e}^{
    -\frac{|\beta|\tau}{4}} \frac{(-1)^{|\beta|}}{\beta!}    D^\beta F(y) z^\beta
    \equiv \sum_{(\beta)} {\mathrm e}^{-\frac{|\beta|\tau}{4}}
\frac{1}{\sqrt{\beta!}} \psi_\beta(y) z^\beta,
\end{equation}
for any $y\in \mathbb{R}^N$, where
\begin{equation*}
    z^{\beta}:=z_1^{\beta_1}\dots z_{N}^{\beta_{N}},
\end{equation*}
and $\psi_{\beta}$ are the normalized eigenfunctions for the operator
$\mathbf{B}$.
The series in \eqref{s10} converges uniformly on compact subsets in
$z \in \mathbb{R}^N$. Indeed, denoting $|\beta|=l$ and estimating the
coefficients
\begin{equation*}
\big|\sum_{|\beta|=l} \frac{(-1)^{l}}{\beta!}  D^\beta F(y) z_1^{\beta_1}
\dots z_{N}^{\beta_{N}}\big| \leq b_l
    |z|^l,
\end{equation*}
by Stirling's formula we have that, for $l\gg 1$,
\begin{equation*}
   b_l = \frac{N^l}{l!} \sup_{y\in \mathbb{R}^N,
    |\beta|=l} |D^\beta F(y)| \approx \frac{N^l}{l!} l^{-l/4}
    {\mathrm e}^{l/4} \approx l^{-3l/4}c^l ={\mathrm e}^{-l \ln 3l/4 +l \ln c}.
\end{equation*}
Note that, the series
$\sum b_l |z|^l$,
has  radius of convergence  $R=\infty$.
 Thus, we obtain the
following representation of the solution:
\begin{equation*}
  w(y,\tau)= \sum_{(\beta)}  {\mathrm e}^{-\frac{|\beta|}{4}\, \tau} M_\beta(u_0) \psi_\beta(y),
    \quad \text{where }
    \lambda_\beta := -\frac{|\beta|}{4},
\end{equation*}
 and $\{\psi_\beta\}$ are the eigenvalues and
eigenfunctions of the operator $\mathbf{B}$, respectively, and
\begin{equation*}
  M_\beta(u_0):= \frac{1}{\sqrt{\beta!}} \int_{\mathbb{R}^N} z_1^{\beta_1}
    \dots z_{N}^{\beta_{N}} u_0(z) \, \mathrm{d}z,
\end{equation*}
are the corresponding momenta of the initial datum $w_0$ defined
by \eqref{s8}.

\subsection{Main spectral properties of the pair $\{\mathbf{B},\,\mathbf{B}^*\}$}

Thus,  the following holds \cite{EGKP}:


\begin{theorem}\label{thms1} \begin{itemize}
\item[(i)] The spectrum of $\mathbf{B}$ comprises real
eigenvalues only with the form
\begin{equation*}
 \sigma(\mathbf{B}):=\big\{\lambda_\beta = -\frac{|\beta|}{4},\,|\beta|=0,1,2,\dots\big\}.
\end{equation*}
Eigenvalues $\lambda_\beta$ have finite multiplicity with eigenfunctions,
\begin{equation} \label{s14}
\psi_\beta(y):= \frac{(-1)^{|\beta|}}{\sqrt{\beta!}} D^\beta F(y) \equiv
\frac{(-1)^{|\beta|}}{\sqrt{\beta!}}
    \big(\frac{\partial}{\partial y_1}\big)^{\beta_1}\dots \big(\frac{\partial}{\partial y_N}\big)^{\beta_N} F(y).
\end{equation}


\item[(ii)] The subset of eigenfunctions
  $  \Phi=\{\psi_\beta\}$
is complete in $L^2(\mathbb{R}^N)$ and in $L_{\rho}^2(\mathbb{R}^N)$.

\item[(iii)] For any $\lambda \notin \sigma(\mathbf{B})$, the resolvent
  $  (\mathbf{B}-\lambda I)^{-1}$
is a compact operator in $L_{\rho}^2(\mathbb{R}^N)$.
\end{itemize}
\end{theorem}

Subsequently, it was also shown in \cite{EGKP} that the adjoint
(in the dual metric of $L^2(\mathbb{R}^N)$) operator of $\mathbf{B}$ given by
\begin{equation*}
\mathbf{B}^*:=-\Delta^2 - \frac{1}{4}\, y \cdot \nabla,
\end{equation*}
 in the weighted space
$L_{{\rho}^*}^2(\mathbb{R}^N)$, with the exponentially decaying weight
function
\begin{equation*}
{\rho}^*(y) \equiv \frac{1}{\rho(y)} = {\mathrm e}^{-a|y|^\alpha} >0,
\end{equation*}
is a bounded linear operator,
$\mathbf{B}^*:H_{{\rho}^*}^4(\mathbb{R}^N) \to L_{{\rho}^*}^2(\mathbb{R}^N)$,
so
\[
    \langle \mathbf{B} v, w\rangle = \langle v, \mathbf{B}^* w\rangle,
  v \in H_{\rho}^4(\mathbb{R}^N), \,\,
    w \in H_{{\rho}^*}^4(\mathbb{R}^N).
\]
Moreover, the following theorem establishes the spectral properties of the
adjoint operator which will be very similar to those ones shown in Theorem
\ref{thms1} for the operator $\mathbf{B}$.

\begin{theorem} \label{thms2}
\begin{itemize}
\item[(i)] The spectrum of $\mathbf{B}^*$ consists of eigenvalues of
finite multiplicity,
\begin{equation*}
\sigma(\mathbf{B}^*)=\sigma(\mathbf{B}):=\big\{\lambda_\beta
= -\frac{|\beta|}{4}\,,\,|\beta|=0,1,2,\dots\big\},
\end{equation*}
and the eigenfunctions    $\psi_\beta^*(y)$
are polynomials of order $|\beta|$.


\item[(ii)]  The subset of eigenfunctions
$    \Phi^*=\{\psi_\beta^*\}$
is complete in $L_{{\rho}^*}^2(\mathbb{R}^N)$.


\item[(iii)] For any $\lambda \notin \sigma(\mathbf{B}^*)$, the resolvent
  $  (\mathbf{B}^*-\lambda I)^{-1}$
is a compact operator in $L_{{\rho}^*}^2(\mathbb{R}^N)$.
\end{itemize}
\end{theorem}


It should be pointed out that, since $\psi_0=F$ and
$\psi_0^* \equiv 1$, we have
\begin{equation*}
 \langle \psi_0, \psi_0^* \rangle= \int_{\mathbb{R}^N} \psi_0\, \mathrm{d}y
      =\int_{\mathbb{R}^N} F (y) \, \mathrm{d}y=1.
\end{equation*}
However, thanks to \eqref{s14}, we have that
\begin{equation*}
\int_{\mathbb{R}^N} \psi_\beta \equiv \langle \psi_\beta, \psi_0^* \rangle =0 \quad
 \text{for }  |\beta|\neq 0.
\end{equation*}
This expresses the orthogonality property to the adjoint
eigenfunctions in terms of the dual inner product.

Note that as shown in \cite{EGKP}, for the eigenfunctions $\{\psi_\beta\}$ of
$\mathbf{B}$ denoted by \eqref{s14}, the corresponding adjoint
eigenfunctions are \emph{generalized Hermite  polynomials} given by
\begin{equation}
\label{s16}
 \psi_\beta^*(y):=\frac{1}{\sqrt{\beta!}}\Big[y^\beta + \sum_{j=1}^{[|\beta|/4]}
    \frac{1}{j!}\, \Delta^{2j} y^\beta\Big].
\end{equation}
Hence, the orthonormality condition holds
\begin{equation*}
    \langle \psi_\beta,\psi_\gamma \rangle=\delta_{\beta\gamma}
    \quad \text{for any} \quad \beta,\;\gamma,
\end{equation*}
where $\langle \cdot,\cdot \rangle$ is the duality product
in $L^2(\mathbb{R}^N)$ and $\delta_{\beta\gamma}$ is  Kronecker's delta. Also,
operators $\mathbf{B}$ and $\mathbf{B}^*$ have zero Morse index (no
eigenvalues with positive real parts are available).

Moreover, the adjoint operator $\mathbf{B}^*$ was used in \cite{TFE4PV}
to analyse blow-up solutions for the TFE-4 \eqref{i1}.

Some key spectral results  can be extended as in \cite{EGKP} to
$2m$th-order linear poly-harmonic flows
 \begin{equation*}
 u_t= - (-\Delta)^m u \quad \text{in} \quad \mathbb{R}^N \times \mathbb{R}_ +,
  \end{equation*}
 where  the elliptic equation for the rescaled kernel $F(y)$ takes the
 form
\begin{equation*}
 \mathbf{B} F \equiv -(-\Delta_y)^m F + \frac{1}{2m} y \cdot \nabla_y F
 +\frac{N}{2m} F=0
    \quad \text{in} \quad \mathbb{R}^N,\quad \int_{\mathbb{R}^N} F(y) \,\mathrm{d}y=1.
\end{equation*}
In particular, for $m=1$, we find the \emph{Hermite operator} and
the \emph{Gaussian kernel} (see \cite{BS} for further information)
\begin{equation*}
  \mathbf{B} F \equiv \Delta F + \frac{1}{2} y\cdot \nabla F  +\frac{N}{2} \,F=0
     \;\Rightarrow\; F(y)= \frac 1{(4 \pi)^{N/2}} \, {\mathrm
     e}^{-|y|^2/4},
\end{equation*}
whose name is connected to fundamental works of Charles Hermite
on orthogonal polynomials $\{H_\beta\}$ about 1870.
 These classic Hermite polynomials are obtained by differentiating
 the Gaussian: up to normalization constants,
 \begin{equation} \label{Ga1}
 D^\beta{\mathrm
     e}^{-\frac{|y|^2}4}= H_\beta(y) \, {\mathrm
     e}^{-\frac{|y|^2}4} \quad \text{for any} \,\,\, \beta.
     \end{equation}
 Note that, for
$N=1$, such operators and polynomial eigenfunctions in 1D were
studied earlier by Jacques C.F.~Sturm in 1836; on this history and
Sturm's main original calculations, see \cite[Ch.~1]{GalGeom}.

The generating formula \eqref{Ga1} for (generalized) Hermite
polynomials is not available if $m \ge 2$, so that \eqref{s16} are
obtained via a different procedure, \cite{EGKP}.

\section{Similarity profiles for the Cauchy problem via  $n$-branching}
\label{S4}

 In general, construction of oscillatory similarity solutions of the
Cauchy problem for the TFE-4 \eqref{i1} is a  difficult nonlinear
problem, which is harder than for the corresponding  FBP one.

On the other hand, for $n=0$, such similarity profiles  exist and are
given by eigenfunctions $\{\psi_\beta\}$. In particular, the first
mass-preserving profile is just the rescaled kernel $F(y)$, so it
is unique, as was shown in Section \ref{S3}.

Hence, somehow, a possibility to visualize such an oscillatory
first ``nonlinear eigenfunction" $f(y)$ of changing sign, which
satisfies the \emph{nonlinear eigenvalue problem} \eqref{sf5},
  at least, for
sufficiently small $n > 0$ can be expected.

This suggests that, via an $n$-branching
approach argument, it is possible to ``connect" $f$ with the rescaled fundamental
profile $F$, satisfying the corresponding linear equation \eqref{s4},
with all the necessary properties of $F$  presented in Section
\ref{S3}.

 Thus, we plan to describe the behaviour of the similarity
 profiles $\{f_\beta\}$, as nonlinear eigenfunctions of \eqref{sf5}
for the TFE performing a ``homotopic" approach when $n\downarrow 0$.

 Homotopic approaches are well-known in the theory of
vector fields, degree and nonlinear operator theory (see
\cite{Deim, KZ} for details). However, we shall be less precise in
order to apply that approach and, here, a ``homotopic path" just
declares existence of a continuous connection (a curve) of
solutions $f \in C_0$ that ends up at $n=0^+$, at the linear
eigenfunction $\psi_0(y)=F(y)$ or further eigenfunctions
$\psi_\beta(y) \sim D^\beta F(y)$, as \eqref{s14} claims.

Using classical branching theory in the case of finite regularity
of nonlinear operators involved, we formally show that the
necessary orthogonality condition holds deriving the corresponding
\emph{Lyapunov--Schmidt branching equation}.
We will try to be as
rigorous as possible in supporting the delivery of the
nonlinear eigenvalues $\{\alpha_k\}$.

Further extensions of solutions
for non-small $n>0$ require a novel essentially non-local
technique of such nonlinear analysis, which remains an open
problem.

\subsection{Nonlinear eigenvalues $\{\alpha_k\}$ and transversality
conditions for the nonlinear eigenfunctions $f$}

In this first part of the section we establish the conditions and terms necessary
for the expansions of the parameter $\alpha$ and the nonlinear eigenfunctions,
as well as the transversality oscillatory conditions for such nonlinear eigenfunctions.

This will allow us to obtain the desired countable number of
solutions \eqref{main} for the similarity equation \eqref{eigp} via
Lyapunov-Schmidt reduction through the subsequent analysis.

The nonlinear eigenvalues $\{\alpha_k\}$ are obtained according to
non-self-adjoint spectral theory from Section \ref{S3}. We then
use  the explicit expressions for the eigenvalues and
eigenfunctions of the linear eigenvalue problem \eqref{LP1} given in
Theorem \ref{thms1}, where we also need the main conclusions of
the ``adjoint" Theorem \ref{thms2}.

 Thus, taking the corresponding linear equation from \eqref{sf5} with $n=0$,
we find, at least, formally, that
\begin{equation*}
 n=0: \quad  \mathcal{L}(\alpha)f:=-\Delta^2 f +\frac{1}{4} y \cdot\nabla f  +\alpha f=0.
\end{equation*}
Moreover, from that equation, combined with the eigenvalues expressions
obtained in the previous section, we ascertain the following
critical values for the parameter $\alpha_k=\alpha_k(n)$,
\begin{equation}
\label{bf4}
  n=0: \quad   \alpha_k(0) := -\lambda_k + \frac{N}{4} \equiv \frac{k+N}4
\quad \text{for } k=0,1,2,\ldots,
\end{equation}
where $\lambda_k$ are the eigenvalues  defined in Theorem \ref{thms1},
so that
\begin{equation*}
  \alpha_0(0)= \frac{N}{4}, \quad \alpha_1(0)= \frac{N+1}{4},\quad
 \alpha_2(0)= \frac{N+2}{4},\ldots ,   \alpha_k(0)= \frac{k+N}{4}\ldots \,.
\end{equation*}
In particular, when $k=0$, we have that $\alpha_0(0)= \frac{N}{4}$ and
the eigenfunction satisfies
\begin{equation*}
    \mathbf{B} F=0, \quad \text{so that }
    \ker \mathcal{L}(\alpha_0) = \operatorname{span}\{\psi_0\} \quad(\psi_0=F),
\end{equation*}
and, hence, since $\lambda_0=0$ is a simple eigenvalue for the operator
$\mathcal{L}(\alpha_0)= \mathbf{B}$, its algebraic multiplicity is 1. In
general, we find that
\begin{equation*}
\ker\big(\mathbf{B} + \frac{k}{4}\, I \big)
= \operatorname{span}\{\psi_\beta, |\beta|=k    \}, \quad \text{for any }
 k=0,1,2,3,\dots\,,
\end{equation*}
where the operator $\mathbf{B} + \frac{k}{4}I$ is Fredholm of
index zero since it is a compact perturbation of the identity of
linear type with respect to $k$. In other words,
$R[\mathcal{L}(\alpha_k)]$ is a closed subspace of $L_{\rho}^2(\mathbb{R}^N)$
and, for each $\alpha_k$,
\begin{equation*}
    \operatorname{dim} \ker(\mathcal{L}(\alpha_k))< \infty
\quad\text{and}\quad \operatorname{codim}R[\mathcal{L}(\alpha_k)]< \infty.
\end{equation*}

Then, for small $n>0$ in \eqref{sf5}, we can assume the following asymptotic
expansions
\begin{gather}\label{br3}
    \alpha_k(n):= \alpha_k+ \mu_{1,k} n+ o(n), \quad\text{and} \\
  \label{br3N}
     |f|^n \equiv  {\mathrm e}^{n\ln |f|}:= 1 +n \ln |f|+o(n).
\end{gather}
 As customary in bifurcation-branching theory \cite{KZ, VainbergTr},
existence of an expansion such as \eqref{br3} will allow one to get
further expansion coefficients in
\begin{equation*}
    \alpha_k(n):= \alpha_k+ \mu_{1,k} n + \mu_{2,k} n^2+ \mu_{3,k} n^3 + \dots\,,
 \end{equation*}
 as the regularity of nonlinearities allows and suggests, though the convergence
 of such an analytic series can be questionable and is not under
 scrutiny here.

Another principle question is that, for oscillatory sign changing
profiles $f(y)$, the last expansion \eqref{br3N} cannot be understood
in the pointwise sense. However, it can be naturally expected to be valid
in other metrics such as weighted $L^2$ or Sobolev spaces, as in
Section \ref{S3}, that used to be appropriate for the functional
setting of the equivalent integral equation and for that with
$n=0$.

Then, since \eqref{br3N} is obviously pointwise violated at the
nodal set $\{f=0\}$ of $f(y)$, this imposes some restrictions on
the behaviour   of corresponding eigenfunctions $\psi_\beta(y)$
($n=0$) close to their zero sets.
 Using well-known asymptotic and other related properties of the
\emph{radial} analytic rescaled kernel $F(y)$ of the fundamental
solutions \eqref{s3}, the generating formula of eigenfunctions
\eqref{s14} confirms that the nodal set of analytic eigenfunctions
$\{\psi_\beta=0\}$
 consists of isolated zero surfaces, which are ``transversal", at least in the a.e. sense,
with the only accumulation point at $y= \infty$.
 Overall, under such conditions, this indicates that
  \begin{equation}   \label{log1}
 \text{Expansion \eqref{br3N} contains not more than
 ``logarithmic" singularities a.e.},
  \end{equation}
which well suited the integral compact operators involved in
   the branching analysis,
  though we are far from claiming this as any rigorous issue.


 Moreover, when $n>0$ is not small enough, such an analogy and statements like
  \eqref{log1} become
unclear, and global extensions of continuous $n$-branches induced by
some compact integral operators, i.e. nonexistence of turning
(saddle-node) points in $n$, require, as usual, some unknown
monotonicity-like results.

Then,  to carry out our homotopic approach
we assume the expansion \eqref{br3N} away from possible zero
surfaces of $f(y)$, which, by transversality, can be localized in
arbitrarily small neighbourhoods.


Indeed, it is clear that when
$$
|f| > \delta >0, \quad \text{for any} \quad \delta>0,
$$
there is no problem in
approximating $|f|^n$ by \eqref{br3N}, i.e.,
$$
|f|^n =1+ O(n) \quad \text{as }  n \to 0^+.
$$
However, when
$$
|f| \leq \delta, \quad \text{for any } \delta>0,
$$
sufficiently small, the proof of such an approximation in weak
topology (as suffices for dealing with equivalent integral
equations) is far from clear unless
\begin{equation*}
 \text{the zeros of the $f$'s are also transversal a.e.},
  \end{equation*}
with a standard accumulating property at the
only interface zero surface. The latter issues have been studied
and described in \cite{EGK2} in the radial setting.
 Hence, we can suppose that such nonlinear eigenfunctions $f(y)$ are oscillatory
and infinitely sign changing close to the interface surface.


Therefore, if we assume that their zero surface is transversal
a.e. with a known geometric-like accumulation at the interface, we
find that, for any $n$ close to zero and any
$\delta= \delta(n) >0$ sufficiently small,
\begin{equation*}
    n| \ln |f| | \gg 1, \quad \text{if } |f| \leq \delta(n),
\end{equation*}
and, hence, on such subsets, $f(y)$ must be exponentially small:
\begin{equation*}
   | \ln |f| | \gg \frac{1}{n}\; \Rightarrow \;\ln |f| \ll -\frac{1}{n}\;
    \Rightarrow \; |f|  \ll {\mathrm e}^{-\frac{1}{n}}.
\end{equation*}
Recall that this happens in  also exponentially small
neighbourhoods of the transversal zero surfaces.

Overall, using the periodic structure of the oscillatory component
at the interface \cite{EGK2} (we must admit that such delicate
properties of oscillatory structures of solutions are known for
the 1D and radial cases only, though we expect that these
phenomena are generic), we can control the singular coefficients
in \eqref{br3N}, and, in particular, to see that
\begin{equation}
 \label{f1loc}
    \ln |f| \in L^1_{\rm loc} (\mathbb{R}^N).
\end{equation}
However, for most general geometric configurations of nonlinear
eigenfunctions $f(y)$, we do not have a proper proof of
\eqref{f1loc} or similar estimates, so our further analysis is
still essentially formal.


\subsection{Derivation of the branching equation}

Under the above-mentioned trans\-versality conditions and assuming the expansions
\eqref{br3}, for the nonlinear eigenvalues $\alpha_k$, and \eqref{br3N}, for the nonlinear
eigenfunctions $f$, we are able to obtain the branching equation applying the
classical Lyapunov-Schmidt method.

It is worth recalling again that our
computations below are to be understood as those dealing with the
equivalent integral equations and operators, so, in particular, we
can use the powerful facts on compactness of the resolvent
$(\mathbf{B}-\lambda I)^{-1}$ and of the adjoint one $(\mathbf{B}^*-\lambda I)^{-1}$ in
the corresponding weighted $L^2$-spaces. Note that, in such an
equivalent integral representation, the singular term in
\eqref{br3N} satisfying \eqref{f1loc} makes no principal
difficulty, so the  expansion  \eqref{br3N} makes rather usual
sense for applying standard nonlinear operator theory.

 Thus, under natural assumptions, substituting \eqref{br3} into \eqref{sf5}, for any
$k=0,1,2,3,\dots$\,, we find that, omitting $o(n)$ terms when
necessary,
\begin{equation*}
    -\nabla \cdot[(1 +n \ln |f|) \nabla \Delta f]
    +\frac{1-\alpha_k n - \mu_{1,k} n^2}{4} y \cdot\nabla  f
    +(\alpha_k+ \mu_{1,k} n) f=0\,,
\end{equation*}
and, rearranging terms,
\begin{equation*}
   -\Delta^2 f -n\nabla \cdot( \ln |f| \nabla \Delta f)
    +\frac{1}{4} y \cdot\nabla  f
    -\frac{\alpha_k n + \mu_{1,k} n^2}{4} y \cdot\nabla  f
    +\alpha_k f+ \mu_{1,k} n f=0\,.
\end{equation*}
Hence, we finally have
\begin{equation*}
   \big(\mathbf{B} + \frac{k}{4}\,I\big) f +n \big[-\nabla \cdot( \ln |f| \nabla \Delta f)
    -\frac{\alpha_k}{4} y \cdot\nabla  f
    + \mu_{1,k} f \big]+o(n)   =0\,,
\end{equation*}
which can be written in the  form
\begin{equation}\label{br5}
 \big(\mathbf{B} + \frac{k}{4}\,I \big) f +n \mathcal{N}_k (f)+o(n)   =0\,,
\end{equation}
with the operator
\begin{equation*}
  \mathcal{N}_k (f):=-\nabla \cdot( \ln |f| \nabla \Delta f)
    -\frac{\alpha_k}{4} y \cdot\nabla  f
    + \mu_{1,k} f\,.
\end{equation*}
Subsequently, as was shown in Section \ref{S3}, we have that
\begin{equation*}
\ker\big(\mathbf{B} + \frac{k}{4}\, I \big)
= \operatorname{span}\{\psi_\beta\}_{|\beta|=k}\quad \text{for  }
 k=0,1,2,3,\dots,
\end{equation*}
where the operator    $ \mathbf{B} + \frac{k}{4} I$
is Fredholm of index zero and
\begin{equation*}
   \operatorname{dim}\ker\big(\mathbf{B} + \frac{k}{4}\, I\big)
=  M_k  \ge 1 \quad \text{for any }k=0,1,2,3,\dots,
\end{equation*}
where $M_k$ stands for the length of the vector $\{D^\beta v, \,
|\beta|=k\}$, so that $M_k>1$ for  $k \ge 1$.

 \textbf{Simple eigenvalue for $k=0$.}
 Since 0 is a simple eigenvalue of $\mathbf{B}$ when $k=0$, i.e.,
\begin{equation*}
    \ker\mathbf{B} \oplus R[\mathbf{B}] = L_{\rho}^2(\mathbb{R}^N),
\end{equation*}
the study of the case $k=0$ seems to be simpler than for other
different $k$'s because the dimension of the eigenspace is
$M_0=1$.


Thus, we shall describe the behaviour of solutions for
small $n>0$ and apply the classical Lyapunov--Schmidt method to
\eqref{br5} (assuming, as usual, some extra necessary regularity
hypothesis to be clarified later on), in order to accomplish the
branching approach as $n\downarrow 0$, in two steps, when $k=0$
and $k$ is different from $0$.


Thus, owing to Section \ref{S3}, we already know that $0$ is
a simple eigenvalue of $\mathbf{B}$, i.e. $\ker\,\mathbf{B}=
\operatorname{span}\{\psi_0\}$ is one-dimensional. Hence, denoting by
$Y_0$ the complementary invariant subspace, orthogonal to
$\psi_0^*$, we set
\begin{equation*}
    f=\psi_0+V_0,
\end{equation*}
where $V_0 \in Y_0$.

Moreover, according to the spectral
properties of the operator $\mathbf{B}$, we define $P_0$ and $P_1$
such that $P_0+P_1=I$, to be the projections onto $\ker\,\mathbf{B}$
and $Y_0$ respectively. Finally, setting
\begin{equation}\label{br8}
 V_0:=n \Phi_{1,0} + o(n),
\end{equation}
substituting the expression for $f$ into \eqref{br5} and passing
to the limit as $n\to 0^+$ leads to a linear inhomogeneous
equation for $\Phi_{1,0}$,
\begin{equation}\label{br9}
    \mathbf{B}\Phi_{1,0}=- \mathcal{N}_0 (\psi_0),
\end{equation}
since $\mathbf{B} \psi_0=0$.

Furthermore, by  Fredholm theory, $V_0 \in
Y_0$ exists if and only if the right-hand side is orthogonal to
the one dimensional kernel of the adjoint operator $\mathbf{B}^*$
with $\psi_0^*=1$, because of \eqref{s16}. Hence, in the topology
of the dual space $L^2$, this requires the standard orthogonality
condition:
\begin{equation}\label{br10}
    \langle \mathcal{N}_0 (\psi_0), 1\rangle=0.
\end{equation}
Then, \eqref{br9} has a unique solution $\Phi_{1,0} \in Y_0$
determining by \eqref{br8} a bifurcation branch for small $n>0$.
In fact, the algebraic equation \eqref{br10} yields  the following
explicit expression for the coefficient $\mu_{1,0}$ of the
expansion \eqref{br3} for the first eigenvalue $\alpha_0(n)$:
\begin{equation*}
    \mu_{1,0} :=
   \frac{\langle \nabla \cdot( \ln |\psi_0| \nabla \Delta \psi_0)
    +\frac{N}{16}  y \cdot\nabla  \psi_0,\psi_0^*\rangle}
    {\langle \psi_0,\psi_0^*\rangle}
= \langle \nabla \cdot( \ln |\psi_0| \nabla \Delta \psi_0)
    +\frac{N}{16} y \cdot\nabla  \psi_0,\psi_0^*\rangle.
\end{equation*}
Consequently, in the particular case of having simple eigenvalues we just
obtain one branch of solutions emanating at $n=0$.


 \textbf{Multiple eigenvalues for $k \ge 1$.}
 Next we ascertain the number of branches
in the case when the eigenvalues of the operator $\mathbf{B}$ are semisimple.
For any
$k\geq 1$, we  know that
\begin{equation*}
\operatorname{dim}\ker\big(\mathbf{B} + \frac{k}{4} \,I\big)=  M_k >1.
\end{equation*}
  Hence, in order to perform a similar analysis to the one done for
simple eigenvalues   we have to use the full eigenspace expansion
\begin{equation}
\label{br12}
 f=\sum_{|\beta|=k} c_\beta \hat{\psi}_\beta +V_k,
\end{equation}
for every $k\geq 1$. Currently, for convenience, we denote
$$\{\hat{\psi}_\beta\}_{|\beta|=k}=\{\hat \psi_1,\dots,\hat \psi_{M_k}\},$$
the natural basis of the $M_k$-dimensional eigenspace
$\ker\big(\mathbf{B} + \frac{k}{4}\, I\big)$ and set
 $$
\psi_k = \sum_{|\beta|=k} c_\beta \hat{\psi}_\beta.
$$
 Moreover,
$$
V_k \in Y_k \quad \text{and} \quad V_k=\sum_{|\beta|>k} c_\beta {\psi}_\beta,
$$
where $Y_k$ is the complementary invariant subspace of
$\ker\big(\mathbf{B} + \frac{k}{4}\, I\big)$.

Furthermore, in the same way, as we did for
the case $k=0$, we define the $P_{0,k}$ and $P_{1,k}$, for every
$k\geq 1$, to be the projections of $\ker\big(\mathbf{B} +
\frac{k}{4}\, I\big)$ and $Y_k$ respectively. We also expand $V_k$
 as
\begin{equation} \label{br13}
    V_k:=n \Phi_{1,k} + o(n).
\end{equation}
Subsequently, substituting \eqref{br12} into \eqref{br5} and
passing to the limit as $n\downarrow 0^+$, we obtain the following
equation:
\begin{equation} \label{br14}
    \big(\mathbf{B}+ \frac{k}{4}\,I\big)\Phi_{1,k}
=- \mathcal{N}_k \big(\sum_{|\beta|=k} c_\beta {\psi}_\beta\big),
\end{equation}
under the natural ``normalizing" constraint
\begin{equation}\label{br15}
 \sum_{|\beta|=k} c_\beta=1 \quad (c_\beta \ge 0).
\end{equation}
Therefore, applying the Fredholm alternative, $V_k \in Y_k$ exists
if and only if the term on the right-hand side of \eqref{br14} is
orthogonal to $\ker\,\big(\mathbf{B}+ \frac{k}{4}\,I\big)$.
Then, multiplying the right-hand side of \eqref{br14} by $\psi_\beta^*$,
for every $|\beta|=k$,  in the topology of the dual space $L^2$, we
obtain an algebraic system of $M_k+1$ equations and the same
number of unknowns, $\{c_\beta, \, |\beta|=k\}$ and $\mu_{1,k}$:
\begin{equation}  \label{alg1}
\langle \mathcal{N}_k (\sum_{|\beta|=k} c_\beta {\psi}_\beta), \psi^*_\beta
\rangle=0 \quad \text{for all }  |\beta|=k,
 \end{equation}
 which is indeed the Lyapunov--Schmidt branching equation
 \cite{VainbergTr}.
In general, such algebraic systems are assumed to allow us to
obtain the branching parameters and, hence, establish the number of
different solutions induced on the given $M_k$-dimensional
eigenspace as the kernel of the operator involved.

However, we must admit and urge that the algebraic system
\eqref{alg1} is a truly difficult issue. One of the main features
of it is as follows:
 \begin{equation}  \label{alg2}
   \text{Equation \eqref{alg1} is not variational.}
\end{equation}
In other words, one cannot use for \eqref{alg1} the classic
category-genus theory of calculus of variation \cite{Berger, KZ},
to claim that the category of the kernel (equal to $M_k$) is the
least number of different critical points and hence of different
solutions.

To see \eqref{alg2}, it suffices to note that, due to \eqref{s14}
and \eqref{s16}, the generalized Hermite polynomials $\psi_\beta^*$
have nothing in common in the algebraic sense with the eigenfunctions
$\psi_\beta$ in the $L^2$-scalar products in \eqref{alg1}.


\subsection{A digression to Hermite classic self-adjoint theory}

 It is worth mentioning that for the classic second-order Hermite operator
 \begin{equation} \label{He1}
 \mathbf{B} = \Delta + \frac 12 y \cdot \nabla + \frac N2  I
 \quad \text{(then, in the $L^2$-metric, 
$\mathbf{B}^*= \Delta- \frac 12 y \cdot \nabla$)},
\end{equation}
statement \eqref{alg2} does not hold. Indeed, by classic theory
 \cite[p.~48]{BS}, these eigenfunctions are related to each other
 by
\begin{equation}   \label{He2}
  \psi_\beta(y) = D^\beta F(y) \equiv H_\beta(y)  F(y) ,\quad \text{where }
F(y)=(4 \pi)^{-\frac N2} {\mathrm e}^{-\frac{|y|^2}4}
   \end{equation}
  is the Gaussian kernel and $H_\beta(y)$ are standard Hermite
  polynomials, which also define the adjoint eigenfunctions:
   \begin{equation}    \label{He3}
   \psi_\beta^*(y)=b_\beta H_\beta(y) \equiv \frac{ b_\beta}{F(y)} \, \psi_\beta(y),
\end{equation}
    where $b_\beta$ are normalization constants.
One knows that this result comes from the symmetry of the operator
\eqref{He1} in the weighted metric of $L^2_\rho(\mathbb{R}^N)$, where
 $$
 \rho(y)= {\mathrm e}^{\frac{|y|^2}4} \sim \frac 1{F(y)}
 \;\Rightarrow\;
\mathbf{B}= \frac 1\rho\,  \nabla \cdot (\rho \nabla) + \frac N2 \, I,
  \quad \text{so } (\mathbf{B})_{L^2_\rho}^*=\mathbf{B}.
 $$
In view of the relations \eqref{He2} and \eqref{He3} of the
bi-orthonormal bases $\{\psi_\beta\}$ and $\{\psi_\beta^*\}$, the
corresponding algebraic systems such as \eqref{alg1} can be
variational. Moreover,  even the original nonlinear elliptic
equation, similar to \eqref{sf5}, where the 4th-order operator is
replaced by a natural 2nd-order, one of the porous medium type:
 $$
 -\nabla(|f|^n \nabla \Delta f) \mapsto \nabla(|f|^n \nabla f),
  $$
then, we are in a situation where it becomes variational.

Thus, in this case, both
  branching (local phenomena) and global extensions of
  $n$-bifurcation branches can be performed on the basis of the
   powerful Lusternik--Schnirel'man category variational theory
from 1920s \cite[\S~56]{KZ}, so that existence and multiplicity
(at least, not less than in the linear case $n=0$) of solutions
are guaranteed.

\subsection{Computations for branching of dipole solutions in 2D}

To avoid excessive computations and as a self-contained example,
we now ascertain some expressions for those coefficients in the
case when $|\beta|=1$, $N=2$, and $M_1=2$, so that, in our notations,
$\{\psi_\beta\}_{|\beta|=1}=\{\hat \psi_1, \hat \psi_2\}$.

Consequently,
in this case, we obtain the following algebraic system: the
expansion coefficients of $\psi_1=c_1 \hat \psi_1+c_2 \hat \psi_2$
satisfy
\begin{equation}\label{br16}
\begin{gathered}
    c_1  \langle \hat \psi_1^*,h_1 \rangle- \frac{c_1 \alpha_1}{4}\,
     \langle \hat \psi_1^*,y \cdot\nabla  \hat{\psi}_1 \rangle +c_1 \mu_{1,1}
    +c_2  \langle \hat \psi_1^*,h_2 \rangle- \frac{c_2 \alpha_1}{4}\,
     \langle \hat \psi_1^*,y \cdot\nabla  \hat{\psi}_2 \rangle = 0,
\\
    c_1  \langle \hat \psi_2^*,h_1 \rangle- \frac{c_1 \alpha_1}{4}\,
     \langle \hat \psi_2^*,y \cdot\nabla  \hat{\psi}_1 \rangle
    + c_2  \langle \hat \psi_2^*,h_2 \rangle- \frac{c_2 \alpha_1}{4}\,
     \langle \hat \psi_2^*,y \cdot\nabla  \hat{\psi}_2 \rangle+ c_2 \mu_{1,1}=0,
\\
    c_1+c_2=1,
    \end{gathered}
\end{equation}
where
\begin{equation*}
   h_1:= -\nabla \cdot  [ \ln (c_1 \hat{\psi}_1+c_2\hat{\psi}_2) \nabla \Delta
   \hat{\psi}_1 ],\quad
 h_2:= -\nabla \cdot
     [ \ln (c_1 \hat{\psi}_1+c_2\hat{\psi}_2) \nabla \Delta \hat{\psi}_2 ],
\end{equation*}
and, $c_1$, $c_2$, and $\mu_{1,1}$ are the coefficients that we
want to calculate, $\alpha_1$ is regarded as the value of the
parameter $\alpha$ denoted by \eqref{bf4} and dependent on the
eigenvalue $\lambda_1$, for which $\hat \psi_{1,2}$ are the associated
eigenfunctions, and $\hat \psi_{1,2}^*$ the corresponding adjoint
eigenfunctions. Hence, substituting the expression $c_2=1-c_1$
from the third equation into the other two, we have the following
nonlinear algebraic system
\begin{equation}\label{br17}
\begin{gathered}
    0=N_1(c_1,\mu_{1,1})- c_1\frac{\alpha_1}{4} \,\big[
     \langle \hat \psi_1^*,y \cdot\nabla  \hat{\psi}_1 \rangle
    -  \langle \hat \psi_1^*,y \cdot\nabla  \hat{\psi}_2 \rangle\big],
\\
    0=N_2(c_1,\mu_{1,1}) - c_1\frac{\alpha_1}{4} \big[
     \langle \hat \psi_2^*,y \cdot\nabla  \hat{\psi}_1 \rangle
    - \langle \hat \psi_2^*,y \cdot\nabla  \hat{\psi}_2 \rangle\big]+ \mu_{1,1},
    \end{gathered}
\end{equation}
where
\begin{gather*}
  N_1(c_1,\mu_{1,1}):= c_1  \langle \hat \psi_1^*,h_1 \rangle
   + \langle \hat \psi_1^*,h_2 \rangle
   - \frac{\alpha_1}{4}\,
     \langle \hat \psi_1^*,y \cdot\nabla  \hat{\psi}_2 \rangle
   -c_1  \langle \hat \psi_1^*,h_2 \rangle +c_1 \mu_{1,1},
   \\
N_2(c_1,\mu_{1,1}):=c_1  \langle \hat \psi_2^*,h_1 \rangle
   +  \langle \hat \psi_2^*,h_2 \rangle
   - \frac{\alpha_1}{4} \,\langle \hat \psi_2^*,y \cdot\nabla  \hat{\psi}_2 \rangle-
   c_1 \langle \hat \psi_2^*,h_2 \rangle-c_1 \mu_{1,1},
\end{gather*}
represent the nonlinear parts of the algebraic system, with $h_0$
and $h_1$  depending on $c_1$.

Subsequently, to guarantee existence of solutions of the system
\eqref{br16}, we apply the Brouwer fixed point theorem to
\eqref{br17} by supposing that the values $c_1$ and $\mu_{1,1}$
are the unknowns, in a disc sufficiently big
$D_R(\hat{c}_1,\hat{\mu}_{1,1})$ centered in a possible
non-degenerate zero $(\hat{c}_1,\hat{\mu}_{1,1})$. Thus, we write
the system \eqref{br17} in the matrix  form
\begin{equation*}
    \binom{0}{0}= \begin{pmatrix}
    -\frac{\alpha_1}{4}\,\big[
   \langle \hat \psi_1^*,y \cdot\nabla  \hat{\psi}_1 \rangle
  -  \langle \hat \psi_1^*,y \cdot\nabla  \hat{\psi}_2 \rangle\big] & 0\\
  -\frac{\alpha_1}{4} \,\big[ \langle \hat \psi_2^*,y \cdot\nabla  \hat{\psi}_1 \rangle
  - \langle \hat \psi_2^*,y \cdot\nabla  \hat{\psi}_2 \rangle\big]
    & 1\end{pmatrix}\binom{c_1}{\mu_{1,1}} + \binom{N_1(c_1,\mu_{1,1})}
    {N_2(c_1,\mu_{1,1})}.
\end{equation*}
Hence, we have that the zeros of the operator
\begin{equation*}
    \mathcal{F}(c_1,\mu_{1,1}):= \mathfrak{M}
\binom{c_1}{\mu_{1,1}} + \binom{N_1(c_1,\mu_{1,1})}
    {N_2(c_1,\mu_{1,1})},
\end{equation*}
are the possible solutions of \eqref{br17}, where $\mathfrak{M}$ is the
matrix corresponding to the linear part of the system, while
 $$
  (N_1(c_1,\mu_{1,1}),N_2(c_1,\mu_{1,1}))^T,
   $$
 corresponds to the
nonlinear part. The application $\mathcal{H}:\mathcal{A} \times
[0,1] \to \mathbb{R}$, defined by
\begin{equation*}
    \mathcal{H}(c_1,\mu_{1,1},t):= \mathfrak{M}
    \binom{c_1}{\mu_{1,1}} + t\binom{N_1(c_1,\mu_{1,1})}
    {N_2(c_1,\mu_{1,1})},
\end{equation*}
provides us with a homotopy transformation from 
$\mathcal{F}(c_1,\mu_{1,1})= \mathcal{H}(c_1,\mu_{1,1},1)$ to its
linearization
\begin{equation}\label{br18}
    \mathcal{H}(c_1,\mu_{1,1},0):= \mathfrak{M} \binom{c_1}{\mu_{1,1}}.
\end{equation}
Thus, the system \eqref{br17} possesses a nontrivial solution if
\eqref{br18} has a non-degenerate zero, in other words, if the next condition
is satisfied
\begin{equation}\label{br19}
     \langle \hat \psi_1^*,y \cdot\nabla  \hat{\psi}_1 \rangle-
     \langle \hat \psi_1^*,y \cdot\nabla  \hat{\psi}_2 \rangle\neq 0.
\end{equation}
Note that, if the substitution would have been $c_1=1-c_2$, the
condition might also be
\begin{equation*}
     \langle \hat \psi_2^*,y \cdot\nabla  \hat{\psi}_2 \rangle-
     \langle \hat \psi_2^*,y \cdot\nabla  \hat{\psi}_1 \rangle\neq 0.
\end{equation*}
Then, under condition \eqref{br19}, the system \eqref{br17} can be
written in the form
\begin{equation*}
    \binom{c_1-\hat{c}_1}{\mu_{1,1}-\hat{\mu}_{1,1}}=-\mathcal{M}^{-1}
    \binom{N_1(c_1,\mu_{1,1})-\hat{c}_1}
    {N_2(c_1,\mu_{1,1})- \hat{\mu}_{1,1}},
\end{equation*}
which can be interpreted as a fixed point equation. Moreover,
applying Brouwer's fixed point theorem, we have that
\begin{align*}
    \text{Ind}((\hat{c}_1,\hat{\mu}_{1,1}),\mathcal{H}(\cdot,\cdot,0))
& =    \mathcal{Q}_{C_R (\hat{c}_1,\hat{\mu}_{1,1})}(\mathcal{H}(\cdot,\cdot,0))\\
&=\deg (\mathcal{H}(\cdot,\cdot,0),D_R (\hat{c}_1,\hat{\mu}_{1,1}))\\
&= \deg (\mathcal{F}(c_1,\mu_{1,1}), D_R (\hat{c}_1,\hat{\mu}_{1,1})),
\end{align*}
where $\mathcal{Q}_{C_R
(\hat{c}_1,\hat{\mu}_{1,1})}(\mathcal{H}(\cdot,\cdot,0))$ defines
the number of rotations of the function
$\mathcal{H}(\cdot,\cdot,0)$ around the curve $C_R
(\hat{c}_1,\hat{\mu}_{1,1})$ and
$\deg (\mathcal{H}(\cdot,\cdot,0),D_R
(\hat{c}_1,\hat{\mu}_{1,1}))$ denotes the topological degree of
$\mathcal{H}(\cdot,\cdot,0)$ in $D_R (\hat{c}_1,\hat{\mu}_{1,1})$.
Owing to classical topological methods, both are equal.

Thus, once we have proved the existence of solutions, we achieve
some expressions for the coefficients required:
\begin{gather*}
\begin{aligned}
    \mu_{1,1} &=c_2 ( \langle \hat \psi_1^*+\hat \psi_2^*,h_1-h_2 \rangle
    - \frac{\alpha_1}{4}\,
     \langle \hat \psi_1^*+\psi_2^*
    ,y \cdot\nabla  \hat{\psi}_1-y \cdot\nabla  \hat{\psi}_2 \rangle)\\
   &\quad  -  \langle \hat \psi_1^*+\hat \psi_2^*,h_1 \rangle
    + \frac{\alpha_1}{4}\,
     \langle \hat \psi_1^*+\hat \psi_2^*
    ,y \cdot\nabla  \hat{\psi}_1  \rangle, \qquad\qquad\qquad\,\,
\end{aligned} \\
    c_1     =1-c_2.
\end{gather*}
The expressions for the coefficients in a general case might be
accomplished after some tedious calculations, otherwise similar to
those performed above.

Note that, in general, those nonlinear
finite-dimensional algebraic problems are rather complicated, and
the problem of an optimal estimate of the number of different
solutions remains open.

Moreover, reliable multiplicity results
are very difficult to obtain. We expect that this number should be
somehow related (and even sometimes coincides) with the dimension
of the corresponding eigenspace of the linear operators $\mathbf{B}+\frac{k}{4}\, I$, for any $k=0,1,2,\ldots\,$. This is a
conjecture only, and may be too illusive; see further supportive
analysis presented below.


 However, we devote the remainder of this section to a possible answer
to that conjecture, which is not totally complete though, since we
are imposing some conditions.

Thus, in order to detect the number of solutions of the nonlinear
algebraic system \eqref{br16}, we proceed to reduce this system to
a single equation for one of the unknowns. As a first step,
integrating by parts in the terms in which $h_1$ and $h_2$ are
involved and rearranging terms in the first two equations of the
system \eqref{br16}, we arrive at
\begin{gather*}
\begin{aligned}
&\int_{\mathbb{R}^N} \nabla \psi_1^* \cdot \ln (c_1 \hat{\psi}_1+c_2\hat{\psi}_2)
\nabla \Delta    (c_1 \hat{\psi}_1    + c_2 \hat{\psi}_2)    \\
& - c_1 \frac{\alpha_1}{4}
    \int_{\mathbb{R}^N} \hat \psi_1^* y \cdot\nabla  \hat{\psi}_1
    +c_1 \mu_{1,1}- c_2  \frac{\alpha_1}{4}
    \int_{\mathbb{R}^N} \hat\psi_1^* y \cdot\nabla  \hat{\psi}_2 = 0,
\end{aligned}    \\
\begin{aligned}
&\int_{\mathbb{R}^N} \nabla \hat \psi_2^*\cdot \ln (c_1 \hat{\psi}_1+c_2\hat{\psi}_2)
\nabla \Delta  (c_2 \hat{\psi}_1+ c_2 \hat{\psi}_2)\\
&  -c_1  \frac{\alpha_1}{4}  \int_{\mathbb{R}^N} \hat \psi_2^* y
\cdot\nabla  \hat{\psi}_1    + c_2 \mu_{1,1}
-c_2  \frac{\alpha_1}{4}  \int_{\mathbb{R}^N} \hat\psi_2^* y \cdot\nabla  \hat{\psi}_2 =0.
\end{aligned}
\end{gather*}
 By the third equation, we have that $c_1 = 1- c_2$, and   hence,
setting
$$
c_1 \hat{\psi}_1+c_2\hat{\psi}_2
= \hat{\psi}_1+(\hat{\psi}_2- \hat{\psi}_1)c_2,
$$
and substituting these into those new expressions for the first two equations of
the system, we find that
 \begin{equation} \label{br59}
\begin{aligned}
&\int_{\mathbb{R}^N} \nabla \hat \psi_1^* \cdot
   \ln (\hat{\psi}_1
+  (\hat{\psi}_2-\hat{\psi}_1)c_2) \nabla \Delta
   (\hat{\psi}_1+(\hat{\psi}_2- \hat{\psi}_1)c_2) +\mu_{1,1}
    -c_2 \mu_{1,1}  \\
&  - \frac{\alpha_1}{4}
    \int_{\mathbb{R}^N} \hat \psi_1^* y \cdot\nabla  \hat{\psi}_1
    + c_2  \frac{\alpha_1}{4}\int_{\mathbb{R}^N} \hat{\psi}_1^* y
    \cdot(\nabla  \hat{\psi}_1-\nabla \hat{\psi}_2)  = 0,
\\
 &\int_{\mathbb{R}^N} \nabla \hat \psi_2^*\cdot   \ln (\hat{\psi}_1+
 (\hat{\psi}_2-\hat{\psi}_1)c_2) \nabla \Delta
   (\hat{\psi}_1+(\hat{\psi}_2- \hat{\psi}_1)c_2)
    + c_2 \mu_{1,1}    \\
& - \frac{\alpha_1}{4}
    \int_{\mathbb{R}^N} \hat \psi_2^* y \cdot\nabla  \hat{\psi}_1
    +c_2  \frac{\alpha_1}{4}
    \int_{\mathbb{R}^N} \hat \psi_2^* y
    \cdot(\nabla  \hat{\psi}_1-\nabla \hat{\psi}_2) =0.
\end{aligned}
\end{equation}
Adding both equations, we have
\begin{align*}
   \mu_{1,1}
&  =    - \int_{\mathbb{R}^N}  (\nabla \hat \psi_1^*+ \nabla \hat \psi_2^*) \cdot
   \ln (\hat{\psi}_1+  (\hat{\psi}_2-\hat{\psi}_1)c_2) \nabla \Delta
   (\hat{\psi}_1+(\hat{\psi}_2- \hat{\psi}_1)c_2)\\
&\quad + \frac{\alpha_1}{4}
    \int_{\mathbb{R}^N} (\psi_1^*+ \psi_2^*) y
    \cdot \nabla  \hat{\psi}_1
    - c_2  \frac{\alpha_1}{4}
    \int_{\mathbb{R}^N} (\hat \psi_1^*+ \hat \psi_2^*) y
    \cdot(\nabla  \hat{\psi}_2-\nabla \hat{\psi}_1).
\end{align*}
Thus, substituting it into the second equation of \eqref{br59}, we
obtain the following equation with the single unknown $c_2$:
\begin{align*}
& -c_2^2 \frac{\alpha_1}{4}\int_{\mathbb{R}^N} (\hat \psi_1^*+ \hat \psi_2^*) y
    \cdot (\nabla  \hat{\psi}_2-\nabla \hat{\psi}_1) \\
&+  c_2 \frac{\alpha_1}{4}\Big(
    \int_{\mathbb{R}^N} (\hat\psi_1^*+ 2 \hat\psi_2^*) y
    \cdot \nabla  \hat{\psi}_1-
    \int_{\mathbb{R}^N} \hat\psi_2^* y \cdot \nabla \hat{\psi}_2\Big) \\
&  - \frac{\alpha_1}{4}
    \int_{\mathbb{R}^N} \hat\psi_2^* y \cdot\nabla  \hat{\psi}_1
    +\int_{\mathbb{R}^N} \nabla \psi_2^*\cdot
   \ln (\hat{\psi}_1+ (\hat{\psi}_2-\hat{\psi}_1)c_2) \nabla \Delta
   (\hat{\psi}_1+(\hat{\psi}_2- \hat{\psi}_1)c_2) \\
&  -c_2 \int_{\mathbb{R}^N} (\nabla \hat \psi_1^* +\nabla \hat\psi_2^*)
   \ln (\hat{\psi}_1+ (\hat{\psi}_2-\hat{\psi}_1)c_2) \nabla \Delta
   (\hat{\psi}_1+(\hat{\psi}_2- \hat{\psi}_1)c_2) =0,
\end{align*}
which can be written as
\begin{equation*}
    c_2^2 A + c_2 B + C + \omega(c_2) \equiv \mathfrak{F} (c_2) + \omega    (c_2)=0.
\end{equation*}
 Here, $\omega(c_2)$ can be considered as a perturbation of the quadratic
form $\mathfrak{F} (c_2)$ with the coefficients
defined by
\begin{gather*}
 A:=  -\frac{\alpha_1}{4}
    \int_{\mathbb{R}^N} (\hat\psi_1^*+ \hat\psi_2^*) y
    \cdot (\nabla  \hat{\psi}_2-\nabla \hat{\psi}_1), \\
  B:= \frac{\alpha_1}{4}\Big(\int_{\mathbb{R}^N} (\hat\psi_1^*+ 2 \hat\psi_2^*) y
    \cdot \nabla  \hat{\psi}_1- \int_{\mathbb{R}^N} \hat\psi_2^* y \cdot \nabla \hat{\psi}_2\Big),
\\
C:= - \frac{\alpha_1}{4}\int_{\mathbb{R}^N} \hat\psi_2^* y \cdot\nabla  \hat{\psi}_1,
\\
\begin{aligned}
\omega(c_2)&:=   \int_{\mathbb{R}^N} \nabla \hat\psi_2^*\cdot
   \ln (\hat{\psi}_1+ (\hat{\psi}_2-\hat{\psi}_1)c_2) \nabla \Delta
   (\hat{\psi}_1+(\hat{\psi}_2- \hat{\psi}_1)c_2)
   \\
&\quad -c_2 \int_{\mathbb{R}^N} (\nabla \hat\psi_1^*+\nabla \hat\psi_2^*)\cdot
   \ln (\hat{\psi}_1+ (\hat{\psi}_2-\hat{\psi}_1)c_2) \nabla \Delta
   (\hat{\psi}_1+(\hat{\psi}_2- \hat{\psi}_1)c_2).
\end{aligned}
\end{gather*}
Since, by the normalizing constraint \eqref{br15}, $c_2 \in
[0,1]$, solving the quadratic equation $ \mathfrak{F} (c_2)$ yields:
\begin{itemize}
\item[(i)] $c_2=0 \Rightarrow \mathfrak{F}(0) = C $;

\item[(ii)] $c_2 =1 \Rightarrow \mathfrak{F}(1)= A+B+C $; and

\item[(iii)] differentiating $\mathfrak{F}$ with respect to $c_2$, we obtain that
$\mathfrak{F}'(c_2) = 2 c_2 A+B$. Then, the critical point of the
function $\mathfrak{F}$ is $c_2^* = -\frac{B}{2A}$ and its image is $
\mathfrak{F}(c_2^*)=- \frac{B}{4A}+C$.
\end{itemize}

Consequently, the conditions to be imposed for having
more than one solution (we already know the existence of at least
one solution) are as follows:
\begin{itemize}
\item[(a)] $C  (A+B+C) >0$;
\item[(b)] $C \big(-\frac{B}{4A}+C\big)<0 $; and
\item[(c)] $0<-\frac{B}{2A}<1$.
\end{itemize}
Note that, for $-\frac{B}{4A}+C= 0$, we have just a single
solution.
Hence, considering the equation again in the form
 $$
 \mathfrak{F} (c_2) + \omega
(c_2)=0,
 $$
  where $\omega (c_2)$ is a perturbation of the quadratic form
$\mathfrak{F} (c_2)$, and bearing in mind that the objective is to
detect  the number of solutions of the system \eqref{br16}, we
need to control  somehow this perturbation.


Under  conditions
(a), (b), and (c), $\mathfrak{F} (c_2)$ possesses exactly two solutions.
Therefore, controlling the possible oscillations of the
perturbation $\omega (c_2)$ in such a way that
\begin{equation*}
    \| \omega (c_2) \|_{L^\infty} \leq \mathfrak{F}(c_2^*),
\end{equation*}
we can assure that the number of solutions for \eqref{br16} is
exactly two. This is  the dimension of the kernel of the operator
$\mathbf{B}+\frac{1}{4}\, I$ (as we expected in our more general
conjecture).



 The above particular example shows how
difficult the questions on existence and multiplicity of
solutions for such non-variational branching problems are.


 Recall that the actual values of the coefficients $A$, $B$, $C$,
 and others, for which the number of solutions crucially depends on, are very
 difficult to estimate, even numerically, in view of the complicated
 nature of the eigenfunctions \eqref{s14} involved. To say nothing of
 the nonlinear perturbation $\omega(c_2)$.

\subsection{Branching computations for $|\beta|=2$}


Overall, the above analysis provides us with some expressions for
the solutions for the self-similar equation \eqref{sf5} depending
on the value of $k$. Actually, we can achieve those expressions
for every critical value $\alpha_k$, but again the calculus gets
rather difficult.

For the sake of completeness, we now analyze the
case $|\beta|=2$ and $M_2=3$, so that $\{\psi_\beta\}_{|\beta|=2}=\{\hat
\psi_1, \hat \psi_2, \hat \psi_3\}$ stands for a basis of the
eigenspace $\ker\big(\mathbf{B} + \frac{1}{2}\, I \big)$, with $k=2$
($\lambda_k=-  k/4$).

Thus, in this case, performing in a similar way as was done for
\eqref{br16} with
$$\psi_2= c_1 \hat \psi_1+c_2 \hat \psi_2+c_3 \hat \psi_3,$$
we arrive at the following algebraic system:
\begin{equation} \label{br61}
\begin{gathered}
\begin{aligned}
 &c_1  \langle \hat \psi_1^*,h_1 \rangle
    +c_2  \langle \hat \psi_1^*,h_2 \rangle
    +c_3  \langle \hat \psi_1^*,h_3 \rangle - \frac{c_1 \alpha_2}{4}\,
     \langle \hat \psi_1^*,y \cdot\nabla  \hat{\psi}_1 \rangle\\
& - \frac{c_2 \alpha_2}{4}
  \langle \hat \psi_1^*,y \cdot\nabla  \hat{\psi}_2 \rangle
  - \frac{c_3 \alpha_2}{4}
     \langle \hat \psi_1^*,y \cdot\nabla  \hat{\psi}_3 \rangle +c_1 \mu_{1,2}= 0,
\end{aligned} \\
\begin{aligned}
&c_1  \langle \hat \psi_2^*,h_1 \rangle
    + c_2  \langle \hat \psi_2^*,h_2 \rangle
    + c_2  \langle \hat \psi_2^*,h_3 \rangle  - \frac{c_1
    \alpha_2}{4}\,
     \langle \hat \psi_2^*,y \cdot\nabla  \hat{\psi}_1 \rangle\\
 &  - \frac{c_2 \alpha_2}{4}\,
     \langle \hat \psi_2^*,y \cdot\nabla  \hat{\psi}_2 \rangle
- \frac{c_3 \alpha_2}{4}\,
     \langle \hat \psi_2^*,y \cdot\nabla  \hat{\psi}_3 \rangle+ c_2 \mu_{1,2}=0,
\end{aligned} \\
\begin{aligned}
&c_1  \langle \hat \psi_3^*,h_1 \rangle
    + c_2  \langle \hat \psi_3^*,h_2 \rangle
    + c_2  \langle \hat \psi_3^*,h_3 \rangle - \frac{c_1
    \alpha_2}{4}
     \langle \hat \psi_3^*,y \cdot\nabla  \hat{\psi}_1 \rangle\\
 &  - \frac{c_2 \alpha_2}{4}
     \langle \hat \psi_3^*,y \cdot\nabla  \hat{\psi}_2 \rangle
    - \frac{c_3 \alpha_2}{4}\,
     \langle \hat \psi_3^*,y \cdot\nabla  \hat{\psi}_3 \rangle+ c_3 \mu_{1,2}=0,
     \end{aligned} \\
    c_1+c_2+c_3=1,
    \end{gathered}
\end{equation}
where
  \begin{gather*}
 h_1:= -\nabla \cdot  [ \ln (c_1 \hat{\psi}_1+c_2\hat{\psi}_2
   +c_3 \hat{\psi}_3) \nabla \Delta    \hat{\psi}_1 ], \\
 h_2:= -\nabla \cdot
     [ \ln (c_1 \hat{\psi}_1+c_2\hat{\psi}_2
     +c_3 \hat{\psi}_3) \nabla \Delta \hat{\psi}_2 ], \\
 h_3:= -\nabla \cdot  [ \ln (c_1 \hat{\psi}_1+c_2\hat{\psi}_2
   +c_3 \hat{\psi}_3) \nabla \Delta
   \hat{\psi}_3 ],
\end{gather*}
and $c_1$, $c_2$, $c_3$, and $\mu_{1,2}$ are the unknowns to be
evaluated. Moreover, $\alpha_2$ is regarded as the value of the parameter
$\alpha$ denoted by \eqref{bf4} and is dependent on the eigenvalue
$\lambda_2$ with $\hat{\psi}_1,\hat{\psi}_2,\hat{\psi}_3$ representing
the associated eigenfunctions and
$\hat{\psi}_1^*,\hat{\psi}_2^*,\hat{\psi}_3^*$ the corresponding
adjoint eigenfunctions.

Subsequently, substituting $c_3=1-c_1-c_2$ into the first three
equations and performing an argument based upon the Brouwer fixed
point theorem and the topological degree as the one done above for
the case $|\beta|=1$, we ascertain the existence of a non-degenerate
solution of the algebraic system if the following condition is
satisfied:
\begin{equation*}
     \langle \hat \psi_1^*,y \cdot\nabla  (\hat{\psi}_3-\hat{\psi}_1) \rangle
     \langle \hat \psi_2^*,y \cdot\nabla  (\hat{\psi}_3-\hat{\psi}_2) \rangle-
     \langle \hat \psi_1^*,y \cdot\nabla  (\hat{\psi}_3-\hat{\psi}_2) \rangle
     \langle \hat \psi_2^*,y \cdot\nabla  (\hat{\psi}_3-\hat{\psi}_1) \rangle\neq 0.
\end{equation*}
Note that, by similar substitutions, other conditions might be
obtained.

Furthermore, once we know the existence of at least one solution,
we proceed now with a possible way of computing the number of
solutions of the nonlinear algebraic system \eqref{br61}.
Obviously, since the dimension of the eigenspace is bigger than
that in the case $|\beta|=1$, the difficulty in obtaining
multiplicity results increases.


Firstly, integrating by parts in the nonlinear terms, in which
$h_1$, $h_2$ and $h_3$ are involved, and rearranging terms in the
first three equations gives
\begin{align*}
&\int_{\mathbb{R}^N} \nabla \psi_1^* \cdot    \ln (c_1 \hat{\psi}_1+c_2\hat{\psi}_2
+ c_3 \hat{\psi}_3)  \nabla \Delta  (c_1 \hat{\psi}_1+ c_2 \hat{\psi}_2
+ c_3 \hat{\psi}_3)\\
& - c_1 \frac{\alpha_2}{4}
    \int_{\mathbb{R}^N} \hat \psi_1^* y \cdot\nabla  \hat{\psi}_1+c_1 \mu_{1,2}
    - c_2  \frac{\alpha_2}{4}
    \int_{\mathbb{R}^N} \hat\psi_1^* y \cdot\nabla  \hat{\psi}_2
    - c_3  \frac{\alpha_2}{4}
    \int_{\mathbb{R}^N} \hat\psi_1^* y \cdot\nabla  \hat{\psi}_3 =    0,
\\
& \int_{\mathbb{R}^N} \nabla \hat \psi_2^*\cdot
   \ln (c_1 \hat{\psi}_1+c_2\hat{\psi}_2+ c_3 \hat{\psi}_3) 
    \nabla \Delta
   (c_2 \hat{\psi}_1+ c_2 \hat{\psi}_2+ c_3 \hat{\psi}_3)\\
 &  -c_1  \frac{\alpha_2}{4}\int_{\mathbb{R}^N} \hat \psi_2^* y \cdot\nabla  \hat{\psi}_1
    + c_2 \mu_{1,2}
   -c_2  \frac{\alpha_2}{4}
    \int_{\mathbb{R}^N} \hat\psi_2^* y \cdot\nabla  \hat{\psi}_2
    -c_3  \frac{\alpha_2}{4}
   \int_{\mathbb{R}^N} \hat \psi_2^* y \cdot\nabla  \hat{\psi}_3=0,
\\
&\int_{\mathbb{R}^N} \nabla \hat \psi_3^*\cdot
   \ln (c_1 \hat{\psi}_1+c_2\hat{\psi}_2+ c_3 \hat{\psi}_3)
    \nabla \Delta
   (c_2 \hat{\psi}_1+ c_2 \hat{\psi}_2+ c_3 \hat{\psi}_3)\\
&  -c_1  \frac{\alpha_2}{4}\int_{\mathbb{R}^N} \hat \psi_3^* y \cdot\nabla  \hat{\psi}_1
    + c_3 \mu_{1,2}
    -c_2  \frac{\alpha_2}{4}\int_{\mathbb{R}^N} \hat\psi_3^* y \cdot\nabla  \hat{\psi}_2
    -c_3  \frac{\alpha_2}{4}\int_{\mathbb{R}^N} \hat \psi_3^* y \cdot\nabla  \hat{\psi}_3=0.
\end{align*}
According to the fourth equation, we have that $c_1 = 1-c_2-c_3$. Then,
setting
\begin{equation*}
    c_1 \hat{\psi}_1+c_2\hat{\psi}_2+ c_3 \hat{\psi}_3 =
    \hat{\psi}_1+c_2(\hat{\psi}_2-\hat{\psi}_1)+ c_3
    (\hat{\psi}_3-\hat{\psi}_1),
\end{equation*}
and substituting it into the expressions obtained above for the
first three equations of the system yield
\begin{equation} \label{br63}
\begin{aligned}
&\int_{\mathbb{R}^N} \nabla \hat \psi_1^* \cdot
   \ln (\hat{\psi}_1+
 (\hat{\psi}_2-\hat{\psi}_1)c_2+(\hat{\psi}_3-\hat{\psi}_1)c_3) \nabla \Delta
  \Big(\hat{\psi}_1+(\hat{\psi}_2- \hat{\psi}_1)c_2\\
& +   (\hat{\psi}_3- \hat{\psi}_1)c_3\Big)
 +\mu_{1,2} -c_2 \mu_{1,2}- c_3 \mu_{1,2} - \frac{\alpha_2}{4}\int_{\mathbb{R}^N} \hat \psi_1^* y \cdot\nabla  \hat{\psi}_1    \\
&  + \frac{\alpha_2}{4} \,\int_{\mathbb{R}^N} \hat{\psi}_1^* y
    \cdot( (\nabla  \hat{\psi}_1-\nabla \hat{\psi}_2)c_2 +
    (\nabla \hat{\psi}_1 - \nabla\hat{\psi}_3)c_3)  = 0,
\\
&\int_{\mathbb{R}^N} \nabla \hat \psi_2^* \cdot
   \ln (\hat{\psi}_1+
 (\hat{\psi}_2-\hat{\psi}_1)c_2+(\hat{\psi}_3-\hat{\psi}_1)c_3) \nabla \Delta
  \Big(\hat{\psi}_1\\
& +(\hat{\psi}_2- \hat{\psi}_1)c_2+
   (\hat{\psi}_3- \hat{\psi}_1)c_3\Big)
 +c_2 \mu_{1,2}
    - \frac{\alpha_2}{4}\int_{\mathbb{R}^N} \hat \psi_2^* y \cdot\nabla  \hat{\psi}_1\\
&  + \frac{\alpha_2}{4}\int_{\mathbb{R}^N} \hat{\psi}_2^* y
    \cdot( (\nabla  \hat{\psi}_1-\nabla \hat{\psi}_2)c_2 +
    (\nabla \hat{\psi}_1 - \nabla\hat{\psi}_3)c_3)  = 0,
\\
&\int_{\mathbb{R}^N} \nabla \hat \psi_3^* \cdot
   \ln (\hat{\psi}_1+
(\hat{\psi}_2-\hat{\psi}_1)c_2+(\hat{\psi}_3-\hat{\psi}_1)c_3) \nabla \Delta
 \Big(\hat{\psi}_1\\
& +(\hat{\psi}_2- \hat{\psi}_1)c_2+
   (\hat{\psi}_3- \hat{\psi}_1)c_3\Big) 
  +c_3\mu_{1,2} - \frac{\alpha_2}{4}
    \int_{\mathbb{R}^N} \hat \psi_3^* y \cdot\nabla  \hat{\psi}_1  \\ 
&  + \frac{\alpha_2}{4} \int_{\mathbb{R}^N} \hat{\psi}_3^* y
    \cdot( (\nabla  \hat{\psi}_1-\nabla \hat{\psi}_2)c_2 +
    (\nabla \hat{\psi}_1 - \nabla\hat{\psi}_3)c_3)  = 0.
 \end{aligned}
\end{equation}
Now, adding the first equation of \eqref{br63} to the other two,
we have 
\begin{align*}
&\int_{\mathbb{R}^N} (\nabla \hat \psi_1^*+ \nabla \hat{\psi}_2^*) \cdot
   \ln (\hat{\psi}_1+
     (\hat{\psi}_2-\hat{\psi}_1)c_2+(\hat{\psi}_3-\hat{\psi}_1)c_3) \nabla \Delta
 \Big(\hat{\psi}_1\\
& +(\hat{\psi}_2- \hat{\psi}_1)c_2+
   (\hat{\psi}_3- \hat{\psi}_1)c_3\Big)
  +\mu_{1,2}    - c_3 \mu_{1,2} - \frac{\alpha_2}{4}
    \int_{\mathbb{R}^N} (\hat \psi_1^*+ \hat{\psi}_2^*) y \cdot\nabla  \hat{\psi}_1 \\ 
&  + \frac{\alpha_2}{4}\int_{\mathbb{R}^N} (\hat \psi_1^*+ \hat{\psi}_2^*) y
    \cdot( (\nabla  \hat{\psi}_1-\nabla \hat{\psi}_2)c_2 +
    (\nabla \hat{\psi}_1 - \nabla \hat{\psi}_3)c_3)  = 0,
\end{align*}
\begin{align*}
&\int_{\mathbb{R}^N} (\nabla \hat \psi_1^*+ \nabla \hat{\psi}_3^*) \cdot
   \ln (\hat{\psi}_1+
 (\hat{\psi}_2-\hat{\psi}_1)c_2+(\hat{\psi}_3-\hat{\psi}_1)c_3) \nabla \Delta
  \Big(\hat{\psi}_1\\
& +(\hat{\psi}_2- \hat{\psi}_1)c_2+
   (\hat{\psi}_3- \hat{\psi}_1)c_3\Big)
 +\mu_{1,2}
    - c_2 \mu_{1,2} - \frac{\alpha_2}{4}
    \int_{\mathbb{R}^N} (\hat \psi_1^*+ \hat{\psi}_3^*) y \cdot\nabla  \hat{\psi}_1 \\ 
& + \frac{\alpha_2}{4}\int_{\mathbb{R}^N} (\hat \psi_1^*+ \hat{\psi}_3^*) y
    \cdot( (\nabla  \hat{\psi}_1-\nabla \hat{\psi}_2)c_2 +
    (\nabla \hat{\psi}_1 - \nabla \hat{\psi}_3)c_3)  = 0.
\end{align*}
Subsequently, subtracting those equations yields
\begin{align*}
\mu_{1,2}
&= \frac{1}{c_3-c_2} \,\big[ \int_{\mathbb{R}^N}
      (\nabla \hat \psi_2^*- \nabla \hat{\psi}_3^*) 
   \ln \Psi \nabla \Delta \Psi 
-\frac{\alpha_2}{4}\int_{\mathbb{R}^N} (\hat \psi_2^*- \hat{\psi}_3^*) y \cdot\nabla  
\hat{\psi}_1\\ 
& +\frac{\alpha_2}{4}\int_{\mathbb{R}^N} (\hat \psi_2^*- \hat{\psi}_3^*) y
    \cdot( (\nabla  \hat{\psi}_1-\nabla \hat{\psi}_2)c_2 +
    (\nabla \hat{\psi}_1 - \nabla \hat{\psi}_3)c_3) \big],
\end{align*}
where $\Psi = \hat{\psi}_1+
(\hat{\psi}_2-\hat{\psi}_1)c_2+(\hat{\psi}_3-\hat{\psi}_1)c_3$.
Thus, substituting it into \eqref{br63} (note that, from the
substitution into one of the last two equations, we obtain the
same equation), we arrive at the following system, with $c_2$ and
$c_3$ as the unknowns:
\begin{align*}
&c_3 \int_{\mathbb{R}^N} (\nabla \hat{\psi}_1^*-
   \nabla \hat{\psi}_2^*  
  +\nabla \hat{\psi}_3^*)    \ln \Psi
     \nabla \Delta \Psi 
- c_2 \int_{\mathbb{R}^N} (\nabla \hat{\psi}_1^*+
   \nabla \hat{\psi}_2^*-\nabla \hat{\psi}_3^*)
   \ln \Psi \nabla \Delta \Psi \\
&+ \int_{\mathbb{R}^N} (
   \nabla \hat{\psi}_2^*-\nabla \hat{\psi}_3^*) 
   \ln \Psi  \nabla \Delta \Psi- \frac{\alpha_2}{4} \int_{\mathbb{R}^N} (\hat{\psi}_2^*-\hat{\psi}_3^*)
   y \cdot\nabla \hat{\psi}_1 \\ 
&+c_2 \frac{\alpha_2}{4}
    [\int_{\mathbb{R}^N} (\hat{\psi}_2^*-\hat{\psi}_3^*) y \cdot\nabla
    (2\hat{\psi}_1 -\hat{\psi}_2)+ \int_{\mathbb{R}^N}
    \hat{\psi}_1^* y \cdot\nabla  \hat{\psi}_1] \\ 
&+ c_3 \frac{\alpha_2}{4} [\int_{\mathbb{R}^N} (\hat{\psi}_2^*-\hat{\psi}_3^*) y \cdot\nabla
    (2\hat{\psi}_1 -\hat{\psi}_3)- \int_{\mathbb{R}^N}
    \hat{\psi}_1^* y \cdot\nabla  \hat{\psi}_1]\\
& + c_2 c_3  \frac{\alpha_2}{4}[ \int_{\mathbb{R}^N} \hat{\psi}_1^* y
    \cdot( \nabla  \hat{\psi}_3-\nabla \hat{\psi}_2) -  \int_{\mathbb{R}^N}
    (\hat{\psi}_2^* - \hat{\psi}_3^*) y \cdot (2\nabla \hat{\psi}_1 -
    \nabla \hat{\psi}_2-\nabla \hat{\psi}_3)] \\ 
&+ c_3^2  \frac{\alpha_2}{4}\int_{\mathbb{R}^N} (\hat{\psi}_1^*-
    \hat{\psi}_2^*+\hat{\psi}_3^*) y
    \cdot( \nabla \hat{\psi}_1 - \nabla\hat{\psi}_3) \\ 
& - c_2^2  \frac{\alpha_2}{4}\int_{\mathbb{R}^N} (\hat{\psi}_1^*+
    \hat{\psi}_2^*-\hat{\psi}_3^*) y
    \cdot(\nabla  \hat{\psi}_1-\nabla \hat{\psi}_2) = 0,
\end{align*}
\begin{align*}
&c_3 \int_{\mathbb{R}^N} \nabla \hat{\psi}_2^* \ln \Psi
 \nabla \Delta \Psi - c_2 \int_{\mathbb{R}^N} \nabla \hat{\psi}_3^* 
   \ln \Psi \nabla \Delta \Psi
 - c_3 \frac{\alpha_2}{4} \int_{\mathbb{R}^N} \hat{\psi}_2^* y 
 \cdot\nabla  \hat{\psi}_1 \\
& +c_2 \frac{\alpha_2}{4}
    \int_{\mathbb{R}^N} \hat{\psi}_3^* y \cdot\nabla  \hat{\psi}_1
 + c_3 \frac{\alpha_2}{4}\int_{\mathbb{R}^N} \hat{\psi}_2^* y
    \cdot( (\nabla  \hat{\psi}_1-\nabla \hat{\psi}_2)c_2 +
    (\nabla \hat{\psi}_1 - \nabla\hat{\psi}_3)c_3) \\ 
&- c_2 \frac{\alpha_2}{4}\int_{\mathbb{R}^N} \hat{\psi}_3^* y
    \cdot( (\nabla  \hat{\psi}_1-\nabla \hat{\psi}_2)c_2 +
    (\nabla \hat{\psi}_1 - \nabla\hat{\psi}_3)c_3)  = 0.
\end{align*}
These can be re-written in the  form
\begin{equation}\label{br64}
    \begin{gathered}
 A_1 c_2^2+ B_1 c_3^2 + C_1 c_2 +D_1 c_3+ E_1 c_2 c_3 +\omega_1 (c_2,c_3)=0,
    \\ 
    A_2 c_2^2+ B_2 c_3^2 + C_2 c_2 +D_2 c_3+ E_2 c_2 c_3 +\omega_2 (c_2,c_3)=0,
    \end{gathered}
\end{equation}
 where
\begin{align*}
\omega_1 (c_2,c_3) 
&   := c_3 \int_{\mathbb{R}^N} (\nabla \hat{\psi}_1^*-
   \nabla \hat{\psi}_2^* +\nabla \hat{\psi}_3^*) 
   \ln \Psi      \nabla \Delta \Psi     \\ 
&\quad - c_2 \int_{\mathbb{R}^N} (\nabla \hat{\psi}_1^*+
   \nabla \hat{\psi}_2^*-\nabla \hat{\psi}_3^*) 
   \ln \Psi \nabla \Delta \Psi \\ 
&\quad  + \int_{\mathbb{R}^N} (
   \nabla \hat{\psi}_2^*-\nabla \hat{\psi}_3^*) \cdot
   \ln \Psi \nabla \Delta \Psi
- \frac{\alpha_2}{4} \int_{\mathbb{R}^N} (\hat{\psi}_2^*-\hat{\psi}_3^*)
   y \cdot\nabla \hat{\psi}_1,
\end{align*}
  and
\[   
\omega_2 (c_2,c_3)
    :=c_3 \int_{\mathbb{R}^N} \nabla \hat{\psi}_2^* \cdot
   \ln \Psi \nabla \Delta \Psi - c_2 \int_{\mathbb{R}^N} \nabla \hat{\psi}_3^* \cdot
   \ln \Psi \nabla \Delta \Psi,
\]
are the perturbations of the quadratic polynomials
 $$
 \mathfrak{F}_i(c_2,c_3) :=
A_i c_2^2+ B_i c_3^2 + C_i c_2 +D_i c_3+ E_i c_2 c_3,
$$
with $i=1,2$.
The coefficients
of those quadratic expressions are 
\begin{gather*}
 A_1:=  -\frac{\alpha_2}{4} \int_{\mathbb{R}^N} (\hat{\psi}_1^*+
    \hat{\psi}_2^*-\hat{\psi}_3^*) y
    \cdot(\nabla  \hat{\psi}_1-\nabla \hat{\psi}_2)  ,\\
  B_1:= \frac{\alpha_2}{4} \int_{\mathbb{R}^N} (\hat{\psi}_1^*-
    \hat{\psi}_2^*+\hat{\psi}_3^*) y
    \cdot( \nabla \hat{\psi}_1 - \nabla\hat{\psi}_3) , \\
  C_1:=
\frac{\alpha_2}{4}
    [\int_{\mathbb{R}^N} (\hat{\psi}_2^*-\hat{\psi}_3^*) y \cdot\nabla
    (2\hat{\psi}_1 -\hat{\psi}_2)+ \int_{\mathbb{R}^N}
    \hat{\psi}_1 y \cdot\nabla  \hat{\psi}_1] ,\\
 D_1:=    \frac{\alpha_2}{4}\, [\int_{\mathbb{R}^N} (\hat{\psi}_2^*-\hat{\psi}_3^*) y 
 \cdot\nabla    (2\hat{\psi}_1 -\hat{\psi}_3)- \int_{\mathbb{R}^N}
    \hat{\psi}_1 y \cdot\nabla  \hat{\psi}_1] , \\
  E_1:=     \frac{\alpha_2}{4}[\int_{\mathbb{R}^N} \hat{\psi}_1^* y
    \cdot( \nabla  \hat{\psi}_3-\nabla \hat{\psi}_2) -  \int_{\mathbb{R}^N}
    (\hat{\psi}_2^* - \hat{\psi}_3^*) y \cdot (2\nabla \hat{\psi}_1 -
    \nabla \hat{\psi}_2-\nabla \hat{\psi}_3)], \\
     A_2:=  -\frac{\alpha_2}{4}\int_{\mathbb{R}^N} \hat{\psi}_3^* y
    \cdot( \nabla  \hat{\psi}_1-\nabla \hat{\psi}_2), \\
 B_2:= \frac{\alpha_2}{4}\int_{\mathbb{R}^N} \hat{\psi}_2^* y
    \cdot(    (\nabla \hat{\psi}_1 - \nabla\hat{\psi}_3) , \\
  C_2:=\frac{\alpha_2}{4}\int_{\mathbb{R}^N} \hat{\psi}_3^* y \cdot\nabla  \hat{\psi}_1, \quad
  D_2:=  -  \frac{\alpha_2}{4}\int_{\mathbb{R}^N} \hat{\psi}_2^* y \cdot\nabla  \hat{\psi}_1,\\
 E_2:=   \frac{\alpha_2}{4}\int_{\mathbb{R}^N} \hat{\psi}_2^* y
    \cdot (\nabla  \hat{\psi}_1-\nabla \hat{\psi}_2) -
    \hat{\psi}_3^* y
    \cdot (\nabla \hat{\psi}_1 - \nabla\hat{\psi}_3).
\end{gather*}


Therefore,  using the conic classification to solve \eqref{br64},
we have the number of solutions through the intersection of
two conics. Then, depending on the type of conic, we shall always
obtain one to four possible solutions for our system. Hence,
somehow, the number of solutions depends on the coefficients we
have for the system and, at the same time, on the eigenfunctions
that generate the subspace $\ker\big(\mathbf{B}+\frac{k}{4}\big)$.


Thus, we have the following conditions, which will provide us with
the conic section of each equation of the system \eqref{br64}:
\begin{itemize}
\item[(i)] If $B_i^2 - 4A_iE_i < 0$, the equation represents an \emph{ellipse},
unless the conic is degenerate, for example $c_2^2 + c_3^2 + k = 0$
for some positive constant $k$. So, if $A_i=B_j$ and $E_i=0$,
the equation represents a \emph{circle};

\item[(ii)] If $B_i^2 - 4A_iE_i = 0$, the equation represents a \emph{parabola};

\item[(iii)] If $B_i^2 - 4A_iE_i > 0$, the equation represents a \emph{hyperbola}.
If we also have $A_i + E_i = 0$ the equation represents a
hyperbola (a rectangular hyperbola).
\end{itemize}
Consequently, the zeros of the system \eqref{br64} and, hence, of
the system \eqref{br61}, adding the ``normalizing" constraint
\eqref{br15}, are ascertained by the intersection of those two
conics in \eqref{br64} providing us with the number of possible
$n$-branches between one and {four}. Note that in case  those
conics are two circles we only have two intersection points at
most. Moreover, due to the dimension of the eigenspaces it looks
like in this  case that we have four possible intersection points
two of them will coincide. However, the justification for this is
far from clear.

Moreover, as it was done for the
previous case when $|\beta|=1$, we need to control the oscillations
of the perturbation functions in order to maintain the number of
solutions. Therefore, imposing that
\begin{equation*}
\| \omega_i (c_2,c_3)\|_{L^\infty} \leq \mathfrak{F}_i(c_2^*,c_3^*),
    \quad \text{with }i=1,2,
\end{equation*}
we ascertain that the number of solutions must be between one and four.
This again gives us an idea of the difficulty of more general
multiplicity results.

\subsection{Further comments on mathematical justification of existence}
 \label{S6}

We return to the self-similar nonlinear eigenvalue  problem
\eqref{sf5}, associated with \eqref{i1},
which can be written in the  form
\begin{equation*}
 \mathcal{L}(\alpha,n)f + \mathcal{N}(n,f)=0  ,\quad \text{where }
    \mathcal{N}(n,f):=\nabla \cdot ((1-|f|^n)\nabla \Delta f)\,.
\end{equation*}
As we have seen, the main difficulty in justifying the
$n$-branching behaviour concerns the distribution and
``transversal topology" of zero surfaces of solutions close to
finite interface hyper-surfaces.

 Recall that,  as in classic nonlinear operator theory \cite{Deim, KZ,
 VainbergTr}, our  analysis above always assumed that we actually
 dealt with and performed computations for the
 integral equation:
\begin{equation}\label{cp8}
f= -\mathcal{L}^{-1}(\alpha,n) \mathcal{N}(n,f) \equiv \mathcal{G}(n,f),\quad
    \mathcal{L}(\alpha,n):= -\Delta^2 +\frac{1-\alpha n}{4} y \cdot \nabla +\alpha I ,
\end{equation}
where $\mathcal{L}(\alpha,n)$ is invertible in $L^2_\rho$ (this is
directly checked via Section \ref{S3}) and, hence compact, for a
fixed $\alpha$, and $f \in C_0(\mathbb{R}^N)$ for small $n>0$. This confirms
that the zeros of the function $\mathcal{F}(n,f)$ are fixed points
of the map $\mathcal{G}(n,f)$.

Note again that \eqref{cp8} is an
eigenvalue problem, where admissible real values of $\alpha$ are
supposed to be defined together with its solvability.
 This makes existence/multiplicity questions for \eqref{cp8}
 extremely difficult.

There are two cases of this problem. The first and simpler one
occurs when the eigenvalue $\alpha$ is determined \emph{a priori},
e.g. in the case $k=0$, where $\alpha_0(0)=  N/4$ denoted as
$\alpha_0(0)=\alpha_0$, and where, for $n>0$, the first nonlinear
eigenvalue is given explicitly (see \eqref{alb1}):
 $$
\alpha_0(n)= \frac N{4+Nn}.
 $$
Then \eqref{cp8} with $\alpha=\alpha_0(n)$ for $n>0$ becomes a standard
nonlinear integral equation with, however,  a quite curious and
hard-to-detect functional setting. Indeed, the right-hand side in
\eqref{cp8}, where the nonlinearity is not in a fully divergent
form, assumes the extra regularity at least such as
 \begin{equation*}
 f \in H_\rho^3.
 \end{equation*}
In view of the known good properties of the compact resolvent
$(\mathcal{L}-\lambda I)^{-1}$, it is clear that the action of the
inverse one $\mathcal{L}^{-1}$ is sufficient to restore the
regularity, since locally in $\mathbb{R}^N$ this acts like $\Delta^{-2}$.
Therefore, it is plausible that
 \begin{equation*}
\mathcal{G}: H^3_\rho \to H^3_\rho,
 \end{equation*}
 and it is not difficult to get an \emph{a priori} bounds at least
 for small enough $f$'s. The accompanying  analysis as $y \to \infty$ (due to the unbounded domain)
  assumes no
 novelties or special difficulties and is standard for such
 weighted $L^2$ and Sobolev spaces.


  Therefore, application of Schauder's fixed point theorem
 (see e.g. \cite[p.~90]{Berger})    to
  \eqref{cp8} is the most powerful tool to imply existence of a
  solution, and moreover a continuous curve of fixed points 
$\Gamma_n=\{f: n>0 \text{ small}\}$.

By scaling invariance of the
  similarity equation, we are obliged to impose the normalization
  condition, say,
   \begin{equation*}
   f(0)= \delta_0>0 \quad \text{sufficiently small}.
    \end{equation*}

On the other hand, uniqueness remains a completely open problem (apart from partial results such as \cite{PVcauchy} when $n$ is sufficiently close to zero). However, studying the
 behaviour of the solution curve $\Gamma_n$ as $n \to 0$ and
 applying (under suitable hypothesis) the branching techniques
 developed above,
  we may conclude that any such continuous
 curve must be originated at a properly scaled eigenfunction
 $\psi_0=F$, so that such a curve is unique due to well-posedness
 of all the asymptotic expansions.

 A possibility of extension of $\Gamma_n$ for larger values of
 $n>0$ represents an essentially more difficult nonlocal open problem. Indeed, via
 compactness of linear operators involved in \eqref{cp8}, one can
 expect that such a curve can end up at a bifurcation
 point only (unless it blows up). However, nonexistence of turning
 saddle-node points at some $n_*>0$ (meaning that the $n$-branch is nonexistent for
 some $n>n_*$) is not that easy to rule
 out. Moreover, such turning points with thin film
 operators involved
 are actually possible, \cite{GalPetII}.


After establishing the existence of such solutions for small $n>0$, we
face the next problem on their asymptotic properties including the
fact that these are compactly supported. On a qualitative level,
these questions were discussed in \cite{EGK1}.

In the case of higher-order nonlinear eigenfunctions of
\eqref{cp8} for $k \ge 1$ including the dipole case $k=1$, the
parameter $\alpha$ becomes an eigenvalue that is essentially involved
into the problem setting. This assumes the consideration of the equation
\eqref{cp8} in the extended space
 \begin{equation}
 \label{ppp4}
 (f,\alpha) \in X= H^3_\rho \times \{\alpha \in \mathbb{R}\} \quad
 \text{and } \mathcal{G}: X \to X,
  \end{equation}
  where  proving the latter  mapping for some compact subsets becomes a
  hard   open problem. Note that here even the necessary convexity issue
  for applying Schauder's Theorem
  can be difficult.
   We still do not know whether representations such
  as \eqref{ppp4} may lead to any rigorous treatment of the
  nonlinear eigenvalue problem \eqref{cp8} for $k \ge 1$.

\subsection*{Acknowledgments}
This work was partially supported by the Ministry of Economy and
 Competitiveness of Spain under research project MTM2012-33258.

\begin{thebibliography}{99}


\bibitem{AEG1}  P. Alvarez-Caudevilla, J. D. Evans,  V. A. Galaktionov;
 \emph{The Cauchy problem for  a tenth-order thin film equation I.
Bifurcation of self-similar oscillatory fundamental solutions}, Mediterranean Journal of Mathematics,
No. \textbf{4}, Vol. 10 (2013),  pp. 1759--1790.


\bibitem{TFE4PV} P. Alvarez-Caudevilla, V. A. Galaktionov;
 \emph{Local bifurcation-branching analysis of global and``blow-up" patterns
 for a fourth-order thin film equation},
 {Nonlinear Differ. Equat.  Appl.},  \textbf{18} (2011), 483--537.

 \bibitem{PVcauchy}
 P. Alvarez-Caudevilla and V.A.~Galaktionov;
 \emph{Well-posedness of the Cauchy
problem for a fourth-order thin film equation via regularization approaches}, 
Nonlinear Analysis, in press.

\bibitem{Ben81}  P. B\'enilan, M. G. Crandall;
 \emph{The continuous dependence on $\varphi$ of solutions
 of $u_t-\Delta \varphi(u)=0$}, {Indiana Univ. Math.~J.,} 
\textbf{30} (1981), 161--177.

\bibitem{BBP} E. Beretta, M. Bertsch,  R. Dal Passo;
 \emph{Nonnegative solutions of a fourth-order nonlinear degenerate parabolic
equation}, {Arch. Rational Mech. Anal.}, \textbf{129} (1995),
175--200.

\bibitem{Berger}   M.~Berger;
\emph{Nonlinearity and Functional Analysis}, Acad.
  Press, New York, 1977.

\bibitem{BF1}  F.~Bernis, A.~Friedman;
 \emph{Higher order nonlinear degenerate
 parabolic equations}, {J.~Differ. Equat.,} \textbf{83} (1990), 179--206.


\bibitem{BPW} F. Bernis, L. A. Peletier,  S. M. Williams;
\emph{Source type solutions of a fourth order nonlinear degenerate
parabolic equation}, {Nonlinear Anal.}, \textbf{18} (1992),
217--234.

\bibitem{BHQ} F. Bernis, J. Hulshof, F. Quir\'os;
\emph{The ``linear" limit of thin film flows as an obstacle-type free 
boundary problem}, {SIAM J. Appl. Math.}, \textbf{61} (2000), 1062--1079.

\bibitem{BS}  M. S. Birman, M. Z. Solomjak;
\emph{Spectral Theory of Self-Adjoint Operators in Hilbert Spaces},
{D. Reidel}, Dordecht/Tokyo (1987).


\bibitem{BW06} M.~Bowen, T. P. Witelski;
\emph{The linear limit of  the dipole problem for the thin film equation}, 
{SIAM J.~Appl.  Math.,} \textbf{66} (2006), 1727--1748.

 \bibitem{Cock99}   B.~Cockburn, G.~Gripenberg;
\emph{Continuous dependence on the nonlinearities of solutions of
  degenerate parabolic equations}, {J.~Differ.
 Equat.,} \textbf{151} (1999), 231--251.

\bibitem {Deim} K. Deimling;
\emph{Nonlinear Functional Analysis},
Springer-Verlag, Berlin/Tokyo, 1985.


\bibitem{EGKP}  Yu. V. Egorov, V. A. Galaktionov, V. A. Kondratiev, S.I. Pohozaev;
\emph{Global solutions of higher-order semilinear parabolic
equations in the supercritical range}, {Adv. Differ. Equat.,}
\textbf{9} (2004),  1009--1038.

\bibitem{EGK1} J. D. Evans, V. A. Galaktionov,  J. R. King,
\emph{Blow-up similarity solutions of the fourth-order unstable
thin film equation}, {European J. Appl. Math.,} \textbf{18}
(2007),  195--231.

\bibitem{EGK2} J. D. Evans, V. A. Galaktionov,  J. R. King;
\emph{Source-type solutions of the fourth-order unstable thin film
equation}, {European J. Appl. Math.,} \textbf{18} (2007), 273--321.

\bibitem{EGK3} J. D. Evans, V. A. Galaktionov,  J. R. King;
\emph{Unstable sixth-order thin film equation I. Blow-up
similarity solutions}, {Nonlinearity}, \textbf{20} (2007),
1799--1841.

\bibitem{EGK4} J. D. Evans, V. A. Galaktionov, J.R. King;
\emph{Unstable sixth-order thin film equation II. Global
similarity patterns}, {Nonlinearity}, \textbf{20} (2007),
1843--1881.


\bibitem{GalGeom} V. A. Galaktionov;
\emph{Geometric Sturmian Theory of Nonlinear
 Parabolic Equations and Applications}, Chapman \& Hall/CRC, Boca Raton,
Florida,  2004.

\bibitem{GalRDE4n} V. A. Galaktionov;
 \emph{Countable branching of similarity solutions of higher-order
 porous medium type equations}, Adv. Differ. Equat., \textbf{13}
 (2008), 641--680.

\bibitem{GalPetII}  V. A. Galaktionov;
 \emph{Very singular solutions for thin film equations with absorption},
Studies in Applied Mathematics, \textbf{124} (2010), 39--63.

\bibitem{GHCo} V. A. Galaktionov, P. J. Harwin;
\emph{On evolution completeness of nonlinear
 eigenfunctions for the porous medium equation in the whole space},
 Adv. Differ. Equat., \textbf{10} (2005), 635--674.


\bibitem{GMPSobI}  V. A. Galaktionov, E. Mitidieri, S. I. Pohozaev;
 \emph{Variational approach to complicated similarity solutions
of higher-order  nonlinear evolution equations of parabolic, hyperbolic,
 and nonlinear dispersion  types},
In: Sobolev Spaces in Mathematics. II, Appl. Anal. and Part.
Differ. Equat., Series: Int. Math. Ser., Vol. \textbf{9}, V. Maz'ya
Ed., Springer, New York, 2009 (an earlier preprint:
arXiv:0902.1425).

\bibitem{Gia08} L. Giacomelli, H. Kn\"upfer,  F. Otto;
\emph{Smooth zero-contact-angle solutions to a  thin film equation
 around the steady state},  {J.~Differ. Equat.,} \textbf{245} (2008), 
1454--1506.

\bibitem{Grun04} G. Gr\"un;
\emph{Droplet spreading under weak slippage - existence for the Cauchy problem}, 
{Commun. Partial Differ. Equat.,} \textbf{29} (2004), 1697--1744.

\bibitem{Kal78}   A. S. Kalashnikov;
\emph{On continuous dependence of generalized solutions of the
equation of unsteady filtration on a function determining the flow
mode}, {J.~Appl. Math. Mech.,} \textbf{42} (1978), 183--185.

\bibitem{Ki01} J. R. King;
\emph{Two generalisations of the thin film equation}, 
{Math. Comput. Modelling,} \textbf{34} (2001), 737--756.

\bibitem{KZ} M. A. Krasnosel'skii, P. P. Zabreiko;
\emph{Geometrical Methods of Nonlinear Analysis}, Springer-Verlag,
Berlin/Tokio, 1984.

\bibitem{Pan08}  J. Pan;
\emph{The approximating character on nonlinearities of solutions of Cauchy problem
   for a singular diffusion equation,}
   {Osaka J. Math.,}  \textbf{45} (2008), 909--919.

\bibitem{Pan07} J. Q. Pan, L. Gang;
 \emph{The linear approach for a nonlinear infiltration equation},
 {Euro. J.~Appl. Math.,}  \textbf{17} (2006), 665--675.

\bibitem{PT} L. A. Peletier, W. C. Troy; Spatial Patterns.
Higher Order Models in Physics and Mechanics, Birkh\"{a}usser,
Boston/Berlin, 2001.

\bibitem{VainbergTr} M. A. Vainberg, V. A. Trenogin;
\emph{Theory of Branching of Solutions of Non-Linear Equations}, 
Noordhoff Int. Publ., Leiden, 1974.

\end{thebibliography}

\end{document}