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