\section{Introduction \label{sec:1}}

The study of Cosmic Microwave Background (CMB) photons using the
\textit{Planck} \cite{Planck} and BICEP2 \cite{Bicep2} data sets has opened up
a new constraint on inflationary models~\cite{encyclo}. Specifically, the
detailed analysis of \textit{Planck} data~\cite{Planck} constrains single
scalar-field models of slow-roll inflation, to have very low tensor-to-scalar
fluctuation ratio $r=n_{t}/n_{s}\ll1$, with a scalar spectral index
$n_{s}\simeq0.96$ possessing no appreciable running. The upper bound set by
the \textit{Planck} Collaboration~\cite{Planck} on tensor-to-scalar
fluctuation ratio, related to the absence of the B-modes of polarization in
the CMB, is $r<0.11$, but their favored regions (higher than 95\% C.L.) point
towards a more stringent bound of $r\leq10^{-3}$. From the theoretical view
point one can see that this is in agreement with the so-called
Starobinsky-type (or $R^{2}$, with $R$ denoting the scalar space-time
curvature) inflationary models~\cite{staro}.

Last year, however, the BICEP2 team~\cite{Bicep2} made an important
claim: to have made the first measurement of B-mode polarization in the CMB
radiation. This measurement was initially interpreted as indication
for gravitational waves at the time of the last scattering, with a
tensor-to-scalar ratio $r=0.16_{-0.05}^{+0.06}$.
If such claims were confirmed, they would constitute the first experimental
observation of (transverse) primordial (possibly quantum) metric fluctuations.
For the scalar spectal index, it has been found that $n_{s}\simeq0.96$
and $dn_{s}/d\mathrm{ln}k\simeq0$, in agreement with the \textit{Planck}
data~\cite{Planck}.

Since then a great effort has been spent in order to compare the BICEP2
tensor-to-scalar ratio against the choice of inflationary paradigm ($R^{2}$
\cite{staro}, chaotic \cite{Linde}, inverse power law \cite{Barr}, hilltop
\cite{Alb}, natural \cite{Free}, supersymmetry \cite{Mav}, D-flation
\cite{Mav1} and the like).

Furthermore, during this period there was an intense debate as to whether the
BICEP2 signal is indubitably due to primordial gravitational waves, or is
polluted by gravitationally-lensed E-modes and Galactic foregrounds.
Recently,
in~ref. \cite{Subir} it was stressed
that magnetized dust associated with radio loops due to supernova remnants
might affect the signal received by BICEP2.
In the case of polarization effects by galactic dust, refs. ~\cite{MS}
and~\cite{Flauger} have shown that the cosmological value of tensor-to-scalar
ratio could be very small $r\ll0.1$. Within this framework, the recent
analysis on the foreground dust in the BICEP2 region released by the
\textit{Planck} collaboration~\cite{Planckrec} points to a significant
foreground pollution which means that the BICEP2 B-mode polarization data
cannot be used as evidence for primordial CMB polarization.

Recently, the joint analysis of BICEP2/Keck Array and \textit{Planck} data
appeared in the literature \cite{joint}. This analysis has placed an upper
bound in the tensor-to-scalar ratio, namely $r<0.12$ at $95\%$ significance
level. In this context, the Planck team has repeated the inflationary analysis
by using the \textit{Planck 2015} data and essentially confirm the
\textit{Planck 2013} results: $n_{s}=0.968\pm0.006$, $dn_{s}/d\mathrm{ln
k=-0.003\pm0.007$ and $r<0.11$ \cite{Planck2015}.

The crux of these studies is that the potential energy of the scalar field, is
not really known and one must introduce it using phenomenological arguments,
starting with the simplest possibilities. There has been intense debate and
speculation about the functional form of the potential energy $V(\phi)$.
Various candidates have been proposed in the literature, such as a power law,
inverse power law, exponential and so on (for review see \cite{Lyth} and
references therein) but exact solutions are not abundant and always possess
special mathematical features.

Some time ago, a special solution for a spatially flat
Friedmann-Lema\^{\i}tre-Robertson-Walker (FLRW) spacetime with a perfect fluid
with a constant equation of state parameter\textbf{ }$P_{m}=(\gamma-1)\rho
_{m}$\textbf{ }(where, for radiation $\gamma=4/3$ and for dust $\gamma=1$) and
a scalar field with a constant equation of state
$P_{\phi}=(\gamma_{\phi}-1)\rho_{\phi}$ was found in \cite{Barrow}. In
particular, Rubano and Barrow \cite{Barrow} (see also \cite{Urena} and
\cite{Saa}) showed that under of specific conditions we can solve the Einstein
equations when the potential $V(\phi)$ has the interesting hyperbolic form:
\begin{equation}
V\left(  \phi\right)  =A\left[  \sinh\left(  \sqrt{3}\frac{(\gamma
-\gamma_{\phi})}{\sqrt{\gamma_{\phi}}}\left(  \phi-\phi_{0}\right)  \right)
\right]  ^{b}, \label{pot
\end{equation}
where the constant $A$ is
\begin{equation}
A=3H_{0}^{2}\left(  1-\Omega_{m0}\right)  \left(  1-\frac{\gamma_{\phi}
{2}\right)  \left(  \frac{1-\Omega_{m0}}{\Omega_{m0}}\right)  ^{-b/2}
\label{AA
\end{equation}
and
\begin{equation}
b=-\frac{2\gamma_{\phi}}{\gamma-\gamma_{\phi}}=\frac{2(1+w_{\phi})}{1+w_{\phi
}-\gamma}. \label{bb
\end{equation}
Note that $H_{0}$ and $\Omega_{m0}$ are the usual cosmological parameters
(although $\Omega_{m0}$ denotes any possible matter content according to the
appropriate choice of $\gamma$). If we trace the late universe (dustlike
matter $P_{m}=0$) we have $\gamma=1$.
We remind the reader that $w_{\phi}$ denotes the equation of state parameter
of the dark energy usually parametrized by $w_{\phi}=\frac{P_{\phi}
{\rho_{\phi}}=\gamma_{\phi}-1$, with $P_{\phi}$ and $\rho_{\phi}$ being the
pressure and density of the dark energy fluid.
Obviously the exponent $b$ plays an important role in the cosmic dynamics,
since it is related with the equation of state parameter, $w_{\phi}$. In this
context, inverting Eq.(\ref{bb}) one can prove
that the equation of state parameter reduces to $w_{\phi}=2/(b-2)$. Note that
the accelerated expansion of the universe poses the restriction $w_{\phi
}<-1/3\Omega_{de0}$ which implies $2(1-3\Omega_{de0})<b<2$, where
$\Omega_{de0}=1-\Omega_{m0}$.

The potential (\ref{pot}) has some interesting geometric characteristics.
Under specific conditions, it behaves either as an exponential or as a
power-law. If $|\lambda\phi|\gg1$ (or $|\lambda\phi|\ll1$),
we find $V\propto e^{-b\lambda\phi}$ (or $V\propto(\lambda\phi)^{b}$), where
$\lambda=\frac{\sqrt{3}(1-\gamma_{\phi})}{\sqrt{\gamma_{\phi}}}$ (see also
\cite{Saa}). The initial motivation in \cite{Barrow} was to use potential
(\ref{pot}) to describe the late-time acceleration of the universe, but we can
also apply Eq.(\ref{pot}) to the very early states of the cosmic evolution
when the cosmic fluid is dominated by the radiation and the inflaton
components.\textbf{ }This might then provide a unified picture of inflation
and dark energy in which both eras are described by the potential of
Eq.(\ref{pot}). It is the purpose of this work to demonstrate compatibility of
the Rubano and Barrow \cite{Barrow} scenario with the \textit{Planck 2015}
data,
taking into account the foreground ambiguities clouding the \textit{Planck
2015} data, as mentioned previously. The structure of the paper is as follows.
The slow-roll inflation and its connection to the hyperbolic potential of ref.
\cite{Barrow} are reviewed in section II. In section III we study the
performance of Eq.(\ref{pot}) against the \textit{Planck 2015} data.
Finally, our conclusions are summarized in section IV.

\begin{figure}[ptb]
\includegraphics[width=0.48\textwidth]{Fig1}\caption{The $n_{s}-r$ diagram for
the hyperbolic potential of \cite{Saa} using $N=55$. The contours borrowed
from \emph{Planck 2015} \cite{Planck2015}. The area which is plotted over the
contours corresponds to the hyperbolic Rubano \& Barrow \cite{Barrow}
potential (see Eq.(\ref{pot1}). The points correspond to chaotic (solid point)
and Starobinsky (star) inflation respectively.
\label{fig:1
\end{figure}

\begin{figure}[ptb]
\caption{The $n_{s}-n_{s}^{\prime}$ diagram. The area corresponds to the
hyperbolic Rubano \& Barrow \cite{Barrow} potential (for more details see the
caption of figure 1).
\label{fig:2
\includegraphics[width=0.48\textwidth]{Fig2}\end{figure}

\section{Slow-roll inflation}

Let us present here the basic ingredients in the context of single-field
inflation. Assuming an inflaton field $\phi$ with potential energy $V(\phi)$,
the slow-roll parameters are given by
\begin{equation}
\epsilon=\frac{M_{pl}^{2}V^{^{\prime}2}}{2V^{2}}, \label{ee
\end{equation
\begin{equation}
\eta=\frac{M_{pl}^{2}V^{^{\prime\prime}}}{V} \label{nn
\end{equation
\begin{equation}
\xi=\frac{M_{pl}^{4}V^{^{\prime}}V^{^{\prime\prime\prime}}}{V^{2}}, \label{xx
\end{equation}
where the prime denotes derivatives with respect to $\phi$ and $M_{pl
^{2}=1/8\pi G$. The corresponding spectral indices are defined in terms of the
slow-roll parameters, as usual \cite{LL}, by
\begin{equation}
n_{s}\simeq1+2\eta-6\epsilon\label{sp1
\end{equation
\begin{equation}
r\simeq16\epsilon\simeq-8n_{t} \label{sp2
\end{equation
\begin{equation}
n_{s}^{^{\prime}}=dn_{s}/d\mathrm{ln}k\simeq16\epsilon\eta-24\epsilon^{2
-2\xi. \label{sp3
\end{equation}
In this framework, the number of e-folds is written as\footnote{In the
literature sometimes we label $\phi$ by $\phi_{\star}$ which denotes the value
at the horizon crossing for which the scalar amplitude is $A_{s}\approx
\Lambda^{4}/24\pi^{2}\epsilon M_{pl}^{4}$, where the energy scale of inflation
is $\Lambda\sim10^{16}$Gev.}
\begin{equation}
N=\int_{t}^{t_{end}}H(t)dt\simeq\frac{1}{M_{pl}^{2}}\int_{\phi_{end}}^{\phi
}\frac{V(\phi)}{V^{^{\prime}}(\phi)}d\phi\;, \label{efold
\end{equation}
where $\phi_{end}$ is the value of the inflaton field at the end of inflation,
namely $\epsilon(\phi_{end})\simeq1$. Now, let us focus on the potential of
Eq.(\ref{pot}) which is written as
\begin{equation}
V(\phi)=A\;\mathrm{sinh}^{b}(\phi/f), \label{pot1
\end{equation}
where $f$ is the scale in units of $M_{pl}$.
Recall that in the dark-energy era the corresponding constants are related
with the cosmological parameters [see Eqs.(\ref{AA}) and (\ref{bb})].
Therefore, in the early universe one may expect that the corresponding
constants in Eq.(\ref{pot1}) are not necessarily equal to those derived by
\cite{Barrow} using arguments from the late universe where pressureless matter
dominates radiation (for comparison with the observational data see
\cite{Saa}). When radiation dominates we have $\gamma=4/3$ and so from the
second equality of Eq.(\ref{bb}) we arrive at
\begin{equation}
\label{bbw}b=\frac{6(1+w_{\phi})}{3w_{\phi}-1}.
\end{equation}
Note that the restriction $b>1$ implies that
the potential contains a critical point.

Now inserting Eq.(\ref{pot1}) into the slow-roll parameters we obtain after
some simple algebra
\begin{equation}
\epsilon=\frac{b^{2}M_{pl}^{2}}{2f^{2}}\mathrm{coth}^{2}(\phi/f), \label{ee1
\end{equation
\begin{equation}
\eta=\frac{bM_{pl}^{2}}{f^{2}}\left[  (b-1)\mathrm{coth}^{2}(\phi/f)+1\right]
, \label{nn1
\end{equation
\begin{equation}
\xi=\frac{b^{2}M_{pl}^{4}}{f^{4}}\mathrm{coth}^{2}(\phi/f)\left[
(b-1)(b-2)\mathrm{coth}^{2}(\phi/f)+(3b-2)\right]  . \label{xx1
\end{equation}
Combining Eq.(\ref{ee1}) with Eq.(\ref{nn1}) and Eq.(\ref{xx1}) we can write
$(\eta,\xi)$ in terms of $\epsilon$:
\begin{equation}
\eta=\frac{bM_{pl}^{2}}{f^{2}}\left[  \frac{2f^{2}\epsilon}{b^{2}M_{pl}^{2
}(b-1)+1\right]  \label{nn2
\end{equation}

\begin{equation}
\label{xx2}\xi=\frac{2M_{pl}^{2} \epsilon}{f^{2}} \left[  \frac{2f^{2
(b-1)(b-2)\epsilon}{b^{2}M_{pl}^{2}}+(3b-2) \right]  .
\end{equation}

In this case, the number of e-folds becomes
\begin{equation}
N\simeq\frac{f^{2}}{bM_{pl}^{2}}\mathrm{ln}\left[  \frac{\mathrm{cosh
(\phi/f)}{\mathrm{cosh}(\phi_{end}/f)}\right]  . \label{efold1
\end{equation}
Inverting the above, we can express the inflaton as $\phi(N)$ by
\begin{equation}
\phi=f\;\mathrm{cosh}^{-1}\left[  e^{NbM_{pl}^{2}/f^{2}}\mathrm{cosh
(\phi_{end}/f)\right]  . \label{efold2
\end{equation}
In order to proceed with the analysis we need to know the values of $N$ and
$\phi_{end}$. Firstly, it is natural to consider that the number of e-folds
lies in the interval $[50,60]$. Here, we set it to 55, for concreteness.
Secondly, using the constraint $\epsilon(\phi_{end})\simeq1$ and
Eq.(\ref{ee1}) we can estimate the value of the scalar field at the end of
inflation to be
\begin{equation}
\phi_{end}\simeq\frac{f}{2}\mathrm{ln}\left(  \frac{\theta+1}{\theta
-1}\right)  ~, \label{efold11
\end{equation}
where $\theta=\sqrt{2}f/bM_{pl}$. Since $\theta>1,$ the scale $f$ obeys the
restriction $f>\sqrt{2}bM_{pl}/2$.

\section{Observational restrictions}

The point of this section is to test the viability of the potential
(\ref{pot1}) at the inflationary level, involving the latest cosmological
data. In particular, the combined analysis between \emph{Planck 2015} and
various data such as WMAP, high-$l$ data, and Baryonic Acoustic Oscillations
(BAO), shows that the scalar spectral index is $n_{s}=0.968\pm0.006$. For the
scalar spectral index in this work we use $n_{s}^{\prime}=-0.003\pm0.007$.
Furthermore, as mentioned in the introduction, the analysis of the
\emph{Planck} collaboration places an upper limit on the tensor-to-scalar
ratio is concerned, $r<0.11,$ which is in agreement with the joint analysis of
BICEP2/Keck Array and \textit{Planck} \cite{joint}.

Let us now briefly present our results. In figures 1 and 2 we show the
confidence contours in the $(n_{s},r)$ and $(n_{s},n_{s}^{\prime})$ planes
which are provided by the Planck team \cite{Planck2015}. On top of that we
present the area for the individual sets of $(n_{s},r)$ which are based on the
potential (\ref{pot1}), whereas in figure 2 we display the corresponding area
in the case of $(n_{s},n_{s}^{\prime})$. Obviously, our results are consistent
with those of \textit{Planck 2015}. Specifically, as it can be seen from
figure 1, the tensor-to-scalar fluctuation ratio could reach the value of
$r\simeq0.075$ which is in a good agreement within $1\sigma$ with that of
BICEP2/Keck Array/\textit{Planck} results ($r\simeq0.05$ see figure 9 in
\cite{joint}). Notice that in order to derive $r\simeq0.07$ the constants in
Eq.(\ref{pot1}) need to obey the following inequalities: $1.02\leq b\leq1.1$
[or inverting Eq.(\ref{bbw}) we obtain $-2.63 \le w_{\phi} \le-2.39$]
and $26M_{pl}\leq f\leq39M_{pl}$. Concerning the running spectral index we
obtain $n_{s}^{\prime}\simeq-0.004$ which is consistent with that of
\emph{Planck 2015}, namely $n_{s}^{\prime}=-0.003\pm0.007$. Moreover, we find
that the running of the scalar spectral index does not change significantly as
a function of $n_{s}$ (see figure 2). Regarding the inflaton field at the
beginning of inflation, we obtain $10.6M_{pl}\leq\phi\leq12.9M_{pl}$; while at
the end of inflation we require $0.7M_{pl}\leq\phi_{end}\leq1.1M_{pl}$.
Furthermore, we find that the allowed region in which our results satisfy the
$2\sigma$ observational restrictions of \emph{Planck 2015} is $f\geq
11.7M_{pl}$ and $1<b\leq1.5$.
Thus the data applied to slow-roll inflation place constraints on $b$,
although the value of scale $f$ has only a lower limit at $2\sigma$ level.

Finally, we would like to compare our results with those found using other
potentials. Specifically, in the case of the chaotic inflation $V(\phi
)=\Lambda^{4}(\phi/M_{pl})^{k}$ \cite{Linde}, the corresponding slow-roll
parameters are written as $\epsilon=k/4N$, $\eta=(k-1)/2N$ and $\xi
=(k-1)(k-2)/4N^{2}$. The latter implies $n_{s}=1-(k+2)/2N$, $r=4k/N$ and
$n_{s}^{\prime}=-(2+k)/2N^{2}$. Using $k=2$ and $N=55$ we obtain $n_{s
\simeq0.964$, $r\simeq0.145$ and $n_{s}^{\prime}\simeq-0.0007$. This also
corresponds to the slow-roll regime of intermediate inflation \cite{int} with
Hubble rate during inflation given by $H\propto t^{k/(4-k)\text{ }}$with
$n_{s}=1-(k+2)r/8k$ and $k=-2$ gives $n_{s}=1$ exactly to first order and
$n_{s}^{^{\prime}}=-2(n_{s}-1)^{2}/(k+2)$. On the other hand the Starobinsky
inflation \cite{staro}, namely $V(\phi)\propto\lbrack1-2\mathrm{e
^{-B\phi/M_{pl}}+\mathcal{O}(\mathrm{e}^{-2B\phi/M_{pl}})]$
, leads to the following slow-roll predictions \cite{Muk81,Ellis13}:
$n_{s}\approx1-2/N$ and $r\approx8/B^{2}N^{2}$, where $B^{2}=2/3$.
Furthermore, for Starobinsky inflation following the notations of
\cite{Ellis13} we find that the running spectral index is given by
$n_{s}^{\prime}\approx-2/N^{2}$. To this end using $N=55$ we obtain
$(n_{s},r,n_{s}^{\prime})\approx(0.963,0.004,-6.6\times10^{-4})$. The above
slow-roll parameters are indicated by the solid points (chaotic inflation) in
figures 1 and 2 while the stars represent values for the Starobinsky
inflation. \emph{\ }

\section{Conclusions}

In the light of the \emph{Planck 2015} results, a debate is taking place in
the literature about the best implementation of the inflationary paradigm. In
the current article we would like to contribute to this discussion. Employing
the hyperbolic dark-energy scalar-field potential of Rubano and Barrow
\cite{Barrow} (see also \cite{Saa}) we study the performance of this model as
a description of inflation. We find that the hyperbolic inflation turns out to
be quite promising in the context of the new data from \emph{Planck 2015}.
Specifically, using the scalar-field potential (\ref{pot1}), we calculate the
slow-roll parameters analytically and then we compare the corresponding
predictions against the observational data. We find that currently hyperbolic
inflation is consistent with the results provided by \emph{Planck 2015} within
$1\sigma$ uncertainties. The combination of ref \cite{Barrow} with
these calculations provide an overall
cosmological investigation of the potential given by Eq.(\ref{pot1}). We find
that the hyperbolic structure of this potential leads to a viable model which
can be used separately to understand the main properties of both inflation and
dark energy in the presence of a single perfect fluid\footnote{If two fluids
are present (for example, matter and radiation) simultaneously then an exact
solution requires the scalar field to have a more complicated potential than
the power of an exponential. Our model gives the limiting form of that more
complicated solution in the limiting cases where one perfect fluid dominates
over the other.}. Finally, following the notations of \cite{Peea} we
can provide a unified picture of dark energy and inflation using the sum of
the potentials (\ref{pot}) and (\ref{pot1})
\[
V_{tot}(\phi)=A_{de}\;\sinh^{b_{de}}\left[  \lambda_{de}\left(  \phi-\phi
_{0}\right)  \right]  +A\;\mathrm{sinh}^{b}(\phi/f)
\]
where $\lambda_{de}=\frac{\sqrt{3}(1-\gamma_{\phi})}{\sqrt{\gamma_{\phi}}}$ and the
constants $A_{de}$, $b_{de}$, $b$, are given by Eqs.(\ref{AA}), (\ref{bb}) and
(\ref{bbw}). Note that similar considerations hold for the early dark-energy
model \cite{Bar00}.

\textbf{Acknowledgments.} SB acknowledges support by the Research Center for
Astronomy of the Academy of Athens in the context of the program
\textquotedblleft\textit{Tracing the Cosmic Acceleration}\textquotedblright.
JDB acknowledges STFC support. Also, we would like to thank N. E. Mavromatos
for useful comments and suggestions.


\section{Introduction}Software Product Line Engineering (SPLE) is a paradigm used to develop software-intensive systems that share common assets \cite{pohl_software_2005, apel_feature-oriented_2013}. In SPLE, a feature is used to represent a characteristic behavior as a unit of functionality of the software product line (SPL) \cite{apel_feature-oriented_2013}. Given a set of features, the configuration of a SPL is constructed by composing the features. Following generative programming mechanisms \cite{barth_generative_2001}, the respective product variant can be derived automatically from a feature selection and  artifacts of each selected feature. The product generation may be realized through Feature-Oriented Programming (FOP) techniques \cite{batory_scaling_2003}, in which each feature is mapped to a feature module containing its artifacts. The product variant is synthesized from the artifacts of the involved feature modules. In this context, we demonstrate our approach that helps the developer to write feature modules containing their artifacts using a dedicated and generic language.

Some authors proposed  approaches for constructing product variants of a product line by reusing and synthesizing artifacts. Using design by contract \cite{meyer_applying_1992} and adhering to FOP techniques, Th\"um proposed a proof composition approach and strategies to reduce the effort of verification  by reusing partial proofs \cite{thum_proof_2011} \cite{thum_family-based_2012}  \cite{thum_classification_2014}. The FEATUREHOUSE framework, described in \cite{apel_featurehouse:_2009}, enables the construction of a new program from existing programs using the structure of the SPL feature tree. Although the above methods can implement software product line efficiently and mostly automatically, we still need a method to generate and prove correctness of all product variants more efficiently and automatically.

In this direction some advances have been proposed in the context of programming language meta-theory within the Coq proof assistant \cite{delaware_product_2011}  \cite{delaware_modular_2013}.
However these tools  are dedicated to a very specific domain and still require an important expertise.

To tackle the above limitations, we depict a methodology that helps the user to automatically produce correct-by-construction product variants. This methodology aims at being independent of any concrete target language but at first we focus and experiment with
\focal{}i{Z}e, a language used to specify, implement and prove (\url{http://focalize.inria.fr}). In this paper, we propose a language, called \textbf{GFML} (for Generic Feature Module Language) allowing the developer to write feature modules which contain three kinds of artifacts:  specifications, code and correctness proofs. GFML is inspired from \focal{}i{Z}e{} but sticks to a FOP approach and tries to reduce the developer's effort. This language is designed so that the artifacts can be easily reused, composed and translated into different languages.

In Section \ref{background}, we first describe the background of our paper, in particular we present  briefly \focal{}i{Z}e. Then in Section \ref{Artifacts}, we illustrate our definition of feature module and a brief description of our language (GFML) used to write feature modules. In Section \ref{sec:GFMLtranslation}, we illustrate the translation of GFML into \focal{}i{Z}e. In Section \ref{sec:composition}, we describe how to compose two feature modules by giving an example, and  depict a composition framework for automatically generating correct-by-construction product variants. In Section \ref{sec:Relate}, we review related work. Finally we conclude the paper and discuss future works in Section \ref{sec:conclusion}.

\section{Background}\label{background}
\subsection{Software Product Line}
In SPLE, a \textbf{feature} is a characteristic behavior specified as a unit of functionality of a product line \cite{apel_feature-oriented_2013}. A set of features, called a \textbf{feature model}, is often graphically depicted as a tree, also called a \textbf{feature diagram}. In addition to the presentation of commonality and variability, a feature diagram also contains various relationships and additional constraints between  features. The main relationships in feature diagrams are optional/mandatory, alternative/or and implies/excludes. Because of visual aid, feature diagrams are considered as standard representations that help the user to choose product configurations of a product line. In this work, we only focus on a simple form of feature diagrams (i.e. Figure \ref{fig:ex-dia}) that can illustrate optional features.

\begin{wrapfigure}{r}{0.35\textwidth}
  \vspace{-15pt}
  \begin{center}
    \includegraphics[width=0.35\textwidth]{ex-dia}
  \end{center}
  \vspace{-15pt}
  \caption{Feature diagram}
  \label{fig:ex-dia}
  \vspace{-10pt}
\end{wrapfigure}

\textbf{Valid Configuration.} A configuration of a product line is a set of features selected in the corresponding feature diagram. All the information such as relations and constraints inferred from the feature diagram is used to check the validity of the configuration. In our work, we only consider valid configurations. However, checking the validity of a configuration is out of the scope of this paper. Many checkers exist (e.g. \cite{batory_feature_2005}, \cite{BenavidesT13}) and we rely on them. As we start from valid configurations, we consider that the configurations contain all the features that are to be composed in order to build the final products. So the distinction between mandatory and optional features is not relevant here.

\subsection{Motivating Example}
\begin{wrapfigure}{r}{0.35\textwidth}
  \vspace{-40pt}
  \begin{center}
    \includegraphics[width=0.35\textwidth, height=0.16\textwidth]{BA-dia}
  \end{center}
  \vspace{-15pt}
  \caption{Feature diagram of bank account product line}
  \vspace{-50pt}
  \label{fig:BA-dia}
\end{wrapfigure}

As a running example, we consider the bank account product line described in \cite{thum_family-based_2012} whose feature diagram is shown in Figure \ref{fig:BA-dia}. It illustrates a family of products allowing the management of bank accounts. The root feature \textit{BankAccount} (BA for short) provides the basic management of an account. It allows the bank storing the current balance and the amount of money added into or withdrawn from the account. A customer can withdraw more money from the account than available if it is within an over limit. The feature BA has three optional child features \textit{DailyLimit} (DL for short), \textit{LowLimit} (LL for short), \textit{Currency}.
The feature DL allows the bank to limit the amount of money withdrawn in a day while the other feature LL indicates that the bank only authorizes a customer to withdraw money from the account only if the amount is greater than a low limit. Related
The feature \textit{Currency} accommodates the management of currency. Finally, an optional feature \textit{CurrencyExchange} is established as a child of \textit{Currency} to enable the calculation of currency exchange.

\subsection{Quick Presentation of \focal{}i{Z}e}\label{focalize}

This subsection presents briefly the technical background necessary to understand the section about the translation of GFML to \focal{}i{Z}e.

The \focal{}i{Z}e{} (\url{http://focalize.inria.fr}) environment provides a set of tools to describe and implement
functions and logical statements together with their proof. A \focal{}i{Z}e{}
source program is analyzed and translated into OCaml sources for execution and
Coq  sources for certification. The \focal{}i{Z}e{} language has an object oriented
flavor allowing inheritance, late binding and redefinition. These characteristics are very helpful to reuse specifications, implementations and proofs.

A \focal{}i{Z}e{} specification can be seen as  a set of algebraic properties describing relations between input and output of the functions implemented in a \focal{}i{Z}e{} program.
For writing code,  \focal{}i{Z}e{} offers a pure functional programming style close to ML, featuring strong typing, recursive functions, data types and pattern-matching.
Proofs written using the \focal{}i{Z}e{} proof language
are sent to the Zenon prover which produces proofs that can be verified by Coq for more confidence \cite{bonichon_zenon_2007}. The \focal{}i{Z}e{}
proof language is a declarative language in which the programmer states a
property and gives hints  to achieve its proof which is performed by Zenon.

\focal{}i{Z}e{} units are called \emph{collections}. They  contain
 entities in a model akin to classes and objects or types and values. Collections have
 functions and  properties which can be called using the ``\texttt{!}''  notation. They are derived from other units called \emph{species} which specify and implement functions.

A species defines a set of entities together with functions and
properties applying to them. At the beginning of a development,  the
representation  of these entities is usually  abstract, it is
precised later in the development. However the type of these entities
is referred as \texttt{Self} in any species. Species may contain specifications, functions and proofs.
More precisely species may specify a function or a property (with resp. \texttt{signature},
\texttt{property} keywords) or implement them
(\texttt{let} keyword when a function is defined, \texttt{proof of} keyword to introduce a proof of a property).
A \texttt{let} defined function must match its signature and similarly a proof
introduced by \texttt{proof of} should prove the statement given by the
\texttt{property} keyword.  Statements belong to first order typed logic.

As said previously, \focal{}i{Z}e{} integrates inheritance, late binding and redefinition to ease reuse and modularity. Inheritance allows the definition of a new species from one or several other species. The new species inherits all the functions, properties and proofs of its parents. Some syntactical mechanisms  are provided to prevent ambiguities. A species may provide a definition for a function that is only specified in its parents. It may also redefine a function when this one is already defined in a parent but in that case the signature is maintained (no overloading). Multiple inheritance comes with a late binding mechanism close to the one found in object oriented languages (even if the resolution is statically  done by the compiler).

A collection is a species where every specified function is defined and every property is proved (or explicitly admitted). Furthermore, in a collection, the concrete \emph{representation} of entities is made private and a programmer using the collection can only use its functions and properties according to the interface. Consequently building a collection from a species provides encapsulation.
Species may have parameters which may either be collections or entities providing
in that way parametric polymorphism.  This parametrization will be intensively used in the translation of GFML to \focal{}i{Z}e{} (see Section~\ref{sec:GFMLtranslation}).

For more details on \focal{}i{Z}e{} please refer to the reference
manual. A more thorough overview  can be found in
\cite{hardin2008}.

\section{Artifacts}\label{Artifacts}
\subsection{Feature Module}
We adhere to the FOP technique that suggests to map each feature to a separate \textbf{feature module} that implements the feature. In our context, each feature module consists of its \textbf{artifacts}: specifications, code and correctness proofs. Specifications are given as a set of expected properties or requirements. Technically these properties are logical formulas relating together some functions described only by their signature. Thus in our setting, a specification is close to an algebraic data-type. Code has to be understood here  as the implementation of the functions introduced/declared in the specifications. We place ourselves in a functional programming setting. The proofs here concern the correctness of the code with respect to the specifications.

\begin{itemize}
\item \textbf{Specification artifact} includes the function declarations  and the properties. A function declaration or signature only describes the name of the function  and the type of its arguments and result. A property is a (first order typed) logical formula. A new property is expressed by a new logical formula while a refining property refers to  existing properties of the $parent$ feature module and may add new premises and conclusions.
\item \textbf{Code artifact} consists of the definition of the concrete representation type $r$  and the function definition/redefinitions $df$. The representation type is the concrete type associated to the abstract data-type specified in the specification artifact. It will be a concrete type (basic or complex types, \textit{à la ML}) or a Cartesian product of concrete types. The representation type is  unique for each feature module and can be also constructed from the existing representation type of $parent$. A function is defined/redefined after the representation type has been defined.
\item \textbf{Proof artifact} contains proofs $pf$. Each proof corresponds to  a property. While writing the proof for a refining property, the mentioned properties of $parent$ can be reused to support the proving process.
\end{itemize}

We define a feature module $fm$ as a 5-tuple: function declarations $d$, properties $p$, representation type $r$, function definitions/redefinitions $df$ and proofs $pf$.
\begin{center}
$fm = (d, p, r, df, pf)$
\end{center}

\subsection{Description of GFML}

We propose a generic language, called \textbf{GFML}, to write these feature modules. Each feature module is embedded into a separate file \textit{.gfm}. This language is suitable for all feature modules possibly containing three kinds of artifacts: specifications, code and correctness proofs.

This language is inspired from \focal{}i{Z}e{}, in particular the styles for writing specifications, code and proofs are common. GFML and  \focal{}i{Z}e{} mainly differ in the way to structure and organize information. Our main objective in defining GFML is to propose a  language close to \focal{}i{Z}e{} that already allows for the three kinds of artifacts in a single setting but closer to the description of commonality and variability we can find in feature diagrams. As we will see in Section~\ref{sec:GFMLtranslation}, the expression of a feature module is much simpler that its translation in \focal{}i{Z}e.
GFML is also inspired from design by contract applied to FOP as in
\cite{ThumSKAS12}. In the same way, GFML syntax allows the programmer to focus in a GFML feature module on the modifications brought to the specifications or implementations of the parent.

\subsection{Examples}
Corresponding to the feature diagram of the bank account product line given in Figure \ref{fig:BA-dia}, the feature module BA implements the root feature BA. Feature modules DL, LL, and \textit{Currency}
are mapped from the features DL, LL, and \textit{Currency}
respectively. They have a common parent, i.e. the root module BA. Feature module \textit{CurrencyExchange} corresponding to feature \textit{CurrencyExchange}, is a child of the module \textit{Currency}.

\begin{figure}
        \begin{subfigure}[b]{0.475\textwidth}
          \begin{center}
                \fbox{\includegraphics[width=1\textwidth]{full_BA}}
                \caption{Module BA}
                \label{fig:BA}
                \end{center}
        \end{subfigure}\qquad
        \begin{subfigure}[b]{0.475\textwidth}
                \fbox{\includegraphics[width=1\textwidth]{full_DL}}
                \caption{Module DL}
                \label{fig:DL}
        \end{subfigure}
         \vspace{-20pt}
        \caption{Feature modules in GFML}\label{fig:BA+DL}
          \vspace{-10pt}
\end{figure}

The root feature module BA, written with GFML, is shown in Figure \ref{fig:BA}. This module includes three signatures: $over$ - is over limit, $balance$ - gets the current balance value  of the account and $update$ - upgrades the new value balance and also returns the status of the operation (\textit{success} or \textit{no\_success}, of type $S$ whose definition is omitted here). The next part of the module BA contains the property $update\_succ\_BA$ (line 5) that specifies a customer can withdraw more money $a$ from the account than available if the balance is within $over$. In that case, its status must be $success$.
The primitive functions $first$ and $second$ used in this property are the usual projections of a Cartesian product.
Then the representation type is defined as $int$ (line 10), it means that an account is only represented by its balance. Then appear the definitions of the functions $balance$ and $update$ (lines 11-15). The proof of property $update\_succ\_BA$ includes two proof hints: by definition of $update$ and $balance$ definitions. This means that the proof must be done by unfolding the definitions of both functions.

Another example of GFML feature module is the module DL defined according to its $parent$ feature module BA (Figure \ref{fig:DL}). Two new declarations $limit\_withdraw$ and  $withdraw$ are added into the module (lines 2-3).  The module introduces the constant $limit\_withdraw$ only declared at that point.  It denotes the limit of withdrawn money in a day. It also introduces  another function  $withdraw$  that returns, for an account, the current amount of withdrawn money in a day. The functions $update$, $balance$ and $over$ defined in the parent are also available in the present feature module. A new property $update\_succ\_DL$ is obtained by modifying the property $update\_succ\_BA$ from the $parent$ (line 6). This modification includes a new premise called H1  (line 7) for short in the figure (but announced below)  and a new conclusion (lines 8-9). The premise H1 is expressed as follows:
\begin{align*}
&r = update(x,a) \rightarrow (a <= 0) \rightarrow
(all \enspace n\_w \enspace w : int, all \enspace n\_s : S, \enspace withdraw (x) = w \enspace \&\&\enspace w + a = \\
&n\_w \enspace \&\& \enspace second(BA!update(x,a)) = n\_s \rightarrow (n\_w >= limit\_withdraw) \enspace \&\&  \enspace (n\_s = success))
\end{align*}
It states that the bank allows a customer to withdraw money only if the amount of withdrawn money in a day is greater than $limit\_withdraw$ (in this case \textit{w}, \textit{n\_w} and $limit\_withdraw$ are negative numbers). The new conclusion states that this operation has to modify the account by updating the amount of withdrawn money. The representation type of module DL is defined as a Cartesian product of $int$ (i.e. the concrete type associated with the amount of money withdrawn in a day) and the representation type of the \textit{parent} (line 11). The functions $withdraw$ and $update$ are defined/redefined (lines 12-23). We can notice that in the redefinition of $update$ the call to the parent function $update$ is done with the parameter $x$, there is here an implicit conversion that will be inserted during the translation into \focal{}i{Z}e.
The proof of the property $update\_succ\_DL$ reuses the property $update\_succ\_BA$ as a proof hint (line 26). All modifications and additions compared to the parent module BA are highlighted or underlined.

\section{From GFML to \focal{}i{Z}e}\label{sec:GFMLtranslation}
In this section, we  illustrate the translation
from GFML to \focal{}i{Z}e. A GFML feature module is mapped into three \focal{}i{Z}e{} separate \textit{species}. The first \textit{species} \textit{\textbf{Inh\_}} contains all desired elements that are inherited in the child module. These elements can be function declarations, properties or predicates. A predicate, which is introduced with the keyword `$logical\ let$' and named $l\_cons$, is associated with each property $p$ found in a module.
For example, the feature module BA written in Figure \ref{fig:BA_foc} is mapped to the species $Inh\_BA$ (line 1). A predicate $l\_cons\_update\_succ\_BA$ (line 5) is associated with the property $update\_succ\_BA$. In Figure \ref{fig:DL_foc}, the species $Inh\_DL$ of the child DL inherits the species $Inh\_BA$ of module BA (line 2).

\begin{figure}
        \begin{subfigure}[b]{0.475\textwidth}
          \begin{center}
                \fbox{\includegraphics[width=1\textwidth]{full_BA_foc}}
                \caption{Module BA in FoCaLize}
                \label{fig:BA_foc}
                \end{center}
        \end{subfigure}\qquad
        \begin{subfigure}[b]{0.475\textwidth}
                \fbox{\includegraphics[width=1\textwidth]{full_DL_foc}}
                \caption{Module DL in FoCaLize}
                \label{fig:DL_foc}
        \end{subfigure}
         \vspace{-20pt}
        \caption{Feature modules in FoCaLize}\label{fig:BA+DL_foc}
          \vspace{-10pt}
\end{figure}

The \focal{}i{Z}e{} inheritance mechanism allows the programmer  to reuse the artifacts without any change (except for functions that can be redefined without changing its signature). But although some properties may be kept in the child module, others may be modified by adding new premises or/and new conclusions. In this last case, the \focal{}i{Z}e{} inheritance mechanism is not sufficient. We have to use the parametrization facility to encode the right meaning. To tackle this limitation, we propose the second \textit{species} \textbf{\textit{Reu\_}} that inherits the first \textit{species} \textit{\textbf{Inh\_}}. It includes all desired specifications that are reused to describe new ones in the child. Each property specified in this \textit{species} must have a corresponding predicate  $l\_cons$ in the first \textit{species} \textit{\textbf{Inh\_}}. Instead of mentioning the property, we use this predicate for expressing new properties in the child. This trick is used because \focal{}i{Z}e{} does not  allow the modification of an inherited  property in a \textit{species}. For example, in Figure \ref{fig:BA_foc} \textit{species} $Reu\_BA$ inherits the \textit{species} $Inh\_BA$. It contains the property $update\_succ\_BA$ (line 13) that is related to the predicate $l\_cons\_update\_succ\_BA$ in \textit{species} $Inh\_BA$. This predicate is reused to specify the refining property $update\_succ\_DL$ of module DL (line 10 of Figure \ref{fig:DL_foc}).

 The parametrization mechanism in FoCaLize
is used here to circumvent the fact that in \focal{}i{Z}e{} when the representation type is fixed in a species P, it cannot be changed in any species which inherits  P. So in order to build a new species that can reuse functions, properties, predicates and proofs that are defined for BA, we have to parametrize this species by an implementation of BA.  Thus,  when the DL feature is selected, the account is implemented as a pair containing a basic account of type BA and an integer corresponding to the new attribute, i.e. the amount of money withdrawn in a day. The implementation of the functions $balance$ and $over$ are, in this context, just a call to the parent's functions combined with the first projection. For example, \textit{species} $Imp\_BA$ inherits \textit{species} $Reu\_BA$ (line 18 of Figure \ref{fig:BA_foc}) while  \textit{species} $Imp\_DL$ inherits \textit{species} $Reu\_DL$ (line 14 of Figure \ref{fig:DL_foc}). A parameter $BA$ encapsulates \textit{species} $Reu\_BA$ (line 13 of Figure \ref{fig:DL_foc}). It is used to call function $update$ of module BA (line 18).
In the proof of the refining property $update\_succ\_DL$, property $update\_succ\_BA$ of the parent BA is considered as a proof hint (line 30). Its associated formula $l\_cons\_update\_succ\_BA$ is also a proof hint (line 29). Let us notice that Zenon, the \focal{}i{Z}e{} prover has done all the proofs with the given hints.

We can notice that the use of a language such as \focal{}i{Z}e{} for implementing and verifying feature modules is quite complex. The proposed language GFML is a solution to write the feature modules together with their artifacts more easily.
We have followed with success this translation scheme for all the features appearing on the feature diagram of Figure \ref{fig:BA-dia}. For the moment, the translation is done manually.

\section{Towards Feature Module Composition}\label{sec:composition}
\subsection{Methodology for composing artifacts}
\begin{figure}
         \centering
         \includegraphics[height=4.3cm, width=16cm]{met.png}
         \vspace{-20pt}
         \caption{Methodology}
         \vspace{-10pt}
         \label{fig:met}
      \end{figure}

The work previously described is the first step in the proposition of a methodology to
help the user to automatically produce correct-by-construction product variants from features selected in a feature diagram.  This methodology is illustrated in Figure \ref{fig:met}. We assume a SPL is described with a feature diagram. First,  the developer writes the different feature modules with the GFML language. He uses the \textbf{GFM Pre-compiler} to translate his/her feature modules into \focal{}i{Z}e{} to verify them (using the Zenon prover). Once this work has been done, the user - he may be different from the developer - chooses  features from the feature diagram. Based on this selection, the corresponding configuration is determined. The related GFML feature modules are then sent to the \textbf{GFM Combiner} that will compose the feature modules of the configuration in order to obtain a composition feature module that will be translated into the desired product variant, which by composition is correct by construction. We expect this method to be independent from the target language which is required to be able to express specifications and code and also proofs (unless these ones can be automatically produced).

In the next subsection we explain the main ideas of our composition mechanism on the running example.

\subsection{Composition of two feature modules: an example}
Some problems may appear when composing two feature modules, in particular  conflicts may appear when synthesizing their artifacts. For example, which properties should be performed first? Which   synthesized representation type will make less complex writing? To solve these problems, we suggest the user gives a composition order of the modules. In our work, we only consider the modified elements such as refining properties, representation types, function redefinitions or proofs of refining properties.

\begin{figure}
        \begin{subfigure}[b]{0.475\textwidth}
          \begin{center}
                \fbox{\includegraphics[width=1\textwidth]{full_LL}}
                \caption{Module LL}
                \label{fig:LL}
                \end{center}
        \end{subfigure}\qquad
        \begin{subfigure}[b]{0.475\textwidth}
                \fbox{\includegraphics[width=1\textwidth]{full_LL_DL}}
                \caption{Module LL\_DL}
                \label{fig:LL_DL}
        \end{subfigure}
         \vspace{-20pt}
        \caption{Composition of feature modules}\label{fig:com}
          \vspace{-10pt}
\end{figure}

Similar to module DL, module LL written in Figure \ref{fig:LL} is extended from the root module BA. A new signature $limit\_low$ declares the low limit of the account. A new property $update\_no\_succ\_LL$ is specified in line 7 and its proof is indicated in line 16. However, we only consider the modifications that can cause conflicts when composing two modules. For example, let us compose  module LL with module DL. Their composition is the module LL\_DL given in Figure \ref{fig:LL_DL}. Property $update\_succ\_LL$ of module LL includes two parts (lines 3-6 of Figure \ref{fig:LL}). Its highlighted part includes a new premise (shortened by H2), expressed as follows:
\begin{align*}
r = update(x,a) \rightarrow
	((a >= 0) \enspace || \enspace (a <= limit\_low))
\end{align*}
The premise specifies that a customer can withdraw money $a$ from the account only if $a$ is less than $limit\_low$ (in this case \textit{a} and $limit\_low$ are negative numbers). The property $update\_succ\_BA$ from the $parent$ BA is its remaining part. Following the composition order, the composite property of two properties $update\_succ\_LL$ and $update\_succ\_DL$ is synthesized by taking H2 first and then mixing with $update\_succ\_DL$. It is named $update\_succ\_LL\_DL$ and  embedded into a composition feature module LL\_DL (lines 2-6 of Figure \ref{fig:LL_DL}). Similarly, the composite representation type is synthesized in line 7. The function $update$ is redefined (lines 8-11). The proof of the composite property $update\_succ\_LL\_DL$ reuses two properties $update\_succ\_DL$ and $update\_succ\_BA$ as its proof hints (lines 13-14) indicate. The highlighted parts of module LL (Figure \ref{fig:LL}) are the parts highlighted in module LL\_DL (Figure \ref{fig:LL_DL}). The other parts which are not highlighted in Figure \ref{fig:LL_DL} are taken from module DL.

\section{Related Work} \label{sec:Relate}
Recently, many authors focused on constructing and verifying program variants in SPL. Together with these, composition methods of feature artifacts are offered. In this subsection, we compare to some of them.

Th\"um et al., in \cite{thum_proof_2011}, presents a composition approach where partial proofs are given in features and then composed to build the complete proofs for an individual product. The specifications are expressed using design by contract \cite{meyer_applying_1992}. The Krakatoa/Why tool \cite{MarchePU04} is used to generate proof obligations that are exported to the Coq proof assistant where proofs are done and verified.

Our work is in line with the previous approach but it is dedicated to a functional setting whereas Th\"um et al's work involves object-oriented programs. Specification artifacts - expressed as logical formulas - found in a GFML feature module are close to contracts. For example, a refining property
is close to a refining contract since it includes a former specification and may add a new premise and a new conclusion. However in our work, each proof is complete. It is extended in a new proof in the result of a composition and thus in the resulting product variant.

Another composition framework FEATUREHOUSE is  offered in \cite{apel_featurehouse:_2009}.
Using the FST (feature structure tree) model, existing artifacts can be composed to construct a new program. The artifacts must have tree structures. In contrast, we presented a framework that allows us to construct products from the artifacts of the features selected by the user. We are interested in expressing feature compositions as  algebraic expressions using composition operations for each kind of artifacts instead of relying on the tree structure. Similar ideas were mentioned in \cite{batory_scaling_2003} \cite{batory_feature_2011}. However, these approaches only focus on constructing products but do not mention how to verify them.

Recent researches \cite{delaware_product_2011}  \cite{delaware_modular_2013}  had proposed some advances in feature composition in the context of meta-theory. However, these tools are dedicated to a very specific domain. Similar to Thum's work, the products in these approaches are verified in Coq. In contrast, implementing feature modules independently of any concrete target language is the purpose of our work. By proposing a generic formal language (GFML) for feature module, the artifacts are easy to reuse and synthesize but don't belong to any concrete language.

\section{Conclusion}\label{sec:conclusion}
In this paper we have described a first step towards the production of a methodology allowing for the development of correct-by-construction product variants according to a FOP paradigm. The contribution is here the definition of the language GFML allowing the developer to write feature modules containing their specifications, code and correctness proofs.

These modules are translated to \focal{}i{Z}e{} for verification and also for obtaining OCaml operational code. Consequently the product variants we are aiming at will be implemented in OCaml.
GFML is here introduced to help the developer when he is writing the different feature modules, allowing him to describe the artifacts of a feature with respect to the artifacts of its parents. Realizing this directly in \focal{}i{Z}e{} is a difficult task as it is exemplified by the translation scheme presented in a previous section. For the moment, nothing is implemented, all the examples found in this paper (all the possible configurations of the case study) have been obtained by a manual translation and composition.

Next step is to provide an implementation for the GFML Pre-compiler that will automatically translate GFML to \focal{}i{Z}e. Then composition of feature modules  will be formally defined and  implemented in the GFML Combiner.

Another important perspective is to make our methodology independent of the target languages. Regarding this point, an intermediate step could be to adapt our methodology and tools to B or EventB where inheritance is not available.

\bibliographystyle{eptcs}


\section{Introduction}
Let $\Omega \subset \mathbb{R}^{2m}$ be a bounded domain. In this context we study the existence of multiple solutions in $W^{m,2}_{\mathcal{N}}(\Om) = \{u \in W^{m,2}(\Omega): \Delta^j u = 0 \mbox{ on } \partial \Omega \text{ for } 0 \leq j < \frac{m}{2}\}$ of the following $2m$-th order problem with Navier boundary condition
$$(P)\hspace{1cm}\left\{
\begin{array}{lllll}\left.
\begin{array}{lllll}
(-\Delta)^m u &= \mu u |u|^pe^{u^2} + \lambda h(x)\\
\hspace{10mm} u &> 0 \end{array}\right\} \;\; \text{in} \;\Om, \\
u=\De u=0=..=\De^{m-1}u \hspace{15mm}\text{ on } \partial \Omega \end{array}\right.
$$
 where $h \geq 0$ in $\Omega$, $\|h\|_{L^2(\Omega)} = 1,$ $ \lambda >0, \mu =1$ if $p>0$ and $\mu \in (0,\lambda_1(\Omega))$ if $p=0$. Also assume that $\lambda_1(\Omega)$ be the first eigenvalue of $(-\Delta)^m $ on $W^{m,2}_{\mathcal{N}}(\Omega)$ with Navier boundary condition and which is strictly positive. The existence of multiple solutions for analogous problems in higher dimension with critical exponent have been studied in \cite{Lu}, \cite{Peral} for the Dirichlet boundary condition and in \cite{Zhang} for Navier boundary condition. The corresponding problem for second order elliptic equations have been studied in \cite{Prashanth} for dimension two, and in \cite{Tarantello} for higher dimensions. The critical growth for the nonlinearity is $u \mapsto |u|^pu$, $p = \frac{4m}{n-2m}$, when $n \geq 2m+1$
 from the Sobolev embedding in $\mathbb{R}^n$. In \cite{Moser} Moser proved the following,
 \begin{thm}          \label{Moser_Thm}
 Let $\Om \subset \mathbb{R}^n, n \geq 2$ be a bounded domain. There exists a constant $C_n >0$ such that for any $u \in W^{1,n}_0(\Omega)$, $n\geq 2$ with $\|\nabla u\|_{L^n(\Omega)} \leq 1$, then
 \begin{equation}
  \int_{\Omega} e^{\alpha |u|^p}dx \leq C_n|\Omega|,              \label{Moser_Sharp}
 \end{equation}
where $$p = \frac{n}{n-1}, \hspace{1mm}\alpha_n:= n w^{\frac{1}{n-1}}_{n-1},$$ and $w_{n-1}$ is the surface measure of the unit sphere $\mathbb{S}^{n-1} \subset \mathbb{R}^n$. Furthermore the integral on the left hand side can be made arbitrarily large if $\alpha > \alpha_n$ by appropriate choice of $u$ with $\|\nabla u\|_{L^n(\Omega)} \leq 1$. The embedding $$W^{1,n}_0(\Omega) \ni u \mapsto e^{\alpha |u|^{\frac{n}{n-1}}} \in L^1(\Omega)$$ is compact for $\alpha < \alpha_n$ and is not compact for $\alpha = \alpha_n$.
\end{thm}
In \cite{Adams} Adams extended the above result of Moser to higher order Sobolev spaces. To state the result of Adams we define the following $m$-th order derivatives of $u \in C^m(\Omega)$:
\[ \nabla^m u = \left\{
            \begin{array}{lllll}
             \Delta ^{\frac{m}{2}} u    \hspace{1cm} \text{for } m \text{ even},\\
             \nabla \Delta^{\frac{m-1}{2}} u \hspace{3mm} \text{for } m \text{ odd.}
            \end{array}\right.\]
            Furthermore, $\| \nabla^m u\|_p$ is the $L^p$ norm of the function $|\nabla^mu|,$ the usual Euclidean length of the vector $\nabla^m u$. We also denote $W^{m,\frac{n}{m}}_0(\Omega)$ to be the completion of $C^{\infty}_0(\Omega)$ under the Sobolev norm
            \begin{equation}
            \|u\|_{W^{m,\frac{n}{m}}(\Omega)}= \left( \|u\|^{\frac{n}{m}}_{\frac{n}{m}} + \sum_{|\alpha|=1}^m \|D^{\alpha} u\|^{\frac{n}{m}}_{\frac{n}{m}}\right)^{\frac{m}{n}}.              \label{full-norm}
            \end{equation}
 \par Adams proved the following embedding:
 \begin{thm}\label{Adams_Thm}
  Let $\Omega \subset \mathbb{R}^n$ be a bounded domain. If $m$ is a positive integer less than $n$, then there exists a constant $C_0=C(n,m)>0$ such that for any $u \in W^{m,\frac{n}{m}}_0(\Omega)$ with $\| \nabla^m u\|_{\frac{n}{m}} \leq 1,$ then
 \begin{equation}
  \frac{1}{|\Omega|} \int_{\Omega} \exp(\beta |u(x)|^{\frac{n}{n-m}})dx \leq C_0,              \label{Adams_Sharp}
  \end{equation}
   for all $\beta \leq \beta_{n,m}$ where
 \[\beta_{n,m}=\left\{
\begin{array}{lllll}
\frac{n}{w_{n-1}}\left[\frac{\pi ^{\frac{n}{2}}2^m \Gamma(\frac{m+1}{2})}{\Gamma(\frac{n-m+1}{2})}\right]^{\frac{n}{n-m}} \hspace{3mm} \text{when } m \text{ is odd,}\\
\frac{n}{w_{n-1}}\left[\frac{\pi ^{\frac{n}{2}}2^m \Gamma(\frac{m}{2})}{\Gamma(\frac{n-m}{2})}\right]^{\frac{n}{n-m}} \hspace{3mm} \text{when } m \text{ is even.}
\end{array}\right.\]
Furthermore, for any $\beta > \beta_{n,m}$, the integral can be made as large as possible by appropriate choice of $u$ with $\|\nabla ^m u\|_{\frac{n}{m}} \leq 1$.
\end{thm}
\begin{Rem}
We remark that for the case $n=2m=4$, Lu-Yang in \cite{Lu-Yang} and in general Zhao-Chang \cite{Zhao-Chang} showed the existence of an explicit sequence for $n=2m$ to prove the sharpness of the constant in $W^{m,\frac{n}{m}}_0(\Omega)$.
\end{Rem}

Let $W^{m,\frac{n}{m}}_{\mathcal{N}}(\Omega)$ denote the following subspace of $W^{m,\frac{n}{m}}(\Omega)$:
$$W^{m,\frac{n}{m}}_{\mathcal{N}}(\Omega) = \left\{u \in W^{m,\frac{n}{m}}(\Omega): \Delta^j u =0, \text{ on } \partial \Omega \text{ for } 0 \leq j\leq [ (m-1)/2]\right\}.$$ Note that $W_0^{m,\frac{n}{m}}(\Omega)$ is strictly contained in $W^{m,\frac{n}{m}}_{\mathcal{N}}(\Omega)$.
Therefore,
 \[ \sup_{u \in W^{m,\frac{n}{m}}_0(\Omega), \|\nabla^m u\|_{\frac{n}{m}} \leq 1} \int_{\Omega} exp(\beta_{n,m} |u|^{\frac{n}{n-m}})dx \\
 \leq\\ \sup_{u \in W^{m,\frac{n}{m}}_{\mathcal{N}}(\Omega), \|\nabla^m u\|_{\frac{n}{m}} \leq 1} \int_{\Omega} exp(\beta_{n,m} |u|^{\frac{n}{n-m}})dx. \]

 Tarsi \cite{Tarsi} later extended Adams' result for the larger space $W^{m,\frac{n}{m}}_{\mathcal{N}}(\Om)$. The key step in her work is to embed $W^{m,\frac{n}{m}}_{\mathcal{N}}(\Om)$ into a Zygmund space. We state her embedding theorem below
 \begin{thm}\label{Tarsi_Thm}
 Let $n >2$ and $\Omega \subset \mathbb{R}^n$ be a bounded domain. Then there is a constant $C_n >0$ such that for all $u \in W^{m,\frac{n}{m}}_{\mathcal{N}}(\Omega)$ with $\|\nabla^m u \|_{\frac{n}{m}} \leq 1$, we have
 \begin{equation}
  \int_{\Om} e^{\beta |u|^{\frac{n}{n-m}}}dx <C_n |\Omega|  \hspace{1cm} \forall \beta \leq \beta_{n,m}               \label{Tarsi}
 \end{equation}
 and the constant $\beta_{n,m}$ appearing in (\ref{Tarsi}) is sharp and $\beta_{n,m}$ is same as in Theorem \ref{Adams_Thm}
 \end{thm}
  \begin{Rem}
  Here we remark that the bilinear form
 \begin{align}
  (u,v) \mapsto \int_{\Omega}\nabla^m u \cdot \nabla^m v= \left\{\begin{array}{lllll}
           \displaystyle\int_{\Omega} \Delta^k u \Delta^k v     \hspace{19mm}   \text{if } m =2k,\\
              \displaystyle \int_{\Omega} \nabla(\Delta^k u) \cdot \nabla (\Delta^k v)  \hspace{3mm} \text{ if } m=2k+1,        \label{poly-norm}
                       \end{array}\right.
 \end{align}
 defines a scalar product on $W^{m,2}_0(\Omega)$ and $W^{m,2}_{\mathcal{N}}(\Omega)$. Furthermore if $\Omega$ is bounded this scalar product induces a norm equivalent to (\ref{full-norm}).
 \end{Rem}
Therefore the above results imply that the problem $(P)$ nonlinearity of critical growth.
\begin{theorem}  \label{Multiplicity}
There exist positive real numbers  $\lambda_* \leq \lambda^*$ with $\lambda_*$ independent of $h$ such that the problem $(P)$ has at least two positive solutions for all $ \lambda \in (0, \lambda_*)$ and no solution for all $\lambda > \lambda^*$.
\end{theorem}
In spite of possible failure of Palais-Smale condition due to the presence of critical exponent we adapt the method of \cite{Tarantello} to prove the existence of the first solution by a decomposition of Nehari manifold into three parts. However for the existence of second solution we rely on the refined version of the Mountain-Pass Lemma introduced by Ghoussoub-Preiss \cite{Preiss}.
\section{Decomposition of Nehari Manifold}
Let $f(u) = \mu |u|^pue^{u^2}.$ The corresponding energy functional to the problem $(P)$ is given by
\begin{equation}
J(u)= \frac{1}{2} \int_{\Omega}|\nabla^m u|^2- \int_{\Omega} F(u) - \lambda \int_{\Om} hu         \label{energy}
\end{equation}
where $F(u) = \int_0^u f(s)ds.$ As the energy functional is not bounded below on $W^{m,2}_{\mathcal{N}}(\Om)$, we need to study $J(u)$ on the Nehari manifold
   \begin{equation}
   \mathcal{M}= \{u \in W^{m,2}_{\mathcal{N}}(\Om) \setminus \{0\}: \langle J'(u),u \rangle =0\},
   \end{equation}
where $J'(u)$ denotes the Frechet derivative of $J$ at $u$, and $\langle .,. \rangle$ is the inner product. Here we note that $\mathcal{M}$ contains every nonzero solution of the problem $(P)$.
We note that for any $u \in W^{m,2}_{\mathcal{N}}(\Om),$
\begin{align*}
\langle J'(u),u \rangle &=  \int_{\Omega} |\nabla^m u|^2 - \int_{\Om} f(u)u - \lambda \int_{\Om} hu, \\
\langle J''(u)u,u \rangle &= \int_{\Omega}|\nabla^m u|^2 - \int_{\Om} f'(u)u^2.
\end{align*}
Similarly to the method used in \cite{Tarantello}, We split $\mathcal{M}$ into three parts:
\begin{align*}
\mathcal{M}^0 = \{u \in \mathcal{M}: \langle J''(u)u,u \rangle = 0\},\\
\mathcal{M}^+ = \{u \in \mathcal{M}: \langle J''(u)u,u \rangle > 0\},\\
\mathcal{M}^- = \{u \in \mathcal{M}: \langle J''(u)u,u \rangle < 0\}.
\end{align*}
\section{Topological Properties of \texorpdfstring{$\mathcal{M}^0,\mathcal{M}^+, \mathcal{M}^-$}{Lg}}
Our first aim is to show, for some small $\lambda, \mathcal{M}^0 = \{0\}$. For this let $\zeta >0$, if $p>0$ and $\zeta < \frac{\lambda_1-\mu}{\mu}$ if $p=0$. Let $ \Lambda = \{u \in W^{m,2}_{\mathcal{N}}(\Om): \int_{\Om} |\nabla^m u|^2 \leq (1+\zeta)\int_{\Omega} f'(u)u^2\}.$
Lemma \ref{nonempty} implies that $\Lambda \neq \{0\}$. We now assume the following important hypothesis:
\begin{align}
 \lambda >0, \|h\|_{L^2(\Om)} = 1, \text{ and}     \notag\\
\inf_{u \in \Lambda \setminus \{0\}}\left(\mu \int_{\Om}(p+2u^2)|u|^{p+2}e^{u^2} - \lambda \int_{\Om}hu\right)>0.   \label{Nehari1}
\end{align}
The condition (\ref{Nehari1}) forces $\lambda$ to be suitably small. Indeed we can prove the following.
\begin{proposition}      \label{lambdastar}
Let
\begin{equation}
\lambda < \mu C_0^{\frac{p+3}{p+4}}|\Om|^{-\left(\frac{p+2}{2p+8}\right)} \label{lambda}
\end{equation}
 where $C_0 = \inf_{u \in \Lambda \setminus \{0\}} \int_{\Om}(p+2u^2)|u|^{p+2}e^{u^2} >0.$ Then (\ref{Nehari1}) holds.
\end{proposition}
\begin{proof}
\textbf{Step 1:} $\inf_{ u \in \Lambda \setminus \{0\}} \|u\|_{W^{m,2}_{\mathcal{N}}(\Om)} >0.$ \\
Assume the contradiction, then there exists a sequence $\{u_n\} \subset \Lambda \setminus \{0\}$ such that $\|u_n\|_{W^{m,2}_{\mathcal{N}}(\Om)} \to 0$ as $n \to \infty$. Let $v_n = \frac{u_n}{\|u_n\|_{W^{m,2}_{\mathcal{N}}(\Om)}}.$ Then $\|v_n\|_{W^{m,2}_{\mathcal{N}}(\Om)} =1$ and $v_n$ satisfies
\begin{equation}
1 \leq (1+\zeta) \int_{\Om} f'(u_n) v_n^2, \quad \forall n.      \label{lambda1}
\end{equation}
Since $u_n \to 0$ in $W^{m,2}_{\mathcal{N}}(\Om)$, by Adams' embedding for the higher order derivative in Theorem \ref{Tarsi_Thm} we get $f'(u_n) \to f'(0)$ in $L^r(\Om)$ for all $r \geq 1$. Since $v_n$ is bounded in $W^{m,2}_{\mathcal{N}}(\Om)$, $v_n$ has a weak limit say $v$ in $W^{m,2}_{\mathcal{N}}(\Om)$. Certainly $\|v\|_{W^{m,2}_{\mathcal{N}}(\Om)} \leq 1$ and up to a subsequence denote it same as $v_n$ which converges strongly to $v$ in $L^r(\Omega)$ for all $r \geq 1$. Hence from (\ref{lambda1}) we get
\begin{equation}
\int_{\Omega} |\nabla^m u|^2 \leq 1 \leq (1+\zeta) f'(0) \int_{\Om} v^2.        \label{lambda2}
\end{equation}
This gives a contradiction if $p>0$ in which case $f'(0)= 0$. If $p=0$, by assumption
$$\int_{\Omega} |\nabla^m u|^2 \geq \lambda_1 \int_{\Om} v^2 > (1+\zeta) \mu \int_{\Om}v^2$$ which gives a contradiction to (\ref{lambda1}) since $f'(0)= \mu$. This proves Step 1.\\
It is easy to check that using Step 1 and the definition of $\Lambda$:
\begin{equation}
 \inf_{u \in \Lambda \setminus \{0\}}\int_{\Omega}(p+2u^2)|u|^{p+2}e^{u^2} =C_0 >0.                   \label{lambda-inf}
\end{equation}

\noindent\textbf{Step 2:} Finally we have,
\begin{align*}
 \lambda \left|\int_{\Omega} hu\right| &\leq \lambda \|u\|_{L^2(\Omega)} \leq \lambda |\Omega|^{\frac{p+2}{2p+8}}\left(\int_{\Omega} |u|^{p+4} \right)^{\frac{1}{p+4}}\\
 &\leq \frac{\lambda |\Omega|^{\frac{p+2}{2p+8}}}{\mu (p+2u^2)|u|^{p+2}e^{u^2})^{\frac{p+3}{p+4}}} (\mu \int_{\Omega} (p+2u^2)|u|^{p+2}e^{u^2})\\
 &\leq \left(\frac{\lambda |\Omega|^{\frac{p+2}{2p+8}}}{\mu C_0^{\frac{p+3}{p+4}}}\right) (\mu \int_{\Omega}(p+2u^2)|u|^{p+2}e^{u^2}).
\end{align*}
Hence from the above inequality together with (\ref{lambda}) and (\ref{lambda-inf}) the proof is complete.
\end{proof}
\begin{lem}   \label{Null Set}
Suppose $\lambda >0$ be such that (\ref{Nehari1}) holds. Then $\mathcal{M}^0 = \{0\}$.
\end{lem}
\begin{proof}
Let $u \in \mathcal{M}^0, u \neq 0$. Then we have
\begin{align}
\int_{\Omega}|\nabla^m u|^2 &= \int_{\Om} f(u) u + \lambda \int_{\Om}hu,           \label{Nehari2}\\
\int_{\Omega}|\nabla^m u|^2 &= \int_{\Om} f'(u)u^2.                          \label{Nehari3}
\end{align}
We note that from (\ref{Nehari3}) $$\int_{\Omega}|\nabla^m u|^2 = \int_{\Om} f'(u)u^2 < (1+\zeta) \int_{\Om} f'(u)u^2$$ it implies that $ u \in \Lambda \setminus \{0\}$. From these two expressions we get
$$\lambda \int_{\Om} hu = \int_{\Om}(f'(u)u-f(u))u = \mu \int_{\Om}(p+2u^2)|u|^{p+2}e^{u^2}$$
which violates the condition (\ref{Nehari1}). Therefore $\mathcal{M}^0 = \{0\}$.
\end{proof}

Next we are going to discuss the topological properties of $\mathcal{M}^+$ and $\mathcal{M}^-$. Given $u \in {W^{m,2}_{\mathcal{N}}(\Om)} \setminus \{0\}$, we define $\xi_u: \mathbb{R}^+ \to \mathbb{R}$ by
\begin{equation}
\xi_u(s) = s \int_{\Om} |\nabla^m u|^2 - \int_{\Om} f(su)u.      \label{xi}
\end{equation}
The choice of the above function is consequence of the following expression,
\begin{equation*}
\langle J'(su), su \rangle = s \left(s \int_{\Om} |\nabla^m u|^2 - \int_{\Om} f(su)u - \lambda \int_{\Om}hu \right).
\end{equation*}
So, $\xi_u(s)=\lambda \int_{\Om}hu$ if and only if $su \in \mathcal{M}$, for $s >0$.

Now we are ready to state the following lemma.
\begin{lem}      \label{nonempty}
For every $u \in W^{m,2}_{\mathcal{N}}(\Om) \setminus \{0\}$ there exists a unique $s_* =s_*(u) >0$ such that $\xi_u(.)$ has its maximum at $s_*$ with $\xi_u(s_*)>0$. Also there holds $s_*u \in \Lambda \setminus \{0\}$.
\end{lem}
\begin{proof}
Differentiating (\ref{xi}) we have,
\begin{equation}
\xi'_u(s) = \int_{\Om} |\nabla^m u|^2- \int_{\Om} f'(su)u^2.     \label{xi1}
\end{equation}
Observe that,
\begin{align}
s^2 \xi'_u(s)&=  \int_{\Om}|\nabla^m (su)|^2 - \int_{\Om}f'(su)(su)^2     \notag \\
             &= \langle J''(su)su,su \rangle.                  \label{xi2}
\end{align}
Now we note that, $\xi_u(.)$ is strictly concave function on $\mathbb{R}^+$, since
\begin{equation}
\xi''_u(s)=-\int_{\Om}f''(su)u^3 <0.               \label{xi3}
\end{equation}
Also from the range of $\mu$ we get
\begin{align*}
\lim_{s \to 0+} \xi'_u(s) >0 \text{ and}\\
\lim_{s \to \infty} \xi_u(s) = -\infty.
\end{align*}
Hence there exists a unique maximum point of $\xi_u(.),$ say $s_* = s_*(u)>0$. Now using (\ref{xi1}) at $s =s_*$ in the definition of $\xi_u$, we deduce,
 \begin{align}
 \xi_u(s_*)&= s_* \int_{\Om} f'(s_*u) u^2 - \int_{\Om} f(s_*u)u,\quad \text{since } \xi'_u(s_*)=0  \notag\\
           &= \frac{1}{s_*} \int_{\Om} (f'(s_*u)s_*u - f(s_*u))s_*u  \notag\\
           &=\frac{\mu}{s_*} \int_{\Om} (p+2(s_*u)^2)|s_*u|^{p+2}e^{(s_*u)^2} > 0.     \label{xi4}
 \end{align}
 here we note that $f'(s)s-f(s) =\mu(p+2s^2)|s|^p s e^{s^2}$. Finally $$s_* \xi'_u(s_*) = \int_{\Om} |\nabla^m(s_*u)|^2 - \int_{\Om}f'(s_*u)(s_*u)^2 =0$$ which implies that $s_*u \in \Lambda \setminus \{0\}$.
\end{proof}

\begin{lem}
Let $\lambda$ be such that (\ref{Nehari1}) holds. Then, for every $u \in W^{m,2}_{\mathcal{N}}(\Om) \setminus \{0\}$, there exists a unique $s_ = s_-(u)>0$ such that $s_-u \in \mathcal{M}^-, s_- > s_*$ and $J(s_-u) = \max_{s \geq s_*} J(su) \hspace{2mm}\forall s \in [s_*, \infty), s \neq s_-.$ Also if we assume $\int_{\Om} hu >0$, then there exists a unique $s_+=s_+(u) >0$ such that $s_+u \in \mathcal{M}^+$. In particular $ s_+ <s_*$ and $J(s_+u) \leq J(su)$ for all $s \in [0, s_-].$      \label{maxima-minima}
\end{lem}
\begin{proof}
Define the functional $\rho_u :[0, \infty) \to \mathbb{R}$ by $\rho_u(s) = J(su)$. Then it is easy to verify that $\rho_u \in C^2([0,\infty],\mathbb{R}) \cap C((0,\infty),\mathbb{R})$. Then we have $$\rho'_u(s) = \xi_u(s)-\lambda \int_{\Om} hu, \hspace{2mm}\rho''_u(s)= \xi'_u(s), \hspace{1mm} \forall t >0.$$
Now from (\ref{Nehari1}) and (\ref{xi4}) we have,
\begin{equation*}
\xi_u(s_*) - \lambda \int_{\Om} hu = \frac{1}{s_*}\left\{ \mu \int_{\Om} (p+2(s_*u)^2)|s_*u|^{p+2}e^{(s_*u)^2} - \lambda \int _{\Om} h(s_*u)\right\} >0
\end{equation*}
Since $\xi_u(.)$ is strictly decreasing in $(s_*,\infty)$ and $\lim_{t \to \infty} \xi_u(s) = -\infty$, there exists a unique $s_- = s_-(u) > s_*$ such that $\xi_u(s_-) = \lambda \int_{\Om} hu.$ That is $s_-u \in \mathcal{M}$. One has $s_+ > s_*$ and $\rho_u'(s_) <0$, we get $s_-u \in \mathcal{M^-}$.

On the other hand when $\int_{\Om} hu >0$ we have $\lim_{s \to 0+} \xi_u(s) <0$ and which gives for $s$ close to $0$, $\xi_u(s)- \lambda \int_{\Om} hu <0$. Hence there exists a unique $s_+$ such that $\xi_u(s_+) =  \lambda \int_{\Om} hu$ which implies $s_+u \in \mathcal{M}.$ From the graph we see that $\xi_u(.)$ is strictly decreasing in $(0,s_*)$. Hence we have $s_+u \in \mathcal{M}^+$.

And the remaining properties of $s_-, s_+$ can be proved by analyzing the identity $\rho_u(s) = \xi_u(s) - \lambda \int_{\Om}hu.$
\end{proof}
\begin{Rem}
If we define the positive cone $\mathcal{P}=\{ u \in W^{m,2}_{\mathcal{N}}(\Om): \int_{\Om}hu >0\}$ in $W^{m,2}_{\mathcal{N}}(\Om)$. Then we note that $\mathcal{M^+} \subset \mathcal{P}$.
\end{Rem}
The next corollary shows some topological properties of $\mathcal{M}^+, \mathcal{M}^-$.
\begin{corl}
Let $S_{W^{m,2}_{\mathcal{N}}(\Om)} = \{u \in W^{m,2}_{\mathcal{N}}(\Om): \|u\|_{W^{m,2}_{\mathcal{N}}(\Om)} = 1\}$. Then there exists a diffeomorphism  $S^+: S_{W^{m,2}_{\mathcal{N}}(\Om)} \to \mathcal{M}^-$ defined by $S^+(u) = s_+(u)u$. Also $\mathcal{M}^+$ is homeomorphic to $S_{W^{m,2}_{\mathcal{N}}(\Om)} \cap \mathcal{P}$.
\end{corl}
\begin{proof}
The function $S^+$ is continuous because $s_+$ is continuous as an application of implicit function theorem applied to $(s,u) \to \xi_u(s) - \lambda \int_{\Om}hu$. And we deduce the continuity of $(S^+)^{-1}$ by the fact that $(S^+)^{-1}(w) = \frac{w}{\|w\|}.$ In a similar manner we can prove that $ \mathcal{M}^+$ is homeomorphic to $S_{W^{m,2}_{\mathcal{N}}(\Om)} \cap \mathcal{P}$.
\end{proof}

Relying on the embedding of $W^{m,2}_{\mathcal{N}}(\Om) \hookrightarrow L^q(\Om)$ for all $1 \leq q <\infty$ and using the estimate $F(s) \leq \frac{\mu|s|^p}{2}(e^{s^2}-1),$ for all $s \in \mathbb{R}$ we have the following lemma on the lower bound and upper bound.
\begin{lem}   \label{upper-lower}
There exists $C_1, C_2>0$ such that $$ -C_2 \lambda^{2p+8} \geq \theta_0 \geq -C_1 \lambda^{\frac{p+3}{p+4}}.$$ Where, $\theta_0 = \inf\{J(u) : u \in \mathcal{M}\}.$
\end{lem}
\begin{proof}
We prove the case of the lower bound.\\
Let $u \in \mathcal{M}$ then from the definition,
\begin{align*}
 J(u)=\frac{1}{2} \int_{\Omega} |\nabla^m u|^2 -\int_{\Omega}F(u) -\lambda \int_{\Omega}hu\\
    =\int_{\Omega} \left[ \frac{1}{2}f(u)u-F(u)\right] -\frac{\lambda}{2} \int_{\Omega}hu.
\end{align*}
We note that a simple calculation gives
\begin{equation}
 F(t) \leq \frac{\mu |t|^p}{2}(e^{t^2}-1), \quad \mbox{ for all } s \in \mathbb{R}.                  \label{F-estimate}
\end{equation}
Using (\ref{F-estimate}) we get
\begin{align}
 J(u) &\geq \frac{\mu}{2} \int_{\Omega} \left((u^2-1)e^{u^2}+1\right) - \frac{\lambda}{2} \int_{\Omega}hu    \notag \\
      &\geq \frac{c\mu}{2} \int_{\Omega} |u|^{p+4} -\frac{\lambda}{2}\int_{\Omega} hu,        \label{lower-bound-lemma}
\end{align}
since $(s^2-1)e^{s^2}+1 \geq cs^4$ for some $c>0$, for all $s \in \mathbb{R}$. By an application of Holder inequality we get
\begin{equation}
 \int_{\Omega} hu \leq |\Omega|^{\frac{p+2}{2(p+4)}} \|u\|_{L^{p+4}(\Omega)}.          \label{lower-bound-holder}
\end{equation}
From (\ref{lower-bound-lemma}) and (\ref{lower-bound-holder}) we get,
\begin{equation}
 J(u) \geq \frac{c \mu}{2} \|u\|^{p+4}_{L^{p+4}(\Omega)} - \left( \frac{\lambda |\Omega|^{\frac{p+2}{2(p+4)}}}{2}\right)\|u\|_{L^{p+4}(\Omega)}.   \label{lower-bound}
\end{equation}
By considering the global minimum of the function
\[ \omega(x)= \left(\frac{c \mu}{2}\right) x^{p+4} - \left( \frac{\lambda |\Omega|^{\frac{p+2}{2(p+4)}}}{2}\right) x,\]
It can be shown that \[J(u) \geq -C_1 \lambda^{\frac{p+4}{p+3}}.\]

      In a similar fashion we can prove the upper bound for $J$.
\end{proof}

As a consequence of Lemma \ref{Null Set} we have:
\begin{lem} \label{Implicit-Inner}
Let $\lambda$ and $h$ satisfy (\ref{Nehari1}). Given $u \in \mathcal{M} \setminus \{0\}$ there exists $\delta >0$ and a differentiable function $s:\{w \in W^{m,2}_{\mathcal{N}}(\Om): \|w\|_{W^{m,2}_{\mathcal{N}}(\Om)} < \delta \}\to \mathbb{R}$, with $$s(0)=1, s(w)(u-w) \in \mathcal{M}, \hspace{1.5mm}\forall \quad \|w\|_{W^{m,2}_{\mathcal{N}}(\Om)} < \delta$$ and
\begin{equation}
\langle s'(0), v \rangle= \frac{2 \int_{\Om} \nabla^m u \cdot \nabla^m v - \int_{\Om} (f'(u)u+f(u))v - \lambda \int_{\Om}hv}{\int_{\Om} |\nabla^m u |^2 - \int_{\Om}f'(u)u^2}                        \label{Nehari4}
\end{equation}
\end{lem}

\begin{proof}
We define the function  $G: \mathbb{R} \times W^{m,2}_{\mathcal{N}}(\Om) \to \mathbb{R}$ by,
$$G(s,w) = s \int_{\Om} |\nabla^m (u-w)|^2 - \int_{\Om} f(s(u-w))(u-w) - \lambda \int_{\Om} h(u-w).$$
Then $ G \in C^1( \mathbb{R} \times W^{m,2}_{\mathcal{N}}(\Om); \mathbb{R})$ and since $u \in \mathcal{M}$ it implies $G(1,0)= \int_{\Om} |\nabla^m u|^2 - \int_{\Om}f(u)u - \lambda \int_{\Om} hu = 0.$ Also $G_s(1,0) \neq 0$, indeed $G_s(1,0)= \int_{\Om} |\nabla^m u|^2 - \int_{\Om} f'(u)u^2 \neq 0$ thanks to Lemma \ref{Null Set}. Then by the Implicit Function Theorem, there exists $\delta >0, s:\{w \in W^{m,2}_{\mathcal{N}}(\Om): \|w\| < \delta\} \to \mathbb{R}$ of class $C^1$ that satisfies:
\begin{align*}
G(s(w),w) &= 0 \text{ for all } w \in W^{m,2}_{\mathcal{N}}(\Om), \|w\|_{W^{m,2}_{\mathcal{N}}(\Om)} < \delta,\\
s(0) &= 1.
\end{align*}
Also
\begin{align*}
0 &= s(w)G(s(w),w)\\
  &= \int_{\Om} (s(w) |\nabla^m (u-w)|)^2 - \int_{\Om}f(s(w)(u-w)(s(w)(u-w)) -\lambda \int_{\Om}h(s(w)(u-w)),
\end{align*}
that is $s(w)(u-w) \in \mathcal{M}$ for all $w \in W^{m,2}_{\mathcal{N}}(\Om)$ with $\|w\| < \delta$.
Now if we differentiate the identity $G(s(w),w)=0$ with respect to $w$, we get
\[ 0 = \langle G_s(s(w),w) + G_w(s(w),w), v \rangle \text{ for all } v \in W^{m,2}_{\mathcal{N}}(\Om). \]
Putting $w=0$ in the above identity
\[ 0 = \langle G_s(1,0) s'(0) + G_w(1,0), v \rangle = G_s(1,0) \langle s'(0), v \rangle + \langle G_w(1,0), v \rangle \] and we deduce from above
\begin{align*}
\langle s'(0), v \rangle &= -\frac{\langle G_w(1,0),v \rangle}{G_s(1,0)}\\
                         &=  \frac{2 \int_{\Om} \nabla^m u \cdot \nabla^m v - \int_{\Om} (f'(u)u+f(u))v - \lambda \int_{\Om}hv}{\int_{\Om} |\nabla^m u|^2 - \int_{\Om}f'(u)u^2}.
  \end{align*}
\end{proof}

\section{Local Minimum of \texorpdfstring{$J$}  \text{ in} \texorpdfstring{$W^{m,2}_{\mathcal{N}}(\Om)$}.}
We are now in a situation to prove the existence of a minimizer for $J$ and hence we guarantee the existence of first solution.
\par Since $\mathcal{M}$ is a closed set of $W^{m,2}_{\mathcal{N}}(\Om)$, hence a complete metric space. Now $J$ is bounded below on $\mathcal{M}$. By the Ekeland's Variational Principle there exists a sequence $\{u_n\} \subset \mathcal{M} \setminus \{0\}$ satisfying:
\begin{equation}
J(u_n) < \theta_0 + \frac{1}{n}, \hspace{2mm} J(v) \geq J(u_n)- \frac{1}{n}\|v-u_n\|_{W^{m,2}_{\mathcal{N}}(\Om)} \hspace{1mm} \forall v \in \mathcal{M}                  \label{Ekeland}
\end{equation}
\begin{proposition}
Let $\lambda$ and $h$ satisfy (\ref{Nehari1}). Then $$\lim _{n \to \infty}\|J'(u_n)\|_{(W^{m,2}_{\mathcal{N}}(\Om))^{-1}} = 0.$$
\end{proposition}
\begin{proof}
We proceed in a few steps. With the help of Lemma \ref{upper-lower} we've $\varliminf_{n \to \infty}\|u_n\|_{W^{m,2}_{\mathcal{N}}}>0$.
\textbf{Claim 1:} $\varliminf_{n \to \infty} \int_{\Omega} (p+2u_n^2)|u_n|^{p+2}e{u_n^2} >0$.\\
If possible let's assume that for a subsequence of $\{u_n\}$, which is still denoted by $\{u_n\}$, we have
\begin{equation}
 \lim_{n \to \infty} (p+2u_n^2)|u_n|^{p+2}e^{u_n^2} \to 0 \text{    as   } n \to \infty      \label{lim-dual}
\end{equation}
Here we note that $u_n \to 0$ in $L^q(\Omega)$ for all $q \in [1,\infty)$ using (\ref{lim-dual}), and if $p>0$,
\[\int_{\Omega} f(u_n)u_n = \mu \int_{\Omega} |u_n|^{p+2} e^{u_n^2} \to 0 \text{  as } n \to \infty.\]
Therefore we have $\int_{\Omega}f(u_n)u_n \to 0,$  $\int_{\Omega} hu_n  \to 0$ as $n \to \infty$. Which imply that $\|u_n\|_{W^{m,2}_{\mathcal{N}}} \to 0$ as $n \to \infty$ because $\{u_n\} \subset \mathcal{M}$ hence a contradiction to the fact that $\varliminf_{n \to \infty}\|u_n\|_{W^{m,2}_{\mathcal{N}}}>0$. Similar argument for $p=0$.\\
\textbf{Claim 2:} $\varliminf_{n \to \infty} \{ |\int_{\Omega}|\nabla^m u|^2 - \int_{\Omega} f'(u_n)u_n^2| >0\}.$\\
Let the claim doesn't hold. Then for a subsequence $\{u_n\}$ we have
\[ \int_{\Omega} |\nabla^m u|^2 - \int_{\Omega} f'(u_n)u_n^2 = o_n(1).\]
From this and the fact $\varliminf_{n \to \infty}\|u_n\|_{W^{m,2}_{\mathcal{N}}}>0$ we deduce that,
\[\varliminf_{n \to \infty}f'(u_n)u_n^2 >0.\]
Therefore we have $u_n \in \Lambda \setminus \{0\}$ for large $n$. Since $\{u_n\} \subset \mathcal{M}$ we get
\begin{align*}
 o_n(1) &= \lambda \int_{\Omega} hu_n + \int_{\Omega} [f(u_n)-f'(u_n)u_n]u_n\\
          &= -\mu \int_{\Omega} (p+2u_n^2)|u_n|^{p+2}e^{u_n^2} + \lambda \int_{\Omega} hu_n,
\end{align*}
which contradicts (\ref{Nehari1}). This completes the proof of the claim.\\
Now we proof the theorem. Let's assume $\|J'(u_n)\|_{(W^{m,2}_{\mathcal{N}})^{-1}}  >0$ for all large $n$ (otherwise obvious). Now we define $u=u_n \in \mathcal{M}$ and $w=\delta \frac{J'(u_n)}{\|J'(u_n)\|}$ (by Riesz representation theorem, we identify $J'(u_n)$ as an element in $W^{m,2}_{\mathcal{N}}(\Omega)$ still denote $J'(u_n)$) for $\delta >0$ small. Therefore we can apply Lemma \ref{Implicit-Inner} for $w$ small we get $s_n(\delta):= s\left[ \delta \frac{J'(u_n)}{\|J'(u_n)\|}\right] >0$ such that \[ w_{\delta}= s_n(\delta) \left[ u_n - \delta \frac{J'(u_n)}{\|J'(u_n)\|}\right] \in \mathcal{M}.\]
Now from (\ref{Ekeland}) and a Taylor expansion we have:
\begin{align*}
 \frac{1}{n} \|w_{\delta} - u_n\| &\geq J(u_n) - J(w_{\delta})\\
                       &=(1-s_n(\delta)) \langle J'(w_{\delta}),u_n \rangle + \delta s_n(\delta) \left\langle J'(w_{\delta}), \frac{J'(u_n)}{\|J'(u_n)\|} \right\rangle +o(\delta)
                       \end{align*}
Dividing by $\delta>0$ and taking limit as $\delta \to 0$ we get:
\[\frac{1}{n}(1+|s_n'(0)| \|u_n\|) \leq -s_n(0) \langle J'(u_n),u_n\rangle +\|J'(u_n)\|=\|J'(u_n)\|.\]
Hence \[ \|J'(u_n)\| \leq \frac{1}{n}(1+s_n'(0)|\|u_n\|).\]
We complete the proof by using, $|s_n'(0)|$ is uniformly bounded on $n$ by (\ref{Nehari4}) and using the Claim 2.
\end{proof}
\begin{theorem}
Let $\lambda,h$ satisfy (\ref{Nehari1}). Then there exists a nonnegative function $ u_0 \in \mathcal{M}^+$ such that $J(u_0)= \inf_{u \in \mathcal{M}\setminus \{0\}} J(u)$. Moreover, $u_0$ is a local minimum for $J$ in $W^{m,2}_{\mathcal{N}}(\Om)$.      \label{1st solution}
\end{theorem}
\begin{proof}
Let $\{u_n\}$ be a sequence which minimizes $J$ on $\mathcal{M} \setminus \{0\}$ as in (\ref{Ekeland}). \\
\noindent\textbf{Step 1:} $\liminf_{n \to \infty} \int_{\Om} hu_n >0$ and hence $u_n \in \mathcal{M}^+$. Indeed $u_n \in \mathcal{M}$ and making some suitable adjustments
 \begin{align}
J(u_n) &= \frac{p}{2(p+2)} \int_{\Om} |\nabla^m u|^2 + \int_{\Om} \left( \frac{1}{p+2} f(u_n)u_n - F(u_n)\right)   \notag\\
        & \quad  -\lambda \frac{p+1}{p+2} \int_{\Om} hu_n < -C \lambda^{2p+8}.    \label{liminf}
\end{align}
Thanks to Lemma \ref{upper-lower} there exists $C>0$.  Now we note that $F(t) < \frac{1}{p+2}f(t)t$ for all $ t \in \mathbb{R}$. Therefore we've from (\ref{liminf}), to make the inequality consistent with sign that \[ \liminf _{n \to \infty} \int_{\Om}hu_n >0.\]

\noindent\textbf{Step 2:} $\limsup _{n \to \infty} \|u\|_{W^{m,2}_{\mathcal{N}}(\Om)} < \infty$. \\
\textit{Case 1.} If $p>0$ then by the means of Sobolev embedding we derive boundedness of  $\{u_n\}$ in $W^{m,2}_{\mathcal{N}}(\Om)$. Using the fact from (\ref{liminf})
$$ \int_{\Om}|\nabla^m u|^2 \leq \lambda \int_{\Om} hu_n.$$
\textit{Case 2.} If $p=0$ by using the fact that $\frac{1}{2}f(t)t - F(t) \geq Ct^4$ for all $t \in \mathbb{R}$ and for some $C>0$ we deduce that $\{u_n\}$ is a bounded sequence in $L^2(\Om)$. And this gives that $\{F(u_n)\}$ is a bounded sequence in $L^1(\Om)$ using (\ref{liminf}) and hence $\{u_n\}$ is a bounded sequence in $W^{m,2}_{\mathcal{N}}(\Om)$.

\noindent \textbf{Step 3:} Existence of $u_0 \in \mathcal{M}^+$.\\
From the previous step up to a subsequence, $u_n \rightharpoonup u_0$ in $W^{m,2}_{\mathcal{N}}(\Om)$. Now from the Proposition 2.2 we note that $\{f(u_n)u_n\}$ is a bounded sequence in $L^1(\Om)$. Therefore from Vitali's convergence theorem (for details see Lemma 8.3 in \cite{Prashanth}), we get that $$\int_{\Om}f(u_n) \phi \to \int_{\Om}f(u_0)\phi, \text{ for all } \phi \in W^{m,2}_{\mathcal{N}}(\Om).$$ Hence $u_0$ will solve $(P)$, in particular $u_0 \in \mathcal{M}$. Here we note that $u_0 \neq 0$ as $h \neq 0$ that is $u_0 \in \mathcal{M} \setminus \{0\}$. We see that $ \theta_0 \leq J(u_0)$. From (\ref{liminf}) we get by using Fatou's Lemma that $\theta_0 = \liminf_{n \to \infty} J(u_n) \geq J(u_0)$. Therefore $u_0$ minimizes $J$ on $\mathcal{M} \setminus \{0\}$. Now we have to show $u_0 \in \mathcal{M}^+$. From the existence of $s_-(u_0)$ and $s_+(u_0)$ in Lemma \ref{maxima-minima} and using the fact $J(s_+(u_0)u_0) < J(s_-(u_0)u_0)$ we get $u_0 \in \mathcal{M}^+$. \\
\noindent \textbf{Step 4:} $u_0$ is a local minimum for for $J$ in $W^{m,2}_{\mathcal{N}}(\Om)$.
\par We see that $s_+(u_0) =1$ because $u_0 \in \mathcal{M}^+$ from Step 3. Also we have from the (\ref{maxima-minima}) we have
\begin{equation*}
s_+(u_0)=1 <s_*(u_0)
\end{equation*}
Now by the continuity of $s_*(u_0)$, for sufficiently small $\delta>0$
\begin{equation}
1 < s_*(u_0-w), \hspace{1.5mm}\forall w \in W^{m,2}_{\mathcal{N}}(\Om), \|w\|_{W^{m,2}_{\mathcal{N}}(\Om)} < \delta.    \label{minimizer}
\end{equation}
Now by the Lemma \ref{Implicit-Inner} for $\delta>0$ small enough if necessary, let $s: \{w \in W^{m,2}_{\mathcal{N}}(\Om): \|w\| < \delta\} \to \mathbb{R}$ such that $s(w)(u_0-w) \in \mathcal{M}$ and $s(0)=1$.
\noindent Whenever $s(w) \to 1 $ when $\|w\| \to 0$, we can assume that
\[ s(w) < s_*(u_0-w), \hspace{1.5mm} \forall w \in W^{m,2}_{\mathcal{N}}(\Om), \|w\| < \delta.\] Hence we get $s(w)(u_0-w) \in \mathcal{M}^+$ using the above inequality and Lemma \ref{maxima-minima}. Again by using the Lemma \ref{maxima-minima} we see,
\[J(s(u_0-w) \geq J(s(w)(u_0-w)) \geq J(u_0), \hspace{1.5mm} \forall s \in [0,s_*(u_0-w)].\]Hence from (\ref{minimizer}) we observe that $J(u_0-w) \geq J(u_0)$ for every $\|w\|_{W^{m,2}_{\mathcal{N}}(\Om)} <\delta$. This shows that $u_0$ is a local minimizer. \\
\noindent{\textbf{Step 5:}} A positive local minimum for $J$. If $u_0 \geq 0$ then we get the positivity by using the strong maximum principle. In case if $u_0 \ngeq 0$ then we consider $\tilde{u}_0 = s_+(u_0)|u_0|>0 \in \mathcal{M}^+$ and also from the definition $\rho_{u_0}(s)=\rho_{|u_0|}(s)$ for all $s>0$. Therefore we get $s_*(|u_0|)=s_*(u_0)$ and from the definition of $s_+$ we deduce $s_+(u_0) \leq s_+(|u_0|)$. Hence from Step 4, $s_+(|u_0|) \geq 1$. Therefore by Lemma \ref{maxima-minima} we get $J(\tilde{u}_0) \leq J(|u_0|)$. Now using the assumption $h \geq 0$ in $\Omega$, we have $J(|u_0|) \leq J(u_0)$ and which implies that $\tilde{u}_0$ minimizes $J$ on $\mathcal{M}\setminus \{0\}$. Hence by repeating the same argument as in Step 4 we get the desired result.
\end{proof}

\section{Existence of The Second Solution}
The existence of the second solution for $(P)$ depends on whether we can apply some version of Mountain Pass Lemma. We wish to look for a solution of the form $u_1 = v + u_0$ where $u_0$ is the local minimum for the functional (\ref{energy}). Then we see that $u_1$ will solve $(P)$ whenever $v$ solves the following equation:
 $$(P_1)\hspace{1cm}\left\{
\begin{array}{lllll}\left.
\begin{array}{lllll}
(-\Delta)^m v &= f(v+u_0)-f(u_0)\\
\hspace{8mm} v &> 0 \end{array}\right\} \;\; \text{in} \;\Om, \\
v = \De v = 0=..=\De^{m-1}u  \hspace{15.5mm}\text{ on } \partial \Omega. \end{array}\right.
$$
We can write the above PDE as following
  $$(\tilde{P}) \hspace{1cm}\left\{
\begin{array}{lllll}\left.
\begin{array}{lllll}
\D2 v &= \tilde{f}(x,v)\\
\hspace{6mm} v &> 0 \end{array}\right\} \;\; \text{in} \;\Om, \\
v = \De v = 0=..=\De^{m-1}v \text{  on } \partial \Omega, \end{array}\right.
$$
by introducing the function $\tilde{f}: \Om \times \mathbb{R} \to \mathbb{R}$ and we define by
\begin{align*}
\tilde{f}(x,s) &= f(s+u_0(x))-f(u_0(x)) \text{ if } s \geq 0,\\
               &= 0 \hspace{2.5cm}\text{   otherwise.}
  \end{align*}
The energy functional corresponding to $(\tilde{P})$ is $J_{u_0}: {W^{m,2}_{\mathcal{N}}(\Om)} \to \mathbb{R}$ defined by
\[ J_ {u_0}(v) = \frac{1}{2}\int_{\Om}|\nabla^m v|^2 - \int_{\Om} \tilde{F}(x,v)dx,\] where $\tilde{F}(x,s) = \int_0^s \tilde{f}(x,t)dt$.
Now onwards, we denote $J_{u_0}$ by $J_0$. These type of functionals were studied by \cite{Zhao-Chang}, \cite{Peral}. We now state the Generalized Mountain Pass Lemma that was introduced by Ghoussoub-Preiss \cite{Preiss}.
\begin{Def}
Let $H$ be a closed subspace of the Banach Space $W^{m,2}_{\mathcal{N}}(\Om)$. We say that a sequence $\{v_n\} \subset W^{m,2}_{\mathcal{N}}(\Om)$ is a Palais-Smale sequence for $J_0$ at the level $c$ around $H$ if:
\begin{itemize}
\item[(i)] $\lim_{n \to \infty} dist(v_n,H) =0$
\item[(ii)] $\lim_{n \to \infty} J_0(v_n) = c$
\item[(iii)] $\lim_{n \to \infty} \|J_0'(v_n)\|_{(W^{m,2}_{\mathcal{N}}(\Om))^{-1}} = 0.$
\end{itemize}
And we say such a sequence a $(PS)_{H,c}$ sequence.
\end{Def}
\begin{Rem}
In case $H =W^{m,2}_{\mathcal{N}}(\Om)$, the above definition coincides with the usual Palais-Smale sequence at the level $c$.
\end{Rem}
\begin{lem} \label{mountain}
Let $H \subset W^{m,2}_{\mathcal{N}}(\Om)$ be a closed set,$c \in \mathbb{R}$. Assume $\{v_n\} \subset W^{m,2}_{\mathcal{N}}(\Om)$ be a $(PS)_{H,c}$ sequence. Then (upto a subsequence), $v_n \rightharpoonup v_0$ in $W^{m,2}_{\mathcal{N}}(\Om)$, and
\begin{equation}
\lim _{n \to \infty} \int_{\Om} \tilde{f}(x,v_n) = \int_{\Om} \tilde{f}(x,v_0), \hspace{1cm} \lim_{n \to \infty} \int_{\Om} \tilde{F}(x,v_n) = \int_{\Om} \tilde{F}(x,v_0).     \label{Vitali}
\end{equation}
\end{lem}
\begin{proof}
From the fact that $\{v_n\}$ is a $(PS)_{H,c}$ sequence we have:
\begin{align}
\frac{1}{2}\int_{\Om} |\nabla^m v_n|^2 - \int_{\Om} \tilde{F}(x,v_n) &= c_0 + o_n(1),             \label{GPS1}\\
\left|\int_{\Om} \nabla^m v_n \cdot \nabla^m \phi - \int_{\Om}\tilde{f}(x,v_n) \phi \right| &\leq o_n(1) \|\phi\|_{W^{m,2}_{\mathcal{N}}(\Om)}, \hspace{5mm} \forall \phi \in W^{m,2}_{\mathcal{N}}(\Om).    \label{GPS2}
\end{align}
Now we claim that,\\
\textbf{Claim:} $\sup_n \|v_n\|_{W^{m,2}_{\mathcal{N}}(\Om)}  < \infty, \hspace{3mm}\sup_n \int_{\Om} \tilde{f}(x,v_n) < \infty.$\\
Given any $\epsilon >0$ there exists $s_{\epsilon} >0$ such that
\begin{equation}
\int_{\Om} \tilde{F}(x,s) \leq \epsilon s \tilde{f}(x,s) \text{ for all } |s| \geq s_{\epsilon}.     \label{GPS3}
\end{equation}
Using (\ref{GPS1}) and (\ref{GPS3}), we see
\begin{align}
\frac{1}{2}\int_{\Om} |\nabla^m v_n|^2  &\leq \int_{\Om \cap \{|v_n| \leq s_{\epsilon}\}} \tilde{F}(x,v_n) + \int_{\Om \cap \{|v_n| \geq s_{\epsilon}\}} \tilde{F}(x,v_n) +c+o_n(1)            \notag \\
                 &\leq \int_{\Om \cap \{|v_n| \leq s_{\epsilon}\}}\tilde{F}(x,v_n) + \epsilon \int_{\Om} \tilde{f}(x,v_n)v_n +c+o_n(1)   \notag\\
                  &\leq C_{\epsilon} + \epsilon \int_{\Om} \tilde{f}(x,v_n)v_n.      \label{GPS4}
\end{align}
Now from (\ref{GPS4}) we obtain,
\begin{align*}
\int_{\Om}\tilde{f}(x,v_n)v_n &\leq \int_{\Om}|\nabla^m v_n|^2 + o_n(1) \|v_n\|_{W^{m,2}_{\mathcal{N}}(\Om)}\\
                             &\leq 2 C_{\epsilon} + 2\epsilon \int_{\Om} \tilde{f}(x,v_n)v_n +o_n(1)\|v_n\|_{W^{m,2}_{\mathcal{N}}(\Om)}
\end{align*}
by substituting $\phi = v_n$ in (\ref{GPS2}).\\
Hence by choosing $\epsilon$ small enough if needed we get
\begin{equation}
 \int_{\Om} \tilde{f}(x,v_n)v_n \leq \frac{2C_{\epsilon}}{1-2\epsilon} + o_n(1)\|v_n\|_{W^{m,2}_{\mathcal{N}}(\Om)}.   \label{GPS5}
\end{equation}
We conclude the claim using (\ref{GPS5}), (\ref{GPS2}) and also $\sup_n\int_{\Omega}\tilde{f}(x,v_n)v_n <\infty$. \\
Since $\{v_n\} \subset W^{m,2}_{\mathcal{N}}(\Om)$ is bounded, up to a subsequence, $v_n \rightharpoonup v_0$ in $W^{m,2}_{\mathcal{N}}(\Om)$, for some $v_0 \in W^{m,2}_{\mathcal{N}}(\Om)$.\\
To prove (\ref{Vitali}) we consider $A$ to be a $2m$ dimensional Lebesgue measure of a set $A \subset \mathbb{R}^{2m}$.
\par Let $C=\sup_n \int_{\Omega}|\tilde{f}(x,v_n)v_n|<\infty$ from the above claim. Given $\epsilon >0$, we define
\[\mu_{\epsilon}= \max_{x \in \bar{\Omega}, |s|\leq \frac{2C}{\epsilon}}|\tilde{f}(x,s)s|.\] Then, for any $A \subset \Omega$ with $|A| \leq \frac{\epsilon}{2C}$, we have
\begin{align*}
 \int_{A} |\tilde{f}(x,v_n)| &\leq \int_{A \cap \{|v_n| \geq \frac{2C}{\epsilon}\}} \frac{|\tilde{f}(x,v_n)v_n|}{|v_n|} + \int_{A \cap \{|v_n| \leq \frac{2C}{\epsilon}\}} |\tilde{f}(x,v_n)|\\
 &\leq \frac{\epsilon}{2} + \mu_{\epsilon}|A| \leq \epsilon.
\end{align*}
Hence $\{\tilde{f}(x,v_n)\}$ is an equi-integrable family in $L^1(\Omega)$ and so is $\{\tilde{F}(x,v_n)\}$(we note that $|\tilde{F}(x,t)| \leq C_1 |\tilde{f}(x,t)|$ for all $x \in \bar{\Omega}, t \in \mathbb{R}$, for some $C_1>0$). By applying the Vitali's convergence theorem we get conclude the lemma.
 \end{proof}

Certainly $J_0(0)=0$ and $v=0$ is a local minimum for $J_0$. Also we have $$ \lim_{s \to \infty}J_0(sv) = -\infty \text{ for any } v \in W^{m,2}_{\mathcal{N}}(\Om)\setminus \{0\}.$$  Hence we can fix $e \in W^{m,2}_{\mathcal{N}}(\Om)\setminus \{0\}$ such that $J_0(e) <0$. Now we define the mountain pass level $$ c_0 = \inf_{\gamma \in \Gamma} \sup_{s \in [0,1]} J_0(\gamma(s)).$$
Where $\Gamma = \{ \gamma \in C([0,1],W^{m,2}_{\mathcal{N}}(\Om)): \gamma(0)=0, \gamma(1) = e\}$. Then from the definition of $c_0$ it follows $c_0 \geq 0$. Define $R_0 =\|e\|_{W^{m,2}_{\mathcal{N}}(\Om)},$ we note that $\inf\{J_0(v):\|v\|_{W^{m,2}_{\mathcal{N}}(\Om)}=R\} = 0$ for all $R \in (0,R_0)$. And we now let $H = W^{m,2}_{\mathcal{N}}(\Om)$ if $c_0 >0$ and $H=\{\|v\|_{W^{m,2}_{\mathcal{N}}(\Om)} =\frac{R_0}{2}\}$ if $c_0=0$. We now state as now the lemma giving upper bound for $c_0$
\begin{lem}The upper bound of the Mountain Pass level is below
\begin{equation}
c_0 < \frac{(4 \pi)^m m!}{2}= \frac{\beta_{2m,m}}{2}. \label{MPL1}
\end{equation}
\end{lem}
\begin{proof}
Without loss of generality we can assume that the unit ball $B_0(1) \subset \Om$. For any $\epsilon >0$ we define
\begin{align}
{\tilde{\tau}}_n(x) := \left\{
\begin{array}{lllll}
\sqrt{\frac{1}{2M}\log n}+\frac{1}{\sqrt{2M \log n}} \sum_{\gamma =1}^{m-1} \frac{(1-k|x|^2)^{\gamma}}{\gamma} \hspace{6mm} |x| \in [0, \frac{1}{\sqrt{n}})\\%{\sqrt{16 \pi^2  \log n}} + \frac{1}{\sqrt{16 \pi^2 \log n }}, \hspace{6mm} |x| \in [0, \frac{1}{\sqrt{n}}),\\
-\sqrt{\frac{2}{M \log n}} \log |x|,                \hspace{3.75cm} |x| \in [\frac{1}{\sqrt{n}},1),\\
\chi_n(x),       \hspace{5.59cm} |x| \in [1,\infty)
\end{array}\right.
\end{align}
where \[ M = \frac{(4 \pi)^m (m-1)!}{2}, \hspace{1cm} \chi_n \in \mathcal{C}^{\infty}_0(\Om), \quad \chi_n|_{\partial B_1(0)}= \chi_n|_{\partial \Omega} =0.\]
Furthermore, for $\gamma=1,2,..,m-1, D^{\gamma} \chi_n|_{\partial B_1(0)}= (-1)^{\gamma}(\gamma-1)! \sqrt{\frac{2}{M \log n}}$,  $\Delta^j \chi_n|_{\partial \Omega} = 0$ for $j =0,1,2,...,[(m-1)/2]$
and $ \chi_n, |\nabla \chi_n|, \De \chi_n$ are all $O\left(\frac{1}{\sqrt{2 \log n}}\right)$. Then, $\tilde{\tau}_n \in W^{m,2}_{\mathcal{N}}(\Om).$
Now we normalize $\tilde{\tau}_n$, setting $$ {\tau}_n := \frac{\tilde{\tau}_n}{\|\tilde{\tau}_n\|_{W^{m,2}_{\mathcal{N}}(\Om)}} \in  W^{m,2}_{\mathcal{N}}(\Om).$$
Suppose (\ref{MPL1}) is not true. This means that, for some $s_n>0$ (see \cite{Lam-Lu}),
\[J_0(s_n \tau_n)= \sup_{s>0}J_0(s \tau_n) \geq \frac{(4 \pi)^mm!}{2}  \hspace{1cm} \forall n.\]
Hence
\begin{equation}
 \frac{s_n^2}{2} - \int_{\Om} \tilde{F}(x,s_n \tau_n) \frac{(4 \pi)^m m!}{2} \hspace{1cm} \forall n.     \label{MPL2}
\end{equation}
 It follows that $\frac{d}{ds}J_0(s \tau_n)=0$ at the point of maximum $s=s_n$ for $J_0$, we get
\begin{equation}
s_n^2 = \int_{\Om} \tilde{f}(x,s_n \tau_n)(s_n \tau_n).     \label{MPL3}
\end{equation}

 Now we note that from the definition of $\tilde{f}$ we see that $\inf_{x \in \bar{\Om}} \tilde{f}(x,s) \geq e^{s^2}$ for $|s|$ large. Then from (\ref{MPL2}) we get for sufficiently large $n$
 \begin{align}
 s_n^2 &\geq \int_{\{|x| \leq \frac{1}{\sqrt{n}}\}} \tilde{f}(x,s_n \tau_n)(s_n \tau_n) \geq \int_{\{|x| \leq \frac{1}{\sqrt{n}}\}} e^{s_n^2 \tau_n ^2}(s_n \tau_n)        \notag\\
        &\geq e^{{s_n^2} \frac{\log n}{2M}} \frac{s_n}{\sqrt{2M}} \sqrt{\log n} \frac{\alpha_{2m}}{n^m}   \notag\\
        &=\frac{\alpha_{2m}}{\sqrt{2M}} e^{\left(\frac{s_n^2}{2M}-m\right)\log n} s_n (\log n)^{\frac{1}{2}},    \label{MPL4}
 \end{align}
  where $\alpha_{2m}$ is the volume of the unit ball in $\mathbb{R}^{2m}$. Using the fact $s_n^2 \geq (4 \pi)^m m!$ from (\ref{MPL2}) and (\ref{MPL4}) it follows that $s_n$ is bounded and also $s_n^2 \to (4 \pi)^m m!$. Also from (\ref{MPL4}) we note \[ s_n \geq \frac{w_{2m}}{\sqrt{2M}} (\log n)^{\frac{1}{2}}, \text{ for all large } n\] which gives the contradiction.
\end{proof}
We now prove the theorem regarding the existence of second solution.
\begin{theorem}  \label{2nd solution}
Given a local minimum $u_0$ of $J$ in $W^{m,2}_{\mathcal{N}}(\Om)$, there exists a point $v_0 \in W^{m,2}_{\mathcal{N}}(\Om)$ with $v_0 >0$ in $\Om$, such that $J_0'(v_0)=0$.
\end{theorem}
\begin{proof}
From Lemma \ref{MPL1} we have $c_0 \in \left[0, \frac{(4 \pi)^ m m!}{2}\right).$ Consider $\{v_n\}$ be a Palais-Smale sequence for $J_0$ at the level $c_0$ around $H$ (such a $(PS)_{H,c_0}$ sequence exists \cite{Preiss}). Then up to a subsequence $v_n \rightharpoonup v_0$ in $W^{m,2}_{\mathcal{N}}(\Om)$ for some $v_0 \in W^{m,2}_{\mathcal{N}}(\Om)$ by Lemma (\ref{mountain}) and (\ref{Vitali}) holds. We can easily check that $v_0$ is a solution of $(\tilde{P})$ and therefore a critical point of $J_0$. It remains to show that $v_0$ is not a trivial solution. \\
\textbf{Case I.} $c_0 = 0,  v_0 = 0$. We note that $H= \{\|v\|_{W^{m,2}_{\mathcal{N}}(\Om)} = \frac{R_0}{2}\}$ in this case. As $\{v_n\}$ is a $(PS)_{H,c}$ sequence we have $v_n \to 0$ strongly in $W^{m,2}_{\mathcal{N}}(\Om)$. From the fact that $dist(v_n,H) = 0$ and $H$ is closed we conclude that $v_n \in H$ and which implies that $v_0 \in H$ and $v_0$ is different from $0$.\\
\textbf{Case II.} $c_0 \in \left(0,\frac{(4 \pi)^m m!}{2}\right), v_0 =0$. Using the fact that $J_0(v_n) \to c_0$ we see that for given any $\epsilon > 0, \|v_n\|^2_{W^{m,2}_{\mathcal{N}}(\Om)} \leq (4 \pi)^m m! -\epsilon$ for all large $n$.  Let $0 < \delta < \frac{\epsilon}{(4 \pi)^m m!}$ and $q = \frac{(4 \pi)^m m!}{(1+\delta)((4 \pi)^m m!-\epsilon)} >1$. We have
$$\int_{\Om} |\tilde{f}(x,v_n)v_n|^q \leq C \int_{\Om} e^{((1+\delta)q\|v_n\|^2)\left(\frac{v_n^2}{\|v_n\|^2}\right)^2},$$
since $\sup_{x \in \bar{\Om}}|\tilde{f}(x,s)s| \leq C e^{(1+\delta)s^2}$, for all $ s \in \mathbb{R}$, for some $C>0$. Now from the Tarsi's embedding \ref{Tarsi} we get that $\sup_{x \in \bar{\Om}} \int_{\Om} |\tilde{f}(x,v_n)v_n|^q  < \infty$ since $(1+\delta)q\|v_n\|^2 \leq (4 \pi)^m m!$. Also by Vitali's convergence theorem we get  $\int_{\Om} \tilde{f}(x,v_n)v_n \to 0$ as $n \to \infty$ since $v_n \to 0$ pointwise almost everywhere in $\Om$. Which implies
\begin{align*}
 o_n(1) \|v_n\|_{W^{m,2}_{\mathcal{N}}(\Om)} &= \langle J_0'(v_n),v_n \rangle = \frac{1}{2} \int_{\Om}|\nabla v_n|^2 - \int_{Om} \tilde{f}(x,v_n)v_n\\
 &= \frac{1}{2} \int_{\Om} |\nabla^m v_n|^2+ o_n(1)
 \end{align*}
which contradicts the fact $\frac{1}{2} \int_{\Om} |\nabla^m v_n|^2 \to c_0$ as $n \to \infty$. Therefore $v_0$ is not identically $ 0$ in $\Om$. And positivity of $v_0$ comes from the fact that $\tilde{f}(x,s) \geq 0$ for all $(x,s) \in \Om \times \mathbb{R}$ and using the maximum principle.
\end{proof}
\section{Proof of Theorem \ref{Multiplicity}}
Define $\lambda_* = \mu C_0^{\frac{p+3}{p+4}}|\Om|^{-\frac{p+2}{2p+8}}$ where $C_0$ is same as in the Proposition (\ref{lambdastar}). Then condition (\ref{Nehari1}) is true whenever $0<\lambda < \lambda_*$. From the Theorem \ref{1st solution} and \ref{2nd solution} we show the existence of at least two positive solutions for $(P)$. \\
Let $\phi_1$ be the eigen function of $(-\Delta)^m$ on $W^{m,2}_{\mathcal{N}}(\Om)$. Define
$$\lambda^* = p \left( \frac{\lambda_1}{p+1}\right)^{\frac{p+1}{p}} \left( \frac{\int_{\Om}\phi_1}{\int_{\Om}h \phi_1}\right).$$
 We prove that there is no solution of $(P)$ when $\lambda > \lambda^*$. Assume that $u_{\lambda}$ be a solution of $(P)$. By multiplying $\phi_1$ with $(P)$ and performing integration by parts over $\Omega$, we get
\begin{align*}
\int_{\Om} (-\Delta)^m u_{\lambda} \phi_1 = \int_{\Om}f(u_{\lambda})\phi_1 + \lambda \int_{\Om}h \phi_1
\end{align*}
implies
\begin{equation}
\lambda \int_{\Om} h \phi_1 = \int_{\Om}(\lambda_1u_{\lambda}-f(u_{\lambda}))\phi_1                 \label{lambdaupper}
\end{equation}
We see that $\lambda_1t-f(t) \leq \lambda_1-\mu t^{p+1}=\Theta(t)$ for all $t >0$. The global maximum for the function $\Theta$ is $p \left(\frac{\lambda_1}{p+1} \right)^{\frac{p+1}{p}}$ on $(0,\infty)$.
Then from (\ref{lambdaupper}) and the definition of $\lambda^*$ we get $\lambda \leq \lambda^*$. This completes Theorem \ref{Multiplicity}.
\section*{Acknowledgement} The author would like to thank Prof. S Prashanth for several helpful discussions.