Temporal derivation operator applied on the historic
and school case of slab waveguides families eigenvalue
equations: another method for computation of
variational expressions

Lucas Garnier, Arthur Doliveira, Fabrice Mahé, Etienne Gaviot, Bruno Bêche

► To cite this version:

Lucas Garnier, Arthur Doliveira, Fabrice Mahé, Etienne Gaviot, Bruno Bêche. Temporal derivation operator applied on the historic and school case of slab waveguides families eigenvalue equations: another method for computation of variational expressions. European Physical Journal: Applied Physics, 2019, 87 (1), pp.10501.1-10501-9. 10.1051/epjap/2019190092 . hal-02183526v2

HAL Id: hal-02183526

https://hal.science/hal-02183526v2

Submitted on 2 Oct 2020

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Temporal derivation operator applied on the historic and school case of slab waveguides families eigenvalue equations: another method for computation of variational expressions

Lucas Garnier<fnref n="1" />, Arthur Doliveira<fnref n="2" />, Fabrice Mahé<fnref n="2" />, Etienne Gaviot<fnref n="3" />, and Bruno Bêche<fnref n="1" />,<fnref n="4" />,*

Received: 19 March 2019 / Received in final form: 1 July 2019 / Accepted: 15 July 2019

**Abstract.** Starting from the well-known and historic eigenvalue equations describing the behavior of 3-layer and 4-layer slab waveguides, this paper presents another specific analytical framework providing time-laws of evolution of the effective propagation constant associated to such structures, in case of temporal variation of its various geometrical features. So as to develop such kind of time-propagator formulation and related principles, a temporal derivation operator is applied on the studied school case equations, considering then time varying values of all the geometrical characteristics together with the effective propagation constant. Relevant calculations are performed on three different cases. For example, we first investigate the variation of the height of the guiding layer for the family of 3-layer slab waveguides: then, considering the 4-layer slab waveguide's family, we successively address the variation of its guiding layer and of its first upper cladding. As regards the family of 4-layer waveguides, calculations are performed for two different families of guided modes and light cones. Such another approach yields rigorous new generic analytical relations, easily implementable and highly valuable to obtain and trace all the family of dispersion curves by one single time-integration and one way.

# 1 Introduction and background

Guided optics has known a considerable evolution since the demonstration of the possibility to guide light by refraction in the mid-nineteenth century by Daniel Colladon [1] and John Tyndall [2]. Both a robust theoretical framework [3,4] and a wide range of technological applications have been developed [5–7]. Concerning the theoretical framework, mainly two kinds of families of waveguide geometries are analytically described by an eigenvalue equation: the pure cylindrical symmetry and the planar symmetry waveguide structures made of several different layers. Most of the other geometries need appropriate approximations (such as the pure Galerkin method or equivalent [8,9], more generally based on spectral method and minimization of residue on approximated functions [10], the Marcatili's method and hybrid one [11,12]), or a step-by-step numerical resolution with a segmentation of the space (multilayer matrix method or finite difference time domain method [13–16] applied on conformal transformation [17,18]). The eigenvalue equations regarding the analytically described geometries are derived by solving the Maxwell's equations in all regions of the global structure, and by properly applying the continuity

conditions of the electro-magnetic fields. The opto-geometrical characteristics of the waveguide entail the emergence of a quantification of the modes (called eigenvectors) that can exist within the structure (TE, TM, HE, EH, LP), with indices highlighting the effective number of space-directions being subjected to quantification. For example, the guided modes belong to the light cone of the guiding structure, i.e., in the geometrical optics framework, these modes are coupled with the structure from a region of space allowing total internal reflection inside the structure. Each one of these modes is associated with an effective propagation constant (eigenvalue); the propagation constant is thus directly correlated to the opto-geometrical characteristics of the overall structure (excitation wavelength $\lambda$, refractive index of each part of the structure $n_i$ and various spatial geometrical dimensions noted $dim_i$, $i$ integers, Fig. 1). This dependency between dimension and propagation constant is entirely summarized by the dispersion curve associated with the structure. The dispersion curve is a tool valuable in numerous fields of physics and represents the interdependency between two different physical magnitudes of the studied system. In case of guided optics, this dispersion curve provides a link between the effective propagation constant and the previously described opto-geometrical

<footnote n="1">¹ Université de Rennes, CNRS, Institut de Physique de Rennes – UMR 6251, 35000 Rennes, France</footnote>

<footnote n="2">² Université de Rennes, CNRS, IRMAR – UMR 6625, 35000 Rennes, France</footnote>

<footnote n="3">³ Université du Maine, CNRS, Laboratoire d'acoustique de l'Université du Maine – UMR 6613, 72000 Le Mans, France</footnote>

<footnote n="4">⁴ Université de Rennes, CNRS, Institut d'Électronique et de Télécommunications de Rennes – UMR S 6164, 35000 Rennes,
France</footnote>

<footnote>* e-mail: bruno.beche@univ-rennes1.fr</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Generic dispersion curves in integrated photonics; each curve is associated with a mode (eigenvector) of a fixed structure at a fixed λ -wavelength.</figcaption></figure>

features of the guiding structure (Fig. 1). The correlation between the effective propagation constant and the opto-geometry is highlighted by the above-mentioned eigenvalue equation (or dispersion relation) noted E.v.eq which mathe-matically corresponds to a functional F of the different parameters of the system: E.v.eq = F ( λ , n$_{i}$ , dim$_{i}$ ) here, the different dim$_{i}$ can be considered either as variables or as fixed parameters). Each set of parameters is associated to a unique family of dispersion curves (Fig. 1); any change of one of these parameters yields a change of the shape of the dispersion curves. The resolution of these eigenvalue equations for a given set of opto-geometrical parameters yields the set of related effective propagation constant β of the optical modes propagating in the structure, which is defined by:

$$
\beta = k \cdot n_{\text{eff}} = \frac{2\pi}{\lambda} \cdot n_{\text{eff}}, \qquad (1)
$$

with k the vacuum wave number, n$_{eff}$ the effective refractive index associated to the mode and λ the vacuum wavelength.

In this paper, we theoretically investigate the correla-
tion between the evolution of time-varying geometrical
features of families of slab guiding structures and their
associated effective propagation constant. To this end, a
temporal derivation operator $\frac{d}{dt}$ is applied onto the generic
eigenvalue equation regarding the whole family of
structures, considering time-dependent values of the
studied geometrical dimension ($dim \equiv dim(t)$) with the
effective propagation constant ($\beta \equiv \beta(t)$), and constant
values of $\lambda$ and optical indices. This operation is
tantamount to determining a kind of evolution propagator
on the eigenvalue equations. Such a way to proceed yields
analytical relations shaped as $\frac{d(dim)}{dt} = f(\frac{d\beta}{dt})$ enabling to
represent the global dynamic displacement along all the
family of dispersion curves as geometrical dimensions vary.
The time derivation operator (TDO) formulation is
performed on the eigenvalue equation of both different

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** (a) Schematic representation of the 3-layer slab waveguides family; the guiding layer yields a refractive index $n_1$ and a height $2h$, the upper cladding $n_3$ and the substrate $n_4$ are both considered as semi-infinite layers. (b) Schematic representation of the family of 4-layer slab waveguides; the index $n_1$ and $n_2$ are associated with respectively the guiding layer and the first upper cladding whose height are respectively $2h$ and $2e$. The upper cladding and the substrate of index $n_3$ and $n_4$ are considered as semi-infinite layers.</figcaption></figure>

essential basic slab-structures: the 3-layer and the 4-layer
waveguides. As the eigenvalue equations are derived for
constant values of β, we assume that the time variation of
the geometrical feature is slow enough to perform an
adiabatic evolution treatment. For the 3-layer waveguide
(Fig. 2a), we consider a temporal variation of the height
noted 2h(t) of the guiding layer. As regards the 4-layer
waveguide (Fig. 2b), two cases are considered: firstly, we
consider a temporal variation of the height of the guiding
layer 2h(t) with a constant height of the first upper
cladding 2e(t) before dealing with the opposite situation.
For both of these 4-layer cases, the two families of guided
modes are investigated [19]. Such an approach should be
quite valuable in applied electromagnetism so as to
describe numerous processes: for example, while monitor-
ing; a direct layer deposition as a growth process, a
processes of sedimentation and soft matter creaming, or
etching processes. The results obtained by this model have
been compared to COMSOL numerical simulations and
show a good agreement.

2 Time derivation formulations and analytical
variational expression of eigenvalue equations

2.1 Family of 3-layer slab waveguides

Figure 2a depicts a schematic representation of the family of 3-layer slab waveguides. We consider the layer of index $n_1$ ($n_1 > n_3, n_4$) with height $2h \equiv 2h(t)$ as the guiding layer of such structures. All the effective propagation constants of guided modes thus verify the light’s cone $n_1 > n_{\text{eff}} > n_3, n_4$. The generic eigenvalue equation ruling all 3-layer slab waveguides is [3,4]:

$$
2h(t)q = \operatorname{Arctg}\left(\frac{p}{q}\right) + \operatorname{Arctg}\left(\frac{r}{q}\right) + m\pi, \quad (2)
$$

with $q = \sqrt{(k_0 n_1)^2 - \beta^2}$, $r = \sqrt{\beta^2 - (k_0 n_3)^2}$, $p = \sqrt{\beta^2 - (k_0 n_4)^2}$,
$h$ the half height of the guiding layer, and $m$ the order of the
considered mode. The optical indices act in the expression of
the projection of the wave vector perpendicular to the plane of the layers (considering Pythagoras with a zig-zag geometric-optical ray vision).

We consider the temporal variation of h and relate it to the temporal variation entailed on the effective propagation constant $\beta \equiv \beta(t)$. To do so, the above equation is rewritten as:

$$h = \frac{\mathrm{Arctg}\left(\frac{p}{q}\right) + \mathrm{Arctg}\left(\frac{r}{q}\right) + m\pi}{2q} \quad (3)$$

Then, a temporal derivation operator $\frac{d}{dt}$ is applied in order to express $\frac{dh}{dt}$ in terms of $\frac{d\beta}{dt}$. The first step of this calculation allows to compute each derivative of parameters p, q and r, considering time dependent values for $\beta$ and constant values of k and each $n_i$ ($i=1,3,4$).

$$\frac{dq}{dt} = \frac{-\beta}{q} \frac{d\beta}{dt}, \quad \frac{dp}{dt} = \frac{\beta}{p} \frac{d\beta}{dt}, \quad \frac{dr}{dt} = \frac{\beta}{r} \frac{d\beta}{dt}. \quad (4)$$

Then, we can compute the derivatives regarding the arguments as quotients of the Arctg function:

$$\frac{d}{dt}\left(\frac{p}{q}\right) = \frac{q^2 + p^2}{pq^3} \beta \frac{d\beta}{dt}, \quad \frac{d}{dt}\left(\frac{r}{q}\right) = \frac{q^2 + r^2}{rq^3} \beta \frac{d\beta}{dt}. \quad (5)$$

Both varying terms of the numerator of equation (3) may be written as:

$$\begin{aligned} \frac{d}{dt}\left(\mathrm{Arctg}\left(\frac{p}{q}\right)\right) &= \frac{q^2 + p^2}{pq^3 + qp^3} \beta \frac{d\beta}{dt}, \\ \frac{d}{dt}\left(\mathrm{Arctg}\left(\frac{r}{q}\right)\right) &= \frac{q^2 + r^2}{rq^3 + qr^3} \beta \frac{d\beta}{dt}. \end{aligned} \quad (6)$$

From these calculations, the derivative of the height of the guiding layer comes up as:

$$\frac{dh}{dt} = \left( \frac{1}{2q} \left( \frac{q^2 + p^2}{pq^3 + qp^3} + \frac{q^2 + r^2}{rq^3 + qr^3} \right) + \frac{1}{2q^3} \left( \mathrm{Arctg}\left(\frac{p}{q}\right) + \mathrm{Arctg}\left(\frac{r}{q}\right) + m\pi \right) \right) \beta \frac{d\beta}{dt}. \quad (7)$$

This result can be re-written and formulated as a generic form,

$$\frac{dh}{dt} = \frac{\alpha_1 + \xi_1 \left( \mathrm{Arctg}\left(\frac{p}{q}\right) + \mathrm{Arctg}\left(\frac{r}{q}\right) + m\pi \right)}{\sigma_1} \beta \frac{d\beta}{dt}, \quad (8)$$

with,

$\alpha_1 = rq^5+q^3r^3+rp^2q^3+qp^2r^3+pq^5+q^3p^3+pr^2q^3+qr^2p^3$ \\
$\xi_1 = prq^4+pr^3q^2+rpq^4+r^3p^3$, \\
$\sigma_1 = 2prq^7+2pr^3q^5+2rp^3q^5+2r^3p^3,q^3$ and $\beta$ defined in equation (1).

## 2.2 Family of 4-layer slab waveguides

Considering 4-layer slab waveguides, two cases are dealt with: first, we consider a time dependent value of the height

of the guiding layer and a constant value of the height of the first upper cladding ($h \equiv h(t)$) and $e = \text{constante}$): conversely, we consider the opposite case ($h = \text{constante}$ and $e \equiv e(t)$). For each one of these configurations, two families of guided modes exist: (i) a strongly confined mode with effective index values verifying $n_1 > n_{\text{eff}} > n_2 > n_3$, $n_4$ and (ii) a weakly confined one for which $n_1 > n_2 > n_{\text{eff}} > n_3$, $n_4$. Both families may be studied in parallel thanks to a double notation in the following equations: the upper notation $|\cdot|$ corresponds to the family of modes verifying condition (i) and the lower one $\perp$ is related to condition (ii). The structure of a 4-layer slab waveguide is depicted in Figure 2b. The eigenvalue equation describing the behavior of such a waveguide is [18]:

$$2hq = \mathrm{Arctg}\left(\frac{p}{q}\right) + \mathrm{Arctg}\left(\frac{s}{q}\middle|\frac{th\left(\mathrm{Arcth}\left(\frac{r}{s}\right) + 2es\right)}{tg\left(\mathrm{Arctg}\left(\frac{r}{s}\right) - 2es\right)}\right) + m\pi, \quad (9)$$

with $q, r$ and $p$ defined the same way as that of the 3-layer waveguide analysis, and $s = \left|\frac{\beta^2 - (k_0n_2)^2}{\sqrt{(k_0n_2)^2 - \beta^2}}\right|$. So as to investigate the effect of any change regarding the height of the guiding layer, equation (9) may be re-written:

$$h = \frac{\mathrm{Arctg}\left(\frac{p}{q}\right) + \mathrm{Arctg}\left(\frac{s}{q}\middle|\frac{th\left(\mathrm{Arcth}\left(\frac{r}{s}\right) + 2es\right)}{tg\left(\mathrm{Arctg}\left(\frac{r}{s}\right) - 2es\right)}\right) + m\pi}{2q}. \quad (10)$$

The time derivative of $q, r$ and $p$ are given in equation (4), and $\frac{ds}{dt}$ verifies:

$$\frac{ds}{dt} = \frac{\left| \frac{+}{-} \beta \right|}{s} \frac{d\beta}{dt}. \quad (11)$$

This allows us to compute the derivatives of the different ratios:

$$\begin{aligned} \frac{d}{dt}\left(\frac{p}{q}\right) &= \frac{q^2 + p^2}{pq^3} \beta \frac{d\beta}{dt}, \quad \frac{d}{dt}\left(\frac{r}{s}\right) \\ &= \frac{s^2}{rs^3} \left| \frac{+}{-} \right| \beta \frac{d\beta}{dt}, \quad \frac{d}{dt}\left(\frac{s}{q}\right) = \frac{s^2}{sq^3} \left| \frac{+}{-} \right| \beta \frac{d\beta}{dt}, \end{aligned} \quad (12)$$

as well as the derivatives of the functions Arctg and Arcth:

$$\begin{aligned} \frac{d}{dt}\left(\mathrm{Arctg}\left(\frac{p}{q}\right)\right) &= \frac{q^2 + p^2}{pq^3 + qp^3} \beta \frac{d\beta}{dt}, \\ \frac{d}{dt}\left(\left|\frac{\mathrm{Arcth}\left(\frac{r}{s}\right)}{\mathrm{Arctg}\left(\frac{r}{s}\right)}\right|\right) &= \frac{s^2}{rs^3} \left| \frac{+}{-} \right| \beta \frac{d\beta}{dt}. \end{aligned} \quad (13)$$
To proceed further with the calculations and so as to
lighten the expressions, we define:

$$
\tau_h = \begin{vmatrix} th(\operatorname{Arcth}(\frac{r}{s}) + 2es) \\ tg(\operatorname{Arctg}(\frac{r}{s}) - 2es) \end{vmatrix}. \quad (14)
$$

Derivating this quantity with respect to time yields:

$$
\frac{d\tau_h}{dt} = \frac{s^2 |_{+}^- r^2 |_{+}^- 2ers^2 + 2er^3 |_{+}^- (s^2 |_{+}^- r^2 |_{+}^- 2ers^2 + 2er^3) \tau_h^2}{rs^3 |_{+}^- sr^3} \times \beta \frac{d\beta}{dt}. \quad (15)
$$

Then, we can compute the derivative of the second term
of the numerator of equation (11):

$$
\frac{d}{dt} \left( \operatorname{Arctg} \left( \frac{s}{q} \tau_h \right) \right) = \frac{\Theta}{\Gamma} \beta \frac{d\beta}{dt}, \quad (16)
$$

with:

$$
\begin{align*}
\Theta &= q^2 s^3 |_{+}^- sr^2 q^2 |_{+}^- 2erq^2 s^3 + 2esq^2 r^3 + (rs^4 |_{+}^- s^2 r^3 \\
&\quad + rs^2 q^2 |_{-}^+ q^2 r^3) \tau_h |_{+}^- (q^2 s^3 |_{+}^- sr^2 q^2 |_{+}^- 2erq^2 s^3 + 2erq^2 r^3) \tau_h^2, \\
\text{and } \Gamma &= rs^3 q^3 |_{+}^- sr^3 q^3 + (rqs^5 |_{+}^- q^3 r^3) \tau_h^2.
\end{align*}
$$

Furthermore, after combining the previous results and
simplifying the expressions, the time derivative of the
height of the guiding layer is given by:

$$
\frac{dh}{dt} = \frac{\alpha_h + \gamma_h \tau_h + \kappa_h \tau_h^2 + \xi_h \left( \operatorname{Arctg}\left(\frac{p}{q}\right) + \operatorname{Arctg}\left(\frac{s}{q}\tau_h\right) + m\pi \right)}{\sigma_h + \chi_h \tau_h^2} \times \beta \frac{d\beta}{dt} \quad (17)
$$

with,

$$
\alpha_h = rs^3 q^5 |_{+}^- sr^3 q^5 + rp^2 s^3 q^3 |_{+}^- sp^2 r^3 q^3 + ps^3 q^5 |_{+}^- spr^2 q^5 \\
\qquad + q^3 p^3 s^3 |_{+}^- sr^2 q^3 p^3 |_{+}^- 2erps^3 q^5 + 2espr^3 q^5 |_{+}^- 2erq^3 p^3 s^3 \\
\qquad + 2esp^3 q^3 r^3,
$$

$$
\gamma_h = prq^3 s^4 |_{+}^- ps^2 r^3 q^3 + qrp^3 s^4 |_{+}^- qs^2 r^3 p^3 - prs^2 q^5 |_{-}^+ pr^3 q^5 \\ - rs^2 p^3 q^3 |_{-}^+ q^3 p^3 r^3,
$$

$$
\begin{align*}
\kappa_h &= r q^3 s^5 |_{+}^- q^3 s^3 r^3 + r q p^2 s^5 |_{+}^- q p^2 s^3 r^3 |_{+}^- p s^3 q^5 \\
&\quad + p s r^2 q^5 |_{+}^- q^3 p^3 s^3 + s r^2 p^3 q^3 + 2 e r p s^3 q^5 |_{+}^- 2 e s p r^3 q^5 \\
&\quad + 2 e r q^3 p^3 s^3 |_{+}^- 2 e s q^3 p^3 r^3,
\end{align*}
$$

$$
\xi_h = rps^3 q^4 |_{+}^- spr^3 q^4 + rp^3 s^3 q^2 |_{+}^- sp^3 r^3 q^2 \\
+ (rpq^2 s^5 |_{+}^- ps^3 r^3 q^2 + rp^3 s^5 |_{+}^- p^3 s^3 r^3) \tau_h^2,
$$

$$
\sigma_h = 2rps^3q^7|_{+}^- 2spr^3q^7 + 2rp^3s^3q^5|_{+}^- 2sp^3r^3q^5,
$$

$$
\chi_h = 2rpq^5s^5|_{+}^- 2ps^3r^3q^5 + 2rq^3p^3s^5|_{+}^- 2q^3p^3s^3r^3 \text{ and } \beta \text{ defined in equation (1).}
$$

Although equation (11) suits well for the study of the
first case ($h \equiv h(t)$) and $e = \text{constante}$), it comes up as
inappropriate for the second one. Then, to study the case
$e \equiv e(t)$ and $h = \text{constante}$, it is re-written:

$$
e = \frac{|-\operatorname{Arcth}\left(\frac{r}{s}\right)| + \operatorname{Arcth}\left(\frac{q}{s} \operatorname{tg}\left(2hq - m\pi - \operatorname{Arctg}\left(\frac{p}{q}\right)\right)\right)}{2s}. \quad (18)
$$

Here, so as to find a relation between $\frac{de}{dt}$ and $\frac{d\beta}{dt}$ the same
methodology as before is applied. We thus compute the
time-derivative of the $\frac{q}{s}$ ratio:

$$
\frac{d}{dt}\left(\frac{q}{s}\right) = \frac{|_{+}^{-}q^2 - s^2}{qs^3} \beta \frac{d\beta}{dt}. \quad (19)
$$

Then, we may define a parameter $\tau_e$ before computing
its time derivative:

$$
\tau_e = \operatorname{tg} \left( 2hq - m\pi - \operatorname{Arctg} \left( \frac{p}{q} \right) \right). \quad (20)
$$

$$
\frac{d\tau_e}{dt} = \frac{-2hpq^2 - 2hp^3 - q^2 - p^2 + (2hpq^2 - 2hp^3 - q^2 - p^2)\tau_e^2}{pq^3 + qp^3} \times \beta \frac{d\beta}{dt}. \quad (21)
$$

Then comes the derivative of the second term of the
numerator of equation (18):

$$
\frac{d}{dt}\left(\left|\operatorname{Arctg}\left(\frac{q}{s}\tau_e\right)\right|\right) = \frac{\Delta}{Z}\beta\frac{d\beta}{dt} \quad (22)
$$
with,

$$ \Delta = -2hps^2q^3 - 2hqs^2p^3 - s^2q^3 - qs^2p^2 + (|\underset{+}{-}pq^4|\underset{+}{-}q^2p^3 - ps^2q^2 - s^2p^3)\tau_e + (-2hps^2q^3 - 2hqs^2p^3 - s^2q^3 - qs^2p^2)\tau_e^2 \text{ and, } Z = pq^3s^3 + qp^3s^3|\underset{+}{-}(psq^5 + sp^3q^3)\tau_e^2. $$

Finally, by combining the previous different results, the time-derivative of e is:

$$ \frac{de}{dt} = \frac{\alpha_e + \gamma_e \tau_e + \kappa_e \tau_e^2|\underset{+}{-}\xi_e \left( \left| \operatorname{Arcth} \left(\frac{r}{s}\right) \right| - \left| \operatorname{Arcth} \left(\frac{q}{s}\tau_e\right) \right| \right)}{\sigma_e + \chi_e \tau_e^2} \beta \frac{d\beta}{dt} \quad (23) $$

with,

$$ \alpha_e = |\underset{+}{-}rq^3s^5|\underset{+}{-}qrp^2s^5 + s^3q^3r^3 + qp^2s^3r^3|\underset{+}{-}pq^3s^5|\underset{+}{-}qp^3s^5 + pr^2q^3s^3 + qr^2p^3s^3|\underset{+}{-}2hprq^3s^5|\underset{+}{-}2hqrp^3s^5 + 2hpq^3r^3s^3 + 2hqp^3r^3s^3, $$

$$ \gamma_e = |\underset{+}{-}prq^2s^5|\underset{+}{-}rp^3s^5 - rps^3q^4 - rq^2p^3s^3 + pq^2s^3r^3 + p^3r^3s^3|\underset{+}{-}psr^3q^4|\underset{+}{-}sq^2p^3r^3, $$

$$ \kappa_e = |\underset{+}{-}rq^3s^5 + q^3s^3r^3|\underset{+}{-}rqp^2s^5 + qp^2s^3r^3 + ps^3q^5 + q^3p^3s^3|\underset{+}{-}spr^2q^5|\underset{+}{-}sr^2p^3q^3|\underset{+}{-}2hprq^3s^5|\underset{+}{-}2hqrp^3s^5 + 2hpq^3r^3s^3 + 2hqp^3r^3s^5, $$

$$ \xi_h = |\underset{+}{-}rpq^3s^4|\underset{+}{-}qrp^3s^4 + pr^3q^3s^2 + qr^3p^3s^2 + (prs^2q^5 + r^3p^3q^3s^2|\underset{+}{-}pr^3q^5|\underset{+}{-}r^3p^3rq^3)\tau_e^2, $$

$$ \sigma_e = 2prq^3s^7 + 2qrp^3s^7|\underset{+}{-}2pr^3q^3s^5|\underset{+}{-}2qr^3p^3s^5, $$

$$ \chi_e = |\underset{+}{-}2prs^5q^5|\underset{+}{-}2rp^3q^3s^5 + 2pr^3s^3q^5 + 2r^3p^3q^3s^3 \text{ and } \beta $$

defined in equation (1).

This analysis provides a valuable analytical link between the variation of geometrical features of the global dimensions structures and the correlated changes of all the effective propagation constants describing their associated modes. Each case we have investigated here yields the same generic architecture of relation, yet some significant differences are remarkable. Differences between 3-layer and 4-layer slab waveguides are easily identifiable; they rely on the presence of the factor $\tau_h$ respectively $\tau_e$ defined with equation (16) (respectively (20)) also present in the final relation (17) (respectively (23)). Considering the expressions of the variation of the guiding layer of the family of 4-layer slab waveguides for strongly and weakly guided modes, the only differences concern their sign and the expression of $\tau_h$. Such sign differences are partially due to the definition of the parameter s whose time derivative only yields unlike signs between the strongly and the

weakly guided mode. The other cause of these changed signs is due to both hyperbolic tangent and hyperbolic arctangent functions in the strongly confined mode configuration instead of the conventional tangent and arctangent functions relevant with the weakly confined mode. Naturally, the time derivation of these terms yields only unlike signs. Eventually, the main difference between relations describing time varying guiding layer and time varying first upper cladding is in their dependence as regards the order m of the considered mode. Indeed, in the case of time varying guiding layer this dependency only occurs with the factor noted $\xi_h$ conversely, considering the time varying first upper cladding, the order m of the mode is present in the parameter $\tau_e$. Hence, the influence of the mode order is more significant in the case of time varying first upper cladding.

## 3 Numerical implementation and discussion

Numerical implementations of the previous relations have been performed through Matlab programs. These implementations aim at plotting the temporal derivative of the effective index ($\frac{dn_{eff}}{dt} = \frac{1}{k} \cdot \frac{d\beta}{dt}$) against time for different sets of opto-geometrical parameters and for different rates of variation of the geometrical features of the structures ($\frac{dh}{dt}$ and $\frac{de}{dt}$). For all the following computations, we consider negative values of $\frac{dh}{dt}$ and $\frac{de}{dt}$, i.e. decreasing values of h and e against time, a wavelength $\lambda = 0.8$ µm and only the first mode ($m=0$) is investigated.

First, we consider the asymptotic case. Figure 3a represents $\frac{dn_{eff}}{dt}$ against time for a 3-layer slab waveguide and a 4-layer slab waveguide for which $2e \to \infty$. These two structures yields the same set of parameters (refractive index, initial height of the guiding layer, rate of variation of this height) in order to compare the results in a relevant way. The relative error between these two curves is plotted in Figure 3b. This graph shows that these two cases are numerically equivalent with a maximum relative error below 0.6%. The maximum of relative error occurs at the end of the decrease of the height of the guiding layer, i.e. when 2h is close to the cutting height of the mode.

Let us now investigate the influence of the rate of variation of the geometrical features. To do so, we plot $\frac{dn_{eff}}{dt}$ against time (Fig. 4a) for a 4-layer structure yielding a variation of the height 2h of the guiding layer, for three different values of the rate $\frac{dh}{dt} = -0.1 \ \mu\text{m.s}^{-1}$; $-0.5 \ \mu\text{m.s}^{-1}$; $-1 \ \mu\text{m.s}^{-1}$. The other parameters are the initial value of the height of the guiding layer $2h_{(t=0s)} = 1 \ \mu\text{m}$, the height of the first upper cladding $2e_{(t=0s)} = 0.4 \ \mu\text{m}$ and the set of indices ($n_1; n_2; n_3; n_4$) = (2; 1.8; 1; 1.5). The plot of Figure 4 clearly yields that the ratio between the different values of $\frac{dn_{eff}}{dt}$ is equal to the ratio between the corresponding $\frac{dh}{dt}$ which is easily predictable from equation (17). Concerning the values of the times at which the mode disappear, their ratio is equal to the inverse of the ratio between the corresponding $\frac{dh}{dt}$. In addition, a simulation under COMSOL on a specific structure corroborates this slope, this time considering two close structures to calculate the rate of change $\Delta n_{eff}$ (Fig. 4b and Tab. 1).
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Asymptotic case. (a) Plot of $\frac{dn_{\text{eff}}}{dt}$ in function of time (time derivation operator (TDO) formulation) for two different structures: the black curve stands for a 3-layer structure (Eq. (8)), and the red curve for a 4-layer structure (Eq. (17)) for which the height *e* of the first upper cladding tends to infinity. The other parameters are $2h_{(t=0s)} = 1$ µm, $\frac{dh}{dt} = -0.1$ µm.s$^{-1}$, and $(n_1; n_3; n_4) = (2; 1.5; 1)$ for the 3-layer structure and $(n_1; n_2; n_3; n_4) = (2; 1.5; 1; 1)$ for the 4-layer one. (b) Relative error between the two cases depicted in (a).</figcaption></figure>

The influence of the set of indices is investigated in Figures 5a and 5b, which represent the evolution of $\frac{dn_{\text{eff}}}{dt}$ against time for a 4-layer structure yielding varying value of the height of the first upper cladding, for strongly and weakly guided mode, respectively. The parameters taken for this simulation are $2h_{(t=0s)} = 0.4$ µm, $2e_{(t=0s)} = 1$ µm, $\frac{de}{dt} = -0.1$ µm.s$^{-1}$ and the investigated set of indices are $(n_1; n_2; n_3; n_4) = (2; 1.5; 1; 1)$; $(2; 1.5; 1.5)$; $(2; 1.5; 1.5)$; $(2; 1.5; 1.5)$. We can see from the graph that the more the mode is confined (great difference between the values of indices) the more abrupt is the slope. Another remarkable property is that the red curve begins to increase before the blue one, which begins to increase before the black one. The reason is that the evanescent part (strongly confined mode) or the oscillating part (weakly confined mode) present in the

$n_2$ layer is more sensitive to the index $n_3$ in the red curve case because the difference $n_2 - n_3$ is greater. Concerning the black curve, as the difference $n_1 - n_2$ is greater than for the other two cases, there is less energy in the $n_2$ layer so the influence of the $n_3$ layer is weaker. Similarly, simulations under COMSOL prove the correspondence in terms of variation of the eigenvalue as the dimensions of the structure change (Figs. 5c and 5d and Tab. 2).

The TDO formulation starts from an initial condition, defining a starting point in the eigenvalue equation (that is a starting opto-geometric structure), then it deploys its evolutionary principle by giving access to the variation of the eigenvalue when the dimensions of the structure changes over time.

## 4 Conclusion

This work presents the derivation of generic and rigorous analytical expressions evaluative in time regarding the propagation constant for various families of slab waveguides in case of temporal variation of their geometrical characteristics.

Since the theoretical framework presented in this article hinges on analytical expressions, it provides a valuable tool, easily implementable for numerical computation in one way. Graphically, these expressions are interpretable as a temporal displacement along all the family of dispersion curves of the investigated versatile structures and modes; computing their integral along time obviously leads to the original and generic eigenvalue equations (2) and (9), which gives all the families of dispersion curves (Fig. 1). Their physical meaning, based on the link between $\frac{d\beta}{dt}$ and $\frac{d(dim)}{dt}$ put them up as specific evolution laws of the propagation constant during any growth or deposition process or during an attack or etching process as regards shaping a given structure. Thus, a direct in situ monitoring becomes possible during thin layer deposition processes in ultra-high-vacuum chambers. The analysis has been performed with 3-layer and 4-layer slab waveguides. In the latter case, the investigation has been carried out for both the variation of the guiding layer and of the first upper cladding, and for the two different families of guided modes. The results computed by the proposed TDO method have been compared to numerical calculations obtained via COMSOL simulations and are in good agreement. Dynamics and therefore evolution are intrinsically contained in the TDO formulation, with the possibility of obtaining a whole family of structures at one time. The advantage brought by this TDO model compared to other numerical simulations, besides that it provides analytical expressions that contain the dynamic, resides in the fact that the other numerical simulations can only resolve static structures, step by step. The simulations have thus to be performed at every time step which take a great amount of time, whereas the proposed TDO formulation (Eqs. (8), (17) and (23)) can be implemented on a regular computer via Matlab or Python, for example, and give results within a few seconds. Furthermore, in a theoretical view, an extended and comprehensive study
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** (a) Plot of $\frac{dn_{\text{eff}}}{dt}$ with time derivation operator (TDO) formulation, against time for the first strongly guided mode of a 4-layer waveguide (Eq. (17)) for three different rates of variation of the height $h$ of the core layer: $\frac{dh}{dt} = -0.1 \ \mu\text{m}\cdot\text{s}^{-1}$; $-0.5 \ \mu\text{m}\cdot\text{s}^{-1}$; $-1 \ \mu\text{m}\cdot\text{s}^{-1}$. The other parameters are $2h_{(t=0s)} = 1 \ \mu\text{m}$, $2e_{\text{fixed}} = 0.2 \ \mu\text{m}$ and $(n_1; n_2; n_3; n_4) = (2; 1.8; 1; 1.5)$. For example, by TDO formulation, the slope at $t=3 \ \text{s}$ (or $2h=0.4 \ \mu\text{m}$) is $0.100 \equiv \Delta n_{\text{eff}}$. (b) The COMSOL simulations associated at the specific point $2h=0.4 \ \mu\text{m}$ (or $t=3 \ \text{s}$) of the previous (a) TDO formulation curves. With COMSOL, the slope is calculated by a simple rate of change $\Delta n_{\text{eff}}$, between two fixed structures or situations $2h=0.4 \ \mu\text{m}$ ($t=3 \ \text{s}$) and $2h=0.2 \ \mu\text{m}$ ($t=4 \ \text{s}$): $\Delta n_{\text{eff}} = 0.098$.</figcaption></figure>
<table><caption>**Table 1.** Compared results of the variation of the eigenvalues $\Delta n_{\text{eff}}$ with our variational method based on time derivation operator (TDO), that is, equation (17) and Figure 4, and COMSOL software. As an example, this opto-geometric guiding structure is chosen at fixed parameters: optical indices ($n_1; n_2; n_3; n_4$) = (2; 1.8; 1; 1.5), wavelength $\lambda$ = 0.8 µm and $2e_{\text{fixed}}$ = 0.2 µm, with $2h$ = 0.4 µm.</caption><thead><tr><th></th><th>TDO (Eq. (17) – Fig. 4a)</th><th>COMSOL (see Fig. 4b)</th></tr></thead><tbody><tr><td>Variation of Eigenvalues Δn<sub>eff</sub></td><td>0.100 (directly by the data of coordinate)</td><td>0.098 (=1.903 – 1.805)</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** (a) Plot of $\frac{dn_{\text{eff}}}{dt}$ with time derivation operator (TDO) formulation, against time for (a) the first strongly guided mode and (b) the first weakly guided mode of a 4-layer structure in case of variation of the height $e$ of the first upper cladding (Eq. (23)) and for different sets of indices ($n_1; n_2; n_3; n_4$) = (2; 1.5; 1; 1); (2; 1.8; 1; 1.5); (2; 1.8; 1.5; 1.5). The other parameters are $2h_{\text{fixed}}$ = 0.4 µm, $2e_{(t=0s)}$ = 1 µm, and $\frac{de}{dt}$ = -0.1 µm.s$^{-1}$. For example, by TDO formulation, the slope at $t$ = 4 s (or $2e$ = 0.2 µm) is respectively $\Delta n_{\text{eff}}$ ≈ 0.007 and 0.048 for both cones of light. (c) and (d) Represent the COMSOL simulations associated at the specific point of the previous (a) and (b) TDO formulation curves. With COMSOL, the slope is calculated by a simple rate of change $\Delta n_{\text{eff}}$, between two fixed structures or situations $2e$ = 0.2 µm ($t$ = 4 s) and $2e$ = 0.1 µm ($t$ = 4.5 s): $\Delta n_{\text{eff}}$ = 0.007 and 0.045 respectively for both cones of light.</figcaption></figure>

could be performed in order to determine the overall dynamics of the effective propagation constants in case of the temporal variation of both the above-mentioned layers at the same time. By considering the temporal variation of both $2h(t)$ and $2e(t)$ in the 4-layer family waveguides, the

time-derivation of equations (10) and (18) should lead to a system of two coupled equations whose resolution should be equivalent to the displacement associated with a set of 2D dispersion maps vision instead of a set of 1D dispersion curves. This model is valid not only for dielectric materials,
<table><caption>**Table 2.** Compared results of the variation of the eigenvalues $\Delta n_{\text{eff}}$ with our variational method based on time derivation operator (TDO), that is equation (23) and Figures 5a and 5b, and COMSOL software. As an example, the opto-geometric guiding structure is chosen at fixed parameters: optical indices $(n_1; n_2; n_3; n_4) = (2; 1.8; 1; 1.5)$, wavelength $\lambda = 0.8 \ \mu\text{m}$ and $2h_{\text{fixed}} = 0.4 \ \mu\text{m}$, with dimension $2e = 0.2 \ \mu\text{m}$.</caption><thead><tr><th></th><th>TDO (Eq. (23) – Figs. 5a and 5b)</th><th>COMSOL (see Figs. 5c and 5d)</th></tr></thead><tbody><tr><td rowspan="2">Variation of Eigenvalues Δn<sub>eff</sub></td><td>0.007 (directly by the data of coordinate)</td><td>0.007 (=1.903 – 1.896)</td></tr><tr><td>0.048 (directly by the data of coordinate)</td><td>0.045 (=1.640 – 1.595)</td></tr></tbody></table>

but also in the presence of a metallic layer, as the eigenvalue equations from which our expressions have been derived are valid whatever the material (see Sect. 2.5 of Ref. [3]). In the case of the presence of a metallic layer, when plasmonic waves can establish, one must take into account the imaginary part of the permittivity $\epsilon = (n + iK)^2$ of the metal. The value of $n^2$ used in the calculations has thus to be modified to $(n^2 - K^2)$ in order to describe the physical behavior of a waveguide comprising a metallic layer, including the propagation of plasmonic waves. Extended studies should also be performed in order to investigate the influence of the temporal variation of the opto-geometrical parameters of more complex photonic structures such as rib waveguide, optical fiber, slab waveguide of an arbitrary large number of layers and metal layers.

## Author contribution statement

Bêche and Gaviot developed the idea of this method (called TDO method) and procedure based on the introduction of a time managing the evolution of electromagnetic structures and thus the course of dispersion curves; they have supervised such project and worked on the manuscript. Garnier performed all the analytical calculation plus the associated simulations; he drafted the manuscript. Doliveira and Mahé simulated under COMSOL to corroborate the previous TDO results.

## References

1. J.D. Colladon, C.R. Acad. Sci. **15**, 800 (1842)

2. J. Hecht, *City of Light: The Story of Fiber Optics* (Oxford University Press, Oxford, 1999)

3. M.J. Adams, *Introduction to Optical Waveguides* (John Wiley & Sons, Chichester, 1981)

4. A.W. Snyder, J.D. Love, *Optical Waveguide Theory* (Kluwer Academic Publishers, Dordrecht, 2000)

5. A. Yariv, *Quantum Electronics*, 2nd edn. (John Wiley & Sons, Chichester, 1957)

6. T. Tamir, *Guided-Wave Optoelectronics*, 2nd edn. (Springer-Verlag, Berlin, 1990)

7. D. Marcuse, *Theory of Dielectric Optical Waveguides* (Academic Press, New York, 1974)

8. D. Marcuse, *IEEE J. Quant. Electron.* **28**, 459 (1992)

9. R.L. Gallawa, I.C. Goyal, Y. Tu, K. Gharak, *IEEE J. Quant. Electron.* **27**, 518 (1991)

10. P.N. Robson, P.C. Kendall, *Rib Waveguide Theory by the Spectral Index Method Research*, (Studies Press Ltd., New York, 1990)

11. E.A.J. Marcatili, *Bell Syst. Techn. J.* **48**, 2071 (1969)

12. T. Begou, B. Bêche, N. Grossard, J. Zyss, A. Goullet, G. Jézéquel, E. Gaviot, J. Opt. A: Pure Appl. Opt. **10**, 053310.1 (2008)

13. M. Koshiba, K. Hayata, M. Suzuki, *Electron. Lett.* **18**, 411 (1982)

14. B. Bêche, J.F. Jouin, N. Grossard, E. Gaviot, J. Zyss, *Sens. Actuators Phys. A* **114**, 59 (2004)

15. P. Yeh, A. Yariv, C.-S. Hong, J. Opt. Soc. Am. **67**, 423 (1977)

16. K.H. Schlereth, M. Tacke, *IEEE J. Quant. Electron.* **26**, 627 (1990)

17. M. Heiblum, J. Harris, *IEEE J. Quantum Electron.* **11**, 75 (1975)

18. L. Garnier, C. Saavedra, R. Castro-Beltran, M.J.L. Lucio, E. Gaviot, B. Bêche, *Optik* **142**, 536 (2017)

19. B. Bêche, E. Gaviot, A. Renault, J. Zyss, F. Artzner, *Optik* **121**, 188 (2010)

**Open Access** This article is distributed under the terms of the Creative Commons Attribution License https://creativecommons.org/licenses/by/4.0 which permits unrestricted use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.

**Cite this article as:** L. Garnier, A. Doliveira, F. Mahé, E. Gaviot, B. Bêche, Temporal derivation operator applied on the historic and school case of slab waveguides families eigenvalue equations: another method for computation of variational expressions, Eur. Phys. J. Appl. Phys. **87**, 10501 (2019)