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

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.


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][6][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][14][15][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 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 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 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, 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 b of the optical modes propagating in the structure, which is defined by: with k the vacuum wave number, n eff 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 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 ≡ 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 analytical relations shaped as d dim À Á 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 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.  We consider the layer of index n 1 (n 1 > n 3 , n 4 ) 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 n 1 > n eff > n 3 , n 4 . The generic eigenvalue equation ruling all 3-layer slab waveguides is [3,4]: q , 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 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: Then, a temporal derivation operator d dt is applied in order to express dh dt in terms of db 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 n i (i = 1,3,4).
Then, we can compute the derivatives regarding the arguments as quotients of the Arctg function: Both varying terms of the numerator of equation (3) may be written as: From these calculations, the derivative of the height of the guiding layer comes up as: This result can be re-written and formulated as a generic form, with, s 1 = 2prq 7 + 2pr 3 q 5 + 2rp 3 q 5 + 2r 3 p 3 ,q 3 and b defined in equation (1).

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 ≡ 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 n 1 > n eff > n 2 > n 3 , n 4 and (ii) a weakly confined one for which n 1 > n 2 > n eff > n 3 , n 4 . 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]: with q, r and p defined the same way as that of the 3-layer waveguide analysis, and s ¼ . So as to investigate the effect of any change regarding the height of the guiding layer, equation (9) may be re-written: The time derivative of q, r and p are given in equation (4), and ds dt verifies: This allows us to compute the derivatives of the different ratios: as well as the derivatives of the functions Arctg and Arcth: To proceed further with the calculations and so as to lighten the expressions, we define: Derivating this quantity with respect to time yields: Then, we can compute the derivative of the second term of the numerator of equation (11): with: and G ¼ rs 3 q 3 À þ sr 3 q 3 þ ðrqs 5 À þ qs 3 r 3 Þt 2 h . Furthermore, after combining the previous results and simplifying the expressions, the time derivative of the height of the guiding layer is given by: with, (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: Here, so as to find a relation between de dt and db dt the same methodology as before is applied. We thus compute the time-derivative of the q s ratio: Then, we may define a parameter t e before computing its time derivative: Then comes the derivative of the second term of the numerator of equation (18): Arcth Arctg with, D ¼ À2hps 2 q 3 À 2hqs 2 p 3 À s 2 q 3 À qs 2 p 2 þ ð À þ pq 4 À þ q 2 p 3 À ps 2 q 2 À s 2 p 3 Þt e þ À2hps 2 q 3 À 2hqs 2 p 3 À s 2 q 3 À ð qs 2 p 2 Þt 2 e and, Z ¼ pq 3 s 3 þ qp 3 s 3 À þ psq 5 þ sp 3 q 3 ð Þ t 2 e . Finally, by combining the previous different results, the time-derivative of e is: with, s e ¼ 2prq 3 s 7 þ 2qrp 3 s 7 À þ 2pr 3 q 3 s 5 À þ 2qr 3 p 3 s 5 ; x e ¼ À þ 2prs 5 q 5 À þ 2rp 3 q 3 s 5 þ 2pr 3 s 3 q 5 þ 2r 3 p 3 q 3 s 3 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 t h respectively t 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 t 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 j h 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. First, we consider the asymptotic case. Figure 3a represents dn eff dt 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.

Numerical implementation and discussion
Let us now investigate the influence of the rate of variation of the geometrical features. To do so, we plot 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 dh dt ¼ À0:1 mm:s À1 ; À 0:5 mm:s À1 ; À1 mm:s À1 . The other parameters 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 (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 dn eff dt 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 Dn eff (Fig. 4b and  The influence of the set of indices is investigated in Figures 5a and 5b, which represent the evolution of dn 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=0 s) = 0.4 mm, 2e (t=0 s) = 1 mm, de dt ¼ À0:1 mm:s À1 and the investigated set of indices are (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). 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.

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 db dt and 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  17)) for which the height e of the first upper cladding tends to infinity. The other parameters are 2h (t=0s) = 1 mm, dh dt ¼ À0:1 mm: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). 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 1. Compared results of the variation of the eigenvalues Dn 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 optogeometric guiding structure is chosen at fixed parameters: optical indices (n 1 ; n 2 ; n 3 ; n 4 ) = (2; 1.8; 1; 1.5), wavelength l = 0.8 mm and 2e fixed = 0.2 mm, with 2h = 0.4 mm.
TDO (Eq. (17) À Fig. 4a) COMSOL (see Fig. 4b) Variation of Eigenvalues Dn eff 0.100 (directly by the data of coordinate) 0.098 (=1.903 À 1.805) 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 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.