Optimized method of lines for non-linear waveguide problems

The Method of Lines is a semi-analytical versatile tool for the solution of partial differential equations. For the analysis of spatial complex linear waveguide structures, this method is combined with impedance/admittance and field transformation, as well as with finite differences. This paper extends this approach to the treatment of structures with non-linear dielectric materials. The non-linear generalized transmission line equations are derived. An iterative algorithm based on the impedance/admittance transformation with the field transformation obtains efficient and self-consistent solutions. For demonstration, a non-linear stripe waveguide is considered. The Kerr non-linearity was investigated, though the general case is treatable. As a criterion for the correctness of the algorithm, a second harmonic generation as well as a bidirectional, spatially, and temporally energy exchange between the harmonics was examined. The specific limiting factors for the algorithm were explored. The approach can be used for any spatial structure, including, for instance, photonic crystal waveguides and metamaterials.


Introduction
The Method of Lines (MoL) in combination with impedance/admittance and field transformation (IAFT) is used to analyze electromagnetic waves, Helfert and Pregla (2002), Pregla and Helfert (2002), Helfert (2003), Pregla (2008). The use cases are linear waveguiding structures of microwave technology and optics. The core of the theory is the solution of generalized transmission line equations (GTL). The GTL equations describe the relationship between the transverse components of the electric and the magnetic field. The GTL equations are derived from Maxwell's equations with specific boundary conditions for a specific structure. In the case of complex structures, a combination with finite differences (FD) can be used, Helfert and Pregla (1996). The conventional procedure is as follows: The calculation area is covered with lines. In the case of the FD, the model of a structure is divided into sufficiently short sections in the direction of the analytical solution. These sections can be assumed to be homogeneous. Their length is set equal to the step length of the FD. A corresponding interpolation between the points of the FD is important for the efficiency and adequacy of the solution. The discretized GTL equations in matrix form are solved analytically for each homogeneous section (or for each step of the FD). The GTL matrices consist of differential operator matrices with corresponding boundary conditions and material parameters. The materials have so far been assumed to be linear. To determine the field distributions, the solution for the entire structure is usually performed in two procedures: First, an impedance/admittance transformation takes place, and then the field transformation.
The impedance/admittance transformation serves to determine all two-port network parameters of all sections and thus also their respective loads. Because the structure consists of different homogeneous sections, the tangential field components have to be matched at interfaces. The impedance/admittance transformation is one such matching over the sections. The two-port network parameters of the sections can be calculated from the conditions for open circuit and short circuit set at the output of the structure. The calculation of the impedance/admittance transformation is done step by step, i.e., section by section, mostly in the direction from the output of the structure to the input. The field transformation serves to determine the electric and magnetic fields on each homogeneous section or (which is identical) in each step of the FD. The procedure is possible if the impedances and/or admittances of all sections (or all FD steps) have already been calculated during the previous run of the impedance/admittance transformation.
The field transformation is also performed recursively, section by section, mostly from the input to the output of the structure, through all sections. The start value of the field at the input of the structure is specified. It can be, i.e., a Floquet fundamental mode transformed into the original domain. The introductory representation of the GTL follows below. A comprehensive treatment of the MoL for the analysis of electromagnetic waves, i.e., various discretization schemes and boundary conditions, linear GTL equations, impedance/admittance, and field transformation, their combination with the FD, periodic problems, and the application for concrete structures is contained in the monograph (Pregla 2008) as well as e.g., in Pregla (2006a or b).
Generally speaking, all materials are non-linear. When dealing with stronger electromagnetic fields, the non-linearity can no longer be neglected. More so, the non-linearity is seen as a new tool to be used, Gómez-Ullate et al. (2010), Agrawal (2007), Boyd and Fischer (2001), Boyd (2020), Knyazyan et al. (2004). A suitable example is the second harmonic generation (SHG). The MoL (or MoL-IAFT-FD) has its characteristic properties, which are particularly interesting for the analysis of non-linear waveguide phenomena: lower computation time compared to the fully numerical methods, stationary behavior, no phenomenon of relative convergence, very accurate calculation of the field distribution because of the relation to the Discrete Fourier Transformation, no spurious modes, etc., Pregla (2008), Jamid and Akram (2002). The following questions arise: How should the non-linear problem be described? Which solutions are possible? What are the limiting factors of the solution? How can the limiting factors be overcome? Is further research in this direction worthwhile?
The first possible description is non-linear generalized transmission line equations (NGTL). These will be derived from Maxwell's equations, taking into account the polarization of the medium. The NGTL equations are the coupled differential equations of the 1-st order. However, here's to do with non-linear processes: The result of the calculation-e.g., a spatial electric field distribution for each FD step-depends in principle on "itself." According to Hermann and Saravi (2016), it is hypothesized that the equations have one and non-singular solutions. This is finally confirmed by means of an iterative algorithm with a self-consistent, convergent solution. The length of the non-linear part of the structure (interaction length) turned out to be the first limiting factor of the numerical solution.
Here, the required accuracy of the calculation and the numerical effort should be weighed for each specific application. It has been found that a suitable way of interpolating the FD steps can significantly overcome this limit: Instead of the linear and quadratic interpolation already established in Pregla (2006a and b), the use of the one-step methods of the 4th order enabled the analysis of a structure up to 25 times longer.
The paper is organized as follows: Sect. 2 introduces the relevant background in the case of linear materials. The NGTL equations are derived in Sect. 3. Here a scheme is followed that was presented earlier by Knyazyan et al. (2004). Section 4 represents the aspects of the discretization and the SHG. The numerical solution mechanism and the specific numerical problems are presented in Sect. 5. Verification of numerical results and recommendations for users can be found in Sect. 6. The conclusions are presented in Sect. 7. The appendix contains an auxiliary derivation.

Method of lines
The MoL-IAFT-FD is integrated into the individual steps of the impedance/admittance and field transformation. Various complex structures, e.g., microwave technology and optics, can be modeled with it, e.g., fiber gratings (Pregla 2004), photonic crystals (Barcz et al. 2002), effects of heat propagation (Conradi et al. 2001), and many others Pregla 2008). The procedure is as follows: The calculation region is covered with lines. The image of a structure is divided into homogeneous sections in the direction of the analytical solution. The numerical analysis generally consists of two parts: -Solving wave equations in homogeneous sections -Field matching at ports between the homogeneous sections The wave equations are mostly derived with the help of generalized transmission line equations (GTL) (Pregla 1999(Pregla , 2002. Then the GTL is discretized. Various discretization schemes are available for an efficient analysis Pregla 2008). They can be set up in different coordinates, equidistant and non-equidistant, in 2D or 3D Pregla 2008). For most practical cases, the discretization of the coordinates perpendicular to the direction of propagation is assumed to be favorable. In the case of 3D, for example, the cross-section of a structure is discretized, and the analytical solution is used in the direction of propagation. All details of the discretization and boundary conditions are shown in Pregla (2008).

Generalized transmission line equations
The GTL describes the relationship between the transverse components of the electric and the magnetic field. The starting point for deriving the GTL is Maxwell's equations, taking into account the boundary conditions of concrete structures. The detailed representation of the general linear GTL, corresponding wave equations, the tensor of the material parameters and the normalization of the linear masses, as well as the magnetic field strength is shown in Pregla (2006aPregla ( , b, 2008. The following relevant aspects are briefly repeated. The H-and E field components are discretized on two shifted discretization line systems. The shift between two adjacent lines, H and E, is one discretization distance across the direction of propagation. The discretized field components are collected in vectors Ê and Ĥ ( F = Ê t ,Ĥ t t ), which corresponds to the spatial distribution of the complex amplitudes of the respective cross-section. The electric and magnetic fields are calculated on two adjacent lines. For details, see in Pregla (2008, p. 15). The GTL has the general form: The GTL equation applies to any homogeneous section, e.g., for each step of the FD (Sect. 2.5). The complex matrix Q contains differential operators with corresponding boundary conditions and diagonal matrices of the discretized material parameters. Only lossless isotropic material parameters are considered in the paper. If an isotropy is assumed ( S E,H = 0 ), then the form of the resulting GTL equations resembles the form of the well-known telegrapher's equations from the transmission line theory (Chen 2004) ch.V. This analogy is characteristic of the GTL equations, and it determines the nature of the underlying calculation processes. In the case of diagonal material tensors, with no losses and lossless boundaries (i.e., Neumann or Dirichlet), the matrices R E,H are real, symmetrical, and in practical applications, indefinite. This indefiniteness only rarely appears to disturb the result because of the stable nature of the impedance/admittance and field transformation, Spiller (2022). The generalized coordinate u = x, y, z is normalized with the free space wave number k 0 = √ 0 0 according to u = k 0 u . In addition, the magnetic field components H u are normalized with the wave impedance 0 = √ 0 ∕ 0 according to H u = 0 H u . All further aspects are shown in detail (Pregla 2008). The GTL solution for the whole multi-sectioned device structure is carried out in two procedures: The impedances or admittances are transformed from the output to the input. With the input impedance/admittance, the transverse electric and magnetic field components at the input are obtained. Then these quantities are transformed back to the output. In contrast to the transmission line theory, more than just one mode is considered (the number of modes that are considered is equal to the number of lines). Each of these procedures is carried out for the entire structure in the FD steps. The field transformation serves to determine the fields in each homogeneous section or (which is identical) in each step of the FD. The concept of impedance and admittance matrix transformation is used to analyze such multi-sectioned structures. These matrices are transformed from one side of a section to the other one and from one side of an interface between two sections to the other one. In many cases it is advantageous to perform the impedance/admittance transformation and the field transformation with the help of the matrix parameters V E , Z H , Y E and V H , see below. (1)

Impedance/admittance transformation
In the case of a structure composed of different homogeneous sections, the tangential field components must be matched to the transitions. The impedance/admittance transformation is one such matching over the homogeneous sections. The basic idea is the follows. The two-port network parameters of the sections with subports A and B can be calculated from the conditions for open circuit and short circuit set at the output of the structure. This calculation is done step by step, section by section, in the direction from the output of the structure to the input. As a final result, all two-port network parameters of all sections and thus also their respective loads are known.
Here is an example of the recursive calculation of the impedance and admittance transformation formulas, respectively: Z A,B and Y A,B are the impedances and admittances, z ij and y ij are the two-port network impedances and admittances, respectively. For the transformation in the opposite direction, one obtains: An alternative calculation is also used. The following expressions for a recursive calculation of the input impedance and admittance can easily be derived from the expressions given in Pregla (2008), 5.3: The above expressions apply to a step n of an impedance/admittance transformation. For example, the input impedance of the following section ( n + 1 ), the discretized matrix Z B , is related to the input impedance/admittance of the current section Z A . The discretized matrix Z 0 is the characteristic impedance of the current homogeneous section. These expressions are similar to the formulas given in Pregla (2008) for the mode domain but apply in this representation to the original domain.

Field transformation
During field transformation, the electric and magnetic fields are determined using the known impedance and admittances. The field transformation runs section by section in the direction from the input to the output of the structure. The start value of the field at the input of the structure is specified. For example, it can be a Floquet fundamental mode transformed into the original domain.
A recursive field transformation can be calculated using the transfer matrices V in the original domain. The calculation of the field can be done both "forward," from the input to the output of the homogeneous section, and "backward," from the output to the input. The side facing the input of the structure is denoted as (subport) A, and the side facing the output of the structure is (subport) B. The corresponding transmission matrices V AB and V BA are set up for each individual section A-B. The field transformation occurs according to The use of the four parameters for the field transformation is already specified in Pregla (2008): The matrix V will serve as a container for the new interpolation methods in the paper. Its four components can easily be extracted and used for the impedance/admittance transformation.
The content of the transmission matrices consists of the material parameters and the differential operators with the corresponding boundary conditions. The transmission matrices depend on length of the respective section. All details are listed in Pregla (2006aPregla ( , b, 2008.

Impedance/admittance transformation and field transformation with finite differences
In the case of structures with a complex distribution of the material parameters, e.g., photonic crystals, the impedance/admittance transformation, and field transformation can be combined with finite differences (Helfert and Pregla 1996;Pregla 2008). The structure is divided into short sections with subports A and B. For sufficiently small distances Δu = between the subports, the method of finite differences with a corresponding interpolation of the differences can be used. In the past, the two methods of interpolation were built into the GTL solution, the linear and the square (Pregla 2006a(Pregla , b, 2008 In Spiller (2022), the possibility of expanding the usability of the MoL-IAFT-FD is taken up. The approach extends the two previously built-in solution methods in principle to all onestep and multistep methods of various orders of accuracy known in numerical mathematics (Bronstein et al. 2005;Zeidler 2004;Samarski 1986). The methods are built into the transmission matrices V BA and/or V AB for each section of the structure, or in other words, for each step of the FD.

Basic formulas for the nonlinear case
This section introduces non-linearity. It is started with time-dependent Maxwell's curl equations, Saleh and Teich (2007), which are given by The polarization dipole moment per unit volume, or polarization → P(t) of a material system can be described as a power series in the field strength → e (t) , according to, e.g., Boyd (2020): → a e and → a e being the unit vector in the direction of electric field → e . Here is known as the susceptibility and 0 is the permittivity of free space. Further on in the paper, a number as a subscript of a quantity denotes the number of the harmonics of the fundamental frequency , i.e., the susceptibility or the complex amplitude of the electric field for the frequency 2 , 2 and E 2 , respectively. Quantities without a subscript or with subscript 1 correspond to the fundamental frequency, i.e., .
So, following (Boyd 2020), the polarization vector → P can be divided into two parts, a linear and a non-linear, In the paper the non-linear part → P NL of second order is assumed.
The propagation of TE modes with the components E y , H x , H z (complex amplitudes) in the 2D-structure of Fig. 1 is investigated. For this purpose, wave equations corresponding to the application are derived. Thus, the general Eq. (12) with consideration of the model conditions in Fig. 1

reduces to
The subscripts x, y, and z identify the corresponding projections on the coordinate axes. After differentiation of the Eq. (16) with respect to the time and substituting Eq. (15b) and P y from Eq. (14a) is obtained where " r " denotes the relative permittivity of the material. Introducing Eq. (15a) into Eq. (17) one obtains the wave equation for the model in Fig. 1

Second harmonic generation
The complex functions of the field components are assumed, (u = x, z): The subscripts "u1" and "u2" identify the corresponding components of the harmonics ("1" denotes the fundamental frequency ); 0 and 0 are the vacuum permittivity and the magnetic permeability of free space. The magnetic field is normalized with the wave impedance in vacuum 0 = √ 0 ∕ 0 , thus, H u = 0 H u . The complex function of the quadratic term can be described by (see Appendix A1) where "*" denotes a complex conjugate. For the first derivative with respect to time, one obtains The complex functions will now be introduced into Eqs. (15)-(16) and are normalized according to z = k 0 z, x = k 0 x with k 0 = ∕v and v = 1∕ √ 0 0 . Instead of the Eqs. (15a) and (16)  The parts with e j t and e j2 t are separated and one obtains Introducing left Eqs. (26) and (27) into right Eqs. (26) and (27), respectively, one obtains the following coupled wave equations or introducing the difference operator D xx = 2 ∕ x 2

SHG with non-linear transmission line equations
An isotropy of the materials was assumed. The H-and E field components are discretized on two shifted discretization line systems, H and E. The shift between two adjacent lines H and E is one discretization distance across the direction of propagation z. The discretized field components are collected in vectors Ê and Ĥ , which corresponds to the spatial distribution of the complex amplitudes of the respective cross-section. The electric and magnetic fields are calculated on two adjacent lines. For details, see in Pregla (2008, p. 15). The analytical solution is performed in the z-direction, whereby the TM and the TE polarization can be analyzed. The following definitions are introduced: In order to avoid the complexity and overlapping of designations, two of the following notations are considered. First, it is further assumed that the magnetic field is normalized as before, see (19), but without the corresponding " ̃ ". Second, the previously announced normalization of the coordinates also remains, e.g., z , but still without an overline, that is, z.
Assuming (26)-(27), the following non-linear transmission line equations can be written: Here P x denotes the matrix differential operator of the 2nd order with respect to normalized z. This contains the boundary conditions on both sides of the structure, Pregla (2008). The matrices R E and R H are the coefficients of the NGTL equations. The material parameters r and r are discretized to ̂ r and ̂ r (further represented as r and r , without " ̂ "). The symbol "diag()", e.g., diag r1 or diag 21 , denotes the diagonal matrices of the discretized material parameters and of the susceptibility, accordingly. The symbol I denotes an adequate unit matrix. The tensor of the material parameters goes into the vectors r1,r2 .
These represent the dielectric permeability in the respective cross-section of the waveguiding structure on the respective homogeneous section. The same arrangement also applies to the vectors 21,22 . The magnetic permeability r is assumed to be equal to I.
Since R E is a function of z and E 1 , the wave equation has to be in terms of Ê Neglecting the dependence of z in E 1 , one obtains Transformation of the wave equation into the mode domain: The transformation matrix T E was determined from the solution of the eigenmode problem, see (41-b). Â and B are square diagonal matrices of corresponding coefficients for forward and backward propagating modes, respectively. 1 and 2 are vectors. When determining from 2 (i.e. when calculating the square root) care must be taken to use the correct sign since is the propagation constant of the forward propagation mode. Field vectors in the transformed domain are marked with an overline, e.g., " E 1 ".
The parameters of the formulas (43)-(44) only apply to a spatial and a temporal point. The formulas show the "nature" of the interaction between the frequency components and 2 and thus a possibly complex course of the field distributions.

Solution of the NGTL equations
A demonstration of the general approach using a simple test structure as an example certainly cannot cover all emergent aspects. Therefore, it seems useful to discuss the most important aspects that the user may encounter in the case of more complex structures and 1 + e − 2 zÂ 2 + X 2 e 1 zB 1 + e 2 zB 2 more complex excitations. The topics to be treated come from the experience of earlier numerical calculations, including, e.g., Spiller (2019Spiller ( , 2022.

MoL and non-linearity
The NGTL equations deal with non-linear processes. In general, the result of a calculation, i.e., an electric field distribution,-depends in principle on "itself", (38-40). An iterative approach to the calculation is, therefore, an option. There is a possibility to use the solutions (43-44), but also the NGTL Eq. (38) or the wave Eqs. (39-40). The iterative approach means that the MoL is a part of an iterative algorithm with self-consistent, convergent solutions. In other words, a MoL-IAFT-FD solution is performed several times in succession for the entire length of the non-linear structure (interaction length of the SHG). Each next iteration uses the results calculated from the previous iteration as input data, in this case the electric fields. However, it should be ensured that the NGTL-as a non-linear differential equation-has an existing and non-singular solution to the initial value problem. According to Hermann and Saravi (2016), the following hypothesis can be considered. It is assumed that Peano's existence theorem and the uniqueness theorem of Picard and Lindelöf are fullfilled, and the Picard iteration converges; see Hermann and Saravi (2016). Above all, this means that left and right parts of a differential equation are continuous and bounded for all points of the argument. However, the formal requirement of continuity cannot be applied directly to the discretized NGTL. The discretization himself, the possible effects of the complex values, and the possible influence of abrupt transitions in the spatial distribution of the material parameters should also be taken into account. However, a rational loosening of the formal requirement for continuity appears as helpful in this case. From practical experience, it seems less likely that instability can occur as a result of sufficiently fine discretization-as opposed to a hypothetical fully "smooth" model. Therefore, the hypothetical assumption of complete "smoothness" of the material parameters seems reasonable. This loosening of the mathematical rigor is hypothetically used as a practiceoriented tool, whereby the core statement of the mathematically strictly formulated context remains physically valid, but the limitations (discontinuity of material parameters) that are less relevant for practice are omitted. It is also known that for technological reasons, there are never perfect interfaces, but a smooth transition between two regions, Pregla (2008) A.1.3. Therefore such an imaginary approximation/assumption of a smooth transition is justified. A convergence of an iterative solution of the NGTL equations to a meaningful field distribution finally serves as a support for the hypothesis set out above.

On the stability and convergence of the MoL-FD
The matrix Q (or the matrices R E and R H ) can no longer be semidefinite depending on the spatial discretization. To assure numerical stability, additional information about the underlying concept should be used: The impedance/admittance transformation. This concept is inherently numerically stable and gives correct results for every length of section Δz . This fact is based on a direct relation to the transmission line theory, Chen (2004) ch.V, which in principle provides exact analytical solutions. Also, in general, the algorithm of the impedance/admittance transformation can be understood as generalized transmission line theory, Pregla (2008). Thus, the GTL (and NGTL) equations, e.g., (38), show the analogy between the MoL and the well-known telegraph equations. Lines used in the MoL can also be represented as transmission lines. If one has only one mode, the recursive calculations of the impedance/admittance from the MoL reduce to to the well-known impedance/ admittance transformation formulas from the transmission line theory. Hence, as shown below, the following applies: As long as individual sections (as steps of the FD) are calculated according to the rules of the transmission line theory, the impedance/admittance transformation remains for Δz → ∞ and Δz → 0 formally exact, Pregla (2008) (e.g., ch. 2). In other words: As long as the two-port parameters of the last section at the output of the structure are determined by short and open circuiting,-one can assume an exact impedance/admittance transformation. Of course, this accuracy relates to the formally calculated mathematical model, which does not necessarily correspond to the original, e.g., because of a too coarse discretization. Here the term "coarse discretization" is understood in the sense of such a low sampling rate that the Nyquist-Shannon sampling theorem is just about fulfilled. This depends on concrete temporal and spatial field distributions and occasionally limits the maximum length of the structure that can be analyzed. For the test structure in the paper, it's approximately 50-60 discretization points per shortest wavelength in the medium (RK4). It is the responsibility of the user to determine whether the discretization is sufficiently dense. If the sampling theorem is disregarded, the results will be incorrect, which, however, has to do with an incorrect application of the model. But even in this case, the impedance/admittance transformation remains formally exact and stable.
If the impedance/admittance transformation is applied iteratively, the additional factors discussed in Sect. 5 should also be taken into account.
Another numerical aspect can also be relevant for complex structures: The mathematical stiffness of the GTL and NGTL equations. In general, differential equations are mathematically "stiff" if they contain some constructs or parameters that cause rapid variations in the solutions. It is generally difficult to integrate the "stiff" equations by ordinary numerical methods. Small errors may rapidly accumulate (Bronstein et al. 2005;Zeidler 2004;Curtiss and Hirschfelder 1952;Hairer and Wanner 1996). With regard to the MoL, the definition of the stiffness in Bronstein et al. (2005) appears to be suitable: An differential equation is stiff if its solutions are made up of different, strongly exponentially decreasing components. In other words, a stiffness occurs when there is a large difference in scale on the same task. It may well happen that certain strongly decreasing components hardly make a contribution to the solution but have a significant influence on the choice of the step size Δz so that the egg flow of the rounding error O(1∕Δz) increases very strongly, Bronstein et al. (2005). In this case, the equation can cause a particularly high computational effort or, in extreme cases, especially with an adaptive choice of the step size, force the user to abort the calculation because of the apparent "standstill" of the calculation. In the sense, the stiff differential equations can certainly pose challenges with regard to the success of the solutions. The stiffness aspect may require a change to another method of integration of the given differential equation, e.g., to an implicit method or to the method of the other order of accuracy.
The aspects of stability and stiffness in MoL-IAFT-FD, especially in the case of waveguiding structures, were discussed and tested in more detail in the paper (Spiller 2022).

Some waveguide specifics
Spatially complex structures can exhibit complex field distributions, which can pose an almost unforeseeable challenge for the solution of the GTL and NGTL, e.g., due to the stiffness. This applies above all to multimode waveguides with abrupt transitions in material parameters and strong permittivity contrasts. As an example, defect waveguides in photonic crystals can be mentioned, which, e.g., in the region of the defect, can show a Gaussian-shaped field in cross-section, but which are accompanied with possible side lobes with relatively abrupt rising and falling sharp peaks.
Another specific aspect occurs relatively rarely but can occasionally falsify the final result of the calculation if handled improperly. This is a change in the composition of the modes caused by certain factors during repeated calculations of the eigenvalue problem. The cause can be an already very small change of certain input parameters, e.g., a slightly different (or differently placed) longitudinal discretization. This can result in an incorrect calculation of the final result because the wrong mode is used-or the previously used mode is no longer capable of propagation under the new conditions. The background of the effect from the point of view of the numerics is the following. The number of eigenvalues computed as a solution to an eigenvalue problem is usually equal to the number of discretization points along the cross-section of the structure. The eigenvalues correspond to the eigenmodes of the waveguide. But only some of them are capable of propagation and are not evanescent (if the real part of the propagation constant is zero). If the user wants to excite one or more modes capable of propagation (such as the excitation in Sect. 6.1), he should select these modes in the matrices appropriately and correctly distinguish them from the other eigenmodes in subsequent calculations. This can be done in two ways, either by the sequence number of the eigenvalues or by the shape of the corresponding field distribution along the cross-section of the structure (the values of this field distribution are contained in the corresponding column of the eigenvector matrix that has the same sequence number). If the calculations of the eigenvalue problem are repeated with slightly changed input parameters (e.g., with a different number of longitudinal discretization points), it can happen that the eigenmode previously assumed to be guiding now has a different sequence number-or even becomes unable to propagate. The consequences of such a "sudden" change, a reallocation, which is usually unexpected for the user-let's call it a "mode jump"-can be fatal for the correctness of the calculation: The software routine recognizes an incorrect mode as the "leading mode"-or no longer finds a suitable mode at all. This phenomenon also occurs in the Floquet domain: A reallocation of the Floquet modes or a loss of their ability to propagate [the criterion of propagation ability is Re G F = 0 , where G F is the Floquet modes phase, Pregla (2008)].
It is reasonable to assume that the effect of the "mode jumps" (among other things) is closely related to the mathematical stiffness of solutions, e.g., with a coarser discretization, abrupt transitions of the material parameters, and their high contrasts.

Modeling of a non-linear stripe waveguide
For demonstration, a simple stripe waveguide (basic model) is considered (Figs. 1 and 2). This test structure was assumed as a finite length waveguide. The infinitely long, structurally identical, but linear waveguides were connected to the input and output. The electric wall and Neumann boundary conditions for the normal component of the electric field were assumed on both sides of the cladding (Pregla 2008). Due to the higher refractive index, the field is mainly confined inside the core, and the spurious reflections from the cladding sides are negligible. The first Floquet mode of the input waveguide was considered as an excitation. This was determined from the following consideration. By using the Floquets theorem, the period structure is transformed into an equivalent transmission line. In this case, the length of the non-linear part (NL) is taken as the period. Thus, the infinitely long input and output waveguides (L) are assumed to be composed of periods of the same length. This makes it possible to calculate the Floquet parameters of the infinite waveguide (L) as a periodic structure. Assume that the fundamental mode is injected with a unit amplitude, then Ẽ A,f = [1, 0, ⋯ , 0] t is the vector of forward propagating Floquet modes. The forward propagating wave in the original domain is E A,f = S EẼA,f , where S E is a Floquet-modal transformation matrix. For more details, see Pregla (2003) or the monograph (Pregla 2008). It is expected that with the material parameters used, the SHG appears as the strongest effect (Knyazyan et al. 2004). Therefore, only the two frequencies are considered for the sake of simplicity.

Introductory considerations
According to the formulas (43)-(44), in the case of a multimode waveguide, a complex nonlinear field superimposition can take place between all spatial modes of the frequency components involved. This superimposition can have different periodicities depending on the respective values of the parameters involved, e.g., phases and amplitudes of the modes involved. For example, a rather favorable or rather unfavorable "phase matching" between the respective modes of the two components ( and 2 ) could significantly influence the resulting field distributions. As a result of this interaction, a bidirectional exchange of energy between the frequency components can take place [see below, also (Boyd 2020)]. There is also the possibility of the emergence of local rapidly increasing and decreasing exponential components of the field profile. From the point of view of the numerical calculations, this would mean a corresponding mathematical stiffness of the solutions and would require an increased number of discretization points. Because of the sequential mechanism of impedance/admittance transformation, where subsequent values are computed from previous ones, accumulation of computational error can occur with larger interaction lengths. This would also require sufficient accuracy of the calculations. In addition, secondary calculated parameters (e.g., propagation constants) can show a certain dependency on the number of discretization points. This-albeit small-dependency can change the periodicities mentioned above and thus also the field distributions in a multimode waveguide. Not all of these effects can be shown within the scope of the paper: Concrete structures can show different combinations of the effects. However, some general recommendations for users are presented in Sect. 6.4.

Numerical results
The hypothetical assumptions regarding the numerical stability and consistency of the iterations appear to have been confirmed in the numerical experiment so far. A barely acceptable immutability of the end result with a further increase in the discretization density was assumed as a criterion for sufficient discretization. The numerical experiments show that an increased number of discretization points in the direction of wave propagation (or the number of FD steps) became necessary, as indicated, for example, by an estimate according to the Nyquist-Shannon sampling theorem. For the assumed model, about 50-60 discretization points per shortest wavelength in the medium ( 2 component) were found to be necessary. However, the resolution in the lateral direction (number of MoL lines per unit length) seemed to have less influence on the final result and only satisfies the general estimates based on the sampling theorem: the selection of 2-3 discretization points per shortest wavelength in the corresponding medium appeared to be sufficient for the model. However, the specific feature of the MoL-IAFT-FD should be noted: The impedance/ admittance transformation is basically stable, even with very large step lengths that are no longer adequate for the application (Spiller 2022). This means that even with a much too coarse discretization (and/or larger FD step length), no instabilities would occur. In this way, the MoL iterations can also converge. However, the calculated field distributions can be noticeably changed: the actually calculated model would no longer be adequate. For this reason, the above criterion of sufficiently dense discretization appears to be particularly useful in practice. With a sufficient discretization density of the test model, a convergence of the solution was given in all cases of the initial conditions so that only a few iterations were necessary until the calculated distribution of the electric field visually showed no more noticeable differences in comparison to the result of the previous iteration. With this accuracy, it was barely necessary for more than 3-5 iterations. Therefore, the iterative algorithm based on , with the discretization points assigned only to them (abscissa axes). Right, b Spatial field distribution in the cross-section at the end of the calculation window the impedance/admittance transformation with the field transformation ensures efficient and self-consistent solutions. Figure 3a, b shows the spatial field distributions of the TE modes with the components Ey, Hx, and Hz. Figure 3a shows the initial time segment of the SHG. With sufficiently strong amplitudes, a secondary harmonic was effectively generated. Figure 3a shows the spatial field distribution in the cross-section at the end of the calculation window. The effect of the SHG can, in this case, serve the purpose of qualitative verification of the results.
In order to explore the further temporal and spatial development of the field distributions-outside the time window in Fig. 3-an extension of the middle nonlinear section of the structure was considered. The question of a maximum possible interaction length for the calculation also arose. Because of the sequential nature of the impedance/admittance transformation, the MoL iterations should be applied to the entire interaction length. It was found that the interaction length can apparently be considered a limiting factor: because of the computer's memory and/or computational time. So there seems to be a practically realizable maximum value N max = N lat N long for each concrete use case and each concrete computer, where N lat and N long are the number of discretization points in lateral and longitudinal directions (in relation to propagation), respectively. The limitation of the interaction length to be calculated should be understood in two ways: First, it is about a practically justifiable numerical effort (e.g., processor time). Second, because of the sequential nature of the impedance/admittance transformation, accumulation of computational errors can occur, such as rounding errors or discretization errors. In the numerical experiment for the test model, no accumulation of the error could be detected, so the computational time became the limit. Considering the reasonable effort, the length of the middle (nonlinear) section of the base model was increased from 22.2 μ m to 222 μ m. The number of discretization points (or FD steps) in the direction of propagation has also been increased by a factor of 10, i.e., to 40,000. Figure 4 shows the corresponding field distributions of the TE modes of the two frequency components, this time for the increased interaction length. , with the discretization points assigned only to them (abscissa axes). Right, (b): view from above, along the "field strength" axis. An exchange of energy between the frequency components in the course of wave propagation is visible. After the decay of the first field strength peak of the 2 component, there follows a characteristic transition to the following increase For the purpose of simplification, an equality of the refractive indices for the two components and 2 was assumed for this example calculation. The effectiveness of the SHG was rapidly diminished with the introduction of a difference in refractive index (more phase mismatch), even by a few percent, which agrees well with (Boyd 2020). Figure 4 also shows a bidirectional exchange of energy between the frequency components in the course of wave propagation. The "period" of the exchange decreases as the total intensity (or the initial field strength of the component) increases. This agrees with the literature, including (Boyd 2020;Dmitriev and Tarasov 1982), and also the results (Knyazyan et al. 2004). Figure 4b also shows an effect that may be of interest in practice: After the onset and decay of a peak in the field strength of the 2 component, there is a characteristic transition to the next increase in field strength. This transition is characterized by a depression curved in the direction of propagation, (Fig. 4b), the upper part. In the numerical experiment for the test model with different input parameters, such a transition was also observed at the fundamental component (not shown here); however, the depression was bent in the opposite direction to the propagation. The spatial form of this transition appears to be characteristic of the effect of energy exchange. The longitudinal profile of the field strengths also corresponds to the literature mentioned above: the deepening is rather pointed, and the maximum is rather flat.
In the paper, the calculation was performed using conventional linear interpolation (Pregla 2008) as well as newly built-in explicit and implicit Runge-Kutta methods of 4th order of accuracy (RK4 and GRK4 accordingly). An interpolation determines the method of calculating the value of a solution at one step of the FD, e.g., using the already calculated value at the previous step (one-step methods). Choosing a specific interpolation method can, e.g., achieve more accuracy in specific cases than the others with the same step size. For more information, see Spiller (2022). The Runge-Kutta methods (especially GRK4) seemed to have the best accuracy with a smaller number of discretization points in all calculated cases, although the computational time was several times higher than when using linear interpolation. It should be noted that the fulfillment of energetic relationships between the frequency components (Manley-Rowe relation) can serve as a criterion for a correct calculation, including sufficient discretization density. Possibly investigating the influence of the methods of interpolation and the required discretization density-especially in the nonlinear case-is the subject of further research.

Recommendations for users
The first limiting factor in practical calculations can be the interaction length-mainly because of the numerical effort. Thus, for each concrete task and computer, there seems to be a practically realizable maximum value N max = N lat N long , where N lat and N long are the numbers of discretization points in the lateral and longitudinal directions (in terms of propagation), respectively. The required accuracy of the calculation and a possible cumulation of the calculation error should be taken into account.
From Sect. 6.3, it follows that an increased number of discretization points in the direction of wave propagation may be necessary. Alternatively, the use of special interpolation methods for finite differences can also be considered (Spiller 2022). These situations can depend on specific waveguiding structures and are difficult to predict. As a criterion for sufficient discretization, a relative invariability of the end result can be assumed with a further increase in the discretization density. The extent to which such a change in the end result is still acceptable is up to the user to decide for each specific application. Because the error sensitivity of the results can hardly be predicted, it seems useful to vary the discretization density exploratively in each concrete case and to observe the change in the end result until this change comes within the range acceptable to the user. Fulfillment of the energetic ratios between the frequency components involved according to the Manley-Rowe relation can serve as a helpful indication of a correct calculation.

Conclusions
As an extension to the conventional linear MoL-IAFT-FD, a treatment of structures with nonlinear dielectric materials is introduced. The non-linear generalized transmission line equations (NGTL) are derived from Maxwell's equations, taking into account the polarization of the medium. For demonstration, a 2D non-linear multimode stripe waveguide was chosen. The Kerr non-linearity was investigated, though the general case is treatable. An iterative algorithm based on the impedance/admittance transformation with the field transformation obtains efficient and self-consistent solutions. This approach was used to calculate wave propagation and second harmonic generation for various parameters. With sufficiently strong amplitudes, a second harmonic was effectively generated, which qualitatively agrees well with the existing literature. Also, as a criterion for the correctness of the algorithm, a bidirectional, spatially, and temporally energy exchange between the harmonics was investigated. The corresponding qualitative character of the field distributions-as well as their functional dependency on certain parameters-fully corresponds to the literature.
So far, the computation time, interaction length, and/or resolution were found to be possible upper limits for practical calculations. It was found that-in order to ensure sufficient accuracy of the solution for the assumed model-an increased discretization density in the longitudinal direction (direction of the solution) compared to the Nyquist-Shannon sampling theorem was found necessary. The necessary discretization density in the lateral direction was found to correspond approximately to the Nyquist-Shannon sampling theorem.
The approach can be used for any spatial structure, including, for instance, photonic crystal waveguides and metamaterials.

A1: Determination of complex function for e 2 y
It is well known that using Euler's formulas, a time harmonic function f = a cos( t + ) with frequency can also be written as f = 1 2 (ae j t + a * e −j t ) with a = ae j . Let us call a c = ae j t the complex function of f.
In this case, e y is a sum of harmonic functions. This sum can be described as where a c is now the complex function of e y . E 1 , E 2 ... are the phasors of the harmonic functions. Now the complex function of e 2 y is to be obtained. Using Eq. (45) one obtains (45) e y = 1 2 (a c + a * c ) a c = E 1 e j t + E 2 e j2 t + E 3 e j3 t + E 4 e j4 t + ...
(46) e 2 y = 1 4 (a 2 c + a * 2 c + 2a c a * c ) = 1 2 ( a 2 c + a * 2 where a c a * c = b c + b * c + C and E q yc is the complex function of e 2 y . As seen E q yc is equal to E q yc = 1 2 a 2 c + b c . For the parts a 2 c one obtains and for the product a c a * c The first row is equal to a constant value which vanishes after differentiation with respect to time. It is evident that the first parts of the other rows correspond to b c and the second part to b * c . Therefore, the complex function E q yc is now given by Neglect all harmonic functions with frequencies larger than 2 one obtains the following expression for the case of second harmonic generation