Using two rigorous electromagnetic approaches, we study plasmon scattering in two-dimensional 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. A resonant system is studied where plasmon roundtrip gain may exceed unity, showing the possibility of plasmon generation.

## I. INTRODUCTION

Two-dimensional materials (such as graphene^{1–4}) and systems (such as GaAs/AlGaAs heterostructures^{5–7}) are able to support plasmons. These are waves that exist due to the collective motion of free charges, and they may propagate in the microwave,^{6,8,9} terahertz,^{1,5,10,11} and infrared frequency^{3,4} ranges. Initial theoretical studies have concentrated on dispersion relations and resonance frequencies of gated and ungated two-dimensional (2D) plasmons,^{12–14} magnetoplasmons,^{15} and edge (magneto-)plasmons.^{16,17} There has recently been an upsurge of interest in plasmon scattering at discontinuities, such as those formed by interfaces between 2D systems with two different carrier densities,^{18–20} between a gated and an ungated system,^{21–23} or by edges^{18,24} and electric contacts.^{25} Plasmon scattering has been studied by adopting methods developed previously for the traditional metallic and dielectric waveguides, such as the mode-matching^{26,27} and the variational methods,^{28,29} transmission-line models,^{28} and the Wiener–Hopf technique.^{26} Two-dimensional plasmons in the absence of dc currents have also been studied extensively experimentally, traditionally in semiconductor heterostructures (e.g., Refs. 5–7, 10, 11, and 30) and more recently in novel 2D materials (e.g., Refs. 1–4, and 31).

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, and 39); evidence of plasmon amplification from interface scattering has come to date from observing radiation from double-grating-gate 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,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 the variational solution of an integral equation for the field at the interface (Sec. III A) and a full-wave numerical model based on the solution of coupled Maxwell’s and hydrodynamic equations (Sec. III B). We compare the 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 $\epsilon d$. The channel is infinitely long in the $y$-direction. We assume that the charge carriers are electrons; the static (dc) electron density in the channel is $n0$. A dc current is flowing in the channel, and the electron drift velocity is $v0$. We assume harmonically varying (ac) signals with the angular frequency $\omega $ 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 of motion (Euler’s equation)

where $m$ is the effective electron mass and $\gamma $ is the relaxation frequency (also expressed as $\gamma =1/\tau $, where $\tau $ is the relaxation time). The third one is the continuity equation

We substitute Eq. (7) in the equation of motion (5) to find an expression for the ac velocity, which we then substitute into Eq. (4) to obtain

### A. 2D plasmons

Maxwell’s equations together with Eq. (8) and the boundary conditions Eqs. (1)–(3) permit eigenmode solutions known as plasmons. These are surface TM waves whose magnetic field component for $x\u22650$ can be written in the form

where $A$ is a constant, $k$ is the longitudinal wavenumber (in the $z$-direction), and $\kappa $ is the decay rate. Both $\kappa $ and $k$ are positive and satisfy the dispersion relation in the dielectric in the form

where $k0=\omega /c$ and $c$ is the light velocity in the dielectric. Substituting Eq. (9) into Maxwell’s equations and Eqs. (1)–(3) and (8) yields a plasmon dispersion relation of the form

where $\Omega p2=e2n0/(2m\epsilon 0\epsilon 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, $v0$, is low so that electron drift provides a perturbation of the two driftless plasmons.

The power carried by a plasmon can be found as^{47}

where $\u2217$ 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 decay 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<\u221e$. The dispersion relation for electromagnetic waves is satisfied, so that

For $0<k\u2264k0$, the values of $q$ are real and correspond to radiation modes. For $k>k0$, 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–50} while in GaAs/AlGaAs heterostructures, lateral patterns of doping with high spatial resolution can be achieved using focused ion beams during molecular beam epitaxy.^{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 $Ex$ and $Hy$ 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 $v0$ is low enough to neglect terms of the order $v02$ 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

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); $kL+$ and $\kappa L+$, $kL\u2212$ and $\kappa L\u2212$ are the wavenumbers of the plasmons propagating, respectively, to the left and to the right in the left channel; $kR+$ and $\kappa 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(\u2212\kappa L\u2212x)$ and integrating along the $x$-coordinate from $x=0$ to $x=\u221e$ yields

Using the plasmon dispersion relation in the form $\Omega pL2\kappa L\u2212=(\omega \u2212kL\u2212v0L)2$, we can write

In the absence of dc current, $v0L=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 $r(q)\u2248r(0)(q)+\delta r(q)\u22c5v0L$, where $r(0)(q)$ is the value in the absence of dc current. Keeping only the terms that are linear in $v0L$ in Eq. (21) implies that $r(q)$ in the integrand should be replaced with $r(0)(q)$. If the value of $r(0)(q)$ is known, the integral can be evaluated, and $R$ can be expressed as

where

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 $\Gamma L\u2212(q~)cos(q~x)+jsin(q~x)$ and integrate $\u222b0\u221edx$. We then ignore higher-order contributions to $R$ and $r(q)$ and obtain

where

and $\u222b\u2212$ 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)$:

where

and

where

and the superscript $(0)$ denotes values in the absence of dc current.

where

and Green’s function $G(x,x\u2032)$ 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 $\psi n(x)$, the unknown field $E(x)$ is expanded into a series of the form

where $\lambda n$ is a constant. Substituting (34) into (31), multiplying by $\psi m(x)$ and integrating over $x$, and then truncating the series transforms the integral equation into a finite matrix equation. Following Refs. 25 and 29, we choose the basis functions as the Laguerre polynomials $Ln(x)$ weighted by an exponential, so that

where $x0$ 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) are equal to zero. They have non-zero values in the presence of dc current because the functions defining the field profiles of the modes in the waveguides are not orthogonal to each other. Equations (24) and (28), the first term in Eq. (30), and the first two terms in Eq. (28) describe contributions due to non-orthogonality between the fields of plasmons and of the continuum modes. On the other hand, the integrals in Eqs. (26) and (30) describe contributions due to non-orthogonality of the fields of different continuum modes. The latter contributions can be expected to be small compared to those involving plasmons, and we neglected these contributions to simplify further calculations.

### B. Full-wave simulations

The other method we developed to study plasmon scattering at interfaces of Fig. 1 is full-wave simulations. Maxwell’s equations were solved using the finite-element method in a commercial solver. Two-dimensional channels were modeled using the boundary condition equations (2) and (3) combined with the equation for the linearized current density equation (8). In the presence of dc current, the wavenumbers and field distributions of plasmons in a 2D channel are different for the opposite propagation directions so that the channel becomes a nonreciprocal waveguide. As a result, we could not rely on such built-in solver techniques as perfectly matched layers, and computation of eigenmode fields and of scattering parameters. We had to design a bespoke computational domain, a method of plasmon excitation, and a method to extract plasmon transmission and reflection coefficients.

The computational domain consisted of four regions as shown in Fig. 2. Regions 2 and 3 were 10 $\mu $m-long, lossless 2D channels with different electron densities. Plasmon was excited at boundary $b1$ and traveled to the right along region 2, reaching the interface between the 2D channels at boundary $b2$. The plasmon then partially reflected back into region 2 and partially transmitted through the interface into region 3. The energy of the reflected and transmitted plasmons had then to be guided away from regions 2 and 3. It was accomplished by the absorbing terminations in regions 1 and 4, each of which was 20 $\mu $m long. These regions were identical to their neighbors except for the loss in the 2D channels. The relaxation frequency [see Eq. (5)] varied along regions 1 and 4 as

where $z0$ is the coordinate of the edge of a region; and $\gamma 0$ and $\alpha $ are constants whose values ($\gamma 0=1013s\u22121$ and $\alpha =6$) were chosen to create the slow variation of loss along the distance. Such a variation prevented plasmons from re-entering 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 top and bottom of the computation domain could then be also terminated by perfectly conducting boundaries placed each 100 $\mu $m away from the 2D channels.

To excite a plasmon, we first solved the plasmon dispersion relation in the left 2D channel, taking electron drift into account, and then calculated, using an analytical expression, the corresponding amplitude of the $Hy$ field component. We then imposed this field distribution at boundary b$1$. The resulting plasmon traveled to the right (into region 2).

To calculate the plasmon reflection coefficient, we extracted the $Hy(x)$ field component immediately to the left of boundary b$1$, projected in onto the analytical expression for the corresponding field component of the incident plasmon, $Hy(inc)$, and normalized the result, leading to

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, $Hy(trans)$, 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) $Hy$ 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 $\mu $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 the results of the two models. The main parameters to explore are the ratio of the electron densities in the two channels, $n0R/n0L$, and the dc electron velocities, for which $n0Lv0L=n0Rv0R$. 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.067m0$ (corresponding to GaAs), and the electron density in the left channel was $n0=5\xd71011cm\u22122$. The electron density in the right channel varied between $0.5n0L$ and $2.5n0L$. The drift velocity in the left channel varied between $\u22123\xd7106$ and $+3\xd7106cm/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 $\xb10.1\pi $ 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, $\delta R$ is the change of the reflection coefficient normalized to the ratio of the drift velocity to the plasmon phase velocity in the absence of drift ($\omega /kL(0)$). Once $R$ and $R(0)$ have been found from the full-wave simulations and the variational method, Eq. (39) can be used to find $\delta R$ for a particular value of the drift velocity. Both models show that the presence of drift does not appreciably change the phase of the reflection coefficient, and $\delta R$ is real for the whole range of the drift velocities used. Figure 4 shows the values of $\delta R$ calculated for different drift velocities and values of $n0R$; Fig. 4(a) is for the full-wave simulations and Fig. 4(b) is for the variational solution.

Both models result in the same behavior of $\delta R$. Once the electron densities in the two channels start to differ, the magnitude of the $\delta R$ becomes nonzero. The change of $\delta R$ is more steep for $n0L>n0R$. Quantitatively, the variational solution results in slightly higher values of $\delta R$ than the full-wave simulations; for example, at $n0R=2.5n0L$, the full-wave simulations give $\delta R\u2248\u22120.35$, whereas the variational model gives $\delta R\u2248\u22120.5$. For most of the range of $n0R$ used, the curves calculated for different drift velocities lie close to each other for both models, which confirms the assumption that the absolute change of the reflection coefficient $R\u2212R(0)$ is proportional to the drift velocity. However, at the lower end of the $n0R$ range, the curves for different drift velocities deviate from each other more noticeably. Because the variational model relies on the assumption that the change of the reflection and transmission coefficients is proportional to the drift velocity, it cannot be expected to yield quantitatively correct values of $R$ for values of $n0R$ lower than the minimum value chosen. The full-wave simulations, however, are free from this assumption. Both models also show that the absolute change of the plasmon transmission coefficient, $|T\u2212T(0)|$ is an order of magnitude lower than the corresponding values of $|R\u2212R(0)|$, which suggests that dc current does not appreciably affect the transmission coefficient.

We have defined the plasmon reflection and transmission coefficients based on the magnetic field component [see Eq. (17)]. However, the magnitude of the coefficients may be different for a different choice of the field component, because the wavenumbers of the counter-propagating plasmons even in the same waveguide are different in the presence of 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, $Pout$, to the power flowing into it, $Pin$, as

where $Pout$ is a sum of the powers of the reflected and the transmitted plasmons, and $Pin$ 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 $n0R/n0L$, calculated by both models. In all four cases, the value of $G$ can be both larger and smaller than unity. When $n0R/n0L<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 $n0R/n0L>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 a finite thickness and is assumed to be perfectly conducting. The electron density in the channel is $n0=5\xd71011cm\u22122$ and the drift velocity is $v0=3\xd7106cm/s$; its direction coincides with the direction of incidence. We varied the contact thickness between 2.5 nm and 200 $\mu $m, and calculated the plasmon reflection coefficients, in the presence and absence of drift, for three values of frequency, 0.5, 0.75, and 1 THz. Figure 6(b) shows how the normalized change of the reflection coefficient $\delta R$ [see Eq. (39)] depends on a normalized contact thickness defined as $\Lambda =hk(0)$, where $k(0)$ is the plasmon wavenumber in the absence of dc current. With an exception of a small dip, $\delta R$ has values around $2$ already for $\Lambda >1$. The value of $\delta R=2$ agrees with the analytical result obtained previously for an infinitely thick contact.^{47} On the other hand, the values of $\delta R$ are close to 4 for $\Lambda <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 $n0L=5\xd71011cm\u22122$, and the right one a density of $n0R=2.5\xd71011cm\u22122$. The length of the left section is 1 $\mu $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 be done in terms of complex frequency derived from a dispersion relation. Assuming the time variation of the form $exp(i\omega t)$, a negative value of the imaginary part of the frequency implies waves that grow in time. The imaginary part of the frequency 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\pi $, and second, the magnitude of the roundtrip amplitude gain should exceed unity.

We have analyzed the configuration of Fig. 7(a) as follows. First, using full-wave simulations, we found all reflection and transmission coefficients at a frequency of 1 THz and several values of the drift velocity, $v0L=0,\xb115,\xb130\xd7106cm/s$. Positive values of the drift velocity correspond to the current flowing to the right. When a plasmon is incident on the interface between the section from the left, it will partially reflect back and partially transmit through the interface. The transmitted plasmon will travel along the right section, reflect from the right contact, and travel back to the interface, where it will undergo further partial reflection and transmission. The total reflection coefficient from the interface, $R$, can then be found using the sum-of-all-paths method.^{22,44,53} The roundtrip gain in the left section is then given by the following expression:

where $Rlc$ is the reflection coefficient of the left contact, and $lL$ is the length of the left section. From Eq. (41), we have determined whether the conditions for plasmon oscillations can be fulfilled. By varying the length of the right section, we changed the argument of $R$ and, therefore, the plasmon roundtrip phase. We then found those values of the right-section length that fulfilled the roundtrip phase condition and calculated the magnitude of the roundtrip gain to see whether it can exceed unity. Figures 7(b)–7(d) show the magnitude of the roundtrip gain (circles) calculated for different values of the relaxation time. Figure 7(b) corresponds to lossless channels. When no dc current flows, the roundtrip gain (hollow circles) is equal to unity, corresponding to a lossless plasmon resonance. When the current is flowing from the left to the right contacts, the roundtrip gain is less than unity, indicating a damped resonance. However, for the opposite direction of current, the roundtrip gain exceeds unity, indicating instability and plasmon generation. The resonant lengths of the left section depend on the drift velocity because dc currents affect plasmon wavenumbers [see Eq. (11)].

As shown in Figs. 7(c) and 7(d), loss in the channels suppresses gain, and this effect is stronger for longer sections. For $\tau =50$ ps, the roundtrip gain still exceeds unity for all resonances and the highest drift velocity of $v0L=\u221230\xd7106cm/s$. For the smaller drift velocity of $v0L=\u221215\xd7106cm/s$, the roundtrip gain exceeds unity only for the first resonance. For $\tau =10$ ps, the roundtrip gain is less than unity for all resonances.

## 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 the 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–43,56} have employed quasi-one-dimensional approaches that can be summarized 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 as the potential in the quasi-static approximation) at 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 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 the excitation of continuum modes (see Fig. 2), in addition to plasmons. One consequence, observed in the absence of drift, is a non-trivial 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 long-wavelength limit.^{22} Using a one-dimensional quasi-static approach, Rejaei and Khavasi also obtained^{19} the non-trivial reflection phase for plasmons incident on junctions between two ungated channels. They suggested that the 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. 43 and 56, 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 the 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, the 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 thicknesses. The effect of dc current on the magnitude of the plasmon reflection coefficient for thin contacts was shown to be twice as large as that 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 the design of future plasmon oscillators in the terahertz range.

## ACKNOWLEDGMENTS

The authors gratefully acknowledge financial support by the EPSRC (Grant Nos. EP/R004994/1 and EP/R00501X/1).

## DATA AVAILABILITY

The data that support the findings of this study are openly available in Imperial College London Data Research Depository at http://doi.org/10.14469/hpc/8088, Ref. 58.