# Whispering gallery modes

Written by I. S. Grudinin

In whispering gallery mode (WGM) resonators, which are typically spheres or disks, electromagnetic field is trapped by total internal reflections from the boundary. Modes of such resonators are distinguished by compact volume and record high quality factors (Q) in a broad range of frequencies.

Modern applications of such resonators include microwave and optical cavities for atomic clocks, cavity optomechanics, nonlinear and quantum optics. Analytical solutions for WG modes are only available for a limited number of idealized geometries, such as sphere or ellipsoid. Since resonator dimensions are typically much larger than optical wavelength, direct application of numerical 3D finite difference time domain (FDTD) or finite element methods (FEM) is not practical. It's possible to solve the vectorial wave equation by reducing it to a two dimensional case by taking axial symmetry into account.

Such reduction leads to a system of 3 equations to be solved in a 2D "$\rho-z$" section of a resonator. Please refer to OXBORROW2007 for a detailed derivation and to GRUDININ2012 for an example of using FreeFem++ to compute WGMs.

## Wave equation for the WGMs#

Since electric field is discontinuous on the surface of a dielectric and magnetic field is typically not, we derive our equations for the magnetic field. The electric field can be easily derived at a later stage from $\vec{E}=\frac{i}{\omega\epsilon_0}\hat{\epsilon}^{-1}\nabla\times\vec{H}$. Following a standard procedure starting with Maxwell equations we derive a wave equation in a single-axis anisotropic medium such as an optical crystal:

$$\label{eqn::wave} \nabla\times\left(\hat{\epsilon}^{-1}\nabla\times\vec{H}\right)-k_0^2\vec{H}-\alpha\nabla\left(\nabla\cdot\vec{H}\right)=0$$

Here $k_0=\omega/c$ is the wavenumber, $\alpha$ is the penalty term added to fight spurious FEM solutions. For anisotropic single-axis medium with $\partial\hat{\epsilon}/\partial\phi=0$ in cylindrical system of coordinates we have:

$$\hat{\epsilon}=\begin{pmatrix} \epsilon_{\rho} & 0 & 0 \\ 0 & \epsilon_{\rho} & 0 \\ 0 & 0 & \epsilon_z \end{pmatrix}. \nonumber$$

We now assume axial symmetry of our electromagnetic fields and insert an imaginary unity in front of the $H_{\phi}$ to allow all field components to be real numbers and also to account for the phase shift of this component $\vec{H}(\rho,\phi,z)=\left\{H_{\rho}(\rho,z),iH_{\phi}(\rho,z),H_z(\rho,z)\right\}\times e^{im\phi}$.

We write the wave equation \eqref{eqn::wave} explicitly in cylindrical coordinates, thus obtaining a set of three differential equations for the domain $\Omega$ given by the resonator's cross section and some space outside:

\begin{eqnarray} \label{eqn::system} A_1\{{H}_{\rho}^t,{H}_{\phi}^t,{H}_{z}^t\}&=&0\\ \nonumber A_2\{{H}_{\rho}^t,{H}_{\phi}^t,{H}_{z}^t\}&=&0\\ \nonumber A_3\{{H}_{\rho}^t,{H}_{\phi}^t,{H}_{z}^t\}&=&0 \end{eqnarray}

The numerical solutions of these equations and boundary conditions can be found with FreeFem++ if we write the system in the weak, or integral form.

## Weak formulation#

In general, to obtain the integral or "weak" statements equivalent to system \eqref{eqn::system} and boundary conditions we form a scalar dot product between an arbitrary magnetic field test function $\mathbf{H}^t=\{{H}_{\rho}^t,{H}_{\phi}^t,{H}_{z}^t\}$ and the components of our vectorial equation $A_1,A_2,A_3$, and integrate over the resonator's cross section domain $\Omega$ (and its boundary for the boundary conditions):

$$\int\limits_{\Omega}(H^t_{\rho}A_1+H^t_{\phi}A_2+H^t_{z}A_3)d\Omega$$

We can reduce the order of partial derivatives in this integral by using the Green's formula for integration by parts. For example:

$$\int\limits_{\Omega}H_z^t \frac{\partial^2 H_z}{\partial \rho^2 }d\Omega= -\int\limits_{\Omega}\frac{\partial H_z^t}{\partial \rho}\frac{\partial H_z}{\partial \rho }d\Omega+\oint H_z^t\frac{\partial H_z}{\partial \rho}n_{\rho}d\Gamma$$

Thus converting equations \eqref{eqn::system} we obtain a large expression for the weak form.

## A dielectric sphere example with FreeFem++#

We now compute the fundamental mode frequency for a fused silica sphere. The sphere is 36 micrometer in diameter, the refractive index is 1.46, the boundary condition is the magnetic wall (which can actually be omitted as it holds automatically).

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 // Parameters real radius = 36; //approximate radius of the cavity real yb = -10, yt = -yb; //window yb=bottom and yt=top coordinates real xl = radius-5, xr = radius+3; //window xl=left and xr=right coordinates real angle = asin((yt)/radius); //angle of the sphere segment to model in radians int Nm = 60; //number of mesh vertices per border real ne = 1.46; //n_e-extraordinary refractive index (root of permittivity parallel to z-axis, epara) real no = 1.46; //n_o-ordinary refractive index (root of permittivity orthogonal to z-axis, eorto) real nm = 1; //refractive index of surrounding medium (air) int nev = 4; // number of eigen values to find int M = 213; //azimuthal mode order ~ 2Pi*n*R/lambda real alpha = 1; //penalty term // Mesh border W1l(t=0, 1){x=xl+(radius*cos(angle)-xl)*(1-t); y=yt; label=1;} border W1r(t=0, 1){x=xr-(xr-radius*cos(angle))*(t); y=yt; label=1;} border W2(t=0, 1){x=xr; y=yb+(yt-yb)*t; label=1;} border W3l(t=0, 1){x=xl+(radius*cos(angle)-xl)*(t); y=yb; label=1;} border W3r(t=0, 1){x=xr-(xr-radius*cos(angle))*(1-t); y=yb; label=1;} border W4(t=0, 1){x=xl; y=yt-(yt-yb)*t; label=1;} border S(t=0, 1){x=radius*cos((t-0.5)*2*angle); y=radius*sin((t-0.5)*2*angle); label=2;} mesh Th = buildmesh(W1r(Nm/4) + W1l(Nm/4) + W4(Nm) + W3l(Nm/4) + W3r(Nm/4) + W2(Nm) + S(Nm)); plot(Th, WindowIndex=0); // Fespace fespace Ph(Th, P0); Ph reg = region; int ncav = reg(xl+1, 0); // cavity int nair = reg(xr-1, 0); //air Ph eorto = no^2*(region==ncav) + nm^2*(region==nair); //subdomains for epsilon values inside and outside the resonators Ph epara = ne^2*(region==ncav) + nm^2*(region==nair); //subdomains for epsilon values inside and outside the resonators //supplementary variables to store eigenvectors, defined on mesh Th with P2 elements - Largange quadratic. fespace Supp(Th, P2); Supp eHsqr; //3d vector FE space fespace Vh(Th, [P2, P2, P2]); Vh [Hr, Hphi, Hz], [vHr, vHphi, vHz]; //magnetic field components on Vh space and test functions vH // Macro //boundary condition macros macro EWall(Hr, Hphi, Hz) ( dy(Hr) - dx(Hz) + Hr*N.x + Hz*N.y - epara*(Hz*M - dy(Hphi)*x)*N.y + eorto*(Hphi - Hr*M+dx(Hphi)*x)*N.x) // macro MWall(Hr, Hphi, Hz) ( Hphi + Hz*N.x - Hr*N.y + epara*(Hz*M - dy(Hphi)*x)*N.x + eorto*(Hphi - Hr*M+dx(Hphi)*x)*N.y ) // // Problem real sigma =(M/(ne*radius))^2+2; // value of the shift (k^2), where the modes will be found varf b ([Hr, Hphi, Hz], [vHr, vHphi, vHz]) = int2d(Th)( x*(Hr*vHr+Hphi*vHphi+Hz*vHz) ) ; // OP = A - sigma B ; // the shifted matrix varf op ([Hr, Hphi, Hz], [vHr, vHphi, vHz])= int2d(Th)( ( (eorto*(vHphi*Hphi - M*(vHphi*Hr + Hphi*vHr) + M^2*vHr*Hr) + epara*M^2*vHz*Hz)/x //A/r + eorto*(dx(vHphi)*(Hphi - M*Hr) + dx(Hphi)*(vHphi - M*vHr)) - epara*M*(vHz*dy(Hphi) + Hz*dy(vHphi)) //B + x*(eorto*dx(vHphi)*dx(Hphi) + epara*((dx(vHz) - dy(vHr))*(dx(Hz) - dy(Hr)) + dy(vHphi)*dy(Hphi))) //C )/(eorto*epara) + alpha*( (vHr*Hr - M*(vHphi*Hr + Hphi*vHr) + M^2*vHphi*Hphi)/x //D/r + (dx(vHr) + dy(vHz))*(Hr - M*Hphi) + (vHr - M*vHphi)*(dx(Hr) + dy(Hz)) //E + x*(dx(vHr) + dy(vHz))*(dx(Hr) + dy(Hz)) //F ) -sigma*x*(vHr*Hr + vHphi*Hphi + vHz*Hz) ) //electric wall boundary condition on the boundary of computation domain +int1d(Th, 1)( EWall(Hr, Hphi, Hz)*EWall(vHr, vHphi, vHz) ) ; //setting sparce matrices and assigning the solver UMFPACK to solve eigenvalue problem matrix B = b(Vh, Vh, solver=UMFPACK); matrix OP = op(Vh, Vh, solver=UMFPACK); // Solve real[int] ev(nev); //to store the nev eigenvalue Vh[int] [eHr, eHphi, eHz](nev); //to store the nev eigenvector //calling ARPACK on sparce matrices with the assigned solver UMFPACK: int k = EigenValue(OP, B, sym=true, sigma=sigma, value=ev, vector=eHr, tol=1e-10, maxit=0, ncv=0); k = min(k, nev); //sometimes the number of converged eigen values //can be greater than nev //file to output mode values ofstream f("modes.txt"); //setting number of digits in the file output int nold = f.precision(11); // Plot & Save for (int i = 0; i < k; i++){ real lambda = 2*pi/sqrt(ev[i]); eHsqr = (sqrt(eHr[i]^2 + eHphi[i]^2 + eHz[i]^2)); //intensity from magnetic field components plot(eHsqr, WindowIndex=i, value=1, nbiso=20,LabelColors=1, aspectratio=1, cmm="Mode "+i+", lambda="+lambda+", F="+(299792.458/lambda)); f << "Mode "<< i << ", ka=" << sqrt(ev[i])*radius << endl; } 

## References#

[OXBORROW2007] OXBORROW, Mark. Traceable 2-D finite-element simulation of the whispering-gallery modes of axisymmetric electromagnetic resonators. IEEE Transactions on Microwave Theory and Techniques, 2007, vol. 55, no 6, p. 1209-1218.

[GRUDININ2012] GRUDININ, Ivan S. et YU, Nan. Finite-element modeling of coupled optical microdisk resonators for displacement sensing. JOSA B, 2012, vol. 29, no 11, p. 3010-3014.