Macroscopic Magnetohydrodynamics

$\epsilon_0 =$ permittivity of free space
${(\chi_e)_i}^j =$ electric susceptability tensor
${(\epsilon_r)_i}^j = \delta_i^j + {(\chi_e)_i}^j =$ relative permittivity of material
${\epsilon_i}^j = \epsilon_0 {(\epsilon_r)_i}^j =$ permittivity of material
$E_i = $ electric field
$P_i = \epsilon_0 {(\chi_e)_i}^j E_j$ = polarization field (is it always linear though?)
$D_i = \epsilon_0 E_i + P_i = \epsilon_0 {(\epsilon_r)_i}^j E_j = {\epsilon_i}^j E_j = $ displacement field

$\mu_0 =$ permeability of free space
${(\chi_m)_i}^j =$ magnetic susceptability
${(\mu_r)_i}^j = \delta_i^j + {(\chi_m)_i}^j = $ relative permeability of material
${\mu_i}^j = \mu_0 {(\mu_r)_i}^j =$ permeability of material
$H_i =$ magnetizing field
$M_i = \mu_0 ({\chi_m)_i}^j H_j =$ magnetization field
$B_i = \mu_0 H_i + M_i = \mu_0 ({\mu_r)_i}^j H_j = {\mu_i}^j H_j =$ magnetic field

$\rho = \rho_{bound} + \rho_{free}$
$\vec{J} = \vec{J}_{bound} + \vec{J}_{free}$

Therefore
$\vec{\nabla} \cdot \vec{E} = \frac{\rho}{\epsilon_0}$ = Microscopic Maxwell Charge:
$\epsilon_0 \vec{\nabla} \cdot \vec{E} = \rho = \rho_{free} + \rho_{bound}$ = Microscopic Maxwell Charge:
Now start with the polarization field divegence:
$\vec{\nabla} \cdot \vec{P} = \vec{\nabla} \cdot (\epsilon \cdot \vec{E})$
$\vec{\nabla} \cdot \vec{P} = (\vec{\nabla} \cdot \epsilon) \vec{E} + \epsilon (\vec{\nabla} \cdot \vec{E})$
...assume $\vec{\nabla} \cdot \epsilon = 0$ and substitute the electric field divergence:
$\vec{\nabla} \cdot \vec{P} = \frac{\epsilon}{\epsilon_0} \rho$
...idk what to do next, but the result needs to be...
$-\vec{\nabla} \cdot \vec{P} = \rho_{bound}$
TODO So where is this definition from? How is it derived?

$\vec{J}_{bound} = \vec{\nabla} \times \vec{M} + \partial_t \vec{P}$ (same, verify plz)

Macroscopic Maxwell Equations (Tells how a charge generates a field around it):
$\vec{\nabla} \cdot \vec{D} = \rho_{free}$
$\vec{\nabla} \cdot \vec{B} = 0$

$\partial_t \vec{B} + \vec{\nabla} \times \vec{E} = 0$
$\partial_t \vec{D} - \vec{\nabla} \times \vec{H} = -\vec{J}_{free}$

aka
$\partial_t \vec{B} + \vec{\nabla} \times (\epsilon^{-1} \cdot \vec{D}) = 0$
$\partial_t \vec{D} - \vec{\nabla} \times (\mu^{-1} \cdot \vec{B}) = -\vec{J}_{free}$

So now if I want to argue that $\partial_t \vec{E} = $ due to electric field steady state, can I still argue that $\partial_t \vec{P} = 0$ such that $\partial_t \vec{D} = 0$?
Only if so do we lose the 2nd time-derivative equation and are left with the 1st time-derivative equation... which is the same found in microscopic ideal MHD, except we would be now using a non-vacuum permittivity.