Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group
Journal Home About Issues in Progress Current Issue All Issues Feature Issues

Open boundary conditions for the simulation of leaky modes

Jonathan Heinz and Miroslav Kolesik

Open Access Open Access

Abstract

We propose an open-boundary method for the simulation of the modes of confining dielectric structures. The technique is inclusive of normal modes, but is especially advantageous for the simulation of quasi-normal, or leaky, modes. The central idea is to utilize the asymptotic form of targeted solutions to eliminate the outer part of the computational domain and bring the numerical boundary close to the simulated structure. While a similar approach was previously demonstrated for scalar quantum models, here we put forward a generalization for fully vectorial fields. Accuracy in this new context is validated using step-index and tube-type hollow core fiber geometries. The method has broad applicability, as quasi-bound modes are intrinsic to many systems of interest in optics and photonics.

© 2023 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

Open systems are widespread in optics and photonics, as well as many other areas of physics. Though perfectly bound solutions are typically more tractable for analytical and numerical treatments, in practice no modes are truly ideal and even high loss resonances can have uniquely desirable properties. As such, the ability to accurately simulate the modal properties and loss of resonances with incomplete confinement is a powerful tool. Resonances of open systems have a wide variety of names even within optics and photonics, including Siegman’s unstable resonator modes [1] and leaky mode resonances (LMR) [2]. The concept of leaky modes has found broad recent application, ranging from sensing [3,4] to holography [5] and nonlinear optics [6,7]. In terms of high energy laser systems, leaky modes are especially useful for enhancing the guidance of high-power pulses using plasma waveguide channels [8–12] as these structures are intrinsically leaky and only support modes with high loss.

Because of the wide range of applications, there has long been interest in mode simulations for open systems. As the modal loss of a resonance is sensitively connected to its character far from the origin, the accuracy of any implementation is highly dependent on the method of truncation for the numerical domain. The introduction of stretched-coordinate perfectly matched layers (PML) [13] and the broad use of finite-element-based multiphysics simulation tools like COMSOL have allowed easier modal analysis for a wide range of open systems. PML are flexible and assume very little about the form of the solution, but they require additional absorbing grid points and a large window to accurately treat systems with long-range resonances. The prevalence of these methods has led to work investigating their numerical properties in the context of open systems [14], such as the tendency of PML to generate a large amount of spurious modes due to the non-physical nature of the boundary layer [15]. Simplistic transparent boundary approaches based on plane waves can enable smaller computational domains [16], but are often inflexible and quickly lose accuracy for more complex systems. Though many numerical approaches are aimed at finding individual modes, there are also cases where knowledge of the full mode expansion is important. Much work has been done to evaluate the characteristics of the full set of so-called quasinormal modes [17–19], as well as to explore its use for the modeling of various types of electromagnetic resonators [20,21].

To address the limitations of existing boundary approaches for the simulation of leaky modes, we present a transparent boundary method that allows both high accuracy and a substantially reduced computational window. The method is based on the separation and analytic solution of the long-range escaping character of individual modes, allowing the boundary to be placed almost arbitrarily close to the furthest extent of the structure. This technique has been previously explored in a scalar quantum mechanics framework, with the aim of calculating ionization rates from Stark resonances [22] under the Metastable Electronic State Approach [23]; while the method showed high accuracy and numerical stability in that work, it has not yet been developed for fully vectorial systems. As vectorial considerations are typically subtle and are necessary for many systems in optics and photonics, it is important to clearly formulate and validate the approach in this new context.

The following section states the problem and gives motivation for the method of solution. Section 2 presents the steps necessary for establishing solutions in the boundary region and matching them to a numerical discretization of the system in the interior region. Implementation details for our numerical experiments are discussed in Section 3, while Section 4 presents validation results using waveguide systems with known vectorial mode solutions.

1. Formulation

Assuming time-harmonic electromagnetic fields, Maxwell’s equations can be used to form the following vectorial wave equation in three dimensions

$$\nabla^2 \vec{\mathrm{E}}(\vec{\mathrm{r}}) + k_0^2 n(\vec{\mathrm{r}})^2 \: \vec{\mathrm{E}}(\vec{\mathrm{r}}) - \nabla \left( \nabla \cdot \vec{\mathrm{E}}(\vec{\mathrm{r}}) \right) = 0,$$
where $k_0 = \frac {2 \pi }{\lambda }$ is the vacuum wavenumber and $n(\vec {\mathrm {r}})$ represents a spatially-varying refractive index. As the resonances found in open systems that derive from this equation have a wide variety of names and conventions, we will use the concept of propagating waveguide modes for the purpose of discussion. This allows a more specific description of the properties of solutions and how we treat them, but the considerations are the same for systems with other formulations and methods of confinement.

Propagating waveguide modes are typically confined in one or both transverse dimensions, and are each characterized by a propagation constant and an associated effective refractive index. While fully bound modes have real effective indices that determine the phase evolution of the mode without loss, leaky modes can have substantial loss that is captured by the imaginary part of the effective index. This loss is due to any outgoing portions of the solution that may exist beyond the limits of the confining structure, so the imaginary part of the effective index is sensitively connected to the transverse profile of the mode at long distances. This means that to accurately simulate leaky modes we must carefully account for both the confinement from the structure and the outgoing behavior of the solution at long ranges.

To accomplish this, the domain is divided by a closed boundary into two parts: an interior region that contains any index variation associated with the structure, and an exterior region where we can assume a constant background index $n_0$. The boundary can in general be of any shape, but considering that this choice only changes the mechanics of the parameterization we will use polar coordinates and a circular boundary of radius $r_B$. A schematic of this arrangement is shown in Fig. 1 for a general dielectric waveguide cross-section.

 figure: Fig. 1.

Fig. 1. Illustration of a waveguide cross-section with the domain divided into interior and exterior regions. The grey dash-dotted line represents the circular boundary at $r_B$, and a hypothetical three-point finite difference stencil is overlaid. The large dot represents a discrete boundary point, the small dots represent neighboring points, and the open circle represents a "ghost" point that must be expressed in terms of interior points.

Download Full Size | PDF

The exterior region is solved analytically for propagating modes, yielding expressions that are parameterized by the propagation constant $\beta$. As the exterior region in the transverse plane has a uniform background, the fields there can be expanded as a sum of solutions with different angular order number $\nu$.

$$\vec{\mathrm{E}}(r,\phi,z) = e^{{-}i \beta z} \,\vec{E}(r,\phi) = \sum_{\nu} e^{i (\nu \phi - \beta z)} \, \vec{\mathcal{E}}^{(\nu)}(r)$$
where $\nu$ is an integer ranging from $-\infty$ to $\infty$ and $\mathrm {E}$, $E$, and $\mathcal {E}$ have been used to delineate fields with three-, two-, and one-dimensional spatial dependence, respectively. Arrows will be used to refer to vectors with components in each spatial direction, while bolded symbols will be used when only the transverse components are being taken into account. Using this expansion, radial solutions are then found for each $\nu$. Though there are six total electric and magnetic field components, only two components must be solved to be able to calculate the rest; for our derivation we will take the two transverse components of the electric field. Once the full form of the terms in Eq. (2) are established, the interior region can be treated numerically and the boundary can be considered; as Eq. (2) is in spectral representation for the angular direction the interior fields must be similarly expanded, matched to the exterior, and transformed back to a real-space representation. Using a finite-difference approach as an example of a discretization of the interior domain, as illustrated in Fig. 1, Eq. (2) can be used to express points past the boundary that are unknown to the interior treatment in terms of points on the boundary. This leads to the following boundary expression for the two transverse components of the electric field, represented by the column vector $\boldsymbol {\mathcal {E}} (r, \phi )$:
$$\boldsymbol{\mathcal{E}} (r_0, \phi) = \frac{1}{2 \pi} \int_{\phi'} \sum_{\nu} e^{i \nu (\phi - \phi')} \mathbf{M}_{\nu}(r_0) \mathbf{M}^{{-}1}_{\nu} (r_B) \boldsymbol{\mathcal{E}}(r_B, \phi') d \, \phi' ,$$
where $\mathbf {M}_{\nu }$ is a matrix representing the external solutions for the transverse field components, which we derive in what follows. This matrix is parameterized both by an individual value of $\nu$ and the propagation constant $\beta$ for the mode as a whole. As such, a guess for the propagation constant or effective index is necessary to be able to evaluate the external solutions. The approximated expressions can then be used as estimated boundary conditions to calculate a new effective index for the interior region. This value is then used to more accurately evaluate the external expressions, and the process is repeated until the two effective indices converge. This results in a full description of the solution for the entire domain. This approach was used to good effect in a scalar context for the simulation of Stark resonances in hydrogen [22], but it has not yet been demonstrated that the same reduction in boundary distance and stable convergence will be seen when simulating multiple simultaneous components of vectorial electromagnetic fields.

2. Boundary region solution

In the region outside of the circular boundary, any derivatives related to the index variation dissapear and (1) reduces to the vector Helmholz equation

$$\nabla^2 \vec{\mathrm{E}}(\vec{\mathrm{r}}) + k_0^2 n_0^2 \: \vec{\mathrm{E}}(\vec{\mathrm{r}}) = 0 .$$

In cylindrical coordinates the z-component of the electric field couples only to itself, so we can write the scalar wave equation for $E_z$ and use the assumptions given by Eq. (2) to simplify angular and longitudinal derivatives.

$$r^2 \frac{\partial^2 \mathcal{E}^{(\nu)}_z}{\partial r^2} + r \frac{\partial \mathcal{E}^{(\nu)}_z}{\partial r} + r^2 \left( k_0^2 n_0^2 - \beta^2 - \frac{\nu^2}{r^2} \right) \mathcal{E}^{(\nu)}_z = 0$$

This equation can be solved by a combination of Bessel functions. For outgoing wave solutions, we choose the Hankel function of the first kind. We can also write the equation for the z-component of the magnetic field in the same way

$$\begin{aligned} \mathcal{E}^{(\nu)}_z(r) &= A_{\nu} (J_\nu (\kappa r) + i Y_\nu (\kappa r) ) = A_{\nu} H_\nu^{(1)}(\kappa r),\\ \mathcal{H}^{(\nu)}_z(r) &= B_{\nu} (J_\nu (\kappa r) + i Y_\nu (\kappa r) ) = B_{\nu} H_\nu^{(1)}(\kappa r), \end{aligned}$$
where $\kappa$ is defined by the relation $\kappa ^2 = k_0^2 n_0^2 - \beta ^2$. Using Maxwell’s equations, we can find expressions for the remaining electric field components in terms of $\mathcal {E}^{(\nu )}_z$ and $\mathcal {H}^{(\nu )}_z$
$$\begin{aligned} \mathcal{E}_r &= \frac{-i \beta}{\kappa^2} \left(\frac{i \nu \omega \mu}{\beta r} \mathcal{H}^{(\nu)}_z + \frac{\partial \mathcal{E}^{(\nu)}_z}{\partial r} \right)\\ \mathcal{E}_{\phi} &= \frac{-i \beta}{\kappa^2} \left(\frac{i \nu}{r} \mathcal{E}^{(\nu)}_z - \frac{\omega \mu}{\beta} \frac{\partial \mathcal{H}^{(\nu)}_z}{\partial r} \right) , \end{aligned}$$
where $\omega$ is the optical frequency and $\mu$ is the vacuum permeability. Expanding this expression using the solutions in Eq. (6) results in a system equations containing Hankel functions and the two unknown coefficients $A_{\nu }$ and $B_{\nu }$. We can write this in condensed matrix form
$$\boldsymbol{\mathcal{E}}^{(\nu)}(r) = \mathbf{M}_{\nu} (r) \mathbf{C}_{\nu}, \quad \mathbf{M}_{\nu}^{{-}1} (r) \boldsymbol{\mathcal{E}}^{(\nu)}(r) = \mathbf{C}_{\nu} ,$$
where
$$\boldsymbol{\mathcal{E}}^{(\nu)}(r) = \begin{bmatrix} \mathcal{E}_r^{(\nu)}(r) \\ \mathcal{E}_{\phi}^{(\nu)}(r) \end{bmatrix}, \quad \mathbf{C}_{\nu} = \begin{bmatrix} A_{\nu} \\ B_{\nu} \end{bmatrix},$$
and $\mathbf {M}_{\nu } (r)$ represents the matrix of Hankel functions. With this expression, we can eliminate the unknown coefficients and find a relationship between the field at two different radii. As we want to apply the exterior solutions as boundary conditions for a general interior solution, we use Eq. (8) to find a relationship between the internal field at boundary point $r_B$ and another point $r_0$. This relationship can be used for the truncation of the interior domain.
$$\boldsymbol{\mathcal{E}}^{(\nu)} (r_0) = \mathbf{M}_{\nu} (r_0) \mathbf{C}_{\nu} = \mathbf{M}_{\nu}(r_0) \left( \mathbf{M}^{{-}1}_{\nu} (r_B) \boldsymbol{\mathcal{E}}^{(\nu)}_{int} (r_B) \right)$$

We can see that when $r_0 \rightarrow r_B$ the two matrices become direct inverses and combine to form the identity matrix. We now have a connection between fields at the boundary, but one that is expressed in terms of $\nu$ rather than $\phi$. To give relations that are in real space for all dimensions, we need to transform back to real space in terms of $\phi$. Explicitly writing out the transformations between real space and spectral representations, and performing a few steps of simplification, results in the final real space equation given in Eq. (3). It is important to note that in principle the $\phi '$ integral in this expression causes each boundary point to be coupled to every other boundary point, though the coupling will still be dominated by the nearest neighbors. We also must truncate the sum over $\nu$ to a finite number of orders for it to be numerically realizable.

With the general form of the boundary equation in place, we can find a specific expression for the Hankel function matrix product. Plugging the $z$ component solutions of Eq. (6) into Eq. (7), we find

$$\mathbf{M}_{\nu} (r) = \begin{bmatrix} - \frac{i \beta}{\kappa} H_{\nu}^{\prime (1)} (\kappa r) & \frac{\omega \mu \nu}{\kappa^2 r} H_{\nu}^{(1)} (\kappa r) \\ \frac{\beta \nu}{\kappa^2 r} H_{\nu}^{(1)} (\kappa r) & \frac{i \omega \mu}{\kappa} H_{\nu}^{\prime (1)} (\kappa r) \end{bmatrix}$$

With this matrix, we can find the matrix product in Eq. (3) and express it as follows:

$$\mathbf{M}_{\nu}(r_0) \mathbf{M}^{{-}1}_{\nu} (r_B) = C(r_B, \nu) \mathbf{M}' (r_0, r_B, \nu) ,$$
where $\mathbf {M}'$ is given by
$$\mathbf{M}' (r_0, r_B, \nu) = \begin{bmatrix} \mathrm{M}'_d (r_0, r_B, \nu) & \mathrm{M}'_o (r_0, r_B, \nu) \\ - \mathrm{M}'_o (r_0, r_B, \nu) & \mathrm{M}'_d (r_0, r_B, \nu) \end{bmatrix}$$
and the scalar factor $C$, diagonal coefficient $\mathrm {M}'_d$, and off-diagonal coefficient $\mathrm {M}'_o$ are defined as follows.
$$C(r_B, \nu) = \frac{1}{\left( \frac{\nu^2}{\kappa^2 r_B^2} H_{\nu}^{(1)} (\kappa r_B)^2 - H_{\nu}^{\prime (1)} (\kappa r_B)^2 \right)} \, ,$$
$$\mathrm{M}'_d (r_0, r_B, \nu) = \frac{\nu^2}{\kappa^2 r_B r_0} H_{\nu}^{(1)} (\kappa r_B) H_{\nu}^{(1)} (\kappa r_0) - H_{\nu}^{\prime (1)} (\kappa r_B) H_{\nu}^{\prime (1)} (\kappa r_0) \, ,$$
$$\mathrm{M}'_o (r_0, r_B, \nu) = \frac{i \nu}{\kappa r_0} H_{\nu}^{\prime (1)} (\kappa r_B) H_{\nu}^{(1)} (\kappa r_0) - \frac{i \nu}{\kappa r_B} H_{\nu}^{(1)} (\kappa r_B) H_{\nu}^{\prime (1)} (\kappa r_0) .$$

The different angular components are coupled through the use of a single guess for the unknown propagation constant $\beta$, which is embedded in these equations through its relationship with $\kappa$. This means that an approximated eigenvalue is necessary to be able to use the exterior solutions to establish boundary conditions for a treatment of the interior region.

3. Implementation

As our approach has not yet been demonstrated for vectorial systems, we created an accompanying implementation using the assumption of radially symmetric waveguide geometries. This was done with the goal of avoiding some of the spatial considerations necessary for a full range of two-dimensional cross sections, as the symmetry reduces the interior domain to a single radial dimension and one choice of the angular mode number $\nu$. This allows the relations shown in (13–16) to be tested and validated for a single set of parameters at a time, and facilitates a more direct comparison between calculated values and known solutions as is discussed in Section 5.

To solve the interior domain allowing for sharp index contrast, we start with the wave equation for transverse components of the field envelope. Only the assumption that the index is constant in $z$ is necessary for the transverse components to remain uncoupled from the propagation direction.

$$\nabla_{{\perp}}^2 \mathbf{E} + k_0^2 n^2 \mathbf{E} = \nabla_{{\perp}} \left( \nabla_{{\perp}} \cdot \mathbf{E} - \frac{1}{n^2} \nabla_{{\perp}} \cdot \left(n^2 \mathbf{E} \right) \right)$$

Here, $\nabla _{\perp }$ is used to represent the transverse components of the vector Laplace operator. Assuming the index has no angular variation, we can reduce this equation to the following coupled equations for the transverse fields:

$$\nabla^2 \mathrm{E}_r - \frac{\mathrm{E}_r}{r^2} - \frac{2}{r^2} \frac{\partial \mathrm{E}_{\phi}}{\partial \phi} + \frac{\partial}{\partial r} \left[ \frac{1}{n^2} \frac{\partial}{\partial r} \left( n^2 \mathrm{E}_r \right) \right] - \frac{\partial^2 \mathrm{E}_r}{\partial r^2} + k_0^2 n^2 \mathrm{E}_r = 0$$
$$\nabla^2 \mathrm{E}_{\phi} - \frac{\mathcal{E}_{\phi}}{r^2} + \frac{2}{r^2} \frac{\partial \mathrm{E}_r}{\partial \phi} + \frac{1}{r} \frac{\partial}{\partial \phi} \left[ \frac{1}{n^2} \frac{\partial}{\partial r} \left( n^2 \mathrm{E}_r \right) \right] - \frac{1}{r} \frac{\partial^2 \mathrm{E}_r}{\partial \phi \partial r} + k_0^2 n^2 \mathrm{E}_{\phi} = 0$$

Selecting an individual solution with a single $\nu$ from Eq. (2) to simplify derivatives in the angular and longitudinal directions, we find

$$\frac{1}{r} \frac{\partial \mathcal{E}_r}{\partial r} - \frac{2 i \nu}{r^2} \mathcal{E}_{\phi} + \frac{\partial}{\partial r} \left[ \frac{1}{n^2} \frac{\partial}{\partial r} \left( n^2 \mathcal{E}_r \right) \right] + \left(k_0^2 n^2 - \frac{1 + \nu^2}{r^2} \right) \mathcal{E}_r = \beta^2 \mathcal{E}_r$$
$$\frac{1}{r} \frac{\partial \mathcal{E}_{\phi}}{\partial r} + \frac{\partial^2 \mathcal{E}_{\phi}}{\partial r^2} + \frac{2 i \nu}{r^2} \mathcal{E}_r + \frac{i \nu}{r} \left[ \frac{1}{n^2} \frac{\partial}{\partial r} \left( n^2 \mathcal{E}_r \right) - \frac{\partial \mathcal{E}_r}{\partial r} \right] + \left(k_0^2 n^2 - \frac{1 + \nu}{r^2} \right) \mathcal{E}_{\phi} = \beta^2 \mathcal{E}_{\phi}$$
which can then be formed into a matrix eigenvalue equation.
$$\boldsymbol{\mathcal{E}} = \begin{bmatrix} \mathcal{E}_r \\ \mathcal{E}_{\phi} \end{bmatrix} , \quad \begin{bmatrix} \mathbf{A}_{rr} & \mathbf{A}_{r \phi}\\ \mathbf{A}_{\phi r} & \mathbf{A}_{\phi \phi} \end{bmatrix} \boldsymbol{\mathcal{E}} = \beta^2 \boldsymbol{\mathcal{E}}$$
where the matrix elements are given by
$$\mathbf{A}_{rr} \mathcal{E}_r = \frac{1}{r} \frac{\partial \mathcal{E}_r}{\partial r} + \frac{\partial}{\partial r} \left[ \frac{1}{n^2} \frac{\partial}{\partial r} \left( n^2 \mathcal{E}_r \right) \right] + \left( k_0^2 n^2 - \frac{1 + \nu}{r^2} \right) \mathcal{E}_r$$
$$\mathbf{A}_{\phi \phi} \mathcal{E}_{\phi} = \frac{1}{r} \frac{\partial \mathcal{E}_{\phi}}{\partial r} + \frac{\partial^2 \mathcal{E}_{\phi}}{\partial r^2} + \left( k_0^2 n^2 - \frac{1 + \nu}{r^2} \right) \mathcal{E}_{\phi}$$
$$\mathbf{A}_{r \phi} \mathcal{E}_{\phi} ={-} \frac{2 i \nu}{r^2} \mathcal{E}_{\phi}$$
$$\mathbf{A}_{\phi r} \mathcal{E}_r = \frac{2 i \nu}{r^2} \mathcal{E}_r + \frac{i \nu}{n^2 r} \frac{\partial n^2}{\partial r} \mathcal{E}_r$$
and the effective index can be calculated from the propagation constant using
$$\beta^2 = n_{eff}^2 k_0^2 .$$

The interior domain is discretized on a radial domain of $N$ points with spacing $\Delta r$. To refer to specific grid points and discrete function values, we use the following notation:

$$r^n = \Delta r \left[ \frac{1}{2} + (n - 1) \right], \quad f^{n} = f(r^{n}),$$
where $n \in [1,N]$ and $r^N = r_B$. Derivatives are approximated using three-point centered finite difference approximations, and the origin is omitted from the grid to avoid the central discontinuity. A more detailed look at the methods of discretization can be found in Appendix A. As we are selecting a mode with well defined angular dependence associated with a single choice of $\nu$, a simplified version of Eq. (10) can be used to find relations between the fields at the boundary point and fields at the "ghost" point past the boundary.
$$\mathcal{E}^{N+1}_r = C^{N, \nu} \left( \mathrm{M}_d^{\prime N, N+1, \nu} \ \mathcal{E}^{N}_r + \mathrm{M}_o^{\prime N, N+1, \nu} \ \mathcal{E}^{N}_{\phi} \right)$$
$$\mathcal{E}^{N+1}_{\phi} = C^{N, \nu} \left( - \mathrm{M}_o^{\prime N, N+1, \nu} \ \mathcal{E}^{N}_r + \mathrm{M}_d^{\prime N, N+1, \nu} \ \mathcal{E}^{N}_{\phi} \right)$$

These boundary conditions illustrate the necessity of a vectorial treatment when $\mathrm {M}'_o$ is non-negligible. It is worth noting that both depend on a single value of $\nu$ as well as the propagation constant for the mode. This makes the linear eigenvalue problem given in Eq. (22) into a nonlinear eigenvalue problem. To find a mode using this approach, we construct a system matrix using both a constant $\nu$ based on the targeted mode profile and a guess for the eigenvalue. The resulting linear system is then solved using FEAST, a subspace iteration eigensolver [24]. The newly calculated eigenvalue is used as an updated guess, and then the process is repeated until the guess and output match to within a given threshold. Only a few iterations are typically necessary to converge to a solution. The iteration procedure performs well in practice, and pairs well with the use of a subspace eigensolver. We use the same guess from the iteration procedure to help define the area in eigenvalue space where the solver will look for solutions.

4. Validation and analysis

4.1 Step-index fiber

Though the boundary conditions were formulated to allow for the solution of leaky modes, it is first important to show that they can be used to accurately simulate more traditional bound modes. To examine the case where the outgoing wave condition becomes evanescent, a sensible starting point is the basic step-index fiber. This structure has a known and easily solvable set of bound modes that still exhibit polarization-dependent behavior. We generated analytic mode results for the lowest order TE, TM, and HE modes following the traditional boundary matching approach [25] for a step-index structure. The profile had a core index of 1.5 and a background index of 1.45 in all cases. Mode eigenvalues were calculated using our finite difference implementation and used to find effective index values through the relationship given in Eq. (27). All of the following numerical eigenvalue results will be presented in this effective index form. Figure 2 shows a comparison of effective index values for the fundamental modes with a range of core radii $R$, assuming a wavelength $\lambda$ of 800 nm and a grid spacing $\Delta r$ of 5 nm. Table 1 shows the convergence of an individual simulation for the fundamental HE mode with a radius of 4 µm, using the same parameters and the core index as the initial guess and requiring that the guess and output index match to eight significant digits. Note that some small deviation from the reference results is to be expected in general due to the finite discretization of the index structure boundaries.

 figure: Fig. 2.

Fig. 2. Calculated effective indices of the fundamental TE, TM, and HE resonances for a step-index fiber with varying radii. $\lambda = 800 \ \mathrm {nm}, \ n_{core} = 1.5, \ n_{clad} = 1.45$. Reference values calculated using the boundary matching approach [25].

Download Full Size | PDF

Tables Icon

Table 1. Convergence, R = 4 $\mu \mathrm {m}$

View Table | View all tables in this article

We see a close agreement between our numerical and reference results, and in this simple system the solution converges in only an iteration or two. As the most extreme indices of the waveguide generally give upper and lower bounds on the effective indices of of modes that are supported by the waveguide, those indices can be used as well-motivated guesses for the effective indices of the highest and lowest order modes.

4.2 Hollow core fiber

The next step is to demonstrate the solution of leaky modes, which we chose to examine using the air-core, air-cladding "tube"-type hollow core fiber. One of the simplest geometries amongst hollow core fibers, the tube fiber has no perfectly bound modes and is often used as a simplified model of other more complicated anti-resonant geometries [26]. For comparison, we use the approximate analytical model presented by Zeisberger and Schmidt [27]. In all cases, the glass tube index is assumed to be 1.45, with a core and background index of 1. Figure 3 shows comparisons of real and imaginary effective index values for the fundamental modes at a range of inner core radii $R$, assuming a wavelength of 1.2 µm, a tube thickness $w$ of 0.7 µm, and a grid spacing $\Delta r$ of 5 nm. The real part of the effective index is shown as a deviation from 1, while the imaginary part is shown in log scale.

 figure: Fig. 3.

Fig. 3. Calculated effective indices of the fundamental TE, TM, and HE resonances for a hollow core fiber of varying radii. $\lambda = 1.2 \ \mathrm {\mu m}, \ w = 0.7 \ \mathrm {\mu m}, \ n_{tube} = 1.45, \ n_{core} = n_{clad} = 1$. Reference values calculated using the model given by Zeisberger and Schmidt [27].

Download Full Size | PDF

The thickness of the tube mostly determines the placement of the resonant and anti-rensonant points in the spectrum; the wavelength for our results was chosen to be near a loss minimum and anti-resonant point pair that lies in the neighborhood of 900 nm [27]. It is also important to note that the fundamental HE mode is the lowest loss mode for this region of parameter space. This is often desirable, as it is easier to match an input beam to this mode’s Gaussian-like profile and achieve high coupling efficiency than it is to match the ring shape of the fundamental TE and TM modes. Table 2 shows the convergence of the iteration procedure to nine significant digits in both real and imaginary parts for the HE mode at a radius of 20 µm.

Tables Icon

Table 2. Convergence, R = 20 $\mu \mathrm {m}$

View Table | View all tables in this article

Even in this more complicated system, we can see five stable digits of the imaginary part after just three iterations. This illustrates just how well the iteration procedure performs in practice. This matches with what we have shown previously for scalar systems, confirming that the vectorial formulation does not substantially impact the performance of the iteration procedure.

For further validation of the numerical characteristics of this implementation, the dependence of calculated values on boundary placement and grid spacing must be investigated. Figure 4(a) shows the impact of varying the spacing between the grid points $\Delta r$ on the more sensitive imaginary part of the effective index, while Fig. 4(b) shows the stability of the imaginary part with respect to the placement of boundary; the stability of the result is measured as a deviation from the average value for all boundary placements, which are specified using the offset $r_B - R$ between the numerical boundary and the ring radius. For different offsets, grid points are added or removed from the end of the domain to keep the placement of the remainder of the grid points constant. All physical parameters were kept constant.

 figure: Fig. 4.

Fig. 4. Numerical studies of the finite difference implementation using the $\mathrm {HE_{11}}$ mode and a 20 µm tube waveguide. (a) Imaginary part of the effective index for a range of grid spacings. (b) Deviation of the imaginary part of the effective index for a range of boundary offsets.

Download Full Size | PDF

The grid spacing shows a well-controlled relationship with the accuracy of the solution, but care must be taken to only directly compare results for grid spacings that are evenly divisible by the tube thickness. The variation of the solution in terms of boundary placement is random and on the order of $10^{-13}$, with no trend for boundaries placed closer to or further away from the structure. Considering that this is essentially at the noise floor of the solver, any dependence between the boundary and the accuracy of the solution is small enough to be hidden by noise. This means that for this system a boundary offset as small as $1 \mu m$ can be used with roughly the same accuracy as a boundary offset of $20 \mu m$, highlighting the potential of this method to dramatically reduce boundary distances and the computational cost of simulations dealing with leaky modes.

5. Conclusion

We have presented a flexible open-boundary approach for the solution of vectorial leaky modes in optics and photonics, and found it to be accurate and stable even for boundaries placed close to the structure of interest. A formulation of the approach for vectorial field envelopes was given in cylindrical coordinates, assuming a radial boundary and an interior region that completely contains the structure of interest. A radially symmetric implementation was detailed, which connects the boundary equations to a finite-difference interior domain and utilizes a simple iteration procedure for convergence. Calculated results were validated using two waveguide systems with known solutions. The implementation demonstrated the substantial capacity of the boundary method to reduce the grid size in simulations of vectorial waveguide modes. The iteration method also exhibited convergence in a small number of iterations, which is similar to previous demonstrations in scalar systems [22].

Further work will incorporate a generalized formulation of the method to be able to simulate a full range of two dimensional index profiles. However, the implementation demonstrated here can already treat many useful systems without modification. The modal analysis of GRIN fibers [28], radial Bragg fibers [29], and many other structures can often benefit greatly from tools that that can quickly simulate both bound and leaky modes. Plasma waveguide channels are of particular interest, as in many cases they are radially symmetric and weakly guiding [8–12].

Appendix A: Radially symmetric Helmholz discretization

First and second derivatives for both the field and index are discretized as using the typical symmetric three-point stencil, involving a point and its two nearest neighbors. The combined term in Eq. (23) that involves both the field and index together is more complicated, and is discretized as follows. Replacing $n^2$ with $\epsilon$ for now for ease of notation, we will approximate the exterior derivative at symmetric points halfway towards the neighbor points.

$$\frac{1}{\Delta r} \left[ \left( \frac{1}{\epsilon} \frac{\partial \epsilon \mathcal{E}_r}{\partial r} \right)_{i+1/2} \! - \: \left( \frac{1}{\epsilon} \frac{\partial \epsilon \mathcal{E}_r}{\partial r} \right)_{i-1/2} \right] ,$$

Then similar differences are performed on the interior derivatives, resulting in a relation that only involves the field at the central and nearest neighbor points.

$$\frac{1}{\Delta r^2} \left[ \left(\frac{1}{\epsilon}\right)^{i+1/2} \left(\epsilon^{i+1} \mathcal{E}_r^{i+1} - \epsilon^i \mathcal{E}_r^i \right) - \left(\frac{1}{\epsilon}\right)^{i-1/2} \left(\epsilon^i \mathcal{E}_r^i - \epsilon^{i-1} \mathcal{E}_r^{i-1} \right) \right]$$

Since the index is typically defined formulaically, this equation could be used directly. However, if the index must be defined on the same grid as the field, we can make the approximation that the function $1/\epsilon$ between points is the average of the values at the two neighboring points, i.e.

$$\left(\frac{1}{\epsilon}\right)^{i+1/2} \approx \frac{1}{2 \epsilon^{i+1}} + \frac{1}{2 \epsilon^i}$$

Gathering coefficients for each field point, this leads to the following discretization for combined term of the eigenvalue equation.

$$\frac{1}{2 \Delta r^2} \left[ \frac{\epsilon^{i+1} + \epsilon^i}{\epsilon^i} \mathcal{E}_r^{i+1} + \frac{\epsilon^i + \epsilon^{i-1}}{\epsilon^i} \mathcal{E}_r^{i-1} - \left( \frac{\epsilon^{i+1} + \epsilon^i}{\epsilon^{i+1}} + \frac{\epsilon^i + \epsilon^{i-1}}{\epsilon^{i-1}} \right) \mathcal{E}_r^i \right]$$

It can be easily seen that if the index is constant in the neighborhood of the given point, this equation reduces to the standard second derivative equation.

For the central boundary, all that is necessary is knowledge of the symmetry across the origin. Since the symmetry about the origin in the transverse plane depends on the angular order number, we have the following boundary conditions

$$\begin{cases} \mathcal{E}_{r,\phi}(0) = 0, \quad & \nu \neq~0\\ \left. \frac{\partial \mathcal{E}_{r,\phi}}{\partial r} \right|_{r = 0} = 0, \quad & \nu = 0 \end{cases}$$

Which lead to relations for the "ghost" point at $-\Delta r/2$.

$$\begin{cases} \mathcal{E}_{r,\phi}^{-\frac{1}{2}} \ = \ \mathcal{E}_{r,\phi}^{\frac{1}{2}} \ , \quad & \nu \neq~0\\ \mathcal{E}_{r,\phi}^{-\frac{1}{2}} \ = \ -\mathcal{E}_{r,\phi}^{\frac{1}{2}} \ , \quad & \nu = 0 \end{cases}$$

The decision between which form to use is made with the same choice of $\nu$ that is used for the inner grid discretization and the outer boundary conditions.

Funding

Air Force Office of Scientific Research (FA9550-19-1-0032).

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. A. E. Siegman, “Unstable optical resonators,” Appl. Opt. 13(2), 353–367 (1974). [CrossRef]  

2. J. Hu and C. R. Menyuk, “Understanding leaky modes: slab waveguide revisited,” Adv. Opt. Photonics 1(1), 58–106 (2009). [CrossRef]  

3. Q. Wang and W.-M. Zhao, “A comprehensive review of lossy mode resonance-based fiber optic sensors,” Opt. Lasers Eng. 100, 47–60 (2018). [CrossRef]  

4. J. Kim, B. Jang, J. Gargiulo, J. Bürger, J. Zhao, S. Upendar, T. Weiss, S. A. Maier, and M. A. Schmidt, “The optofluidic light cage - on-chip integrated spectroscopy using an antiresonance hollow core waveguide,” Anal. Chem. 93(2), 752–760 (2021). [CrossRef]  

5. D. E. Smalley, S. Jolly, G. E. Favalora, and M. G. Moebius, “Status of leaky mode holography,” Photonics 8(8), 292 (2021). [CrossRef]  

6. O. Peleg, Y. Plotnik, N. Moiseyev, O. Cohen, and M. Segev, “Self-trapped leaky waves and their interactions,” Phys. Rev. A 80(4), 041801 (2009). [CrossRef]  

7. M. Cassataro, D. Novoa, M. C. Günendi, N. N. Edavalath, M. H. Frosz, J. C. Travers, and P. S. Russell, “Generation of broadband mid-ir and uv light in gas-filled single-ring hollow-core pcf,” Opt. Express 25(7), 7637–7644 (2017). [CrossRef]  

8. P. Sprangle, E. Esarey, J. Krall, and G. Joyce, “Propagation and guiding of intense laser pulses in plasmas,” Phys. Rev. Lett. 69(15), 2200–2203 (1992). [CrossRef]  

9. T. R. Clark and H. M. Milchberg, “Optical mode structure of the plasma waveguide,” Phys. Rev. E 61(2), 1954–1965 (2000). [CrossRef]  

10. L. Feder, B. Miao, J. E. Shrock, A. Goffin, and H. M. Milchberg, “Self-waveguiding of relativistic laser pulses in neutral gas channels,” Phys. Rev. Res. 2(4), 043173 (2020). [CrossRef]  

11. B. Miao, L. Feder, J. Shrock, A. Goffin, and H. Milchberg, “Optical guiding in meter-scale plasma waveguides,” Phys. Rev. Lett. 125(7), 074801 (2020). [CrossRef]  

12. P. Panagiotopoulos, M. Kolesik, S. Tochitsky, and J. V. Moloney, “Generation of long homogeneous plasma channels with high power long-wave ir pulsed bessel beams,” Opt. Lett. 46(21), 5457–5460 (2021). [CrossRef]  

13. W. C. Chew and W. H. Weedon, “A 3D perfectly matched medium from modified Maxwell’s equations with stretched coordinates,” Microw. Opt. Technol. Lett. 7(13), 599–604 (1994). [CrossRef]  

14. S. Kim and J. Pasciak, “The computation of resonances in open systems using a perfectly matched layer,” Math. Comp. 78(267), 1375–1398 (2009). [CrossRef]  

15. J. C. Araujo-Cabarcas and C. Engström, “On spurious solutions in finite element approximations of resonances in open systems,” Comput. Math. with Appl. 74(10), 2385–2402 (2017). [CrossRef]  

16. G. Mur, “Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic-field equations,” IEEE Trans. Electromagn. Compat. EMC-23(4), 377–382 (1981). [CrossRef]  

17. P. T. Leung, S. Y. Liu, and K. Young, “Completeness and time-independent perturbation of the quasinormal modes of an absorptive and leaky cavity,” Phys. Rev. A 49(5), 3982–3989 (1994). [CrossRef]  

18. P. T. Leung, S. Y. Liu, and K. Young, “Completeness and orthogonality of quasinormal modes in leaky optical cavities,” Phys. Rev. A 49(4), 3057–3067 (1994). [CrossRef]  

19. M. Mansuripur, M. Kolesik, and P. Jakobsen, “Leaky modes of dielectric cavities,” in Spintronics IX, vol. 9931H.-J. Drouhin, J.-E. Wegrowe, and M. Razeghi, eds., International Society for Optics and Photonics (SPIE, 2016), pp. 14–33.

20. P. Lalanne, W. Yan, K. Vynck, C. Sauvan, and J.-P. Hugonin, “Light interaction with photonic and plasmonic resonances,” Laser Photonics Rev. 12(5), 1700113 (2018). [CrossRef]  

21. P. T. Kristensen, K. Herrmann, F. Intravaia, and K. Busch, “Modeling electromagnetic resonators using quasinormal modes,” Adv. Opt. Photonics 12(3), 612–708 (2020). [CrossRef]  

22. J. Heinz and M. Kolesik, “Transparent boundary conditions for discretized non-hermitian eigenvalue problems,” Int. J. Mod. Phys. C 32(10), 2150138 (2021). [CrossRef]  

23. A. Bahl, J. K. Wahlstrand, S. Zahedpour, H. M. Milchberg, and M. Kolesik, “Nonlinear optical polarization response and plasma generation in noble gases: Comparison of metastable-electronic-state-approach models to experiments,” Phys. Rev. A 96(4), 043867 (2017). [CrossRef]  

24. E. Polizzi, “Density-matrix-based algorithm for solving eigenvalue problems,” Phys. Rev. B 79(11), 115112 (2009). [CrossRef]  

25. C. R. Pollock and M. Lipson, Step-Index Circular Waveguides (Springer US, Boston, MA, 2003), pp. 73–97.

26. A. Hartung, J. Kobelke, A. Schwuchow, K. Wondraczek, J. Bierlich, J. Popp, T. Frosch, and M. A. Schmidt, “Double antiresonant hollow core fiber – guidance in the deep ultraviolet by modified tunneling leaky modes,” Opt. Express 22(16), 19131–19140 (2014). [CrossRef]  

27. M. Zeisberger and M. A. Schmidt, “Analytic model for the complex effective index of the leaky modes of tube-type anti-resonant hollow core fibers,” Sci. Rep. 7(1), 11761 (2017). [CrossRef]  

28. R. Olshansky, “Leaky modes in graded index optical fibers,” Appl. Opt. 15(11), 2773–2777 (1976). [CrossRef]  

29. P. Yeh, A. Yariv, and E. Marom, “Theory of bragg fiber,” J. Opt. Soc. Am. 68(9), 1196–1201 (1978). [CrossRef]  

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (4)
Fig. 1.
Fig. 1. Illustration of a waveguide cross-section with the domain divided into interior and exterior regions. The grey dash-dotted line represents the circular boundary at $r_B$, and a hypothetical three-point finite difference stencil is overlaid. The large dot represents a discrete boundary point, the small dots represent neighboring points, and the open circle represents a "ghost" point that must be expressed in terms of interior points.

View in Article | Download Full Size | PDF

Fig. 2.
Fig. 2. Calculated effective indices of the fundamental TE, TM, and HE resonances for a step-index fiber with varying radii. $\lambda = 800 \ \mathrm {nm}, \ n_{core} = 1.5, \ n_{clad} = 1.45$. Reference values calculated using the boundary matching approach [25].

View in Article | Download Full Size | PDF

Fig. 3.
Fig. 3. Calculated effective indices of the fundamental TE, TM, and HE resonances for a hollow core fiber of varying radii. $\lambda = 1.2 \ \mathrm {\mu m}, \ w = 0.7 \ \mathrm {\mu m}, \ n_{tube} = 1.45, \ n_{core} = n_{clad} = 1$. Reference values calculated using the model given by Zeisberger and Schmidt [27].

View in Article | Download Full Size | PDF

Fig. 4.
Fig. 4. Numerical studies of the finite difference implementation using the $\mathrm {HE_{11}}$ mode and a 20 µm tube waveguide. (a) Imaginary part of the effective index for a range of grid spacings. (b) Deviation of the imaginary part of the effective index for a range of boundary offsets.

View in Article | Download Full Size | PDF

Tables (2)

Tables Icon

Table 1. Convergence, R = 4 μ m

View Table | View all tables in this article
Tables Icon

Table 2. Convergence, R = 20 μ m

View Table | View all tables in this article

Equations (36)

Equations on this page are rendered with MathJax. Learn more.

∇ 2 E → ( r → ) + k 0 2 n ( r → ) 2 E → ( r → ) − ∇ ( ∇ ⋅ E → ( r → ) ) = 0 ,
E → ( r , ϕ , z ) = e − i β z E → ( r , ϕ ) = ∑ ν e i ( ν ϕ − β z ) E → ( ν ) ( r )
E ( r 0 , ϕ ) = 1 2 π ∫ ϕ ′ ∑ ν e i ν ( ϕ − ϕ ′ ) M ν ( r 0 ) M ν − 1 ( r B ) E ( r B , ϕ ′ ) d ϕ ′ ,
∇ 2 E → ( r → ) + k 0 2 n 0 2 E → ( r → ) = 0 .
r 2 ∂ 2 E z ( ν ) ∂ r 2 + r ∂ E z ( ν ) ∂ r + r 2 ( k 0 2 n 0 2 − β 2 − ν 2 r 2 ) E z ( ν ) = 0
E z ( ν ) ( r ) = A ν ( J ν ( κ r ) + i Y ν ( κ r ) ) = A ν H ν ( 1 ) ( κ r ) , H z ( ν ) ( r ) = B ν ( J ν ( κ r ) + i Y ν ( κ r ) ) = B ν H ν ( 1 ) ( κ r ) ,
E r = − i β κ 2 ( i ν ω μ β r H z ( ν ) + ∂ E z ( ν ) ∂ r ) E ϕ = − i β κ 2 ( i ν r E z ( ν ) − ω μ β ∂ H z ( ν ) ∂ r ) ,
E ( ν ) ( r ) = M ν ( r ) C ν , M ν − 1 ( r ) E ( ν ) ( r ) = C ν ,
E ( ν ) ( r ) = [ E r ( ν ) ( r ) E ϕ ( ν ) ( r ) ] , C ν = [ A ν B ν ] ,
E ( ν ) ( r 0 ) = M ν ( r 0 ) C ν = M ν ( r 0 ) ( M ν − 1 ( r B ) E i n t ( ν ) ( r B ) )
M ν ( r ) = [ − i β κ H ν ′ ( 1 ) ( κ r ) ω μ ν κ 2 r H ν ( 1 ) ( κ r ) β ν κ 2 r H ν ( 1 ) ( κ r ) i ω μ κ H ν ′ ( 1 ) ( κ r ) ]
M ν ( r 0 ) M ν − 1 ( r B ) = C ( r B , ν ) M ′ ( r 0 , r B , ν ) ,
M ′ ( r 0 , r B , ν ) = [ M d ′ ( r 0 , r B , ν ) M o ′ ( r 0 , r B , ν ) − M o ′ ( r 0 , r B , ν ) M d ′ ( r 0 , r B , ν ) ]
C ( r B , ν ) = 1 ( ν 2 κ 2 r B 2 H ν ( 1 ) ( κ r B ) 2 − H ν ′ ( 1 ) ( κ r B ) 2 ) ,
M d ′ ( r 0 , r B , ν ) = ν 2 κ 2 r B r 0 H ν ( 1 ) ( κ r B ) H ν ( 1 ) ( κ r 0 ) − H ν ′ ( 1 ) ( κ r B ) H ν ′ ( 1 ) ( κ r 0 ) ,
M o ′ ( r 0 , r B , ν ) = i ν κ r 0 H ν ′ ( 1 ) ( κ r B ) H ν ( 1 ) ( κ r 0 ) − i ν κ r B H ν ( 1 ) ( κ r B ) H ν ′ ( 1 ) ( κ r 0 ) .
∇ ⊥ 2 E + k 0 2 n 2 E = ∇ ⊥ ( ∇ ⊥ ⋅ E − 1 n 2 ∇ ⊥ ⋅ ( n 2 E ) )
∇ 2 E r − E r r 2 − 2 r 2 ∂ E ϕ ∂ ϕ + ∂ ∂ r [ 1 n 2 ∂ ∂ r ( n 2 E r ) ] − ∂ 2 E r ∂ r 2 + k 0 2 n 2 E r = 0
∇ 2 E ϕ − E ϕ r 2 + 2 r 2 ∂ E r ∂ ϕ + 1 r ∂ ∂ ϕ [ 1 n 2 ∂ ∂ r ( n 2 E r ) ] − 1 r ∂ 2 E r ∂ ϕ ∂ r + k 0 2 n 2 E ϕ = 0
1 r ∂ E r ∂ r − 2 i ν r 2 E ϕ + ∂ ∂ r [ 1 n 2 ∂ ∂ r ( n 2 E r ) ] + ( k 0 2 n 2 − 1 + ν 2 r 2 ) E r = β 2 E r
1 r ∂ E ϕ ∂ r + ∂ 2 E ϕ ∂ r 2 + 2 i ν r 2 E r + i ν r [ 1 n 2 ∂ ∂ r ( n 2 E r ) − ∂ E r ∂ r ] + ( k 0 2 n 2 − 1 + ν r 2 ) E ϕ = β 2 E ϕ
E = [ E r E ϕ ] , [ A r r A r ϕ A ϕ r A ϕ ϕ ] E = β 2 E
A r r E r = 1 r ∂ E r ∂ r + ∂ ∂ r [ 1 n 2 ∂ ∂ r ( n 2 E r ) ] + ( k 0 2 n 2 − 1 + ν r 2 ) E r
A ϕ ϕ E ϕ = 1 r ∂ E ϕ ∂ r + ∂ 2 E ϕ ∂ r 2 + ( k 0 2 n 2 − 1 + ν r 2 ) E ϕ
A r ϕ E ϕ = − 2 i ν r 2 E ϕ
A ϕ r E r = 2 i ν r 2 E r + i ν n 2 r ∂ n 2 ∂ r E r
β 2 = n e f f 2 k 0 2 .
r n = Δ r [ 1 2 + ( n − 1 ) ] , f n = f ( r n ) ,
E r N + 1 = C N , ν ( M d ′ N , N + 1 , ν   E r N + M o ′ N , N + 1 , ν   E ϕ N )
E ϕ N + 1 = C N , ν ( − M o ′ N , N + 1 , ν   E r N + M d ′ N , N + 1 , ν   E ϕ N )
1 Δ r [ ( 1 ϵ ∂ ϵ E r ∂ r ) i + 1 / 2 − ( 1 ϵ ∂ ϵ E r ∂ r ) i − 1 / 2 ] ,
1 Δ r 2 [ ( 1 ϵ ) i + 1 / 2 ( ϵ i + 1 E r i + 1 − ϵ i E r i ) − ( 1 ϵ ) i − 1 / 2 ( ϵ i E r i − ϵ i − 1 E r i − 1 ) ]
( 1 ϵ ) i + 1 / 2 ≈ 1 2 ϵ i + 1 + 1 2 ϵ i
1 2 Δ r 2 [ ϵ i + 1 + ϵ i ϵ i E r i + 1 + ϵ i + ϵ i − 1 ϵ i E r i − 1 − ( ϵ i + 1 + ϵ i ϵ i + 1 + ϵ i + ϵ i − 1 ϵ i − 1 ) E r i ]
{ E r , ϕ ( 0 ) = 0 , ν ≠   0 ∂ E r , ϕ ∂ r | r = 0 = 0 , ν = 0
{ E r , ϕ − 1 2   =   E r , ϕ 1 2   , ν ≠   0 E r , ϕ − 1 2   =   − E r , ϕ 1 2   , ν = 0
Select as filters
Select Topics Cancel
Top
       
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.
Privacy | Terms of Use