---
abstract: |
    The linearized operator for non-radial oscillations of spherically symmetric self-gravitating gaseous stars is analyzed in view of the functional analysis. The evolution of the star is supposed to be governed by the Euler-Poisson equations under the equation of state of the ideal gas, and the motion is supposed to be adiabatic. We consider the case of not necessarily isentropic, that is, not barotropic motions. Basic theory of self-adjoint realization of the linearized operator is established. Some problems in the investigation of the concrete properties of the spectrum of the linearized operator are proposed.\
    Key Words and Phrases. Gaseous star. Adiabatic Oscillation. Self-adjoint operator. Friedrichs extension. Spectrum of Sturm-Liouville type. Brunt-Va̋isa̋la̋ frequency. g-Mode.

    2010 Mathematical Subject Classification Numbers. 35P05, 35L51, 35Q31, 35Q85, 46N20, 76N15.
author:
- 'Tetu Makino [^1]'
title: 'On Linear Adiabatic Perturbations of Spherically Symmetric Gaseous Stars Governed by the Euler-Poisson Equations '
---

Introduction
============

We consider the adiabatic hydrodynamic evolution of a self-gravitating gaseous star governed by the Euler-Poisson equations

$$\begin{aligned}
&\frac{\partial\rho}{\partial t}+\sum_{k=1}^3\frac{\partial}{\partial x^k}(\rho v^k)=0, \label{EPa} \\
&\rho\Big(\frac{\partial v^j}{\partial t}+\sum_{k=1}^3 v^k\frac{\partial v^j}{\partial x^k}\Big)
+\frac{\partial P}{\partial x^j}+\rho\frac{\partial\Phi}{\partial x^j}=0, \quad j=1,2,3,
 \label{EPb} \\
&\rho\Big(\frac{\partial S}{\partial t}+\sum_{k=1}^3 v^k\frac{\partial S}{\partial x^k}\Big)=0, \label{EPc} \\
&\triangle \Phi =4\pi\mathsf{G}\rho. \label{EPd}\end{aligned}$$

Here $t\geq 0, \mathbf{x}=(x^1,x^2,x^3) \in \mathbb{R}^3$. The unknowns $\rho \geq 0, P, S, \Phi \in \mathbb{R}$ are the density, the pressure, the specific entropy, the gravitational potential, and $\mathbf{v}=(v^1,v^2,v^3) \in \mathbb{R}^3$ is the velocity fields. $\mathsf{G}$ is a positive constant, the gravitation constant.

In this article the pressure $P$ is supposed to be a prescribed function of $\rho, S$. But for the sake of simplicity, we assume the equation of state of the ideal fluid, that is, we assume

$P$ is the function of $(\rho, S) \in [0,+\infty[ \times \mathbb{R}$ given by $$P=\rho^{\gamma}\exp\Big(\frac{S}{\mathsf{C}_V}\Big), \label{DefP}$$ where $\gamma$ and $\mathsf{C}_V$ are positive constants such that $$1<\gamma <2.$$

The constant $\gamma$ is the adiabatic exponent and $\mathsf{C}_V$ is the specific heat per unit mass at constant volume.

Since we are concerned with compactly supported density distribution $\rho$ in this article, the Poisson equation will be replaced by the Newtonian potential $$\Phi(t,\mathbf{x})=-4\pi\mathsf{G}\mathcal{K}\rho(t,\cdot)(\mathbf{x}), \label{NPot}$$ where $$\mathcal{K}f(\mathbf{x}):=\frac{1}{4\pi}\int\frac{f(\mathbf{x}')}{|\mathbf{x}-\mathbf{x}'|}d\mathbf{x}'. \label{defK}$$

We suppose that there is fixed a spherically symmetric equilibrium $\bar{\rho}, \bar{P},
\bar{S}, \bar{\Phi}$, which satisfy , , , , such that $\bar{\rho}(\mathbf{x})>0 \Leftrightarrow r=|\mathbf{x}| <R$ with a finite positive number $R$, the radius of the equilibrium.

We consider the perturbation $\mbox{\boldmath$\xi$}=\delta\mathbf{x}, \delta\rho, \delta P, \delta S, \delta\Phi$ at this fixed equilibrium. We use the Lagrangian co-ordinate which will be dented by the diversion of the letter $\mathbf{x}$ of the Eulerian co-ordinate. So, $\mathbf{x}$ runs on the fixed domain $B_R:=\{ \mathbf{x} \in \mathbb{R}^3 | r=|\mathbf{x}| < R\}$, while $\{\rho >0\} $ described by the Eulerian co-ordinate may move along $t$.

Then the linearized equation which governs the perturbations turns out to be $$\frac{\partial^2\mbox{\boldmath$\xi$}}{\partial t^2}+\mathbf{L}\mbox{\boldmath$\xi$}=0,
\label{Eqxi}$$ where $$\begin{aligned}
\mathbf{L}\mbox{\boldmath$\xi$}&=\frac{1}{\bar{\rho}}\mathrm{grad}\delta P+\frac{\delta\rho}{\bar{\rho}}\mathrm{grad}\bar{\Phi}+
\mathrm{grad}\delta\Phi \nonumber \\
&=\frac{1}{\bar{\rho}}\mathrm{grad}\delta P
-\frac{\delta\rho}{\bar{\rho}^2}\mathrm{grad}\bar{P}+
\mathrm{grad}\delta\Phi. \label{L}\end{aligned}$$

We have $$\begin{aligned}
&\delta\rho=-\mathrm{div}(\bar{\rho}\mbox{\boldmath$\xi$}), \label{drho} \\
&\delta\Phi=-4\pi\mathsf{G}\mathcal{K}(\delta\rho). \label{dPhi}\end{aligned}$$

Here $\delta$ denotes the Eulerian perturbation, while $\Delta$ will denote the Lagrangian perturbation. Recall their definition $$\begin{aligned}
&\Delta Q(t,\mathbf{x})=Q(t, \mbox{\boldmath$\varphi$}(t,\mathbf{x}))-\bar{Q}(\mathbf{x}), \\
&\delta Q(t,\mathbf{x})=Q(t,\mbox{\boldmath$\varphi$}(t,\mathbf{x}))-\bar{Q}
(\mbox{\boldmath$\varphi$}(t,\mathbf{x})),\end{aligned}$$ where $\mbox{\boldmath$\varphi$}(t,\mathbf{x})=\mathbf{x}+\mbox{\boldmath$\xi$}(t,\mathbf{x})$ is the steam line given by $$\frac{\partial}{\partial t}\mbox{\boldmath$\varphi$}(t,\mathbf{x})=\mathbf{v}(t,\mbox{\boldmath$\varphi$}(t,\mathbf{x})),
\quad \mbox{\boldmath$\varphi$}(0,\mathbf{x})=\mathbf{x}.$$ Thus it holds that $$\Delta Q=\delta Q+(\mbox{\boldmath$\xi$}|\mathrm{grad}\bar{Q})$$ in the linearized approximation, for any quantity $Q$.

Supposing that the initial perturbation of the density vanishes, that is, $\Delta\rho|_{t=0}=0$, the equation implies $\Delta\rho=0$ always, which is . Supposing $\Delta S|_{t=0} =0$, the equation implies $\Delta S=0$ always, therefore $$\Delta P=\overline{\frac{\gamma P}{\rho}}\Delta\rho.$$ This implies $$\delta P=\overline{\frac{\gamma P}{\rho}}\delta\rho+
\gamma\mathscr{A}\bar{P}(\mbox{\boldmath$\xi$}|\mathbf{e}_r) \label{dP}.$$\
Here we denote by $\mathbf{e}_r$ the unit vector $\partial/\partial r$ and we define the ‘Schwarzschild’s discriminant of convective stability’ $\mathscr{A}$ with the ‘Brunt-Va̋isa̋la̋ frequency’ $\mathscr{N}$ by

We put $$\mathscr{A}:=\overline{\frac{1}{\rho}\frac{d\rho}{dr}}-\overline{\frac{1}{\gamma P}\frac{dP}{dr}}\quad \Big(=-\frac{1}{\gamma \mathsf{C}_V}\frac{d\bar{S}}{dr}\Big)$$ and $$\mathscr{N}^2:=\mathscr{A}\overline{\frac{1}{\rho}\frac{dP}{dr}}
=-\mathscr{A}\overline{\frac{d\Phi}{dr}}.$$

For the physical meaning of these quantities, see [@LedouxW] or [@Cox Chapter III, Section 17]. When $\overline{dP/dr}<0$, $\mathscr{N}$ is real if and only if $\mathscr{A} \leq 0$. The condition $\mathscr{A}<0$ is that of the convective stability.\
Thus by , , we can see the right-hand side of is an integro-differential operator acting on the unknown $\mbox{\boldmath$\xi$}$, provided that the spherically symmetric equilibrium $(\bar{\rho}, \bar{S}, \bar{P})$ is fixed.

For the derivation of $\mathbf{L}$, see e.g., [@LedouxW], [@Cox] or [@Unno].\
The purpose of this article is to clarify the functional analysis properties of this integro-differential operator $\mathbf{L}$.\
Nonlinear evolution of spherically symmetric perturbations has been investigated sufficiently well in [@TM.OJM] and [@JJ.APDE]. In these studies spectral properties of the linearized operator for spherically symmetric perturbations, which was established by [@Beyer1995] and independently by [@Lin], are fully presupposed. Its spectrum was proved to be actually of the Sturm-Liouville type, and it was not obvious because of the singularity of the coefficients, caused by the physical vacuum boundary of the equilibrium. Therefore if we want to study nonlinear evolution of not necessarily spherically symmetric perturbations around a spherically symmetric equilibrium, we should prepare a sufficiently strong functional analysis study of spectral properties of the linearized operator for general, not necessarily spherically symmetric, perturbations. As for barotropic case, we have attacked this task, and have gotten sufficiently strong results in [@JJTM]. Thus here we consider the case of not necessarily barotropic motions. Unfortunately the results which we have established is little bit weaker than the barotropic case. There remains some open problems. But mathematically rigorous treatment of the problem is quiet new.\
This article is organized as follows.

In Section 2, we discuss on the existence of spherically symmetric equilibrium for prescribed entropy distribution. The concept of the ‘admissible’ equilibrium will play a crucial rôle throughout the mathematically rigorous investigations of this article. In Section 3, we prove the self-adjoint realization of the operator $\mathbf{L}$ as the Friedrichs extension. Astrophysical texts lacked mathematically rigorous proof. But such a strong assertion on the concrete form of the spectrum as that of the barotropic case given in [@JJTM] is not yet obtained. In order to investigate the specified concrete form of the spectrum we investigate eigenfunctions represented by spherical harmonics in Section 4. The situation is clarified to be quite different from the barotropic case. That is, it may be impossible to reduce the problem to that of Sturm-Liouville type. But the justification of the self-adjoint realization of the associated operator $\vec{L}_l$ for each degree $l$ of the harmonics $Y_{lm}$ can be done with success. A strong guess that the form of the spectrum of $\vec{L}_l$ is quite different from that of the barotropic case is suggested by the so-called ‘g-modes’. Section 5 is devoted to the study of the ‘g-modes’ and ‘p-modes’ proposed by astrophysicists. Justification of the Sturm-Liouville properties of the operators associated to ‘g-modes’, ‘p-modes’ can be done with success. But it may be difficult to prove that the sequence of eigenvalues which accumulates both to 0 and to $+\infty$ give a good approximation of the real eigenvalues of the original operators $\vec{L}_l$. This is a task to be done in a future work. The last Section 6 is devoted to examination of the arguments in the work [@Eisenfeld] by J. Eisenfeld. Besides the lack of the explicit announcement of the assumption that $\displaystyle \frac{1}{\gamma-1}$ be a rational number, because of the serious lack of proof of the identification of the spectra and the eigenvalues, the argument of [@Eisenfeld] does not succeed to prove the completeness of eigenfunctions as claimed. Therefore there is an open problem here concerning the proof of the absence of continuous spectrum.\
We shall use the following notations:

We denote

$$\begin{aligned}
&B_R:=\{ \mathbf{x}\in\mathbb{R}^3 \  |\  r=|\mathbf{x}| <R \}, \\
& \overline{B_R}:=\{ \mathbf{x}\in\mathbb{R}^3 \  |\  r=|\mathbf{x}| \leq R \}.\end{aligned}$$

\[Not.1\] 1) A function $F$ on a subset of $\mathbb{R}^3$ is said to be spherically symmetric if there exists a function $f$ on a subset of $[0,+\infty[$ such that $F(\mathbf{x})=f(|\mathbf{x}|)$ for $\forall \mathbf{x}$ in the domain of $F$. Then we shall denote $f=F^{\sharp}$.

2\) For a function $f$ on a subset of $[0,+\infty[$, we shall denote by $f^{\flat}$ the function on a subset of $\mathbb{R}^3$ such that $f^{\flat}(\mathbf{x})=
f(|\mathbf{x}|)$ for $\forall \mathbf{x}$ such that $r=|\mathbf{x}|$ is in the domain of $f$.

3\) When it is expected that no confusion may occur, we shall divert the symbols $f$ or $F$ instead of $f^{\flat}$ or $F^{\sharp}$.

Here let us note the following lemma, which can be verified easily:

\[Lem.0\] If a function $f$ defined on $[0,R[$ satisfies $f^{\flat}\in C^{k+2}(B_R), k=0,1,2,\cdots$, then $f \in C^{k+2}([0,R[), \quad \displaystyle \frac{df}{dr}\Big|_{r=+0}=0$,   and $\displaystyle \frac{1}{r}\frac{df}{dr} \in C^k([0,R[)$.

Proof. We can show inductively that $$D^kh(r)=\frac{1}{r^{k+1}}\int_0^rD^{k+2}f(s)s^kds$$ for $\displaystyle h=\frac{1}{r}\frac{df}{dr}$. Therefore $$D^kh(r) \rightarrow \frac{D^{k+2}f(0)}{k+1}\quad\mbox{as}\quad r \rightarrow +0.$$ $\square$\

We denote the unit vectors $$\mathbf{e}_r=\frac{\partial}{\partial r},\quad \mathbf{e}_{\vartheta}=\frac{1}{r}\frac{\partial}{\partial\vartheta},
\quad \mathbf{e}_{\phi}=\frac{1}{r\sin\phi}\frac{\partial}{\partial\phi}$$ for the spherical polar co-ordinates $$\begin{aligned}
&x^1=r\sin\vartheta\cos\phi, \nonumber \\
&x^2=r\sin\vartheta\sin\phi, \nonumber \\
&x^3=r\cos\vartheta.\end{aligned}$$

Existence of spherically symmetric equilibrium for prescribed entropy distribution
==================================================================================

In this section we establish the existence of spherically symmetric equilibria which enjoy good properties used in the following consideration on $\mathbf{L}$.\
Let us put the following

\[Def.2\] A pair of $t$-independent spherically symmetric functions $(\bar{\rho}, \bar{S})
\in C_0^1(\mathbb{R}^3; [0,+\infty[)\times C^1(\mathbb{R}^3; \mathbb{R})$ which satisfies with $\mathbf{v}=0$ and $\Phi$, $P$ determined by , is called an admissible spherically symmetric equilibrium, if there is a finite positive number $R$ such that

1\) $\{\bar{\rho} >0\}=B_R $;

2\) $\bar{\rho}^{\gamma-1}, \bar{S} \in C^{\infty}(B_R) \cap
C^{2,\alpha}(\overline{B_R})$, $\alpha$ being a positive number such that $0<\alpha <\min(\frac{\gamma}{\gamma-1}-2, 1)$ ;

3\) $\displaystyle \frac{d\bar{\rho}^{\sharp}}{dr}, \frac{dP^{\sharp}}{dr} <0$ for $0<r<R$ and $$\frac{1}{r}\frac{d\bar{\rho}^{\sharp}}{dr}\Big|_{r=+0}<0, 
\quad \frac{1}{r}\frac{dP^{\sharp}}{dr}\Big|_{r=+0} <0; \label{PON}$$

4\) The boundary $\partial B_R$, on which $\bar{\rho}=0$, is a physical vacuum boundary, that is, $$-\infty <\frac{d}{dr}(\bar{\rho}^{\gamma-1})^{\sharp}\Big|_{r=R-0}<0, \label{PhysVacBd}$$ which means $$-\infty < \frac{d}{dr}\overline{\frac{\gamma P}{\rho}}^{\sharp}\Big|_{r=R-0} <0,$$ where $\overline{\gamma P/\rho}=\overline{(\partial P/\partial\rho)_{S=\mbox{Const}}}$ is the square of the sound speed.

Note that $\gamma <2$ implies $2<\frac{\gamma}{\gamma-1}$, therefore such an $\alpha$ exists. Moreover we see that $\bar{\rho}^{\gamma-1}\in C^{2,\alpha}(\overline{B_R})$ implies $\bar{\rho}^{\gamma}, \bar{P}, \bar{S} \in 
C^{2,\alpha}(\overline{B_R})$. However $\bar{\rho}\in C^{1,\alpha}(\overline{B_R})$ and $\bar{\rho} \not\in C^2(\overline{B_R})$ unless $\frac{1}{\gamma-1}\geq 2 \Leftrightarrow \gamma \leq \frac{3}{2}$.\
Note that 4) of Definition \[Def.2\] is equivalent to $$-\infty <\frac{1}{\bar{\rho}^{\sharp}}\frac{d\bar{P}^{\sharp}}{dr}\Big|_{r=R-0} <0.$$ In fact, we have $$\frac{1}{\rho}\frac{dP}{dr}=\Big[\frac{\gamma}{\gamma-1}\frac{d}{dr}\rho^{\gamma-1}+
\frac{\rho^{\gamma-1}}{\mathsf{C}_V}\frac{dS}{dr}\Big]\exp\Big(\frac{S}{\mathsf{C}_V}\Big),$$ and $\rho=\bar{\rho}^{\sharp}\rightarrow 0$ as $r\rightarrow R-0$ and $S=\bar{S}\in C^1(\bar{B}_R)$.\
We claim

\[Th.1\] Let a smooth function $\Sigma$ on $\mathbb{R}$ and a positive number $\rho_{\mathsf{O}}$ be given. Assume that it holds, for $\eta>0$, that $$\gamma +\frac{\gamma-1}{\mathsf{C}_V}\eta\frac{d}{d\eta}\Sigma(\eta) >0. \label{pDP}$$ Either if $\displaystyle\frac{4}{3}<\gamma <2$ or if $\displaystyle \frac{6}{5}<\gamma \leq \frac{4}{3}$ and $\rho_{\mathsf{O}}$ is sufficiently small, then there exists an admissible spherically symmetric equilibrium $(\bar{\rho}, \bar{S})$ such that $\bar{S}=\Sigma(\bar{\rho}^{\gamma-1})$ and $\rho(O)=\rho_{\mathsf{O}}$.

Proof . Consider the functions $f^P, f^u$ defined by $$\begin{aligned}
&f^P(\rho):=\rho^{\gamma}\exp\Big[\frac{\Sigma(\rho^{\gamma-1})}{\mathsf{C}_V}\Big], \\
&f^u(\rho):=\int_0^{\rho}
\frac{Df^P(\rho')}{\rho'}d\rho'\end{aligned}$$ for $\rho >0$. Thanks to the assumption we have $$Df^P(\rho) >0$$ for $\rho >0$, and there exists a smooth function $\Lambda$ on $\mathbb{R}$ such that $\Lambda(0)=0$ and $$f^P(\rho)=\mathsf{A}\rho^{\gamma}(1+\Lambda(\rho^{\gamma-1}))$$ for $\rho >0$. Here $\mathsf{A}:=\exp(\Sigma(0)/\mathsf{C}_V)$ is a positive constant. Then $$u=f^u(\rho)=\frac{\gamma\mathsf{A}}{\gamma-1}\rho^{\gamma-1}(1+\Lambda_u(\rho^{\gamma-1}))$$ and the inverse function $f^{\rho}$ of $f^u$ $$\rho=f^{\rho}(u)=\Big(\frac{\gamma-1}{\gamma\mathsf{A}}\Big)^{\frac{1}{\gamma-1}}
u^{\frac{1}{\gamma-1}}(1+\Lambda_{\rho}(u))$$ are given. Here $\Lambda_u, \Lambda_{\rho}$ are smooth functions on $\mathbb{R}$ such that $\Lambda_u(0)=0, \Lambda_{\rho}(0)=0$.

Therefore the problem is reduced to that for barotropic case to solve $$-\frac{1}{r^2}\frac{d}{dr}r^2\frac{du}{dr}=4\pi\mathsf{G}f^{\rho}(u), \quad
u=u_{\mathsf{O}}+O(r^2) \qquad ( r \rightarrow +0 )$$ by the shooting method. Here $u_{\mathsf{O}}=f^u(\rho_{\mathsf{O}})$ is given. The proof of the asserted result can be found in [@TM1984] and [@RendallS]. $\square$

In the barotropic case, the quantity $u$ means the specific enthalpy. But in the general baroclinic case, $u$ is not the specific enthalpy $\chi$ which should be defined as $$\chi:=\mathsf{C}_VT+\frac{P}{\rho}=\frac{\gamma}{\gamma-1}\frac{P}{\rho},$$ $T=P/((\gamma-1)\mathsf{C}_V\rho)$ being the absolute temperature. In fact we have $$\frac{d}{dr}(u-\chi)=-\frac{1}{\mathsf{C}_V}\eta\frac{d}{d\eta}\Sigma(\eta)
\Big|_{\eta=\bar{\rho}^{\gamma-1}}\overline{\frac{P}{\rho^2}\frac{d\rho}{dr} }$$ does not vanish if $\Sigma$ is not constant, that is, $\bar{S}$ is not constant, for $d\bar{\rho}/dr <0$.

Hereafter in this article we fix an admissible spherically symmetric equilibrium $(\bar{\rho}, \bar{S})$, and denote $$\rho_{\mathsf{O}}=\bar{\rho}(O)(>0),
\quad
\rho_{\mathsf{O}1}=-\lim_{r\rightarrow +0}\frac{1}{r}\frac{d\bar{\rho}}{dr}(>0), 
\quad
P_{\mathsf{O}1}=-\lim_{r\rightarrow +0}\frac{1}{r}\frac{dP}{dr}(>0).$$\
As for the Schwarzschild’s discriminant, since $\displaystyle \mathscr{A}=-\frac{1}{\gamma \mathsf{C}_V}\frac{d\bar{S}}{dr} $ and we are assuming $ \bar{S} \in 
C^{2,\alpha}(\overline{B_R})$, we have $$\mathscr{A}=O(r)$$ as $r\rightarrow +0$. (Recall Lemma \[Lem.0\].)

If $\bar{S}=\Sigma(\bar{\rho}^{\gamma-1})$, then we have $$\mathscr{A}=-\frac{1}{\gamma \mathsf{C}_V}\frac{d\bar{S}}{dr}
=-\frac{1}{\gamma \mathsf{C}_V}\frac{d}{d\eta}\Sigma(\eta)\Big|_{\eta=\bar{\rho}^{\gamma-1}}
\frac{d}{dr}\bar{\rho}^{\gamma-1}.$$ Therefore we have the following

\[Prop.1\] Suppose $\bar{S}=\Sigma(\bar{\rho}^{\gamma-1})$. If $$\frac{d}{d\eta}\Sigma(\eta) <0 \quad [\!( \frac{d}{d\eta}\Sigma(\eta) \leq 0 )\!]$$ for $\eta \geq 0$, then it holds $$\mathscr{A} <0 \quad [\!( \mathscr{A} \leq 0 )\!]$$ for $0<r\leq R$.

Self-adjoint realization of $\mathbf{L}$
========================================

We are considering the integro-differential operator $$\mathbf{L}\mbox{\boldmath$\xi$}=\frac{1}{\rho}\mathrm{grad}\delta P-
\frac{\delta\rho}{\rho^2}\mathrm{grad} P +\mathrm{grad}\delta\Phi,$$ where

$$\begin{aligned}
\delta\rho&=-\mathrm{div}(\rho\mbox{\boldmath$\xi$}), \\ 
\delta P&=\frac{\gamma P}{\rho}\delta \rho+\gamma \mathscr{A}P(\mbox{\boldmath$\xi$}|\mathbf{e}_r) ,\\ 
\delta \Phi&=-4\pi\mathsf{G}\mathcal{K}(\delta\rho). \end{aligned}$$

Here and hereafter the bars to denote the quantities evaluated at the fixed equilibrium are omitted, that is, $\rho, P$ etc stand for $\bar{\rho}, \bar{P}$ etc.\
Let us consider the operator $\mathbf{L}$ in the Hilbert space $\mathfrak{H}=L^2((B_R, \rho dx), \mathbb{C}^3)$ endowed with the norm $\|\mbox{\boldmath$\xi$}\|_{\mathfrak{H}}$ defined by $$\|\mbox{\boldmath$\xi$}\|_{\mathfrak{H}}^2=\int_{B_R}|\mbox{\boldmath$\xi$}(\mathbf{x})|^2\rho(\mathbf{x})d\mathbf{x}$$

We shall use

For complex number $z=x+\sqrt{-1}y, x, y \in \mathbb{R}$, the complex conjugate is denoted by $z^*=x-\sqrt{-1}y$. Thus, for $\displaystyle \mbox{\boldmath$\xi$}=\sum \xi^k\frac{\partial}{\partial x^k}, \xi^k \in
\mathbb{C}$, we denote $\displaystyle \mbox{\boldmath$\xi$}^*=\sum (\xi^k)^*\frac{\partial}{\partial x^k}$.

First we observe $\mathbf{L}$ restricted on $C_0^{\infty}(B_R, \mathbb{C}^3)$. Let us write $$\begin{aligned}
\mathbf{L}\mbox{\boldmath$\xi$}&=\mathbf{L}_0\mbox{\boldmath$\xi$} +4\pi\mathsf{G}\mathbf{L}_1\mbox{\boldmath$\xi$}, \\
\mathbf{L}_0\mbox{\boldmath$\xi$}&=
\mathrm{grad}\Big(-\frac{\gamma P}{\rho^2}\mathrm{div}(\rho\mbox{\boldmath$\xi$})+
\frac{\gamma\mathscr{A}P}{\rho}(\mbox{\boldmath$\xi$}|\mathbf{e}_r)\Big) + \nonumber \\
&+\frac{\gamma\mathscr{A}P}{\rho^2}\Big(-\mathrm{div}(\rho\mbox{\boldmath$\xi$})+
\frac{d\rho}{dr}(\mbox{\boldmath$\xi$}|\mathbf{e}_r)\Big)\mathbf{e}_r, \\
\mathbf{L}_1\mbox{\boldmath$\xi$}&=\mathrm{grad}\mathcal{K}(\mathrm{div}(\rho\mbox{\boldmath$\xi$})).\end{aligned}$$

Using this expression for $\mbox{\boldmath$\xi$}_{(\mu)} \in C_0^{\infty}(B_R), \mu=1,2$, we have the following formula by integration by parts: $$\begin{aligned}
(\mathbf{L}_0\mbox{\boldmath$\xi$}_{(1)}|\mbox{\boldmath$\xi$}_{(2)})_{\mathfrak{H}}&=
\int\frac{\gamma P}{\rho^2}\mathrm{div}(\rho\mbox{\boldmath$\xi$}_{(1)})\mathrm{div}(\rho\mbox{\boldmath$\xi$}_{(2)}^*) + \\
&+\int\frac{\gamma\mathscr{A}P}{\rho}
\Big[(\mbox{\boldmath$\xi$}_{(1)}|\mathbf{e}_r)\cdot\mathrm{div}(\rho\mbox{\boldmath$\xi$}_{(2)}^*)
-\mathrm{div}(\rho\mbox{\boldmath$\xi$}_{(1)})\cdot (\mbox{\boldmath$\xi$}_{(2)}|\mathbf{e}_r)^*\Big] + \\
&+\int\frac{\gamma\mathscr{A}P}{\rho}\frac{d\rho}{dr}(\mbox{\boldmath$\xi$}_{(1)}|\mathbf{e}_r)
(\mbox{\boldmath$\xi$}_{(2)}|\mathbf{e}_r)^*, \\
(\mathbf{L}_1\mbox{\boldmath$\xi$}_{(1)}|\mbox{\boldmath$\xi$}_{(2)})_{\mathfrak{H}}&=
-\int\mathcal{K}(\mathrm{div}(\rho\mbox{\boldmath$\xi$}_{(1)})\cdot \mathrm{div}(\rho\mbox{\boldmath$\xi$}_{(2)}^*).\end{aligned}$$

Thanks to the symmetry of $\mathcal{K}$, we have $$(\mathbf{L}\mbox{\boldmath$\xi$}_{(1)}|\mbox{\boldmath$\xi$}_{(2)})_{\mathfrak{H}}=
(\mbox{\boldmath$\xi$}_{(1)}|\mathbf{L}\mbox{\boldmath$\xi$}_{(2)})_{\mathfrak{H}},$$ that is, $\mathbf{L}$ restricted on $C_0^{\infty}(B_R)$ is a symmetric operator. Of course $C_0^{\infty}(B_R)$ is dense in $\mathfrak{H}$.

Moreover we have $$\begin{aligned}
(\mathbf{L}_0\mbox{\boldmath$\xi$}|\mbox{\boldmath$\xi$})&=
\int\frac{\gamma P}{\rho^2}|\mathrm{div}(\rho\mbox{\boldmath$\xi$})|^2 + \\
&+2\sqrt{-1}\mathfrak{Im}\Big[
\int\frac{\gamma\mathscr{A}P}{\rho}(\mbox{\boldmath$\xi$}|\mathbf{e}_r)\cdot\mathrm{div}(\rho\mbox{\boldmath$\xi$}^*)\Big] +
\int\frac{\gamma\mathscr{A}P}{\rho}\frac{d\rho}{dr}|(\mbox{\boldmath$\xi$}|\mathbf{e}_r)|^2.\end{aligned}$$

Since $\mathscr{A} \in C^{1,\alpha}(\overline{B_R})$, we have $$|\mathscr{A}|\sqrt{\frac{\gamma P}{\rho}}\leq C_1$$ on $0<r\leq R$, for $P/\rho =O(R-r)$. Therefore $$\begin{aligned}
\Big|\int
\frac{\gamma\mathscr{A}P}{\rho}(\mbox{\boldmath$\xi$}|\mathbf{e}_r)
\mathrm{div}(\rho\mbox{\boldmath$\xi$})^*\Big|
&\leq C_1\int
\sqrt{ \frac{\gamma P}{\rho} }|(\mbox{\boldmath$\xi$}|\mathbf{e}_r)|
|\mathrm{div}(\rho\mbox{\boldmath$\xi$})| \\
&\leq \frac{C_1}{2}\Big[
\frac{1}{\epsilon}\int \rho |(\mbox{\boldmath$\xi$}|\mathbf{e}_r)|^2+
\epsilon\int
\frac{\gamma P}{\rho^2}|\mathrm{div}(\rho\mbox{\boldmath$\xi$})|^2\Big] \\
&\leq \frac{C_1}{2}\Big[
\frac{1}{\epsilon}\|\mbox{\boldmath$\xi$}\|_{\mathfrak{H}}^2+
\epsilon\int\frac{\gamma P}{\rho^2}|\mathrm{div}\rho\mbox{\boldmath$\xi$})|^2\Big].\end{aligned}$$

Since $\displaystyle \frac{P}{\rho}\frac{d\rho }{dr}=O(\rho)$, we have $$\Big|\frac{\gamma\mathscr{A}P}{\rho}\frac{d\rho}{dr}\Big|\leq C \rho$$ Therefore we have $$\Big|\int\frac{\gamma\mathscr{A}P}{\rho}\frac{d\rho}{dr}|(\mbox{\boldmath$\xi$}|\mathbf{e}_r)|^2\Big|\leq C_2\|\mbox{\boldmath$\xi$}\|_{\mathfrak{H}}^2.$$

Thus $$(\mathbf{L}_0\mbox{\boldmath$\xi$}|\mbox{\boldmath$\xi$})_{\mathfrak{H}}\geq
\Big(1-\frac{\epsilon C_1}{2}\Big)\int
\frac{\gamma P}{\rho^2}|\mathrm{div}(\rho\mbox{\boldmath$\xi$})|^2 
-\Big(\frac{C_1}{2\epsilon}+C_2\Big)\|\mbox{\boldmath$\xi$}\|_{\mathfrak{H}}^2.$$ Taking $\epsilon$ so small that $1-\frac{\epsilon C_1}{2} \geq 0$, we get $$(\mathbf{L}_0\mbox{\boldmath$\xi$}|\mbox{\boldmath$\xi$})_{\mathfrak{H}}\geq
-\Big(\frac{C_1}{2\epsilon}+C_2\Big)\|\mbox{\boldmath$\xi$}\|_{\mathfrak{H}}^2.$$

On the other hand, it is known that $$(\mathbf{L}_1\mbox{\boldmath$\xi$}|\mbox{\boldmath$\xi$})_{\mathfrak{H}} \geq -\rho_{\mathsf{O}}\|\mbox{\boldmath$\xi$}\|_{\mathfrak{H}}
^2.$$ For a proof , see [@JJTM Proof of Proposition 2].

Summing up, $\mathbf{L}$ is bounded from below in $\mathfrak{H}$. Therefore, thanks to [@Kato Chapter VI, Section 2.3], we have

The integro-differential operator $\mathbf{L}$ on $C_0^{\infty}(B_R, \mathbb{C}^3)$ admits the Friedrichs extension, which is a self-adjoint operator, in $\mathfrak{H}$.

We want to clarify the spectral property of the self-adjoint operator $\mathbf{L}$. But this task has not yet been completely done.

At least we can claim that the spectrum of $\mathbf{L}$ cannot be of the Sturm-Liouville type in the sense defined in [@JJTM], (that is, the spectrum $\sigma(T)$ of a self-adjoint operator $T$ in a Hilbert space $\mathsf{X}$ is said to be of the Sturm-Liouville type if $\sigma(T)$ consists of isolated eigenvalues with finite multiplicities,) since $\mathrm{dim}\mathsf{N}(\mathbf{L})=\infty$, where $\mathsf{N}(\mathbf{L})$ denotes the kernel of $\mathbf{L}$. In fact, if we consider a scalar field $a$ on $B_R$ given by a function $a^{\sharp}: (r,\vartheta,\phi) \mapsto a(\mathbf{x})$ which belongs to $C_0^{\infty}(]0,R[\times ]0, \pi[\times \mathbb{S}^1)$, then the field $$\mbox{\boldmath$\xi$}=\frac{1}{\sin\vartheta}\frac{\partial a}{\partial \phi}\mathbf{e}_r
-\frac{\partial a}{\partial\vartheta}\mathbf{e}_{\phi}$$ belongs to $C_0^{\infty}(B_R)$ and satisfies $\mathrm{div}(\rho\mbox{\boldmath$\xi$})=0$ and $(\mbox{\boldmath$\xi$}|\mathbf{e}_r)=0$, therefore it belongs to $\mathsf{N}(\mathbf{L})$. Since the dimension of spaces of such $a^{\sharp}$ is infinite, we see $\mathrm{dim}\mathsf{N}(\mathbf{L})=\infty$.

In the work [@JJTM], we proved that, when $S$ is constant so that $\mathscr{A}=0$, then the spectrum of the operator $\mathbf{L}$ is $\{0\}\cup \{\lambda_n , n=1,2,\cdots \}$, where 0 is an essential spectrum and $\lambda_n$ are eigenvalues of finite multiplicities, $\lambda_n \rightarrow +\infty$ as $n \rightarrow \infty$, provided that $\mathbf{L}$ is considered in the Hilbert space $$\mathfrak{F}=\mathfrak{H}\cap\{ \mbox{\boldmath$\xi$} | 
\mathrm{div}(\rho\mbox{\boldmath$\xi$})\in \mathfrak{G}\},$$ while $$\mathfrak{G}=
L^2(B_R, \frac{1}{\rho}\frac{dP}{d\rho}d\mathbf{x})\cap\{g | \int gd\mathbf{x}=0\}.$$

So, we have the following

\[Q.1\] Even if $\mathscr{A}$ does not vanish somewhere, the spectrum of the operator $\mathbf{L}$ considered in a suitable Hilbert space, e.g. $\mathfrak{F}$, or something like it, which is dense in $\mathfrak{H}$ and in which $\mathbf{L}$ is self-adjoint, is of the form $\{0\}\cup \{\lambda_j\}, \lambda_1\leq \lambda_2 \leq \cdots \leq \lambda_n \leq \rightarrow +\infty$, where the multiplicities of $\lambda_j$ are finite?

We guess that the answer is ‘No’ when $S$ is not constant and $\mathscr{A}<0$ for $0<r<R$. In such a situation, it is expected that the so called g-modes may appear, that is, there may exist a sequence of eigenvalues which accumulates to 0. See Section 6. But a rigorous mathematical justification of this observation has not been found.

Solutions represented by spherical harmonics
============================================

In this section we consider the perturbation $\mbox{\boldmath$\xi$}$ of the particular form $$\mbox{\boldmath$\xi$}=V^r(r)Y_{lm}(\vartheta,\phi)\mathbf{e}_r+
V^h(r)\nabla_sY_{lm}+V^v(r)\nabla_s^{\perp}Y_{lm}.$$ Here $l, m \in \mathbb{Z}, 0 \leq l, |m|\leq l$, and $Y_{lm}$ is the spherical harmonics: $$\begin{aligned}
&Y_{lm}(\vartheta,\phi)=\sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}}
P_l^m(\cos\vartheta)e^{\sqrt{-1}m\phi}, \\
&Y_{l,-m}=(-1)^mY_{lm}^*\end{aligned}$$ for $m \geq 0$, while $P_l^m$ is the associated Legendre function given by $$P_l^m(\zeta)=\frac{(-1)^m}{2^ll!}(1-\zeta^2)^{m/2}\Big(\frac{d}{d\zeta}\Big)^{m+l}(\zeta^2-1)^l.$$ See [@Jackson]. We use

We denote

$$\begin{aligned}
&\nabla_sf:=\frac{\partial f}{\partial\vartheta}\mathbf{e}_{\vartheta}+
\frac{1}{\sin\vartheta}
\frac{\partial f}{\partial\phi}\mathbf{e}_{\phi}, \\
&\nabla_s^{\perp}f:=\frac{1}{\sin\vartheta}
\frac{\partial f}{\partial\phi}\mathbf{e}_{\vartheta}-\frac{\partial f}{\partial\vartheta}\mathbf{e}_{\phi}.\end{aligned}$$

Note that

$$\begin{aligned}
&\mathrm{div}(\psi(r)Y_{lm}\mathbf{e}_r)=\Big(
\frac{1}{r^2}\frac{d}{dr}r^2\frac{d\psi}{dr}\Big)Y_{lm}\mathbf{e}_r, \\
&\mathrm{div}(\psi(r)\nabla_sY_{lm})=-\frac{l(l+1)}{r}\psi(r)Y_{lm}\mathbf{e}_r, \\
&\mathrm{div}(\psi(r)\nabla_s^{\perp}Y_{lm})=0.\end{aligned}$$

\
Then , , read

$$\begin{aligned}
\delta\rho&=\delta\check{\rho}(r)Y_{lm}(\vartheta,\phi), \\
\delta P&=\delta\check{P}(r)Y_{lm}(\vartheta,\phi), \\
\delta\Phi&=\delta\check{\Phi}(r)Y_{lm}(\vartheta, \phi),\end{aligned}$$

where

$$\begin{aligned}
\delta\check{\rho}&=-\frac{1}{r^2}\frac{d}{dr}(r^2\rho V^r)+\frac{l(l+1)}{r}\rho V^h, \label{4.6a}\\
\delta\check{P}&=\frac{\gamma P}{\rho}\delta\check{\rho}+\gamma
\mathscr{A}P V^r \nonumber \\
&=-\frac{\gamma P}{r^2\rho}\frac{d}{dr}(r^2\rho V^r)+\gamma \mathscr{A}PV^r+
l(l+1)\frac{\gamma P}{r}V^h, \label{4.6b}\\
\delta\check{\Phi}&=-4\pi\mathsf{G}\mathcal{H}_l(\delta\check{\rho}) \nonumber \\
&=4\pi\mathsf{G}\mathcal{H}_l\Big(
\frac{1}{r^2}\frac{d}{dr}(r^2\rho V^r)-\frac{l(l+1)}{r}\rho V^h\Big).\label{4.6c}\end{aligned}$$

Here the integral operator $\mathcal{H}_l$ is defined by

We put $$\mathcal{H}_lf(r)=\frac{1}{2l+1}\Big[
\int_r^{\infty}f(r')\Big(\frac{r}{r'}\Big)^lr'dr'+
\int_0^rf(r')\Big(\frac{r}{r'}\Big)^{-(l+1)}r'dr'\Big],$$ provided that $f \in L^2([0,+\infty[, r^2dr)$ and $f(r)=0$ for $r\geq R$.

The equation reads $$\frac{\partial^2V^r}{\partial t^2}+L_l^r=0, \quad \frac{\partial^2V^h}{\partial t^2}+L_l^h=0, \quad \frac{\partial^2V^v}{\partial t^2} =0, \label{EqVl}$$ where

$$\begin{aligned}
&L_l^r=\frac{1}{\rho}\frac{d}{dr}\delta\check{P}-\frac{1}{\rho^2}\frac{dP}{dr}\delta\check{\rho}+\frac{d}{dr}\delta\check{\Phi}, \\
&L_l^h=\frac{1}{r}\Big(\frac{\delta\check{P}}{\rho}+\delta\check{\Phi}\Big).\end{aligned}$$

We mean $$\mathbf{L}(V^rY_{lm}\mathbf{e}_r+V^h\nabla_sY_{lm}
+V^v\nabla_s^{\perp}Y_{lm})=L_l^rY_{lm}\mathbf{e}_r+
L_l^h\nabla_sY_{lm}.$$

We are going to analyze the operator $$\vec{L}_l=\begin{bmatrix}
L_l^r \\
\\
L_l^h
\end{bmatrix}$$ which acts on $$\vec{V}=\begin{bmatrix}
V^r \\
V^h
\end{bmatrix}.$$

Actually we can neglect the component $V^v(r)$, which should be an arbitrary affine function of $t$ in order to satisfy , and we are going to consider the eigenvalue problem $$\vec{L}_l\vec{V}=\lambda\vec{V}.$$\
As for the integral operator $\mathcal{H}_l$, we shall keep in mind the following lemma, which is easy to prove:

\[Lem.2\] 1) Let $l =0$ and $f \in L^2([0,+\infty[, r^2dr), f(r)=0$ for $r \geq R$. Then $H=\mathcal{H}_0f \in C^1(]0,+\infty[)$ and satisfies $$\begin{aligned}
&H(r)=O(1), \quad \frac{d}{dr}H(r) =O(r^{-\frac{1}{2}})\quad \mbox{as}\quad r \rightarrow +0, \\
&H(r)=\frac{C}{r}\quad\mbox{as}\quad r \geq R,\end{aligned}$$ $$-\frac{1}{r^2}\frac{d}{dr}\Big(r^2\frac{dH}{dr}\Big)=f. \label{EqH0}$$ Here $C$ is the constant given by $$C=\int_0^Rf(r)r^{2}dr.$$

Conversely, if $H$ is absolutely continuous and satisfies on $]0,+\infty[$, there exist constants $C_1,C_2$ such that $$H=\mathcal{H}_0f +C_1  +\frac{C_2}{r},$$ therefore, $H\in L^{\infty}([0,+\infty[)$and $ H \rightarrow  0 $ as $ r \rightarrow +\infty$ if and only if $H=\mathcal{H}_0f$.

2\) Let $l \geq 1$ and $f \in L^2([0,+\infty[, r^2dr), f(r)=0$ for $r \geq R$. Then $H=\mathcal{H}_lf \in C^1(]0,+\infty[)$ and satisfies $$\begin{aligned}
&H(r)=O(r^{\frac{1}{2}}), \quad \frac{d}{dr}H(r)=O(r^{-\frac{1}{2}})\quad \mbox{as}\quad r \rightarrow +0, \\
&H(r)=Cr^{-(l+1)}\quad\mbox{as}\quad r \geq R,\end{aligned}$$ $$-\frac{1}{r^2}\frac{d}{dr}\Big(r^2\frac{dH}{dr}\Big)+
\frac{l(l+1)}{r^2}H=f. \label{EqHl}$$ Here $C$ is the constant given by $$C=\frac{1}{2l+1}\int_0^Rf(r)r^{l+2}dr.$$

Conversely, if $H$ is absolutely continuous and satisfies on $]0,+\infty[$, there exists constants $C_1,C_2$ such that $$H=\mathcal{H}_lf +C_1 r^l +C_2 r^{-(l+1)},$$ and $H\in L^2([0,+\infty[)$ if and only if $H=\mathcal{H}_lf$.

Case $l=0$
----------

First let us consider the case $l=0$, when only $m=0$ is possible, and $\displaystyle Y_{00}=\frac{1}{\sqrt{4\pi}}$. We are considering $$\mbox{\boldmath$\xi$}=\frac{V(r)}{\sqrt{4\pi}}\mathbf{e}_r,$$ where we write $V$ instead of $V^r$, while we need not consider $V^h$.

We are concerned with the operator $$L_0=\frac{1}{\rho}\frac{d}{dr}\delta\check{P}-\frac{1}{\rho^2}\frac{dP}{dr}\delta\check{\rho}+\frac{d}{dr}\delta\check{\Phi}$$ with

$$\begin{aligned}
\delta\check{\rho}&=-\frac{1}{r^2}\frac{d}{dr}(r^2\rho V), \\
\delta\check{P}&=\frac{\gamma P}{\rho}\delta\check{\rho}+\gamma
\mathscr{A}P V \nonumber \\
&=-\frac{\gamma P}{r^2\rho}\frac{d}{dr}(r^2\rho V)+\gamma \mathscr{A}PV, \\
\delta\check{\Phi}&=-4\pi\mathsf{G}\mathcal{H}_0(\delta\check{\rho}) \nonumber \\
&=4\pi\mathsf{G}\mathcal{H}_0\Big(
\frac{1}{r^2}\frac{d}{dr}(r^2\rho V)\Big)\end{aligned}$$

so that $$\frac{d}{dr}\delta\check{\Phi}=-4\pi\mathsf{G}\rho V.$$ We mean $$\mathbf{L}\Big(\frac{V}{\sqrt{4\pi}}\mathbf{e}_r\Big)=\frac{L_0V}{\sqrt{4\pi}}\mathbf{e}_r.$$\
Introducing the variable $\psi$ by $$V=r\psi,$$ and putting $$L^{\mathsf{ss}}\psi=\frac{1}{r}L_0(r\psi),$$ we analyze the differential operator operator $L^{\mathsf{ss}}$ in the Hilbert space $L^2([0,R], \rho r^4dr), \mathbb{C})$, since $$\|\mbox{\boldmath$\xi$}\|_{\mathfrak{H}}^2=\int_0^R|\psi(r)|^2\rho(r)r^4dr$$ for $\displaystyle \mbox{\boldmath$\xi$}=\frac{r\psi}{\sqrt{4\pi}}\mathbf{e}_r$.\
We claim

The operator $L^{\mathsf{ss}}$ on $C_0^{\infty}(]0,R[)$ admits the Friedrichs extension, a self-adjoint operator bounded from below in $L^2([0,R], \rho r^4dr)$, and its spectrum consists of simple eigenvalues $\lambda_1^{\mathsf{ss}}<\lambda_2^{\mathsf{ss}}<\cdots <
\lambda_n^{\mathsf{ss}}<\cdots \rightarrow +\infty$.

Proof. First we write $L^{\mathsf{ss}}$ as $$L^{\mathsf{ss}}\psi=-\frac{1}{r^4\rho}
\frac{d}{dr}\Big(\gamma r^4P\frac{d\psi}{dr}\Big)+q_{00}(r)\psi,$$ where $$\begin{aligned}
q_{00}(r)&=-\frac{\gamma P}{\rho^2}\frac{d^2\rho}{dr^2}
-\frac{\gamma-1}{\rho^2}\frac{dP}{dr}\frac{d\rho}{dr}
+\frac{\gamma P}{\rho^3}\Big(\frac{d\rho}{dr}\Big)^2 \nonumber \\
&-\frac{\gamma P}{r\rho^2}\frac{d\rho}{dr}
-\frac{3(\gamma-1)}{r\rho}\frac{dP}{dr}
+\frac{1}{r\rho}\frac{d}{dr}(\gamma r\mathscr{A}P)-4\pi\mathsf{G}\rho.
\label{q00}\end{aligned}$$

We see that $$|q_{00}(r)|\leq C\quad\mbox{for}\quad 0<r<R.$$ In fact, although each term in the first line of the right-hand side of is of order $(R-r)^{-1}$, these singularities are canceled after the summation, which turns out to be $$-\frac{\gamma}{\gamma-1}\Big[
\frac{d^2\eta}{dr^2}+\frac{\gamma^1}{\gamma \mathsf{C}_V}\eta^{\frac{2-\gamma}{\gamma-1}}
\frac{dS}{dr}\frac{d\eta}{dr}\Big]\exp\Big(\frac{S}{\mathsf{C}_V}\Big)=O(1),$$ where $\eta:=\rho^{\gamma-1} \in C^{2,\alpha}(\overline{B_R})$.

So, we perform the Liouville transformation of $L^{\mathsf{ss}}\psi=\lambda\psi$ to $$-\frac{d^2y}{dx^2}+q_0(x)y=\lambda y,$$ where $$\begin{aligned}
x&=\int\sqrt{\frac{c}{a}}dr,\qquad y=(ac)^{\frac{1}{4}}\psi, \\
q_0&=q_{00}+\frac{1}{4}\frac{a}{c}Q, \\
Q&=\frac{d^2}{dr^2}\log(ac)-\frac{1}{4}\Big(\frac{d}{dr}
\log(ac)\Big)^2+\Big(\frac{d}{dr}\log a\Big)\Big(
\frac{d}{dr}\log(ac)\Big)\end{aligned}$$ with $$a=\gamma r^4P,\qquad c=r^4\rho.$$ (See [@BirkhoffR p. 275, Theorem 6] or [@Yosida p.110].)

Since $$\sqrt{\frac{c}{a}}=\sqrt{\frac{\rho}{\gamma P}},$$ we can put $$x=\int_0^r\sqrt{\frac{\rho}{\gamma P}}(r')dr'$$ so that $x$ runs on the interval $[0,x_+]$, where $x_+=\int_0^R\sqrt{\rho/\gamma P}dr$. We have $$\begin{aligned}
&x \sim \sqrt{\frac{\rho}{\gamma P}}\Big|_Or\quad\mbox{as}\quad r \rightarrow +0 \\
&x_+-x \sim C_R(R-r)^{\frac{1}{2}}\quad\mbox{as}\quad r \rightarrow R-0\end{aligned}$$ with a positive constant $C_R$.

It can be verified that $$\begin{aligned}
&q_0 \sim \frac{2}{x^2} \quad\mbox{as}\quad x\rightarrow +0 \\
&q_0 \sim
\frac{(\gamma+1)(3-\gamma)}{4(\gamma-1)^2}\frac{1}{(x_+-x)^2}
\quad\mbox{as}\quad x \rightarrow x_+-0.\end{aligned}$$ Hence the assertion follows from [@ReedS p.159, Theorem X.10]. $\square$

Case $l \geq 1$
---------------

Suppose $\l \geq1$.\
Let us consider the Hilbert space $\mathfrak{X}_l$ of functions $\vec{f}=(f^r, f^h)^{\top}$ defined on $[0,R[$ endowed with the norm $\|\vec{f}\|_{\mathfrak{X}_l}$ given by

$$\|\vec{f}\|_{\mathfrak{X}_l}^2=\int_0^R
(|f^r(r)|^2+l(l+1)|f^h(r)|^2)\rho r^2dr.$$

Of course $\vec{f} =(f^r, f^h)^{\top} \in \mathfrak{X}_l$ if and only if $$\mbox{\boldmath$\xi$}=\vec{f}^{\flat/lm}:=
f^r(r)Y_{lm}\mathbf{e}_r+f^h(r)\nabla_sY_{lm} \in \mathfrak{H}$$ for $|m|\leq l$, and $$\|\vec{f}^{\flat/lm}\|_{\mathfrak{H}}=\sqrt{4\pi}\|\vec{f}\|_{\mathfrak{X}_l}.$$ we consider the operator $\vec{L}_l$ in $\mathfrak{X}_l$. We claim

\[Th.3\] The integro-differential operator $\vec{L}_l$ on $C_0^{\infty}([0,R[, \mathbb{C}^2)$ admits the Friedrichs extension, which is a self-adjoint operator bounded from below, in $\mathfrak{X}_l$.

Proof. First we look at $L_l^r, L_l^h$ by writing them as

$$\begin{aligned}
L_l^r&=-\frac{\gamma P^{\frac{1}{\gamma}}}{\rho}\frac{d}{dr}\Big(\frac{P^{1-\frac{2}{\gamma}}}{r^2}
\frac{d}{dr}(r^2P^{\frac{1}{\gamma}}V^r)\Big)
+\frac{\mathscr{A}}{\rho}\frac{dP}{dr}V^r+ \nonumber \\
&+l(l+1)\frac{\gamma P^{\frac{1}{\gamma}}}{\rho}\frac{d}{dr}
\Big(\frac{P^{1-\frac{1}{\gamma}}}{r}V^h\Big)+\frac{d\Psi}{dr}, \\
L_l^h&=-\frac{\gamma P^{1-\frac{1}{\gamma}}}{r^3\rho}
\frac{d}{dr}(r^2P^{\frac{1}{\gamma}}V^r)+
\frac{l(l+1)}{r^2}\frac{\gamma P}{\rho}V^h+\frac{\Psi}{r},\end{aligned}$$

where $$\begin{aligned}
\Psi=\delta\check{\Phi}&=-4\pi\mathsf{G}\mathcal{H}_l(\delta\check{\rho}) \nonumber \\
&=4\pi\mathsf{G}\mathcal{H}_l
\Big(\frac{1}{r^2}\frac{d}{dr}(r^2\rho V^r)-\frac{l(l+1)}{r}\rho V^h\Big).\end{aligned}$$

Using this expression, we see that the operator $\vec{L}_l$ restricted on $C_0^{\infty}([0,R[, \mathbb{C}^2)$ is symmetric and bounded from below in $\mathfrak{X}_l$.

In fact, if $\vec{V}_{(\mu)} \in C_0^{\infty}, \mu=1,2$, then the integration by parts leads us to $$\begin{aligned}
(\vec{L}_l\vec{V}_{(1)}|\vec{V}_{(2)})_{\mathfrak{X}_l} &=\gamma \int W_{(1)}W_{(2)}^*dr +\int\mathscr{A}\frac{dP}{dr}V_{(1)}^rV_{(2)}^{r*}r^2dr + \nonumber \\
&-4\pi\mathsf{G}\int\mathcal{H}_l(\delta\check{\rho}_{(1)})(\delta\check{\rho}_{(2)})^*r^2dr,\end{aligned}$$ where $$W_{(\mu)}:=\frac{P^{\frac{1}{2}-\frac{1}{\gamma}}}{r}
\frac{d}{dr}(r^2P^{\frac{1}{\gamma}}V_{(\mu)}^r)
-l(l+1)P^{\frac{1}{2}}V_{(\mu)}^h.$$ Since the integral operator $\mathcal{H}_l$ is symmetric, we see that the restriction of $\vec{L}_l$ onto $C_0^{\infty}$ is symmetric.

Let us estimate $$(\vec{L}_l\vec{V}|\vec{V})_{\mathfrak{X}_l}
=\gamma\int |W|^2dr +\int\mathscr{A}\frac{dP}{dr}|V^r|^2r^2dr
-4\pi\mathsf{G}\int\mathcal{H}_l(\delta\check{\rho})(\delta\check{\rho})^*r^2dr$$ from below for $\vec{V} \in C_0^{\infty}$.

Since $$\Big|\mathscr{A}\frac{dP}{dr}\Big|\leq C_1 \rho ,$$ we have $$\Big|\int
\mathscr{A}\frac{dP}{dr}|V^r|^2r^2dr\Big|
\leq C_1\|\vec{V}\|_{\mathfrak{X}_l}^2.$$ On the other hand, we know $$\Big|\int\mathcal{H}_l(\delta\check{\rho})(\delta\check{\rho})^*r^2dr\Big|\leq
\rho_{\mathsf{O}}\|\vec{V}\|_{\mathfrak{X}_l}^2.$$ For a proof, see [@JJTM Section 5.2]. Therefore we have $$(\vec{L}_l\vec{V}|\vec{V})_{\mathfrak{X}_l}\geq
\gamma \int |W|^2dr -(C_1+4\pi\mathsf{G}\rho_{\mathsf{O}})\|\vec{V}\|_{\mathfrak{X}_l}^2,$$ that is, $\vec{L}_l$ is bounded from below.

Therefore, thanks to [@Kato Chapter VI, Section 2.3], the restriction of $\vec{L}_l$ onto $C_0^{\infty}$ admits the Friedrichs extension. This completes the proof of Theorem \[Th.3\]. $\square$\
We see $$W=-\frac{1}{\gamma}rP^{-\frac{1}{2}}\delta\check{P}$$ by a direct calculation.\
Hereafter we denote by $\vec{L}_l$ the self-adjoint operator in $\mathfrak{X}_l$.

Note that the domain $\mathsf{D}(\vec{L}_l)$ of the Friedrichs extension is given by $$\mathsf{D}(\vec{L}_l)=\overset{\circ}{\mathfrak{W}}_l\cap\{\vec{V} \  | \  \vec{L}_l\vec{V} \in \mathfrak{X}_l\quad\mbox{in distribution sense}\quad \}.$$ Here $\overset{\circ}{\mathfrak{W}}_l$ is the closure of $C_0^{\infty}([0,R[,\mathbb{C}^2)$ in the Hilbert space $\mathfrak{W}_l$ endowed with the norm $\|\cdot\|_{\mathfrak{W}_l}$ given by $$\|\vec{V}\|_{\mathfrak{W}_l}^2=
\|\vec{V}\|_{\mathfrak{X}_l}^2+\|
\delta\check{\rho}\|_{L^2(\frac{\gamma P}{\rho^2}r^2dr)}^2, \label{NormWl}$$ where $$\|\delta\check{\rho}\|_{L^2(\frac{\gamma P}{\rho^2}r^2dr)}^2=
\int_0^R
\Big|-\frac{1}{r^2}\frac{d}{dr}(r^2\rho V^r)+
\frac{l(l+1)}{r}\rho V^h\Big|^2\frac{\gamma P}{\rho^2}r^2dr.$$

Here let us note that we can claim

\[Prop.2\] If $\vec{V}\in \mathfrak{W}_l$ satisfies $|V^r| \leq C$, then $\vec{V}$ belongs to $\overset{\circ}{\mathfrak{W}}_l$.

Proof can be done by taking $\varphi_n\in C^{\infty}([0,R])$ such that $\varphi_n(r)=1$ for $0 \leq r \leq R-\frac{1}{n},
\varphi_n(r)=0$ for $R-\frac{1}{2n} \leq r \leq R$, and $0 \leq \varphi_n\leq 1, |d\varphi_n/dr|\leq Cn$, and considering $\varphi_n\cdot \vec{V}$ for $\vec{V} \in \mathfrak{W}_l$. Let us omit the details.\
Now, actually is equivalent to $$\|\vec{V}\|_{\mathfrak{X}_l}^2+\gamma \int_0^R|W|^2dr$$ with $$W=\frac{P^{\frac{1}{2}-\frac{1}{\gamma}}}{r}
\frac{d}{dr}(r^2P^{\frac{1}{\gamma}}V^r)
-l(l+1)P^{\frac{1}{2}}V^h 
=-\frac{1}{\gamma}rP^{-\frac{1}{2}}\delta\check{P},$$ that is, is equivalent to $((\vec{L}_l+\kappa)\vec{V}|\vec{V})_{\mathfrak{X}_l}$ with a sufficiently large constant $\kappa$.\
On a parallel with the Question \[Q.1\] we have

\[Q.2\] Let $l \geq 1$. Even if $S$ is not constant, the spectral properties of $\vec{L}_l$ can be reduced to the study of an operator with spectrum of the Sturm-Liouville type as discussion of [@JJTM Section 5.2]?

We guess that the answer is ‘No’.\
As for the dimension of the kernel of $\vec{L}_l$, we have the following

1\) Let $l \geq1$. Suppose that $\mathscr{A}=0$ identically on $]0,R[$. Then $\mathrm{dim}\mathsf{N}(\vec{L}_l)=\infty$.

2\) Suppose that $\mathscr{A}<0$ everywhere on $]0,R[$. Then $\mathrm{dim}\mathsf{N}(\vec{L}_l)=0$ when $l \geq 2$ and $\mathrm{dim}\mathsf{N}(\vec{L}_1)=1$ when $l=1$.

Proof. 1) Let $l \geq 1$ and suppose that $\mathscr{A}=0$ identically. Then, if $$-\frac{1}{r^2}\frac{d}{dr}(r^2\rho V^r)+\frac{l(l+1)}{r}\rho V^h=0, \label{4.29}$$ then $\delta\check{\rho}=\delta\check{P}=\delta\check{\Phi}=0$ by , , , therefore $L_l^r=L_l^h=0$, that is, $\vec{V}\in \mathsf{N}(\vec{L}_l)$. But, if $a \in C_0^{\infty}(]0,R[)$, then $\vec{V}=(V^r,V^h)^T$ given by $$V^r(r)=\frac{1}{\rho}\cdot\frac{l(l+1)}{r^2}\int_0^ra(r')r'dr',\quad
V^h(r)=\frac{1}{\rho}a(r)$$ belongs to $C_0^{\infty}(]0,R[)$ and satisfies so that $\vec{L}_l\vec{V}=0$. Since $a \in C_0^{\infty}$ is arbitrary, we see $\mathrm{dim}\mathsf{N}(\vec{L}_l)=\infty$.

2\) Suppose that $\mathscr{A} <0$ everywhere. Let us consider $\vec{V}\in \mathsf{N}(\vec{L}_l)$, $l\geq 1$. Of course $\vec{V} \in \mathfrak{X}_l$, and moreover, since $\vec{V} \in
\mathsf{D}(\vec{L}_l)$, we have $$W=-\frac{1}{\gamma}rP^{-\frac{1}{2}}\delta\check{P}\in L^2([0,R],dr),$$ therefore $\delta\check{P}\in L^2([0,R], P^{-1}r^2dr)$. It follows that $$\frac{\gamma P}{\rho}\delta\check{\rho}=\delta\check{P}-\gamma\mathscr{A}PV^r
\in L^2([0,R], P^{-1}r^2dr),$$ therefore $$\delta\check{\rho} \in L^2([0,R], P\rho^{-2}r^2dr) \subset L^2([0,R], r^2dr),$$ since $P\rho^{-2}\geq 1/C$. Thus $\delta\check{\Phi}=-4\pi\mathsf{G}\mathcal{H}_l(\delta\check{\rho})$ enjoys the properties listed in Lemma \[Lem.2\]. Now $\vec{V}\in \mathsf{N}(\vec{L}_l)$ means

$$\begin{aligned}
&\frac{1}{\rho}\frac{d}{dr}\delta\check{P}-\frac{1}{\rho^2}\frac{dP}{dr}\check{\rho}+
\frac{d}{dr}\delta\check{\Phi}=0, \label{N4.33a}\\
&\frac{\delta\check{P}}{\rho}+\delta\check{\Phi}=0. \label{N4.33b}\end{aligned}$$

It follows from that $$\Phi'\delta\check{\rho}=\frac{d\rho}{dr}\delta\check{\Phi}, \label{4.31}$$ where we put $$\Phi':=\frac{d\Phi}{dr}=-\frac{1}{\rho}\frac{dP}{dr}.$$ On the other hand, $\delta\check{\Phi}=-4\pi\mathsf{G}
\mathcal{H}_l(\delta\check{\rho})$ implies $$\Big[\frac{1}{r^2}\frac{d}{dr}r^2\frac{d}{dr}-\frac{l(l+1)}{r^2}\Big]\delta\check{\Phi}=
4\pi\mathsf{G}\delta\check{\rho}. \label{4.33}$$ Let us consider $$X:=\frac{\delta\check{\Phi}}{\Phi'},$$ keeping in mind that $$\Phi'
=\frac{d\Phi}{dr}=\frac{4\pi\mathsf{G}}{r^2}\int_0^r\rho(r')(r')^2dr'>0
\quad\mbox{for}\quad r>0.$$ Then $\delta\check{\Phi}=\Phi'X$ and reads $\displaystyle \delta\check{\rho}=\frac{d\rho}{dr}X$. Thus, eliminating $\delta\check{\rho}$ from and , we can derive the equation $$-\frac{d}{dr}\Big(r^2(\Phi')^2\frac{dX}{dr}\Big)+(l+2)(l-1)(\Phi')^2X=0, \label{4.35}$$ using $\displaystyle \frac{1}{r^2}\frac{d}{dr}\Big(r^2\Phi'\Big)=4\pi\mathsf{G}\rho$. We owe this trick to N. R. Lebovitz, [@Lebovitz], but we are considering that the equation holds for $0<r<+\infty$ in view of $\delta\check{\rho}=\frac{d\rho}{dr}=0$ for $r \geq R$ so that we do not care the ‘boundary condition’ of $\delta\check{\Phi}$ at $r=R$.

Let us multiply by $X$ and integrate it on $[0,+\infty[$. We are going to perform the integration by parts using the following observations: For $r \geq R$ we have $$\begin{aligned}
&\Phi'=\frac{C_0}{r^2},\quad C_0=\mathrm{Const.}>0, \\
&\delta\check{\Phi}=\frac{C_1}{r^{l+1}},\quad C_1=\mathrm{Const.}, \\
&X=\frac{C_1}{C_0}\frac{1}{r^{l-1}},\end{aligned}$$ and, on the other hand, as $r\rightarrow +0$, we have $$\begin{aligned}
&\Phi'=O(r), \quad\frac{d\Phi'}{dr}=\frac{d^2\Phi}{dr^2}=O(1), \\
&\delta\check{\Phi}=O(r^{\frac{1}{2}}), \quad \frac{d}{dr}\delta\check{\Phi}=O(r^{-\frac{1}{2}}), \\
\mbox{therefore}& \\
&\frac{dX}{dr}=O(r^{-\frac{3}{2}}),
\quad
r^2(\Phi')^2\frac{dX}{dr}X=O(r^2).\end{aligned}$$ Thus the contributions from the boundaries vanish, that is, $$r^2(\Phi')^2\frac{dX}{dr}X \rightarrow 0$$ both as $r \rightarrow +0$ and as $r \rightarrow +\infty$. So the integration by parts gives $$\int_0^{+\infty}r^2(\Phi')^2\Big(\frac{dX}{dr}\Big)^2dr+
(l+2)(l-1)\int_0^{+\infty}(\Phi')^2X^2dr=0. \label{4.36}$$

Let $l \geq 2$, that is, $(l+2)(l-1)>0$. Then implies $X=0$, therefore, $\delta\check{\Phi}=\delta\check{\rho}=\delta\check{P}=0$. Then we have $$\gamma\mathscr{A}PV^r=\delta\check{P}-\frac{\gamma P}{\rho}\delta\check{\rho}=0$$ by . Since $\mathscr{A}\not=0$ everywhere, we have $V^r=0$. Since $\delta\check{\rho}=0$, this implies $V^h=0$ in view of . Thus $\vec{V}=\vec{0}$, and $\mathrm{dim}\mathsf{N}(\vec{L}_l)=0$.

Let $l=1$. Then $\displaystyle \frac{dX}{dr}=0$, that is, $X$ is a constant $\kappa$. Then we have $$\delta\check{\Phi}=-\kappa\frac{1}{\rho}\frac{dP}{dr},\quad
\delta\check{P}=\kappa\frac{dP}{dr}, \quad
\delta\check{\rho}=\kappa\frac{d\rho}{dr},$$ which imply $V^r=V^h=-\kappa$. Conversely, if $V^r=V^h=1$, then we have $$\delta\check{\rho}=-\frac{d\rho}{dr}, \quad
\delta\check{P}=-\frac{dP}{dr}, \quad
\delta\check{\Phi}=4\pi\mathsf{G}\mathcal{H}_1\Big(-\frac{d\rho}{dr}\Big)=-\frac{d\Phi}{dr},$$ and $L_1^r=L_1^h=0$. But $(1,1)^{\top} \in \mathsf{D}(\vec{L}_1)$, since $(1,1)^{\top} \in \mathfrak{W}_l$. Recall Proposition \[Prop.2\].

Summing up, we can claim $\mathrm{dim}\mathsf{N}(\vec{L}_1)=1$. This completes the proof. $\square$\
As noted in [@Lebovitz], the identity $$Y_{10}\mathbf{e}_r+\nabla_sY_{10}=\sqrt{\frac{3}{4\pi}}\frac{\partial}{\partial x^3}$$ leads us to the interpretation that the eigenfunction $\vec{V}=(1,1)^{\top}$ for $l=1$ means a uniform translation, and it can be eliminated by requiring that the center of mass remains fixed in space, or, by requiring $\displaystyle \int_{B_R}(\delta\rho)x^3d\mathbf{x}=0$.

Cowling approximation, g-modes, p-modes
=======================================

We are considering the eigenvalue problem $$\vec{L}_l\vec{V}=\lambda\vec{V}.$$ In this section we suppose $ l \geq 1 $ and consider the case $\lambda \not=0$.

Cowling approximation
---------------------

The Cowling approximation is done by neglecting the gravitational perturbations. Thus we consider the eigenvalue problem $$\vec{L}_{[\mathrm{C}]l}\vec{V}=\lambda \vec{V}, \label{5.2}$$ where $$\vec{L}_{[\mathrm{C}]l}=
\begin{bmatrix}
L_{[\mathrm{C}]l}^r \\
\\
\\
L_{[\mathrm{C}]l}^h
\end{bmatrix}
=
\begin{bmatrix}
\displaystyle \frac{1}{\rho}\frac{d}{dr}\delta\check{P} -\frac{1}{\rho^2}\frac{dP}{dr}\delta\check{\rho} \\
\\
\displaystyle \frac{1}{r\rho}\delta\check{P}
\end{bmatrix}.$$

We are going to rewrite the problem by introducing the variables $$v=r^2P^{\frac{1}{\gamma}}V^r, \qquad w=P^{-\frac{1}{\gamma}}\delta\check{P}.$$

Then we have the following expressions:

$$\begin{aligned}
V^r&=\frac{1}{r^2}P^{-\frac{1}{\gamma}}v, \label{5.5a}\\
V^h&=\frac{1}{l(l+1)}\Big[
\frac{1}{r}P^{-\frac{1}{\gamma}}\frac{dv}{dr}+\frac{1}{\gamma}rP^{-1+\frac{1}{\gamma}}
w\Big], \\
\delta\check{P}&=P^{\frac{1}{\gamma}}w, \\
\delta\check{\rho}&=\frac{1}{\gamma}\rho P^{-1+\frac{1}{\gamma}}w
-\frac{1}{r^2}\mathscr{A}\rho P^{-\frac{1}{\gamma}}v, \\
\frac{d}{dr}\delta\check{P}-\frac{1}{\rho}\frac{dP}{dr}\delta\check{\rho}&=
P^{\frac{1}{\gamma}}\frac{dw}{dr}+\frac{1}{r^2}\mathscr{A}P^{-\frac{1}{\gamma}}
\frac{dP}{dr}v. \label{5.5e}\end{aligned}$$

\
Since we are considering the case in which $\lambda \not=0$, using the expressions - , we can write the problem as

$$\begin{aligned}
&\frac{dv}{dr}=\Big(\frac{\mathscr{S}_l^2}{\lambda}-1\Big)\frac{r^2\mathscr{B}}{\mathsf{c}^2}w, \label{5.7a} \\
&\frac{dw}{dr}=(\lambda-\mathscr{N}^2)\frac{1}{r^2\mathscr{B}}v. \label{5.7b}\end{aligned}$$

Here we denote $$\mathscr{B}:=\frac{P^{\frac{2}{\gamma}}}{\rho},$$ and we use

We put $$\mathscr{S}_l^2:=\frac{l(l+1)}{r^2}\frac{\gamma P}{\rho}=\frac{l(l+1)\mathsf{c}^2}{r^2},
\qquad \mathsf{c}^2:=\frac{\gamma P}{\rho}.$$

Recall $$\mathscr{N}^2=\frac{\mathscr{A}}{\rho}\frac{dP}{dr}.$$

$\mathscr{S}_l$ is called the Lamb frequency, and $\displaystyle \mathsf{c}=\sqrt{\frac{\gamma P}{\rho}}=\sqrt{\frac{dP}{d\rho}}$ is the speed of sound when $P \propto \rho^{\gamma}$.\

Eigenvalue problem for g-modes
------------------------------

Now, so called ‘g-modes’ are given by putting $\lambda=0$ in . Then we get the approximation $$\frac{dw}{dr}+\frac{\mathscr{N}^2}{r^2\mathscr{B}}v=0 \label{5.12}$$

Now let us assume

\[A.3\] There exists a positive numbers $\delta$ such that $$\frac{1}{r}\frac{d\bar{S}}{dr}
\geq \delta$$ for $0<r <R$.

Note that, when $\bar{S}=\Sigma(\bar{\rho}^{\gamma-1})$, this Assumption \[A.3\] is fulfilled if

$$\frac{d\Sigma}{d\eta}(\rho_{\mathsf{O}}^{\gamma-1})<0$$ for $ 0\leq \eta \leq \rho_{\mathsf{O}}^{\gamma-1}$, thanks to and .\
The Assumption \[A.3\] implies $\displaystyle \mathscr{A} =-\frac{1}{\mathsf{C}_V}\frac{d\bar{S}}{dr} < 0$ on $\overline{B_{R}} \setminus \{ O\}$. So, we can take $\mathscr{N} >0$ on $\overline{B_{R}}\setminus\{ O\}$.\
Under this assumption, which guarantees $\mathscr{N}^2>0$ for $0<r<R$, inserting this into , we get the following eigenvalue problem for eigenvalue $1/\lambda$: $$-\frac{d}{dr}\Big(\frac{r^2\mathscr{B}}{\mathscr{N}^2}\frac{dw}{dr}\Big)+
\frac{r^2\rho}{\gamma P}\mathscr{B}w=
\frac{1}{\lambda}l(l+1)\mathscr{B}w. \label{5.13}$$\
Let us consider the eigenvalue problem : $$L^{\mathsf{g}}w=\frac{1}{\lambda}w,$$ where $$L^{\mathsf{g}}w=
-\frac{1}{l(l+1)\mathscr{B}}\frac{d}{dr}\Big(\frac{r^2\mathscr{B}}{\mathscr{N}^2}\frac{dw}{dr}\Big)+
\frac{1}{\mathscr{S}_l^2}w.$$\
First we note

That there is a constant $C$ such that $$\frac{\rho}{C}\leq \mathscr{B} \leq C\rho, \label{ineqB}$$ and $$\frac{r^2}{C}\leq \mathscr{N}^2 \leq C r^2. \label{ineqN}$$ Moreover $\mathscr{B}/\rho$ and $r^2/\mathscr{N}^2$ as functions of $r^2$ belong to $C^2([0, R])$.

In fact, the estimate is clear, since $$\mathscr{B}=\rho\exp\Big[\frac{2S}{\gamma \mathsf{C}_V}\Big],$$ and the estimate can be shown by the expression $$\mathscr{N}^2=-\frac{1}{\gamma \mathsf{C}_V}\frac{dS}{dr}\frac{1}{\rho}\frac{dP}{dr}
=\frac{1}{\gamma \mathsf{C}_V}\frac{dS}{dr}\frac{d\Phi}{dr}$$ under the Assumption $2^{\sharp}$ and the physical vacuum boundary condition . Moreover recall Lemma \[Lem.0\].\
We perform the Liouville transformation of to $$-\frac{d^2y}{dx^2}+q^{\mathsf{g}}y=\frac{1}{\lambda}y,$$ where $$\begin{aligned}
&x=\int \sqrt{\frac{c}{a}}dr,\qquad y=(ac)^{\frac{1}{4}}w, \\
&q^{\mathsf{g}}=\frac{b}{c}+\frac{1}{4}\frac{a}{c}Q, \\
&Q=\frac{d^2}{dr^2}\log(ac) -\frac{1}{4}\Big(\frac{d}{dr}\log(ac)\Big)^2
+\Big(\frac{d}{dr}\log a\Big)\Big(\frac{d}{dr}\log(ac)\Big),\end{aligned}$$ with $$a=\frac{r^2\mathscr{B}}{\mathscr{N}^2},\quad
b=\frac{r^2\rho \mathscr{B}}{\gamma P},
\quad c=l(l+1)\mathscr{B}.$$ (See [@BirkhoffR p. 275, Theorem 6] or [@Yosida p.110].)

Noting that $$\frac{1}{C}\leq \sqrt{\frac{c}{a}}=\sqrt{l(l+1)}\frac{\mathscr{N}}{r}\leq C,$$ we can put $$x=\int_0^r\sqrt{\frac{c}{a}}$$ so that $x$ runs on the finite interval $[0, x_+^{\mathsf{g}}]$, where $x_+^{\mathsf{g}}:=\int_0^R\sqrt{c/a}$.

Look at $$\begin{aligned}
q^{\mathsf{g}}&=\frac{r^2}{l(l+1)}\Big(\frac{\rho}{\gamma P} +\frac{Q}{4\mathscr{N}^2}\Big),\\
Q&=\frac{d^2}{dr^2}\log\frac{r^2\mathscr{B}^2}{\mathscr{N}^2}
-\frac{1}{4}\Big(\frac{d}{dr}\log\frac{r^2\mathscr{B}^2}{\mathscr{N}^2}\Big)^2
+\Big(\frac{d}{dr}\log\frac{r^2\mathscr{B}}{\mathscr{N}^2}\Big)
\Big(\frac{d}{dr}\log\frac{r^2\mathscr{B}^2}{\mathscr{N}^2}\Big).\end{aligned}$$

It can be shown that $$q^{\mathsf{g}}=O(1)\quad\mbox{as}\quad x\rightarrow +0,$$ keeping in mind that $r^2/\mathscr{N}^2 \in C^2([0,R[)$.

On the other hand, we see $$q^{\mathsf{g}} \sim \frac{2\gamma-1}{4(\gamma-1)^2}\frac{1}{(x_+^{\mathsf{g}}-x)^2}
\quad\mbox{as}\quad x \rightarrow x_+^{\mathsf{g}}-0.$$ Here, under the assumption $1<\gamma <2$, we have $\displaystyle \frac{2\gamma-1}{4(\gamma-1)^2}
>0$, and $\displaystyle \frac{2\gamma-1}{4(\gamma-1)^2}\gtrless\frac{3}{4}$ if and only if $\displaystyle \gamma \gtrless \frac{3}{2}$, when the boundary point $x_+^{\mathsf{g}}$ is of limit point/circle type. See [@ReedS p.159, Theorem X.10].

Anyway there exists a constant $C$ such that $$q^{\mathsf{g}} \geq -C \quad  \mbox{for}\quad 0<x<x_+^{\mathsf{g}}.$$

Hence we can claim that the operator $\displaystyle -\frac{d^2}{dx^2}+q^{\mathsf{g}}$ defined on $C_0^{\infty}(]0,x_+^{\mathsf{g}}[)$ admits the Friedrichs extension, a self-adjoint operator in $L^2([0,x_+^{\mathsf{g}}])$, whose spectrum consists of simple or double eigenvalues. In other words we can claim the following

The operator $L^{\mathsf{g}}$ defined on $C_0^{\infty}(]0,R[)$ admits the Friedrichs extension, a self-adjoint operator in $L^2([0,R], l(l+1)\mathscr{B}dr)$, whose spectrum consists of simple or double eigenvalues $\displaystyle \frac{1}{\lambda_n^{\mathsf{g}}},
n=1,2,\cdots$ : $\displaystyle \frac{1}{\lambda_1^{\mathsf{g}}}\leq \frac{1}{\lambda_2^{\mathsf{g}}} \leq \cdots
\leq\frac{1}{\lambda_n^{\mathsf{g}}} \leq \rightarrow +\infty$. The eigenvalues are simple if $\gamma >3/2$.

Note that $1/\lambda_1^{\mathsf{g}} >0$. In fact, if $y_1$ is an associated eigenfunction such that $\|y_1\|_{L^2([0,x_+^{\mathsf{g}}])}=1$, then $w_1=(ac)^{-\frac{1}{4}}y_1$ satisfies $$\int_0^R\Big(a(r)\Big|\frac{dw_1}{dr}\Big|^2+b(r)|w_1|^2\Big)dr=\frac{1}{\lambda_1^{\mathsf{g}}}
\int_0^R|w_1|^2c(r)dr
=\frac{1}{\lambda_1^{\mathsf{g}}}.$$ Since $ a=r^2\mathscr{B}/\mathscr{N}^2>0, b=r^2\rho \mathscr{B}/\gamma P>0$, we see $1/\lambda_1^{\mathsf{g}} >0$. Thus $$\lambda_1^{\mathsf{g}} \geq \lambda_2^{\mathsf{g}} \geq \cdots \geq \lambda_n^{\mathsf{g}} \geq \cdots
\rightarrow +0.$$\

Eigenvalue problem for p-modes
------------------------------

So called ‘p-modes’ are given by putting $1/\lambda =0$ in . Then we get the approximation $$\frac{dv}{dr}+\frac{r^2\rho}{\gamma P}\mathscr{B}w=0. \label{5.14}$$ Inserting this into , we get the following eigenvalue problem for eigenvalue $\lambda$: $$-\frac{d}{dr}\Big(\frac{\gamma P}{r^2\rho \mathscr{B}}\frac{dv}{dr}\Big)
+\frac{\mathscr{N}^2}{r^2\mathscr{B}}v=
\lambda\frac{1}{r^2\mathscr{B}}v. \label{5.15}$$\
Let us consider the eigenvalue problem : $$L^{\mathsf{p}}v=\lambda v,$$ where $$L^{\mathsf{p}}v=
-r^2\mathscr{B}\frac{d}{dr}\Big(\frac{\gamma P}{r^2\rho \mathscr{B}}\frac{dv}{dr}\Big)
+\mathscr{N}^2v.$$

We perform the Liouville transformation of to $$-\frac{d^2y}{dx^2}+q^{\mathsf{p}}y=\lambda y,$$ where $$\begin{aligned}
&x=\int\sqrt{\frac{c}{a}}dr, \qquad y=(ac)^{\frac{1}{4}}v, \\
&q^{\mathsf{p}}=\frac{b}{c}+\frac{1}{4}\frac{a}{c}Q, \\
&Q=\frac{d^2}{dr^2}\log(ac) -\frac{1}{4}\Big(\frac{d}{dr}\log(ac)\Big)^2
+\Big(\frac{d}{dr}\log a\Big)\Big(\frac{d}{dr}\log(ac)\Big),\end{aligned}$$ with $$\begin{aligned}
&a=\frac{\gamma P}{r^2\rho \mathscr{B}}=\frac{\gamma P^{1-\frac{2}{\gamma}}}{r^2\rho},\qquad
b=\frac{\mathscr{N}^2}{r^2\mathscr{B}}=\frac{\mathscr{N}^2}\rho P^{-\frac{2}{\gamma}}, \\
&c=l(l+1)\frac{1}{r^2\mathscr{B}}=l(l+1)\frac{\rho P^{-\frac{2}{\gamma}}}{r^2}.\end{aligned}$$ Noting that $$\sqrt{\frac{c}{a}}=\sqrt{\frac{\rho}{\gamma P}}\qquad
\begin{cases}
\quad\rightarrow C_0=\sqrt{\frac{\rho}{\gamma P}}(O) \quad\mbox{as}\quad r\rightarrow +0 \\
\quad \sim C_1(R-r)^{-\frac{1}{2}}\quad\mbox{as}\quad r\rightarrow R-0,
\end{cases}$$ where $C_0, C_1$ are positive constants, we can put $$x=\int_0^r\sqrt{\frac{c}{a}}$$ so that $x$ runs over the finite interval $[0,x_+^{\mathsf{p}}]$, where $x_+^{\mathsf{p}}=\int_0^R\sqrt{c/a}$, and $$\begin{aligned}
&x \sim C_0r \quad \mbox{as}\quad r \rightarrow +0 \\
&x_+^{\mathsf{p}}-x \sim 2C_1(R-r)^{\frac{1}{2}} \quad\mbox{as}\quad r \rightarrow R-0.\end{aligned}$$

Look at $$\begin{aligned}
&q^{\mathsf{p}}=\mathscr{N}^2+\frac{1}{4}\frac{\gamma P}{\rho}Q, \\
&Q=\frac{d^2}{dr^2}\log\frac{\rho P^{1-\frac{4}{\gamma}}}{r^4} 
-\frac{1}{4}\Big(\frac{d}{dr}\log\frac{\rho P^{1-\frac{4}{\gamma}}}{r^4}\Big)^2 +\Big(\frac{d}{dr}\log\frac{P^{1-\frac{2}{\gamma}}}{r^2}\Big)
\Big(\frac{d}{dr}\log\frac{\rho P^{1-\frac{4}{\gamma}}}{r^4}\Big).\end{aligned}$$ It can be shown that $$q^{\mathsf{p}} \sim \frac{2}{x^2}\quad\mbox{as}\quad x \rightarrow +0$$ and $$q^{\mathsf{p}} \sim \frac{(3-\gamma)(\gamma+1)}{4(\gamma-1)^2}\frac{1}{(x_+^{\mathsf{p}}-x)^2}
\quad\mbox{as}\quad x \rightarrow x_+^{\mathsf{p}} -0.$$ Here note that $$\frac{(3-\gamma)(\gamma+1)}{4(\gamma-1)^2}>\frac{3}{4}$$ for $1<\gamma <2$. Anyway there is a constant $C$ such that $$q^{\mathsf{p}} \geq -C \quad\mbox{for}\quad 0<x<x_+^{\mathsf{p}}.$$ The boundary points $0, x_+^{\mathsf{p}}$ are of limit point type. Hence, thanks to [@ReedS p.159, Theorem X.10], we can claim that the operator $\displaystyle -\frac{d^2}{dx^2}+q^{\mathsf{p}}$ defined on $C_0^{\infty}(]0,x_+^{\mathsf{p}}[)$ admits the Friedrichs extension, a self-adjoint operator in $L^2([0,x_+^{\mathsf{p}}])$, whose spectrum consists of simple eigenvalues. In other words we can claim the following

The operator $L^{\mathsf{p}}$ defined on $C_0^{\infty}(]0,R[)$ admits the Friedrichs extension, a self-adjoint operator in $\displaystyle L^2([0, R], \frac{1}{r^2\mathscr{B}}dr)$, whose spectrum consists of simple eigenvalues $\lambda_n^{\mathsf{p}}, n=1,2,\cdots,$ : $\lambda_1^{\mathsf{p}}<\lambda_2^{\mathsf{p}} < \cdots <\lambda_n^{\mathsf{p}}<\cdots \rightarrow +\infty$.

Note that $\lambda_1^{\mathsf{p}} >0$, since $a>0, b>0$.

Behavior of $\mathscr{S}_l^2$ and $\mathscr{N}^2$ as functions of $r$ and related speculations
----------------------------------------------------------------------------------------------

Summing up the above discussions, we have two sequences of eigenvalues $\lambda_n^{\mathsf{g}}$ and $\lambda_n^{\mathsf{p}}$, of which the former accumulates to 0 and the latter accumulates to $+\infty$. Astrophysicists believe that they give, at least approximately, eigenvalues of the problem : $\vec{L}_{[\mathrm{C}]l}\vec{V}=\lambda\vec{V}$. See, e.g., [@LedouxW], [@Cox], and so on. The priority of this observation may go back to the work by T. G. Cowling, [@Cowling], on November 3, 1941. So, we have the following

\[Q.3\] Can we perform a mathematically rigorous justification of the observation that this sequence of eigenvalues which runs to two directions and accumulates both to 0 and to $+\infty$ gives a good landscape of the real eigenvalues of $\vec{L}_l$ ?

Mathematically rigorous affirmative answer to the Question \[Q.3\] has not yet been found. But we guess that the Cowling approximation by $\vec{L}_{[\mathrm{C}]l}$ of $\vec{L}_l$ by neglecting the gravitational perturbation may be relatively easily justified. However it may be not so easy to justify the approximation of the eigenvalue problem by the g-modes / p-modes argument.

Let us reconsider the way of approximations leading the g-modes and p-modes.\

### Behavior of $\mathscr{S}_l^2$ and $\mathscr{N}^2$ 

First let us consider the behavior of the quantity $\mathscr{N}^2$ as a function of $r$.\
Here let us weaken the Assumption \[A.3\] as the following\
[**Assumption $\mathbf{2}^{\flat}$** ]{}  [*There exist positive numbers $\delta $ and $r_0 (<R)$ such that $$\frac{1}{r}\frac{d\bar{S}}{dr} \geq \delta$$ for $0<r\leq r_0$.* ]{}\
Namely we do not require that $d\bar{S}/dr > 0$ throughout $0<r<R$ but admit that $d\bar{S}/dr < 0$ so that $ \mathscr{N}^2 < 0$ for $0< R-r \ll 1$ in accordance with the more realistic view due to the helioseismology.\
Let us recall how we derived the equation of g-modes. We put $\lambda=0$ in . If $\mathscr{N}^2$ were a positive constant, $\mathscr{N}^2-\lambda$ could be replaced by $\mathscr{N}^2$ for very small $\lambda$. But now $\mathscr{N}^2$ is a function of $r$. The limit $\displaystyle \mathscr{N}_{\mathsf{O}1}^2=\lim_{r\rightarrow +0}\frac{1}{r}
\frac{d\mathscr{N}^2}{dr}$ exists. Since $$\mathscr{N}^2=-\frac{1}{\gamma \mathsf{C}_V}\frac{dS}{dr}\frac{1}{\rho}\frac{dP}{dr},$$ the Assumption $\mathbf{2}^{\flat}$ and $P_{\mathsf{O}1}>0$ imply that $\mathscr{N}_{\mathsf{O}1}^2$ is positive. Thus $\mathscr{N}^2(r)$ is strictly increasing on an interval $]0, \delta_0], \delta_0 \ll 1$ and $\mathscr{N}^2 \sim \frac{1}{2}\mathscr{N}_{\mathsf{O}1}^2r^2$ as $r\rightarrow +0$. Therefore we have positive numbers $\lambda_{\mathscr{N}}$ and $r_{\mathscr{N}}^*$ with the following property: For any $\lambda \in ]0, \lambda_{\mathscr{N}}]$ there is a unique $r_{\mathscr{N}}(\lambda) \in ]0, r_{\mathscr{N}}^*[$ such that, for $r \in ]0, r_{\mathscr{N}}^*]$, $$\lambda < \mathscr{N}^2(r) \quad \Leftrightarrow \quad
r_{\mathscr{N}}(\lambda)<r.$$

Then $\mathscr{N}^2-\lambda$ changes the sign near $0$. In other words, no matter how small $\lambda$ is, we still have $\mathscr{N}^2(r)
<\lambda$ for sufficiently small $r$. Thus the way of approximation is strongly doubtful.\
However, since $\mathscr{N}^2$ is bounded on $]0, R[$, we can generally claim at least that there exists sufficiently large $\lambda_{\mathscr{N}}^+$ such that, for any $\lambda \geq \lambda_{\mathscr{N}}^+$, it holds that $\mathscr{N}^2(r) <\lambda$ for $0<r<R$.\
On the other hand, let us consider the behavior of $\mathscr{S}_l^2$ as a function of $r$.

We see that, for $\lambda > 0$ fixed, $\displaystyle \frac{\mathscr{S}_l^2}{\lambda} -1 > 0$ for $0<r\ll 1$ and $\displaystyle \frac{\mathscr{S}_l^2}{\lambda} -1<0$ for $0<R-r \ll 1$, since $$\mathscr{S}_l^2 \sim \frac{l(l+1)\gamma P(O)}{\rho_{\mathsf{O}}}\frac{1}{r^2}
\quad\mbox{as}\quad r\rightarrow +0,$$ and $$\mathscr{S}_l^2 \sim l(l+1)C(R-r)\quad\mbox{as}\quad
r\rightarrow R-0,$$ where $C$ is a positive constant. Looking at $$\frac{1}{\mathscr{S}_l^2}\frac{d}{dr}\mathscr{S}_l^2 =\varphi(r)
+\frac{1}{\mathsf{C}_V}\frac{dS}{dr},$$ where $$\varphi(r):=-\frac{2}{r}+(\gamma-1)\frac{1}{\rho}\frac{d\rho}{dr},$$ we see that $ \varphi(r) <0$ for $0<r<R$, $\displaystyle \varphi(r) \sim -\frac{2}{r}$ as $r\rightarrow +0$, and $\displaystyle \varphi(r)\sim -\frac{\gamma(\gamma-1)}{R-r}$ as $r \rightarrow R-0$. Therefore we have $\displaystyle \frac{d}{dr}\mathscr{S}_l^2 <0$ for $0<R-r\ll 1$, since $dS/dr$ is bounded. Hence there exists $\lambda_{\mathscr{S}}>0$ such that, for any $0<\lambda \leq \lambda_{\mathscr{S}}$, there exists a unique $r_{\mathscr{S}}(\lambda) \in ]0, R[$ such that $$\mathscr{S}_l^2(r)<\lambda \quad \Leftrightarrow \quad
r_{\mathscr{S}}(\lambda)<r.$$ Then $\displaystyle \frac{\mathscr{S}_l^2}{\lambda}-1$ changes the sign near the surface $r=R$.

On the other hand, since $\displaystyle \frac{d}{dr}\mathscr{S}_l^2 <0$ for $0<r\ll 1$ and $\mathscr{S}_l^2 \rightarrow +\infty$ as $r\rightarrow +0$, there exists sufficiently large $\lambda_{\mathscr{S}}^+$ such that, for any $\lambda \geq \lambda_{\mathscr{S}}^+$, there exists unique $r_{\mathscr{S}}(\lambda) \in ]0,R[$ such that $$\mathscr{S}_l^2(r) >\lambda\quad \Leftrightarrow r< r_{\mathscr{S}}(\lambda).$$ Then $\displaystyle \frac{\mathscr{S}_l^2}{\lambda}-1$ changes the sign near the center $r=0$.

Moreover, we can claim that $\mathscr{S}_l^2$ is monotone decreasing from $+\infty$ to $0$ as $r$ varies from $0$ to $R$, provided that $dS/dr$ is sufficiently small uniformly. In fact, we see $\displaystyle \varphi(r) <-\frac{2}{R}$ for $0<r<R$, and, if $\displaystyle \sup_{0<r<R}\frac{dS}{dr} \leq \frac{2\mathsf{C}_V}{R}$, then $\displaystyle \frac{d}{dr}\mathscr{S}_l^2 <0$ for $0<r<R$. Actually this can be the case if, $4/3 < \gamma <2$ and $\rho_{\mathsf{O}}$ being fixed, we take $S=\Sigma(\rho^{\gamma-1})=-\varepsilon\rho^{\gamma-1}$ with $0<\varepsilon \ll 1$.\
Summing up, we can claim the following

\[Prop.NS\] There exist sufficiently small positive numbers $\lambda_0$ and $r_{\mathscr{N}}^*$ and a sufficiently large positive number $\lambda^+$, say, $0<\lambda_0<\lambda^+<+\infty$, such that

1\) if $0<\lambda \leq \lambda_0$, then there exist $r_{\mathscr{N}}(\lambda)
$ and $
r_{\mathscr{S}}(\lambda) $, while $0<r_{\mathscr{N}}(\lambda) <r_{\mathscr{N}}^* <
r_{\mathscr{S}}(\lambda) <R $, such that, for $0<r \leq r_{\mathscr{N}}^*$, $$\lambda -\mathscr{N}^2 (r)>0\quad\Leftrightarrow 
\quad 0<r<r_{\mathscr{N}}(\lambda),$$ and $$\frac{\mathscr{S}_l^2(r)}{\lambda}-1<0\quad\Leftrightarrow 
\quad r_{\mathscr{S}}(\lambda) < r < R;$$

2\) if $\lambda^+<\lambda$, then $$\lambda -\mathscr{N}^2(r) >0 \quad \mbox{for}\quad 0<r<R,$$ and there exists $r_{\mathscr{S}}^+(\lambda) \in
]0,R[$ such that $$\frac{\mathscr{S}_l^2(r)}{\lambda}-1>0\quad\Leftrightarrow
\quad
0<r<r_{\mathscr{S}}^+(\lambda).$$

Moreover $r_{\mathscr{N}}(\lambda) \rightarrow +0, r_{\mathscr{S}}(\lambda)\rightarrow R-0$ as $\lambda \rightarrow +0$, $r_{\mathscr{S}}^+(\lambda) \rightarrow +0$ as $ \lambda \rightarrow +\infty $, and we can suppose $$\frac{d}{dr}\mathscr{N}^2 \Big|_{r=r_{\mathscr{N}}(\lambda)}>0,\qquad
\frac{d}{dr}\mathscr{S}_l^2\Big|_{r=r_{\mathscr{S}}(\lambda)} <0 \quad\mbox{for}\quad 0<\lambda<\lambda_0$$ and $$\frac{d}{dr}\mathscr{S}_l^2\Big|_{r=r_{\mathscr{S}}^+(\lambda)} <0.
\quad\mbox{for} \quad \lambda >\lambda^+.$$

### Analysis according to H. Shibahashi

Section 16 of [@Unno2] discusses the properties of oscillations precisely for various modes by the asymptotic methods using Airy functions. Let us examine the argument found there.\
We are studying the system of equations ,, which we write as

$$\begin{aligned}
\frac{dv}{dr}&=G(r, \lambda)w, \\
\frac{dw}{dr}&=H(r, \lambda)v,\end{aligned}$$

where

$$\begin{aligned}
G&:=\Big(\frac{\mathscr{S}_l^2}{\lambda}-1\Big)\frac{r^2\mathscr{B}}{\mathsf{c}^2}, \\
H&:=(\lambda-\mathscr{N}^2)\frac{1}{r^2\mathscr{B}}.\end{aligned}$$

We put $$\mathsf{k}^2(r, \lambda):=-GH=\frac{1}{\mathsf{c}^2\lambda}(\mathscr{S}_l^2-\lambda)(\mathscr{N}^2-\lambda). \label{DefK}$$

By Proposition \[Prop.NS\], when $0<\lambda \ll 1$, $G(r)$ changes the sign at $r_{\mathscr{S}}(\lambda)$ and $H(r)$ changes the sign at $r_{\mathscr{N}}(\lambda)$, and $\mathsf{k}^2(r)<0$ for $ 0<r<r_{\mathscr{N}}(\lambda)$, $\mathsf{k}^2(r)>0$ for $r_{\mathscr{N}}(\lambda)<r
<r_{\mathscr{N}}^*$ and $\mathsf{k}^2(r) <0$ for $r_{\mathscr{S}}(\lambda)<r<r_{\mathscr{N}}^*$. On the other hand, when $1 \ll \lambda$, $\mathsf{k}^2(r)<0$ for $0<r<r_{\mathscr{S}}^+(\lambda)$ and $\mathsf{k}^2(r) >0$ for $r_{\mathscr{S}}^+(\lambda) <r <R$.

The change of variables $$v=|G|^{\frac{1}{2}}y, \qquad w=|H|^{\frac{1}{2}}z$$ leads us to the equations

$$\begin{aligned}
&-\frac{d^2y}{dr^2}+(\mathsf{F}[G]-\mathsf{k}^2)y=0, \\
&-\frac{d^2z}{dr^2}+(\mathsf{F}[H]-\mathsf{k}^2)z=0. \label{524b}\end{aligned}$$

Here the nonlinear differential operator $\mathsf{F}$ is defined by the following

For any function $Q$ of $r$ in an open interval, we put $$\begin{aligned}
\mathsf{F}[Q](r)&=|Q|^{\frac{1}{2}}\frac{d^2}{dr^2}|Q|^{-\frac{1}{2}} \nonumber \\
&=\frac{3}{4}\Big(\frac{1}{Q}\frac{dQ}{dr}\Big)^2-\frac{1}{2}\frac{1}{Q}\frac{d^2Q}{dr^2}\end{aligned}$$ for $r \in I$ such that $Q(r)\not=0$.

By Proposition \[Prop.NS\], we see $$\begin{aligned}
&\mathsf{F}[G](r) \quad\sim\quad \frac{3}{4}\frac{1}{(r-r_{\mathscr{S}}(\lambda))^2}
\quad\mbox{as}\quad r \rightarrow r_{\mathscr{S}}(\lambda), \\
&\mathsf{F}[H](r)\quad\sim\quad \frac{3}{4}\frac{1}{(r-r_{\mathscr{N}}(\lambda))^2}
\quad\mbox{as}\quad r \rightarrow r_{\mathscr{N}}(\lambda)\end{aligned}$$ when $ 0<\lambda \ll 1$, and $$\mathsf{F}[G](r)\quad\sim\quad
\frac{3}{4}\frac{1}{(r-r_{\mathscr{S}}^+(\lambda))^2}
\quad\mbox{as}\quad r\rightarrow r_{\mathscr{S}}^+(\lambda)$$ when $1 \ll \lambda$.

From Line 8 from bottom of Page 133, [@Unno2] considers the case of $\lambda \gg 1$, and, to obtain the function $z$ in the region containing the turning point $r_{\mathscr{S}}^+(\lambda)$, introduces the transformation of variables $(r, z) \mapsto (\zeta, Z)$ defined by

$$\begin{aligned}
\zeta&=\Big(\frac{dr}{d\zeta}\Big)^2(\mathsf{k}^2-\mathsf{F}[H]), \label{536a} \\
Z&=\Big|\frac{dr}{d\zeta}\Big|^{-\frac{1}{2}}z.\end{aligned}$$

In view of [@Unno2 (16.22)], it seems that is replaced by $$\zeta=\Big(\frac{dr}{d\zeta}\Big)^2\mathsf{k}^2 \label{537}$$ as a plausible approximation, in which $ \mathsf{F}[H]$ is neglected in comparison to $\mathsf{k}^2$. If we adopt this approximation, the equation for $Z$ turns out to be $$\frac{d^2Z}{d\zeta^2}+\Big(\zeta -\mathsf{F}\Big[\frac{dr}{d\zeta}\Big]\Big) Z=0,\label{528}$$ provided that the equation is replaced by $$-\frac{d^2z}{dr^2}-\mathsf{k}^2z=0 \label{529}$$ in view of the approximation of neglecting $\mathsf{F}[H]$ in comparison to $\mathsf{k}^2$. Moreover again $\displaystyle \mathsf{F}\Big[\frac{dr}{d\zeta}\Big] $ is neglected in comparison to $\zeta$ in . Then we are lead to the equation $$\frac{d^2Z}{d\zeta^2}+\zeta Z=0$$ and its solutions can be written by the Airy functions. Thus [@Unno2] continues the discussion of WBKJ asymptotic analysis.

However, it seems difficult to justify the approximation by $\mathsf{k}^2$ of $\mathsf{k}^2-\mathsf{F}[H]$ near the turning point $r_{\mathscr{S}}^+(\lambda)$ by the following reason.

Note that Proposition \[Prop.NS\] says $$\mathsf{k}^2 \quad \sim\quad C(r-r_{\mathscr{S}}^+(\lambda))$$ as $ r\rightarrow r_{\mathscr{S}}^+(\lambda)$, $C$ being a positive constant, since $\displaystyle \frac{d}{dr}\mathscr{S}_l^2 <0$ and $\mathscr{N}^2-\lambda >0$ near the point $r=r_{\mathscr{S}}^+(\lambda)$. On the other hand, since $H <0$ near the turning point, we have, at $r=r_{\mathscr{S}}^+(\lambda)$, $$\begin{aligned}
\mathsf{F}[H]&=\frac{1}{2}(\mathscr{N}^2-\lambda)^{-2}
\Big(\frac{d}{dr}\mathscr{N}^2\Big)^2
-\frac{1}{2}(\mathscr{N}^2-\lambda)^{-1}\frac{d^2}{dr^2}\mathscr{N}^2 + \\
&-\frac{1}{2}\frac{1}{r^2} +\frac{1}{2}\frac{d}{dr}\Big(\frac{1}{\mathscr{B}}\frac{d\mathscr{B}}{dr}\Big) + \\
&+\Big\{
-\frac{1}{2}(\mathscr{N}^2-\lambda)^{-1}\frac{d}{dr}\mathscr{N}^2
+\frac{1}{2}\frac{1}{r}+
\frac{1}{2}\frac{1}{\mathscr{B}}\frac{d\mathscr{B}}{dr}\Big\}^2 \\
&=\frac{1}{2}(\mathscr{N}^2-\mathscr{S}_l^2)^{-2}
\Big(\frac{d}{dr}\mathscr{N}^2\Big)^2
-\frac{1}{2}(\mathscr{N}^2-\mathscr{S}_l^2)^{-1}\frac{d^2}{dr^2}\mathscr{N}^2 + \\
&-\frac{1}{2}\frac{1}{r^2} +\frac{1}{2}\frac{d}{dr}\Big(\frac{1}{\mathscr{B}}\frac{d\mathscr{B}}{dr}\Big) + \\
&+\Big\{
-\frac{1}{2}(\mathscr{N}^2-\mathscr{S}_l^2)^{-1}\frac{d}{dr}\mathscr{N}^2
+\frac{1}{2}\frac{1}{r}+
\frac{1}{2}\frac{1}{\mathscr{B}}\frac{d\mathscr{B}}{dr}\Big\}^2 \\
&\sim\quad -\frac{1}{4r^2}\end{aligned}$$ as $r=r_{\mathscr{S}}^+(\lambda)\rightarrow +0$ for $\lambda \rightarrow +\infty$, while $\mathscr{N}^2 \rightarrow +0, \mathscr{S}_l^2\rightarrow +\infty$. Thus we see that $-\mathsf{F}[H](r_{\mathscr{S}}^+(\lambda)) \gg 1$, and it could not be neglected in comparison to $\mathsf{k}^2$ which vanishes at $r=r_{\mathscr{S}}^+(\lambda)$. In other words, the approximation [@Unno2 (16.22)], that is, $$\zeta=\frac{\mathsf{k}^2(r)}{|\mathsf{k}^2(r)|}\Big|\frac{3}{2}\int_{r_{\mathscr{S}}^+(\lambda)}^r
|\mathsf{k}^2(r')|^{\frac{1}{2}}dr'\Big|^{\frac{2}{3}}, \label{531}$$ is strongly doubtful.

Can repeating neglects of $\mathsf{F}[H]$, two times, and $\mathsf{F}[dr/d\zeta]$ produce a justifiable good approximation as the result by canceling errors in combination?

For trial let us observe the situation when we avoid the neglect of $\mathsf{F}[H]$ in comparison to $\mathsf{k}^2$. Even if we do not replace by and if we do not replace by , we get the same equation , that is, is the exact equation without approximation by neglecting $\mathsf{F}[H]$ in comparison to $\mathsf{k}^2$. But, if we use the exact , we should use it only for $r -r_{\mathscr{S}}^+>0$ which should be mapped to $\zeta >0$, since, as we know, $(\mathsf{k}^2-\mathsf{F}[H])(r) > 0$ for $|r-r_{\mathscr{S}}^+| \ll 1$. Then instead of , we have, for $0 < r- r_{\mathscr{S}}^+ \ll 1$, $$\zeta=
\Big(\frac{3}{2}\int_{r_{\mathscr{S}}^+}^r
\sqrt{\mathsf{k}^2-\mathsf{F}[H]}(r')dr'\Big)^{\frac{2}{3}}.$$ Therefore we can claim $$\begin{aligned}
& \zeta \sim \Big(\frac{3}{2}\Big)^{\frac{2}{3}}
|\mathsf{F}[H](r_{\mathscr{S}}^+)|^{\frac{1}{3}}
(r-r_{\mathscr{S}}^+)^{\frac{2}{3}}, \\
&\frac{dr}{d\zeta} \sim |\mathsf{F}[H](r_{\mathscr{S}}^+)|^{-\frac{1}{2}}
\zeta^{\frac{1}{2}}, \\
& \mathsf{F}\Big[\frac{dr}{d\zeta}\Big](\zeta) \sim \frac{5}{16}\frac{1}{\zeta^2}\end{aligned}$$ as $\zeta \rightarrow +0$. Thus the quantity $\mathsf{F}[dr/d\zeta]$ could not be neglected in comparison to $\zeta$. In other words, the approximation by for the equation is still strongly doubtful.\
Here let us remark that the corresponding argument of the first edition [@Unno] in its Section 15, which is the same as that of [@Shibahashi], considers the case of $0<\lambda \ll 1$ to investigate the g-modes near the turning point $r_{\mathscr{N}}(\lambda)$. However the same difficulty can be observed also in this discussion. The transformation $(r, y) \mapsto (\xi, Y)$ is defined by $$\begin{aligned}
\xi&=\Big(\frac{dr}{d\xi}\Big)^2(\mathsf{k}^2-\mathsf{F}[G]), \label{5.33}\\
Y&=\Big|\frac{dr}{d\xi}\Big|^{-\frac{1}{2}}y. \label{5.34}\end{aligned}$$ The term $\mathsf{F}[G]$ is neglected in comparison to $\mathsf{k}^2$ at $r=r_{\mathscr{N}}(\lambda)$ to deduce [@Unno (15.22)] or [@Shibahashi (15)], say, $$\xi=\frac{\mathsf{k}^2(r)}{|\mathsf{k}^2(r)|}
\Big|\frac{3}{2}\int_{r_{\mathscr{N}}(\lambda)}^r
|\mathsf{k}^2(r')|^{\frac{1}{2}}dr'\Big|^{\frac{2}{3}}.$$ But $\mathsf{k}^2$ vanishes at $r=r_{\mathscr{N}}(\lambda)$, while $\mathsf{F}[G](r_{\mathscr{N}}(\lambda)) \gg 1$, since we can show that, at $r=r_{\mathscr{N}}(\lambda)$, $$\mathsf{F}[G] \quad \sim\quad \frac{3}{4}\frac{1}{r^2}$$ as $r=r_{\mathscr{N}}(\lambda)\rightarrow +0$ for $\lambda \rightarrow +0$. Thus the approximation proposed here is doubtful, too. Let us observe the situation by exact use of the equation . Now, since $\mathsf{k}^2-\mathsf{F}[G] < 0$ near $r_{\mathscr{N}}$, used for $r > r_{\mathscr{N}}$ should give the mapping of $ r >r _{\mathscr{N}}$ to $\xi <0$. Therefore we have, for $0< r-r_{\mathscr{N}} \ll 1$, $$\xi =-\Big(\frac{3}{2}\int_{r_{\mathscr{N}}}^r
\sqrt{\mathsf{F}[G]-\mathsf{k}^2}(r')dr'\Big)^{\frac{2}{3}}.$$ We have $$\begin{aligned}
& \xi \sim -\Big(\frac{3}{2}\Big)^{\frac{2}{3}}(\mathsf{F}[G](r_{\mathscr{N}}))^{\frac{1}{3}}
(r-r_{\mathscr{N}})^{\frac{2}{3}}, \\
& \frac{dr}{d\xi} \sim -(\mathsf{F}[G](r_{\mathscr{N}}))^{-\frac{1}{2}}|\xi|^{\frac{1}{2}}, \\
&\mathsf{F}\Big[\frac{dr}{d\xi}\Big]=\mathsf{F}\Big[-\frac{dr}{d\xi}\Big] 
\sim \frac{5}{16}\frac{1}{|\xi|^2}\end{aligned}$$ as $\xi \nearrow 0 \quad (\Leftrightarrow r \searrow r_{\mathscr{N}})$. Thus the approximation by $\displaystyle \frac{d^2Y}{d\xi^2}+\xi Y=0$ for the equation $$\frac{d^2Y}{d\xi^2}+\Big(\xi-\mathsf{F}\Big[\frac{dr}{d\xi}\Big]\Big)Y=0$$ is doubtful, too.\

### Analysis according to D. O. Gough

Another formulation for the asymptotic analysis of the eigenvalue problem is done by D. O. Gough [@Gough]. Let us introduce it.\
We start with the eigenvalue problem under the Cowling approximation. Let us consider $\lambda \not=0$. Then the system to be considered is

$$\begin{aligned}
&\frac{1}{r^2}\frac{d}{dr}(r^2V^r)-\frac{g}{\mathsf{c}^2}V^r+
\Big(1-\frac{l(l+1)\mathsf{c}^2}{\lambda r^2}\Big)\frac{\delta\check{P}}{\rho \mathsf{c}^2}=0 \label{543.1a} \\
&\frac{d}{dr}\delta\check{P}+\frac{g^2}{\mathsf{c}^2}\delta\check{P}+
(\mathscr{N}^2-\lambda)\rho V^r=0. \label{543.1b}\end{aligned}$$

Here we put $$g:=-\frac{1}{\rho}\frac{dP}{dr}=\frac{d\Phi}{dr},\qquad \mathsf{c}^2=\frac{\gamma P}{\rho}.$$ Let us recall $$\mathscr{N}^2=-g\mathscr{A}=\frac{1}{\mathsf{C}_{V}}g\frac{dS}{dr}=g\Big(
\frac{1}{\mathsf{H}[\rho]}-\frac{g}{\mathsf{c}^2}\Big),$$ where and hereafter we use the following

For any function $Q$ of $r$ in an open interval $I$, we put $$\frac{1}{\mathsf{H}[g]}:=-\frac{1}{Q}\frac{dQ}{dr}$$ for $r \in I$ such that $Q(r)\not= 0$.

$\mathsf{H}[Q]$ is called the ‘scale hight of the quantity $Q$’. We note the relation $$\mathsf{F}[Q]=\frac{1}{4\mathsf{H}[Q]^2}\Big(1-2\frac{d}{dr}\mathsf{H}[Q]\Big).$$

Anyway the system for $(\xi, \eta)$, where $$\xi= V^r,\qquad \eta=\Delta\check{P}=\delta\check{P}+\frac{dP}{dr}V^r,$$ reads

$$\begin{aligned}
&\frac{d\xi}{dr}+A_{11}\xi+A_{12}\eta =0, \label{543.6a} \\
&\frac{d\eta}{dr}+A_{21}\xi+A_{22}\eta=0 \label{543.6b}\end{aligned}$$

with $$\begin{aligned}
&A_{11}=\frac{2}{r}-\frac{l(l+1)g}{\lambda r^2},\qquad A_{12}=
\Big(1-\frac{l(l+1)\mathsf{c}^2}{\lambda r^2}\Big)\frac{1}{\rho \mathsf{c}^2}, \nonumber \\
&A_{21}=\frac{g\rho F}{\lambda r},\qquad A_{22}\frac{l(l+1)g}{\lambda r^2},\end{aligned}$$ where we put $$F=F(\lambda):=l(l+1)\frac{g}{r}
-\Big(2+\frac{r}{\mathsf{H}[g]}\Big)\lambda
-\frac{r}{g}\lambda^2.$$

Now the system turns out to be the single second order equation for $\eta$: $$\frac{d^2\eta}{dr^2}+\frac{1}{\mathsf{H}[\mathfrak{Q}]}\frac{d\eta}{dr}+
\Big[\frac{1}{\mathsf{c}^2}\Big(\lambda+g\Big(\frac{1}{\mathsf{H}[g]}+
\frac{2}{r}\Big)\Big)-\frac{l(l+1)}{r^2}\Big(1-\frac{\mathfrak{N}^2}{\lambda}\Big)\Big]\eta=0, \label{543.9}$$ where we put $$\mathfrak{Q}:=\frac{g\rho F}{r^3},$$ and $$\mathfrak{N}^2:=g\Big(\frac{1}{\mathsf{H}[\mathfrak{Q}]}-\frac{g}{\mathsf{c}^2}-\frac{2}{\mathsf{H}[g]}
-\frac{4}{r}\Big) \quad
\Big( =\mathscr{N}^2+\frac{g}{\mathsf{H}[rF/g]}\Big). \label{543.47}$$

The second order equation is transformed to the standard form $$-\frac{d^2\Psi}{dr^2}+\mathfrak{V}(r)\Psi=0 \label{543.12}$$ by the change of variable $$\eta=\mathfrak{Q}^{\frac{1}{2}}\Psi.$$ It can be calculated that $$\mathfrak{V}(r)=\frac{1}{\mathsf{c}^2}(\lambda_c-\lambda)+\frac{l(l+1)}{r^2}
\Big(1-\frac{\mathfrak{N}^2}{\lambda}\Big),$$ where we put $$\lambda_c:=
\frac{\mathsf{c}^2}{4\mathsf{H}[\mathfrak{Q}]^2}
\Big(1-2\frac{d}{dr}\mathsf{H}[\mathfrak{Q}]\Big)
-g\Big(\frac{1}{\mathsf{H}[g]}+\frac{2}{r}\Big). \label{543.51}$$

So, D. O. Gough [@Gough] performs asymptotic analysis of this equation , which is [@Gough (5.4.7)], while $K^2$ of [@Gough] is nothing but $-\mathfrak{V}(r)$ here.

The reduction of the problem to the above standard form is exact, and, in order to perform the WKB analysis, D. O. Gough rewrites $\mathfrak{V}(r)$ as $$\begin{aligned}
-\lambda\mathsf{c}^2\mathfrak{V}(r)&=\lambda(\lambda-\lambda_c)-\mathscr{S}_l^2(\lambda-\mathfrak{N}^2) \nonumber \\
&=(\lambda-\lambda_+)(\lambda-\lambda_-),\end{aligned}$$ where $$\lambda_{\pm}=\frac{1}{2}\Big(
(\mathscr{S}_l^2+\lambda_c)\pm\sqrt{(\mathscr{S}_l^2+\lambda_c)^2-4\mathfrak{N}^2\mathscr{S}_l^2}\Big).$$ See [@Gough (5.5.1)(5.5.2)]. But we should note that $\lambda_c$ and $\mathfrak{N}$ depends not only on $r$ but also on $\lambda$. In fact, using the quantity $F(\lambda)=F(r,\lambda)$, we can write $$\begin{aligned}
\frac{1}{\mathsf{H}[\mathfrak{Q}]}&=\frac{1}{\mathsf{H}[g]}+
\frac{1}{\mathsf{H}[\rho]}+\frac{1}{\mathsf{H}[F(\lambda)]}+\frac{3}{r} \nonumber \\
&=\frac{1}{\mathsf{H}[\rho]}+\frac{1}{\mathsf{H}[gF(\lambda)/r^3]}, \\
\mathfrak{N}^2&=g\Big(-\frac{1}{\mathsf{H}[g]}+
\frac{1}{\mathsf{H}[\rho]}+\frac{1}{\mathsf{H}[F(\lambda)]}-\frac{1}{r}-\frac{g}{\mathsf{c}^2}\Big) \nonumber \\
&=\mathscr{N}^2+\frac{g}{\mathsf{H}[rF(\lambda)/g]}.\end{aligned}$$ (Here let us note that $g/r$ and $r/g$ are bounded on $0<r<R$. ) So, it is not clear whether the term $1/\mathsf{H}[F(\lambda)]$ can be neglected or, at least, can be replaced by a quantity independent of $\lambda$, or not, in order to perform the asymptotic analysis to justify the existence of g-modes.

In fact, in the preceding review [@DeubnerG] by D. O. Gough himself and in the later book [@Aerts] by C. Aerts, J. Christensen-Dalsgaard and D. W. Kurtz, the approximation $${\mathsf{H}[\mathfrak{Q}]}  \approx {\mathsf{H}[\rho]}, \quad
\mathfrak{N}^2  \approx \mathscr{N}^2, \quad
\lambda_c  \approx \frac{\mathsf{c}^2}{4\mathsf{H}[\rho]}\Big(1-2\frac{d}{dr}\mathsf{H}[\rho]\Big)$$ is adopted. If we use this approximation, we can forget the dependence on $\lambda$ of $\mathfrak{N}^2$ and $\lambda_c$. As D. O. Gough clarifies in [@Gough p. 440], this approximation is done by ‘not taking the spherical geometry fully into account’. In other words, , ‘reduce to the above approximation if $\displaystyle \frac{1}{\mathsf{H}[g]}+\frac{2}{r} \rightarrow 0$ and $\displaystyle \frac{1}{\mathsf{H}[\mathfrak{Q}]}\rightarrow\frac{1}{\mathsf{H}[\rho]}$’. It is said that ‘ $\mathfrak{N}^2$ and $\lambda_c$ are approximated well except very close to the center of the star’. So D. O.Gough calls this limits of $\lambda_c$ and $\mathfrak{N}^2$ the ‘planer values’.

Formulation as a first order system of four ordinary differential equations
===========================================================================

Let us consider the eigenvalue problem $$\vec{L}_l\vec{V}=\lambda \vec{V} \label{6.1}$$ for $ l\geq 1$. We consider $\lambda \not=0$.\
According to [@Dziembowski], we introduce the variables

$$\begin{aligned}
y_1&=\frac{V^r}{r}, \\
y_2&=-\lambda\Big(\frac{1}{\rho}\frac{dP}{dr}\Big)^{-1}V^h, \\
y_3&=-\frac{1}{r}\Big(\frac{1}{\rho}\frac{dP}{dr}\Big)^{-1}\delta\check{\Phi}, \\
y_4&=-\Big(\frac{1}{\rho}\frac{dP}{dr}\Big)^{-1}\frac{d}{dr}
\delta\check{\Phi}.\end{aligned}$$

Then the eigenvalue problem reads $$r\frac{d\vec{y}}{dr}=A(r, \lambda)\vec{y}, \label{6.3}$$ where $$\vec{y}=\begin{bmatrix}
y_1 \\
y_2 \\
y_3 \\
y_4
\end{bmatrix},
\quad
A(r,\lambda)=
\begin{bmatrix}
a_{11} & a_{12} & a_{13} & a_{14} \\
a_{21} & a_{22} & a_{23} & a_{24} \\
a_{31} & a_{32} & a_{33} & a_{34} \\
a_{41} & a_{42} & a_{43} & a_{44} 
\end{bmatrix}$$ with $$\begin{aligned}
&a_{11}=-3-\frac{r}{\gamma P}\frac{dP}{dr}, \quad a_{12}=
\Big(1-\frac{\mathscr{S}_l^2}{\lambda}\Big)\frac{r}{\gamma P}\frac{dP}{dr}, 
\quad a_{13}=-\frac{r}{\gamma P}\frac{dP}{dr}, \quad  a_{14}=0 \nonumber \\
&a_{21}=(\mathscr{N}^2-\lambda)\Big(\frac{1}{r\rho}\frac{dP}{dr}\Big)^{-1},
\quad a_{22}=-\Big(\frac{dP}{dr}\Big)^{-1}\frac{d}{dr}\Big(r\frac{dP}{dr}\Big)+
\frac{r}{\gamma P}\frac{dP}{dr}, \nonumber \\
&a_{23}=r\mathscr{A}, \quad a_{24}=0 \nonumber \\
&a_{31}=0, \quad a_{32}=0, \quad a_{33}=\frac{r^2}{\rho}\frac{dP}{dr}\frac{d}{dr}\Big(\frac{r}{\rho}\frac{dP}{dr}\Big)^{-1} , \quad a_{34}=1, \nonumber  \\
&a_{41}=4\pi\mathsf{G}r^2\rho^2\Big(\frac{dP}{dr}\Big)^{-1}\mathscr{A},
\quad a_{42}=
4\pi\mathsf{G}\frac{r^2\rho^2}{\gamma P},  \nonumber \\
&a_{43}=-4\pi\mathsf{G}\frac{r^2\rho^2}{\gamma P}+l(l+1), 
\quad a_{44}=-\Big(\frac{r}{\rho}\frac{dP}{dr}\Big)^{-1}\frac{d}{dr}
\Big(\frac{r^2}{\rho}\frac{dP}{dr}\Big).\end{aligned}$$

Here recall $$\mathscr{A}=\frac{1}{\rho}\frac{d\rho}{dr}-\frac{1}{\gamma P}\frac{dP}{dr}, \quad
\mathscr{N}^2=\mathscr{A}\frac{1}{\rho}\frac{dP}{dr}, \quad
\mathscr{S}_l^2=\frac{l(l+1)}{r^2}\frac{\gamma P}{\rho}.$$

The boundary conditions at $r=+0$, which correspond to $\vec{V}\in {\mathfrak{W}}_l$ and $\delta\check{\Phi}=-4\pi\mathsf{G}\mathcal{H}_l(\delta\check{\rho})$, read

$$\begin{aligned}
&\int_0^{\epsilon}|y_1|^2\rho r^4dr +
\frac{l(l+1)}{|\lambda|^2}\int_0^{\epsilon}
|y_2|^2\frac{1}{\rho}\Big(\frac{dP}{dr}\Big)^2r^2dr  + \nonumber \\
&+\frac{1}{\gamma}\int_0^{\epsilon}
|y_2-y_3|^2\frac{1}{P}\Big(\frac{dP}{dr}\Big)^2r^4dr < \infty, \\
&y_3=O(r^{-\frac{3}{2}}),\quad y_4=O(r^{-\frac{3}{2}})\quad\mbox{as}\quad r \rightarrow +0.\end{aligned}$$

Here $0< \epsilon \ll 1$.

On the other hand the boundary conditions at $r=R-0$ read

$$\begin{aligned}
&\int_{R-\epsilon}^{R}|y_1|^2\rho r^4dr +
\frac{l(1+1)}{|\lambda|^2}\int_{R-\epsilon}^{R}
|y_2|^2\frac{1}{\rho}\Big(\frac{dP}{dr}\Big)^2r^2dr + \nonumber \\
&+\frac{1}{\gamma}\int_{R-\epsilon}^{R}
|y_2-y_3|^2\frac{1}{P}\Big(\frac{dP}{dr}\Big)^2r^4dr < \infty, \\
& y_3(R)=\lim_{r\rightarrow R-0}y_3(r),
y_4(R)=\lim_{r\rightarrow R-0}y_4(r)\quad\mbox{exist and} \nonumber \\
&y_4(R)+(l+1)y_3(R)=0. \label{6.9b}\end{aligned}$$

Note that the condition , which means $$r\frac{d}{dr}\delta\check{\Phi}=-(l+1)\delta\check{\Phi}\quad\mbox{at}\quad r=R,$$ comes from that $\delta\check{\Phi}=-4\pi\mathsf{G}\mathcal{H}_l(\delta\check{\rho})$ should satisfy $$\delta\check{\Phi}(r)=\frac{C}{r^{l+1}} \quad\mbox{for}\quad r \geq R.$$ Otherwise, the corresponding $\delta\check{\Phi}$ for which $\displaystyle -r\frac{1}{\rho}\frac{dP}{d\rho}y_3$ is the candidate might be equal to $-4\pi\mathsf{G}\mathcal{H}_l(\delta\check{\rho})+C'r^{l}$ with $C'\not=0$.\
Now let us suppose

\[Ass.3\] $\rho$ and $S$ are analytic functions of $r^2$ near $r=0$.

This assumption holds if $S=\Sigma(\rho^{\gamma-1})$ with a function $\Sigma(\eta)$ which is analytic near $\eta=\rho_{\mathsf{O}}^{\gamma-1}$.\
Let us consider the system at $r=+0$. Under the Assumption \[Ass.3\], we see reads $$r\frac{d\vec{y}}{dr}=\Big(K_0+\sum_{m\geq 1}r^{2m}K_m\Big)\vec{y}, \label{6.10}$$ where $$K_0=
\begin{bmatrix}
-3 & \frac{l(l+1)\beta}{\lambda} & 0 & 0 \\
\frac{\lambda}{\beta} & -2 & 0 & 0 \\
0 & 0 & -2 & 1 \\
0 & 0 & l(l+1) & -3
\end{bmatrix}$$ and $\sum r^{2m}K_m$ is a convergent matrix-valued power series in $r^2$ with positive radius of convergence. Here $\beta=P_{\mathsf{O}1}/\rho_{\mathsf{O}}$ with a positive number $$P_{\mathsf{O}1}:=-\lim_{r\rightarrow +0}\frac{1}{r}\frac{dP}{dr}.$$ The eigenvalues of $K_0$ are $l-2, -(l+3)$, which are double.

Thus, putting $$z=r^2,\quad \vec{y}=
\begin{bmatrix}
l\beta & 0 & (l+1)\beta & 0 \\
\lambda & 0 & -\lambda & 0 \\
0 & 1 & 0 & 1 \\
0 & l & 0 & -(l+1)
\end{bmatrix}
\vec{w},$$ we have a system $$z\frac{d\vec{w}}{dz}=\Big(\mathbf{R}+
\sum_{m\geq 0}z^{m+1}A_m\Big)\vec{w} \label{6.14}$$ with $$\mathbf{R}=\mathrm{diag}(\rho_1, \rho_1, \rho_2, \rho_2),$$ where $$\rho_1=\frac{l-2}{2}, \qquad \rho_2=-\frac{l+3}{2}.$$ Since $\rho_1-\rho_2=l+\frac{1}{2}$ is not an integer, [@CoddingtonL Chapter 4, Theorem 4.1] gives a fundamental matrix $\Phi_w$ of the system of the form $$\Phi_w(z)=\Big(I_4+\sum_{m\geq 1}z^mP_m\Big)z^{\mathbf{R}},$$ where $\sum z^mP_m $ is a convergent matrix-valued power series. Note that $$z^{\mathbf{R}}=\mathrm{diag}(z^{\rho_1}, z^{\rho_1}, z^{\rho_2}, z^{\rho_2}).$$

As result, we have a fundamental system of solutions $\vec{y}=\mbox{\boldmath$\varphi$}_{0j}(r), j=1,2,3,4,$ of the system of the form

$$\begin{aligned}
& \mbox{\boldmath$\varphi$}_{01}=
\begin{bmatrix}
l{\beta} \\
\lambda \\
0 \\
0
\end{bmatrix}
r^{l-2}(1+O(r^2)), \\
&\mbox{\boldmath$\varphi$}_{02}=
\begin{bmatrix}
0 \\
0 \\
1 \\
l
\end{bmatrix}
r^{l-2}(1+O(r^2)), \\
&\mbox{\boldmath$\varphi$}_{03}=
\begin{bmatrix}
(l+1){\beta} \\
-\lambda  \\
0 \\
0
\end{bmatrix}
r^{-(l+3)}(1+O(r^2)), \\
&\mbox{\boldmath$\varphi$}_{04}=
\begin{bmatrix}
0 \\
0 \\
1 \\
-(l+1)
\end{bmatrix}
r^{-(l+3)}(1+O(r^2)).\end{aligned}$$

It is easy to see that only $\mbox{\boldmath$\varphi$}_{01}, \mbox{\boldmath$\varphi$}_{02}$ satisfy the boundary conditions. Therefore we have $$\vec{y}=C_{01}\mbox{\boldmath$\varphi$}_{01}+C_{02}\mbox{\boldmath$\varphi$}_{02},$$ with constants $C_{01}, C_{02}$ in order that $\vec{y}$ gives $\vec{V} \in \mathfrak{W}_l$.\
Let us consider the system at the boundary point $r=R-0$. Here we suppose

\[Ass.4\] 1) As $r \rightarrow R-0$, it holds that

$$\begin{aligned}
&\rho=C(R-r)^{\nu}\Big(1+
\sum_{k_1+k_2\geq 1}a_{k_1k_2}
(R-r)^{k_1}(R-r)^{(\nu+1)k_2}\Big), \\
&S=\sum_{k_1+k_2\geq 0}b_{k_1k_2}(R-r)^{k_1}(R-r)^{(\nu+1)k_2},\end{aligned}$$

where $C$ is a positive constant, $$\nu=\frac{1}{\gamma-1},$$ and $\sum a_{k_1k_2}X_1^{k_1}X_2^{k_2}, \sum b_{k_1k_2}X_1^{k_1}X_2^{k_2}$ are double power series with positive radii of convergence.

2\) The index $\displaystyle \nu=\frac{1}{\gamma-1}$ is a rational number.

The first part 1) of the Assumption \[Ass.4\] holds if $S=\Sigma(\rho^{\gamma-1})$ with a function $\Sigma(\eta)$ which is analytic near near $\eta=0$. For a proof, see Appendix.\
Let $\nu=\nu_N/\nu_D, \nu_N, \nu_D $ being mutually prime natural numbers.

Putting $$s=(R-r)^{1/\nu_D},$$ the Assumptions \[Ass.4\] gives that the system reads $$s\frac{d\vec{y}}{ds}=\Big(K_0^R+\sum_{m\geq 1}s^m
K_m^R\Big)\vec{y}, \label{6.22}$$ where $$K_0^R=
\begin{bmatrix}
-\nu_N & \nu_N & \nu_N & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0
\end{bmatrix}.$$ So, putting $$z=s,\qquad \vec{y}=
\begin{bmatrix}
1 & 1 & 0 & 1 \\
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0
\end{bmatrix}
\vec{w},$$ we have a system $$z\frac{d\vec{w}}{dz}=\Big(\mathbf{R}^R+\sum_{m\geq 0}z^{m+1}A_m^R\Big)\vec{w},$$ where $$\mathbf{R}^R=\mathrm{diag}(0, 0, 0, -\nu_N).$$ Here recall that $\nu_N$ is supposed to be an integer $\geq 2$. Hence we should apply the recipe prescribed in the proof of [@CoddingtonL Chapter 4, Theorem 4.2]. The result is a fundamental matrix $\Phi_y^R(s)$ of of the form $$\Phi_y^R(s)=
\begin{bmatrix}
1 & 1 & 0 & 1 \\
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0
\end{bmatrix}
\begin{bmatrix}
I_3 & \sum_{j=1}^{\nu_N}s^{j-\nu_N}\mathbf{c}_j \\
\mathbf{0}^{\top} & s^{-\nu_N}
\end{bmatrix}
(I_4+\sum_{m\geq 1}s^mP_m)
\begin{bmatrix}
I_3 & (\log s)\mathbf{b} \\
\mathbf{0}^{\top} & 1
\end{bmatrix},$$ where $\sum s^mP_m$ is a convergent matrix-valued power series, $\mathbf{c}_j$ and $\mathbf{b}$ are 3-vectors.

Thus we have solutions $\vec{y}=\mbox{\boldmath$\varphi$}_{Rj}, j=1,2,3,4$ of the form

$$\begin{aligned}
&\mbox{\boldmath$\varphi$}_{R1}=
\begin{bmatrix}
1 \\
1 \\
0 \\
0
\end{bmatrix}
+O(s), \\
&\mbox{\boldmath$\varphi$}_{R2}=
\begin{bmatrix}
1 \\
0 \\
1 \\
0
\end{bmatrix} 
+O(s), \\
&\mbox{\boldmath$\varphi$}_{R3}=
\begin{bmatrix}
0 \\
0 \\
0 \\
1
\end{bmatrix}
+O(s), \\
&\mbox{\boldmath$\varphi$}_{R4}=
\begin{bmatrix}
1 \\
0 \\
0 \\
0
\end{bmatrix}
s^{-\nu_N}(1+O(s\log s)).\end{aligned}$$

By the boundary conditions, we have $$\vec{y}=C_{R1}\mbox{\boldmath$\varphi$}_{R1}+C_{R2}(\mbox{\boldmath$\varphi$}_{R2}-(l+1)\mbox{\boldmath$\varphi$}_{R3}).$$\
Summing up, $\vec{y}$ should be $$\vec{y}=C_{01}\mbox{\boldmath$\varphi$}_{01}+C_{02}\mbox{\boldmath$\varphi$}_{02}=
C_{R1}\mbox{\boldmath$\varphi$}_{R1}+C_{R2}
(\mbox{\boldmath$\varphi$}_{R2}-(l+1)\mbox{\boldmath$\varphi$}_{R3}) \label{6.30}$$ in order that $\vec{y}$ gives $\vec{V}\in \mathfrak{W}_l$. Conversely, if so, then the corresponding $\vec{V} $ belongs to $\mathfrak{W}_l$ and we see that $V^r= ry_1$ is bounded, therefore, thanks to Proposition \[Prop.2\], $\vec{V}$ belongs to $\overset{\circ}{\mathfrak{W}}_l$ so to $\mathsf{D}(\vec{L}_l)$ as an eigenfunction.\
The condition reads $$D(r,\lambda):=\mathrm{det}(\mbox{\boldmath$\varphi$}_{01}(r,\lambda),
\mbox{\boldmath$\varphi$}_{02}(r,\lambda),
\mbox{\boldmath$\varphi$}_{R1}(r,\lambda),
\mbox{\boldmath$\varphi$}_{R2}
(r,\lambda)-(l+1)\mbox{\boldmath$\varphi$}_{R3}(r,\lambda)) =0.$$ But the condition $D(r,\lambda)=0$ is independent of $r$. So, fixing $r_0 \in ]0,R[$, we can consider $f(\lambda)=D(r_0,\lambda)$, which is a holomorphic function of $\lambda \in \mathbb{C}\setminus \{0\}$. Now we have that $f(\lambda)=0$ if and only if $\lambda$ is an eigenvalue of . Since $f(\lambda)\not=0$ for $ \lambda \notin [ -C, +\infty[ $, which $\lambda$ belongs to the resolvent set of $\vec{L}_l$, $f$ is not an identical $0$. Therefore the zeros of $f$ cannot accumulate to a value in $\mathbb{C}\setminus \{0\}$. Thus we can claim

Suppose the Assumptions \[Ass.3\],\[Ass.4\]. If there exist eigenvalues to the eigenvalue problem , then they are at most countably many eigenvalues located on the real axis, and cannot accumulate to a value $\not=0$.

Note that the above argument does not give the answer to the following

\[Q.ex\] Do exist actually eigenvalues to the problem ? In other words, does the function $f=D(r_0,\cdot)$ actually attains zeros on $\mathbb{R}\setminus \{0\}$ ?

In order to answer to this Question \[Q.ex\] we should examine the connection formula between the regular singular points $r=0, r=R$ as follows.

Since $(\mbox{\boldmath$\varphi$}_{01}, \mbox{\boldmath$\varphi$}_{02}, \mbox{\boldmath$\varphi$}_{03}, \mbox{\boldmath$\varphi$}_{04} )$ is a fundamental matrix of solutions, there are coefficients $a_{1j}(\lambda), a_{2j}(\lambda), j=1,2,3,4,$ of the connection formula such that $$\begin{aligned}
&\mbox{\boldmath$\varphi$}_{R1}=\sum_{j=1}^4a_{1j}\mbox{\boldmath$\varphi$}_{0j}, \\
&\mbox{\boldmath$\varphi$}_{R2}-(l+1)\mbox{\boldmath$\varphi$}_{R3}=
\sum_{j=1}^4a_{2j}\mbox{\boldmath$\varphi$}_{0j}.\end{aligned}$$ Then we have $$D(r,\lambda)=
\mathrm{det}
\begin{bmatrix}
a_{13}(\lambda) & a_{23}(\lambda) \\
a_{14}(\lambda) & a_{24}(\lambda)
\end{bmatrix}
\mathrm{det}(\mbox{\boldmath$\varphi$}_{01},\mbox{\boldmath$\varphi$}_{02},\mbox{\boldmath$\varphi$}_{03}, \mbox{\boldmath$\varphi$}_{04}).$$ Since $\mathrm{det}(\mbox{\boldmath$\varphi$}_{01},\mbox{\boldmath$\varphi$}_{02},\mbox{\boldmath$\varphi$}_{03}, \mbox{\boldmath$\varphi$}_{04})\not=0$, we have $$D(r,\lambda)=0 \quad \Leftrightarrow \quad
\mathrm{det}
\begin{bmatrix}
a_{13}(\lambda) & a_{23}(\lambda) \\
a_{14}(\lambda) & a_{24}(\lambda)
\end{bmatrix}
=0,$$ which is the condition for that $\lambda$ is an eigenvalue. Thus we should examine the dependence of the coefficients $a_{1j}, a_{2j}, j=3,4,$ of the global connection formula on $\lambda$. But this task seems not to be so easy.\
Anyway, suppose the answer to the Question \[Q.ex\] is ‘Yes’, and let $\{\lambda_j\}$ be the eigenvalues, which are actually infinitely many. J. Eisenfeld [@Eisenfeld] claims that the spectrum $\sigma(\vec{L}_l)$ of the self-adjoint operator $\vec{L}_l$ can be exhausted by them, that is, $\sigma(\vec{L}_l)\setminus \{0\}=\{\lambda_j\}$. However the proof of [@Eisenfeld] seems to be incomplete. Let us examine it.

Let $\lambda\not=0 $ be not an eigenvalue. Then $(\vec{L}_l-\lambda)^{-1}$ exists and the range $\mathsf{R}(\vec{L}_l-\lambda)$ is dense in $\mathfrak{X}_l$, since the residual spectrum is empty for the self-adjoint $\vec{L}_l$. The problem is:\

\[Q.4\] Does it hold that $\mathsf{R}(\vec{L}_l-\lambda) =\mathfrak{X}_l$? Or, is $(\vec{L}_l-\lambda)^{-1}$ bounded with respect to the norm $\|\cdot\|_{\mathfrak{X}_l}$ ?

If the answer is always ‘Yes’, then $\lambda \in \varrho(\vec{L}_l)$ for $\lambda$ which is neither $0$ nor an eigenvalue,, and we can claim $\sigma(\vec{L}_l)\setminus \{0\}=\{\lambda_j\}$ as J. Eisenfeld claims. Then the answer to the following question is ‘Yes’:

\[Q.5\] Do the eigenfunctions associated with the eigenvalues $\{\lambda_j\}$ form a complete orthonormal system of $\mathfrak{X}_l$ ?

Actually J. Eisenfeld [@Eisenfeld] claims the answer ‘Yes’ to the Question \[Q.5\] by referring [@DunfordS p.905, Theorem X.3.4]. However, in order to apply (the unbounded version of ) [@DunfordS Theorem X.3.4], it is necessary to prove that $\sigma(\vec{L}_l)\setminus\{0\}=
\{\lambda_j\}$, that is, to prove the answer ‘Yes’ to the Question \[Q.4\]. So let us consider the Question \[Q.4\]

Of course, if the Question \[Q.5\] is affirmatively answered, it guarantees the very existence of the eigenvalues, say affirmative answer to the Question \[Q.ex\]. But, for us, this is a kind of circular argument.

However, we should note that, maybe we should consider the spectral property of the operator $\vec{L}_l$ in a suitable Hilbert space $\mathfrak{Y}$, e.g., $\mathfrak{W}_l$, which is not equal to but dense in $\mathfrak{X}_l$ such that the restriction of $\vec{L}_l$ to $\mathfrak{Y}$ is self-adjoint. Then we should consider the Questions \[Q.4\], \[Q.5\] in this situation. In fact, for the barotropic case, according to [@JJTM], the spectral property of $\mathbf{L}$ has been clarified not in $\mathfrak{H}=L^2(B_R, \rho d\mathbf{x}) $ but in $\mathfrak{F}=\mathfrak{H}\cap
\{ \mbox{\boldmath$\xi$} |  \mathrm{div}(\rho\mbox{\boldmath$\xi$}) \in \mathfrak{G}\} $ with $\mathfrak{G}=
L^2(B_R,\frac{1}{\rho}\frac{dP}{d\rho}d\mathbf{x})\cap\{g|\int gd\mathbf{x}=0\}$ and the situation of $\vec{L}_l$ is the same. Thus maybe we should modify the Questions \[Q.4\], \[Q.5\] by replacing $\mathfrak{X}_l$ by this suitable $\mathfrak{Y}$.

In order to try to answer the Question \[Q.4\], maybe we should look at a concrete representation of the solution of the equation $$\vec{L}_l\vec{V}-\lambda\vec{V}=\vec{f}, \label{6.31}$$ where $\vec{f}=(f^r, f^h)^{\top} $ is given in $\mathsf{R}(\vec{L}_l-\lambda)$. Here we suppose that $\lambda\not=0$ and $\lambda$ is not an eigenvalue.

It can be shown that reads $$r\frac{d\vec{y}}{dr}=A(r, \lambda)\vec{y}+\vec{h}, \label{6.32}$$ where $\vec{h}=(h_1, h_2, h_3, h_4)^{\top}$ with

$$\begin{aligned}
h_1=&\frac{r\rho}{\gamma P}f^h, \\
h_2=&\Big(\frac{1}{\rho}\frac{dP}{dr}\Big)^{-1}
\Big[-f^r+
r\frac{d}{dr}f^h + (1+r\mathscr{A})f^h \Big], \\
h_3=&0, \\
h_4=&4\pi\mathsf{G}r^2\frac{\rho^2}{P}\Big(\frac{1}{\rho}\frac{dP}{d\rho}\Big)^{-1}f^h.\end{aligned}$$

Put $$\mbox{\boldmath$\Phi$}(r)=(\mbox{\boldmath$\varphi$}_{01}, \mbox{\boldmath$\varphi$}_{02},\mbox{\boldmath$\varphi$}_{R1}, \mbox{\boldmath$\varphi$}_{R2}
-(l+1)\mbox{\boldmath$\varphi$}_{R3}).$$

Since $\lambda\not=0$ is not an eigenvalue, $\mbox{\boldmath$\Phi$}(r)$ is non-singular, that is, $\mathrm{det}\mbox{\boldmath$\Phi$}(r)\not=0$. The solution of should be of the form $$\vec{y}(r)=\mbox{\boldmath$\Phi$}(r)\int_{R/2}^r\mbox{\boldmath$\Phi$}(r')^{-1}
\vec{h}(r')\frac{dr'}{r'}+\mbox{\boldmath$\Phi$}(r)\mathbf{c},$$ where $\mathbf{c}$ is a suitable constant vector, which should be determined from $\vec{h}$ of $\vec{f}$ so that the corresponding $\vec{V}$ to $\vec{y}$ belong to $\mathsf{D}(\vec{L}_l)$.

In [@Eisenfeld p.365] J. Eisenfeld claims that the constant vector $\mathbf{c}=(c_1,c_2, c_3, c_4)^{\top}$ should be chosen as $$\begin{aligned}
& c_1=-\int_{R/2}^Rk_1(r)\frac{dr}{r},\qquad 
c_2=-\int_{R/2}^Rk_2(r)\frac{dr}{r}, \\
&c_3=\int_0^Rk_3(r)\frac{dr}{r},\qquad
c_4=\int_0^Rk_4(r)\frac{dr}{r},\end{aligned}$$ where $k_j=k_i(r), i=1,2,3,4,$ is the $i$-th component of the vector $\mbox{\boldmath$\Phi$}(r)^{-1}\vec{h}(r)$, that is, $$\mbox{\boldmath$\Phi$}(r)^{-1}\vec{h}(r)=
\begin{bmatrix}
k_1 \\
k_2 \\
k_3 \\
k_4
\end{bmatrix}.$$ In other words, it is claimed that the solution $\vec{y}$ should be given as $$\begin{aligned}
\vec{y}(r)&=\Big(-\int_r^Rk_1(r')\frac{dr'}{r'}\Big)\mbox{\boldmath$\varphi$}_{01}(r) +
\Big(-\int_r^Rk_2(r')\frac{dr'}{r'}\Big)\mbox{\boldmath$\varphi$}_{02}(r) + \\
&+\Big(\int_0^rk_3(r')\frac{dr'}{r'}\Big)\mbox{\boldmath$\varphi$}_{R1}(r) +
\Big(\int_0^rk_4(r')\frac{dr'}{r'}\Big)(\mbox{\boldmath$\varphi$}_{R2}(r)-(l+1)\mbox{\boldmath$\varphi$}_{R3}(r)).\end{aligned}$$ However there we can find no persuasive argument on why the constants $c_j$ should be chosen so ?, and why $c_j$ can be chosen so?, that is, why the definite integrals $\displaystyle \int_{R/2}^Rk_1(r)\frac{dr}{r},\int_{R/2}^Rk_2(r)\frac{dr}{r},\int_0^Rk_3(r)\frac{dr}{r},
\int_0^Rk_4(r)\frac{dr}{r} $ are well-determined as finite numbers?

In order to determine $\mathbf{c}$, maybe we need asymptotic behaviors of $\mbox{\boldmath$\Phi$}(r)$ and $\mbox{\boldmath$\Phi$}(r)^{-1}$ as $r \rightarrow +0$ and as $r\rightarrow R-0$. Let us consider $r \rightarrow R-0$. Since $\mbox{\boldmath$\Phi$}(r)$ and $$\mbox{\boldmath$\Phi$}_R=(\mbox{\boldmath$\varphi$}_{R1}, \mbox{\boldmath$\varphi$}_{R2}, \mbox{\boldmath$\varphi$}_{R3}, \mbox{\boldmath$\varphi$}_{R4})$$ are fundamental matrices of the homogeneous system, there exists a constant matrix $C=(c_{ij})_{ij}$ such that $$\mbox{\boldmath$\Phi$}(r)=\mbox{\boldmath$\Phi$}_R(r)C.$$ See [@CoddingtonL Chapter 1, Theorem 7.3]. Of course $C$ is of the form $$C=
\begin{bmatrix}
c_{11} & c_{12} & 1 & 0 \\
c_{21} & c_{22} & 0 & 1 \\
c_{31} & c_{32} & 0 &-(l+1) \\
c_{41} & c_{42} & 0 & 0
\end{bmatrix}.$$ We should determine $c_{ij}, i=1,2,3,4, j=1,2$. But J. Eisenfeld claims a relation which can be interpreted as $$\begin{aligned}
\mbox{\boldmath$\varphi$}_{01}&=c_{31}\mbox{\boldmath$\varphi$}_{R3}+c_{41}\mbox{\boldmath$\varphi$}_{R4}, \\
\mbox{\boldmath$\varphi$}_{02}&=c_{32}\mbox{\boldmath$\varphi$}_{R3}+c_{42}\mbox{\boldmath$\varphi$}_{R4}\end{aligned}$$ at [@Eisenfeld p. 366, line 7]. In other words it is claimed that $c_{ij}=0$ for $i,j=1,2$. However we cannot find its proof there. It seems to be doubtful. Thus the asymptotic behaviors of the matrix given by [@Eisenfeld (5.11), (5.12)], which can be interpreted as the asymptotic behaviors of $\mbox{\boldmath$\Phi$}(r)^{-1}$ as $r \rightarrow +0$ and as $r\rightarrow R-0$, seem to be not well-grounded.

As conclusion, the argument by J. Eisenfeld [@Eisenfeld] is too weak to determine the answer to the Question \[Q.4\]. The problem is still open.\
[**Supplementary Remark** ]{}\
The same speculation can be done for the system of equations , under the Cowling approximation. Namely let us consider the system

$$\begin{aligned}
&\frac{dv}{dr}=\frac{1}{\lambda}(\mathscr{S}_l^2-\lambda)
\frac{r^2\mathscr{B}}{\mathsf{c}^2}w \label{SR201a} \\
&\frac{dw}{dr}=(\lambda-\mathscr{N}^2)\frac{1}{r^2\mathscr{B}}v \label{SR201b}\end{aligned}$$

for the variables $$v=r^2P^{\frac{1}{\gamma}}V^r,\quad w=P^{-\frac{1}{\gamma}}\delta\check{P}.$$ Recall $$\begin{aligned}
&\mathscr{B}=\frac{P^{\frac{2}{\gamma}}}{\rho},\quad
\mathscr{S}_l^2=\frac{l(l+1)}{r^2}\frac{\gamma P}{\rho},
\quad \mathsf{c}^2=\frac{\gamma P}{\rho}, \\
&\mathscr{N}^2=\frac{\mathscr{A}}{\rho}\frac{dP}{d\rho}=
-\frac{1}{\gamma \mathsf{C}_V}\frac{dS}{dr}\frac{1}{\rho}\frac{dP}{d\rho}=
\frac{1}{\gamma\mathsf{C}_V}\frac{dS}{dr}\frac{d\Phi}{dr}.\end{aligned}$$

The boundary conditions for $(v,w)$ at the boundaries $r=0, r=R$ can be derived from the condition that the corresponding $\vec{V}=(V^r,V^h)^{\top}$ belong to $\overset{\circ}{\mathfrak{W}}_l$. It reads $$\|\vec{V}\|_{\mathfrak{X}_l}^2=
\int_0^R|v|^2\frac{dr}{\mathscr{B}}+
\frac{1}{l(l+1)}\int_0^R\Big|\frac{dv}{dr}+
\frac{r^2\mathscr{B}}{\mathsf{c}^2}w\Big|^2\frac{dr}{\mathscr{B}} <\infty,$$ and $$\gamma\int_0^R|W|^2dr=
\int_0^R|w|^2\frac{r^2\mathscr{B}}{\mathsf{c}^2}dr <\infty,$$ provided that $r^2P^{-\frac{1}{\gamma}}|v|$ is bounded on $]0,R[$.

Let us consider the boundary $r=0$. The system , can be written $$r\frac{d\vec{u}}{dr}=
\begin{bmatrix}
-1 & \frac{1}{\lambda}(\mathscr{S}_l^2-\lambda)\frac{r^2\mathscr{B}}{\mathsf{c}^2} \\
\frac{\lambda-\mathscr{N}^2}{\mathscr{B}} & 0
\end{bmatrix}
\vec{u},$$ where $\displaystyle \vec{u}=(u,w)^{\top}, u=\frac{v}{r}$. Under the Assumption \[Ass.3\], we see $$\begin{bmatrix}
-1 & \frac{1}{\lambda}(\mathscr{S}_l^2-\lambda)\frac{r^2\mathscr{B}}{\mathsf{c}^2} \\
\frac{\lambda-\mathscr{N}^2}{\mathscr{B}} & 0
\end{bmatrix}
=
\begin{bmatrix}
-1 & \frac{l(l+1)B_{\mathsf{O}}}{\lambda} \\
\frac{\lambda}{B_{\mathsf{O}}} & 0
\end{bmatrix}
+[r^2]_1,$$ where $B_{\mathsf{O}}=\mathscr{B}\Big|_{r=0}=P_{\mathsf{O}}^{\frac{2}{\gamma}}/\rho_{\mathsf{O}}$. Therefore there is a fundamental system of solutions $\vec{u}_{01}=(u_{01},w_{01})^{\top},
\vec{u}_{02}=(u_{02}, w_{02})^{\top}$ such that $$\vec{u}_{01}\sim 
\begin{bmatrix}
lB_{\mathsf{O}} \\
\lambda
\end{bmatrix}
r^l,\qquad
\vec{u}_{02}\sim
\begin{bmatrix}
(l+1)B_{\mathsf{O}} \\
-\lambda
\end{bmatrix}
r^{-(l+1)}$$ as $r\rightarrow +0$. Put $\mbox{\boldmath$\varphi$}_{0j}(r)=(ru_{0j}(r), w_{0j}(r))^{\top}, j=1,2$. Clearly only $\mbox{\boldmath$\varphi$}_{01}$ is the solution of which satisfies the boundary condition.

As for the boundary $r=R$, we consider the system $$s\frac{d\vec{U}}{ds}=
\begin{bmatrix}
-\nu_N          &
-\nu_D\Big(\frac{\mathscr{S}_l^2}{\lambda}-1\Big)s^{-\nu_N+\nu_D}\frac{r^2\mathscr{B}}{\mathsf{c}^2} \\
-\nu_Ds^{\nu_N+\nu_D}(\lambda-\mathscr{N}^2)\frac{1}{r^2\mathscr{B}}& 0
\end{bmatrix}
\vec{U}$$ for $s=(R-r)^{1/\nu_D}, \vec{U}=(U, w)^{\top}, U=s^{-\nu_N}v$. Under the Assumptions \[Ass.3\] and \[Ass.4\], we see $$\begin{bmatrix}
-\nu_N          &
-\nu_D\Big(\frac{\mathscr{S}_l^2}{\lambda}-1\Big)s^{-\nu_N+\nu_D}\frac{r^2\mathscr{B}}{\mathsf{c}^2} \\
-\nu_Ds^{\nu_N+\nu_D}(\lambda-\mathscr{N}^2)\frac{1}{r^2\mathscr{B}}& 0
\end{bmatrix}
=\begin{bmatrix}
-\nu_N & b \\
0 & 0
\end{bmatrix}
+[s]_1,$$ where $$b:=\lim_{s\rightarrow +0}\nu_Ds^{-\nu_N+\nu_D}
\frac{r^2\mathscr{B}}{\mathsf{c}^2}\Big|_{r=R-s} >0$$ Thus there is a fundamental system of solutions $\vec{U}_{Rj}=(U_{Rj}, w_{Rj})^{\top}, j=1,2$ such that $$\vec{U}_{R1}(r) \sim
\begin{bmatrix}
b \\
\nu_N
\end{bmatrix}
,
\quad
\vec{U}_{R2}\sim
\begin{bmatrix}
1 \\
0
\end{bmatrix}
s^{-\nu_N}$$ as $s \rightarrow +0$. Put $\mbox{\boldmath$\varphi$}_{Rj}=((R-r)^{\nu_N}U_{Rj}, w_{Rj})^{\top}, j=1,2$. Clearly only $\mbox{\boldmath$\varphi$}_{R1}$ is the solution of which satisfies the boundary condition.

Therefore the eigenfunction, if exists, should be $$\mbox{\boldmath$\varphi$}=C_1\mbox{\boldmath$\varphi$}_{01}=C_2\mbox{\boldmath$\varphi$}_{R1}.$$ That is, if and only if $$D(r,\lambda)=\mathrm{det}(\mbox{\boldmath$\varphi$}_{01}(r,\lambda),\mbox{\boldmath$\varphi$}_{R1}(r,\lambda))=0 
\quad\mbox{for}\quad \forall / \exists r \in ]0,R[,$$ $\lambda$ is an eigenvalue. So, the question is

Do there exist a sequence of infinitely many zeros of $D(r_0,\cdot)$ which accumulates to $0$ and $+\infty$? To begin with, does $D(r_0,\cdot)$ admit at least one zero?

Of course there are determined the coefficients of connection of the two regular singular points $c_j(\lambda), j=1,2$ such that $$\mbox{\boldmath$\varphi$}_{01}(r,\lambda)=c_1(\lambda)\mbox{\boldmath$\varphi$}_{R1}(r,\lambda)+
c_2(\lambda)\mbox{\boldmath$\varphi$}_{R2}(r,\lambda),$$ and then we have $$D(r_0,\lambda)=0\quad\Leftrightarrow\quad c_2(\lambda)=0.$$ The dependence of the connection coefficient $c_2(\lambda)$ on $\lambda$ is problematic. Moreover, if $\lambda\not=0$ is not an eigenvalue, then

Does it hold $\mathsf{R}(\vec{L}_{[\mathrm{C}]l}-\lambda)=\mathfrak{X}_l$?

[**Acknowledgment**]{}\
This work is financially supported by JSPS KAKENHI Grant Number JP18K03371.\
[**Appendix**]{}\
Let us consider a function $f$ of the form $$f(u)=Ku^{\nu}(1+\sum_{k\geq 1}f_ku^k)$$ as $u \rightarrow +0$, where $K$ is a positive constant, $1<\nu <+\infty$, $\sum f_ku^k$ is a convergent power series, while $f(u)>0$ for $u >0$. Suppose $u=u(r), 0<r <R,$ satisfies $u >0, du/dr <0,$ $$\frac{d^2u}{dr^2}+\frac{2}{r}\frac{du}{dr}+f(u)=0$$ on $0<r<R$, $u \rightarrow +0$ as $r\rightarrow R-0$, and the limit $\displaystyle  \lim_{r\rightarrow R-0}\frac{du}{dr} $ exists to be finite and strictly negative.

Then we have the expansion $$u=C\frac{R-r}{R}
\Big[
1+
\sum_{k_1+k_2+k_3\geq 1}b_{k_1k_2k_3}
\Big(\frac{R-r}{R}\Big)^{k_1}
\Big(C'\Big(\frac{R-r}{R}\Big)^{\nu+1}\Big)^{k_2}
\Big(C\frac{R-r}{R}\Big)^{k_3}\Big],$$ where $C$ is a positive constant, $C'=R^2KC^{\nu-1}$, and $\sum b_{k_1k_2k_3}X_1^{k_1}X_2^{k_2}X_3^{k_3}$ is a convergent triple power series.\
Let us sketch the proof. First note that there is a convergent power series $$G(u)=\sum_{k\geq 1}G_ku^k$$ such that $$u\frac{d}{du}\log f(u)=\nu+G(u).$$

Putting $$x_1=-\Big(\frac{r}{u}\frac{du}{dr}\Big)^{-1},
\quad
x_2=\frac{r^2f(u)}{u}\Big(\frac{r}{u}\frac{du}{dr}\Big)^{-2},
\quad x_3=u,$$ we get the autonomous system $$\begin{aligned}
&u\frac{dx_1}{du}=(1-x_1+x_2)x_1, \\
&u\frac{dx_2}{du}=(\nu+1-4x_1+2x_2+G(x_3))x_2, \\
&u\frac{dx_3}{du}=x_3.\end{aligned}$$ Here $(x_1,x_2, x_3)$ is considered as functions of $u >0$, and we have $(x_1,x_2,x_3) \rightarrow (0,0,0)$ as $u \rightarrow +0$.

Then there is a transformation of variables $(x_1,x_2,x_3) \leftrightarrow (\xi_1,\xi_2,\xi_3)$ of the form $$x_j=\xi_j(1+[\xi_1,\xi_2,\xi_3]_1),\quad j=1,2,3,$$ which reduce the system to $$u\frac{d\xi_1}{du}=\xi_1,\quad
u\frac{d\xi_2}{du}=(\nu+1)\xi_2 \quad u\frac{d\xi_3}{du}=\xi_3.$$ Here and hereafter $[X_1, X_2, X_3]_1$ generally stands for various convergent triple power series of the form $\sum_{k_1+k_2+k_3\geq 1}a_{k_1k_2k_3}X_1^{k_1}X_2^{k_2}X_3^{k_3}$. Take a general solution $$\xi_1=C_1u,\quad \xi_2=C_2u^{\nu+1},\quad \xi_3=u.$$ putting $C=1/C_1$, we get the desired expansion, since the integration of $$-\frac{u}{r}\frac{dr}{du}=x^1$$ gives $$\begin{aligned}
\log \frac{R}{r}&=-\log \Big(1-\frac{R-r}{R}\Big)=\frac{R-r}{R}\Big(1+
\Big[\frac{R-r}{R}\Big]_1\Big) \\
&=C_1u(1+[C_1u,C_2u^{\nu+1},u]_1).\end{aligned}$$

[99]{} C. Aerts, J. Christensen-Dalsgaard and D. W. Kurtz, Astroseismology, Springer, Dordrecht-Heidelberg-London-New York, 2010. H. R. Beyer, The spectrum of radial adiabatic stellar oscillations, J. Math. Physics, 36(1995), 4815-4825. G. Birkhoff and G.-C. Rota, Ordinary Differential Equations, 3rd Ed., John Wiley and Sons, New York, 1959. E. A. Coddington and N. Levinson, Theory of Ordinary Differential Equations, McGrawhill, New York-Toronto-London, 1955. T. G. Cowling, The non-radial oscillations of polytropic stars, Monthly Notices Roy. Astronom. Soc. London, 101(1941), 367-375 . J. Cox, Theory of Stellar Pulsation, Princeton University Press, Princeton, 1980. F. -L. Deubner and D. O. Gough, Helioseismology: Oscillations as a diagnostic of the star interior, Ann. Rev. Astron. Astrophys., 22(1984), 594-619. N. Dunford and J. T. Schwartz, Linear Operators, Part II: Spectral Theory, Interscience, New York, 1963. W. Dziembowski, Nonradial oscillations of evolved stars. I. Quasiadiabatic approximation, Acta Astronomica, 21(1971), 289-306. J. Eisenfeld, A complete theorem for an integro-differential operator, J. Math. Anal. and Appl., 26(1969), 357-375. D. O. Gough, Linear adiabatic stellar pulsation, in J. -P. Zahn et al eds, Les Houches Session XLVII 1987, North-Holland, Amsterdam-London-New York-Tokyo, 1993, 1-32. J. D. Jackson, Classical Electrodynamics, Wiley and Sons, New York, 1962. Juhi Jang, Time periodic approximations of the Euler-Poisson system near Lane-Emden stars, Analysis and PDE., 9(2016), 1043-1078. Juhi Jang and T. Makino, Linearized analysis of barotropic perturbations aroud spherically symmetric gaseous stars governed by the Euler-Poisson equations, arXiv:1810.08294. T. Kato, Perturbation Theory for Linear Operators, Springer, Berlin-Heidelberg-New York, 1980. N. R. Lebovitz, On the onset of convective instability, Astrophys. J., 142(1965), 1257-1260. P. Ledoux and Th. Walraven, Variable stars, in S. Flügge ed., Handbuch der Physik, Band LI, Springer, Berlin-Göttingen-Heidelberg, 1958, 353-604. S.-S. Lin, Stability of gaseous stars in spherically symmetric motions, SIAM J. Math. Anal., 28(1997), 539-569. T. Makino, On the existence of positive solutions at infinity for ordinary differential equations of Emden type, Funkcialaj Ekvacioj, 27(1984), 319-329. T. Makino, On spherically symmetric motions of a gaseous star governed by the Euler-Poisson equations, Osaka J. Math., 52(2015), 545-580. M. Reed and B. Simon, Methods of Modern Mathematical Physics, II, Fourier Analysis, Self-Adjointness, AcademicPress, New York, 1975. A. D. Rendall and B. G. Schmidt, Existence and properties of spherically symmetric static fluid bodies with a given equation of state, Classical Quantum Gravity, 8(1991), 985-1000. H. Shibahashi, Modal analysis of stellar nonradial oscillations by an asymptotic method, Publ. Astron. Soc. Japan, 31(1979), 87-104. W. Unno, Y. Osaki, H. Ando and H. Shibahashi, Nonradial Oscillations of Stars, University of Tokyo Press, Tokyo, 1979. W. Unno, Y. Osaki, H. Ando and H. Shibahashi, Nonradial Oscillations of Stars, 2nd Ed., University of Tokyo Press, Tokyo, 1989. K. Yosida, Lectures on Differential and Integral Equations, Interscience, New York, 1960.

[^1]: Professor Emeritus at Yamaguchi University, Japan E-mail: makino@yamaguchi-u.ac.jp
