Imaging of THz waves in 2D photonic crystal structures embedded in a slab waveguide

We present space- and time-resolved simulations and measurements of single-cycle terahertz (THz) waves propagating through two-dimensional (2D) photonic crystal structures embedded in a slab waveguide. Specifically, we use a plane wave expansion technique to calculate the band structure and a time-dependent finite-element method to simulate the temporal evolution of the THz waves. Experimentally, we measure the space–time evolution of the THz waves through a coherent time-resolved imaging method. Three different structures are laser machined in LiNbO3 crystal slabs and analyzing the transmitted as well as the reflected THz waveforms allows determination of the bandgaps. Comparing the results with the calculated band diagrams and the time-dependent simulations shows that the experiments are consistent with 3D simulations, which include the slab waveguide geometry, the birefringence of the material, and a careful analysis of the excited modes within the band diagrams.


Introduction
The interaction of electromagnetic waves with randomly or periodically structured media has gained considerable interest, especially when the structure sizes become comparable with the wavelength. Important examples are photonic crystals (PhCs), which are composite materials with periodic variations of the electromagnetic properties in one dimension (1D), 2D or 3D [1]. In the past, 1D PhC structures have been used extensively, for example as dielectric multilayer mirrors [2]. The propagation of electromagnetic waves through PhCs may be described by applying the band theory known from solid state physics [3]- [5] and Yablonovitch et al proposed that the dispersion relation splits into several photonic bands analogous to the electronic bands in semiconductors. In between the bands, photonic bandgaps may exist in 1D, 2D or 3D depending on the wavelength, the propagation direction, the polarization of the wave and the dielectric properties of the periodic structure [6,7]. So far, artificially created PhCs of dielectric materials have been employed to observe a variety of fascinating phenomena. On the most fundamental level, PhCs offer the possibility to control the linear [6] as well as the nonlinear behavior [8,9] of electromagnetic waves propagating within them; of special interest in this context is the tailoring of linear and nonlinear properties of PhC fibers [10]. In integrated optics, PhCs have often led to better performance of existing devices, such as bends [11] or splitters, or have even facilitated devices with novel functionalities. PhCs have also been used as multidimensional resonators, which in turn are employed for example to localize photons within defects of the photonic structure [4], to inhibit the spontaneous emission of excited atoms [3], to construct a micro-cavity within a periodic dielectric structure [12], or to realize a nano-scale laser [13]. A variety of powerful numerical techniques have been developed to simulate wave propagation through those media [14]- [18]. However, a few experimental techniques are appropriate to measure the electromagnetic field distribution with sufficient spatial and temporal resolution. In the visible part of the spectrum, near-field microscopy combined with interferometry is often employed [19]- [24]. Ideally, the spatial resolution should 3 be a fraction of a wavelength and the temporal resolution a fraction of an oscillation period. A spectral regime where these requirements can be readily fulfilled is the terahertz (THz) part of the spectrum. Moreover, the fabrication process of the appropriate PhC structures is greatly simplified as the relevant length scales are of the order of tens of microns [25,26]. In the recent past, significant progress has been made in THz science and applications, such as high-bandwidth signal processing, THz imaging, THz sensing or THz spectroscopy [27]. Also PhCs were studied in this frequency range. For example, Kitahara et al analyzed the dispersion of THz radiation in 2D PhCs; they measured the far-field amplitude and phase of the transmitted electric field by THz time-domain spectroscopy in order to reconstruct the dispersion relation of the PhC structure [28]. Liu et al numerically calculated the band structures of PhCs and the field distribution within functional devices, such as waveguides or beam splitters [29]. Recently, an integrated platform for THz waveform generation, guidance, processing and characterization has been demonstrated and has led to THz polaritonics [30]. In the THz polaritonics platform, the signal carriers are phonon-polaritons, coupled admixtures of electromagnetic waves and polar lattice vibrations of comparable frequency and wavevector. They are generated by femtosecond optical pulses through a second-order nonlinear optical process [31,32] in suitable crystals, such as LiNbO 3 . They are fully coherent, they can be almost single cycle and they propagate with light-like speeds. Their wavelengths are of the order of tens of micrometers and their oscillation periods are about picoseconds. The shape and the frequency distribution of these THz waves can be controlled using appropriate beam shapes in space and time [33,34]. The THz waves propagate primarily in lateral directions away from the optical excitation beam. This key property enables interactions of the THz waveforms with further optical pulses, functional elements, or with structured media. The THz waves induce a change in the crystal's refractive index, whose sign and magnitude reflect the THz electric field amplitude. A time-delayed probe beam illuminating the whole crystal will therefore experience a phase modulation where it overlaps with the THz wave. This property allows for spatiotemporal imaging of THz waves. That is, snapshots showing the spatial distribution of the electric field can be taken at any time after excitation. Since the probe pulses have visible wavelengths and durations of approximately 100 fs, the THz waves can be visualized with very high spatial and temporal resolutions.
Here, we employ the integrated THz polaritonics platform to characterize 2D PhCs embedded in slab waveguides for THz applications [35]. This paper is organized as follows. Section 2 summarizes the basic properties of the THz polaritonics platform and briefly introduces the simulation tools used to calculate the band diagrams and the time-dependent THz wave propagation through the structured materials. We recapitulate that the band structure critically depends on the geometry and material properties. We show that the shape and the wavevector spectrum of the single-cycle THz wave determine which of all possible bands are excited and, thus, they influence the positions and widths of the bandgaps observed. In order to emphasize this point, we exemplarily calculate the band structure and the THz wave propagation for a hexagonal array of circular air voids in LiNbO 3 , first in 2D and then for the same structure in a slab waveguide, and finally we incorporate the birefringence of LiNbO 3 . In section 3, we present the experimental methods, and in section 4, the experimental results for three different structures. Specifically, we measure the reflected and the transmitted THz waveform through spatiotemporal imaging. A comparison with the simulations shows that good agreement is found between simulated and measured bandgaps only if, firstly, the waveguide properties, secondly, the birefringence and, thirdly, the excited mode spectrum are considered. Phonon-polariton excitation in a thin, structured LiNbO 3 crystal. Inset: the focused pump beam generates two counter-propagating THz waves. The wavevector k THz is split into an in-plane k x and a forward k y component.

Theory and simulations
In the following sections, we address the theoretical aspects required to model the experimental results presented towards the end of the paper. We start with a brief summary of THz wave generation, propagation and detection in suitable crystals followed by a summary of the simulations. Specifically, we define the problem space and increase the complexity step by step from a 2D to a 3D problem and, lastly, include the birefringent nature of the material. We use a plane wave method (BandSolve, RSoft) to calculate the band diagram of the structures and we use a finite-element method (FEM) to simulate the time-dependent propagation of THz waves through the structured materials. Although we measure the transmitted and the reflected THz waveforms, in the theory section we compare the in-structure THz waveform of the FEM simulations with the band structure of the plane wave expansion method because of the higher sensitivity in identifying discrepancies between the two. The comparison of the plane wave expansion with the FEM results is used to determine which of the bands in the band diagram are actually excited by the incident THz waveform; this in turn influences the bandgaps observed in reflection or transmission.

THz phonon-polaritons
2.1.1. Generation, propagation and detection. Phonon-polaritons travel at light-like speeds with partly mechanical and partly electromagnetic energy. At the low wavevectors and THz frequencies of interest here, the phonon-polaritons exhibit an almost linear dispersion described by ω = ck/n(ω), where the refractive index n(ω) is essentially a constant.
In our experiments phonon-polaritons are excited through optical rectification by focusing a short laser pulse to a line as shown in figure 1. From the line focus two almost plane singlecycle THz waves propagate in opposite directions with a non-negligible forward wavevector component k y [31]. The polarization of the THz waves is parallel to the line focus and parallel to the optical axis of the crystal. When the THz waves reach the back surface they are totally reflected as the angle of incidence exceeds the critical angle for total internal reflection. As a consequence, the forward wavevector component reverses its sign at each crystal-air interface.
Spatiotemporal probing of the THz waves relies on the electro-optic effect. The THz electric field induces a change of the crystal's refractive index at optical wavelengths according to n opt (x, y, z, t) ∝ E z (x, y, z, t). A time-delayed probe pulse experiences a phase modulation integrated along its direction of propagation y which appears to be stationary as long as the probe pulse is phase-matched to the forward wavevector component [30]. The integrated phase modulation can be converted to an amplitude modulation using a Sagnac interferometer [36]. In our setup, the probe beam typically illuminates the whole crystal and images the complete spatial electric field distribution at a given time delay. That is, the intensity measured by a CCD camera is directly proportional to the electric field of the THz radiation, i.e. I CCD (x, z, t) ∝ E z (x, y = 0, z, t). To analyze the images, we integrate them along one spatial direction, here along the z-direction, thereby creating a 1D plot I CCD (x, t). Recording a sequence of those plots for increasing time delays and stacking them on top of each other result in a 2D space-time distribution that shows the electric field amplitude as a function of x and t. A 2D Fourier transformation converts the space-time distribution to a wavevector-frequency distribution and all nonzero amplitudes should coincide with allowed bands in the band diagram of the material.

Time-dependent simulations.
To simulate the time-dependent excitation and propagation of the THz electromagnetic field, we use the FEM toolbox COMSOL; specifically, we use the transient propagation solver with two or three spatial dimensions depending on the problem. COMSOL requires four essential inputs: firstly, the spatially varying electromagnetic properties of the medium; secondly, the appropriate boundary conditions; thirdly, the definition of the discrete mesh; and, lastly, the solver parameters, such as time-step size. The THz wave generation process was approximated by defining an appropriate source at the left boundary such that the emerging THz wave matches the experimentally observed one. All simulations in three spatial dimensions include the forward wavevector component.
To demonstrate that it is essential to include the forward wavevector component in 3D simulations, we analyze the simple case of a slab waveguide. The thickness of the waveguide is 50 µm, it is 4 mm long, and waveguide theory predicts the appearance of several waveguide modes. Here, the electric field is s-polarized with respect to the crystal-air interfaces and the corresponding modes split into two subsets, namely even and odd modes with the symmetry of the lowest, zeroth-order mode being even. Figure 2 shows that all modes are in between the two light lines; that is, their effective index of refraction is between 1.0 (air) and 5.1 (LiNbO 3 ). Figure 2(a) shows that only the lowest order waveguide mode is excited when the THz source at the left boundary has no forward wavevector component. To obtain the results shown in figure 2(b), the forward wavevector component is included and in that case also the first-order mode is excited. In other words, the wavevector spectrum of the THz waveform with forward wavevector component is such that it excites all accessible waveguide modes irrespective of symmetry, whereas the THz source without forward wavevector component can only excite the symmetric modes within reach. Consequently, a correct interpretation of the experiments requires simulations that include the forward wavevector component.
Further, we analyze, again with the help of the simple slab waveguide, to what extent the probing process must be simulated in detail. As mentioned above, probing uses a second-time-delayed visible pulse that propagates through the entire slab and integrates the THz-induced phase variations along the y-axis. In the simulations, probing can be incorporated by an appropriate postprocessing of the stored 3D electric field distribution. For example, to obtain the measured 2D electric field distribution at a time delay of t, we extract N x z-slices from the 3D data that are stored at times t + j y/v g , i.e. E z (x, j y, z, t + j y/v g ), with the group velocity of the probe pulse being v g , j ∈ [0 . . . N − 1] and the thickness of the waveguide being N y. Summing all N x z-slices yields a 2D x zelectric field distribution that is proportional to the measured images. It turns out that the resulting 2D electric field distribution is approximately the same as in any x z-slice not too close to the center or the crystal-air interfaces of the slab waveguide. This is because phase-matching plays almost no role in such thin crystal slabs. The result of a simulation that includes the probing process is shown in figure 2(c). It is practically identical to the field distribution in figure 2(b), which displays the field distribution in a single x z-slice. Since all the experiments following this were recorded in 50 µm thick waveguides, we will henceforth always analyze the electric field distribution in a single x z-slice rather than simulating the entire probing process in order to minimize memory usage and computation time.
The Bloch-Floquet theorem for periodic eigenvalue problems states that H ( r ) = e i k r H n, k ( r ) is a solution to equation (1) and the functions H n, k satisfy The discrete eigenvalues ω n ( k) are continuous functions of k forming the bands with band number n when plotted versus the latter. The solutions are periodic functions of k as well, i.e. the solution at k is the same as the solution at k + G j , with R i · G j = 2πδ i j . It can be shown that a certain periodicity of the permittivity can lead to a photonic bandgap, i.e. a range of frequencies ω for which no waves with real wavevectors k exist. A detailed discussion can be found for example in [1,7]. Mid-gap frequencies and widths of the bandgaps of the 1D structure can be calculated using analytical models. The band diagrams of 2D and 3D periodic structures have to be calculated with special software tools, for example making use of the plane wave expansion technique [37,38]. In the frequency domain, the fields H k ( r ) are expanded in a complete basis, i.e. H k ( r ) = n h n b n ( r ), leading to a matrix eigenvalue problem for all h n . Here, we use a software package (BandSolve) from RSoft to calculate the eigenvalues and the band diagrams.  To calculate the band diagram we apply periodic boundary conditions, i.e., the structure is expanded to an infinite 2D array of elementary cells. The polarization is TM, i.e. the magnetic field vector is normal to the plane. The calculations are performed for a selected range of wavevectors within the Brillouin zone, i.e. those connecting the points of high symmetry, , M and K. The wavevector path is indicated by red lines in figure 3(b). BandSolve calculates the frequencies f n ( k) for a user defined number of bands as a function of the wavevectors specified. The output of a band calculation consists of a frequency versus wavevector plot with the symmetry points of the PhC highlighted and shows the eigenvalues of the different bands including the bandgaps between them. We performed a parametric scan of the period and/or the diameter of the air voids in order to shift the bandgaps in the accessible frequency range of our experiment, i.e. between 0.1 and 1.5 THz. The band diagram for an optimized hexagonal structure is shown in figure 4; the hole diameter is d = 50 µm and the periodicity is x p = 100 µm. We calculate 20 different bands and find two bandgaps within the frequency range of interest.

2D hexagonal lattice.
While the band calculations are based on an infinitely large PhC, the time-dependent FEM simulations calculate the THz waveform propagation through a PhC with 10 × 10 elementary cells (as used in the experiments). From the space-time plots E z (x, t) we obtain E z (k x , f ) through a 2D Fourier transformation. The THz waveform is generated at the left boundary, propagates first through an unstructured part of the LiNbO 3 slab and impinges on the structured area where parts of the waveform are reflected and other parts are transmitted through the PhC and appear to the right of the PhC. The resulting space-time plot of E z (x, t) is shown in figure 5.
The plot may be divided into three distinct regions. To the left of the PhC, we find the incident and the reflected waveforms, within the PhC a complex structured waveform, and to the right of the PhC the transmitted waveform. A 2D Fourier transform of the region within the PhC yields the wavevector-frequency distribution E z (k x , f ), and its amplitude |E z (k x , f )| can be compared with the band structure obtained with BandSolve. In the simulations as well as in the experiments, we rotate the structure by 90 • rather than changing the wavevector from k x to k z in order to extract E z (k z , f ). In figure 6, the in-structure electric field amplitudes are PhC refl. trans.
x (mm) t (ps)  Some of the weak features of the in-structure amplitude |E z (k x , f )| are a result of the space-time window function used in combination with the Fourier transformation; we employ either a Gaussian or a Kaiser-Bessel window function. The finite width of all features when compared with the bands (lines) is determined mostly by the finite THz spectrum. The asymmetry of the in-structure amplitude |E z (k x , f )| comes from the fact that the incident THz waveform originates to the left of the PhC. As a consequence, the amplitudes of those parts of the band diagram with positive velocities are more intense than the corresponding parts with negative velocities. A striking feature is that two of the bands (dashed red lines in figure 6(a) are not excited at all even though they are within the accessible frequency range. Figure 6(c) shows the corresponding scan along the k z direction and its direction and range are indicated in figure 6(d). Again, features of those parts traveling towards positive k z are more intense than those traveling in the opposite direction and again two bands are not excited at all. Ochiai et al [39] and Robertson et al [40] have shown that further symmetry aspects determine whether the incoming THz wave can couple to a specific mode or not. Those modes that are not excited here have additional symmetry properties either with respect to x or z that prohibit an interaction with the incoming THz wave. An important consequence of this effect is that the upper limit of the lowest bandgap in the k x direction shifts towards higher frequencies. In other words, a band calculation alone is not always sufficient to determine the positions and widths of the bandgaps observed in the experiment.

2D hexagonal lattice in an isotropic slab waveguide.
A more sophisticated simulation of the experiments has to account for the finite thickness of the slab waveguide in which the 2D PhCs are embedded, in our case 50 µm. The additional confinement of the THz radiation in the third dimension leads to a modified wave propagation and, consequently, to a modified band diagram. To account for the additional boundaries, we have to extend the simulation to three dimensions. The whole slab that includes the PhC volume is embedded in a rectangular simulation box with a refractive index of 1 (air). In BandSolve the unit cell is extended in the y-direction with a total height of four times the slab thickness. The Brillouin zone is 3D but the wavevector paths considered here remain in the x z-plane.
The calculations yield a significantly different band diagram as shown in figure 7 when compared with the purely 2D situation in figure 4. We show only those bands with even symmetry with respect to the y-direction and only those parts of the bands lying in between the light lines and thus corresponding to guided modes. We find two bandgaps and they range from 0.51 to 0.61 THz and from 0.81 to 0.89 THz. Most importantly, the bandgaps shift to different frequencies, which is caused by the modified effective index of refraction of the slab waveguide embedded in air. Again we overlay the results of 3D FEM simulations. The simulation of a structure consisting of 10 × 10 unit cells in a 50 µm thick slab waveguide would require a huge amount of grid points and of computer resources. These can be reduced quite substantially when we exploit the symmetries of the system by defining a simulation box as shown in figure 8.
The PhC volume is reduced along the z-direction to one unit cell assuming an infinite periodicity along this axis. Along the x-axis the PhC is sandwiched between two unstructured pieces of LiNbO 3 . THz waves with their electric field vector along the z-axis are generated   (2). The 3D time-dependent solver calculates the 3D electric field distribution E z (x, y, z, t) for progressively later times t. From these we extract an x z-slice (see section 2.1.2) and integrate over the z-direction resulting in E z (x, t). A Fourier transform of the in-structure field distribution then yields E z (k x , f ) and its amplitude is included in the band diagram. Figures 9(a) and (c) show the band diagram together with |E z (k x , f )| and |E z (k z , f )|, respectively. The corresponding wavevector paths are indicated in figures 9(b) and (d).
Again we find excellent agreement between the time-dependent simulations and the even bands of the band diagram. Most of the energy of the THz wave couples to the even modes, whereas very little energy couples to the odd modes. Although some bands are not excited at all (drawn in red and for reasons explained above [39,40]) the bandgap widths are not affected by this. The total k x range covers more than one Brillouin zone as indicated in (b). (c) Band diagram and in-structure amplitude |E z (k z , f )|. The total k z range covers more than one Brillouin zone as indicated in (d).

2D hexagonal lattice in an anisotropic slab waveguide.
The final refinement is to include the birefringent nature of the host material. With the optic axis of LiNbO 3 parallel to the z-axis, the permittivity tensor of LiNbO 3 is diagonal, with the elements x x = 41, yy = 41 and zz = 26 [30]. The orientation of the PhC with respect to the anisotropy of LiNbO 3 is different for the k x and the k z scans as we rotate the crystal rather than the incident THz wave. That is, for a k z scan the diagonal elements must be interchanged, i.e. x x = 26, yy = 41 and zz = 41. A complete wavevector scan connecting the points of high symmetry is shown in figure 10.
The most striking features are a shift of the two lowest bandgaps along the frequency axis and the appearance of a third bandgap. The frequency shift arises from the strong birefringence of LiNbO 3 at THz frequencies and the bandgaps are now found between 0.47 and 0.54 THz, between 0.78 and 0.82 THz and between 1.08 and 1.10 THz.
Firstly, we focus on the in-structure region as was done in the previous sections. The amplitude of the electric field in wavevector-frequency space is compared with the band diagram as shown in figure 11.
The wavevector paths continue through more than a single Brillouin zone and we find excellent agreement between the two simulations. Several bands are not excited at all and are highlighted by red dashed lines. In contrast to the previous simulations, the fact that some bands  cannot be excited strongly influences the bandgaps. For example, the upper boundary of the first bandgap in figure 11(a) depends on whether the band (red dashed curve) is excited or not. Also, the excitation of the second red dashed band in figure 11(a) determines whether an additional bandgap in this frequency range will appear or not. Secondly, we examine the reflected and the transmitted parts of the THz waveform as these will eventually be compared with the experimental results. Some of the PhC properties can be extracted from these regions as well, especially the positions and widths of the bandgaps. Since the reflected as well as the transmitted parts of the THz waveform propagate through unstructured parts of the crystal before being measured, we expect that the amplitude |E z (k x , ω)| contains information on the PhC as well as on the waveguide modes [41]. To show that this is true we crop the relevant parts of the space-time distribution before Fourier transformation. The amplitude |E z (k x , ω)| of the transmitted THz waveform is shown in figure 12(a) and that of the reflected THz waveform is shown in figure 12(b).
In both cases the amplitude is nonzero only along the allowed waveguide modes. The dashed red lines denote the edges of the bandgaps derived from the bands having eigenmodes with even parity in the y-direction. The bandgaps are best seen in the transmitted part of the lowest order waveguide mode; only this mode has the same symmetry as the bands. We omit plotting the odd bandgaps, because of the small amount of energy transported on bands with this symmetry. The bandgaps are much harder to identify in the reflected part of the THz waveform because of the reduced contrast. While the transmission coefficient drops to zero for frequencies within the bandgaps, the reflection coefficient increases from the Fresnel value, here approximately 0.6, to 1.

Summary and comparison.
The last few subsections have illustrated that the bandgaps of a periodic hexagonal lattice of air voids in an LiNbO 3 slab waveguide depend on the geometry, the birefringence of LiNbO 3 , and on the wavevector spectrum of the incident THz wave. For comparison, figure 13 summarizes the lowest order bandgaps within the relevant frequency range for the three different types of simulations, i.e. purely 2D, 3D in an isotropic and in an anisotropic waveguide. The wavevector is parallel to k x .
The bandgaps are derived from the band diagrams of even order bands and consider only those bands that are excited by the incident THz wave. The differences are substantial and a comparison to the experimental results should indicate whether the most sophisticated  simulations are based on the correct assumptions. Similar results are found when analyzing the bandgaps in the perpendicular direction.

THz generation and probing
The intense laser pulses used for THz wave generation originate from a kHz Ti:sapphire multipass amplifier and have a pulse duration of approximately 100 fs. A beam splitter is used to split the laser beam into a probe and a pump beam. The pump beam is focused with a cylindrical lens to a line into the LiNbO 3 crystal near the PhC structure (see figure 1) where two counterpropagating near plane wave THz wavelets are generated. Their center frequency is around 0.75 THz and their bandwidth ranges from 0.1 to 1.2 THz. The less intense probe beam propagates through a Sagnac interferometer. The LiNbO 3 crystal is located close to the beam splitter of the Sagnac interferometer and is imaged onto a CCD camera. Due to this asymmetry, the reference pulse in the interferometer passes through the crystal about 1 ns before the probe pulse. The time delay between the probe and the pump is adjusted with a motorized delay stage. Details of the detection setup can be found in [36].

Li N bO 3 micro-machining
The 2D structures in the slab waveguides were prepared by laser micro-machining of LiNbO 3 crystals having dimensions of 5 mm × 5 mm × 0.05 mm. The laser source was an amplified Ti:sapphire laser system generating pulses with a maximum energy of 1 mJ at a repetition rate of 1 kHz. A polarizer and a half-wave plate were used to reduce the pulse energy to the optimal value and the compressor was deliberately 'misaligned' to achieve a pulse duration between 1 and 2 ps. An external Pockels cell and a pair of polarizers were used to pass the required number of pulses to the samples. The pulses were focused through an objective with a numerical aperture of 0.15 and a focal length of 21 mm, resulting in a spot size of 5 µm at the sample. The LiNbO 3 samples were placed on a computer-controlled, motorized three-axis translation stage. The setup allowed for fabrication of arrays of uniform circular holes with a minimum diameter of approximately 5 µm. Larger holes were produced by milling and the difference between entrance and exit diameters was found to be of the order of a few microns only. A selection of samples is shown in figure 14.

Results
We start by analyzing THz wave propagation in 50 µm thin slab waveguides. Then we proceed with hexagonal and cubic 2D photonic structures integrated in parts of the same waveguides. All the LiNbO 3 slabs were clamped in a custom-built mount providing protection and high stability.

Planar waveguide
First, we analyze the propagation of THz waves through unstructured thin crystals, i.e. thin slab waveguides [42]. The measured space-time plots for slab thicknesses of 50 and 150 µm are shown in figures 15(a) and (c).
For small pump-probe time delays, a nearly plane single-cycle wave emerges from the excitation region. Its small forward wavevector component leads to a reflection of the THz waveform at the back surface of the crystal. That is, the THz waveform is reflected back and forth between the lower and the upper crystal-air interface and follows a zigzag path. Phase matching between the THz waveform and probe beam is best if the former propagates in the  figure 2(b)) and shows that the tilted THz waveform indeed excites all even and odd order modes, which are within its frequency spectrum.

Hexagonal structure I
The geometry is identical to that in section 2.1.2 and is shown in figure 14(b). Figure 16 shows the measured space-time distribution of the THz electric field. The THz wave is excited at t = 0 ps and x = 0 mm to the left of the PhC structure, which is indicated by the two vertical dashed black lines. One waveform propagates towards the left and is eventually reflected at the left crystal edge. The other waveform propagates towards the right Trans.

Refl. from edge
THz field amplitude (a.u.) Figure 16. Space-time THz electric field distribution for a hexagonal PhC structure, which is located between x = 0.2 mm and 1.2 mm as indicated by the two vertical dashed black lines. and impinges at the PhC at t = 4.5 ps and at x = 0.28 mm. Part of the waveform is reflected and reverses its direction of propagation. Its shape suggests that there are multiple scattering events originating from different rows of air holes. There is no signal visible within the PhC, because there is some debris deposited on the remaining LiNbO 3 during laser-machining, effectively making the area between the air holes opaque. To the right of the PhC structure the transmitted THz waveform appears at around 15 ps. Most of the waveform is located below the dashed white line that corresponds to the THz velocity within bulk LiNbO 3 . That is, the velocity within the PhC structure must be considerably higher, i.e. the corresponding effective index of refraction must be lower than in bulk, which is no surprise because the structured region is composed 23% of air. Both the transmitted and the reflected THz waves are no longer single-cycle waves but rather complex wave trains reflecting the frequency filter characteristics of the PhC. Specifically, the transmitted signal shows a clear signature of a chirped instantaneous frequency, the long wavelength components appear first and are followed by the high-frequency components. A 2D Fourier transform of the relevant parts of the space-time distribution helps analyzing the properties of the PhC structure. To that end we crop the region of interest, i.e. the transmitted or reflected waveform, and calculate E z (k x , f ). The amplitudes of the transmitted and reflected waveforms are shown in figure 17; the gray scale is logarithmic. The FEM simulations suggest that we almost exclusively excite even order modes, which is why the red dotted lines here and in all following figures of the same kind indicate the boundaries of the bandgaps with even symmetry only. While the first two bandgaps result from the band diagram, the highest bandgap appears only because the FEM simulations indicate that some of the modes cannot be excited (see figure 11(a)).
The wavevector-frequency distribution agrees well with the simulations in figure 12 and the maxima coincide closely with the slab waveguide modes. The transmitted signal is shown in figure 17(a). While most of the energy is concentrated in the first waveguide mode (even), there is a small amount of energy in the second waveguide mode (odd). The first mode shows clear signatures of the PhC's bandgaps and the first and the second bandgap are very pronounced. Integrating the wavevector-frequency distribution with respect to the wavevector yields the transmitted power spectrum, which is shown to the right of the wavevector-frequency distributions in figure 17 and which shows clear signatures of the bandgaps. The reflected signal reveals a strong contribution on the second (odd) waveguide mode, which is indicative of a strong reflection of the energy contained in this mode, whereas the better part of the lowest order (even) mode is transmitted. The only signature of the (even) bandgaps in the integrated reflected power spectrum is that the amplitude is somewhat higher than in the other parts of the power spectrum. Below the first bandgap the two parts, i.e. transmitted and reflected signals, show a pronounced complementarity. Almost all of the signal in this frequency range is transmitted, and close to nothing is reflected. In order to measure a k z wavevector scan, another sample was fabricated where the structure was rotated by 90 • and the wavevector-frequency distributions are shown in figure 18. The simulations predict two bandgaps for even symmetry within the frequency range of interest and their borders are indicated by the dotted red lines.
While the transmission data in figure 18(a) indicate that most of the energy is concentrated on the lowest order mode, the reflected signal as shown in figure 18(b) is almost equally distributed among the two lowest order modes. Analyzing the transmitted signal in figure 18(a) and the corresponding wavevector integrated power spectrum clearly shows the existence of two bandgaps. The majority of the lowest bandgap forms a complete bandgap, i.e. it appears in the k x scan from to K (see above) as well as in the k z scan from ro M. It ranges from roughly 0.5 to 0.6 THz and agrees well with the measurements. Figure 19 compares the measured transmitted power spectrum along k x with the calculated bandgaps resulting from three different models, namely 2D, 3D and 3D plus anisotropy. The best agreement when analyzing the position and width of the lowest and broadest bandgap is found with the most appropriate model (labeled '3d n eo /n o , sym.') that accounts for the correct dimensionality, the anisotropy of the crystalline material, and the excited modes. Sym.

Hexagonal structure II
A further experiment was performed using a hexagonal structure with a different geometry. The air-hole diameter and its period were increased to 90 and 150 µm, respectively (see figure 14(d)). The array consists of 10 × 10 unit cells. The measured wavevector-frequency distributions of the transmitted and the reflected signals together with the waveguide modes and bandgaps are shown in figure 20.
Simulations predict that for this geometry the bandgaps should shift towards lower frequencies and their widths should change, in excellent agreement with the experimental observations as shown in figure 20(a). The overall strength of the transmitted signal is much lower when compared with the previous structure, especially at higher frequencies. The band diagram hints that the higher lying bands are rather flat, which means that the higher frequency components travel at much reduced velocities through the structure. This process not only slows down but may even hinder the propagation of higher frequency components, causing reduced high frequency transmission. Nevertheless, we find the lowest bandgap located between 0.4 and 0.5 THz and the second bandgap roughly between 0.7 and 0.75 THz. The calculations indicate that the second bandgap exists for both parities even and odd with respect to y and it is located just below the lowest allowed frequency of the second waveguide mode. The third bandgap cannot be determined from the data shown in figure 20(a); however, we find enhanced reflection within the third bandgap as shown in figure 20(b). Similar to the previous structure the reflected power spectrum shows peaks in the reflectivity curve at frequencies where bandgaps are located, but clear signatures of those bandgaps are obscured.

Cubic structure
The last structure examined is a cubic array of air holes in LiNbO 3 , as shown in figure 14(a). BandSolve calculations show that such symmetry exhibits no complete bandgap. Nevertheless, there exist bandgaps for specific directions, for example along the k x -direction. The measured wavevector-frequency distributions are shown in figure 21.
In the transmitted signal ( figure 21(a)) the electric field amplitude is almost exclusively concentrated on the lowest order waveguide mode and the three lowest bandgaps can unambiguously be identified. The bandgaps are also clearly visible in the power spectrum shown to the right of figure 21(a). The reflected power spectrum is rather structured with a high reflectivity in the bandgap regions and a lower reflectivity for frequencies that are allowed to propagate through the PhC structure. Almost all frequency components below the first bandgap are completely transmitted through the structure, suggesting that the effective refractive index of the corresponding band is almost identical to the surrounding LiNbO 3 crystal.

Conclusion
Within the polaritonics platform we have demonstrated that it is possible to coherently measure the reflected and the transmitted components of a THz single-cycle plane wave impinging on a finite size 2D PhC structure embedded in a slab waveguide. The photonic structures were fabricated through laser micro-machining of thin crystals of LiNbO 3 . The measurements yield the space-time distribution of the THz waveform before, in (in principle) and after the PhC and can be used to calculate the wavevector-frequency distributions. Those can then be compared to calculations resulting from a plane wave expansion method. The complementary timedependent FEM simulations help to determine which of the many bands of a given symmetry are actually excited in an experiment and thus they facilitate identification of the various bandgaps. The experimental results, i.e. the positions and widths of the different bandgaps, can only be explained through 3D simulations including the birefringence of the LiNbO 3 crystalline material, and the agreement found is excellent. Pure 2D or 3D simulations in an isotropic crystal fail to explain the experimental results. The in-structure propagation was simulated and analyzed, but up to now it could not be observed, in experiments, with sufficient quality owing to debris in the structured areas. So far, the analysis yields only a qualitative and pictorial understanding of the electric field distribution within the PhC structures, as shown in figure 22.
Although the quality of the data is insufficient to extract quantitative information, it shows that in principle such an in-structure field measurement with the possibility to extract almost the entire band diagram is within reach.