Purcell effect in nonlinear photonic structures : A coupled mode theory analysis

We develop a coupled mode theory (CMT) model of the behavior of a polarization source in a general photonic structure, and obtain an analytical expression for the resulting generated electric field; loss, gain and/or nonlinearities can also be modeled. Based on this treatment, we investigate the criteria needed to achieve an enhancement in various nonlinear effects, and to produce efficient sources of terahertz radiation, in particular. Our results agree well with exact finite-difference time-domain (FDTD) results. Therefore, this approach can also in certain circumstances be used as a potential substitute for the more numerically intensive FDTD method. © 2008 Optical Society of America OCIS codes: (350.4238) Nanophotonics and photonic crystals; (190.4720) Optical nonlinearities of condensed matter; (260.2110) Electromagnetic optics References and links 1. P. A. Franken, A. E. Hill, C. W. Peters, and G. Weinreich, “Generation of Optical Harmonics,” Phys. Rev. Lett. 6, 118-119 (1961). 2. M. Bass, P. A. Franken, A. E. Hill, C. W. Peters, and G. Weinreich, “Optical Mixing,” Phys. Rev. Lett. 8, 18-18 (1962). 3. J. A. Giordmaine, “Mixing of Light Beams in Crystals,” Phys. Rev. Lett. 8, 19-20 (1962). 4. P. D. Maker, R. W. Terhune, M. Nisenhoff, and C. M. Savage, “Effects of Dispersion and Focusing on the Production of Optical Harmonics,” Phys. Rev. Lett. 8, 21-22 (1962). 5. J. A. Armstrong, N. Bloembergen, J. Ducuing, and P. S. Pershan, “Interactions between Light Waves in a Nonlinear Dielectric,” Phys. Rev. 127, 1918-1939 (1962). 6. N. Bloembergen and Y. R. Shen, “Quantum-Theoretical Comparison of Nonlinear Susceptibilities in Parametric Media, Lasers, and Raman Lasers,” Phys. Rev. 133, A37-A49 (1964). 7. G. Eckardt, R. W. Hellwarth, F. J. McClung, S. E. Schwartz, D. Weiner, and E. J. Woodbury, “Stimulated Raman Scattering From Organic Liquids,” Phys. Rev. Lett. 9, 455 (1962). 8. G. P. Agrawal, A. Hasegawa, Y. Kivshar, and S. Wabnitz, “Introduction to the Special Issue on Nonlinear Optics,” IEEE J. Sel. Top. In Quantum Electron. 8, 405-407 (2002). 9. J. Squier and M. Müller, “High resolution nonlinear microscopy: a review of sources and methods for achieving optimal imaging,” Rev. Sci. Instrum. 72, 2855-2867 (2001). 10. Yuri Kivshar, “Spatial solitons: Bending light at will,” Nature Physics 2, 729-730 (2006). 11. Günter Steinmeyer, “Laser physics: Terahertz meets attoscience,” Nature Physics 2, 305-306 (2006). 12. J. M. Dudley, C. Finot, D. J. Richardson, and G. Millot, “Self-similarity in ultrafast nonlinear optics,” Nature Phys. 3, 597-603 (2007). #97787 $15.00 USD Received 24 Jun 2008; revised 17 Jul 2008; accepted 18 Jul 2008; published 4 Aug 2008 (C) 2008 OSA 18 August 2008 / Vol. 16, No. 17 / OPTICS EXPRESS 12523 13. B. Ferguson and X.-Cheng Zhang, “Materials for terahertz science and technology,” Nature Materials 1, 26-33 (2002). 14. E. Mueller, “Terahertz radiation sources for imaging and sensing applications,” Photonics Spectra 40, 60 (2006). 15. E. Purcell, “Spontaneous emission probabilities at radio frequencies,” Phys. Rev. 69, 681 (1946). 16. M. Soljac̆ić, and J. D. Joannopoulos, “Enhancement of non-linear effects using photonic crystals,” Nature Materials 3, 211-219 (2004). 17. G. DAguanno, M. Centini, M. Scalora, C. Sibilia, Y. Dumeige, P. Vidakovic, J. A. Levenson, M. J. Bloemer, C. M. Bowden, J. W. Haus, and M. Bertolotti, “Photonic band edge effects in finite structures and applications to χ(2) interactions,” Phys. Rev. E. 64, 016609 (2001). 18. A. Taflove, Computational Electrodynamics: The Finite-Difference Time-Domain Method (Artech House, Norwood, 1995). 19. J. Bravo-Abad, A. Rodriguez, P. Bermel, S. G. Johnson, J. D. Joannopoulos, and M. Soljac̆ić, “Enhanced nonlinear optics in photonic-crystal microcavities,” Opt. Express 15, 16161-16176 (2007). 20. T. Laroche, F. I. Baida, and D. Van Labeke, “Three-dimensional finite-difference time-domain study of enhanced second-harmonic generation at the end of a apertureless scanning near-field optical microscope metal tip,” J. Opt. Soc. Am. B 22, 1045-1051 (2005). 21. K. Sakoda and K. Ohtaka, “Optical response of three-dimensional photonic lattices: Solutions of inhomogeneous Maxwells equations and their applications,” Phys. Rev. B 54, 5732-5741 (1996). 22. K. Sakoda and K. Ohtaka, “Sum-frequency generation in a two-dimensional photonic lattice,” Phys. Rev. B 54, 5742-5749 (1996). 23. H. A. Haus, Waves and Fields in Optoelectronics (Prentice-Hall, New Jersey, 1984). 24. J. D. Joannopoulos, R. D. Meade, J. N. Winn, Photonic Crystals: Molding the Flow of Light (Princeton University Press, 1995). 25. C. Luo, M. Ibanescu, E. J. Reed, S. G. Johnson, and J. D. Joannopoulos, “Doppler Radiation Emitted by an Oscillating Dipole Moving inside a Photonic Band-Gap Crystal,” Phys. Rev. Lett. 96, 043903 (2006). 26. D. H. Auston, K. P. Cheung, J. A. Valdmanis, and D. A. Kleinman, “Cherenkov Radiation from Femtosecond Optical Pulses in Electro-Optic Media,” Phys. Rev. Lett. 53, 1555-1558 (1984). 27. M. Ibanescu, E. J. Reed, and J. D. Joannopoulos, “Enhanced Photonic Band-Gap Confinement via Van Hove Saddle Point Singularities,” Phys. Rev. Lett. 96, 033904 (2006). 28. 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). 29. T. Tanabe, K. Suto, J.-ichi Nishizawa, K. Saito, and T. Kimura, “Frequency-tunable terahertz wave generation via excitation of phonon-polaritons in GaP,” J. Phys. D: Applied Physics 36, 953-957 (2003). 30. S. G. Johnson and J. D. Joannopoulos, “Block-iterative frequency-domain methods for Maxwell’s equations in a planewave basis ,” Opt. Express 8, 173-190 (2001). 31. R. W. Boyd, Nonlinear Optics (Academic Press, 2002). 32. A. Farjadpour, D. Roundy, A. Rodriguez, M. Ibanescu, P. Bermel, J. D. Joannopoulos, S. G. Johnson, and G. Burr, “Improving accuracy by subpixel smoothing in the finite-difference time domain,” Opt. Lett. 31, 29722974 (2006).


Introduction
Since the emergence of nonlinear optics [1] in 1961, and the major breakthroughs [4,3,2,5,6,7] that subsequently marked its growth, the field of nonlinear optics has been producing continuous scientific excitement: Numerous nonlinear optical phenomena have been discovered, and have significantly impacted scientific progress in many respects [9,12,10,11,8].In particular, the generation of terahertz radiation via nonlinear optical techniques, such as optical rectification, difference frequency generation and parametric generation, presents unique features [13,14] that are also proving to be very useful tools for biomedical imaging, sensing, spectroscopy, sample characterization, etc.This wealth of applications, made possible by nonlinear phenomena, has necessitated the perpetual quest for improving the efficiency of nonlinear optical processes.In addition to phase matching considerations, the Purcell effect [15], well known for its ability to control spontaneous emission rates, provides a hint on how to boost nonlinear conversion efficiencies, by modifying (increasing) the local density of photonic states (LDOS) at the location of the nonlinear polarization source.A light source embedded in a photonic structure is known to emit faster when the LDOS at the source frequency is larger, and this applies for light sources originating from nonlinear polarization as well.Photonic crystals constitute a versatile tool for tailoring LDOS, and hence they offer extraordinary opportunities to enhance nonlinear effects [16], also by exploiting the Purcell effect.
So far, a few proposals have been made to exploit the Purcell effect in order to improve the nonlinear optical response of photonic crystals.For instance, the effective medium approach has been used to investigate how the density of states affects the conversion efficiency for second harmonic generation in a one dimensional photonic band gap (PBG) structure [17].The finite-difference time-domain (FDTD) method [18] could also be used to solve the problem numerically [19,16,20].However, a better comprehension of the problem requires a more elaborate treatment, that elucidates the implications of the Purcell effect on efficient nonlinear conversion.A method using Green's function has been developed [21,22] to calculate the nonlinear optical response of a photonic crystal to an external polarization source.In this respect, temporal coupled mode theory [23] (CMT) seems to be of great promise, since it provides a physically intuitive framework that easily addresses the problem, where the photonic structure can also exhibit any sort of loss or gain.
In this paper, we use CMT to study the response of a general photonic structure to a traveling polarization source of arbitrary extent.As long as CMT is still valid (i.e. as long as the quality factors are large enough), absorptive and/or radiative losses, as well as gain, can be modeled.Given the current challenge of making high-power terahertz sources, we illustrate our general approach, by proposing, for the first time, a method for efficient generation of terahertz waves by optical rectification in a two-dimensional (2D) photonic crystal.For this purpose, we theoretically investigate sending an optical beam into an appropriately designed 2D nonlinear photonic crystal, and we use our CMT-based result to calculate the total energy radiated at terahertz frequencies.Next, we repeat the same calculation for the terahertz generated energy, but now using the exact FDTD method instead.The obtained close agreement between the CMT and FDTD results, not only validates our CMT approach, but it also suggests that our CMT treatment could serve as a substantially less numerically intensive alternative to the FDTD method, for certain systems of interest.For instance, in problems involving frequencies and wavelengths that range over many orders of magnitude, FDTD calculations are very difficult to perform.However, it is far less intricate to solve such problems by using our CMT-based formalism.A specific example of such cases is the generation of terahertz radiation by optical rectification, where a proper FDTD calculation should involve simultaneously wavelengths in both the optical and terahertz regimes, resulting in a computational burden that could go beyond the capabilities of currently available computers.Such complications, however, are not present if one uses our CMT formalism, instead.
The structure of the paper is as follows: In Sec.II, we develop the CMT model and calculate the generated electric field in the most general case.Then, in Sec.III, we discuss how the Purcell effect follows from our general result, by calculating the power emitted by a point dipole source.Finally, in Sec.IV, we apply our work to the specific problem of THz generation and compare our results with FDTD results.

Coupled-mode-theory model
To start with, we suppose that there exists a polarization source P( r,t) inside a photonic structure.Let ( E ν , H ν ) label a mode of the source-free solutions to Maxwell's equations obtained by using linear real indices of refraction of the photonic structure; while calculating the modes ( E ν , H ν ), we assume that the photonic structure does not involve any loss, gain or nonlinearities.The effects of loss/gain will be addressed perturbatively below through the use of CMT.Nonlinear effects that are of interest for THz generation are included in the sense that the polarization source itself is generated through nonlinear effects, starting from some external electric fields of different frequencies; this polarization source, in turn excites the modes of the structure.{ D ν ( r)/ ε( r)}, where D ν ( r) = ε( r) E ν ( r), form a complete set [24], since they are eigenmodes of the hermitian operator 1  √ . The polarization source P( r,t) induces in the structure electromagnetic fields E( r,t) and H( r,t).If we denote the electric displacement vector by D( r,t), then D( r,t)/ ε( r) = ε( r) E( r,t) can be expressed as a linear superposition [24] of D ν ( r)/ ε( r) = ε( r) E ν ( r).Hence, we can write where a ν is the amplitude of the mode labeled by ν, normalized such that |a ν | 2 is the total energy in that particular mode.A note about the mode labeling is in place.For a uniform medium, the modes are labeled by the wavevector k and polarization σ , hence ν is to be identified with ( k, σ ) in this case.For a photonic crystal, it is most convenient to label the modes by a band index n, a polarization σ , and a Bloch wavevector k that lies in the first Brillouin zone.Therefore, we identify ν ≡ (n, k, σ ) for photonic crystals.For a finite photonic crystal structure, the allowed values of k consistent with boundary conditions, are discrete.However, for infinite photonic crystal structures, there is a continuum of allowed values of k, and the sum in Eq. ( 1) transforms into a discrete sum over n and an integral over k.
Treating each mode as a CMT port, the CMT equation for the mode amplitude a ν becomes [23] da where ω ν is the frequency of the mode labeled by ν, Γ ν rad and Γ ν abs are the rates of radiative (out of the structure) and absorptive decay, respectively, and Γ ν g is the rate of gain.κ ν s ν + is the square root per unit time of the portion of the polarization source's energy, that couples to the photonic structure; i.e. this term models the excitation of the mode of the structure by the polarization source.From Poynting's theorem, κ ν s ν + is given by where the current density source J(r,t) is related to the polarization P( r,t) by J(r,t) ≡ ∂ P( r,t) ∂t.The general solution to Eq. ( 2) can be easily obtained by multiplying both sides by an integrating factor Hence where t o is any reference time preceding the turn-on of the polarization source, and t > t o .Plugging this into Eq.( 1), the electric field induced in the photonic structure becomes The total power radiated out of the photonic structure, at time t, is

and the flux through a surface A at time t is
where E( r,t) is given by Eq. ( 6), and H( r,t) can be obtained from Faraday's law: ∂t, μ 0 being the magnetic permeability.

Connection with Purcell effect
To get more intuition on the above result, let us consider an oscillating point dipole of moment p, embedded at position r o in the structure.If we denote the dipole's angular frequency by ω s , then the polarization is P( r,t) = pδ (3) ( r − r o )e −iω s t .Let us further assume that Γ ν g = 0 ∀ ν (no gain).In this case, the resulting mode amplitude a ν becomes Denoting Γ ν = Γ ν rad + Γ ν abs , the total radiated power In the limit of zero loss If we further assume that there is perfect polarization match between the source and all the modes ν, then p ∀ ν, and where is the local density of photonic states at r o .Therefore, the total radiated power (even from a nonlinear source) is proportional to the local density of states, as expected from the Purcell effect.It also increases quadratically with the dipole moment of the polarization source.

Connection with Doppler radiation in a PhC crystal
To verify further the validity if our approach, we apply our general result (Eq.( 6)) to another well-known special case: the radiation emitted by an oscillating dipole moving with a fixed velocity v in a PhC [25].In this case, the current density can be simply expressed as Plugging this into Eq.( 6), and assuming that Γ ν = 0 ∀ ν, we obtain As we mentioned before, the modes of a PhC are labeled by ν ≡ (n, k, σ ), and from Bloch theorem E n k ( r) = e i k. r u n k ( r), where u n k ( r) has the periodicity of the PhC, and hence only reciprocal lattice vectors G appear in its Fourier expansion u n k ( r) = ∑ G e n k, G e i G. r .Hence Eq. ( 14) becomes which simply evaluates to which is consistent with Eq.(1) in ref (25) in the limit ω n kσ ω s + ( k + G). v. Thus, while our formalism properly reproduces the well-known predictions of the Purcell effect (Eq. ( 11)) and of Doppler radiation in PhCs (Eq.( 16)), its scope extends further, because it analyzes in a systematic and quantitative way the behavior of an arbitrary polarization in any photonic structure.Note that the polarization source can have complicated spatiotemporal variations, and is not restricted to static collections of point dipoles.
A glance at Eq. ( 6) enables us to recognize the various possibilities for boosting the magnitude of the electric electric field E( r,t) induced in the structure by a given external polarization source.Clearly, we need to have the largest possible density of modes that yield a significant value of the summand on the right hand side of Eq. (6).To achieve such large values of the summand, the phase-matched modes ought to have their resonance frequency ω ν close enough to the frequency ω s of the source.Their polarization should also be matched to the source's polarization, and their spatial overlap with the source's extent should be considerable.This latter requirement suggests that the modes should be highly concentrated at the position of the source.And of course, we should attempt to reduce losses to the lowest possible level.

Optimizing the structure
Having identified the essential features required to enhance the response of a photonic structure to an external polarization source, and motivated by the existing need for efficient high-power terahertz sources, we apply our formalism to the specific problem of terahertz generation by optical rectification.To this end, we select a nonlinear photonic crystal (PhC) structure that features a χ (2) nonlinearity, and that satisfies as many as possible of the criteria stated in section III.Before getting into the details of the structure, let's outline briefly the fundamental principles.We propose to excite our PhC structure with a beam centered around an optical frequency; we denote the optical electric field by E opt .By optical rectification, this optical beam generates a polarization P T Hz ∼ χ (2) E opt E * opt in the terahertz frequency range [26].According to our formalism, the rate at which this polarization P THz excites terahertz waves depends on several considerations pertaining to the PhC structure.First, one of course strives to use a material with χ (2) as large as possible.Next, another option is to increase the magnitude of E opt as much as possible; ultimately there is a limit to which one can pursue this approach, imposed by the optical breakdown threshold.Finally one can adjust the properties of the PhC structure at terahertz frequencies, to achieve high density of states and good spatial overlap between the terahertz modes and the source P T Hz .
A promising PhC structure, in this respect, has been proposed by some of the current authors [27].Here we consider a very similar structure consisting of a two-dimensional (2D) square lattice of rods in air.The spacing between the rods is denoted by a, the nonlinear susceptibility by χ (2) , and the linear refractive index at terahertz frequencies by n r .The rods of radius r are connected along the x-direction by thin waveguides, each of thickness t and consisting of the same material as the rods.These waveguides are needed, in order to provide guiding for the optical light.A unit cell of the 2D periodic PhC structure is depicted in Fig. 1(a) as a color contour plot of the linear dielectric function ε(x, y).A similar contour plot is presented in Fig. 1(b), showing 5 units of the structure in each direction, together with the optical beam.If we choose r = 0.13a, t = 0.04a and n r = 3.5 (close to the refractive index of GaP [28,29], commonly used for terahertz generation by optical rectification), the second band (n = 2) of the transverse magnetic (σ =TM) modes is characterized by a saddle point where the band is ultraflat, and consequently the density of states is enhanced.The band structure was computed by preconditioned conjugate-gradient minimization of the block Rayleigh quotient in a planewave basis, using a freely available software package [30].The projected band diagram for σ =TM, shown in Fig. 2 along with a color contour plot of the second band, indicates that the second TM band is narrowest at k x = 0.1559(2π/a), where it has a width of only 1.27% around the central frequency f = 0.509(c/a); thereby if the targeted frequency is e.g.f = 1 THz, one needs to choose a = 152.7μm.So, the PhC structure we have chosen indeed has a considerable number of modes polarized along z (parallel to the rods), with k x = 0.1559(2π/a), and by choosing a appropriately, we can tune the structure to be optimized for any frequeny in the THz regime.Since we are interested in the very small frequency range around the saddle point of the second band, we consider a narrow-bandwidth excitation.To this end, we assume that the optical beam is sent through the waveguide centered at the origin, and is particularly chosen such that the current density J T Hz associated with the polarization source P THz , has the form: where k s x = 0.1559(2π/a) to ensure phase matching with the modes at the ultraflat portion of the second band.The angular frequency of the terahertz polarization source is ω s = 0.509(2πc/a), and the values of the remaining parameters are: ζ = 0.02a and τ = 100(a/c).The time t is expressed in units of a/c.Since J T Hz points along ẑ, all the TM modes have their polarization perfectly matched to that of the source.For convenience in subsequent calculations, and because J T Hz (x, y;t) is separable in space and time, we write it as where and We also set J T Hz (x, y;t) = 0 for t < 0 and t ≥ 1000a/c.Up to this point most of the criteria for efficient THz generation have been met, and we are left with the issue of the spatial overlap between the modes (n = 2; k x = 0.1559(2π/a), k y ; σ =TM) and the terahertz polarization source.The degree of localization of the modes E n=2;(k x=0.1559(2π/a) ,k y ;σ =T M) ( r) at the source's extent can be checked by looking at the color contour plot of w k y (x, y) ≡ ε(x, y) E n=2;(k x=0.1559(2π/a) ,k y ;σ =T M) (x, y) 2 (21) for several k y values.We computed the fields E n=2;(k x=0.1559(2π/a) ,k y ;σ =T M) ( r) by using the same software package [30], and calculated w k y (x, y) for different k y values.A representative plot is shown in Fig. 3 for k y = 0. Despite the good localization of the modes of interest close to the source's extent, this does not guarantee an optimum for the overlap integral for n = 2. Indeed, for most of the different possible k y in the first brillouin zone, the integrand of the overlap integral O n=2;(k x =k s x ,k y ) oscillates along x, and hence integrates to a negligible value.That is, for n = 2 and for many k y , the integrand is observed to flip the sign at x ±0.25a.One way to circumvent this problem is to periodically pole the structure every half period.More specifically, we propose to flip the sign of χ (2) at x = ±(0.25a+ 0.5a) where is a zero or positive integer, as illustrated in Fig. 4. Mathematically, this corresponds to multiplying the  (2) is positive, while in the blue portions, the crystal's orientation is flipped resulting in a negative χ (2) .integrand in Eq. ( 22) by a factor q(x ), which is +1 in regions where χ (2) is positive, and −1 where χ (2) is negative.In this way, the integrand of the overlap integral O n=2; k preserves the same sign for most of x, and therefore the integral evaluates to a substantially higher value than without any poling.Although this looks similar to quasi-phase matching [31], one should keep in mind that our motivation behind poling was to prevent the overlap integral of the 'second' band from vanishing.For the first band, the integrand of the overlap integral O n=1; k does not change sign along x, and so we wouldn't have to pole the structure if we were interested in modes of the first band.Similarly, if we were interested in modes of the third band, we would have to use a poling configuration different from that used in the second band.Also note that although periodic poling leads to much more efficient coupling between the nonlinear polarization source and the second-band modes of the PhC structure, all the linear properties on which the numerical calculations are based, remain intact.

Calculation of generated energy
Let us now calculate the total emitted energy at THz frequencies.First, we need to calculate the electric field E( r,t) induced by the THz current source (Eq.( 17)) in the periodically poled PhC structure, using Eq. ( 6).We consider a time of t = 1010(a/c) (which is 10a/c later than source turn-off), and we assume that a n k = 0 ∀ (n, k) prior to source turn-on at t = 0; this corresponds to the case where none of the modes of the PhC were excited before t = 0. We use only modes with k x = 0.1559(2π/a) = k s x , and we further assume Γ n k g = Γ n k abs = 0. From Eq. ( 6), it is evident that we cannot proceed before computing the following two functions of k y O poled all space (k y ) ≡ all space and  The integral in Eq. ( 23) can be simplified by noting that in a PhC, the electric fields can be written as where u n k ( r) has the same periodicity as the PhC.Consequently, O poled all space (k y ) becomes The integrand of this equation is periodic function of x .Hence the integral over all x simplifies to the number N x of periods in the x-direction times the integral over a single period, with x ranging from −a/2 to a/2.The integration over y can also be taken to range between −a/2 and a/2, because the optical beam is sent through the central waveguide only, and thus J T Hz space is zero for y outside the interval To evaluate this integral, we discretize space with a resolution of 256 gridpoints/a, and calculate the TM fields u n=2;(k s x ,k y ) ( r) based on the TM modes E n=2;(k s x ,k y ) ( r) computed by using the software package [30].Note that the TE modes have their electric field polarized in the xy plane, and hence they don't get excited by the polarization source.So, in all what follows, when we refer to PhC modes, we consider only the TM modes, and we omit the notation σ = T M. While computing the modes, we make sure that all of them have their phases fixed relative to each other, and that they are all normalized in the same way, e.g.such that A plot of the overlap integral O poled one period (k y ), as a function of k y , is shown in Fig. 5 (after poling is performed).Clearly, the integral takes values close to maximum for most of the k y 's.The reason for which it vanishes as k y → ±0.5(2π/a) is that those modes have an extended node over the whole source, and consequently the integrand in Eq. ( 28) vanishes almost everywhere, either because the source is zero or because the electric field of the modes is zero.So whether the structure is poled or not, the overlap integral vanishes for modes with k y = ±0.5(2π/a).
Next, we calculate T t (k y ) at t = 1010(a/c), using the previously calculated band structure.Eq. ( 6) becomes Since the PhC is of infinite extent in the x and y directions, the k y values consistent with periodic boundary conditions, become dense enough that to a good approximation, the discrete sum over k y can be converted into an integral x(a) To calculate the total emitted energy, we consider calculating the electric field E( r,t) at time t = 1010(a/c), over a region of space large enough in the y direction, hoping that no energy would have left it by t = 1010(a/c).So, we take our computational domain to be a 2D box of size 1a along x and 80a along y.To save on computational memory, we use a spatial resolution of 64 gridpoints/a only, and calculate the TM fields E n=2;(k s x ,k y ) for 201 equally spaced values of k y ranging between −π/a and π/a, again using the software package [30].Next we multiply each field labeled with k y , by O poled one period (k y ) • T 1010(a/c) (k y ), and sum the resulting fields over all values of k y .Finally, we multiply the result by a/π and obtain the THz electric field induced by the optical beam in the PhC at time t = 1010(a/c).Note that we attempted to use a denser grid of k y values (334 equidistant k y points), without change in the final result.The electric energy density profile ε(x, y)|E(x, y,t)| 2 , at t = 1010(a/c), can now be simply calculated over the computational domain, and is presented in Fig. 6.
As we intended, all the emitted energy is still confined inside our large computational domain; more specifically, the energy density is non-vanishing only in the regions y ∈ ± [8a, 19a].Thus, to calculate the total emitted energy, it is sufficient to integrate the energy density at t = 1010(a/c) over a 2D box of size 1a along x, and large enough in the y direction to enclose the regions y ∈ ± [8a, 19a].Setting a = 1, we obtain a value of 0.7108 for the total THz energy emitted in the PBG.To enable comparison with the FDTD results, we normalize the total emitted energy and express it in dimensionless units; that is, we divide the above-mentioned integral of the energy density in the poled PhC case, by the same quantity (0.1839) similarly calculated in an unpoled bulk of the same nonlinear material as that used in the PhC structure.We obtain a dimensionless value of 3.86.
To confirm the validity of our analytical model in general, and of the result for the total emitted energy in the specific example of terahertz generation, we repeat exactly the same calculation for the terahertz emitted energy, but now instead of using our analytical formalism, we employ the finite-difference time-domain (FDTD) method [18], using a freely available software package with subpixel smoothing for increased accuracy [32].We consider a computational cell (centered at the origin) of size 1a along x, and 90a along y, and we discretize the structure with a resolution of 128 pixels/a.We send the same source described before through the waveguide centered at y = 0, and we impose Bloch-periodic boundary conditions along x, with a wavevector k x = 0.1559(2π/a).At the boundaries in the y direction, we set up perfectly matched layers (PML), each of thickness 3a.We simulate the effect of periodic poling by explicitely reversing the sign of the polarization source at x = ±0.25a.To compute the total emitted energy at terahertz, we record the time evolution of the energy E in a box of size 1a (along x)×80a (along y), centered at the origin.We also place flux calculation planes at y = ±40a, and compute the flux that leaves through each of the flux planes as a function of time.The total emitted energy at a particular time t is then given by the sum of the integral up to time t of the net flux through the flux planes, and the energy remaining in the box surrounded by the flux planes, at time t.The time development of the energy E , and of the time integral of the net flux, is shown in Fig. 7, together with the total emitted energy, which is given by their sum.Again, we normalize the total energy emitted in the poled PhC to that emitted in the corresponding unpoled bulk, and we get a value of 3.94 in dimensionless units.Thus, our result for the normalized emitted energy at terahertz, differs only by 2% from the exact result computed by FDTD.To gain even more confidence in the validity our analytical result, we perform FDTD calculations of the energy density profile in the poled PhC, at t = 1010(a/c), and we show the result as a color contour plot in Fig. 8.The agreement between Fig. 6 and Fig. 8 is indeed remarkable; not only do we get a coincidence of the intervals y ∈ ± [8a, 19a] in which the energy density is nonvanishing, but also the waveguides at which the energy density is maximum occur at exactly the same position, namely y = ±13a, according to both methods (FDTD and CMTbased analytics).In addition to validating our analytical CMT-based formalism, the agreement of the results obtained from our analytical model with the exact FDTD results, suggests that our approach would work as a simpler alternative to the numerically intensive FDTD method.Our procedure has the advantage of being far less demanding than the brute-force FDTD technique, in terms of the computational time and resources, especially in problems involving frequencies that range over many orders of magnitude, such as terahertz generation by optical rectification.To get a concrete estimate, we mention that, although we assumed the terahertz polarization source to be given from the beginning, and dealt with terahertz frequencies only, the calculation of energy density profile shown in Fig. 6 took around 15 minutes on a single processor, while the FDTD calculations resulting in Fig. 7 took around 10 hours using 360 processors on a supercomputer.Finally, the FDTD calculations of Fig. 8 took 10 minutes using 32 processors on a supercomputer when a spatial resolution of 64 gridpoints/a was used, and 2 hours using 360 processors when the spatial resolution was 128 gridpoints/a.Note that the only assumption on which our analytical formalism is based is that CMT be valid, meaning that the rates Γ ν be much smaller than the frequency ω ν for each mode ν of the photonic structure [23].

Conclusion
In conclusion, we developed an analytical model for calculating the electric field produced by an external polarization source in a photonic structure.Then, we investigated how the implications of this formalism can provide an insight on enhancing the efficiency of nonlinear effects.We illustrated the procedure, by applying it to the specific case of terahertz generation by optical rectification in a 2D nonlinear photonic crystal, and finally we checked our results against exact calculations based on the FDTD method.The results are in a good agreement with each other, thus validating our approach, and motivating us to propose our technique as a potentially simpler alternative to FDTD, as long as CMT is applicable.

Fig. 2 .
Fig. 2. (Color online) (a) Projected band diagram of the first three bands.(b) Color contour plot of the second band, showing the saddle point where the band is narrowest.

Fig. 4 .
Fig. 4. (Color online) Periodically poled PhC structure.The red portions correspond to the original orientation of the nonlinear crystal when χ(2) is positive, while in the blue portions, the crystal's orientation is flipped resulting in a negative χ(2) .

Fig. 5 .
Fig. 5. (Color online) Overlap integral O poled one period as a function of k y for modes at the narrowest portion (k x = 0.1559(2π/a)) of the second band.

Fig. 6 .
Fig. 6. (Color online) CMT calculations: (a) Color contour plot of the unnormalized terahertz energy density ε(x, y)|E(x, y,t)| 2 at t = 1010(a/c), in a 2D box of size 1a × 80a.(b) A zoom-in version of the plot in Fig. a, showing more details of the interval y ∈ ±[8a, 19a].Note that the optical beam was originally sent through the waveguide at y = 0.

Fig. 7 .
Fig. 7. (Color online) FDTD calculation for the terahertz emitted energy in the PhC structure (not normalized to the bulk).

Fig. 8 .
Fig. 8. (Color online) FDTD calculations: (a) Color contour plot of the terahertz energy density ε(x, y)|E(x, y,t)| 2 at t = 1010(a/c), in a 2D box of size 1a × 80a.(b) A zoom-in version of the plot in Fig. a, showing more details of the interval y ∈ ±[8a, 19a].Note that the optical beam was originally sent through the waveguide at y = 0.