Oblique incidence over a stratified medium - the transfer matrix is not unitary?

As you’ve read on the title, my problem is about an oblique incidence problem over a stratified electrical medium.

I’ve made a quick sketch with the incidence on the first layer. Let’s suppose that I have a TE incidence, even if it is not important to my question.

I tried to solve the problem using the transfer matrices. For the first interface between layer 1 and 2, for example, I have:

$$
left[begin{matrix}
E_1^+ E_1^-
end{matrix}right]
= underbrace{frac{1}{t_{12}}left[ begin{matrix} 1 & r_{12} r_{12} & 1 end{matrix} right]}_{M_{12}} left[begin{matrix}
E_2^+ E_2^-
end{matrix}right]
$$

with

$$r_{12}, t_{12} in mathbb{C} quad mathrm{and} quad M_{12}in mathbb{C}^{2 times 2}$$

The coefficients $t_{12} $ and $r_{12}$ are computed using the Fresnel equations.

Obviously before I can add a second interface matrix, I should add a propagation matrix:

$$
P_2 = left[
begin{matrix} e^{mathrm{i} beta_2 ell_2} & 0 0 & e^{-mathrm{i} beta_2 ell_2} end{matrix}
right]
$$

Where $beta_2 $is the propagation constant in layer 2 in the axis orthogonal to the interfaces and $l_2$ is the thickness of layer 2.
For N layers I have:

$$
left[begin{matrix}
E_1^+ E_1^-
end{matrix}right] = underbrace{M_{12} P_2 M_{23} P_3 cdots M_{N-1, N}}_{M_mathrm{tot}}
left[begin{matrix}
E_N^+ E_N^-
end{matrix}right]
$$

I numerically computed:

$$
M_mathrm{tot} = left[ begin{matrix} m_{11} & m_{12} m_{21} & m_{22} end{matrix} right] = left[ begin{matrix} 1/t & r/t r/t & 1/r end{matrix} right] qquad mathrm{and I plotted how} |r|^2= left|frac{m_{21}}{m_{11}}right|^2 mathrm{ varies with} theta
$$

Where $|r|^2$ is the reflection coefficient.

This gave me nonsense results: I have for some angles a reflection coefficient greater than one.

I think that the main problem is that the propagation matrices are not unitary.

If I compute the determinant for $M_{12}$, I get:

$$
det(M_{12}) = frac{1}{t_{12}^2}(1-r_{12}^2),
$$

Using the definition $t = 1+r$, then:
$$
det(M_{12}) = frac{1-r_{12}^2}{(1+r_{12})^2} = frac{1-r_{12}}{1+r_{12}} = frac{1-(t_{12}-1)}{t_{12}} = frac{2}{t_{12}}-1.
$$

This means that the determinant of these matrices can be also greater than one with certain angles.

I am doing something wrong, but I can’t see it.

I’ve also thought to use the scattering coefficient S21 instead of transmissivity t, but the transfer matrix requires the use of Fresnel coefficients.