[proofplan]
The proof rewrites the numerator as an inner product for the positive definite quadratic form induced by $\Sigma$. A direct Cauchy--Schwarz argument for this inner product gives the upper bound $J(a) \leq d^\top \Sigma^{-1}d$. The equality condition in that Cauchy--Schwarz inequality identifies exactly when $a$ is a non-zero scalar multiple of $\Sigma^{-1}d$, which proves both the maximisers and the maximum value.
[/proofplan]
[step:Define the inner product generated by $\Sigma$ and rewrite the criterion]
Since $\Sigma$ is real symmetric positive definite, the map
\begin{align*}
(\cdot,\cdot)_\Sigma: \mathbb{R}^p \times \mathbb{R}^p &\to \mathbb{R} \\
(x,y) &\mapsto x^\top \Sigma y
\end{align*}
is an inner product on $\mathbb{R}^p$. Define
\begin{align*}
b := \Sigma^{-1}d \in \mathbb{R}^p.
\end{align*}
Because $d \neq 0$ and $\Sigma$ is invertible, $b \neq 0$. Moreover $\Sigma b = d$, so for every $a \in \mathbb{R}^p_0$,
\begin{align*}
a^\top d
=
a^\top \Sigma b
=
(a,b)_\Sigma,
\qquad
a^\top \Sigma a
=
(a,a)_\Sigma.
\end{align*}
Thus
\begin{align*}
J(a)
=
\frac{(a,b)_\Sigma^2}{(a,a)_\Sigma}.
\end{align*}
[/step]
[step:Prove the Cauchy--Schwarz bound in the $\Sigma$ inner product]
[claim:Cauchy--Schwarz inequality for $(\cdot,\cdot)_\Sigma$]
For all $x,y \in \mathbb{R}^p$,
\begin{align*}
(x,y)_\Sigma^2 \leq (x,x)_\Sigma (y,y)_\Sigma.
\end{align*}
If $y \neq 0$, equality holds if and only if $x$ is a scalar multiple of $y$.
[/claim]
[proof]
If $y = 0$, then both sides are $0$. Assume $y \neq 0$. For every $t \in \mathbb{R}$, positive definiteness gives
\begin{align*}
0
\leq
(x - ty, x - ty)_\Sigma
=
(x,x)_\Sigma - 2t(x,y)_\Sigma + t^2(y,y)_\Sigma.
\end{align*}
Since $(y,y)_\Sigma > 0$, substitute
\begin{align*}
t := \frac{(x,y)_\Sigma}{(y,y)_\Sigma}.
\end{align*}
Then
\begin{align*}
0
\leq
(x,x)_\Sigma
-
\frac{(x,y)_\Sigma^2}{(y,y)_\Sigma}.
\end{align*}
Multiplying by $(y,y)_\Sigma > 0$ gives
\begin{align*}
(x,y)_\Sigma^2 \leq (x,x)_\Sigma (y,y)_\Sigma.
\end{align*}
Equality holds exactly when
\begin{align*}
(x - ty, x - ty)_\Sigma = 0.
\end{align*}
By positive definiteness, this is equivalent to $x - ty = 0$, hence to $x$ being a scalar multiple of $y$.
[/proof]
Applying the claim with $x := a$ and $y := b$ gives
\begin{align*}
(a,b)_\Sigma^2
\leq
(a,a)_\Sigma (b,b)_\Sigma.
\end{align*}
Since $a \in \mathbb{R}^p_0$, positive definiteness gives $(a,a)_\Sigma > 0$, so division by $(a,a)_\Sigma$ yields
\begin{align*}
J(a)
=
\frac{(a,b)_\Sigma^2}{(a,a)_\Sigma}
\leq
(b,b)_\Sigma.
\end{align*}
Finally,
\begin{align*}
(b,b)_\Sigma
=
b^\top \Sigma b
=
(\Sigma^{-1}d)^\top \Sigma(\Sigma^{-1}d)
=
d^\top \Sigma^{-1}d,
\end{align*}
where the last equality uses $\Sigma^\top = \Sigma$.
[/step]
[step:Identify exactly when equality occurs]
The previous step shows that every $a \in \mathbb{R}^p_0$ satisfies
\begin{align*}
J(a) \leq d^\top \Sigma^{-1}d.
\end{align*}
Equality holds exactly when equality holds in the Cauchy--Schwarz inequality for $x := a$ and $y := b$. Since $b = \Sigma^{-1}d \neq 0$, this is equivalent to the existence of $\lambda \in \mathbb{R}$ such that
\begin{align*}
a = \lambda b = \lambda \Sigma^{-1}d.
\end{align*}
The condition $a \in \mathbb{R}^p_0$ forces $\lambda \neq 0$.
Conversely, if $a = \lambda \Sigma^{-1}d$ with $\lambda \in \mathbb{R}\setminus\{0\}$, then
\begin{align*}
J(a)
=
\frac{(\lambda b^\top d)^2}{\lambda^2 b^\top \Sigma b}
=
\frac{(b^\top d)^2}{b^\top \Sigma b}.
\end{align*}
Since $d = \Sigma b$, we have $b^\top d = b^\top \Sigma b$, and therefore
\begin{align*}
J(a)
=
b^\top \Sigma b
=
d^\top \Sigma^{-1}d.
\end{align*}
Thus the maximisers are precisely the non-zero scalar multiples of $\Sigma^{-1}d$, and the maximum value is $d^\top \Sigma^{-1}d$.
[/step]
