Scattering-induced ampliﬁcation of two-dimensional plasmons: electromagnetic modeling

Using two rigorous electromagnetic approaches, we study plasmon scattering in two-dimensional systems and show that plasmon ampliﬁcation is possible in the presence of dc currents. Two scenarios are considered: plasmon scattering from an interface between diﬀerent two-dimensional channels and plasmon reﬂection from electric contacts of arbitrary thickness. In each case, the eﬀect of a dc current of the plasmon reﬂection and transmission coeﬃcients, and the plasmon power are both quantiﬁed. A resonant system is studied where plasmon roundtrip gain may exceed unity, showing the possibility of plasmon generation.

Using two rigorous electromagnetic approaches, we study plasmon scattering in twodimensional systems and show that plasmon amplification is possible in the presence of dc currents. Two scenarios are considered: plasmon scattering from an interface between different two-dimensional channels and plasmon reflection from electric contacts of arbitrary thickness. In each case, the effect of a dc current of the plasmon reflection and transmission coefficients, and the plasmon power are both quantified.
It was also suggested that, in the presence of a dc current, plasmon scattering by discontinuities in 2D systems may lead to plasmon amplification. The first theoretical study was performed by Dyakonov and Shur 32 who showed that plasmon oscillations in the channel of a field-effect transistor may become unstable provided specific boundary conditions are realized at the source and the drain. Their model was developed further for different boundary conditions 33,34 and geometries 35,36 . Amplification of plasmons scattering from interfaces between 2D systems (see Fig. 1) has received less attention. In contrast to transistor source and drain, however, such interfaces are partially transparent to plasmons, and may serve as building blocks for more sophisticated geometries with higher gain, such as plasmonic crystals and Bragg mirrors, for example. Experimental studies of 2D plasmons in the presence of dc currents are scarce (e.g. Refs. [11,[37][38][39]); evidence of plasmon amplification from interface scattering has come to date from observing radiation from double-gratinggate transistors 40 . Several theoretical studies 41-43 concentrated on analytical models for the plasmon transmission and reflection coefficients. A typical approximation is to reduce the three-dimensional interface to a single line of contact between two-dimensional systems. This quasi one-dimensional approach tends to disregard that plasmons reflecting from interfaces acquire additional phase, whose existence has been revealed, in the absence of dc current, by models that take into account fields across entire interfaces 18,19,21,22,44 . On the other hand, prior numerical studies 45,46 of two-dimensional plasmons in the presence of dc current mostly concentrated on excitation of entire geometries by a free-space electromagnetic wave or a filament. Plasmon scattering at single interfaces is then obscured by the excitation, multiple reflections, losses, and matching to free space radiation.
The recent advances in analytical and numerical electromagnetic modeling of passive interfaces motivate further efforts in studying such interfaces in the presence of dc current.
In this paper, we first consider amplification of plasmons incident on an interface formed by two different two-dimensional systems, see Fig. 1. To this end, we develop two models: a numerical model based on variational solution of an integral equation for the field at the interface (Sec. III A), and a full-wave numerical model based on solution of coupled Maxwell's and hydrodynamic equations (Sec. III B). We compare results of the three models to each other in Sec. III C. We then apply the full-wave model to plasmon reflection from electric contacts in Sec. III D. We discuss plasmon generation in Sec. IV. In Sec. V, we compare full electromagnetic and quasi one-dimensional approaches. We draw conclusions in Sec. VI.

II. BASIC EQUATIONS
This section discusses known results that underpin a theoretical description of 2D plasmons in the system shown schematically in Fig. 1. A 2D channel occupies the plane x = 0 and is embedded in a uniform dielectric with a relative permittivity ε d . The channel is in-finitely long in the y-direction. We assume that the charge carriers are electrons; the static (dc) electron density in the channel is n 0 . A dc current is flowing in the channel, and the electron drift velocity is v 0 . We assume harmonically-varying (ac) signals with the angular frequency ω and will write all equations in terms of their amplitudes. Since our interest is in plasmons, which are TM waves, we assume that the electric field has x-and z-components, and the magnetic field only a y-component. Consequently, the ac current density only has a z-component.
We use the full system of Maxwell's equations to describe the fields in the dielectric. At the channel, the fields obey the standard field boundary conditions of the form and Here, n is the amplitude of the ac electron density, J is the amplitude of the ac current density, e is the electron charge, and superscripts (1) and (2) denote the regions above and below the channel, respectively.
The electron and current densities are related to each other through the following three equations. The first is the standard linearized expression for the current density where v is the amplitude of the ac electron velocity. The second is the linearized equation where m is the effective electron mass, and γ is the relaxation frequency (also expressed as γ = 1/τ , where τ is the relaxation time). The third one is the continuity equation Substituting Eq. (4) into Eq. (6) gives velocity, which we then substitute into Eq. (4) to obtain Equation (8) can be recast in terms of field quantities using Eqs. (2) and (3).

A. 2D plasmons
Maxwell's equations together with Eq.
where A is a constant; k is the longitudinal wavenumber (in the z-direction) and κ is the decay rate. Both κ and k are positive and satisfy the dispersion relation in the dielectric in the form where k 0 = ω/c and c is the light velocity in the dielectric. Substituting Eq.
where Ω 2 p = e 2 n 0 /(2mε 0 ε d ). In the absence of dc current and loss, Eq. (11) permits two real solutions for k, which correspond to identical counter-propagating plasmons. In the presence of a dc current, there can be up to four real-valued solutions. The full solution of the dispersion relation was discussed in Ref. [47]. However, we will assume that the drift velocity, v 0 , is low, so that electron drift provides a perturbation of the two drift-less plasmons.
The power carried by a plasmon can be found as 47 where * denotes complex conjugation. The first term in Eq. (12) corresponds to electromagnetic power, and the second one to kinetic power due to electron motion. Using Eq. (9) and ignoring loss, the expression for plasmon power can be written in the form

B. Continuum modes
Aside from plasmons, whose fields decays away from the channel, Maxwell's equations permit other eigenmode solutions whose fields do not grow away from the channel. Our interest is in TM modes that have the same symmetry as plasmons. Their y-component of the magnetic field can be written in the form 47 where k is the wavenumber in the z-direction, and q is the wavenumber in the x-direction and In contrast to plasmons, the values of k are unrestricted, 0 < k < ∞. The dispersion relation for electromagnetic waves is satisfied, so that For 0 < k ≤ k 0 , the values of q are real and correspond to radiation modes. For k > k 0 , the values of q are imaginary and correspond to evanescent modes.
In the absence of electron drift, the kinetic power of all continuum modes is zero. As expected, the electromagnetic power of the evanescent modes is also zero. As shown in Ref. [47], however, the electromagnetic power carried by the evanescent modes may be non-zero in the presence of drift. On the other hand, their total power (as defined by Eq. (12)) is zero, because the electromagnetic and the kinetic powers are carried in the opposite directions. Non-zero total power can only be carried by the radiation modes and the plasmons.

III. PLASMON SCATTERING BY INTERFACES
Having considered plasmons in uniform waveguides, we now proceed with discussing plasmons incident on an interface formed by two waveguides with different electron densities, see Fig. 1. A plasmon is incident on the interface (placed at z = 0) from the left waveguide. It will partially reflect back to the left waveguide and partially transmit through into the right waveguide. We are interested in finding the plasmon transmission and reflection coefficients. To this end, we employ two different approaches: a variational solution and full-wave simulations.
We assume a step-like variation of the dc electron density at the interface and do not consider conducting gates placed above channels, which are often used to control the electron density electrostatically. We also assume the dielectric surrounding the channels to be uniform. As a result, the model is symmetric, and the spectrum of the continuum modes, Eq. (14), is simple, allowing for a relatively compact variational formulation. On the other hand, the presence of conducting gates and non-uniform dielectrics results in a more complicated spectrum of continuum modes 21 even in the absence of dc current. This model of abrupt junctions between two ungated sections is the same as developed for graphene in Refs. [18, 19, and 24]. In graphene, the electron density in a channel can be controlled without a gate, for example, by chemical doping [48][49][50] . In GaAs/AlGaAs heterostructures, on the other hand, doping patterns can be written, with high spatial resolution, using a focused ion gun 51 . Two-dimensional channels can then be integrated with coplanar waveguides and photoconductive switches for terahertz-frequency characterization using the time-domain technique reported in Refs. [7 and 52].

A. Variational solution
Here, we present a variational solution of the interface problem. We adopt a method 29 originally developed for open dielectric waveguides supporting surface waves, in which an integral equation is formulated for the field at a waveguide interface. The equation is then solved using the variational Ritz-Galerkin procedure, by expanding the field into a series of orthonormal functions with unknown coefficients. Truncating the series, substituting it into the boundary conditions, and using the mode orthogonality conditions reduces the boundary conditions into a matrix equation for the unknown coefficients. Once these are calculated, mode reflection and transmission coefficients can be found.
In a previous study 25 , we adopted this variational method to the problem of Fig. 1 but in the absence of dc current. In this case, the mode orthogonality condition may be written in terms of solely the E x and H y field components, and it can be applied directly to the two field boundary conditions. In the presence of dc current, however, the mode orthogonality also involves ac electron current and velocity 53 . As a result, the original variational method can no longer be applied.
As we shall show, however, the variational method can be modified if one assumes that the drift velocity v 0 is low enough to neglect terms of the order v 2 0 and higher. A plasmon incident upon the interface from the left channel will excite a reflected and a transmitted plasmon and continuum modes in the both channels, so that the boundary condition for the continuity of the y-component of the magnetic field takes the form (17) Similarly, the boundary condition for the continuity of the x-component of the electric field can be written as Here, as before, R and T denote the plasmon reflection and transmission coefficients defined for the magnetic field of the form of Eq. (9); k L+ and κ L+ , k L− and κ L− are the wavenumbers of the plasmons propagating, respectively, to the left and to the right in the left channel; k R+ and κ R+ are the wavenumbers of the plasmon in the right channel; r(q) and t(q) are the amplitudes of the continuum modes in, respectively, the left and right channels.
We denote the unknown x-component of the electric field at the interface z = 0 as E(x), so that Eq. (18) can be rewritten as the following two equations and Taking Eq. (19), multiplying by exp(−κ L− x) and integrating along the x-coordinate from Using the plasmon dispersion relation in the form Ω 2 In the absence of dc current, v 0L = 0, so that Eq. (22) and the integral in Eq. (21) are zero.
Then R would be the only unknown in Eq. (21) and could have been expressed directly.
It is a consequence of the mode orthogonality, on which the original variational method relies. The next steps of the original method would be to find, by similar manipulations, expressions for T , r(q), and t(q) and substitute them into the other boundary condition for the magnetic field, Eq. (17).
In the presence of dc current, however, Eq. (21) contains two unknowns R and r(q), and is an integral equation for the latter. Neither R nor r(q) can be found from Eq. (21) and the original solution procedure cannot be applied.
However, if the electron velocity is small, we can expand r(q) into a power series with respect to the drift velocity and ignore all but the first two terms, so that is the value in the absence of dc current. Keeping only the terms that are linear in v 0L in Eq. (21) implies that r(q) in the integrand should be replaced with . If the value of r (0) (q) is known, the integral can be evaluated, and R can be expressed The values of r (0) (q) can be found, for example, by following the method introduced in Ref. [25].
The next step is to find an expression for r(q) using a similar procedure. We multiply Eq. (19) by Γ L− (q) cos(qx) + j sin(qx) and integrate ∞ 0 dx. We then ignore higher-order contributions to R and r(q) and obtain where and − denotes the Cauchy principal value. The plasmon reflection coefficient in the absence of drift R (0) can also be found using the variational method 25 ; analytical expressions are also available that are applicable in the absence of retardation 19 .
Similar manipulations with Eq. (20) give the following expressions for T and t(q) and where and the superscript (0) denotes values in the absence of dc current.
Substituting now Eqs. (23), (25), (27) and (29) into Eq. (17) and rearranging yields where and Green's function G(x, x ′ ) is of the form Equation (31) is an integral equation for the unknown field E(x), and it can be solved by the Ritz-Galerkin procedure. Choosing a complete set of basis orthogonal functions ψ n (x), the unknown field E(x) is expanded into a series of the form where λ n is a constant. Substituting (34) into (31), multiplying by ψ m (x) and integrating over x, and then truncating the series transforms the integral equation into a finite matrix equation. Following Refs. [29] and [25], we choose the basis functions as the Laguerre polynomials L n (x) weighted by an exponential, so that where x 0 is a constant.
In the absence of dc current, the quantities A, B(q), C and D(q) defined by Eqs. (24), (26), (28) and (30) where z 0 is the coordinate of the edge of a region; and γ 0 and α are constants whose values (γ 0 = 10 13 s −1 and α = 6) were chosen to create slow variation of loss along the distance.
Such a variation prevented plasmons from re-entring regions 2 and 3 while gradually absorbing their power. The field amplitudes at the ends of regions 1 and 4 were negligible, and these regions could be terminated with perfectly conducting boundaries. Since plasmon fields decay exponentially away from the 2D channels, radiation at the interface can be expected to be negligible 25 The plasmon transmission coefficient was found in a similar fashion, by extracting the field at boundary b 3 , projecting it onto the analytical expression for the field of the transmitted plasmon, H (trans) y , and normalizing, so that To demonstrate that the above approach indeed allows us to excite a single plasmon incident on the interface and then extract the reflected and transmitted plasmons without reentry, Fig. 2 also shows an example of the amplitudes of calculated total (black solid line) and reflected (orange dashed line) H y field components. The transmitted field coincides with the total field to the right of the interface. In region 2, the total field amplitude has a standing-wave pattern created by interference of the incident and reflected plasmons. Away from the interface, the transmitted and reflected field amplitudes are constant, showing propagation of a single wave in each region. However, close to the interface, field amplitudes vary significantly, indicating excitation of evanescent modes. The length of regions 2 and 3 (10 µm) was chosen specifically to allow the evanescent modes to decay. In regions 1 and 4, the field amplitudes remain initially almost constant but then quickly fall off.

C. Results
In this section, we compare results of the two models. The main parameters to explore are the ratio of the electron densities in the two channels, n 0R /n 0L , and the dc electron velocities, for which n 0L v 0L = n 0R v 0R . In all our calculations made for the interface geometry, the chosen frequency was 1 THz, the relative dielectric permittivity was 12.4, the effective electron mass was 0.067m 0 (corresponding to GaAs), and the electron density in the left channel was n 0 = 5 × 10 11 cm −2 . The electron density in the right channel varied between 0.5n 0L and 2.5n 0L . The drift velocity in the left channel varied between −3 × 10 6 and +3 × 10 6 cm/s, where negative values correspond to the dc current flowing to the left (opposite to the direction of plasmon incidence).
As a first step, to validate our two models and to calculate the zero-order transmission and reflection coefficients needed for the variational approach, we apply our models to plasmon scattering at the interface in the absence of dc current and compare their results in Fig. 3 with the model of Ref. [19]. All three models give almost identical results for the amplitude and phase of the reflection, R (0) , and the transmission, T (0) , coefficients. The nontrivial contribution to the phase of the reflection reaches around ±0.1π at the extremes of the density ranges, but the phase of the transmission coefficient is zero. The agreement between the models suggests that all three capture the essential physical mechanism of plasmon scattering in the absence of dc current. However, the expressions for plasmon reflection and transmission coefficients in Ref. [19] cannot be applied directly in the presence of a dc current in the 2D channels.
We then write the reflection coefficient in the presence of drift in the form Here, δR is the change of the reflection coefficient normalized to the ratio of the drift velocity  Figure 4 shows the values of δR calculated for different drift velocities and values of n 0R ; Fig. 4(a) is for the full-wave simulations and Fig. 4(b) is for the variational solution.
Both models result in the same behaviour of δR. Once the electron densities in the two channels start to differ, the magnitude of the δR becomes nonzero. The change of δR is more steep for n 0L > n 0R . Quantitatively, the variational solution results in slightly higher drift. This ambiguity may be removed by using power relationships. As mentioned in Sec. II B and discussed in more detail in Ref. [47], the evanescent modes may carry electromagnetic power in the presence of drift, but do not carry total power (the sum of the electromagnetic and the kinetic terms is zero). Therefore, a coefficient based on the total plasmon power, Eq. (13), appears to be an appropriate choice. We define a coefficient based on the ratio between the power flowing out of the interface, P out , to the power flowing into it, P in , as where P out is a sum of the powers of the reflected and the transmitted plasmons, and P in is the power of the incident plasmon. All powers are calculated from Eq. (12) and take into account both the electromagnetic and the kinetic terms. When G > 1, the total power flowing out of the interface exceeds the power incident upon it. In this sense, the total ac power is amplified at the interface. The value of G cannot, however, be interpreted as gain in resonators aiming to achieve plasmon oscillations (see Sec. IV). Figure 5 shows the dependence of G on the drift velocity for four values of n 0R /n 0L , calculated by both models. In all four cases, the value of G can be both larger and smaller than unity. When n 0R /n 0L < 1, G > 1 for positive values of the drift velocity (the direction of dc current coincides with the direction of incidence). On the other hand G < 1 for the opposite direction of the dc current. The situation reverses for n 0R /n 0L > 1. The effect of dc current is larger for greater differences between the electron densities in the two channels; the maximum relative increase of power is around 7%.

D. Reflection from electric contacts
Current is supplied to 2D channels through electric contacts, on which plasmons may also scatter. Using the full-wave model, we have studied how plasmons reflect from contacts in the configuration shown in Fig. 6(a). The contact has finite thickness and is assumed to be perfectly conducting. The electron density in the channel is n 0 = 5 × 10 11 cm −2 and the drift velocity is v 0 = 3 × 10 6 cm/s; its direction coincides with the direction of incidence.
We varied the contact thickness between 2. of δR are close to 4 for Λ < 0.1. Therefore, the thin contacts affect the plasmon reflection coefficient more strongly than thick ones. The behavior is the same for all three frequencies.

IV. PLASMON GENERATION
Combining plasmon amplification at interfaces with feedback may lead to plasmon generation. Feedback can be realized by multiple reflectors, such as interfaces and electric contacts. A basic system is shown in Fig. 7(a). It consists of two sections of 2D channels of different lengths and with different carrier densities. DC current is supplied through electric contacts, which we assumed to be thick. The left channel has an electron density of n 0L = 5 × 10 11 cm −2 , and the right one a density of n 0R = 2.5 × 10 11 cm −2 . The length of the left section is 1 µm.
Plasmons propagate along the two sections, and they reflect from the contacts as well as partially transmit through and partially reflect from the interface between the channels.
The acts of reflection and transmission may, depending on the relative directions of incidence and of dc current, lead to plasmon amplification or de-amplification. In addition, plasmons are losing their energy while propagating along the lossy channels. In the original study by Dyakonov and Shur as well as several subsequent studies where analytical expressions for the plasmon coefficients in terms of plasmon wavenumbers were available, the analysis of plasmon instability could been done in terms of complex frequency derived from a dispersion has then the meaning of the instability increment. However, this approach is unsuitable for our study, because the plasmon coefficients are calculated numerically assuming real-valued frequency. To demonstrate plasmon instability, therefore, we have used a somewhat different approach, adopted from analysis of threshold conditions in lasers 54 . Instability requires two conditions to be fulfilled. First, the roundtrip change of the plasmon phase should be equal to a multiple of 2π, and second, the magnitude of the roundtrip amplitude gain should exceed unity.
We have analyzed the configuration of Fig. 7(a) where R lc is the reflection coefficient of the left contact, and l L is the length of the left section.

V. COMPARISON WITH OTHER APPROACHES
The effects we have considered here are similar to those discussed in several previous theoretical studies, starting with the pioneering paper by Dyakonov and Shur 32 . Different types of instabilities have been investigated, but the common idea is to use discontinuities in 2D channels to realize plasmon amplification in the presence of a dc current, and to use reflections from the discontinuities to provide feedback needed for oscillations. Such discontinuities can be realized by electric contacts 32,33,47,55 , by interfaces between gated and ungated channels 42,44,53 , or by a variation in the channel geometry 36,41,43,56 . However, physical interpretation of the effects may differ depending on the approach chosen.
A number of studies 23,32,33,36,[41][42][43]56 have employed quasi one-dimensional approaches that can be summarised as follows. The properties of plasmons (for example, their dispersion relation) in a uniform channel are derived taking into account the variation of the relevant quantities (such the potential in the quasi-static approximation) at the both sides of the channel. However, plasmon scattering at interfaces is then treated using a one-dimensional model. Different studies have used different boundary conditions at interfaces, but a common essential feature is that these boundary conditions can be satisfied solely by the incident and the scattered plasmons.
The fully-electromagnetic approach we have adopted here differs by requiring that the field boundary conditions are satisfied across the entire plane of an interface. These boundary conditions do not reduce to the continuity of the potential and current at the channels because of the different decay rates of the incident, reflected, and transmitted plasmons.
Matching the fields across the entire interface leads to excitation of continuum modes (see Fig. 2), in addition to plasmons. One consequence, observed in the absence of drift, is a nontrivial phase of the reflected plasmons, see Fig. 3(b) and Refs. [18, 21, and 22]. This phase may not vanish for a junction between a gated and an ungated channel even in the longwavelength limit 22 . Using a one-dimensional quasi-static approach, Rejaei and Khavasi also obtained 19 non-trivial reflection phase for plasmons incident on junctions between two ungated channels. They suggested that excitation of the evanescent modes at the junction was responsible for the non-trivial phase, which agrees with the physical interpretation offered by our model. The continuum modes are also excited at the interface in the presence of drift, but as we have shown, the drift does not change the reflection phase for the configuration of Fig. 1. On the other hand, the complexity of the fully-electromagnetic approaches means that, for example for the three-dimensional structures considered in Refs. [56] and [43], a variational solution similar to the one developed in Sec. III A is unlikely to be practical, and full-wave simulations 57 would require significant computational resources.

VI. CONCLUSIONS
Using two rigorous electromagnetic approaches, we showed that plasmons may be amplified when scattering from interfaces between 2D channels in the presence of a dc current.
For moderate values of the carrier drift velocity and ratio of carrier densities in the two channels, the change of the magnitude of the plasmon reflection coefficient is proportional to the drift velocity. However, dc current does not affect considerably the phase of the reflection coefficient and both the amplitude and the phase of the transmission coefficient.
We further studied plasmon reflection from electric contacts of different thickness. The effect of dc current on the magnitude of the plasmon reflection coefficient for thin contacts was shown to be twice as large than for thick ones. Finally, we showed that conditions for plasmon generation can be met in a resonant system comprising two different channels and two electric contacts, and we quantified the effects of loss on the roundtrip gain. The results contribute to the understanding of plasmonic gain in two-dimensional systems, and may aid design of future plasmon oscillators in the terahertz range.