How does one derive these equations for vacuum polarization and magnetization?
$\vec P = \frac{2 \epsilon_0 \alpha}{E_\text{c}^2} \bigg ( 2 ( E^2 - c^2 B^2) \vec E + 7 c^2 (\vec E \cdot \vec B) \vec B\bigg )$
$\vec M = - \frac{2 \alpha}{\mu_0 E_\text{c}^2} \bigg ( 2 ( E^2 - c^2 B^2) \vec E + 7 c^2 (\vec E \cdot \vec B) \vec B\bigg )$
Start from the QED Lagrangian $$\mathcal L_{\mathrm{QED}} = -\frac14 F_{\mu\nu} F^{\mu\nu} + \bar{\psi}(i {\not} D -m)\psi $$ where $m$ is the electron mass and $D$ is the covariant derivative. The path integral is $$Z = \int DA\, D\bar\psi\, D\psi\, e^{i S_{\mathrm{QED}}}$$ We can integrate out the fermions, obtaining a fermion determinant $$Z = \int DA\, \det{M}\, e^{i S_{\mathrm{EM}}}= \int DA\, e^{i \left(S_{\mathrm{EM}}-i\log\det{M}\right)}$$ where $M=i {\not} D -m$ and $S_{\mathrm{EM}}$ is the action for the free electromagnetic field, $S_{\mathrm{EM}} = \int d^4x\,\left(-\frac14 F_{\mu\nu} F^{\mu\nu}\right)$. So instead of considering QED (a theory of photons and electrons) we can consider a theory of photons only with effective action $$S_{\mathrm{eff}}=S_{\mathrm{EM}}-i\log\det{M}$$ This procedure is useful for example if one wants to consider photon processes with energies lower than the mass of the electron. In this case we can compute the correction $\log\det{M}$ to leading order in $m^{-1}$.
The calculation can be done through standard methods, but is somewhat complicated. It is done in the original Euler-Heisenberg papers, or in this pdf. In the appropriate limit the end result is $$\mathcal{L}_{\mathrm{eff}} = -\frac14 F_{\mu\nu} F^{\mu\nu}+\frac{\alpha^2}{90 m^4}\left[\left(F_{\mu\nu} F^{\mu\nu}\right)^2 +\frac{7}{4}\left(F_{\mu\nu} \tilde{F}^{\mu\nu}\right)^2\right] +\mathcal{O}(m^{-6})$$ Noting that $$-F_{\mu\nu} F^{\mu\nu} = 2 \left(\mathbf{E}^2 - \mathbf{B}^2\right)\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,F_{\mu\nu} \tilde{F}^{\mu\nu}=4 \mathbf{E} \cdot \mathbf{B}$$ So in terms of the electric and magnetic fields, the effective Lagrangian is $$\mathcal{L}_{\mathrm{eff}} = \frac12 \left(\mathbf{E}^2 - \mathbf{B}^2\right)+\frac{2\alpha^2}{45 m^4}\left[\left(\mathbf{E}^2 - \mathbf{B}^2\right)^2 +7\left(\mathbf{E} \cdot \mathbf{B}\right)^2\right]$$ From here, you can derive the Euler-Lagrange equations associated to this effective Lagrangian, which are $$\nabla \cdot \mathbf{D} = 0\,\,\,\,\,\,\,\,\,\,\,\,\nabla \times \mathbf{H} -\frac{\partial \mathbf{D}}{\partial t}=0\\ \nabla \cdot \mathbf{B} = 0\,\,\,\,\,\,\,\,\,\,\,\, \nabla \times \mathbf{E} +\frac{\partial \mathbf{B}}{\partial t}=0$$ where $$\mathbf{D} = \frac{\partial \mathcal{L}_{\mathrm{eff}}}{\partial \mathbf{E}} = \mathbf{E} + \mathbf{P}\,\,\,\,\,\,\,\,\,\,\,\,\mathbf{H} = -\frac{\partial \mathcal{L}_{\mathrm{eff}}}{\partial \mathbf{B}}=\mathbf{B} - \mathbf{M}$$ apart from possibly some constants that I haven't kept track of, where $\mathbf{P}$ and $\mathbf{M}$ are as given by the OP. These are Maxwell's equations in a medium with polarization $\mathbf{P}$ and magnetisation $\mathbf{M}$. The medium is interpreted to be the vacuum itself.
