Hamiltonian for a Lagrangian with coupling 
I am dealing with the following Lagrangian density
  $$\mathscr{L}_{em}=
-\frac{1}{2}\rho\omega^2 u^2
+\frac{1}{2}\nabla u:\Sigma :\nabla u-\frac{1}{2}\nabla\phi\cdot\epsilon\cdot\nabla\phi+\nabla\phi\cdot P:\nabla u$$
  where $\rho,\omega\in\mathbb{R}^+$, $\Sigma_{ij,kl}=\Sigma_{ji,kl}=\Sigma_{ij,lk}=\Sigma_{kl,ij}$, $\epsilon_{ij}=\epsilon_{ji}$, $P_{ijk}=P_{ikj}$, $\phi$ is a scalar field and $ u \in\mathbb{R}^3$.
I need to compute the associated Hamiltonian density.

If it was just
$$\mathscr{L}_m=-\frac{1}{2}\rho\omega^2 u^2
+\frac{1}{2}\nabla u:\Sigma :\nabla u$$
by defining the momentum $\sigma_m=\Sigma:\nabla u$, and using the Legendre transform
$$\mathscr{H}=p\cdot\nabla q(q,p)-\mathscr{L}(q,p),$$
where $q$ are the field variables and $p$ the momentum, I obtain
$$\mathscr{H}_m=\frac{1}{2}\rho\omega^2 u^2+\frac{1}{2}\sigma_m:\Sigma^{-1} :\sigma_m.$$
Also for
$$\mathscr{L}_e=-\frac{1}{2}\nabla\phi\cdot\epsilon\cdot\nabla\phi$$
I can obtain
$$\mathscr{H}_e=-\frac{1}{2} d_e\cdot\epsilon^{-1}\cdot d_e$$
with $d_e=-\varepsilon\cdot\nabla\phi$.
But now, what about the Hamiltonian density for $\mathscr{L}_{em}$?
Can I write something like
$$\mathscr{H}_{em}=\frac{1}{2}\rho\omega^2 u^2
+\frac{1}{2}\sigma_m:\Sigma^{-1}:\sigma_m
-\frac{1}{2}d_e\cdot\epsilon^{-1}\cdot d_e\pm d_e\cdot Q:\sigma_m~?$$
Or must I rely on the introduction of the momentum
$$\sigma_{em}=\Sigma:\nabla u+P^T\cdot\nabla\phi.$$
$$d_{em}=-\varepsilon\cdot\nabla\phi+P:\nabla u~ ?$$
Who is the matrix $Q$?
Is something related to this Phys.SE post: Lagrangian and hamiltonian of interaction ?
I am new in the argument, but every suggestion is appreciated.
 A: In order to perform the (possibly singular) Legendre transformation, it is necessary to have information about pertinent rank conditions of the structure constants $\rho$, $\omega$, $\Sigma_{ij,k\ell}$, $\epsilon_{ij}$ and $P_{ijk}$. 
In this answer, we will sketch how the (possibly singular) Legendre transformation is performed in principle:


*

*We will use DeWitt's condensed notation to hide all spatial derivatives for simplicity.

*Assume that the Lagrangian density 
$$\tag{1}\mathcal{L}~=~\mathcal{L}_2+\mathcal{L}_1+\mathcal{L}_0$$ 
is a quadratic function of the velocities $\dot{\Phi}^A$ (=temporal derivatives of the fields).

*We may for later convenience redefine the fields 
$$\tag{2}\Phi^A~\longrightarrow ~\Phi^{\prime A}~=~R^A{}_B~\Phi^B.$$

*After possible redefinition (2) of the fields 
$$\tag{3}\Phi^A~=~\{\phi^{\alpha}, \ldots \},$$
we may assume that $\mathcal{L}_2$ is of the form 
$$\tag{4} \mathcal{L}_2~=~\frac{1}{2}\dot{\phi}^{\alpha}m_{\alpha\beta}\dot{\phi}^{\beta},$$
where the symmetric matrix $m_{\alpha\beta}$ is invertible.

*Let us for simplicity assume that $\mathcal{L}_1$ is quadratic in the variables $\Phi^B$.

*After possible redefinition (2) of the fields 
$$\tag{5}\Phi^A~=~\{\phi^{\alpha};z^I; \lambda^a\},$$
we may assume that $\mathcal{L}_1$ is of the form
$$\tag{6}\mathcal{L}_1~=~A_{\alpha}\dot{\phi}^{\alpha}
+\frac{1}{2}z^I\omega_{IJ}\dot{z}^J,$$
where the matrix $\omega_{IJ}$ is invertible, and $A_{\alpha}$ depends linearly on the fields $\Phi^B$.

*Define momenta
$$\tag{7}\pi_{\alpha}
~:=~\frac{\partial\mathcal{L}}{\partial \dot{\phi}^{\alpha}}
~=~m_{\alpha\beta}\dot{\phi}^{\beta}+A_{\alpha}.$$

*By possibly redefining the fields
$$\tag{8}z^I~\longrightarrow ~z^{\prime I}~=~r^I{}_J~z^J,$$
with
$$\tag{9}z^I~=~\{q^i;p_j\},$$
and possibly throwing away total time derivative terms, we may assume that
$$\tag{10}\frac{1}{2}z^I\omega_{IJ}\dot{z}^J~=~p_i \dot{q}^i.$$

*Introduce momenta $\rho_b$ to the auxiliary variables $\lambda^a$.

*The Hamiltonian density becomes
$$\tag{11} {\cal H}~=~\frac{1}{2}(\pi_{\alpha}-A_{\alpha})(m^{-1})^{\alpha\beta}(\pi_{\beta}-A_{\beta})-\mathcal{L}_0, $$
cf. the Faddeev-Jackiw method. See also e.g. this Phys.SE post.

*In the Hamiltonian formulation, $\{\phi^{\alpha};\pi_{\beta}\}$, $z^I =\{q^i;p_j\}$ and $\{\lambda^a;\rho_b\}$ are canonical variables.
A: Thanks to the procedure suggested by Qmechanic, I have clarified myself. I need just to invert the matrix $m$, since it has for the momenta
$\sigma_{em}=\Sigma :\nabla u+P^T\nabla \phi$
$d_{em}=P:\nabla u-\epsilon\cdot\nabla\phi$
or
$\left(\begin{array}{c}\sigma_{em}\\d\end{array}\right)=
\left(\begin{array}{cc}\Sigma & P^T\\ P & -\epsilon\end{array}\right)
\cdot\left(\begin{array}{c}\nabla u\\\nabla\phi\end{array}\right)$.
In my case, the matrix $m$, has tensor entries, but considering the symmetries for $\Sigma$, $\epsilon$ and $P$, it can be reduced to a 9x9 matrix $m'$ with scalar entries and the vector for the "velocities" has now just 9 entries. Then is matter to do the computation.
The procedure suggested by Qmechanic is quite general and I've very appreciated his/her suggestion, Thanks!
