You are absolutely correct. This is because the real and imaginary parts of the field must satisfy the same wave equation. What are the ways to get coupling between real and imaginary parts of the field?
- If there are lossy materials present. This is easy to see. But since we are dealing with the lossless problem, we need to look for other means.
- Through boundary conditions. A port or radiation boundary condition looks like a resistive load, hence will couple the real and imaginary field components. (I have a post in this in this discussion showing this boundary condition). This is the case for scattering that you show in your last post
- The other way is in situations where hybrid modes are present (not pure TE or TM). This is your waveguide case. This happens because the boundary condition between materials sometimes requires the existence of complex modes.
Your waveguide eigenvalue system is indefinite, which means that under certain circumstances it may yield complex conjugate pair eigenvalues. This will happen in certain frequency ranges for certain geometries. Also, \gamma^2 is complex. Your assembled matrices and DoF vectors must reflect this fact when you construct them.