text:
\section{\textbf{Introduction}}

Let $(M,g)$ be a closed (compact without boundary) Riemannian manifold of dimension $n\geq 3$. We will consider the Yamabe type equations
\begin{equation}\label{Eq:Main}
-\Delta_{g}u+\lambda u=\mu\left\vert u\right\vert ^{p-1}u\qquad\text{ on }M
\end{equation}
where $\lambda\in\mathcal{C}^\infty(M)$, $\mu\in\mathbb R$  and $p>1$. In case $\lambda=R_g$ is the scalar curvature and $p=p_n:=\frac{n+2}{n-2}$ is the critical Sobolev exponent, equation \eqref{Eq:Main} is  the well known Yamabe equation, widely studied in the last 50 years (see, for example, \cite{Au,BrMa,cfer2} and the references therein). When $p<\frac{n+2}{n-2}$, we will say that the equation \eqref{Eq:Main} is subcritical and we will call it supercritical if $p>\frac{n+2}{n-2}$. In the subcritical case, as the Sobolev embedding $H^1(M,g)\hookrightarrow L^p(M,g)$ is compact, the existence of positive and sign-changing solutions can be obtained using standard variational methods \cite{St,Wi}.  When $M=\Omega$ is a bounded domain of $\mathbb R^{n+1}$ with smooth boundary, there has been recent progress in handling supercritical exponent problems like \eqref{Eq:Main}. A fruitful approach consists in reducing the supercritical problem to a more general elliptic critical or subcritical problem, either by considering rotational symmetries or by means of maps preserving the Laplace operator or by a combination of both, see \cite{ClPi} and the references therein. In case of closed Riemannian manifolds, these reduction methods also apply and have been combined with the Lyapunov-Schmidt reduction method in order to obtain sequences of positive and sign-changing solutions to similar supercritical problems, such that they blow-up or concentrate at minimal submanifolds of $M$ \cite{ClGhMi,DeMuPi,GhMiPi,MiPiVe,PiVa}.

The main interest of this paper is to seek for sign-changing solutions (also called \emph{nodal solutions}) to the problem
\begin{equation}\label{Eq:Yamabe sphere}
-\Delta_{g_0}u+\lambda u=\lambda\left\vert u\right\vert ^{p-1}u\qquad\text{ on }(\mathbb{S}^n,g_0).
\end{equation}
when $p$ is either subcritical or supercritical. Here $g_0$ denotes the round metric and we will assume from now on that $\lambda>0$ is constant. When $p=p_n$, \eqref{Eq:Yamabe sphere} is a renormalization of the Yamabe problem \eqref{Eq:Main} and this kind of solutions have been studied in \cite{Cl,cfer2,dpmpp,dpmpp2,MuWe} and more recently in \cite{fp,MeMuWe,pv}. The slightly subcritical case has been studied in \cite{RoVe}, where the authors obtained multiplicity of nodal solutions blowing-up at points, while the general subcritical case has been studied in \cite{hp} and in \cite{BrLi}. The method introduced in \cite{hp} allowed the authors to obtain more information about the qualitative behavior of the solutions, for they showed the existence of an infinite number of non constant positive solutions having prescribed level sets in terms of isoparametric hypersurfaces. This method has been further generalized to supercritical exponents in \cite{BeJuPe} to prove a similar result on general closed Riemannian manifolds, including the round sphere. Other results concerning the existence and concentration of positive solutions along minimal submanifolds for the supercritical and slightly supercritical can be found in \cite{DeMuPi,MiPiVe}. However, very little is known about the existence, multiplicity and blow-up of nodal solutions for the supercritical problem on the sphere \eqref{Eq:Yamabe sphere}. One of the few results known by the authors is given in \cite{h}, where the existence of at least one sign-changing solution to the supercritical problem was settled. In this direction we will follow and generalize the ideas introduced in \cite{fp} to obtain an infinite number of nodal solutions to the supercritical and subcritical problem \eqref{Eq:Yamabe sphere}, having as level and critical sets isoparametric hypersurfaces and its focal submanifolds.

To state our main result and to describe the method, we briefly recall some aspects of the theory of isoparametric functions and hypersurfaces. For the details, we refer the reader to \cite{bco,cr}. A smooth function $f:(M,g)\rightarrow\mathbb R$ is isoparametric if there exist smooth functions $a,b:\mathbb R\rightarrow\mathbb R$ such that
\begin{equation}\label{Eq:Isoparametric}
\vert \nabla f\vert_g^2=b(f)\quad\text{and}\quad\Delta f=a(f).
\end{equation}
The regular level sets of $f$ are called isoparametric hypersurfaces.

The theory of isoparametric hypersurfaces in the round sphere $(\mathbb S^n,g_0)$ is very rich and it is a vast research topic. In this case, isoparametric hypersurfaces coincide with the hypersurfaces of constant principal curvatures. Its classification began with E. Cartan \cite{c} and it is still an open problem, see \cite{ch,miy2,miy3} and the references therein. Some major progresses in the theory were made by Cartan himself, H. F. M\"{u}nzner \cite{m1,m2} and D. Ferus, H. Karcher and M\"{u}nzner \cite{fkm}. Given an isoparametric hypersurface, there exist a huge number of isoparametric functions having it as level hypersurface, for if $f:\mathbb S^n\rightarrow\mathbb R$ is isoparametric, $\nu:Im(f)\rightarrow\mathbb R$ is monotone and $\alpha\in\mathbb R\smallsetminus\{0\}$, then $\alpha(\nu\circ f)$ is again isoparametric.  However, there are ``canonical'' isoparametric functions, which are obtained by restricting Cartan-M\"{u}nzner polynomials to the sphere \cite[Section 3.5]{cr}. They are well understood and have some nice properties. For instance, if $f:\mathbb S^n\rightarrow\mathbb R$ is obtained in this way, then $\text{Im}f=[-1,1]$, the inverse image of a regular value is a connected isoparametric hypersurface, its only critical values are $t=\pm 1$ and the functions $a$ and $b$ defined in \eqref{Eq:Isoparametric} can be written explicitly. To give these explicit expressions, let $\ell$ be the number of distinct principal curvatures of the level sets of $f$. M\"{u}nzner showed that $\ell\in\{1,2,3,4,6\}$ and if $\ell$ is odd, then all the multiplicities of the principal curvatures are the same, while if $\ell$ is even, there exist, at most, two different multiplicities $m_{-}$ and $m_{+}$ with $1\leq m_{-},m_{+}\leq n-1$. With this notation, if $\Delta_{g_0} f=a(f)$ and $\vert\nabla f\vert_{g_0}=b(f)$, then
\[
a(t)=-\ell(n+\ell-1)t+\frac{\ell^2(m_{+}-m_{-})}{2}\quad\text{and}\quad b(t)=-\ell^2t^2+\ell^2.
\]
The sets $M_{-}:=f^{-1}(-1)$ and $M_{+}:=f^{-1}(1)$ are smooth submanifolds of $\mathbb S^n$ of dimension $n_{-}=(n-1)-m_{-}$ and $n_{+}:=(n-1)-m_{+}$, called focal submanifolds, see \cite[Section 2.4]{cr}. The main feature of these submanifolds is that every isoparametric hypersurface is a tube around $M_{-}$ and $M_{+}$.

If $u$ denotes a sign-changing smooth function defined on a Riemannian manifold, we define the nodal and the critical sets of $u$ to be the sets $\{u=0\}$ and $\{\nabla u = 0\}$, respectively. We state the main result of this paper, which generalizes Theorem 1.3 in \cite{fp}.

\begin{theorem}\label{Th:supercritical} Let $S\subset \mathbb S^n$ be an isoparametric hypersurface and let $n_{-}\leq n_{+}$ be the dimensions of its corresponding focal submanifolds. Then, for any  $\lambda>0$, any $k\in\mathbb N$ and any $p\in (1,\frac{n-n_{+}+2}{n-n_{+}-2})$, equation \eqref{Eq:Yamabe sphere} admits a  nodal solution $u_k$ such that its nodal set has at least $k$ connected components, each of them being  an isoparametric hypersurface diffeomorphic to $S$. The critical set of $u_k$ consists in the focal submanifolds $M_{-}$ and $M_{+}$ and, at least, $k-1$ isoparametric hypersurfaces diffeomorphic to $S$. Moreover, the solutions $u_k$ satisfy
\begin{equation}\label{Eq:Blow-up}
	\lim_{k\rightarrow\infty}\vert u_k(x)\vert=\infty,
\end{equation}
for every $x\in M_-$ or for every $x\in M_+$.
\end{theorem}

Here the numbers $\frac{n-n_{\pm}+2}{n-n_{\pm}-2}\geq p_n$ are just the critical Sobolev exponents in dimensions $n-n_\pm$. Our Theorem improves the existence result stated by Henry in \cite{h}, giving an infinite number of distinct solutions instead of one. It also extends the multiplicity result in \cite{fp} to the subcritical and supercritical exponents. However, this last result gives a better description of the nodal set of the solutions. We strongly believe that a refinement of our methods may give a prescribed number of connected components for the nodal sets of the sign-changing solutions to problem \eqref{Eq:Main}, as in Theorem 1.2 in \cite{fp}.

The last assumption of Theorem \ref{Th:supercritical} says that the sequence $(u_k)$ is not compact with the $C^0$ topology and that the blow-up occurs on one of the focal submanifolds, which are minimal submanifolds of the sphere \cite{cr}. Other noncompactness phenomenon of the same nature appears in the solutions to the critical Yamabe problem obtained in  \cite{dpmpp}, where the blow-up occurs at a single point. However, it was recently proved by Premoselli and V\'{e}tois that this sequence of solutions is uniformly bounded from below, but not from above \cite{pv}. This does not holds true in general, as we state next.

\begin{corollary}\label{Cor:Blowup}
Let $S$ be an isoparametric hypersurface with focal submanifolds $M_-$ and $M_+$ satisfying $\dim M_{\pm}>0$. Then there exists a sequence $(u_k)$ of sign-changing solutions to the Yamabe problem on the sphere
\begin{equation}\label{Eq:CritYamabe sphere}
-\Delta_{g_0}u+\frac{n(n-2)}{4} u=\frac{n(n-2)}{4}\left\vert u\right\vert ^{p_n-1}u\qquad\text{ on }\mathbb{S}^n,
\end{equation}
satisfying that
\[
\lim_{k\rightarrow\infty}u_k(x)=\infty \text{ for every }x\in M_-\ \text{ and } \ \lim_{k\rightarrow\infty}u_k(x)=-\infty \text{ for every }x\in M_+.
\]
\end{corollary}

As another consequence of Theorem \ref{Th:supercritical}, we obtain a multiplicity result for the Yamabe problem on projective spaces.

\begin{corollary}\label{Cor:YamabeProjective}
Let $(M,g)$ be the complex projective space $\mathbb C P^m$ or the quaternionic projective space $\mathbb{H}P^m$ endowed with their canonical metric. Then, if  $j=2$ in case of $\mathbb C P^m$ and $j=4$ in case of $\mathbb H P^m$, for every $k\in\mathbb N$, the Yamabe equation
\begin{equation}\label{Eq:Yamabe projective}
-\frac{4(jm-1)}{jm-2}\Delta_{g}v+R_g v=\left\vert v\right\vert ^{p_{jm}-1}v\qquad\text{ on }(M,g),
\end{equation}
admits a sequence of sign-changing solutions $(u_k)$ such that the regular level sets of $u_k$ consist of isoparametric hypersurfaces in $(M,g)$ and
\begin{equation}\label{Eq:Blow-upProjective}
\lim_{k\rightarrow\infty}\max_{x\in M}\vert u_k(x)\vert=\infty.
\end{equation}
\end{corollary}

We describe briefly the method we shall use in order to prove Theorem \ref{Th:supercritical}. The details will be given in Section \ref{Sec:Double Shooting} and Section \ref{Sec:Projective}.

Let $f:\mathbb S^n\rightarrow\mathbb R$ be an isoparametric function obtained as the restriction of a Cartan-M\"{u}nzner polynomial, and let $\ell$, $m_{-}$ and $m_{+}$ be the number of principal curvatures and the multiplicities associated to the isoparametric hypersurfaces that $f$ defines, as it was explained before. Then, it is easy to see that $z:[-1,1]\rightarrow\mathbb R$ is a solution to the problem
\begin{equation}\label{Eq:Subcritical Yamabe ODE}
b(t)z'' + a(t)z' + \lambda [\vert z\vert^{p-1}z-z] = 0 \quad \text{ on }\  [-1,1],
\end{equation}
with $a(t):=-\ell(n+\ell-1)t+\frac{\ell^2(m_{+}-m_{-})}{2}$ and $b(t):=-\ell^2t^2+\ell^2$, if and only if $u=z\circ f$ is a solution to the problem \eqref{Eq:Yamabe sphere} (Cf. \cite{fp}). Therefore, if $u=z\circ f$ is a solution to \eqref{Eq:Yamabe sphere}, its regular level sets and the set of its critical points are conformed by isoparametric hypersurfaces and focal submanifolds. We can simplify equation \eqref{Eq:Subcritical Yamabe ODE} even more  by considering the new variable $w(r)=z(\cos r)$, and, in this way,  solving \eqref{Eq:Subcritical Yamabe ODE} is equivalent to solving the singular ODE
\begin{equation}\label{Eq:Subcritical radial Yamabe ODE}
w'' + \frac{h(r)}{\sin r} w' + \frac{\lambda}{\ell^2}\left(\vert w\vert^{p-1}w - w\right)=0 \text{ on } [0,\pi],
\end{equation}
where $h(r)=\frac{n-1}{\ell}\cos r - \frac{m_{+}-m_{-}}{2}=\frac{m_{-}+m_{+}}{\ell}\cos r - \frac{m_{+}-m_{-}}{2}$.

Observe that the natural boundary conditions associated to this problem are given by $w'(0)=w'(\pi)=0$. Theorem \ref{Th:supercritical} will be a consequence of the following one.

\begin{theorem}\label{Th:Nodal Yamabe ODE}	For any $p\in (1,\frac{n-n_{-}+2}{n-n_{-}-2})$ and any $k\in\mathbb N$, the equation \eqref{Eq:Subcritical radial Yamabe ODE} with boundary conditions $w'(0)=w'(\pi)=0$ admits a sign changing solution $w_k$ having at least $k$ isolated zeroes in $[0,\pi]$ and at least $k+1$ isolated critical points.
\end{theorem}

We will prove this theorem in Section \ref{Sec:Projective}.

The function $h$ appearing in Equation \eqref{Eq:Subcritical radial Yamabe ODE} has a unique zero $a_0 \in (0,\pi )$. To prove Theorem \ref{Th:Nodal Yamabe ODE} we will use the double shooting method developed in \cite{fp}, which consists in considering the solutions $w_d$, $\widetilde{w}_c$ of
Equation \eqref{Eq:Subcritical radial Yamabe ODE} with initial conditions $w_d'(0) = \widetilde{w}_c '(\pi )=0$,
$w_d (0)=d$, $\widetilde{w}_c (\pi )=c$ and consider the maps $I(d) = (w_d (a_0 ), w_d ' (a_0 ))$ and
$J(c) =(\widetilde{w}_c (a_0 ), \widetilde{w}_c ' (a_0 ))$. If $I(c) = J(d)$, then
$w_d = \widetilde{w}_c$ is a solution of Equation \eqref{Eq:Subcritical radial Yamabe ODE} with $w_d'(0) = w_d '(\pi )=0$, as one can readily see.
To understand the intersections of the curves $I, J$ one needs information of the functions $w_d$, $\widetilde{w}_c$. In the next section we will prove that, for large $d$ and $c$, these functions have many zeroes close to $0$ and $\pi$ (respectively) and then, generalizing an argument based on a Pohozaev-type identity and presented in \cite{cf}, we will prove that $\vert I(d)\vert,\vert J(c)\vert\rightarrow\infty$ as $c,d\rightarrow\infty$. These two results will allow us to conclude that the curves $I$ and $J$ behave as spirals rotating in opposite directions and from this we will obtain the intersections needed to solve the double shooting problem.

\section{\textbf{Double shooting and the proof of Theorem \ref{Th:Nodal Yamabe ODE}.}}\label{Sec:Double Shooting}

We now develop the double shooting method used to prove Theorem \ref{Th:Nodal Yamabe ODE}. First, observe that the function $h$ defined in \eqref{Eq:Subcritical radial Yamabe ODE} satisfies $h(0)=m_{-}$, $h(\pi ) = -m_{+}$, it is strictly decreasing, has a unique zero $a_0\in(0,\pi)$ and $h(r)>0$ in $[0,a_0)$, while $h(r)<0$ in $(a_0,\pi]$. Moreover, the function $\widetilde{h}(r):=-h(\pi-r)=\frac{m_{-} + m_{+}}{2} \cos r +\frac{m_{+} - m_{-}}{2}$ has the same properties with $m_{-}$ and $m_{+}$ interchanged and a unique zero at $\pi-a_0$. To handle both singularities in \eqref{Eq:Subcritical radial Yamabe ODE} at the same time, the strategy is to shoot solutions from each of them and expect that, for some suitable initial and final conditions, the solutions coincide. That is, we consider the initial value problem
\begin{equation}\label{Eq:Singular forward}
\left\{\begin{tabular}{cc}
$w_{i}''(r) + \frac{h(r)}{\sin r} w_i'(r) + \frac{\lambda}{ \ell^2 }( |w_i(r)|^{p-1} w_i -w_i) =0$& in $[0,a_0]$,\\
$w_i(0)=d,\ \  w_i'(0)=0,$ &
\end{tabular}\right.
\end{equation}
and the ``final'' value problem
\begin{equation}\label{Eq:Singular backward}
\left\{\begin{tabular}{cc}
$w_f''(r) + \frac{h(r)}{\sin r} w_f'(r) + \frac{\lambda}{ \ell^2}( |w_f(r)|^{p-1} w_f -w_f) =0$& in $[a_0,\pi]$,\\
$w_f(\pi)=c, \ \  w_f'(\pi)=0,$ &
\end{tabular}\right.
\end{equation}
looking for initial and final conditions $d$ and $c$ such that $w_i(a_0,d)=w_f(a_0,c)$ and $w'_i(a_0,d)=w'_f(a_0,c)$. Hence, by uniqueness of the solution, we would have a well defined solution to problem \eqref{Eq:Subcritical radial Yamabe ODE} given by $w(r)=w_i(r,d)$ if $r\in[0,a_0]$ and $w(r)=w_f(r,c)$ if $r\in[a_0,\pi]$. To construct the solutions with an arbitrarily large number of zeroes, we will need to use that the number of zeroes before and after $a_0$ goes to infinity as $\vert d\vert,\vert c\vert\rightarrow\infty$.

Actually, problem \eqref{Eq:Singular backward} can be written as an initial condition problem having the form of \eqref{Eq:Singular forward}. Indeed, if we consider the function $\widetilde{h}(r)$ defined above, then $w_f$ solves \eqref{Eq:Singular backward} if and only if $\omega(r)=w_f(\pi-r)$ solves the initial value problem
\begin{equation}\label{Eq:Singular backward equivalent}
\left\{\begin{tabular}{cc}
$\omega''(r) + \frac{\widetilde{h}(r)}{\sin r} \omega'(r) + \frac{\lambda}{ \ell^2}( |\omega(r)|^{p-1} \omega -\omega) =0$& in $[0,\pi-a_0]$,\\
$\omega(0)=c,\ \  \omega'(0)=0,$ &
\end{tabular}\right.
\end{equation}
So, in order to understand problem \eqref{Eq:Singular forward} it is enough to consider problem \eqref{Eq:Singular backward}.

In what follows, we will consider the more general equation

\begin{equation}\label{Eq:General singular}
\left\{\begin{tabular}{cc}
$w''(r) + \frac{H(r)}{r} w'(r) + \mu( |w(r)|^{p-1} w -w) =0$ & in $[0,A]$\\
$w(0)=d, \ \ w'(0)=0,$ &
\end{tabular}\right.
\end{equation}
where $A>0$, $\mu>0$ and $H$ is a non negative $C^1$ function defined in the interval $[0,A]$. Notice that equations \eqref{Eq:Singular forward} and \eqref{Eq:Singular backward equivalent} are special cases of the former by taking $\mu=\frac{\lambda}{\ell^2}$, $H(r)=\frac{h(r)r}{\sin r}$ in $[0,A]$ with $A<a_0$ and $H(r)=\frac{\widetilde{h}(r)r}{\sin r}$ in $[0,A]$ with $A<\pi-a_0$. Observe that now we are just dealing with a single singularity at $r=0$.

A standard contraction map argument (Cf. \cite{fp,k}) yields the existence and uniqueness of the solutions to equation \eqref{Eq:General singular} with initial conditions $w(0)=d\in\mathbb R$ and $w'(0)=0$, depending continuously on $d$. For $d> 0$, let $w_d:=w(\cdot,d)$ be the solution with initial values $w_{d } (0) =d$ and  $w_{d } ' (0)=0$. To assure the existence of an arbitrarily large number of zeroes, we use the following result, proven in \cite{fp,h}.

\begin{theorem}\label{Prop:prescribed zeroes}
	Suppose  that $H(0) >0$, $p>1$ and that the following inequality
	\begin{equation}\label{Eq: Subcriticality ODE}
	\frac{H(0) +1 }{2} < \frac{p+1}{p-1}
	\end{equation}
	holds true. Then, for any $\varepsilon >0$ and  any positive integer $k$ there exists $D_k > 0$ so that the solution $w_{d }$ of \eqref{Eq:General singular} has at least $k$ zeroes in $(0,\varepsilon )$ for any $d\geq D_k$.
\end{theorem}

In case of equations \eqref{Eq:Singular forward} and \eqref{Eq:Singular backward equivalent}, we have that $H(0)=h(0)=m_{-}$ and $H(0)=\widetilde{h}(0)=m_{+}$ respectively. Taking $1\leq m_{-}\leq m_{+}\leq n-1$, and recalling that $n_{-}=(n-1)-m_{-}$ and $n_{+}=(n-1)-m_{+}$,  inequality
\begin{equation}\label{Eq:Key inequality}
1<p<\frac{n-n_{+}+2}{n-n_{+}-2}=\frac{m_{+}+3}{m_{+}-1}\leq \frac{m_{-}+3}{m_{-}-1}=\frac{n-n_{-}+2}{n-n_{-}-2}
\end{equation}
implies the validity of inequality \eqref{Eq: Subcriticality ODE} for both equations \eqref{Eq:Singular forward} and \eqref{Eq:Singular backward equivalent}.

Since $p_n\leq \frac{m_{+}+3}{m_{+}-1}$ is true for every $1\leq m_{+}\leq n-1$, we may guarantee the existence of an arbitrary large number of zeroes for equations \eqref{Eq:Singular forward} and \eqref{Eq:Singular backward equivalent} when $p<p_n$, corresponding to the subcritical  Yamabe problem, or when $p_n<p<\frac{m_{+}+3}{m_{+}-1}$, corresponding to the supercritical one. Also observe that inequality \eqref{Eq:Key inequality} is true when $p=p_n$ and $1\leq m_{-},m_{+}<n-1$, a fact used in \cite{fp} to assure the existence of a prescribed number of zeroes to equation \eqref{Eq:Subcritical radial Yamabe ODE} in the critical case. For the rest of this section, we will suppose that $p>1$ satisfies inequality \eqref{Eq:Key inequality}.

Even if the number of zeroes is arbitrarily large, it can  not be infinite as we next show.

\begin{lemma}\label{Cor:No zeroes concentration}
	For $d_\ast>0$, $d_\ast\neq 1$ fixed, the zeroes and critical points of $w_{d_\ast}$ are isolated in $[0,A]$.
\end{lemma}

\begin{proof}
First observe that if $r_0\in[0,A]$ is a zero of $w_{d_\ast}$, then uniqueness of the solution implies that $w_{d_\ast}(r_0)\neq 0$. Therefore $w_{d_\ast}$ is monotone in a neighborhood of $r_0$ and it is an isolated zero. Now, to see that the critical points are isolated, suppose that $w'_{d_\ast}(r_0)=0$ for some $r_0\in[0,A]$. As $d\neq 0,1$, by uniqueness of the solutions, $w_{d_\ast}(r_0)\neq 0$ and this together with equation \eqref{Eq:General singular} implies that $w_{d_\ast}''(r_0)\neq 0$. Without loss of generality, suppose $w''_{d_\ast}>0$. By continuity, there exists $\varepsilon>0$ such that $w''_{d_\ast}(r)>0$ in $(r_0-\varepsilon,r_0+\varepsilon)\cap[0,A]$. Hence $w'_{d_\ast}$ is monotone in $(r_0-\varepsilon,r_0+\varepsilon)\cap[0,A]$ and, therefore, $w_{d_\ast}'(r)\neq 0$ in $(r_0-\varepsilon,r_0+\varepsilon)\cap[0,A]\smallsetminus\{r_0\}$ and $r_0$ is an isolated critical point.
\end{proof}

For $\varepsilon<A$, take $D_k>0$ as given in Theorem \ref{Prop:prescribed zeroes}. Then, for  $d \geq D_k>0$, $w_d$ has at least $k$ zeroes before $A$. Denote them by $r_1(d)<r_2(d)<\ldots < r_k(d)$.  The following statement holds true.

\begin{lemma}\label{Lemma:Limit zeroes}
	For each $j=1,\ldots,k$,
	\[
	\lim_{d\rightarrow\infty}r_j(d)=0
	\]
\end{lemma}

\begin{proof}
Let
\[
z_{d } (r): = d ^{-\frac{2}{p-1}} \   w_{d ^{\frac{2}{p-1}}}\left( \frac{r}{d   \sqrt{\lambda} } \right).
\]
For any compact subset $K$ of $[0,A]$, the functions $z_d$ converge $C^1$-uniformly on $K$ to the unique solution of the limit Cauchy problem
\begin{equation}\label{Eq:Limit problem}
\left\{\begin{tabular}{cc}$v''(r) + \frac{H(0)}{r} v'(r) +  |v(r)|^{p-1} v  =0$ & \text{ in } $[0,\infty)$,\\
$v(0)=1,\ \ v'( 0)=0$ & \end{tabular}\right.
\end{equation}
and as $p>1$ satisfies inequality \eqref{Eq:Key inequality}, $v$ has an infinite number of isolated zeroes in $(0,\infty)$
(see Lemma 3.2 and Theorem 3.3 in \cite{fp} and Proposition 3.6 in \cite{hw}).

Now observe that
	\[
	w_d(r)=dz_{d^{\frac{p-1}{2}}}\left(\sqrt{\lambda}d^{\frac{p-1}{2}r}\right)\qquad\text{and}\qquad w'_d(r)=\sqrt{\lambda}d^{\frac{p+1}{2}}z'_{d^{\frac{p-1}{2}}}\left( \sqrt{\lambda}d^{\frac{p-1}{2}}r \right).
	\]
Therefore $r$ is a zero of $w_d$ if and only if $\sqrt{\lambda}d^{\frac{p-1}{2}}r$ is a zero of $z_{d^{\frac{p-1}{2}}}$. For each $j=1,\ldots,k$, denote by $r^\ast_j(d)$ the $j$-th zero of $z_{d^{\frac{p-1}{2}}}$ and by $a_j$ the $j$-th zero of the solution $v_0$ to the limit problem \eqref{Eq:Limit problem}. Then $r_j^\ast(d)=\lambda d^{\frac{p-1}{2}}r_j(d)$ and the $C^1$ uniform convergence of $z_d$ to $v_0$ implies that $r_j^\ast(d)\rightarrow a_j$ as $d\rightarrow\infty$ for each $j=1,\ldots,k$. Hence $r_j(d)\rightarrow 0$ as $d\rightarrow\infty$.
\end{proof}

Now we focus on equation \eqref{Eq:Singular forward} in order to do a phase plane analysis. Let $a_0$ be the unique zero of $h$ in $(0,\pi )$. Let $w_d$ be the solution of \eqref{Eq:Singular forward}
with initial conditions $w_d (0)=d$, $w_d' (0)=0$.

Consider the curve $I: \mathbb{R} \rightarrow \mathbb{R}^2$ given by $I(d) = (w_d (a_0 ) , w_d ' (a_0 ) )$. Note that
$I(0) = (0,0)$, $I(-d) =-I(d)$ and $I(d) \neq (0,0)$ if $d\neq 0$. It is then easy to see that we have
a well defined continuous function $\theta : (0, \infty ) \rightarrow \mathbb{R}$ such that $\theta (1)=0$
and $\theta (d)$ gives an angle between $I(d)$ and the positive $x$-axis for any $d>0$. Note that, in a similar way, there is a unique continuous function  $\theta : (-\infty ,0) \rightarrow \mathbb{R}$ such that $\theta (-1) = -\pi$
and $\theta (d)$ gives an angle between $I(d)$ and the positive $x$-axis. Thus, we have that for any
$d>0$, $\theta (-d) = \theta (d) - \pi $. Also notice that $w_d (a_0 )=0$ if and only if $\theta (d) = -\frac{\pi}{2} - k \pi$ for some  integer $k$.

Next we proceed to define a second curve in the phase space, corresponding to the solutions to problem \eqref{Eq:Subcritical radial Yamabe ODE} in $[0,\pi]$ with condition $w'(\pi)=0$. Let $\widetilde{h}(r)=-h(\pi-r)=\frac{m_{-}+m_{+}}{2}\cos r + \frac{m_{+}-m_{-}}{2}$ and consider the  problem \eqref{Eq:Singular backward equivalent}. If $\omega$ is a solution to this problem, then $\widetilde{w}(r):=\omega(\pi-r)$ solves the ``final'' conditions problem \eqref{Eq:Singular backward}. For $c\in\mathbb R$, denote by $\widetilde{w}_c$ the solution to the problem \eqref{Eq:Singular backward} and define the map $J(c):=(\widetilde{w}_c(a_0),\widetilde{w}_c'(a_0))$.

In an entirely similar way, $J(1)=(1,0)$, $J(0)=(0,0)$, $J(c) \neq (0,0)$ if $c\neq 0$ and $J(-c)=-J(c)$. So, there is a well defined angle $\vartheta:\mathbb R\smallsetminus\{0\}\rightarrow\mathbb R$ such that $\vartheta(1)=0$ and
\[
J(c)=(\vert J(c)\vert\cos(\vartheta(c)),\vert J(c)\vert\sin(\vartheta(c)).
\]
Note that $\theta$ and $\vartheta$ run in opposite directions. If $n(d)$ denotes the number of zeroes of $w_d$ before $a_0$ and $N(c)$ the number of zeroes of $\widetilde{w}_c$ after $a_0$, then the angles $\theta$ and $\vartheta$ are related with these numbers by the following formulas, proved in \cite{fp},
\begin{equation}\label{Eq:FormulaZeroes}
n (d) = -\left\lfloor\frac{\theta(d) - \pi /2}{\pi} \right\rfloor -1\quad\text{ and }\quad N(c)= - \left\lfloor \frac{-\vartheta(c) - \pi /2}{\pi}\right\rfloor  -1,
\end{equation}
where, as usual for $x\in\mathbb R$, $\lfloor x\rfloor$ denotes the maximum integer such that $\lfloor x\rfloor \leq x$. As a consequence of these formulas and Theorem \ref{Prop:prescribed zeroes}, we have the following limits
\begin{equation}\label{Eq:limit argument}
\lim_{d\rightarrow\infty}\theta(d)=-\infty\quad\text{ and }\quad\lim_{c\rightarrow\infty}\vartheta(c)=\infty.
\end{equation}

We will show that the curves $I$ and $J$ behave like spirals turning in opposite directions. The above limits show that the spirals actually turn. Next we see that the spirals are not bounded. If $w_d$ is a solution to the initial value problem \eqref{Eq:Singular forward} and $w_c$ is a solution to \eqref{Eq:Singular backward}, define
$\rho(r,d):=\sqrt{w_d^2(r)+(w_d')^2(r)}$ and $\varrho(r,c)=\sqrt{w_c^2(r)+(w_c')^2(r)}.$

\begin{lemma}\label{Lemma:Radius}
We have that
	\[
	\lim_{d\rightarrow\infty}\rho(r,d)=\infty,\quad\text{and}\quad 	\lim_{c\rightarrow\infty}\varrho(r,c)=\infty,
	\]
uniformly in $[0,a_0]$ and in $[a_0,\pi]$, respectively.
	\end{lemma}

The proof of this Lemma is technical and we postpone it to the Appendix \ref{Sec:Energy Analysis}.

Consider the functions $\rho(d):=\rho(d,a_0)$, $\varrho(c):=\varrho(c,a_0)$, $\theta(d):=\theta(d,a_0)$ and $\vartheta(c):=\vartheta(c,a_0)$.

Now, define in the radius-argument plane the curves $R:[1,\infty)\rightarrow \mathbb R\times\mathbb R_{>0}$ and $S:[1,\infty)\rightarrow\mathbb R\times\mathbb R_{>0}$ given by
\[
R(d):=(\theta(d),\rho(d)),\quad\text{ and }\quad S(c):=(\vartheta(c),\varrho(c))
\]
From uniqueness of the solutions to problems \eqref{Eq:Singular forward} and \eqref{Eq:Singular backward}, the curves $R$ and $S$ are simple and they intersect at the point $(0,1)$.

For each $i,j\in\mathbb N$, define
\[
d_i:=\max\{d\;:\; \theta(d)=-i\pi\}\quad \text{ and }\quad c_j:=\max\{c\;:\;\vartheta(c)=j\pi\},
\]
and
\[
\widehat{d}_i:=\min\{d\;:\; \theta(d)=-i\pi\}\quad \text{ and }\quad \widehat{c}_j:=\min\{c\;:\;\vartheta(c)=j\pi\},
\]
These numbers are well defined by \eqref{Eq:limit argument} and they form unbounded sequences by the same limit.  Observe that $\widehat{d}_i$ and $d_i$ are first and last time that the curve $R$ hits the line  $\theta=-i\pi$, respectively, while $\widehat{c}_i$ and $c_j$ correspond to the first and last time that $S$ hits the line $\vartheta=j\pi$. It was also shown in \cite{fp} that for any $c,d>0$, $\theta(d)<\pi/2$ and $\vartheta(c)>-\pi/2$. So, it follows that the curve $R$ is completely contained in $(-\infty,\pi/2)\times\mathbb R_{>0}$, while $S$ is contained in $(-\pi/2,\infty)\times\mathbb R_{>0}$, and that $R$ restricted to $[d_1,\infty)$ and $S$ restricted to $[c_1,\infty)$ do not intersect.

We can now prove Theorem \ref{Th:Nodal Yamabe ODE}.

\begin{proof}[Proof of Theorem \ref{Th:Nodal Yamabe ODE}.]
Let $k\in\mathbb N$ and for each $i,j\in\mathbb N$, set $x_i:=\rho(d_i)$, $\widehat{x}_i:=\rho(\widehat{d}_i),$ $y_j:=\varrho(c_j)$ and $\widehat{y}_j:=\varrho(\widehat{c}_j)$.
By Lemma \ref{Lemma:Radius}, these sequences are unbounded. Therefore, we can find $i,j>k$ and $\alpha>j$ such that $y_j<\min\{x_i,\widehat{x}_i,x_{\alpha+1-j}\}<\widehat{y}_\alpha$. Observe that the curves $R$ and $S-((\alpha+j)\pi,0)$ restricted to the intervals $[d_i,\widehat{d}_{\alpha+i-j}]$ and $[c_j,\widehat{c}_\alpha]$, respectively, are both contained in $\mathcal{A}_k:=[-(\alpha+i-j)\pi,-i\pi]\times[y_j,\widehat{y}_\alpha]$. As $\widehat{x}_{\alpha+i-j}>y_j$ and $x_i<\widehat{y}_{\alpha}$ and as $R$ restricted to $[d_j,\widehat{d}_{\alpha+j-i}]$ intersects $\mathcal{A}_k$ only at the points $(-(\alpha+j-i)\pi,\widehat{x}_{\alpha+j-i})$ and $(-j\pi,x_j)$, the intermediate value Theorem implies that the curve $R$ must intersect $S-((\alpha+j)\pi,0)$. Let $d_R> 1$ and $c_S>1$ be the points such that $R(d_R)=S(c_S)-((\alpha+j)\pi,0)$. Using the formulas \eqref{Eq:FormulaZeroes}, we can argue as in the proof of Theorem 1.2 in \cite{fp} to conclude that $w_{d_R}=\widetilde{w}_{c_S}$ is a solution to the problem \eqref{Eq:Subcritical radial Yamabe ODE} with exactly $\alpha+j>k$ zeroes and, since $w'_{d_R}(0)=w'_{d_R}(\pi)=0$, it has, at least, $k+1$ critical points by Rolle's Theorem. The fact that the zeroes and that the critical points are isolated follows from Lemma \ref{Cor:No zeroes concentration}.
\end{proof}

This theorem implies \ref{Th:supercritical} as follows.

\begin{proof}[Proof of Theorem \ref{Th:supercritical}.]
Associated to $S$, there is a Cartan-M\"{u}nzner polynomial such that its restriction to the sphere is an isoparametric function $f:\mathbb S^n\rightarrow[-1,1]$. Then equation \eqref{Eq:Yamabe sphere} can be reduced into equation \eqref{Eq:Subcritical radial Yamabe ODE}. By Theorem \ref{Th:Nodal Yamabe ODE}, this equation admits a sign-changing solution $w_k$ having at least $k$ isolated zeroes and at least $k+1$ isolated critical points in $[0,\pi]$. Therefore, $u_k=w_k(\arccos f)$ is a solution to the problem \eqref{Eq:Yamabe sphere} having as regular level sets a disjoint union of connected isoparametric hypersurfaces diffeomorphic to $S$. As $w_k$ has at least $k$ isolated zeroes, then the nodal set $u_k^{-1}(0)$ has at least $k$ connected components, all of them being isoparametric hypersurfaces diffeomorphic to $S$. As $w_k$ has at least $k-1$ isolated critical points in $(0,\pi)$ and as $w_k'(0)=0=w_k'(\pi)$, the critical set of $u_k$ consists in, at least, $k-1$ connected isoparametric submanifolds diffeomorphic to $S$, together with the focal submanifolds $M_{-}=f^{-1}(-1)$ and $M_{+}=f^{-1}(1)$. Finally, by construction of $w_k$, we have that $w_k(0)\rightarrow\infty$ or $\vert w_k(\pi)\vert \rightarrow\infty$ as $k\rightarrow\infty$, for the number of zeroes increases as the initial or final conditions increase. Suppose, without loss of generality, that $w_k(0)\rightarrow\infty$. In this case, as $\arccos\left( f(x)\right)=0$ for all $x\in M_-$ we concluding the limit \eqref{Eq:Blow-up} for every $x\in M_-$.
\end{proof}

Even if Theorem \ref{Th:supercritical} holds for subcritical, critical and supercritical values of $p>1$, the methods developed here do not allow us to prove the existence of solutions such that the final value $w_d(\pi)=c$ is negative. However, in case of the critical exponent $p_n$, we can use the refined version of this theorem given in \cite{fp} to construct a sequence of solutions which is not uniformly bounded from below.

\begin{proof}[Proof of Corollary \ref{Cor:Blowup}]
As the focal manifolds that generate $S$ have positive dimensional, the number of distinct principal curvatures $\ell$ must be bigger that $1$. Let $f$ be the Cartan-M\"{u}nzner polynomial associated to $S$. Then, for any $k\in\mathbb N$, Theorem 1.2 in \cite{fp} gives a solution to the Yamabe problem on the sphere \eqref{Eq:CritYamabe sphere} having the form $u_k=\arctan(f(w_{k}))$, where $w_{k}$ is a solution to the problem \eqref{Eq:Subcritical radial Yamabe ODE}, with $p=p_n$, having exactly $k$ zeroes in $[0,\pi]$. Lemma \ref{Lemma:Radius} together with Lemma 4.6 in \cite{fp} imply that the sequences $(x_i)$ and $(y_j)$ are both increasing and unbounded. In this way, the situation in Lemma 4.7 in \cite{fp} can not happen and the number of zeroes before and after $a_0$ must increase as the initial and final conditions, $w_k(0)$ and $\vert w_k(\pi)\vert$ respectively, increase. For $k$ odd, by Lemma \ref{Cor:No zeroes concentration} and its proof, necessarily $w_k(0)>0$ and $w_k(\pi)<0$ and $w_k$ can not have an infinite number of zeroes before and after $a_0$. Therefore there exists  a subsequence $(w_{k_j})$, with each $k_j$ odd, such that $w_{k_j}(0)\rightarrow\infty$ and $w_{k_j}(\pi)\rightarrow-\infty$ as $j\rightarrow\infty$. Since $\arctan(f(x))=0$ for every $x\in M_+$ and $\arctan(f(x))=\pi$ for every $x\in M_-$, the sequence $(u_{k_j})$ satisfies the desired limits as $j\rightarrow\infty$.
\end{proof}

\section{Submersion with minimal fibers and the proof of Corollary \ref{Cor:YamabeProjective}.}\label{Sec:Projective}

We shall study Riemannian submersions of the form $\pi:(\mathbb S^n,g_0)\rightarrow (M^m,g)$, where $(M,g)$ is a closed Riemannian manifold with dimension $m<n$, and such that the fibers are minimal. It is well known that under these conditions, $\pi$ is a harmonic morphism with dilation $\varphi\equiv 1$, see \cite{bw} or \cite[Section 11.2]{f}. In this case, a function $v:M\rightarrow\mathbb R$ solves the equation
\begin{equation}\label{Eq:Reduced equation}
-\Delta_g v + \lambda v = \mu \vert v\vert^{p-1}v,\quad\text{on }M
\end{equation}
with $\lambda,\mu>0$, if and only if $u=\left[ \frac{\mu}{\lambda}\right]^{\frac{1}{p-1}}v\circ \pi$ is a solution to the Yamabe problem on the sphere \eqref{Eq:Yamabe sphere} (Cf. \cite{cfp}). Observe that if we are considering the critical exponent $p_n$ for the equation \eqref{Eq:Yamabe sphere}, then this same exponent is subcritical in equation \eqref{Eq:Reduced equation} for, in this case, $m<n$ implies $p_n<p_m$.

We next show that isoparametric functions are preserved by Riemannian submersions with minimal fibers.

\begin{lemma}\label{Lemma:Isoparam functions under Riemannian submersions}
	Let $\pi:(N,h)\rightarrow(M,g)$ be a Riemannian submersion with minimal fibers. If a smooth function $\widehat{f}:M\rightarrow\mathbb R$ is isoparametric, then $f:=\widehat{f}\circ\pi:N\rightarrow\mathbb R$ is isoparametric. Conversely, if $f:N\rightarrow \mathbb R$ is isoparametric and $f(x)=f(y)$ whenever $\pi(x)=\pi(y)$, then there exists a smooth function $\widehat{f}:M\rightarrow \mathbb R$ such that $f=\widehat{f}\circ\pi$ and $\widehat{f}$ is isoparametric.
\end{lemma}

\begin{proof}
	Suppose $\widehat{f}:M\rightarrow\mathbb R$ is isoparametric with $\vert \nabla \widehat{f}\vert_g^2=\widehat{a}(\widehat{f})$ and $\Delta \widehat{f}=\widehat{b}(\widehat{f})$. On the one hand, as $\pi:N\rightarrow M$ is a harmonic morphism with dilation $\varphi\equiv 1$ (see \cite{bw}), we have that
	\[
	\Delta_h f=\Delta_h (\widehat{f}\circ\pi)=(\Delta_g \widehat{f})\circ\pi=(\widehat{b}\circ \widehat{f})\circ\pi=\widehat{b}(\widehat{f}\circ\pi)=\widehat{b}(f).
	\]

	Next, observe that
	\[
	\nabla \widehat{f} = \pi_\ast\nabla f
	\]
	because $\pi:N\rightarrow M$ is a Riemannian submersion.

	Therefore,
	\begin{align*}
	&\vert \nabla f(x)\vert_h^2 = h(\nabla f(x),\nabla f(x))=g(\pi_\ast\nabla f(x),\pi_\ast\nabla f(x))\\
	&=g(\nabla \widehat{f}(\pi(x)),\nabla \widehat{f}(\pi(x)))=\vert \nabla \widehat{f}(\pi(x))\vert_g^2=a(\widehat{f}(\pi(x)))=\widehat{a}(f(x)),
	\end{align*}
	and $f$ is an isoparametric function.

	Now suppose that $f$ is an isoparametric function such that $f(x)=f(y)$ if $\pi(x)=\pi(y)$. Let $a,b:\mathbb R\rightarrow\mathbb R$ be such that $\vert \nabla f\vert_h^2=a(f)$ and $\Delta_h f=b(f)$. As $\pi:N\rightarrow M$ is a submersion, it is a quotient map \cite{lee} and the function $f$ passes to the quotient as an smooth function $\widehat{f}:M\rightarrow \mathbb R$ such that $f=\widehat{f}\circ\pi$. Then, as before
	\[
	a\circ \widehat{f}\circ\pi(x)=a(f(x))=\vert \nabla f(x)\vert_n^2=\vert \nabla \widehat{f}(\pi(x))\vert_g^2
	\]
	and
	\[
	b\circ \widehat{f}\circ\pi=b(f)=\Delta_h f= \Delta_h(\widehat{f}\circ\pi)=(\Delta_g \widehat{f})\circ\pi.
	\]
	Since $\pi$ is surjective, it has a right inverse and we conclude that $\vert \nabla \widehat{f}\vert_g^2=a(\widehat{f})$ and $\Delta_g \widehat{f}=b(\widehat{f})$.
\end{proof}

Let $(M,g)$ denote $(\mathbb C P^m,\tilde{g}_{0})$ or $(\mathbb H P^m,\tilde{g}_0)$, the complex and quaternionic spaces with their corresponding canonical metrics. Both of them are Einstein manifolds with constant positive scalar curvature (see \cite[Theorem 14.39]{b} and Sections 8.1 and 8.3 in \cite{p}). We will denote it by $R_{\tilde{g}}=\Lambda_m>0$. Recall that $\mathbb{C}P^m=\mathbb S^{2m+1}/\mathbb S^1$ and that $\mathbb H^m=\mathbb S^{4m+3}/SU(2)$, where $\mathbb S^1$ is the circle group and $SU(2)$ is the group of unit quaternions. We need the following lemma.

\begin{lemma}\label{Lemma:ExistenceIsoparametric}
We have the following:
\begin{enumerate}
	\item For each $m\geq 3$, there exists an $\mathbb S^1$-invariant Cartan-M\"{u}nzner polynomial $f:\mathbb S^{2m+1}\rightarrow\mathbb R$ such that the associated focal manifolds $M_{-}$ and $M_{+}$ satisfy $\dim M_{\pm}\geq 2$.
	\item For each $m\geq 3$ there exists an $SU(2)$-invariant Cartan-M\"{u}nzner polynomial $f:\mathbb S^{4m+3}\rightarrow\mathbb R$ such that the associated focal manifolds $M_{-}$ and $M_{+}$ satisfy $\dim M_\pm\geq 4$.
\end{enumerate}
\end{lemma}

\begin{proof}
We begin with the proof of (1). In this case, we consider $\mathbb S^{2m+1}\subset\mathbb R^{2m+2}$. As $m\geq 3$, we can write $m=2+k$ with $k\in\mathbb N\smallsetminus\{0\}$ and we can decompose $\mathbb R^{2m+2}\equiv\mathbb R^{2\alpha}\times\mathbb R^{2\beta}\equiv\mathbb C^\alpha\times\mathbb C^\beta$ where $\alpha=\beta=2$ if $k=1$ and $\alpha=k$, $\beta=3$ if $k\geq 2$. Hence, we have an action of $\mathbb S^1$ by isometries given as
$\zeta (x,y):=(\zeta x,\zeta y)$, where $\zeta(x,y):=(\zeta x_1,\ldots,\zeta x_{\alpha},\zeta y_1,\ldots,\zeta y_{\beta})$, with $\zeta\in\mathbb S^1\subset\mathbb C$ and $x_i,y_j\in\mathbb C$. We can then consider the degree two Cartan-M\"{u}nzner polynomial $f:\mathbb R^{2m+2}\equiv\mathbb C^\alpha\times\mathbb C^\beta\rightarrow\mathbb R$, $f(x,y)=\vert x\vert^2 - \vert y\vert^2$, which is clearly $\mathbb S^1$-invariant. The restriction of $f$ to $\mathbb S^{2m+1}$ is an isoparametric function with focal submanifolds $M_{-}=\mathbb S^{2\alpha-1}\times \{0\}$ and $M_{+}=\{0\}\times\mathbb S^{2\beta-1}$, see \cite{cr}. If $k=1$, $\dim M_\pm=3$ while if $k\geq 2$, $\dim M_-=2k-1\geq 3$ and $\dim M_+=5$, and the lemma follows.

Now we prove (2) in a similar way. For each $m\geq 4$, write $m=4+k$ with $k\in\mathbb N\cup\{0\}$. Hence $\mathbb S^{4m+3}\subset \mathbb R^{4m+4}\equiv\mathbb R^{4\alpha}\times\mathbb R^{4\beta}\equiv\mathbb H^{\alpha}\times\mathbb H^{\beta}$, where $\alpha=\beta=4$ if $m=3$ and $\alpha=k+2$ and $\beta=3$ if $m\geq 4$. Hence we have a natural action of $SU(2)$ on $\mathbb R^{4m+4}$ given by $\zeta(x,y)=(\zeta x,\zeta y)=(\zeta x_1,\ldots,\zeta x_{\alpha},\zeta y_1,\ldots,\zeta y_{\beta})$ where $\zeta\in SU(2)$ and $x_i,y_j\in \mathbb H$. Since $\vert \zeta q\vert=\vert \zeta\vert\vert q\vert=\vert q\vert$ for each $q\in \mathbb H$ and $\zeta\in SU(2)$, we have that the Cartan-M\"{u}nzner polynomial $f:\mathbb R^{2m+2}\equiv\mathbb H^\alpha\times\mathbb H^\beta\times\mathbb R^{2m+2}\rightarrow\mathbb R$ given by $f(x,y)=\vert x\vert^2 - \vert y\vert^2$ is $SU(2)$-invariant. As before, the corresponding focal submanifolds are given by $M_{-}=\mathbb S^{4\alpha-1}\times\{0\}$ and $M_{+}=\{0\}\times\mathbb S^{\beta - 1}$. Thus, $\dim M_\pm= 7$ if $m=3$, while $\dim M_-=4(k+2)-1\geq 11$ and $\dim M_+=11$ for every $m\geq 4$, where we conclude the lemma in this case.
\end{proof}

\begin{proof}[Proof of Corollary \ref{Cor:YamabeProjective}]
We analize the case of the quaternionic projective space, being the case for $\mathbb C P^m$ completely analogous. First we write \eqref{Eq:Yamabe projective} in the form \eqref{Eq:Reduced equation},
\begin{equation}\label{Eq:YamabeQuaternionic}
-\Delta_{g}v+\lambda_m v =\mu_m \vert v\vert^{p_{4m}-1}\qquad \text{ on }\mathbb H P^m,
\end{equation}
where $\lambda_m:=\frac{\Lambda_m(4m-2)}{16m-4}$ and $\mu_m:=\frac{4m-2}{16m-4}$. Let $\pi:\mathbb S^{4m+3}\rightarrow \mathbb H P^m$ be the natural projection. This map is a Riemannian submersion with minimal fibers and, so, it is a harmonic morphism \cite{bw}.

Therefore, $v$ is a solution to the Yamabe problem
\eqref{Eq:Yamabe projective} if and only if $u=\left[\frac{\mu_m}{\lambda_m}\right]^{\frac{1}{p_{4m}-1}}v\circ\pi$ is a solution to the supercritical problem on the sphere
\begin{equation}\label{Eq:SupercriticalSphereQuaternionic}
-\Delta_{g_0}u+\lambda_{m}u=\lambda_m\vert u\vert^{p_{4m}-1}u\qquad\text{on }\mathbb S^{4m+3}.
\end{equation}
By Lemma \ref{Lemma:ExistenceIsoparametric}, there exists an $SU(2)$-invariant isoparametric function $f:\mathbb S^{4m+3}\rightarrow\mathbb R$, which is the restriction of a Cartan-M\"{u}nzner polynomial and such that its corresponding focal submanifold have dimension at least $4$. Hence, if we take $\kappa = \min\{\dim M_-,\dim M_+\}\geq 4$, then
\[
p_{4m}<\frac{(4m+3-\kappa)+2}{(4m+3-\kappa)-2}
\]
and the supercritical problem \eqref{Eq:SupercriticalSphereQuaternionic} admits a sequence of nodal solutions $(u_k)$ of the form $u_k=z_k\circ f$ where $z_k$ is a solution to the problem \eqref{Eq:Subcritical Yamabe ODE} with at least $k$ zeroes. As $f$ is $SU(2)$-invariant, there exists $\widehat{f}:\mathbb H P^{m}\rightarrow\mathbb R$ such that $f=\widehat{f}\circ \pi$. Therefore $u_k=z_k\circ f=(z_k\circ\widehat{f})\circ\pi$ is a solution to \eqref{Eq:SupercriticalSphereQuaternionic}, implying that $v_k:=\left[\frac{\mu_m}{\lambda_m}\right]^{-\frac{1}{p_{4m}-1}}z_k\circ\widehat{f}$ is a solution to the Yamabe problem \eqref{Eq:Yamabe projective}. Moreover, by Lemma \ref{Lemma:Isoparam functions under Riemannian submersions}, $\widehat{f}$ is an isoparametric function and, thus, $v_k$ has isoparametric hypersurfaces in $\mathbb H P^{m}$ as level sets.
Finally, the blow-up of the sequence $(v_k)$ is a consequence of the limit \eqref{Eq:Blow-up} for the sequence $(u_k)$ at $M_-$ or at $M_+$.
\end{proof}


text:
\section{Introduction}
The electrodynamic lattice-Boltzmann method (ELBM) has been implemented for modeling interactions of electromagnetic waves with uniform dielectric \cite{Mendoza2010,Hanasoge2011,Hauser2017,Hauser2019,hauser2022simulation}, conductive \cite{Mendoza2010,Hauser2017}, non-linear and dispersive media \cite{Mendoza2011,Chen2013}. These methods can have superior accuracy and/or computational efficiency compared with the finite-difference-time-domain (FDTD) method \cite{KaneYee1966,Mendoza2010,Hauser2019}. The ELBM is advantageous to the FDTD method in its use of a single lattice, as opposed to an offset Yee-lattice \cite{KaneYee1966}. The Hauser \& Verhey (HV) ELBM has the additional capability of representing multi-dimensional propagation of electromagnetic waves which are free of numerical dispersion in free space \cite{Hauser2019}. Previously, both Chen \emph{et al.} and Hauser \& Verhey have used piece-wise linear recursive convolution (PLRC) methods to model the frequency dispersion of electromagnetic waves \cite{Chen2013,Hauser2021,hauser2022simulation}. The HV ELBM is second-order accurate, demonstrating an improvement over Chen \emph{et al.}'s ELBM which is first-order accurate \cite{Chen2013}, while remaining amenable to an ADE approach \cite{Warner2021}. Compared with the PLRC approach, which is typically used in the ELBM, an auxiliary differential equation (ADE) approach has demonstrated greater stability and efficiency for evaluating second-order accurate frequency dispersion of electromagnetic waves in the FDTD method \cite{Okoniewski1997,Chun2013}. The ADE with multiple complex-conjugate pole-residue pairs (CCPRPs)\cite{MinghuiHan2006} may also describe Debye and Lorentz media in a unified manner, reducing its implementation complexity from the PLRC method which treats these two types of media in an independent manner\cite{Chun2013,Chen2013}. In a previous conference paper \cite{Warner2021}, we integrated a second-order central difference of the ADE using multiple CCPRPs\cite{Okoniewski1997,MinghuiHan2006} within a one-dimensional (1D) HV ELBM\cite{Hauser2019} to model the broadband frequency-dispersion of arbitrary electromagnetic waves transmitted through dispersive silver \cite{MinghuiHan2006}. This study elaborates on our conference paper by describing how a central difference CCPRP ADE dispersive model designed for the FDTD method \cite{Okoniewski1997,MinghuiHan2006} can be incorporated into the HV ELBM \cite{Hauser2019} without modifying the ADE model or the HV ELBM. It further compares the 1D HV ELBM and the 1D FDTD method's relative errors, order-of-accuracy, computational time, and computational resources for evaluating the analytical transmittance of evaporated silver from Palik's measurements \cite{Lynch1985,MinghuiHan2006,Chen2013}. The relative error and computational time comparisons are drawn for the simplest case of steady-state continuous wave transmission through an evaporated silver slab of 100 \si{nm} thickness \cite{Lynch1985,MinghuiHan2006,Warner2021}. Finally, a Dirac-delta wave-function is used to demonstrate the stability of the 1D HV ELBM model, while recovering a broadband transmittance spectrum for the same evaporated silver slab of 100 \si{nm} thickness within a single simulation, which reduces the computational time by 2-3 orders of magnitude compared with the continuous-wave simulation.  Since electric and magnetic fields are synchronous in the 1D HV ELBM, we find it requires 40\% less computational time and only 10\% more memory for obtaining similar numerical solutions as the 1D FDTD method \cite{KaneYee1966,Hauser2019}. Therefore, including an ADE in the ELBM is a stable and accurate alternative for modeling broadband frequency dispersion of electromagnetic waves at optical frequencies.

\section{Theory}
This section describes the auxiliary differential equation (ADE) with added complex-conjugate pole-residue pairs (CCPRPs) for modeling dispersion of unified Debye and Lorentz media.

\subsection{Auxiliary Differential Equation}
The ADE\cite{Okoniewski1997} can be implemented using CCPRPs to model dispersion of unified Debye and Lorentz media \cite{MinghuiHan2006}.  The ADE with CCPRPs uses a dispersive permittivity,
\begin{equation}\tag{1}\label{eq:1}
    \epsilon(\omega)=\epsilon_0\epsilon_\infty + \epsilon_0\sum_{p=1}^{N_p} \left(\frac{c_p}{j\omega-a_p}+\frac{c_p^*}{j\omega-a_p^*}\right),
\end{equation}
where $\omega$ is the angular frequency, $a_p$ and $c_p$ are the complex poles and residues, while $c_p^*$ and  $a_p^*$ are their conjugates \cite{MinghuiHan2006}. The dispersive permittivity $\epsilon(\omega\rightarrow\infty)=\epsilon_0\epsilon_\infty$, where $\epsilon_\infty$ is the relative permittivity of the medium as $\omega\rightarrow\infty$. The corresponding representations of $c_p$ and $a_p$ in terms of the pole relaxation time, $\tau_p$, plasmon frequency, $\omega_p$, and damping constant, $\delta_p$, are listed in Tab. \ref{tab:1}.
\begin{table}[!t]
    \centering
   \caption{Representation of pole-pairs for Debye and Lorentz media according to Eq. (\ref{eq:2})}
    \begin{tabularx}{\textwidth}{| X | X |}
    \hline
    Debye medium parameters & Lorentz medium parameters\\
    \hline
       $c_p=\Delta \epsilon_p/(2\tau_p)$  &  $c_p=j\Delta\epsilon_p\omega_p^2/(2\sqrt{\omega_p^2-\delta_p^2})$ \\
        $a_p=-1/\tau_p$ &  $a_p=-\delta_p - j\sqrt{\omega_p^2-\delta_p^2}$ \\
        \hline
    \end{tabularx}
    \label{tab:1}
\end{table}
Treating Debye and Lorentz media in a unified manner generalizes the implementation of dispersion \cite{MinghuiHan2006}. In the frequency domain, the corresponding polarization current density for each CCPRP is defined as
\begin{equation}\tag{2}\label{eq:2}
    \begin{split}
    \mathbf{J}_p(\omega) = \epsilon_0\frac{c_p}{j\omega-a_p}j\omega \mathbf{E}(\omega), \\
        \mathbf{J}_p^*(\omega) = \epsilon_0\frac{c_p^*}{j\omega-a_p^*}j\omega \mathbf{E}(\omega),
    \end{split}
\end{equation}
where $\mathbf{E}(\omega)$ is the electric field intensity. In the time-domain, only real polarization current density is required, so the conjugate terms are redundant \cite{MinghuiHan2006}.
Therefore, the auxiliary differential equation (ADE) for each CCPRP's current density in the time-domain becomes
\begin{equation}\tag{3}\label{eq:3}
    \frac{d }{dt} \mathbf{J}_p(\mathbf{x},t) - a_p\mathbf{J}_p(\mathbf{x},t) - \epsilon_0c_p\frac{d}{dt}\mathbf{E}(\mathbf{x},t)=0.
\end{equation}
This ADE can be implemented for multiple CCPRPs in both the finite-difference time-domain (FDTD) method and the Hauser \& Verhey electrodynamic lattice-Boltzmann method (HV ELBM).
\section{Methodology}
We compare the HV ELBM to the FDTD method for modeling one-dimensional frequency dispersion of electromagnetic waves transmitted through a 100 \si{nm} silver (Ag) slab \cite{Hauser2019,MinghuiHan2006}. This section first discusses the implementation of a central-difference ADE to describe current density in the HV ELBM \cite{MinghuiHan2006,Okoniewski1997,Hauser2019}. Then this section describes the CCPRP model used to represent silver \cite{MinghuiHan2006}, and concludes with a description of the methodology used to compare the HV ELBM with the same CCPR model used in the FDTD method \cite{MinghuiHan2006,Okoniewski1997}.
\subsection{Electrodynamic lattice-Boltzmann method}
The HV ELBM is implemented on a one-dimensional three-velocity (D1Q3) lattice which can incorporate $N=4$ independent non-zero left-handed and right-handed orthogonal bases with velocity, $\mathbf{c}_n$, electric field, $\mathbf{e}_n$, and magnetic field, $\mathbf{h}_n$, lattice vectors that satisfy a Chapman-Enskog analysis \cite{Hauser2019}.  In one dimension there is no orthogonal propagation and the sets of lattice vectors reduce to scalar sequences constrained in their respective dimensions. The sequences for the D1Q3 lattice used in this work are
\begin{align}
     &  \mathbf{c}_{n} = c_{n}\hat{\mathbf{x}} = \{1 \quad \mbox{-}1 \quad \mbox{-}1 \quad 1\}\hat{\mathbf{x}}, \tag{4a}\label{eq:4a} \\
     & \mathbf{e}_n =  e_{n}\hat{\mathbf{y}} = \{1 \quad 1 \quad \mbox{-}1 \quad \mbox{-}1\}\hat{\mathbf{y}}, \:\:\textrm{and} \tag{4b}\label{eq:4b} \\
     &  \mathbf{h}_n = h_{n}\hat{\mathbf{z}} = \{1 \quad \mbox{-}1 \quad 1 \quad \mbox{-}1\}\hat{\mathbf{z}}. \tag{4c}\label{eq:4c}
\end{align}
Given that there is no orthogonal propagation in one dimension, we find that the electric and magnetic fields propagate one lattice unit cell per each iteration using the sequences (\ref{eq:4a}-\ref{eq:4c}), such that the lattice unit is $\delta x=1$ and the time unit is $\delta t=1$, to represent the lattice velocity $c_0=1$. With a greater number of dimensions, orthogonal outcomes are required; therefore $\delta t=1/2$ for $\delta x=1$ and $c_0=1$ \cite{Hauser2019}.
 In one-dimension, we find that the 1/2 scaling factor for Hauser \& Verhey's Chapman-Enskog analysis for Maxwell's equations manifests exclusively within the initial condition while expected wave behaviour occurs at each consecutive time-step. The corresponding equilibrium pseudo-particle population is
\begin{equation} \tag{5}\label{eq:5}
    f_n^{eq}(\mathbf{x},t) = \frac{1}{4}(\mathbf{E}(\mathbf{x},t)\cdot\mathbf{e}_n+\mathbf{H}(\mathbf{x},t)\cdot\mathbf{h}_n),
\end{equation}
while the relaxation time is $\tau=1/2$. Therefore, the collision operation is
\begin{equation} \tag{6}\label{eq:6}
   f_n(\mathbf{x}+\mathbf{c}_n\delta t,t+\delta t) = 2f_n^{eq}(\mathbf{x},t)-f_n(\mathbf{x},t),
\end{equation}
where $f_n(\mathbf{x},t)$ is outcome $n$ of the pseudo-particle population.
\\
\subsubsection{Non-dispersive dielectric media}
The equilibrium polarization density, $\mathbf{P}^{eq}(\mathbf{x},t)$, is the outcome of the electric field remaining at its present location. This outcome is not included within $f^{eq}(\mathbf{x},t)$, but is crucial for electric and magnetic field interactions with dielectric media \cite{Hauser2017}. In a uniform non-dispersive medium, $\mathbf{P}^{eq}(\mathbf{x},t)$ can be defined as \cite{Hauser2017,Hauser2019}
\begin{equation}\tag{7}\label{eq:7}
    \mathbf{P}^{eq}(\mathbf{x},t)=(\epsilon_r(\mathbf{x})-1)\mathbf{E}(\mathbf{x},t),
\end{equation}
while the dynamic polarization density, $\mathbf{P}(\mathbf{x},t)$, is calculated with a collision operation,
\begin{equation}\tag{8}\label{eq:8}
    \mathbf{P}(\mathbf{x},t+\delta t)=2\mathbf{P}^{eq}(\mathbf{x},t)-\mathbf{P}(\mathbf{x},t).
\end{equation}
The polarization density is included within a re-normalized electric field, at the same location and time instant as the magnetic field. In this way, the electric and magnetic fields are calculated as synchronous orthogonal moments of $f(\mathbf{x},t)$,
\begin{equation}\tag{9}\label{eq:9}
\begin{split}
    & \mathbf{E}(\mathbf{x},t) = \frac{1}{\epsilon_r(\mathbf{x})}\left(\sum_{n=1}^Nf_n(\mathbf{x},t)\mathbf{e}_n+\mathbf{P}(\mathbf{x},t)\right) \\
   & \mathbf{H}(\mathbf{x},t) = \sum_{n=1}^Nf_n(\mathbf{x},t)\mathbf{h}_n.
    \end{split}
\end{equation}
\subsubsection{Dispersive dielectric media}
To model dispersive media using the HV ELBM, we include loss in the equilibrium condition for the non-dispersive polarization density, $\mathbf{P}^{eq}(\mathbf{x},t)$, in (\ref{eq:7}) via the Bhatnaghar-Gross-Krook (BGK) collision operator \cite{Bhatnagar1954}. This approach is analogous to the temporal implicit central difference scheme of the ADE used in the FDTD method \cite{Okoniewski1997,MinghuiHan2006}, and the Shan-Chen forcing scheme used in the kinetic lattice-Boltzmann method \cite{Shan1993}. It is also similar to the finite-element scheme used by Hauser \& Verhey for a uniform current density \cite{Hauser2019}. First the CCPRP parameters are converted to dimensionless forms as $\alpha_p=a_p(\Delta t e/\hbar)$, and $\chi_p = c_p(\Delta t e/\hbar)$, using  $\Delta t$ as the time-step, $e$ as the elementary unit of charge, and $\hbar$ as the reduced Planck constant for conversion from angular- to time-frequency \cite{MinghuiHan2006}. We interpret each CCPRP polarization current density in the ADE of  (\ref{eq:3}) as a population of the collision operator for the polarization density \cite{Hauser2021,hauser2022simulation}. Written in dimensionless form,
\begin{equation}\tag{10}\label{eq:10}
\begin{split}
    \mathbf{\Omega}_p(\mathbf{P}_p) & = -\mathbf{j}_p(\mathbf{x},t)=\frac{\chi_p}{\alpha_p}\frac{d}{dt} \mathbf{E}(\mathbf{x},t)-\frac{1}{\alpha_p}\frac{d}{dt} \mathbf{j}_p(\mathbf{x},t), \\
    & = -\frac{d}{dt}\mathbf{P}_p(\mathbf{x},t)=-\frac{\delta t}{\tau}(\mathbf{P}_p(\mathbf{x},t)-\mathbf{P}_p^{eq}(\mathbf{x},t)),
     \end{split}
\end{equation}
where $\mathbf{\Omega}_p(\mathbf{P}_p)$ is a population of the collision operator, which we propose takes a BGK form consistent with non-dispersive media \cite{Hanasoge2011,Hauser2017,Hauser2019}, as a function of the dimensionless polarization density for each CCPRP, $\mathbf{P}_p(\mathbf{x},t)$. We use a conventional backward difference in time to estimate each dimensionless CCPRP polarization current density in the present iteration
\begin{equation}\tag{11}\label{eq:11}
    \mathbf{j}_p( \mathbf{x},t)= \kappa_p  \mathbf{j}_p(\mathbf{x},t-\delta t)+\frac{\beta_p}{\delta t} (\mathbf{E}(\mathbf{x},t)-\mathbf{E}(\mathbf{x},t-\delta t)),
\end{equation}
where $k_p$ and $\beta_p$ are defined as the dimensionless parameters
\begin{equation}\tag{12}\label{eq:12}
    \kappa_p = \frac{1+\alpha_p\tau}{1-\alpha_p\tau}\qquad \beta_p = \frac{\chi_p}{1-\alpha_p\tau},
\end{equation}
which are identical to their formulation in the FDTD method \cite{MinghuiHan2006}. In higher dimensions, additional considerations must be made to preserve Gauss's law \cite{Hauser2019}, and more pole-pairs to preserve the model's stability \cite{Hauser2021,hauser2022simulation}.  Sixteen-pole Debye and Lorentz models for dispersive media have been described elsewhere \cite{Hauser2021,hauser2022simulation}.
Now a dispersive equilibrium polarization density, $\mathbf{j}^{eq}(\mathbf{x},t)$, is constructed as the sum of $p=\{1,...,N_p\}$ finite-difference approximations for each $\mathbf{j}_p(\mathbf{x},t)$ in a similar manner to its evaluation in the FDTD method \cite{MinghuiHan2006}. However, we re-interpret their sum as a moment of the BGK CCPRP collision operator, $\mathbf{\Omega}_p(\mathbf{P}_p)$, from (\ref{eq:10})
 \begin{equation}\tag{13}\label{eq:13}
     \begin{split}
\mathbf{j}^{eq}(\mathbf{x},t) & =\textrm{Re}\left\{\sum_{p=1}^{N_p}\mathbf{j}_p(\mathbf{x},t)(1+\kappa_p)\right\}, \\
& = \textrm{Re}\left\{\sum_{p=1}^{N_p}\mathbf{\Omega}_p(\mathbf{P}_p)(1+\kappa_p)\right\}, \\
& = \frac{\delta t}{\tau}(\mathbf{P}(\mathbf{x},t)-\mathbf{P}^{eq}(\mathbf{x},t)).
 \end{split}
 \end{equation}
 In an earlier work by Hauser \& Verhey \cite{Hauser2017}, the update equation for the polarization density was described using the BGK collision operator
 \begin{equation}\tag{14}\label{eq:14}
    \begin{split}
     \mathbf{P}(\mathbf{x},t+\delta t) & =\mathbf{P}(\mathbf{x},t)-\frac{\delta t}{\tau}(\mathbf{P}(\mathbf{x},t)-\mathbf{P}^{eq}(\mathbf{x},t)), \\
     & = \mathbf{P}(\mathbf{x},t)-\mathbf{j}^{eq}(\mathbf{x},t).
     \end{split}
 \end{equation}
 Therefore, it must also be true that
 \begin{equation}\tag{15}\label{eq:15}
    \mathbf{j}^{eq}(\mathbf{x},t)= -\left(\mathbf{P}(\mathbf{x},t+\delta t) -\mathbf{P}(\mathbf{x},t)\right).
 \end{equation}
 We can rearrange the collision operator from (\ref{eq:14}) to reveal the second-order equilibrium polarization density in terms of the equilibrium current density
 \begin{equation}\tag{16}\label{eq:16}
 \begin{split}
     \mathbf{P}^{eq}(\mathbf{x},t) & =\mathbf{P}(\mathbf{x},t)+\frac{\tau}{\delta t}(\mathbf{P}(\mathbf{x},t+\delta t)-\mathbf{P}(\mathbf{x},t)) \\
     & =\mathbf{P}(\mathbf{x},t)-\frac{\tau}{\delta t}\mathbf{j}^{eq}(\mathbf{x},t).
     \end{split}
 \end{equation}
This new equilibrium population for the polarization density, $\mathbf{P}^{eq}(\mathbf{x},t)$, can be adopted without modifying the original HV ELBM formulation \cite{Hauser2019}. However, to enforce its stability we further subject the dynamic polarization density, $\mathbf{P}(\mathbf{x},t)$, to the constraint that it must relax to a state of non-dispersive equilibrium from (\ref{eq:7}), i.e. $\mathbf{P}(\mathbf{x},t)\rightarrow\mathbf{P}^{eq}(\mathbf{x},t)=(\epsilon_r(\mathbf{x})-1)\mathbf{E}(\mathbf{x},t)$, within one time unit.
 Given that the relaxation time is constrained to $\tau=1/2$ for elastic scattering \cite{Hanasoge2011}, we simply subtract the effective equilibrium polarization current density from (\ref{eq:7}) to represent dispersion of the equilibrium polarization density,
\begin{equation}\tag{17}\label{eq:17}
    \mathbf{P}^{eq}(\mathbf{x},t)=(\epsilon_r(\mathbf{x})-1)\mathbf{E}(\mathbf{x},t)-\frac{\tau}{\delta t}\mathbf{j}^{eq}(\mathbf{x},t),
\end{equation}
where
\begin{equation}\tag{18}\label{eq:18}
 \epsilon_r = \epsilon_\infty + \sum_{p=1}^{N_p} \textrm{Re} \{ \beta_p \},
\end{equation}
represents the non-dispersive relative permittivity of the CCPRPs at equilibrium \cite{MinghuiHan2006}. Following this calculation of the equilibrium polarization density, $\mathbf{P}(\mathbf{x},t)$ can be evaluated according to (\ref{eq:8}) and $\mathbf{E}(\mathbf{x},t)$ and $\mathbf{H}(\mathbf{x},t)$ according to (\ref{eq:9}) without modifying the HV ELBM formulation \cite{Hauser2019}. Further work is required to show that this ADE approach agrees with analytical solutions in higher dimensions. This work employs a one dimensional analysis to demonstrate the model's unconditional stability and recovery of broadband transmittance, which is useful for quick frequency domain simulations.

\subsection{The discrete Dirac-delta wave-function}
The 1D HV ELBM is unconditionally stable, and may therefore represent the discrete Dirac-delta wave-function as its initial condition. The discrete Dirac-delta wave-function has a uniform frequency spectrum up to the spatiotemporal Nyquist limit in one dimension, making it the ideal source condition for modeling broadband electromagnetic frequency dispersion. Its initial condition is an electric and magnetic field represented at a single location in space, $\mathbf{x}_0$, for a single instant in time, $t_0$,
\begin{equation}\tag{19}\label{eq:19}
    \mathbf{E}(\mathbf{x}_0,t_0)=\frac{1}{2\delta x}\mathbf{e}_n, \qquad
    \mathbf{H}(\mathbf{x}_0,t_0) = \frac{1}{2\eta \delta x}\mathbf{h}_n.
\end{equation}
The $1/2$ factor manifests exclusively within the initial condition as a scaling factor of Hauser \& Verhey's Chapman-Enskog analysis to Maxwell's equations \cite{Hauser2019}. This is expected, given the duplicate elements in $\mathbf{e}_n$ and $\mathbf{h}_n$ that are required to describe all configurations of the electric field, the magnetic field and velocity on the lattice. Initializing one outcome with a half-amplitude is necessary to represent a unit amplitude, which is maintained after the first iteration in the 1D HV ELBM. Alternatively, the 1D FDTD method is conditionally stable up to the Courant-Friedrichs-Lewy (CFL) limit \cite{Courant1928,KaneYee1966,Sarris2011}, such that it is not possible for the 1D FDTD method to represent the discrete Dirac-delta wave-function and its frequency components at the Nyquist limit. The 1D FDTD method will instead become unstable and exhibit unbounded growth in numerical error.
\subsection{Power spectral density}
The power spectral density is recovered from the Poynting vector of electromagnetic sine waves fully transmitted through the 100 \si{nm} silver slab at multiple frequencies. We believe this is the most intuitive way to evaluate and compare the relative error of the simulated dispersion in both the HV ELBM and the FDTD method. One advantage of the 1D HV ELBM over the 1D FDTD method is that the 1D HV ELBM can represent a Dirac-delta wave-function with a uniform power spectral density approaching the Nyquist limit. Given this input, we demonstrate the accurate recovery of the 100 \si{nm} silver slab's transmittance as the multiplication of fast-Fourier transforms (FFTs) of electric and magnetic fields fully transmitted through the slab. Their multiplication in the frequency domain is equivalent to their convolution in the time-domain,
\begin{equation}\tag{20}\label{eq:20}
    \mathbf{\hat{S}}(\nu)=\mathbf{\hat{E}}(\nu)\times\mathbf{\hat{H}}(\nu),
\end{equation}
where $\nu=\Delta t/T'$ represents the dimensionless time-frequency for a dimensional (physical) time unit $\Delta t$ and time period $T'$. The power spectral density of the Dirac-delta wave-function's fully transmitted electric and magnetic fields approximates the transmittance spectrum of the 100 \si{nm} silver slab as an impulse response function with an $l$\textsuperscript{th} order-of-accuracy, $\mathcal{O}(\Delta \nu^l)$,
\begin{equation}\tag{21}\label{eq:21}
    T_{slab}(\nu)=|\mathbf{\hat{S}}(\nu)|+\mathcal{O}(\Delta \nu^l).
\end{equation}
In this manner, one simulation can be used to retrieve a broadband transmittance spectrum representing the frequency-dispersion of electromagnetic waves through a dispersive medium, as opposed to the inherently band-limited spectra that represent all continuous wave-functions. The error order-of-accuracy is quantified using the first-order, $L_1$, and second-order, $L_2$, relative error norms as a function of the number of photon energy bins, $N_\nu$, across the spectrum. The first-order and second-order norms are defined as
\begin{equation}\label{eq:22}\tag{22}
\begin{aligned}
   & L_1=\frac{1}{N_\nu}\sum_{\nu}^{N_\nu}\left|\frac{T_{slab}(\nu)-|\hat{\mathbf{S}}(\nu)|}{T_{slab}(\nu)}\right|, \\
   & L_2=\left(\frac{\sum_{\nu}^{N_\nu}|T_{slab}(\nu)-|\hat{\mathbf{S}}(\nu)||^2}{\sum_{\nu}^{N_\nu}|T_{slab}(\nu)|^2}\right)^{1/2}.
\end{aligned}
\end{equation}
\subsection{Material Model}
We choose a CCPRP model which represents dispersive silver (Ag) using six CCPRPs \cite{MinghuiHan2006}. We compare this model to measurements of the refractive index, $n$, and extinction coefficient, $k$, for evaporated silver \cite{Lynch1985}, and a more recent fit obtained using four CCPRPs \cite{Gharbi2020}. We find that six CCPRPs, including $a_p$ and $b_p$, can more accurately represent the measurement data for evaporated silver \cite{MinghuiHan2006,Lynch1985}. Tab. \ref{tab:2} lists these parameters, and Fig. \ref{fig:1} demonstrates their fit to empirical data. Nonetheless, we find that this CCPRP model has insufficient pole-pairs to accurately represent the available bandwidth of the 1D HV ELBM, so we only compare its error with the FDTD method in the photon energy range of interest, 1-5 \si{eV}.
\begin{figure}[!t]
    \centering
    \includegraphics[width=\textwidth]{fig1}
    \caption{Comparison of Han \emph{et al.} \cite{MinghuiHan2006}, and Gharbi \emph{et al.} \cite{Gharbi2020} fitting to measured index of refraction, $n$, and extinction coefficient, $k$, for evaporated silver (Ag) \cite{Lynch1985}. Han \emph{et al.}'s model provides a better fit.}
    \label{fig:1}
\end{figure}
\begin{table}[!t]
    \centering
        \caption{Values of the Complex-Conjugate Pole-Residue Pairs for the permittivity of Ag\cite{MinghuiHan2006}}
    \begin{tabularx}{\textwidth}{|c X|c X|}
    \hline
    $c_p$, & \si{eV} & $a_p$, & \si{eV} \\
    \hline
       \scriptsize $c_1=$ & \scriptsize $5.987\cdot10^{\mbox{-}1}+\textrm{j}4.195\cdot 10^3$ & \scriptsize $a_1=$ & \scriptsize $\mbox{-}2.502\cdot10^{\mbox{-}2} - \textrm{j}8.626\cdot10^{\mbox{-}3}$   \\
        \scriptsize$c_2=$ & \scriptsize$\mbox{-}2.211\cdot10^{\mbox{-}1}+\textrm{j}2.680\cdot10^{\mbox{-}1}$ &\scriptsize $a_2=$ &\scriptsize $\mbox{-}2.021\cdot10^{\mbox{-}1} - \textrm{j}9.407\cdot10^{\mbox{-}1}$  \\
        \scriptsize$c_3=$ &\scriptsize  $\mbox{-}4.240+\textrm{j}7.324\cdot10^{\mbox{-}2}$ &\scriptsize $a_3=$ &\scriptsize $\mbox{-}1.467\cdot10^1 - \textrm{j}1.338$ \\
        \scriptsize$c_4=$ &\scriptsize $6.391\cdot10^{\mbox{-}1}-\textrm{j}7.186\cdot10^{\mbox{-}2}$ &\scriptsize $a_4=$ &\scriptsize $\mbox{-}2.997\cdot10^{\mbox{-}1}-\textrm{j}4.034$ \\
       \scriptsize $c_5=$ &\scriptsize $1.806 + \textrm{j}4.563$ &\scriptsize $a_5=$ &\scriptsize $\mbox{-}1.896-\textrm{j}4.808$ \\
       \scriptsize $c_6=$ &\scriptsize $1.443-\textrm{j}8.129\cdot10^1$ &\scriptsize $a_6=$ &\scriptsize $\mbox{-}9.396-\textrm{j}6.477$ \\
        \hline
    \end{tabularx}
    \label{tab:2}
\end{table}
\subsection{Simulation}\label{sec:sim}
 Three different source conditions are used to analyze the accuracy of electromagnetic wave frequency-dispersion from a 100 \si{nm} silver slab in the photon energy range of interest (1-5 \si{eV}): (1) a Ricker wavelet with a peak photon energy of $\hbar\omega_p/e=3.8735$ \si{eV} and a half-breadth of $\tau_b=145.32$ \si{as} (7.02 time units); (2) a sequence of time-harmonic sine-waves of variable period; and (3) a Dirac-delta wave-function.
 Source condition (1) is used to compare the numerically solved electric and magnetic fields in the 1D FDTD method and the 1D HV ELBM. Source condition (2) is used to compare their numerical accuracy and computational time. Source condition (3) is used to demonstrate the versatility of the 1D HV ELBM for modeling single-simulation broadband frequency dispersion of electromagnetic waves. \\ \\
 Source condition (1) is situated at node $x_0=125$ on a lattice consisting of $N_x=500$ spatial unit cells. The silver slab is situated at the center of the lattice, and its thickness is 99.2 \si{nm}, to represent an integer number of dimensional spatial units each of 1/40th the smallest wavelength of interest ($\Delta x=6.2$ \si{nm}). This 3100 \si{nm} spatial domain is identical to that of a former FDTD simulation \cite{MinghuiHan2006}, to ensure a faithful reproduction of the original numerical solution. \\ \\
 Source condition (2) is situated at position $x_0=N_{x_k}/4$ within a 620 \si{nm}  spatial domain. The spatial domain's dimensional length is kept constant, but the number of its lattice units is linearly increased for each consecutive simulation indexed by $k$. The $k^{\textrm{th}}$ simulation has $N_{x_k}=100k$ spatial lattice units and represents $N_{\nu_k}=16k$ frequencies and $16k$ points of the transmittance spectrum. Each simulation achieves a steady state and the time-averaged Poynting vector is measured to estimate the numerical transmittance. As the number of spatial units increases, the silver slab thickness approaches 100 \si{nm}. The silver slab's exact thickness (from the number of spatial lattice units) is used in the analytical solution to assess the relative error norms. The 1D HV ELBM and 1D FDTD method relative error norms are plotted with respect to $N_\nu$, and given a similar accuracy, their computational time and memory are compared via linear regression.\\ \\
Source condition (3) is first initialized at dimensionless time $t_0=0$ and position $x_0=N_{x}/4$ in a 3100 \si{nm} spatial domain to illustrate the stable transmission of electromagnetic waves through the 100 \si{nm} slab in the 1D HV ELBM with the added ADE. Source condition (3) is then initialized at dimensionless time $t_0=0$ and position $x_0=N_{x_k}/4$ in a 620 \si{nm} spatial domain used for a sequence of $k$ simulations with $N_{x_k}=100k$ lattice units and $N_{t_k}=4000k$ time units. The transmitted power spectral density from source condition (3) is compared with the analytical transmittance spectrum of the 100 \si{nm} silver slab. The power spectral density's relative error (relative difference from the analytical transmittance spectrum), for source (3) is plotted with respect to the number of unit cells to determine the order-of-accuracy of the 1D HV ELBM solutions from the slope of a proportional error norm.

\subsection{Boundary conditions}
Absorbing boundary conditions are used on each end of the one dimensional spatial domain described in subsection \ref{sec:sim}, but with a distinct implementation in the FDTD and the HV ELBM. We choose the Mur boundary condition for the FDTD simulation, since it perfectly absorbs electric and magnetic fields in one dimension \cite{mur1981absorbing}. In the HV ELBM simulation we implement similar absorbing boundary conditions, called ``free'' boundaries, which assign the  boundary particle populations to their nearest adjacent nodes inside the domain,
\begin{equation}\tag{23}\label{eq:23}
    f_n(\mathbf{x}_b-\delta x\hat{\mathbf{n}}_b,t)=f_n(\mathbf{x}_b,t),
\end{equation}
where $\hat{\mathbf{n}}_b$ is the unit vector normal to boundary $b$ at position $\mathbf{x}_b$.
\newpage

\section{Results and Discussion}
Fig. \ref{fig:2} illustrates a Ricker wavelet with a peak photon energy of $\hbar\omega_p/e=3.8735$ \si{eV} and a half-breadth of $\tau_b=145.32$ \si{as} (7.02 time units) as it transmits through a 100 \si{nm} thick silver slab in both the 1D FDTD method and the 1D HV ELBM. The 1D HV ELBM and 1D FDTD method predict similar incident, reflected, and transmitted electric and magnetic fields. However, there is a slight difference in phase between the magnetic field evaluated in the 1D HV ELBM and 1D FDTD method. This slight difference in phase may be attributed to the staggered Yee-lattice in the FDTD method, which will evaluate the magnetic field at a different time and location in space than the fully synchronous 1D HV ELBM.

\begin{figure} [H]
    \centering
    \includegraphics[width=\textwidth]{fig2}
    \caption{[Color online]: The electric field (top) and magnetic field (bottom) for the transmission of a Ricker wavelet with a peak photon energy of $\hbar\omega_p/e=3.8735$ \si{eV} and a half-breadth of $\tau_b=145.32$ \si{as} (7.02 time units) through a silver slab of 100 \si{nm} thickness (boundaries represented by black-dashed lines). Both the 1D FDTD method and the 1D HV ELBM recover similar electric and magnetic fields, apart from a phase delay in the magnetic field as a consequence of the Yee-lattice. Axis units are in nanometers, \si{nm}.}
    \label{fig:2}
\end{figure}
Fig. \ref{fig:3} compares the relative errors, while Fig. \ref{fig:4} compares the computational time and memory, of the 1D HV ELBM and 1D FDTD method for source condition (2). While the two methods predict identical error norms from 1-5 \si{eV}, the 1D HV ELBM requires 60\% of the computational time for only 10\% more memory. Producing a propagating Dirac-delta wave-function in the 1D HV ELBM allows it to recover the broadband transmittance spectrum of a dispersive medium within a single simulation, which reduces the computational time by 2-3 orders of magnitude compared with the continuous-wave simulation.
\begin{figure}[!t]
    \centering
    \subfloat[]{\includegraphics[width=0.49\textwidth]{Transmittance_source_2}\label{fig:3a}}
     \subfloat[]{\includegraphics[width=0.49\textwidth]{Error_norms_source_2}\label{fig:3b}}
    \caption{[Color online]: (a) The transmittance spectrum from 160 continuous waves [source condition (2)] in the photon energy range 1-5 \si{eV} transmitted through a 100 \si{nm} silver slab. (b) The error norms ($L_1$) and ($L_2$) from the transmittance spectrum measured for an increasing number of spatial lattice units ($N_x$) and frequencies ($N_{\nu}$) between 1-5 \si{eV}. Red symbols are beneath the black symbols, both solutions have the same numerical error.}
    \label{fig:3}
\end{figure}
\begin{figure}[!t]
    \centering
    \subfloat[]{\includegraphics[width=0.49\textwidth]{Computational_time_source_2}\label{fig:4a}}
     \subfloat[]{\includegraphics[width=0.49\textwidth]{Memory_source_2}\label{fig:4b}}
    \caption{[Color online]: (a) The computational time and (b) memory for the simulation of multiple continuous waves [source condition (2)]. See \ed{Supplementary material} for MATLAB scripts. The 1D HV ELBM requires 60\% of the computation time, but 10\% greater memory to achieve the same numerical solutions from 1-5 \si{eV}. Simulations were performed on an 11\textsuperscript{th} Gen Intel\textsuperscript{\textregistered} Core\textsuperscript{TM} i7-11700K @ 3.60 GHz processor with 32 GB of RAM and 8 cores. }
    \label{fig:4}
\end{figure}
Fig. \ref{fig:5} demonstrates the transmission of a Dirac-delta wave-function through the 100 \si{nm} thick silver slab in the 1D HV ELBM. The numerical stability of the simulation is maintained as the Dirac-delta wave-function transmits through the silver slab to represent its entire transmittance spectrum within the photon energy range of interest (1-5 \si{eV}).
\newpage
\begin{figure} [t!]
    \centering
    \includegraphics[width=\textwidth]{Dirac_delta_fields.pdf}
    \caption{Transmission of electric and magnetic fields from a Dirac-Delta wave-function [source condition (3)] as it interacts with a silver slab of 100 \si{nm} thickness (boundaries represented by black-dashed lines). Axis units are in nanometers, \si{nm}.}
    \label{fig:5}
\end{figure}
\newpage
Fig. \ed{\ref{fig:6b}} shows a plot of the transmittance's relative error, which has a mean of 0.28\%  from 1-5 \si{eV}. Since the simulation was performed using $N_x=1000$ lattice units and $N_t=40000$ time units, the spectrum has a much greater bandwidth (0-1000 \si{eV}). While the 1D HV ELBM may represent frequencies up to the Nyquist limit, this requires developing a model with a greater number of CCPRPs representing the transmittance of silver over a broader spectrum. Nonetheless, a model fitted with six CCPRPs provides a reasonably accurate solution in the photon energy range of interest (1-5 \si{eV}), while remaining stable up to the Nyquist limit. Stability of the ADE allows calculating the relative error norms from the entire spectrum, for an accurate evaluation of the 1D HV ELBM order of accuracy.
\begin{figure}[!t]
    \centering
    \subfloat[]{\includegraphics[width=0.49\textwidth]{Transmittance_Dirac}\label{fig:6a}}
    \subfloat[]{\includegraphics[width=0.49\textwidth]{Transmittance_Dirac_error}\label{fig:6b}}
    \caption{(a)  The 1D HV ELBM Dirac-delta wave-function's transmitted power spectral density compared with the analytical transmittance for a 100 \si{nm} thickness silver slab. The two are equivalent in this case. (b) The relative error in transmittance, with a mean of 0.28\% in the range 1-5 \si{eV} for $N_x=1000$ spatial lattice units, after $N_t=40000$ iterations for a computational time of 7.97 s. Simulation was performed on an 11\textsuperscript{th} Gen Intel\textsuperscript{\textregistered} Core\textsuperscript{TM} i7-11700K @ 3.60 GHz processor with 32 GB of RAM and 8 cores.}
    \label{fig:6}
\end{figure}
The transmittance spectrum, evaluated from (\ref{eq:20}) for the case of the Dirac-delta wave-function transmitted through the 100 \si{nm} silver slab, is used to evaluate the 1D HV ELBM's order of accuracy across its entire bandwidth. The number of unit cells, $N_{x_k}$, is varied from $N_{x_1}=100$ ($\Delta x=6.2$ \si{nm}) to $N_{x_{50}}=5000$ ($\Delta x=0.124$ \si{nm}), while the number of time units, $N_{t_k}$, is varied from $N_{t_1}=4000$ ($\Delta t=20.7$ \si{as}), to $N_{t_{50}}=200000$ ($\Delta t=0.414$ \si{as}) to measure the first- and second-order relative error norms, $L_1$ and $L_2$, each with respect to both the number of lattice units, $N_x$, and the number of photon energy bins (frequency bins), $N_\nu$. The first- and second-order error norms have slopes that are parallel to the second-order relative error, as demonstrated in Fig. \ref{fig:7}. Therefore, the 1D HV ELBM in this study obtains a second-order accuracy and corroborates an earlier study in which the 1D HV ELBM's energy density achieved an absolute error of second-order accuracy for non-dispersive media \cite{Hauser2019}. Unlike absolute error \cite{Chen2013}, relative error is independent of the transmittance amplitude and is a preferable measure for the order-of-accuracy in this application. The ADE CCPRP model utilized by Han \emph{et al.} is applicable to photon energies ranging from 0.125-5 \si{eV} \cite{MinghuiHan2006}. Analyzing this model's transmittance over the 1D HV ELBM's complete spectrum (ranging from 0-100 \si{eV} at $N_{x_1}=100$ and $N_{\nu_1}=1986$, to 0-5000 \si{eV} at $N_{x_{50}}=5000$ and $N_{\nu_{50}}=99176$), we find its analytical transmittance amplitude increases significantly at photon energies exceeding 25 \si{eV}.  The relative error in Fig. \ref{fig:6} accounts for this variation in transmittance amplitude, and continues to decrease with a second-order slope. The error is significant for lower spatial discretizations, demonstrating that even a six-pole model is ill-suited for the broad spectrum represented by the 1D HV ELBM. Nonetheless, the relative error decreases to less than 1\% for $N_x>1100$, and demonstrates a second order-of-accuracy for the decreasing dimensional unit cell and time unit. Second order-of-accuracy is expected for the 1D HV ELBM model \cite{Hauser2019}, and the central difference equation used to incorporate the ADE \cite{MinghuiHan2006}. The second-order accurate 1D HV ELBM is an improvement from an earlier first-order accurate 1D ELBM for dispersive media \cite{Chen2013}.
\begin{figure}[!t]
    \centering
    \subfloat[]{\includegraphics[width=0.49\textwidth]{Error_norms_DD_Nx.pdf}}
    \subfloat[]{\includegraphics[width=0.49\textwidth]{Error_norms_DD_Nnu.pdf}}
    \caption{ The relative error is plotted with a parallel slope to a second order of accuracy with respect to both (a) number of lattice units and (b) number of frequency bins. This CCPRP model requires $N_x>1100$ in order to obtain solutions within 1\% accuracy across the entire available spectrum (0-100 \si{eV} at $N_{x_1}=100$, 0-5000 \si{eV} at $N_{x_{50}}=5000$). This demonstrates that even a six pole CCPRP model is ill-suited for the HV ELBM's broad spectrum. The relative error is a more reliable measure for order-of-accuracy, given the large variation in the transmittance spectrum of silver. }
    \label{fig:7}
\end{figure}
\section{Conclusion}
The auxiliary-differential-equation with complex-conjugate pole-residue pairs can be added to the one-dimensional electrodynamic Hauser \& Verhey electrodynamic lattice-Boltzmann method (1D HV ELBM) to model the electromagnetic frequency dispersion of an evaporated silver slab. The implementation treats Debye and Lorentz media in a unified manner while recovering electric and magnetic fields similar to the one-dimensional finite-difference-time-domain (1D FDTD) method. Compared with the 1D FDTD method for the same lattice unit and accuracy, the 1D HV ELBM requires 60\% the computational time and only 10\% greater memory. However, this is not the main advantage of the 1D HV ELBM. Since the 1D HV ELBM predicts the stable transmission of a Dirac-delta wave-function through silver, it can describe its transmittance spectrum within a single simulation. Therefore, given a sufficiently stable complex-conjugate pole-pair model for a dispersive medium, we find the 1D HV ELBM is a powerful alternative to the 1D FDTD method for modeling the broadband frequency dispersion of electromagnetic waves.

\bibliographystyle{IEEEtran}

redpajama_set_name: RedPajamaArXiv