# Scalable numerical approach for the steady-state ab initio laser theory

###### Abstract

We present an efficient and flexible method for solving the non-linear lasing equations of the steady-state ab initio laser theory. Our strategy is to solve the underlying system of partial differential equations directly, without the need of setting up a parametrized basis of constant flux states. We validate this approach in one-dimensional as well as in cylindrical systems, and demonstrate its scalability to full-vector three-dimensional calculations in photonic-crystal slabs. Our method paves the way for efficient and accurate simulations of microlasers which were previously inaccessible.

###### pacs:

42.55.Sa, 42.55.Ah, 42.25.Bs## I Introduction

As lasers become increasingly complicated, especially in nanophotonic systems with wavelength-scale features Painter et al. (1999); Noda et al. (2000); Altug et al. (2006); Ellis et al. (2011), there has been a corresponding increase in the computational difficulty of solving for their nonlinear behavior, as described by the Maxwell–Bloch (MB) equations Haken (1985). To address this key challenge in the design and understanding of lasers, a highly efficient approach to finding the non-linear steady-state properties of complex laser systems has recently been introduced, known by the acronym SALT (steady-state ab initio laser theory) ^{1}^{1}1This is the same theory as the ab initio self consistent laser theory abbreviated as AISC, but the name has been changed.. In this paper, we present a technique to directly solve the SALT formulation Türeci et al. (2006, 2008); Ge et al. (2010) of the steady-state MB equations (using finite-difference frequency-domain (FDFD) Christ and Hartnagel (1987); Champagne et al. (2001) or finite-element methods (FEM) Chew et al. (2001)), and we demonstrate that, unlike previous approaches to the SALT equations Türeci et al. (2008); Ge et al. (2010), our technique scales to full three-dimensional (3D) low-symmetry geometries (such as photonic-crystal slabs Joannopoulos et al. (2001)).

The SALT equations (reviewed in Sec. II) simplify the general MB equations by removing the time dependence for steady-state modes, which allows SALT solvers to be potentially far more efficient than previous time-domain approaches Ziolkowski et al. (1995); Nagra and York (1998), while providing comparable accuracy Ge et al. (2008); Cerjan et al. (2012). However, all earlier approaches to SALT required the intermediate construction of a specialized constant-flux (CF) basis for the laser modes. While efficient and yielding numerous insights in highly symmetric geometries where it can be constructed semi-analytically, the CF basis becomes unwieldy and numerically expensive for complex low-symmetry laser geometries, especially in three dimensions. In our approach, we solve the SALT equations directly as a set of coupled nonlinear partial differential equations (PDEs), using a combination of Newton-Raphson Press et al. (2007), sparse-matrix solver Hénon et al. (2002), and nonlinear eigenproblem Guillaume (1999) algorithms in standard FDFD or FEM discretizations. In Sec. IV, we validate our solver against previous CF solutions for one-dimensional (1D) and cylindrical systems, while demonstrating that even in one dimension the CF basis rapidly becomes large and expensive as the system is brought farther above threshold. Furthermore, we show in Sec. III.5 that analytical outgoing-radiation boundary conditions, which are difficult to generalize to three dimensions Taflove and Hagness (2005), can be substituted by the standard PML (perfectly matched layer) method Taflove and Hagness (2005); Champagne et al. (2001); Shin and Fan (2012) which is equally effective at modeling open systems. We also demonstrate multi-mode laser solutions (Secs. IV.2 and IV.3), and reproduce the non-trivial avoided crossing interaction between lasing and non-lasing modes found in Ref. Liertzer et al. (2012).

We conclude in Sec. IV.3 with full 3D vectorial laser-mode solutions for a photonic-crystal slab microcavity Joannopoulos et al. (2001). The appendixes provide further details on the computational techniques we use in this paper, but in general any standard computational method in electromagnetism could be combined with our nonlinear solver algorithms. We believe that this computational approach provides a powerful tool to design and explore laser phenomena in the complex geometries accessible to modern nanofabrication, which were previously intractable for accurate modeling.

The Maxwell-Bloch (MB) equations provide the most basic formulation of semi-classical laser theory. The propagation of the electromagnetic field is given by the classical Maxwell equations and only the interaction of the field with the gain medium, represented by an ensemble of two-level atoms embedded in a cavity or background linear medium, is treated quantum mechanically. The MB equations are a set of time-dependent coupled nonlinear equations that are typically hard to solve analytically, except by using many approximations and idealizations. In the generic case of laser systems where such approximations are not valid, the MB equations have typically been solved using numerically expensive time-domain simulations Ziolkowski et al. (1995); Nagra and York (1998). For the case of steady-state lasing, as noted above, a much more efficient theory for calculating the multi-periodic solutions of the MB equations is the steady-state ab-initio lasing theory (SALT) Türeci et al. (2006); Ge et al. (2010); Türeci et al. (2009). This theory has proven to be a viable tool for describing laser systems ranging from random lasers Türeci et al. (2008); Stano and Jacquod (2013); Hisch et al. (2013) to coupled laser systems Liertzer et al. (2012) and photonic-crystal lasers Chua et al. (2011a). It makes no a priori assumptions about the geometry of the laser system, treats the open (non-Hermitian) character of the laser system exactly, and the non-linear hole-burning interactions between the laser modes to infinite order. More realistic and quantitative laser modeling typically requires treating a gain medium with three, four or more relevant atomic levels, but it has been shown that for the steady-state properties, under the same assumptions as SALT, the semiclassical equations can be reduced to an effective two-level (MB) system with renormalized parameters and solved with essentially the same efficiency as two-level SALT Ge (2010); Cerjan et al. (2012). SALT can also be used to describe quantum properties of lasers by combining the non-linear scattering matrix of SALT with input-output theory, leading specifically to a general formula for the linewidth of each mode in the non-linear steady-state Chong and Stone (2012).

For readers familiar with linear resonant cavities in photonics, which essentially trap light for a long time in a small volume, a laser can be semiclassically understood via the introduction of nonlinear gain (amplification) whose strength is determined by an input-energy “pump” Yariv (1989). As the pump strength is increased, one eventually reaches a “threshold” at which the gain balances the cavity loss and a steady-state real-frequency lasing (“active”) mode comes into existence. A key element is that the gain is nonlinear: increasing the laser-mode amplitude depletes the excited states of the gain medium (via a “hole-burning” term in the gain), and so at a given pump strength above threshold there is a self-consistent stable laser amplitude. At higher pump strengths, however, this picture is complicated by the introduction of additional lasing modes, which interact nonlinearly and whose individual gains and losses are balanced simultaneously by the SALT equations. Also, while a linear “resonant mode” technically refers to a pole in the Green’s function (or scattering matrix) at a complex frequency lying slightly below the real axis, a lasing mode can arise from any pole that is pushed up to the real axis by the gain, even poles that start out far from the real axis and do not resemble traditional resonant-cavity modes (for example, in random lasers Türeci et al. (2008)).

A strategy for efficiently solving the SALT equations was introduced in Türeci et al. (2008, 2009) and significantly extended in Ge et al. (2010). These existing methods can be viewed as a spectral integral-equation method Boyd (2001): they solve the nonlinear problem by first parametrizing each laser mode in terms of a specialized “spectral” basis, called the “constant-flux (CF) states”, that solve a linear non-Hermitian Maxwell eigenproblem parametrized by its (unknown) real lasing frequency. Because the frequency is required to be real outside the cavity, the photon flux outside the laser cavity is conserved, unlike the well-known quasi-bound states of the system, which are also purely outgoing, but do not conserve flux. This basis is defined so that at the lasing threshold for each mode, where the non-linear hole-burning interaction term is zero, one member of the basis set is the lasing solution. Hence, by construction, the basis expansion for the SALT solution above but near threshold converges rapidly even when the non-linear terms are taken into account, and the SALT equations reduce to finding a relatively small number of expansion coefficients for each mode. In highly symmetric geometries such as 1D or cylindrical systems with uniform pumping, the CF states can be found semi-analytically in terms of known solutions of the Helmholtz equation in each homogeneous region (e.g., in terms of sinusoid or Bessel functions), and such a basis will typically converge exponentially quickly Boyd (2001) to the SALT solutions. Furthermore, the CF basis can be used as a starting point for other analyses of laser systems, such as to identify the cause of mode suppression due to modal interactions Türeci et al. (2008); Ge et al. (2010) and exceptional points Liertzer et al. (2012); Brandstetter et al. (2014). However, the CF basis also has some disadvantages for complex geometries or for lasers operating far above threshold where the nonlinearities are strong and the convergence is not so rapid. In complex geometries where Helmholtz solutions are not known analytically, the CF basis itself must be found numerically by a generic discretization (e.g., FDFD or FEM) for many real frequencies (since the lasing frequency is not known a priori above threshold) and for multiple CF eigenvalues at each frequency in order to ensure convergence. The lack of separable solutions in low-symmetry two-dimensional (2D) and 3D geometries also increases the number of basis functions that are required (in contrast to cylindrical systems, for example, where the solutions can be solved one at a time). In three dimensions, where the discretization might have millions of points (e.g., on a grid), even storing a CF basis consisting of hundreds or thousands of modes becomes a challenge, not to mention the expense of computing this many 3D eigenfunctions numerically or of computing the resulting SALT equation terms. As a consequence, our approach in this paper is to abandon the construction of the intermediate CF basis and instead to directly discretize and solve the nonlinear SALT PDEs. This approach enables us to solve even low-symmetry 3D systems, and greatly enhances the power of the SALT approach for modeling and for the design of realistic laser structures.

## Ii Review of SALT

The origin of the SALT equations are the MB equations, which nonlinearly couple an ensemble of two-level atoms with transition frequency to the electric field Haken (1985); Sargent et al. (1977):

(1) | |||

(2) | |||

(3) |

Here, and are the positive-frequency components of the electric field and polarization, respectively. The coupling to the negative-frequency components is neglected in terms of a rotating wave approximation (RWA) which is both very useful for simplifying the equations and very accurate under general conditions. Note that at no point did we or will we assume the standard slowly-varying envelope approximation, which, if used, reduces the accuracy of the MB solutions. The population inversion of the medium is given by in the absence of lasing, which is roughly proportional to the external pumping rate and thus generally referred to as the pump strength. One of the useful features of SALT is that this pump strength can have an arbitrary spatial profile in addition to a varying global amplitude, such that one can represent different experimental pumping protocols by evolving along a “pump trajectory” which we parametrize here by , following Ref. Liertzer et al. (2012). Note that if there are gain atoms in unpumped regions of the laser, then the pump strength will be negative in these regions and thus the SALT equations will automatically take into account absorption due to unexcited gain atoms. and are the relaxation rates of the polarization and inversion, respectively. The linear cavity dielectric function is homogeneous outside the cavity region, and consequently a finite spatial domain can be used for the laser system with an outgoing boundary condition. We have assumed a scalar and dipole matrix element although in anisotropic gain media they can be generalized to tensors.

The attractive feature of SALT is that it provides access to the spatial profiles of the lasing modes as well as to the lasing frequencies of a multi-mode microlaser at very low computational costs. To achieve this high performance, SALT makes two essential assumptions. First, it assumes that for a fixed pump strength the electric field and polarization eventually reach a multi-periodic steady state,

(4) | |||||

(5) |

with unknown lasing modes and real lasing frequencies . Second, SALT makes the stationary inversion approximation (SIA), i.e, . In the single-mode regime the SIA is not necessary, as the average inversion in steady-state is exactly zero, but in the multimode regime the inversion is in general not stationary and only under certain conditions is . However, the development of SALT was specifically oriented towards describing novel solid state microlasers and the necessary conditions are typically satisfied for such lasers, as we discuss in the following.

If the laser is operating in the multimode regime, then the term in Eq. (3) above will drive the inversion at all beat frequencies of active modes, which is of order , the free spectral range of the laser. In addition, the polarization can respond at the rate and could additionally drive time variation in the inversion. However, if the condition holds, then the inversion is being driven non-resonantly and responds quite weakly, except to the dc part of the drive which represents static gain saturation. The effects of the residual four-wave mixing can be included perturbatively if desired, as was done in Ref. Ge et al. (2008), but are neglected in standard SALT. The condition is satisfied in essentially all solid state lasers due to strong dephasing, but the condition depends on the linear dimensions and geometry of the laser cavity and is typically not satisfied for macro scale tabletop lasers. However for a linear cavity it typically would be satisfied for and hence the SIA tends to be a good approximation for multimode lasing in micro lasers. This general argument was made by Fu and Haken Fu and Haken (1991) in 1991 and was applied to Fabry-Perot lasers, for which they provided a stability proof for the multimode state under these conditions. These assumptions leading to the SIA allow the derivation of the much more general SALT equations, which were then tested extensively in comparison to full FDTD simulations for many multimode lasing structures in Refs. Ge et al. (2008); Cerjan et al. (2012); Liertzer et al. (2012). A general linear stability analysis in the SALT framework is challenging due to the necessity of testing stability against all possible spatial fluctuations, something not ever done in standard analyses, where the spatial degrees of freedom are frozen. However, work in this direction is in progress and partial results have been obtained. A condition relating to the stability of multimode solutions is discussed immediately below.

For completeness we note that this analysis of the validity of the SIA differs from the well-established classification of lasers into categories denoted class A, B, and C, depending on whether two, one, or zero of the fields can be adiabatically eliminated, meaning that the rapidly responding field instantaneously follows the slowly varying field(s). By far the most important case is class B, in which adiabatically follows (even in the transient dynamics), and the three MB equations are reduced to two equations for the field and the inversion. The condition for class B is expressed by the inequality , where is the cavity decay rate. This condition is neither necessary nor sufficient for the validity of the SIA in the multimode regime.

The class B condition is not sufficient, as is well-known, because once two or more modes lase, the beat frequency can drive complex and even chaotic dynamics of the inversion and field. One needs the further condition as just noted. This laser classification was introduced before the advent of the micro laser, for which this inequality holds, and hence it was assumed that multimode lasing would never be stable for class B. However, if the SIA condition holds then we do not need full adiabatic elimination of . The SIA and SALT can still describe the multimode steady-state which is eventually reached. The condition (“good cavity” limit) is not necessary to have a stable multimode solution. The magnitude of only affects the steady-state in terms of its stability to fluctuations. A noise driven fluctuation will oscillate as it decays at the relaxation frequency, ; if the beat frequency, , then the multimode interference can drive the fluctuations resonantly and destabilize the multimode solution. This mechanism was analyzed carefully in a number of works on the approach to chaos in lasers with injection Arecchi et al. (1984); Oppo et al. (1986) or multiple modes. This yields a third implicit condition on the validity of the SIA and SALT, i.e., . If one has a good cavity with then this is easily satisfied; but if one has a “bad cavity” with (which can be achieved) then the condition can still be satisfied if is sufficiently small. Thus we do expect SIA to hold in the multimode regime even for bad cavity lasers which are not standard class B, as long as this inequality holds, and SALT should describe multimode lasing in the bad-cavity limit. Comparisons of SALT with FDTD for bad cavities confirm this expectation, as well as recent work by Pillay et al. which uses SALT to compute the laser linewidth in the bad cavity limit Pillay et al. (2014).

Using these well-motivated approximations, Eq. (1) can then be written for each lasing mode as

(6) |

where the two-level active gain material is described by the non-linear susceptibility . Here, , is the Lorentzian gain curve, where

(7) |

and the population inversion. The latter contains the spatial hole-burning term that nonlinearly couples all lasing modes,

(8) |

where the are in their natural unit .

The non-linear SALT equations, Eq. (6), for the electric field of the lasing modes, , and for the associated lasing frequencies can be conceived of as the limit of an amplifying scattering process in which the input goes to zero, corresponding to purely outgoing solutions with real frequency or, equivalently, to a pole in the relevant scattering matrix on the real axis. Until the external pump is strong enough for the gain to balance the loss there will be no solution of this type, i.e., . However, when increasing the pump strength, non-trivial solutions appear at a sequence of thresholds and at different frequencies . The non-linear interaction between these solutions is through the spatial hole-burning and depletion of the gain medium, Eq.(8): each lasing mode extracts energy from the pump in a space-dependent manner which in general makes it more difficult for subsequent modes to reach threshold, and also effectively changes the index of refraction of the gain medium.

As already noted, Eq. (6) has been solved in 1D and 2D geometries, where either the electric or the magnetic field can be treated as scalar, for diverse systems such as random, microdisk or photonic crystal lasers using an algorithm based on expansion of the solutions in the CF basis Ge et al. (2010). In the most recent and most efficient formulation, the linear non-Hermitian eigenvalue problem,

(9) |

is used to define the optimal set of threshold CF states and eigenvalues .

The function adapts the basis to the spatial pump profile of the experiment of interest and is nonzero only inside the gain medium. The form a complete basis and satisfy a biorthogonality relation at any frequency . Equation (6) is solved by projecting the lasing modes into the CF basis. The resulting non-linear eigenvalue equation can only be satisfied at discrete frequencies which hence determine the lasing frequencies, . In principle one does not need to pre-calculate and store the CF basis at different real values of but it is numerically favorable to do so in general. However, the wider the Lorentzian gain curve, Eq. (7), is compared to the free spectral range, the more memory intensive the storage of the CF basis becomes, which makes calculations problematic in two and three dimensions. Moreover, if the pump profile is fixed and only its amplitude is varied experimentally, then CF states need only be calculated for various values, but if the pump profile also varies along a pump trajectory then one has to calculate new CF states also for many values of Liertzer et al. (2012). For a limited set of highly symmetric cavities, including piecewise-homogeneous 1D slabs and uniform cylinders, the solution of Eq. (9) is known semi-analytically at any . However, for all other geometries, Eq. (9) must be solved numerically for all relevant needed to build a basis. Consequently, for a fully-vectorial treatment of SALT in arbitrary cavities, CF bases cannot be used without significant computational costs. Our direct solution method eliminates the computation and storage of CF bases and scales easily to 3D geometries.

## Iii Solution method

### iii.1 Overview

The basic idea of our new solution method to obtain the lasing modes in the SALT is as follows: We discretize Eq. (6), using standard discretization techniques like FEM (see Appendix B) or FDFD (see Appendix C), and iteratively solve for the lasing modes and their frequencies at successively increasing values of the pump parameter . This nonlinear coupled problem is most conveniently solved by using the Newton-Raphson method. For initial guesses, we use the modes at threshold when we are close above threshold, and the modes at the previous pump step when we are far above threshold. In order to find the first threshold and the corresponding solution, Eq. (6) is initially solved for as an eigenvalue problem (EVP). The solutions are the resonances or quasi bound states of the passive cavity, corresponding to the poles of the passive scattering matrix ( matrix) Ge et al. (2010) with frequencies lying in the negative imaginary half plane (note that we will label all quantities below threshold with overbars throughout the paper). While increasing the pump , Eq. (6) is solved without the nonlinearity in Eq. (8) and the nonlasing modes near the gain frequency are tracked until the first reaches the real axis and turns the corresponding mode into an active lasing mode, . Once we have crossed the first threshold, we use the solutions for and of the eigenvalue problem at threshold as a first guess for the solution of and in the non-linear Newton solver slightly above threshold. The latter already includes the non-linearity which, once the Newton solver has converged, we treat as a fixed function like to examine the remaining non-lasing modes at the current pump strength . This has to be done in order to verify if further modes cross the lasing threshold. For the non-lasing modes, Eq. (6) is thus only nonlinear in and linear in , such that this problem can be cast into a nonlinear EVP Guillaume (1999). The procedure of increasing the pump is now continued by tracking the lasing mode solving the nonlinear coupled SALT system, while the non-lasing modes are evaluated from the corresponding nonlinear EVP until a second mode reaches threshold. At this point the number of lasing modes is increased by 1 and the procedure continues with two and more lasing modes in essentially the same way.

To illustrate this approach in more detail, we apply it to the simple one-dimensional edge-emitting laser shown in Fig. 1(a) which already captures all the main features. We pump the 1D slab cavity uniformly along its length with a pump strength which, above the first threshold, leads to emission to the right. Starting with , where the SALT system reduces to a simple resonance problem, we increase and observe that the resonance poles move upwards in the complex plane; see Fig. 1(b) where the starting point is marked by circles and the pump value at the first threshold, , is marked by triangles. Below this first threshold no mode is lasing, such that the non-linear spatial hole-burning term is zero, resulting in the following PDE for all non-lasing modes,

(10) |

which is linear with respect to , but into which the resonance values enter non-linearly. Starting at the first threshold, the terms and of the first lasing mode enter the spatial hole-burning denominator in Eq. (8) (where ), resulting in the following equation for the first lasing mode and its wavenumber ,

(11) |

which is now nonlinear with respect to both and . When continuing to increase the pump, the frequencies corresponding to the active modes are forced to stick to the real axis, while the eigenvalues associated to all other inactive modes continue moving upwards, see Fig. 1(b). To detect the activation of further modes, the inactive modes have to be recalculated again, however, this time by additionally taking into account the spatial hole burning contribution of the currently lasing active mode at a given pump strength . For this, we insert the currently active mode into the denominator of Eq. (11) which turns the above nonlinear problem into another nonlinear (in ) eigenvalue problem,

(12) |

which, however, has the same structure as Eq. (10). As soon as the imaginary part of another eigenvalue reaches the real axis, a new laser mode becomes active which increases the size of the nonlinear problem by 1. For even higher pump strength and a larger number of lasing modes this procedure continues accordingly. Note also, that the case when a mode shuts down during the pumping process can be incorporated without major effort.

To summarize, the solution of the SALT equations reduces essentially to computing the full nonlinear (in and ) system of PDEs through a Newton-Raphson method and the computation of an EVP which is linear in but which still remains nonlinear in . Details of how to obtain the active or lasing solutions of the Newton problem as well as the inactive or non-lasing solutions through the nonlinear EVP are provided in the following two sections.

### iii.2 Lasing modes

For modes that are lasing, Eq. (6) is nonlinear in the unknowns . As theses modes are all coupled together through the spatial hole-burning interaction, they must be solved simultaneously. In general, such systems of nonlinear equations can be written in the form

(13) |

where the vector of equations is an analytic nonlinear function of the unknown solution vector which again gathers all unknowns . This nonlinear problem can generally be solved by using the Newton-Raphson method Press et al. (2007). The basic idea is that for a guess for the solution , one can write

(14) |

where is the Jacobian matrix of partial derivatives of with respect to . A solution can usually be obtained by iterating Eq. (14) using only the linear terms. This iterative algorithm converges “quadratically” (squaring the errors on each step Press et al. (2007)) if is small. Further, we use an analytic evaluation for the Jacobian from Eq. (6), as described in Appendix A, and do not need to compute it using numerical differentiation schemes. Since is then a sparse matrix each iteration can exploit fast algorithms for sparse linear equations Trefethen and Bau (1997); Davis (2006).

To solve Eq. (13) on a discrete level, we project the complex fields of each lasing mode onto a discrete -component basis as explained in Appendixes B (for a FEM approach) and C (for an FDFD approach). Unlike the CF basis, we use a localized basis generated once from a grid or mesh. This is the key to producing sparse matrices and hence makes the method scalable to the larger bases required in two and three dimensions. The discretizations on such a basis turn the fields into complex coefficient vectors , while is required to be purely real. Because the SALT equations are not differentiable in the complex fields (due to the complex conjugation), we split our unknown coefficient vectors (and the vector function accordingly) into their real and imaginary parts. The discretized version of then consists of real unknowns (fields and frequencies). However, we only obtain real-valued equations from . The underspecification comes from the fact that the hole-burning term happens to be invariant under global phase rotations . In addition to the problem of underspecification, there is also a problem of stability: for lasing modes slightly above threshold, the amplitude is nearly zero, which would result in problems distinguishing between the solution we want and the trivial solution . We resolve both issues by normalizing the amplitude and fixing the phase of all lasing modes while keeping track of their amplitudes using a separate variable. This procedure results in both the number of real unknowns and the number of real equations being . Further details for our method for lasing modes are given in Appendix A.

Note that for the Newton-Raphson iteration to be scalable to higher dimensions and to high-resolution meshes, it is also important to use a scalable solver (in our case, the sparse direct solver Davis (2006) PaStiX Hénon et al. (2002) was called from the PETSc library Balay et al. (1997) because the Jacobian is sparse). For very large-scale 3D systems, it may become necessary to use iterative linear solvers Trefethen and Bau (1997) for each Newton step instead, in which case it is important to select certain PML formulations Shin and Fan (2012).

### iii.3 Non-lasing modes

In order to find the first pump threshold and the corresponding lasing solution as well as to verify when a new mode activates, the non-lasing modes have to be monitored while changing the pump. These non-lasing modes are defined as complex-frequency solutions to Eq. (6) that do not enter into the nonlinear hole burning term in , see Eq. (8). Due to causality constraints, the complex eigenvalues associated with non-lasing modes, always feature , and usually approach the real axis as is increased (interesting exceptions are discussed in Sec. IV.2). When all lasing modes have been determined for a particular , the function is known and can be treated as a fixed function like , see Eq. (12). As outlined in Sec. III.1, this reduces Eq. (6) to a non-Hermitian, nonlinear eigenvalue problem (NEVP) which is linear in the eigenvectors , but nonlinear in the complex eigenvalues .

For situations where we are only interested in the behavior of a few lasing modes in a small range of the pump parameter , Newton’s method is still a convenient approach to determine the non-lasing modes and, in fact, the only viable method in terms of computational cost for high resolution 2D or 3D computations. In this case, we typically use standard EVP algorithms to solve Eq. (6) first for (which is usually either linear or quadratic in , depending on the method for implementing the outgoing radiation condition). This provides us all the modes of interest which we then track to threshold with Newton’s method as is increased. As in Sec. III.2, convergence is “quadratic”, but, unlike for the lasing modes, Eq. (14) can be used with complex unknowns and equations since Eq. (6) is differentiable in all unknowns once is fixed. The downside of Newton’s method is that, in the absence of a good initial guess, it can be very unpredictable and slow to converge. Such a situation arises, e.g., when the modes that can lase are not known a priori as in the case where a large number of near-threshold modes are clustered together, all with frequencies close to the gain center . In this instance, a more general and comprehensive method for evaluating the non-lasing modes is required.

Such more general techniques exist in terms of NEVP solvers Guillaume (1999). One conceptually simple method for our problem is to divide Eqs. (6) and 12) by , turning the rational EVP into a cubic EVP which can then be linearized at the expense of making the problem three times as large and possibly also very ill-conditioned. Other, more sophisticated solution methods include “trimmed” linearization Su and Bai (2011), Newton Peters and Wilkinson (1979), Jacobi-Davidson Sleijpen et al. (1996), rational Krylov Ruhe (2003), and nonlinear Arnoldi Voss (2004). Independently of the chosen solution strategy, we can take into account that only modes which have a spectral overlap with the gain curve near its center frequency are expected to be candidates for active laser modes. In addition, the Lorentzian gain curve of width produces a singularity in the NEVP at which may result in spurious numerical solutions. Combining these observations, we restrict our attention to those eigenvalues that are in the following cropped subpart of the complex plane: . A suitable method that allows us to conveniently include such auxiliary restrictions is the contour integral method presented recently in Asakura et al. (2010); Beyn (2012) and reviewed in Appendix D. There, the search for eigenvalues is restricted to a region within a smooth contour such as a circle or an ellipse; see Fig. 2(a). By using the residue theorem, all poles of the inverse of the differential operator, which are equivalent to the eigenvalues of the same operator, are obtained within the specified contour. This feature is not only useful for employing this method as a stand-alone solver for non-lasing modes, but also as a complementary tool to check if, in addition to the limited set of non-lasing modes that are tracked with a Newton solver, no new modes have entered the region of interest within the chosen contour.

### iii.4 Alternative strategy for a single pump

Similar to the single-pole approximation in the CF-expansion method Ge et al. (2010), it is possible to speed up the calculations of the direct solver when the intensity of the laser is only desired at or starting from a specific pump strength . For this, we first solve the SALT equations, Eq. (6), only at this desired pump strength by neglecting any spatial hole-burning interactions in . If happens to be above the first threshold, the corresponding NEVP will yield complex frequencies that partly lie in the non physical region above the real axis in the complex plane. This is shown in Fig. 2(a), again for the simple 1D edge-emitting laser introduced above. Next, the most non-physical mode, i.e., the one which has the eigenvalue with the highest imaginary part, is selected and the corresponding solution vector as well as the real part of the corresponding eigenvalue are used as initial guesses in the nonlinear SALT solver. After the nonlinear iteration converges, the corresponding solution is then included in the spatial hole burning term, which effectively reduces the pump within the system and pulls down all inactive modes in the complex plane; see Fig. 2(b). If some of the remaining inactive modes are still located above the real axis, this procedure is repeated by increasing the number of active lasing modes until all modes lie on or below the real axis. The latter are then the true lasing modes of the SALT at the desired pump strength ; see Fig. 2(c). Hence, as long as the nonlinear solver manages to converge, a solution to the SALT can be obtained rather quickly.

### iii.5 Outgoing radiation condition

For numerical computations, the outgoing radiation condition must be implemented within a truncated, finite domain. In one dimension, the radiation condition can be expressed exactly Tolstikhin et al. (1998). This also allows us to shift the boundary of the domain right to the border of the cavity, which decreases the computational cost. This method is, however, not easily generalizable to two and three dimensions Taflove and Hagness (2005).

An efficient and robust alternative is to use the standard perfectly matched layer (PML) technique Berenger (1994); Bramble and Pasciak (2008) in which an artificial material is placed at the boundaries. The material has a certain complex permittivity and permeability such that it is absorbing and analytically reflectionless. In one dimension, the PML technique can be tested against an exact outgoing boundary condition, and the two methods yield results that are nearly indistinguishable, as shown in Fig. 3. Also shown in Fig. 3 is a comparison with conventional methods of solving the MB equations using finite difference time domain (FDTD) simulations demonstrating the validity of the stationary inversion approximation used in the derivation of the SALT equations. Both the quantitative agreement between SALT and FDTD solutions as well as the former’s substantial numerical efficiency over the latter have been previously documented Ge et al. (2008); Cerjan et al. (2012). Of course, the precise computation times depend on many factors, including hardware details, parameter choices in the algorithms, and software implementation quality, but the magnitude of the difference here makes it unlikely that any FDTD implementation could be competitive with the SALT approach.

## Iv Assessment and application of the solution method

In this section we will validate our solution strategy against the traditional method based on CF states and we will show first results for prototypical laser cavities.

### iv.1 1D slab laser as test case

We demonstrate here the accuracy of the presented direct solver method by studying in more detail the 1D edge-emitting slab laser introduced in Sec. III.1. One of the advantages of the direct solver, as compared to the CF method, is

the accuracy of its solutions far above the threshold. In this regime the CF basis becomes a poorer match for the lasing modes and, as explained in Sec. I, a large number of basis functions is required for convergence compared to near threshold. This is especially relevant for low- (short-lifetime) laser resonators such as random lasers or cavities featuring gain-induced states, as considered, e.g., in Ge et al. (2011) . In Fig. 4 the intensity of such a low- cavity is plotted with respect to an overall pump strength for a constant spatial pump profile. The figure contains both the results of the direct and of the CF state solver. For the latter the solution for different numbers of CF states are depicted, demonstrating that for a larger basis the solution converges towards the solution of the direct solver. Our solution method thus leads to an accuracy far above threshold which can only be achieved by the traditional approach with a considerably large number of CF states.

### iv.2 Non-uniform pump and avoided resonance crossings

In the following we consider an example for a laser for which the overall spatial profile of the applied pump, evolves non uniformly as a function of the pump parameter . As recently pointed out in Liertzer et al. (2012) such a spatially varying pump function can strongly influence the laser output in a counterintuitive way, an effect which has meanwhile also been verified experimentally Brandstetter et al. (2014). The system we consider to realize such a behavior consists of two coupled one-dimensional ridge cavities (see inset Fig. 5) which feature strong loss in the absence of pump. The pump function is defined as follows: For values of the pump parameter between and , only the left cavity of the system is pumped uniformly with an amplitude that is linearly increasing from zero (at ) to a value where the laser is close above threshold (at ). For between and the pump in the left cavity is kept constant (at the value for ), while the pump in the right cavity is linearly increased from zero (at ) to the same pump strength as in the left cavity (at ). Since the overall pump strength in the cavity steadily increases, one would expect that also the overall intensity of the laser should increase along this “pump trajectory” from to . Instead, the laser displays an intermittent shut down within a whole interval of around , as shown in Fig. 5.

In Ref. [23] this shutdown as obtained with SALT has been quantitatively verified against a traditional FDTD method to show that the solutions are stable and not an artifact of SALT. Furthermore, the shutdown was attributed to the occurrence of an exceptional point, corresponding to a non-Hermitian degeneracy in the TCF eigenvalues [see Eq. (9)] when parametrized over both the outside frequency and the pump parameter . In the direct solver, there no longer exists such a two-dimensional parameter space since the frequency can no longer be freely adjusted outside the cavity. Instead, the frequency is already obtained simultaneously with the corresponding lasing mode. We can thus expect that the poles associated with the (non)lasing modes reflect, in some form, their vicinity to the exceptional point through a nontrivial behavior along this pump trajectory. Indeed, our calculations show that the intermittent laser shut down is realized in terms of an avoided crossing between a lasing pole and a non-lasing pole in the complex plane (see Fig. 6). Here, the solid lines represent the solutions of the full SALT while the dashed lines show the movement of the complex eigenvalues when spatial hole burning is neglected.

In fact, we observe two avoided crossings in this plot. The first one occurs in the range between (marked as circles in Fig. 6) and (marked as squares). In this case the poles associated with the blue and the red mode first attract each other and then undergo an avoided crossing which pushes the red mode towards and, ultimately, beyond the real axis, i.e., the lasing threshold.

The second avoided crossing occurs in the interval between and , where we observe that the blue pole moves towards the real axis and interacts with the red pole such as to pull it below the real axis, corresponding to switching this mode off. In a corresponding experiment Brandstetter et al. (2014), only the second pole interaction can be directly observed in terms of an intermittent laser shut down, followed by a re-emergence of the laser modes at slightly detuned lasing frequencies.

Figure 6 also illustrates a crucial point touched on earlier: If one neglects the non-linear spatial hole burning interaction (dashed lines) one obtains poles in the upper half of the complex plane which violate the causality principle for the dielectric response. Including spatial hole burning (solid lines), keeps all poles below or on the real axis, as required [see Fig. 6(a)]. Note, that one also observes how the hole-burning interaction influences the movement of the non-lasing modes in terms of a delayed turn on of the blue mode [see the line between the two triangles in Fig. 6(a)].

### iv.3 Scalability to full-vector 2D and 3D calculations

In this section we briefly explore the applicability of our solution strategy to 2D and 3D setups by considering the following prototypical examples: In the 2D case we investigate a circular dielectric resonator and in the 3D case a photonic-crystal slab.

In the former situation we study a circular disk with uniform index, which is routinely used in the experiment due to its long-lived resonances associated with “whispering gallery modes” Vahala (2003). For this system we study lasing based on TM polarized modes and compare the Newton method presented here (based on FDFD) with the previously developed CF-state method Türeci et al. (2006); Ge (2010); Ge et al. (2010). Due to the azimuthal symmetry, the resonant TM modes Türeci et al. (2006); Dettmann et al. (2009) are exact solutions of the Bessel equation characterized by an azimuthal phase (with being an integer angular-momentum quantum number) and subject to outgoing boundary conditions. Due to the circular symmetry, each of the modes with a given value of comes with a degenerate partner mode, characterized by the quantum number . In the presence of the lasing nonlinearities, a preferred superposition will typically be selected as the stable solution, e.g., the circulating modes , rather than the and standing waves. The determination of this stable solution in a degenerate lasing cavity is a complex problem that we plan to address in future work. For validation and demonstration purposes in this paper, we simply select a priori a single solution from each degenerate pair by imposing corresponding symmetry boundary conditions. In the case of the circular cavity, we choose the circulating modes with a phase for comparison with the CF solutions. We obtain these by solving for both the sine and cosine modes (using the appropriate boundary conditions at the and symmetry planes) and by combining them to construct the exponentially circulating mode.

Under these premises, we find that for uniform pump the first mode turns on at and increases linearly in intensity, as seen in Fig. 7. The second mode turns on at about twice the pump strength as the first threshold. As the intensity of the second mode increases, we observe a reduction and ultimately a complete suppression of the first mode intensity. This mode competition can be attributed to the following two effects: The two modes have a significant spatial overlap, such that they compete for the same gain through non-linear spatial hole burning which is fully incorporated in SALT. In addition, as being spectrally closer to the peak of the gain curve , the second mode can profit more strongly from the gain in the disk than the first mode. As a result, the second mode prevails against the first mode in this non-linear competition. This behavior of interaction-induced mode switching is general and can be found in other laser configurations and nonlinear media as well Ge et al. (2014). In Fig. 7 we show that this behavior is faithfully reproduced with our approach, not only in terms of the modal intensities as a function of the applied pump (see top panel), but also in terms of the corresponding lasing modes which mirror those obtained with the CF-state technique very accurately (see bottom panel).

The second example we consider is a photonic crystal slab with a “defect” (see inset Fig. 8) engineered to efficiently trap a mode Lin et al. (2001). The photonic crystal is formed in a dielectric slab by holes which are arranged in a hexagonal lattice and the defect is created by decreasing the radius of seven of the holes in the center. In our study, we focus on a TE-like lasing mode, situated at the defect (spatially) and in the bandgap of the lattice (spectrally). To select one of the degenerate standing-wave solutions, we impose even and odd symmetry at and , respectively, as well as an even symmetry at . Staying close to a potential experimental realization, we choose the pump profile to be uniform inside the slab material’s defect region but zero outside and in the air holes. Increasing the overall amplitude of this pump profile, we find the lasing behavior shown in Fig. 8 (main panel). This calculation was performed with 16 nodes (using one CPU per node) of the Kraken Cray XT5 at the University of Tennessee. With pixels (the mirror conditions effectively halve these), the total wall-clock time for the computation, from passive resonance at to lasing above threshold at , was 5.9 min. Pump steps of were taken, with three to four Newton iterations per pump value.

## V Concluding Remarks

In this paper, we have presented an algorithm for solving the SALT equations which describe the steady-state lasing modes and frequencies of lasers with a free spectral range and a dephasing rate that are both large as compared to the population decay rate and the relaxation oscillation frequency. These conditions are typically satisfied by microlasers with a linear dimension that does not exceed a few hundred wavelengths. Our solution strategy proceeds by a direct discretization using standard methods as FEM or FDFD, without the need for an intermediate CF basis. The resulting increase in efficiency lets our approach scale to complex 2D and 3D lasing structures, which paves the way for future work in a number of directions.

First, it is now possible to study lasing in much more complex geometries than could previously be readily simulated, offering the possibility of discovering geometries that induce unexpected new lasing phenomena. Going one step further, future computations could search a huge space of lasing structures via large-scale optimization (“inverse design”), which has already been applied to the design of linear microcavities Frei et al. (2008); Kao and Santosa (2008); Lu and Vuckovic (2010). Since our approach is only more expensive than the solution of linear cavity modes by a small constant factor (e.g., the number of modes and the number of Newton iterations) it will be the ideal tool for this purpose. More complicated gain profiles, lineshapes, and other material properties can easily be incorporated into our approach as well. SALT can, e.g., be coupled to a diffusion equation in order to model the migration of excited atoms in molecular-gas lasers Chua et al. (2011b); Cerjan et al. (2014). Based on the mathematical relation of the multimode lasing equations to incoherent vector solitons (Appendix E), we believe that numerical methods commonly used in soliton theory can also be adopted to efficiently solve the multimode SALT equations. Another intriguing direction of research is the development of a more systematic approach to modeling lasers with degenerate linear modes, which requires a technique to evaluate the stability of the solution and evolve an unstable mode to a stable mode. Finally, many refinements are possible to the numerical methods, such as efficient iterative solvers and preconditioners for the Newton iterations of the lasing modes or criteria to alternate between systematic contour-integral evaluation and simpler Newton-inverse tracking of the non-lasing modes. In this sense our approach has more in common with standard sparse discretization methods used to solve other nonlinear PDEs than the CF-basis approach (which is specialized to the SALT problem) and thus opens the door for more outside researchers and numerical specialists to study lasing problems.

###### Acknowledgements.

The authors would like to thank the following colleagues for very fruitful discussions: S. Burkhardt, D. Krimer, X. Liang, Z. Musslimani, L. Nannen, H. Reid, and J. Schöberl. S.E., M.L., K.G.M., J.M.M., and S.R. acknowledge financial support by the Vienna Science and Technology Fund (WWTF) through Project No. MA09-030, by the Austrian Science Fund (FWF) through Projects No. SFB IR-ON F25-14 and No. SFB NextLite F49-P10, by the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme (FP7/2007-2013) under REA Grant Agreement No. PIOF-GA-2011-303228 (Project NOLACOME) as well as computational resources on the Vienna Scientific Cluster (VSC). D.L. and S.G.J. were supported in part by the AFOSR MURI for Complex and Robust On-chip Nanophotonics (Dr. Gernot Pomrenke), Grant No. FA9550-09-1-0704 and by the Army Research Office through the Institute for Soldier Nanotechnologies (ISN), Grant No. W911NF-07-D-0004. A.D.S. and A.C. acknowledge the support by the NSF Grants No. DMR-0908437 and No. DMR-1307632 and the Yale HPC for their computing resources.## Appendix A Details for lasing mode solution method

In this appendix, we provide further details on setting up the Newton-Raphson iteration for lasing modes. First, we describe how to fix the phase and normalization for each mode (as mentioned in Sec. III.2). We choose a point and a constant unit vector such that is nonzero for all lasing modes. This condition is usually satisfied provided that is neither far outside the cavity nor a point of high symmetry. We then define the quantity and rescale the field such that the physical field becomes and the rescaled field satisfies

(15) |

With this redefinition, the rescaled field has a fixed phase and a normalization that distinguishes it from the trivial solution . Further we treat the quantity as a separate unknown that contains the mode’s amplitude. The spatial hole-burning [Eq. (8)] then becomes

Now, we describe how to construct the vector of unknowns which, after rescaling, should contain , , and . First, the discretized fields are described by -component complex vectors . The real unknowns for each mode can then be written in block form as

(16) |

The vector we use for the Newton-Raphson method contains all in sequence, since the lasing modes are all coupled together through the spatial hole-burning interaction and thus must be solved simultaneously.

Next, we construct the equation vector by discretizing the operator into a sparse complex matrix . In the discrete basis, Eq. (6) becomes which gives complex scalar equations, and the normalization condition that fixes the phase [Eq. (15)] becomes the complex scalar equation , where is the discrete-basis representation of the vector function consisting of the unit vector at point and zero everywhere else. The real and imaginary parts of these complex equations can be written in block form as

(17) |

The vector we use for the Newton-Raphson method contains all in sequence, due to the intermodal coupling.

Finally, we describe how to construct the (real) Jacobian matrix (which is real), which consists of blocks that each have size and have the block form

We explicitly construct these blocks by taking derivatives of column blocks of with respect to row blocks of , as defined in Eqs. (16) and (17). First, we see that , while all other blocks of with are zero. Second, we have the columns

and

where the derivatives of are diagonal complex matrices that can be obtained straightforwardly by discretizing the same derivatives of the complex scalar function . (In the case that exact outgoing radiation conditions are used for , the matrix for may also depend on and this dependence must also be included in the derivative.)

Finally, the remaining blocks are given by

where is the Kronecker , and is the matrix discretization of the real tensor function

with the outer product taken over the real and imaginary parts of the vector components of .

## Appendix B FEM formalism

In this appendix we provide details on how to implement our SALT solution strategy with a high order finite element method (-FEM) Demkowicz (2007); Solín (2005). Most importantly, our approach does not depend on this specific discretization method, but it entails several significant advantages. Specifically, -FEM can handle highly complex, irregularly shaped geometries and the higher order discretizations tend to exponentially increase the accuracy of the computations (if all of the boundary discontinuities and corner singularities are properly taken into account).

In order to obtain the discretized formulation, we truncate the open problem to a bounded computational domain . Then we multiply Eq. (6) with an arbitrary test function and integrate over the domain of the cavity. Using the Green’s formula this leads to the weak formulation

(18) | ||||

where denotes the outer normal vector at the boundary . The boundary term has to be replaced by a term incorporating the radiation condition at infinity. Formally, this can be done by , where is the Dirichtlet-to-Neumann operator Nédélec (2001). In one dimension and in two dimensions for TE modes the Maxwell equations reduce to the well-known scalar Helmholtz equation. Furthermore, in the one-dimensional case the boundary integral can be simply replaced by

(19) |

In higher dimensions an appropriate representation of the open boundary becomes more sophisticated, as mentioned in Sec. III.5. For a detailed discussion see Monk (2003), but Eq. (19) can also be simply used as the first-order approximation of the operator.

The unknown laser modes are sought as a linear combination of element basis functions such that

(20) |

where are piecewise polynomials with local support and is the number of degrees of freedom of the system. For more details on the choice of such a basis, based on high-order elements, we refer to Demkowicz (2007). We use the notation with a complex FEM-coefficient vector . Inserting the ansatz Eq. (20) into Eq. (18), extracting sums out of the integrals, and assembling the contributions for all elements we arrive at the following finite-element scheme in matrix form:

(21) |

where the sparse -matrices are the stiffness matrix

corresponding to the Laplacian or curl-curl term, respectively, the mass matrix

containing the passive dielectric function, the matrix

which only involves the boundary elements and incorporates the outgoing boundary condition (19), and the nonlinear contribution

which accounts for the nonlinear coupling including the spatial hole burning effect.

## Appendix C FDFD formalism

For finite-difference calculations shown in the main text, the discretization code implemented in Liang and Johnson (2013); Liang (2012) was used. The complex electric fields were discretized on an pixel grid of equally spaced points with the operator being conveniently discretized using second-order centered differences on a Yee lattice Taflove and Hagness (2005). To impose outgoing boundary conditions, additional pixels of PML were added at the boundaries with the appropriate absorption, as explained in Sec. III.5. For each mirror symmetry in a geometry, we were able to halve the computational domain by replacing the PML at the lower walls with the corresponding boundary conditions of the mirror plane. Furthermore, for the cases of TM () and TE polarization (), the problem size can be reduced by factors of and , respectively, by projecting and into the nonzero field components only. Additionally, 2D calculations were performed by setting and the boundary condition in the direction to be periodic. 1D calculations were performed by doing so for both the and directions.

## Appendix D Contour integral method

This section reviews the algorithm for solving a nonlinear EVP of the form as discussed in Ref. [47]. The corresponding inverse matrix can generically be written as

where are the desired complex eigenvalues and are the corresponding left and right eigenvectors. The residual term is holomorphic. Using the residue theorem the following two matrices can be defined:

where is a diagonal matrix with the diagonal entries corresponding to all the poles of the inverse matrix inside the contour which in turn are the eigenvalues of . Before we show how the desired matrix is computed from the two matrices and , we discuss the realization of the contour integration which is obtained by numerical quadrature. Very fast (i.e., exponential) convergence is achieved with the trapezoidal rule (Kress, 1998, Theorem 9.28), if the contour is an analytic curve such as a circle or an ellipse. Moreover, inverting the large matrix for each quadrature point is numerically expensive and may even be infeasible given that the inverses are fully populated. This can be remedied by an approximation scheme that exploits the fact that the rank of the matrices and is given by the number of eigenvalues inside the contour and is thus very small compared to : For a random matrix , we merely evaluate at each quadrature point on the contour, i.e., we reduce the computational cost to the solution of linear systems for each quadrature point. The parameter has to be selected slightly larger than the expected size of the number of eigenvalues inside the contour. To obtain the matrix , we first compute the (reduced) singular value decomposition (SVD) of

Then we define the matrix and observe that and are similar, i.e., for some matrix . Therefore their eigenvalues are the same such that the desired eigenvalues can be obtained from the reduced eigenvalue problem

In this short sketch we have assumed that is exactly the number of eigenvalues inside the contour so that exists; if is larger, the SVD of has to be replaced with a rank-revealing variant. In total, the algorithm involves a (dense) eigenvalue problem, an SVD of an matrix, and sparse linear solves, where is the number of quadrature points for the contour integration. Consequently, the bottleneck of a large scale eigenvalue problem is essentially shifted to the solution of perfectly parallelizable linear systems. The contour integration method assumes that the contour does not pass through eigenvalues; since eigenvalues close to the contour could affect the convergence of the quadrature scheme, we compute the residuum of the computed eigenvalues found by the algorithm before proceeding.

Our computations are performed with ellipsoidal contours and the trapezoidal rule. If contours are desired that are no longer analytic but only piecewise analytic, then the trapezoidal rule should be replaced by other exponentially convergent schemes such as Gaussian or Clenshaw Curtis quadrature. However, with quadrature error control in place, this method guarantees to find all eigenvalues inside the contour, but avoids spurious eigenvalues at . Other eigenvalue solvers either compute all eigenvalues which is neither realistic nor necessary or rely on local convergence properties.

## Appendix E Multimode Lasing and Vector Solitons

In this appendix we discuss an interesting mathematical connection between the SALT lasing equations in the multimode regime and the nonlinear incoherent “vector solitons” in photorefractive media Christodoulides et al. (1996). The noninstantaneous nonlinearity in such media allows more than two components of solitons to be self-trapped and thus to form vector solitons based on the mutual incoherence between their various components (the term “vector” refers here to the locked components that propagate together). Following Christodoulides et al. (1996), the normalized multimode soliton equations for the scalar electric-field envelopes of interacting beams are:

(22) |

where is the total intensity at infinity and is the peak nonlinear index. Equation (22) describes coupled beams in a saturable optical photorefractive medium. Such an interaction can form vector solitons that consist of two or more components mutually self-trapped in the nonlinear medium. In the small-intensity regime (Kerr limit), and when , the resulting governing equations describe the so-called Manakov solitons Manakov (1974); Kang et al. (1998), which in 1+1 dimensions are known to be integrable. The stationary soliton solutions have the form with the soliton eigenvalues being real numbers. Substituting this ansatz into Eq. (22) leads to the following nonlinear eigenvalue problem:

(23) |

subject to the boundary conditions . By comparing with Eq. (6), we stress that Eq. (23) is an eigenvalue problem defined on the whole real line , while the laser problem is restricted to the domain of a finite cavity length, . Hence, as a result of the different boundary conditions Eq. (23) admits a continuum family of soliton solutions for any nonzero positive value of the soliton eigenvalue , whereas Eq. (6) admits a discrete family of solutions corresponding to different lasing frequencies that are determined self-consistently. Note also that in extension of the above soliton equations, in SALT not only the eigenvectors but also the eigenvalues appear non-linearly. In spite of such characteristic differences we can envision applying various numerical methods developed in the field of nonlinear optics and soliton theory to solve efficiently the multimode SALT equations. In particular, the most commonly used techniques used to numerically determine vector soliton solutions are the multi-dimensional Newton-Raphson method Efremidis et al. (2003), the self-consistent method Cohen et al. (2003) and the recently developed spectral renormalization method Ablowitz and Musslimani (2005).

## References

- Painter et al. (1999) O. Painter, R. K. Lee, A. Scherer, A. Yariv, J. D. O’Brien, P. D. Dapkus, and I. Kim, Science 284, 1819 (1999).
- Noda et al. (2000) S. Noda, A. Chutinan, and M. Imada, Nature (London) 407, 608 (2000).
- Altug et al. (2006) H. Altug, D. Englund, and J. Vuckovic, Nat. Phys. 2, 484 (2006).
- Ellis et al. (2011) B. Ellis, M. A. Mayer, G. Shambat, T. Sarmiento, J. Harris, E. E. Haller, and J. Vuckovic, Nat. Photon. 5, 297 (2011).
- Haken (1985) H. Haken, Light: Laser Light Dynamics (North-Holland, Amsterdam, 1985).
- (6) This is the same theory as the ab initio self consistent laser theory abbreviated as AISC, but the name has been changed.
- Türeci et al. (2006) H. E. Türeci, A. D. Stone, and B. Collier, Phys. Rev. A 74, 043822 (2006).
- Türeci et al. (2008) H. E. Türeci, L. Ge, S. Rotter, and A. D. Stone, Science 320, 643 (2008).
- Ge et al. (2010) L. Ge, Y. D. Chong, and A. D. Stone, Phys. Rev. A 82, 063824 (2010).
- Christ and Hartnagel (1987) A. Christ and H. L. Hartnagel, IEEE Trans. Microwave. Theory and Tech. 35, 688 (1987).
- Champagne et al. (2001) N. J. Champagne, J. G. Berryman, and H. M. Buettner, J. Comput. Phys. 170, 830 (2001).
- Chew et al. (2001) W. C. Chew, J. M. Jin, E. Michielssen, and J. M. Song, eds., Fast and Efficient Algorithms in Computational Electromagnetics (Artech House, 2001).
- Joannopoulos et al. (2001) J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals: Molding the Flow of Light (, Princeton Univeristy Press, 2001).
- Ziolkowski et al. (1995) R. W. Ziolkowski, J. M. Arnold, and D. M. Gogny, Phys. Rev. A 52, 3082 (1995).
- Nagra and York (1998) A. S. Nagra and R. A. York, Antennas Propag. 46, 334 (1998).
- Ge et al. (2008) L. Ge, R. J. Tandy, A. D. Stone, and H. E. Türeci, Opt. Express 16, 16895 (2008).
- Cerjan et al. (2012) A. Cerjan, Y. Chong, L. Ge, and A. D. Stone, Opt. Express 20, 474 (2012).
- Press et al. (2007) W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes: The Art of Scientific Computing, 3rd ed. (Cambridge University Press, 2007).
- Hénon et al. (2002) P. Hénon, P. Ramet, and J. Roman, Parallel Comput. 28, 301 (2002).
- Guillaume (1999) P. Guillaume, SIAM J. Matrix. Anal. Appl. 20, 575 (1999).
- Taflove and Hagness (2005) A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed. (Artech House, London, 2005).
- Shin and Fan (2012) W. Shin and S. Fan, J. Comp. Phys. 231, 3406 (2012).
- Liertzer et al. (2012) M. Liertzer, L. Ge, A. Cerjan, A. D. Stone, H. E. Türeci, and S. Rotter, Phys. Rev. Lett. 108, 173901 (2012).
- Türeci et al. (2009) H. E. Türeci, A. D. Stone, L. Ge, S. Rotter, and R. J. Tandy, Nonlinearity 22, C1 (2009).
- Stano and Jacquod (2013) P. Stano and P. Jacquod, Nat. Photon. 7, 66 (2013).
- Hisch et al. (2013) T. Hisch, M. Liertzer, D. Pogany, F. Mintert, and S. Rotter, Phys. Rev. Lett. 111, 023902 (2013).
- Chua et al. (2011a) S.-L. Chua, Y. Chong, A. D. Stone, M. Soljacic, and J. Bravo-Abad, Opt. Express 19, 1539 (2011a).
- Ge (2010) L. Ge, Steady-state Ab Initio Laser Theory and its Applications in Random and Complex Media, PhD thesis, Yale University (2010).
- Chong and Stone (2012) Y. D. Chong and A. D. Stone, Phys. Rev. Lett. 109, 063902 (2012).
- Yariv (1989) A. Yariv, Quantum Electronics, 3rd ed. (Wiley, New York, 1989).
- Boyd (2001) J. P. Boyd, Chebyshev and Fourier Spectral Methods, 2nd ed. (Dover, New York, 2001).
- Brandstetter et al. (2014) M. Brandstetter, M. Liertzer, C. Deutsch, P. Klang, J. Schöberl, H. E. Türeci, G. Strasser, K. Unterrainer, and S. Rotter, Nat. Commun. 5, 4034 (2014).
- Sargent et al. (1977) M. Sargent, M. O. Scully, and W. E. Lamb, Laser Physics (Westview, Boulder, CO, 1977).
- Fu and Haken (1991) H. Fu and H. Haken, Phys. Rev. A 43, 2446 (1991).
- Arecchi et al. (1984) F. T. Arecchi, G. L. Lippi, G. P. Puccioni, and J. R. Tredicce, Opt. Commun. 51, 308 (1984).
- Oppo et al. (1986) G. L. Oppo, A. Politi, G. L. Lippi, and F. T. Arecchi, Phys. Rev. A 34, 4000 (1986).
- Pillay et al. (2014) J. C. Pillay, Y. Natsume, A. D. Stone, and Y. D. Chong, Phys. Rev. A 89, 033840 (2014).
- Trefethen and Bau (1997) L. N. Trefethen and D. Bau, Numerical Linear Algebra (SIAM, 1997).
- Davis (2006) T. A. Davis, Direct Methods for Sparse Linear Systems (SIAM, 2006).
- Balay et al. (1997) S. Balay, W. D. Gropp, L. C. McInnes, and B. F. Smith, in Modern Software Tools in Scientific Computing, edited by E. Arge, A. M. Bruaset, and H. P. Langtangen (Birkhäuser, 1997) pp. 163–202.
- Su and Bai (2011) Y. Su and Z. Bai, SIAM J. Matrix. Anal. Appl. 32, 201 (2011).
- Peters and Wilkinson (1979) G. Peters and J. H. Wilkinson, SIAM Review 21, 339 (1979).
- Sleijpen et al. (1996) G. L. Sleijpen, A. G. Booten, D. R. Fokkema, and H. A. van der Vorst, BIT Numer. Math. 36, 595 (1996).
- Ruhe (2003) A. Ruhe, J. Math. Sci. 114, 1854 (2003).
- Voss (2004) H. Voss, BIT Numer. Math. 44, 387 (2004).
- Asakura et al. (2010) J. Asakura, T. Sakurai, H. Tadano, T. Ikegami, and K. Kimura, Japan J. Appl. Math. 27, 73 (2010).
- Beyn (2012) W.-J. Beyn, Lin. Alg. Appl. 436, 3839 (2012).
- Tolstikhin et al. (1998) O. I. Tolstikhin, V. N. Ostrovsky, and H. Nakamura, Physical Review A 58, 2077 (1998).
- Berenger (1994) J.-P. Berenger, J. Comput. Phys. 114, 185 (1994).
- Bramble and Pasciak (2008) J. H. Bramble and J. E. Pasciak, Math. Comp. 77, 1 (2008).
- Ge et al. (2011) L. Ge, Y. D. Chong, S. Rotter, H. E. Türeci, and A. D. Stone, Phys. Rev. A 84, 023820 (2011).
- Vahala (2003) K. J. Vahala, Nature 424, 839 (2003).
- Dettmann et al. (2009) C. P. Dettmann, G. V. Morozov, M. Sieber, and H. Waalkens, Europhys. Lett. 87, 34003 (2009).
- Ge et al. (2014) L. Ge, D. Liu, S. G. Johnson, S. Rotter, H. Cao, A. Cerjan, A. D. Stone, and H. E. Tureci, (2014), in preparation.
- Lin et al. (2001) S. Y. Lin, E. Chow, S. G. Johnson, and J. D. Joannopoulos, Opt. Lett. 26, 1903 (2001).
- Frei et al. (2008) W. R. Frei, H. T. Johnson, and K. D. Choquette, J. App. Phys. 103, 033102 (2008).
- Kao and Santosa (2008) C.-Y. Kao and F. Santosa, Wave Motion 45, 412 (2008).
- Lu and Vuckovic (2010) J. Lu and J. Vuckovic, Opt. Express 18, 3793 (2010).
- Chua et al. (2011b) S.-L. Chua, C. A. Caccamise, D. J. Phillips, J. D. Joannopoulos, M. Soljačić, H. O. Everitt, and J. Bravo-Abad, Opt. Express 19, 7513 (2011b).
- Cerjan et al. (2014) A. Cerjan, Y. D. Chong, and A. D. Stone, arXiv:1406.6659 [physics] (2014), arXiv: 1406.6659.
- Demkowicz (2007) L. Demkowicz, Computing With Hp-adaptive Finite Elements: Frontiers, Chapman and Hall/CRC Applied Mathematics and Nonlinear Science Series (Chapman & Hall/CRC, 2007).
- Solín (2005) P. Solín, Partial Differential Equations and the Finite Element Method, 1st ed. (John Wiley & Sons, 2005).
- Nédélec (2001) J. C. Nédélec, Acoustic and Electromagnetic Equations (Springer, New York, 2001).
- Monk (2003) P. Monk, Finite Element Methods for Maxwell’s Equations, Numerical Mathematics and Scientific Computation (Clarendon Press, 2003).
- Liang and Johnson (2013) X. Liang and S. G. Johnson, Opt. Express 21, 30812 (2013).
- Liang (2012) X. Liang, Modeling of fluids and waves with analytics and numerics., PhD thesis, MIT (2012).
- Kress (1998) R. Kress, Numerical Analysis (Springer, 1998).
- Christodoulides et al. (1996) D. N. Christodoulides, S. R. Singh, M. I. Carvalho, and M. Segev, App. Phys. Lett. 68, 1763 (1996).
- Manakov (1974) S. Manakov, Sov. Phys. JETP 38, 692 (1974).
- Kang et al. (1998) J. Kang, J. Aitchison, G. Stegeman, and N. Akhmediev, Optical and Quantum Electronics 30, 649 (1998).
- Efremidis et al. (2003) N. K. Efremidis, J. Hudock, D. N. Christodoulides, J. W. Fleischer, O. Cohen, and M. Segev, Phys. Rev. Lett. 91, 213906 (2003).
- Cohen et al. (2003) O. Cohen, T. Schwartz, J. W. Fleischer, M. Segev, and D. N. Christodoulides, Phys. Rev. Lett. 91, 113901 (2003).
- Ablowitz and Musslimani (2005) M. J. Ablowitz and Z. H. Musslimani, Opt. Lett. 30, 2140 (2005).