Let's compute this step by step.
Let $\xi^a$ be a Killing vector. We can lower the index to make it a one-form, so that we can apply the covariant derivative,
$$(d\xi)_{ab} = 2 \nabla_{[a}\xi_{b]} = (\nabla_{a}\xi_{b}-\nabla_{b}\xi_{a}) = 2 \nabla_a\xi_b$$
where we used Killing's equation, $\nabla_a \xi_b+\nabla_b \xi_a=0$. Now we can apply the Hodge dual. In general, the Hodge dual of a $p-$form $\alpha$ is an $(n-p)-$form $\star \alpha$, with components given as follows:
$$(\star\alpha)_{\mu_{1} \cdots \mu_{n-p}} = \frac{1}{p!} \alpha^{\nu_1\cdots \nu_p}\epsilon_{\nu_1 \cdots \nu_p \mu_1 \cdots \mu_{n-p}}$$
where $\epsilon$ is the Levi-Civita tensor, defined as
$$\epsilon_{\mu_1 \cdots \mu_n}=\sqrt{|\det{g}|}\, \widetilde{\epsilon}_{\mu_1 \cdots \mu_n}$$
where $\widetilde{\epsilon}$ is the Levi-Civita symbol, whose components are $\pm 1$. Then, applying this to the two-form $d\xi$ gives
$$(\star d\xi)_{ab} = \frac{1}{2}\epsilon_{abcd} (d\xi)^{cd} = \epsilon_{abcd} \nabla^c\xi^d$$
We can take an exterior derivative of this. In general, if $\alpha$ is a two-form, then
$$(d\alpha)_{abc} = 3 \nabla_{[a}\alpha_{bc]}$$
which in our case implies
$$(d\star d \xi)_{abc} = \frac{1}{4}\bigg[\nabla_a ((\star d\xi)_{bc})-\nabla_b ((\star d\xi)_{ac})+\nabla_c ((\star d\xi)_{ab})\bigg]=\\
=\frac{1}{4}\bigg[ \epsilon_{bcde} \nabla_a \nabla^d \xi^e -\epsilon_{acde} \nabla_b \nabla^d \xi^e+\epsilon_{abde} \nabla_c \nabla^d \xi^e\bigg]$$
where we used the fact that the Levi-Civita tensor is covariantly constant, see for example eq.s $(1.2.36-37)$ in these lecture notes.
Now we will compute the dual of this. According to the above definition,
$$(\star d\star d \xi)_f =\epsilon_{abcf} (d\star d \xi)^{abc}$$
Now we can use the following identity for $\epsilon$, which you can find on Wikipedia:
$$\epsilon^{\mu_1 \cdots \mu_p \alpha_1 \cdots \alpha_{n-p}}\epsilon_{\mu_1 \cdots \mu_n \beta_1 \cdots \beta_{n-p}} = s\, p! (n-p)!\, \delta^{[\alpha_1}_{\beta_1} \delta^{\alpha_2}_{\beta_2} \cdots \delta^{\alpha_{n-p}]}_{\beta_{n-p}}$$
where $s$ is the sign of the determinant of the metric, in our case $s=-1$. Specialising to the case of interest for us, with $n=4$, $p=2$, we find
$$\epsilon^{abcd}\epsilon_{abef} = - 4\, \delta^{[c}_{e} \delta^{d]}_{f} = - 2\, (\delta^{c}_{e} \delta^{d}_{f}-\delta^{d}_{e} \delta^{c}_{f})$$
Therefore, we find
\begin{align*}
(\star d\star d \xi)_f &=\frac{1}{4}\bigg[\epsilon_{abcf} \epsilon^{bcde} \nabla^a \nabla_d \xi_e -\epsilon_{abcf} \epsilon^{acde} \nabla^b \nabla_d \xi_e+\epsilon_{abcf} \epsilon^{abde} \nabla^c \nabla_d \xi_e\bigg]=\\
&=\frac{1}{4}\bigg[\epsilon_{bcfa} \epsilon^{bcde} \nabla^a \nabla_d \xi_e +\epsilon_{acbf} \epsilon^{acde} \nabla^b \nabla_d \xi_e+\epsilon_{abcf} \epsilon^{abde} \nabla^c \nabla_d \xi_e \bigg]= \\
&=-\frac{1}{2}\bigg[\nabla^a \nabla_f \xi_a - \nabla^a \nabla_a \xi_f + \nabla^b \nabla_b \xi_f-\nabla^b \nabla_f \xi_b+ \nabla^c \nabla_c \xi_f-\nabla^c \nabla_f \xi_c\bigg]=\\
&=-\frac{1}{2}\bigg[\nabla^a \nabla_a \xi_f-\nabla^a \nabla_f \xi_a\bigg]
\end{align*}
Now we note the Ricci identity for Killing vectors,
$$\nabla_a \nabla_b \xi^c = R^{c}_{\,\,bad} \xi^d$$
which implies
$$\nabla^a \nabla_b \xi_a = R^{a}_{\,\,bad} \xi^d = R_{bd} \xi^d = J_b\\
\nabla_a \nabla^a \xi_c = R^{\,\,a}_{c\,\,ad} \xi^d = - R_{cd} \xi^d = -J_c$$
This means
$$(\star d\star d \xi)_f = J_f$$
or simply $\star d\star d \xi = J$. Now apply $\star$ on both sides, and note that $\star^2 = s (-1)^{p(n-p)}$ on a $p-$form in $n-$dimensional spacetime. So in this case $\star^2 = 1$ and we conclude that $\star J = d\star d \xi$.
