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, (ω₀ + Ω/2) - (ω₀ - Ω/2), of the two main frequency components of a fs-IR pulse with center frequency, ω₀. We derive a closed expression for the power density spectra, dP_THz, for situations in which OR occurs between the 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.

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, α_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.

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⁻⁴ Watts⁻¹ 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⁻⁴ W⁻¹.

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₀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]

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,

The 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_c for completely phase-matched structures, which is unphysical. In the other limiting scenario with completely phase-matched structures, Δβ_R = 0, and we find $L_c=2\ln ({\alpha }_m/\Delta {\beta }_{IR})/({\alpha }_m-\Delta {\beta }_{IR})$ which agrees with the result in [4].

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.

 

Fig. 1 (a) THz power density spectra, dP_THz, vs. THz frequencies, Ω, computed with the three-wave-approximation for the devices shown in Fig. 5(a) of [1]; the emitter length, L, is indicated with arrows. (b) Ratio of dP_THz obtained through both theories, dP_Full/ dP_Approx. vs Ω,; the results correspond to the same devices in (a). (c) Nonlinear THz conversion efficiencies, η_THz vs. L, obtained for the results in (a) and in Fig. 5(a) of [1].

Download Full Size | PDF

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₃₃ ~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).

 

Fig. 2 (a) Symmetric five-layer slab waveguide, the core is made out of the guest-host EO polymer composite AJTB203/APC, the arrows represent the dipole moment of the chromophores. (b) Second order susceptibility vs. wavelength for AJTB203 in the telecom range. Inset shows the chemical structure of the chromophore.

Download Full Size | PDF

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_g,eff, vs. core half-thickness, 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 n_g,eff matches N_eff,m and, since decreasing the cladding thickness (bigger 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, α_eff. 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.

 

Fig. 3 For structures in Fig. 2(a) with t = 10 μm, and distinct cores sizes (d = 5-10 μm): (a) Effective index Neff,0 (inset respective mode attenuation α₀) vs THz frequencies (Ω), for the fundamental, m = 0, TEM-like THz mode. The dispersion for AJTB203 (attenuation in inset) is indicated with arrows. PS index is ~1.59 remaining flat for the THz frequencies. (b) Effective group indices, ng,eff vs core half-thickness (d), for the fundamental, l = 0, and first excited, l = 1, even IR modes for the pumping wavelength λ₀, indicated with arrows.

Download Full Size | PDF

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_c, computed using both the traditional approach (no loss) and our numerical approach (including loss) solving Eq. (4) for the structures that showed the higher L_c values when pumping at 1567 nm and 1900 nm. The numerical approach always yields a finite L_c 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 L_c for completely phase-matched structures. In Fig. 4(a) we observe how both methods to compute L_c differ significantly for the structure with 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.

 

Fig. 4 For structures in Fig. 2(a) with t = 10 µm: (a) Coherence length, L_c vs THz frequencies (Ω), computed neglecting losses (solid black and red), and with our numerical method including losses (dotted black and red) for structures pumped at 1567 nm (d = 9.5 µm) and 1900 nm (d = 9.125 µm). (b) and (c) coherence length contour plots, L_c vs (Ω, d) with our numerical method pumping at 1567 nm (b) and 1900 nm (c). The plot color range in L_c at 1567 nm is almost twice that at 1900 nm and, the blank regions near the origin are out of range, L_c~Ω⁻¹.

Download Full Size | PDF

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_c for different waveguide parameters can also be used for design purposes since structures that present high L_c are either well phase-matched or have low effective mode attenuation, or both. Thus, from Fig. 4(b) we would expect the structure with 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_c values in 0-4 THz range, for example.

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).

 

Fig. 5 (a), (b) and (c), THz power density spectrum, dP_THz vs. THz frequencies, Ω, for OR between the l = 0 IR mode and the TEM-like (m = 0) THz mode for structures in Fig. 2(a) with t = 10 μm, d (indicated in figure), different emitter lengths, L (indicated with arrows), and pumped at λ₀ (indicated in figure). (d) Nonlinear THz conversion efficiency, η_THz, vs emitter length, L, for the structures in (a), (b) and (c). (e) and (f) efficiency contour plots, η_THz vs. (L, d), for OR mediated respectively by the l = 0 and l = 1 even IR modes pumped at 1567 nm.

Download Full Size | PDF

The 35 x 10⁻⁴ W⁻¹ 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⁻⁴ W⁻¹ 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⁻⁴ Watts⁻¹ 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_g,eff is reduced for the 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)).

 

Fig. 6 For structures in Fig. 2(a) with t = 10 μm: (a) Mode profiles for the l = 0 and l = 1 even IR modes for a structure with d = t (i.e., no cladding) at 1567 nm. (b) OR effective THz attenuation, α_eff = α_THz + Δβ_i, and its different contributions for OR between the l = 0 IR mode and the TEM-like THz mode for a structure d = 10 µm, pumped at 1900 nm. Δβ_IR, is defined by Δβ_i = Δβ_IR - α_THz /2. (c) IR contribution to the imaginary part of the phase mismatch, Δβ_IR vs core half-size, d, for the (l = 1) even IR mode evaluated at 7.5 THz for different pumping wavelenghts. Inset: Δβ_i,IR vs. d, for the l = 0 IR mode, the units are the same as in (c).

Download Full Size | PDF

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 – α_m/2, and since α_eff = α_m + Δβ_i, α_eff = α_m/2 + Δβ_IR. Thus, the two main contributions to α_eff result from the THz effective mode attenuation α_m 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 Δβ_IR is relatively high for the l = 0 IR mode since the high plasmonic IR losses make α_eff significantly larger than α_m. In Fig. 6(c) we show Δβ_IR, as the cladding layers are thinned, $d\to t$,for the 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, α₀(ω₀)< α₁(ω₀), and α_eff ~α_m/2 for both modes. However, as the claddings are thinned Δβ_IR increases rapidly for the l = 0 IR mode as it turns plasmonic when α₀(ω₀) 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 α_eff are much larger for the 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_opt). In Fig. 7(a) we show η_max for different core thickness (or cladding layers) when the OR is mediated through the l = 0 or l = 1 IR mode when pumping either at 1567 nm or 1900 nm. We further define the overall maximum efficiency, η_MAX, as the maximum of η_max(d), indicated by the black stars and the red dots in Fig. 7(a). Note that for OR with the l = 1 IR η_MAX occurs when there are no cladding layers, i.e. η_MAX = η_max(d = t) or simply a metal/core/metal structure. In Fig. 7(b) we plot the overall maximum efficiency, η_MAX, vs. pumping wavelength, λ₀, for OR mediated through either the l = 0 or l = 1 IR mode.

 

Fig. 7 For OR between the l = 0 and l = 1 IR mode (indicated with arrows) and the (m = 0) THz mode for structures with t = 10 μm: (a) Nonlinear THz conversion efficiency, η_max vs core half-thickness, d, for structures with optimal emitter length, L_opt, pumped at 1567 nm (solid lines) and 1900 nm (dotted lines); η_max = η_THz(L_opt,d) is maximum for any L value. (b) Maximum nonlinear THz conversion efficiency, η_MAX vs pumping wavelength, λ. η_MAX is the maximum of η_max for any d, in (a) η_MAX is indicated with asterisks (l = 0) and red dots (l = 1). (c) Overall optimal emitter length, L_OPT, and L_80%, vs λ_0, η_THz decreases to 80% of η_MAX at L_80%.

Download Full Size | PDF

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, λ₀ > 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_g,IR – n_THz ~0.02, near 7 THz. For longer pumping wavelengths, Δ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_OPT, the emitter length for which η_MAX occurs, and, L_80%, the emitter length for which η_THz has decreased to 80% of η_MAX. The gap between L_OPT and L_80% is a measure of how stable the efficiency is to variations in emitter lengths. Similar L_OPT values are obtained by pumping either the l = 0 or the l = 1 but L_80% is considerably larger for structures with cladding layers where OR is mediated by the 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⁻⁴ W⁻¹, 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.

 

Fig. 8 For structures in Fig. 2(a) with t = 10 μm, d = 9.85 μm (thin cladding layers ~150 nm) pumped at 1567 nm: (a) Mode profiles for the l = 0 and l = 1 even IR modes, Gaussian input profile with a 6.4 µm beam waist and corresponding mode expansion with the l = 0 and l = 1 modes. (b) Nonlinear THz conversion efficiency, η_THz, vs emitter length, L, computed with integral DFG CMT [1] using the l = 0 and l = 1 IR modes and, with three-wave-approximation for OR mediated by the l = 0 and l = 1 IR modes. (c) THz power density spectrum, dP_THz vs. THz frequencies, Ω, for different emitter lengths, L computed with integral DFG CMT [1]. The THz power density spectrum computed with the three wave approximation have a similar behavior as observed (c) but reduced by a factor of 1/3 and 1/10 for the l = 0 and l = 1 modes respectively.

Download Full Size | PDF

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, η_THz(L,d), 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.

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⁻⁴ W⁻¹ 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⁻⁴ W⁻¹ 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]

(5)$$d{P}_{THz}(z,\Omega )=4\pi f_{rep}{\displaystyle \sum _{m,m\text{'}}[}({\displaystyle \int d{r}_{\perp }\widehat{z}\cdot ({\widehat{e}}_m\times {\widehat{h}}_{m\text{'}}{}^{\ast }))}{A}_m(z,\Omega )A_{m\text{'}}{(z,\Omega )}^{*}{e}^{i({\beta }_m-{\beta }_{m\text{'}}^{*})z}+c.c.].$$

m and m’, are the THz modes and f_rep, 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 (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/{ϵ}^2)/ϵ\sqrt{\pi }\to \delta (x)$ when $ϵ\to 0$. Then,

(6)$$A_m^{THz}(z,\Omega )=\frac{\sqrt{\pi }i{ϵ}_0\Omega }{{\tau }_0{P}_m}{{\rm K}}_{l,l}^m({\omega }_0)e^{i\Delta \beta ({\omega }_0,\Omega )z/2}\frac{\sin (\Delta \beta ({\omega }_0,\Omega )z/2)}{\Delta \beta ({\omega }_0,\Omega )/2i}{A}_l^{IR}(z,{\omega }_0+\frac{\Omega }{2})A_{l^{\prime }}^{IR}{(z,{\omega }_0-\frac{\Omega }{2})}^{*}.$$

Thus, the THz spectral power density, dP_THz, describing OR between a single-IR-mode l and a single-THz-mode m is,

(7)$$d{P}_{THz}(z,\Omega )=C(\Omega )e^{-({\alpha }_m(\Omega )+\Delta {\beta }_i)z}\frac{(\cos h(\Delta {\beta }_{i}z)-\cos (\Delta {\beta }_{R}z))}{\Delta {\beta }_R{}^2+\Delta {\beta }_i{}^2}.$$

Here $C(\Omega )=4\pi f_{rep}|{ϵ}_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 ${\displaystyle \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

(8)$$\begin{array}{c}c\Delta \beta =[({\omega }_0+\Omega /2)N_{eff,l}({\omega }_0+\Omega /2)-({\omega }_0-\Omega /2)N_{eff,l}({\omega }_0-\Omega /2)-\Omega N_{eff,m}(\Omega )]\\ +i[(\omega +\Omega /2){\kappa }_{eff,l}({\omega }_0+\Omega /2)+({\omega }_0-\Omega /2){\kappa }_{eff,l}({\omega }_0-\Omega /2)-{\kappa }_{eff,m}\Omega (\Omega )].\end{array}$$

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 ω₀

(9)$$({\omega }_0+\Omega /2)f({\omega }_0+\Omega /2)\pm ({\omega }_0-\Omega /2)f({\omega }_0-\Omega /2))=\{\begin{array}{c}2{\omega }_{0}f({\omega }_0)+\frac{{\Omega }^2}{2}{f}^{\prime }({\omega }_0)+O({\Omega }^2)\\ \Omega (f({\omega }_0)+{\omega }_0{f}^{\prime }({\omega }_0))+O({\Omega }^3).\end{array}$$

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

(10)$$\Delta \beta =\frac{\Omega }{c}[n_{g,IR}({\omega }_0)-N_{eff,m}(\Omega )]+i[{\alpha }_l({\omega }_0)(1-\frac{1}{4}\frac{{\Omega }^2}{{\omega }_0{}^2})+\frac{1}{4}\frac{{\Omega }^2}{{\omega }_0}{{\alpha }^{\prime }}_l({\omega }_0)-\frac{{\alpha }_m(\Omega )}{2}].$$

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

(11)$$ϵ(\Omega )={ϵ}_0\left(n_b{}^2+{\displaystyle \sum _{i=1}^6\frac{a_i{\Omega }_i^2}{{({\Omega }_i^2-\Omega )}^2-2i{\Gamma }_i\Omega }}\right).$$

Here n_b is the background refractive index, Ω_i is the i-th resonance frequency and, Γ_i and a_i 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.

Table 1. Fit parameters for the six-oscillator-Lorentz model for AJTB203

View Table

It was necessary to include a term k₀ω 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₀ = 251.5 cm⁻¹ THz⁻¹.

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]

(12)$${(n(\Omega )+ik(\Omega ))}^2={\epsilon }_b+\frac{i}{{\epsilon }_0\Omega }(\frac{\delta e^2\tau /m_e}{1-i\Omega \tau })$$

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 e and m_e are the charge and mass of the electron respectively. The fitting parameters for this global fit were ε_b = 1754.4, δ = 9.85 × 10²² cm⁻³ 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]