What is the matter/energy flux on the cosmological horizons? Consider the usual standard FLRW cosmology.  Depending on the model, there may be two distinct horizons around any local comoving observer:  the particles (or causality) horizon and the event horizon.  In general, these spherical horizons are evolving (i.e expanding or contracting) with the cosmic time.  Matter and energy are also expanding, at a different rate.  What should be the general expressions defining the matter/energy flux through both causal and event horizons, in the standard FLRW cosmology?  What should we except, from the Friedmann-Lemaître equations?
For reference, the proper distance at cosmic time $t_0$ ("today") to the causal and event horizons are respectively the following (assuming models that start by a Big Bang at cosmic time $t = 0$):
\begin{align}
D_{\mathcal{C}}( t_0) &= a(t_0) \int_0^{t_0} \frac{1}{a(t)} \, dt, \tag{1} \\[2ex]
D_{\mathcal{E}}( t_0) &= a(t_0) \int_{t_0}^{\infty} \frac{1}{a(t)} \, dt. \tag{2}
\end{align}
Friedmann-Lemaître equations are the following (I'm using units such that $c \equiv 1$):
\begin{gather}
\frac{\dot{a}^2}{a^2} + \frac{k}{a^2} = \frac{8 \pi G}{3} \, \rho + \frac{\Lambda}{3}, \tag{3} \\[2ex]
\frac{\ddot{a}}{a} = -\, \frac{4 \pi G}{3}(\rho + 3 p) + \frac{\Lambda}{3}, \tag{4}
\end{gather}
The cosmological fluid particles (or "galaxies") global velocity at time $t_0$ is
$$\tag{5}
v_{\mathcal{P}}(t_0) \equiv \frac{d D_P}{d t} = H(t_0) \, D_{\mathcal{P}}(t_0),
$$
where $D_{\mathcal{P}}(t_0)$ is the particle's proper distance and $H(t_0) = \dot{a}/a$ is the expansion rate today (i.e. Hubble's constant).
Using the previous equations, I would like to know the matter flux through both horizons, if possible in a fluid model independent way.  In other words: How many "galaxies" (i.e mass and energy) are appearing or disappearing from sight, per unit of cosmological time?
Take note that an event horizon exists only if space is accelerated, i.e when there is dark energy in the cosmological model.  But all models that start with a Big Bang have a causality horizon.  I'm interested in both horizons, but mostly the causality horizon.  I'm not interested in speculations about the cosmological horizon's entropy and its temperature, or Hawking particles flux through the horizons.
I suspect that the total energy/matter flux through an horizon should be proportional to $\rho + p$, since the deSitter universe (with $p = -\, \rho$) has a static event horizon, so its vacuum energy flux should be 0.
 A: For the particle horizon
$${\rm r_P=c \ a(t) \ \int_0^t} \frac{{\rm d}t}{{\rm a}(t)} = {\rm c \ a \ \int_0^a} \frac{{\rm d}a}{a^2 \ {\rm H}(a)}$$
the time derivative is
$$\rm \dot{r}_P = \frac{d r_P}{d t}=c+r_P \ H $$
and the matter density thins out with
$$\rm \rho_M = \frac{\rho_{M0}}{a^3}$$
so that time derivative is
$$\rm \dot{\rho}_M = \frac{d \rho_M}{d t}=- \frac{3 \ \rho_{M0} \ H}{a^4}$$
The mass of the matter inside the volume
$$\rm V_P=\frac{4 \ \pi \ r_P^3}{3}$$
with the time derivative
$$\rm \dot{V}_P = \frac{d V_P}{d t} = 4 \ \pi \ r_P^2 \ \dot{r}_P$$
of the particle horizon is
$$\rm M_P = \rho_M \ V_P$$
so the time derivative of that is
$$\rm \dot{M}_P = \frac{d M_P}{d t} = V_P \ \dot{\rho}_M + \dot{V}_P \ \rho_M$$
which today at about $\rm t=13.8 \ Glyr$ is around
$$\rm \dot{M}_P= +2 \cdot 10^{36} \ kg/s$$
for matter and dark matter combined with
$$\rm \rho_{M0}=3  \ H_0^2 \ \Omega_{M0}/( 8 \ \pi \ G), \ \Omega_{M0}=0.315, \ H_0=67150 \ m/s/Mpc$$
so there is matter (bright and dark) worth around a million solar masses per second coming into our growing particle horizon $\rm r_P$ and therefore the sign of $\rm \dot{M}_P$ is positive .
The same method is used for the event horizon $\rm r_E$, just with other integration limits. In that case there are galaxies leaving the volume, and the sign of $\rm \dot{M}_E$ is negative.
That the matter flux is positive for the particle horizon, but negative for the event horizon can also be seen in comoving coordinates where the matter per comoving volume stays constant and $\rm R_P=r_P/a$ grows with time converging to $\rm 63 \ Glyr$, while $\rm R_E=r_E/a$ asymptotically shrinks to $0$.
Here I used the $\rm H(a)$ and $\rm a(t)$ for the matter and dark energy dominated universe we live in today, if you want to look into the radiation dominated era where $\rm \dot{r}_P=2 c$ or construct some other models or density proportions see here for the required functions.
I know that is not a perfect answer to your question since you wanted to use your own equations given in the question, but that's the best I can do for now.
