Hermitian Hamiltonians with real eigenvalues and orthogonal eigenstates are generally used to describe closed (isolated) systems, while for an open system, there is the additional interaction of the states with the environment, resulting in a non-Hermitian Hamiltonian with complex eigenvalues and nonorthogonal eigenstates ^[1–3]. The extension of a Hermitian to a non-Hermitian Hamiltonian is generally achieved by introducing complex potentials or gain and loss ^[4–7]. Note that if the complex potentials have parity–time (PT) symmetry ^[4,8], the system experiences a phase transition from the unbroken phase to the broken one across an exceptional point (EP), where the eigenvalues and eigenstates coalesce. All these features lead to many interesting applications in various optical systems, including PT-symmetric lasers ^[9,10], coherent absorption ^[11], Bragg reflectors ^[12,13], nonreciprocal light transmission ^[14,15], and others ^[16–18]. Also, for passive PT-symmetric systems, through gauge transformation, the structure can behave in a PT-symmetric-like fashion ^[19,20]. When the nonlinearity included in these PT systems has been taken into account, it gives rise to a wide array of new phenomena. In particular, the gain saturation is unavoidable in amplifying optical waveguides ^[21,22].

Among various photonic structures, coupled optical waveguides provide one of the most appealing and experimentally accessible platforms for implementing the idea of PT symmetry ^[23,24], while for asymmetric coupling systems arising from asymmetric geometry or arbitrary gain and loss, the associated properties of PT symmetry are lost ^[24,25]. Yet with some subtle approaches, PT-symmetric-like properties can be recovered to a certain extent ^[26]. In Ref. ^[27], we revealed the absence of EP, even for a situation where the system apparently fulfills the PT symmetry condition with uniform gain and loss. This is ascribed to the different modes coupling together with unequal propagation constants. In addition to the different modes coupling, the detuning of the coupled modes arising from nonuniform gain and loss can also eliminate the EP, but more apparently yields lossless modes, which are of great importance for photonic signal processing.

In this paper, we present such a study and examine the evolution of propagation constants for four coupled waveguides in a square array with nonuniform gain and loss. The non-Hermitian Hamiltonian obtained from the classical coupled-mode theory (CMT) ^[28–30] can be used to describe the propagation dynamics; however, it fails to capture the influence of the gain and loss, particularly for heavy gain/loss involvement and for identical modes coupling. In order to deal with this problem, another well-established approach, referred to as general CMT, from the first-principle point of view, has been developed ^[31]. Based on this approach, a detailed description of the effect of the gain and/or loss factor on the propagation constants is illustrated. In particular, we present the conditions for the emergence of a purely real propagation constant via properly designed gain and loss distribution, utilizing the classical CMT and general CMT simultaneously. All of these are fundamentally related to the energy balance of the coupled modes in their spatial propagation channels (locations in the cross sections of the waveguide array).

We begin with a simple two-states model to describe the system that has two coupled propagation modes. The Hamiltonian is written as where β1 is the propagation constant, κ12 and κ21 are the coupling strengths of the two modes, γ1(2) defines the intrinsic gain (γ1(2)>0) or loss (γ1(2)<0), and Δβ represents the additional and tunable propagation constant introduced in the second mode. In Eq. (1), j is the imaginary unit.

The eigenvalues of the Hamiltonian, i.e., Eq. (1), take the form Note that the gain/loss of the propagation modes are characterized by the imaginary parts of the propagation constants (obtained either from theoretical eigenvalues of Eq. (1) or from numerical full-wave simulation, such as COMSOL Multiphysics ^[32]). In experiments, the gain/loss may be achieved via erbium doping optically pumped with a laser beam ^[27]. In the nonlinear regime, the presence of gain saturation results in a mismatch with the loss. So a rational value of the gain/loss coefficient should be carefully adopted in the experiment in the linear regime, and also in the nonlinear regime, which is not discussed in this paper.

If Δβ=0 and γ1=−γ2, the system reduces to the general PT-symmetric case, and experiences a transition from a fully real spectrum to a complex conjugate spectrum ^[4,8]. If Δβ=0, γ1=0, and γ2<0, the system corresponds to the passive PT symmetry case ^[33]. If Δβ≠0, the system no longer has real eigenvalues, and the EP will not emerge. In this work, we mainly focus on the deviation of the propagation constants, stemming from the additional nonuniform gain and loss for the coupling of the originally identical modes, in addition to the coupling of different modes supported in the four coupled waveguides in a square array [see Fig. 1(a)].

To investigate and compare the two different cases (e.g., the identical mode coupling and different mode coupling cases), we consider a usual and unique structure: four coupled waveguides in a square array ^[27]. The Hamiltonian of this system is not transparent, as each waveguide supports two degenerate fundamental modes. To circumvent that, the system is considered as a coupled structure with supermodes supported in diagonal waveguides and the supermodes supported in off-diagonal waveguides. The supermode here refers to a hybridized mode from the fundamental ones in each individual circular waveguide. It turns out that the total structure can be represented by four isolated 2×2 Hamiltonians (the modes are located in diagonal and off-diagonal waveguides), as shown in Figs. 1(b)–1(e) (respectively corresponding to coupling combinations 1–4, labeled as C1–C4). Coupling schemes C1–C4 are specific combinations of two decoupled modes (from the four modes, labeled as M1, M2, M3, and M4), supported in the diagonal and off-diagonal waveguides as shown in Figs. 1(b)–1(e), respectively.

Here, a more rigorous and systematic approach, (i.e., the nonorthogonal vector CMT based on the linear superposition of the modes for individual waveguides ^[30]) can be utilized to theoretically analyze the mode interactions in a two-state model. The coupled-mode equation reads where βA(B) is the propagation constant for each supermode supported in diagonal (off-diagonal) waveguides, κA(B) is the coupling strength, and aA(B) is the corresponding normalized amplitude. γA(B) denotes the gain (γA(B)>0) or loss (γA(B)<0) coefficient. Here, the wavelength is set to 1550 nm. We note that κA(B) and γA(B) can be calculated based on the field distribution of the decoupled modes and their overlapping integral ^[28,29]. The cross-interacting term is essentially treated as a complex inner product between the fields: Here, uA(B) represents the normalized electric field distribution of the individual waveguide system, and SA(B) is the integral cross section where waveguides A(B) are located. In this work they were obtained numerically by the finite element method (FEM) ^[32]. The associated coefficients, especially κA(B) and γA(B), for different coupling schemes (C1–C4) are shown in Table 1, calculated by Eqs. (4) and (5). Note that γA(B) are assumed to be proportional to ϵA(B),I in the conventional CMT.

According to the eigenvalues of Eq. (3) and the coefficients in Table 1, if the permittivity of the system obeys ϵ(x)=ϵ*(−x) (* stands for complex conjugate), it possesses general PT symmetry for the combinations C2 and C3, while for the other two combination schemes, the system loses the EP ^[27]. The absence of the EP is ascribed to the different propagation constants between the two decoupled supermodes, which originate in the intrinsic properties of the structure.

Up to now, we have considered the cases of uniform gain and loss distributions briefly. But for nonuniform gain and loss, particularly for ϵB,I=0 and ϵA,I=ϵI<0, which are analogous to passive PT symmetry ^[33], two of the total coupling cases (C2 and C3) should experience a similar transition from an unbroken phase with identical eigenvalues to a broken phase with unequal ones, according to Eq. (3), based on the classical CMT.

The corresponding numerical results obtained from FEM simulation are shown in Fig. 2 for varying amounts of gain and loss. Although they are not consistent with the above prediction, we will carefully demonstrate the reason. Note that the EPs of all four of these combinations are lost. However, after a critical amount of loss is added, the total transmission of the waveguide array increases, even though the loss amount is increased. This is similar to the passive PT symmetry, leading to loss-induced transmission ^[34].

The CMT presented in Section 3 is valid for a definite and conserved optical power of the whole system ^[30], while for non-Hermitian systems with arbitrary gain and loss, and especially for systems with strong gain/loss, the total integrated power is not a conserved quantity ^[31]. Thus, a general CMT analysis is necessary to deal with these non-Hermitian systems. This method is based on the variational principles, in which the scalar inner product for non-Hermitian systems is used ^[31].

In this section, we briefly outline the theoretical procedure for constructing the general CMT; for further details, see ^[31,35]. We assume that the forward (+z) propagation modes supported in the unperturbed systems are expressed as φ=[Em+Hm+]T=[em+hm+]Texp[j(ωt−βmz)]. The corresponding adjoint fields, i.e., the backward (−z) propagation modes are ϕ=[Em−Hm−]T=[em−hm−]Texp[j(ωt+βmz)]. The relationships between the fields are as follows: em+=[exeyez]m, hm+=[hxhyhz]m, em−=[exey−ez]m, and hm−=[−hx−hyhz]m. Here em± (hm±) is the vector electric (magnetic) field in the transverse cross section for forward (+) and backward (−) propagation. Importantly, the counterpropagating modes (em± and hm±) have a definite relationship, so one can be deduced from the other. The subscript m denotes the index of the supported supermodes.

The normalized forward fields satisfy the Maxwell equations where ∇t=x∂/∂x+y∂/∂y and x, y, and z are unit vectors. When a small perturbation (Δϵ) is taken into account, the fields of the perturbed system can be approximately written as a linear combination of the unperturbed adjoint fields, i.e., Here, they also satisfy the Maxwell equations By considering only the effect of the perturbation, we can take certain operations of Eq. (6) multiplied by Σnhn− minus Eq. (10) multiplied by hm+, and Eq. (7) multiplied by Σnen− minus Eq. (11) multiplied by em+(instead of multiplying by the complex conjugate of the fields); then, taking the integral operator over the transverse cross section, the equation finally can be simplified to The coefficients in Eq. (12) are as follows: For the unperturbed situation (Δϵ=0), it is easy to have bmn=−j(βn−βm)pmn. As a result, Eq. (12) reduces to Equation (16) can be rewritten as an eigenvalue problem, where a=(a1,a2,…,an)T, (H1)mn=βnpmn−kmn, and (H2)mn=pmn. The eigenvalue β of the Hamiltonian H is the propagation constant of the propagating modes, and the corresponding eigenvector a denotes the proportion of the decoupled modes (i.e., the unperturbed modes) written as Eqs. (8) and (9). As a consequence, the propagation constant β is the result of the original unperturbed modes coupling due to the involved perturbation. This is the key equation referred to as the general CMT.

In view of Eq. (16), a small perturbation Δϵ has a great effect on the final propagation constants. We would utilize the general CMT to analyze the evolution of the propagation constants for the cases shown in Figs. 2(c) and 2(d). Intuitively, it comes from the detuned Δβ due to the nonuniform gain and loss for the coupling modes. In the preceding section, it was mentioned that the Hamiltonian was considered as the coupled modes formed in the diagonal and off-diagonal waveguides, and therefore we would particularly examine the evolution of each mode supported in the diagonal waveguides when the imaginary part of the permittivity is introduced. Here, the four modes (labeled as M1, M2, M3, and M4) are exactly orthogonal for the corresponding Hermitian system. As such, they are treated as single-mode interactions in the general CMT, i.e., m=n in Eq. (16). In this case, Eq. (16) becomes Here m corresponds to modes M1, M2, M3, and M4, respectively. The propagation constant after perturbation can be calculated step by step through the fields obtained by numerical simulation based on the FEM. The results are shown in Fig. 3. It is obvious that as the imaginary part of the permittivity (ϵA,I=ϵI) increases from zero, the real part of the propagation constants increases, and the increasing speed accelerates, while the imaginary part keeps growing linearly. We note that the real part is even and the imaginary part is odd with respect to ϵI. The results of the FEM agree well with the general CMT. Figure 3 further shows that the imaginary part of the permittivity has great influence on the real part of the propagation constants, especially for large ϵI values. This is the main cause of the introduction of Δβ for the coupled modes, which is responsible for the results in Figs. 2(c) and 2(d). The emergence of zero Δβ is also the reason for the formation of PT symmetry with uniform gain and loss.

According to the above analysis and the results presented in Fig. 3, it is possible to achieve a stable propagation mode (without any attenuation and amplification during the propagation) in systems with nonuniform gain and loss distributions. Lossless propagation modes are ubiquitous in PT-symmetric systems, for example, in coupled waveguides with balanced gain and loss below the EP, and attenuation is frequently encountered in the absence of PT symmetry. In this section, we show the possibility of achieving purely real eigenvalues in waveguide systems with nonuniform gain and loss, combining the classical CMT and the general CMT. Here, we use the coupling combination C1 as the example, and the results can be applied to the other cases (e.g., C2, C3, and C4).

To simplify the discussion, we introduce the degree of freedom α (ϵA,I=−ϵI, ϵB,I=αϵI). A stable propagation mode requires that Im(β±)=0. By inserting into Eq. (3) the parameters (listed in Table 1) obtained from the classical CMT, we get a phase diagram describing the properties of the mode coupling (M3 and M4) as shown in Fig. 4, where the black solid line corresponds to the case of Im(β+)=0, representing one lossless eigenstate and one lossy one. Notice that the red dashed line corresponds to Im(β−)=0, for which the system has a gain eigenstate and a lossless one. For parameters in the black and the red lines, the system is tuned to a point (let us call it the ‘α- point’ for convenience) where one of the two supermodes we considered conserves its energy during propagation. For this reason, the whole parameter space is separated into three regions labeled I, II, and III in Fig. 4. It is easy to see that region I consists of double lossy eigenstates, region II has a gain eigenstate and a lossy one, and region III contains two gain eigenstates. Importantly, we can see that there exist two independent solutions for Im(β−)=0 (for a certain ϵI value range) and a single solution for Im(β+)=0. Also, there exists a limit for the choice of ϵI to fulfill Im(β−)=0. Thus, the properties of the entire system can be distinguished by the specific parameter pair (ϵI, α). The propagation behaviors of the system are then discussed in some detail, essentially using Fig. 4 as a rough guideline. The numerical results for a purely real propagation constant are also shown in Fig. 4 with different symbols (circles and squares).

To check the detailed modal properties, we first focus on a specific α-point of (ϵI,α)=(0.6,0.9459), which falls right on the black solid curve in Fig. 4. With these settings, the system supports two supermodes: one of effective index β+/k0=2.53 and the other of β−/k0=2.44−0.005j. To see the consequence of such a nonuniform combination of gain and loss, we use the fourth-order Runge–Kutta method to solve Eq. (3). Figure 5 shows the evolution of the slowly varying intensity |ai(z)|2 (i=A,B) as a function of the propagation distance z. This is for the case when the lossy waveguide mode aA is initially excited. It is seen that the optical field is able to be coupled into the gain waveguide mode aB. There is also obvious beam oscillation for both modes, resulting from the interference of the asymmetric mode and their reciprocal coupling. The oscillation period Lo≈17.2μm is determined by the mismatch of the real parts of the supermode refractive indices as Lo=2π/Re(Δβ′), where Δβ′=β+−β−. Alternatively, for the case of the initial excitation of the gain waveguide mode aB, there is similar propagation behavior.

For another particular pair (ϵI,α)=(0.5,0.9459), which falls into the phase region I in Fig. 4, similar beam oscillation can be observed [see Figs. 5(c) and 5(d)]. However, the optical field propagation suffers from overall attenuation, decaying approximately to zero after a distance z≈400μm. For this system, the two supermodes are both lossy, in agreement with the phase diagram in Fig. 4. It is noted that the attenuation at relatively small ϵI values in region I is quite small, which indicates that light can propagate over long distances (in the order of millimeters).

We emphasis that the differences between the results of classical CMT and FEM in Fig. 4 are partially due to the fact that the coefficients in Table 1 are obtained by the original field distributions of the corresponding Hermitian system. This treatment neglects the perturbation of mode profile and coupling by including the imaginary part of the permittivity. To have a clearer and more exact theoretical modeling, the abovementioned general CMT in Section 4 is applied here. Our system can be represented by a two-mode coupling model, so based on Eq. (16), the corresponding general coupled-mode equation can be written as Here, the perturbation is Δϵ=jϵA,I (diagonal waveguides) and Δϵ=jϵB,I (off-diagonal waveguides).

Figure 6 shows the corresponding results from the FEM, the classical CMT, and the general CMT. For different α values, the results of the FEM and general CMT all agree quite well. When α increases from the loss-dominated value (α<1) to the gain-dominated value (α>1), the real parts of the two bands approach, and the imaginary parts experience a transition, meaning that the system goes from two lossy eigenstates to two gain ones. This exactly illustrates and corroborates the phase region partition in Fig. 4. For Im(β−)=0 at ϵI≈0.6, there exist two choices of α, and the corresponding dispersion relations are shown in Figs. 6(c)–6(f). This clearly demonstrates the difference between the selection of Im(β−)=0 and Im(β+)=0 in order to have a stable mode. Despite the very good agreement with respect to the FEM results, we note that the general CMT could not directly provide a phase boundary prediction as shown in Fig. 4, while the classical CMT can give a rough estimation of these phase boundaries. Therefore, the classical CMT and the general CMT can be complementary to each other in terms of verification and prediction, apart from brute force numerical methods.

In conclusion, we have revealed the eigenvalue dynamics for a system of four geometry identical waveguides with nonuniform gain and loss. Based on the first principle of the Maxwell equations, the general CMT is used to understand the origin of the offset of originally equal propagation constants in the presence of permittivity perturbation, especially for the real part. The offset leads to asymmetric mode coupling and the absence of EP. Also, the system can be tuned to support supermodes with real propagation constants, via a carefully designed gain/loss ratio in the individual waveguides. We show that such a system has a phase transition that interfaces two regimes: one with both coupled propagation modes being lossy, the other with one gain and one lossy mode. It can also be optimized to demonstrate a different transition from the phase with all gain states to a phase with one gain and one lossy. The required tuning parameter is the gain/loss ratio α: for the former, α<1, while for the latter, α>1. The additional degree of freedom α is introduced to modify the energy distribution, in order to compensate the mode asymmetry-induced imbalance. It is noted that α closely relies on the propagation constant of the decoupled modes. Furthermore, our results may be extended to plasmonic or plasmonic–dielectric hybridized waveguide systems.