Lasing condition for trapped modes in subwavelength–wired PT–symmetric resonators

The ability to control the laser modes within a subwavelength resonator is of key relevance in modern optoelectronics. This work deals with the theoretical research on optical properties of a PT–symmetric nano–scaled dimer formed by two dielectric wires, one is with loss and the other with gain, wrapped with graphene sheets. We show the existence of two non–radiating trapped modes which transform into radiating modes by increasing the gain–loss parameter. Moreover, these modes reach the lasing condition for suitable values of this parameter, a fact that makes these modes to achieve an ultra high quality factor that is manifested on the response of the structure when it is excited by a plane wave. Unlike other mechanism that transform trapped modes into radiating modes, we show that the variation of gain–loss parameter in the balanced loss–gain structure here studied leads to a variation in the phase difference between induced dipole moments on each wires, without appreciable variation in the modulus of these dipole moments. We provide an approximated method that reproduces the main results provided by the rigorous calculation. Our theoretical findings reveal the possibility to develop unconventional optical devices and structures with enhanced functionality. PACS numbers: 81.05.ue,73.20.Mf,78.68.+m,42.50.Pq


Introduction
It is known that, under certain conditions, a plasmonic system can support trapped electromagnetics modes, which are electromagnetic non radiating oscillations that stay localized inside the structure [1]. These plasmon modes are characterized by an ultrahigh quality factor provided that the structural symmetry remains unbroken. In fact, the narrow resonances characterizing these mode excitations rely on the breaking degree of the symmetry which allows the trapped mode to couples with free photons [2]. In this sense, plasmon trapped modes can be considered as a kind of symmetry protected states. Due to Ohmic loss in the plasmonic material, the eigenfrequencies of the trapped states are complex valued. However, by introducing gain material elements into the structure its possible to compensate this material loss and, as a consequense, to reach the lasing condition for which the eigenfrequency associated to an eigenmode is real valued [3].
Because of their fundamental properties as well as their potential applications, the study of hybrid systems composed of gain material media and plasmonic materials, which are loss media, is a topic of continuous increasing interest [4,5,6,7,8]. In areas such as condensed matter and surface optics, the amplification of eigenmodes by stimulated emission of radiation has played a key role in the interpretation of a wide variety of experiments, the understanding of various fundamental properties of solids and the engineering of nanolaser devices [9,10,11]. In particular, new phenomena associated with parity-time (PT) symmetry have been observed in optical systems with balanced loss and gain [12,13]. These optical systems belong to a large family of non-Hermitian systems, which can have a real spectrum provided that the system be invariant under combined operations of parity (P) and time-reversal (T) symmetry [14,15].
Possibilities have been widened towards PT symmetric structures that incorporate graphene as plasmonic material. Doped graphene allows the propagation of surface plasmons with low Ohmic losses, i.e., with high quality factor, from terahertz to nearinfrared range [16]. Moreover, the plasmon resonance spectrum, and consequently the optical responses of the system, can be tuned by varying the doped level on graphene. In this context, several works have focused on the influence that long living and tunable surface plasmons on graphene have on the optical responses of a PT-symmetric system [17,18,19,20].
The present work contains all the ingredients above entered and focusing on the study of a sub-wavelength graphene plasmonic system with balanced volume losses. In particular, our system is composed of two parallel dielectric cylinders, one is with loss and the other with the same level of gain, both wrapped with graphene shets. Here, the sole effect of graphene is to provide LSP resonances in sub-wavelength wires. Although other more complex structures such as an array of many or infinite cylinders can be designed to present PT-symmetry, our motivation is based in the simplicity of the dimer structure that we have chosen, since it allows us to understand the underlying physics behind the PT-symmetry. We use a rigorous formalism based on Mie theory to calculate the eigenmodes dependence with the gain-loss parameter for both trapped and radiating eigenmodes. To do this, we solve the boundary value problem with the corresponding boundary conditions without an incident wave (homogeneous problem). The procedure requires the analytic continuation of the eigenvalue, the frequency in our case, in the complex plane. Analytic continuation is inevitable, even in the case of media with no intrinsic losses, since the open nature of our resonator generates non-null imaginary parts due to radiation losses.
In addition, as it is well established nowadays, the modal lasing analysis in an open plasmonic resonator, as our PT-symmetric structure, can be carried out by applying the lasing eigenvalue problem (see [24] and references therein), in which the modal eigenfrequencies are assumed to be real valued at the lasing condition (or lasing threshold). Following this concept, we solve the homogeneous problem to find the real valued eigenfrequencies, gain-loss parameters and chemical potentials for each modal lasing condition. Moreover, we study the eigenmode influence on the optical response of the system when it is excited by a plane wave near the lasing condition.
By using the quasistatic approximation valid in the long wavelength limit, we show that despite the graphene ohmic losses, a fact that gives rise to a complex spectrum, the structure exhibits a set of properties in common with a PT-symmetric system. For instance, two eigenmode branches coalesce at an exceptional point and, for the gain-loss parameter above certain threshold, these both branches are repelled in the direction of the imaginary part of frequency, making that one of these branches achieves the lasing condition. In addition, we demonstrate that branches eigenmode, which are trapped modes in case of null gain-loss parameter, transform into radiating modes that can be excited by plane wave incidence when the gain-loss parameter is increased. Unlike previous works [21,23], where transitions from trapped to radiating eigenmodes are achieved by producing a difference between the modulus of individual dipole moments on each particle forming the dimer, i.e., by producing a weakly asymmetry with respect to the center of the dimer modifying the dipole moment amplitudes, here, we show that the increment of the gain-loss parameter leads to a change in the phase difference between these individual dipole moments, maintaining their modulus values constant. This paper is organized as follows. In section 2 we present a brief description of the rigorous method used in this work to calculate the scattering of a dimmer composed of two graphene wires. From this method, using the quasistatic approximation, we deduce analytical expressions for eigenfrequencies and eigenvectors as a function of geometrical and constitutive parameters that explain the main features calculated with the rigorous method. In section 3 we present results of two parallel dielectric cylinders tightly coated with a graphene layer, one of then with small inner losses and the other one with the same level of gain. Concluding remarks are provided in Section 4. The Gaussian system of units is used and an exp(−i ω t) time-dependence is implicit throughout the paper, where ω is the angular frequency, t is the time coordinate, and i = √ −1. The symbols Re and Im are used for denoting the real and imaginary parts of a complex quantity, respectively.

Rigorous description of the fields and scattering efficiencies
We consider a cluster consisting of two parallel and non overlapping cylindrical dielectric wires, one with gain, ε a = ε 1 −iε i (ε i > 0), and the other with equal loss, ε b = ε 1 +iε i , as shown in Figure 1. Both wires have the same radius R a = R b = R and are wrapped with a graphene sheet. The system is embedded in a lossless and non-magnetic dielectric denoted as medium v with permittivity ε v . In this case, a PT symmetry around the central axis, denoted by O, is fulfilled. We assume that the radius R is sufficiently large to describe the optical properties of the wires as characterized by the same surface conductivity as planar graphene (see appendix Appendix A). We denote by r j (r), φ j (r) (j = a, b) the polar coordinates of a point at position r with respect to the local origin O j . A plane wave radiation impinges on the wires with an angle of incidence φ inc with respect to the y axis. Although some enhanced optical effects related to invisibility modes are observed for s polarization (electric field along the z axis) [22], this work focus on p polarization (magnetic fields along the z axis) for which the electric field in the graphene coating induces electric currents directed along the azimuthal direction φ j (r) and LSPs exist in the graphene circular cylinder. In this way, the incident magnetic field (along the z axis) can be written in a system linked to the j-cylinder as [23], where r j , φ j are the polar coordinates of the j-cylinder, k v = √ ε v ω c , is the modulus of the photon wave vector in medium v, ω is the angular frequency, c is the vacuum speed of light and J m (x) is the nth Bessel function. The scattered magnetic field in medium v (r j (r) > R) can be written as a superposition of the field scattered by each of the cylinders, where H m (x) is the nth Hankel functions of the first kind. Note that the jth term of the summation corresponds to the field scattered by the jth cylinder linked to the local system with origin O j . In the region inside the cylinders, r j (r) < R, the transmitted field is written as where j = a, b. To find the unknown complex amplitudes of the reflected b jm and transmitted a jm fields (2) and (3), we use the usual boundary conditions and the addition theorem for Bessel and Hankel functions [25]. This theorem allow us to write one of the terms in (2), associated to the scattered field of one of the cylinders (for example, the j = b cylinder) in the other local coordinates (the j = a cylinder). In this way, the scattered field (2) will be represented in the form of expansions in Hankel functions written in the local coordinate j = a. By replacing this expression and Eq. (3) with j = a into the boundary conditions along the surface r a = R of the j = a cylinder, one obtain a set of 2 equations for the 2 × 2 unknown amplitudes. Similarly, we can write the scattered field (2) in the local coordinate j = b and use the boundary conditions on the surface of the j = b cylinder to obtain other set of 2 equations for the 2 × 2 unknown amplitudes. However, we closely follow a variant of the method, developed in [26], that allows to reduce to half the dimension of the system of equations. The detailed of this implementation has been given in [23], and leads to the following system of equations for the amplitudes b jm where b l and Q l are vectors whose coordinates are the elements b lm and respectively, T lj is the matrix with elements and S j is the matrix with elements S jmn = s jm δ mn , where s jm are the elements of the scattering matrix associated to the jth cylinder [27] where x = k j R, y = k v R and the prime denotes the derivative with respect to the argument. Knowing the total electromagnetic field allows us to calculate the scattering cross sections. The time-averaged scattered power is calculated from the integral of the radial component of the complex Poynting vector flux through an imaginary cylinder of length L of radius ρ 0 which envelops the graphene wire system (see Fig. 1), It is convenient to express each of the terms in Eq. (2) in a same coordinate O [23]. Substituting the obtained expression into Eq. (9) we obtain where Inserting Eq. (10) into (8) and taking into account the wronskian , after some algebraic manipulations, we obtain the power scattered by the cylinders, The scattering cross section is defined as the ratio between the total power scattered by the cylinders, given by Eq. (12), and the incident power P inc intersected by the area of all the cylinders. It is known that the scattering cross section and the near to the cylinders field are strongly affected by complex singularities in the field amplitudes b jm . Singularities occur at complex locations and they represent the frequency of the eigenmodes supported by the cylinders system. Complex frequencies of these modes are obtained by solving the homogeneous problem, i.e., by imposing the vectors Q a and Q b in Eq. (4) be zero. Then, a condition to determine eigenfrequencies is to require the determinant of the matrix in Eq. (4) to be zero, This condition corresponds to the full retarded dispersion equation (FR) for eigenmodes and it determines the complex frequencies ω = ω R + iω I (ω I < 0) in terms of all the geometrical and constitutive parameters of the system.

Quasistatic approximation. A simple model based on two coupled electric dipoles
Although the rigorous treatment represented by Eq. (13) gives us all the kinematics and dynamics characteristics of the PT-symmetric eigenmodes, the method lacks of analytical expressions that explain the main dependencies with both geometrical and constitutive parameters. For the purpose of showing this dependencies, by applying the quasistatic method, here, we reduce the full treatment provided by the homogeneous part of Eq. (4) to a simple 2x2 matrix description as follows.
Assuming that the radius R of cylinders is much smaller than the wavelength λ = 2π/ω, the problem can be treated using the dipole approximation where the dimer eigenfunctions calculated with (4) are, as a good approximation, a superposition of single plasmons with angular momentum m = ±1 linked to each cylinder. In this way, only four coordinates, b j±1 (j = a, b), define the amplitudes vector.
We first consider the case in which the induced dipole moments are along the ±x directions, i.e., the local magnetic field associated to each of the cylinders has a dependence ≈ sin φ j (j = a, b). As a consequence, the amplitudes Here, the subscript j − 1 stand for cylinder j and angular momentum m = −1. In this way, the matrix equation for the modal amplitudes reduces to a 2×2 matrix system for amplitudes b a1 and b b1 , where we have used . Using the small argument asymptotic expansions for Bessel and Hankel functions [25], it follows that s j can be written as are the dipolar polarizability of the cylinders j = a or j = b, respectively, g(ω) = − ω 2 g ω 2 +iωγg , ω 2 g = 4e 2 µ 2 R is the effective plasma frequency for the dipolar mode [28]. Note that the value of α a differs from α b in the sign of the gain-loss parameter ε i , signs + and − correspond to j = a and j = b respectively. Near the resonance frequency, the polarizability α j (j = a, b) can be written as where ω j = ω jR + iω jI (ω jI < 0) is the complex pole of α j , i.e., the eigenfrequency of the single graphene cylinder j, g ′ (ω) is the derivative of g(ω) and The single cylinder eigenfrequencies for cylinders a and b are written as [28] where the signs + and − corresponds to j = a and j = b, respectively, and are the damping rates corresponding to the ohmic loss in graphene covers and dielectric cores, respectively. We have used the fact that ε i << ε 1 + ε v in the last equality in Eqs. (19). From Eq. (20) we see that the real part of the eigenfrequencies do not depends on ε i , a fact that also results by solving the fully retarded dispersion equation for a single graphene cylinder [29]. By taking into account the small argument z = k v R ab in the Hankel functions, the off diagonal elements in the matrix in Eq. (14) are written as Therefore, the matrix equation (14) takes the form where A j j = a, b are given by Eq. (18). Eq. (24) gives us a simple description for the dimer dynamic. For large separations, A j /R 2 ab << ω a (or R/R ab << 1), the matrix (24) is diagonal and thus the eigenfrequencies correspond to that of each individual graphene wires composing the dimer, i.e., the plasmonic wires do not interact between them. Taking into account the PT parameters, the real parts of both frequencies are the same whereas their imaginary parts differ in 2Γ core . For small enough values of R ab , extra diagonal terms take appreciable values and a splitting between the real parts of eigenfrequencies occur. From Eq. (B8) we can see that this splitting is proportional to R 2 /R 2 ab ω 0 . For system (24) to have a non-trivial solution, its determinant must be equal to zero, a condition which can be written as A detailed developed of the right hand side of this equation can be seen in appendix Appendix B. By replacing expression (B8) into Eq. (25) and solving for ω eigenfrequencies, where in the last equality we have used Eq.
ab ω 0 (ε 1 +εv) ε v while their imaginary parts remain degenerated at the value −γ g /2. From equation (24), we see that the eigenvector associated to the upper branch verify b a1 = −b b1 pointing out that the induced dipole moments on each of the cylinders move out of phase, i.e., the state corresponds to a trapped mode, while the eigenvector associated to the lower branch verify b a1 = b a1 pointing out that this mode corresponds to a radiating mode that can be excited, for example, by an incident plane wave [23]. We use the notation +x−x and +x+x to refer to the upper and lower branches, respectively. Since the terms inside the root have opposed signs between them, the splitting between upper (+x−x) and lower (+x+x) branches monotonically decreases as ε i is increased until reaching the value for which the eigenfrequencies coalesce at the exceptional point ω = −i γg 2 + ωg √ ε 1 +εv . At the same time, the imaginary parts of both eigenfrequencies remain equal. Beyond the exceptional point and for the gain-loss parameter ε i > ε ep , the imaginary parts of the eigenfrequencies bifurcate, while their real parts remain degenerated. In Figure 2 we illustrate in the complex plane the above described trajectories for the upper and lower branches (26) as parametric functions of the gain-loss parameter ε i . We have taken the straight segment z = 0, with ℜ z > 0 as the cut line for the square root function in Eq. (26) so that √ −w 2 = iw for w real and positive. We now consider the case of polarization along the y axis, in which each of the induced dipole moments are along ±y direction, i.e., the local magnetic field associated to each of the cylinders has a dependence ≈ cos φ j (j = a, b). As a consequence, the amplitudes b j1 = −b j−1 (j = a, b). In this case, the matrix equation for these amplitudes is where we have used 2H ′ 1 (z)/z = H 0 (z) − H 2 (z), and H ′ 1 (z) is the derivative of H 1 (z). Following the same steps as in the case of x polarization, the matrix equation (28) takes the form where A j are given by Eq. (18). Note that the matrix in Eq. (29) differs from that in Eq. (24) only in the signs of the non diagonal terms. Therefore, the eigenfrequencies corresponding to induced dipole moments along y direction are formally given by Eq. (26). Unlike the x polarization case, for which the upper frequency branch for ε i = 0 corresponds to a trapped mode, From Eq. (29) and taking ε i = 0, we see that the eigenvector associated to the upper branch verify b a1 = b b1 , i.e., it corresponds to a radiating mode. Conversely, it is straightforward to verify that the lower frequency branch corresponds to a trapped mode for which b a1 = −b b1 . It is worth noting that in this case, ±y oscillations, it is convenient to take the cut of the complex square root function √ z as the straight line z = iw (w real and positive) so that √ −w 2 = −iw. In this way, beyond the exceptional point the upper branch +y+y moves away from the real axis whereas the lower branch +y−y reaches the real axis, as shown in Figure 2b.

Results
We consider a system of two dielectric cylinders with permittivities ε a = 2.13 + iε i (ε i > 0), ε b = 2.13 − iε i for lossy and gain cylinders, respectively. The radii R a = R b = R = 0.03µm and both cylinders are coated with a graphene monolayer and immersed in vacuum (ε v = 1). The graphene parameters are: temperature T = 300K, chemical potentials µ 1 = µ 2 = 0.5eV and the carriers scattering rates γ 1 = γ 2 = 0.1meV. The positions of the cylinders are r a = −0.05µmx, r b = 0.05µmx (center to center distance R ab = 0.1µm) and the gap between them is ∆ = 0.04µm. Figure 3a shows the trajectory of the eigenfrequencies in the complex plane as a parametric function of the gain-loss parameter ε i calculated by solving the full retarded dispersion equation (FR). To solve this equation, we use a Newton-Raphson method adapted to treat complex variable. Four branches are observed, two of them (+x+x and +x−x) corresponds to both cylinders polarized along the x axis (dipole moments along the x axis) while the other two branches (+y−y and +y+y) corresponds to the case in which both cylinders are polarized in the y axis (dipole moments along y axis). We are using the notation of the asymptotic case when ε i = 0 for naming the dimer surface plasmons, so that +y+y (+x+x) configuration corresponds to the two dipole moments moving in phase on the y axis (x axis), and the +y−y (+x−x) configuration corresponds Dashed lines correspond to QA branches plotted in Figure 2. (b) Scattering cross section for p-polarized incident waves for illumination direction along (φ inc = 90 • ) and perpendicular (φ inc = 0) to the axis joining the cylinder centers. The permittivity of the cylinders are ε a = 2.13 + iε i , ε b = 2.13 − iε i , the radius R a = R b = R = 0.03µm and the gap ∆ = 0.04µm. The graphene parameters are T = 300K, γ a = γ b = 0.1meV and µ a = µ b = 0.5eV.
to the dipole moments oscillating in opposite phase on the y axis (x axis). For instance, the branch +x+x start at frequency ω/c = 0.8589226 − i0.5789274 10 −3 µm −1 for ε i = 0, where both dipole moments are in phase, and it moves to the right side leaving away from the real axis. Moreover, the real part of this trajectory approaches asymptotically to the value ω/c = 0.8868µm −1 corresponding to the real part of the eigenfrequency of a single graphene cylinder [28]. On the contrary, the branch +x−x start at ω/c = 0.9107941−i0.253260710 10 −3 µm −1 for ε i = 0, with both dipole moments in opposite phase, and it approaches to the real axis where the lasing threshold at ω crit /c ≈ 0.8954615m −1 is reached for ε i = ε crit = 0.166 (point A in Figure 3a).
On the other hand, Figure 3a also shows the trajectory of the branch corresponding to the polarization along the y axis. The branch +y−y start at frequency ω/c = 0.8594583−i0.25349871 10 −3 µm −1 for ε i = 0, where both dipole moments are in opposed phase, and it moves toward the right side, reaching the real axis at point C where the critical eigenfrequency ω crit /c = 0.8755407µm −1 for a critical value of the gain-loss parameter ε crit = 0.1715275. On the other hand, the +y+y branch start at frequency ω/c = 0.9102954 − i0.5560595 10 −3 µm −1 for ε i = 0 and it moves to the left side leaving away from the real axis. As in the x polarized case, both branches are repelled in the direction of the imaginary axis allowing the +y−y branch to achieve the lasing threshold at point C.
In order to understand the gain-loss compensation near the critical points at which the eigenfrequencies are almost real, in Figure 3b we plot the scattering cross sections for a plane wave impinging at an angle φ = 0 (electric field along the x axis) and φ = 90 • (electric field along the y axis). The corresponding gain-loss parameter is ε i = 0.166 for φ = 0, and ε i = 0.1715 for φ = 90 • , i.e., it values are near to the critical values at which lasing conditions are achieved. In the φ = 0 case, we observe that the scattering cross section is enhanced at frequency near ω/c ≈ 0.895m −1 that agree well with the lasing frequency for the +x−x branch calculated by solving the eigenmode problem (point A in Figure 3a). Moreover, we observe another peak (less intense) near ω/c = 0.873µm −1 , a value that falls near the real part of the eigenfrequency of an state with ε i = 0.166 but corresponding to the +x+x branch (point B in Figure 3a). A similar behavior presents the scattering curve for φ = 90 • and ε i = 0.1715. From Figure 3b we observe a very sharp peak at frequency that coincides with the lasing frequency ω crit /c ≈ 0.8755407µm −1 calculated by solving the eigenmode problem (point C in Figure 3a). Moreover, a second peak is observed at a frequency ≈ 0.894µm −1 associated to the excitation of the state on the +y+y branch for ε i = 0.1715 (point D in Figure 3a).
The question that arises from the above results is how eigenmodes on the +x−x and +y−y branches, which correspond to trapped modes for ε i = 0, can be excited with a plane wave by varying the gain-loss parameter ε i . Furthermore, these branches reach their lasing threshold for a critical value of the gain-loss parameter. To find a response, we have calculated the eigenvectors, which contains all the field amplitudes of the eigenmodes. In particular, we have verified that coefficients with |m| = 1 are orders of magnitude less than those corresponding to m = ±1, suggesting that the dimer plasmons can be considered, as a good approximation, as a superposition of single plasmon with m = ±1 linked to each cylinder. As a consequence, to gain further insight into the underlying physics of these branches excitations, we applied the QA as follows. Without loss of generality, we consider the case for which the induced dipole moments on cylinders are in ±x direction. By replacing Eq. (26) into Eq. (24) and using Eq. (19), we find the following relation between the dimer amplitudes where the − and + signs correspond to the +x−x and +x+x branches, respectively. From this equation we see that the modulus of the b a1 and b b1 amplitudes are equal providing that the gain-loss parameter be less than ε ep . Taking into account the fact that b a−1 = b a1 and b b−1 = b b1 and using Eq. (2), we can write the field scattered by the cylinders as Comparing this expression with that corresponding to a single dipole moment p along the x axis and centered at the origin (H ≈ pH 1 (k v r) sin(φ)), we deduce that Eq. (31) corresponds to a superposition of two fields, one of them due to a dipole moment of amplitude p a ≈ b a1 centered at the cylinder a and other due to a dipole moment of amplitude p b ≈ b b1 centered at the cylinder b. Since b a1 = b b1 regardless of ε i (ε i < ε ep ), the induced dipole moments p a = p b and, as a consequence, both emitted fields for each of the dipoles have the same intensity. In addition, from Eq. (30) we see that the phase difference Φ between b a1 and b b1 amplitudes, and thus between p a and p b , is shifted from π to π/2 when ε i increases from 0 up to above the value ε ep . This fact implies that the eigenmodes on the +x−x trajectory pass from trapped to bright by increasing the gain-loss parameter. Our calculation confirm this expectation, as can be seen in Figure 4 where we show plots of the phase Φ as a function of the gain-loss parameter ε i by applying the FR dispersion rigorous method (continuous line) and using Eq. (30) (dashed line). From this Figure we see that the curve calculated with the QA agree well with that calculated by using the FR dispersion method. In particular, we see that the phase Φ = π for ε i = 0, which means that the eigenmodes are dark states, and monotonically decreases for until reach the lasing modes A (for +x−x brach) and B (for +y−y branch). It is worth noting some similarities and differences between the eigenmode branches calculation by using the FR dispersion equation and those calculated by using the QA. On the one hand, the +x+x and +x−x branches in Figure 3a are repelled in the direction of the imaginary axis, as predicted by the QA for a PT-symmetric system (Figure 2), a feature that allows the +x−x branch to achieve the lasing threshold at point A. Furthermore, the critical value (27) calculated by using the QA results ε ep = 0.18, a value that agree well with the gain-loss parameter for which the lasing  threshold is achieved in Figure 3a. Moreover, the phase difference between the induced dipole moments on each cylinders calculated by FR and QA methods matches quite well. On the other hand, +x+x and +x−x branches in Figure 3a do not start at points with the same value of their imaginary parts as occur for the branches in Figure 2 and, as a consequence, these trajectories does not coalesce at an exceptional point as in Figure 2. This is true because the lack of a radiation losses term in the QA, i.e., the radiation losses would prevent the system from having all the properties of a full loss compensated PT-symmetric structure, shown in Figure 2, such as the existence of an exceptional point. In order to study the behavior with the chemical potential on graphene, we set µ a = µ b = µ and vary the values of µ. In Figure 5 we have plotted the lasing frequency ω crit and the corresponding gain-loss parameter ε crit as functions of µ, calculated with the rigorous FR method. From Figure 5a we can see that the lasing frequency for both +x−x and +y−y branches are increasing functions of µ. This fact can be understood by taking into account that the real part of the eigenfrequency corresponding to a single graphene cylinder, which falls between the lasing frequencies for +y−y and +x−x branches, is proportional to √ µ (see Eq. (20)). On the other hand, from Figure 5b we see that the critical value of the gain-loss parameter for which the lasing condition is achieved is a decreasing function of the chemical potential µ. This behaviour has a similarity with that presented by a single graphene cylinder for which has been demonstrated that the gain level to achieve the lasing condition decreases with the chemical potential value [29]. µ (eV) Figure 5. Lasing frequency (a) and critical gain-loss parameter (b) as functions of the chemical potential µ (µ a = µ b = µ) for branches +x−x and +y−y. The calculations have been carried out by applying the FR dispersion method. Geometrical and other constitutive parameters are the same as in Figure 3.

Conclusions
In conclusion, we have analytically studied the scattering and the eigenmode problems for a dimer composed of two graphene coated dielectric cylinders, one of them with loss and the other with the same level of gain. We have demonstrated the existence of two branches, corresponding to trapped modes when ε i = 0, that reach the lasing conditions for suitable values of the gain-loss parameter. While the phase difference between the induced dipole moments on individual cylinders changes from π to a value near π/2 when the gain-loss parameter is incremented, a fact that provides the mode transformation from trapped to radiating modes, the modulus of each individual dipole moments maintains equals in between.
Other mechanisms to transform a trapped mode into a resonant observable which can be excited by a plane wave have been reported in other works. All these methods are based on the introduction of a small asymmetry with respect to the center of the dimer by producing dissimilar dipole moment modulus on each of the cylinders. Interestingly, here, we found that in the transformation from trapped to radiating eigenmodes on both +x−x and +y−y branches, the modulus of the individual dipole moments does not change, while it is changed the phase between them. A distinction between these two kind of mechanism, the phase variation and modulus variation mechanisms, to transform a trapped mode into a resonant observable has not been reported before to our knowledge. We believe that our results will be usefull for a deeper understanding of the PT-symmetric LSP characteristics and that the LSP effects we have demonstrated opens up possibilities for practical applications involving sub-wavelength laser structures.