What's the most general form of Maxwell's stress tensor for EM fields in matter? In Griffith's Introduction to Electrodynamics we find a derivation of the relation
$$\mathbf f + \frac{\partial \mathbf g}{\partial t} = \nabla \cdot \overleftrightarrow{\mathbf T} ,$$
where $\mathbf f$ is the Lorentz force per unit volume, $\mathbf g$ is the electromagnetic linear momentum per unit volume, and $\overleftrightarrow{\mathbf T}$ is Maxwell's stress tensor. However, this derivation employs Maxwell's equations in a vacuum, or at least in media that exhibit no polarization or magnetization.
My aim is to derive a similar relation for Maxwell's tensor, but starting from the more general set of equations that include the fields $\mathbf D$, $\mathbf E$, $\mathbf B$, and $\mathbf H$, with the general relations $\mathbf D = \varepsilon_0 \mathbf E + \mathbf P$ and $\mathbf B = \mu_0(\mathbf H + \mathbf M)$. Starting again from the local expression of the Lorentz force on free charges and current densities
$$\mathbf f = \rho_0 \mathbf E + \mathbf J_0 \times\mathbf B, $$
I've managed to show that
$$ \mathbf f = -\frac{\partial \mathbf g}{\partial t} + (\nabla \cdot \mathbf D)\mathbf E + (\nabla \cdot \mathbf B) \mathbf H - \mathbf D \times (\nabla \times \mathbf E) - \mathbf B \times (\nabla \times \mathbf H), $$
but I cannot seem to be able to continue. Intuitively, I think I should get something of the form
$$T_{ij} = E_i D_j + H_i B_j - \frac 1 2 \delta_{ij} (\mathbf E \cdot \mathbf D + \mathbf H \cdot \mathbf B), $$
which does revert to Griffith's
$$T_{ij} = \varepsilon_0 E_i E_j + \frac 1 {\mu_0} B_i B_j - \frac 1 2 \delta_{ij} \left( \varepsilon_0 E^2 + \frac 1 {\mu_0} B^2\right)   $$
as soon as $\mathbf P$ and $\mathbf M$ both vanish, and also happens to contain the expression for the electromagnetic energy density
$$u = \frac 1 2 \mathbf E \cdot \mathbf D + \frac 1 2 \mathbf H \cdot \mathbf B. $$
Maybe I'm missing some relevant vector calculus identities, or relations between the fields. Otherwise, I guess the only way out is to express $\mathbf D$ and $\mathbf B$ in terms of $\mathbf E$ and $\mathbf H$ respectively, through the usual relations (valid in all homogeneous isotropic linear media)
$$\mathbf D = \varepsilon \mathbf E, \qquad \mathbf B = \mu \mathbf H, $$
with $\varepsilon = \varepsilon_0 \varepsilon_r = \varepsilon_0 (1+\chi_e)$ and $\mu = \mu_0 \mu_r = \mu_0(1+\chi_m)$, which would yield
$$T_{ij} = \varepsilon E_i E_j + \frac 1 {\mu} B_i B_j - \frac 1 2 \delta_{ij} \left( \varepsilon E^2 + \frac 1 {\mu} B^2\right), $$
but this means we lose generality, and we won't be able to treat e.g. ferroelectric or ferromagnetic materials, or anisotropic linear media.
Any thoughts?
 A: Your guessed answer is almost correct.  The right order must be
$T_{ij}=D_iE_j+B_iH_j$+...
A derivation with the necessary assumptions is in pp. 283-284 of my EM textbook.
A: The incomparable Griffiths answers your question himself (indirectly because he doesn't use tensors) in his paper with Hnizdo What's the Use of Bound Charge (arXiv:1506.02590).
They derive two results that are relevant to your question, which in the form of your equation are:
$$\mathbf f + \frac \partial {\partial t} \left( \varepsilon_0 \mathbf E \times \mathbf B \right) = \nabla \cdot \overleftrightarrow{\mathbf T}$$
and
$$\mathbf f_o + \frac \partial {\partial t} \left( \frac 1 {c^2} \mathbf E \times \mathbf H \right) = \nabla \cdot \overleftrightarrow{\mathbf T}$$
where
$$T_{ij} = E_i D_j + H_i B_j - \frac 1 2 \delta_{ij} \left[ \varepsilon_0 E^2 + \mu_0 \left( H^2 - M^2 \right) \right]$$
$\mathbf f$ is the rate of change of total momentum density and $\mathbf f_o$ is the ‘overt’ force density, which is equal to $\mathbf f$ less the rate of change of the ‘hidden’ momentum density:
$$\mathbf f_o = \mathbf f - \frac \partial {\partial t} \left( \frac 1 {c^2} \mathbf M \times \mathbf E \right)$$
A: Further to my answer above, I thought it would be interesting to answer the question for the stress-energy tensor also.
This turns out to be very controversial, and the answer depends critically on how the tensor is partitioned between matter and field components. I'll take the matter component to correspond to Griffiths and Hnizdo's overt force density and break the field component into three to make things clearer: the (standard) electromagnetic stress-energy tensor, a polarization stress-energy tensor, and a magnetization stress-energy tensor. They are made up as follows.
Electromagnetic stress-energy tensor
$$T^{00} = \frac 1 2 \left(\varepsilon_0 E^2 + \frac 1 {\mu_0} B^2 \right)$$
$$T^{i0} = T^{0j} = \frac 1 {\eta_0} \mathbf E \times \mathbf B$$
$$T^{ij} = - \varepsilon_0 E_i E_j - \frac 1 {\mu_0} B_i B_j + \frac 1 2 \left( \varepsilon_0 E^2 + \frac 1 {\mu_0} B^2 \right) \delta_{ij}$$
Polarization stress-energy tensor
$$T_P^{00} = 0$$
$$T_P^{i0} = T_P^{0j} = \mathbf 0$$
$$T_P^{ij} = - E_i P_j$$
Magnetization stress-energy tensor
$$T_M^{00} = - \mathbf M \cdot \mathbf B$$
$$T_M^{i0} = T_M^{0j} = \frac 1 c \mathbf M \times \mathbf E$$
$$T_M^{ij} = M_i B_j - \mathbf M \cdot \mathbf B \delta_{ij}$$
The overt force density is then given by
$$f_o^\mu = - \partial_\nu \left(T^{\mu \nu} + T_P^{\mu \nu} + T_M^{\mu \nu} \right)$$
which is conistent with my previous answer.
This is almost the same as de Groot and Suttorp's result (https://staff.fnwi.uva.nl/l.g.suttorp/articles/energy1.pdf), except that their $T_M^{00} = 0$, whereas mine relates to the energy stored in the current loops that make up the magnetic dipoles (in the same way that $T_M^{i0}$ and $T_M^{0j}$ relate to the hidden momentum of the current loops). However, other authors present significantly different results: caveat emptor!
A: Using the formula (https://en.wikipedia.org/wiki/Vector_calculus_identities#Cross_product_rule)
$$\mathbf{V}\times \left(\nabla \times \mathbf{W} \right)= \left(\nabla \mathbf{W} \right) \cdot \mathbf{V}-\left(\mathbf{V} \cdot \nabla \right) \mathbf{W}$$
That, put in the formula gives
$$ \mathbf{f}+ \frac{\partial \mathbf{g}}{\partial t}= \left(\nabla \cdot \mathbf{D} \right) \mathbf{E}+\left(\nabla \cdot \mathbf{B} \right) \mathbf{H} - \left(\nabla \mathbf{E} \right) \cdot \mathbf{D}+\left(\mathbf{D} \cdot \nabla \right) \mathbf{E} - \left(\nabla \mathbf{H} \right) \cdot \mathbf{B}+\left(\mathbf{B} \cdot \nabla \right) \mathbf{H}$$
Expressed in terms of components, limited just to the terms with $\mathbf{E}$ and $\mathbf{D}$,
$$f_i+\frac{\partial g_i}{\partial t}=\left(\frac{\partial D_j}{\partial x_j} \right) E_i+ D_j \frac{\partial E_i}{\partial x_j}-\left(\frac{\partial E_i}{\partial x_j} \right) D_i+...$$
Now, the first 2 terms can be expressed as
$$\frac{\partial}{\partial x_j} (E_i \, D_j)$$
but, to add something like $-\delta_{ij} \mathbf{E}\cdot \mathbf{D}$, we need to also have a term $-E_i \left(\frac{\partial D_i}{\partial x_j} \right)$. 
Honestly, I don't know how to fix this. 
By the way, $\mathbf{g}=\mathbf{D}\times \mathbf{B}$.
