Fix $\delta$ satisfying Assumption <ref> and define $f^{(3)}(\boldsymbol{\phi})$ and $f(\boldsymbol{\phi}) $ by (<ref>) and (<ref>) respectively. Then, there exists a positive constant $c$ so that

\begin{align*}
f(\boldsymbol{\phi})  \leq c \left( h_{\Omega}(\boldsymbol{\phi}) \right)^2, \quad \| f^{(3)} (\boldsymbol{\phi}) \|_{\mathcal{P}} \leq c \| \boldsymbol{\phi} \|_{\mathcal{P}}^3,
\end{align*}

for all $\boldsymbol{\phi} \in \mathcal{P} $.

Fix $\delta$ satisfying Assumption <ref> and define $f^{(3)}(\boldsymbol{\phi})$ and $f(\boldsymbol{\phi}) $ by (<ref>) and (<ref>) respectively for all $\boldsymbol{\phi} =(\phi,\partial_{t}\phi) \in \mathcal{P}$. By density (Section <ref>), we can assume that $\boldsymbol{\phi} \in C^\infty( 0, \pi/2)$. Firstly, we focus on $f(\boldsymbol{\phi})  \leq c \left( h_{\Omega}(\boldsymbol{\phi}) \right)^2$ and prove that

\begin{align} \label{InequalitySobolevHardyGoalKG}
\left( \int_{0}^{\pi/2} \mathsf{W} \phi^4 d\mu  \right)^{1/4}  \lesssim \left( \int_{0}^{\pi/2}   \phi \mathsf{L}\phi  d\mu  \right)^{1/2},
\end{align}

\begin{align*}
\mathsf{W}(x) &= \cos^{2(\delta-1)}(x), \quad  \mathsf{w}(x)=\sin^2(x) \cos^{2(\delta-1)}(x), \\ d\mu(x)&=\mathsf{w}(x) dx, \quad \mathsf{L} \phi   =-\frac{1}{\mathsf{w}} \partial_{x} \left(\mathsf{w} \partial_{x}\phi \right)+\delta^2 \phi.
\end{align*}

We set $\psi=\cos^{\delta-1}(x) \phi $ so that $\psi(\pi/2)=0$ and compute

\begin{align*}
\mathsf{w} \partial_{x} \phi &=   \mathsf{w} \partial_{x} ( \cos^{-(\delta-1)}(x) \psi) \\
&=  \sin^2(x) \cos^{2\delta-2}(x)  [ (\delta-1)\cos^{-\delta}(x) \sin(x) \psi + \cos^{-\delta+1}(x) \partial_{x}\psi] \\
&=   (\delta-1)   \sin^3(x) \cos^{\delta-2}(x) \psi
+  \sin^2(x) \cos^{\delta-1}(x) \partial_{x}\psi, \\
\partial_{x} \left(  \mathsf{w} \partial_{x} \phi \right)&=
(\delta-1)  \partial_{x} \left( \sin^3(x) \cos^{\delta-2}(x) \psi \right)\\
&+  \cos^{\delta-1}(x) \partial_{x} \left( \sin^2(x)\partial_{x}\psi \right)
+  \sin^2(x)\partial_{x}\psi \partial_{x} \left( \cos^{\delta-1}(x) \right) \\
3 (\delta-1)     \sin^2(x) \cos^{\delta-1}(x) \psi  -
(\delta-1) (\delta-2)     \sin^4(x) \cos^{\delta-3}(x) \psi   \\
& +
(\delta-1)     \sin^3(x) \cos^{\delta-2}(x) \partial_{x}\psi+  \cos^{\delta-1}(x) \partial_{x} \left( \sin^2(x)\partial_{x}\psi \right)\\
&- (\delta-1) \sin^3(x)  \cos^{\delta-2}(x) \partial_{x}\psi   \\
3 (\delta-1)     \sin^2(x) \cos^{\delta-1}(x) \psi  -
(\delta-1) (\delta-2)     \sin^4(x) \cos^{\delta-3}(x) \psi   \\
& +  \cos^{\delta-1}(x) \partial_{x} \left( \sin^2(x)\partial_{x}\psi \right),   \\
\frac{1}{ \mathsf{w}} \partial_{x} \left(  \mathsf{w} \partial_{x} \phi \right)&=
3 (\delta-1)   \cos^{-\delta+1}(x)     \psi  -
(\delta-1) (\delta-2)  \cos^{-\delta-1}(x)
\sin^2(x)   \psi   \\
& +  \cos^{-\delta+1  }(x) \sin^{-2}(x) \partial_{x} \left( \sin^2(x)\partial_{x}\psi \right),   \\
\mathsf{L} \phi &= -3 (\delta-1)   \cos^{-\delta+1}(x)     \psi  +
(\delta-1) (\delta-2)  \cos^{-\delta-1}(x)
\sin^2(x)   \psi \\
& + \delta^2 \cos^{-\delta+1}(x) \psi
-  \cos^{-\delta+1  }(x) \sin^{-2}(x) \partial_{x} \left( \sin^2(x)\partial_{x}\psi \right), \\
\mathsf{w}  \mathsf{L} \phi &= -3 (\delta-1)  \sin^2(x) \cos^{\delta-1}(x)      \psi  +
(\delta-1) (\delta-2) \sin^4(x) \cos^{\delta-3}(x)    \psi \\
& + \delta^2  \sin^2(x) \cos^{\delta-1}(x) \psi
-  \cos^{\delta-1}(x)  \partial_{x} \left( \sin^2(x)\partial_{x}\psi \right)
\end{align*}

and hence

\begin{align*}
\mathsf{w} \phi \mathsf{L} \phi  &= -3 (\delta-1)  \sin^2(x)     \psi  ^2    +
(\delta-1) (\delta-2) \sin^4(x) \cos^{ -2}(x)    \psi^2  \\
& + \delta^2  \sin^2(x)  \psi   ^2
-     \psi  \partial_{x} \left( \sin^2(x)\partial_{x}\psi \right)   .
\end{align*}

Now, using the Dirichlet boundary condition $\psi(\pi/2)=0$, integration by parts yields

\begin{align*}
\int_{0}^{\pi/2} \psi  \partial_{x} \left( \sin^2(x)\partial_{x}\psi \right)       dx=   - \int_{0}^{\pi/2}     \left( \partial_{x}\psi \right) ^2   \sin^2(x)   dx
\end{align*}

and hence we deduce

\begin{align*}
\int_{0}^{\pi/2} \phi \mathsf{L} \phi d \mu &=
\int_{0}^{\pi/2} \mathsf{w} \phi \mathsf{L} \phi  d x  \\
&= -3 (\delta-1) \int_{0}^{\pi/2} \sin^2(x)     \psi  ^2  dx  +
(\delta-1) (\delta-2)\int_{0}^{\pi/2} \sin^4(x) \cos^{ -2}(x)    \psi^2  dx \\
& + \delta^2  \int_{0}^{\pi/2} \sin^2(x) \psi   ^2dx
-    \int_{0}^{\pi/2} \psi  \partial_{x} \left( \sin^2(x)\partial_{x}\psi \right)       dx \\
&=\left[  \delta^2 -3 (\delta-1) \right] \int_{0}^{\pi/2}     \psi  ^2  \sin^2(x) dx  +
(\delta-1) (\delta-2)\int_{0}^{\pi/2}   \tan^2(x)  \psi^2  \sin^2(x)  dx  \\
& + \int_{0}^{\pi/2}    \left( \partial_{x}\psi \right) ^2   \sin^2(x)    dx.
\end{align*}

Notice that

\begin{align*}
\delta^2 -3 (\delta-1) \geq 1, \quad (\delta-1) (\delta-2) \geq 0,
\end{align*}

for all integers $\delta \geq 2$, and consequently we arrive at

\begin{align*}
\int_{0}^{\pi/2} \phi  \mathsf{L} \phi d \mu  &  \geq   \int_{0}^{\pi/2}     \psi  ^2  \sin^2(x) dx   + \int_{0}^{\pi/2}    \left( \partial_{x}\psi \right) ^2   \sin^2(x)    dx  .
\end{align*}

We use once again the Dirichlet boundary condition $\psi(\pi/2)=0$ to extend $\psi:(0,\pi/2) \rightarrow \mathbb{R}$ to a smooth function $\psi:(0,\pi)\rightarrow \mathbb{R}$ imposing the reflection symmetry $\psi(x)=-\psi(\pi-x)$, for all $x \in (0,\pi)$. Furthermore, let $f: \mathbb{S}^3= (0,\pi) \times (0,\pi) \times (0,2\pi)\rightarrow \mathbb{R}$ defined by $f(x,\vartheta,\varphi)=\psi(x)$ and note that $d\text{vol}_{\mathbb{S}^3}(x,\vartheta,\varphi) = \sin^2 (x) \sin(\vartheta) dx d\vartheta d\varphi$ is the volume form on the entire $3$-sphere. On the one hand, we have that

\begin{align} \label{Hardy1}
\left( 	\int_{0}^{\pi/2} \phi  \mathsf{L} \phi d \mu \right)^{1/2} &  \geq \left(  \int_{0}^{\pi/2}   \left[  \psi  ^2 +\left( \partial_{x}\psi \right) ^2 \right] \sin^2(x) dx \right)^{1/2}   \simeq \left(  \int_{0}^{\pi}   \left[  \psi  ^2 +\left( \partial_{x}\psi \right) ^2 \right] \sin^2(x) dx \right)^{1/2}  \nonumber \\
& \simeq \left(  \int_{\mathbb{S}^3}  \left[   f ^2 +   \left( \nabla_{\mathbb{S}^3}f\right) ^2 \right] d\text{vol}_{\mathbb{S}^3}  \right)^{1/2}   =  \| f\|_{H^1(\mathbb{S}^3)}.
\end{align}

On the other hand, we also have that

\begin{align} \label{Hardy2}
\left(	\int_{0}^{\pi/2}  \mathsf{W} \phi^4 d\mu \right)^{1/4} & =
\left( \int_{0}^{\pi/2}  \psi^4 \sin^2(x) dx \right)^{1/4}  \simeq
\left( 	\int_{0}^{\pi}  \psi^4 \sin^2(x) dx   \right)^{1/4} \nonumber \\
& \simeq \left( \int_{\mathbb{S}^3}  f^4   d\text{vol}_{\mathbb{S}^3} \right)^{1/4}    \simeq  \| f\|_{L^4(\mathbb{S}^3)} .
\end{align}

In view of (<ref>) and (<ref>), in order to prove (<ref>), it suffices to show that there exists a positive constant so that

\begin{align*}
\| f\|_{L^4(\mathbb{S}^3)} \lesssim  \| f\|_{H^1(\mathbb{S}^3)}.
\end{align*}

The latter is a direct consequence of the standard Sobolev inequality on the 3-sphere (Theorem 5.1, page 121 in [31]),

\begin{align*}
\| f\|_{L^6(\mathbb{S}^3)} \lesssim  \| f\|_{H^1(\mathbb{S}^3)},
\end{align*}

together with the embedding $L^6(\mathbb{S}^3) \hookrightarrow L^4(\mathbb{S}^3) $ since $ \mathbb{S}^3  $ is compact. Finally, we focus on the second estimate and prove that $\| f^{(3)} (\boldsymbol{\phi}) \|_{\mathcal{P}} \lesssim \| \boldsymbol{\phi} \|_{\mathcal{P}}^3$. To this end, we proceed as above and obtain

\begin{align*}
\| f^{(3)} (\boldsymbol{\phi}) \|^{2}_{\mathcal{P}} &=
\|   \mathsf{W} \phi^3  \|^{2}_{L^2 (  (0,\pi/2);d\mu  )}
=  \int_{0}^{\pi/2}  \mathsf{W}^2(x) \phi^6(x)   d\mu(x)   =  \int_{0}^{\pi/2} 	\psi^6 (x) \sin^2(x)   dx \\
\simeq 	\int_{0}^{\pi}  \psi^6 \sin^2(x) dx
\simeq  \int_{\mathbb{S}^3}  f^6   d\text{vol}_{\mathbb{S}^3}     \simeq  \| f\|_{L^6(\mathbb{S}^3)}^6 \lesssim  \| f\|_{H^1(\mathbb{S}^3)}^6 \\
& \lesssim \|\phi\|_{H^1 ( (0,\pi/2);d\mu )}^6 \leq \|  \boldsymbol{\phi}  \| _{\mathcal{P}}^6 ,
\end{align*}

that completes the proof.

Next, we compute the harmonic energy in the Fourier space.

Let $\boldsymbol{\phi}=(\phi_1,\phi_2) \in \mathcal{P}$ be expanded in terms of the eigenfunctions $\{e_m:m\geq 0 \}$ to the linearized operator $L$,

\begin{align*}
\phi_{ \alpha} (t)=\sum_{m=0}^{\infty} \phi_{ \alpha}^{(m)} (t)e_m, \quad
\alpha \in \{1,2\}.
\end{align*}

\begin{align*}
h_{\Omega}(\boldsymbol{\phi}) &= \sum_{m=0}^{\infty}
( \phi_2^{(m)} )^2  +\sum_{m=0}^{\infty} \omega_{m}^2  ( \phi_1^{(m)} )^2  .
\end{align*}

Let $\boldsymbol{\phi}=(\phi_1,\phi_2) \in \mathcal{P}$ and expand its coefficients   in terms of the eigenfunctions $\{e_m:m\geq 0 \}$ to the linearized operator $L$,

\begin{align*}
\phi_{ \alpha} (t)=\sum_{m=0}^{\infty} \phi_{ \alpha}^{(m)} (t)e_m, \quad
\alpha \in \{1,2\}.
\end{align*}

Then, we use the fact that the set of eigenfunctions $\{e_m:m\geq 0\}$ is orthonormal with respect to the inner product in $L^2((0,\pi/2);d\mu)$, to compute that

\begin{align*}
h_{\Omega}(\boldsymbol{\phi}) &= \frac{1}{2} \Omega(A\boldsymbol{\phi},\boldsymbol{\phi})
\phi_2 ^2 + \left(L^{1/2}\phi_{1}\right)^2 \right)  d\mu  \\
& =\sum_{i,j=0}^{\infty}  \int_{0}^{\pi/2}\left(
\phi_2^{(i)}\phi_2^{(j)} e_i (x)e_j(x) +  \phi_1^{(i)}\phi_1^{(j)} L^{1/2}e_i (x)L^{1/2}e_j(x) \right)  d\mu (x)  \\
& =\sum_{i,j=0}^{\infty}  \left(
\phi_2^{(i)}\phi_2^{(j)}   + \omega_{i}\omega_{j} \phi_1^{(i)}\phi_1^{(j)}   \right)\int_{0}^{\pi/2}e_i (x) e_j(x)  d\mu (x)  \\
& =\sum_{i,j=0}^{\infty}  \left(
\phi_2^{(i)}\phi_2^{(j)}   + \omega_{i}\omega_{j} \phi_1^{(i)}\phi_1^{(j)}   \right) \delta_{ij}  = \sum_{m=0}^{\infty}
( \phi_2^{(m)} )^2  +\sum_{m=0}^{\infty} \omega_{m}^2  ( \phi_1^{(m)})^2  ,
\end{align*}

that completes the proof.

Next, we compute the linear flow in the Fourier space. To do so, we denote by $ e^{At}$ the unitary group generated by $A$.

Let $\boldsymbol{\zeta}=(\zeta_1,\zeta_2) \in \mathcal{P}$ be expanded in terms of the eigenfunctions $	\{ e_m: m\geq 0 \}$ to the linearized operator $L$,

\begin{align*}
\zeta_\alpha =\sum_{m=0}^{\infty} \zeta_\alpha^{(m)} e_m,\quad  \alpha \in \{1,2\}.
\end{align*}

\begin{align*}
\boldsymbol{\Phi}^{t}(\boldsymbol{\zeta}) := e^{At}\boldsymbol{\zeta}=
\begin{pmatrix}
\sum_{m=0}^{\infty} \left(\zeta_{1}^{(m)}  \cos(\omega_{m}t) + \frac{\zeta_2^{(m)}}{\omega_{m}} \sin (\omega_{m}t) \right) e_m \\
\sum_{m=0}^{\infty} \left(-\omega_{m}\zeta_{1}^{(m)}  \sin(\omega_{m}t) +  \zeta_2^{(m)}  \cos (\omega_{m}t) \right) e_m
\end{pmatrix}.
\end{align*}

Let $\boldsymbol{\zeta}=(\zeta_1,\zeta_2) \in \mathcal{P}$ and recall that $\boldsymbol{\phi}_{\text{linear}}(t)=\boldsymbol{\Phi}^{t}(\boldsymbol{\zeta})=e^{tA}\boldsymbol{\zeta}$  is the solution to the initial value problem (<ref>). For $\boldsymbol{\phi}_{\text{linear}}(t)=(\phi_{\text{linear},1}(t),\phi_{\text{linear},2}(t)) \in \mathcal{P}$, we compute

\begin{align*}
\begin{pmatrix}
\partial_{t}\phi_{\text{linear},1}(t) \\
\partial_{t}\phi_{\text{linear},2}(t)
\end{pmatrix}
\begin{pmatrix}
0 & 1 \\
-L & 0
\end{pmatrix}
\begin{pmatrix}
\phi_{\text{linear},1}(t) \\
\phi_{\text{linear},2}(t)
\end{pmatrix}=
\begin{pmatrix}
\phi_{\text{linear},2}(t) \\
-L \phi_{\text{linear},1}(t)
\end{pmatrix}
\end{align*}

that leads to

\begin{align}\label{AuxiliarlySecondOrderODE}
\partial_{t}^2 \phi_{\text{linear},1}(t) +L  \phi_{\text{linear},1}(t) =0.
\end{align}

Expanding $\phi_{\text{linear},1}(t)$, $\phi_{\text{linear},2}(t)$, $\zeta_1$ and $\zeta_2$ in terms of the eigenfunctions of $L$,

\begin{align*}
\phi_{\text{linear},\alpha} (t)=\sum_{m=0}^{\infty} \phi_{\text{linear},\alpha}^{(m)} (t)e_m, \quad
\zeta_\alpha =\sum_{m=0}^{\infty} \zeta_\alpha^{(m)} e_m,\quad  \alpha \in \{1,2\},
\end{align*}

we infer

\begin{align*}
\phi_{\text{linear},1}^{(m)}(0) &= \left( \phi_{\text{linear},1}(0) | e_m \right)  = \left( \zeta_{1} | e_m \right)= \zeta_1^{(m)}, \\
\partial_{t}  \phi_{\text{linear,1}}^{(m)}(0)&= \left(\partial_{t}\phi_{\text{linear,1}}(0) | e_m \right)  = \left( (A\zeta)_{1} | e_m \right)=\left( \zeta_{2} | e_m \right)= \zeta_2^{(m)},
\end{align*}

for all integers $m\geq0$, and hence (<ref>) yields

\begin{align*}
\frac{d^2}{dt^2} \ddot{\phi}_{\text{linear,1}}^{(m)}(t)+\omega_{m}^2 \phi_{\text{linear},1}^{(m)}(t)=0
\Longrightarrow   \phi_{\text{linear},1}^{(m)}(t)= \zeta_{1}^{(m)}  \cos(\omega_{m}t) + \frac{\zeta_2^{(m)}}{\omega_{m}} \sin (\omega_{m}t),
\end{align*}

for all integers $m\geq0$. Consequently, we obtain

\begin{align*}
\boldsymbol{\Phi}^{t}(\boldsymbol{\zeta}) = e^{At}\boldsymbol{\zeta}=
\begin{pmatrix}
\sum_{m=0}^{\infty} \left(\zeta_{1}^{(m)}  \cos(\omega_{m}t) + \frac{\zeta_2^{(m)}}{\omega_{m}} \sin (\omega_{m}t) \right) e_m \\
\sum_{m=0}^{\infty} \left(-\omega_{m}\zeta_{1}^{(m)}  \sin(\omega_{m}t) +  \zeta_2^{(m)}  \cos (\omega_{m}t) \right) e_m
\end{pmatrix},
\end{align*}

that completes the proof.

Notice that $\boldsymbol{\Phi}^{t} (\boldsymbol{\zeta})$ is $2\pi$-periodic in time. The leading order term $f(\boldsymbol{\phi})$ in the Hamiltonian is given by (<ref>) and we define its time-average over the linear flow  by

\begin{align*}
\langle f  \rangle (\boldsymbol{\zeta})= \frac{1}{2\pi} \int_{0}^{2\pi} f (\boldsymbol{\Phi}^{t}(\boldsymbol{\zeta}) )	dt,
\end{align*}

where $\boldsymbol{\Phi}^{t}(\boldsymbol{\zeta}) $ stands for the linear flow, given by Lemma <ref>, associated with $\boldsymbol{\zeta}$. Next, we compute this time-averaging operator in the Fourier space.

Let $\boldsymbol{\zeta}=(\zeta_1,\zeta_2) \in \mathcal{P}$ be expanded in terms of the eigenfunctions $\{e_m:m\geq 0 \}$ to the linearized operator $L$,

\begin{align*}
\zeta_{ \alpha} =\sum_{m=0}^{\infty} \zeta_{ \alpha}^{(m)} e_m,
\quad
\alpha \in \{1,2\}.
\end{align*}

\begin{align*}
\langle f \rangle (\boldsymbol{\zeta})
&=\frac{3}{32}\sum_{i,j,k,m=0}^{\infty} C_{ijkm}
P_{ijkm}(\boldsymbol{\zeta})  \mathds{1}(m=i+j-k \geq 0) ,
\end{align*}

\begin{align*}
P_{ijkm}(\boldsymbol{\zeta})  &= -\frac{\zeta_2^{  ( i )} \zeta_2^{  ( j )} \zeta_1^{  ( k )} \zeta_1^{  ( m )}}{\omega_{ i } \omega_{ j }}+\frac{\zeta_2^{  ( i )} \zeta_1^{  ( j )} \zeta_2^{  ( k )} \zeta_1^{  ( m )}}{\omega_{ i } \omega_{ k }}+\frac{\zeta_1^{  ( i )} \zeta_2^{  ( j )} \zeta_2^{  ( k )} \zeta_1^{  ( m )}}{\omega_{ j } \omega_{ k }} +\frac{\zeta_2^{  ( i )} \zeta_1^{  ( j )} \zeta_1^{  ( k )} \zeta_2^{  ( m )}}{\omega_{ i } \omega_{ m }} \\
&+\frac{\zeta_1^{  ( i )} \zeta_2^{  ( j )} \zeta_1^{  ( k )} \zeta_2^{  ( m )}}{\omega_{ j } \omega_{ m }}-\frac{\zeta_1^{  ( i )} \zeta_1^{  ( j )} \zeta_2^{  ( k )} \zeta_2^{  ( m )}}{\omega_{ k } \omega_{ m }}+\zeta_1^{  ( i )} \zeta_1^{  ( j )} \zeta_1^{  ( k )} \zeta_1^{  ( m )}+\frac{\zeta_2^{  ( i )} \zeta_2^{  ( j )} \zeta_2^{  ( k )} \zeta_2^{  ( m )}}{\omega_{ i } \omega_{ j } \omega_{ k } \omega_{ m }}.
\end{align*}

Let $\boldsymbol{\phi} =(\phi_{ 1},\phi_{ 2}) \in \mathcal{P}$ and expand its component in terms of the eigenfunctions to the lienarized operator $L$, namely

\begin{align*}
\phi_{ \alpha} (t)=\sum_{m=0}^{\infty} \phi_\alpha^{(m)} (t)e_m, \quad  \alpha \in \{1,2\}.
\end{align*}

Then, we also expand

\begin{align}\label{fphi}
f(\boldsymbol{\phi}) &=  \frac{1}{4} \int_{0}^{\pi/2} W \phi_1^4 d\mu \nonumber \\
&=  \frac{1}{4}\sum_{i,j,k,m=0}^{\infty} \phi_{1}^{(i)} (t)\phi_{1}^{(j)} (t)\phi_{1}^{(k)} (t)\phi_{1}^{(m)} (t) \int_{0}^{\pi/2} W(x)e_i(x)e_j(x)e_k(x)e_m(x)d\mu (x)\nonumber \\
&=  \frac{1}{4}\sum_{i,j,k,m=0}^{\infty} C_{ijkm}   \phi_{1}^{(i)} (t)\phi_{1}^{(j)} (t)\phi_{1}^{(k)} (t)\phi_{1}^{(m)} (t) ,
\end{align}

where we also used the definition of the Fourier coefficient from (<ref>). Let $\boldsymbol{\zeta}=(\zeta_1,\zeta_2) \in \mathcal{P}$ and consider the linear flow $\boldsymbol{\phi}_{\text{linear}}(t)=\boldsymbol{\Phi}^{t}(\boldsymbol{\zeta})=(\phi_{\text{linear},1}(t),\phi_{\text{linear},2}(t)) \in \mathcal{P}$ associated with $\boldsymbol{\zeta}$ given by Lemma <ref>. In particular, we have that

\begin{align}\label{philinear1}
\phi_{\text{linear},1}^{(m)}(t) = \zeta_{1}^{(m)}  \cos(\omega_{m}t) + \frac{\zeta_2^{(m)}}{\omega_{m}} \sin (\omega_{m}t) ,
\end{align}

for all integers $m \geq 0$. To compute the time average

\begin{align*}
\langle f  \rangle (\boldsymbol{\zeta})= \frac{1}{2\pi} \int_{0}^{2\pi} f (\boldsymbol{\Phi}^{t}(\boldsymbol{\zeta}))	dt,
\end{align*}

we firstly plug (<ref>) into (<ref>) and use trigonometric identities to expand all products $\cos(\omega_{l_1}t)\sin(\omega_{l_2}t)$, $\sin(\omega_{l_1}t)\sin(\omega_{l_2}t)$ and $\cos(\omega_{l_1}t)\cos(\omega_{l_2}t)$ with $l_1,l_2\in \{i,j,k,m\}$ in terms of $\cos((\omega_{i}\pm \omega_{j}\pm \omega_{k}\pm \omega_{m})t)$ and $\sin((\omega_{i}\pm \omega_{j}\pm \omega_{k}\pm \omega_{m})t)$. We obtain a long formula of the form

\begin{align}\label{ComplicatedResult}
f(\boldsymbol{\Phi}^{t}(\boldsymbol{\zeta}))	=\sum_{T\in \{\sin,\cos\}} \sum_{(i,j,k,m)\in I_0 \cup I_1 \cup I_2}Z_{ijkm} C_{ijkm} T ( (\omega_{i}\pm \omega_{j}\pm \omega_{k}\pm \omega_{m})t ),
\end{align}

where $Z_{ijkm}$ depends on $\zeta_1^{(l)}$ and $\zeta_2^{(l)}$ for all $l\in\{i,j,k,m \}$ and

\begin{align*}
I_0 &= \left\{(i,j,k,m):~ \text{no minus signs in } \omega_{i}\pm \omega_{j}\pm \omega_{k}\pm \omega_{m} \right\}, \\
I_1 &= \left\{(i,j,k,m):~ \text{only one minus sign in } \omega_{i}\pm \omega_{j}\pm \omega_{k}\pm \omega_{m} \right\}, \\
I_2 &= \left\{(i,j,k,m):~ \text{only two minus signs in } \omega_{i}\pm \omega_{j}\pm \omega_{k}\pm \omega_{m} \right\}.
\end{align*}

According to Lemma <ref>, we have that $ C_{`ijkm} =0$, for all $(i,j,k,m) \in I_1$, hence we only need to compute the sum with respect to $(i,j,k,m) \in I_0 \cup I_2$. Then, we apply the integral with respect to $t\in(0,2\pi)$ to each of the summands in (<ref>) and use the fact that

\begin{align*}
\frac{1}{2\pi }\int_{0}^{2\pi } \cos(\omega t)dt = \mathds{1}(\omega=0), \quad
\frac{1}{2\pi }\int_{0}^{2\pi } \sin(\omega t)dt=0
\end{align*}

to compute each integral obtaining

\begin{align*}
\langle f  \rangle (\boldsymbol{\zeta})=	\frac{1}{2\pi }\int_{0}^{2\pi }f(\boldsymbol{\Phi}^{t}(\boldsymbol{\zeta}))dt=\sum_{(i,j,k,m)\in  I_1 \cup I_2}Z_{ijkm} C_{ijkm} \mathds{1} ( \omega_{i}\pm \omega_{j}\pm \omega_{k}\pm \omega_{m} =0 ).
\end{align*}

In conclusion, we infer

\begin{align*}
\langle f  \rangle (\boldsymbol{\zeta})
&=\frac{1}{32}\sum_{i,j,k,m=0}^{\infty} C_{ijkm}  \Big[
P_{ijkm}(\boldsymbol{\zeta})  \mathds{1}(m=i+j-k \geq 0)\nonumber \\
&+Q_{ijkm}(\boldsymbol{\zeta})  \mathds{1}(m=i-j+k \geq 0)
+R_{ijkm}(\boldsymbol{\zeta}) \mathds{1}(m=-i+j+k \geq 0)
\Big],
\end{align*}

\begin{align*}
P_{ijkm}(\boldsymbol{\zeta})  &= -\frac{\zeta_2^{  ( i )} \zeta_2^{  ( j )} \zeta_1^{  ( k )} \zeta_1^{  ( m )}}{\omega_{ i } \omega_{ j }}+\frac{\zeta_2^{  ( i )} \zeta_1^{  ( j )} \zeta_2^{  ( k )} \zeta_1^{  ( m )}}{\omega_{ i } \omega_{ k }}+\frac{\zeta_1^{  ( i )} \zeta_2^{  ( j )} \zeta_2^{  ( k )} \zeta_1^{  ( m )}}{\omega_{ j } \omega_{ k }} +\frac{\zeta_2^{  ( i )} \zeta_1^{  ( j )} \zeta_1^{  ( k )} \zeta_2^{  ( m )}}{\omega_{ i } \omega_{ m }} \\
&+\frac{\zeta_1^{  ( i )} \zeta_2^{  ( j )} \zeta_1^{  ( k )} \zeta_2^{  ( m )}}{\omega_{ j } \omega_{ m }}-\frac{\zeta_1^{  ( i )} \zeta_1^{  ( j )} \zeta_2^{  ( k )} \zeta_2^{  ( m )}}{\omega_{ k } \omega_{ m }}+\zeta_1^{  ( i )} \zeta_1^{  ( j )} \zeta_1^{  ( k )} \zeta_1^{  ( m )}+\frac{\zeta_2^{  ( i )} \zeta_2^{  ( j )} \zeta_2^{  ( k )} \zeta_2^{  ( m )}}{\omega_{ i } \omega_{ j } \omega_{ k } \omega_{ m }}, \\
Q_{ijkm}(\boldsymbol{\zeta})  &=  \frac{\zeta_2^{  ( i )} \zeta_2^{  ( j )} \zeta_1^{  ( k )} \zeta_1^{  ( m )}}{\omega_{ i } \omega_{ j }}-\frac{\zeta_2^{  ( i )} \zeta_1^{  ( j )} \zeta_2^{  ( k )} \zeta_1^{  ( m )}}{\omega_{ i } \omega_{ k }}+\frac{\zeta_1^{  ( i )} \zeta_2^{  ( j )} \zeta_2^{  ( k )} \zeta_1^{  ( m )}}{\omega_{ j } \omega_{ k }}+\frac{\zeta_2^{  ( i )} \zeta_1^{  ( j )} \zeta_1^{  ( k )} \zeta_2^{  ( m )}}{\omega_{ i } \omega_{ m }} \\
&+\zeta_1^{  ( i )} \zeta_1^{  ( j )} \zeta_1^{  ( k )} \zeta_1^{  ( m )} -\frac{\zeta_1^{  ( i )} \zeta_2^{  ( j )} \zeta_1^{  ( k )} \zeta_2^{  ( m )}}{\omega_{ j } \omega_{ m }}+\frac{\zeta_1^{  ( i )} \zeta_1^{  ( j )} \zeta_2^{  ( k )} \zeta_2^{  ( m )}}{\omega_{ k } \omega_{ m }}+\frac{\zeta_2^{  ( i )} \zeta_2^{  ( j )} \zeta_2^{  ( k )} \zeta_2^{  ( m )}}{\omega_{ i } \omega_{ j } \omega_{ k } \omega_{ m }} , \\
R_{ijkm}(\boldsymbol{\zeta})  &=  \frac{\zeta_2^{  ( i )} \zeta_2^{  ( j )} \zeta_1^{  ( k )} \zeta_1^{  ( m )}}{\omega_{ i } \omega_{ j }}+\frac{\zeta_2^{  ( i )} \zeta_1^{  ( j )} \zeta_2^{  ( k )} \zeta_1^{  ( m )}}{\omega_{ i } \omega_{ k }}-\frac{\zeta_1^{  ( i )} \zeta_2^{  ( j )} \zeta_2^{  ( k )} \zeta_1^{  ( m )}}{\omega_{ j } \omega_{ k }}-\frac{\zeta_2^{  ( i )} \zeta_1^{  ( j )} \zeta_1^{  ( k )} \zeta_2^{  ( m )}}{\omega_{ i } \omega_{ m }} \\
&+\zeta_1^{  ( i )} \zeta_1^{  ( j )} \zeta_1^{  ( k )} \zeta_1^{  ( m )} +\frac{\zeta_1^{  ( i )} \zeta_1^{  ( j )} \zeta_2^{  ( k )} \zeta_2^{  ( m )}}{\omega_{ k } \omega_{ m }}+\frac{\zeta_1^{  ( i )} \zeta_2^{  ( j )} \zeta_1^{  ( k )} \zeta_2^{  ( m )}}{\omega_{ j } \omega_{ m }}+\frac{\zeta_2^{  ( i )} \zeta_2^{  ( j )} \zeta_2^{  ( k )} \zeta_2^{  ( m )}}{\omega_{ i } \omega_{ j } \omega_{ k } \omega_{ m }} .
\end{align*}

Finally, notice that from the formulas above one can immediately verify that

\begin{align*}
Q_{ijkm}(\boldsymbol{\zeta})=P_{ikjm}(\boldsymbol{\zeta}), \quad
\end{align*}

for all integers $i,j,k,m \geq 0$. Moreover, we also have the symmetries of the Fourier coefficients $C_{ijkm}=C_{ikjm}=C_{kjim}$, for all integers $i,j,k,m \geq 0$, from (<ref>) and (<ref>). Putting these together, we exchange $k$ with $j$ in the second sum,

\begin{align*}
\sum_{i,j,k,m=0}^{\infty} C_{ijkm}Q_{ijkm}(\boldsymbol{\zeta})  \mathds{1}(m=i-j+k \geq 0)  &=
\sum_{i,j,k,m=0}^{\infty} C_{ikjm}Q_{ikjm}(\boldsymbol{\zeta})  \mathds{1}(m=i-k+j \geq 0) \\
\sum_{i,j,k,m=0}^{\infty} C_{ijkm}P_{ijkm}(\boldsymbol{\zeta})  \mathds{1}(m=i-k+j \geq 0),
\end{align*}

and $k$ with $i$ in the third sum,

\begin{align*}
\sum_{i,j,k,m=0}^{\infty} C_{ijkm}R_{ijkm}(\boldsymbol{\zeta}) \mathds{1}(m=-i+j+k \geq 0) &=
\sum_{i,j,k,m=0}^{\infty} C_{kjim}R_{kjim}(\boldsymbol{\zeta}) \mathds{1}(m=-k+j+i \geq 0) \\
\sum_{i,j,k,m=0}^{\infty} C_{ijkm}P_{ijkm}(\boldsymbol{\zeta}) \mathds{1}(m=-k+j+i \geq 0),
\end{align*}

to infer

\begin{align*}
\langle f  \rangle (\boldsymbol{\zeta})
&=\frac{3}{32}\sum_{i,j,k,m=0}^{\infty} C_{ijkm}
P_{ijkm}(\boldsymbol{\zeta})  \mathds{1}(m=i+j-k \geq 0) ,
\end{align*}

that completes the proof.

Next, we verify Condition (<ref>) in Theorem <ref>. To do so, we use the method of Lagrange multipliers and show that the rescaled first 1-mode defined in (<ref>) is a critical point for the  function  $\langle f\rangle $ restricted on a hypersurface on constant harmonic energy.

Let $\boldsymbol{\xi}=(\xi_1,\xi_2)$ be the rescaled first 1-mode  given by (<ref>) and let

\begin{align*}
\Sigma= \left \{\boldsymbol{\zeta}\in \mathcal{P}:~  h_{\Omega}(\boldsymbol{\zeta})=\kappa_{0}^2 \omega_{0}^2 \right\}.
\end{align*}

be the hypersurface of constant harmonic energy $\kappa_{0}^2 \omega_{0}^2$. Then, $\boldsymbol{\xi} \in \Sigma$ is a critical point for the function $\langle f\rangle |_{\Sigma}$.

Pick any $\boldsymbol{\zeta}=(\zeta_1,\zeta_2) \in \mathcal{P}$, expand its coefficients in terms of the eigenfunctions $\{e_m:m\geq 0 \}$ to the linearized operator $L$,

\begin{align*}
\zeta_\alpha =\sum_{m=0}^{\infty} \zeta_\alpha^{(m)} e_m,\quad  \alpha \in \{1,2\},
\end{align*}

and consider the functions

\begin{align*}
h_{\Omega}(\boldsymbol{\zeta}) &=
h_{\Omega} (\zeta_1^{(0)},\zeta_1^{(1)},\zeta_1^{(2)},\dots, \zeta_2^{(0)}, \zeta_2^{(1)}, \zeta_2^{(2)},\dots), \\
\langle f  \rangle (\boldsymbol{\zeta})&=
\langle f  \rangle (\zeta_1,\zeta_2)=
\langle f  \rangle (\zeta_1^{(0)},\zeta_1^{(1)},\zeta_1^{(2)},\dots, \zeta_2^{(0)}, \zeta_2^{(1)}, \zeta_2^{(2)},\dots).
\end{align*}

These are given in the Fourier space by Lemmata <ref> and <ref> respectively. According to the method of Lagrange multipliers, the point $\boldsymbol{\zeta}$ is a critical point for $\langle f\rangle |_{\Sigma}$ if it satisfies the equations

\begin{align}\label{LagrangeMultipliers}
\nabla_{\boldsymbol{\zeta}} \left[ \lambda  h_{\Omega}(\boldsymbol{\zeta}) +   \langle f  \rangle (\boldsymbol{\zeta}) \right] =0 \Longleftrightarrow
\begin{dcases}
\lambda \frac{\partial}{\partial \zeta_1^{(l)}} h_{\Omega}(\boldsymbol{\zeta})+
\frac{\partial}{\partial \zeta_1^{(l)}}  \langle f  \rangle (\boldsymbol{\zeta})=0, \\
\lambda \frac{\partial}{\partial \zeta_2^{(l)}}  h_{\Omega}(\boldsymbol{\zeta})+
\frac{\partial}{\partial \zeta_2^{(l)}}  \langle f  \rangle (\boldsymbol{\zeta})=0,
\end{dcases}
\end{align}

for all integers $l \geq 0$ and some $\lambda \in \mathbb{R}$.  On the one hand, we use Lemma <ref> to compute

\begin{align*}
\frac{\partial}{\partial \zeta_1^{(l)}} h_{\Omega}(\boldsymbol{\zeta}) = 2\omega_{l}^2  \zeta_1^{(l)}, \quad
\frac{\partial}{\partial \zeta_2^{(l)}} h_{\Omega}(\boldsymbol{\zeta}) = 2   \zeta_2^{(l)},
\end{align*}

for all integers $l \geq 0$. On the other hand, we use Lemma <ref> together with the differentiation rule

\begin{align*}
\frac{\partial}{\partial \zeta_{\alpha}^{(l)}} (\zeta_{\alpha}^{(i)} \zeta_{\alpha}^{(j)} ) =
\zeta_{\alpha}^{(i)}\frac{\partial}{\partial \zeta_{\alpha}^{(l)}} \zeta_{\alpha}^{(j)}+
\zeta_{\alpha}^{(j)}\frac{\partial}{\partial \zeta_{\alpha}^{(l)}}\zeta_{\alpha}^{(i)} =
\zeta_{\alpha}^{(i)} \mathds{1}(l=j) +
\zeta_{\alpha}^{(j)}\mathds{1}(l=i),
\end{align*}

for all integers $\alpha \in \{1,2\}$ and $i,j,l \geq 0$, to obtain that

\begin{align*}
\frac{\partial}{\partial \zeta_1^{(l)}} 	P_{ijkm}(\boldsymbol{\zeta})
\end{align*}

is given by

\begin{align*}
% \frac{\partial}{\partial \zeta_1^{(l)}} 	P_{ijkm}(\boldsymbol{\zeta})  &=
% -\frac{\zeta_2^{  ( i )} \zeta_2^{  ( j )}  }{\omega_{ i } \omega_{ j }}
%  \frac{\partial}{\partial \zeta_1^{(l)}} \left[ \zeta_1^{  ( k )} \zeta_1^{  ( m )}\right]
% %
% %
% +\frac{\zeta_2^{  ( i )}  \zeta_2^{  ( k )}  }{\omega_{ i } \omega_{ k }}
%  \frac{\partial}{\partial \zeta_1^{(l)}} \left[\zeta_1^{  ( j )}\zeta_1^{  ( m )} \right]
% %
% %
% +\frac{ \zeta_2^{  ( j )} \zeta_2^{  ( k )} }{\omega_{ j } \omega_{ k }}
%  \frac{\partial}{\partial \zeta_1^{(l)}} \left[ \zeta_1^{  ( i )}\zeta_1^{  ( m )}\right]\\
% %
% %
% & +\frac{\zeta_2^{  ( i )}   \zeta_2^{  ( m )}}{\omega_{ i } \omega_{ m }}
%  \frac{\partial}{\partial \zeta_1^{(l)}} \left[\zeta_1^{  ( j )}\zeta_1^{  ( k )} \right]
%+\frac{ \zeta_2^{  ( j )}  \zeta_2^{  ( m )}}{\omega_{ j } \omega_{ m }}
% \frac{\partial}{\partial \zeta_1^{(l)}} \left[\zeta_1^{  ( i )} \zeta_1^{  ( k )} \right]
%-\frac{ \zeta_2^{  ( k )} \zeta_2^{  ( m )}}{\omega_{ k } \omega_{ m }}
% \frac{\partial}{\partial \zeta_1^{(l)}} \left[\zeta_1^{  ( i )} \zeta_1^{  ( j )} \right] \\
%& +
% \frac{\partial}{\partial \zeta_1^{(l)}} \left[ \zeta_1^{  ( i )} \zeta_1^{  ( j )} \zeta_1^{  ( k )} \zeta_1^{  ( m )}\right]  \\
& -\frac{\zeta_2^{  ( i )} \zeta_2^{  ( j )}  }{\omega_{ i } \omega_{ j }}
\zeta_1^{  ( m )}\mathds{1}(l=k)
-\frac{\zeta_2^{  ( i )} \zeta_2^{  ( j )}  }{\omega_{ i } \omega_{ j }} \zeta_1^{  ( k )} \mathds{1}(l=m)
+\frac{\zeta_2^{  ( i )}  \zeta_2^{  ( k )}  }{\omega_{ i } \omega_{ k }}
\zeta_1^{  ( m )}\mathds{1}(l=j)
+\frac{\zeta_2^{  ( i )}  \zeta_2^{  ( k )}  }{\omega_{ i } \omega_{ k }}\zeta_1^{  ( j )} \mathds{1}(l=m)   \\
& +\frac{ \zeta_2^{  ( j )} \zeta_2^{  ( k )} }{\omega_{ j } \omega_{ k }}
\zeta_1^{  ( m )}  \mathds{1}(l=i)
+\frac{ \zeta_2^{  ( j )} \zeta_2^{  ( k )} }{\omega_{ j } \omega_{ k }}   \zeta_1^{  ( i )} \mathds{1}(l=m)
+\frac{\zeta_2^{  ( i )}   \zeta_2^{  ( m )}}{\omega_{ i } \omega_{ m }}
\zeta_1^{  ( k )} \mathds{1}(l=j)
+\frac{\zeta_2^{  ( i )}   \zeta_2^{  ( m )}}{\omega_{ i } \omega_{ m }}   \zeta_1^{  ( j )} \mathds{1}(l=k)   \\
& +\frac{ \zeta_2^{  ( j )}  \zeta_2^{  ( m )}}{\omega_{ j } \omega_{ m }}
\zeta_1^{  ( k )} \mathds{1}(l=i)
+\frac{ \zeta_2^{  ( j )}  \zeta_2^{  ( m )}}{\omega_{ j } \omega_{ m }} \zeta_1^{  ( i )}  \mathds{1}(l=k)
-\frac{ \zeta_2^{  ( k )} \zeta_2^{  ( m )}}{\omega_{ k } \omega_{ m }}
\zeta_1^{  ( j )}  \mathds{1}(l=i)
-\frac{ \zeta_2^{  ( k )} \zeta_2^{  ( m )}}{\omega_{ k } \omega_{ m }}  \zeta_1^{  ( i )} \mathds{1}(l=j)    \\
& +
\zeta_1^{  ( j )} \zeta_1^{  ( k )} \zeta_1^{  ( m )}\mathds{1}(l=i)
+\zeta_1^{  ( i )}  \zeta_1^{  ( k )} \zeta_1^{  ( m )}\mathds{1}(l=j)
+\zeta_1^{  ( i )} \zeta_1^{  ( j )} \zeta_1^{  ( m )}\mathds{1}(l=k)
+\zeta_1^{  ( i )} \zeta_1^{  ( j )} \zeta_1^{  ( k )} \mathds{1}(l=m) , \\
\end{align*}

\begin{align*}
\frac{\partial}{\partial \zeta_2^{(l)}} 	P_{ijkm}(\boldsymbol{\zeta})
\end{align*}

is given by

\begin{align*}
% \frac{\partial}{\partial \zeta_2^{(l)}} 	P_{ijkm}(\boldsymbol{\zeta})  &=
%-\frac{ \zeta_1^{  ( k )} \zeta_1^{  ( m )}}{\omega_{ i } \omega_{ j }}
% \frac{\partial}{\partial \zeta_2^{(l)}}  \left[ \zeta_2^{  ( i )} \zeta_2^{  ( j )}\right]
%+\frac{ \zeta_1^{  ( j )}  \zeta_1^{  ( m )}}{\omega_{ i } \omega_{ k }}
%\frac{\partial}{\partial \zeta_2^{(l)}}  \left[\zeta_2^{  ( i )} \zeta_2^{  ( k )}\right]
%+\frac{\zeta_1^{  ( i )}  \zeta_1^{  ( m )}}{\omega_{ j } \omega_{ k }}
%\frac{\partial}{\partial \zeta_2^{(l)}}  \left[\zeta_2^{  ( j )} \zeta_2^{  ( k )} \right] \\
%& +\frac{ \zeta_1^{  ( j )} \zeta_1^{  ( k )} }{\omega_{ i } \omega_{ m }}
%\frac{\partial}{\partial \zeta_2^{(l)}}  \left[\zeta_2^{  ( i )}\zeta_2^{  ( m )} \right]
%+\frac{\zeta_1^{  ( i )}  \zeta_1^{  ( k )} }{\omega_{ j } \omega_{ m }}
%\frac{\partial}{\partial \zeta_2^{(l)}}  \left[\zeta_2^{  ( j )}\zeta_2^{  ( m )} \right]
%-\frac{\zeta_1^{  ( i )} \zeta_1^{  ( j )} }{\omega_{ k } \omega_{ m }}
%\frac{\partial}{\partial \zeta_2^{(l)}}  \left[ \zeta_2^{  ( k )} \zeta_2^{  ( m )} \right]\\
%&+\frac{1}{\omega_{ i } \omega_{ j } \omega_{ k } \omega_{ m }}
%\frac{\partial}{\partial \zeta_2^{(l)}}  \left[\zeta_2^{  ( i )} \zeta_2^{  ( j )} \zeta_2^{  ( k )} \zeta_2^{  ( m )} \right] \\
% %%%%%%%%%
% &=
& -\frac{ \zeta_1^{  ( k )} \zeta_1^{  ( m )}}{\omega_{ i } \omega_{ j }}
\zeta_2^{  ( j )} \mathds{1}(l=i)
-\frac{ \zeta_1^{  ( k )} \zeta_1^{  ( m )}}{\omega_{ i } \omega_{ j }} \zeta_2^{  ( i )}\mathds{1}(l=j)
+\frac{ \zeta_1^{  ( j )}  \zeta_1^{  ( m )}}{\omega_{ i } \omega_{ k }}
\zeta_2^{  ( k )}  \mathds{1}(l=i)
+\frac{ \zeta_1^{  ( j )}  \zeta_1^{  ( m )}}{\omega_{ i } \omega_{ k }}  \zeta_2^{  ( i )}\mathds{1}(l=k)   \\
&+\frac{\zeta_1^{  ( i )}  \zeta_1^{  ( m )}}{\omega_{ j } \omega_{ k }}
\zeta_2^{  ( k )}  \mathds{1}(l=j)
+\frac{\zeta_1^{  ( i )}  \zeta_1^{  ( m )}}{\omega_{ j } \omega_{ k }}   \zeta_2^{  ( j )} \mathds{1}(l=k)
+\frac{ \zeta_1^{  ( j )} \zeta_1^{  ( k )} }{\omega_{ i } \omega_{ m }}
\zeta_2^{  ( m )}  \mathds{1}(l=i)
+\frac{ \zeta_1^{  ( j )} \zeta_1^{  ( k )} }{\omega_{ i } \omega_{ m }} \zeta_2^{  ( i )}\mathds{1}(l=m)    \\
&+\frac{\zeta_1^{  ( i )}  \zeta_1^{  ( k )} }{\omega_{ j } \omega_{ m }}
\zeta_2^{  ( m )}  \mathds{1}(l=j)
+\frac{\zeta_1^{  ( i )}  \zeta_1^{  ( k )} }{\omega_{ j } \omega_{ m }}  \zeta_2^{  ( j )}\mathds{1}(l=m)
-\frac{\zeta_1^{  ( i )} \zeta_1^{  ( j )} }{\omega_{ k } \omega_{ m }}
\zeta_2^{  ( m )}  \mathds{1}(l=k)
-\frac{\zeta_1^{  ( i )} \zeta_1^{  ( j )} }{\omega_{ k } \omega_{ m }}  \zeta_2^{  ( k )} \mathds{1}(l=m)   \\
&+\frac{  \zeta_2^{  ( j )} \zeta_2^{  ( k )} \zeta_2^{  ( m )}}{\omega_{ i } \omega_{ j } \omega_{ k } \omega_{ m }}\mathds{1}(l=i)
+\frac{\zeta_2^{  ( i )}  \zeta_2^{  ( k )} \zeta_2^{  ( m )}}{\omega_{ i } \omega_{ j } \omega_{ k } \omega_{ m }}\mathds{1}(l=j)
+\frac{ \zeta_2^{  ( i )} \zeta_2^{  ( j )}  \zeta_2^{  ( m )}}{\omega_{ i } \omega_{ j } \omega_{ k } \omega_{ m }} \mathds{1}(l=k)
+\frac{\zeta_2^{  ( i )} \zeta_2^{  ( j )} \zeta_2^{  ( k )}}{\omega_{ i } \omega_{ j } \omega_{ k } \omega_{ m }}  \mathds{1}(l=m),
\end{align*}

for all integers $i,j,k,m,l \geq 0$. Next, we restrict our attention to the rescaled first 1-mode $\boldsymbol{\xi}=(\xi_1,\xi_2)$ defined by (<ref>). In this case, the latter yields

\begin{align*}
\frac{\partial}{\partial \zeta_1^{(l)}} 	P_{ijkm}(\boldsymbol{\xi})  &= \kappa_{0}^3 \Big[
\mathds{1}(k=j=m=0)  \mathds{1}(l=i)
+ \mathds{1}(i=k=m=0)\mathds{1}(l=j) \\
&+ \mathds{1}(i=j=m=0)\mathds{1}(l=k)
+ \mathds{1}(i=j=k=0) \mathds{1}(l=m) \Big] , \\
\frac{\partial}{\partial \zeta_2^{(l)}} 	P_{ijkm}(\boldsymbol{\xi})  &=
\end{align*}

for all integers $i,j,k,m,l \geq 0$. Now, one can see that there exists a choice of $\lambda$ such that both equations in (<ref>) hold true for the rescaled first 1-mode. Indeed, the second equation in (<ref>) is trivially satisfied for all $\lambda$
whereas the first equation in (<ref>) is satisfied provided that

\begin{align}\label{ValueLagrangeMultiplierLambda}
\lambda=  - \frac{1}{2}.
\end{align}

Indeed, for all integers $l \geq 0$, we have that

\begin{align*}
-2\lambda \omega_{0}^2 \kappa_{0} \mathds{1}(l=0) &=	-2\lambda \omega_{l}^2  \xi_1^{(l)} =\frac{\partial  \langle f  \rangle (\boldsymbol{\xi}) }{\partial \zeta_1^{(l)}}  \\
&=\frac{3}{32}\sum_{i,j,k,m=0}^{\infty} C_{ijkm}
\frac{\partial  P_{ijkm}(\boldsymbol{\xi}) }{\partial \zeta_1^{(l)}} \mathds{1}(m=i+j-k \geq 0) \\
&=\frac{3}{32}\kappa_{0} ^3 \sum_{i=0}^{\infty} C_{i0 0 0}
\mathds{1}(l=i =0) +
\frac{3}{32}\kappa_{0} ^3\sum_{j=0}^{\infty} C_{0 j 0 0}
\mathds{1}( l= j=0  )\\
&+\frac{3}{32}\kappa_{0} ^3\sum_{k=0}^{\infty} C_{0 0 k 0}
\mathds{1}(l=k=0 )
+\frac{3}{32}\kappa_{0} ^3\sum_{ m=0}^{\infty} C_{0 0 0 m}
\mathds{1}(l=m=0  ) \\
&= \frac{3}{8} \kappa_{0} ^3\mathds{1}(l =0  )   C_{0 0 0 0} \Longleftrightarrow  \\
-2\lambda \omega_{0}^2  \mathds{1}(l=0) &= \frac{3}{8} \kappa_{0} ^2 \mathds{1}(l =0  )   C_{0 0 0 0} =\frac{3}{8}     \frac{8\omega_{0}^2}{3C_{00 0 0}   }\mathds{1}(l=0) C_{0 0 0 0} \Longleftrightarrow \lambda=-\frac{1}{2},
\end{align*}

where we also used the definition of the rescaled constant $\kappa_{0}$ from (<ref>). Finally, we use Lemma <ref> to compute the harmonic energy of the rescaled 1–modes and show that $\boldsymbol{\xi} \in \Sigma$. We obtain

\begin{align*}
& = \sum_{m=0}^{\infty}
( \xi_2^{(m)} )^2  +\sum_{m=0}^{\infty} \omega_{m}^2  ( \xi_1^{(m)})^2
= \sum_{m=0}^{\infty} \omega_{m}^2  \left( \kappa_{0}\mathds{1}(m=0) \right)^2 = \kappa_{0}^2 \omega_{0}^2,
\end{align*}

that completes the proof.

Next, we show that the rescaled first 1-mode generates a curve of critical points.

Let $\boldsymbol{\xi}=(\xi_1,\xi_2)$ be the rescaled first 1-mode  given by (<ref>),

\begin{align*}
\Sigma= \left \{\boldsymbol{\zeta}\in \mathcal{P}:~  h_{\Omega}(\boldsymbol{\zeta})=\kappa_{0}^2 \omega_{0}^2 \right\}
\end{align*}

the hypersurface of constant harmonic energy and

\begin{align*}
\Gamma_{\text{linear}}(\boldsymbol{\xi}) = \left \{\boldsymbol{\Phi}^{t}(\boldsymbol{\xi}):t \in \mathbb{R} \right\}
=\left \{
\begin{pmatrix}
\kappa_{0}  \cos(\omega_{0}t)   e_{0} \\
- \omega_{0}\kappa_{0}  \sin(\omega_{0}t) e_{0}
\end{pmatrix}
:t \in \mathbb{R} \right\}
\end{align*}

the trajectory of the solution to the linearized equation with $\boldsymbol{\xi} $ as initial data. Then,  $ \Gamma_{\text{linear}}(\boldsymbol{\xi}) \subset \Sigma$ and every point $\boldsymbol{\eta}\in \Gamma_{\text{linear}}(\boldsymbol{\xi}) $ is a critical point for the function $\langle f\rangle |_{\Sigma}$.

It follows from the fact that both $h_{\Omega}$ and $\langle f\rangle |_{\Sigma}$ are invariant under the linear flow.

Next, we study the tangent space of the hypersurface of constant  harmonic energy.

Let $\boldsymbol{\xi}=(\xi_1,\xi_2)$ be the rescaled first 1-mode  given by (<ref>),

\begin{align*}
\Sigma= \left \{\boldsymbol{\zeta}\in \mathcal{P}:~  h_{\Omega}(\boldsymbol{\zeta})=\kappa_{0}^2 \omega_{0}^2 \right\}
\end{align*}

the hypersurface of constant  harmonic energy and

\begin{align*}
\Gamma_{\text{linear}}(\boldsymbol{\xi}) = \left \{\boldsymbol{\Phi}^{t}(\boldsymbol{\xi}):t \in \mathbb{R} \right\}
=\left \{
\begin{pmatrix}
\kappa_{0}  \cos(\omega_{0}t)   e_{0} \\
- \omega_{0}\kappa_{0}  \sin(\omega_{0}t) e_{0}
\end{pmatrix}
:t \in  \mathbb{R}  \right\} \subset \Sigma
\end{align*}

the trajectory of the solution to the linearized equation with $\boldsymbol{\xi} $ as initial data. Then, for any point $\boldsymbol{\eta}\in \Gamma_{\text{linear}}(\boldsymbol{\xi})$, the tangent space of $  \Sigma$ is given by

\begin{align*}
T_{\boldsymbol{\eta}}\Sigma   =
\left \{\boldsymbol{X} =
\begin{pmatrix}
\sum_{\substack{ m=0 \ }}^{\infty} X_{1}^{(m)} e_{m} \\
\sum_{\substack{ m=0  }}^{\infty} X_{2}^{(m)} e_{m}
\end{pmatrix} :~ \sin (\omega_{0}t)
X_{2}^{(0)}  =\omega_{0} \cos (\omega_{0}t)X_{1}^{(0)} \right \}.
\end{align*}

Moreover, it splits into

\begin{align*}
T_{\boldsymbol{\eta}}\Sigma =T_{\boldsymbol{\eta}}\Gamma_{\text{linear}}(\boldsymbol{\xi}) \oplus \left(T_{\boldsymbol{\eta}}\Gamma_{\text{linear}}  (\boldsymbol{\xi})\right)^{\perp},
\end{align*}

\begin{align*}
T_{\boldsymbol{\eta}}\Gamma_{\text{linear}}(\boldsymbol{\xi}) & = \left \{\boldsymbol{X} \in T_{\boldsymbol{\eta}} \Sigma:~\boldsymbol{X}=
\mu_{0}
\begin{pmatrix}
\sin(\omega_{0}t)   e_{0} \\
\omega_{0} \cos(\omega_{0}t) e_{0}
\end{pmatrix},~\mu_{0} \in \mathbb{R}
\right\}, \\
\left(T_{\boldsymbol{\eta}}\Gamma_{\text{linear}} (\boldsymbol{\xi}) \right)^{\perp}  &  =\left \{\boldsymbol{X} \in T_{\boldsymbol{\eta}} \Sigma:~ X_{1}^{(0)}=X_{2}^{(0)}=0 \right \}  .
\end{align*}

Let $t \in \mathbb{R} $, set  $\boldsymbol{\eta}=\boldsymbol{\Phi}^{t}(\boldsymbol{\xi}) $ and consider the function

\begin{align*}
\eta_2 ^2 + \left(L^{1/2}\eta_{1}\right)^2 \right)  d\mu
\end{align*}

studied in Lemma <ref>. Lemma <ref> yields that $\boldsymbol{\eta}=\boldsymbol{\Phi}^{t}(\boldsymbol{\xi})\in \Gamma_{\text{linear}}(\boldsymbol{\xi}) \subset \Sigma$. A normal vector to $\Sigma$ at the point $\boldsymbol{\eta}$ is given by the gradient $\nabla_{\boldsymbol{\eta}} h_{\Omega}(\boldsymbol{\eta}) $ and hence

\begin{align*}
T_{\boldsymbol{\eta}}\Sigma   &=
\left \{\boldsymbol{X} \in T_{\boldsymbol{\eta}} \mathcal{P}:~ \langle \nabla_{\boldsymbol{\eta}} h_{\Omega}(\boldsymbol{\eta}),\boldsymbol{X} \rangle =0 \right \}
\left \{\boldsymbol{X} \in T_{\boldsymbol{\eta}} \mathcal{P}:~ d h_{\Omega}(\boldsymbol{\eta})(\boldsymbol{X})   =0 \right \}.
\end{align*}

Let $\boldsymbol{X}=(X_1,X_2) \in T_{\boldsymbol{\eta}} \mathcal{P}$ be expanded in terms of the eigenfunctions $\{e_m: m \geq 0\}$ to the linearization operator $L$,

\begin{align*}
X_{\alpha}=\sum_{m=0}^{\infty} X_{\alpha}^{(m)} e_{m}, \quad \alpha \in \{1,2\}.
\end{align*}

We compute

\begin{align*}
d h_{\Omega}(\boldsymbol{\eta})(\boldsymbol{X}) &=
\int_{0}^{\pi/2}\left(
2\eta_2 X_2  +2  L^{1/2}\eta_{1}  L^{1/2}X_{1}  \right)  d\mu \\
& =2
\int_{0}^{\pi/2}\left(
\eta_2 X_2  +  X_{1} L \eta_{1}   \right)  d\mu \\
&= -2\omega_{0}\kappa_{0}\sin (\omega_{0}t)
\int_{0}^{\pi/2}
e_{0} X_2 d\mu +2\kappa_{0}\cos (\omega_{0}t)  \int_{0}^{\pi/2}  X_{1} L  e_{0}     d\mu \\
&= -2\omega_{0}\kappa_{0}\sin (\omega_{0}t)
\sum_{m=0}^{\infty} X_{2}^{(m)}\int_{0}^{\pi/2}
e_{0}  e_{m} d\mu +2\kappa_{0}\cos (\omega_{0}t)\sum_{m=0}^{\infty} X_{1}^{(m)}  \int_{0}^{\pi/2} e_{m} L  e_{0}     d\mu \\
&= -2\omega_{0}\kappa_{0}\sin (\omega_{0}t)
\sum_{m=0}^{\infty} X_{2}^{(m)} \mathds{1}(m=0) +2\omega_{0}^2 \kappa_{0}\cos (\omega_{0}t)\sum_{m=0}^{\infty} X_{1}^{(m)}  \mathds{1}(m=0) \\
&= -2\omega_{0}\kappa_{0}\sin (\omega_{0}t)
X_{2}^{(0)}  +2\omega_{0}^2\kappa_{0}\cos (\omega_{0}t)X_{1}^{(0)} \\
&=2\omega_{0}\kappa_{0} \left[ -\sin (\omega_{0}t)
X_{2}^{(0)}  +\omega_{0} \cos (\omega_{0}t)X_{1}^{(0)} \right]
\end{align*}

and hence

\begin{align*}
T_{\boldsymbol{\eta}}\Sigma   =
\left \{\boldsymbol{X} =
\begin{pmatrix}
\sum_{\substack{ m=0 \ }}^{\infty} X_{1}^{(m)} e_{m} \\
\sum_{\substack{ m=0  }}^{\infty} X_{2}^{(m)} e_{m}
\end{pmatrix} :~ \sin (\omega_{0}t)
X_{2}^{(0)}  =\omega_{0} \cos (\omega_{0}t)X_{1}^{(0)} \right \}.
\end{align*}

Moreover, the vectors tangent to the curve $\Gamma_{\text{linear}}(\boldsymbol{\xi}) = \left \{\boldsymbol{\Phi}^{t}(\boldsymbol{\xi}):t \in \mathbb{R} \right\}$ are proportional to

\begin{align*}
\partial_{t}\boldsymbol{\Phi}^{t}(\boldsymbol{\xi})=\partial_{t}\begin{pmatrix}
\kappa_{0}  \cos(\omega_{0}t)   e_{0} \\
- \omega_{0}\kappa_{0}  \sin(\omega_{0}t) e_{0}
\end{pmatrix}=
\begin{pmatrix}
- \kappa_{0} \omega_{0} \sin(\omega_{0}t)   e_{0} \\
- \kappa_{0} \omega_{0}^2 \cos(\omega_{0}t) e_{0}
\end{pmatrix} \in T_{\boldsymbol{\eta}}\Gamma_{\text{linear}} (\boldsymbol{\xi}) \subset   T_{\boldsymbol{\eta}}\Sigma
\end{align*}

and hence

\begin{align*}
T_{\boldsymbol{\eta}}\Gamma_{\text{linear}}(\boldsymbol{\xi}) = \left \{\boldsymbol{X} \in T_{\boldsymbol{\eta}} \Sigma:~\boldsymbol{X}=
\mu_{0}
\begin{pmatrix}
\sin(\omega_{0}t)   e_{0} \\
\omega_{0} \cos(\omega_{0}t) e_{0}
\end{pmatrix},~\mu_{0} \in \mathbb{R}
\right\}.
\end{align*}

Finally, we compute

\begin{align*}
\left(T_{\boldsymbol{\eta}}\Gamma_{\text{linear}} (\boldsymbol{\xi}) \right)^{\perp}
&= \left\{\boldsymbol{X} \in T_{\boldsymbol{\eta}} \Sigma:~ \langle \boldsymbol{X},\boldsymbol{Y}\rangle =0,~\forall \boldsymbol{Y} \in T_{\boldsymbol{\eta}}\Gamma_{\text{linear}}(\boldsymbol{\xi}) \subset   T_{\boldsymbol{\eta}}\Sigma
\right \} \\
&= \left\{\boldsymbol{X} \in T_{\boldsymbol{\eta}} \Sigma:~ \langle \boldsymbol{X},\partial_{t}\boldsymbol{\Phi}^{t}(\boldsymbol{\xi})\rangle =0 \right \} \\
&= \left\{\boldsymbol{X} \in T_{\boldsymbol{\eta}} \Sigma:~\left  \langle \boldsymbol{X},
\begin{pmatrix}
\sin(\omega_{0}t)   e_{0} \\
\omega_{0} \cos(\omega_{0}t) e_{0}
\end{pmatrix} \right \rangle =0 \right \}.
\end{align*}

For every $\boldsymbol{X}=(X_1,X_2) \in T_{\boldsymbol{\eta}} \Sigma$, we use the orthogonality condition

\begin{align*}
\left  \langle e_m,   e_{0}\right \rangle_{L^2 ((0,\pi/2);d\mu )} &=(e _{m}  |e _{0}  )=  \mathds{1}(m=0),
\end{align*}

for all integers $n,m \geq 0$, derived in (<ref>), as well as the fact that

\begin{align*}
\left  \langle e_m,   e_{0}\right \rangle_{H^1((0,\pi/2);d\mu )} &=(e^{\prime}_{m}  |e^{\prime}_{0}  )+(e _{m}  |e _{0}  )=  \mathds{1}(m=0),
\end{align*}

since $e_0$ is constant in the KG model, to compute

\begin{align*}
& \left \langle \boldsymbol{X},
\begin{pmatrix}
\sin(\omega_{0}t)   e_{0} \\
\omega_{0} \cos(\omega_{0}t) e_{0}
\end{pmatrix} \right \rangle =\sin(\omega_{0}t)
\left  \langle X_1,   e_{0}\right \rangle_{H^1((0,\pi/2);d\mu )}+
\omega_{0} \cos(\omega_{0}t)
\left  \langle X_2,e_{0} \right \rangle_{L^2((0,\pi/2);d\mu )} \\
& =\sin(\omega_{0}t)
\sum_{m=0}^{\infty}X_{1}^{(m)} \left  \langle e_m,   e_{0}\right \rangle_{H^1((0,\pi/2);d\mu )}+
\omega_{0} \cos(\omega_{0}t)
\sum_{m=0}^{\infty}X_{2}^{(m)}
\left  \langle e_m ,e_{0} \right \rangle_{L^2((0,\pi/2);d\mu )} \\
&= \sin(\omega_{0}t)
X_{1}^{(0)} +
\omega_{0} \cos(\omega_{0}t)
X_{2}^{(0)}  .
\end{align*}

Hence, the tangent vectors  $\boldsymbol{X}  \in   T_{\boldsymbol{\eta}} \Sigma \cap \left(T_{\boldsymbol{\eta}}\Gamma_{\text{linear}}(\boldsymbol{\xi})  \right)^{\perp} $ are given by the linear system

\begin{align*}
\begin{pmatrix}
\omega_{0} \cos(\omega_{0}t)  &
- \sin(\omega_{0}t) \\
\sin  (\omega_{0}t) &
\omega_{0} \cos(\omega_{0}t)
\end{pmatrix}
\begin{pmatrix}
X_{1}^{(0)} \\ X_{2}^{(0)}
\end{pmatrix}=
\begin{pmatrix}
0 \\ 0
\end{pmatrix}.
\end{align*}

The determinant is given by $\omega_{0}^2 \cos^2 (\omega_{0}t)+\sin^2 (\omega_{0}t) $ and it is non-zero for all $t \in \mathbb{R}$ that leads us to the unique solution $X_{1}^{(0)}=X_{2}^{(0)}=0$. We conclude that

\begin{align*}
\left(T_{\boldsymbol{\zeta}}\Gamma_{\text{linear}} (\boldsymbol{\xi}) \right)^{\perp} =\left \{\boldsymbol{X} =
\begin{pmatrix}
\sum_{\substack{ m=0 \ }}^{\infty} X_{1}^{(m)} e_{m} \\
\sum_{\substack{ m=0  }}^{\infty} X_{2}^{(m)} e_{m}
\end{pmatrix} :~ X_{1}^{(0)}=X_{2}^{(0)}=0 \right \},
\end{align*}

that completes the proof.

Next, we compute the second differential of $ \langle f\rangle $ restricted on a hypersurface of constant hormonic energy and evaluated at the rescaled first 1–mode.

Let $\boldsymbol{\xi}=(\xi_1,\xi_2)$ be the rescaled first 1-mode  given by (<ref>) and let

\begin{align*}
\Sigma= \left \{\boldsymbol{\zeta}\in \mathcal{P}:~  h_{\Omega}(\boldsymbol{\zeta})=\kappa_{0}^2 \omega_{0}^2 \right\}
\end{align*}

be the hypersurface of constant  harmonic energy. Then, for any tangent vector

\begin{align*}
\boldsymbol{X}=(V,W)=\left(\sum_{m=0}^{\infty}V^{(m)}e_{m},\sum_{m=0}^{\infty}W^{(m)}e_{m} \right)  \in T_{\boldsymbol{\xi}}\Sigma,
\end{align*}

we have that

\begin{align*}
d^2  ( \langle f\rangle |_{\Sigma} )(\boldsymbol{\xi})(\boldsymbol{X},\boldsymbol{X}) & =  - \frac{\omega_{0}^2}{C_{0 0 0 0}} \Bigg[  \sum_{l=0}^{\infty} \Bigg(  \frac{C_{0 0 0 0}}{\omega_{0}^2} - 2 \frac{ C_{l l0   0       } }{ \omega_{l}^2} \Bigg) ( \omega_{l}V^{(l)})^2
+\sum_{l=0}^{\infty} \Bigg( \frac{C_{0 0 0 0}}{\omega_{0}^2} -  2\frac{  C_{ l l 0   0  } }{\omega_{l}^2     }    \Bigg)(W^{(l)})^2   \Bigg] \\
&+\omega_0^2 (V^{(0)})^2
-     \left(   W^{(0)}  \right)^2.
\end{align*}

Pick any $\boldsymbol{\zeta}=(\zeta_1,\zeta_2) \in \mathcal{P}$, expand its coefficients in terms of the eigenfunctions $\{e_m:m\geq 0 \}$ to the linearized operator $L$,

\begin{align*}
\zeta_\alpha =\sum_{m=0}^{\infty} \zeta_\alpha^{(m)} e_m,\quad  \alpha \in \{1,2\},
\end{align*}

and consider the functions

\begin{align*}
h_{\Omega}(\boldsymbol{\zeta}) &=
h_{\Omega} (\zeta_1^{(0)},\zeta_1^{(1)},\zeta_1^{(2)},\dots, \zeta_2^{(0)}, \zeta_2^{(1)}, \zeta_2^{(2)},\dots), \\
\langle f  \rangle (\boldsymbol{\zeta})&=
\langle f  \rangle (\zeta_1,\zeta_2)=
\langle f  \rangle (\zeta_1^{(0)},\zeta_1^{(1)},\zeta_1^{(2)},\dots, \zeta_2^{(0)}, \zeta_2^{(1)}, \zeta_2^{(2)},\dots),
\end{align*}

given by Lemma <ref> and <ref> respectively. Firstly, we use the computations for

\begin{align*}
\frac{\partial}{\partial \zeta_1^{(l)}} 	P_{ijkm}(\boldsymbol{\zeta}) , \quad \frac{\partial}{\partial \zeta_2^{(l)}} 	P_{ijkm}(\boldsymbol{\zeta})
\end{align*}

for all integers $l \geq 0$, derived in Lemma <ref>, to compute

\begin{align*}
\frac{\partial^2}{\partial \zeta_1^{(l)} \zeta_1^{(n)}  } 	P_{ijkm}(\boldsymbol{\zeta}) , \quad
\frac{\partial^2}{\partial \zeta_1^{(l)} \zeta_2^{(n)}   } 	P_{ijkm}(\boldsymbol{\zeta}) , \quad
\frac{\partial^2}{\partial \zeta_2^{(l)} \zeta_2^{(n)}  } 	P_{ijkm}(\boldsymbol{\zeta}),
\end{align*}

for all integers $l,n \geq 0$, and for generic $\boldsymbol{\zeta}=(\zeta_1,\zeta_2)$. We obtain that

\begin{align*}
\frac{\partial^2}{\partial \zeta_1^{(l)}\zeta_1^{(n)} } 	P_{ijkm}(\boldsymbol{\zeta})
\end{align*}

is given by

\begin{align*}
&-\frac{\zeta_2^{  ( i )} \zeta_2^{  ( j )}  }{\omega_{ i } \omega_{ j }}
\mathds{1}(n=m) \mathds{1}(l=k)
-\frac{\zeta_2^{  ( i )} \zeta_2^{  ( j )}  }{\omega_{ i } \omega_{ j }}  \mathds{1}(n=k) \mathds{1}(l=m)
+\frac{\zeta_2^{  ( i )}  \zeta_2^{  ( k )}  }{\omega_{ i } \omega_{ k }}
\mathds{1}(n=m)\mathds{1}(l=j) \\
& +\frac{\zeta_2^{  ( i )}  \zeta_2^{  ( k )}  }{\omega_{ i } \omega_{ k }} \mathds{1}(n=j) \mathds{1}(l=m)
+\frac{ \zeta_2^{  ( j )} \zeta_2^{  ( k )} }{\omega_{ j } \omega_{ k }}
\mathds{1}(n=m) \mathds{1}(l=i)
+\frac{ \zeta_2^{  ( j )} \zeta_2^{  ( k )} }{\omega_{ j } \omega_{ k }}    \mathds{1}(n=i) \mathds{1}(l=m)   \\
&  +\frac{\zeta_2^{  ( i )}   \zeta_2^{  ( m )}}{\omega_{ i } \omega_{ m }}
\mathds{1}(n=k) \mathds{1}(l=j)
+\frac{\zeta_2^{  ( i )}   \zeta_2^{  ( m )}}{\omega_{ i } \omega_{ m }}    \mathds{1}(n= j) \mathds{1}(l=k)
+\frac{ \zeta_2^{  ( j )}  \zeta_2^{  ( m )}}{\omega_{ j } \omega_{ m }}
\mathds{1}(n= k) \mathds{1}(l=i)  \\
&  +\frac{ \zeta_2^{  ( j )}  \zeta_2^{  ( m )}}{\omega_{ j } \omega_{ m }} \mathds{1}(n= i) \mathds{1}(l=k)
-\frac{ \zeta_2^{  ( k )} \zeta_2^{  ( m )}}{\omega_{ k } \omega_{ m }}
\mathds{1}(n= j)  \mathds{1}(l=i)
-\frac{ \zeta_2^{  ( k )} \zeta_2^{  ( m )}}{\omega_{ k } \omega_{ m }}  \mathds{1}(n= i) \mathds{1}(l=j)    \\
& +  \mathds{1}(n= j) \zeta_1^{  ( k )} \zeta_1^{  ( m )}\mathds{1}(l=i)
+  \zeta_1^{  ( j )} \mathds{1}(n= k) \zeta_1^{  ( m )}\mathds{1}(l=i)
+  \zeta_1^{  ( j )} \zeta_1^{  ( k )} \mathds{1}(n= m)\mathds{1}(l=i) \\
&+\mathds{1}(n= i)  \zeta_1^{  ( k )} \zeta_1^{  ( m )}\mathds{1}(l=j)
+\zeta_1^{  ( i )} \mathds{1}(n= k) \zeta_1^{  ( m )}\mathds{1}(l=j)
+\zeta_1^{  ( i )}  \zeta_1^{  ( k )} \mathds{1}(n= m)\mathds{1}(l=j) \\
& +\mathds{1}(n= i) \zeta_1^{  ( j )} \zeta_1^{  ( m )}\mathds{1}(l=k)
+\zeta_1^{  ( i )} \mathds{1}(n= j) \zeta_1^{  ( m )}\mathds{1}(l=k)
+\zeta_1^{  ( i )} \zeta_1^{  ( j )} \mathds{1}(n= m)\mathds{1}(l=k) \\
& +\mathds{1}(n= i) \zeta_1^{  ( j )} \zeta_1^{  ( k )} \mathds{1}(l=m)
+\zeta_1^{  ( i )} \mathds{1}(n= j) \zeta_1^{  ( k )} \mathds{1}(l=m)
+\zeta_1^{  ( i )} \zeta_1^{  ( j )} \mathds{1}(n= k)\mathds{1}(l=m) ,
\end{align*}

\begin{align*}
\frac{\partial^2}{\partial \zeta_1^{(l)} \partial \zeta_2^{(n)}  } 	P_{ijkm}(\boldsymbol{\zeta})
\end{align*}

is given by

\begin{align*}
& -\frac{ \zeta_2^{  ( j )}  }{\omega_{ i } \omega_{ j }}
\zeta_1^{  ( m )}\mathds{1}(n= i)\mathds{1}(l=k)
-\frac{\zeta_2^{  ( i )}  }{\omega_{ i } \omega_{ j }}
\zeta_1^{  ( m )}\mathds{1}(n= j)\mathds{1}(l=k)
-\frac{ \zeta_2^{  ( j )}  }{\omega_{ i } \omega_{ j }} \zeta_1^{  ( k )} \mathds{1}(n= i) \mathds{1}(l=m) \\
& -\frac{\zeta_2^{  ( i )}   }{\omega_{ i } \omega_{ j }} \zeta_1^{  ( k )} \mathds{1}(n= j)\mathds{1}(l=m)
+\frac{ \zeta_2^{  ( k )}  }{\omega_{ i } \omega_{ k }}
\zeta_1^{  ( m )} \mathds{1}(n= i) \mathds{1}(l=j)
+\frac{\zeta_2^{  ( i )}    }{\omega_{ i } \omega_{ k }}
\zeta_1^{  ( m )} \mathds{1}(n= k) \mathds{1}(l=j) \\
& +\frac{ \zeta_2^{  ( k )}  }{\omega_{ i } \omega_{ k }}\zeta_1^{  ( j )} \mathds{1}(n= i)  \mathds{1}(l=m)
+\frac{\zeta_2^{  ( i )}    }{\omega_{ i } \omega_{ k }}\zeta_1^{  ( j )}\mathds{1}(n= k) \mathds{1}(l=m)
+\frac{   \zeta_2^{  ( k )} }{\omega_{ j } \omega_{ k }}
\zeta_1^{  ( m )}\mathds{1}(n= j)  \mathds{1}(l=i)\\
&  +\frac{ \zeta_2^{  ( j )}  }{\omega_{ j } \omega_{ k }}
\zeta_1^{  ( m )}\mathds{1}(n= k)  \mathds{1}(l=i)
+\frac{   \zeta_2^{  ( k )} }{\omega_{ j } \omega_{ k }}   \zeta_1^{  ( i )} \mathds{1}(n= j) \mathds{1}(l=m)
+\frac{ \zeta_2^{  ( j )}   }{\omega_{ j } \omega_{ k }}   \zeta_1^{  ( i )}
\mathds{1}(n= k) \mathds{1}(l=m)   \\
&  +\frac{   \zeta_2^{  ( m )}}{\omega_{ i } \omega_{ m }}
\zeta_1^{  ( k )}\mathds{1}(n= i)  \mathds{1}(l=j)
+\frac{\zeta_2^{  ( i )}  }{\omega_{ i } \omega_{ m }}
\zeta_1^{  ( k )} \mathds{1}(n= m) \mathds{1}(l=j)
+\frac{   \zeta_2^{  ( m )}}{\omega_{ i } \omega_{ m }}   \zeta_1^{  ( j )}  \mathds{1}(n= i)\mathds{1}(l=k)  \\
& +\frac{\zeta_2^{  ( i )}   }{\omega_{ i } \omega_{ m }}   \zeta_1^{  ( j )} \mathds{1}(n= m)\mathds{1}(l=k)
+\frac{  \zeta_2^{  ( m )}}{\omega_{ j } \omega_{ m }}
\zeta_1^{  ( k )} \mathds{1}(n= j) \mathds{1}(l=i)
+\frac{ \zeta_2^{  ( j )}}{\omega_{ j } \omega_{ m }}
\zeta_1^{  ( k )}  \mathds{1}(n= m) \mathds{1}(l=i) \\
&  +\frac{    \zeta_2^{  ( m )}}{\omega_{ j } \omega_{ m }} \zeta_1^{  ( i )}  \mathds{1}(n= j) \mathds{1}(l=k)
+\frac{ \zeta_2^{  ( j )}   }{\omega_{ j } \omega_{ m }} \zeta_1^{  ( i )}  \mathds{1}(n= m) \mathds{1}(l=k)
-\frac{  \zeta_2^{  ( m )}}{\omega_{ k } \omega_{ m }}
\zeta_1^{  ( j )}  \mathds{1}(n= k) \mathds{1}(l=i)\\
&  -\frac{ \zeta_2^{  ( k )}  }{\omega_{ k } \omega_{ m }}
\zeta_1^{  ( j )} \mathds{1}(n= m) \mathds{1}(l=i)
-\frac{   \zeta_2^{  ( m )}}{\omega_{ k } \omega_{ m }}  \zeta_1^{  ( i )} \mathds{1}(n= k) \mathds{1}(l=j)
-\frac{ \zeta_2^{  ( k )} }{\omega_{ k } \omega_{ m }}  \zeta_1^{  ( i )} \mathds{1}(n= m) \mathds{1}(l=j)    ,\\
\end{align*}

\begin{align*}
\frac{\partial^2}{\partial \zeta_2^{(n)} \partial \zeta_2^{(l)}} 	P_{ijkm}(\boldsymbol{\zeta})
\end{align*}

is given by

\begin{align*}
%  ws pros \zeta_2^{(n)}    \mathds{1}(n=  )
& -\frac{ \zeta_1^{  ( k )} \zeta_1^{  ( m )}}{\omega_{ i } \omega_{ j }}
\mathds{1}(n=  j)\mathds{1}(l=i)
-\frac{ \zeta_1^{  ( k )} \zeta_1^{  ( m )}}{\omega_{ i } \omega_{ j }} \mathds{1}(n=  i) \mathds{1}(l=j)
+\frac{ \zeta_1^{  ( j )}  \zeta_1^{  ( m )}}{\omega_{ i } \omega_{ k }}
\mathds{1}(n=  k) \mathds{1}(l=i) \\
& +\frac{ \zeta_1^{  ( j )}  \zeta_1^{  ( m )}}{\omega_{ i } \omega_{ k }}  \mathds{1}(n=  i)\mathds{1}(l=k)
+\frac{\zeta_1^{  ( i )}  \zeta_1^{  ( m )}}{\omega_{ j } \omega_{ k }}
\mathds{1}(n= k ) \mathds{1}(l=j)
+\frac{\zeta_1^{  ( i )}  \zeta_1^{  ( m )}}{\omega_{ j } \omega_{ k }}   \mathds{1}(n=  j) \mathds{1}(l=k)  \\
& +\frac{ \zeta_1^{  ( j )} \zeta_1^{  ( k )} }{\omega_{ i } \omega_{ m }}
\mathds{1}(n=  m) \mathds{1}(l=i)
+\frac{ \zeta_1^{  ( j )} \zeta_1^{  ( k )} }{\omega_{ i } \omega_{ m }} \mathds{1}(n=  i)\mathds{1}(l=m)
+\frac{\zeta_1^{  ( i )}  \zeta_1^{  ( k )} }{\omega_{ j } \omega_{ m }}
\mathds{1}(n=  m) \mathds{1}(l=j) \\
&  +\frac{\zeta_1^{  ( i )}  \zeta_1^{  ( k )} }{\omega_{ j } \omega_{ m }}  \mathds{1}(n=  j)\mathds{1}(l=m)
-\frac{\zeta_1^{  ( i )} \zeta_1^{  ( j )} }{\omega_{ k } \omega_{ m }}
\mathds{1}(n=  m) \mathds{1}(l=k)
-\frac{\zeta_1^{  ( i )} \zeta_1^{  ( j )} }{\omega_{ k } \omega_{ m }}  \mathds{1}(n=  k) \mathds{1}(l=m)   \\
&+\frac{ \zeta_2^{  ( k )} \zeta_2^{  ( m )}}{\omega_{ i } \omega_{ j } \omega_{ k } \omega_{ m }}\mathds{1}(l=i)\mathds{1}(n= j)
+\frac{  \zeta_2^{  ( j )}  \zeta_2^{  ( m )}}{\omega_{ i } \omega_{ j } \omega_{ k } \omega_{ m }}\mathds{1}(l=i)\mathds{1}(n= k)
+\frac{  \zeta_2^{  ( j )} \zeta_2^{  ( k )} }{\omega_{ i } \omega_{ j } \omega_{ k } \omega_{ m }}\mathds{1}(l=i)\mathds{1}(n= m) \\
& +\frac{  \zeta_2^{  ( k )} \zeta_2^{  ( m )}}{\omega_{ i } \omega_{ j } \omega_{ k } \omega_{ m }}\mathds{1}(n= i) \mathds{1}(l=j)
+\frac{\zeta_2^{  ( i )} \zeta_2^{  ( m )}}{\omega_{ i } \omega_{ j } \omega_{ k } \omega_{ m }} \mathds{1}(l=j) \mathds{1}(n= k)
+\frac{\zeta_2^{  ( i )}  \zeta_2^{  ( k )}}{\omega_{ i } \omega_{ j } \omega_{ k } \omega_{ m }} \mathds{1}(n= m) \mathds{1}(l=j) \\
& +\frac{ \zeta_2^{  ( j )}  \zeta_2^{  ( m )}}{\omega_{ i } \omega_{ j } \omega_{ k } \omega_{ m }} \mathds{1}(n= i) \mathds{1}(l=k)
+\frac{ \zeta_2^{  ( i )} \zeta_2^{  ( m )}}{\omega_{ i } \omega_{ j } \omega_{ k } \omega_{ m }} \mathds{1}(l=k) \mathds{1}(n= j)
+\frac{ \zeta_2^{  ( i )} \zeta_2^{  ( j )} }{\omega_{ i } \omega_{ j } \omega_{ k } \omega_{ m }} \mathds{1}(l=k) \mathds{1}(n= m) \\
&+\frac{ \zeta_2^{  ( j )} \zeta_2^{  ( k )}}{\omega_{ i } \omega_{ j } \omega_{ k } \omega_{ m }}  \mathds{1}(l=m)\mathds{1}(n= i)
+\frac{\zeta_2^{  ( i )}  \zeta_2^{  ( k )}}{\omega_{ i } \omega_{ j } \omega_{ k } \omega_{ m }}  \mathds{1}(l=m)\mathds{1}(n= j)
+\frac{\zeta_2^{  ( i )} \zeta_2^{  ( j )}}{\omega_{ i } \omega_{ j } \omega_{ k } \omega_{ m }}  \mathds{1}(l=m) \mathds{1}(n= k),
\end{align*}

for all integers $i,j,k,m,n,l \geq 0$. Secondly, we restrict our attention to the rescaled first 1-mode $\boldsymbol{\xi}=(\xi_1,\xi_2)$ defined by (<ref>). In this case, the latter yields that

\begin{align*}
\frac{\partial^2}{\partial \zeta_1^{(l)}\zeta_1^{(n)} } 	P_{ijkm}(\boldsymbol{\xi})
\end{align*}

is given by

\begin{align*}
& + \kappa_{0}^2 \mathds{1}(k=m=0)  \mathds{1}(n= j) \mathds{1}(l=i)
+  \kappa_{0}^2 \mathds{1}(j=m=0)  \mathds{1}(n= k)  \mathds{1}(l=i) \\
&  + \kappa_{0}^2 \mathds{1}(j=k=0)  \mathds{1}(n= m)\mathds{1}(l=i)
+\kappa_{0}^2 \mathds{1}(k=m=0) \mathds{1}(n= i)  \mathds{1}(l=j) \\
& +\kappa_{0}^2 \mathds{1}(i=m=0)  \mathds{1}(n= k)  \mathds{1}(l=j)
+\kappa_{0}^2 \mathds{1}(i=k=0) \mathds{1}(n= m)\mathds{1}(l=j) \\
& +\kappa_{0}^2 \mathds{1}(j=m=0) \mathds{1}(n= i) \mathds{1}(l=k)
+\kappa_{0}^2 \mathds{1}(i=m=0)  \mathds{1}(n= j) \mathds{1}(l=k) \\
& +\kappa_{0}^2 \mathds{1}(i=j=0)  \mathds{1}(n= m)\mathds{1}(l=k)
+\kappa_{0}^2 \mathds{1}(j=k=0) \mathds{1}(n= i)   \mathds{1}(l=m)\\
& + \kappa_{0}^2 \mathds{1}(i=k=0)  \mathds{1}(n= j)   \mathds{1}(l=m)
+\kappa_{0}^2 \mathds{1}(i=j=0)  \mathds{1}(n= k)\mathds{1}(l=m) ,
\end{align*}

\begin{align*}
\frac{\partial^2}{\partial \zeta_1^{(l)} \partial \zeta_2^{(n)}  } 	P_{ijkm}(\boldsymbol{\xi})    =0
\end{align*}

\begin{align*}
\frac{\partial^2}{\partial \zeta_2^{(n)} \partial \zeta_2^{(l)}} 	P_{ijkm}(\boldsymbol{\xi})
\end{align*}

is given by

\begin{align*}
& -    \frac{  \kappa_{0}^2}{\omega_{ i } \omega_{ j }}
\mathds{1}( k= m= 0) \mathds{1}(n=  j)\mathds{1}(l=i)
-\frac{ \kappa_{0}^2 }{\omega_{ i } \omega_{ j }} \mathds{1}( k= m= 0)\mathds{1}(n=  i) \mathds{1}(l=j) \\
& +\frac{ \kappa_{0}^2 }{\omega_{ i } \omega_{ k }}
\mathds{1}( j= m= 0) \mathds{1}(n=  k) \mathds{1}(l=i)
+\frac{ \kappa_{0}^2 }{\omega_{ i } \omega_{ k }} \mathds{1}( j= m= 0) \mathds{1}(n=  i)\mathds{1}(l=k) \\
&+\frac{\kappa_{0}^2  }{\omega_{ j } \omega_{ k }}
\mathds{1}( i= m= 0) \mathds{1}(n= k ) \mathds{1}(l=j)
+\frac{\kappa_{0}^2 }{\omega_{ j } \omega_{ k }} \mathds{1}(i = m= 0)   \mathds{1}(n=  j) \mathds{1}(l=k)  \\
& +\frac{ \kappa_{0}^2  }{\omega_{ i } \omega_{ m }}
\mathds{1}( j= k= 0)\mathds{1}(n=  m) \mathds{1}(l=i)
+\frac{ \kappa_{0}^2  }{\omega_{ i } \omega_{ m }}\mathds{1}( j= k= 0) \mathds{1}(n=  i)\mathds{1}(l=m)  \\
& +\frac{\kappa_{0}^2 }{\omega_{ j } \omega_{ m }}
\mathds{1}( i= k= 0)\mathds{1}(n=  m) \mathds{1}(l=j)
+\frac{ \kappa_{0}^2  }{\omega_{ j } \omega_{ m }}  \mathds{1}(n=  j)\mathds{1}( i= k= 0)\mathds{1}(l=m)  \\
& -\frac{\kappa_{0}^2  }{\omega_{ k } \omega_{ m }}
\mathds{1}(i =j = 0)\mathds{1}(n=  m) \mathds{1}(l=k)
-\frac{\kappa_{0}^2 }{\omega_{ k } \omega_{ m }} \mathds{1}(i =j = 0)  \mathds{1}(n=  k) \mathds{1}(l=m)  ,
\end{align*}

for all integers $i,j,k,m,n,l \geq 0$. Now, restrict our attention to the rescaled first 1-mode, combine the three previous results together with Lemma <ref> and the symmetries for the Fourier coefficients from (<ref>) and (<ref>) to obtain

\begin{align*}
\frac{\partial^2 \langle f\rangle (\boldsymbol{\xi})}{\partial \zeta_1^{(l)} \partial \zeta_1^{(n)}}  &= \frac{3}{32}\sum_{i,j,k,m=0}^{\infty} C_{ijkm}
\frac{\partial^2 P_{ijkm}(\boldsymbol{\xi})}{\partial \zeta_1^{(l)} \partial \zeta_1^{(n)}}  \mathds{1}(m=i+j-k \geq 0)  \\
&=  \frac{3}{32} \kappa_{0}^2\Big[
C_{ln 00 }    \mathds{1}(0=l+n \geq 0)
+  C_{l0 n 0  }     \mathds{1}(n=l \geq 0)
+  C_{l 0 0 n }  \mathds{1}(n=l \geq 0) \\
&+ C_{ n l 00 }   \mathds{1}(0=n+l \geq 0)
+ C_{0 l n 0 }    \mathds{1}(n=l \geq 0)
+ C_{0 l 0 n  }   \mathds{1}(n=l \geq 0)   \\
+ C_{ n 0 l 0  }  \mathds{1}(l=n \geq 0)
+ C_{ 0 n l 0  }  \mathds{1}(l=n \geq 0)
+ C_{ 00 l n }    \mathds{1}(n=-l \geq 0)   \\
+ C_{ n 00 l  }   \mathds{1}(l=n \geq 0)
+  C_{0 n 0 l}     \mathds{1}(l=n \geq 0)
+ C_{00 n l }    \mathds{1}(l=-n \geq 0)
\Big]    \\
&=  \frac{3}{32} \kappa_{0}^2\Big[
C_{ln 00 }    \mathds{1}(l=n=0 )
+  C_{l0 n 0  }     \mathds{1}( n=l )  \\
& +  C_{l 0 0 n }  \mathds{1}(n=l )
+ C_{ n l 00 }   \mathds{1}(l=n=0)
+ C_{0 l n 0 }    \mathds{1}(  n=l  )
+ C_{0 l 0 n  }   \mathds{1}(n= l   )   \\
& + C_{ n 0 l 0  }  \mathds{1}(n=l  )
+ C_{ 0 n l 0  }  \mathds{1}(n=l  )
+ C_{ 00 l n }    \mathds{1}(n=-l   )
+ C_{ n 00 l  }   \mathds{1}(l=n  ) \\
& +  C_{0 n 0 l}     \mathds{1}(l= n  )
+ C_{00 n l }    \mathds{1}(l=n=0)	\Big]    \\
&=  \frac{3}{32} \kappa_{0}^2\Big[
4 C_{ln 00 }    \mathds{1}(l=n=0 )
+ 8 C_{l0 n 0  }     \mathds{1}( n=l )  	\Big] \\
&=  \frac{3}{8} \kappa_{0}^2\Big[
C_{ln 00 }    \mathds{1}(l=n=0 )
+ 2 C_{l0 n 0  }     \mathds{1}( n=l )  	\Big]  ,  \\
\frac{\partial^2 \langle f\rangle (\boldsymbol{\xi})}{\partial \zeta_2^{(n)} \partial \zeta_1^{(l)}}  &=  \frac{3}{32}\sum_{i,j,k,m=0}^{\infty} C_{ijkm}
\frac{\partial^2  P_{ijkm}(\boldsymbol{\xi})}{\partial \zeta_2^{(n)} \partial \zeta_1^{(l)}}   \mathds{1}(m=i+j-k \geq 0)
=  0,\\
\frac{\partial^2 \langle f\rangle (\boldsymbol{\xi})}{\partial \zeta_2^{(n)} \partial \zeta_2^{(l)}}   &=   \frac{3}{32}\sum_{i,j,k,m=0}^{\infty} C_{ijkm}
\frac{\partial^2 P_{ijkm}(\boldsymbol{\xi}) }{\partial \zeta_2^{(n)} \partial \zeta_2^{(l)}}   \mathds{1}(m=i+j-k \geq 0)  \\
&=  \frac{3}{32}\kappa_{0}^2 \Big[
-    \frac{  C_{ ln 0 0 }  }{\omega_{ l } \omega_{ n }}
\mathds{1}(0 =l+n-0  \geq 0)
-\frac{C_{n l 0 0   }  }{\omega_{ n } \omega_{ l }}
\mathds{1}(-l =n  \geq 0)
+\frac{ C_{ l 0  n 0  }  }{\omega_{ l } \omega_{ n }}
\mathds{1}(n =l  \geq 0)  \\
+\frac{ C_{ n 0  l 0  } }{\omega_{ n} \omega_{ l }}
\mathds{1}(l =n  \geq 0)
+\frac{ C_{0  l n 0   } }{\omega_{ l } \omega_{ n}}
\mathds{1}(n = l \geq 0)
+\frac{ C_{ 0  n l 0  } }{\omega_{ n } \omega_{ l }}
\mathds{1}(l = n \geq 0)   \\
& +\frac{  C_{l 0  0  n }  }{\omega_{ l } \omega_{ n }}
\mathds{1}(n=l  \geq 0)
+\frac{  C_{ n 0 0 l }  }{\omega_{ n } \omega_{ l }}
\mathds{1}(l=n  \geq 0)
+\frac{ C_{ 0  l 0  n } }{\omega_{ l } \omega_{ n}}
\mathds{1}(n=l  \geq 0)
+\frac{  C_{ 0 n 0  l }  }{\omega_{ n } \omega_{ l }}
\mathds{1}(l=n  \geq 0)    \\
& -\frac{ C_{ 0  0  l n  } }{\omega_{ l } \omega_{ n }}
\mathds{1}(n=-l \geq 0)
-\frac{ C_{ 0  0  n l } }{\omega_{ n } \omega_{ l }}
\mathds{1}(l= -n \geq 0)  \Big] \\
&=  \frac{3}{32}\kappa_{0}^2 \Big[
-    \frac{  C_{ ln 0 0 }  }{\omega_{ l } \omega_{ n }}
\mathds{1}( l=n=0 )
-\frac{C_{n l 0 0   }  }{\omega_{ n } \omega_{ l }}
\mathds{1}(l=n=0)
+\frac{ C_{ l 0  n 0  }  }{\omega_{ l } \omega_{ n }}
\mathds{1}(l=n) \\
&+\frac{ C_{ n 0  l 0  } }{\omega_{ n} \omega_{ l }}
\mathds{1}(n=l)  +\frac{ C_{0  l n 0   } }{\omega_{ l } \omega_{ n}}
\mathds{1}(l=n)
+\frac{ C_{ 0  n l 0  } }{\omega_{ n } \omega_{ l }}
\mathds{1}(n=l)
+\frac{  C_{l 0  0  n }  }{\omega_{ l } \omega_{ n }}
\mathds{1}(n=l ) \\
& +\frac{  C_{ n 0 0 l }  }{\omega_{ n } \omega_{ l }}
\mathds{1}(l=n )
+\frac{ C_{ 0  l 0  n } }{\omega_{ l } \omega_{ n}}
\mathds{1}(n= l )
+\frac{  C_{ 0 n 0  l }  }{\omega_{ n } \omega_{ l }}
\mathds{1}(l= n )    \\
& -\frac{ C_{ 0  0  l n  } }{\omega_{ l } \omega_{ n }}
\mathds{1}(l=n=0)
-\frac{ C_{ 0  0  n l } }{\omega_{ n } \omega_{ l }}
\mathds{1}(l=n=0)  \Big] \\
&=  \frac{3}{32}\kappa_{0}^2 \Big[
-  4  \frac{  C_{ ln 0 0 }  }{\omega_{ l } \omega_{ n }}
\mathds{1}( l=n=0 )
+8\frac{ C_{ l 0  n 0  }  }{\omega_{ l } \omega_{ n }}
\mathds{1}(l=n)  \Big] \\
&=  \frac{3}{8 } \frac{\kappa_{0}^2}{\omega_{ l } \omega_{ n }} \Big[
-     C_{ ln 0 0 }
\mathds{1}( l=n=0 )
+2 C_{ l 0  n 0  }    \mathds{1}(l=n)  \Big] ,
\end{align*}

for all integers $n,l \geq 0$. Moreover, for any tangent vector

\begin{align*}
\boldsymbol{X}=(V,W)=\left(\sum_{m=0}^{\infty}V^{(m)}e_{m},\sum_{m=0}^{\infty}W^{(m)}e_{m} \right)  \in T_{\boldsymbol{\xi}}\Sigma,
\end{align*}

given by Lemma <ref>, the latter yields

\begin{align*}
d^2  &( \langle f\rangle   )(\boldsymbol{\xi})(\boldsymbol{X},\boldsymbol{X}) = \sum_{ \alpha,\beta=0}^{\infty}  \frac{\partial^2  \langle f\rangle (\boldsymbol{\xi})}{ \partial \boldsymbol{\zeta}_{\alpha} \partial \boldsymbol{\zeta}_{\beta} } \boldsymbol{X}^{\alpha} \boldsymbol{X}^{\beta} \\
&=\sum_{n,l=0}^{\infty} \Bigg[
\frac{\partial^2 \langle f\rangle (\boldsymbol{\xi})}{\partial \zeta_1^{(n)} \partial \zeta_1^{(l)}}  V^{(n)}V^{(l)}+
2\frac{\partial^2 \langle f\rangle (\boldsymbol{\xi})}{\partial \zeta_2^{(n)} \partial \zeta_1^{(l)}}  W^{(n)}V^{(l)}+
\frac{\partial^2 \langle f\rangle (\boldsymbol{\xi})}{\partial \zeta_2^{(n)} \partial \zeta_2^{(l)}}  W^{(n)}W^{(l)}
\Bigg] \\
&=\frac{3}{8} \kappa_{0}^2 \sum_{n,l=0}^{\infty} \Bigg[  \big(
C_{ln 00 }    \mathds{1}(l=n=0 )
+ 2 C_{l0 n 0  }     \mathds{1}( n=l )  	\big) V^{(n)}V^{(l)} \\
&+    \big(
-     C_{ ln 0 0 }
\mathds{1}( l=n=0 )
+2 C_{ l 0  n 0  }    \mathds{1}(l=n)  \big)\frac{W^{(n)}W^{(l)}}{\omega_{ n } \omega_{ l }}  \Bigg] \\
&=\frac{3}{8} \kappa_{0}^2 \Bigg[
C_{0000}   (V^{(0)})^2
+ 2 \sum_{l=0}^{\infty} C_{l0 l 0  } (V^{(l)})^2
-       C_{ 0000 }  \left(\frac{W^{(0)} }{\omega_{ 0 } }   \right)^2
+2   \sum_{l=0}^{\infty} C_{ l 0  l 0  }   \left( \frac{W^{(l)}}{\omega_{ l }  } \right)^2 \Bigg] \\
&=  \frac{\omega_{0}^2}{C_{00 0 0}   } \Bigg[
C_{0000 }  \left(V^{(0)}\right)^2
+ 2 \sum_{l=0}^{\infty} C_{l0 l 0  } (V^{(l)})^2
-         C_{ 0000} \left( \frac{W^{(0)} }{\omega_{0 } }  \right)^2
+2   \sum_{l=0}^{\infty} C_{ l 0  l 0  }   \left( \frac{W^{(l)}}{\omega_{ l }  } \right)^2 \Bigg] \\
&= \omega_{0}^2 \left(V^{(0)}\right)^2
+ 2 \omega_{0}^2\sum_{l=0}^{\infty} \frac{C_{l0 l 0  }}{C_{0000}} (V^{(l)})^2
-         \left( W^{(0)}  \right)^2
+2   \omega_{0}^2 \sum_{l=0}^{\infty}  \frac{C_{l0 l 0  }}{C_{0000}}    \left( \frac{W^{(l)}}{\omega_{ l }  } \right)^2  ,
\end{align*}

where we also used the definition of the rescaled constant $\kappa_{0}$ given by (<ref>). Thirdly, we use Lemma <ref> to compute

\begin{align*}
\frac{\partial^2}{\partial \zeta_1^{(l)} \zeta_1^{(n)}  } h_{\Omega} (\boldsymbol{\zeta}) , \quad
\frac{\partial^2}{\partial \zeta_1^{(l)} \zeta_2^{(n)}   } h_{\Omega}(\boldsymbol{\zeta}) , \quad
\frac{\partial^2}{\partial \zeta_2^{(l)} \zeta_2^{(n)}  } 	h_{\Omega}(\boldsymbol{\zeta}) ,
\end{align*}

for all integers $l,n \geq 0$, and for any $\boldsymbol{\zeta}=(\zeta_1,\zeta_2)$. We obtain that

\begin{align*}
\frac{\partial^2}{\partial \zeta_1^{(l)} \zeta_1^{(n)}  } h_{\Omega} (\boldsymbol{\zeta}) = 2\omega_{l}^2 \mathds{1}(l=n)  ,  \quad
\frac{\partial^2}{\partial \zeta_1^{(l)} \zeta_2^{(n)}   } h_{\Omega}(\boldsymbol{\zeta})  = 0,  \quad
\frac{\partial^2}{\partial \zeta_2^{(l)} \zeta_2^{(n)}  } 	h_{\Omega}(\boldsymbol{\zeta}) = 2  \mathds{1}(l=n) ,
\end{align*}

for all integers $l,n \geq 0$. Moreover, restricting our attention to the rescaled first 1-mode $\boldsymbol{\xi}=(\xi_1,\xi_2)$ given by (<ref>), for any tangent vector

\begin{align*}
\boldsymbol{X}=(V,W)=\left(\sum_{m=0}^{\infty}V^{(m)}e_{m},\sum_{m=0}^{\infty}W^{(m)}e_{m} \right)  \in T_{\boldsymbol{\xi}}\Sigma,
\end{align*}

given by Lemma <ref>, the latter yields

\begin{align*}
d^2  h_{\Omega}(\boldsymbol{\xi})(\boldsymbol{X},\boldsymbol{X}) &=
\sum_{ \alpha,\beta=0}^{\infty}  \frac{\partial^2  h_{\Omega}(\boldsymbol{\xi})}{ \partial \boldsymbol{\zeta}_{\alpha} \partial \boldsymbol{\zeta}_{\beta} } \boldsymbol{X}^{\alpha} \boldsymbol{X}^{\beta} \\
&=\sum_{n,l=0}^{\infty} \Bigg[
\frac{\partial^2 h_{\Omega}(\boldsymbol{\xi})}{\partial \zeta_1^{(n)} \partial \zeta_1^{(l)}}  V^{(n)}V^{(l)}+
2\frac{\partial^2 h_{\Omega}(\boldsymbol{\xi})}{\partial \zeta_2^{(n)} \partial \zeta_1^{(l)}}  W^{(n)}V^{(l)}+
\frac{\partial^2 h_{\Omega}(\boldsymbol{\xi})}{\partial \zeta_2^{(n)} \partial \zeta_2^{(l)}}  W^{(n)}W^{(l)}
\Bigg] \\
&=\sum_{n,l=0}^{\infty} \Bigg[
2\omega_{l}^2 \mathds{1}(l=n)  V^{(n)}V^{(l)}
+   2  \mathds{1}(l=n) W^{(n)}W^{(l)}
\Bigg] \\
\omega_{l}^2  (  V^{(l)}  )^2
+   2 \sum_{l=0}^{\infty}   (  W^{(l)} )^2   .
\end{align*}

Finally, according to Remark page 78-79 in [5], we have that

\begin{align*}
d^2  ( \langle f\rangle |_{\Sigma} )(\boldsymbol{\xi})(\boldsymbol{X},\boldsymbol{X})  = \lambda  d^2  h_{\Omega}(\boldsymbol{\xi})(\boldsymbol{X},\boldsymbol{X}) + d^2  ( \langle f\rangle   )(\boldsymbol{\xi})(\boldsymbol{X},\boldsymbol{X}) ,
\end{align*}

with $\lambda=-1/2$ due to (<ref>), and hence we conclude that

\begin{align*}
d^2 &  ( \langle f\rangle |_{\Sigma} )(\boldsymbol{\xi})(\boldsymbol{X},\boldsymbol{X})  = -\frac{1}{2}  d^2  h_{\Omega}(\boldsymbol{\xi})(\boldsymbol{X},\boldsymbol{X}) + d^2  ( \langle f\rangle   )(\boldsymbol{\xi})(\boldsymbol{X},\boldsymbol{X}) \\
& = -\sum_{l=0}^{\infty}
\omega_{l}^2  (  V^{(l)}  )^2
-   \sum_{l=0}^{\infty}   (  W^{(l)} )^2  +
\omega_{0}^2      (V^{(0)} )^2
+  2\omega_{0}^2\sum_{l=0}^{\infty} \frac{C_{l0 l 0  } }{C_{00 0 0}   } (V^{(l)})^2 \\
& 	 -      \left(  W^{(0)}     \right)^2
+  2\omega_{0}^2 \sum_{l=0}^{\infty} \frac{ C_{ l 0  l 0  } }{C_{00 0 0}   }   \left( \frac{W^{(l)}}{\omega_{ l }  } \right)^2  \\
&  =  - \frac{\omega_{0}^2}{C_{0 0 0 0}} \Bigg[  \sum_{l=0}^{\infty} \Bigg(  \frac{C_{0 0 0 0}}{\omega_{0}^2} - 2 \frac{ C_{l l0   0       } }{ \omega_{l}^2} \Bigg) ( \omega_{l}V^{(l)})^2
+\sum_{l=0}^{\infty} \Bigg( \frac{C_{0 0 0 0}}{\omega_{0}^2} -  2\frac{  C_{ l l 0   0  } }{\omega_{l}^2     }    \Bigg)(W^{(l)})^2   \Bigg]\\
&  +  \omega_0^2 (V^{(0)})^2
-     \left(   W^{(0)}  \right)^2    ,
\end{align*}

that completes the proof.

Finally, we verify Condition (<ref>) in Theorem <ref> which states that the rescaled first 1-mode must also be an extremum point for the function $\langle f\rangle |_{\Sigma}$, namely there exists a positive constant $c_{\perp}$ such that

\begin{align*}
\pm  d^2  ( \langle f\rangle |_{\Sigma} )(\boldsymbol{\xi})(\boldsymbol{X},\boldsymbol{X}) \geq c_{\perp} \| \boldsymbol{X}\|^2, \quad \forall \boldsymbol{X} \in \left(T_{\boldsymbol{\xi}}\Gamma_{\text{linear}}(\boldsymbol{\xi})  \right)^{\perp},
\end{align*}

with the plus and minus sign if $\boldsymbol{\xi} \in \Sigma$ is a minimum or maximum for the function $\langle f\rangle |_{\Sigma}$ respectively. Recall that the orthogonal complement $\left(T_{\boldsymbol{\xi}}\Gamma_{\text{linear}} (\boldsymbol{\xi}) \right)^{\perp}$ is given by Lemma <ref>. As we will see next, in our case, the rescaled first 1-mode maximizes $\langle f\rangle |_{\Sigma}$ for both KG and WM if and only if $\delta$ satisfies Assumption <ref>.

Let  $\boldsymbol{\xi}=(\xi_1,\xi_2)$ be the rescaled first 1-mode  given by (<ref>) and let

\begin{align*}
\Sigma= \left \{\boldsymbol{\zeta}\in \mathcal{P}:~  h_{\Omega}(\boldsymbol{\zeta})=\kappa_{0}^2 \omega_{0}^2 \right\}
\end{align*}

be the hypersurface of constant  harmonic energy. Then, $\delta$ satisfies Assumption <ref> if and only if there exists a positive constant $c_{\perp}=c_{\perp}(\delta)$, depending only on $\delta$, such that

\begin{align*}
-d^2  ( \langle f\rangle |_{\Sigma} )(\boldsymbol{\xi})(\boldsymbol{X},\boldsymbol{X}) \geq c_{\perp} \| \boldsymbol{X}\|^2,  \quad \forall \boldsymbol{X} \in \left(T_{\boldsymbol{\xi}}\Gamma_{\text{linear}}  (\boldsymbol{\xi}) \right)^{\perp}.
\end{align*}

Let  $\boldsymbol{\xi}=(\xi_1,\xi_2)$ be the rescaled first 1-mode  given by (<ref>), pick any $\boldsymbol{X}=(V,W)\in  \left(T_{\boldsymbol{\xi}}\Gamma_{\text{linear}}  (\boldsymbol{\xi})\right)^{\perp} $ and expand its coefficients terms of the eigenfunctions $\{e_m:m\geq 0 \}$ to the linearized operator $L$,

\begin{align*}
V=\sum_{m=0}^{\infty} V^{(m)} e_{m}, \quad W=\sum_{m=0}^{\infty} W^{(m)} e_{m}.
\end{align*}

Then, according to Lemma <ref>, we have that

\begin{align*}
V^{(0)}=W^{(0)}=0 .
\end{align*}

For convenience, we rescale $V$ and set

\begin{align*}
Z^n = \omega_{n}V^n,
\end{align*}

for all integers $n \geq 0$. On the one hand, according to Lemma <ref>, we have that

\begin{align*}
- \frac{C_{0 0 0 0}}{\omega_{0}^2}  d^2  ( \langle f\rangle |_{\Sigma} )(\boldsymbol{\xi})(\boldsymbol{X},\boldsymbol{X}) & =     \sum_{l=1}^{\infty} \Bigg(   \frac{C_{0 0 0 0}}{\omega_{0}^2} -2 \frac{ C_{l l0   0       }  }{ \omega_{l}^2} \Bigg) (Z^{(l)})^2
+\sum_{l=1}^{\infty} \Bigg(   \frac{C_{0 0 0 0}}{\omega_{0}^2} - 2 \frac{  C_{ l l 0   0  } }{\omega_{l}^2     }  \Bigg)(W^{(l)})^2 .
\end{align*}

On the other hand, we have that

\begin{align*}
\| \boldsymbol{X}\|_{\mathcal{P}}^2 &=
% \| \boldsymbol{X}\|_{ H^1((0,\pi/2) ;d\mu) \times L^2((0,\pi/2) ;d\mu) }^2 =  \sum_{\substack{ l=1}}^{\infty}\left(  (\omega_{l} V^{(l)}))^2 +  (W^{(l)})^2 \right)=
\sum_{\substack{ l=1   }}^{\infty}\left(  (Z^{(l)})^2 +  (W^{(l)})^2 \right).
\end{align*}

Let $\{U^{(m)} : m \geq 1\}$ be either $\{Z^{(m)} : m \geq 1\}$ or $\{W^{(m)} : m \geq 1\}$. Then,  we deduce that

\begin{align*}
\exists c_{\perp}>0:\Hquad  -d^2  ( \langle f\rangle |_{\Sigma} )(\boldsymbol{\xi})(\boldsymbol{X},\boldsymbol{X}) &\geq c_{\perp} \| \boldsymbol{X}\|^2, \quad \forall	\boldsymbol{X} \in \left(T_{\boldsymbol{\xi}}\Gamma_{\text{linear}} (\boldsymbol{\xi}) \right)^{\perp}  \Longleftrightarrow \\
\exists \tilde{c}_{\perp}>0:\Hquad  \sum_{\substack{ l=1   }}^{\infty} \Bigg(   \frac{C_{0 0 0 0}}{\omega_{0}^2} - 2 \frac{ C_{l l0   0       } }{ \omega_{l}^2}\Bigg) ( U^{(l)})^2
& \geq \tilde{c}_{\perp} \sum_{\substack{ l=1   }}^{\infty}  ( U^{(l)})^2, \quad \forall	 \{U^{(m)}: m \geq 1 \} \subset \mathbb{R} \Longleftrightarrow \\
\exists \tilde{c}_{\perp}>0:\Hquad     \frac{C_{0 0 0 0}}{\omega_{0}^2} - 2 \frac{ C_{l l0   0       } }{ \omega_{l}^2}
& \geq \tilde{c}_{\perp} ,
\end{align*}

which, according to Proposition <ref>, is equivalent to the Assumption <ref> for $\delta$, that completes the proof.

§ ZEILBERGER'S ALGORITHM

In this section, we discuss Zeilberger's algorithm. In Section <ref>, we follow Sections 5 and 6 in [47] to describe how it works in general. In Section <ref>, we fix $\delta$ and run it through by hand producing step-by-step the results of Lemma <ref>.

§.§ How the algorithm works

Let $F=F(n,k)$ be a given hypergeometric function with respect to both arguments $n$ and $k$. This means that both $F(n+1,k)/F(n,k)$ and $F(n,k+1)/F(n,k)$ are rational functions of $n$ and $k$ and hence there exist polynomials $r_1,r_2,s_1,s_2$ so that

\begin{align}\label{Fhypergeometricfunction}
\frac{F(n+1,k)}{F(n,k)}=\frac{s_1(n+1,k)}{s_2(n+1,k)} ,\quad \frac{F(n,k+1)}{F(n,k)}=\frac{r_1(n,k)}{r_2(n,k)}.
\end{align}

We are interested in finding a recurrence relation for the sum

\begin{align}\label{DefinitionSumf}
f(n)=\sum_{k=0}^{K} F(n,k).
\end{align}

For simplicity, we assume that $K$ is a finite integer. However, as it will become clear later, this assumption can be relaxed. We claim that it suffices to find a recurrence relation for the summand $F(n,k)$. In other words, given $F$, we would like to find a function $G=G(n,k)$ in the form

\begin{align}\label{AnsatzforG}
\end{align}

for some function $R=R(n,k)$, as well as coefficients $\{\alpha_{j}(n):j =0,1,\dots,J\}$, for some (possibly large) integer $J$, all in closed forms, so that

\begin{align}\label{RecurrenceRelationF}
\sum_{j=0}^{J} \alpha_{j}(n)F(n+j,k)=G(n,k+1)-G(n,k).
\end{align}

Then, if (<ref>) is possible, a recurrence relation for $f(n)$, given by (<ref>), follows immediately by summing both sides of (<ref>) with respect to $k$. Indeed, since the coefficients $\{\alpha_{j}(n):j =0,1,\dots,J\}$ are all independent of $k$, summing (<ref>) with respect to $k$ yields

\begin{align}\label{RecurrenceRelationSmallf}
\sum_{j=0}^{J} \alpha_{j}(n)F(n+j,k) & =G(n,k+1)-G(n,k) \Longrightarrow \nonumber \\
\sum_{k=0}^{K} \left(	\sum_{j=0}^{J} \alpha_{j}(n)F(n+j,k) \right) &=\sum_{k=0}^{K} \left(G(n,k+1)-G(n,k) \right) \Longrightarrow \nonumber \\
\sum_{j=0}^{J}\alpha_{j}(n)\left(	\sum_{k=0}^{K} F(n+j,k) \right) &= G(n,K+1)-G(n,0)  \Longrightarrow \nonumber \\
\sum_{j=0}^{J}\alpha_{j}(n)f(n+j) &=G(n,K+1)-G(n,0),
\end{align}

which is the desired recurrence relation for $f(n)$. Consequently, we turn our attention to how one can find the function $G=G(n,k)$ and coefficients $\{\alpha_{j}(n):j =0,1,\dots,J\}$ so that (<ref>) holds[We note that, for all proper hypergeometric functions $F$, a recurrence relation of the form (<ref>) always exists (Theorem 6.2.1, page 105 in [47]).]. To this end, we fix $J$ and denote the left-hand-side of (<ref>) by $t_k$,

\begin{align}\label{DefinitiontkStart}
t_k=\sum_{j=0}^{J} \alpha_{j}(n)F(n+j,k).
\end{align}

Using (<ref>), we infer

\begin{align*}
\frac{t_{k+1}}{t_{k}} & = \frac{\sum_{j=0}^{J} \alpha_{j}(n)F(n+j,k+1)}{\sum_{j=0}^{J} \alpha_{j}(n)F(n+j,k)}
= \frac{\sum_{j=0}^{J} \alpha_{j}(n)\frac{F(n+j,k+1)}{F(n,k+1)}}{\sum_{j=0}^{J} \alpha_{j}(n)
\frac{F(n+j,k)}{F(n,k)}
\frac{F(n,k+1)}{F(n,k)}  \\
& = \frac{\sum_{j=0}^{J} \alpha_{j}(n)\prod_{i=0}^{j-1} \frac{F(n+j-i,k+1)}{F(n+j-i-1,k+1)}}{\sum_{j=0}^{J} \alpha_{j}(n)\prod_{i=0}^{j-1} \frac{F(n+j-i,k)}{F(n+j-i-1,k)}
} \frac{F(n,k+1)}{F(n,k)}
= \frac{\sum_{j=0}^{J} \alpha_{j}(n)\prod_{i=0}^{j-1} \frac{s_1(n+j-i,k+1)}{s_2(n+j-i,k+1)}}{\sum_{j=0}^{J} \alpha_{j}(n)\prod_{i=0}^{j-1} \frac{s_1(n+j-i,k)}{s_2(n+j-i,k)}
} \frac{r_1(n,k)}{r_2(n,k)}  \\
& = \frac{\sum_{j=0}^{J} \alpha_{j}(n) \left[
\prod_{i=0}^{j-1}  s_1(n+j-i,k+1)
\prod_{r=j+1}^{J}  s_2(n+r,k+1)
\right]}{\sum_{j=0}^{J} \alpha_{j}(n) \left[
\prod_{i=0}^{j-1}  s_1(n+j-i,k)
\prod_{r=j+1}^{J}  s_2(n+r,k)
\right]
} \frac{r_1(n,k)}{r_2(n,k)} \frac{\prod_{r=1}^{J}s_2(n+r,k) }{\prod_{r=1}^{J}s_2(n+r,k+1)} .
\end{align*}

Hence, we have that

\begin{align*}
\frac{t_{k+1}}{t_{k}} & = \frac{p_0(k+1)}{p_0(k)} \frac{r(k)}{s(k)},
\end{align*}

\begin{align}
p_0( k) &= \sum_{j=0}^{J} \left[ \prod_{i=0}^{j-1}  s_1(n+j-i,k) \prod_{r=j+1}^{J}  s_2(n+r,k) \right],\label{Defintionp0}  \\
r( k) &=r_1(n,k)  \prod_{r=1}^{J}  s_2(n+r,k)  ,\label{Defintionr} \\
s( k) &=r_2(n,k)  \prod_{r=1}^{J}  s_2(n+r,k+1)  .\label{Defintions}
\end{align}

Next, we use Theorem 5.3.1, page 82 in [47]
to find polynomials $p_1,p_2,p_3$ with $\text{gcd}(p_2(k),p_3(k+j))=1$, for all $j \in \{0,1,2,\dots\}$ so that

\begin{align} \label{Dividerands}
\frac{r(k)}{s(k)} = \frac{p_1(k+1)}{p_1(k)} \frac{p_2(k)}{p_3(k)}.
\end{align}

Then, we set $p(k)=p_0(k)p_1(k)$ and have that

\begin{align}\label{FinalEquationtk}
\frac{t_{k+1}}{t_{k}} & = \frac{p_0(k+1)}{p_0(k)} \frac{r(k)}{s(k)}
= \frac{p_0(k+1)}{p_0(k)} \frac{p_1(k+1)}{p_1(k)} \frac{p_2(k)}{p_3(k)}
= \frac{p(k+1)}{p(k)} \frac{p_2(k)}{p_3(k)}.
\end{align}

Now, one can show that (<ref>) has a indefinitely summable hypergeometric solution $t_k$ if and only if the recurrence relation

\begin{align}\label{Recurrencebk}
\end{align}

has a polynomial solution $b(k)$. For a detailed proof of this result, see Section 5.2 and Theorem 5.2.1 on page 76 in [47]
. We look for a polynomial solution to (<ref>) of degree $\Delta$ and set

\begin{align*}
b(k)=\sum_{l=0}^{\Delta} \beta_{l}k^l,
\end{align*}

where the coefficients $\{\beta_l:l =0,1,\dots,\Delta\}$ are to be determined. By substituting this expression for $b(k)$ into (<ref>), one finds a system of simultaneous linear equations for $\Delta+J+2$ unknowns $\alpha_{0},\dots,\alpha_{J}$ and $\beta_{0},\dots,\beta_{\Delta}$. Finally, we solve this linear system, find the coefficients $\alpha_{0},\dots,\alpha_{J}$ and $\beta_{0},\dots,\beta_{\Delta}$ and set

\begin{align}\label{DefinitionofG}
G(n,k)=\frac{p_3(k-1)}{p(k)} b(k) t_k.
\end{align}

For more details, see Sections 5 and 6 in [47].

§.§ The case of Lemma <ref>

Next, we produce step-by-step the results of Lemma <ref> by running Zeilberger's algorithm through by hand. Due to the complexity of the formulas appearing, we fix the conformal mass and set $\delta=3$ and we are looking for a recurrence relation with $J=2$ steps for the sum

\begin{align*}
\mathsf{f}_{m}=	 \frac{ \mathsf{C}_{00 mm}}{\somega_{m}^2} =\sum_{k=0}^{2m+1} \mathsf{F}(m,k)
\end{align*}

that is in the form of (<ref>). Here, $\mathsf{F}(m,k)$ is given by a closed formula as stated in Lemma <ref>. Firstly, one can use the fact that $\Gamma(x+1)=x\Gamma(x)$, valid for all $x \in \mathbb{R}$, to compute

\begin{align*}
\frac{\mathsf{F}(m,k)}{\mathsf{F}(m-1,k)} & = \frac{(2 m+1) (2 k-4 m-1) (2 k-4 m+1) (k+2 m+3) (k+2 m+4)}{(2 m+3) (k-2 m-1) (k-2 m) (2 k+4 m+7) (2 k+4 m+9)}, \\
\frac{\mathsf{F}(m,k+1)}{\mathsf{F(m,k)}} & =\frac{(k+2) (2 k-1) (2 k+1) (2 k+7) (4 k+11) (k-2 m-1) (k+2 m+5)}{(k+1) (k+4) (k+5) (2 k+5) (4 k+7) (2 k-4 m-1) (2 k+4 m+11)}.
\end{align*}

Since are both rational functions we conclude that $\mathsf{F}(m,k)$ is a hypergeometric function both respect to $k$ and $m$ and hence Zeilberger's algorithm applies. According to the computations above and (<ref>), we define

\begin{align*}
s_1(m,k) &= (2 m+1) (2 k-4 m-1) (2 k-4 m+1) (k+2 m+3) (k+2 m+4),\\
s_2 (m,k)& =(2 m+3) (k-2 m-1) (k-2 m) (2 k+4 m+7) (2 k+4 m+9), \\
r_1 (m,k)& = (k+2) (2 k-1) (2 k+1) (2 k+7) (4 k+11) (k-2 m-1) (k+2 m+5)  ,\\
r_2 (m,k)& = (k+1) (k+4) (k+5) (2 k+5) (4 k+7) (2 k-4 m-1) (2 k+4 m+11) .
\end{align*}

According to (<ref>), we set

\begin{align*}
t_{k}=   \alpha_{0}(m)F(m,k)+\alpha_{1}(m)F(m+1,k)+\alpha_{2}(m)F(m+2,k).
\end{align*}

From now, we simplify the notation and write $\alpha_j$ instead of $\alpha_j(m)$, for all $j \in \{0,1,2\}$. Using the definitions (<ref>), (<ref>) and (<ref>), we compute

\begin{align*}
p_0(k)  &= \alpha_0 (2 m+5) (2 m+7) (-k+2 m+2) (-k+2 m+3) (-k+2 m+4) (-k+2 m+5)\cdot \\
& \cdot (2 k+4 m+11) (2 k+4 m+13) (2 k+4 m+15) (2 k+4 m+17)
\\
+(2 m+3) (k+2 m+5) (k+2 m+6) (2 k-4 m-5) (2 k-4 m-3) (\alpha_1 (2 m+7)\cdot \\
& \cdot  (-k+2 m+4) (-k+2 m+5) (2 k+4 m+15) (2 k+4 m+17)   \\
& +\alpha_2 (2 m+5) (k+2 m+7) (k+2 m+8) (-2 k+4 m+7) (-2 k+4 m+9)), \\
r(k) &= (k+2) (2 k-1) (2 k+1) (2 k+7) (4 k+11) (2 m+5) (2 m+7) (k-2 m-5) (k-2 m-3)\cdot \\
& \cdot  (k-2 m-1) (k+2 m+5) (2 k+4 m+11) (2 k+4 m+13) (2 k+4 m+15) (2 k+4 m+17) \cdot \\
& \cdot (k-2 (m+1)) (k-2 (m+2)), \\
s(k) &=(k+1) (k+4) (k+5) (2 k+5) (4 k+7) (2 m+5) (2 m+7) (2 k-4 m-1) (k-2 m-3)\cdot \\
& \cdot  (k-2 m-1) (2 k+4 m+11) (2 k+4 m+13) (2 k+4 m+15) (2 k+4 m+17)\cdot \\
& \cdot  (2 k+4 m+19) (k-2 (m+1)) (k-2 (m+2)).
\end{align*}

Notice that the only trace of the unknowns $\alpha_0$, $\alpha_1$ and $\alpha_2$ is in the polynomial $p_0(k)$ and not in $r(k)$ nor in $s(k)$. In addition, observe that $r(k)$ and $s(k)$ have many common factors. By dividing these two, the common terms are canceled and one obtains the canonical form

\begin{align*}
\frac{r(k)}{s(k)} &=\frac{(k+2) (2 k-1) (2 k+1) (2 k+7) (4 k+11) (k-2 m-5) (k+2 m+5)}{(k+1) (k+4) (k+5) (2 k+5) (4 k+7) (2 k-4 m-1) (2 k+4 m+19)}.
\end{align*}

This is of the form (<ref>) with

\begin{align}
p_1  (k) &= (1 + k) (5 + 2 k) (7 + 4 k)  ,\nonumber  \\
p_2  (k) &=  (-1 + 2 k) (1 + 2 k) (-5 + k - 2 m) (5 + k + 2 m),\nonumber \\
p_3  (k) &=  (k+2) (k+3) (2 k-4 m-1) (2 k+4 m+15) . \label{Definitionp3}
\end{align}

We set

\begin{align*}
\end{align*}

so that (<ref>) holds, and compute

\begin{align}
p(k) &=(k+1) (2 k+5) (4 k+7) \big[\alpha_0 (2 m+5) (2 m+7) (-k+2 m+2) (-k+2 m+3)\cdot \nonumber \\
& \cdot (-k+2 m+4) (-k+2 m+5) (2 k+4 m+11) (2 k+4 m+13) (2 k+4 m+15)\cdot \nonumber \\
& \cdot  (2 k+4 m+17)+(2 m+3) (k+2 m+5) (k+2 m+6) (2 k-4 m-5) (2 k-4 m-3) \cdot \nonumber \\
& \cdot \big[\alpha_1 (2 m+7) (-k+2 m+4) (-k+2 m+5) (2 k+4 m+15) (2 k+4 m+17)
\nonumber \\
&  +
\alpha_2 (2 m+5) (k+2 m+7) (k+2 m+8) (-2 k+4 m+7) (-2 k+4 m+9) \big]  \Big]. \label{Definitionp}
\end{align}

Notice that the only trace of the unknowns $\alpha_0$, $\alpha_1$ and $\alpha_2$ is in the polynomial $p(k)$ and not in $p_1(k)$, $p_2(k)$ nor in $p_3(k)$. Now, we focus on whether the recurrence relation (<ref>) has a polynomial solution $b(k)$. On the one hand, notice that both $p_2(k)$ and $p_3(k)$ are polynomials of degree $4$ with respect to $k$. On the other hand, $p(k)$ is a polynomial of degree $11$ with respect to $k$. Hence, $b(k)$ should be a polynomial of degree at least $7$ with respect to $k$. However, we relax the algorithm[This is due to some cancellations in (<ref>) for these particular choices of $p_2$, $p_3$ and $p$.] and search for a solution $b(k)$ that is a polynomial of degree $\Delta=8$ with respect to $k$ and define

\begin{align}\label{Definitionb}
b(k)=\beta_8 k^8+\beta_7 k^7+ \beta_6 k^6+\beta_5 k^5+\beta_4 k^4+\beta_3 k^3+\beta_2 k^2+\beta_1 k + \beta_0,
\end{align}

where the coefficients $\{\beta_j: j \in \{0,1,\dots,8\}\}$ are all to be determined (together with the $\{\alpha_0,\alpha_1,\alpha_2\}$). Next, we plug (<ref>) into the recurrence relation (<ref>) and match the coefficients of same powers of $k$ on both sides. After a long computation, the result is a linear system consisting of $12$ equations for the $12$ unknowns  $\{\alpha_0,\alpha_1,\alpha_2,\beta_0,\beta_1,\dots,\beta_8\}$. This is given by

Equation 1:

\begin{align*}
& \left[-140 (1 + m) (2 + m) (3 + 2 m) (5 + 2 m)^2 (7 + 2 m) (11 + 4 m) (13 +
4 m) (15 + 4 m) (17 + 4 m)  \right] \alpha_0 + \\
&	\left[ -140 (2 + m) (3 + m) (3 + 2 m) (5 + 2 m)^2 (7 + 2 m) (3 + 4 m) (5 +
4 m) (15 + 4 m) (17 + 4 m) \right] \alpha_1 + \\
&	\left [ -140 (3 + m) (4 + m) (3 + 2 m) (5 + 2 m)^2 (7 + 2 m) (3 + 4 m) (5 +
4 m) (7 + 4 m) (9 + 4 m) \right] \alpha_2 + \\
&	\left[49 (13 + 4 m (5 + m))  \right] \beta_0 +
\left[ (5 + 2 m)^2 \right] \beta_1 +
\left[(5 + 2 m)^2 \right] \beta_2 +
\left[(5 + 2 m)^2 \right] \beta_3 +\\
\left[(5 + 2 m)^2  \right] \beta_4 +
\left[(5 + 2 m)^2  \right] \beta_5 +
\left[ (5 + 2 m)^2 \right] \beta_6 +
\left[ (5 + 2 m)^2 \right] \beta_7 +
\left[ (5 + 2 m)^2 \right] \beta_8 =0,
\end{align*}

Equation 2:

\begin{align*}
& \big[ -2 (5 + 2 m) (7 + 2 m) (97615125 +
2 m (21 + 4 m) (10060700 +
m (21 + 4 m) (744557 \\
& +
4 m (21 + 4 m) (5927 + 69 m (21 + 4 m))))) \big] \alpha_0 + \\
&	\big[ -2 (3 + 2 m) (5 + 2 m)^2 (7 + 2 m) (1722825 +
2 m (5 + m) (792789 + 32 m (5 + m) (6679 \\
& + 552 m (5 + m)))) \big] \alpha_1 + \\
&	\big [ -2 (3 + 2 m) (5 + 2 m) (28653975 +
2 m (19 + 4 m) (4255230 +
m (19 + 4 m) (430337 + \\
4 m (19 + 4 m) (4547 + 69 m (19 + 4 m))))) \big] \alpha_2 + \\
&	\left[ 7 (3 + 16 m (5 + m)) \right] \beta_0 +
\left[49 (13 + 4 m (5 + m))  \right] \beta_1 +
\left[2 (5 + 2 m)^2 \right] \beta_2 +
\left[ 3 (5 + 2 m)^2 \right] \beta_3 +\\
\left[4 (5 + 2 m)^2  \right] \beta_4 +
\left[ 5 (5 + 2 m)^2\right] \beta_5 +
\left[ 6 (5 + 2 m)^2 \right] \beta_6 +
\left[ 7 (5 + 2 m)^2  \right] \beta_7 +
\left[ 8 (5 + 2 m)^2 \right] \beta_8 =0,
\end{align*}

Equation 3:

\begin{align*}
& \big[ -7 (5 + 2 m) (7 + 2 m) (-4523755 +
8 m (21 + 4 m) (75805 +
m (21 + 4 m) (16433 + 4 m (21 + 4 m)\cdot \\
& \cdot (200 + 3 m (21 + 4 m))))) \big] \alpha_0 + \\
&	\big[ -7 (3 + 2 m) (7 + 2 m) (-6556875 +
8 m (5 + m) \cdot \\
& \cdot (-96085 +
4 m (5 + m) (36407 + 32 m (5 + m) (325 + 24 m (5 + m))))) \big] \alpha_1 + \\
&	\big [ -7 (3 + 2 m) (5 + 2 m) (-5009355 +
8 m (19 + 4 m) (-34525 +
m (19 + 4 m) (6233 + 4 m (19 + 4 m)\cdot \\
& \cdot  (140 + 3 m (19 + 4 m))))) \big] \alpha_2 + \\
&	\left[-294  \right] \beta_0 +
\left[ -80 + 96 m (5 + m) \right] \beta_1 +
\left[536 + 180 m (5 + m) \right] \beta_2 +
\left[-26 - 4 m (5 + m) \right] \beta_3 +\\
\left[49 + 8 m (5 + m)  \right] \beta_4 +
\left[ 149 + 24 m (5 + m) \right] \beta_5 +
\left[274 + 44 m (5 + m)  \right] \beta_6  \\
&  +
\left[ 424 + 68 m (5 + m) \right] \beta_7 +
\left[599 + 96 m (5 + m)  \right] \beta_8 =0,
\end{align*}

Equation 4:

\begin{align*}
& \big[-(5 + 2 m) (7 + 2 m) (-91051067 +
4 m (21 + 4 m) (-3125231 +
2 m (21 + 4 m) (-60393 \\
& + 16 m (21 + 4 m) (-15 + m (21 + 4 m)))))  \big] \alpha_0 + \\
&	\big[ -(3 + 2 m) (7 + 2 m) (-58593875 +
8 m (5 + m) (-5015259 +
64 m (5 + m) (-14847 \\
& + 8 m (5 + m) (-55 + 8 m (5 + m))))) \big] \alpha_1 + \\
&	\big [ -(3 + 2 m) (5 + 2 m) (-40305047 +
4 m (19 + 4 m) (-1969371 +
2 m (19 + 4 m) (-54393 \\
& + 16 m (19 + 4 m) (-35 + m (19 + 4 m))))) \big] \alpha_2 + \\
&	\left[-56  \right] \beta_0 +
\left[ -294 \right] \beta_1 +
\left[ -181 + 80 m (5 + m) \right] \beta_2 +
\left[334 + 148 m (5 + m) \right] \beta_3 +\\
\left[ -16 (19 + 3 m (5 + m)) \right] \beta_4 +
\left[ -5 (51 + 8 m (5 + m)) \right] \beta_5 +
\left[ -2 (53 + 8 m (5 + m)) \right] \beta_6  \\
& +
\left[28 (2 + m) (3 + m)  \right] \beta_7 +
\left[ 592 + 96 m (5 + m) \right] \beta_8 =0,
\end{align*}

Equation 5:

\begin{align*}
& \big[ 7 (5 + 2 m) (7 + 2 m) (963731 +
16 m (21 + 4 m) (27348 + 5 m (21 + 4 m) (401 + 8 m (21 + 4 m)))) \big] \alpha_0 + \\
&	\big[ 7 (3 + 2 m) (7 + 2 m) (-319565 +
16 m (5 + m) (65057 + 80 m (5 + m) (335 + 32 m (5 + m)))) \big] \alpha_1 + \\
&	\big [ 7 (3 + 2 m) (5 + 2 m) (-502109 +
16 m (19 + 4 m) (10298 + 5 m (19 + 4 m) (281 + 8 m (19 + 4 m)))) \big] \alpha_2 + \\
\left[-52  \right] \beta_1 +
\left[ -290 \right] \beta_2 +
\left[ -278 + 64 m (5 + m)\right] \beta_3 +\\
\left[ 5 (7 + 20 m (5 + m)) \right] \beta_4 +
\left[ -881 - 140 m (5 + m) \right] \beta_5 +
\left[-4 (284 + 45 m (5 + m))  \right] \beta_6 \\
\left[ -2 (621 + 98 m (5 + m)) \right] \beta_7 +
\left[ -6 (179 + 28 m (5 + m)) \right] \beta_8 =0,
\end{align*}

Equation 6:

\begin{align*}
& \big[ -(-5 - 2 m) (7 + 2 m) (-12838089 +
16 m (21 + 4 m) (-63803 + 2 m (21 + 4 m) (-227 \\
&+ 16 m (21 + 4 m)))) \big] \alpha_0 + \\
&	\big[ (3 + 2 m) (7 + 2 m) (-10661929 +
32 m (5 + m) (-125789 + 16 m (5 + m) (-359 + 64 m (5 + m)))) \big] \alpha_1 + \\
&	\big [ (3 + 2 m) (5 + 2 m) (-7979449 +
16 m (19 + 4 m) (-56863 + 2 m (19 + 4 m) (-467 + 16 m (19 + 4 m)))) \big] \alpha_2 + \\
\left[-48 \right] \beta_2 +
\left[-282 \right] \beta_3 +
\left[-367 + 48 m (5 + m)  \right] \beta_4 +
\left[ -353 + 36 m (5 + m) \right] \beta_5  \\
\left[-2 (923 + 148 m (5 + m))  \right] \beta_6
\left[-14 (213 + 34 m (5 + m)) \right] \beta_7 +
\left[ -96 (44 + 7 m (5 + m)) \right] \beta_8 =0,
\end{align*}

Equation 7:

\begin{align*}
& \big[ -98 (5 + 2 m) (7 + 2 m) (18097 +
16 m (21 + 4 m) (191 + 6 m (21 + 4 m))) \big] \alpha_0 + \\
&	\big[ -98 (3 + 2 m) (7 + 2 m) (10641 + 32 m (5 + m) (325 + 48 m (5 + m))) \big] \alpha_1 + \\
&	\big [-98 (3 + 2 m) (5 + 2 m) (5217 +
16 m (19 + 4 m) (131 + 6 m (19 + 4 m)))  \big] \alpha_2 + \\
\left[-44 \right] \beta_3 +
\left[ -270 \right] \beta_4 +
\left[ -444 + 32 m (5 + m) \right] \beta_5 +
\left[ -818 - 44 m (5 + m) \right] \beta_6 \\
\left[ -28 (117 + 19 m (5 + m)) \right] \beta_7 +
\left[ -42 (149 + 24 m (5 + m)) \right] \beta_8 =0,
\end{align*}

Equation 8:

\begin{align*}
& \big[ -16 (5 + 2 m) (7 + 2 m) (-39241 +
8 m (21 + 4 m) (-103 + 6 m (21 + 4 m))) \big] \alpha_0 + \\
&	\big[ -16 (3 + 2 m) (7 + 2 m) (-38265 + 16 m (5 + m) (-263 + 48 m (5 + m))) \big] \alpha_1 + \\
&	\big [ -16 (3 + 2 m) (5 + 2 m) (-33921 +
8 m (19 + 4 m) (-163 + 6 m (19 + 4 m))) \big] \alpha_2 + \\
\left[-40  \right] \beta_4 +
\left[ -254 \right] \beta_5 +
\left[-505 + 16 m (5 + m)  \right] \beta_6 +
\left[ -28 (48 + 5 m (5 + m)) \right] \beta_7 \\
& +
\left[ -48 (109 + 18 m (5 + m)) \right] \beta_8 =0,
\end{align*}

Equation 9:

\begin{align*}
& \big[ 672 (5 + 2 m) (7 + 2 m) (179 + 12 m (21 + 4 m)) \big] \alpha_0 + \\
&	\big[672 (3 + 2 m) (7 + 2 m) (155 + 48 m (5 + m))  \big] \alpha_1 + \\
&	\big [ 672 (3 + 2 m) (5 + 2 m) (119 + 12 m (19 + 4 m)) \big] \alpha_2 + \\
& \left[ -36 \right] \beta_5 +
\left[ -234 \right] \beta_6 +
\left[ -546 \right] \beta_7 +
\left[ -21 (91 + 12 m (5 + m)) \right] \beta_8 =0,
\end{align*}

Equation 10:

\begin{align*}
& \big[16 (5 + 2 m) (7 + 2 m) (-421 + 32 m (21 + 4 m))  \big] \alpha_0 + \\
&	\big[ 16 (3 + 2 m) (7 + 2 m) (-485 + 128 m (5 + m)) \big] \alpha_1 + \\
&	\big [ 16 (3 + 2 m) (5 + 2 m) (-581 + 32 m (19 + 4 m)) \big] \alpha_2 + \\
\left[ -32 \right] \beta_6 +
\left[ -210 \right] \beta_7 +
\left[ -563 - 16 m (5 + m) \right] \beta_8 =0,
\end{align*}

Equation 11:

\begin{align*}
& \big[-2464 (5 + 2 m) (7 + 2 m)  \big] \alpha_0 +
\big[-2464 (3 + 2 m) (7 + 2 m)  \big] \alpha_1 + \\
&	\big [ -2464 (3 + 2 m) (5 + 2 m) \big] \alpha_2 +
\left[-28  \right] \beta_7 +
\left[ -182 \right] \beta_8 =0,
\end{align*}

Equation 12:

\begin{align*}
& \big[-128 (5 + 2 m) (7 + 2 m)  \big] \alpha_0 +
\big[-128 (3 + 2 m) (7 + 2 m)  \big] \alpha_1 + \\
&	\big [ -128 (3 + 2 m) (5 + 2 m) \big] \alpha_2 +
\left[-24  \right] \beta_8 =0.
\end{align*}

We note that, this system has the matrix form

\begin{align*}
\left[\begin{array}{@{}ccc|cccccccccc@{}}
x_{00} & x_{01} & x_{02} & y_{00} & y_{01}  & y_{02} & y_{03}  & y_{04} & y_{05}  & y_{06}  & y_{07} & y_{08}  \\
x_{10} & x_{11} & x_{12} & y_{10} & y_{11}  & y_{12} & y_{13}  & y_{14} & y_{15}  & y_{16}  & y_{17} & y_{18} \\
x_{20} & x_{21} & x_{22} & y_{20} & y_{21}  & y_{22} & y_{23}  & y_{24} & y_{25}  & y_{26}  & y_{27} & y_{28} \\
x_{30} & x_{31} & x_{32} & y_{30} & y_{31}  & y_{32} & y_{33}  & y_{34} & y_{35}  & y_{36}  & y_{37} & y_{38} \\
x_{40} & x_{41} & x_{42} &0 & y_{41}  & y_{42} & y_{43}  & y_{44} & y_{45}  & y_{46}  & y_{47} & y_{48} \\
x_{50} & x_{51} & x_{52} & 0 & 0  & y_{52} & y_{53}  & y_{54} & y_{55}  & y_{56}  & y_{57} & y_{58} \\
x_{60} & x_{61} & x_{62} & 0 & 0 &0 & y_{63}  & y_{64} & y_{65}  & y_{66}  & y_{67} & y_{68} \\
x_{70} & x_{71} & x_{72} & 0 & 0  & 0 & 0  & y_{74} & y_{75}  & y_{76}  & y_{77} & y_{78} \\
x_{80} & x_{81} & x_{82} &0 & 0   & 0  & 0   & 0 & y_{85}  & y_{86}  & y_{87} & y_{88} \\
x_{90} & x_{91} & x_{92} & 0 & 0  &0 & 0  & 0 & 0  & y_{96}  & y_{97} & y_{98} \\
x_{10,0} & x_{10,1} & x_{10,2} & 0 & 0  &0 & 0  &0 & 0 & 0  & y_{10,7} & y_{10,8} \\
x_{11,0} & x_{11,1} & x_{11,2} & 0 & 0  &0 & 0  & 0 & 0  &0  &0 & y_{11,8}
\end{array}\right]
\begin{bmatrix}
\alpha_0  \\
\alpha_1 \\
\alpha_2 \\
\hline
\beta_0 \\
\beta_1 \\
\beta_2 \\
\beta_3 \\
\beta_4 \\
\beta_5 \\
\beta_6 \\
\beta_7 \\
\beta_8
\end{bmatrix} =
\begin{bmatrix}
0  \\
0  \\0  \\0  \\0  \\0  \\0  \\0  \\0  \\0  \\0  \\0
\end{bmatrix} ,
\end{align*}

where the coefficients of the matrix are given as above. One can solve this linear system with respect to the unknowns  $\{\alpha_0,\alpha_1,\alpha_2,\beta_0,\dots,\beta_6\}$ so that $\alpha_0$ and $\beta_0$ are free variables. Then, we choose the consistent values

\begin{align}\label{Choiceofbeta0}
\beta_0=0
\end{align}

\begin{align}\label{Definitionalpha0}
\alpha_0 = 2 (1 + m)^2 (3 + m) (1 + 2 m) (3 + 2 m)^2 (5 + 2 m) (6 + 2 m) (34 +
24 m + 4 m^2).
\end{align}

With these choices, one can solve the linear system above, obtaining the solutions

\begin{align}
\alpha_1 &=  -8 (5 + 2 m)^4 (819 + 4 m (5 + m) (106 + m (5 + m) (18 + m (5 + m)))),\label{Definitionalpha1} \\
\alpha_2 &= 8 (2 + m)^2 (4 + m)^2 (5 + 2 m) (7 + 2 m)^2 (9 + 2 m) (7 +
2 m (4 + m)),   \label{Definitionalpha2}
\end{align}

as well as

\begin{align*}
\beta_1 &=-240 (2 + m) (3 + m) (3 + 2 m) (5 + 2 m)^2 (7 + 2 m) (424785 +
2 m (5 + m) (159981 \\
& +
4 m (5 + m) (11215 + 4 m (5 + m) (347 + 16 m (5 + m))))), \\
\beta_2 &=-32 (2 + m) (3 + m) (3 + 2 m) (5 + 2 m)^2 (7 + 2 m) (5795205 +
2 m (5 + m) (2262318 \\
m (5 + m) (655655 + 8 m (5 + m) (10457 + 496 m (5 + m))))), \\
\beta_3 & = -80 (2 + m) (3 + m) (3 + 2 m) (5 + 2 m)^2 (7 + 2 m) (937797 +
4 m (5 + m) (214330 \\
m (5 + m) (70297 + 64 m (5 + m) (155 + 8 m (5 + m))))), \\
\beta_4 &=-32 (2 + m) (3 + m) (3 + 2 m) (5 + 2 m)^2 (7 + 2 m) (-732039 +
m (5 + m) (-252965 \\
8 m (5 + m) (1543 + 8 m (5 + m) (185 + 16 m (5 + m))))), \\
\beta_5 &= 320 (2 + m) (3 + m) (3 + 2 m) (5 + 2 m)^2 (7 + 2 m) (49665 +
4 m (5 + m) (7427 \\
&+ m (5 + m) (1469 + 96 m (5 + m)))), \\
\beta_6 &= 128 (2 + m) (3 + m) (3 + 2 m) (5 + 2 m)^2 (7 + 2 m) (3255 +
2 m (5 + m) (2443 \\
&+ 2 m (5 + m) (373 + 32 m (5 + m)))), \\
\beta_7 & = -5120 (2 + m)^2 (3 + m)^2 (3 + 2 m) (5 + 2 m)^2 (7 + 2 m) (19 +
4 m (5 + m)) , \\
\beta_8 &=  -512 (2 + m)^2 (3 + m)^2 (3 + 2 m) (5 + 2 m)^2 (7 + 2 m) (19 +
4 m (5 + m))  .
\end{align*}

On the one hand, using (<ref>), we compute the left-hand-side of the recurrence relation (<ref>), that is

\begin{align}\label{Defintitionoftk}
t_{k} &= \frac{6720}{\pi ^2} \frac{(m+3) (2 m+5) (m+2) }{
\prod_{i=2}^{5} (2m-k+i)
\prod_{j=5}^{8} (4m+2k+2j+1) }\frac{(k+1) (4 k+7)}{(2 k+1) (2 k+3)} \cdot \nonumber \\
& \cdot \frac{\Gamma \left(k-\frac{1}{2}\right) \Gamma \left(k+\frac{7}{2}\right)}{\Gamma (k+4) \Gamma (k+5)}\frac{\Gamma \left(-k+2 m+\frac{3}{2}\right) \Gamma (k+2 m+5)}{\Gamma (-k+2 m+2) \Gamma \left(k+2 m+\frac{11}{2}\right)}T_m(k) ,
\end{align}

where $T_m(k)$ is a polynomial with respect to $k$ of degree $8$,

\begin{align*}
T_m(k) &= \sum_{j=0}^{8} \tilde{T}_m (j) k^j,
\end{align*}

and the coefficients are given by

\begin{align*}
\tilde{T}_m (0)  &=  2 (2 + m) (3 + m) (5 + 2 m)^2 (24057 +
2 m (5 + m) (8121 + 8 m (5 + m) (223 + 16 m (5 + m))))  ,\\
\tilde{T}_m (1)  &= 21 (2535165 +
m (5 + m) (2762223 +
2 m (5 + m) (588501 +
4 m (5 + m) (30749 + 8 m (5 + m)\cdot \\
& \cdot (395 + 16 m (5 + m))))))  ,\\
\tilde{T}_m (2)  &= -36879489 +
2 m (5 + m) (-12352033 +
m (5 + m) (-2533381 +
12 m (5 + m)\cdot \\
& \cdot  (-3819 + 16 m (5 + m)  (199 + 16 m (5 + m)))))   ,\\
\tilde{T}_m (3)  &=  -42 (406041 +
2 m (5 + m) (195702 +
m (5 + m) (66649 + 4 m (5 + m) (2421 + 128 m (5 + m))))) ,\\
\tilde{T}_m (4)  &=  5822802 -
m (5 + m) (-2447647 +
4 m (5 + m) (-30745 + 48 m (5 + m) (299 + 32 m (5 + m))))  ,\\
\tilde{T}_m (5)  &=  504 (4315 + m (5 + m) (2771 + 2 m (5 + m) (291 + 20 m (5 + m)))) ,\\
\tilde{T}_m (6)  &=  8 (-13212 + m (5 + m) (1877 + 40 m (5 + m) (53 + 6 m (5 + m)))) ,\\
\tilde{T}_m (7)  &=-672 (2 + m) (3 + m) (19 + 4 m (5 + m))   ,\\
\tilde{T}_m (8)  &=  -48 (2 + m) (3 + m) (19 + 4 m (5 + m))  .
\end{align*}

On the other hand, using (<ref>) together with the definitions of of $t_k$ from (<ref>), $p_3(k)$ from (<ref>), $p(k)$ from (<ref>) and $b(k)$ from (<ref>), we find $G=G (m,k)$ in the form of (<ref>), that is

\begin{align*}
G (m,k)= R(m,k) \mathsf{F}(m,k),
\end{align*}

for some function $R(m,k)$ which we compute explicitly. This has the form

\begin{align*}
R(m,k)=k \tilde{R}(m,k),
\end{align*}

for some complicated function $ \tilde{R}(m,k)$ which we leave out. From the latter, we get immediately $R(m,0)=0$ and hence[The fact that $G(m,0)=0$, for all integers $m \geq 1$, also follows from the choice (<ref>). Indeed, according to (<ref>), we have that $\beta_0=0$ and hence, based on the definition (<ref>), we get  $b(0)=\beta_0=0$. Then, (<ref>) yields $G(m,0)=0$. ]

\begin{align*}
\end{align*}

Since $K=2m+1$, we also compute $\mathsf{F}_{m} (2m+2)=0$ and hence we get

\begin{align*}
G_{m} (2m+2)=0 .
\end{align*}

Putting these together, we infer

\begin{align*}
G (m,2m+2)-G (m,0)=0.
\end{align*}

Consequently, the right-hand-side of (<ref>) vanishes and the desired recurrence relation for $\mathsf{f}(m) $ boils down to

\begin{align*}
\alpha_{0}(m)\mathsf{f}_{m} +\alpha_{1}(m)\mathsf{f}_{m+1}+\alpha_{2}(m)\mathsf{f}_{m+2}=0,
\end{align*}

where the coefficients $\{\alpha_{0},\alpha_{1},\alpha_{2}\}$ are given by (<ref>), (<ref>) and (<ref>) respectively. Finally, we note that these coeffiicents coincide with the corresponding coefficients $\mathsf{P}_m$, $-\mathsf{Q}_m$ and $\mathsf{R}_m$ from Lemma <ref> for $\delta=3$.

§ AUXILIARY COMPUTATIONS FOR THE KG

Fix a real number $\delta  \geq 3/2$ and let $\mathsf{R}_m $, $\mathsf{A}_m$ and $\mathsf{B}_m$ be the functions defined in Lemma <ref>. Then,

\begin{align*}
\mathsf{R}_m > 0, \quad \mathsf{B}_m >0, \quad \mathsf{A}_m -\mathsf{B}_m <1,
\end{align*}

for all integers $m \geq 1$

Fix a real number $\delta  \geq 3/2$. Firstly, notice that $\mathsf{R}_m$ is given in a factorized form. From this, it is clear that $\mathsf{R}_m >0$, for all integers $m \geq 1$ and real $\delta \geq 3/2$. Then,

\begin{align*}
\mathsf{A}_m -\mathsf{B}_m -1 <0 \Longleftrightarrow \mathsf{R}_m (\mathsf{A}_m -\mathsf{B}_m -1) <0 \Longleftrightarrow  \mathsf{Q}_m -\mathsf{P}_m -\mathsf{R}_m <0.
\end{align*}

Now, a long computation shows that

\begin{align*}
\mathsf{Q}_m -\mathsf{P}_m -\mathsf{R}_m = \sum_{j=0}^{8} \mathsf{c}_{j}  m^j,
\end{align*}

where the coefficients are given by

\begin{align*}
\mathsf{c}_0 &= -85 \delta ^7-1532 \delta ^6-6202 \delta ^5-10288 \delta ^4-7401 \delta ^3-2196 \delta ^2-64 \delta +120   ,\\
\mathsf{c}_1 &=  -108 \delta ^7-3140 \delta ^6-21004 \delta ^5-51892 \delta ^4-54816 \delta ^3-25736 \delta ^2-4760 \delta +176  ,\\
\mathsf{c}_2 &=  -32 \delta ^7-2092 \delta ^6-24908 \delta ^5-96308 \delta ^4-149484 \delta ^3-101032 \delta ^2-29816 \delta -2600  ,\\
\mathsf{c}_3 &= -448 \delta ^6-12416 \delta ^5-83392 \delta ^4-198400 \delta ^3-191808 \delta ^2-78784 \delta -11904   ,\\
\mathsf{c}_4 &= -2208 \delta ^5-33984 \delta ^4-137568 \delta ^3-198144 \delta ^2-111136 \delta -22688   ,\\
\mathsf{c}_5 &= -5248 \delta ^4-47744 \delta ^3-113664 \delta ^2-89984 \delta -23296   ,\\
\mathsf{c}_6 &= -6528 \delta ^3-33920 \delta ^2-41728 \delta -13440   ,\\
\mathsf{c}_7 &=  -4096 \delta ^2-10240 \delta -4096  ,\\
\mathsf{c}_8 &=  -1024 \delta -512  .
\end{align*}

One can easily see that $\mathsf{c}_j<0$, for all $j \in \{0,1,\dots,8\}$ and all real $\delta \geq 3/2$, that proves the claim.

Fix a real number $\delta  \geq 3/2$ and let $\mathsf{x}_1$ be the initial value defined in Lemma <ref>. Then, $\mathsf{x}_1<0$.

Fix a real number $\delta  \geq 3/2$ and recall that

\begin{align*}
\mathsf{x}_1 &  = \frac{63 \delta ^5+196 \delta ^4+111 \delta ^3-78 \delta ^2-52 \delta -24}{80 \delta ^5+448 \delta ^4+492 \delta ^3-136 \delta ^2-236 \delta -48}.
\end{align*}

One can easily show that the denominator is strictly positive for all real $\delta \geq 3/2$. Then, a direct computation shows that $\mathsf{x}_1 <1$ is equivalent to

\begin{align*}
% x_1 <1 & \Longleftrightarrow 63 \delta ^5+196 \delta ^4+111 \delta ^3-78 \delta ^2-52 \delta -24 < 80 \delta ^5+448 \delta ^4+492 \delta ^3-136 \delta ^2-236 \delta -48& \Longleftrightarrow   \\
-17 \delta ^5-252 \delta ^4-381 \delta ^3+58 \delta ^2+184 \delta +24 <0.
\end{align*}

Finally, a straight forward computation verifies the latter and completes the proof.

§ AUXILIARY COMPUTATIONS FOR THE WM

Fix a real numbabusier $\delta  >0$ and let $\mathfrak{R}_m $, $\mathfrak{A}$ and $\mathfrak{B}$ be the functions defined in Lemma <ref>. Then,

\begin{align*}
\mathfrak{R}_m >0, \quad \mathfrak{B}_m >0, \quad \mathfrak{A}_m- \mathfrak{B}_m <1,
\end{align*}

for all integers $m \geq 1$.

Fix a real number $\delta  >0$. Firstly, notice that $\mathfrak{R}_m$ is given in a factorized form. From this, it is clear that $\mathfrak{R}_m >0$, for all integers $m \geq 1$ and real $\delta  >0$. Then,

\begin{align*}
\mathfrak{A}_m -\mathfrak{B}_m -1 <0 \Longleftrightarrow \mathfrak{R}_m (\mathfrak{A}_m -\mathfrak{B}_m -1) <0 \Longleftrightarrow  \mathfrak{Q}_m -\mathfrak{R}_m -\mathfrak{R}_m <0.
\end{align*}

Now, a long computation shows that

\begin{align*}
\mathfrak{Q}_m -\mathfrak{R}_m -\mathfrak{R}_m &=
\frac{-\mathfrak{c}(m)}{(m+2) (\delta +m+2) (\delta +m+3) (\delta +2 m+3) (\delta +2 m+6) }\cdot \\
\frac{1}{(2 \delta +2 m+7) (\delta +2 m (\delta +m+3)+3)} ,
\end{align*}

where $\mathfrak{c}(m)$ is a polynomial with respect to $m$ of degree $6$,

\begin{align*}
\mathfrak{c}(m)=\sum_{j=0}^{6}\mathfrak{c}_j m^j,
\end{align*}

and its coefficients are given by

\begin{align*}
\mathsf{c}_0 &= \delta ^6+29 \delta ^5+291 \delta ^4+1379 \delta ^3+3260 \delta ^2+3472 \delta +1072  ,\\
\mathsf{c}_1 &=  2 \delta ^6+68 \delta ^5+784 \delta ^4+4200 \delta ^3+11054 \delta ^2+12868 \delta +4144   ,\\
\mathsf{c}_2 &=  26 \delta ^5+544 \delta ^4+4044 \delta ^3+13256 \delta ^2+17998 \delta +6156  ,\\
\mathsf{c}_3 &=  112 \delta ^4+1576 \delta ^3+7352 \delta ^2+12512 \delta +4608 ,\\
\mathsf{c}_4 &= 216 \delta ^3+1924 \delta ^2+4620 \delta +1856   ,\\
\mathsf{c}_5 &= 192 \delta ^2+864 \delta +384   ,\\
\mathsf{c}_6 &= 64 \delta +32 .
\end{align*}

One can easily see that $\mathfrak{c}_j>0$, for all $j \in \{0,1,\dots,6\}$ and all real $\delta  >0$, that proves the claim.

Fix a real number $\delta  >0$	and let $\mathfrak{x}_1 $  be the initial value defined in Lemma  <ref>. Then, $ \mathfrak{x}_1<1$.

Fix a real number $\delta  >0$ and recall that

\begin{align*}
\mathfrak{x}_1 &  = \frac{3 \delta ^5+38 \delta ^4+189 \delta ^3+482 \delta ^2+708 \delta +560}{4 \delta ^5+62 \delta ^4+364 \delta ^3+1038 \delta ^2+1472 \delta +840}.
\end{align*}

One can easily show that the denominator is strictly positive for all real $\delta >0$. Then, a direct computation shows that $\mathfrak{x}_1 <1$ is equivalent to

\begin{align*}
-\delta ^5-24 \delta ^4-175 \delta ^3-556 \delta ^2-764 \delta -280<0.
\end{align*}

Finally, a straight forward computation verifies the latter and completes the proof.

[1]
J. Ablinger, J. Blümlein, C. G. Raab, and C. Schneider.
Iterated binomial sums and their associated iterated integrals.
J. Math. Phys., 55(11):112301, 57, 2014.

[2]
S. A. Abramov.
When does Zeilberger's algorithm succeed?
Adv. in Appl. Math., 30(3):424–441, 2003.

[3]
George E. Andrews and Peter Paule.
Some questions concerning computer-generated proofs of a binomial
double-sum identity.
Journal of Symbolic Computation, 16(2):147–151, 1993.

[4]
Venkat Balasubramanian, Alex Buchel, Stephen R. Green, Luis Lehner, and
Steven L. Liebling.
Holographic thermalization, stability of anti–de sitter space, and
the fermi-pasta-ulam paradox.
Phys. Rev. Lett., 113:071601, Aug 2014.

[5]
D. Bambusi and N.N. Nekhoroshev.
A property of exponential stability in nonlinear wave equations near
the fundamental linear mode.
Physica D: Nonlinear Phenomena, 122(1):73–104, 1998.

[6]
D. Bambusi and S. Paleari.
Families of periodic solutions of resonant PDEs.
J. Nonlinear Sci., 11(1):69–87, 2001.

[7]
Hans Bantilan, Pau Figueras, Markus Kunesch, and Paul Romatschke.
Nonspherically symmetric collapse in asymptotically ads spacetimes.
Phys. Rev. Lett., 119:191103, Nov 2017.

[8]
Farid Bencherif, Rachid Boumahdi, and Tarek Garici.
Symmetric identity for polynomial sequences satisfying
Commun. Math., 29(3):343–355, 2021.

[9]
Massimiliano Berti and Philippe Bolle.
Periodic solutions of nonlinear wave equations with general
Comm. Math. Phys., 243(2):315–328, 2003.

[10]
Anxo Biasi, Ben Craps, and Oleg Evnin.
Energy returns in global $\rm AdS_4$.
Phys. Rev. D, 100(2):024008, 17, 2019.

[11]
Piotr Bizo ńń and Joanna
Jałmu żżna.
Globally regular instability of 3-dimensional anti–de sitter
Phys. Rev. Lett., 111:041102, Jul 2013.

[12]
Piotr Bizo ńń, Maciej Maliborski, and Andrzej
Resonant dynamics and the instability of anti–de sitter spacetime.
Phys. Rev. Lett., 115:081103, Aug 2015.

[13]
Piotr Bizo ńń and Andrzej Rostworowski.
Weakly turbulent instability of anti–de sitter spacetime.
Phys. Rev. Lett., 107:031102, Jul 2011.

[14]
Piotr Bizoń.
Is AdS stable?
Gen. Relativity Gravitation, 46(5):Art. 1724, 14, 2014.

[15]
Piotr Bizoń and Andrzej Rostworowski.
Instability of anti-de Sitter spacetime.
In General relativity, cosmology and astrophysics, volume 177
of Fundam. Theor. Phys., pages 239–244. Springer, Cham, 2014.

[16]
Piotr Bizoń and Andrzej Rostworowski.
Gravitational turbulent instability of $\rm AdS_5$.
Acta Phys. Polon. B, 48(8):1375–1381, 2017.

[17]
John M. Campbell.
Applications of Zeilberger's algorithm to Ramanujan-inspired
series involving harmonic-type numbers.
Discrete Math. Lett., 11:7–13, 2023.

[18]
Athanasios Chatzikaleas.
On the Fourier analysis of the Einstein-Klein-Gordon system: Growth
and Decay of the Fourier constants.
arXiv e-prints, page arXiv:2004.11049, April 2020.

[19]
Athanasios Chatzikaleas and Jacques Smulevici.
Non-linear periodic waves on the Einstein cylinder.
arXiv e-prints, page arXiv:2201.05447, January 2022.

[20]
Frédéric Chyzak.
An extension of Zeilberger's fast algorithm to general holonomic
volume 217, pages 115–134. 2000.
Formal power series and algebraic combinatorics (Vienna, 1997).

[21]
Ben Craps and Oleg Evnin.
AdS (in)stability: an analytic approach.
Fortschr. Phys., 64(4-5):336–344, 2016.

[22]
Ben Craps, Oleg Evnin, Puttarak Jai-akson, and Joris Vanhoof.
Ultraviolet asymptotics for quasiperiodic $\rm Ads_4$
J. High Energy Phys., (10):080, front matter+18, 2015.

[23]
Ben Craps, Oleg Evnin, and Joris Vanhoof.
Renormalization group, secular term resummation and AdS
J. High Energy Phys., (10):048, front matter+30, 2014.

[24]
Ben Craps, Oleg Evnin, and Joris Vanhoof.
Ultraviolet asymptotics and singular dynamics of AdS
J. High Energy Phys., (10):079, front matter+16, 2015.

[25]
M Dafermos.
The black hole stability problem.
Talk at the Newton Institute, Cambridge, http://www-
old.newton.ac.uk/webseminars/pg+ws/2006/gmx/1010/dafermos/, 2006.

[26]
M Dafermos and G Holzegel.
Dynamic instability of solitons in $4+1$ dimesnional gravity with
negative cosmological constant.
https://www.dpmms.cam.ac.uk/ md384/ADSinstability.pdf, 2006.

[27]
Nils Deppe, Allison Kolly, Andrew Frey, and Gabor Kunstatter.
Stability of anti–de sitter space in einstein-gauss-bonnet gravity.
Phys. Rev. Lett., 114:071102, Feb 2015.

[28]
Óscar J. C. Dias, Gary T. Horowitz, and Jorge E. Santos.
Gravitational turbulent instability of anti-de Sitter space.
Classical Quantum Gravity, 29(19):194002, 7, 2012.

[29]
Oleg Evnin and Puttarak Jai-akson.
Detailed ultraviolet asymptotics for AdS scalar field
J. High Energy Phys., (4):054, front matter+19, 2016.

[30]
Stephen R. Green, Antoine Maillard, Luis Lehner, and Steven L. Liebling.
Islands of stability and recurrence times in ads.
Phys. Rev. D, 92:084001, Oct 2015.

[31]
Emmanuel Hebey.
Nonlinear analysis on manifolds: Sobolev spaces and
inequalities, volume 5 of Courant Lecture Notes in Mathematics.
New York University, Courant Institute of Mathematical Sciences, New
York; American Mathematical Society, Providence, RI, 1999.

[32]
Dominika Hunik-Kostyra and Andrzej Rostworowski.
AdS instability: resonant system for gravitational perturbations of
\(\mathrm{AdS}_5\) in the cohomogeneity-two biaxial Bianchi IX
J. High Energy Phys., 2020(6):40, 2020.
Id/No 2.

[33]
Wolfram Research, Inc.
Mathematica, Version 12.2.

[34]
Wolfram Koepf.
Hypergeometric summation.
Advanced Lectures in Mathematics. Friedr. Vieweg & Sohn,
Braunschweig, 1998.
An algorithmic approach to summation and special function identities.

[35]
Alois Kufner.
Weighted Sobolev spaces.
A Wiley-Interscience Publication. John Wiley & Sons, Inc., New York,
Translated from the Czech.

[36]
Maciej Maliborski and Andrzej Rostworowski.
Time-periodic solutions in an einstein ads–massless-scalar-field
Phys. Rev. Lett., 111:051102, Aug 2013.

[37]
Georgios Moschidis.
The Einstein–null dust system in spherical symmetry with an inner
mirror: structure of the maximal development and Cauchy stability.
arXiv e-prints, page arXiv:1704.08685, April 2017.

[38]
Georgios Moschidis.
A proof of the instability of AdS for the Einstein–massless Vlasov
arXiv e-prints, page arXiv:1812.04268, December 2018.

[39]
Georgios Moschidis.
Two Instability Results in General Relativity.
ProQuest LLC, Ann Arbor, MI, 2018.
Thesis (Ph.D.)–Princeton University.

[40]
Georgios Moschidis.
A proof of the instability of AdS for the Einstein-null dust
system with an inner mirror.
Anal. PDE, 13(6):1671–1754, 2020.

[41]
Georgios Moschidis.
A proof of the instability of AdS for the Einstein-massless Vlasov
Invent. math., 2022.

[42]
Georgios Moschidis.
The instability of anti-de sitter spacetime for the einstein-scalar
field system.
Princeton University, Seminar, Nov. 2021.

[43]
Frank W. J. Olver, Daniel W. Lozier, Ronald F. Boisvert, and Charles W. Clark,
NIST handbook of mathematical functions.
U.S. Department of Commerce, National Institute of Standards and
Technology, Washington, DC; Cambridge University Press, Cambridge, 2010.
With 1 CD-ROM (Windows, Macintosh and UNIX).

[44]
Peter Paule and Axel Riese.
A Mathematica $q$-analogue of Zeilberger's algorithm based on
an algebraically motivated approach to $q$-hypergeometric telescoping.
In Special functions, $q$-series and related topics
(Toronto, ON, 1995), volume 14 of Fields Inst. Commun., pages
179–210. Amer. Math. Soc., Providence, RI, 1997.

[45]
Peter Paule and Markus Schorn.
A Mathematica version of Zeilberger's algorithm for proving
binomial coefficient identities.
volume 20, pages 673–698. 1995.
Symbolic computation in combinatorics $\Delta_1$ (Ithaca, NY, 1993).

[46]
Peter Paule, Markus Schorn, and Axel Riese.
Fast zeilberger package, Version 3.61.
Copyright Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria.

[47]
Marko Petkovšek, Herbert S. Wilf, and Doron Zeilberger.
A K Peters, Ltd., Wellesley, MA, 1996.
With a foreword by Donald E. Knuth, With a separately available
computer disk.

[48]
J. Sánchez-Ruiz.
Linearization and connection formulae involving squares of
Gegenbauer polynomials.
Appl. Math. Lett., 14(3):261–267, 2001.

[49]
Gábor Szegő.
Orthogonal polynomials.
American Mathematical Society, Providence, R.I., fourth edition,
American Mathematical Society, Colloquium Publications, Vol. XXIII.

[50]
Bengt Ove Turesson.
Nonlinear potential theory and weighted Sobolev spaces,
volume 1736 of Lecture Notes in Mathematics.
Springer-Verlag, Berlin, 2000.

[51]
C. M. Warnick.
The massive wave equation in asymptotically AdS spacetimes.
Comm. Math. Phys., 321(1):85–111, 2013.

[52]
Kurt Wegschaider, Axel Riese, and Burkhard Zimmermann.
Multisum package, Version 2.3.
Copyright Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria.