Solution (source code)

= Solution

Well before <matter-radiation equality>, <radiation domination> makes the integrand approximately constant:
$$
r_p(a)\simeq\frac{ca}{H_0\sqrt{\Omega_{r,0}}}.
$$
At equality, $\Omega_{r,0}a_{\rm eq}^{-4}=\Omega_{m,0}a_{\rm eq}^{-3}$, so $a_{\rm eq}\simeq1/3000$ and $\Omega_{r,0}\simeq10^{-4}$. Extrapolating the radiation-only result gives $r_p(a_{\rm eq})\simeq133\,\mathrm{Mpc}$. This overestimates the horizon because matter is no longer negligible at equality and increases the expansion rate.

Retain both matter and radiation, still neglecting the insignificant early <cosmological constant>. The <particle horizon through radiation-matter equality> is
$$
r_p(a)=\frac{2c}{H_0\Omega_{m,0}}\left(\sqrt{\Omega_{r,0}+\Omega_{m,0}a}-\sqrt{\Omega_{r,0}}\right).
$$
Therefore
$$
\boxed{r_p(a_{\rm eq})\simeq2(\sqrt2-1)(133\,\mathrm{Mpc})\simeq110\,\mathrm{Mpc}.}
$$
This is a comoving distance; the proper horizon at equality is approximately $110/3000\,\mathrm{Mpc}\simeq37$ kiloparsecs.

For today's horizon, extrapolating <matter domination> gives
$$
r_p(1)\simeq\frac c{H_0}\int_0^1\frac{da}{\sqrt{\Omega_{m,0}a}}=\frac{2c}{H_0\sqrt{\Omega_{m,0}}}\simeq1.46\times10^4\,\mathrm{Mpc}.
$$
This is an order-of-magnitude estimate: <dark energy> is significant late, and radiation matters early. Numerical integration of the complete expression with the supplied density parameters and $\Omega_{r,0}\simeq10^{-4}$ gives
$$
\boxed{r_p(t_0)\simeq1.30\times10^4\,\mathrm{Mpc}.}
$$
Both neglected positive-density terms increase $H$ at fixed $a$ and reduce the horizon relative to the matter-only extrapolation.