Coupling of scattered TE mode to TM mode using vector Helmholtz equation

Question

I solve a scattered field computation problem using the frequency domain Maxwell's/Helmholtz PDEs. Particularly, I'm studying the behavior of light, which is essentially a plane wave propagating in a 3D domain, as it interacts with a scatterer characterized by a complex refractive index. This interaction can result in diffraction, scattering, absorption, or reflection of the incident light. In ideal conditions for a TE mode, we consider only the $E_y$ and there should not be any electric field component ($E_x=0$ and $E_z=0$) along the direction of propagation (z-direction). However, in practical situations, there may be imperfections/irregularities/resonances on the scatterer boundaries that can cause some field components to scatter in unintended directions. As a result, I have a leakage of the electric field $E_y$ and scattering of it into other modes or directions. So for TE mode, $E_x$ and $E_z$ are not zero anymore, especially on the material interfaces (see Figure).

Before solving the scattered field computation problem, expanding Maxwell's equations

$$ \begin{bmatrix} \frac{\partial}{\partial y} E_z - \frac{\partial}{\partial z} E_y\\ \frac{\partial}{\partial z} E_x - \frac{\partial}{\partial x} E_z\\ \frac{\partial}{\partial x} E_y - \frac{\partial}{\partial y} E_x\\ \end{bmatrix} = j \omega \mu \begin{bmatrix} H_x\\ H_y\\ H_z\\ \end{bmatrix} $$

$$ \begin{bmatrix} \frac{\partial}{\partial y} H_z - \frac{\partial}{\partial z} H_y\\ \frac{\partial}{\partial z} H_x - \frac{\partial}{\partial x} H_z\\ \frac{\partial}{\partial x} H_y - \frac{\partial}{\partial y} H_x\\ \end{bmatrix} = -j \omega \epsilon(x,y,z) \begin{bmatrix} E_x\\ E_y\\ E_z\\ \end{bmatrix} $$

I derive the set of equations for E-field:

$$ \frac{\partial^2}{\partial y^2} E_x + \frac{\partial^2}{\partial z^2} E_x - \frac{\partial}{\partial y}\frac{\partial}{\partial x} E_y - \frac{\partial}{\partial z}\frac{\partial}{\partial x} E_z + k_{0}^2 \varepsilon E_x= 0 $$

$$ \frac{\partial^2}{\partial z^2} E_y + \frac{\partial^2}{\partial x^2} E_y - \frac{\partial}{\partial z}\frac{\partial}{\partial y} E_z - \frac{\partial}{\partial x}\frac{\partial}{\partial y} E_x + k_{0}^2 \varepsilon E_y= 0 $$

$$ \frac{\partial^2}{\partial x^2} E_z + \frac{\partial^2}{\partial y^2} E_z - \frac{\partial}{\partial x}\frac{\partial}{\partial z} E_x - \frac{\partial}{\partial y}\frac{\partial}{\partial z} E_y + k_{0}^2 \varepsilon E_z = 0 $$

The residual is the result of the summation of all terms comprising the PDE. For the problem, specified by initial/boundary conditions, the residuals are to be zero at the points where the approximate solution is obtained. By minimizing the residual using solvers, I solve the PDE numerically and I could easily compute the light at different modes: $E_y$ at TE, and $E_x$, $E_z$ at TM. My current objective is to compute how much $E_y$ light was rescattered to $E_x$, $E_z$.

Problem: by setting $E_x$, $E_z$ initially to zero for TE mode, I exclude these variables from the PDEs and cannot compute them later. To address how much of the scattered TE mode light is converted into TM mode light, I should consider the vector Helmholtz equation and somehow determine the coupling of scattered TE mode to TM mode.

Your assistance in providing a proper set of coupled equations would be greatly appreciated. I hope this clarifies the problem, but if you have any further questions, please don't hesitate to inquire.

Answer 1

This is not an answer to your question, which I do not really understand, but it is too long for a comment.

Permittivity is position dependent: $\epsilon = \epsilon(x,y,z)$, and can be complex for lossy medium. You seem to have assumed no free charges $\nabla \cdot \mathbf D = 0.$

$$\nabla \times (\nabla \times \mathbf {E}) =\nabla \times (-j\mu \omega \mathbf H)=\omega ^2\mu\epsilon \mathbf E \tag{1}$$ For any $\mathbf U $ vector field:$$\nabla \times (\nabla \times \mathbf {U}) = \nabla(\nabla \cdot \mathbf U)-\nabla^2 \mathbf U\tag{2}$$ and from $$0=\nabla \cdot\mathbf D =\nabla (\epsilon \mathbf E) =\mathbf E\cdot \nabla \epsilon +\epsilon \nabla \cdot\mathbf E,$$ and then $$\nabla \cdot \mathbf E = -\mathbf E\cdot \tfrac{1}{\epsilon}\nabla\epsilon ,\tag{3}$$ also $$\nabla(\nabla \cdot \mathbf E) = - \nabla(\mathbf E\cdot \tfrac{1}{\epsilon}\nabla \epsilon)\tag{4}.$$ You get $\nabla^2 \mathbf E+\nabla(\mathbf E\cdot \tfrac{1}{\epsilon} \nabla \epsilon)+\omega ^2\mu\epsilon \mathbf E=0 \tag{5}.$ It seems to me that the 2nd term in Eq (5) is missing from your setup and instead you have cross derivatives. Note too that if your medium is piecewise continuous then in every "piece" $\nabla \epsilon = 0$ and (5) simplifies to the standard Helmholtz equation $\nabla^2 \mathbf E+\omega ^2\mu\epsilon \mathbf E=0$ as it should. This is easier to solve by matching boundary conditions directly at the interface between two such homogeneous pieces.