Determinant for a coupled fluctuation Lagrangian Lets consider a bosonic physical system in variables $t, x$ and $y(x)$ ($x$ dependent) with a classical Lagrangian $L$. To first order in fluctuations $x \to x+\xi_1$ and $y \to y+\xi_2$ the fluctuated action reads
$$S_{fl}=\int dt dx \sqrt{g} ~{\vec \xi}^T~{\tilde O}~{\vec \xi}$$
Where the differential operator ${\tilde O}$ is given by
$${\tilde O}=\begin{pmatrix}
\nabla^2+M_1(x) & -\frac{1}{\sqrt{g}}\partial_x +M_{12}  \\
\frac{1}{\sqrt{g}}\partial_x +M_{12} & \nabla^2+M_2(x) 
\end{pmatrix}$$
Following the standard procedure, one has to compute the functional determinant of  ${\tilde O}$ to obtain the effective action for the fluctuations. Now, since the system is coupled (off-diagonal elements are present), the computation of $\det[{\tilde O}]$ will be rather hard, if not impossible:
$$\det[{\tilde O}]=\det\begin{pmatrix}
\nabla^2+M_1(x) & -\frac{1}{\sqrt{g}}\partial_x +M_{12}  \\
\frac{1}{\sqrt{g}}\partial_x +M_{12} & \nabla^2+M_2(x) 
\end{pmatrix}\\=\det\left(\left(\nabla^2+M_1(x)\right)\left(\nabla^2+M_2(x) \right)-\left(-\frac{1}{\sqrt{g}}\partial_x +M_{12}\right)\left(\frac{1}{\sqrt{g}}\partial_x +M_{12}\right)\right)$$
However, formally, before attempting to compute the determinant, one can multiply out all the terms in the fluctuation Lagrangian $L_{fl}={\vec \xi}^T~{\tilde O}~{\vec \xi}$, partially integrate once in the off-diagonal terms and rearrange the fields $\xi_1 \xi_2$ such that if one puts all the terms back into matrix notation, one finds:
$$S_{fl}=\int dt dx \sqrt{g} ~{\vec \xi}^T~{\tilde O}~'~{\vec \xi}$$
where now
$${\tilde O}~'=\begin{pmatrix}
\nabla^2+M_1(x) & 0  \\
\frac{2}{\sqrt{g}}\partial_x +2M_{12} & \nabla^2+M_2(x) 
\end{pmatrix}$$
If one proceeds to compute the determinant now:
$$\det[{\tilde O}~']=\det\left(\left(\nabla^2+M_1(x)\right)\left(\nabla^2+M_2(x) \right)\right)\\=\det\left(\nabla^2+M_1(x)\right)\det\left(\nabla^2+M_2(x) \right)$$
one finds a much simpler problem to solve.
Now my question - is this "simplification" actually allowed, or did I miss some subtlety which forbids such a rearrangement of bosonic fields at the Lagrangian level?
 A: I asked my professor and in a discussion we came up with the 
following.
The process of establishing the effective action for a 
fluctuation Lagrangian to consist of the functional determinant 
of the initial differential operator involved, relies on the 
equality:
$$\det(A)=e^{Tr(\log(A))}$$
for a matrix $A$, which is only true for diagonalizable 
matrices. A generic matrix of the type
$$B=\begin{pmatrix}
b_{11} & 0  \\
b_{21} & b_{22}
\end{pmatrix}$$
cannot be diagonalized and therefore the procedure does not hold for these kind of matrices.
One must always choose a diagonalizable matrix structure for the differential operators to be able to apply the QFT machinery.
EDIT
Due to questions in the comments I decided to mention one further idea which forbids matrices $B$ as given above. Remember the basics of QFT on a discretized lattice. Usually the following relation is found for real scalar fields (use them as an example) by explicit evaluation of Gaussian integrals:
$$\int D\varphi ~e^{-\frac{1}{2}\varphi_i A_{ij}\varphi_j}=\frac{1}{\sqrt{\det A}}$$
Now, remember why we involve the notation $\det A$ and not some other symbol in the result. It stems from the fact that we recognize by inspection that the result is exactly a determinant of a symmetric matrix $A$. Since the Gaussian integrals are evaluated separately, in components, it is even unavoidable that if we use
$$B=\begin{pmatrix}
b_{11} & 0  \\
b_{21} & b_{22}
\end{pmatrix}$$
as our matrix, the result will come out to be $1/\sqrt{\det B'}$ with
$$B'=\begin{pmatrix}
b_{11} & \frac{1}{2}b_{21}  \\
\frac{1}{2}b_{21} & b_{22}
\end{pmatrix}$$
Therefore, if we want to skip the part of explicitly evaluating Gaussian integrals and proceed directly to evaluate a determinant, we have to be careful to involve a matrix as it actually emerges from the formalism, namely symmetric.
Hope this different argument helps.
