I have been working on waveguide mode analysis using FEM in Mathematica for a week, but I haven't succeeded until now.

The optical fiber-like waveguide is featured with different refractive index in core and in clad, and the interface between the core and the clad should have the boundary condition of $D_⊥$ (the normal component of $\mathbf{D}$) and $\mathbf{E}_∥$ (the tangential component of $\mathbf{E}$) are continuous. But I don't know how to express this kind of boundary condition in Mma. I think this is of course different in Neumann, Dirichlet and Robin conditions.

The physical model is described below.

For Helmholtz equation for optical waveguide:

$$\nabla ^2 E(x,y,z)+\epsilon (\frac{2 \pi} {\lambda})^2 E(x, y, z)=0$$

Assuming that $$E(x,y,z)=E(x,y) e^{i \beta z}$$

We have

$$\nabla ^2 E(x,y)+\epsilon \left(\frac{2 \pi }{\lambda }\right)^2- \beta^2 E(x,y) = 0$$

The $\beta$ in this equation is to be solved and $\beta ^2$ can be considered as the eigen value of this eigen equation

$$\nabla ^2 E(x,y)+\epsilon \left(\frac{2 \pi }{\lambda }\right)^2 = \beta^2 E(x,y)$$

$\epsilon$ is different for core and cladding, i.e. , $\epsilon_\text{core}$ and $\epsilon_\text{clad}$, respectively.

The boundary conditions at the interface should be : (1) the tangential component of the $\mathbf{E}$, i.e. $\mathbf{E}_∥$, is continuous. (2) the normal component of the $\mathbf{D}$, i.e. $D_⊥$, is continuous, in which $\mathbf{D}=\epsilon \mathbf{E}$. In the cylindrical coordinates $(r, \theta, z)$, the boundary conditions at interface should be, $E_z$ and $E_\theta$ is continuous, and $D_r$ is continuous.

These conditions are my main concern when using FEM for the analysis of the eigenmode. Although they can be formulated easily in some special cases such as in rectangular or circular waveguide, but I'd like to try a more general form.

Here is my unsuccessful try.( Mma 12.0, Win 10)

To make the mesh points on the boundary, it can be used like this,

<< NDSolve`FEM`

r = 0.8;
outerCirclePoints = 
    With[{r = 2.}, 
      Table[{r Cos[θ], r Sin[θ]}, {θ, Range[0, 2 π, 0.05 π] // Most}]]; (* the outer circle  *)
innerCirclePoints = 
    With[{r = r}, 
      Table[{r Cos[θ], r Sin[θ]}, {θ, Range[0, 2 π, 0.08 π] // Most}]]; (* the inner circle *)

bmesh = ToBoundaryMesh[
      "Coordinates" -> Join[outerCirclePoints, innerCirclePoints], 
      "BoundaryElements" -> {LineElement[
                RotateLeft[Range[Length@outerCirclePoints], 1]] // 
              Partition[#, 2] &], 
            Riffle[Range[Length@outerCirclePoints + 1, 
                  Length@Join[outerCirclePoints, innerCirclePoints]], 
                  Range[Length@outerCirclePoints + 1, 
                   Length@Join[outerCirclePoints,innerCirclePoints]],1]] //Partition[#,2] &]}];                                                     
    mesh = ToElementMesh[bmesh];
{bmesh["Wireframe"], mesh["Wireframe"]}
 (* generate the boundary and element mesh, to make the mesh points \
on the outer and inner circles   *)

glass = 1.45^2; air = 1.; k0 = (2 π)/1.55;
ϵ[x_, y_] := If[x^2 + y^2 <= r^2, glass, air]

helm = \!\(\*SubsuperscriptBox[\(∇\), \({x, y}\), \(2\)]\(u[x,y]\)\) + ϵ[x, y]*k0^2*u[x, y];
boundary = DirichletCondition[u[x, y] == 0., True];


{vals, funs} = NDEigensystem[{helm, boundary}, u[x, y], {x, y} ∈ mesh, 1,Method -> {"Eigensystem" -> {"FEAST","Interval" -> {k0^2, glass* k0^2}}}];

 Table[Plot3D[funs[[i]], {x, y} ∈ mesh, PlotRange -> All, 
    PlotLabel -> vals[[i]]], {i, Length[vals]}]

Although the profile in the figure seems right, but the eigenvalue is not right, since I can check it using analytical solutions.

Edit 1

I notice their is a very closely related post here, where PML is employed. However, there are some bug there, and it couldnot properly run.

Are there some more examples? Thank you in advance.

Edit 2:

For analytical solution, I have to mention that it can be obtained by solving a transdental equation shown in the figure for a circular waveguide. The derivation of this equation can be found in papers, for example, https://www.osapublishing.org/oe/abstract.cfm?uri=oe-12-6-1025.

In addition, I can check the FEM results with other FEM tools like COMSOL.

Edit 3

I'm thankful to user21 for his kindness and patience. In this post, I will first give more information on the analytical solution, then I would present more clearly what I want to do with FEM.

Since there are too many equations, I'd like to post the revelant contents from textbooks (the appendix of “Photonics”, written by A. Yariv) about the analytical approch to the eigenmodes of optical fibers. I have underlined the most important equations in the pictures to make these stuffs to be understand more easily.

enter image description here

enter image description here

enter image description here

enter image description here

enter image description here

Therefore, for circular optical fibers, I just need to solve the transcedental equation of (B-11), that was shown in the figure in my original post.

Since FEM is a more general way to cope with eigensystems with diverse operators and boundary conditions in arbitrary calculation region, now I'd like to use the FEM in mathematica to get the eigenvalue of $\beta ^2$ of $\nabla ^2 E(x,y)+\epsilon \left(\frac{2 \pi }{\lambda }\right)^2 = \beta^2 E(x,y)$, corresponding to the equation (A-7) underlined. Note that the eigenvalue is the square of the transcedental equation root in (B-11) underlined.

The codes in the newest post by user21 improves a lot compared in my original one, but I think it it still not right, since the calcutation result from COMSOL is exactly equal to the analytical solution.

Edit 4:

The following is the code for analytical solution.

λ = 1.55; ρ = 0.8; (* fiber core radius*)
ncore = 1.45;(*fiber core index*)
nclad = 1;(*fiber cladding index*) 
s = 1;(*β for Subscript[HE, sm] or Subscript[EH, sm] mode*)

k = (2 π)/λ;  
V = k ρ (ncore^2 - nclad^2)^(1/2); (*Print["V=",V];*)
U = ρ (k^2 ncore^2 - β^2)^(1/2); W = ρ (-k^2 nclad^2 + β^2)^(1/2);

Subscript[L, 1] = D[BesselJ[s, y], y]/(y BesselJ[s, y]) + 
       D[BesselK[s, z], z]/(z BesselK[s, z]); 
Subscript[L, 2] = D[BesselJ[s, y], y]/(y BesselJ[s, y]) + 
       (nclad^2/ncore^2) (D[BesselK[s, z], z]/(z BesselK[s, z])); 

lft = Subscript[L, 1 ] Subscript[L, 2] /. {y -> U, z -> W};
rght = ((s β)/(k ncore))^2 (V/(U W))^4;
bv = FindRoot[
   lft == rght, {β, 
    Boole[V <= 2.405]*(k nclad + 10^-10) + 
     Boole[V > 2.405]*(k ncore - 10^-10), k nclad + 10^-10, 
    k ncore - 10^-10}];

Edit 5

Actually, some of the FEM model for fiber eigenmode analysis uses this expression to calculate the eigenfrequency of the modes. I'm sorry that I cann't understand this.

$\omega^2=\frac{\int\left[(\nabla\times H)^*{\varepsilon}(\nabla \times H)+\rho (\nabla \times H)^*(\nabla \times H)\right]\text{dx}\text{dy}}{\int \text{dxdy}H\mu H^*}$

I'm not sure it is an answer. But I'll post it here and waiting for elegant Mathematica codes.

All the above answers are very helpful. However, I'm afraid none of them is correct.

Here I find something that may be useful on the github. I hope this can work as a hint for elegant Mma codes.

I copied this theory part from here. It also included python codes.

I think models like this would be very helpful for non-expert users like me, and I don't know whether it is possible to include such a model in Mathematica, like in COMSOL. Actually, this is the real motivation for this question.

Theory (taken from the fenics book)

We are solving the Helmholtz equation:

$$\nabla \times \frac{1}{\mu_r}\nabla \times \boldsymbol{E} - k_0^2 \epsilon_r \boldsymbol{E}=0,$$ with the boundary condition of electric mirrors. $$k_0 = \frac{2\pi f_0}{c_0}$$

$$\epsilon_r = (n+\mathrm{i}k)^2$$

The functional is: $$F(\boldsymbol{E}) = \frac{1}{2} \int_\Omega \left[\frac{1}{\mu_r} \left(\nabla \times \boldsymbol{E}\right)\cdot \left(\nabla \times \boldsymbol{E}\right) - k^2_0 \epsilon_r \boldsymbol{E}\boldsymbol{E}\right]dx$$

In order to find it the vector is split in to the transverse and axial components.

$$\nabla = \nabla_t + \nabla_z,$$ with $$\nabla_t = \hat{e_x}\frac{\partial}{\partial x} + \hat{e_y}\frac{\partial}{\partial y}$$ and $$\nabla_z = \hat{e_z}\frac{\partial}{\partial z}$$ We also assume that the electric field can be split like: $$\boldsymbol{E}(x,y,z) = \left[E_t(x,y) + \hat{e_z} E_z(x,y) \right] e^{-\beta z}$$ where $\beta$ is the complex propagation constant, $$\boldsymbol{E}_t = \hat{e_x} E_x + \hat{e_y} E_y$$ and $$ \boldsymbol{E}_{axial} = \hat{e_z} E_z$$

By splitting the $\nabla$ and by substituting the above for the electric field the functional becomes: $$ F_{\beta}(\boldsymbol{E}) = \int_{\Omega} \frac{1}{\mu_r} \left(\nabla_t \times \boldsymbol{E}_{t,\beta}\right) \cdot \left(\nabla_t \times \boldsymbol{E}_{t,\beta}\right) -k^2_0 \epsilon_r E_{t,\beta}E_{t,\beta} \\ + \beta^2 \left[\frac{1}{\mu_r} \left(\nabla_t E_z+\beta E_{t,\gamma}\right)\left(\nabla_t E_z+\beta E_{t,\gamma}\right) - k^2_0 \epsilon_r E_z E_z dx \right]$$

Since the functional is built now the electric field components need to be discetized in order to be calculated. The axial component can be discretised by nodal basis functions (Lagrandge) but the transverse need to be by two dimensional curl-conforming vector basis functions(Nedelec).

The electric field in the nedelec space is described by: $$E_{t} = \sum_{i=1}^{N_{N}} (e_t)_i N_i$$

and the axial in lagrange elements as:

$$E_z = \sum_{i=1}^{N_{L}} (e_z)_i L_i$$ Where $N_N$ and $N_L$ are the number of nedelec and lagrange elements respectively. $N_i$, $L_i$ are the basis functions and $e_t$, $e_z$ are the equivalent coefficients.

While the cuttoff method is easier to implement it only prety much finds the cuttoff wavenumber (shown in fenics book and in other papers). This is not what is of interest. I am interested in calculating the propagation constants $\beta$ and by extension the effective indices and the Electric field that they equate to.

The matrix eqution to be solved is

$$\begin{bmatrix} A_{tt} &amp; 0 \\ 0 &amp; 0 \end{bmatrix}\begin{bmatrix} e_t \\ e_z \end{bmatrix} = -\beta^2 \begin{bmatrix} B_{tt} &amp; B_{tz} \\ B_{zt} &amp; B_{zz} \end{bmatrix}\begin{bmatrix} e_t \\ e_z \end{bmatrix} $$ Where $$A_{tt} = S_{tt} - k^2_0 T_{tt} = \int_{\Omega} \frac{1}{\mu_r} \left(\nabla_t \times N_i \right) \left(\nabla_t \times N_j \right)dx \\ - k^2_0 \int_{\omega} \epsilon_r N_i N_jdx$$

$$B_{zz} = S_{zz} - k^2_0 T_{zz} = \int_{\Omega} \frac{1}{\mu_r} \left(\nabla_t L_i \right) \left(\nabla_t L_j \right)dx \\ - k^2_0 \int_{\omega} \epsilon_r L_i L_j dx$$ $$ B_{tt} = \int_{\Omega} \frac{1}{\mu_r} N_i \cdot N_j dx$$$$ B_{tz} = \int_{\Omega} \frac{1}{\mu_r} N_i \cdot \nabla_t L_j dx$$$$ B_{zt} = \int_{\Omega} \frac{1}{\mu_r} \nabla_t L_i \cdot N_j dx$$


  • $\begingroup$ Comments are not for extended discussion; this conversation has been moved to chat. $\endgroup$
    – Kuba
    Commented Jun 3, 2019 at 6:31

3 Answers 3


There are three conditions when we want to get eigenfunctions in Cartesian coordinates, similar to eigenfunctions in cylindrical coordinates. The first is the correspondence of boundaries. The second is the azimuth number match, eg $l_1=l_2=0$.Third, the radius of the outer circle must meet the boundary condition. All three conditions are violated in the author’s code. I show how to find the eigenfunction with $\beta = 5.336$.

    << NDSolve`FEM`

r = 0.8; ne = 10; om = 0.0; kap = 1000;
reg = ImplicitRegion[x^2 + y^2 <= 2.0928^2, {x, y}]; f = 
 Function[{vertices, area}, 
  Block[{x, y}, {x, y} = Mean[vertices]; 
   If[x^2 + y^2 <= r^2, area > 0.001, area > 0.01]]];
mesh = ToElementMesh[reg, MeshRefinementFunction -> f];

glass = 1.45; air = 1.; k0 = (2 \[Pi])/1.55; b = 5;
n[R_] := ( .5*(1 - Tanh[kap*(R - r)])*(glass^2 - air^2) + air^2)*k0^2

helm = -Laplacian[u[x, y], {x, y}] - (b^2 + n[Sqrt[x^2 + y^2]])*
    u[x, y] + I*om*(x*D[u[x, y], y] - y*D[u[x, y], x]);
boundary = DirichletCondition[u[x, y] == 0, True];
{vals, funs} = 
    NDEigensystem[{helm, boundary}, u[x, y], {x, y} \[Element] mesh, 

Sqrt[Re[vals] + b^2]

(* {5.01271, 5.01285, 5.03804, 5.03825, 4.92974, 4.92969, \
5.28019, 5.28066, 5.33615, 5.60683}*)

Here we can see that the 9-th eigenvalue is equal to 5.33615, which corresponds to the desired $\beta = 5.336$. Figure 1 shows the mesh and eigenfunction along with the cylinder bounding the fiberglass.

{Show[ mesh["Wireframe"], 
  ContourPlot[x^2 + y^2 == r^2, {x, -1, 1}, {y, -1, 1}, 
   ColorFunction -> Hue]], 
 Show[Plot3D[Re[funs[[9]]], {x, y} \[Element] mesh, PlotRange -> All, 
     PlotLabel -> Sqrt[vals[[9]] + b^2], Mesh -> None, 
   ColorFunction -> Hue], 
  Graphics3D[{Gray, Opacity[.4], 
    Cylinder[{{0, 0, -1}, {0, 0, 1.}}, r]}]]}

Figure 1

Figure 2 shows the remaining functions with $l\ne 0$ and desired eigenfunction with $l=0$ Figure 2

To isolate monotone solutions in the clad with l = 1, we add to the Helmholtz operator (b^2 + l^2/(x^2 + y^2))*u[x, y] and choose eigenfunctions that fade out in the outer region what is achieved when b = I*Sqrt[glass]*k0. Figure 3 shows one of the eigenfunctions. In this case, the desired value $\beta = 5.336$ is achieved with increasing size of the clad. In fig. 4 shows the same eigenfunction with a 2-fold increase in the size of the region of integration.

<< NDSolve`FEM`
r = 0.8; ne = 10;  kap = 1000; l = 1;
reg = ImplicitRegion[x^2 + y^2 <= 2^2, {x, y}]; f = 
 Function[{vertices, area}, 
  Block[{x, y}, {x, y} = Mean[vertices]; 
   If[x^2 + y^2 <= r^2, area > 0.001, area > 0.01]]];
mesh = ToElementMesh[reg, MeshRefinementFunction -> f];

glass = 1.45; air = 1.; k0 = (2 \[Pi])/1.55; b = I*Sqrt[glass]*k0;
n[R_] := ( .5*(1 - Tanh[kap*(R - r)])*(glass^2 - air^2) + air^2)*k0^2

helm = -Laplacian[
     u[x, y], {x, y}] - (b^2 + n[Sqrt[x^2 + y^2]] + l^2/(x^2 + y^2))*
    u[x, y];
boundary = DirichletCondition[u[x, y] == 0, True];

{vals, funs} = 
    NDEigensystem[{helm, boundary}, u[x, y], {x, y} \[Element] mesh, 

Sqrt[vals + b^2]

(*{0. + 4.93777 I, 0. + 5.29335 I, 0. + 5.29463 I, 
 0. + 3.9743 I, 0. + 3.97351 I, 0. + 3.51044 I, 0. + 3.50924 I, 
 0. + 3.23389 I, 0. + 2.86891 I, 0. + 2.86774 I}*)
{Show[ mesh["Wireframe"], 
  ContourPlot[x^2 + y^2 == r^2, {x, -1, 1}, {y, -1, 1}, 
   ColorFunction -> Hue]], 
 Show[Plot3D[Im[funs[[3]]], {x, y} \[Element] mesh, PlotRange -> All, 
     PlotLabel -> Row[{"\[Beta] = ", Im[Sqrt[vals[[3]] + b^2]]}], 
   Mesh -> None, ColorFunction -> Hue], 
  Graphics3D[{Gray, Opacity[.4], 
    Cylinder[{{0, 0, -1}, {0, 0, 1.}}, r]}]]}

Table[Plot3D[Im[funs[[i]]], {x, y} \[Element] mesh, PlotRange -> All, 
    PlotLabel -> Sqrt[vals[[i]] + b^2], Mesh -> None, 
  ColorFunction -> Hue], {i, Length[vals]}]

Figure 3 Figure 4 And finally I give the best result that was obtained in this model with l = 1:

<< NDSolve`FEM`
r = 0.8; ne = 1; kap = 10000; l = 1;
reg = ImplicitRegion[x^2 + y^2 <= 5.3^2, {x, y}];
mesh = ToElementMesh[reg, 
   MeshRefinementFunction -> 
    Function[{vertices, area}, 
     area > 0.0004 (1 + 9 Norm[Mean[vertices]])]];

glass = 1.45; air = 1.; k0 = (2 \[Pi])/1.55; b = I*Sqrt[glass]*k0*1.1;
n[R_] := ( .5*(1 - Tanh[kap*(R - r)])*(glass^2 - air^2) + air^2)*k0^2

helm = -Laplacian[
     u[x, y], {x, y}] - (b^2 + n[Sqrt[x^2 + y^2]] + l^2/(x^2 + y^2))*
    u[x, y];
boundary = DirichletCondition[u[x, y] == 0, True];

{vals, funs} = 
    NDEigensystem[{helm, boundary}, u[x, y], {x, y} \[Element] mesh, 

{Show[ mesh["Wireframe"], 
  ContourPlot[x^2 + y^2 == r^2, {x, -1, 1}, {y, -1, 1}, 
   ColorFunction -> Hue]], 
 Show[Plot3D[Im[funs[[1]]], {x, y} \[Element] mesh, PlotRange -> All, 
     PlotLabel -> Row[{"\[Beta] = ", Im[Sqrt[vals[[1]] + b^2]]}], 
   Mesh -> None, ColorFunction -> Hue], 
  Graphics3D[{Gray, Opacity[.4], 
    Cylinder[{{0, 0, -1}, {0, 0, 1.}}, r]}]]}

Figure 5

  • 1
    $\begingroup$ @yulinlinyu Virtually nothing has changed. b is a calibration constant that is included in the definition of $\beta $. b is used instead of option Method -> {"Eigensystem" -> {"FEAST","Interval" -> {k0^2, glass* k0^2}}}. n[Sqrt[x^2 + y^2]] practically is the same as yours \[Epsilon][x, y]*k0^2 . $\endgroup$ Commented Jun 9, 2019 at 5:20
  • 1
    $\begingroup$ In addtion to the oscillatory behavior in the clad, the fundamental mode of optical fiber corresponding to l=1. Is it possible to show this? In expression (A-7), the Helmholtz equation is Laplacian[u[x,y],{x,y}]+(k^2-\beta ^2)u[x,y]=0. Therefore, I'm just wondering why you use b^2 in front of k^2, and how can you find the value 2.0928or 1.9615 in your former post? Thank you very much. $\endgroup$
    – yulinlinyu
    Commented Jun 9, 2019 at 14:50
  • 1
    $\begingroup$ can you explain it more? I get lost here. In eq(A-7)I post, there is no such term. I thought it came from the term of d/d[fai]. But all the first three terms equals to Laplacian[u[x,y],{x,y}]. I know you are right. In fact, I have tried to exclude this term, but the eigenvalue is far from the correct value although the eigenvector seems right. $\endgroup$
    – yulinlinyu
    Commented Jun 11, 2019 at 12:53
  • 1
    $\begingroup$ I also try to use a larger clad in your codes, but the eigenvalue persists be about 5.30 or something like that. It is a little bit smaller than the exact value of 5.33. However, the COMSOL can get 5.33 exactly. Can you comment this issue? $\endgroup$
    – yulinlinyu
    Commented Jun 11, 2019 at 12:57
  • 1
    $\begingroup$ I'm sorry, but you've forgotten to add the minus back when calculating $\beta$ i.e. it should be Sqrt[-(vals + b^2)] rather than Sqrt[vals + b^2]. And, I'm afraid it's just impossible to obtain the correct result only with Helmholtz equation. Please see my answer below for more details. $\endgroup$
    – xzczd
    Commented Jun 14, 2019 at 16:36

This is not a complete answer but a few more things to think about:

Generate the mesh:

glass = 1.45^2;
air = 1.;
k0 = (2 Pi)/1.55;
\[Epsilon][x_, y_] := 
 If[ElementMarker == 1, Evaluate[glass], Evaluate[air]]
mesh = ToElementMesh[Annulus[{0, 0}, {0.8, 2}], "RegionHoles" -> None,
    "RegionMarker" -> {{{0, 0}, 1}, {{3/2, 0}, 2}}(*,

enter image description here

Set up the equation and boundary condition only at the outer boundary (Using True will also make use of the internal boundary)

helm = Laplacian[u[x, y], {x, y}] + \[Epsilon][x, y]*k0^2*u[x, y];
boundary = DirichletCondition[u[x, y] == 0, x^2 + y^2 >= 2];

Note that there are more eigenvalues in the interval you request:

{vals, funs} = 
 NDEigensystem[{helm, boundary}, u[x, y], {x, y} \[Element] mesh, 3, 
  Method -> {"Eigensystem" -> {"FEAST", 
      "Interval" -> {k0^2, glass*k0^2}}}]

(* {{21.8177, 21.8179, 29.2875},... } *)

But you requested NDEigensystem to only return one.

Note there are also negative eigenvalues:

{vals, funs} = 
  NDEigensystem[{helm, boundary}, u[x, y], {x, y} \[Element] mesh, 40];

{-2.09566, -2.10176, 2.30094, 2.30241, -2.74458, -2.74508, -3.84583, \
-3.85689, 4.79376, 5.27311, 5.27699, 7.27767, 7.2782, -8.27121, \
-8.27829, 9.51751, 9.51881, -11.3267, -11.3395, -11.6234, -11.6335, \
-11.6822, 13.2792, 13.6627, 13.6649, -14.3816, -14.3887, -19.6736, \
-19.6883, -20.4762, -20.4798, -20.9379, -20.9583, -21.1027, -21.1095, \
21.8177, 21.8179, -28.4373, -28.4643, 29.2875}

You could not catch those with the FEAST interval specified.

Last, I am not sure if you are looking for lambda or lambda^2:


{4.39178, 4.41741, 5.29434, 5.30111, 7.5327, 7.53544, 14.7904, \
14.8756, 22.9802, 27.8057, 27.8467, 52.9645, 52.9722, 68.4129, \
68.5301, 90.5831, 90.6077, 128.295, 128.584, 135.102, 135.339, \
136.475, 176.338, 186.669, 186.73, 206.831, 207.035, 387.05, 387.627, \
419.275, 419.424, 438.397, 439.251, 445.322, 445.609, 476.013, \
476.022, 808.683, 810.215, 857.756}

The 27.8 value will get closer to 28 when you refine the mesh.

  • $\begingroup$ Thank you very much for your kindness and patience. But I'm afraid it is still problematic. Please see the new Edit 3 . $\endgroup$
    – yulinlinyu
    Commented Jun 4, 2019 at 8:32
  • $\begingroup$ @yulinlinyu, what boundary conditions did you apply in COMSOL? $\endgroup$
    – user21
    Commented Jun 4, 2019 at 8:54
  • $\begingroup$ @yulinlinyu, what exact refractive index did you use? $\endgroup$
    – user21
    Commented Jun 4, 2019 at 9:04
  • $\begingroup$ @yulinlinyu, what eigenvalue do you expect? $\endgroup$
    – user21
    Commented Jun 4, 2019 at 9:06
  • $\begingroup$ @yulinlinyu, can you post the analytical solution as Mathematica code? $\endgroup$
    – user21
    Commented Jun 4, 2019 at 9:07

Not a perfect answer, but I believe I've found the correct direction.


First of all, I'd like to emphasize the following:

  1. Helmholtz equation for electric field is deduced from Maxwell's equations in frequency domain assuming permittivity $\epsilon$ is constant.

  2. Piecewise constant is not constant.

In other words, simply solving Helmholtz equation in the whole domain of definition won't produce the correct result, even if the piecewise constant is approximated with a smooth function. Linking the 2 sub-domain of definition with proper interface condition (as shown in the screenshot given by OP) is a possible solution, but that's just troublesome. A better approach is to turn to a more general equation that's valid even at the interface (in the sense of limit). Then which equation should we use? Well, I'm not an expert of electromagnetism, but I found the one mentioned in Full-vectorial mode calculations by finite difference method by C.L.Xu seems to be a possible choice.

I'll repeat the deduction in this post for completeness. Starting from Maxwell's equations in frequency domain

$$\nabla \times \mathbf{E}=- j \omega \mu_0 \mathbf{H} \tag{1}$$ $$\nabla \times \mathbf{H}= j \omega n^2 \epsilon_0 \mathbf{E} \tag{2}$$

Take the curl of $(1)$ and substitute it into $(2)$, we have

$$\nabla \times \nabla \times \mathbf{E}- n^2 k^2 \mathbf{E}=0 \tag{3}$$

where $k=\omega/c$ and $c=1/\sqrt{\epsilon_0 \mu_0}$.

With the vector identity

$$ \nabla \times \left( \nabla \times \mathbf{A} \right) = \nabla(\nabla \cdot \mathbf{A}) - \nabla^{2}\mathbf{A} \tag{4}$$

$(3)$ becomes

$$\nabla^2 \mathbf{E}+n^2 k^2 \mathbf{E}=\nabla(\nabla \cdot \mathbf{E}) \tag{5}$$

Notice if $\epsilon$ is constant, $\nabla \cdot \mathbf{E}$ will be $0$ based on Gauss's law, thus $(5)$ will simplify to Helmholtz equation, but of course we can't do this here.

Then let's eliminate $E_z$. The transverse components of $(5)$ are

$$\nabla_t^2\mathbf{E}_t+n^2 k^2 \mathbf{E}_t=\nabla_t(\nabla_t \cdot \mathbf{E}_t+\frac{\partial E_z}{\partial z}) \tag{6}$$

where $\nabla_t=\left(\frac{\partial}{\partial x},\frac{\partial}{\partial y}\right)$ and $\mathbf{E}_t=(E_x,E_y)$ in Cartesian coordinate.

Since $n$ is $z$-invariant in waveguide i.e. $\partial n/\partial z=0$, Gauss's law

$$\nabla \cdot( n^2 \mathbf{E})=0 \tag{7}$$

can be transformed to

$$\frac{\partial E_z}{\partial z}=-\frac{1}{n^2} \nabla_t \cdot (n^2 \mathbf{E}_t) \tag{8}$$

Substitute $(8)$ and

$$\mathbf{E}(x,y,z)=\mathbf{E}(x,y)e^{-j \beta z} \tag{9}$$

into $(6)$, we have

$$\boxed{\nabla_t^2\mathbf{E}_t+n^2 k^2 \mathbf{E}_t - \nabla_t\left(\nabla_t \cdot \mathbf{E}_t -\frac{1}{n^2} \nabla_t \cdot \left(n^2 \mathbf{E}_t\right)\right) = \beta^2 \mathbf{E}_t} \tag{10}$$

The remaining work is just to solve the equation.

"FiniteElement"-based Approach

It's straightforward to solve equation $(10)$ with the built-in NDEigensystem. First, code the equations:

lap = Laplacian[#, {x, y}] &;
grad = Grad[#, {x, y}] &;
div = Div[#, {x, y}] &;
elst = e[#][x, y] & /@ Range[2];

lhs = With[{n2 = n2[x, y]}, 
   lap@elst + (n2 k^2) elst - grad[div@elst - 1/n2 div[n2 elst]]];

Next step is to approximate the piecewise constant with a smooth function:

r = 8/10;
glass = (145/100)^2; air = 1;
k = (2 π)/155 100;
appro = With[{k = 100}, ArcTan[k #]/Pi + 1/2 &];
n2 = (Function[{x, y}, #] &@(Simplify`PWToUnitStep@
       PiecewiseExpand@If[x^2 + y^2 <= r^2, glass, air] /. UnitStep -> appro))

Notice this step is necessary, or the solution won't converge to the desired result and you'll see the eigenvalue close to 29 again.

Finally, generate the mesh and solve for the eigenfunctions:

outer = 2;
mesh = ToElementMesh[Annulus[{0, 0}, {r, outer}], "RegionHoles" -> None, 
   "MaxBoundaryCellMeasure" -> 0.01];

{val, vec} = 
   NDEigensystem[{lhs /. para /. {e[1] -> e1, e[2] -> e2}}, {e1, e2}, {x, y} ∈ mesh, 6, 
     Method -> {"Eigensystem" -> {"Arnoldi", "Shift" -> k^2 glass}}]; // AbsoluteTiming
(* 13.5384 seconds, in a 8-core machine. *)
(* {19.8019, 19.8068, 20.1241, 21.7348, 28.351, 28.3528} *)

DensityPlot[#[[1]][x, y], {x, y} ∈ mesh, PlotRange -> All, PlotPoints -> 50] &@vec[[-1]]

enter image description here

Still, the eigenvalue isn't 5.336^2 == 28.4729, and the resulting graphic is slightly but apparently skrew, but this time we've at least been close to the analytic solution. Sadly I fail to improve the result further.

FDM-based Approach

In the "FiniteElement"-based approach, the derivative of $\epsilon$ is calculated symbolically so a very dense grid is needed to capture the sudden variation of $\epsilon$, and this may be a reason for the inaccuracy of the result, so I decide to turn to FDM, with which the symbolic derivation is avoided naturally.

We first introduce 3 intermediate variables to avoid symbolic derivation:

r = .8;
glass = (1.45)^2; air = 1.;
k = (2 π)/1.55;
n2 = Function[{x, y}, #] &@(Simplify`PWToUnitStep@
     PiecewiseExpand@If[x^2 + y^2 <= r^2, glass, air]);
With[{n2 = n2[x, y]}, lhs = lap@elst + (n2 k^2) elst - grad@term[x, y];
  termrhs = div@elst - 1/n2 div@{termx[x, y], termy[x, y]};
  {termxrhs, termyrhs} = n2 elst];

Then discretize the left hand side (lhs) based on FDM. I'll use pdetoae, and discretizing the system in Cartesian coordinate for simplicity:

difforder = 1; points = 400; L = 2; domain = {-L, L}; grid = Array[# &, points, domain];

n2[#, grid] & /@ grid // ArrayPlot

enter image description here

As we can see, the approximation for the waveguide isn't bad with the dense enough grid.

ptoafunc = pdetoae[Flatten[{e /@ {1, 2}, term, termx, termy}][x, y], {grid, grid}, 
del = #[[2 ;; -2]] &;

ae = del /@ del@# & /@ ptoafunc@lhs; // AbsoluteTiming
(* {26.967, Null} *)
{aetermrhs, aetermxrhs, aetermyrhs} = 
   ptoafunc@{termrhs, termxrhs, termyrhs}; // AbsoluteTiming
(* {24.0939, Null} *)
vare = Outer[e[#][#2, #3] &, Range@2, del@grid, del@grid, 1] // Flatten;

Block[{term, termx, termy},
  Evaluate@ptoafunc@Through[{term, termx, termy}[x, y]] = {aetermrhs, aetermxrhs, 
  {barray, marray} = CoefficientArrays[ae // Flatten, vare]; // AbsoluteTiming]
(* {58.6109, Null} *)
{val, vec} = 
   Eigensystem[marray, -6, Method -> {"Arnoldi", "Shift" -> k^2 glass}]; // AbsoluteTiming
(* {23.9645, Null}, in a 8-core machine. *)    
mat = ArrayReshape[#, {2, points - 2, points - 2}] & /@ vec;
 MapThread[ArrayPlot[#[[1]], PlotLabel -> Sqrt@#2, PlotRange -> All, 
    ColorFunction -> "AvocadoColors"] &, {mat, val}]

enter image description here

As we can see, the result is closer to 5.336, but once again, I fail to improve the result further. Simply make the grid denser or make L larger won't help. Perhaps the automatic discretization by pdetoae is too naive in this case and a better difference scheme is necessary.


  1. NDEigenSystem won't work without the replacement {e[1] -> e1, e[2] -> e2}. (e[1] and e[2] causes the warning NDEigensystem::baddep.) I'm not sure about the reason.

  2. "FEAST" method can't be used, otherwise the warning Eigensystem::nosymh pops up. I'm not sure about the reason.

  3. The utilization of Gauss's law in the deduction of $(10)$ seems to be critical. Actually one can still eliminate $E_z$ with $(1)$ and $(2)$ only, but the deduced equation just can't converge to the desired result. I guess this might be related to the observation that numeric algorithm that doesn't obey Gauss's law can be inaccurate. (Check Introduction section of this paper for more information. )

  4. If you insists on solving the problem with Helmholtz equation with interface condition, notice the interface conditions in the question aren't sufficient. One still needs another 3 conditions for $\mathbf{H}$ i.e. the tangential component of $\mathbf{H}$ is continuous across the surface if there's no surface current present.

  5. Dirichlet b.c.s are essentially set by del@grid in vare. Since the elements at the boundary of domain aren't included in vare, they will be treated as constant and moved to barray by CoefficientArray, which is equivalent to setting zero Dirichlet b.c.s.

  • $\begingroup$ I have not been able to follow the discussion, but I am curious how the correct weak formulation of these Helmholtz equations would look like. Actually, I expect the equation should read $\nabla \cdot \big( \varepsilon^{-1} \, \nabla E \big) + \lambda\, E = 0$, where $\varepsilon$ has a jump-type discontinuity along a spherical surface. This way, all the fuzz with virtual boundary conditions on an interior surface might become obsolete. One has only to take care that the stiffness matrix is assembled correctly. No? $\endgroup$ Commented Jun 14, 2019 at 16:30
  • 1
    $\begingroup$ @HenrikSchumacher I thought so, too, but this equation won't lead to the correct solution, and after playing with Maxwell's equations for quite a while, I noticed this seemingly reasonable equation is simply baseless…The continuity conditions for Maxwell equations are just not as simple as equation of heat conduction. $\endgroup$
    – xzczd
    Commented Jun 14, 2019 at 16:40
  • $\begingroup$ I see. Very interesting. $\endgroup$ Commented Jun 14, 2019 at 16:48
  • $\begingroup$ @xzczd Here you found a solution with l = 0. Now we need to find a solution with l = 1 and with $\beta =5.336$. $\endgroup$ Commented Jun 14, 2019 at 18:14
  • $\begingroup$ The evaluation can be faster by searching for a fewer eigenvalues if you adopt the trick of shifting the eigenvalues shown by Alex. $\endgroup$
    – yulinlinyu
    Commented Jun 17, 2019 at 6:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.