for what it's worth, Mathematica gives a closed-form expression in terms of a hypergeometric function: $$\int_{0}^{\infty} k^{\ell+1} e^{-p^2k^2}(K_{\mu}(k))^2\,dk=$$ $$\mu^{-1}2^{-2 \mu-3} p^{-\ell-2 \mu-2} \left[\mu \Gamma (-\mu)^2 \Gamma \left(\frac{\ell}{2}+\mu+1\right) \, _2F_2\left(\mu+\frac{1}{2},\frac{\ell}{2}+\mu+1;\mu+1,2 \mu+1;\frac{1}{p^2}\right)+16^{\mu} \mu p^{4 \mu} \Gamma (\mu)^2 \Gamma \left(\frac{\ell}{2}-\mu+1\right) \, _2F_2\left(\frac{1}{2}-\mu,\frac{\ell}{2}-\mu+1;1-2 \mu,1-\mu;\frac{1}{p^2}\right)+\pi l \left(-4^{\mu}\right) \Gamma \left(\frac{\ell}{2}\right) p^{2 \mu}\frac{1}{ \sin (\pi \mu)} \, _2F_2\left(\frac{1}{2},\frac{\ell}{2}+1;1-\mu,\mu+1;\frac{1}{p^2}\right)\right],$$$$\mu^{-1}2^{-2 \mu-3} p^{-\ell-2 \mu-2} \left[\mu \Gamma (-\mu)^2 \Gamma \left(\frac{\ell}{2}+\mu+1\right) \, _2F_2\left(\mu+\frac{1}{2},\frac{\ell}{2}+\mu+1;\mu+1,2 \mu+1;\frac{1}{p^2}\right)+16^{\mu} \mu p^{4 \mu} \Gamma (\mu)^2 \Gamma \left(\frac{\ell}{2}-\mu+1\right) \, _2F_2\left(\frac{1}{2}-\mu,\frac{\ell}{2}-\mu+1;1-2 \mu,1-\mu;\frac{1}{p^2}\right)-4^\mu \pi l \Gamma \left(\frac{\ell}{2}\right) p^{2 \mu}\frac{1}{ \sin \pi \mu} \, _2F_2\left(\frac{1}{2},\frac{\ell}{2}+1;1-\mu,\mu+1;\frac{1}{p^2}\right)\right],$$ for $p>0$, $\ell>2|\mu|-2$.
for what it's worth, Mathematica gives a closed-form expression in terms of a hypergeometric function: $$\int_{0}^{\infty} k^{\ell+1} e^{-p^2k^2}(K_{\mu}(k))^2\,dk=$$ $$\mu^{-1}2^{-2 \mu-3} p^{-\ell-2 \mu-2} \left[\mu \Gamma (-\mu)^2 \Gamma \left(\frac{\ell}{2}+\mu+1\right) \, _2F_2\left(\mu+\frac{1}{2},\frac{\ell}{2}+\mu+1;\mu+1,2 \mu+1;\frac{1}{p^2}\right)+16^{\mu} \mu p^{4 \mu} \Gamma (\mu)^2 \Gamma \left(\frac{\ell}{2}-\mu+1\right) \, _2F_2\left(\frac{1}{2}-\mu,\frac{\ell}{2}-\mu+1;1-2 \mu,1-\mu;\frac{1}{p^2}\right)+\pi l \left(-4^{\mu}\right) \Gamma \left(\frac{\ell}{2}\right) p^{2 \mu}\frac{1}{ \sin (\pi \mu)} \, _2F_2\left(\frac{1}{2},\frac{\ell}{2}+1;1-\mu,\mu+1;\frac{1}{p^2}\right)\right],$$ for $p>0$, $\ell>2|\mu|-2$.
for what it's worth, Mathematica gives a closed-form expression in terms of a hypergeometric function: $$\int_{0}^{\infty} k^{\ell+1} e^{-p^2k^2}(K_{\mu}(k))^2\,dk=$$ $$\mu^{-1}2^{-2 \mu-3} p^{-\ell-2 \mu-2} \left[\mu \Gamma (-\mu)^2 \Gamma \left(\frac{\ell}{2}+\mu+1\right) \, _2F_2\left(\mu+\frac{1}{2},\frac{\ell}{2}+\mu+1;\mu+1,2 \mu+1;\frac{1}{p^2}\right)+16^{\mu} \mu p^{4 \mu} \Gamma (\mu)^2 \Gamma \left(\frac{\ell}{2}-\mu+1\right) \, _2F_2\left(\frac{1}{2}-\mu,\frac{\ell}{2}-\mu+1;1-2 \mu,1-\mu;\frac{1}{p^2}\right)-4^\mu \pi l \Gamma \left(\frac{\ell}{2}\right) p^{2 \mu}\frac{1}{ \sin \pi \mu} \, _2F_2\left(\frac{1}{2},\frac{\ell}{2}+1;1-\mu,\mu+1;\frac{1}{p^2}\right)\right],$$ for $p>0$, $\ell>2|\mu|-2$.
for what it's worth, Mathematica gives a closed-form expression in terms of a hypergeometric function: $$\int_{0}^{\infty} k^{\ell+1} e^{-p^2k^2}(K_{\mu}(k))^2\,dk=$$ $$\mu^{-1}2^{-2 \mu-3} p^{-\ell-2 \mu-2} \left[\mu \Gamma (-\mu)^2 \Gamma \left(\frac{\ell}{2}+\mu+1\right) \, _2F_2\left(\mu+\frac{1}{2},\frac{\ell}{2}+\mu+1;\mu+1,2 \mu+1;\frac{1}{p^2}\right)+16^{\mu} \mu p^{4 \mu} \Gamma (\mu)^2 \Gamma \left(\frac{\ell}{2}-\mu+1\right) \, _2F_2\left(\frac{1}{2}-\mu,\frac{\ell}{2}-\mu+1;1-2 \mu,1-\mu;\frac{1}{p^2}\right)+\pi l \left(-4^{\mu}\right) \Gamma \left(\frac{\ell}{2}\right) p^{2 \mu}\frac{1}{ \sin (\pi \mu)} \, _2F_2\left(\frac{1}{2},\frac{\ell}{2}+1;1-\mu,\mu+1;\frac{1}{p^2}\right)\right],$$ for $p>0$, $\ell>2|\mu|-2$.