The homogenization limit and waveguide gradient index devices demonstrated through direct visualization of THz fields

Electromagnetic homogenization approximation calculates an effective refractive index of a composite material as a weighted average of its components, and has found uses in gradient refractive index and transformation optics devices. However, the utility of the homogenization approximation is hindered by uncertainty in its range of applicability. Harnessing the capability of time-resolved imaging provided by the terahertz polaritonics platform, we determined the dispersion curves of slab waveguides with periodic arrays of holes, and we quantified the breakdown of the homogenization approximation as the period approached the terahertz wavelength and the structure approached the photonic bandgap regime. We found that if the propagation wavelength in the dielectric waveguide was at least two times as large as the Bragg condition wavelength, the homogenization approximation held independent of the detailed geometry, propagation direction, or fill fraction. This value is much less demanding than the estimate of 10:1 often assumed for homogenization. We further used the experimental capabilities to extract the effective refractive index of the photonic crystals in the homogenization approximation limit, and we used this to analyze the predictive strength of analytical formulas. These formulas enabled rapid design of a Luneburg lens and a bi-directional cloak in a waveguide platform without the need for numerical simulations. Movies of terahertz waves interacting with these structures, which were fabricated using femtosecond laser machining, reveal excellent performance. The combination of an analytical formula and confidence in the homogenization approximation will aid in fast design and prototyping of gradient index devices.


Introduction
When appropriate conditions are met, the homogenization approximation can simplify a complex optical system comprised of multiple elements, such as an array of holes or inclusions in a host dielectric, by approximating it as a homogenous composite material with an effective refractive index [1][2][3][4]. The effective index can be continuously tuned by selecting the 'concentration' of the inclusion material. By controlling the spatially varying size, density, or type of inclusions, one can specify the three-dimensional (3D) index of refraction, n x y z ( , , ). This approach has found widespread application in the design of graded index (GRIN) devices. When the inclusions are dielectrics or metamaterials, GRIN devices can achieve effects such as negative refraction [5] and can function as elements such as flat surface lenses [6] or diffraction gratings [7]. When combined with transformation optics [8,9], the homogenization approximation can also be used to realize spatial permittivity and permeability maps with diverse capabilities including cloaking [10,11], wave bending [12], and electromagnetic wavefront rotation [13]. The great simplification of the homogenization approximation is critical for the understanding and design of such versatile photonic components.
Despite its widespread application, the range of applicability of the homogenization approximation has not been well quantified. A necessary condition for treatment of a material as homogenous is that the wavelength λ of light be sufficiently large compared to the characteristic size Λ of the structure or inclusion or period in a repeating pattern [14,15], with experiments often conducted with λ Λ ratios of 5-10 [16][17][18]. Yet simulation has shown that the homogenization regime may extend to ratios as small as 2.3 [19]. The condition becomes further unclear in waveguides where the wavelength of light can be on the order of the structure or inclusion thickness, requiring full 3D simulations to treat the geometry. With an understanding of the range of applicability of the homogenization approximation, analytical formulas can be used instead of simulations and can enable rapid device design and prototyping. This may enable successful designs in cases where fabrication technology limits the feature size such that high λ Λ cannot be achieved. For metamaterials (e.g. single resonance structures), it is also important to understand the limits of the homogenization approximation since the size of the structure must be specified to achieve a desired resonance frequency and it is not possible to make the unit cell any smaller. Recently, it was observed that breakdown of homogenization occurs when the ratio is 1.78 in a 45°-rotated square lattice of air-holes embedded in a dielectric host [20], but the influence of geometry, volume fill fraction (FF), and host material were not investigated. Here we systematically analyze the breakdown of homogenization in a dielectric slab, and compare direct experimental measurements of the effective index to the predictions made using analytical theory to rapidly design GRIN devices.
We experimentally study the breakdown of the homogenization regime using the terahertz (THz) polaritonics platform. In polaritonics, THz-frequency electromagnetic waves are generated, waveguided, and detected with full spatiotemporal resolution, typically on a 30-100 μm thick slab of lithium niobate (LN) or lithium tantalate (LT) [21][22][23][24]. In the present work, the approach is used to accurately map out dispersion curves of a diverse set of composite waveguides composed of lattices of air-holes in 50 μm thick LN or LT waveguides (see figure 1). The dispersion curve fully specifies the frequency-dependent refractive index of the composite slab by relating the wavelength (or wave vector) of the propagating THz light to the corresponding eigenfrequency. The facile ability to determine the dispersion curve is a distinctive capability of the polaritonics platform and is fundamental to this study, since the curve shows the response of the composite slab to a broad range of wavelengths. The range spans the transition from long wavelengths where the homogenization approximation is excellent to short wavelengths where it is quite poor and where the wavelength-dependent error can be determined. We chose to study homogenization in a lattice of air-holes i.e. a photonic crystal (PhC), because we could vary parameters such as lattice geometry, material, and the air-hole volume FF. This allowed us to examine different analytical theories for homogenization and determine their ranges of applicability. To further validate the applicability of the analytical theories, we designed two gradient index devices without any simulations: a Luneburg lens and a bi-directional cloak. We fabricated these devices on the polaritonics platforms, and recorded time-resolved movies of each.

Experiment
The lattice of air-holes, Luneburg lens, and bi-directional cloak were all fabricated in slab waveguides using a chemically-assisted femtosecond laser machining process [25]. A bare 50 μm thick slab of LN or LT was first coated with a protective layer of silicon dioxide (SiO 2 ) using plasma enhanced chemical vapor deposition. It was then machined by focusing the output of an ultrafast Ti:sapphire laser (800 nm center wavelength, 1 kHz repetition rate, 100 fs pulse duration) through a 0.25 NA objective, each laser shot ablating material from the sample surface. Repeated shots were used to drill holes all the way through the slabs. The samples were then etched to remove the SiO 2 , amorphous or damaged LN or LT around the hole edges, and re-deposited material.
Measurements of time-resolved THz wave propagation were performed using a self-compensating polarization gating imaging system [24]. First, a linearly polarized, 800 nm, cylindrically focused femtosecond pulse generated a transverse electric (TE) polarized THz wave through optical rectification, that was waveguided down the LN or LT slab (see figure 1(a) inset, bottom-right). The propagating THz electric field (E-field) induced a time-and position-dependent change in the LN or LT index of refraction through the electro-optic effect, imparting a position-dependent phase shift, ϕ x y ( , ), to a time-delayed, spatially expanded probe beam that was passed through the sample. The polarization gating imaging system converted the phase information into a spatially varying intensity change I x y ( , )that was detected by a camera. By changing the relative time delay between the THz-generating excitation pulse and the probe pulse, a series of images was recorded with the THz wave at different stages of propagation. The images can be played back in time-sequence to present a movie of THz wave propagation. For the air-hole lattice, the optical pump line source generated a broadband THz wave and spanned the 0.1-2 THz frequency range. For the gradient index devices, a titled intensity front [26] was used to generate a multicycle THz wave centered at 0.26 THz (Luneburg lens) or 0.30 THz (bi-directional cloak).
The generated THz waves are essentially plane waves as shown in the top-left inset of figure 1(a) where the THz wave is shown 12 ps after generation propagating along the +x direction in a bare LN slab of 50 μm thickness. Uniformity of the THz E-field along the y direction was used to collapse the intensity information from one frame in the movie into a 1D row vector representing I(x). Repeating this procedure with the image recorded at each time delay and displaying the collapsed row vectors I x t ( , )in time order along the vertical axis generated a 2D matrix or space-time plot as shown in figure 1(a). Figure 1(b) is the corresponding dispersion diagram for the bare LN waveguide, which was obtained by performing a 2D Fourier transform of the spacetime plot [27] and shows the first two TE waveguide modes. In figure 1(b), the horizontal axis is the propagation constant, β π λ = 2 , where λ is the propagation wavelength, and the vertical axis is the corresponding eigenfrequency. In this article, we will refer to λ as the propagation wavelength in the slab waveguide rather than the TEM wavelength λ ω = c n TEM because the former is the longitudinal wavelength component most relevant for interactions with the structures or inclusions. We overlay on figure 1(b) a waveguide dispersion curve that Top left inset shows a THz wave recorded by the camera, with uniformity along the ydimension exploited to collapse the image and create the space-time plot. Bottom right inset shows the experimental geometry, with cylindrically focused 800 nm pump pulse passing through the LN slab and generating counterpropagating THz waves, one of which may encounter the composite structure. (b) The wave amplitude at each frequency and wave vector is obtained from the 2D Fourier transform of (a) and shows the experimental TE waveguide dispersion diagram. The overlaid numerically calculated waveguide dispersion curves demonstrate agreement at all wavelengths as expected for this truly homogenous waveguide. (c) Space-time plot for the LN slab with a square lattice of air-holes, showing coherent reflections at 100 μm intervals corresponding to the lattice constant a. Inset shows an experimentally recorded image of a THz wave propagating inside the PhC slab. (d) PhC waveguide dispersion curve obtained by 2D Fourier transform of (c), showing good agreement at long wavelengths (small β) to a homogeneous slab waveguide curve dispersion curve (green) with a best-fit effective index value n eff . Inset shows magnified view of the lowest order mode where deviation between the homogenous (green) and experimental (Expʼt, red) curve demonstrates breakdown of the homogenization approximation. was calculated numerically [28] with the refractive index as an adjustable parameters optimized to reproduce the experimental bare LN slab waveguide curve accurately. In figure 1(c), we show the space-time plot for the 50 μm thick LN waveguide that was patterned by cutting a square lattice of air-holes through the slab (see figure 1(c) inset) with 26 μm diameter and 100 μm period. The spatially periodic intensity modulation observed in the xdirection is reduced at the locations of the air holes because the signal from the THz field is only detected in the electro-optic crystalline material, so the vertically integrated intensity is lower where some of the crystalline material is replaced by air. This DC spatial variation does not affect the dispersion curves which are only calculated for finite frequencies. THz wave propagation through and reflection from each row of air-holes gives rise to a frequency-dependent contribution. Figure 1(d) shows the corresponding dispersion diagram. The upper part of the lowest-order dispersion curve shows that for wavevector values above approximately 20 rad mm −2 , the patterned LN slab can no longer be treated as a bare slab with a constant index value. The dispersion curves are heavily modified because the patterned slab corresponds to a PhC waveguide, with several modes visible. Despite the non-uniformity of some PhCs along the y-direction (e.g. hexagonal lattice), the intensity can be collapsed along the y-dimension because of the use of a line-source excitation, which most significantly projects onto modes that are spatially uniform along y. The validity of this approach is proven by strong agreement of the experimental dispersion curves with simulations [29].
According to the homogenization approximation, for all wavelengths considerably larger than the size of the air-holes, the PhC waveguide is equivalent to a slab waveguide with a single effective index (neglecting material dispersion in the applicable long-wavelength range). To demonstrate this, we overlay on figure 1(d) a numerical waveguide dispersion curve (green solid) calculated for 50 μm thickness with the refractive index determined by fitting to the experimental curve at long wavelengths. We focus specifically on the lowest order TE mode since this is the one that can satisfy the long wavelength regime in our experiments. Agreement between experimental data and the fitted numerical curve at long wavelengths suggests that the PhC waveguide can be treated as a homogenous slab with an effective index of = n 4.15 eff in this regime. However, at shorter wavelengths the homogenization approximation breaks down (see figure 1(d) inset); the two dispersion curves begin to diverge as the THz wavelength approaches the size of the air-holes and coherent in-phase reflections begin to heavily modify the dispersion from that of a bulk material to that characteristic of a PhC with bandgaps. Since the experimental dispersion curve includes the range of wavelengths in which the deviation begins to become significant, we can determine the wavelength-dependent percent difference in frequencies between the fitted and experimental curves (see figure 2).

Results
We determine the influence of lattice geometry and orientation (relative to the wave propagation direction), host material, and areal FF of the holes by comparing the cutoff wavelengths in a number of diverse composite waveguides. The influence of various parameters on λ cutoff can be used as an expedient measure of their influence on the wavelength range over which the PhC structure can be treated as homogeneous. These parameters also the most readily modifiable in creating GRIN or transformation optics devices, and understanding their effects will guide the design of future devices. In PhC structures, a fundamental limit for homogenization is given by the Bragg effect in which the reflections from alternating layers add in phase to provide a strong net reflection. This typically belongs to a wave at the edge of the Brillouin zone with wave vector k Br or Bragg wavelength Λ π = k 2 Br Br , and will differ depending on geometry and orientation. To account for this, we normalized the propagation wavelength (λ π β = 2 ) to the Bragg wavelength. Further, we define λ cutoff as the wavelength at which a 5% difference occurs between the homogenized and experimental frequencies, which is often accurate enough for device design (see applications) and so can be taken as a lower limit for the applicability of the homogenization treatment. The results are summarized in table 1, while figure 2 shows the wavelength-dependent percent difference between the experimental and homogenized dispersion curves. As expected, the percent change is considerable when λ Λ Br is near unity, but decreases toward zero as λ Λ Br increases and the homogenization approximation becomes valid.

Lattice geometry and orientation
We first studied the influence of lattice geometry and orientation while holding the fill factor and lattice periodicity constant at = ± FF 0.21 0.1 and μ = a 100 m in a 50 μm thick LT slab (see figure 2(a)). We recorded the dispersion curves and compared the cutoff wavelengths in square lattices oriented at 0 o (Γ to X) and 45 o (Γ to M) and in hexagonal lattices oriented at 0 o (Γ to M) and 30 o (Γ to K). The optical and THz wave polarizations were along the optic axis of the LT crystal in all cases. In all four lattices, the cutoff wavelength was found to be λ Λ = ± 1.3 0.1 cutoff Br (see figure 2(a), table 1(a)), indicating consistency between different lattice geometries and orientations. The consistency shows that Λ Br better represents characteristic length scale than the lattice constant a, i.e. despite all four lattices having different Λ Br values (see table 1(a)), homogenization breaks down in all cases at the same λ Λ Br ratio. The role of Λ Br demonstrates the influence of the in-phase reflections (i.e. Bragg reflections) and indicates that the approach toward the PhC regime serves as the fundamental limiting factor to the homogenization approximation. For some applications, it would be sensible to use the geometry with the smallest values of Λ Br because it allows for the use of a smaller λ cutoff and thereby relaxes fabrication. Alternatively, using an aperiodic lattice of inclusions can remove coherent in-phase reflections and also extend λ cutoff to smaller wavelengths. when all other parameters are kept constant (see figure 2(b), table 1(a)). The result is unsurprising due to the difference in Table 1. Influence of geometry and material on the cutoff wavelength (±0.1), defined as the wavelength at which a 5% difference between the homogenized and experimental frequencies occurs. Bulk extraordinary indices are (±0.1 dependent on sample) = n 6.44 LT , = n 5.0 LN . Lattice constant μ = a 100 m in all structures listed. There is consistency between geometries at the same FF, but differences between materials due to index differences. bulk index between LT (n = 6.4) and LN (n = 5.0). Material clearly plays a role in the form of index contrast [30], with higher index contrasts resulting in stronger Fresnel reflections and stronger modification of the dispersion relation of the composite material, thereby increasing λ cutoff . When designing homogenization approximation based devices, it may be useful to try to minimize the index contrast to allow operation at smaller wavelengths, while still allowing sufficient contrast to meet the required index tuning range. In section 4, we build our devices in LT instead of LN since it offers a wider index tuning range while only having a slightly larger λ Λ cutoff Br . In general, other performance characteristics of a material that should be considered are the intrinsic material losses and dispersion.

Volume FF
Although λ Λ cutoff Br is minimally influenced by lattice geometry and orientation, it depends strongly on the volume fill factor (see figure 2(c)). Here, we analyze FF of 0.12, 0.2 and 0.40 by varying the radii of the air-holes in a hexagonal lattice with results summarized in table 1(b). The host material is LN and Λ Br is preserved at 173 m Br in all three structures. The plot shows that larger FF increases λ cutoff . It might be suggested that larger air-holes produce a smaller effective index, thereby naturally increasing λ cutoff . However, our results show that λ cutoff increases by more than the inverse of the effective index. The results indicate that lattice periodicity cannot be the only metric used to evaluate the homogenization approximation, as the size of the inclusions must be taken into account. Larger inclusions result in more scattering and produce a larger band gap. This requires that the light depart from homogenous behavior and begin PhC behavior at longer wavelengths, giving rise to a larger λ cutoff value.
Our experimental methods make it difficult to discern the behavior of λ cutoff at FF far past 0.6, since the signal comes entirely from the host material and is reduced when the air-holes take up most of the volume. Further complicating matters is that as FF becomes large, the waveguiding condition is lost since the material becomes largely air. In any specific system, the actual FF at which λ cutoff reaches a maximum will depend on many factors including the nature of the inclusions and the lattice geometry. Our simulations for square and hexagonal lattices indicate that the largest λ cutoff occurs when the perturbation to a bulk host medium is largest, which is in the vicinity of ∼ − FF 0.5 0.8. At even larger FF the width of the band gap, and therefore λ cutoff , decreases as the perturbation shrinks and host and inclusions reverse roles.

Length scales
In a bulk substrate (neglecting material dispersion), λ Λ cutoff Br remains constant as the entire system is scaled up, i.e. λ cutoff and Λ Br both increase proportionally. This is due to the scale invariance of Maxwellʼs equations. But in a waveguide with a fixed slab thickness, λ Λ cutoff Br will decrease as the system is scaled up because λ cutoff does not increase proportional to Λ Br . At larger wavelengths, the light is more heavily extended into the surrounding air and less affected by structures or inclusions, thereby reducing λ cutoff compared to the bulk substrate value. However, this was not a significant factor for the length scales investigated, where increasing the size of the system ( μ = r 21 m, μ = a 60 m to μ = r 35 m, μ = a 100 m) did not show a change in λ Λ cutoff Br . Nonetheless, it is desirable to use structures on the order of the slab thickness when designing waveguide devices because this allows for both localization in the slab and strong interaction with the structures.

Parameter trends
Several FFs, lattice geometries and orientations, length scales, and two high index contrast materials were tested in order to determine trends in λ Λ cutoff Br . Despite the dependence of λ Λ cutoff Br on these parameters, we found that satisfying λ Λ ⩾ 2 Br allowed the homogenization approximation to hold regardless of the parameters used.

Analytical formulae
An analytical formula that can closely predict the effective index of a composite material in the long wavelength limit is a useful tool for rapid development of GRIN devices, and can prove useful as a starting point for simulations. Here we compare the effective index calculated by analytical formulae to that obtained by fitting to the experimental data. Comparison of these values can be used as a metric to determine the predictive power of the formulas. Specifically, we examine the well-established Maxwell-Garnett mixing formula (MGF) [31] that is often used to estimate the effective index for cylindrical inclusions in a host material when the E-field is perpendicular to the height of the cylinder [14,19]  In deriving both JSF and MGF, the E-field within the inclusions is related to the local external field in order to calculate n eff . Figure 3 compares the fitted effective index (n fit ) with the analytical effective index (n MGF and n JSF ) for a number of lattices, materials, and FFs. Error bars are a result of an upper and lower bound on the air-hole radius; the black annuli around the air-holes (see figure 1(c) inset) makes it difficult to report an absolute air-hole diameter. Points closer to the y = x line demonstrate an analytical effective index that more closely matches the fitted effective index, and shows that MGF has slightly stronger predictive power. This is particularly interesting at larger FF (smaller n eff ), where MGF is expected to suffer by not accurately factoring in the effect of neighboring inclusions, with exceptions in the case of random configurations [33]. The validity of MGF at high FF in the periodic arrays of inclusions studied here ( ⩽ FF 0.5-0.7 in all data) suggests that the dipole field from neighboring inclusions may not be significant in our experiments, and JSF is overestimating their contributions to n eff . Based on previous studies, it is expected that JSF will show stronger predictive power beyond a threshold FF, when the inclusions become connected [32]. But in our PhC waveguides, MGF proves it can be very useful in predicting the effective index.

Applications
To test the applicability of the homogenization condition and analytical formulae, we designed and fabricated a Luneburg lens [34] and a cloak in the polaritonics platform. The time-resolved imaging of the THz fields inside the devices offers direct insight into device operation. In addition to testing homogenization theories, transformation optics and GRIN devices expand the capabilities of the polaritonics toolkit and offer great potential for THz telecommunications or all-optical processing systems. Designing devices in a waveguide is more complicated than in bulk due to the finite lateral dimension and electromagnetic field profile. In a bare slab waveguide, higher frequency waves are more highly concentrated in the slab, and this can be translated to a frequency-dependent waveguide index ω n ( ) wg [28]. In designing devices for a specific frequency ω o , we used ω n ( ) wg o as the value for the bulk index. We fabricated the devices in a 50 μm thick LT slab waveguide, and tuned the index by varying the radius of air-holes and calculating the effective index using MGF. The analytical formula proves extremely useful for determining the air-hole size required to achieve the desired index, which would otherwise require an extensive database of fitted curves or time-consuming simulations.
The first device demonstrated here is a rotationally symmetric lens (see figure 4(a), media 1), known as a Luneburg lens. Unlike conventional lenses, the Luneburg lens is capable of focusing a plane wave incident from any direction to a point on the other side of its surface. The lens was fabricated to have an index distribution where r is the radius from the center, = R a 10 , and n o = 4. Preserving a unit cell at μ = a 100 m and a center wavelength of λ μ = 300 m, this translates to an increasing index of 2.6 at the center where there are no air-holes, to an index of 4.3 at the edge with a hole radius of 37 μm. Figure 4 and media 1 show the time evolution of a THz wave as it passes through the lens. As it propagates, the curvature of the THz wave changes (compare figures 4(a) and (b)) before finally focusing within the lens ( figure 4(c)). The operation of the device is encouraging, with a focal spot size of μ ∼200 m that is slightly larger than the diffraction limit of λ μ = 2 150 m in an ideal Luneburg lens [14]. The incomplete focusing can be attributed to the curvature and finite bandwidth of the incoming THz field that result from the multicycle THz wave generation process. The amount of light coupled in to the lens is reduced due to the impedance mismatch at the edge, but can be overcome by gradually varying the effective index from n host to n edge in order to create an anti-reflection coating. The highest FF in this device is 0.44 where it is estimated that λ μ ≃ 380 m cutoff (λ Λ ≃ 1.9 cutoff Br , square lattice Λ = a 2 Br , μ = a 100 m). The successful qualitative operation of this device at λ λ < cutoff is a potential indication that it is possible to operate beyond the cutoff wavelength by removing the periodicity that creates the Bragg diffraction, or that a 5% difference between homogenized and experimental frequencies used in our analysis is an overly strict condition for the cutoff wavelength.  A bi-directional cloak, a device that routes the path of incoming light around an object, was also studied in order to demonstrate further control over the THz waveform (see figure 5, media 2, and media 3). A THz wave of λ μ = 260 m was used with a square lattice constant of μ = a 73 m. The device has an inverted index profile of a lens so that it excludes light from the focal point [19]. Additional air trenches above and below the cloak (not shown in figure 5) allow the redirected light to be collected at the output face. Successful operation of the device can be seen by comparing the interaction of a THz wave with a structure in the absence and presence of the cloak. Figures 5(a) and (b) show before and after images of a THz wave as it encounters a large triangular air-hole, after which the THz wave forms a distinct diffraction pattern (media 2). Figures 5(c) and (d) show the results when the triangular air-hole is embedded in the cloak (media 3). Here the THz wave is routed around the object, and the resulting output retains the curvature of the input wave. The cloak operates independent of the material placed inside because the light is diverted around the center of the cloak.

Conclusion
In this article, we first analyzed composite waveguides to quantitatively determine the wavelength range over which the homogenization approximation is valid. The onset of homogenization, when appropriately normalized to the Bragg wavelength, was nearly independent of lattice geometry and propagation direction. However, the homogenization did depend strongly on FF and was weakest when the sample was roughly 50-70% 'inclusion'. Fundamentally, homogenization broke down as the composite material entered the PhC regime of operation where coherent in-phase reflections lead to a photonic band gap. Also, note that the THz index contrast between air and LN or LT is much larger than the contrast at optical wavelengths where typical waveguides have refractive indices in the range of 1-3.5. Combining the (i) larger index contrasts studied here, (ii) the cutoff condition of 5%, and (iii) having probed FF up to range of 0.5-0.7 near the perturbation to a bulk slab is largest, suggest that λ Λ ⩾ 2 Br should serve as an adequate condition for most conventional dielectric effective index devices.
We subsequently used the effective index extracted from our experimental curves to compare the predictive strength of analytical formulae, with the MGF showing excellent predictive quality even at higher FF. Equipped with an analytical formula, we designed and performed time-resolved studies of GRIN devices. These devices seamlessly interface with existing polaritonics capabilities such as focusing [21], field enhancement [35], and pulse shaping of the driving femtosecond pulse [23], and could serve practical uses in optical processing on the polaritonics platform. Specifically, a Luneburg lens could prove useful for coupling light into an on-chip THz detector, collimating light from an on-chip THz source, or for applications requiring high field strengths in tightly confined dimensions. The cloak may benefit systems requiring planar integration of electrical and optical components, since it is possible to cloak objects such that propagating light is minimally affected by their presence. Other platforms with similar fabrication constraints may also benefit from the results shown here, since the conditions are dimensionless and the results should be valid at all wavelengths. Future work may include analyzing the effect of more complicated inclusion shapes and symmetries.