Although I believe that the optimal constant is $1$. I was only able to prove the inequality up to a multiplicative factor of $\color{blue}{\frac{3}{2}}$,
$$ \sum_{k=1}^{n} \frac{k^2}{a_1^2 + a_2^2 + \cdots + a_k^2} \leq \color{blue}{\frac{3}{2}} \Biggl( \sum_{k=1}^{n} \frac{1}{a_k} \Biggr)^2, $$
a slight improvement from other answers.
Step 0 - Integral Analogue. Its integral analog is much easier to prove. For positive increasing function $f$,
\begin{align*}
\left( \int_{0}^{\infty} \frac{1}{f(x)} \, \mathrm{d}x \right)^2
&= \int_{0}^{\infty}\frac{2}{f(r)} \biggl( \int_{0}^{r} \frac{1}{f(t)} \, \mathrm{d}t \biggr) \mathrm{d}y \\
&\geq \int_{0}^{\infty} \frac{2r}{f(r)^2} \, \mathrm{d}r \\
&= \int_{0}^{\infty} \frac{6r^4}{f(r)^2} \left( \int_{r}^{\infty} \frac{\mathrm{d}x}{x^4} \right) \, \mathrm{d}r \\
&= \int_{0}^{\infty} \frac{6}{x^4} \left( \int_{0}^{x} \frac{r^4}{f(r)^2} \, \mathrm{d}r \right) \mathrm{d}x \\
&\geq \int_{0}^{\infty} \frac{6}{x^4} \frac{\left( \int_{0}^{x} r^2 \, \mathrm{d}r \right)^2}{\int_{0}^{x} f(r)^2 \, \mathrm{d}r} \mathrm{d}x \tag{$\because$ C–S} \\
&= \frac{2}{3} \int_{0}^{\infty} \frac{x^2}{\int_{0}^{x} f(r)^2 \, \mathrm{d}r} \mathrm{d}x.
\end{align*}
Step 1. First, we extend the sequence by letting $a_{n+k} = +\infty$ for each $k \geq 1$. This allows us to convert the finite sums to infinite series, thereby lessening some notational load. Next, we rearrange $(a_k)_{k=1}^{\infty}$ in increasing order to obtain an increasing sequence $(b_k)_{k\in\mathbb{N}}$. Then we rearrange the sum as
\begin{align*}
\left( \sum_{k=1}^{n} \frac{1}{a_k} \right)^2
&= \left( \sum_{k=1}^{\infty} \frac{1}{b_k} \right)^2
= \sum_{k,l=1}^{\infty} \frac{1}{b_k b_l}
= \sum_{r=1}^{\infty} \sum_{\substack{j, k \geq 1 \\ \max\{j,k\}=r}} \frac{1}{b_k b_l} .
\end{align*}
Step 2. Now for each $r \geq 1$, there are precisely $2r-1$ pairs $(j, k)$ for which $\max\{j,k\} = r$ holds, and for each of such pairs, we have $b_k b_l \leq b_r^2$. So
\begin{align*}
\left( \sum_{k=1}^{n} \frac{1}{a_k} \right)^2
\geq \sum_{r=1}^{\infty} \frac{2r-1}{b_r^2}
&= \sum_{r=1}^{\infty} \frac{(2r-1)r^3}{b_r^2} \sum_{j=r}^{\infty} \left( \frac{1}{j^3} - \frac{1}{(j+1)^3} \right) \tag{1} \\
&= \sum_{j=1}^{\infty} \left( \frac{1}{j^3} - \frac{1}{(j+1)^3} \right) \sum_{r=1}^{j} \frac{(2r-1)r^3}{b_r^2}.
\end{align*}
Step 3. We now get into the dirty part of analyzing the sum in the last step. By the Cauchy–Schwarz inequality,
$$ \Biggl(\sum_{r=1}^{j} \frac{(2r-1)r^3}{b_r^2} \Biggr)\Biggl( \sum_{r=1}^{j} b_r^2 \Biggr) \geq \Biggl( \sum_{r=1}^{j} \sqrt{(2r-1)r^3} \Biggr)^2. \tag{2} $$
Plugging this back and noting that $ \sum_{r=1}^{j} a_r^2 \geq \sum_{r=1}^{j} b_r^2 $, we are led to
$$ \left( \sum_{k=1}^{n} \frac{1}{a_k} \right)^2
\geq \sum_{j=1}^{\infty} \frac{j^2} {\sum_{r=1}^{j} a_r^2} \cdot \underbrace{ \frac{1}{j^2} \left( \frac{1}{j^3} - \frac{1}{(j+1)^3} \right) \Biggl( \sum_{r=1}^{j} \sqrt{(2r-1)r^3} \Biggr)^2 }_{=(*)} $$
Step 4. So it suffices to prove that $(*)$ is bounded below by $\frac{2}{3}$. By noting that $f(x) = \sqrt{2 - x}$ is concave on $[0, 1]$, we get $f(x) \geq x f(1) + (1-x) f(0) $. Then plugging $x = 1/r$ gives
$$ \sqrt{(2r-1)r^3} \geq \sqrt{2} r^2 - (\sqrt{2} - 1) r $$
for $r \geq 1$. Summing both sides for $r = 1, 2, \ldots, j$,
$$ \sum_{j=1}^{r} \sqrt{(2r-1)r^3}
\geq \sum_{r=1}^{j} \left( \sqrt{2} r^2 - (\sqrt{2} - 1) r \right)
= \frac{1}{6} j (j + 1) \bigl( 2\sqrt{2} j + (\sqrt{2}-1)^2 \bigr). $$
So, squaring both sides, multiplying $\frac{1}{j^2}\Bigl( \frac{1}{j^3} - \frac{1}{(j+1)^3} \Bigr)$, and substituting $x = 1/j \in [0, 1]$ gives
\begin{align*}
(*)
&\quad \geq \frac{1}{36} \biggl( \frac{3j^2 + 3j + 1}{j^3(j+1)} \biggr) \bigl( 2\sqrt{2} j + (\sqrt{2}-1)^2 \bigr)^2 \\
&\quad = \frac{1}{36} \biggl( 3 + \frac{x^2}{1+x} \biggr) \bigl( 2\sqrt{2} + (\sqrt{2}-1)^2 x \bigr)^2
\end{align*}
This is increasing in $x$, hence is further bounded below by
$$ \frac{1}{36} \cdot 3 \cdot \bigl( 2\sqrt{2} \bigr)^2 = \frac{2}{3}. $$
Therefore the claim follows.
Remarks.
Here is my thought as to why we pick up an extra multiplicative factor.
Note that $(1)$ is rather conservative, becoming precise when $b_1 = b_2 = \cdots = b_n$. On the other hand, the Cauchy–Schwarz inequality, $(2)$ will become closer to an equality if $b_r \asymp r$. So, the steps $(1)$ and $(2)$ are optimized for different types of sequences, thereby rendering the entire argument unoptimized.
In $\text{(1)}$, we utilize the equality
$$ 1 = r^p \sum_{j=r}^{\infty} \left( \frac{1}{j^p} - \frac{1}{(j+1)^p} \right) $$
with $p = 3$. Although one may play with the values of $p$, the choice $p = 3$ optimizes the resulting multiplicative constant.