Constructing a partially transparent computational boundary for UPPE using leaky modes

TL;DR

This paper presents a novel method for creating a partially transparent computational boundary in optical simulations using leaky modes and artificial-index materials, reducing reflections more effectively than traditional PMLs.

Contribution

The authors introduce a new boundary construction technique employing leaky modes and artificial materials, enhancing spectral propagator simulations of optical wave propagation.

Findings

Significantly reduced boundary reflections in optical simulations.

Effective for spectral propagators and TE electromagnetic waves.

Outperforms traditional PML methods in specific scenarios.

Abstract

In this paper we introduce a method for creating a transparent computational boundary for the simulation of unidirectional propagation of optical beams and pulses using leaky modes. The key element of the method is the introduction of an artificial-index material outside a chosen computational domain and utilization of the quasi-normal modes associated with such artificial structure. The method is tested on the free space propagation of TE electromagnetic waves. By choosing the material to have appropriate optical properties one can greatly reduce the reflection at the computational boundary. In contrast to the well-known approach based on a perfectly matched layer, our method is especially well suited for spectral propagators.

Figures

Equations

$$\nabla \times \text{E}$$

$$\nabla \times \text{B}$$

$$\nabla \cdot \text{E}$$

$$\nabla \cdot \text{B}$$

$$\text{E}(\text{r},t)={\int }_{-\infty }^{\infty }d\omega \text{E}(\text{r},\omega )e^{-i\omega t}.$$

$${\text{P}}_L(\text{r},\omega )={\epsilon }_0\chi (\omega )\text{E}(\text{r},\omega ),$$

$$\text{D}(\text{r},\omega )$$

$$n^2(\omega ,\text{x}))=1+\chi (\omega ,\text{x}).$$

$$\text{E}(\text{r},\omega )$$

$$\text{P}(\text{r},\omega )$$

$${\partial }_{zz}e(x,z,\omega )+{\partial }_{xx}e(x,z,\omega )+$$

$$\Downarrow \text{linearization}$$

$${\partial }_{zz}e(x,z,\omega )+{\partial }_{xx}e(x,z,\omega )+$$

$$e(\pm a_{-},z,\omega )$$

$${\partial }_{x}e(\pm a_{-},z,\omega )$$

$$e(x,z,\omega )=D{e}^{i\beta z}{e}^{i\xi x},$$

$$e(x,z,\omega )=e^{i{\beta }_{0}z}\left(B{e}^{i{\xi }_{0}x}+C{e}^{-i{\xi }_{0}x}\right),$$

$$e(x,z,\omega )=A{e}^{i\beta z}{e}^{-i\xi x},$$

$$\beta$$

$${\beta }_0$$

$${\xi }^2={\left(\frac{\omega }{c}\right)}^2(n^2-1)+{\xi }_0^2.$$

$$\left(\begin{array}{cccc}{e}^{ia\xi }& -e^{-ia{\xi }_0}& -e^{ia{\xi }_0}& 0\\ -i{e}^{ia\xi }\xi & -i{e}^{-ia{\xi }_0}{\xi }_0& i{e}^{ia{\xi }_0}{\xi }_0& 0\\ 0& e^{ia{\xi }_0}& e^{-ia{\xi }_0}& -e^{ia\xi }\\ 0& i{e}^{ia{\xi }_0}{\xi }_0& -i{e}^{-ia{\xi }_0}{\xi }_0& -i{e}^{ia\xi }\xi \end{array}\right)\left(\begin{array}{c}A\\ B\\ C\\ D\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right),$$

$$\tan (2a{\xi }_0)+i\frac{2\xi {\xi }_0}{{\xi }^2+{\xi }_0^2}=0.$$

$$\left\{\left\{\xi ,{\xi }_0\right\},\left(A,B,C,D\right)\right\}$$

$$\left(\begin{array}{c}A\\ B\\ C\\ D\end{array}\right)=\left(\begin{array}{c}0\\ -1\\ 1\\ 0\end{array}\right).$$

$$\tan (2a{\xi }_0)$$

$$\tan (2a{\xi }_0)=-i\frac{2{\xi }_0|{\xi }_0|\sqrt{1+\frac{\alpha }{{\xi }_0^2}}}{2{\xi }_0^2\left(1+\frac{\alpha }{2{\xi }_0^2}\right)}=i\sqrt{1+\frac{\alpha }{{\xi }_0^2}}\left(\frac{1}{1+\frac{\alpha }{2{\xi }_0^2}}\right)\approx i\left(1-\frac{{\alpha }^2}{8x^4}\right),$$

$$-i\frac{e^{i(2ax+i2ay)}-e^{-i(2ax+i2ay)}}{e^{i(2ax+i2ay)}+e^{-i(2ax+i2ay)}}\approx i\left(1-\frac{{\alpha }^2}{8x^4}\right),$$

$$-\frac{e^{i(4ax+i4ay)}-1}{e^{i(4ax+i4ay)}+1}\approx \left(1-\frac{{\alpha }^2}{8x^4}\right),$$

$$-\frac{r{e}^{i\theta }-1}{r{e}^{i\theta }+1}\approx \left(1-\frac{{\alpha }^2}{8x^4}\right),$$

$$\frac{1-r^2}{r^2+2r\cos \theta +1}-i\frac{2r\sin \theta }{r^2+2r\cos \theta +1}\approx \left(1-\frac{{\alpha }^2}{8x^4}\right).$$

$$\theta$$

$$\frac{1-r^2}{r^2+2r\cos \theta +1}$$

$$\frac{1-r}{1+r}$$

$$r$$

Full text

1 Introduction

Treatment of domain boundaries in numerical simulations, especially in the solution of partial differential equations, presents a long-standing problem. While powerful methods have been developed for certain situations, they often introduce significant additional complexity and computational overhead. The perfectly matched layer (PML) [1] approach stands as a prime example of methods that work extremely well in situations where a transparent boundary is meant to mimic a connection of the given “computational box” to an infinite outside space. Indeed, PML-based methods are routinely employed for wave-propagation simulation, for example in finite-difference Maxwell solvers [2] and in beam-propagation simulation [3].

Nevertheless, there are applications for which good boundary treatments are still lacking. For example, in extreme nonlinear optics, characterized by high intensity, few cycle pulses, which, through their interaction with material degrees of freedom, display very broad and complex spatio-temporal spectra, spectral beam and pulse propagators [4] are the preferred methods of choice. Unfortunately they do not mesh well with the boundary treatments developed for the finite-difference solvers such as PML.

While spectral propagators applied to pulses and/or beams shine in many situations that are next to impossible to handle with finite-difference approaches, the boundary treatment can be a significant problem. For example, long-distance propagation of highly nonlinear optical pulses [5] is often connected with light-matter interactions that send significant energy propagating toward the boundaries of computational domains where it must be “absorbed” as if propagating into infinite space. In connection to spectral-based numerical simulation of beam and pulse propagation, this is a difficult problem that we aim to address in this work.

The method we put forward can be understood as an extension of an approximation that is sometimes used to simulate beam and pulse propagation in leaky waveguides [6, 7] such as hollow-core fibers or capillaries [8]. In such a context, the propagating modes are approximated [9, 10] by real parts of the true leaky modes for the given waveguide [11], while their propagation constants are redefined by inclusion of the imaginary parts that reflect the propagation loss of a leaky mode. Such an approach can be interpreted as a first-order perturbation theory where eigenvalues are corrected while the wave functions are kept unchanged. Needless to say, this only works when the physics dictates that the propagation is dominated by a relatively small number of modes that have small propagation losses.

We propose to utilize the true leaky modes, without approximations, as the basis for both the numerical representation of the optical field and for the realization of transparent boundary conditions. We introduce an artificial structure outside of the given computational domain in order to introduce an infinite set of quasi-normal modes, and construct an expansion of an arbitrary beam profile. While we present the treatment for a fixed frequency, the generalization to pulsed waveform is straightforward.

Leaky modes has had a long history in the field of electromagnetics. They were used already as early as in 1884 by J. J. Thomson[12] in his study of decay phenomena in electromagnetics. Since then, they have been of enduring interest in electromagnetics, for resonator cavities [13, 14], optical waveguides [15], photonic [16] and plasmonic  [17, 18] structures, and are often used for numerical simulations, which is also what we propose to do in the current paper. Leaky modes are decaying eigenstates and as such have played an important role in quantum theory from its very inception until today. In this setting they describe unstable states. Such states were first defined in terms of the absence of incoming waves by Siegert[19] for the nuclear scattering matrix. Siegert’s definition of unstable states was taken up by Peierls[20],Couteur[21] and Humblet[22] and by them refined into an important tool for nuclear scattering theory. The wave functions satisfying the Siegert outgoing-wave conditions are known as resonant states, and their properties has been of interest for many years [23, 24, 25, 26, 27, 28].

As is evident from the previous paragraph, leaky modes and unstable states have a long history and has been, and is, of great utility [29, 30, 31, 32] in various fields. However, the fact that they are decaying eigenstates means that the corresponding eigenvalue problems are not self adjoint. Consequently, the matter of projecting general field configurations into sum of leaky modes or resonant states, and the question of completeness of the resulting expansions, are not backed up by any general theory, like for the self adjoint case. In fact the leaky modes and resonant states are invariably growing exponentially in space, and thus can not be placed in some well known inner product spaces. The lack of a general theory for non-self adjoint operators is challenging, and it means that questions of projection and completeness have to be handled in a case by case basis. In this paper we will introduce a projection method for leaky modes based on a naturally occurring complex non-Hermitian inner product, but will not present a convergence proof for our leaky mode expansions.

The paper is organized as follows. In section two we introduce the model which we will use to access the feasibility of our proposed approach to setting up a partially reflective boundary for UPPE. The model describes the propagation of TE electromagnetic waves in a homogeneous medium, that we for convenience assume is a vacuum. We then proceed to set up and solve the eigenvalue problem for the complex transverse wave numbers that define the leaky modes. In this section we also derive a very accurate explicit asymptotic formula for the location of leaky mode wave numbers in the complex plane. In section three we introduce the leaky modes and show that, by using the technique, well know from the study of resonant states in quantum theory, of shifting them over to a complex spatial contour outside the transverse computational domain, the leaky modes can be identified with vectors in a vector space of functions on the real line that is endowed with a complex non-Hermitian inner product[33][34]. The leaky modes are orthogonal with respect to this product and we can thus write down generalized Fourier series for any given function based on the orthogonal leaky modes and this non-Hermitian inner product. This solves the projection problem for our leaky modes.

We have done extensive numerical experiments using our leaky mode expansions and in section four we presents some examples, and the conclusions we draw from these examples, with regards to their suitability for representing initial data for UPPE. We argue that the leaky mode expansions converge point wise for all sufficiently smooth functions in our space, but that they do not always converge to the function used to generate the expansion. The point wise convergence only becomes problematic in the limit when the index, of the artificial material introduced outside the computational domain, tends towards the same value as the index inside the domain. One would expect that problems with the leaky mode expansions would appear in this limit, since if the limit is reached, there is no index difference between the inside and the outside of the computational domain and leaky modes cease to exist. However, in order to minimize the reflection from the boundary of the computational domain we want to choose the difference between the inside index and the outside, artificial index, as small as possible. It thus becomes a trade off between making it small in order to minimize reflections, and not making it so small that the leaky mode expansions stops giving a good representation of the functions used to generate the expansions. At the end of section four we argue, using a dimensionless quantity that appear from our theory, that there is an acceptable trade off that can be made.

In this paper we do not present a proof that the leaky mode expansion converge to the function used to generate the series. The chief reason for this is that we believe that they never really do converge point wise to the function used to generate them. This is what our numerical results from section four indicated. In section five we present analytical arguments that points to the same conclusion. However, the conjectured lack of point wise convergence to the desired function does not make the leaky mode expansions useless from a more practical point of view. This what we argued in section four where we used a certain dimensionless quantity to specify what we mean by a practical point of view in this context.

2 The model

In a situation where there are no free charges or currents, Maxwell’s equations in the frequency domain can be written in the form

[TABLE]

Here, we use the sign convention for the inverse Fourier transform with respect to time that is standard in optics,

[TABLE]

The polarization is a sum of a linear and a non-linear part. The linear part, which is the focus of the current paper, takes in frequency domain, for a material that is dispersive and possibly spatially inhomogeneous, the form

[TABLE]

and thus the electric displacement takes the form

[TABLE]

where $n=n(\omega,\textbf{x})$ is the refractive index of the material, defined as usual by the identity

[TABLE]

We will now assume that the spatial inhomogeneity of the refractive index takes the form of a straight channel oriented along the z-axis of our coordinate system, of uniform width $2a$ in the transverse direction, which is oriented along the x-axis of our coordinate system. The geometry of the channel is illustrated in figure 1.

Consistent with the geometry we assume that the electromagnetic field is transverse electric. Thus we have

[TABLE]

Using Maxwell’s equations we find that $e(x,z,\omega)$ is a solution to the following model equation

[TABLE]

In addition to the model equation (10), the electric field $e(x,z,\omega)$ must satisfy the conditions

[TABLE]

which follows from the electromagnetic interface conditions for transverse electric fields at $x=\pm a$.

The goal is now to find leaky modes for the linearized equation. These modes can then be used to write down a UPPE version of the nonlinear equation (10), where a leaky mode transform takes the place of the usual transverse Fourier transform. The rest of the paper is focused on constructing the leaky modes and evaluating for which transverse field configurations they form a suitable basis.

3 Leaky modes

Leaky modes are solutions to the linearized model equation (11) that are propagating in the direction of the the positive z-axis, satisfy the electromagnetic interface conditions (12),(13), and are outgoing at positive and negative infinity.

Such functions must be of the form

[TABLE]

where $\beta_{0},\xi_{0}$ and $\beta,\xi$ are the propagation constants and transverse wave numbers inside and outside the channel, respectively. The propagation constants are determined by the transverse wave numbers by the identities

[TABLE]

From the physical point of view, the modes represents electromagnetic disturbances that propagate in the direction of the positive z-axis while they are partially reflected and transmitted at the lateral boundaries defining the index channel. This is illustrated in figure 2.

In order for the functions (14) to satisfy the electromagnetic boundary conditions (13), and thus be leaky modes, the two propagation constants $\beta$ and $\beta_{0}$ must be equal, which is only true if the following identity holds

[TABLE]

This is Snell’s law. In addition, the following linear algebraic system

[TABLE]

must have a unique solution. This can only happen if the the determinant of the matrix defining the system is zero. One can show that the determinant is zero if and only if the transverse wavenumber satisfy following equation

[TABLE]

Equation (18), together with Snell’s law (16), will determine the dispersion law pertaining to each separate leaky mode.

Note that the system of (16)-(18) has two symmetries connecting solutions. If we denote solutions using the notation $\left\{\left\{\xi,\xi_{0}\right\},\left(A,B,C,D\right)\right\}$, the two symmetries are of the form

[TABLE]

Let us start by observing that using the symmetries (19) and (20), it is enough to find all the solutions

3.1 Dispersion laws

In this section we will design asymptotic formulas for all solutions $\xi,\xi_{0}$ to equations (16),(18), and thus determine all modes for the system and their respective dispersion laws.

Let us start by observing that $\xi_{0}=0$ is a solution to equation (18) and that the corresponding solution vector to the linear system (16) is given by

[TABLE]

However, we observe that if we insert $\xi_{0}=0$ and the vector (21) into the formula for the modes (14) we find that the corresponding mode is identically zero. Thus the solution $\xi_{0}=0$ only gives us a trivial mode which can be disregarded when we use the modes for expanding electric field configurations.

Observe that because of the symmetries (19) and (20), it is sufficient to consider the case when $\xi_{0}$ is in the second quadrant. Any solution in one of the other quadrants can be generated from a solution in the second quadrant by using the symmetries. In the second quadrant, we can split the system (16),(18) into to separate systems depending on which square root we take when equation (16) is used to express $\xi$ as a function of $\xi_{0}$

[TABLE]

where we have defined $\alpha=(\omega/c)^{2}(n^{2}-1)$. In figure (3) we display the solutions of the first of the two systems, (22). In the figure, the solutions are defined by the intersection of the zero contours for the real and imaginary part of the equation for $\xi_{0}$ in (22). There clearly exists an infinite set of solutions, each one corresponding to a distinct mode with its associated dispersion law. A similar plot for the second of the two systems, (23), gives convincing numerical evidence that it has no solutions in the second quadrant and thus this system does not give us any additional modes in the second quadrant.

It is evident that for most solutions displayed in figure (3), the real part strongly dominate the imaginary part. This fact can be used to find an asymptotic formula for the solutions to equation (18).

Assuming that $|\xi_{0}|\gg\sqrt{\alpha}$, equation (18) can be approximated by

[TABLE]

because $|\xi_{0}|\equiv\sqrt{\xi_{0}^{2}}$, is equal to $-\xi_{0}$ when $\xi_{0}$ is in the second quadrant.

Judging from the locations of the zeros in figure (3), most of them will be found in regions of the complex plane where $|\text{Re}[\xi_{0}]|\gg|\text{Im}[\xi_{0}]|$. We therefore write $\xi_{0}=x+iy$, where $|x|\gg|y|$, and use this to simplify the first equation from (24) as follows

[TABLE]

where $r=\text{exp}(-4ay)$ and $\theta=4ax$. We find the real and imaginary part of the left-hand side of equation (25) to be

[TABLE]

The right-hand side in equation (26) is real, so we must have

[TABLE]

where $k$ is a whole number. We know from figure (3), that for this choice of $\xi$ the determinant has solutions in the $2^{nd}$ quadrant, so $k$ is a negative whole number. Also, we assume that $y>0$. According to equation (27), $\theta$ is a whole number multiple of $\pi$. Let $k$ be first an odd number $k=2p+1$. Then the imaginary part of the left-hand side in equation (26) we get

[TABLE]

Notice that this result is $>1$, while the real part of the right-hand side in equation (26) is $<1$. This is a contradiction which tells is that $k$ cannot be odd. On the other hand, with $k$ being even, $k=2p$, we get

[TABLE]

where both sides are less than one. Solving equation (29) for $r$ and consequently for $y$ we get

[TABLE]

This gives us the approximative solutions to equation (18) in the $2^{nd}$ quadrant.

[TABLE]

where $p=1,2,\cdots$. In order to get an even better approximation, we can write eq. (29) as

[TABLE]

where the fourth power on the right-hand side allows us to write $\xi_{0}$ instead of $x$ , because we work with asymptotic expression under the assumption that $|x|\gg|y|$. This leads to the following iteration scheme for the solutions $\xi_{0p}$

[TABLE]

where $n$ is the index in the recursive formula with $\xi_{0p}^{0}=-p\pi/(2a)$. As it turns out, we get a very good approximation already for $n=2$ and thus an asymptotic approximation to the solutions of (16),(18) in the 2st quadrant is

[TABLE]

Using the symmetries (19),(20) we get the formulas for modes residing in the other three quadrants in the form

[TABLE]

Writing the formula defining $\xi$ in terms of $\xi_{0}$ from (22) in the form $\xi=\xi^{\prime}+i\xi^{\prime\prime}$, it is evident that if $\xi_{0}$ is in the second quadrant, then $\xi^{\prime}>0$ and $\xi^{\prime\prime}<0$. Using the formula (14) and our convention for the inverse Fourier transform (5) we can conclude that the modes in the second quadrant, determined by formula (34), are outgoing and exponentially growing in the transverse direction. From the formula for the propagation constants (15) it is also evident that they are decaying in the propagation direction. These are thus leaky modes. In a similar way the modes determined by formula (37) are also outgoing and decaying in the propagation direction, and thus are also leaky modes. We find however that the modes determined by formulas (35) and (36) are incoming and growing in the propagation direction. These modes are thus not leaky modes, but gaining modes.

Even if we assumed $p$ being large, the formulas for all four quadrants give surprisingly good results, even when $p$ is of order 1. However, it is exactly in this region where the formulas can break down. Observe that the inner logarithm in formulas (35-37) must be positive in order to stay in the same quadrant. Therefore, these formulas become invalid if

[TABLE]

Let us assume that $\sqrt{\alpha}\gg|\xi_{0p}|$. Applying this assumption to the equation for $\xi_{0}$ from (22), which determine the leaky modes in the 2st quadrant, gives us

[TABLE]

where $\xi_{0}=x+iy$ and $r=\text{exp}(-4ay)$, $\theta=4ax$. Using the same approach as before, we find that

[TABLE]

Under the assumed condition $\sqrt{\alpha}\gg|\xi_{0p}|$, the second term in equation (40) is a small correction to the first term. This allows us to look at equation (40) as a recursion formula for the solution $\xi_{0p}$. Starting from the leading term $\xi_{0p}^{0}=-\pi p/(2a)$, we obtain from the first iteration

[TABLE]

We have thus obtained a different asymptotic formula for the solutions $\xi_{0p}$ in the 2st quadrant, a formula where we know that the imaginary part is a small correction to the real part, $\pi p/(2a)\gg\pi p/(2a^{2}\sqrt{\alpha})$, or equivalently, $p\ll a\sqrt{\alpha}$. This condition implies that condition (38) holds. Thus, we can conclude that the asymptotic formula (41) holds exactly when the asymptotic formula (35) breaks down. Formulas similar to (41) can be derived for the other quadrants.

Some of the leaky modes are paraxial whereas others are not. In order to be more precise about which modes are paraxial, note that the propagation vector for the light beam is of the form $(\xi_{0},\beta_{0})$. This is clear from equation (14). This allows us to calculate the propagation angle of the beam with respect to the $z$-axis. This angle is

[TABLE]

Clearly, for each $\xi_{0p}$, we get a different angle. In order for a mode to be paraxial, the angle $\theta_{p}$ must be small, and this holds only if

[TABLE]

where we have used the fact to leading order $\xi_{0p}\approx-\frac{\pi p}{2a}$. Formula (43) determine which leaky modes are paraxial.

We now investigate if there are zeros in parts of the complex $\xi_{0}$-plane that are not covered by the asymptotic formulas we have found so far. We will focus on the second quadrant, the other quadrants can be treated in a similar way with corresponding results. These investigations are necessary, because the exponential smallness of the equation determining $\xi_{0}$, in the part of the second quadrant well away from the real axis, makes a direct numerical search for solutions, like the one in figure 3, very challenging.

Let us first look for zeros in the part of the second quadrant where $y=\text{Im}[\xi_{0}]$ is much larger than $x=\text{Re}[\xi_{0}]$ and $\text{Im}[\xi_{0}]\gg 1\gg\alpha$. Under these conditions on $\xi_{0}$ we have $\tan(2a\xi_{0})\approx i$ and equation for $\xi_{0}$ in (22) takes the simplified form

[TABLE]

where we included the next term of the Taylor series for the right-hand side in equation (22). We thus end up with a solution $\xi_{0}=i\sqrt{\alpha}$ that contradicting the assumptions imposed on $\xi_{0}$. Hence, no zeros can exist in this part of the second quadrant.

Let us next look at the region where $y=\text{Im}[\xi_{0}]\sim\text{Re}[\xi_{0}]=x$ and $y,x\gg 1\gg\alpha$. At this point, observe that the left hand side of equation (18) is not actually equal to the determinant of the system(17). Imposing equation (18) only implies that the determinant is zero. If we rather equates the full determinant of (17) to zero we get

[TABLE]

The above equation can be simplified using the assumptions $|\xi_{0}|\gg\alpha$. Writing $\xi_{0}=-x+iy$, where $x,y>0$ we get

[TABLE]

The only possible solutions to equation (46) are $x=y(\sqrt{2}-1),x=y(\sqrt{2}+1)$. Substituting these solutions into (47), and solving for $y$, we find in both cases t $y=0$. This contradict our assumptions and thus there are no solutions in this region of the second quadrant either. We have now covered all possible regions of the second quadrant and thus conclude that there are no other zeros of the determinant, and thus leaky modes, than the ones we have already found and that is covered by our asymptotic formulas.

There is however one remaining issue related to the zeros, and thus leaky modes, that needs to be discussed. As we have already noted, the asymptotic formulas for the zeros, which, by design, are expected to be accurate only in the limit when the index $p$ is very large, in fact works surprisingly well even for $p$ as small as 2. However, the very first zero, the one corresponding to $p=1$, is never very accurate. The first zero also behave differently when the parameter $\alpha$ is varied. Recall that the value of this parameter is proportional to the size of the index step defining the channel where the waves will be propagating. We are interested in minimizing reflections from the edges of the channel and therefore would want to make the index step, and hence the parameter $\alpha$, as small as possible. When we let alpha decrease, we observe that all the zeros in the second quadrant, except the first, move slowly up, and even more slowly towards the imaginary axis. This behavior is to be expected from of the logarithmic dependence of the imaginary part of the zero on the parameter $\alpha$. The first zero approaches the imaginary axis at at fast rate when $\alpha$ is decreased, and for a finite value of $\alpha=\alpha_{c}$, it simply vanishes. For $\alpha<\alpha_{c}$ we observe that for one value of the index $p$, the formula (35) indicate the presence of a double zero. The index for which this occurs increase when $\alpha$ keeps decreasing towards zero. These double zeros are however spurious, careful numerical investigations show that there are no double zeros. However, this abrupt change in the prediction derived from formula (35), when $\alpha$ vary smoothly, alerted us to the possibility that the root cause to why our formula predicted both the vanishing of the first zero and the existence of double zeros, is the crossing of a branch cut. Observe that the argument of the logarithm in formula (35) is an expression with complex values, so there is indeed a branch cut implied by the formula and thus the argument crossing this branch cut when $\alpha$ vary smoothly is a real possibility.

Formula (34) can be written in the form

[TABLE]

We use the standard branch of the logarithm in our calculations, and thus there is a branch cut along the negative real axis. The real part of $\xi_{0p}$ is negative, so we will have a crossing of the branch cut whenever the imaginary part of $z$ vanish. Expanding the polynomial expression defining $z$ in equation (48), and taking the imaginary part, we find that there is a crossing of the branch cut whenever $\alpha$ solves the equation

[TABLE]

Denoting $x=\ln\left(\frac{p\pi}{a\sqrt{\alpha}}\right)$ we have

[TABLE]

Solutions to eq. (50) are $x=0,\pm p\pi/2$. Investigating all three solutions we find that the one we are looking for is $+p\pi/2$ which yields

[TABLE]

Further numerical investigations show that $\alpha_{c}=\alpha_{1}$ and that $\alpha_{p}$ for $p>1$ correspond to the values of $\alpha$ where formula (35) predicts a double zero for the value of the index equal to $p$. The impact of the disappearing of the first zero on our leaky mode expansions will be discussed later, at the end of section four.

3.2 Mode shapes, normalization and projection

For $\{\xi,\xi_{0}\}$, solving equations (16),(18) with $\xi_{0}$ in the second quadrant and $\xi$ in the fourth quadrant, this is the case specified in (22), we have a leaky mode whose formula which, according to (14), is given by

[TABLE]

Using the symmetries (19),(20) we can conclude that there is a corresponding incoming, gaining mode, in the first quadrant whose formula is given by

[TABLE]

Observe that we have $(u_{p}^{+})^{*}=u_{p}^{-}$. In figure (4) we see an outgoing mode corresponding to the index $p=20$. The mode is evidently exponentially growing in $x$. This holds true for all modes, both incoming and outgoing.

Since the modes are exponentially growing in $x$, they are clearly not normalizable. We can make the modes normalizable by analytically continuing them into a complexified spatial domain, and restricting the analytically continued modes to carefully chosen complex contours. The contours will be different depending on whether the modes are incoming or outgoing. The contours we will be using are of the form

[TABLE]

where $z^{+}$ is used for the incoming modes and $z^{-}$ is used for the outgoing modes.

Evaluating eq. (59),(55) on these contours we find that they exponentially decay in both directions on the real axis. Define functions $\psi^{+}_{p}(x)$ and $\psi^{-}_{p}(x)$ on the positive real axis as

[TABLE]

The formulas for these functions are

[TABLE]

Note that for these complexified modes we also have the relation $\left(\psi^{+}_{p}(x)\right)^{*}=\psi^{-}_{p}(x)$. Figure (5) shows the functions $u^{-}_{20}(x)$ and $\psi^{-}_{20}(x)$ in the same picture. We observe that $\psi^{-}_{20}(x)$ decay exponentially outside the channel, which is confined to the interval $[-a,a]$. Also note that the complexified modes are not continuously differentiable at the points $x=\pm a$. This is because we restricted the analytically continued modes to a contour that is singular at $x=\pm a$. We made this choice in order to get fastest possible decay of the complexified modes and the simplest possible expressions for certain key differential operators acting on the modes.

Using the analyticity of the complexified modes at the two points $z=\pm a$ ,and the formulas for the two singular contours (67), it is easy to verify that the following boundary conditions holds for $\psi^{+}_{p}(x)$ and $\psi_{p}^{-}(x)$ at the two points $x=\pm a$

[TABLE]

This fact tells us that complexified modes $\psi^{+},\psi^{-}$ belong to two different spaces of functions, $V^{+}$ and $V^{-}$. Here $V^{-}$ is the space of smooth functions on real line which satisfies the boundary conditions for $\psi_{p}^{-}(x)$ (78), and similarly for $V^{+}$. We evidently have

[TABLE]

It is easy to verify that the complexified modes are in fact eigenfunctions to the differential operator

[TABLE]

We have

[TABLE]

In order use the complexified modes as a tool for expanding functions in $V^{-}$, functions that are in the span of $\left\{\psi^{-}_{p}(x)\right\}_{p=1}^{\infty}$, we need an inner product on the space. Furthermore, with respect to this inner product the leaky modes must be normalizable and orthogonal. Orthogonality would be assured if the operator $\mathcal{L}_{x}$, defined in (83), is self-adjoint with respect to the chosen inner product. This, however, seems like an impossible task, since we know that the eigenvalues $\lambda_{p}$, defined in (85), are in fact complex.

Nevertheless, an inner product that satisfy all the requirements can be constructed. In order to do this, note that for any contour $\mathcal{C}$ in the complex plane we can define a complex values scalar product on the space of functions analytic in an open set containing the contour

[TABLE]

where $\overline{\Psi}(z)$ is an analytic function defined by $\overline{\Psi}(z)=\Psi^{*}\left(z^{*}\right)$. Applying this definition of scalar product of analytic functions to the contour $z^{-}$, we get the following complex valued scalar product on the space $V^{-}$, defined for any pair of function $\psi,\phi\in V^{-}$ by the expression

[TABLE]

It is now straight forward to show that the differential operator $\mathcal{L}_{x}$ is self-adjoint with respect to the inner product (87) on the space $V^{-}$. The orthogonality of the leaky modes then follows by the familiar argument

[TABLE]

Any function in $f\in V^{-}$ which is in the span of the leaky modes $\left\{\psi^{-}_{p}(x)\right\}_{p=1}^{\infty}$ can now be expanded in terms of a generalized Fourier series of the form

[TABLE]

In a similar way an inner product can be introduced on the space of gaining modes $V^{+}$

[TABLE]

which can be used to expand gaining modes in a generalized Fourier series of the form

[TABLE]

Observe that the boundary conditions (78) implies that $\phi\in V^{-}\Leftrightarrow\psi^{*}\in V^{+}$. Thus the complex conjugate maps between these two spaces. In a similar way, the complex conjugate maps between the inner products on the two spaces

[TABLE]

The spaces of leaky modes $V^{-}$ and gaining modes $V^{+}$ are not only linear spaces, but also complex algebras. This holds because products of functions preserve the boundary conditions at $x=\pm a$. For any pair of functions in $\psi,\phi\in V^{-}$ we have for example

[TABLE]

Thus, we can conclude that $(\psi\phi)\in V^{-}$.

4 Numerical results

In this paper we will not present a formal proof specifying precisely which space of functions are in the span of the set of leaky modes, and thus for which space of functions the expansions (89) converge point wise. We will however present some arguments in section five that addresses the question of convergence of the leaky mode expansions (89).

In this section we present some numerical tests of the leaky mode expansions that will indicate strongly that they are indeed useful for the optical beam propagation context we have designed them for. In order for the leaky mode expansions to be useful for modelling (semi) transparent computational boundaries for UPPE there are two conditions that must be met.

Firstly, physically reasonable initial data must be in the span of the leaky modes. Second, products of functions in the span must also be in the span. In figure (6) we display an expansion of a Gaussian wave packet using only 30 terms in the leaky mode expansion (89). In these plots, the refractive index outside the slab is $n=1+10^{-12}$. As we can see, the Gaussian wave packet and its leaky mode expansion are indistinguishable.

Secondly, since UPPE must be able to handle nonlinear interactions, products of functions in the span must also be in the span. In order to investigate this we expanded squares of the leaky modes, namely, $\left(\psi^{-}_{11}(x)\right)^{2}$ and $\left(\psi^{-}_{50}(x)\right)^{2}$. The first propagate at an angle of $2^{\circ}$ with respect to the $z$-axis, while the second one propagate at an angle of $9^{\circ}$. The results are clearly very satisfying. Note that for the second mode we needed more terms in the leaky mode expansion because of its highly oscillatory nature.

In this paper we are not going to implement our leaky mode expansions in a fully nonlinear UPPE propagation algorithm. Before this can be done, more work has to be put into ensuring the accuracy and efficiency of the transformation from a function to its leaky mode expansion and back again. Here we will show a linear propagation example, where we compare the approach using the leaky modes, to one using regular Fourier modes, which corresponds to imposing perfectly reflecting boundary condition at at $x=\pm a$. Both are compared to the exact, infinite domain solution, which, for any $z$, can be approximated arbitrarily well by using a regular Fourier series on a much larger transverse domain. To appreciate how well our leaky modes expansion does, we demonstrate a numerical experiment where model a CW Gaussian beam propagation using Fourier expansion for finite as well as for infinite domain and compare it to the leaky modes expansion.

For the regular Fourier solution with perfectly reflecting boundary conditions at $x=\pm a$ we have

[TABLE]

where $\beta_{k}=\left(\left(\pi kx/(2a)\right)^{2}-(\omega/x)^{2}\right)^{1/2}$, and where the domain is $x\in(-a,a),z\in(0,\infty)$. As the boundary required by UPPE at $z=0$, we use the Gaussian $e(x,0)=f(x)=\text{exp}(-cx^{2})$ for some parameter $c>0$. Notice that the expansion functions in (96) form an orthogonal set

[TABLE]

for $k,l$ both being either even or odd. Imposing the Gaussian as a boundary condition at $z=0$ determine the coefficients (96), of the expansion (93)

[TABLE]

As indicated earlier, we also express the infinite domain solution using Fourier modes, now on a larger domain, say $3a$. The exact solution and this numerical solution will not deviate until the diffracting Gaussian hit the boundary of the extended domain.

Figure (8) depicts the solutions for all three approaches. The width of the domain for the finite Fourier method as well as the slab for the leaky modes is $a=10$ and the refractive index outside the slab is $n=1+10^{-15}$ in the optical regime. On (8(a)) we see that the leaky modes solution and infinite domain Fourier overlap perfectly, but approximately at this point ($z\approx 1300$), the wave hits the boundary of the slab and the finite Fourier solution starts to deviate from the other two solutions, as expected. Propagating the wave further in the slab, we observe that around $z\approx 4500$ the leaky modes solution starts to deviate from the infinite domain Fourier solution. Up to this moment, the leaky modes and infinite Fourier solutions were very close to each other. We can therefore say, that with the leaky modes method we were able to propagate the wave approximately $3.4$ times longer than with finite Fourier method. The reason why the leaky modes eventually collapsed is, that the slab is not perfectly transparent. In other words, the difference in the refractive indices for the slab and the outside domain is non-zero. This leads to reflections that gradually build up as the wave propagates in $z$ causing it to interfere with itself.

To be able to propagate the wave using leaky modes even further, we could lower the index $n$ even more, to say $n=1+10^{-18}$. However, it turns out that here we come across some serious issues. Let us first look at an expansions for a Gaussian function using $n=1+10^{-18}$. Looking at figure (9) we see that the width of the domain is the same as in the propagation example, however we made refractive index $n$ outside the slab closer to 1. The number of terms used in this expansion was 200. The badness of this expansion suggests that one should use perhaps more terms to make it better. But the truth is, the expansion does not change with more terms. Thus, the numerics indicate that the leaky mode expansion for the Gaussian does converge point wise, but unfortunately to some other function than the target Gaussian. There are two possible explanations for what happens here. The first is that the series actually diverges, but so slowly that we cannot detect it numerically. The second is that the series does converge point wise, but not to the function used to generate it.

Extensive numerical investigations, using very high numerical precision, leads us to conjecture that it is the second explanation that is correct. In fact, we suspect that the leaky mode never converge point wise to the function used to generate it. We will look more into these issues in the next section using asymptotic methods. Here we just note that even though we very likely do not have point wise convergence for the leaky mode expansion, the expansion is nevertheless for the task it was designed for. We find that the deviation between a function and its leaky mode expansion is only noticeable when the dimensionless number

[TABLE]

is not too small. For the series to give a, practically speaking, faithful representation of functions, we need at least $\eta\gtrsim 10^{-4}$. We find that for a pulse in the optical regime ($\omega/c\approx 1.58153\times 10^{7}$) we have a good representation of Gaussian initial data if

[TABLE]

An important requirement for using the leaky mode expansion is that the main part of the pulse, where the bulk of the nonlinear interactions takes place, is well inside the domain $[-a,a]$. The choices for the transverse width $a$ of the domain in the above list are chosen because they corresponds to actual dimensions used in high energy, long distance, propagation of optical pulses in air, using the UPPE code developed at the Center for Mathematical Sciences at the University of Arizona.

In the previous chapter we mentioned a phenomenon that occurs when one manipulates with the value of $\alpha$. In particular, if $\alpha$ becomes less than (51) for $p=1$, the first zero disappears. In other words, we loose the first eigenfunction completely, an eigenfunction which determine the first term in the leaky mode expansion. Since one usually expects that the first terms in the expansion are the most important ones for generic functions, the loss off the first eigenfunction is ominous. We conjecture that this loss, at least partly, explains why the leaky mode expansion loses its ability to accurately represent important boundary data like a Gaussian, when the parameter $\alpha$ become small enough.

In support of this conjecture, note that the eigenfunctions $\psi^{-}_{p}(x)$ alternate between being odd and even functions depending on the index $p$. Before the disappearing of the first zero, $\psi_{1}^{-}(x)$ is an even function. Let us denote (51) for $p=1$ as $\alpha^{*}$. Then for $\alpha\gtrsim\alpha^{*}$ the expansion is a good representation of the Gaussian, which is even. However, for $\alpha<\alpha^{*}$ the expansion goes bad because we have lost the first term in the sum. Now we understand why the expansion goes wrong. Because the most important first term in the expansion is an odd function trying to represent an even Gaussian. With this in mind, let’s expand an odd function instead of the even Gaussian for $\alpha<\alpha^{*}$. Let us for example use derivative of the Gaussian. And indeed, as can be seen in figure (10), the leaky mode expansion represent the odd functions much better than the even Gaussian. While the Gaussian was badly represented by its leaky mode expansion for $n=1+10^{-18}$, for the derivative of the Gaussian, which is an odd function, we have a very precise leaky mode expansion, even for an index step as small as $n=1+10^{-24}$.

5 Asymptotic series

In the previous section we conjectured that for small index steps, the leaky mode expansions does converge, but not to the functions used to generate the expansion. We used high precision numerical calculations to support this conjecture. In this section we will lend additional support to the conjecture by proving that in the limit of small index step, the leaky mode expansion does indeed converge point wise, but to the wrong function. The asymptotic regime we are exploring are more here conveniently defined in terms of the parameter $\alpha$. The requirement of the analysis in this section is that

[TABLE]

The validity of this inequality is what we in this section mean by the asymptotic limit. Note that (103) is in fact also the requirement for the formula (104) to be an accurate approximation to the locating of the zeros of the determinant $\xi_{0p}$. Recall that the leaky mode expansion for some function $f(x)$ is given by

[TABLE]

We will be interested in finding an asymptotic approximation to the terms in this sum for two sample functions. However, before we proceed to the actual terms for our sample functions, we first need to know, how the coefficients $A,B,C,D$ in eq. (77) depend on the index $p$ in the asymptotic limit. First of all, we realize that the vector $(A,B,C,D)^{T}$ is the null space and thus eigenvector belonging to the eigenvalue 0. To compute this eigenvector we can proceed, as we would normally do, by row-reducing the matrix equation (17), which gives us the matrix

[TABLE]

where

[TABLE]

Observe that $\eta(\xi_{0})$ contains the determinant of M in the numerator, so we get all zeros in the last row of the matrix, if $\xi_{0}=\xi_{0p}$. Using this simplification it is easy to find a basis for the one dimensional null space in the form

[TABLE]

Using the asymptotic expression (30), for the location of the zeros of the determinant $\xi_{0}=\xi_{0p}$, we find the asymptotic formula for the term $\text{exp}(\pm ia\xi_{0p})$ in the following form

[TABLE]

So the term $\text{exp}(-ia\xi_{0p})$ grows linearly in $p$, while $\text{exp}(+ia\xi_{0p})$ decays. This helps us write approximate (107) by the expression

[TABLE]

In the asymptotic limit we also clearly have

[TABLE]

so that equation (109) can be further simplified into

[TABLE]

For any given function $f(x)$ we can write the leaky mode expansion in the form

[TABLE]

where according to eq. (87)

[TABLE]

We now turn our attention to the normalization terms $N_{p}$ in equation (116). In the limit of small $\alpha$ since the zeros $\xi_{0p}$ in eq. (104) are from the 2st quadrant, we can assume $\xi_{p}\approx-\xi_{0p}$. Using this assumption, and doing the integrals in (116) exactly, we obtain

[TABLE]

where $k=A^{2}+D^{2}=B^{2}+C^{2}$. We know from eq. (108) that $\text{exp}(2ia\xi_{0p})$ decays as $p^{-2}$, thus the whole expression (117) simplifies into

[TABLE]

Let us next make a general statement about the decay rate, as a function of the index $p$, of the projection (115) of a given function $f(x)$ onto the leaky mode $\psi^{-}_{p}(x)$. Let us assume that $f(x)$ is a function that is zero at $x=\pm a$ and is $n$-times continuously differentiable. It is clear that each time we perform integration by parts in (115), i.e. differentiating $f(x)$ and integrating $\psi^{-}_{p}(x)$, we get an extra factor $i\xi_{0p}$ in the denominator. After $n$ consecutive integrations by parts we get

[TABLE]

where

[TABLE]

The term $1/\xi_{0p}^{n}$ can be approximated as $\approx(-2a)^{n}/(p\pi)^{n}$, so the asymptotic expression for eq. (119) becomes

[TABLE]

Next, we will find asymptotic approximations to (123) for our two chosen sample functions. This will give us an asymptotic approximation to the terms of the leaky mode expansion for the two sample functions.

As our first sample we choose the following triangle function

[TABLE]

whose derivative is 1 for $-d<x<0$ and -1 for $0<x<d$. In this case, the approximative coefficients $c_{p}(x)=b_{p}(x)/N_{p}$ after one integration by parts are

[TABLE]

Using the approximations from (108), we get

[TABLE]

In figures (11) and (12) we compare the asymptotic expressions for the terms in the leaky mode expansion with the exact terms calculated using high precision numerics. As we can see, there is a remarkable agreement between the values predicted by the asymptotic formulas and the exact values, even for small values of the mode index $p$. Numerically, the terms appear to approach zero fairly quickly, indicating the the series itself converge.

In order to see if, and for which values of $x$ the series converge or diverge, we write the terms in the series (128) into the following form

[TABLE]

The sum over $p$ of each term in (129) can be expressed using the polylogarithm function $Li(n,z)$, which is defined by the expression

[TABLE]

Looking at the form of the exponents of the terms in $c_{p}(x)$, the values of $n$ occurring in the polylogarithms needed to sum all the terms, are $2\pm x/a,2\pm(d+x)/a$ and $2\pm(d-x)a$. If $n$ is strictly larger than 1, the series defining the polylogarithm converge absolutely. After analyzing the various inequalities we find that we get absolute convergence of the leaky mode expansion for the triangle function only if $-a+d<x<a-d$. For the case in figure (12) this region is $-1/2<x<1/2$. However, we also have a convergence in the entire channel. In the region outside the region of absolute convergence, we have also have convergence. The convergence here is ensured by cancellations among terms spiraling towards the origin in the complex plane. In the region of the channel outside the domain $-a+d<x<a-d$, the amplitude of the terms does not decay fast enough ensure absolute convergence and the cancellation among the spiraling terms are needed for convergence. The resulting convergence is evidently only conditional.

As our second sample function we pick a Gaussian wave packet

[TABLE]

for some real numbers $m,k>0$. The asymptotic terms $c_{p}(x)$ for this case are

[TABLE]

Let the parameters $m,k$ be such that $f(x)$ has its support well inside the slab and $f(\pm a)\approx 0$. Then we can evaluate the integral analytically as

[TABLE]

where $\xi_{0p}$ is defined in (31). The terms $c_{p}(x)$ in this case decay exponentially and thus ensure that the leaky mode expansion converge for all $x$ in the channel. From (13) we see that the leaky mode expansion and the exact numerical expansion both are very close to the original Gaussian wave packet for all $x$ in the channel.

We have seen that the expansion we introduced in eq. (104) can represent a function very well as long as the parameter values are not exceeded outside their bounds. These bounds are sufficient for all practical purposes. However, looking at the expansion under such circumstances where the value of $\alpha$ is small enough, we see that the expansion is a very bad representation of the target function. Although we can’t state that we know the reason for this, we have done some preliminary investigations that points to a likely explanation.

Recall that we do have completeness for the scattering modes. Formally this is expressed by the identity

[TABLE]

Here, $\varphi_{\xi_{0}}(x)$, can be any linear combination of scattering modes. The usual way to get from the completeness for scattering modes to the completeness for the leaky modes is to analytically extend the scattering modes into the complex frequency space and then use the Cauchy theorem. This allows us to write it as a discrete sum of residues evaluated at the poles which are $\xi_{0j}$. Thus the scattering states get converted into resonant leaky modes at these points leaving us with a sum similar to the one in to (104).

In order to be more precise about this, we introduce an integration contour $\mathcal{C}$ in figure (14) that contains the zeros $\xi_{0j}$ in the second and fourth quadrant. In the scattering states, the continuity coefficients contain the determinant of the matrix (17), which contains the variable $\xi=\sqrt{\alpha+\xi_{0}^{2}}$. This is a complex square root that has a branch cut on the negative real axis. Figure (14) depicts one possible complex contour. We indicated the branch points, where the branch cut begins.

Integrating the integrand in (134) over the contour $\mathcal{C}$, we get

[TABLE]

where $\mathcal{C}_{R}$ denotes the circular part of $\mathcal{C}$, $\mathcal{C}_{r}$ is the contour part along the real axis and $\mathcal{C}_{i}$ is the one along the imaginary axis. Letting $R\to\infty$ and assuming that the contribution from the integrals over $\mathcal{C}_{R}$ and $\mathcal{C}_{i}$ vanish in the limit, we have

[TABLE]

With the aid of Cauchy theorem we write the left-hand side in (136) as

[TABLE]

Let us now obtain the expressions for the scattering states. We can write them in the form

[TABLE]

With the usual boundary conditions we end up with a system where we have one free parameter. Solving the system leaves us with the following solution

[TABLE]

where M is the matrix found in (17). Once we have obtained the scattering states, we use the relation (137) to expand the function $f(x)$ by multiplying both sides with it and integrate wrt $x^{\prime}$. This gives us the identity

[TABLE]

The goal is to match the two sums (104) and (143). Comparing them, we observe that the inner product defined as in (87) does not include any complex conjugate and neither does the sum (104). In order to match (104) and (143), we therefore need $\varphi_{\xi_{0}}(x)$ to be real. One possibility for the integrand in (134) would be

[TABLE]

where $\overline{\psi^{-}}_{\xi_{0}}(x)=\left(\psi^{-}_{\xi_{0}^{*}}(x)\right)^{*}$ to ensure the analyticity of the state as a function of complex wavenumber $\xi_{0}$.

The free parameter in (142) controls which term in the expression (144) are going to have poles at $\xi_{0j}$. In other words, it controls which terms vanish after evaluating the residues. Notice that after evaluating the residues at $\xi_{0j}$, the only survivor we want is $\psi^{-}_{\xi_{0}}(x)$. A suitable choice for achieving this is

[TABLE]

With this choice, the coefficients in (142) have common factor $\overline{\det\textbf{M}}^{\frac{1}{2}}/\det\textbf{M}^{\frac{1}{2}}$. The cross terms from $\varphi_{\xi_{0}}(x)\left(\varphi\right)^{*}_{\xi_{0}}(x^{\prime})$ vanish leaving us with one term of the form $\overline{\psi^{-}}_{\xi_{0}}(x)\overline{\psi^{-}}_{\xi_{0}}(x^{\prime})$ whose factor is $\det\textbf{M}/\overline{\det\textbf{M}}$. Since $\overline{\det\textbf{M}}$ has no zeros at the outgoing $\xi_{0j}$, this terms vanishes too upon taking the residues, leaving only the desired term $\psi^{-}_{\xi_{0}}(x)\psi^{-}_{\xi_{0}}(x^{\prime})$.

Thus, with these choices made, the completeness of the scattering states implies that the following identity holds

[TABLE]

We now equate the sums (104) and (146).

[TABLE]

where $\psi^{-}_{\xi_{0j}}(x)$ is the function defined in (141) with $\xi_{0}\to\xi_{0j}$. In this limit, the right hand side goes to zero and the matrix becomes singular and the coefficients $(A^{-},B,C,D)$ approach a basis for the null space of the matrix which leads the solution to be the nullspace of M, so (141) becomes the resonant states with the coefficients $(A^{-},B,C,D)$ whose definitions are in (142) except of the common factor $1/\det\textbf{M}$ and $a^{+}$ instead of $A^{+}$. This factor was modified by (145) and put into limit in (147). We wrote the inner product of $f(x)$ and the resonant state on the left-hand side as an integral from $-a$ to $a$ because of the compact support of $f(x)$ in this area.

It is now evident that we can make the expressions on the two sides of (147) identical only if the free parameter $a^{+}(\xi_{0})$ is chosen to be

[TABLE]

However, this choice for $a^{+}(\xi_{0j})$ is not an analytic function because both functions $\det^{\prime}\textbf{M}(\xi_{0j})$ and $\left(\psi^{-}_{\xi_{0j}}(x),\psi^{-}_{\xi_{0j}}(x)\right)^{-}$ are zero inside the integration contour. Each of these two families of countably many zeros give rise to equally many branch cuts. The parametric formulas for these branch cuts are possible to find but while applying the Cauchy theorem we must now include terms representing integrals around all these additional branch cuts. Thus what we get from the Cauchy theorem is that any function with compact support inside that channel is equal to its leaky mode expansion, plus additional terms that includes integrals along the branch cuts on and off the imaginary axis as described above. What we know is that, unless alpha is smaller than the critical value $\alpha^{*}$, which we introduced in section four, the function is well represented by the leaky mode expansion alone. This means that the contribution from all the other terms for such values of alpha are negligible. For smaller values of alpha the contributions from the rest of the terms are not negligible and as $\eta=a^{2}\alpha$ approaches zero, these terms will come to dominate. For such values of $\alpha$ the leaky mode expansion still converge, but it does not converge to the function used to construct the expansion. By deriving asymptotic formulas for the all the terms defined by integrals around branch cuts, in the limit when $\eta$ approaches zero, one could compare their sizes and identify the dominant ones. If, say, one term dominate, then this term could be added to the leaky mode expansion resulting in an expansion that represents the function to be expanded in a much better way than the leaky mode expansion is able to do on its own. We believe that the asymptotic expressions for the terms could be found, but there might not be a dominant term, and even if there is, extending the leaky mode expansion by adding this term might easily make the expansion too hard to use for practical calculations.

6 Conclusion

In this paper we have presented an new approach to minimizing the reflections from finite computational boundaries for wave equations formulated in the spectral domain. This approach is based on representing the field in the transverse spatial direction using leaky mode expansions supported by an artificial index channel. We have shown that at the linear level, our approach makes it possible propagate the waves much further than what is possible if a regular Fourier expansion is used. The leaky modes are not reflectionless at the boundary, and eventually the small but finite reflections build up, and the computed solutions starts to deviate from the infinite domain solutions. This reflection can be minimized by reducing the index step, but at the price of getting a progressively worse representation of the solution to the wave equation. In section four and five we have argued, using both numerical and analytical approaches, that a practically useful trade off can be made between minimizing reflections from the boundary and maximizing the accuracy of the representation of solutions of the wave equations using leaky modes.

We have illustrated our approach using the case of a TE electromagnetic wave in vacuum, but the approach can clearly be generalized to much more general wave propagation problems than this. In the optical context the obvious next step would be to consider waves with cylinder symmetry. An important issue that we have not discussed in this paper is how to compute the transformation from fields to leaky mode amplitudes and back in an accurate, stable and efficient way.

7 Acknowledgments

The authors are thankful for support from the Department of mathematics and statistics at the Arctic University of Norway, from the Arizona Center for Mathematical Sciences at the University of Arizona, and for the support from the Air Force Office for Scientific Research under Grant No. FA9550-19-1-0032