Design of ultra-broadband terahertz polymer waveguide emitters for telecom wavelengths using coupled mode theory

We use coupled mode theory, adequately incorporating optical losses, to model ultra-broadband terahertz (THz) waveguide emitters (0.120 THz) based on difference frequency generation of femtosecond infrared (IR) optical pulses. We apply the model to a generic, symmetric, five-layer, metal/cladding/core waveguide structure using transfer matrix theory. We provide a design strategy for an efficient ultra-broadband THz emitter and apply it to polymer waveguides with a nonlinear core composed of a poled guest-host electro-optic polymer composite and pumped by a pulsed fiber laser system operating at 1567 nm. The predicted bandwidths are greater than 15 THz and we find a high conversion efficiency of 1.2×10 W by balancing both the modal phase-matching and effective mode attenuation. ©2013 Optical Society of America OCIS codes: (190.4390) Nonlinear optics, integrated optics; (320.7110) Ultrafast nonlinear optics; (300.6495) Spectroscopy, terahertz. References and links 1. X. Zheng, C. V. McLaughlin, P. D. Cunningham, and L. M. Hayden, "Organic broadband terahertz sources and sensors," J. Nanoelectron. Optoelectron. 2, 58-76 (2007). 2. K. Liu, J. Xu, and X.-C. Zhang, "GaSe crystals for broadband terahertz wave detection," Appl. Phys. Lett. 85, 863-865 (2004). 3. U. Peschel, K. Bubke, D. C. Hutchings, J. S. Aitchison, and J. M. Arnold, "Optical rectification in a traveling-wave geometry," Phys. Rev. A 60, 4918-4926 (1999). 4. Y. J. Ding, "Terahertz parametric converters by use of novel metallic-dielectric hybrid waveguides," J. Opt. Soc. Am. B 23, 1354-1359 (2006). 5. A. Marandi, T. E. Darcie, and P. P. M. So, "Design of a continuous-wave tunable terahertz source using waveguide-phase-matched GaAs," Opt. Express 16, 10427-10433 (2008). 6. M. Cherchi, A. Taormina, A. C. Busacca, R. L. Oliveri, S. Bivona, A. C. Cino, S. Stivala, S. R. Sanseverino, and C. Leone, "Exploiting the optical quadratic nonlinearity of zinc-blende semiconductors for guided-wave terahertz generation: A material comparison," IEEE J. Quantum Electron. 46, 368-376 (2010). 7. C. Staus, T. Kuech, and L. McCaughan, "Continuously phase-matched terahertz difference frequency generation in an embedded-waveguide structure supporting only fundamental modes," Opt. Express 16, 13296-13303 (2008). 8. 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, 13502-13515 (2009). 9. C. Tao, S. Junqiang, L. Linsen, and T. Jianguan, "Proposal for efficient terahertz-wave difference frequency generation in an AlGaAs photonic crystal waveguide," J. Lightw. Technol. 30, 2156-2162 (2012). 10. Y. Li, X. Hu, F. Liu, J. Li, Q. Xing, M. Hu, C. Lu, and C. Wang, "Terahertz waveguide emitters in photonic crystal fiber form," J. Opt. Soc. Am. B 29, 3114-3118 (2012). 11. Z. Wang, H. Liu, N. Huang, Q. Sun, and J. Wen, "Efficient terahertz-wave generation via four-wave mixing in silicon membrane waveguides," Opt. Express 20, 8920-8928 (2012). 12. F. F. Lu, T. Li, J. Xu, Z. D. Xie, L. Li, S. N. Zhu, and Y. Y. Zhu, "Surface plasmon polariton enhanced by optical parametric amplification in nonlinear hybrid waveguide," Opt. Express 19, 2858-2865 (2011). 13. S. B. Hasan, C. Rockstuhl, T. Pertsch, and F. Lederer, "Second-order nonlinear frequency conversion processes in plasmonic slot waveguides," J. Opt. Soc. Am. B 29, 1606-1611 (2012). 14. F. M. Pigozzo, D. Modotto, and S. Wabnitz, "Second harmonic generation by modal phase matching involving optical and plasmonic modes," Opt. Lett. 37, 2244-2246 (2012). 15. S. B. Bodrov, I. E. Ilyakov, B. V. Shishkin, and A. N. Stepanov, "Efficient terahertz generation by optical rectification in Si-LiNbO3-air-metal sandwich structure with variable air gap," Appl. Phys. Lett. 100, 201114 (2012). 16. M. I. Bakunov, S. B. Bodrov, A. V. Maslov, and M. Hangyo, "Theory of terahertz generation in a slab of electro-optic material using an ultrashort laser pulse focused to a line," Phys. Rev. B 76, 085346 (2007). 17. Y.-C. Huang, T.-D. Wang, Y.-H. Lin, C.-H. Lee, M.-Y. Chuang, Y.-Y. Lin, and F.-Y. Lin, "Forward and backward THz-wave difference frequency generations from a rectangular nonlinear waveguide," Opt. Express 19, 24577-24582 (2011). 18. V. Berger and C. Sirtori, "Nonlinear phase matching in THz semiconductor waveguides," Semicond. Sci. Technol. 19, 964-970 (2004). 19. H. Cao, R. A. Linke, and A. Nahata, "Broadband generation of terahertz radiation in a waveguide," Opt. Lett. 29, 1751-1753 (2004). 20. G. Chang, C. J. Divin, J. Yang, M. A. Musheinish, S. L. Williamson, A. Galvanauskas, and T. B. Norris, "GaP waveguide emitters for high power broadband THz generation pumped by Yb-doped fiber lasers," Opt. Express 15, 16308-16315 (2007). 21. F. Peter, S. Winnerl, H. Schneider, and M. Helm, "Excitation wavelength dependence of phase matched terahertz emission from a GaAs slab," Opt. Express 18, 19574-19580 (2010). 22. K. Saito, T. Tanabe, Y. Oyama, K. Suto, and J.-i. Nishizawa, "Terahertz-wave generation by GaP rib waveguides via collinear phase-matched difference-frequency mixing of near-infrared lasers," J. Appl. Phys. 105, 063102 (2009). 23. V. A. Kukushkin, "Generation of terahertz pulses from tightly focused single near-infrared pulses in double-plasmon waveguides," J. Opt. Soc. Am. B 25, 818-824 (2008). 24. K. Saito, T. Tanabe, and Y. Oyama, "Elliptically polarized THz-wave generation from GaP-THz planar waveguide via collinear phase-matched difference frequency mixing," Opt. Express 20, 26082-26088 (2012). 25. P. Bienstman, "Rigorous and efficient modeling of wavelength scale photonic components," (Universiteit Gent 2000-2001). 26. 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, 671-677 (1987). 27. X. Ying and I. Katz, "A simple reliable solver for all the roots of a nonlinear function in a given domain," Computing 41, 317-333 (1989). 28. R. E. Smith, G. W. Forbes, and S. N. Houde-Walter, "Unfolding the multivalued planar waveguide dispersion relation," IEEE J. Quantum Electron. 29, 1031-1034 (1993). 29. 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, 968-976 (1987). 30. A. Nahata, J. Shan, J. T. Yardley, and C. Wu, "Electro-optic determination of the nonlinear-optical properties of a covalently functionalized Disperse Red 1 copolymer," J. Opt. Soc. Am. B 10, 1553-1564 (1993). 31. Ş. E. Kocabaş, G. Veronis, D. A. B. Miller, and S. Fan, "Modal analysis and coupling in metal-insulatormetal waveguides," Phys. Rev. B 79, 035120 (2009). 32. L. Mandel and E. Wolf, Optical Coherence and Quantum Optics (Cambridge University, 1995). 33. 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, 043505 (2011). 34. E. D. Palik, ed., Handbook of Optical Constants of Solids (Academic, 1998) Vol. I III. 35. M. A. Ordal, R. J. Bell, J. R. W. Alexander, L. A. Newquist, and M. R. Querry, "Optical properties of Al, Fe, Ti, Ta, W, and Mo at submillimeter wavelengths," Appl. Opt. 27, 1203-1209 (1988). 36. S. N. Kasarova, N. G. Sultanova, C. D. Ivanov, and I. D. Nikolov, "Analysis of the dispersion of optical plastic materials," Opt. Mater. 29, 1481-1490 (2007). 37. M. Dellnitz, O. Schütze, and Q. Zheng, "Locating all the zeros of an analytic function in one complex variable," J. Comput. Appl. Math. 138, 325-333 (2002). 38. P. N. Robson and P. C. Kendall, Rib Waveguide Theory by the Spectral Index Method, Electronic & electrical engineering research studies: Optoelectronics series (Research Studies, 1990).


Introduction
A widely used method to produce ultra-broadband terahertz (THz) pulses is to mix the frequency components of infrared (IR) femtosecond (fs) pulses through difference frequency generation (DFG) in organic and inorganic crystals [1,2].Pulsed fiber laser systems operating at telecom wavelengths are compact and have the potential to provide portable ultrabroadband THz sources outside the laboratory.However, these systems are power limited and only a few materials are well phase-matched at these wavelengths for generation in a broadband THz range (DAST, electro-optic (EO) polymers [1]).Poled EO guest-host polymer films are attractive since their material properties are tunable and provide an inexpensive alternative to crystal based ultra-broadband THz emitters.Unfortunately, these EO composites are somewhat absorptive in the THz range.In addition, to achieve maximum conversion when using them as bulk emitters the IR radiation must be polarized parallel to the poling direction, normal to the film surface.This slanted geometry results in a non-optimal overlap of the input IR polarization with the emitter nonlinearity.
In this work, we report a model for ultra-broadband THz DFG with a frequency domain (FD) coupled mode theory (CMT), similar to [23].The model is specially equipped to describe broadband applications since it is formulated in the FD and our expressions account for the contributions of all the possible frequency pairs in the pump spectrum that drive the DFG of single THz component.A novelty of our approach, different from [23], is that our CMT formulation is equipped to deal with lossy modes [8,13,25] and we account for higher order dispersion effects.Loss effects cannot be neglected in THz applications particularly if they include metallic structures.Group velocity dispersion effects must also be accounted for in structures with nonlinear cores whose dispersion curves present considerable curvature at the pumping wavelength; pump pulses coupled to the fundamental IR mode would broaden affecting the generated THz power spectra.Also, the CMT expressions that we derived are valid for an interaction between an arbitrary number of IR modes and THz modes.In previous works, references [3,23], the generation of THz pulses by self-frequency mixing of an IR fs pulse was described by single IR-to-single THz-mode interactions.Furthermore, the time-domain formulations [3,9], are not adequate to describe broadband THz generation since they assume that the THz mode profiles remain constant across the THz spectrum, and higher order dispersion effects are usually neglected.We specifically give expressions for symmetric, five-layer, slab waveguides whose nonlinear core belongs to a !mm point symmetry class.We apply the transfer matrix theory and mode matching methods to solve for the complex propagation constants [26][27][28].
A symmetric, five-layer, slab waveguide, Fig. 1(a), similar to the one reported in [19], has the ideal geometry for exploiting the EO polymer nonlinearity, achieving good phasematching conditions and managing THz attenuation.The two metal capping layers help confine the modes involved in the interaction.The IR pump is coupled to the fundamental TM mode that is bounded inside the EO core and entirely polarized parallel to the poling direction.The THz radiation is collected by the 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 metal capping layers.
Using a design strategy where the conversion between a single-IR-mode and a single-THz -mode is optimized, we demonstrate an efficient broadband THz emitter (>15 THz) composed of a symmetric five-layer slab waveguide with the core composed of the guest-host EO composite DAPC (in the appendix we demonstrate how our theory allows for modeling the conversion of two IR pump modes into a single THz mode and use that result to verify the output from a previous experimental demonstration by Cao et. al.,[19]).In DAPC, guest EO chromophores DCDHF-6-V (inset Fig. 1(b)) are embedded in an amorphous polycarbonate (APC) matrix host at 40% mass ratio.For our computation we assumed a pulsed fiber laser pump operating at 1567 nm with a pulse width of 47-fs, 100 MHz repetition rate and with a bandwidth extending from (1504 to 1630) nm.Our results show wide bandwidths greater than 15 THz and a maximum conversion efficiency of !THz = P THz / P pump 2 !1.2 " 10 #4 W #1 .This efficiency is higher than previously proposed narrow band DFG THz waveguide emitters with !" 8 # 10 $5 W $1 [5, 8, 9], and those broadband THz waveguide emitters previously demonstrated, !< 7 " 10 #7 W #1 [19-22, 24].Modifying the waveguide dimensions allows tuning both modal effective indices to achieve better phase-matching and lower loss conditions for this single mode to single mode DFG interaction.

Frequency domain coupled mode theory
Consider a nonlinear medium interacting with electromagnetic radiation and denote !" E( " r,! ) the Fourier transform of the electric field !E( !r,t) .Using the CMT approach the solution to the perturbed problem !" E can be expanded as Here êl (r !," )e i# l (" )z are the modes of the waveguide structure when the nonlinearity is not present, i.e. !" P NL = 0 .The ! l 's are the mode constants, the êl 's are the mode profile functions, z is the longitudinal direction of propagation and r ! is the transverse position component.Using the Lorentz reciprocity theorem together with a mirror symmetry operation it can be proved that the amplitude functions A l (z,! ) satisfy the following relation [8,25] !A l (z," ) Here P l (! ) = 2 d S " r # ẑ $( êl,# % ĥl,# ) , êl,! and êl,z are the transverse and longitudinal components of the mode profile, and S is the cross-sectional area of the nonlinear section in the waveguide.Different from other coupled mode theories it is worth noting that this version holds even in the case of lossy materials where both the mode constants !l = k 0 N eff ,l + i" l / 2 and the mode profiles are complex.The time-domain electric field is given by  If we consider only second order processes and assume that the nonlinearity acts instantaneously we find from Eqs. ( 1) and (3) that the nonlinear polarization is given by, P !" NL (t) = !0 ! (2): The generation of THz pulses through difference frequency mixing is due to the second term in the last equation, !P DFG (2) .The first term is responsible for SHG and we will neglect it here.
To extract the !P DFG (2) we perform the following change in variable, (!, " ! ) # ($ % !& " !,! ) , the Jacobian of this transformation is 1 and !P DFG (2) If there are N modes involved in the interaction, N coupled first order partial differential equations can be derived describing the interaction.From Eqs. ( 2) and ( 5), we find that the THz mode amplitudes are given by, Where the phase-mismatch is given by !"(#,$) Equation ( 7) measures the coupling strength since it is proportional to χ (2) and to the overlap integral among the modes involved in the interaction.Note that Eq. ( 6) accounts for all the contributions of all the possible frequency pairs (! " #,! ) in the pump spectrum that drive the DFG of single THz component, even when multiple IR modes are involved, which is a noted difference to our approach compared to previous descriptions in the literature.Up to now the theory is completely general.However, Eq. ( 7) varies depending on the nonlinear material used in the DFG interaction and the specific geometry of the waveguide.

TM to TM conversion for !mm point symmetry class nonlinearities
For the symmetric five-layer slab waveguide (Fig. 1(a)) we choose the nonlinear core to be made of a poled EO guest-host polymer composite.This material falls into the !mm point symmetry class.We take the x-axis to be the polar axis, which coincides with the polinginduced macroscopic dipole moment (arrows in Fig. 1(a)).The nonzero tensor components under this convention are !xxx (2) , ! yyx (2) = !zzx (2) , ! yxy (2) = !zxz (2) and !xyy (2) = !xzz (2) [29].To access the maximum nonlinearity we will assume that the input IR beam is injected with its polarization parallel to the poling direction, i.e. the IR beam is x polarized.Since the birefringence effects for EO composites are generally small for IR wavelengths we can assume IR radiation will couple to a TM mode and we set E y = 0 , then ! P DFG (!) has only x and z components (in the plane of incidence), hence the DFG interaction favors a TM to TM conversion.Typically for EO materials what is measured is the value of the EO tensor r ij,k ; the !ijk (2) tensor is related to r ij,k by !ijk (2) (";" ,0) = #n i (" ) 2 n j (" ) 2 r ij,k (" ) / 2 [29, 30].However, these ! (2)values correspond to the EO process where != !+ 0 .We must compute the values of ! (2)that correspond to difference frequency generation ! (2)(";#," $ # ) ; here THz frequencies are denoted by ! and the IR frequencies are denoted by ! .Since !<< " we can set ! " 0 , thus we need only to determine the tensor for the optical rectification (OR) effect ! (2) (0;",#" ) .Imposing Kleinman symmetry the values for OR can be obtained from the ones from the EO process by !ijk (2) (0;" ,#" ) = !jki (2) (";" ,0) * .For well poled guest-host EO composites with chromophores with high electric dipole moments !xxx (2) " 3! xyy (2) " 3! yxy (2) [29].Therefore, it is only necessary to know the value of r 33 (! ) to compute all the values of !ijk (2) needed in Eqs. ( 7) and ( 6).

Input beam conditions, completeness and energetic considerations
We assume a Gaussian input beam, both in time and in the x coordinate, injected at z = 0 and uniformly distributed across the y direction along a strip of length Here !0 is a measure of the pulse duration, !0 is the central frequency and w 0 is the beam waist of the IR pump laser in the x direction at the input face of the waveguide.The value of E 0 depends on the average power of the pump laser P pump and the repetition rate 1/T rep of the train of pulses produced in the laser system.To extract E 0 we compute the energy per pulse Here L y is twice the beam waist in the y direction and we assume that T rep >> !0 and !o " 0 >> 1 .Under these conditions E 0 is given entirely from experimental quantities, and !E in (! ) is given by the Fourier transform of Eq. ( 9) Note the extra factor of 1/2 that arises from the convention established through Eq. ( 4).We proceed to find a valid expansion for !E in (! , x) in terms of the waveguide modes.Care must be taken here since when considering materials with losses the differential operator that gives the mode decomposition is non-Hermitian (non-self-adjoint) and its spectrum is complex, see Appendix A Eqs. ( 16) and (18).The general classification of the spectrum for non-selfadjoint operators is still an open problem and completeness is difficult to prove [31].Also, the standard inner product < i, j >= !" (x)E i (x)E j (x) * dx = # ij , with !(x) = n(x) 2 for TM and !(x) = 1 for TE modes, is not necessarily valid for lossy materials.In fact, for TM modes when !(x) "! the solution space for the differential operator will lack a well-defined metric < i,i >!! and the solution space is not a Hilbert space.Nevertheless, for metal-insulatormetal waveguides it has been shown that its solution space forms a complete basis even when the modes are complex and result from a non-self-adjoint operator [31].Fortunately, using the Lorentz reciprocity theorem combined with symmetry operations yields [8,25] Here S is the waveguide cross section that lies on the xy plane.Eq. ( 12) holds independently if the modes are lossy or not.This orthogonality relation provides means to normalize by dx as a pseudo inner product.It follows that: Contrary to usual conventions |A l | 2 is not the power in Watts transported by the l-th mode.The total power needs to be calculated using Eq. ( 14).

Energetic considerations
The ensemble formed by the train of IR pump pulses and the corresponding THz pulses generated through DFG form a stationary random process.The average THz power is ( ) Here, the sum runs only on the THz modes.To compute the integral in !' we used sin(x / !) / (! x) " # (x) as ! !0 where !(x) is the Dirac delta function, note 2 / T rep !MHz and (! ± !) " THz .The first term in Eq. ( 14) cancels while computing the integral since E(!) = 0 for !< 0 .Note that the quantity dP THz (!) has units of Watt/Hz and it represents the power spectral density of the generated THz radiation.The fact that we can indeed define this "power spectrum" is not surprising due to the stochastic nature of the THz generation process.This is explicitly stated in the Wiener-Khintchine theorem: the autocorrelation function of a stationary random process and the spectral density form a Fourier transform pair, i.e. !(" ) = S # ($ )e % i$" [32].In this realm, the average THz power takes the role of the same-time autocorrelation function, !(" = 0) = P THz , and S(! ) = dP THz .For this TM to TM conversion the output THz power density generated by the structure in Fig. 1(a) is

Results and discussion
To benchmark the model we applied the theory to the waveguide device fabricated by Cao et.al., in 2004 that served as a proof of concept for broadband THz generation in waveguides [19].Our predictions for this structure, reported in detail in Appendix B, are in excellent agreement with the experimental results of the peak-to-peak ratio of the waveforms produced by emitters with different propagation lengths, the only experimental result reported in that work.We consider this to be a satisfactory benchmark since we had no knowledge of the values of extinction coefficients for the cladding and core in that experiment nor did we have access to the exact detector response function of their photoconductive antenna.
In this section, we proceed to study structures composed of DAPC, an inexpensive polymeric EO guest-host composite with a high EO coefficient at IR frequencies.The DAPC films we produce through contact poling typically have EO coefficients r 33 !40 pm / V at 830 nm [1].
The corresponding susceptibility value !xxx (2) " 53pm / V at 1567 nm is obtained from the measured values at 830 nm assuming a two level model [29], since DAPC is known to have a single resonance at 610 nm (see Fig. 1(b)).The choice of cladding materials can lead to better performance based on their optical properties at both THz and IR frequency ranges.In particular, we report here results for two polymeric cladding materials that were shown [33] to have low THz losses in a broadband range, (0.1-10 THz), polyethylene cyclic olefin copolymer (TOPAS ® ) and polystyrene (PS).To calculate both the THz and the IR modes, we must know both the refractive index and absorption coefficients of the waveguide constituents at both frequency ranges.We choose Al for the metal capping layers since its properties are well characterized both at the THz and IR ranges.We used a global Drude fit to the experimental values measured on thin films of vapor deposited Al in the THz range obtained from [34,35], and we used interpolating fits in the IR range for experimental values from [34].In Fig. 2, we show fits and experimental values.The refractive index and attenuation at THz frequencies for DAPC, TOPAS ® and PS were measured with a broadband THz spectrometer in the (0.1-10) THz range [33].PS and TOPAS ® showed flat indices, varying only 0.005 for TOPAS ® and 0.01 for PS in the (0.1-10) THz range.Over this range TOPAS ® loss is less than 3 cm -1 and PS is less than 15 cm -1 .We fit these data using a permittivity model with constant index and extinction coefficient.DAPC on the other hand showed several absorptions features in the THz range and was fit to a five-oscillator-Lorentz model.In Figs.3(c) and 3(d) we show the fits used in these models for the THz.For the IR wavelengths, we measured indices for TOPAS ® and DAPC using a prism coupler device Metricon model 2010/M at 1064 nm, 1311 nm and 1537 nm, and the index values over the (430-1052) nm range for PS were taken from [36].Using these values, Sellmeier fits, n 2 = a + b! 2 / (! 2 " !0 2 ) , were computed for DAPC, TOPAS ® and PS, see Fig.

3(a)
. The fit parameters a,b,! 0 ( ) are respectively (2.35, 0.28, 610 nm), (2.2,0.093, 358nm) and (2.26, 0.19, 302 nm).The losses in the IR are small for DAPC, TOPAS ® and PS and were neglected.In Fig. 3 we show the calculated modal dispersion, β m (ω) = k 0 N eff,m (ω) + iα m (ω)/2, for even modes at both the THz frequencies (Ω) and the IR wavelengths (λ) for DAPC structures with either PS or TOPAS ® as the cladding material.Since we assumed symmetric input conditions, only even modes are involved in the DFG interaction.Note that for a waveguide of sub-wavelength dimensions in which the wavelength of the guided radiation is greater than the waveguide size, such as the one studied here at the THz frequency range (Ω), there is a single symmetric TM guided mode accepted by the waveguide.Figures 3(a and IR wavelength range !" (1.3-1.9)µm.We choose a broad range for the IR frequencies to correctly calculate the integral appearing in Eq. ( 6).The main challenge we faced in our implementation results from including losses in the model.The search for the solutions to the dispersion equation must be carried out in the complex plane since β m is complex.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.The details of this procedure are out of the scope of this work and are explained elsewhere [26-28, 31, 37].
As part of our design strategy we wanted our structures to be single-core-mode across the entire IR pump beam bandwidth centered at 1567 nm, delimited by the vertical dashed lines in Fig. 3(a).Coupling the input beam into a single mode inside the nonlinear core guarantees that most of the pump power contributes coherently towards the generation process.This determines the value of the core half thickness d.In Figs.3(a) and 3(b) we show the modal dispersion at IR wavelengths for a DAPC structure with d = 1.579 µm, t = 20 µm and PS cladding.The fundamental mode (m = 0) is a core mode since it is well confined in the core and its effective index lies between the core and cladding indices.The (m = 1) mode is plasmonic in nature and will not be excited in focused scenarios since its mode profile is highly confined to the cladding layers at the interface with the metal capping layers.Note in Fig. 3(b) that the attenuation for this plasmon mode is an order of magnitude larger than for the other modes.Higher order modes (m = 2, 3, 4, ...) are cladding modes since their effective indices are lower than the cladding index across the pump beam bandwidth.For thicker cores the (m = 2) mode will become a core mode and this establishes an upper limit for the core half thickness, d = 1.274 µm for structures with TOPAS ® cladding and d = 1.579 µm with PS cladding.
To have single mode structures in the (0.1-20) THz range, the half thickness of the waveguide t is also limited to half of the minimum THz wavelength, t < 15 µm.For THz frequencies the fundamental (m = 0) TM THz mode is a TEM-like mode since the mode profile for the THz frequencies is uniform in the core and cladding sections and the effective index dispersion resembles the dispersion of the core.In Figs.3(c) and 3(d), we show the calculated modal dispersion, for the fundamental (m = 0) TEM-like mode in the THz range for the most efficient structures we found for t ![d, 15 µm] with both TOPAS ® and PS claddings.The phase-matching concept for THz generation in bulk geometries can be extended to waveguide emitters [1,18,5].Better phase-matching (i.e., more efficient mode conversion from IR to THz) is obtained when the effective group index of the IR mode, n g,IR = N eff ,0 − λ(dN eff ,0 / dλ) , matches the effective index of the THz mode N eff,0 .The structures with TOPAS ® cladding are less efficient than the PS structures due to a higher effective index mismatch caused by the lower THz index in TOPAS ® , see Fig. 4 (a).Varying t has a very significant effect on the effective index of the TEM-like THz mode see Fig. 3(c).Whereas for the IR range, the modal dispersion is relatively insensitive to variations in t as long as the cladding thickness is comparable to the pump's center wavelength, t -d ≥ λ 0 .This provides a useful tool to design a THz waveguide emitter achieving better phase matching by changing t.For thinner claddings, N eff,m for the TEM-like mode increases approaching the index of the core since most of the radiation will be guided in the higher index core and that is desired for phase-matching purposes, see Fig. 4(a).Also, thinner claddings result in better mode overlap, which increases the coupling Eqs. ( 6) and (7).
Unfortunately, the mode attenuation increases too since the core is also more absorptive, Fig. 3(d).Surprisingly, since the conversion efficiency depends on both the effective index matching and on the mode attenuation of the THz modes, the most efficient structure is not necessarily the one that is better phase-matched at 1567 nm.THz losses are extremely important and in an optimized structure a perfect balance between both modal phase-matching and mode losses effects must be found.This is seen in Fig. 4(b) where the nonlinear efficiency, !THz = P THz / P IR 2 , is shown for the most optimal structures we found.For example, even though the PS structure with t = 2.3 µm has a lower index mismatch (Fig. 4(a)), structures with thicker PS claddings (t = 2.5 µm and 2.7 µm) and with lower THz mode attenuation Fig. 3(d)) have higher conversion efficiencies (Fig. 4(b)).
In our computations, we obtain a closed solution to Eq. ( 6) by decoupling the amplitudes for the THz modes by: (i) neglecting pump beam depletion of the IR modes (set A l IR (z) = A l IR (0) ), (ii) assuming that the THz modes are zero valued at the input face of the waveguide A m THz (0) = 0 and (iii) performing the straightforward integral in z .We assumed a Gaussian pulse with 47-fs pulse width centered at 1567nm, 10 mW of input pump power, and repetition rate of 100 MHz.In addition, we assumed symmetric input conditions and, an input Gaussian beam focused into the core with the aid of an elliptical lens to a beam waist (in the x-axis) w 0 = 1.4 µm for the TOPAS structure and 1.8 µm for those with PS.To couple most of the input pump power to the fundamental IR mode, the beam waist of the input profile is chosen such that the expansion in Eq. ( 1) is dominated by the l = 0 mode.For these values of w 0 about 81% of the pump power is coupled into the fundamental IR mode.To determine the power coupled to a certain IR mode we use Eq. ( 15) limited to the l-th IR mode using Eq. ( 13) to compute A l IR ; integrating across the pump spectral range yields the total power coupled into the l-th mode.In our calculations we did not include coupling losses resulting from reflections at the input face of the waveguide or at the cylindrical coupling lens.Therefore, we assume 10 mW of power crossing the input interface of the waveguide emitter just after the z ≥ 0, and only around 8 mW couple into the fundamental IR mode.For particular values of t and d, there exists an optimal emitter length, L opt , for which the nonlinear efficiency η THz is a maximum.In Fig. 5(a) we plot the spectral power density for our most optimal structure with d = 1.597 µm, t = 2.5 µm and PS cladding, for different emitter lengths L. The maximum η THz = 1.2 × 10 -4 W -1 for this structure occurs when L = 550 µm.In Fig. 5(b) we show the pulse produced by such a structure with corresponding bandwidth >15 THz.For longer structures the higher attenuation observed for higher frequencies degrades the bandwidth of the device.This attenuation effect is less significant for the structure with t = 2.7 µm which has a lower mode attenuation but due to the greater phasemismatch results in a slightly lower maximum efficiency than the t = 2.5 µm structure.Unfortunately efficiency values were not reported in [19].In our benchmark (see Appendix B) we predicted nonlinear conversion efficiencies ~ 10 -6 W -1 for the structure in [19].These values are 100 times smaller than the efficiency of our optimized structure.An enhancement factor 25 results just from using an EO polymer with five times higher χ (2)  values, P TH ~ (χ (2) ) 2 .Input conditions are also crucial since this effect is quadratic in η THz ; a proper choice of w 0 allows coupling a higher fraction of the pump power into the fundamental IR mode.Another source of enhancement is the access to power at higher THz frequencies due mainly to our shorter pump pulse duration that provides broader IR bandwidth than the one used in [19].Structures with optimal emitter lengths, L opt , can effectively collect the power at all accessible THz frequencies.The structures in [19] had longer propagation lengths (a few millimeters) and the mode attenuation we calculate for these structures are bigger than 50 cm -1 past 3 THz.The THz power generated above 3 THz by the structures in [19] is considerably reduced for longer emitter lengths, similarly as in Fig. 5(a).Finally, our optimized structure also has a slightly higher overlap coefficient Eq. ( 7) than those in [19] due to higher mode confinement.
We proceed to study the effect of group velocity dispersion (GVD) in our results.A pulse coupled to the fundamental IR mode of one of our structures will be confined within the DAPC core.Due to the low curvature observed in the dispersion curve of the DAPC index (quasi-linear) in Fig 3(a) the GVD effects should be small; in addition the emitter lengths are short (550 µm).We confirmed this by using an approximate solution to Eq. ( 6) where we compute analytically the integral in frequency but we ignore GVD and higher order dispersion effects.This simplified 3-wave approach produces almost identical results as those observed in Fig. 5(a) and Fig. 4(b) for such short interaction lengths.
As a final remark, we point out that from our results we expect that laser systems with longer central wavelengths would yield better performances for the DAPC waveguide emitters we propose here.Although DAPC χ (2) values are slightly lower for longer wavelengths (see Fig. 1(b)), the phase-mismatch for these materials is greater for shorter wavelengths (Fig. 4(a)).In fact, above 1.7 µm there exists a wavelength for which each PS structure is completely phase-matched.However, since current commercial fiber laser systems have central wavelengths near 1567 nm, the structures we propose here are more suited for stand alone broadband waveguide THz emitters.An optimal set of materials is one for which the crossing between n g,IR and N eff,0 shifts to lower wavelengths (or closer to the resonance) and occurs for thicker claddings (which lowers the modal losses).

Conclusions
We presented a semi-analytical model that is able to predict the complete bandwidth generated by a broadband THz waveguide emitter based on the difference frequency mixing of a pulsed fs IR pump laser.This frequency domain CMT model is suited to accurately incorporate both THz and IR losses.For the ∞mm point class type of second order nonlinearity, we applied the transfer matrix theory and mode-matching methods to model the symmetric five-layered metal/cladding/nonlinear-core waveguide.We applied a design strategy based in single IR mode to single THz mode conversion, to design an efficient five-layered waveguide structure made out of a DAPC core and PS cladding layers.The most efficient THz emitter we found, η THz = 1.2 × 10 -4 W -1 , yields an efficiency comparable to state of the art CW THz waveguide emitters but the wide bandwidth ~ 15 THz constitutes a record high for broadband THz waveguide emitters.We found that for an optimized structure there must be a perfect balance between both modal phase-matching and mode losses effects.This can be achieved by careful selection of the cladding, cladding thickness and emitter length.The theoretical method presented here is entirely general and can be applied to any material system provided the appropriate material refractive index and absorption dispersion information is known.

Appendix A: Waveguide modes obtained through transfer matrix theory
We assume that the fields inside the structure shown in Fig. 1(a) are propagating in the z direction with a mode constant . Since metal layers are present in the structure, the TM mode profiles will present an abrupt discontinuity at the interface between metal and cladding layers due to the high index contrast, see below Eq. ( 18).Mode solvers based on finite difference (FD) methods have to be carefully modified and implemented to account for this effect.Luckily for layered structures !/ !y = 0 and the mode profiles depend only on the x component, !E(x, y) = !E(x) .In the absence of any nonlinearity an analytical expression for the modes is: with k 0 = ω/c for each i-th layer in the structure (i = 0,...,4).Note that for the symmetric fivelayer waveguide n 0 = n 4 and n 1 = n 3 .It must be pointed out that Eq. ( 16) yields the modes only when combined with appropriate boundary conditions since all of Maxwell's equations must be simultaneously fulfilled.In the linear regime and for this type of waveguide, the evolution of any input radiation propagating longitudinally can be expanded as a superposition of two independent sets of solutions to Eq. ( 16), the TE {H x , E y , H z } and TM {E x , H y , E z } modes.Using an electric field formulation in which we solve Eq. ( 16) only for the principal electric component, E = E y for TE modes or E = E x for TM modes, the remaining components of the mode for both TM and TE are given in terms of the principal field component E [38].For example, TM modes are given by H y = n i 2 E / (n eff Z 0 ) and , and E is given by at the i-th layer, here k i Here !i+1 = 1 for TE modes and !i+1 = n i+1 2 / n i 2 for TM modes.This results in a linear system of equations that relate A i , B i { } to A i+1 , B i+1 { } through matrix relations.To guarantee that the radiation in the mode is confined in the waveguide and decays exponentially in the outer layers we impose the guiding mode condition A 0 = B 4 = 0 [26].By sequential matrix multiplication the coefficients of both outer layers are related by:

The dispersion relation
Bound modes must satisfy m 1,1 = 0 , this requirement is known as the dispersion relation for both TE and TM modes.Its solutions yield the allowed mode constants {! n } 's.Bound modes can be further classified m 2,1 = 1 for even and m 2,1 = !1 for odd modes.The dispersion relation for the TM modes is explicitly [26]: ). (20)   Here The upper term in the brackets applies to even modes and the lower term to odd modes.

Mode profiles
The A i , B i { } coefficients can be solved in terms of B 0 = m 2,1 A 4 .Setting x = 0 at the center of the waveguide we find for the even modes, Here B 0 is a normalization constant and is related to the power carried by the mode.

Appendix B.
To benchmark the model, we applied the theory to the waveguide device fabricated by Cao, Linke, and Nahata in 2004 that served as a proof of concept for broadband THz generation in waveguides [19].The THz output generated from this device was well documented in their paper and most of the experimental parameters necessary to apply our theory were available.The THz emitter used in this demonstration consisted of a symmetric five-layer waveguide as the one in Fig. 1(a) on top of a silicon substrate.The outer metal layers consisted of 200 nm layers of Al but we assumed that the metal capping layers were infinitely thick since the absorption depth in the metal is considerably thinner than the thickness of the metal layers.The nonlinear core was made out of a 3 µm layer (d = 1.5 µm) of the copolymer Disperse Red 1 covalently functionalized at a 10%-mol ratio to a methyl-methacrylate backbone (DR1-MMA 10%).The structure was corona poled to achieve a second order nonlinearity.We used r 33 = 3.6 pm/V which was the value measured from corona poled films of DR1-MMA 16.5% measured at 810 nm [30], this corresponds to r 33 = 2.21 pm/V at 820 nm for DR1-MMA 10%.We obtain this value scaling by 10/16 and using a two level system model on the susceptibility with a resonance at 470 nm [30].With the 3.5 µm thick cladding layers of UVcurable acrylate (Norland optical adhesive NOA 81) the cladding/core/cladding section was 10 µm thick, (t = 5 µm).The waveguide was pumped with a mode-locked Ti:Sapphire laser oscillator with a pulse width of !FWHM = 100 fs operating at 820 nm, !0 = 365.6THz .The other pump beam parameters used in the model were P pump = 3 mW and L y = 6 mm which correspond to the reported values.We assumed that the input pump beam was focused to a beam waist w 0 = 1.5 µm in the x axis at the input interface with a cylindrical lens used to couple the pump beam into the waveguide.We also assumed that the repetition rate of the laser system was 100 MHz.The THz indexes for both DR1-MMA 10% and NOA81, 1.725 and 1.560 respectively, were measured with a THz spectrometer based on a photoconductive antenna whose bandwidth typically covers the (0.1 -1.5) THz range [19].In our calculations we assumed that these values remained constant through the entire (0.1-10) THz range.We must account for THz losses, especially for broadband applications and devices made out of acrylates.Unfortunately, such information is not available for the UV-curable acrylates used for the core and the cladding.However, from a previous work using a broadband THz spectrometer we reported the index and attenuation for an acrylate with a similar composition, poly-methyl methacrylate (PMMA), over the entire (0.1-10) THz range [33].In our computations, we assumed that both DR1-MMA and NOA81 had the same linear attenuation PMMA exhibits over the (0.1-1.5) THz range.We limit to this range since the attenuation is linear, that is α =2κω/2 with κ = 0.034 constant, which is consistent with our assumption that these materials exhibit constant indices in the entire (0.1-10) THz.For the IR range we used the function n(! ) = 1.5364 + 2.0392 " 10 #2 / ! 2 for DR1-MMA 10%, where λ is given in µm.This function fits the experimental values from 500 nm to 1.2 µm for DR1-MMA 10% and for NOA81 we used 1.557 for the index since it was found that this material is dispersionless at 820 nm [19].The IR losses were neglected for both DR1-MMA and NOA81 since this effect is insignificant for the millimeter long devices modeled here.
For the Gaussian input beam in Eq. ( 9) with w 0 = 1.5 µm that we used in our computations, we need about 5 modes in the mode expansion to completely describe the input mode profile.However, in our implementation only two IR modes were necessary to solve for Eq. ( 7) since the expansion with two terms provides a reasonable approximation, see Fig. 6(a).About 70% of the input IR pump couples to the fundamental IR mode and 7% couples to the first order mode.The waveguide structure has a single guided THz mode, this fundamental (m = 0) TM THz mode is a TEM like mode since the mode profile for the THz frequencies is uniform in the core, see Fig.The THz output produced by our model for Cao's structure for the 3 waveguide devices they fabricated, with propagation length L = 1, 2 and 3 mm are shown in Fig. 7.The time-domain waveforms are obtained by computing the Fourier transform in Eq. ( 3) after being convolved with a detector response function for a mode expansion with one THz mode and two IR modes.The fields are evaluated at the center of the waveguide (x = 0).To reproduce the experimental conditions we had to account for the response of the detector used in their experimental setup.In Cao's setup a photo-conductive dipole antenna limited the detected signal to approximately 2 THz, yet our model predicts a bandwidth up to 6 THz.Since we did not have access to the response function of this antenna we assumed a response function R(!) = 1 / 2Erfc(r(!" !0 )) where, Erfc(x) is the complementary error function, Ω 0 /2π = 1.5 THz and r = 1/(1 THz) is a constant that provides a soft and continuous decay from 1 to zero.With this response function the bandwidth of the convolved THz spectra is limited to approximately 2 THz.In Fig. 7(b) we show the response function R(Ω), along with the complete generated and corresponding convolved bandwidths for the 3 devices.The waveforms shown in Fig. 7(a) correspond to the convolved bandwidths and since the response function is not precisely the one used for the detector the exact structure of the pulse does not correspond to the experimental results.However, the ratio of peak-to-peak field amplitudes obtained by our model for the 1 mm: 2 mm: 3 mm waveguides in the Cao devices, are 1:1.57:1.75 which are in very good agreement with the experimental results 1:1.56:1.77.The nonlinear conversion efficiency we predict for these devices is on the order of ~ 10 -6 W -1 .

Fig. 1 .
Fig. 1.(a) Symmetric five-layer slab waveguide.The arrows represent the dipole moment of the chromophores.(b) Second order susceptibility, χ xxx (2) , vs. wavelength for DAPC in the telecom range.Inset shows the chemical structure of the chromophore used in DAPC.

Fig. 2 .
Fig. 2. (a) Refractive index and extinction coefficient for Al at THz frequencies (Ω).(b) Extinction coefficient and (c) refractive index for Al at IR wavelengths (λ).The dots represent experimental data and the lines correspond to the fits in the model.A global Drude fit is used in the THz range whereas interpolation fits are used for the IR range.
) and 3(c) show the modal effective index N eff,m and Figs.3(b) and 3(d) show the mode attenuation α m .These results are obtained by finding all the solutions (β m 's) to the dispersion equation (Eq.(20) in Appendix A) repeatedly for 250 points at both the THz frequency range !" (0.1-20) THz

Fig. 3 .
Fig. 3. (a) Effective index Neff,m and (b) respective mode attenuation αm, for even modes in the IR wavelength range (λ) for a DAPC structure with d = 1.579 µm, t = 20 µm and PS cladding.The horizontal dashed lines enclose the bandwidth of the IR pump used in our calculations.For structures with DAPC cores (d = 1.579 µm for PS cladding and d = 1.274 µm for TOPAS claddings), we show (c) the effective index Neff,m and (d) the respective mode attenuation αm, for the TEM-like mode for THz frequencies (Ω) for a range of PS cladding thicknesses (t=2.3-2.7 µm) and for a TOPAS cladding where t=1.9 µm.The index dispersion for the core and cladding are shown for comparison.Core and cladding losses are neglected for the IR modes and αm for m =1 in (b) is scaled for comparison, the scaling factor is indicated next to the curve.

Fig. 4 ( 4 PS :-- 2 . 3 9
Fig. 4 (a) Group index n g,IR for the m = 0 IR mode and the value of Neff,0 for the TEM-like mode at 5 THz for structures with different half thickness values t.For the TOPAS structure d = 1.274 µm and d = 1.579 µm for the PS structure.(b) Nonlinear efficiency ηTHz vs emitter length L for the different structures.

Fig. 5
Fig. 5 For a structure with d = 1.579 µm and t = 2.5 µm and PS cladding: (a) Spectral power densities dPTHz vs. THz frequencies for different emitter lengths, L. (b) Output electric field at the center of the waveguide produced by the most efficient structure we found. 6(b).

Fig. 6 .
Fig. 6.(a) Gaussian input profile (w 0 = 1.5 µm) for the IR pump beam at 820 nm and corresponding mode decompositions using 1 and 2 modes and (b) mode profile for the fundamental TEM-like mode at 10 THz for a structure with (t = 5 µm, d = 1.5 µm)

Fig. 7 .
Fig. 7. Output THz electric field for a device with (d = 1.5 µm, t = 5 µm) for different propagation lengths L. (a) THz waveforms after being convolved with the detector response function R(Ω).(b) The corresponding THz spectral fields before and after being convolved with R(Ω), shown for comparison.
2. To solve for the unknown coefficients A i , B i H must be continuous at each one of the boundaries at x i , E y , H z { } for TE modes and {H y , E z } for the TM modes.This is equivalent to