PT -symmetry enabled stable modes in multi-core fiber

Open systems with balanced gain and loss, described by parity-time ( PT -symmetric) Hamiltonians have been deeply explored over the past decade. Most explorations are limited to finite discrete models (in real or reciprocal spaces) or continuum problems in one dimension. As a result, these models do not leverage the complexity and variability of two-dimensional continuum problems on a compact support. Here, we investigate eigenvalues of the Schr¨odinger equation on a disk with zero boundary condition, in the presence of constant, PT -symmetric, gain-loss potential that is confined to two mirror-symmetric disks. We find a rich variety of exceptional points, re-entrant PT -symmetric phases, and a non-monotonic dependence of the PT -symmetry breaking threshold on the system parameters. By comparing results of two model variations, we show that this simple model of a multi-core fiber supports propagating modes in the presence of gain and loss.

After their experimental realizations in numerous platforms, it has become clear that PT -symmetric Hamiltonians accurately model open systems with balanced, spatially separated gain (V I > 0) and loss (V I < 0) [31].Their standard phenomenology is as follows: starting from the Hermitian Hamiltonian H 0 with real spectrum and Dirac-orthogonal eigenfunctions, as the imaginary part of the potential V I (x) is increased, two or more real eigenvalues undergo level attraction, become degenerate, and then develop into complex-conjugate pairs.This eigenvalue degeneracy, called exceptional point (EP) degeneracy [32][33][34], is characterized by the coalescence * tgrat2@pdx.edu of corresponding eigenfunctions and lowering of the rank of the Hamiltonian operator.Due to the antilinearity of the PT -operator, an eigenfunction f n (x) is simultaneously an eigenfunction of the PT operator with unit eigenvalue if and only if the corresponding eigenvalue λ n is real; if λ n is complex, then it follows that PT f n (x) is an eigenfunction with complex-conjugate eigenvalue λ * n .The transition across the EP from a real spectrum to one with complex-conjugate eigenvalues is called PTsymmetry breaking transition, since the corresponding eigenfunctions lose that symmetry, PT f n (x) ̸ = f n (x).
Here, we investigate a two-dimensional continuum model on a compact domain subject to hardwall (vanishing eigenfunctions) boundary condition in the presence of constant PT -symmetric complexvalued potentials.In one dimension, such potential leads to a single PT -symmetry breaking transition when the strength of the imaginary part of the potential, γ, exceeds a threshold γ PT set by the Hermitian Hamiltonian H 0 .The two-dimensional case, we will show, differs dramatically.It leads to multiple transitions where stable modes (real spectra) change into amplifying and leaky modes (complex conjugate eigenvalues) as γ is increased.More surprisingly, we also find PT -restoring transitions where, as the pure gain-loss potential V I is increased, amplifying and leaky modes are pairwise stabilized.
The plan of the paper is as follows.In Sec.I we introduce the model and recall the Hermitianlimit results for a cylindrical waveguide.Section III contains the outline of the numerical procedure we use for discretization.Results for eigenspectra and eigenfunctions across multiple PT -breaking and restoring transitions are shown in Sec.IV.Section V concludes the paper.As a physical example, we consider a lengthwise uniform, multi-core fiber with circular crosssection of radius R = 1 (purple) centered at the origin in the x 1 -x 2 plane, a lossy core of radius ρ centered at distance d/2 from the origin (green), and a gain-medium core of the same radius ρ centered at the mirror-symmetric location (pink).The Maxwell's equation for a TM-mode electric field E(x) = E 3 (x 1 , x 2 ) exp[i(kx 3 − ωt)]x 3 propagating along the fiber gives a Schrödinger-like eigenvalue problem with a potential V (x 1 , x 2 ) proportional to the local index of refraction.When expressed in terms of suitable dimensionless variables, this problem can be cast as follows. Let denote the disk of dimensionless radius r centered at p. Our problem is set in the domain Ω = B 1 (0, 0) where a purely imaginary gainloss potential is introduced into non-intersecting left and right subdomains D L = B ρ (−d/2, 0) and D R = B ρ (d/2, 0) where ρ < R and 2ρ ≤ d ≤ 2(R − ρ) ensures that the two domains do not intersect (Fig. 1).The eigenvalue problem is to find complex-valued functions f n (x 1 , x 2 ) on Ω that vanish on the boundary ∂Ω and are square-integrable (f n ∈ L 2 (Ω)), and complex numbers λ n such that ).The dimensionless gain-loss potential V (x) is given by ( We define the parity operator P : ), i.e., P mirrors functions about the second axis.It is easy to see that P is a linear, self-adjoint, and unitary operator in L 2 (Ω).The antilinear time-reversal operator T : For unbounded operators H defined on a proper subspace dom(H) ⊂ L 2 (Ω) rather than all of L 2 (Ω), namely H : dom(H) → L 2 (Ω), Eq.( 3) means f, PT f ∈ dom H and the equality (3) holds.The operator of interest to us, A = −∆ + V (x 1 , x 2 ), is unbounded, and its domain is given by dom(A) = H 2 (Ω) ∩ ∂H 1 0 (Ω).Here H k (Ω) denotes the Sobolev space of square-integrable functions all of whose derivatives of order at most k ≥ 1 are also square integrable and ∂H 1 (Ω) denotes the subspace of H 1 (Ω)-functions that vanish on the boundary ∂Ω.It is straightforward to check that A is PT -symmetric.
To investigate the eigenvalues of A(γ), we start with the Hermitian limit of Eq.( 2), γ = 0.In this case, the cylindrical symmetry in the x 1 -x 2 plane gives unnormalized eigenfunctions in polar coordinates, where the corresponding eigenvalue λ mp is determined by the p th zero of the m th Bessel function, J m ( λ mp ) = 0, which enforces the hard-wall boundary condition f | ∂Ω = 0. Except for m = 0, these solutions with exp(±imθ) are degenerate, and represent positive and negative angular momentum states respectively.The semi-analytically obtained eigenvalues λ mp of A(γ = 0) are the starting point for computing eigenvalue trajectories λ mp (γ).They also serve to verify our numerical methodology by bench-marking it against the γ = 0 case.

III. NUMERICAL DISCRETIZATION IN 2D
In one dimension, the discretization of the Schrödinger operator leads to tridiagonal matrix with no corner elements.Since the boundary of the connected domain in this case consists of its two end points, no special considerations are needed to ensure that the discretization scheme conforms to it.Two-dimensional domains, on the other hand, require more care.Let us denote the complex L 2 (Ω)inner product by ⟨•|•⟩.For any smooth function g vanishing on ∂Ω, the eigenvalue equation ( 1) implies The finite element method imposes the same equation on the Lagrange finite element [35] space X h consisting of continuous functions, vanishing on the boundary ∂Ω, which are polynomials of degree at most p in each mesh element.Here the mesh is a geometrically conforming mesh of triangles subdividing the domain, respecting the material interfaces, with curved elements near the circular boundaries and interfaces.The subscript h indicates the maximal diameter of all elements in the mesh.As h becomes smaller, or as p becomes larger, the discretization becomes finer and dim X h becomes larger.
Our numerical method computes the eigenvalues of a discretization A h : X h → X h of the infinitedimensional operator A. It is defined by for all f h , g h ∈ X h .Namely, we compute an eigenvalue approximation λ h,n and right eigenfunction f h,n satisfying Standard finite element theory [36] can be used to show that the approximate eigenpairs (f h,n , λ h,n ) converge to the exact ones under suitable assumptions as h → 0; the symmetry of the mesh is immaterial for such convergence.The right eigenfunction f h ∈ X h in ( 6) is equivalently given by Using a (non-orthogal) basis ψ i of finite element shape functions, Eq. ( 7) can be converted to a matrix eigenvalue problem Ax = λBx where A ij = ⟨ψ i |A h ψ j ⟩ and B ij = ⟨ψ i |ψ j ⟩.This generalized eigenproblem is then solved for a cluster of selected eigenvalue using a contour integral method called the FEAST algorithm [37,38].While much of our ensuing analysis use meshes without symmetry, we have also experimented with meshes with parity symmetry that are invariant under reflection by the vertical axis (x 1 = 0).On such meshes, the discretized A h is exactly P T -symmetric; specifically, (5) implies that for all f h ∈ X h , recovering the perfect analogue of Eq.( 3) on the discrete space X h .Figure 2: Flow of eigenvalues λ mp (γ) for the first seven eigenvalues.All of them except m = 0 cases are doubly degenerate at γ = 0 and this degeneracy is lifted with increasing γ.The first PT -symmetry breaking transition occurs at γ ≈ 97 immediately followed by PT -restoring transition near γ = 102 (detailed in the second plot).This is followed by a significantly broad PT -broken region, and another small PT -breaking and restoring transition.This reentrant PT -symmetric phase in a model with single gain-loss parameter is uncommon.
In Figure 2a we track real (pink) and imaginary (blue) parts of the lowest seven eigenvalues as a function of gain-loss strength γ.All but two of them are doubly degenerate, and the gain-loss potential lifts this degeneracy.As γ is increased, we see level attraction, leading to degeneracy and emergence of complex conjugate pair, indicated by equal and opposite imaginary parts.The first such transition occurs near γ = 97, shown in detail in the second plot of Fig. 2 (with rescaled axes), and is followed by a PT -restoring transition where the spectrum becomes purely real again near γ = 102.The third plot of Figure 2 shows the zoomed-in and rescaled view of another such small window near γ = 290.

V. CONCLUSION
In this work, we have investigated the rich diversity of PT -symmetry breaking transitions that arise in a two-dimensional, continuum, bounded domain with uniform gain or loss potentials confined to parity-symmetric regions.We have found multiple PT -symmetry breaking and restoring transitions, and analyzed the behavior of eigenfunctions and transition thresholds as a function of system parameters.Our results suggest that stable, propagating modes are supported in multi-core fibers with gain and loss regions.

Figure 1 :
Figure 1: The P T -symmetric model