Search

Photonic structures with high-$Q$ resonances are essential for many practical applications, and they can be relatively easily realized by modifying ideal structures with bound states in the continuum (BICs). When an ideal photonic structure with a BIC is perturbed, the BIC may be destroyed (becomes a resonant state) or may continue to exist with a slightly different frequency and a slightly different wavevector (if appropriate). Some BICs are robust against certain structural perturbations, but most BICs are nonrobust. Recent studies suggest that a nonnegative integer $n$ can be defined for any generic nondegenerate BIC with respect to a properly defined set of structural perturbations. The integer $n$ is the minimum number of tunable parameters needed to preserve the BIC for perturbations arbitrarily chosen from the set. Robust and nonrobust BICs have $n=0$ and $n\ge 1$, respectively. A larger $n$ implies that the BIC is more difficult to find. If a structure is given by $m$ real parameters, the integer $n$ is the codimension of a geometric object formed by the parameter values at which the BIC exists in the $m$-dimensional parameter space. In this paper, we suggest a formula for $n$, give some justification for the general case, calculate $n$ for different types of BICs in two-dimensional structures with a single periodic direction, and illustrate the results by numerical examples. Our study improves the theoretical understanding on BICs and provides useful guidance to their practical applications.

Numerical mode matching (NMM) methods are widely used for analyzing wave propagation and scattering in structures that are piece-wise uniform along one spatial direction. For open structures that are unbounded in transverse directions (perpendicular to the uniform direction), the NMM methods use the perfectly matched layer (PML) technique to truncate the transverse variables. When incident waves are specified in homogeneous media surrounding the main structure, the total field is not always outgoing, and the NMM methods rely on reference solutions for each uniform segment. Existing NMM methods have difficulty handing gracing incident waves and special incident waves related to the onset of total internal reflection, and are not very efficient at computing reference solutions for non-plane incident waves. In this paper, a new NMM method is developed to overcome these limitations. A Robin-type boundary condition is proposed to ensure that non-propagating and non-decaying wave field components are not reflected by truncated PMLs. Exponential convergence of the PML solutions based on the hybrid Dirichlet-Robin boundary condition is established theoretically. A fast method is developed for computing reference solutions for cylindrical incident waves. The new NMM is implemented for two-dimensional structures and polarized electromagnetic waves. Numerical experiments are carried out to validate the new NMM method and to demonstrate its performance., 26 pages, 15 figures

Starting from Hadamard's method, we extend Babich's ansatz to the frequency-domain point-source (FDPS) Maxwell's equations in an inhomogeneous medium in the high-frequency regime. First, we develop...

We have developed a multiple level-set method for inverting magnetic data produced by weak induced magnetization only. The method is designed to deal with a specific class of 3D magnetic inverse problems in which the magnetic susceptibility is known and the objective of the inversion is to find the boundary or geometric shape of the causative bodies. We adopt the conceptual representation of the subsurface geologic structure by a set of magnetic bodies, each having a uniform magnetic susceptibility embedded in a nonmagnetic background. This representation enables us to reformulate the magnetic inverse problem into a domain inverse problem for those unknown domains defining the supports of the magnetic causative bodies. Because each body may take on a variety of shapes, and we may not know the number of bodies a priori either, we use multiple level-set functions to parameterize these domains so that the domain inverse problem can be further reduced to an optimization problem of multiple level-set functions. To efficiently compute gradients of the nonlinear functional arising from the multiple level-set formulation, we take advantage of the rapid decay of the magnetic kernels with distance to significantly speed up the matrix-vector multiplications in the minimization process. We apply the new method to the synthetic and field data sets and determine its effectiveness.

We have developed a level-set method for the inverse gravimetry problem of imaging salt structures with density contrast reversal. Under such a circumstance, a part of the salt structure contributes two completely opposite anomalies that counteract with each other, making it unobservable to the gravity data. As a consequence, this amplifies the inherent nonuniqueness of the inverse gravimetry problem so that it is much more challenging to recover the whole salt structure from the gravity data. To alleviate the severe nonuniqueness, it is reasonable to assume that the density contrast between the salt structure and the surrounding sedimentary host depends upon the depth only and is known a priori. Consequently, the original inverse gravity problem reduces to a domain inverse problem, where the supporting domain of the salt body becomes the only unknown. We have used a level-set function to parametrize the boundary of the salt body so that we reformulated the domain inverse problem into a nonlinear optimization problem for the level-set function, which was further solved for by a gradient descent method. Both 2D and 3D experiments on the SEG/EAGE salt model were carried out to demonstrate the effectiveness and efficiency of the new method. The algorithm was able to recover dipping flanks of the salt model, and it only took 40 min in a 2.5 GHz CPU to invert for a 3D model of 97,000 unknowns.

We propose a novel Babich-like ansatz consisting of an infinite series of dyadic coefficients (three-by-three matrices) and spherical Hankel functions for solving point-source Maxwell's equations in an inhomogeneous medium so as to produce the so-called dyadic Green's function. Using properties of spherical Hankel functions, we derive governing equations for the unknown asymptotics of the ansatz including the traveltime function and dyadic coefficients. By proposing matching conditions at the point source, we rigorously derive asymptotic behaviors of these geometrical-optics ingredients near the source so that their initial data at the source point are well-defined. To verify the feasibility of the proposed ansatz, we truncate the ansatz to keep only the first two terms, and we further develop partial differential equation--based Eulerian approaches to compute the resulting asymptotic solutions. Since the system of governing equations for each dyadic coefficient is strongly coupled, we introduce auxiliary variables to transform these strongly coupled systems into decoupled scalar equations. Furthermore, we develop high-order Lax--Friedrichs weighted essentially nonoscillatory schemes for computing these auxiliary variables so that the Green's function can be constructed. Numerical examples demonstrate that our new ansatz yields a uniform asymptotic solution in the region of space containing a point source but no other caustics.

In some applications, it is reasonable to assume that geodesics (rays) have a consistent orientation so that Maxwell's equations may be viewed as an evolution equation in one of the spatial directions. With such applications in mind, we propose a new Eulerian geometrical-optics method, dubbed the fast Huygens sweeping method, for computing Green's functions of Maxwell's equations in inhomogeneous media in the high-frequency regime and in the presence of caustics. The first novelty of the fast Huygens sweeping method is that a new dyadic-tensor-type geometrical-optics ansatz is proposed for Green's functions which is able to utilize some unique features of Maxwell's equations. The second novelty is that the Huygens--Kirchhoff secondary source principle is used to integrate many locally valid asymptotic solutions to yield a globally valid asymptotic solution so that caustics associated with the usual geometrical-optics ansatz can be treated automatically. The third novelty is that a butterfly algorithm is adapted to carry out the matrix-vector products induced by the Huygens--Kirchhoff integration in $O(N\log N)$ operations, where $N$ is the total number of mesh points, and the proportionality constant depends on the desired accuracy and is independent of the frequency parameter. To reduce the storage of the resulting traveltime and amplitude tables, we compress each table into a linear combination of tensor-product based multivariate Chebyshev polynomials so that the information of each table is encoded into a small number of Chebyshev coefficients. The new method enjoys the following desired features: (1) it precomputes a set of local traveltime and amplitude tables; (2) it automatically takes care of caustics; (3) it constructs Green's functions of Maxwell's equations for arbitrary frequencies and for many point sources; (4) for a specified number of points per wavelength it constructs each Green's function in nearly optimal complexity $O(N\log N)$ in terms of the total number of mesh points $N$, where the prefactor of the complexity depends only on the specified accuracy and is independent of the frequency parameter. Three-dimensional numerical experiments are presented to demonstrate the performance and accuracy of the new method.

We have developed a local level-set method for inverting 3D gravity-gradient data. To alleviate the inherent nonuniqueness of the inverse gradiometry problem, we assumed that a homogeneous density contrast distribution with the value of the density contrast specified a priori was supported on an unknown bounded domain [Formula: see text] so that we may convert the original inverse problem into a domain inverse problem. Because the unknown domain [Formula: see text] may take a variety of shapes, we parametrized the domain [Formula: see text] by a level-set function implicitly so that the domain inverse problem was reduced to a nonlinear optimization problem for the level-set function. Because the convergence of the level-set algorithm relied heavily on initializing the level-set function to enclose the gravity center of a source body, we applied a weighted [Formula: see text]-regularization method to locate such a gravity center so that the level-set function can be properly initialized. To rapidly compute the gradient of the nonlinear functional arising in the level-set formulation, we made use of the fact that the Laplacian kernel in the gravity force relation decayed rapidly off the diagonal so that matrix-vector multiplications for evaluating the gradient can be accelerated significantly. We conducted extensive numerical experiments to test the performance and effectiveness of the new method.

We propose a multiple level set method for inverting three-dimensional magnetic data induced by magnetization only. To alleviate inherent non-uniqueness of the inverse magnetic problem, we assume that the subsurface geological structure consists of several uniform magnetic mass distributions surrounded by homogeneous non-magnetic background such as soil, where each magnetic mass distribution has a known constant susceptibility and is supported on an unknown sub-domain. This assumption enables us to reformulate the original inverse magnetic problem into a domain inverse problem for those unknown domains defining the supports of those magnetic mass distributions. Since each uniform mass distribution may take a variety of shapes, we use multiple level-set functions to parameterize these domains so that the domain inverse problem can be further reduced to an optimization problem for multiple level-set functions. To compute rapidly gradients of the nonlinear functional arising in the multiple-level-set formulation, we utilize the fact that the kernel function in the field-susceptibility relation decays rapidly off the diagonal so that matrix-vector multiplications for evaluating the gradients can be speeded up significantly. Numerical experiments are carried out to illustrate the effectiveness of the new method.

Photonic crystal fibers (PCFs) with many air holes and complicated geometries can be difficult to analyze using conventional waveguide mode solvers such as the finite element method. Boundary integral equation (BIE) methods are suitable for PCFs, since they formulate eigenvalue problems only on the interfaces and are capable of computing leaky modes accurately. To improve the efficiency, it is desirable to have high-order BIE methods that calculate the minimum number of functions on the interfaces. Existing BIE methods calculate two or four functions on the interfaces, but high-order implementations are only available for those with four functions. In this paper, a new high-order BIE method is developed and it calculates two functions on the interfaces. Numerical results indicate that the new BIE method achieves exponential convergence and extremely high accuracy.

Full-vectorial beam propagation methods (FVBPMs) are widely used to model light waves propagating in high-index-contrast optical waveguides. We show that the paraxial FVBPM and wide-angle FVBPMs based on diagonal Pade approximants are analytically unstable for waveguides with complex modes. The instability cannot be removed by enlarging the computational domain, increasing the numerical resolution, or using perfectly matched layers, because the complex modes are highly confined around the waveguide core.

A new integral equation mode solver is developed for dielectric waveguides with sharp corners and high index-contrast. The method makes use of a Nystrom method for discretizing integral operators with logarithmic singularities, a graded mesh technique for handling corners, and computes boundary operators for domains of constant refractive indices. The method is illustrated by a number of numerical examples.

An accurate boundary integral equation method is developed for analyzing 2D photonic crystals where the cylinders in the unit cells have arbitrary shapes and corners. It first calculates the so-called Neumann-to-Dirichlet map for unit cells.

Conventional integral equation methods for diffraction gratings require lattice sum techniques to evaluate quasi-periodic Green's functions. The boundary integral equation Neumann-to-Dirichlet map (BIE-NtD) method in Wu and Lu [J. Opt. Soc. Am. A 26, 2444 (2009)], [J. Opt. Soc. Am. A 28, 1191 (2011)] is a recently developed integral equation method that avoids the quasi-periodic Green's functions and is relatively easy to implement. In this paper, we present a number of improvements for this method, including a revised formulation that is more stable numerically, and more accurate methods for computing tangential derivatives along material interfaces and for matching boundary conditions with the homogeneous top and bottom regions. Numerical examples indicate that the improved BIE-NtD map method achieves a high order of accuracy for in-plane and conical diffractions of dielectric gratings.