# Perturbation theory for Maxwell’s equations in anisotropic materials with
shifting boundaries

Di Yu  Institute of Fundamental and Frontier Sciences, University of
Electronic Science and Technology of China, Chengdu 610054, China  Yingcai
Honors College, University of Electronic Science and Technology of China,
Chengdu 611731, China  Xiaomin Lv  Institute of Fundamental and Frontier
Sciences, University of Electronic Science and Technology of China, Chengdu
610054, China  Boyu Fan  Institute of Fundamental and Frontier Sciences,
University of Electronic Science and Technology of China, Chengdu 610054,
China  Ju Gao  Institute of Fundamental and Frontier Sciences, University of
Electronic Science and Technology of China, Chengdu 610054, China  School of
Physics, University of Electronic Science and Technology of China, Chengdu
611731, China  Jingdao Tang  Institute of Fundamental and Frontier Sciences,
University of Electronic Science and Technology of China, Chengdu 610054,
China  Nan Xu  Institute of Fundamental and Frontier Sciences, University of
Electronic Science and Technology of China, Chengdu 610054, China  You Wang
Institute of Fundamental and Frontier Sciences, University of Electronic
Science and Technology of China, Chengdu 610054, China  Southwest Institute of
Technical Physics, Chengdu 610041, China  Haizhi Song  Institute of
Fundamental and Frontier Sciences, University of Electronic Science and
Technology of China, Chengdu 610054, China  Southwest Institute of Technical
Physics, Chengdu 610041, China  Qiang Zhou  Institute of Fundamental and
Frontier Sciences, University of Electronic Science and Technology of China,
Chengdu 610054, China  School of Optoelectronic Science and Engineering,
University of Electronic Science and Technology of China, Chengdu 610054,
China  CAS Key Laboratory of Quantum Information, University of Science and
Technology of China, Hefei 230026, China  Guangwei Deng<EMAIL_ADDRESS>Institute of Fundamental and Frontier Sciences, University of Electronic
Science and Technology of China, Chengdu 610054, China  CAS Key Laboratory of
Quantum Information, University of Science and Technology of China, Hefei
230026, China

###### Abstract

Perturbation theory is a kind of estimation method based on theorem of Taylor
expansion, and is useful to investigate electromagnetic solutions of small
changes. By considering a sharp boundary as a limit of smoothed systems,
previous study has solved the problem when applying standard perturbation
theory to Maxwell’s equations for small shifts in isotropic dielectric
interfaces[1]. However, when dealing with anisotropic materials, an
approximation is conducted and leads to an unsatisfactory error. Here we
develop a modified perturbation theory for small shifts in anisotropically
dielectric interfaces. By using optimized smoothing function for each
component of permittivity, we obtain a method to calculate the intrinsic
frequency shifts of anisotropic permittivity field when boundaries shift,
without approximation. Our method shows accurate results when calculating
eigenfrequency’s shifts in strong-anisotropy materials, and can be widely used
for small shifts in anisotropically dielectric interfaces.

Perturbation theory and Anisotropy and Shifting boundary

##  I Introduction

Perturbation theory is an approximation theory for solving eigenvalues and
eigenvectors of a characteristic equation based on theorem of Taylor
expansion. It is widely used in calculation in quantum mechanics,
electrodynamics and other physical fields. The classical electrodynamics
perturbation theory has great limitations[3, 4, 5], including the inability to
deal with material boundary movement[1, 6, 7, 8]. In fact, the formula for
calculating eigenfrequency shifts given by the classical electrodynamics
perturbation theory is only suitable for smooth permittivity constant field
perturbations and moving metal boundaries[2, 9]. Here, a smooth perturbation
on a field means that the overall field shifts smoothly with respect to the
change of arguments parameterizing the perturbation. If we use the calculation
formula of the intrinsic frequency shifts given by the electrodynamic
perturbation theory as a direct generalization without deliberateness, that
formula will yield a result including a ill-defined term, which should have
been defined at the boundary but unfortunately discontinuous near the
boundary[1]. If we simply take the value of the term on either side of the
boundary and substitute it into the result aforementioned, we will get a wrong
result, especially when there is a big dielectric contrast between the two
sides of the boundary[7, 8, 10]. This problem has been solved by Johnson, et
al., when the medium in the two sides are both isotropic[1]. They treat the
value of the permittivity at the boundary (as described by a step function) as
the limit of some continuous function. When dealing with anisotropic
materials, a common method is to approximately treat them as isotropic
materials in order to apply Johnson’s formula, such as the approximation on
lithium niobate photonic crystal nanocavities[12, 11]. However, there will be
non-negligible errors when using this approximation method, especially when
dealing with more complicated calculations, such as the opto-mechanical
coupling coefficients.

In this paper, we report a method for calculating the intrinsic frequency
shift of anisotropic materials with shifting boundaries. By comparing our
method with the previous one in some specific computing tasks, we find that
our method is more accurate. What’s more, we find Johnson’s formula can be
derived from our perturbation theory when restricted to the isotropic case.
Our method extends the application of Maxwell’s equations to strongly
anisotropic materials.

##  II Perturbation theory for the shift of isotropic material boundaries

There are two requirements to obtain the eigenvalue for the first order
perturbation theory: first, the perturbation is smooth; second, the operation
in the characteristic equation is Hermitian. The corresponding procedures and
requirements for the generalized eigenvalue equation are exactly the same.
With these two conditions the perturbation theory can be used to solve the
frequency shift of resonant electromagnetic field when the dielectric field is
perturbed — because the time-harmonic electromagnetic field satisfies the
following equation[2]:

$\nabla\times\nabla\times\left|E\right\rangle=\left(\frac{\omega}{c}\right)^{2}\left(\frac{\varepsilon}{\varepsilon_{0}}\right)\left|E\right\rangle,$
(1)

which is a generalized eigenvalue equation. Consequently, one should assume
that the perturbation of the dielectric field is smooth to satisfy the first
requirement mentioned above. In order to meet the second requirement, we also
need to assume that there is no dielectric loss, which leads to a real
permittivity. This procedure can also be performed in the case of anisotropic
materials, where the dielectric tensor is a real symmetric tensor under the
assumption of no loss. In this case, by using a Hermitian operation, we can
obtain similar results with isotropic materials, where the only change is that
the $\varepsilon$ here is a tensor.

Assuming the perturbation on a given permittivity field is parameterized with
an argument $q$, the derivative of intrinsic frequency with respect to $q$
writes

$\frac{d\omega}{dq}=-\frac{\omega^{(0)}}{2}\frac{\left\langle
E^{(0)}\vphantom{\frac{d\varepsilon}{dq}E^{(0)}}\right|\frac{d\varepsilon}{dq}\left|E^{(0)}\vphantom{E^{(0)}\frac{d\varepsilon}{dq}}\right\rangle}{\left\langle
E^{(0)}\vphantom{\varepsilon
E^{(0)}}\right|\varepsilon\left|E^{(0)}\vphantom{E^{(0)}\varepsilon}\right\rangle},$
(2)

where superscript $(0)$ indicates corresponding initial value without
perturbation.

We then consider the boundary of two materials that are denoted by their
permittivities $\varepsilon_{1}$ and $\varepsilon_{2}$ respectively, see
Fig.1. Obviously, there is an abrupt change of permittivity with the nearby
boundary moving, so the perturbation of permittivity is not smooth in this
case, which means it is illegal to apply general perturbation here.

$\varepsilon_{1}$$\varepsilon_{2}$$x$$y$$h$ Figure 1: A 2-dimensional
schematic diagram showing shifts of boundaries. The permittivities of the
materials on each side of the boundary are $\varepsilon_{1}$ and
$\varepsilon_{2}$ respectively. The local coordinate system is defined to be
an orthogonal coordinate system with $x$ axis perpendicular to the boundary.
$h$ is the displacement of the boundary in $x$ orientation. These definitions
are still available in corresponding 3-dimensional case.

To solve this problem, Johnson et al. [1] smooth the permittivity field around
the boundary to eliminate the abrupt change of the permittivity field
accompanying the movement of boundaries so that the change of the permittivity
field $\Delta\varepsilon_{q}(x,y,z)$ is smooth. Then the procedure of
perturbation method can be carried out and the approximated shift of
eigenfrequency can be obtained. Hence, to get the final result we need to
calculate the limitation of Eq. (2), namely to calculate the numerator
$\left\langle
E^{(0)}\vphantom{\frac{d\varepsilon}{dq}E^{(0)}}\right|\frac{d\varepsilon}{dq}\left|E^{(0)}\vphantom{E^{(0)}\frac{d\varepsilon}{dq}}\right\rangle$
and the denominator $\left\langle E^{(0)}\vphantom{\varepsilon
E^{(0)}}\right|\varepsilon\left|E^{(0)}\vphantom{E^{(0)}\varepsilon}\right\rangle$,
respectively. Typically, the denominator $\left\langle
E^{(0)}\vphantom{\varepsilon
E^{(0)}}\right|\varepsilon\left|E^{(0)}\vphantom{E^{(0)}\varepsilon}\right\rangle$
is easy to calculate, while it is much more difficult to calculate the
numerator $\left\langle
E^{(0)}\vphantom{\frac{d\varepsilon}{dq}E^{(0)}}\right|\frac{d\varepsilon}{dq}\left|E^{(0)}\vphantom{E^{(0)}\frac{d\varepsilon}{dq}}\right\rangle$.
Johnson et al. [1] calculated this limitation for isotropic materials to be

$\left\langle
E^{(0)}\vphantom{\frac{d\varepsilon}{dq}E^{(0)}}\right|\frac{d\varepsilon}{dq}\left|E^{(0)}\vphantom{E^{(0)}\frac{d\varepsilon}{dq}}\right\rangle=\int
dA\frac{dh}{dq}[\Delta\varepsilon_{12}|E_{||}^{(0)}|^{2}-\Delta(\varepsilon_{12}^{-1})|D_{\perp}^{(0)}|^{2}],$
(3)

where the superscript $(0)$ indicates the corresponding initial values without
perturbation, the integration domain is the whole boundary. Here
$\Delta\varepsilon_{12}=\varepsilon_{2}-\varepsilon_{1}$,
$\Delta(\varepsilon_{12}^{-1})=\varepsilon_{2}^{-1}-\varepsilon_{1}^{-1}$,
$|E_{||}^{(0)}|$ means the module of the parallel component for electric field
around the border, and $|D_{\perp}^{(0)}|$ denotes the perpendicular component
of electric displacement field at the boundary, and $h$ indicates the distance
of the movement of infinitesimal boundary patches.

However, this procedure, when conducted on anisotropic materials, will
encounter many difficulties because of the complexity of the algebra structure
of matrix. To our knowledge, there is no perturbation theory for the case of
moving boundaries of anisotropic materials up to now.

##  III The perturbation theory for the shifts of anisotropic material
boundaries

Considering an infinitesimal area of the boundary, we only need to consider
its movement in perpendicular orientation to the plane, because the parallel
movements of an infinitesimal area does not influence the distribution of the
permittivity field. The coordinate system is depicted in Fig.1, where we
define $x$ as the coordinate in the perpendicular orientation. $x=0$ at the
boundary, and $x$ is positive in the $\varepsilon_{2}$ side. A given
perturbation is imposed on the boundary. $h(q)$ is defined to be the change of
$x$ coordinate of an infinitesimal area during the perturbation parameterized
by argument $q$, and it depends on the location of the infinitesimal area. The
unperturbed permittivity distribution is described by
$\varepsilon(x)=\varepsilon_{1}+(\varepsilon_{2}-\varepsilon_{1})\Theta(x)$,
while the perturbed one by
$\varepsilon(x)=\varepsilon_{1}+(\varepsilon_{2}-\varepsilon_{1})\Theta(x-h)$,
where $\Theta(x)=0$ when $x<0$, and $\Theta(x)=1$ when $x>0$. As a result, the
perturbation of the permittivity field is
$(\varepsilon_{2}-\varepsilon_{1})(\Theta(x-h)-\Theta(x))$, as shown by the
red line in Fig.2. Because it is not smooth at $x=0$, there will be problems
in subsequent calculations.

$x$$\varepsilon$O$\varepsilon(x)$$\varepsilon_{s}(x)$ Figure 2: The smoothing
from $\varepsilon(x)$ to $\varepsilon_{s}(x)$, where $s$ parameterizes the
smoothing. Obviously, if $x$ changes a little, the change of
$\varepsilon_{s}(x)$ on each position is small, satisfying our smoothing
condition.

To solve this problem, we replace the original permittivity field
$\varepsilon(x)$ with a smoothed permittivity distribution function
$\varepsilon_{s}(x)$ and replace $\varepsilon(x-h)$ with
$\varepsilon_{s}(x-h)$, where the argument $s$ parameterizes the smoothing
method and $\varepsilon_{s}(x)\rightarrow\varepsilon(x)$ when $s\rightarrow
0$. In this way the perturbation of the permittivity field becomes smooth. The
specific form of $\varepsilon_{s}(x)$ is undetermined, but we can apply some
restrictions on it. Owing to that $\varepsilon_{1}$ and $\varepsilon_{2}$ are
both real, symmetric, positive definite tensor under lossless assumption, we
naturally stipulate the $\varepsilon_{s}(x)$ to be a real, symmetric, positive
definite tensor, which means $\varepsilon_{s}(x)$ is a Hermitian operation in
our calculation. Hence perturbation theory can be applied on the
characteristic equation

$\nabla\times\nabla\times\left|E\right\rangle=\left(\frac{\omega}{c}\right)^{2}\frac{\varepsilon_{s}}{\varepsilon_{0}}\left|E\right\rangle.$
(4)

The result of the perturbation theory is

$\frac{d\omega}{dq}=-\frac{\omega^{(0)}}{2}\frac{\left\langle
E^{(0)}\vphantom{\frac{d\varepsilon_{s}}{dq}E^{(0)}}\right|\frac{d\varepsilon_{s}}{dq}\left|E^{(0)}\vphantom{E^{(0)}\frac{d\varepsilon_{s}}{dq}}\right\rangle}{\left\langle
E^{(0)}\vphantom{\varepsilon_{s}E^{(0)}}\right|\varepsilon_{s}\left|E^{(0)}\vphantom{E^{(0)}\varepsilon_{s}}\right\rangle}.$
(5)

In order to get the true value of $\frac{d\omega}{dq}$ with no smoothing, we
need to calculate the limitation of $\left\langle
E^{(0)}\vphantom{\frac{d\varepsilon_{s}}{dq}E^{(0)}}\right|\frac{d\varepsilon_{s}}{dq}\left|E^{(0)}\vphantom{E^{(0)}\frac{d\varepsilon_{s}}{dq}}\right\rangle$
and $\left\langle
E^{(0)}\vphantom{\varepsilon_{s}E^{(0)}}\right|\varepsilon_{s}\left|E^{(0)}\vphantom{E^{(0)}\varepsilon_{s}}\right\rangle$
when $s\rightarrow 0$ according to Eq. 5. While $\left\langle
E^{(0)}\vphantom{\varepsilon_{s}E^{(0)}}\right|\varepsilon_{s}\left|E^{(0)}\vphantom{E^{(0)}\varepsilon_{s}}\right\rangle$
verges $\left\langle E^{(0)}\vphantom{\varepsilon
E^{(0)}}\right|\varepsilon\left|E^{(0)}\vphantom{E^{(0)}\varepsilon}\right\rangle$
with $s\rightarrow 0$, the limitation of $\left\langle
E^{(0)}\vphantom{\frac{d\varepsilon_{s}}{dq}E^{(0)}}\right|\frac{d\varepsilon_{s}}{dq}\left|E^{(0)}\vphantom{E^{(0)}\frac{d\varepsilon_{s}}{dq}}\right\rangle$
is elusive. To transform this term to an well-defined form at boundaries for
the calculation of the limitation, we now transform it to include only terms
continuous at the border. Because $D_{x},E_{y},E_{z}$ are all continuous
there, it is assumed that the limitation is a function of $D_{x},E_{y},E_{z}$,
where $x$ axis indicates the orientation of the shifts of the interface.

For convenience, $k_{ij}$ is employed to denote
$\frac{d\varepsilon_{s,ij}}{dq}\ (i,j=1,2,3)$, and $E_{ij}$ is employed to
denote the matrix with only $(i,j)$ element equal to 1 while others are $0$.
In the following discussion, we will omit the field quantity superscript
$(0)$.

To begin with, we have

$[D_{x},E_{y},E_{z}]^{\prime}=E_{11}\varepsilon_{s}\left|E\right\rangle+(E_{22}+E_{33})\left|E\right\rangle=A\left|E\right\rangle,$
(6)

where $A=E_{11}\varepsilon_{s}+E_{22}+E_{33}$, and ’′’ indicates
transposition. Here $A$ is reversible, because $\varepsilon_{s}$ is positive
definite tensor with a positive $(1,1)$ element. Hence

$A^{-1}[D_{x},E_{y},E_{z}]^{\prime}=\left|E\right\rangle.$ (7)

Here $\frac{d\varepsilon}{dh}\neq 0$ only in the domain where the permittivity
is smoothed, which is a small area around the boundary. $\Delta h_{s}$ is
defined so that the $x$ coordinate ranges from $-\Delta h_{s}$ to $\Delta
h_{s}$ in this area. $\Delta h_{s}$ describe the width of the area with
smoothed permittivity, so $\mathop{lim}\limits_{s\rightarrow 0}\Delta
h_{s}=0$. Hence, we have

$\left\langle
E\vphantom{\frac{d\varepsilon_{s}}{dq}E}\right|\frac{d\varepsilon_{s}}{dq}\left|E\vphantom{E\frac{d\varepsilon_{s}}{dq}}\right\rangle=\int\limits_{R^{3}}dV[E_{x},E_{y},E_{z}]^{*}\frac{d\varepsilon_{s}}{dq}[E_{x},E_{y},E_{z}]^{\prime}=\int\limits_{S}dA\frac{dh}{dq}\int_{-\Delta
h_{s}}^{\Delta
h_{s}}[E_{x},E_{y},E_{z}]^{*}\frac{d\varepsilon_{s}}{dh}[E_{x},E_{y},E_{z}]^{\prime}dx.$
(8)

Substituting Eq. (7) into Eq. (8), we have

$\left\langle
E\vphantom{\frac{d\varepsilon_{s}}{dq}E}\right|\frac{d\varepsilon_{s}}{dq}\left|E\vphantom{E\frac{d\varepsilon_{s}}{dq}}\right\rangle=\int\limits_{S}dA\frac{dh}{dq}\int_{-\Delta
h_{s}}^{\Delta
h_{s}}[D_{x},E_{y},E_{z}]^{*}A^{-1}\frac{d\varepsilon_{s}}{dh}A^{-1}[D_{x},E_{y},E_{z}]^{\prime}dx.$
(9)

When $s\rightarrow 0$,
$\frac{d\varepsilon_{s}}{dh}=(\varepsilon_{1}-\varepsilon_{2})\delta(x-h),$
where $\delta(x)$ is Dirac function. Because $[D_{x},E_{y},E_{z}]$ is
continuous at $x=0$, $[D_{x},E_{y},E_{z}]^{*}$ and
$[D_{x},E_{y},E_{z}]^{\prime}$ can be taken out of the inner integral, so

$\left\langle
E\vphantom{\frac{d\varepsilon_{s}}{dq}E}\right|\frac{d\varepsilon_{s}}{dq}\left|E\vphantom{E\frac{d\varepsilon_{s}}{dq}}\right\rangle=\int\limits_{S}dA\frac{dh}{dq}[D_{x},E_{y},E_{z}]^{*}\int_{-\Delta
h_{s}}^{\Delta
h_{s}}A^{-1}\frac{d\varepsilon_{s}}{dh}A^{-1}dx[D_{x},E_{y},E_{z}]^{\prime}.$
(10)

Now we can only need to ponder how to calculate the limitation of the
following term:

$\int_{-\Delta h_{s}}^{\Delta h_{s}}Xdx:=\int_{-\Delta h_{s}}^{\Delta
h_{s}}(A^{-1})^{\prime}\frac{d\varepsilon_{s}}{dh}A^{-1}dx,$ (11)

as other parts in Eq. 10 are already well-defined at the boundary.

Applying Cramer’s rule, $A^{-1}$ can be written as

$A^{-1}=\begin{bmatrix}\varepsilon_{s,11}^{-1}&-\frac{\varepsilon_{s,12}}{\varepsilon_{s,11}}&-\frac{\varepsilon_{s,13}}{\varepsilon_{s,11}}\\\
0&1&0\\\ 0&0&1\end{bmatrix}.$ (12)

By substituting Eq. (12) into Eq. (11), we have

$\displaystyle\int_{-\Delta h_{s}}^{\Delta h_{s}}Xdx=$ (13)
$\displaystyle\int_{-\Delta h_{s}}^{\Delta
h_{s}}\begin{bmatrix}\varepsilon_{s,11}^{-2}k_{11}&\varepsilon_{s,11}^{-1}\left(-\frac{\varepsilon_{s,12}}{\varepsilon_{s,11}}k_{11}+k_{12}\right)&\varepsilon_{s,11}^{-1}\left(-\frac{\varepsilon_{s,13}}{\varepsilon_{s,11}}k_{11}+k_{13}\right)\\\
\varepsilon_{s,11}^{-1}\left(-\frac{\varepsilon_{s,12}}{\varepsilon_{s,11}}k_{11}+k_{12}\right)&\varepsilon_{s,12}^{2}\varepsilon_{s,11}^{-2}k_{11}-2\frac{\varepsilon_{s,12}}{\varepsilon_{s,11}}k_{12}+k_{22}&\varepsilon_{s,12}\varepsilon_{s,13}\varepsilon_{s,11}^{-2}k_{11}-\frac{\varepsilon_{s,13}}{\varepsilon_{s,11}}k_{12}+k_{23}\\\
\varepsilon_{s,11}^{-1}\left(-\frac{\varepsilon_{s,13}}{\varepsilon_{s,11}}k_{11}+k_{13}\right)&\varepsilon_{s,12}\varepsilon_{s,13}\varepsilon_{s,11}^{-2}k_{11}-\frac{\varepsilon_{s,13}}{\varepsilon_{s,11}}k_{12}+k_{23}&\varepsilon_{s,13}^{2}\varepsilon_{s,11}^{-2}k_{11}-2\frac{\varepsilon_{s,13}}{\varepsilon_{s,11}}k_{13}+k_{33}\end{bmatrix}dx.$

This is a symmetric matrix, so we only need to consider 6 components of it. We
will calculate these components respectively.

In order to calculate the $(1,1)$ element, we consider the smoothing method s.
t.

$\varepsilon_{s,11}^{-1}=\int_{-\infty}^{-\infty}g_{s}(x-x^{\prime})\varepsilon_{11}^{-1}dx^{\prime},$
(14)

where $g_{s}(x)$ is an even smoothing function, and it is non-zero only within
a small range with $x$ in $(-\Delta h_{gs},\Delta h_{gs})$. A smoothing
function $S(x)$ is a real function that intensively distributes around $x=0$
and has a unit integral on the whold $x$ axis
$\int_{-\infty}^{-\infty}S(x)dx=1$ and verges Dirac function when the related
parameter verges 0. According to Eq. (14), we have
$k_{11}=-\varepsilon_{s,11}^{2}\Delta(\varepsilon_{11}^{-1})g_{s}(x-h)$, and
$\int_{-\Delta h_{s}}^{\Delta h_{s}}X_{11}dx=-\Delta(\varepsilon_{11}^{-1})$
for any value of $s$, where
$\Delta(\varepsilon_{11}^{-1})=\varepsilon_{1,11}^{-1}-\varepsilon_{2,11}^{-1}$,
so $\mathop{lim}\limits_{s\rightarrow 0}\int_{-\Delta h_{s}}^{\Delta
h_{s}}X_{11}dx=-\Delta(\varepsilon_{11}^{-1})$.

$f_{s}(x)$ is another smoothing function that is non-zero only within a small
range with x in $(-\Delta h_{fs},\Delta h_{fs})$. An extra limitation imposed
the $f_{s}(x)$ is that $f_{s}(x)$ verges Dirac function at a lower speed than
$g_{s}(x)$ does when $s\rightarrow 0$, which means $\Delta h_{gs}<<\Delta
h_{fs}$ for the same $s$. Taking advantage of $f_{s}(x)$, we construct the
following smoothing method

$x$$y$$O$$g_{s}(x)$$f_{s}(x)$$f_{s0}(x)$ Figure 3: A schematic for different
smoothing functions. A smoothing function $S(x)$ is a real function that
intensively distributes around $x=0$, has a unit integral on the whold $x$
axis $\int_{-\infty}^{-\infty}S(x)dx=1$, and verges Dirac function when the
related parameter verges 0. A possible form of $f_{s}(x)$ and $g_{s}(x)$ for
the same $s$ is shown in this graph, where $g_{s}(x)$ is much sharper than
$f_{s}(x)$, satisfying the condition we impose. $f_{s0}(x)$ is another
possible form of $f_{s}(x)$ corresponding to rectangular smoothing method, and
the parameter s here may represent the width of the rectangle.

$\varepsilon_{s,ij}=\int_{-\infty}^{\infty}f_{s}(x-x^{\prime})\varepsilon_{ij}dx^{\prime}.$
(15)

For convenience, We define $k_{ij}=\Delta\varepsilon_{ij}f_{s}(x-h)$, where
$i,j\in{1,2,3}$, and
$\Delta\varepsilon_{ij}=\varepsilon_{1,ij}-\varepsilon_{2,ij}$. Substituting
these equations into $X_{12},X_{13}$ and considering the parity of $f_{s}(x)$
and $g_{s}(x)$, we have

$\displaystyle\mathop{lim}\limits_{s\rightarrow 0}\int_{-\Delta h_{s}}^{\Delta
h_{s}}X_{1k}dx=\mathop{lim}\limits_{s\rightarrow 0}\int_{-\Delta
h_{s}}^{\Delta
h_{s}}X_{k1}dx=\frac{\varepsilon_{1,1k}}{\varepsilon_{1,11}}-\frac{\varepsilon_{2,1k}}{\varepsilon_{2,11}}\
(k=2,3).$ (16)

For the remaining parts $\int_{-\Delta h_{s}}^{\Delta h_{s}}X_{ij},\ i,j=2,3$,
each one consists of three terms in terms of Eq. 13. The last term o which is
easy to calculate: $\int_{-\Delta h_{s}}^{\Delta
h_{s}}k_{ij}dx=\Delta\varepsilon_{ij}$. As for the first term $\int_{-\Delta
h_{s}}^{\Delta
h_{s}}\varepsilon_{s,1i}\varepsilon_{s,1j}\varepsilon_{s,11}^{-2}k_{11}dx$, we
can replace $\varepsilon_{s,12}$ and $\varepsilon_{s,13}$ by
$\frac{\varepsilon_{1,12}+\varepsilon_{2,12}}{2}$ and
$\frac{\varepsilon_{1,13}+\varepsilon_{2,13}}{2}$, respectively, because the
non-zero region of $\varepsilon_{s,11}^{-2}k_{11}$ is much smaller than that
of $\varepsilon_{s,12},\varepsilon_{s,13}$. Finally we turn to the second term

$\frac{\varepsilon_{s,1i}}{\varepsilon_{s,11}}k_{1j}\ (i,j=2,3).$ (17)

Owing to the high convergence rate of $g_{s}(x)$, $\varepsilon_{s,11}$ can be
seen as a piecewise function taking a uniform value when $x<0$, while taking
the other uniform value when $x>0$. Whereupon

$\int_{-\Delta h_{s}}^{\Delta
h_{s}}\frac{\varepsilon_{s,12}}{\varepsilon_{s,11}}k_{12}dx=\frac{1}{\varepsilon_{1,11}}\int_{-\Delta
h_{s}}^{0}\varepsilon_{s,12}k_{12}dx+\frac{1}{\varepsilon_{2,11}}\int_{0}^{\Delta
h_{s}}\varepsilon_{s,12}k_{12}dx.$ (18)

Noting that $(-\Delta h_{fs},\Delta h_{fs})=(-\Delta h_{s},\Delta h_{s})$, we
have $\varepsilon_{s,12}|_{x=-\Delta h_{s}}=\varepsilon_{1,12}$,
$\varepsilon_{s,12}|_{x=0}=\frac{\varepsilon_{1,12}+\varepsilon_{2,12}}{2}$,
and $\varepsilon_{s,12}|_{x=\Delta h_{s}}=\varepsilon_{2,12}$. Then it is
reasonable for us to approximate
$\varepsilon_{s,12}\approx\frac{3}{4}\varepsilon_{1,12}+\frac{1}{4}\varepsilon_{2,12}$,
when$x<h$; $\frac{1}{4}\varepsilon_{1,12}+\frac{3}{4}\varepsilon_{2,12}$, when
$x>h$. Substituting these into Eq. (18), we have

$\mathop{lim}\limits_{s\rightarrow 0}\int_{-\Delta h_{s}}^{\Delta
h_{s}}\frac{\varepsilon_{s,12}}{\varepsilon_{s,11}}k_{12}dx\approx\frac{1}{8}\left[\frac{3\varepsilon_{1,12}+\varepsilon_{2,12}}{\varepsilon_{1,11}}+\frac{\varepsilon_{1,12}+3\varepsilon_{2,12}}{\varepsilon_{2,11}}\right](\varepsilon_{1,12}-\varepsilon_{2,12}).$
(19)

Many specific choice of the smoothing function can satisfy the equation above
accurately, so this step would not introduce errors. For example, we can
choose rectangular impulse function $f_{s0}(x)$, then we have

$\mathop{lim}\limits_{s\rightarrow 0}\int_{-\Delta h_{s}}^{\Delta
h_{s}}\frac{\varepsilon_{s,12}}{\varepsilon_{s,11}}k_{12}dx=\frac{1}{8}\left[\frac{3\varepsilon_{1,12}+\varepsilon_{2,12}}{\varepsilon_{1,11}}+\frac{\varepsilon_{1,12}+3\varepsilon_{2,12}}{\varepsilon_{2,11}}\right](\varepsilon_{1,12}-\varepsilon_{2,12}).$
(20)

Similarly, for other three terms in Eq. (17), we have

$\displaystyle\mathop{lim}\limits_{s\rightarrow 0}\int_{-\Delta h_{s}}^{\Delta
h_{s}}\frac{\varepsilon_{s,1i}}{\varepsilon_{s,11}}k_{1j}dx=\frac{1}{8}\left[\frac{3\varepsilon_{1,1i}+\varepsilon_{2,1i}}{\varepsilon_{1,11}}+\frac{\varepsilon_{1,1i}+3\varepsilon_{2,1i}}{\varepsilon_{2,11}}\right](\varepsilon_{1,1j}-\varepsilon_{2,1j})\
(i,j=2,3).$ (21)

By substituting the four formulas above into $\int_{-\Delta h_{s}}^{\Delta
h_{s}}Xdx$, the final components are calculated:

$\displaystyle\mathop{lim}\limits_{s\rightarrow 0}\int_{-\Delta h_{s}}^{\Delta
h_{s}}X_{ij}dx=\Delta\left(\varepsilon_{ij}-\frac{\varepsilon_{1i}\varepsilon_{1j}}{\varepsilon_{11}}\right)\
(i,j=2,3).$ (22)

To summarize, we transform Eq .13 to be

$\mathop{lim}\limits_{s\rightarrow 0}\int_{-\Delta h_{s}}^{\Delta
h_{s}}Xdx=\Delta\varepsilon_{11}^{-1}\begin{bmatrix}-1&\varepsilon_{12}&\varepsilon_{13}\\\
\varepsilon_{12}&\varepsilon_{11}\varepsilon_{22}-\varepsilon_{12}^{2}&\varepsilon_{11}\varepsilon_{23}-\varepsilon_{12}\varepsilon_{13}\\\
\varepsilon_{13}&\varepsilon_{11}\varepsilon_{23}-\varepsilon_{12}\varepsilon_{13}&\varepsilon_{11}\varepsilon_{33}-\varepsilon_{13}^{2}\end{bmatrix}=\xi$
(23)

Combining this equation with Eq. (10), we derive a perturbation theory for
calculating the eigenfrequency of anisotropic materials when shifts of
boundaries occur.

##  IV Simulation Verification

In this section, we use the Finite Element Method (FEM) to analyze the
accuracy of our formula, in comparison with the previous isotropic-material
approximation scheme.

xyz100 nm100 nm100 nmO

Figure 4: (a) The cube resonant cavity with the length of each edge L=100 nm.
The shifts of boundaries are imposed on the blue color side. The red arrow
indicates the orientation of the shifts. Definition of Euler angles $\alpha,\
\beta,$ and $\gamma$ is shown in (b), where the axes $X$, $Y$, and $Z$
indicate three principal axes of the index ellipsoid of the cavity material.

We consider a resonant cavity and calculate its eigenfrequency derivative as a
function of deformation through both our and previous methods. Here we choose
a cubic cavity with the length of each edge $L=100\ nm$ surrounded by air, as
depicted in Fig.4. The computational domain is a cube with the length of each
edge to be $5L$, and all sides of the computational domain cube are parallel
to the corresponding sides of the resonant cavity cube, while the cavity
locates in the center of the cube. We set the first-order scattering boundary
condition on all interfaces of the computational domain. Dielectric loss is
ignored and the refractive indexes of the cavity material are set to be
$n_{o}=10,\ n_{e}=9$. With such high refractive indexes, we can tighten the
distribution of electromagnetic field so as to reduce the amount of
calculation. As a comparison, we set the refractive index of the cavity to be
$n=9.5$ for the isotropic approximation situation, following a recent work on
Lithium Niobate [12]. The Euler angles $\alpha$, $\beta$, $\gamma$
corresponding to the rotation angles of three spindles $X$, $Y$, and $Z$ of
the index ellipsoid are defined with basis vectors in the space coordinate
system (see Fig.4). The displacement imposed on the blue side $\Delta
x_{q}(y,z)=\frac{q}{2\pi{\sigma}^{2}}exp\left(-\frac{y^{\prime 2}+z^{\prime
2}}{2\sigma^{2}}\right)$ is perpendicular to that side of the cavity, where
$\sigma=20,\ y^{\prime}=y-50,\ z^{\prime}=z-50$, and the orientation of the
displacement is indicated by the red arrow in Fig.4. (Here we omit the common
unit of ’nm’.) To simplify, here we investigate the case where
$\beta=\gamma=0$, $\alpha\in[0,\frac{\pi}{2}]$ and the frequency of the
investigated mode is about 360 THz. The parameter of the perturbation is set
as $q=1000$, when the maximum shift of 0.4 nm occurs at the center of the blue
side.

In order to compare the accuracy between our method and previous isotropic-
approximation scheme, we first calculate the eigenfrequency derivative as a
function of deformation $\frac{d\omega}{dq}$ by numerical method, which will
be considered as the ‘true’ value. Then we calculate $\frac{d\omega}{dq}$ by
our ’anisotropic’ method and previous isotropic-approximation scheme,
respectively. Compared with the ’true’ value, the estimated relative errors
are shown in Fig.5, as a function of $\alpha$. It is shown that our method
yields much smaller errors in most cases.

$0$$20$$40$$60$$80$$-10$$-5$$0$$5$$10$$\alpha\ $(deg)relative error of
$-\frac{d\omega}{dq}$ (%)

Figure 5: Estimated relative errors for $\frac{d\omega}{dq}$, as a function of
$\alpha$. Here $\beta=\gamma=0$. The results of our anisotropic perturbation
theory and previous isotropic method are in red and blue, respectively.

s

##  V Discussion and conclusion

In Fig.3, we have shown three smoothing functions, and $g_{s}(x)$ and
$f_{s0}(x)$ are used for smoothing different components. Why not directly
choose a rectangular smoothing function $f_{s0}(x)$ to simplify the
derivation? Actually, by using a rectangular smoothing function, one can
derive the same result. However, it seems that there is an assumption that all
physical smoothing methods will lead to the same result. Conversely, our
derivation does not rely on such an assumption.

Here, ’physical’ means that the smoothing should have physical meanings. For
example, the imposed restrictions on the smoothing that the smoothed
permittivity tensor field should be symmetric positive-defined means that the
smoothed permittivity field may correspond to some real matter.

A few factors that would cause frequency shifts are not taken into account,
such as contribution from the edges. Johnson et al. [1] have demonstrated that
the influence from factors except the shifts of sides are negligible for
1-order approximation, which guarantees the accuracy of our results.

Our perturbation theory incorporates Johnson’s scheme. In the case of
isotropic permittivity, if we denote $\varepsilon_{1}-\varepsilon_{2}$ by
$\Delta\varepsilon$, then it will yield $\xi=\Delta
Diag\\{-\varepsilon^{-1},\varepsilon,\varepsilon\\}$. In this way our formula
is the same as Johnson’s result[1] under isotropic conditions.

Our anisotropic perturbation theory can help in many situations. For example,
Lithium Niobate optomechanical crystals have several promising applications
like ultra-low-power modulators and quantum information processing, and a key
step of designing these crystals is to calculate the single-phonon
optomechanical coupling rate in anisotropic materials[12]. The optomechanical
coupling rate includes two terms, one term comes from the change of the
crystal permittivity by pressure and the other term originates from the shifts
of crystal interfaces. The connection is, the calculation of the second term
involves the calculation of eigenfrequency’s derivative with respect to
deformation, which needs Johnson’s or our perturbation method. While previous
isotropic-approximation method causes errors in the second term for
anisotropic situations, our method can almost eliminate the errors. Using our
scheme to calculate the second term, the coupling rate is

$g_{0,MB}=-\frac{\omega^{(0)}}{2}\frac{\int\limits_{S}dA[D_{x},E_{y},E_{z}]^{*}\xi[D_{x},E_{y},E_{z}]^{\prime}h}{\int
E\cdot D^{*}dV},$ (24)

where $h$ is the displacement of the interface in $x$-axis orientation for
mechanical vibration. The coordinate system is defined in Fig.1, and $\xi$ is
defined in Eq. (23). The modification of eigenfrequency derivative based on
our anisotropic theory is typically $\Delta g_{0,MB}=10$ kHz, while the total
coupling rate is several hundred kHz.

Nevertheless, the calculation of perturbed field modes is still elusive as the
expansion of $E^{(1)}$ by $E^{(0)}$ (or similar expansion of magnetic field)
may fail as the shifts of boundaries[1]. This is, in fact, a fundamental
theoretical problem of perturbation theory. Another problem is how to develop
the perturbation theory of Maxwell’s equations dielectric loss is taken into
account. A possible method is to transform the form of Maxwell’s equations to
construct a real and symmetric coefficient tensor. More theoretical work is
still needed in this field.

In conclusion, we report a perturbation theory of Maxwell’s equations for
anisotropic dielectric interfaces. This theory shows a better accuracy
compared with previous isotropic-approximation method. Our method could play
an important role in the future applications of anisotropic materials.

###### Acknowledgements.

This work is partially supported by National Key Research and Development
Program of China (2018YFA0306102, 2018YFA0307400, 2017YFA0304000); National
Natural Science Foundation of China (91836102, 61704164, 12074058, 61775025,
61705033, 61405030, 61308041).

## References

  * [1] Steven G. Johnson, M. Ibanescu, M. A. Skorobogatiy, O. Weisberg, J. D. Joannopoulos, and Y. Fink, Phys. Rev. E 65, 066611 (2002).
  * [2] J. D. Joannopoulos, S. G. Johnson, J. N. Winn, R. D. Meade, Photonic Crystals: Electromagnetism in Mixed Dielectric Media (Princeton University, Princeton, NJ, 1995).
  * [3] D. Marcuse, Theory of Dielectric Optical Waveguides, 2nd ed. (Academic Press, San Diego, 1991).
  * [4] A. W. Snyder and J. D. Love, Optical Waveguide Theory (Chapman and Hall, London, 1983).
  * [5] C. Vassallo, Optical Waveguide Concepts (Elsevier, Amsterdam, 1911).
  * [6] N. R. Hill, Phys. Rev. B 24, 7112 (1981).
  * [7] S. G. Jhonson, M. Ibanescu, M. Skorobogatiy, O. Weisberg, T. D. Engeness, M. Solja$\check{c}$i$\acute{c}$, S. A. Jacobs, J. D. Joannopoulos, and Y. Fink, Opt. Express 9, 748 (2001).
  * [8] M. Skorobogatiy, M. Ibanescu, S. G. Jhonson, O. Weisberg, T. D. Engeness, M. Solja$\check{c}$i$\acute{c}$, S. A. Jacobs, and T. Fink, J. Opt. Soc. Am. B 12, 2867 (2002).
  * [9] C. Cohen-Tannoudji, B. Din, and F. Lalo$\ddot{e}$, Quantum Mechanics (Hermann, Paris, 1977), Vol. 1, Chap. 2; Vol. 2, Chaps. 11 and 13.
  * [10] D. Q. Chowdhury and D. A. Nolan, Opt. Lett. 20, 1973 (1995).
  * [11] H. X. Liang, R. Luo, Y. He, H. W. Jiang and Q. Lin, Optica 4, 1251 (2017).
  * [12] W. T. Jiang, R. N. Patel, F. M. Mayor, T. P. Mckenna, P. Arrangoiz-Arriola, C. J. Sarabalis, J. D. Witmer, P. V. Laer, and A. H. Safavi-Naeini, Optica 6, 845 (2019).