## Abstract

We analyze, theoretically and experimentally, the dynamics of the wavepackets in plasmonic beaming devices, and show that the beam evolution in this class of structures is a multiscale phenomenon, which initiates in the near-field proximity of the device, develops continuously over a new length scale many times the wavelength of the light, and is completed well into the far-field of the system. We develop a quantitative description of the light evolution in the beaming structures and verify our theoretical predictions with experiments. Our analytical results are utilized to develop plasmonic geometries for shaping the mid-field beam evolution, and experimental results from these structures are demonstrated.

©2011 Optical Society of America

## 1. Introduction

Recent advances in fabrication, theory, and computational techniques have enabled a large class of composite materials [1–4] where light molding is an implicitly multiscale process. Here we consider the problem of beaming in plasmonic structures whose extreme far-field properties have been described in [5–9], which link the far-field beam pattern to near-field propagation and scattering of excited surface waves. We show, experimentally and theoretically, that while beam formation in these unique devices initiates in the near-field, the beam develops continuously over a length scale many times the wavelength of the light (through the device “mid” field), and is completed well into the far-field of the system. This work allows for a new understanding of complex beam formation over this new lengthscale, opening the door to the integration of beam -steering and -shaping structures with miniaturized optical circuits, on-chip photonics, and lithographic and imaging systems.

The evolution of optical beams in simple structures can be attributed to either deep subwavelength- or extreme far-field phenomena [10]. As technology has progressed, the development of optical structures with complex geometries and material properties has necessitated an understanding of optical materials at a multiscale level. As was first discovered in Ref. [5], the diffraction of light from a subwavelength slit fundamentally changes when the slit is surrounded by a periodically corrugated plasmonic surface. In contrast to the broad, cylindrical wave, emitted by the isolated slit, emission by the plasmonic structure is dominated by two highly directional beams. This difference stems from the unique TM-polarized optical mode, guided by the interface between the dielectric and plasmonic media, known as a surface plasmon polariton (SPP) [11]. SPPs, often used for compact optical and opto-electronic structures [12–14], can be coupled to conventional, free-space optics by surface inhomogeneities. This scattering of SPPs from surface defects was experimentally demonstrated and theoretically described in earlier work demonstrating scattering from surface defects at short wavelengths [15–17]. Schematics of the plasmonic beam steering devices fabricated on GaAs substrates, and characterized in this work, are shown in Fig. 1. In a typical beam steering device such as the ones depicted in Fig. 1, SPPs, excited by a subwavelength aperture, propagate away from the slit along the metal-dielectric interface. The periodic surface corrugation results in partial out-coupling of SPPs into free-space. As result, the light emitted by the device represents the interference of the cylindrical wave emitted by the central hole and the portion of the original SPP scattered by the periodic surface corrugations [5–8]. Since the emission of the central hole lacks directionality, the direction of the beams is determined solely by the interference of SPP-scattered waves.

The direction corresponding to the interference maximum is given by the condition that, at this ‘beaming’ angle, the phase of the light scattered by two neighboring surface corrugations differs by 2*π* [6]:

*k*

_{0}and

*k*

_{SPP}correspond to the wavevector of the propagating waves above the grating and that of the SPP, respectively, and the parameter

*δ*describes the retardation of the phase of the SPP caused by an individual corrugation (see below). This relationship describes the extreme far-field response of the system, and is typically measured experimentally, with a detector positioned at a centimeter- (or larger) distance from the sample, well into the far-field of the plasmonic structure.

The situation is completely different in the relative proximity of the plasmonic structure. As evident from Fig. 2, beam evolution in this region is dominated by the prolonged high-intensity region, the “focal range” (labeled as “f” in Fig. 2), extending from the slit, across near-, mid-, and far-field range. As optical devices and systems become smaller and more compact, understanding the evolution of light on the scale of 1 to 100 wavelengths becomes critical for the correct design and development of these systems. This motivates the *necessity* to gain insight into the origin of the focal range formation in the beamsteering devices.

In this manuscript we provide the comprehensive analysis of this exciting phenomenon. The manuscript is structured as follows: Section 2 presents our main analytical results, describes the numerical approach for their validation, and summarizes the numerical studies. Section 3 describes the experimental data supporting the phenomenon of the focal range formation. In Section 4, we use the developed formalism to design and realize plasmonic-grating mid-field focusing structures. Section 5 concludes the paper.

## 2. Analytical Results and their Numerical Verification

#### 2.1. SPP transmission and scattering from corrugations

In the limit when the size of the corrugation is small when compared to the confinement scale of the SPP, the effect of an individual corrugation on the propagation of the SPP can be described by a (complex) transmission amplitude that characterizes the ratio of the amplitude of a SPP that propagates through a corrugation to the amplitude of a SPP that propagates on the uncorrugated surface (Fig.3). The value of the transmission amplitude can be estimated using the formalism described in [17]. In this work the calculation of the SPP transmission coefficients was performed using numerical solutions of Maxwell equations on commercial FEM software [20]. The transmission of the SPP through *n* = 1 to 5 inhomogeneities was calculated and compared to unconstrained propagation of the SPP; The transmission amplitude (defined as *t _{n}* = (

*A*/

_{n}*A*

_{0})

^{1/}

*), with*

^{n}*A*being the SPP amplitude at the exit side of the simulation domain upon the scattering over

_{n}*n*inhomogeneities, was calculated. Our analysis demonstrates that for the class of beaming structures used in this work, the amplitude

*t*only weakly depends on the number of scatterers, with the transmission of a SPP through an individual corrugation exceeding 90%, and relatively small reflection and phase retardation at each corrugation.

_{n}Since the effect of the corrugation is relatively weak, we can make the assumption that the scattering of SPPs by different inhomogeneities within the same system is independent of each other, and that the collective effect of the beaming structure can be considered as a sum of the effects due to individual scattering elements. At a height *y* directly above the slit, the phases of the light scattered by two neighboring corrugations differ by

*x*= Λ

_{i}*i*corresponds to the distance from the central slit to the

*i*

^{th}corrugation.

Obviously the condition Δ* _{i}* = 2

*π*corresponds to the interference maximum at the

*i*focal point of the device,

^{th}*f*:

_{i}*f*corresponding to different corrugations do not overlap with each other. As result, different focal points coalesce into a single focal range stemming from the device boundary into the far field. The extent of the focal range can be estimated as where

_{i}*i*∼ −1/ln |

_{max}*t*| with the parameter

*t*= |

*t*|exp(

*iδ*) describing the transmission of the SPP through the corrugation, as described above.

The existence of the focal zone was first noticed in numerical simulations in Ref. [7]. However, the phenomenon was suggested to appear only in the limited frequency range under normal beaming conditions (*θ* ≃ 0). One of the main results of this work is the demonstration (both analytically and experimentally) of the fact that light focusing is a universal phenomenon, appearing in the majority of plasmonic gratings and other grating-outcoupled waveguides [9, 22].

In all these structures, the end of the focal region corresponds to the transition between the few-corrugation interference (dominating the near- to mid-field zone of the device) and the collective interference dominating the extreme far-field behavior. In particular, this transition results in the formation of the beams, which are seen to be originating in the proximity to *f*
_{max}. For shorter wavelengths *f _{i}* < 0, indicating that the focal spot is “virtual”, so that the beams appear to be emanating from behind the corrugated surface.

As pointed out above, the beams themselves are the result of constructive interference of light scattered by the corrugated plasmonic surface. However, depending on the geometry of the device [11], this light can be in- or out-of-phase with the light emitted by the central slit. The phase relation between the slit- and surface-originated beams determines whether the beams appear at maxima or minima in the far-field response.

In the mid-field, however, the phase difference between the light emitted from the central slit and the light scattered by corrugations strongly depends on the coordinate *y*. Thus, the focal range always manifests itself as set of interference maxima.

#### 2.2. Calculation of field distribution in the system

In order to gain insight into the wave propagation in multiscale plasmonic structures and to verify the validity of the analytical equations describing the far- (Eq. 1) and mid- (Eq. 3) field beam profiles we calculate the intensity distribution in the plasmonic beaming structures via numerical solutions of Maxwell equations. Two alternative techniques have been used to calculate the field distribution in the system. First, the beam evolution has been calculated via commercially-available finite-element-method (FEM) PDE solver [20]. The simulation space in our models was at least 300*μm*-long and at least 450*μm*-high. The maximum size of the mesh was 300*nm*, corresponding to approximately
${\lambda}_{0}/10\sqrt{{\varepsilon}_{\text{GaAs}}}$, ensuring the convergence of the field profile. The simulation space was terminated by perfectly matched layers to adequately emulate the unbounded space in experiment. This class of simulations provides us with capability to analyze real-space intensity distribution, validating Eq.(3), and it was primarily used in this work to analyze the mid-field dynamics of model- and experimentally-realized structures.

Alternatively, to assess the interplay between real-space and wavevector-space dynamics and to verify extreme far-field behavior of the structures [Eq.(1)], the field profile was calculated using the modal expansion technique outlined in Ref. [21]. In these semi-analytical calculations, we have the metal film was considered to be relatively thick, so that the subwavelength slit was represented by a subwavelength planar waveguide with highly-conductive walls. The field inside this guide was represented as a set of waveguide modes, with the incident field being in *TM*
_{0} mode propagating upwards. At the slit-GaAs boundary the incident field was partially reflected, partially transmitted into the field in GaAs layer. The reflected field can be analytically written as

The field inside the GaAs layer was represented as a set of plane waves propagating away from the slit. Explicitly, this field is represented by

*w*are selected to ensure that the sum adequately approximates the improper integral.

_{n}The unknown modal amplitudes are then found using boundary conditions.

*ε*(

_{b}*x*) describing the permittivity of a substrate at a particular location.

The above equations are solved for
$a({k}_{x}^{n})$ with a technique introduced in Ref. [21], specifically developed to calculate the modal amplitudes when the number of modes in the region *x* < 0 differs from that in the region *x* > 0.

Since the mode-matching formalism assumes substrate to be flat (although having position-dependent permittivity), the developed technique is accurate in the absence of resonant-cavity-like modes [6] in the structure (low-scattering approximation implicitly assumes the same regime).

To verify the validity of the mode-matching technique we compare the intensity distribution calculated with this method against the distribution generated by FEM simulation. As an example, we choose the system with deep (*h* = −0.5*μm*), *L* = 1.4*μm*-long inserts of dielectric with *ε* = 10.89(1 + *i*) (see Fig. 2 for notations). The agreement between two fundamentally different types of numerical solutions of Maxwell equations is shown in Fig. 5. It is clearly seen that the mode-matching-based solutions adequately represent the results of FEM-based simulations, with exception of some minor interference-like fringing. The fringing is a common artifact originating from discrete Fourier expansion used as the basis for the field representation in Eq. 6. Since the field pattern produced by the discrete Fourier spectrum is always periodic (with period given by the smallest spacing in the spectrum), the real solution is effectively accompanied by light generated by infinite number of virtual copies of the system. The fringing is produced when light emitted by the central structure interferes with light from these virtual copies. The fringing can be significantly reduced when the signal from
${k}_{x}\hspace{0.17em}\sim \hspace{0.17em}\sqrt{{\varepsilon}_{\text{GaAs}}\omega /c}$, propagating parallel to the metal interface, is eliminated from the plot.

#### 2.3. Numerical Verification of Analytical Results

In this work, we select the operating wavelength in the mid-IR range (*λ*
_{0} ∼ 10*μm*) that is of crucial importance for sensing and security applications [8]. Another advantage of working in the mid-IR frequency range are the extremely long propagation length of SPPs, which allow a systematic study of focal range dependence on the number of corrugations. To adequately compare the response of several structures, we fix the period of corrugations Λ = 2.8*μm*, thereby fixing the extreme far-field direction of the beams emanating from all these structures. To compare our analytical and numerical results to our experimental results, we also assume that the formation of the beams takes place in GaAs (*ε* = 10.89) [18], that serves as a substrate in our experimental device (see below).

We first consider the response of the structure with deep dielectric inserts, discussed above. Note that due to the relatively weak confinement of the mid-IR SPPs, this geometry, despite the deep grooves, still exhibits a significant (∼ 90%) transmission of the SPP wave across each groove (Fig.3). The profile of the field across the system, calculated with FEM techniques, and its comparison to the predictions of Eqs.(1, 4) are shown in Fig.6. It is clearly seen that the developed analytical expressions adequately describe the response of the system.

Alternatively, the mode-matching calculations yield wavevector-space representation of the field that unambiguously demonstrates the dominant role of the interference of light scattered by the SPPs in the optical response of the system.

To demonstrate universality of the observed phenomenon, we simulate the formation of the beams in a similar system, where the relatively deep grooves are replaced with *δ*-like scatterers, modeled as 250nm-high and 100-nm-thick bumps. In contrast to the previously described structure, about 99% of SPP energy is transmitted through the individual *δ*-scatterer, so the *collective* effect of the SPP bumps on the formation of the beams can be adequately evaluated. Fig.7 clearly reveals that the size of the focal zone is determined by the number of surface inhomogeneities contributing to the beam formation.

## 3. Experimental Set-Up and Results

Mid-IR plasmonic beam steering structures were fabricated in order to experimentally verify the analytical and numerical simulations described above. The devices were fabricated using standard optical lithography and metallization techniques on semi-insulating GaAs wafers. The first set of samples were fabricated to support mid-IR SPPs on the metal/GaAs interface [Fig. 1(a)]. For these samples, the pattern for two periodic groove arrays, separated by 4.2*μm* were exposed by photolithography in photoresist, which then served as an etch mask to define the period grooves (Λ = 2.8*μm*), etched a depth d=0.5*μm* into the GaAs surface. A second photoresist pattern, consisting of a single, 1.4*μm* slit, centered between the two groove arrays, was then exposed. A metal layer consisting of 10nm Ti and 200nm Au was deposited, and the central slit lifted off in acetone. The thin Ti layer is used to ensure adhesion of the metal layer to the semiconductor surface, and does not significantly change the metal layer optical properties, as both Ti and Au, in the mid-IR wavelength range, have large and negative values of permittivity [19]. The resulting structure consisted of the 1.4*μm* central slit, flanked by the two groove arrays of period 2.8*μm*. The GaAs wafer was then lapped down to a thickness of ∼ 200*μm*, and the back surface roughened, to allow field measurements as close to the slit as possible and to prevent multiple reflections within the GaAs substrate, respectively (good agreement between experimental results and theoretical calculations indicate that such roughening does not introduce major scattering in the system). A second set of samples was also fabricated, in this case to support mid-IR SPPs on the air/metal interface [Fig. 1(b)]. Because the period of the patterned corrugations for our structures scale as 1/*ε _{d}* (where

*ε*is the permittivity of the dielectric at the metal/dielectric interface), the air/metal structures afford us greater flexibility in geometry, allowing for the fabrication of smaller features (relative to

_{d}*λ*) than the metal/GaAs structures. For our metal/air structures, a narrow ridge of photoresist is patterned on the GaAs substrate, the surface metallized, and the slit lifted off. A second layer of metal is then deposited over the first, consisting of the two periodic ridge arrays flanking the central slit. The air/metal structures developed for this work, also allow a probing of the beam fields close to the plasmonic surface (relative to the metal/GaAs structures) and were utilized for the focusing samples to be discussed in the following section.

_{SPP}Two experimental set-ups were developed to characterize both the mid- and far-field beam characteristics of our structures. First, light transmitted through the beaming structures is captured by a lens pair (collection optics) and focused onto a HgCdTe detector. For this set-up, both the collection optics and detector are mounted on a rotational stage, set to rotate about the central slit of our structures (and the focal point of the beam incident upon the sample), giving a picture of the far-field beam intensity of the light transmitted through the structure. We compare the prediction of the developed analytical and numerical simulation techniques to the experimental results. As expected, device performance was strongly polarization-dependent. Thus, the transmission of TE-polarized light, which does not couple to SPPs at the central slit, was negligibly small, while the transmission of TM-polarized light was significantly larger and its far-field directionality (measured at ∼ 5*cm*) was in perfect agreement with predictions of the generalized diffraction theory explained above and represented by Eq.1. The far-field imaging set-up and beaming patterns can be seen in our previous work [8].

A second experimental set-up, shown in Fig. 8(a), was developed to investigate the mid- to far- field beam evolution in our plasmonic steering structures. Light from a Daylight Solutions external cavity tunable quantum cascade laser (tunable from *λ* = 8.8 – 9.4*μm*) is focused onto the back side of the sample, after passing through a polarizer pair (to allow us the switch between TE and TM excitation of the structure). A narrow (75*μm*-wide) vertical slit was mounted on a computer-controlled moving stage behind the sample. The light transmitted through the slit is collected by a ZnSe lens pair and focused onto a liquid-nitrogen cooled HgCdTe detector with 1*mm* × 1*mm* element size. Finally a beam pick-off is placed before the focusing lens and a small fraction of the incident laser light is directed to a second HgCdTe detector. The second detector allows us to normalize the transmitted intensity to the incident laser power, which prevents fluctuations in laser power over the duration of the scan from distorting the collected field intensity data. The dependence of the mid/far-field intensity transmitted through the vertical slit, as a function of wavelength and the position of the slit, was analyzed for each of the structures fabricated [Fig.8(b,c)]. The small scan area, coupled with the small f-number of our collection lens (*f* /1.5) and the relatively large detector element, allows us to maintain relatively consistent collection efficiency across our scan area.

Field intensity data was collected on the backside (from the GaAs substrate surface) of our steering structures. The measurements demonstrate the existence of the prolonged focal range behind the sample [Fig. 8(b)], with the double-beam forming behind the focal range, in complete agreement with the theory described above. Note that in these last calculations, Eqs.(2,4) were modified to explicitly account for the finite thickness of the GaAs substrate.

## 4. Manipulating the Mid-Field

Having demonstrated both analytical and numerical methods for modeling the beam evolution in our beam steering structures, the concept of beam interference can now be used to reshape the profile of the light in the mid and far-field of the device for various applications. In Fig. 9 we illustrate this process by refocusing the radiation to a particular spot above the corrugated surface.

The single focal spots in Fig. 9 appear when light scattered by all corrugations collectively interferes at a pre-designed point in space. Obviously, this condition is not achievable with a periodically corrugated structure (see discussion after Eq. 3). Starting from Eq. 2, and requiring that Δ = 2*π*, we obtain the following expression describing the position of the next corrugation based on the position of the current corrugation:

*f*

_{0}being the desired distance from the slit to the focal spot, and

*κ*=

_{i}*k*

_{SPP}+

*k*

_{0}

*x*/

_{i}*f*

_{0}.

Samples were fabricated to demonstrate the focusing ability of the non-periodic plasmonic structures. For the focusing structures, the light is incident upon the GaAs substrate, and SPPs are generated on the patterned air/metal interface. Multiple structures were designed and fabricated, intended to focus *λ* = 9*μm* light to points *y* = 250*μm* and *y* = 500*μm* above the structures’ central slit. Fig. 9 shows both the numerical simulations and the experimental results for the focusing structures designed to focus to points *y* = 250*μm* and *y* = 500*μm* above the central slit. Clear focal points are seen in the experimental data at *y* = 250*μm*, shown in Fig. 9(c) and at *y* = 500*μm*, shown in Fig. 9(f), as expected for the designed structures. It should be noted that our experiment, being a far-field measurement technique, cannot reproduce the high-intensity regions close the metal/dielectric surface that correspond to evanescent and SPP components. Thus, all experimental results will only allow for imaging of the propagating component of the fields in the dielectric, which explains the only significant deviation of our experiment from the theoretical models. The results from the focusing samples clearly indicate the potential for utilizing patterned plasmonic interfaces for control of the mid-field intensity for a variety of applications in sensing, lithography, and on-chip photonic technologies.

## 5. Conclusions

To conclude, we have demonstrated, experimentally, analytically, and numerically, that beam formation in plasmonic systems is an inherently multiscale phenomenon and developed a quantitative description of this effect. Although the present work focused on mid-IR focal range, the developed formalism is directly applicable to visible, near-IR, and and especially THz frequencies (where surface waves also propagate for long distances). The formalism can be further extended to predict the beam formation in new types of quantum cascade lasers [9], and in other structures where the emission pattern is formed by partial scattering of guided modes. A wide class of materials and devices where similar phenomena are expected to take place includes photonic crystals, and other metamaterials- as well as composite structures with typical inhomogeneity lengths scales comparable to wavelength.

## Acknowledgments

This work was supported by ONR ( # N00014-07-1-0457), NSF ( #ECCS-0724763), AFRL ( #FA8650-08-C-1445), and the AFOSR Young Investigator Program ( #FA9550-10-1-0226).

## References and links

**1. **J. D. Joannopoulos, R. D. Meade, and J. N. Winn, *Photonic Crystals: Molding the flow of light* (Princeton U. Press, 1995).

**2. **D. R. Smith, W. J. Padilla, D. C. Vier, S. C. Nemat-Nasser, and S. Schultz, “Composite medium with simultaneously negative permeability and permittivity,” Phys. Rev. Lett. **84**, 4184–4187 (2000). [CrossRef] [PubMed]

**3. **U. Leonhardt, “Optical conformal mapping,” Science **312**, 1777–1780 (2006). [CrossRef] [PubMed]

**4. **W. Cai, U. K. Chettiar, A. V. Kildishev, and V. M. Shalaev, “Optical cloaking with metamaterials,” Nature Photon. **1**, 224–227 (2007). [CrossRef]

**5. **H. J. Lezec, A. Degiron, E. Devaux, R. A. Linke, L. Martín-Moreno, F. J. García-Vidal, and T. W. Ebbesen, “Beaming light from a subwavelength aperture,” Science **297**, 820–822 (2002). [CrossRef] [PubMed]

**6. **L. Martín-Moreno, F. J. García-Vidal, H. J. Lezec, A. Degiron, and T. W. Ebbesen, “Theory of highly directional emission from a single subwavelength aperture surrounded by surface corrugations,” Phys. Rev. Lett. **90**, 167401 (2003). [CrossRef] [PubMed]

**7. **F. J. García-Vidal, L. Martín-Moreno, H. J. Lezec, and T. W. Ebbesen, “Focusing light with a single subwavelength aperture flanked by surface corrugations,” Appl. Phys. Lett. **83**, 4500–4502 (2003). [CrossRef]

**8. **D. C. Adams, S. Thongrattanasiri, T. Ribaudo, V. A. Podolskiy, and D. Wasserman, “Plasmonic mid-infrared beam steering,” Appl. Phys. Lett. **96**, 201112 (2010).

**9. **N. Yu, J. Fan, Q. J. Wang, C. Pflügl, L. Diehl, T. Edamura, M. Yamanishi, H. Kan, and F. Capasso, “Small-divergence semiconductor lasers by plasmonic collimation,” Nature Photon. **2**, 564–570 (2008). [CrossRef]

**10. **M. Born and E. Wolf, *Principles of Optics* (Cambridge Univ., 1999).

**11. **H. Raether, *Surface Plasmons on Smooth and Rough Surfaces and on Gratings* (Springer-Verlag, 1988).

**12. **S. I. Bozhevolnyi, V. S. Volkov, E. Devaux, J.-Y. Laluet, and T. W. Ebbesen, “Channel plasmon subwavelength waveguide components including interferometers and ring resonators,” Nature **440**, 508–511 (2006). [CrossRef] [PubMed]

**13. **W. L. Barnes, A. Dereux, and T. W. Ebbesen, “Surface plasmon subwavelength optics,” Nature **424**, 824–830 (2003). [CrossRef] [PubMed]

**14. **A. Krier, *Mid-infrared Semiconductor Optoelectronics* (Springer, 2006). [CrossRef]

**15. **I. I. Smolyaninov, D. L. Mazzoni, and C. D. Davis, “Imaging of surface plasmon scattering by lithographically created individual surface defects,” Phys. Rev. Lett. **77**, 3877–3880 (2003). [CrossRef]

**16. **I. I. Smolyaninov, D. L. Mazzoni, J. Mait, and T. W. Ebbesen, “Experimental study of surface-plasmon scattering by individual surface defects,” Phys. Rev. B **56**, 1601–1611 (2003). [CrossRef]

**17. **A. V. Shchegrov, I. V. Novikov, and A. A. Maradudin, “Scattering of surface plasmon polaritons by a circularly symmetric surface defect,” Phys. Rev. Lett. **78**, 4269–4272 (2003). [CrossRef]

**18. **E.D. Palik, *The Handbook of Optical Constants of Solids* (Academic, 1997).

**19. **M. A. Ordal, L. L. Long, R. J. Bell, R. R. Bell, R. W. Alexander Jr., and C. A. Ward, “Optical properties of the metals Al, Co, Cu, Au, Fe, Pb, Ni, Pd, Pt, Ag, Ti, and W in the infrared and far infrared,” Appl. Optics **22**, 1099–1119 (1983). [CrossRef]

**20. **For details see COMSOL Multiphysics 3.5a User’s Guide and RF Module User’s Guide, www.comsol.com.

**21. **S. Thongrattanasiri, J. Elser, and V. A. Podolskiy, “Quasi-planar optics: computing light propagation and scattering in planar waveguide arrays,” J. Opt. Soc. Am. B **26**, B102–B110 (2009). [CrossRef]

**22. **M. I. Amanti, M. Fischer, G. Scalari, M. Beck, and J. Faist, “Low-divergence single-mode terahertz quantum cascade lase,” Nature Photon. **3**, 586–590 (2009) [CrossRef]