Electromagnetic irradiation of a dielectric: Transforming the striction force equation

Consider the simple case of electromagnetic irradiation of a homogeneous isotropic dielectric, neglecting the dispersion of the refractive index. Assuming a transparent medium, the spatial density of forces acting on the dielectric in a static external electromagnetic field can be given as

$$mathbf{f} = – nabla p – nabla epsilon dfrac{langle mathbf{E}^2 rangle}{8 pi} – nabla mu dfrac{langle mathbf{H}^2 rangle}{8 pi} + nabla left[ left( rho dfrac{partial{epsilon}}{partial{rho}} right)_T dfrac{langle mathbf{E}^2 rangle}{8 pi} + left( rho dfrac{partial{mu}}{partial{rho}} right)_T dfrac{langle mathbf{H}^2 rangle}{8 pi} right] + dfrac{epsilon mu – 1}{4 pi c} dfrac{partial}{partial{t}}langle [ mathbf{E} times mathbf{H}] rangle.$$

$p$ is the pressure in the medium (for a given density $rho$ and temperature $T$ in zero field.
$epsilon$ and $mu$ are the permittivity and magnetic permeability.
$c$ is the speed of light.
The angular brackets denote averaging over a time period far greater than the characteristic alternation period of light.

It is said that, by expressing $langle E^2 rangle$ through $I$ (the light intensity) and introducing the refractive index $n = sqrt{epsilon}$, we can transform the striction force equation to

$$mathbf{f}_{text{str}} = nabla left[ left( rho dfrac{partial{epsilon}}{partial{rho}} right)_T dfrac{langle mathbf{E}^2 rangle}{8 pi} right] = nabla left[ left( rho dfrac{partial{n}}{partial{rho}} right)_T dfrac{I}{c} right].$$

I’m trying to understand how exactly we get $nabla left[ left( rho dfrac{partial{epsilon}}{partial{rho}} right)_T dfrac{langle mathbf{E}^2 rangle}{8 pi} right] = nabla left[ left( rho dfrac{partial{n}}{partial{rho}} right)_T dfrac{I}{c} right]$. I’ve been doing a lot of research to try and understand this, but I’m stuck.

My best attempt is as follows. As said here, in optics, the time-averaged value of the radiated flux is technically known as the irradiance, more often simply referred to as the intensity. The Wikipedia article for intensity says that, if $I$ is the local intensity (I’m not completely sure if this is the correct assumption for our case), then we have that $I = dfrac{cn epsilon_0}{2}|E|^2$, where $epsilon_0$ is the vacuum permittivity. And so, if we assume that $langle mathbf{E}^2 rangle = |E|^2$ (which seems to be true, given the answer here), then we get that $|E|^2 = dfrac{2I}{cn epsilon_0}$, and so $nabla left[ left( rho dfrac{partial{epsilon}}{partial{rho}} right)_T dfrac{langle mathbf{E}^2 rangle}{8 pi} right] = nabla left[ left( rho dfrac{partial{n^2}}{partial{rho}} right)_T dfrac{I}{4 pi c n epsilon_0} right]$. But it is not clear how one proceeds from here.

Some other potentially relevant facts that I found during my research are as follows:

• According to the article on irradiance (different from the article on intensity), $E_{{mathrm {e}}}={frac {n}{2mu _{0}{mathrm {c}}}}E_{{mathrm {m}}}^{2}cos alpha ={frac {nvarepsilon _{0}{mathrm {c}}}{2}}E_{{mathrm {m}}}^{2}cos alpha$. If we let that $cos(alpha) = 1$ for our case, then this might be relevant.
• The article on vacuum permittivity states that $varepsilon _{0}={frac {1}{mu _{0}c^{2}}}$, where $mu_0$ is the vacuum permeability.
• This page on "energy density, flux and power" has numerous relevant-looking facts that include $E$ and time-averaged values, and look like they could potentially cancel the necessary factors, such as $4pi$ or $8pi$, somehow.

I would greatly appreciate it if people would please take the time to explain exactly how we get from $nabla left[ left( rho dfrac{partial{epsilon}}{partial{rho}} right)_T dfrac{langle mathbf{E}^2 rangle}{8 pi} right]$ to $nabla left[ left( rho dfrac{partial{n}}{partial{rho}} right)_T dfrac{I}{c} right]$.