I'm attempting to recreate some plots from this paper on neutrino opacity calculations for interacting matter at supra-nuclear densities. Namely, I'm trying to write a Python script to perform the integration in Eqn. $(15)$, so as to obtain Fig. 3. Eqn. $(15)$ gives the absorption mean free path as:
$$\frac{\sigma(E_1)}{V} = \frac{G_F^2}{4\pi^2}(\mathcal{V}^2 + 3\mathcal{A}^2)\int_{-\infty}^{E_1}\frac{E_3}{E_1}[1-f_3(E_3)]dq_0 \;\times\;\int_{|q_0|}^{2E_1 - q_0}qS(q_0,q)dq$$
Where $G_F, \mathcal{V}$ and $\mathcal{A}$ are constants, $f_i$ the Fermi-Dirac distribution, $E_i$ the single particle energies, $q_0 = E_1 - E_3$ the energy transfer, $q = |\vec{p_1} - \vec{p_3}|$, and $S(q_0,q)$ the structure function, which characterizes the isospin response of the system, and which depends upon the chemical potentials.
The paper notes that "the single particle energies and chemical potentials depend on the ambient matter conditions (i.e. the density and temperature), and also on the interactions among the various other particles". It isn't clear to me, unfortunately, how I might go about determining the single particle energies $E_i$ such that I can perform the integrations above. The figure I wish to recreate simply notes that the neutrino energy is taken to be the Fermi energy $E_\nu = \mu_\nu$.
I'm thinking that, since the particles are relativistic, I can re-write their energies in terms of their masses and momenta (e.g. $E_1^2 = m_1^2 - p_1^2$). The integral then becomes:
$$\frac{\sigma(E_1)}{V} = \frac{G_F^2}{4\pi^2}(\mathcal{V}^2 + 3\mathcal{A}^2)\int_{-\infty}^{\sqrt{m_1^2 + p_1^2}}\frac{\sqrt{m_3^2 + p_3^2}}{\sqrt{m_1^2 + p_1^2}}[1-f_3(\sqrt{m_3^2 + p_3^2})]dq_0 \;\times\;\int_{|q_0|}^{2(\sqrt{m_1^2 + p_1^2}) - q_0}qS(q_0,q)dq$$
$$ = \frac{G_F^2}{4\pi^2}(\mathcal{V}^2 + 3\mathcal{A}^2)\int_{-\infty}^{\sqrt{m_1^2 + p_1^2}}\frac{\sqrt{m_3^2 + p_3^2}}{\sqrt{m_1^2 + p_1^2}}[1-f_3(\sqrt{m_3^2 + p_3^2})]d(\sqrt{m_1^2 + p_1^2} - \sqrt{m_3^2 + p_3^2}) \;\times\;\int_{|\sqrt{m_1^2 + p_1^2} - \sqrt{m_3^2 + p_3^2}|}^{2(\sqrt{m_1^2 + p_1^2}) - (\sqrt{m_1^2 + p_1^2} - \sqrt{m_3^2 + p_3^2})}qS(q_0,q)dq$$
It isn't clear to me, however, that this can be massaged into a product of integrals written e.g. explicitly in terms of $q, q_0$.
How do I determine or compute the single particle energies, or re-express the integral in terms of $q,q_0$?
Edit. It occurs to me that, if I let $E_\nu = \mu_\nu$ (as in the figure I wish to recreate), I can express the original integral as:
$$\frac{\sigma(\mu_\nu)}{V} = \frac{G_F^2}{4\pi^2}(\mathcal{V}^2 + 3\mathcal{A}^2)\int_{-\infty}^{\mu_\nu}\frac{E_3}{\mu_\nu}[1-f_3(E_3)]d(\mu_\nu-E_3) \;\times\;\int_{|(\mu_\nu-E_3)|}^{2\mu_\nu - (\mu_\nu-E_3)}qS((\mu_\nu-E_3),q)dq$$
I think this is the same as:
$$\frac{\sigma(\mu_\nu)}{V} = \frac{G_F^2}{4\pi^2}(\mathcal{V}^2 + 3\mathcal{A}^2)\int_{-\infty}^{\mu_\nu}\int_{|(\mu_\nu-E_3)|}^{\mu_\nu+E_3}\frac{E_3}{\mu_\nu}[1-f_3(E_3)] qS((\mu_\nu-E_3),q)dE_3dq$$
---I'm not quite sure that this is correct, and it still isn't clear to me how one might go about determining the chemical potentials (upon which the function $S(q_0, q)$ depends) in this case.
