© L. Garnier et al., EDP Sciences, 2019

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 λ, 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 opto-geometry is highlighted by the above-mentioned eigenvalue equation (or dispersion relation) noted E.v.eq which mathematically corresponds to a functional ℱ of the different parameters of the system: E.v.eq = ℱ (λ, 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}},$$(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 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 $\frac{d}{\mathit{d}t}$ 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 (β ≡ β(t)), and constant values of λ 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\left(dim\right)}{\mathit{d}t}=f\left(\frac{d\beta }{\mathit{d}t}\right)$ 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 β, 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.

Fig. 1 Generic dispersion curves in integrated photonics; each curve is associated with a mode (eigenvector) of a fixed structure at a fixed λ-wavelength.

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.

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 ₁ (n ₁ > n ₃, n ₄) 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 ₁ > n _eff  > n ₃, n ₄. The generic eigenvalue equation ruling all 3-layer slab waveguides is [3,4]:$$2h\left(t\right)q=\text{Arc}\text{t}\text{g}\left(\frac{p}{q}\right)+\text{Arc}\text{t}\text{g}\left(\frac{r}{q}\right)+m\pi \text{,}$$(2) with $q=\sqrt{{\left(k_0{n}_1\right)}^2-{\beta }^2}$, $r=\sqrt{{\beta }^2-{\left(k_0{n}_3\right)}^2}$, $p=\sqrt{{\beta }^2-{\left(k_0{n}_4\right)}^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 β ≡ β (t). To do so, the above equation is re-written as:$$h=\frac{\text{Arc}\text{t}\text{g}\left(\frac{p}{q}\right)+\text{Arc}\text{t}\text{g}\left(\frac{r}{q}\right)+m\pi }{2q}\text{.}$$(3)

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

Then, we can compute the derivatives regarding the arguments as quotients of the Arctg function:$$\frac{d}{\mathit{d}t}\left(\frac{p}{q}\right)=\frac{q^2+p^2}{p{q}^3}\beta \frac{d\beta }{\mathit{d}t},\frac{d}{\mathit{d}t}\left(\frac{r}{q}\right)=\frac{q^2+r^2}{r{q}^3}\beta \frac{d\beta }{\mathit{d}t}\text{.}$$(5)

Both varying terms of the numerator of equation (3) may be written as:$$\frac{d}{\mathit{d}t}\left(\text{A}\text{r}\text{c}\text{t}\text{g}\left(\frac{p}{q}\right)\right)=\frac{q^2+p^2}{p{q}^3+q{p}^3}\beta \frac{d\beta }{\mathit{d}t},\frac{d}{\mathit{d}t}\left(\text{A}\text{r}\text{c}\text{t}\text{g}\left(\frac{r}{q}\right)\right)=\frac{q^2+r^2}{r{q}^3+q{r}^3}\beta \frac{d\beta }{\mathit{d}t}\text{.}$$(6)

From these calculations, the derivative of the height of the guiding layer comes up as:$$\frac{dh}{\mathit{d}t}=\left(\frac{1}{2q}\left(\frac{q^2+p^2}{p{q}^3+q{p}^3}+\frac{q^2+r^2}{r{q}^3+q{r}^3}\right)+\frac{1}{2q^3}\left(\text{Arc}\text{t}\text{g}\left(\frac{p}{q}\right)+\text{Arc}\text{t}\text{g}\left(\frac{r}{q}\right)+m\pi \right)\right)\beta \frac{d\beta }{\mathit{d}t}\text{.}$$(7)

This result can be re-written and formulated as a generic form,$$\frac{dh}{\mathit{d}t}=\frac{{\alpha }_1+{\xi }_1\left(\text{Arc}\text{t}\text{g}\left(\frac{p}{q}\right)+\text{Arc}\text{t}\text{g}\left(\frac{r}{q}\right)+m\pi \right)}{{\sigma }_1}\beta \frac{d\beta }{\mathit{d}t}\text{,}$$(8)with,$${\alpha }_1=r{q}^5+q^3{r}^3+r{p}^2{q}^3+q{p}^2{r}^3+p{q}^5+q^3{p}^3+p{r}^2{q}^3+q{r}^2{p}^3$$

ξ₁ = prq⁴ + pr³q² + rpq⁴ + r³p³,

σ₁ = 2prq⁷ + 2pr³q⁵ + 2rp³q⁵ + 2r³p³,q³ and β 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 ≡ 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 ₁ > n _eff  > n ₂ > n ₃, n ₄ and (ii) a weakly confined one for which n ₁ > n ₂ > n_eff > n ₃, n ₄. Both families may be studied in parallel thanks to a double notation in the following equations: the upper notation |^▪ corresponds to the family of modes verifying condition (i) and the lower one |_▪ 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=\text{Arc}\text{t}\text{g}\left(\frac{p}{q}\right)+\text{Arc}\text{t}\text{g}\left(\frac{s}{q}|\begin{array}{c}\hfill th\left(\text{Arc}\text{t}\text{h}\left(\frac{r}{s}\right)+2es\right)\hfill \\ \hfill tg\left(\text{Arc}\text{t}\text{g}\left(\frac{r}{s}\right)-2es\right)\hfill \end{array}\right)+m\pi \text{,}$$(9)with q, r and p defined the same way as that of the 3-layer waveguide analysis, and $s=|\begin{array}{c}\hfill \sqrt{{\beta }^2-{\left(k_0{n}_2\right)}^2}\hfill \\ \hfill \sqrt{{\left(k_0{n}_2\right)}^2-{\beta }^2}\hfill \end{array}$. So as to investigate the effect of any change regarding the height of the guiding layer, equation (9) may be re-written:$$h=\frac{\text{Arc}\text{t}\text{g}\left(\frac{p}{q}\right)+\text{Arc}\text{t}\text{g}\left(|\begin{array}{c}\hfill \frac{s}{q}th\left(\text{Arc}\text{t}\text{h}\left(\frac{r}{s}\right)+2es\right)\hfill \\ \hfill \frac{s}{q}tg\left(\text{Arc}\text{t}\text{g}\left(\frac{r}{s}\right)-2es\right)\hfill \end{array}\right)+m\pi }{2q}\text{.}$$(10)

The time derivative of q, r and p are given in equation (4), and $\frac{ds}{\mathit{d}t}$ verifies:$$\frac{ds}{\mathit{d}t}=\frac{|\begin{array}{c}\hfill +\hfill \\ \hfill -\hfill \end{array}\beta }{s}\frac{d\beta }{\mathit{d}t}\text{.}$$(11)

This allows us to compute the derivatives of the different ratios:$$\frac{d}{\mathit{d}t}\left(\frac{p}{q}\right)=\frac{q^2+p^2}{p{q}^3}\beta \frac{d\beta }{\mathit{d}t},\frac{d}{\mathit{d}t}\left(\frac{r}{s}\right)=\frac{s^2|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}{r}^2}{r{s}^3}\beta \frac{d\beta }{\mathit{d}t},\frac{d}{\mathit{d}t}\left(\frac{s}{q}\right)=\frac{s^2|\begin{array}{c}\hfill +\hfill \\ \hfill -\hfill \end{array}{q}^2}{s{q}^3}\beta \frac{d\beta }{\mathit{d}t}\text{,}$$(12)as well as the derivatives of the functions Arctg and Arcth:$$\frac{d}{\mathit{d}t}\left(\text{Arc}\text{t}\text{g}\left(\frac{p}{q}\right)\right)=\frac{q^2+p^2}{p{q}^3+q{p}^3}\beta \frac{d\beta }{\mathit{d}t},\frac{d}{\mathit{d}t}\left(|\begin{array}{c}\hfill \text{Arc}\text{t}\text{h}\left(\frac{r}{s}\right)\hfill \\ \hfill \text{Arc}\text{t}\text{g}\left(\frac{r}{s}\right)\hfill \end{array}\right)=\frac{s^2|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}{r}^2}{r{s}^3|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}s{r}^3}\beta \frac{d\beta }{\mathit{d}t}\text{.}$$(13)

To proceed further with the calculations and so as to lighten the expressions, we define:$${\tau }_h=|\begin{array}{c}\hfill th\left(\text{Arc}\text{t}\text{h}\left(\frac{r}{s}\right)+2es\right)\hfill \\ \hfill tg\left(\text{Arc}\text{t}\text{g}\left(\frac{r}{s}\right)-2es\right)\hfill \end{array}\text{.}$$(14)

Derivating this quantity with respect to time yields:$$\frac{d{\tau }_h}{\mathit{d}t}=\frac{s^2|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}{r}^2|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}2er{s}^2+2e{r}^3|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}\left(s^2|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}{r}^2|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}2er{s}^2+2e{r}^3\right){\tau }_h^2}{r{s}^3|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}s{r}^3}\times \beta \frac{d\beta }{\mathit{d}t}\text{.}$$(15)

Then, we can compute the derivative of the second term of the numerator of equation (11):$$\frac{d}{\mathit{d}t}\left(\text{Arc}\text{t}\text{g}\left(\frac{s}{q}{\tau }_h\right)\right)=\frac{\mathrm{\Theta }}{\mathrm{\Gamma }}\beta \frac{d\beta }{\mathit{d}t}\text{,}$$(16)with: $\mathrm{\Theta }=q^2{s}^3|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}s{r}^2{q}^2|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}2er{q}^2{s}^3+2es{q}^2{r}^3+\left(r{s}^4|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}{s}^2{r}^3+r{s}^2{q}^2|\begin{array}{c}\hfill +\hfill \\ \hfill -\hfill \end{array}{q}^2{r}^3\right){\tau }_h|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}\left(q^2{s}^3|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}s{r}^2{q}^2|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}2er{q}^2{s}^3+2er{q}^2{r}^3\right){\tau }_h^2\text{,}$ and $\mathrm{\Gamma }=r{s}^3{q}^3|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}s{r}^3{q}^3+\left(rq{s}^5|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}q{s}^3{r}^3\right){\tau }_h^2$.

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}{\mathit{d}t}=\frac{{\alpha }_h+{\gamma }_h{\tau }_h+{\kappa }_h{\tau }_h^2+{\xi }_h\left(\text{Arc}\text{t}\text{g}\left(\frac{p}{q}\right)+\text{Arc}\text{t}\text{g}\left(\frac{s}{q}{\tau }_h\right)+m\pi \right)}{{\sigma }_h+{\chi }_h{\tau }_h^2}\times \beta \frac{d\beta }{\mathit{d}t}$$(17)with,$${\alpha }_h=r{s}^3{q}^5|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}s{r}^3{q}^5+r{p}^2{s}^3{q}^3|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}s{p}^2{r}^3{q}^3+p{s}^3{q}^5|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}sp{r}^2{q}^5+q^3{p}^3{s}^3|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}s{r}^2{q}^3{p}^3|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}2erp{s}^3{q}^5+2esp{r}^3{q}^5|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}2er{q}^3{p}^3{s}^3+2es{p}^3{q}^3{r}^3,$$ $${\gamma }_h=pr{q}^3{s}^4|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}p{s}^2{r}^3{q}^3+qr{p}^3{s}^4|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}q{s}^2{r}^3{p}^3-pr{s}^2{q}^5|\begin{array}{c}\hfill +\hfill \\ \hfill -\hfill \end{array}p{r}^3{q}^5-r{s}^2{p}^3{q}^3|\begin{array}{c}\hfill +\hfill \\ \hfill -\hfill \end{array}{q}^3{p}^3{r}^3,$$ $${\kappa }_h=r{q}^3{s}^5|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}{q}^3{s}^3{r}^3+rq{p}^2{s}^5|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}q{p}^2{s}^3{r}^3|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}p{s}^3{q}^5+ps{r}^2{q}^5|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}{q}^3{p}^3{s}^3+s{r}^2{p}^3{q}^3+2erp{s}^3{q}^5|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}2esp{r}^3{q}^5+2er{q}^3{p}^3{s}^3|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}2es{q}^3{p}^3{r}^3,$$ $${\xi }_h=rp{s}^3{q}^4|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}sp{r}^3{q}^4+r{p}^3{s}^3{q}^2|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}s{p}^3{r}^3{q}^2+\left(rp{q}^2{s}^5|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}p{s}^3{r}^3{q}^2+r{p}^3{s}^5|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}{p}^3{s}^3{r}^3\right){\tau }_h^2,$$ ${\sigma }_h=2rp{s}^3{q}^7|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}2sp{r}^3{q}^7+2r{p}^3{s}^3{q}^5|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}2s{p}^3{r}^3{q}^5,$ ${\chi }_{\text{h}}=2{\mathrm{rpq}}^5{s}^5|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}2{\mathrm{ps}}^3{r}^3{q}^5+2{\mathrm{rq}}^3{p}^3{s}^5|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}2q^3{p}^3{s}^3{r}^3$ and β defined 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:$$e=\frac{|\begin{array}{c}\hfill -\text{Arc}\text{t}\text{h}\hfill \\ \hfill \text{Arc}\text{t}\text{g}\hfill \end{array}\left(\frac{r}{s}\right)|\begin{array}{c}\hfill +\text{Arc}\text{t}\text{h}\hfill \\ \hfill -\text{Arc}\text{t}\text{g}\hfill \end{array}\left(\frac{q}{s}tg\left(2hq-m\pi -\text{Arc}\text{t}\text{g}\left(\frac{p}{q}\right)\right)\right)}{2s}\text{.}$$(18)

Here, so as to find a relation between $\frac{de}{\mathit{d}t}$ and $\frac{d\beta }{\mathit{d}t}$ the same methodology as before is applied. We thus compute the time-derivative of the $\frac{q}{s}$ ratio:$$\frac{d}{\mathit{d}t}\left(\frac{q}{s}\right)=\frac{|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}{q}^2-s^2}{q{s}^3}\beta \frac{d\beta }{\mathit{d}t}\text{.}$$(19)

Then, we may define a parameter τ _e before computing its time derivative:$${\tau }_e=tg\left(2hq-m\pi -\text{Arc}\text{t}\text{g}\left(\frac{p}{q}\right)\right)\text{.}$$(20) $$\frac{d{\tau }_e}{\mathit{d}t}=\frac{-2hp{q}^2-2h{p}^3-q^2-p^2+\left(2hp{q}^2-2h{p}^3-q^2-p^2\right){\tau }_e^2}{p{q}^3+q{p}^3}\times \beta \frac{d\beta }{\mathit{d}t}\text{.}$$(21)

Then comes the derivative of the second term of the numerator of equation (18):$$\frac{d}{\mathit{d}t}\left(|\begin{array}{c}\hfill \text{Arc}\text{t}\text{h}\hfill \\ \hfill \text{Arc}\text{t}\text{g}\hfill \end{array}\left(\frac{q}{s}{\tau }_e\right)\right)=\frac{\mathrm{\Delta }}{Z}\beta \frac{d\beta }{\mathit{d}t}$$(22)with,

$\mathrm{\Delta }=-2hp{s}^2{q}^3-2hq{s}^2{p}^3-s^2{q}^3-q{s}^2{p}^2+\left(|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}p{q}^4|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}{q}^2{p}^3-p{s}^2{q}^2-s^2{p}^3\right){\tau }_e+\left(-2hp{s}^2{q}^3-2hq{s}^2{p}^3-s^2{q}^3-q{s}^2{p}^2\right){\tau }_e^2\text{and},Z=p{q}^3{s}^3+q{p}^3{s}^3|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}\left(ps{q}^5+s{p}^3{q}^3\right){\tau }_e^2$

Finally, by combining the previous different results, the time-derivative of e is:$$\frac{de}{\mathit{d}t}=\frac{{\alpha }_e+{\gamma }_e{\tau }_e+{\kappa }_e{\tau }_e^2|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}{\xi }_e\left(|\begin{array}{c}\hfill \text{Arc}\text{t}\text{h}\hfill \\ \hfill \text{Arc}\text{t}\text{g}\hfill \end{array}\left(\frac{r}{s}\right)-|\begin{array}{c}\hfill \text{Arc}\text{t}\text{h}\hfill \\ \hfill \text{Arc}\text{t}\text{g}\hfill \end{array}\left(\frac{q}{s}{\tau }_e\right)\right)}{{\sigma }_e+{\chi }_e{\tau }_e^2}\beta \frac{d\beta }{\mathit{d}t}$$(23)with,$${\alpha }_e=|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}r{q}^3{s}^5|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}qr{p}^2{s}^5+s^3{q}^3{r}^3+q{p}^2{s}^3{r}^3|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}p{q}^3{s}^5|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}q{p}^3{s}^5+p{r}^2{q}^3{s}^3+q{r}^2{p}^3{s}^3|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}2hpr{q}^3{s}^5|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}2{\mathrm{hqrp}}^3{s}^5+2{\mathrm{hpq}}^3{r}^3{s}^3+2h{\mathrm{qp}}^3{r}^3{s}^3,{\gamma }_e=|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}pr{q}^2{s}^5|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}r{p}^3{s}^5-{\mathrm{rps}}^3{q}^4-r{q}^2{p}^3{s}^3+p{q}^2{s}^3{r}^3+p^3{r}^3{s}^3|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}\mathrm{ps}{r}^3{q}^4|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}{\mathrm{sq}}^2{p}^3{r}^3\text{,}$$ $${\kappa }_e=|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}r{q}^3{s}^5+q^3{s}^3{r}^3|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}rq{p}^2{s}^5+q{p}^2{s}^3{r}^3+p{s}^3{q}^5+q^3{p}^3{s}^3|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}sp{r}^2{q}^5|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}s{r}^2{p}^3{q}^3|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}2hpr{q}^3{s}^5|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}2hqr{p}^3{s}^5+2hp{q}^3{r}^3{s}^3+2hq{p}^3{r}^3{s}^5,$$ $${\xi }_h=|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}r{p}^{\mathrm{}}{q}^3{s}^4|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}{\mathrm{qrp}}^3{s}^4+p{r}^3{q}^3{s}^2+q{r}^3{p}^3{s}^2+{\mathrm{(prs}}^2{q}^5+r^p{q}^3{s}^2|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}p{r}^3{q}^5|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}{r}^3{p}^{3}r{q}^3){\tau }_e^2,$$ $${\sigma }_e=2pr{q}^3{s}^7+2qr{p}^3{s}^7|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}2p{r}^3{q}^3{s}^5|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}2q{r}^3{p}^3{s}^5,$$ ${\chi }_e=|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}2pr{s}^5{q}^5|\begin{array}{c}\hfill -\hfill \\ \hfill +\hfill \end{array}2r{p}^3{q}^3{s}^5+2p{r}^3{s}^3{q}^5+2r^3{p}^3{q}^3{s}^3$ and β 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 τ _h respectively τ _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 τ _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 ξ _h conversely, considering the time varying first upper cladding, the order m of the mode is present in the parameter τ _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{d{n}_{\text{eff}}}{\mathit{d}t}=\frac{1}{k}.\frac{d\beta }{\mathit{d}t})$ against time for different sets of opto-geometrical parameters and for different rates of variation of the geometrical features of the structures $(\frac{dh}{\mathit{d}t}$ and $\frac{de}{\mathit{d}t}$). For all the following computations, we consider negative values of $\frac{dh}{\mathit{d}t}$ and $\frac{de}{\mathit{d}t}$, i.e. decreasing values of h and e against time, a wavelength λ = 0.8 µm and only the first mode (m = 0) is investigated.

First, we consider the asymptotic case. Figure 3a represents $\frac{d{n}_{\text{eff}}}{\mathit{d}t}$ 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 $\frac{d{n}_{\text{eff}}}{\mathit{d}t}$ 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}{\mathit{d}t}=-0.1\text{µ}m.{\text{s}}^{-1}\text{;}-0.5\text{µ}m.{\text{s}}^{-1};-1\text{µ}m.{\text{s}}^{-1}$. The other parameters are the initial value of the height of the guiding layer 2h _(t=0 s) = 1 µm, the height of the first upper cladding 2e _(t=0 s) = 0.4 µm and the set of indices (n ₁; n ₂; n ₃; n ₄) = (2; 1.8; 1; 1.5). The plot of Figure 4 clearly yields that the ratio between the different values of $\frac{d{n}_{\text{eff}}}{\mathit{d}t}$ is equal to the ratio between the corresponding $\frac{dh}{\mathit{d}t}$ 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}{\mathit{d}t}$. 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 Δn _eff (Fig. 4b and Tab. 1).

The influence of the set of indices is investigated in Figures 5a and 5b, which represent the evolution of $\frac{d{n}_{\text{eff}}}{\mathit{d}t}$ 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 µm, 2e _(t=0 s) = 1 µm, $\frac{de}{\mathit{d}t}=-0.1\text{µ}m.{\text{s}}^{-1}$ and the investigated set of indices are (n ₁; n ₂; n ₃; n ₄) = (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 ₂ layer is more sensitive to the index n ₃ in the red curve case because the difference n ₂ − n ₃ is greater. Concerning the black curve, as the difference n ₁ − n ₂ is greater than for the other two cases, there is less energy in the n ₂ layer so the influence of the n ₃ 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.

Fig. 3 Asymptotic case. (a) Plot of  $\frac{d{n}_{\text{eff}}}{\mathit{d}t}$  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}{\mathit{d}t}=-0.1\text{µ}m.{\text{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).

Fig. 4 (a) Plot of $\frac{d{n}_{\text{eff}}}{\mathit{d}t}$ 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}{\mathit{d}t}=-0.1\text{µ}m.{\text{s}}^{-1}\text{;}-0.5\mu \text{m}.{\text{s}}^{-1};-1\text{µ}m.{\text{s}}^{-1}$. The other parameters are 2h  (t=0s) = 1 µm, 2e fixed = 0.2 µ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 s (or 2h = 0.4 µm) is 0.100 ≡ Δ n eff. (b) The COMSOL simulations associated at the specific point 2h = 0.4 μm (or t = 3 s) of the previous (a) TDO formulation curves. With COMSOL, the slope is calculated by a simple rate of change Δn eff, between two fixed structures or situations 2h = 0.4 μm (t = 3 s) and 2h = 0.2 μm (t = 4 s): Δ n eff = 0.098.

Fig. 5 (a) Plot of $\frac{d{n}_{\text{eff}}}{\mathit{d}t}$ 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 fixed = 0.4 µm, 2e  (t=0s) = 1 µm, and  $\frac{de}{\mathit{d}t}=-0.1\text{µ}m.{\text{s}}^{-1}\text{.}$ For example, by TDO formulation, the slope at t = 4 s (or 2e = 0.2 µm) is respectively Δ n 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 Δn eff, between two fixed structures or situations 2e = 0.2 μm (t = 4 s) and 2e = 0.1 μm (t = 4.5 s): Δ n eff  = 0.007 and 0.045 respectively for both cones of light.

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 }{\mathit{d}t}$ and $\frac{d\left(dim\right)}{\mathit{d}t}$ 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 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, 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 ε = (n + iK)² of the metal. The value of n ² used in the calculations has thus to be modified to (n ²–K ²) 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

All Tables

All Figures

	Fig. 1 Generic dispersion curves in integrated photonics; each curve is associated with a mode (eigenvector) of a fixed structure at a fixed λ-wavelength.
In the text

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.

In the text

Fig. 3 Asymptotic case. (a) Plot of  $\frac{d{n}_{\text{eff}}}{\mathit{d}t}$  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}{\mathit{d}t}=-0.1\text{µ}m.{\text{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).

In the text

Fig. 4 (a) Plot of $\frac{d{n}_{\text{eff}}}{\mathit{d}t}$ 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}{\mathit{d}t}=-0.1\text{µ}m.{\text{s}}^{-1}\text{;}-0.5\mu \text{m}.{\text{s}}^{-1};-1\text{µ}m.{\text{s}}^{-1}$. The other parameters are 2h  (t=0s) = 1 µm, 2e fixed = 0.2 µ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 s (or 2h = 0.4 µm) is 0.100 ≡ Δ n eff. (b) The COMSOL simulations associated at the specific point 2h = 0.4 μm (or t = 3 s) of the previous (a) TDO formulation curves. With COMSOL, the slope is calculated by a simple rate of change Δn eff, between two fixed structures or situations 2h = 0.4 μm (t = 3 s) and 2h = 0.2 μm (t = 4 s): Δ n eff = 0.098.

In the text

Fig. 5 (a) Plot of $\frac{d{n}_{\text{eff}}}{\mathit{d}t}$ 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 fixed = 0.4 µm, 2e  (t=0s) = 1 µm, and  $\frac{de}{\mathit{d}t}=-0.1\text{µ}m.{\text{s}}^{-1}\text{.}$ For example, by TDO formulation, the slope at t = 4 s (or 2e = 0.2 µm) is respectively Δ n 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 Δn eff, between two fixed structures or situations 2e = 0.2 μm (t = 4 s) and 2e = 0.1 μm (t = 4.5 s): Δ n eff  = 0.007 and 0.045 respectively for both cones of light.

In the text