Eur. Phys. J. Appl. Phys. 87, 10501 (2019)
© L. Garnier et al., EDP Sciences, 2019
https://doi.org/10.1051/epjap/2019190092
THE EUROPEAN
PHYSICAL JOURNAL
APPLIED PHYSICS
Regular Article
Temporal derivation operator applied on the historic and school
case of slab waveguides families eigenvalue equations: another
method for computation of variational expressions
Lucas Garnier1, Arthur Doliveira2, Fabrice Mahé2, Etienne Gaviot3, and Bruno Bêche1,4,*
1
2
3
4
Université
Université
Université
Université
France
de Rennes, CNRS, Institut de Physique de Rennes UMR 6251, 35000 Rennes, France
de Rennes, CNRS, IRMAR – UMR 6625, 35000 Rennes, France
du Maine, CNRS, Laboratoire d’acoustique de l’Université du Maine UMR 6613, 72000 Le Mans, France
de Rennes, CNRS, Institut d’Électronique et de Télécommunications de Rennes UMR S 6164, 35000 Rennes,
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
* e-mail: bruno.beche@univ-rennes1.fr
conditions of the electro-magnetic fields. The optogeometrical 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 l, refractive index of each part of
the structure ni and various spatial geometrical dimensions noted dimi, 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
10501-p1
2
L. Garnier et al.: Eur. Phys. J. Appl. Phys. 87, 10501 (2019)
Fig. 1. Generic dispersion curves in integrated photonics; each
curve is associated with a mode (eigenvector) of a fixed structure
at a fixed l-wavelength.
features of the guiding structure (Fig. 1). The correlation
between the effective propagation constant and the optogeometry is highlighted by the above-mentioned eigenvalue
equation (or dispersion relation) noted E.v.eq which mathematically corresponds to a functional F of the different
parameters of the system: E.v.eq = F (l, ni, dimi) here, the
different dimi 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 b of the optical modes
propagating in the structure, which is defined by:
b ¼ k⋅neff ¼
2p
⋅neff ;
l
ð1Þ
with k the vacuum wave number, neff the effective
refractive index associated to the mode and l the vacuum
wavelength.
In this paper, we theoretically investigate the correlation 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 dtd is applied onto the generic
eigenvalue equation regarding the whole family of
structures, considering time-dependent values of the
studied geometrical dimension (dim ≡ dim(t)) with the
effective propagation constant (b ≡ b(t)), and constant
values of l and optical indices. This operation is
tantamount to determining a kind of evolution propagator
on the eigenvalue equations. Such a way to proceed yields
db
Þ
analytical relations shaped as dðdim
dt ¼ f 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
Fig. 2. (a) Schematic representation of the 3-layer slab
waveguides family; the guiding layer yields a refractive index
n1 and a height 2h, the upper cladding n3 and the substrate n4 are
both considered as semi-infinite layers. (b) Schematic representation of the family of 4-layer slab waveguides; the index n1 and n2
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 n3 and n4 are considered
as semi-infinite layers.
essential basic slab-structures: the 3-layer and the 4-layer
waveguides. As the eigenvalue equations are derived for
constant values of b, 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 monitoring; 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
n1 (n1 > n3, n4) with height 2h ≡ 2h (t) as the guiding layer
of such structures. All the effective propagation constants
of guided modes thus verify the light’s cone n1 > neff > n3,
n4. The generic eigenvalue equation ruling all 3-layer slab
waveguides is [3,4]:
p
r
2hðtÞq ¼ Arctg
þ Arctg
þ mp;
ð2Þ
q
q
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
with q ¼ ðk0 n1 Þ2 b2 , r ¼ b2 ðk0 n3 Þ2 , p ¼ b2 ðk0 n4 Þ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
10501-p2
L. Garnier et al.: Eur. Phys. J. Appl. Phys. 87, 10501 (2019)
the projection of the wave vector perpendicular to the plane of
the layers (considering Pythagoras with a zig-zag geometricoptical ray vision).
We consider the temporal variation of h and relate it to
the temporal variation entailed on the effective propagation constant b ≡ b (t). To do so, the above equation is rewritten as:
Arctg pq þ Arctg rq þ mp
:
ð3Þ
h¼
2q
Then, a temporal derivation operator dtd is applied in
db
order to express dh
dt in terms of dt . The first step of this
calculation allows to compute each derivative of parameters p, q and r, considering time dependent values for b and
constant values of k and each ni (i = 1,3,4).
dq
b db dp b db dr b db
¼
;
¼
;
¼
:
dt
q dt dt p dt dt r dt
ð4Þ
Then, we can compute the derivatives regarding the
arguments as quotients of the Arctg function:
d p
q2 þ p2 db d r
q2 þ r2 db
;
b
b : ð5Þ
¼
¼
dt q
dt dt q
dt
pq3
rq3
Both varying terms of the numerator of equation (3)
may be written as:
d
p
q2 þ p2 db
Arctg
b ;
¼ 3
dt
q
pq þ qp3 dt
d
r
q2 þ r2 db
Arctg
ð6Þ
b :
¼ 3
dt
q
rq þ qr3 dt
From these calculations, the derivative of the height of
the guiding layer comes up as:
2
dh
1
q þ p2
q2 þ r2
1
p
¼
þ
þ 3 Arctg
dt
2q pq3 þ qp3 rq3 þ qr3
2q
q
r
db
ð7Þ
þArctg
þ mp b :
q
dt
of the guiding layer and a constant value of the height of the
first upper cladding (h ≡ h (t)) and e = constante): conversely,
we consider the opposite case (h = constante and e ≡ e (t)).
For each one of these configurations, two families of guided
modes exist: (i) a strongly confined mode with effective index
values verifying n1 > neff > n2 > n3, n4 and (ii) a weakly
confined one for which n1 > n2 > neff > n3, n4. Both families
may be studied in parallel thanks to a double notation in the
following equations: the upper notation j ▪ corresponds to the
family of modes verifying condition (i) and the lower one j ▪ 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]:
r
1
0
th Arcth
þ 2es
p
s
Bs
C
þ Arctg@
2hq ¼ Arctg
Aþ mp;
q
q tg Arctg r
2es
s
ð9Þ
with q, r and p defined the sameq
way
as that offfi the 3-layer
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
b2 ðk0 n2 Þ2
waveguide analysis, and s ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi . So as to
ðk0 n2 Þ2 b2
investigate the effect of any change regarding the height of
the guiding layer, equation (9) may be re-written:
r
1
0 s
th Arcth
þ 2es
B q
rs
C
Arctg pq þ Arctg@ s
A þ mp
tg Arctg
2es
q
s
h¼
:
2q
ð10Þ
The time derivative of q, r and p are given in
equation (4), and ds
dt verifies:
þ
ds
¼
dt
b
db
:
s dt
ð11Þ
This allows us to compute the derivatives of the
different ratios:
This result can be re-written and formulated as a
generic form,
p
r
a
þ
j
Arctg
þ
Arctg
þ
mp
1
1
q
q
dh
db
¼
b ; ð8Þ
dt
dt
s1
d p
q2 þ p2 db d r
b ;
¼
dt q
dt dt s
pq3
s2
¼
with,
a1 ¼ rq5 þq3 r3 þ rp2 q3 þ qp2 r3 þ pq5 þ q3 p3 þ pr2 q3 þ qr2 p3
j1 = prq4 + pr3q2 + rpq4 + r3p3,
s 1 = 2prq7 + 2pr3q5 + 2rp3q5 + 2r3p3,q3 and b defined in
equation (1).
3
þ
rs3
r2
s2 þ q2
db d s
db
b ;
b ;
¼
3
dt dt q
dt
sq
ð12Þ
as well as the derivatives of the functions Arctg and Arcth:
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
10501-p3
d
p
q2 þ p2 db
Arctg
b ;
¼ 3
dt
q
pq þ qp3 dt
r 1
0
s2 r2
d @ Arcth s A
db
þ
r ¼
b :
dt
dt
Arctg
rs3 sr3
s
þ
ð13Þ
4
L. Garnier et al.: Eur. Phys. J. Appl. Phys. 87, 10501 (2019)
To proceed further with the calculations and so as to
lighten the expressions, we define:
r
th Arcth
þ 2es
s
ð14Þ
th ¼
r
:
2es
tg Arctg
s
Derivating this quantity with respect to time yields:
dth
dt
2
s
þ
¼
r
2
2
þ
2ers þ 2er
3
þ
rs3
s
þ
þ
r
2ers þ 2er t2h
2
2
2
þ
3
sr3
db
b :
dt
ð15Þ
Then, we can compute the derivative of the second term
of the numerator of equation (11):
d
s
Q db
Arctg t h
¼ b ;
dt
q
G dt
ð16Þ
kh ¼ rq3 s5
þ
þpsr2 q5
q3 s3 r3 þ rqp2 s5
þ
þ2erq3 p3 s3
þ
qp2 s3 r3
þ
ps3 q5
q3 p3 s3 þ sr2 p3 q3 þ 2erps3 q5
þ
þ
2espr3 q5
2esq3 p3 r3 ;
jh ¼ rps3 q4 spr3 q4 þ rp3 s3 q2 sp3 r3 q2
þ
þ
3 3 2
3 5
2 5
þ rpq s
p3 s3 r3 t2h ;
ps r q þ rp s
þ
þ
s h ¼ 2rps3 q7
þ
2spr3 q7 þ 2rp3 s3 q5
þ
2sp3 r3 q5 ;
2ps3 r3 q5 þ 2rq3 p3 s5 2q3 p3 s3 r3 and b deþ
þ
fined in equation (1).
Although equation (11) suits well for the study of the
first case (h ≡ h (t)) and e = constante), it comes up as
inappropriate for the second one. Then, to study the case
e ≡ e (t) and h = constante, it is re-written:
xh ¼ 2rpq5 s5
e¼
Arcth r þArcth
Arctg s
Arctg
q
s tg
2hq
mp
Arctg
2s
p
q
:
ð18Þ
with:
Q ¼ q2 s3
þrs2 q2
þ
þ
sr2 q2
q2 r3 Þt h
þ
þ
2erq2 s3 þ 2esq2 r3 þ ðrs4
ðq2 s3
þ
sr2 q2
þ
þ
s2 r3
2erq2 s3 þ 2erq2 r3 Þt2h ;
sr3 q3 þ ðrqs5 qs3 r3 Þt2h .
þ
þ
Furthermore, after combining the previous results and
simplifying the expressions, the time derivative of the
height of the guiding layer is given by:
db
Here, so as to find a relation between de
dt and dt the same
methodology as before is applied. We thus compute the
time-derivative of the qs ratio:
and G ¼ rs3 q3
dh
dt
¼
ah þ g h th þ kh t 2h þ jh
q2 s2
d q
db
þ
b :
¼
dt s
dt
qs3
ð19Þ
Then, we may define a parameter t e before computing
its time derivative:
te ¼ tg 2hq
Arctg pq þ Arctg sq th þ mp
mp
Arctg
p
:
q
ð20Þ
s h þ xh t 2h
b
db
dt
ð17Þ
with,
ah ¼ rs3 q5
þ q3 p3 s3
¼
þ
sr3 q5 þ rp2 s3 q3
sr2 q3 p3
þ
þ 2esp3 q3 r3 ;
g h ¼ prq3 s4
dt e
dt
þ
þ
sp2 r3 q3 þ ps3 q5
2erps3 q5 þ 2espr3 q5
ps2 r3 q3 þ qrp3 s4 qs2 r3 p3
þ
þ
þ
rs2 p3 q3 q3 p3 r3 ;
þ
þ
spr2 q5
b
2erq3 p3 s3
prs2 q5
þ
2hpq2
2hp3
q2
p2 þ ð2hpq2
pq3 þ qp3
2hp3
db
:
dt
q2
p2 Þt2e
ð21Þ
Then comes the derivative of the second term of the
numerator of equation (18):
pr3 q5
d
dt
10501-p4
D db
Arcth q
t
¼ b
Arctg s e
Z dt
ð22Þ
L. Garnier et al.: Eur. Phys. J. Appl. Phys. 87, 10501 (2019)
with,
2 3
D¼
2hps q
q2 p3
ps2 q2
2 3
2 3
2hqs p
sq
2 2
4
qs p þ ð
pq
þ
þ
2hqs2 p3 s2 q3
s2 p3 Þt e þ ð 2hps2 q3
ðpsq5 þ sp3 q3 Þt2e .
þ
Finally, by combining the previous different results, the
time-derivative of e is:
qs2 p2 Þt2e and, Z ¼ pq3 s3 þ qp3 s3
de
¼
dt
ke t2e
ae þ g e t e þ
Arcth r
j
þ e Arctg s
s e þ xe t2e
Arcth q
t
Arctg s e
db
b
dt
ð23Þ
with,
ae ¼
þ
rq3 s5
þ
qrp2 s5 þ s3 q3 r3 þ qp2 s3 r3
þ pr2 q3 s3 þ qr2 p3 s3
þ
2hprq3 s5
þ
þ
pq3 s5
þ
qp3 s5
2hqrp3 s5 þ2hpq3 r3 s3
þ 2hqp3 r3 s3 ;
ge ¼
þ
prq2 s5
þ p3 r3 s3
ke ¼
þ
þ
þ
rp3 s5
psr3 q4
þ
rq3 s5 þ q3 s3 r3
rps3 q4
rq2 p3 s3 þ pq2 s3 r3
sq2 p3 r3 ;
þ
rqp2 s5 þ qp2 s3 r3 þ ps3 q5
spr2 q5 sr2 p3 q3 2hprq3 s5 2hqrp3 s5
þ
þ
þ
þ
þ 2hpq3 r3 s3 þ 2hqp3 r3 s5 ;
þ q3 p3 s3
jh ¼
þ
rpq3 s4
þ rp3 q3 s2
þ
þ
qrp3 s4 þ pr3 q3 s2 þ qr3 p3 s2 þ ðprs2 q5
pr3 q5
þ
s e ¼ 2prq3 s7 þ 2qrp3 s7
xe ¼
þ
2prs5 q5
þ
r3 p3 rq3 Þt2e ;
þ
2pr3 q3 s5
þ
2qr3 p3 s5 ;
2rp3 q3 s5 þ 2pr3 s3 q5 þ 2r3 p3 q3 s3 and b
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 threspectively tedefined
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 th. 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
5
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 jh conversely, considering the
time varying first upper cladding, the order m of the mode
is present in the parameter t 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 ðdndteff ¼ k1 : db
dt Þ against time for different sets of
opto-geometrical parameters and for different rates of
variation of the geometrical features of the structures ðdh
dt
and de
dt ). For all the following computations, we consider
dh
de
negative values of dt and dt , i.e. decreasing values of h and e
against time, a wavelength l = 0.8 mm and only the first
mode (m = 0) is investigated.
First, we consider the asymptotic case. Figure 3a
represents dndteff against time for a 3-layer slab waveguide and
a 4-layer slab waveguide for which 2e ! ∞. 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
dneff
against time (Fig. 4a) for a 4-layer structure
dt
yielding a variation of the height 2h of the guiding
layer, for three different values of the rate
dh
0:1 mm:s 1 ; 0:5 mm:s 1 ; 1 mm:s 1 . The other padt ¼
rameters are the initial value of the height of the guiding
layer 2h (t=0 s) = 1 mm, the height of the first upper
cladding 2e (t=0 s) = 0.4 mm and the set of indices (n1; n2;
n3; n4) = (2; 1.8; 1; 1.5). The plot of Figure 4 clearly yields
that the ratio between the different values of dndteff is equal
to the ratio between the corresponding 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 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 Dneff
(Fig. 4b and Tab. 1).
10501-p5
6
L. Garnier et al.: Eur. Phys. J. Appl. Phys. 87, 10501 (2019)
n2 layer is more sensitive to the index n3 in the red curve case
because the difference n2 n3 is greater. Concerning the
black curve, as the difference n1 n2 is greater than for the
other two cases, there is less energy in the n2 layer so the
influence of the n3 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
Fig. 3. Asymptotic case. (a) Plot of dndteff 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 mm, dh
0:1 mm:s 1 , and
dt ¼
(n1; n3; n4) = (2; 1.5; 1) for the 3-layer structure and (n1; n2; n3;
n4) =(2; 1.5; 1; 1) for the 4-layer one. (b) Relative error between
the two cases depicted in (a).
The influence of the set of indices is investigated in
Figures 5a and 5b, which represent the evolution of dndteff
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=0 s) = 0.4 mm, 2e (t=0 s) = 1 mm,
de
0:1 mm:s 1 and the investigated set of indices are
dt ¼
(n1; n2; n3; n4) = (2; 1.5; 1; 1) ; (2; 1.8; 1; 1.5); (2; 1.8; 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
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 db
dt and
dðdimÞ
put them up as specific evolution laws of the
dt
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
10501-p6
L. Garnier et al.: Eur. Phys. J. Appl. Phys. 87, 10501 (2019)
7
Fig. 4. (a) Plot of dndteff with time derivation operator (TDO) formulation, against time for the first strongly guided mode of a 4-layer
1
waveguide (Eq. (17)) for three different rates of variation of the height h of the core layer: dh
0:5 mm:s 1 ; 1 mm:s 1 .
dt ¼ 0:1 mm:s ;
The other parameters are 2h (t=0s) = 1 mm, 2efixed = 0.2 mm and (n1; n2; n3; n4) = (2 ; 1.8 ; 1 ;1.5) . For example, by TDO formulation, the
slope at t = 3 s (or 2h = 0.4 mm) is 0.100 ≡ Dneff. (b) The COMSOL simulations associated at the specific point 2h = 0.4 mm (or t = 3 s) of
the previous (a) TDO formulation curves. With COMSOL, the slope is calculated by a simple rate of change Dneff, between two fixed
structures or situations 2h = 0.4 mm (t = 3 s) and 2h = 0.2 mm (t = 4 s): Dneff = 0.098.
10501-p7
8
L. Garnier et al.: Eur. Phys. J. Appl. Phys. 87, 10501 (2019)
Table 1. Compared results of the variation of the eigenvalues Dneff with our variational method based on time
derivation operator (TDO), that is, equation (17) and Figure 4, and COMSOL software. As an example, this optogeometric guiding structure is chosen at fixed parameters: optical indices (n1; n2; n3; n4) = (2; 1.8; 1; 1.5), wavelength
l = 0.8 mm and 2efixed = 0.2 mm, with 2h = 0.4 mm.
TDO (Eq. (17)
Variation of Eigenvalues Dneff
Fig. 4a)
0.100 (directly by the data of coordinate)
COMSOL (see Fig. 4b)
0.098 (=1.903
1.805)
Fig. 5. (a) Plot of dndteff 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 (n1; n2; n3; n4) = (2; 1.5; 1; 1) ; (2; 1.8; 1; 1.5); (2; 1.8; 1.5; 1.5). The other parameters are 2hfixed = 0.4 mm, 2e (t=0s) = 1 mm,
0:1 mm:s 1 : For example, by TDO formulation, the slope at t = 4 s (or 2e = 0.2 mm) is respectively Dneff ≡ 0.007 and 0.048 for
and de
dt ¼
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 Dneff, between two fixed structures or situations
2e = 0.2 mm (t = 4 s) and 2e = 0.1 mm (t = 4.5 s): Dneff = 0.007 and 0.045 respectively for both cones of light.
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,
10501-p8
L. Garnier et al.: Eur. Phys. J. Appl. Phys. 87, 10501 (2019)
9
Table 2. Compared results of the variation of the eigenvalues Dneff 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 (n1; n2 ; n3; n4) = (2 ; 1.8 ; 1; 1.5),
wavelength l = 0.8 mm and 2hfixed = 0.4 mm, with dimension 2e = 0.2 mm.
TDO (Eq. (23)
Variation of Eigenvalues Dneff
Figs. 5a and 5b)
0.007 (directly by the data of coordinate)
0.048 (directly by the data of coordinate)
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 e = (n +
iK)2 of the metal. The value of n2 used in the calculations
has thus to be modified to (n2–K2) 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)
COMSOL (see Figs. 5c and 5d)
0.007 (=1.903
0.045 (=1.640
1.896)
1.595)
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. (SpringerVerlag, 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://creativecom
mons.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)
10501-p9