## Abstract

Femtosecond optical pulses were used to generate THz-frequency phonon polariton waves in a 50 micrometer lithium niobate slab, which acts as a subwavelength, anisotropic planar waveguide. The spatial and temporal electric field profiles of the THz waves were recorded for different propagation directions using a polarization gating imaging system, and experimental dispersion curves were determined via a two-dimensional Fourier transform. Dispersion relations for an anisotropic slab waveguide were derived via analytical analysis and found to be in excellent agreement with all observed experimental modes. From the dispersion relations, we analyze the propagation-direction-dependent behavior, effective refractive index values, and generation efficiencies for THz-frequency modes in the subwavelength, anisotropic slab waveguide.

©2010 Optical Society of America

## 1. Introduction

Terahertz-frequency phonon polariton generation, control and detection have received extensive attention in recent years due to their outstanding capabilities in terahertz (THz) spectroscopy, imaging and advanced signal processing [1–5]. Phonon polariton waves result from the coupling of lattice vibrational waves and electromagnetic waves, and can be generated in ferroelectric crystals such as LiNbO_{3} (LN) via impulsive stimulated Raman scattering (ISRS) using femtosecond optical pulses [6,7]. The electromagnetic component of the phonon polariton wave can be coupled into free space and is a source for intense THz pulses [8–12]. THz waves generated in the sample do not propagate collinearly with the pump beam due to the large index-mismatch between optical and THz frequencies. Instead they generate a Cherenkov radiation pattern and propagate primarily in the lateral direction [13, 14]. This lateral propagation facilitates coherent control of the THz wave, which can easily be made to interact with subsequent optical pulses, other THz waves, or patterned structures all in the same small crystal of LN. As a result, a LN slab can serve as a platform for THz processing because generation, propagation, detection, and control can be fully integrated in one sample [5, 15]. Furthermore, when the sample thickness becomes comparable to or less than the THz wavelength, the strong evanescent field of the THz wave can interact with material deposited on the crystal surface. This opens the door for spectroscopic analysis and interfacing the LN slab with other optical or photoelectric devices.

Because the THz wave propagates almost perpendicular to the optical pump beam, it is possible to obtain time-resolved images of the electric field in the LN slab. As the THz wave propagates through the crystal, its electric field changes the refractive index through the electro-optic effect. The time-delayed probe pulses, which can be expanded to illuminate the whole crystal, experience a spatially dependant phase shift proportional to the refractive index change. Four methods have been introduced to convert this phase pattern to an amplitude image: Talbot imaging [2], Sagnac interferometry [16], polarization gating [16,18], and phase contrast imaging [17]. In a recent comparison [18], an improved geometry for polarization gating was found to offer the best sensitivity and most reliable field quantification, while phase contrast imaging was best in situations requiring high spatial resolution. In this paper, we used the polarization gating system similar to that shown in [18] to record a sequence of images. The full spatio-temporal evolution was extracted from the image sequence and double Fourier transformed to obtain the wave vector vs. frequency dispersion curves [e.g 17.]. The data collection and analysis were performed as a function of wave propagation direction to study the complex mode structure present in an anisotropic slab waveguide, which was found to be in excellent agreement with theory. From the dispersion relations we extract the mode and propagation-angle dependent effective refractive index (ERI) and discuss pumping efficiencies for THz phonon polariton waves in a LN waveguide.

## 2. Experimental section

The experiments were performed with a Ti:sapphire regenerative amplifier whose pulse duration was 120 fs, central wavelength was 800 nm, and repetition rate was 1 KHz. The laser pulses were divided into a pump beam (370 μJ per pulse) and probe beam (35 μJ per pulse). The vertically polarized pump beam was routed through a mechanical delay stage and then focused to a line on the sample by a 200 mm focal length cylindrical lens (about 1 TW/cm^{2}). The probe was frequency-doubled to 400 nm in a BBO crystal and expanded to be larger than the sample. The probe beam is nearly collinear with the pump by using a dichroic mirror, so the second harmonic wave of the pump on the sample, whose wavelength is the same as probe, can be blocked with a razor blade on the focal plane of the imaging lens. Figure 1(a)
shows a sketch of the experimental setup and the coordinate system. A quarter-wave plate (QW1) and a retroreflective mirror were used in a 4-f system. The mirror and lenses imaged the sample precisely back onto itself without magnification or inversion. The axis of QW1, which was the same as the first Glan-Taylor polarizer (GTP1), was at + 45° so it exchanged the ordinary and the extraordinary polarization components of the probe. In this way the spatially varying phase shift between the vertical and horizontal polarization components accumulated from the probe’s first pass through the sample was compensated after the second pass. The phase shift after the first pass resulted from the intrinsic birefringence of the LN slab, and self-compensation was necessary to correct for spatial inhomogeneities in the phase shift due to thickness variation, strain, or other imperfections in the slab. The phase shift electro-optically induced by the THz wave, however, was not compensated because the THz wave was launched only after the probe pulse had passed through the sample the first time. The THz-induced phase information was converted to amplitude information prior to detection with the camera by QW2 (oriented vertically) and GTP2 (oriented at −45°). In this geometry a positive field results in a positive amplitude change and vice versa [18].

The pump geometry is shown in Fig. 1(b). Red lines represent the 800 nm pump beam and green the broadband THz waves generated when the pump is focused into the 50 μm thick LiNbO_{3} crystal slab. Because the center wavelength of the THz phonon polariton wave is about 100 μm, the slab acts as a sub-wavelength waveguide. As Fig. 1(c) shows, the THz wave propagation direction was changed by rotating the cylindrical lens. Because of the strong anisotropy of LN at THz frequencies, the nature and behavior of the waveguide modes change drastically as the propagation direction rotates relative to the optic axis.

## 3. Results

By changing the delay between the pump and probe pulses, a series of images can be obtained. The image sequence can be compiled to form a movie showing THz propagation [2–5, 17, 18]. Because the line focus launches plane-wave THz transients that propagate laterally away from the focal region (i.e. the origin), in each image recorded at a different probe time delay the signal was uniform along the direction of the line focus. Therefore at each propagation distance from the origin, we averaged the signal over this direction, collapsing each 2D matrix of values and the corresponding image to 1D. We then displayed each 1D image corresponding to a selected time delay as a horizontal line, and we displayed the 1D images one above the other in time order, showing the full temporal and spatial information in a single graphic (Fig. 2(a)
). From Fig. 2(a) we can see dispersion, reflection, and different waveguide modes clearly and the THz electric field E(t) also can be acquired using the same method as [18]. A 2-dimensional Fourier transformation of Fig. 2(a) yields the THz dispersion curves (Fig. 2(b)). Along the vertical axis, time is transformed to frequency, and along the horizontal axis, space is transformed to the wave vector, *k _{x}*, which is often called the waveguide propagation constant,

*β*. In Fig. 2,

*θ*= 0, so the crystal’s

*c*-axis is parallel to the 800 nm pump polarization. In this geometry, which has been most used in previous work [2–5,15–18], only

*z*-polarized THz is generated, and true transverse electric (TE) modes are launched in the slab. Overlaid on the experimental data are the dispersion curves for air (white line), bulk LN (magenta line), and the calculated TE mode dispersion curves (see e. g [19].) up to a frequency of 2 THz for an isotropic slab waveguide with

*n*=

*n*(dotted blue lines). The curves show four TE waveguide modes, which propagate at different group velocities, ${v}_{g}=d\omega /dk$, and phase velocities, ${v}_{p}=\omega /k$. Cutoff frequencies can be seen for the all but the first mode as expected. Although the isotropic waveguide analysis is predictive in this simple geometry, a more complete analysis is required when

_{e}*θ*≠ 0, as will be shown below.

In an anisotropic waveguide, constraints relating to propagation in bulk anisotropic material and constraints relating to propagation in a waveguide both come into play. In bulk anisotropic material waves are divided into two normal modes, ordinary waves and extraordinary waves, which propagate through the material at different velocities [20]. In an isotropic waveguide there are also two uncoupled eigenmodes, the transverse electric (TE) and transverse magnetic (TM) modes, which propagate through the waveguide at different velocities [19]. When *θ* = 0 or 90°, these modes map directly onto one another. For *θ* = 0° the TE mode is an extraordinary wave and the TM mode is an ordinary wave while for *θ* = 90° the opposite pairing holds. When *θ* ≠ 0, however, the high-symmetry configuration is broken and all the modes couple together. The new eigenmodes of the system are neither purely TE nor TM and also not purely ordinary or extraordinary. The coupling effects the mode profiles, dispersion curves, and effective refractive indices in a fundamental and significant way, as will be demonstrated experimentally (presented immediately below) and theoretically (the full analysis can be found in the appendix) in the remainder of this paper.

With the experimental system mentioned above, we measured the dispersion curves for different propagation directions by rotating the cylindrical lens and CCD camera together, which kept the THz wavefront aligned vertically in the images. In this manner the THz wave propagation direction was varied from 0 to 90 degrees relative to the *c*-axis as shown in Fig. 1(c). The polarization of the 800 nm pump light was not rotated and thus was parallel to the *c*-axis in all measurements. Because of the strong *r*
_{33} electro-optic coefficient in LN, this ensured efficient pumping of THz waves with a large component polarized along the optic axis [21–23]. Using the same data collection and analysis procedure as was used to generate Fig. 2(b), the dispersion curves were measured for different angles *θ*, some examples of which are shown in Fig. 3
. Overlaid on the experimental data are the theoretical solutions for a bound mode propagating in an anisotropic, dielectric slab waveguide with a thickness of 50 μm, an extraordinary index of 5.11, and an ordinary index of 6.8. The full derivation is presented in appendix A. For all modes and all angles, the data agree very well with theoretical predictions.

In Fig. 3(a) where *θ* = 20°, one set of modes is very TE-like, and one is strongly TM-like. Because the TE-like modes have their primary polarization component along the *c*-axis, they were pumped much more strongly than the TM-like modes, which were too weak to be observed clearly. Blue dotted lines are calculated TE-like modes and green dashed lines are TM-like modes. An interesting effect resulting from propagation in the anisotropic waveguide when *θ* ≠ 0 or 90° is visible in the region containing the white lines, a magnified view of which is shown in the lower right corner. Although the TM-like modes are not visible, we still see an avoided crossing when two modes with the same symmetry (symmetric or antisymmetric) cross. The avoided crossing, visible in experiment and predicted by theory, results from coupling between TE- and TM-like modes in the anisotropic waveguide. Here the two lowest symmetric modes, the lowest TE-like mode and the second TM-like mode, avoid each other.

Figure 3(b) shows dispersion curves for the case of *θ* = 50°. The three TE-like modes predominate, although their strength is reduced, and TM modes still cannot be observed. No avoided crossings occur between modes of the same symmetry within the bandwidth of the experiment. As *θ* increases, the velocity of the extraordinary wave approaches that of the ordinary wave [30], which means that the slopes of TE-like and TM-like modes tend to be more similar at higher frequencies. Figure 3(c) shows the results for *θ* = 70°, where both TM-like and TE-like modes can be seen clearly. The strength of the TE-like modes continues to decrease with increased *θ* and TM-like modes are finally pumped strongly enough to detect. Although some of the modes are weak, the first seven modes can be observed in the experiment, all of which agree with theoretical predictions. Continuing the trend, at high frequencies the slopes of the TM-like and TE-like modes become even more similar. Finally, Fig. 3(d) shows results for *θ* = 90°, where only TM modes can be observed.

Using the derivation in appendix A, we can calculate the E-field profile of THz waves as shown in Fig. 4
. Figures 4(a) and (b) show field profiles for TE and TM modes respectively at *θ =* 0°. The coordinate system in Fig. 4 is the same as in the appendix (see Fig. 8
), where the axes are defined by the propagation direction of the wave and not by the lab frame as in Fig. 1. Blue, green and red lines represent electric field along *x-*axis, *y*-axis and *z*-axis respectively. The electric field along the *y*-axis, whose polarization is perpendicular to the surface of the slab, changes drastically at the slab surface ( ± 25 μm). As mentioned above, pure TE and TM modes only exist at 0 and 90 degrees. At any other angle the eigenmodes are superpositions of TE and TM modes and contain all three polarization components, as shown in Fig. 4 (c) and (d) where *θ* = 50°.

One can extract the group and phase effective refractive index (ERI) from dispersion curves like those shown in Fig. 3. The phase ERI can be retrieved directly from the dispersion curve, ${n}_{p}=c/{v}_{p}=ck/(2\pi f)$, and the group ERI can be retrieved from the slope of the dispersion curve, ${n}_{g}=c/{v}_{g}=c\Delta k/(2\pi \Delta f)$. Here *v _{p}* and

*v*are the phase and group velocities, and the wave vector,

_{g}*k*, and frequency,

*f*, correspond to the axes in Fig. 3. From the derivation in the appendix, we can calculate both phase and group ERI for all propagation directions of TE-like and TM-like modes in the LN slab waveguide. Based on the agreement between experimental data and theoretical predictions shown in Fig. 3, Fig. 5 gives the theoretically calculated phase and group ERI for different angles, modes and frequencies. The ERI values are important for phase-matching in THz generation and for many nonlinear as well as linear optical processes.

In Fig. 5 dotted blue lines are the calculated ERI for TE-like modes and dashed green lines are the ERI for TM-like modes. From Fig. 5(a) and (b), we can see that the phase ERI for both TE-like and TM-like modes transitions from 1 (the index of air) to the bulk effective index. For the TM-like waves the bulk index is always the ordinary index of refraction, *n _{o} ≈* 6.8, while for the TE-like waves the bulk index is that for the extraordinary wave in the anisotropic material, and changes from ~5.1 at 0° to ~6.8 at 90°. At all angles, the low-frequency TM-like modes have most of their energy in the evanescent field in the air and have ERI values near unity. The ERI then transitions rapidly to bulk-like values at higher frequencies. Much like the phase ERI, the group ERI transitions from 1 to the bulk effective index (see Fig. 5 (c) and (d)). In contrast to the phase ERI, however, the group ERI rises well above the bulk values before approaching them asymptotically at high frequencies. In contrast to the phase index, where higher modes always have lower ERIs, the group ERI is usually higher for higher modes. Another difference is that the peak group index changes drastically with

*θ*for both TE-like and TM-like modes, while the peak phase ERI for the TM-like modes is just the bulk value and insensitive to angle.

A useful way to display the ERI is with an index ellipse, which highlights the angle-dependant behavior. Figure 6
follows the phase ERI for the first three TE-like modes at wave vector magnitude *β* = 50 rad/mm as a function of angle, tracing out the phase ERI ellipse. We measured data in the first quadrant, and because of the symmetry these results can also be used for 90° to 360°. Values over 70° were not recorded because the TE-like modes were too weak to be observed. Figure 6 shows the measured values for the first three TE-like modes as open symbols and the calculated values predicted by the derivation in the appendix as solid lines. The experimental data can be fit to an ellipse, where the long and short axes are 5.44 and 4.18 for the first mode, 3.36 and 2.59 for the second mode and 2.01 and 1.74 for the third mode. The value of the long axis represents the ERI for an ordinary wave (the TE mode is purely ordinary at 90°) and the short axis represents the ERI for an extraordinary wave (the TE mode is purely extraordinary at 0°).

Because the 800 nm light was always polarized along the *c*-axis, we only pumped through the *r*
_{33} electro-optic coefficient, which generated THz polarized along the optic axis [5]. The THz generated by the pump can be represented by a linear combination of waveguide modes, and the magnitude of the contribution from a given mode is related to the projection of its polarization along the *c*-axis. Thus, for a given mode and frequency, one can make a rough estimate of the relative pumping efficiency $\eta (\theta )$ by looking at the fraction of the mode energy corresponding to a field inside the crystal oriented along the optic axis:

Figure 7
shows *η* as a function of *θ*. When *θ* = 0°, the TE mode is polarized purely along the optic axis and is pumped most efficiently. As *θ* increases, the component of the TE wave along the optic axis slowly decreases. In contrast, the component of the TM wave along the optic axis increases, especially after 60°, and at 90°only the TM mode is pumped. At 70°, both modes are pumped with similar efficiencies. The qualitative trends in *η* explain the mode amplitudes observed in Fig. 3. For *θ* less than about 60°, TM-like modes are too weak to be observed, while for *θ* more than about 85°, the TE-like modes are not visible. As predicted, both kinds of modes are visible at 70° as shown in Fig. 3(c). Using different pump polarizations and reflective elements integrated into the waveguides, it will be possible to generate modes not observed in this study.

## 4. Conclusions

We have measured the propagation properties of THz waves in a 50 μm LiNbO_{3} anisotropic slab waveguide using a self-compensating polarization gating imaging system. This system can detect the THz electric fields both temporally and spatially over a wide wavelength range. Using the system, we studied the propagation-direction-dependent behavior of waveguide modes and determined the dispersion curves and effective refractive index for THz waves. A general solution for waveguide modes in a uniaxial slab waveguide was derived and found to agree with the experimental data.

Dispersion is integral to many processes in THz science and generally in linear and nonlinear optics, including broadening of ultrashort pulses, walk-off between pump and probe pulses, phase-matching of parametric processes, and generation of optical solitons. Because dispersion in a waveguide is determined by both the intrinsic material dispersion and geometric dispersion, it is essential to understand waveguiding effects. The results presented here will facilitate the design of functional devices with new capabilities in the LiNbO_{3} platform for integrated THz experiments and processing.

**Appendix A:** The general solution to a uniaxial slab waveguide with isotropic cladding

Anisotropic slab waveguides were extensively studied in the 1970’s [24–30]. In many cases, attention was focused on anisotropic films deposited on a substrate that was itself anisotropic because mode converters, polarization mode filters, and other devices of that time had such a geometry [26]. The case in this paper is somewhat simpler because the geometry is symmetric (see Fig. 8). An additional simplification is that the anisotropic core (the slab) is embedded within an isotropic cladding (air in our experiment). In the derivation of the waveguide dispersion curves and mode profiles presented below, we assume the experimentally relevant conditions that the crystal is uniaxial (like LiNbO_{3}) and its optic axis is parallel to the slab surface. The slab is assumed to extend infinitely along *x* and *z* and both core and cladding have no magnetic response. The wave is assumed to propagate along the *x*-direction and extend infinitely along the *z*-direction. In the experiment the cylindrical lens generating the wave was rotated instead of the sample, so the derivation here is performed in the coordinate frame of the lens. Finally, to simplify the analysis we assume that the waves are harmonic in space and along the propagation direction: $\stackrel{\rightharpoonup}{E}(x,y,z,t)=\stackrel{\rightharpoonup}{E}(y)\mathrm{exp}[i(\beta x-\omega t)]$, where $\beta ={k}_{x}$ is the propagation constant.

The derivation presented below will loosely follow the analysis of Marcuse and Kaminow [30] where the more complicated symmetric geometry of an anisotropic slab with anisotropic cladding is studied. For the sake of brevity our analysis will skip some intermediate steps, many of which can be found in [30]. The first step in the derivation is to determine the characteristics of waves in bulk material, i.e. the dispersion curves and polarizations, in both core and cladding. Linear combinations of these bulk waves, constrained by system symmetry, are used to build the waveguide modes and lay out the general functional form of the solution. The boundary conditions at the waveguide surface generate a homogeneous system of equations which can be used to solve for the coefficients in the linear combination. Solutions exist for this system of equations, i.e. the determinant of the corresponding matrix is zero, only for certain pairs of frequency and propagation constant. These allowed solutions correspond to the waveguide dispersion curves.

To simplify notation we define several important variables. The propagation constant, $\beta \equiv {k}_{x}$, was defined above, and the wave vector orthogonal to the slab surface is defined both outside the crystal, $i\alpha \equiv {k}_{y}^{\text{out}}$, and for the ordinary and extraordinary waves inside the crystal, ${\kappa}_{o},\text{\hspace{0.17em}}{\kappa}_{e}\equiv {k}_{yo}^{\text{in}},\text{\hspace{0.17em}}{k}_{ye}^{\text{in}}$. *α* is defined as imaginary because bound modes will have evanescent, decaying fields in the cladding. There are three relevant bulk dispersion curves which define the relationships between wave vector, frequency, and index, one for the cladding and one each for the ordinary and extraordinary waves in the uniaxial core. They are:

where $k=\omega /c$ is the wave vector in free space and ${n}_{c},\text{\hspace{0.17em}}{n}_{o},\text{\hspace{0.17em} and}{n}_{e}$ are the cladding index, ordinary index in the slab, and extraordinary index in the slab respectively. These relations are used to eliminate $\alpha ,\text{\hspace{0.17em}}{\kappa}_{o},\text{\hspace{0.17em} and}{\kappa}_{e}$from the equations which follow, so everything is expressed in terms of *β* and *k*.

For a specific pair of *β* and *k*, there are four possible plane waves in each region, two signs for ${k}_{y}$ and two polarizations. In the anisotropic medium, the polarizations correspond to the ordinary (later represented by $\stackrel{\rightharpoonup}{o}$) and extraordinary (later represented by $\stackrel{\rightharpoonup}{e}$) waves. In the cladding, any pair of orthogonal polarizations can be chosen, so for convenience we choose the TE polarization (represented below by $\stackrel{\rightharpoonup}{v}$ for vertical polarization) and the TM polarization (represented later by $\stackrel{\rightharpoonup}{h}$for horizontal polarization). We can write out the most general form of the waveguide mode solution as:

where *A*
_{i}, *B*
_{i}, and *C*
_{i} are scalar constants and the +/− superscripts correspond to the sign of ${k}_{y}$.

The polarizations in the expression above can be determined from the appropriate vector constraints. In the cladding:

where ${h}_{x}$ and ${h}_{y}$ are the magnitudes of the components of the normalized polarization vector. In contrast to the isotropic cladding, where any orthogonal polarizations could be chosen, in the slab the polarizations are uniquely determined as the ordinary and extraordinary wave polarizations in bulk material. The ordinary wave will be orthogonal to the plane containing the crystal axis and the wave vector: $\stackrel{\rightharpoonup}{o}\propto \stackrel{\rightharpoonup}{k}\times \stackrel{\rightharpoonup}{c}$. The extraordinary wave will be located in the plane of $\stackrel{\rightharpoonup}{k}$ and $\stackrel{\rightharpoonup}{c}$. The displacement field will be given by$\stackrel{\rightharpoonup}{D}\propto \stackrel{\rightharpoonup}{k}\times (\stackrel{\rightharpoonup}{k}\times \stackrel{\rightharpoonup}{c})$, and the electric field is given through the constitutive relation: $\stackrel{\rightharpoonup}{E}=\overline{\overline{R}}(-\theta ){\overline{\overline{\epsilon}}}^{-1}\overline{\overline{R}}(\theta )\stackrel{\rightharpoonup}{D}$ where $\overline{\overline{R}}$is the rotation matrix for rotation around the *y*-axis. This yields:

where ${o}_{x},\text{\hspace{0.17em}}{o}_{y},\text{\hspace{0.17em}}{o}_{z},\text{\hspace{0.17em}}{e}_{x},\text{\hspace{0.17em}}{e}_{y},\text{and}{e}_{z}$ are the magnitudes of the components of the normalized polarization vectors.

With the dispersion curves (Eq. 2) and polarizations (Eqs. 4 & 5) of waves in the bulk material in hand, we can simplify the expressions in Eq. 3 for $\stackrel{\rightharpoonup}{E}(y)$. For bound solutions, we require that the electric field decays to zero as $y\to \pm \infty $, so the terms in the cladding that are exponentially growing can be discarded. We now apply the symmetry condition that there is a reflection plane down the center of the sample, which eliminates half of the coefficients. In this situation, the solution must be made of symmetric and antisymmetric modes. Absorbing some constant factors into the coefficients, we have:

Symmetric:

Antisymmetric:

Applying the symmetry conditions eliminated half the unknowns, so now we need only apply boundary conditions at one interface to solve for the coefficients. The boundary condition is that the tangential *E* and *H* fields must be continuous across the boundary [20]. Using Faraday’s law and the fact that $\partial /\partial z=0,$
$\partial /\partial x=i\beta ,$ and $\partial /\partial t=-i\omega $ for our functional form, $\stackrel{\rightharpoonup}{E}(x,y,z,t)=\stackrel{\rightharpoonup}{E}(y)\mathrm{exp}[i(\beta x-\omega t)]$, we can express all the boundary conditions in terms of the electric field components:

The constant coefficients in functional form of the solutions (Eqs. 6 & 7) must be chosen so the above boundary conditions are satisfied at the interface ($y=\ell $). They must be solved independently for the symmetric and antisymmetric modes. The four expressions above yield a set of homogeneous equations which can be recast in matrix notation.

Symmetric:

Antisymmetric:

The polarizations (Eqs. 4 & 5) and bulk dispersion curves (Eq. 2) can be used to remove all dependence on $\alpha ,\text{\hspace{0.17em}}{\kappa}_{o},\text{\hspace{0.17em} and}{\kappa}_{e}$, so for a given angle *θ*, the only variables are *β* and *k*. The determinant will be zero, i.e. the set of equations has a solution, only for *β* and *k* pairs that are on the waveguide dispersion curve, and finding all allowed pairs traces out these curves. Using the allowed pairs, the bulk dispersion curves, and one additional “normalization condition” such as ${B}_{1}+{B}_{2}=1$, all wave vectors and coefficients can be completely determined. The theoretical dispersion curves for several angles are plotted along with the experimental data in Fig. 3, selected electric field profiles are shown in Fig. 4, and the effective indices of refraction for two angles are shown in Fig. 5.

## Acknowledgements

The authors would like to thank Alexei Maznev for his invaluable insights regarding anisotropic waveguide behavior and solutions. This work was supported by the National Basic Research Program of China (2010CB933801), the 111 Project (B07013), the National Natural Science Foundation of China (10604033), the Tianjin Natural Science Foundation (09JCYBJC15100), U.S. National Science Foundation Grant No. ECCS-0824185, and a National Science Foundation Graduate Research Fellowship (C.A.W.).

## References and links

**1. **Y.-S. Lee, *Principles of Terahertz Science and Technology* (Springer, New York, 2009).

**2. **R. M. Koehl, S. Adachi, and K. A. Nelson, “Real-Space Polariton Wave Packet Imaging,” J. Chem. Phys. **110**(3), 1317–1320 (1999). [CrossRef]

**3. **N. S. Stoyanov, D. W. Ward, T. Feurer, and K. A. Nelson, “Terahertz polariton propagation in patterned materials,” Nat. Mater. **1**(2), 95–98 (2002). [CrossRef]

**4. **T. Feurer, J. C. Vaughan, and K. A. Nelson, “Spatiotemporal coherent control of lattice vibrational waves,” Science **299**(5605), 374–377 (2003). [CrossRef] [PubMed]

**5. **T. Feurer, N. S. Stoyanov, D. W. Ward, J. C. Vaughan, E. R. Statz, and K. A. Nelson, “Terahertz Polaritonics,” Annu. Rev. Mater. Res. **37**(1), 317–350 (2007). [CrossRef]

**6. **T. P. Dougherty, G. P. Wiederrecht, and K. A. Nelson, “Impulsive Stimulated Raman Scattering Experiments in The Polariton Regime,” J. Opt. Soc. Am. **9**(12), 2179–2189 (1992). [CrossRef]

**7. **Y. X. Yan, E. B. Gamble, and K. A. Nelson, “Impulsive Stimulated Scattering: General Importance in Femtosecond Laser Pulse Interactions with Matter, and Spectroscopic Applications,” J. Chem. Phys. **83**(11), 5391–5399 (1985). [CrossRef]

**8. **D. Auston and M. Nuss, “Electrooptic Generation and Detection of Femtosecond Electrical Transients,” IEEE J. Quantum Electron. **24**(2), 184–197 (1988). [CrossRef]

**9. **Y. S. Lee, T. Meade, V. Perlin, H. Winful, T. B. Norris, and A. Galvanauskas, “Generation of narrow-band terahertz radiation via optical rectification of femtosecond pulses in periodically poled lithium niobate,” Appl. Phys. Lett. **76**(18), 2505–2507 (2000). [CrossRef]

**10. **J. Hebling, A. G. Stepanov, G. Almási, B. Bartal, and J. Kuhl, “Tunable THz pulse generation by optical rectification of ultrashort laser pulses with tilted pulse fronts,” Appl. Phys. B **78**, 593–599 (2004). [CrossRef]

**11. **K. L. Yeh, M. C. Hoffmann, J. Hebling, and A. Keith, “Generation of 10 μJ ultrashort terahertz pulses by optical rectification,” Appl. Phys. Lett. **90**(17), 171121 (2007). [CrossRef]

**12. **K.-H. Lin, C. A. Werley, and K. A. Nelson, “Generation of multicycle THz phonon-polariton waves in a planar waveguide by tilted optical pulse fronts,” Appl. Phys. Lett. **95**(10), 103304 (2009). [CrossRef]

**13. **D. H. Auston, K. P. Cheung, J. A. Valdmains, and D. A. Kleinman, “Cherenkov Radiation from Femtosecond Optical Pulses in Electro-Optic Media,” Phys. Rev. Lett. **53**(16), 1555–1558 (1984). [CrossRef]

**14. **J. K. Wahlstrand and R. Merlin, “Cherenkov radiation emitted by ultrafast laser pulses and the generation of coherent polaritons,” Phys. Rev. B **68**(5), 054301 (2003). [CrossRef]

**15. **N. S. Stoyanov, T. Feurer, D. W. Ward, and K. A. Nelson, “Integrated diffractive terahertz elements,” Appl. Phys. Lett. **82**(5), 674–676 (2003). [CrossRef]

**16. **P. Peier, S. Pilz, F. Müller, K. A. Nelson, and T. Feurer, “Analysis of phase contrast imaging of terahertz phonon-polaritons,” J. Opt. Soc. Am. B **25**(7), B70–B75 (2008). [CrossRef]

**17. **Q. Wu, C. A. Werley, K. H. Lin, A. Dorn, M. G. Bawendi, and K. A. Nelson, “Quantitative phase contrast imaging of THz electric fields in a dielectric waveguide,” Opt. Express **17**(11), 9219–9225 (2009). [CrossRef] [PubMed]

**18. **C. A. Werley, Q. Wu, K. H. Lin, C. R. Tait, A. Dorn, and K. A. Nelson, “A comparison of phase sensitive imaging techniques for studying THz waves in structured LiNbO_{3},” J. Opt. Soc. Am. B **27**(11) 2350-2359 (2010). [CrossRef]

**19. **E. A. Bahaa, B. Saleh and M. C. Teich, *Fundamentals of photonics* (JOHN WILEY&SONS, 1991)

**20. **M. Born and E. Wolf, *Principles of Optics (Cambridge University Press, Cambridge**,**1999**).*

**21. **D. H. Auston and M. C. Nuss, “Electrooptic generation and detection of femtosecond electrical transients,” IEEE J. Quantum Electron. **24**(2), 184–197 (1988). [CrossRef]

**22. **A. S. Barker Jr and R. Loudon, “Dielectric properties and optical phonons in LiNbO_{3},” Phys. Rev. **158**(2), 433–445 (1967). [CrossRef]

**23. **N. S. Stoyanov, T. Feurer, D. Ward, E. Statz, and K. Nelson, “Direct visualization of a polariton resonator in the THz regime,” Opt. Express **12**(11), 2387–2396 (2004). [CrossRef] [PubMed]

**24. **S. Wang, M. L. Shah, and J. D. Crow, “Wave propagation in thin film optical waveguides using gyrotropic and anisotropic materials as substrates,” IEEE J. Quantum Electron. **8**(2), 212–216 (1972). [CrossRef]

**25. **D. P. G. Russo and J. H. Harris, “Wave propagation in anisotropic thin-film optical waveguides,” J. Opt. Soc. Am. **63**(2), 138–145 (1973). [CrossRef]

**26. **W. K. Burns and J. Warner, “Mode dispersion in uniaxial optical waveguides,” J. Opt. Soc. Am. **64**(4), 441–446 (1974). [CrossRef]

**27. **V. Ramaswamy, “Propagation in asymmetrical anisotropic film waveguides,” Appl. Opt. **13**(6), 1363–1371 (1974). [CrossRef] [PubMed]

**28. **S. Nemoto and T. Makimoto, “Further discussion of the relationship between phase and group indices in anisotropic inhomogeneous guiding media,” J. Opt. Soc. Am. **67**(9), 1281–1283 (1977). [CrossRef]

**29. **D. Marcuse, “Modes of a symmetric slab optical waveguide in birefringent media-Part I: Optical axis not in plane of slab,” IEEE J. Quantum Electron. **14**(10), 736–741 (1978). [CrossRef]

**30. **D. Marcuse and I. P. Kaminow, “Modes of a Symmetric Slab Optical Waveguide in Birefringent Media Part II: Slab with Coplanar Optical Axis,” IEEE J. Quantum Electron. **15**(2), 92–101 (1979). [CrossRef]