The electromagnetic field momentum density in a region changes for two reasons.
One is a flux of momentum density (for instance if an electromagnetic wave enters the region), and the other is an exchange of momentum with electric charges.
Both are important. The link you cite actually covers both. If you take the time rate of change of the electromagnetic momentum density plus the force density you get the flux.
As an analogy you can look at conservation of charge. The divergence of $\vec J$ is proportional to the time rate of change of the charge density. You even have the continuity equation $\frac{\partial \rho}{\partial t}=-\vec \nabla \cdot \vec J.$ So the current density is the flux of the electric charge.
I'll say more later on to keep this answer self contained, but right now note if you look at equation 1068 of the link you cited, then yes the $\rho \vec E$ term appears in the derivation at the same time as $ \epsilon_0 \vec E \vec {\nabla} \cdot \vec E$ appears. But you should have started with the total change of the total momentum to see clearly what is going on.
Let's do that now to be self contained.
Start with $\rho \vec E + \vec J \times \vec B+\frac{\partial}{\partial t}\left(\epsilon_0 \vec E \times \vec B\right)$ as the time rate at which total momentum (field and mechanical) changes in a region. We will want to show that it can be represented as a total divergence because then we can see the flux of momentum which tells us how momentum flows. And yes at some point you have to turn $\rho \vec E$ into $ \epsilon_0 \vec E \vec {\nabla} \cdot \vec E$ and yes that represents mechanical momentum exchanging with electromagnetic momentum, but when you start with $\rho \vec E + \vec J \times \vec B+\frac{\partial}{\partial t}\left(\epsilon_0 \vec E \times \vec B\right)$ what you are computing is the flux, which is really an expression of how the momentum moves and doesn't care if momentum at that point is switching from one type to another.
We can look at both things, first how the momentum is moving, and second how it is changing right there (though the latter you might understand just fine from the Lorentz force). The only reason the total momentum can change is if it flows in. And so the field momentum in a region can change for two reasons, from the flow of momentum (flux) and from exchanging with mechanical momentum right there. So just like you can write $\frac{ \partial \rho}{\partial t}$ as a divergence of the current; similarly, you want to write $\rho \vec E + \vec J \times \vec B+\frac{\partial}{\partial t}\left(\epsilon_0 \vec E \times \vec B\right)$ as a total divergence because then that thing you are taking the divergence of can be the flux of momentum.
So now that we know what we are computing $\rho \vec E + \vec J \times \vec B+\frac{\partial}{\partial t}\left(\epsilon_0 \vec E \times \vec B\right)$ (because it is the time rate of change of total momentum density, both mechanical and field momentum) and we know that we are aiming to write it as a total divergence so that the thing being diverged is the flux of momentum. So given that, let's look at the parts. Actually I'll just repeat the derivation so we don't need your link.
So field momentum density is $
\epsilon_0 \vec E \times \vec B$ so the rate of change of field momentum is $\frac{\partial}{\partial t}\left(\epsilon_0 \vec E \times \vec B\right)$. The rate of change of mechanical momentum density is $ \rho \vec E + \vec J \times \vec B$, so the total change of total momentum density is $\rho \vec E + \vec J \times \vec B+\frac{\partial}{\partial t}\left(\epsilon_0 \vec E \times \vec B\right)$. Now we will do math and will need the identity $\vec \nabla \left(E^2/2\right)=\vec E \times \left(\vec \nabla \times \vec E\right) + \left(\vec E \cdot \vec \nabla\right)\vec E$ as well as product rules and Maxwell:
$$\rho \vec E + \vec J \times \vec B+\frac{\partial}{\partial t}\left(\epsilon_0 \vec E \times \vec B\right)=$$
$$\rho \vec E + \vec J \times \vec B+\epsilon_0 \frac{\partial \vec E}{\partial t} \times \vec B+\epsilon_0 \vec E \times \frac{\partial \vec B}{\partial t}=$$
$$\rho \vec E + \vec J \times \vec B+\left(\frac{1}{\mu_0} \vec \nabla \times \vec B-\vec J\right)\times \vec B+\epsilon_0 \vec E \times \frac{\partial \vec B}{\partial t}.$$
Now we recall that our goal is to write this as a total divergence. But right now we see two terms $ \vec J \times \vec B$ and $-\vec J \times \vec B$ and mathematically it isn't obvious that cancelling them will help us get a total divergence but it doesn't hurt and then we have less to write. And I want to point out that you didn't seem to object to this part. However there is more than simple mathematics here. The term $\vec J \times \vec B$ tells us how the mechanical momentum changes and we know that the change in mechanical momentum just goes into the fields. But this tells us that the mechanical momentum exchanged because of the magnetic force is only one of the two things that contributes to the change in $\vec E$. Maxwell should really be written as evolution equations for $\vec E$ and $\vec B$ and really its just that two things contribute to the change in $ \vec E$, circulating $\vec B$ and current. So both (circulating $\vec B$ and current) contribute to the change in field momentum there. It's just that one (the current) contributes to the field momentum by giving up mechanical momentum whereas the other changes the momentum here by having it flow in from nearby. And only that flow of momentum contributes to the flux of total momentum. So cancelling them is physically correct, not just mathematically allowed and typographically convenient. The same exact thing will happen with $\rho \vec E$ later because we only care about the flux of total momentum because that is what we are computing. So we have:
$$\rho \vec E + \left(\frac{1}{\mu_0} \vec \nabla \times \vec B\right)\times \vec B+\epsilon_0 \vec E \times \frac{\partial \vec B}{\partial t}=$$
$$\rho \vec E + \frac{-1}{\mu_0} \vec B \times \left(\vec \nabla \times \vec B\right)+\epsilon_0 \vec E \times \frac{\partial \vec B}{\partial t}=$$
$$\rho \vec E + \frac{-1}{\mu_0} \left(\vec \nabla \left(B^2/2\right)-\left(\vec B \cdot \vec \nabla\right)\vec B\right)+\epsilon_0 \vec E \times \frac{\partial \vec B}{\partial t}=$$
$$\rho \vec E + \frac{-1}{\mu_0} \left(\vec \nabla \left(B^2/2\right)-\left( \vec \nabla \cdot \vec B\right)\vec B-\left(\vec B \cdot \vec \nabla\right)\vec B\right)+\epsilon_0 \vec E \times \frac{\partial \vec B}{\partial t}=$$
$$\rho \vec E + \frac{-1}{\mu_0} \left(\vec \nabla \left(B^2/2\right)-\left( \vec \nabla \cdot \vec B\right)\vec B-\left(\vec B \cdot \vec \nabla\right)\vec B\right)+\epsilon_0 \vec E \times \left(-\vec\nabla \times \vec E\right).$$
Now we notice that every term has either magnetic fields or electric fields and the terms with magnetic fields are actually a total divergence. Let's see if we can make the electric part look like a total divergence as well.
$$\rho \vec E + \frac{-1}{\mu_0} \left(\vec \nabla \left(B^2/2\right)-\left( \vec \nabla \cdot \vec B\right)\vec B-\left(\vec B \cdot \vec \nabla\right)\vec B\right)+\epsilon_0 \vec E \times \left(-\vec\nabla \times \vec E\right)=$$
$$\frac{-1}{\mu_0} \left(\vec \nabla \left(B^2/2\right)-\left( \vec \nabla \cdot \vec B\right)\vec B-\left(\vec B \cdot \vec \nabla\right)\vec B\right)+\rho \vec E -\epsilon_0 \vec E \times \left(\vec\nabla \times \vec E\right)=$$
$$ \frac{-1}{\mu_0} \left(\vec \nabla \left(B^2/2\right)-\left( \vec \nabla \cdot \vec B\right)\vec B-\left(\vec B \cdot \vec \nabla\right)\vec B\right)+\rho \vec E -\epsilon_0\left(\vec \nabla \left(E^2/2\right)-\left(\vec E \cdot \vec \nabla\right)\vec E\right)=$$
$$ \frac{-1}{\mu_0} \left(\vec \nabla \left(B^2/2\right)-\left( \vec \nabla \cdot \vec B\right)\vec B-\left(\vec B \cdot \vec \nabla\right)\vec B\right)+\epsilon_0\left(\vec \nabla \cdot\vec E\right) \vec E -\epsilon_0\left(\vec \nabla \left(E^2/2\right)-\left(\vec E \cdot \vec \nabla\right)\vec E\right)=$$
$$ \frac{-1}{\mu_0} \left(\vec \nabla \left(B^2/2\right)-\left( \vec \nabla \cdot \vec B\right)\vec B-\left(\vec B \cdot \vec \nabla\right)\vec B\right) -\epsilon_0\left(\vec \nabla \left(E^2/2\right)-\left(\vec \nabla \cdot\vec E\right) \vec E-\left(\vec E \cdot \vec \nabla\right)\vec E\right).$$
Now both collections of terms are a total divergence. We could have just taken the curl and checked that we got $\vec 0$ if that's what we wanted to know, but the point is that we know how to get a specific thing whose divergence is what we want.
Specifically $\sum_j\frac{\partial}{\partial x_j}\left(\epsilon_0E^2\delta_{ij}/2-\epsilon_0E_iE_j+B^2\delta_{ij}/2\mu_0-B_iB_j/\mu_0\right)$ is a divergence that gives the negative of the $i$th component of $$ \frac{-1}{\mu_0} \left(\vec \nabla \left(B^2/2\right)-\left( \vec \nabla \cdot \vec B\right)\vec B-\left(\vec B \cdot \vec \nabla\right)\vec B\right) -\epsilon_0\left(\vec \nabla \left(E^2/2\right)-\left(\vec \nabla \cdot\vec E\right) \vec E-\left(\vec E \cdot \vec \nabla\right)\vec E\right).$$
So conceptually we wanted an equation like $\frac{\partial \rho }{\partial t}=-\vec \nabla \cdot \vec J$ where the flux of the charge tells us the time rate of change of the charge. We have it:
$$\rho E_i + \left[\vec J \times \vec B\right]_i + \frac{\partial}{\partial t}\left[\epsilon_0\vec E \times \vec B\right]_i=-\sum_j\frac{\partial}{\partial x_j}\left(\epsilon_0E^2\delta_{ij}/2-\epsilon_0E_iE_j+B^2\delta_{ij}/2\mu_0-B_iB_j/\mu_0\right)$$
Two conceptual issues remain. One, you mentioned, was about the exchange of momentum between the field and the charge. If there weren't any charge density the term would have already been a total divergence, just like the magnetic terms were. When there is charge, then the electric field has a divergence. This just means that field momentum can actually flow in or out from where the field lines terminate. This is not a hack. If you look at total momentum, there isn't a problem, it just flows, so there is a momentum flux. And the time rate of change of the total momentum is just the negative of the divergence of the momentum flux. Nothing mysterious in the slightest. Sure, if you look at just the field momentum now there are two ways for it to change, the momentum flux can flow some momentum in from nearby and it can exchange with the mechanical momentum.
Now, for the magnetic force it exchanged momentum with some of the changing field momentum from the changing electric field. That's actually one of the ways field momentum changes, from magnetic forces acting on moving charges and the moving charges changing the electric field above and beyond how the circulating magnetic field changes the electric field. And the other contributor to changing electric fields, circulating $\vec B$ fields, happens even when there are no moving charges, so it is about the flow of momentum through the space between charges. However for electric forces, things might at first appear different. The $\vec{B}$ field changes solely because of circulating electric fields. However, the same thing happens anyway.
Neither the mechanical momentum by itself nor the the field momentum neither by itself has its time rate of change be a total divergence. So every intuition about conserved quantities fails. Field momentum isn't conserved, it isn't a thing that just flows. Same with mechanical momentum. But total momentum is conserved, it does just flow from here to there. What you can imagine is writing each as a total divergence plus some other term, where the non divergence terms are equal an opposite. This is not unique, and it isn't physically meaningful since field momentum and mechanical momentum by themselves (in regions where there are both) aren't conserved, so we are free to mathematically pick whichever we want, and unfortunately for conceptual understanding the mathematically easiest is to have the $\rho\vec{E}$ be 100% the other (so the divergence part is zero). Meanwhile the field momentum part is $$-\epsilon_0\left(\vec \nabla \left(E^2/2\right)-\left(\vec E \cdot \vec \nabla\right)\vec E\right)=$$
$$-\epsilon_0\left(\vec \nabla \left(E^2/2\right)-\left(\vec \nabla \cdot\vec E\right) \vec E-\left(\vec E \cdot \vec \nabla\right)\vec E+\left(\vec \nabla \cdot\vec E\right) \vec E\right)=$$
$$-\epsilon_0\left(\vec \nabla \left(E^2/2\right)-\left(\vec \nabla \cdot\vec E\right) \vec E-\left(\vec E \cdot \vec \nabla\right)\vec E\right)-\epsilon_0\left(\vec \nabla \cdot\vec E\right) \vec E=$$
$$-\epsilon_0\left(\vec \nabla \left(E^2/2\right)-\left(\vec \nabla \cdot\vec E\right) \vec E-\left(\vec E \cdot \vec \nabla\right)\vec E\right)-\rho \vec E.$$
So a total divergence minus a term like $\rho\vec E$. But this is not a hack, it is just saying when the divergence of $\vec E$ is nonzero, the field momentum is not conserved and the rate at which it is lost (rather than just flowing around, since that is tracked by the total divergence term) is equal and opposite to what the mechanical momentum gains. This is exactly what we want, and therefore the opposite of a hack. So we can see how momentum flows (because the flux is the thing we take a total divergence of) and yeah it interconverts between different kinds of momentum because that's what momentum does, it is no more mysterious than when one particle exchanges momentum with a different particle. We could call it particle one momentum and particle two momentum and notice that neither by itself is conserved.
Just to be totally clear about the magnetic field, the magnetic field changes solely because of circulating electric fields. It does not change because of charges, or because of currents. However, only in the absence of charges does the electric field part flow in a way where you have a total divergence for the field momentum by itself. And this is just like how the only in the absence of current does the magnetic field part flow in a way where you have a total divergence. That's just because field momentum isn't conserved in the other cases. So its flux is not a total divergence it is a divergence plus a term just like the Lorentz Force Density. Some of it is being exchanged between the two types of momentum, and some is flowing around from here to there. And total momentum is conserved.