## Abstract

We present a simplified coupled mode theory (CMT), suited for high losses, to describe ultra-broadband THz generation through optical rectification (OR) of *fs* infrared pulses in waveguides. We derive a new expression that incorporates loss effects into the coherence length for OR. The simplified approach reproduces the results of a computationally rigorous integral CMT that must be used for broadband THz generation. With the new model we perform a parametric study to establish the optimal conditions for OR in symmetric, five-layer, metal/cladding/core structures with electro optic polymer cores. We find conversion efficiencies as high as 35 × 10^{−4} W^{−1} and bandwidths up to 20 THz when pumping at 1900 nm. We find that low-loss-cladding layers enhance the efficiency for phase-matched structures, increase the interaction length, and improve the stability of the efficiency with respect to variations in waveguide parameters.

© 2013 Optical Society of America

## 1. Introduction

Recently we published a frequency domain (FD) coupled mode theory (CMT) [1] to describe the generation of broadband THz pulses by mixing the frequency components of infrared (IR) femtosecond (*fs*) pulses through difference frequency generation (DFG). Since *fs* pulses have a broadband spectrum, all the possible frequency pairs ($\omega -\Omega /2$,$\omega +\Omega /2$) in the pump spectrum drive the DFG of a single THz component, *Ω*, and all these contributions must be accounted when modeling broadband THz pulses. This formulation works even in the situation where the DFG occurs between several IR and THz modes. Since it is a FD theory, it accounts for distinct absorption features and higher order dispersion effects observed in the waveguide composites at both the THz and IR frequencies. These features must be incorporated in the refractive index and absorption coefficient of the waveguide composites when solving for the waveguide modes. Unfortunately, the model requires an *integral approach* that is computationally expensive and is tedious to implement.

In this paper, we present a simplified, yet accurate, FD *three-wave-approximation* approach to describe the ultra-broadband generation of THz pulses through optical rectification (OR) of *fs*-IR pulses in waveguide emitters. This simplified approach is suited for highly lossy waveguides and maps the convoluted DFG integral interaction to a simpler three-wave interaction. The generation of the THz frequency, *Ω*, can be thought of as resulting from the difference frequency interaction, (*ω _{0} + Ω/2*)

*-*(

*ω*), of the two main frequency components of a

_{0}- Ω/2*fs-*IR pulse with center frequency,

*ω*. We derive a closed expression for the power density spectra,

_{0}*dP*, for situations in which OR occurs between the

_{THz}*l-*th IR mode and the

*m-*th THz mode. We also provide, up to first order in dispersion effects, a very simple and intuitive expression for the loss-dependent complex phase-mismatch term, Δ

*β*, for OR between these two modes. The resulting equations are simple and depend only on the THz frequency

*Ω*, the pump central frequency

*ω*, and, the IR and THz modes involved in the interaction.

_{0}The *three-wave-approximation* presented here can be seen as an extension of the plane wave formulation presented by Vodopyanov [2] to model OR of *fs* IR pulses in crystals but is specifically suited for waveguide architectures [3]. A novelty of our approach, different from [2, 3], is that our CMT formulation is equipped to deal with lossy modes [1, 4], unavoidable for most THz devices. Phase matching ultra-broadband THz emitters is more complicated than its narrowband counterparts [4]. The whole spectrum has to be considered and several THz frequencies must be simultaneously phased-matched to achieve high THz conversion efficiency. For broadband THz applications the total power, *P*_{THz}, is non-trivially distributed across the THz frequencies and *dP _{THz}* is an important quantity containing all the spectral information of the generated THz pulse. However, the expressions derived in the simplified model presented here can only be used when a single-IR-mode and a single-THz-mode interact. Care must be taken since situations where multiple IR modes (or THz modes) are required, are not necessarily rare. For these cases, the full CMT is required [1].

Using our new approach, we derive an expression for the coherence length that *incorporates loss effects* to the OR process, but is also relevant to second harmonic generation (SHG). The coherence length calculated through this method is *always finite* and is valid for structures that are not completely phased matched. Importantly, we introduce the concept of effective attenuation for the OR process, *α _{eff}*(

*Ω,ω*), which depends on the THz effective-mode-attenuation,

_{0}*α*(

_{m}*Ω*), and also on the properties of the IR effective-mode-attenuation,

*α*(

_{l}*ω*). This effective loss concept could be generalized to other nonlinear processes (e.g., SHG) and are potentially useful for quasi-phase matched architectures in which the optimal interaction length is increased by surpassing phase-matching limitations. In those situations it is important to correctly account for the effect of guiding-related-losses on the nonlinear conversion.

_{0}We apply our technique to find the optimal pumping conditions for a symmetric, five layer, metal/cladding/core waveguide emitter with a core composed of a poled electro-optic (EO) polymer with a high EO coefficient of 130 pm/V AJTB203 [5], low-THz-loss polystyrene (PS) cladding layers, and Al capping layers. For this geometry OR is mediated through even TM IR modes and even TM THz modes. We compute the modes of the slab waveguide using transfer matrix theory and implement a simple numerical scheme to solve for the transcendental equation that yields the loss dependent coherence length. We then perform a parametric study for different waveguide parameters and pumping wavelengths, in the telecom and SWIR region, to establish under what conditions this structure results in high conversion efficiencies. Coupling conditions are optimized to guarantee that most of the incident pump power (~%80) is utilized in the interaction. We find a maximum nonlinear conversion efficiency of 35 x 10^{−4} Watts^{−1} for a structure pumped at 1900 nm, seven times higher than efficiencies reported in other models [6]. In addition, the bandwidth is almost 20 THz, for a pump laser with a 47*-fs* pulsewidth. Several structures pumped at the commercial telecom wavelength near 1567 nm also show broad bandwidth (~8 THz) and high efficiencies peaking at 7 x 10^{−4} W^{−1}.

We discovered that the fundamental TM IR mode hybridizes to a plasmon mode as the cladding layers are removed and the THz generation occurs through OR of the first excited even TM IR mode (*l = 1)* that behaves as a quasi-fundamental mode. Even this simple metal/EO-core/metal structure shows high conversion efficiencies with this mode. We predict an enhancement in efficiency (up to a factor of two at 1900 nm) for structures with cladding layers when the phase-mismatch of AJTB203, *n _{g,IR} – n_{THz}*< 0.02, which occurs for pumping wavelengths > 1700 nm. In the limiting case of thin cladding layers (~1

*µm*) the

*l = 0*and

*l = 1*IR modes must both be accounted for to correctly describe the OR process. The cladding layers help improve the stability to variations of the waveguide parameters.

## 2. Theory

In FD coupled mode theories, the solution to the nonlinear interaction is expanded in terms of the unperturbed waveguide modes of the structure, $\tilde{E}(\overrightarrow{r},\omega )={\displaystyle {\sum}_{l}{A}_{l}(z){\widehat{e}}_{l}({r}_{\perp},\omega )exp(i{\beta}_{l}(\omega )z)}$. Here *z* and ${r}_{\perp}$ are respectively the longitudinal and transverse directions, *ω* is the frequency and, ${\widehat{e}}_{l}$, ${\widehat{h}}_{l}$ and *β _{l}(ω) = k_{0}N_{eff,l}(ω) + iα_{l}(ω)/2* are respectively the complex electric and magnetic mode profiles, and the mode propagation constants of the lossy modes in our model. Using a theory based on the Lorentz reciprocity theorem, valid for systems with high losses, coupled first order partial differential equations can be derived describing the interaction [1]

*l*and

*l’*IR modes, ${{\rm K}}_{l,l\text{'}}^{m}={\displaystyle \int dA\left({\widehat{e}}_{m,\perp}-{\widehat{e}}_{m,z}\right)\cdot {\chi}^{(2)}:{\widehat{e}}_{l}{\widehat{e}}_{{l}^{\prime}}{}^{*}}$ measures the overlap of the modes involved in the nonlinear interaction and depends on the symmetry class and geometry of the nonlinearity, ${\widehat{e}}_{l,\perp}$ and ${\widehat{e}}_{l,z}$ are the respective transverse and longitudinal components of the mode profile, ${P}_{l}(\omega )=2{\displaystyle \int d{r}_{\perp}\widehat{z}\cdot ({\widehat{e}}_{l,\perp}\times {\widehat{h}}_{l,\perp})}$, and ${A}_{l}^{IR}(\omega )={\displaystyle \int d{r}_{\perp}\widehat{z}\cdot ({\tilde{E}}_{in}(\omega )\times {\widehat{h}}_{l})/2}$is the coupling coefficient of the input IR pump to the guiding IR modes [1].

Equation (1) holds *even when multiple IR modes are involved in the DFG interaction* and is rigorous provided the waveguide modes are complete [1]. The resulting expressions for *A _{m}*(

*z,Ω*) involve

*integral equations*over the IR pump frequencies,

*ω*, and have to be evaluated for each THz frequency,

*Ω*, in the output spectrum. This is required for broadband THz generation since it is necessary to account for the contributions of all the possible frequency pairs $(\omega -\Omega /2,\omega +\Omega /2)$ in the spectrum of the

*fs*-pulse pump that drives the DFG of a single THz component

*Ω*. However, this integral approach is computationally demanding.

For OR of *fs* IR pulses, it is often common to model systems where the IR pump beam is coupled to a single IR mode, *l = l’*, and a single THz mode, *m*. In this case, the main contribution in Eq. (1) results from the DFG interaction of the pair of IR frequencies centered at, ${\omega}_{0}$, the center frequency of the pump beam. Here, the complex phase-mismatch term associated with OR becomes $\Delta \beta ({\omega}_{0},\Omega )={\beta}_{l}({\omega}_{0}+\Omega /2)-{\beta}_{{l}^{\prime}}{({\omega}_{0}-\Omega /2)}^{*}-{\beta}_{m}(\Omega )$. Expanding this about ${\omega}_{0}$ to first order in dispersion effects, gives,

Here, ${n}_{}{}_{g,eff}={N}_{eff,l}+\omega (d{N}_{eff,l}/d\omega )$, is the group effective index of the *l*-th IR mode and ${N}_{eff,m}$is the effective index of the THz mode. In the Appendix, we provide a derivation of Eq. (2) and for the THz spectral power density, *dP _{THz}*, for OR between the

*l*-th IR mode and the

*m*-th THz mode,

*β =*Δ

*β*Δ

_{R}+ i*β*.

_{i}*C*(Ω) depends on the pump beam optical properties and coupling condition, an expression for which is given in the Appendix. All information of the spectral content of the generated radiation is contained in

*dP*and${P}_{THz}={\displaystyle \int d{P}_{THz}(z,\Omega )d\Omega}$ is the total THz power. A close inspection of Eq. (3) reveals some further insight into the role of both the real and imaginary parts of the phase-mismatch, Δ

_{THz}*β*

_{R}and Δ

*β*

_{i}, in the OR process. We can define, ${\alpha}_{eff}(\Omega )={\alpha}_{m}(\Omega )+\Delta {\beta}_{i}(\Omega )$, as the effective attenuation for the OR process. If the IR losses are negligible, ${\alpha}_{eff}(\Omega )\sim {\alpha}_{m}(\Omega )/2$ and the effective OR attenuation is dominated by THz mode attenuation; the THz radiation is continuously generated along the entire waveguide, thus the effective THz losses are reduced by a factor of two. The OR effective attenuation,

*α*does not only depend on

_{eff}(Ω),*α*, but also depends on the mode attenuation of the IR mode,

_{m}(Ω)*α*(

_{l}*ω*), through Δ

_{0}*β*

_{i}. This is important in situations where

*α*(

_{l}*ω*) is high and

_{0}*α*has a big contribution resulting from the IR modes. Note that Δ

_{eff}(Ω)*β*can also be reduced if

_{i}*α*(

_{l}’*ω*) is negative and has a high magnitude, then

_{0}*α*can be

_{eff}*< α*.

_{m}/2The coherence length, *L _{c}*, is an important figure of merit in the design of emitters based on any type of nonlinear process since it gives the optimal emitter length for maximum conversion. The coherence length results from $d(d{P}_{THz}(z,\Omega ))/dz|{}_{{L}_{c}}=0$ [7], which yields

In general, since Eq. (4) is transcendental it must be solved it numerically. In simple cases, an analytic expression for the coherence length can be derived from Eq. (4). For example, in the case where no losses are considered, Δ*β _{i} =* 0 and we find that $\Delta {\beta}_{R}{L}_{c}=\pi $ which yields, ${L}_{c}=c/(2\Omega |{n}_{g,eff}-{N}_{eff,m}|)$ a generally accepted expression [8]. Although very useful, this expression does not account for losses and yields infinite

*L*for completely phase-matched structures, which is unphysical. In the other limiting scenario with completely phase-matched structures, Δ

_{c}*β*= 0, and we find ${L}_{c}=2\mathrm{ln}({\alpha}_{m}/\Delta {\beta}_{IR})/({\alpha}_{m}-\Delta {\beta}_{IR})$ which agrees with the result in [4].

_{R}To test the simplified three-wave-approximation described here in Eqs. (2) and (3), we apply the model to the optimized structures reported in [1] and computed with the full integral DFG approach that solves Eq. (1). The results for dP_{THz} in Fig. 1(a), obtained through the three-wave-approximation approach and computed with same parameters used in [1], appear to be identical to the ones Fig. 5(a) of [1]. Nevertheless, some slight differences can be observed, mostly for devices with longer emitter lengths *L* where the values produced by the three-wave-approximation overestimate dP_{THz}. To study the difference in both methods we show in Fig. 1(b) the ratio of the THz power density spectra, dP_{Full}/ dP_{Approx.}, obtained by dividing dP_{THz} from the full integral DFG approach by the value obtained through the three-wave-approach. Note that the ratios in Fig. 1(b) go to zero for increasing THz frequencies since dP_{THz} obtained through the full integral DFG theory goes to zero faster than for the three wave approximation approach. Since we only account for first order dispersion effects in the three-wave-approach, *comparing both methods provides a way to study high order dispersion effects* such as group velocity dispersion (GDV) and others, which are naturally included in the full CMT of [1]. In Fig. 1(b) we see how GVD and higher order dispersion effects degrade the THz generation for higher THz frequencies and for longer emitter lengths. GVD effects cause pulse broadening of the IR pump in the waveguide, this reduces the pump’s bandwidth, which in turn reduces the number of frequency pairs $(\omega -\Omega /2,\omega +\Omega /2)$ in the spectrum that contribute to a single THz frequency. This is more detrimental for higher THz frequencies as seen in Fig. 1(b). Finally, in Fig. 1(c) we show the difference between both methods in calculating the nonlinear THz conversion efficiencies, ${\eta}_{THz}(L)={P}_{THz}/{P}_{pump}^{2}$. The results differ considerably (~20%) for devices with emitter lengths are greater than ~1 mm. Larger effects will be seen for materials with higher dispersion curvature.

## 3. Results and discussion

Poled EO guest-host polymer composites have been synthesized to possess very high EO coefficients at telecom wavelengths and they constitute an inexpensive alternative for efficient emitters based on second order nonlinearities. A symmetric, five-layer, slab waveguide Fig. 2(a) has the ideal geometry for exploiting the EO polymer nonlinearity [1, 9]. In this configuration, the IR pump is coupled to a TM mode that is bounded inside the EO core and the THz radiation is collected by the single fundamental TEM-like mode that extends spatially in the entire core/cladding section of the waveguide. The low-loss-cladding layers help mitigate both the high THz losses introduced by the lossy core and the metal capping layers that help confine the modes involved in the interaction, which increases the mode overlap. For the cladding layers we chose polystyrene (PS) since it presents low THz losses [10] that help reduce the THz mode attenuation and has adequate index properties suited for phase-matching. For the metal capping layers we used Al. A poled polymer guest-host EO composite with guest chromophore AJTB203 mixed at 24% mass concentration in an amorphous polycarbonate (APC) host has been measured to have an *r _{33}* ~130 pm/V at 1330 nm [5]. The corresponding susceptibility value ${\chi}_{xxx}^{(2)}$ at IR wavelengths is modeled assuming a two-level model [11]; AJTB203 has a single absorption peak at 803 nm, Fig. 2(b).

We followed a strategy to design waveguides with single guided TM modes in the (0.1-10) THz frequency range but allowed the device to be multi-mode in IR frequencies. That sets an upper limit to the waveguide sizes of about 20 µm thick (*t = 10 µm*). There is a single symmetric guided mode at frequencies for which the wavelength is greater than the waveguide size. Special care has to be taken to maximize the pump coupling to the IR mode that drives the OR. This differs from the approach in [1,9] that resulted in optimized devices with single IR modes with thin layers (~1 µm) that are difficult to build, pump, and handle. Thus, fixing *t = 10 µm* in Fig. 2(a) results in waveguides with thicker cores that guarantee that more chromophores contribute to the THz generation process. In addition, it is now possible to couple higher input powers (> 50 mW) since we are not required to focus tightly into the core while maintaining IR intensities below the damage threshold.

In Fig. 3(a) and its inset we show the respective real and imaginary part of the calculated modal dispersion, *β _{m}(ω) = ω/c N_{eff,m}(ω) + iα_{m}(ω)/2*, for the fundamental (

*m = 0*) TEM-like mode at THz frequencies, Ω, for 7 structures (out of the 200 we computed) with

*t = 10 µm*and $d\in [5\mu m,t]$. In the Appendix, we explain how the modes were calculated. In Fig. 3(b) we show the IR effective group index,

*n*, vs. core half-thickness,

_{g,eff}*d*, for structures with

*t = 10 µm*, for the fundamental (

*l = 0*) and first excited (

*l = 1*) TM even IR modes for three different pumping wavelengths,

*λ*(1330, 1567, 1900) nm. We assumed a symmetric input beam injected in the middle of the waveguide so only even modes can get excited. Better phase-matching is obtained when

_{0}=*n*matches

_{g,eff}*N*and, since decreasing the cladding thickness (bigger

_{eff,m}*d*’s) increases the modal effective index, thinner cladding layers result in a better phase-matched structure. Unfortunately, as seen in the inset of Fig. 3(a), the THz mode attenuation increases too since the core is also more absorptive and that is also reflected in the OR effective attenuation,

*α*. Comparing Figs. 3(a) and 3(b) we see that, for the materials in this study, better phase matching is achieved for higher pumping wavelengths. It is impossible to perfectly phase-match more than a single THz frequency at 1567 nm and none at 1330 nm, while excellent broadband phase-matching occurs for 1900 nm.

_{eff}In Fig. 4, we show the effect of losses on the coherence length, *L _{c}*, for OR between the

*l = 0*IR mode and the TEM-like (

*m = 0*) THz mode for structures with

*t = 10 µm*. In Fig. 4(a), we show

*L*, computed using both the traditional approach (no loss) and our numerical approach (including loss) solving Eq. (4) for the structures that showed the higher

_{c}*L*values when pumping at 1567 nm and 1900 nm. The numerical approach always yields a finite

_{c}*L*since it accounts for losses. On the other hand, the traditional method ${L}_{c}=c/(2\Omega |{n}_{g,eff}-{N}_{eff,m}|)$ neglects loss effects and yields overestimated optimal emitter lengths that result in unphysical singularities in

_{c}*L*for completely phase-matched structures. In Fig. 4(a) we observe how both methods to compute

_{c}*L*differ significantly for the structure with

_{c}*d ~9.1 µm*when pumping at 1900 nm, this occurs since several distinct THz frequencies get phase-matched in the (0.1-20) THz range. In contrast, both methods have a fair agreement for poor phase-matching scenarios at high THz frequencies, differing significantly only around 1 THz where the structure with

*d ~9.5 µm*gets phase-matched when pumping at the commercial telecom wavelength 1567 nm.

Including the effect of losses when computing coherence length, *L _{C}*, provides additional criteria for choosing the right waveguide parameters that maximize the THz output. The coherence length is then strictly the emitter length that maximizes the THz OR conversion at a given THz frequency. The behavior of

*L*for different waveguide parameters can also be used for design purposes since structures that present high

_{c}*L*are either well phase-matched or have low effective mode attenuation, or both. Thus, from Fig. 4(b) we would expect the structure with

_{c}*d ~9.5 µm*to show relatively efficient and broadband THz output (>5 THz) for emitter lengths

*L*~1 mm when pumped at 1567 nm. For the structures in Fig. 4(b) the coherence length decreases monotonically for increasing frequencies but, it peaks when

*d ~9.5 µm*and

*L*

_{C}> 1 mm for THz frequencies <5 THz. Similarly, from Fig. 4(c) we would expect that structures with

*L ~2*mm and

*d*~9.1 µm to have an efficient and ultra-broadband (>10 THz) THz output when pumped at 1900 nm. Note also that at 1900 nm

*d*can be chosen to enhance a certain THz frequency range, in Fig. 4(c) emitters with

*d ~8 µm*have high

*L*values in 0-4 THz range, for example.

_{c}In Figs. 5(a) and 5(c) we show the THz power density spectrum, *dP _{THz}*, for the structures mentioned in the previous paragraph for different emitter lengths,

*L*, pumped either at 1567 nm or 1900 nm. These results, computed using Eqs. (2) and (3) for OR between the

*l = 0*IR mode and the TEM-like THz mode, confirm our observations about the coherence lengths plots in Fig. 4. The most optimal configuration is indeed the 1.5 mm long structure with

*d = 9.1 µm*pumped at 1900 nm, shown in Fig. 5(c). This structure shows about 3 times more spectral THz signal, in the (0.1-20) THz range, than the 600 µm long structure with

*d = 9.5 µm*in the (0.1-8) THz range at 1567 nm, shown in Fig. 5(a). This translates to a factor of six in the nonlinear THz conversion efficiency, ${\eta}_{THz}(z)={P}_{THz}/{P}_{pump}^{2}$with ${P}_{THz}={\displaystyle \int d{P}_{THz}(z,\Omega )d\Omega}$, in Fig. 5(d). The dips observed in the THz power spectra for all structures are caused by the absorption features of AJTB203, Fig. 3(a).

The 35 x 10^{−4} W^{−1} efficiency results mainly due to the high nonlinearity of the EO polymer AJTB203 and the much broader bandwidth attained using the highly phase-matched structure with *d = 9.1 µm* at 1900 nm. Several structures pumped at 1567 nm also present high efficiencies peaking at 7 x 10^{−4} W^{−1} and it is even possible to enhance narrow bands for lower THz frequencies by increasing the emitter length, notice the 1 mm long structure with *d = 9.5 µm* in Fig. 5(a). To optimize the THz output spectra around a narrow THz frequency range we can increase the cladding layers without losing much in efficiency (Fig. 5(d)). In Fig. 5(b), the structures with *d = 8.6 µm* can be tuned to have highly peaked narrow bands of THz output tunable in the (0.1-2) THz range by changing *L* from 1 mm to 2 cm.

In Fig. 5(e), we show how the efficiency, *η _{THz}*, is affected by varying both the cladding layers (

*t-d*) and the emitter length,

*L*, for structures pumped at 1567 nm. We observe a broad zone of high efficiencies > 5 x 10

^{−4}Watts

^{−1}for structures with

*d ~8.9 ± 1 µm*and

*L*$\in $(

*0.5-4*) mm. It is interesting to note that OR from the

*l = 0*IR mode is suppressed for structures with very thin claddings. This does not imply that a metal/EO-core/metal will not efficiently generate THz, however. Indeed, as seen in Fig. 5(f), OR between the

*l = 1*IR mode and the TEM-like THz mode for structures

*with thin claddings*and shorter emitters length (

*L*~1 mm) is almost as efficient as for the optimized structures for OR with

*l = 0*IR mode in Fig. 5(e), but, is null for structures with thick claddings. The coupling conditions were optimized for each structure with different

*d,*for either the

*l = 0*or

*l = 1*IR mode, choosing the beam waist of the input Gaussian pump that maximizes the coupling efficiency to the EO core up to ~80%. The resulting beam waists vary in the (4-10)

*µm*range depending on core sizes. This is experimentally achievable with an elliptical focusing lens. For the calculations we assumed a train of 47

*fs*Gaussian pulses at a 100 MHz repetition rate and 10 mW of input power.

Different from the five-layered-slab waveguide in Fig. 2(a), the metal/EO-core/metal structure (no cladding) has a plasmon mode in its mode decomposition see Fig. 6(a). As the cladding layers are thinned, the fundamental *l = 0* mode hybridizes to the plasmon mode while the first excited *l = 1* even IR mode assumes the role of a quasi-fundamental symmetric state still with two crossings at zero (Fig. 6(a)). The plasmon mode is not suited to mediate OR under the conditions we assume here. For that reason, OR for the metal/EO-core/metal structure is mediated only through the *l* = 1 even IR mode. The efficiency for the plasmon is low since: (i) the maximum coupling efficiency is reduced by an order of magnitude in the hybridization process since the pump beam is assumed to be launched at the center of the core and the radiation is confined near the metal-core boundary for the plasmon, (ii) the overlap integral with the TEM-like THz mode, ${K}_{l,l}^{m}$, is also reduced for the plasmon since the THz mode profile is quasi-uniform inside the core and (iii) the *l =* 0 IR mode becomes less phase-matched as the cladding layers are removed since the effective group index, *n _{g,eff}*, increases to its highest value. On the other hand,

*n*is reduced for the

_{g,eff}*l*= 1 even IR mode, which makes the

*l = 1*mode better phased matched for OR for structures with no cladding layers(Fig. 3(b)).

OR through the plasmon mode is also heavily affected by the high IR losses resulting from proximity to the metal. In Fig. 6(b), we plot the OR effective THz attenuation, *α _{eff}*, and its different contributions for the plasmon (

*l = 0*) IR mode for a metal/EO-core/metal structure pumped at 1900 nm. From Eq. (2) the imaginary part of the phase-mismatch is Δ

*β*Δ

_{i}=*β*

_{IR}

*– α*, and since

_{m}/2*α*=

_{eff}*α*,

_{m}+ Δβ_{i}*α*=

_{eff}*α*Δ

_{m}/2 +*β*

_{IR}. Thus, the two main contributions to

*α*result from the THz effective mode attenuation

_{eff}*α*and the IR mode contribution to imaginary part of the phase mismatch, $\Delta {\beta}_{IR}={\alpha}_{l}({\omega}_{0})[1-{\Omega}^{2}/(4{\omega}_{0}{}^{2})]+{\Omega}^{2}{{\alpha}^{\prime}}_{l}({\omega}_{0})/(4{\omega}_{0})$; which depends on the THz frequency but whose value remains relatively constant across the (0-20) THz range since ${\Omega}^{2}/(4{\omega}_{0}{}^{2})$<<1. For metal/EO-core/metal structures Δ

_{m}*β*

_{IR}is relatively high for the

*l = 0*IR mode since the high plasmonic IR losses make

*α*significantly larger than

_{eff}*α*. In Fig. 6(c) we show

_{m}*Δβ*, as the cladding layers are thinned, $d\to t$,for the

_{IR}*l*= 1 even IR mode (for the

*l =*0 in the inset) for different pumping wavelengths,

*λ*. For cladding layers > 1µm the IR losses for the IR modes are small,

_{0}*α*

_{0}(

*ω*

_{0})<

*α*

_{1}(

*ω*

_{0}), and

*α*for both modes. However, as the claddings are thinned Δ

_{eff}~α_{m}/2*β*

_{IR}increases rapidly for the

*l = 0*IR mode as it turns plasmonic when

*α*

_{0}(

*ω*

_{0}) reaches its maximum, inset Fig. 6(c). On the other hand for the

*l = 1*mode, Δ

*β*

_{IR}increases to a maximum for a structure with thin but non-zero cladding layers,

*d ~9.8µm*, where the

*l = 1*mode starts assuming the role of the fundamental (Fig. 6(c)); After which Δ

*β*

_{IR}decreases to a small but non-zero value as the cladding layers are completely removed. In the absence of cladding layers, Δ

*β*

_{IR}and thus

*α*are much larger for the

_{eff}*l*= 0 IR mode than for the

*l*= 1 even IR mode. Thus, the

*l*= 1 even IR mode is more suited to mediate the OR.

The cladding layers reduce the loss effects added by the confining metal layers to the OR effective attenuation, *α _{eff}*, for both the

*l = 0*and

*l = 1*IR modes. However, we see in Figs. 5(e) and 5(f), similar maximum efficiencies for structures with or without cladding layers. Is it possible to design structures with cladding layers that enhanced OR efficiencies?

We performed a parametric study, computing the efficiency contour plots Figs. 5(e) and 5(f) for both the *l = 0* and *l = 1* IR mode, for 10 wavelengths across the telecom and mid IR ranges in (1300-1900) nm; the results of the study are summarized in Fig. 7. For each wavelength 800 structures with 40 distinct *d* $\in $(5-t) µm values and 20 distinct *L*$\in $(0.1-10) mm values were modeled. Similar trends as the ones observed in Figs. 5(e) and 5(f) are seen for all wavelengths. However, the maximum efficiency, the position of this maximum in the *(d, L)* plane, and the range of parameters for which the efficiency remains relatively high or *stable* are different for each wavelength. For every *d* there is always an *optimal emitter length, L _{opt}*, for which the efficiency peaks to a value

*η*(

_{max}= η_{THz}*L*). In Fig. 7(a) we show

_{opt}*η*for different core thickness (or cladding layers) when the OR is mediated through the

_{max}*l =*0 or

*l*=

*1*IR mode when pumping either at 1567 nm or 1900 nm. We further define the

*overall maximum efficiency*,

*η*, as the maximum of

_{MAX}*η*(

_{max}*d*), indicated by the black stars and the red dots in Fig. 7(a). Note that for OR with the

*l*= 1 IR

*η*occurs when there are no cladding layers, i.e.

_{MAX}*η*or simply a metal/core/metal structure. In Fig. 7(b) we plot the overall maximum efficiency,

_{MAX}= η_{max}(d = t)*η*, vs. pumping wavelength,

_{MAX}*λ*for OR mediated through either the

_{0},*l =*0 or

*l*= 1 IR mode.

The result is that adding low-loss-cladding layers to a metal/core/metal structure results in an efficiency enhancement for THz generation through OR in highly phase-matched scenarios, *λ _{0}* > 1700 nm see Figs. 3 and 7(b). For those wavelengths the phase-mismatch for the material used in the EO core (AJTB203) Δ

*n*≅

*n*~0.02, near 7 THz. For longer pumping wavelengths, Δ

_{g,IR}– n_{THz}*n*decreases resulting in enhanced efficiencies for OR through the (

*l*= 0) IR mode for structures with cladding layers. Thus, a minimal condition to enhance OR through the addition of cladding layers is to have core mismatches, Δ

*n*<0.02. We also observe that adding the cladding layers results in higher stability of THz conversion efficiency to variations in emitter length,

*L*, and core half thickness,

*d*. This is an important design consideration since it is directly related to fabrication tolerances. In Fig. 7(c) we plot the overall optimal emitter length,

*L*, the emitter length for which

_{OPT}*η*occurs, and,

_{MAX}*L*, the emitter length for which

_{80%}*η*has decreased to 80% of

_{THz}*η*. The gap between

_{MAX}*L*and

_{OPT}*L*is a measure of how stable the efficiency is to variations in emitter lengths. Similar

_{80%}*L*values are obtained by pumping either the

_{OPT}*l =*0 or the

*l =*1 but

*L*is considerably larger for structures with cladding layers where OR is mediated by the

_{80%}*l =*0 IR mode.

Finally, we briefly study OR for structures in the *thin cladding regime t-d* < 1µm. In Fig. 7(a) we observe that the *three-wave-approximation* predicts very low efficiencies for structures with thin claddings (~200 nm or *d*~9.8 µm) when the OR is mediated through either the *l* = 0 or the *l =* 1 IR mode. However, for structures in the thin cladding regime both the *l* = 0 and *l* = 1 IR modes are involved in the DFG interaction and must both be accounted for to completely describe OR.

In the *thin cladding regime* the modal characteristics of the structure are a hybrid of the mode decompositions of both the metal/core/metal and the five-layer structures. Thus, in the transition zone of the hybridization of the fundamental IR mode to a plasmon mode, neither the *l = 0* mode nor the *l* = 1 modes behave like a symmetrical bound state. In Fig. 8(a) we observe that both modes must be accounted for to obtain a correct expansion of a Gaussian input beam centered at the core of the waveguide. As a result, the three-wave-approximation is unsuited to model the OR for structures with thin claddings since it assumes a single IR mode interacting with a single THz mode. Indeed, in Fig. 8(b) we see that the efficiencies predicted with the three-wave-approximation for the *l = 0* and *l = 1* IR modes do not match the results predicted by the more accurate integral DGF theory [1]. Moreover, adding the efficiencies predicted by three-wave-approximation with both IR modes stills results in lower efficiencies than those predicted by the integral DFG method. Interestingly, the maximum efficiency predicted by the integral method for structures with thin claddings, 7 x 10^{−4} W^{−1}, is the same as those predicted with both modes in Fig. 7(b) at 1567 nm. In Fig. 8(c) we also see that a structure *L*_{opt} ~200 µm shows a broad bandwidth of 15 THz.

## 4. Conclusions

We present a simplified FD *three-wave-approximation* approach to describe the ultra-broadband generation of THz pulses through OR of *fs*-IR pulses in waveguide emitters. This simplified approach is suited for highly lossy waveguides and agrees with the rigorous DFG integral theory [1] for situations in which the interaction occurs between a single IR mode and a single THz mode. Comparing the results produced by both methods we are able to quantify GVD and higher order dispersion effects on OR in waveguide emitters. A great advantage of using the three-wave-approximation approach is the considerable savings in computer time for obtaining the results. With the three-wave-approximation approach it takes about half an hour to compute the 9 traces in Fig. 1(a), with 500 points per trace distributed across the 20 THz bandwidth. This is to be compared with 4.5 days that are required to compute the same results using the integral method [1] using a modern, multi-core desktop computer.

Of considerable importance is the development of a new expression for the coherence length that *incorporates loss effects* into the OR process and we show how analyzing the coherence length plots for different waveguide parameters would be helpful in designing efficient THz waveguide emitters. In general, structures with longer coherence lengths correspond to structures with higher efficiencies. However, to determine the optimal emitter length we need to compute the THz power density spectra, *dP _{THz}* , for different emitters lengths to find the structure with the higher performance. Studying the efficiency contour plots,

*η*like those in Figs. 5(e) and 5(f) allows one to isolate the regions of the waveguide parameters that produce the most efficient THz generation through OR by exciting a particular IR mode and pumping at a given wavelength. We also find that it is possible to have situations where the output THz pulses are enhanced for very narrow THz frequency ranges, which is useful for CW THz sources.

_{THz}(L,d),We performed a parametric study for symmetric, five layer, Al/PS/AJTB203 structures to find the most optimal conditions for high efficiency for 10 pumping wavelengths over the range, 1300-1900 nm. The 1.5 mm long structure pumped at 1900 nm showed the maximum nonlinear conversion efficiency of 35 x 10^{−4} W^{−1} and had a bandwidth of almost 20 THz, when pumped with a 47-*fs* pulse. Several structures pumped at 1567 nm also show high efficiencies, peaking at 7 x 10^{−4} W^{−1} and it is even possible to enhance narrow bands for lower THz frequencies.

Finally, we discovered that the fundamental TM IR mode hybridizes into a plasmon mode as the cladding layers are thinned and the THz generation occurs through OR of the first excited even TM IR mode (*l* = 1*)* that acts as a quasi-fundamental mode. This simpler Al/AJTB203/Al structure shows high conversion efficiencies for OR with the *l* = 1 mode. We observe an enhancement in efficiency (up to a factor of two at 1900 nm) for structures with cladding layers when the phase-mismatch of the material used in the core, *n _{g,IR} – n_{THz}* is < 0.02, which occurs for pumping wavelengths > 1700 nm. In the limiting case of thin cladding layers (< 1

*µm*) the

*l*= 0 and

*l*= 1 IR modes must both be accounted for to correctly describe the OR process. The cladding layers help improve the stability to variations of the waveguide parameters and adding these layers results in increased optimal emitter lengths (~2 mm).

## 5. Appendix

#### 5.1 Theory

The general expression for the THz power spectrum in the case of multiple-mode THz generation is [1]

*m*and

*m’,*are the THz modes and

*f*, is the repetition of the laser system. In the general case in which multiple IR modes are involved in OR, each term in the sum in Eq. (1) represents the interaction of at most three modes, two IR modes (

_{rep}*l*and

*l’*) interacting via DFG, and the

*m*THz mode collecting the generated radiation. If only the

*l*-th IR mode is excited, the sum in Eq. (1) collapses to a single term with

*l = l’*. A closed form for

*A*(

_{m}*z,Ω*) results from Eq. (1) by performing the integral in

*z*, neglecting pump beam depletion of the IR modes and assuming no THz at the input, assuming an input Gaussian beam as in [1], with a pulse duration, ${\tau}_{0}$, and central frequency, ${\omega}_{0}$. Then, ${A}_{l}^{IR}(\omega )\propto exp(-{\tau}_{0}{}^{2}{(\omega -{\omega}_{0})}^{2}/2)$and ${A}_{l}^{IR}(\omega +\Omega /2){A}_{l\text{'}}^{IR}{(\omega -\Omega /2)}^{\ast}\propto exp[-{((\omega -{\omega}_{0}){\tau}_{0})}^{2}+{(\Omega {\tau}_{0}/2)}^{2})]$. Thus, the integral in ω in Eq. (1) can be solved analytically noting that $exp(-{x}^{2}/{\u03f5}^{2})/\u03f5\sqrt{\pi}\to \delta (x)$ when $\u03f5\to 0$. Then,

Thus, the THz spectral power density, *dP _{THz}*, describing OR between a single-IR-mode

*l*and a single-THz-mode

*m*is,

Here $C(\Omega )=4\pi {f}_{rep}|{\u03f5}_{0}\Omega {K}_{l,l}^{m}{A}_{l}^{IR}({\omega}_{0}+\Omega /2){A}_{l}^{IR}{({\omega}_{0}-\Omega /2)}^{\ast}/{P}_{m}{|}^{2}\mathrm{Re}[{\displaystyle \int d{r}_{\perp}\widehat{z}\cdot ({\widehat{e}}_{m}\times {\widehat{h}}_{m}{}^{\ast})}]$and we split the phase-mismatch into both real and imaginary parts through,$\Delta \beta =\Delta {\beta}_{R}+i\Delta {\beta}_{i}$. The standard inner product used for the mode decomposition is not necessarily valid for lossy waveguides, the differential operator that gives the mode decomposition is non-Hermitian (non-self-adjoint) and its spectrum and modes are complex [12]. Thus, the power carried by the *m*-th THz mode accounted in the surface integral with *m’ = m* in Eqs. (5) and (7) must be strictly computed since we normalize the modes such that $\int d{r}_{\perp}\widehat{z}\cdot ({\widehat{e}}_{n}\times {\widehat{h}}_{n\text{'}})}=2{\delta}_{n,n\text{'}$ [4]. Finally, the phase-mismatch term, $\Delta \beta ({\omega}_{0},\Omega )={\beta}_{l}({\omega}_{0}+\Omega /2)-{\beta}_{{l}^{\prime}}{({\omega}_{0}-\Omega /2)}^{*}-{\beta}_{m}(\Omega )$, for OR between a single IR mode *l* and a single THz mode *m* is

Here, ${\kappa}_{eff,l}$ is the effective absorption coefficient of the *l*-th mode, ${\alpha}_{l}(\omega )/2={\kappa}_{eff,l}\omega /c$. It is possible to simplify Eq. (8) using this identity from a Taylor series of $f({\omega}_{0}\pm \Omega /2)$ around *ω _{0}*

The upper term in brackets corresponds to + in the left side of the equation and the lower to -. Replacing *f* by either ${N}_{eff,l}$ or ${\kappa}_{eff,l}$ results in an expression for the complex phase-mismatch for OR between one IR mode and a THz mode up to first order in dispersion effects

Here, ${n}_{}{}_{g,IR}={N}_{eff,l}+\omega (d{N}_{eff,l}/d\omega )$, is the group effective index of the *l*-th IR mode.

#### 5.2 Material optical properties and mode solving technique

The refractive index and attenuation at THz frequencies for AJTB203 were measured with a broadband THz spectrometer in the (0.1-10) THz range [10]. To model the absorptions features observed in AJT203 in the THz range the dielectric function in the THz range is modeled by a six-oscillator-Lorentz model

Here *n _{b}* is the background refractive index, Ω

*is the*

_{i}*i*-th resonance frequency and, Γ

*and*

_{i}*a*are the corresponding line widths and oscillator strengths respectively. A global fit was performed on the experimental data in [10] for THz index of refraction and attenuation. In Table 1 we report the parameters for the global fit and the curves are shown in Fig. 3(a) and its inset with dashed lines.

_{i}It was necessary to include a term *k*_{0}*ω* when fitting the THz attenuation due to the linear scattering effects that are typically observed in the extracted attenuation for EO polymer composites [10]. The fitting parameter found was *k*_{0} = 251.5 cm^{−1} THz^{−1}.

At IR wavelengths, we measured the refractive index of AJTB203 using the prism coupling technique and modeled it with a Sellmeier fit, ${n}^{2}=2.21122+0.346966\times {\lambda}^{2}/({\lambda}^{2}-{\lambda}_{0}{}^{2})$ with ${\lambda}_{0}=803\text{nm}$. We found that it was crucial for our calculation to include IR losses to avoid obtaining oversized optimal emitter lengths. Thus, we assumed a constant absorption coefficient for AJTB203, $\kappa =c\alpha /(2\omega )$, taking α = 1.8 dB/cm at 1500 nm. This attenuation value was extracted from the one measured for the EO polymer composite from Pacific Wave’s CLD-1, noting that both AJTB203 and CLD-1 have polycarbonate hosts [13]. To model the refractive index and absorption coefficients for PS and Al, both at THz and IR frequencies, we use the fits and values used in [1]. Since IR losses for PS were neglected in [1], here, we assumed a constant absorption coefficient for PS of 0.8 dB/cm at 1064 nm [13].

To model Al’s refractive index, *n*(Ω), and absorption coefficient, *k*(Ω), in the THz range we use the same Drude global fit as in [1]

Here *ε _{b}* is the background dielectric constant,

*δ*is the electron density in the metal,

*τ*is the scattering time,

*ε*is the permittivity of free space, and

_{0}*e*and

*m*are the charge and mass of the electron respectively. The fitting parameters for this global fit were

_{e}*ε*1754.4,

_{b}=*δ = 9.85 × 10*cm

^{22}^{−3}and

*τ = 12.6 fs.*

Due to the symmetry in the *y* direction the mode profiles accept an analytic expression for the symmetric, five-layer, slab waveguide in Fig. 2(a). We apply transfer matrix theory and mode matching methods to solve for the complex propagation constants [1, 14, 15]. The main challenge results from including losses in our model, requiring that the search for the solutions to the dispersion equation be carried out in the complex plane. Furthermore, to avoid missing any solutions and facilitating the search in the complex plane we implemented an algorithm designed to account for the ill-defined dispersion equation whose entire set of solutions lies on a two-sheeted Riemannian surface [12, 15, 16].

## Acknowledgments

This work was partially supported by the National Institute for Standards and Technology under grant #60NANB10D137 and by the STC program of the National Science Foundation Grant No. DMR-0120967. FAV would especially like to thank J. R. Knab for many illuminating discussions during the preparation of this manuscript.

## References and links

**1. **F. A. Vallejo and L. M. Hayden, “Design of ultra-broadband terahertz polymer waveguide emitters for telecom wavelengths using coupled mode theory,” Opt. Express **21**(5), 5842–5858 (2013). [CrossRef] [PubMed]

**2. **K. L. Vodopyanov, “Optical generation of narrow-band terahertz packets in periodically inverted electro-optic crystals: conversion efficiency and optimal laser pulse format,” Opt. Express **14**(6), 2263–2276 (2006). [CrossRef] [PubMed]

**3. **V. A. Kukushkin, “Efficient generation of terahertz pulses from single infrared beams in C/GaAs/C waveguiding heterostructures,” J. Opt. Soc. Am. B **23**(12), 2528–2534 (2006). [CrossRef]

**4. **Z. Ruan, G. Veronis, K. L. Vodopyanov, M. M. Fejer, and S. Fan, “Enhancement of optics-to-THz conversion efficiency by metallic slot waveguides,” Opt. Express **17**(16), 13502–13515 (2009). [CrossRef] [PubMed]

**5. **C. V. McLaughlin, L. M. Hayden, B. Polishak, S. Huang, J. Luo, T.-D. Kim, and A. K.-Y. Jen, “Wideband 15 THz response using organic electrooptic polymer emitter-sensor pairs at communications wavelengths,” Appl. Phys. Lett. **92**(15), 151107 (2008). [CrossRef]

**6. **J. B. Khurgin and G. Sun, “The case for using gap plasmon-polaritons in second-order optical nonlinear processes,” Opt. Express **20**(27), 28717–28723 (2012). [CrossRef] [PubMed]

**7. **A. Schneider, M. Neis, M. Stillhart, B. Ruiz, R. U. A. Khan, and P. Günter, “Generation of terahertz pulses through optical rectification in organic DAST crystals: theory and experiment,” J. Opt. Soc. Am. B **23**(9), 1822–1835 (2006). [CrossRef]

**8. **A. Nahata, A. S. Weling, and T. F. Heinz, “A wideband coherent terahertz spectroscopy system using optical rectification and electro-optic sampling,” Appl. Phys. Lett. **69**(16), 2321–2323 (1996). [CrossRef]

**9. **H. Cao, R. A. Linke, and A. Nahata, “Broadband generation of terahertz radiation in a waveguide,” Opt. Lett. **29**(15), 1751–1753 (2004). [CrossRef] [PubMed]

**10. **P. D. Cunningham, N. N. Valdes, F. A. Vallejo, L. M. Hayden, B. Polishak, X.-H. Zhou, J. Luo, A. K.-Y. Jen, J. C. Williams, and R. J. Twieg, “Broadband terahertz characterization of the refractive index and absorption of some important polymeric and organic electro-optic materials,” J. Appl. Phys. **109**(4), 043505 (2011). [CrossRef]

**11. **K. D. Singer, M. G. Kuzyk, and J. E. Sohn, “Second-order nonlinear-optical processes in orientationally ordered materials: relationship between molecular and macroscopic properties,” J. Opt. Soc. Am. B **4**(6), 968–976 (1987). [CrossRef]

**12. **Ş. E. Kocabaş, G. Veronis, D. A. B. Miller, and S. Fan, “Modal analysis and coupling in metal-insulator-metal waveguides,” Phys. Rev. B **79**(3), 035120 (2009). [CrossRef]

**13. **H. Ma, A. K. Y. Jen, and L. R. Dalton, “Polymer-based optical waveguides: materials, processing, and devices,” Adv. Mater. **14**(19), 1339–1365 (2002). [CrossRef]

**14. **Y.-F. Li and J. W. Y. Lit, “General formulas for the guiding properties of a multilayer slab waveguide,” J. Opt. Soc. Am. A **4**(4), 671–677 (1987). [CrossRef]

**15. **R. E. Smith, G. W. Forbes, and S. N. Houde-Walter, “Unfolding the multivalued planar waveguide dispersion relation,” IEEE J. Quantum Electron. **29**(4), 1031–1034 (1993). [CrossRef]

**16. **X. Ying and I. Katz, “A simple reliable solver for all the roots of a nonlinear function in a given domain,” Computing **41**(4), 317–333 (1989). [CrossRef]