[proofplan]
We identify $C_{p,q}^{\mathrm{reg}}$ as the cone over the product link $L_{p,q}=S^p(\sqrt{p/(p+q)})\times S^q(\sqrt{q/(p+q)})\subset \mathbb{S}^{p+q+1}$. The chosen radii make the two families of principal curvatures cancel, so the link is minimal in the unit sphere, and therefore the cone is minimal in Euclidean space. For $p=q=3$, the link has $|A_L|^2=6$, so the second variation quadratic form on the cone reduces to a radial Hardy inequality with weight $r^6$. The sharp one-dimensional Hardy constant is $25/4$, which dominates $6$, giving nonnegativity for every compactly supported variation away from the vertex.
[/proofplan]
[step:Describe the cone as the radial extension of its product link]
Set $m:=p+q$. Define
\begin{align*}
a&:=\sqrt{\frac{p}{m}},&
b&:=\sqrt{\frac{q}{m}},
\end{align*}
and define the product submanifold
\begin{align*}
L_{p,q}:=
S^p(a)\times S^q(b)
\subset
\mathbb{R}^{p+1}\times\mathbb{R}^{q+1}.
\end{align*}
Since $a^2+b^2=1$, we have $L_{p,q}\subset \mathbb{S}^{m+1}$. Define the smooth map $\Psi:(0,\infty)\times L_{p,q} \to C_{p,q}^{\mathrm{reg}}$ by
\begin{align*}
\Psi(r,z):=rz \quad \text{for } (r,z)\in (0,\infty)\times L_{p,q}.
\end{align*}
Then $\Psi$ is a smooth diffeomorphism.
Indeed, if $z=(u,v)\in L_{p,q}$, then
\begin{align*}
q|u|^2-p|v|^2
=
q\frac{p}{m}-p\frac{q}{m}
=
0,
\end{align*}
so $rz\in C_{p,q}^{\mathrm{reg}}$. Conversely, if $(x,y)\in C_{p,q}^{\mathrm{reg}}$, set
\begin{align*}
r:=\sqrt{|x|^2+|y|^2}.
\end{align*}
Then $r>0$ and $z:=(x/r,y/r)\in \mathbb{S}^{m+1}$. The defining equation gives
\begin{align*}
q|x|^2=p|y|^2,
\qquad
|x|^2+|y|^2=r^2,
\end{align*}
hence
\begin{align*}
|x/r|^2=\frac{p}{m},
\qquad
|y/r|^2=\frac{q}{m}.
\end{align*}
Thus $z\in L_{p,q}$ and $(x,y)=rz$.
[/step]
[step:Compute the principal curvatures of the link in the unit sphere]
Define the smooth unit normal field $\nu_L:L_{p,q}\to T\mathbb{S}^{m+1}\big|_{L_{p,q}}$ by
\begin{align*}
\nu_L(u,v):=\left(\frac{b}{a}u,-\frac{a}{b}v\right) \quad \text{for } (u,v)\in L_{p,q}.
\end{align*}
For $(u,v)\in L_{p,q}$,
\begin{align*}
|\nu_L(u,v)|^2
=
\frac{b^2}{a^2}|u|^2+\frac{a^2}{b^2}|v|^2
=
b^2+a^2
=
1,
\end{align*}
and
\begin{align*}
\nu_L(u,v)\cdot (u,v)
=
\frac{b}{a}|u|^2-\frac{a}{b}|v|^2
=
ab-ab
=
0.
\end{align*}
Thus $\nu_L$ is tangent to $\mathbb{S}^{m+1}$ and normal to $L_{p,q}$ inside $\mathbb{S}^{m+1}$.
Let $A_L$ denote the shape operator of $L_{p,q}\subset \mathbb{S}^{m+1}$ with respect to $\nu_L$, using the convention $A_L(W)=-D_W\nu_L$, where $D$ is the Euclidean derivative. If $W=(\xi,0)$ is tangent to the $S^p(a)$ factor, then
\begin{align*}
D_W\nu_L=\frac{b}{a}(\xi,0),
\qquad
A_L(W)=-\frac{b}{a}W.
\end{align*}
If $W=(0,\eta)$ is tangent to the $S^q(b)$ factor, then
\begin{align*}
D_W\nu_L=-\frac{a}{b}(0,\eta),
\qquad
A_L(W)=\frac{a}{b}W.
\end{align*}
Therefore the principal curvatures of $L_{p,q}$ in $\mathbb{S}^{m+1}$ are
\begin{align*}
-\frac{b}{a}=-\sqrt{\frac{q}{p}}
\quad\text{with multiplicity }p,
\qquad
\frac{a}{b}=\sqrt{\frac{p}{q}}
\quad\text{with multiplicity }q.
\end{align*}
The mean curvature scalar of $L_{p,q}$ in $\mathbb{S}^{m+1}$ is the trace of $A_L$:
\begin{align*}
\operatorname{tr}A_L
=
-p\sqrt{\frac{q}{p}}+q\sqrt{\frac{p}{q}}
=
-\sqrt{pq}+\sqrt{pq}
=
0.
\end{align*}
Thus $L_{p,q}$ is minimal in $\mathbb{S}^{m+1}$.
[/step]
[step:Pass minimality from the link to the cone]
Under the parametrization $\Psi:(0,\infty)\times L_{p,q}\to C_{p,q}^{\mathrm{reg}}$, the induced metric on $C_{p,q}^{\mathrm{reg}}$ is
\begin{align*}
g_C
=
d\mathcal{L}^1(r)^2+r^2 g_L,
\end{align*}
where $g_L$ is the induced metric on $L_{p,q}$. The radial direction has zero normal curvature, and each principal curvature $\kappa$ of $L_{p,q}\subset \mathbb{S}^{m+1}$ contributes the cone principal curvature $\kappa/r$ on $C_{p,q}^{\mathrm{reg}}$. Hence the Euclidean mean curvature scalar of the cone is
\begin{align*}
H_C(r,z)
=
\frac{1}{r}\operatorname{tr}A_L(z)
=
0
\end{align*}
for every $(r,z)\in (0,\infty)\times L_{p,q}$. Therefore $C_{p,q}^{\mathrm{reg}}$ is a minimal hypersurface of $\mathbb{R}^{m+2}$.
[/step]
[step:Write the second variation of the Simons cone in radial and angular variables]
Now assume $p=q=3$. Then $m=6$, $a=b=1/\sqrt{2}$, and
\begin{align*}
L:=L_{3,3}
=
S^3(1/\sqrt{2})\times S^3(1/\sqrt{2})
\subset \mathbb{S}^7.
\end{align*}
From the principal curvature computation above, the principal curvatures of $L\subset\mathbb{S}^7$ are $-1$ with multiplicity $3$ and $1$ with multiplicity $3$. Hence
\begin{align*}
|A_L|^2=3\cdot 1^2+3\cdot 1^2=6.
\end{align*}
The cone $C:=C_{3,3}^{\mathrm{reg}}$ has dimension $7$, and its squared second fundamental form is
\begin{align*}
|A_C(r,z)|^2=\frac{|A_L(z)|^2}{r^2}=\frac{6}{r^2}.
\end{align*}
Let $\phi\in C_c^\infty(C)$, and define its pullback $\widetilde{\phi}:(0,\infty)\times L\to \mathbb{R}$ under $\Psi$ by
\begin{align*}
\widetilde{\phi}(r,z):=\phi(rz) \quad \text{for } (r,z)\in (0,\infty)\times L.
\end{align*}
The [second variation formula](/theorems/2729) for a Euclidean minimal hypersurface gives
\begin{align*}
\delta^2 A_C[\phi\nu]
=
\int_C
\left(
|\nabla_C\phi|^2-|A_C|^2\phi^2
\right)
\,d\mathcal{H}^7.
\end{align*}
Here $\nabla_C\phi$ is the intrinsic gradient on $C$, and $\mathcal{H}^7$ is the $7$-dimensional [Hausdorff measure](/page/Hausdorff%20Measure) on $C$. In the warped product coordinates $(r,z)$,
\begin{align*}
|\nabla_C\phi|^2
=
|\partial_r\widetilde{\phi}|^2+\frac{1}{r^2}|\nabla_L\widetilde{\phi}|^2,
\end{align*}
and the Hausdorff measure decomposes as
\begin{align*}
d\mathcal{H}^7(rz)
=
r^6\,d\mathcal{L}^1(r)\,d\mathcal{H}^6_L(z),
\end{align*}
where $\mathcal{H}^6_L$ is the induced $6$-dimensional Hausdorff measure on $L$. Therefore
\begin{align*}
\delta^2 A_C[\phi\nu]
=
\int_0^\infty\int_L
\left(
|\partial_r\widetilde{\phi}(r,z)|^2
+
\frac{1}{r^2}|\nabla_L\widetilde{\phi}(r,z)|^2
-
\frac{6}{r^2}\widetilde{\phi}(r,z)^2
\right)
r^6
\,d\mathcal{H}^6_L(z)\,d\mathcal{L}^1(r).
\end{align*}
Since $|\nabla_L\widetilde{\phi}|^2\geq 0$, it is enough to prove
\begin{align*}
\int_0^\infty\int_L
|\partial_r\widetilde{\phi}(r,z)|^2 r^6
\,d\mathcal{H}^6_L(z)\,d\mathcal{L}^1(r)
\geq
6
\int_0^\infty\int_L
\widetilde{\phi}(r,z)^2 r^4
\,d\mathcal{H}^6_L(z)\,d\mathcal{L}^1(r).
\end{align*}
[guided]
The purpose of this step is to translate geometric stability into an explicit inequality. The cone is minimal, so the second variation of area for a compactly supported normal variation $V=\phi\nu$ is governed by the Jacobi quadratic form
\begin{align*}
\delta^2 A_C[\phi\nu]
=
\int_C
\left(
|\nabla_C\phi|^2-|A_C|^2\phi^2
\right)
\,d\mathcal{H}^7.
\end{align*}
The negative term is the possible source of instability: curvature tries to make the quadratic form negative.
For the Simons cone, the link
\begin{align*}
L=
S^3(1/\sqrt{2})\times S^3(1/\sqrt{2})
\end{align*}
has six nonzero principal curvatures, namely three equal to $1$ and three equal to $-1$. Thus
\begin{align*}
|A_L|^2=3\cdot 1^2+3\cdot (-1)^2=6.
\end{align*}
When passing from the link to the cone, each angular principal curvature scales by $1/r$, while the radial direction contributes principal curvature $0$. Hence
\begin{align*}
|A_C(r,z)|^2=\frac{|A_L(z)|^2}{r^2}=\frac{6}{r^2}.
\end{align*}
Now write the variation function in cone coordinates. Define the map $\widetilde{\phi}:(0,\infty)\times L\to \mathbb{R}$ by
\begin{align*}
\widetilde{\phi}(r,z):=\phi(rz) \quad \text{for } (r,z)\in (0,\infty)\times L.
\end{align*}
Because $\phi$ is compactly supported away from the vertex, $\widetilde{\phi}$ is smooth and compactly supported in $[\varepsilon,R]\times L$ for some numbers $0<\varepsilon<R<\infty$. This compact support justifies all integrations by parts and all weighted integrals below.
The metric on the cone is
\begin{align*}
g_C=d\mathcal{L}^1(r)^2+r^2g_L.
\end{align*}
Therefore the intrinsic gradient decomposes into radial and angular parts:
\begin{align*}
|\nabla_C\phi|^2
=
|\partial_r\widetilde{\phi}|^2+\frac{1}{r^2}|\nabla_L\widetilde{\phi}|^2.
\end{align*}
The $7$-dimensional Hausdorff measure on the cone decomposes as
\begin{align*}
d\mathcal{H}^7(rz)
=
r^6\,d\mathcal{L}^1(r)\,d\mathcal{H}^6_L(z),
\end{align*}
because the link has dimension $6$. Substituting these identities into the second variation formula gives
\begin{align*}
\delta^2 A_C[\phi\nu]
=
\int_0^\infty\int_L
\left(
|\partial_r\widetilde{\phi}(r,z)|^2
+
\frac{1}{r^2}|\nabla_L\widetilde{\phi}(r,z)|^2
-
\frac{6}{r^2}\widetilde{\phi}(r,z)^2
\right)
r^6
\,d\mathcal{H}^6_L(z)\,d\mathcal{L}^1(r).
\end{align*}
The angular term is nonnegative, so it can only help stability. Thus the entire problem reduces to proving that the radial energy controls the negative curvature term:
\begin{align*}
\int_0^\infty\int_L
|\partial_r\widetilde{\phi}(r,z)|^2 r^6
\,d\mathcal{H}^6_L(z)\,d\mathcal{L}^1(r)
\geq
6
\int_0^\infty\int_L
\widetilde{\phi}(r,z)^2 r^4
\,d\mathcal{H}^6_L(z)\,d\mathcal{L}^1(r).
\end{align*}
This is exactly where the dimension of the Simons cone enters: the radial weight is $r^6$, and the corresponding Hardy constant is large enough to dominate the curvature constant $6$.
[/guided]
[/step]
[step:Apply the weighted Hardy inequality to dominate the curvature term]
For each fixed $z\in L$, define the map $f_z:(0,\infty)\to\mathbb{R}$ by
\begin{align*}
f_z(r):=\widetilde{\phi}(r,z) \quad \text{for } r\in (0,\infty).
\end{align*}
Since $\widetilde{\phi}$ is compactly supported away from $r=0$, each $f_z$ belongs to $C_c^\infty(0,\infty)$. We prove the weighted Hardy inequality
\begin{align*}
\int_0^\infty |f_z'(r)|^2 r^6\,d\mathcal{L}^1(r)
\geq
\frac{25}{4}
\int_0^\infty f_z(r)^2 r^4\,d\mathcal{L}^1(r).
\end{align*}
Integrating by parts with no boundary term, because $f_z$ is compactly supported in $(0,\infty)$, gives
\begin{align*}
\int_0^\infty f_z(r)^2 r^4\,d\mathcal{L}^1(r)
&=
-\frac{2}{5}
\int_0^\infty f_z(r)f_z'(r) r^5\,d\mathcal{L}^1(r).
\end{align*}
Taking absolute values and applying the [Cauchy-Schwarz inequality](/theorems/432) in the [measure space](/page/Measure%20Space) $((0,\infty),\mathcal{B}(0,\infty),\mathcal{L}^1)$ to the functions $r\mapsto f_z'(r)r^3$ and $r\mapsto f_z(r)r^2$ yields
\begin{align*}
\int_0^\infty f_z(r)^2 r^4\,d\mathcal{L}^1(r)
&\leq
\frac{2}{5}
\left(
\int_0^\infty |f_z'(r)|^2 r^6\,d\mathcal{L}^1(r)
\right)^{1/2}
\left(
\int_0^\infty f_z(r)^2 r^4\,d\mathcal{L}^1(r)
\right)^{1/2}.
\end{align*}
If the rightmost weighted $L^2$ norm of $f_z$ is zero, the desired inequality is immediate. Otherwise, dividing by its square root gives
\begin{align*}
\left(
\int_0^\infty f_z(r)^2 r^4\,d\mathcal{L}^1(r)
\right)^{1/2}
&\leq
\frac{2}{5}
\left(
\int_0^\infty |f_z'(r)|^2 r^6\,d\mathcal{L}^1(r)
\right)^{1/2}.
\end{align*}
Squaring proves the claimed Hardy inequality.
Integrating this inequality over $L$ with respect to $\mathcal{H}^6_L$ gives
\begin{align*}
\int_0^\infty\int_L
|\partial_r\widetilde{\phi}(r,z)|^2 r^6
\,d\mathcal{H}^6_L(z)\,d\mathcal{L}^1(r)
\geq
\frac{25}{4}
\int_0^\infty\int_L
\widetilde{\phi}(r,z)^2 r^4
\,d\mathcal{H}^6_L(z)\,d\mathcal{L}^1(r).
\end{align*}
Since $25/4>6$, we have
\begin{align*}
\int_0^\infty\int_L
|\partial_r\widetilde{\phi}(r,z)|^2 r^6
\,d\mathcal{H}^6_L(z)\,d\mathcal{L}^1(r)
-
6
\int_0^\infty\int_L
\widetilde{\phi}(r,z)^2 r^4
\,d\mathcal{H}^6_L(z)\,d\mathcal{L}^1(r)
\geq
\frac{1}{4}
\int_0^\infty\int_L
\widetilde{\phi}(r,z)^2 r^4
\,d\mathcal{H}^6_L(z)\,d\mathcal{L}^1(r)
\geq 0.
\end{align*}
Returning to the second variation formula and keeping the nonnegative angular term, we obtain
\begin{align*}
\delta^2 A_C[\phi\nu]\geq 0.
\end{align*}
Thus $C_{3,3}$ is stable for every compactly supported smooth normal variation away from the origin.
[/step]