1. Introduction

In a uniform refractive index χ⁽²⁾ photonic crystal, pairs of waves interact via the quadratic susceptibility of the material to produce harmonics with sum or difference frequencies. If a wave vector of such a harmonic is phase matched by a reciprocal lattice vector of the photonic crystal then waves generated at different source points throughout the crystal interfere constructively. In a 2-dimensional photonic crystal [1, 2, 3, 4] two of these quasi-phase-matching (QPM) processes [5] can be supported, enabling efficient third harmonic generation (THG) or fourth harmonic generation (FHG).

The efficiency of a QPM harmonic generation scheme depends on the χ⁽²⁾ poling pattern of the photonic crystal. The poling function p(x, y)=±1 describes the direction of the optical axis of the material at position (x, y) to be aligned (p=+1) or anti-aligned (p=-1) with the x ³≡z coordinate axis. Since the quadratic susceptibility is a rank 3 tensor field, those components ${\mathrm{\chi }}_{\mathit{\text{ijk}}}^{\left(2\right)}$ for which an odd number of the indices i, j, k are equal to 3 have their sign determined by p(x, y). Thus, χ^{(2) 333}(x, y)=constant × p(x, y). Since ${\mathrm{\chi }}_{333}^{\left(2\right)}$ is the largest component, this is the one generally used in practice, and the case studied here. For trigonal 3m symmetry materials (eg., lithium niobate) and for an incident field propagating in the x-direction with z-polarization, the quadratic processes involve only ${\mathrm{\chi }}_{333}^{\left(2\right)}$.

A poling function that is periodic in the {x, y}-plane defines a photonic crystal. The photonic crystal will support the required QPM processes provided its unit cell geometry (size and shape) satisfies the corresponding QPM relations (Eq. (2) below). These relations are solved by a 2-parameter family of unit cell geometries. Each such QPM unit cell determines a corresponding photonic crystal reciprocal lattice. Since the poling function can be defined by its Fourier coefficients with respect to the reciprocal lattice, the problem of choosing a poling function splits naturally into two parts: (1) selection of a reciprocal lattice basis (or QPM unit cell); and (2) prescription of the poling function Fourier coefficients. Problem (1) has no obvious formulation in terms of harmonic generation efficiency, and is not considered here (we believe the choice of QPM unit cell should be based on fabrication constraints). On the other hand, problem (2) does have a formulation as a THG/FHG efficiency optimization problem and, remarkably, can be solved in full generality for an arbitrary QPM unit cell. For different choices of QPM unit cell, the resulting optimal poling patterns differ only by a linear transformation of the {x, y}-plane.

2. Phase matching in two dimensions

Two independent QPM processes can be accommodated in a 2-dimensional photonic crystal. Since only a single wave of frequency ω is incident, the harmonic generated in the first process is necessarily of frequency 2ω. For the second process one may phase match either a third harmonic wave (3ω=ω+2ω) or a fourth harmonic wave (4ω=2ω+2ω).

A QPM scheme is generic if at each harmonic only one wave vector is phase matched. There exist symmetric non-generic schemes [6]. We consider generic schemes with phase matching diagrams as in Fig. 1. The QPM harmonic with frequency qω has wave vector k _q=(cosθ_q, sin θ_q, 0) k_q , where k_q =n_qqω/c and n_q =n_e (qω) is the extraordinary refractive index of the crystal at frequency qω. The fundamental propagates in the x-direction, so θ₁=0. The propagation angles θ_q, q=2, 3, 4 are free parameters that may be negative or zero, provided the resulting phase matching vectors G ₁ and G ₂ of Fig. 1 remain linearly independent.

 

Fig. 1. Generic QPM schemes for (a) THG, (b) FHG. Reciprocal lattice vectors G ₁ and G ₂ phase match quadratic interactions between jth harmonic wave vectors k _j.

Download Full Size | PDF

The electric fields are z-polarized and have amplitudes varying slowly with position,

Equations for the amplitudes A_q (x) follow from a perturbation analysis of Maxwell’s equations [7]. Each term in the equations is associated with a QPM relation,

where q=m+n, and the phase matching vector G belongs to the reciprocal lattice of the photonic crystal, so appears in the Fourier expansion of the poling function,

Here x=(x, y), G _ab=a G ₁₀+b G ₀₁, and {G ₁₀, G ₀₁} is the reciprocal lattice basis. The equations for A_q (x) depend on the poling pattern only through a few of its Fourier coefficients, each being associated with a phase matching vector. In the following we suppose the two phase matching vectors are ${\mathbf{G}}_1={\mathbf{G}}_{a_1{b}_1}$ and ${\mathbf{G}}_2={\mathbf{G}}_{a_2{b}_2}$ , and we denote the corresponding poling pattern Fourier coefficients by ${\sigma }_1=p_{a_1{b}_1}$ and ${\sigma }_2=p_{a_2{b}_2}$ .

3. Standard THG solutions

The slowly varying amplitude equations for THG with QPM diagram Fig. 1(a), are

Here χ=$\widehat{\chi }$ ⁽²⁾ ₃₃₃=2d ₃₃ denotes the quadratic nonlinearity, assumed here to be independent of frequency (Kleinman symmetry). For lithium niobate we use χ=82×10^-12m/V [8]. The perturbation theory leading to Eqs. (4)–(6) relies on expansion in the dimensionless small parameter ∊=χE _max, where E _max is a typical strong electric field.

THG solutions solve the system (4)–(6) with initial conditions A ₂(0)=A ₃(0)=0 and A ₁(0) arbitrary. Our task is to find, amongst all THG solutions of all possible THG systems (differing by system parameters σ_j), a solution and corresponding poling pattern, that are in some sense optimal. To do this we first establish the following result: Every THG solution can be obtained from some member of a 1-parameter space of standard solutions, using symmetry transformations of the THG system.

The 1-parameter space of THG standard solutions is defined as follows. The standard solution at parameter value ϕ has standard initial data A ₂(0)=A ₃(0)=0 and $A_1\left(0\right)=\sqrt{\frac{2c{\mu }_{0}u}{n_1}}$ for an incident energy flux of U=30 MW/cm², and solves the THG system for which the Fourier parameters are given by

where the real valued function σ(ϕ) is defined to be as large as possible (subject to σ₁, σ₂ being Fourier coefficients of some poling function). The function σ(ϕ) clearly exists, since by Parseval’s theorem σ(ϕ)²≤∫ p(x)²=1. An efficient numerical method for evaluating σ(ϕ) is given in [6].

THG solutions are physically equivalent for our purpose if they predict the same energy fluxes. The energy flux of the qth harmonic through a plane perpendicular to the x direction is U_q =n_q cos(θ_q)A*_qA_q /(2cµ₀). By Eqs. (4)–(6) the total energy flux U≡U ₁+U ₂+U ₃, satisfies dU/dx=0, so is uniform throughout the crystal.

The THG system (4)–(6) is invariant under the transformations,

Solutions that are related by the phase symmetry (8) have the same energy fluxes. Those for which A ₁(0) are real may therefore be taken as representative. Similarly, we need only consider solutions for which |A ₁(0)| is some conveniently chosen constant because the scale symmetry (9) shows that the effect of increasing |A ₁(0)| by a factor μ is to compress the evolution into a length of crystal μ times shorter.

We may take σ_j real. To see this, note that a translation of the poling pattern introduces phase factors into its Fourier coefficients. If p(x) ↦ p(x-x₀) then

Since ${\mathbf{G}}_{a_1{b}_1},{\mathbf{G}}_{a_2{b}_2}$ are linearly independent, there exits x ₀ such that σ_j are real. THG systems with translated poling patterns are physically equivalent. The associated phase symmetry is (10), by which solutions of complex σ_j systems with standard initial data can be transformed into solutions of real σ_j systems with standard initial data. For the latter, the solutions are such that A ₁(x) and A ₃(x) are real, and A ₂(x) is imaginary.

The scaling symmetry (11) shows that for each direction ϕ in the {σ₁, σ₂}-plane we need only consider one value of σ(ϕ)≡(${\mathrm{\sigma }}_1^2$+${\mathrm{\sigma }}_2^2$)^1/2. Solutions for different σ(ϕ) can be obtained by rescaling in x. Since we want the frequency conversion process to take place over the shortest length of crystal, we suppose σ(ϕ) is chosen as large as possible.

4. Σ-patterns

A poling pattern with the property that some vector of its Fourier coefficients (as in Eq. (7)) is as long as possible, is referred to here as a Σ-pattern. Poling patterns for THG standard solutions are thus (by definition, Section 3) examples of Σ-patterns.

 

Fig. 2. (a) Each point on the curve Σ corresponds to a Σ-pattern. (b) (2.4MB) Video of the Σ-patterns defined by Eq. (14) for varying angle ψ.

Download Full Size | PDF

A construction for Σ-patterns [6] has been used to plot Fig. 2(a), which shows the curve Σ defined by Eq. (7). The explicit general form for Σ-patterns is [6],

Here n is the unit normal to the hypersurface Σ, and δ_x is a Fourier coefficient vector for the periodic extension of the Dirac delta function translated to position x (see [6] for details). In the present case Σ is 1-dimensional and we let its unit normal in the {σ₁, σ₂}-plane be n=(cosψ, sin ψ). Then Eq. (13) evaluates to,

Since p(x) is clearly 2π-periodic in the arguments $\mathbf{x}·{\mathbf{G}}_{a_1{b}_1}$ and $\mathbf{x}·{\mathbf{G}}_{a_2{b}_2}$ we may, without loss of generality, take the phase matching vectors to be G ₁₀ and G ₀₁.

The video linked to Fig. 2(b) shows 4×4 unit cells of the Σ-patterns defined by Eq. (14) as the parameter ψ varies through [0, 2π]. The unit cell has been mapped onto the unit square in the {t ¹, t ²}-plane, where t ¹=x·G ₁₀/(2π) and t ²=x·G ₀₁/(2π).

5. Optimal THG

To examine the space of standard THG solutions we graph the energy fluxes, animated by angle ψ (it is more convenient to parametrize Σ by ψ in Eq. (14) rather than ϕ in Eq. (7)). A video linked to Fig. 3, shows THG energy fluxes for lithium niobate at 140°C and λ=1.536 µm (refractive indices n ₁=2.1434, n ₂=2.1857 and n ₃=2.2510 [9], and we used the approximation cos θ _q =1). One finds a critical value ψ_crit ≈ 41.42° for which lim_x→∞ U_q (x)=0 for q=1, 2. The critical solution is clearly optimal for THG. For a crystal of length x=0.5 cm almost all (99.995%) input power is converted to the 3rd harmonic. Fig. 4 shows the corresponding amplitudes. The Σ-pattern for ψ=41.42° is shown in Fig. 2(b), and has Fourier coefficients σ₁=0.4834 and σ₂=0.3228.

 

Fig. 3. (1.9MB) Energy fluxes for standard THG solutions plotted against distance into the photonic crystal. The video shows how the energy fluxes change as the poling parameter ψ is varied. The red, green and blue curves are U_q (x), q=1, 2, 3 (1st, 2nd and 3rd harmonics) respectively. For the critical value ψ_crit ≈ 41.42°, 100% conversion is attained at x=∞. (a) ψ=41.41°. (b) ψ=41.43°.

Download Full Size | PDF

The critical THG solution exists independently of material and operating conditions. This can be seen by making an appropriate change of variables in Eqs. (4)–(6). Numerically, we find that the critical THG solution is characterized by the condition $\frac{{\sigma }_2}{{\sigma }_1}\approx 0.651780083\sqrt{n_3\cos \frac{{\theta }_3}{n_1}}$ . This critical Fourier coefficient ratio for THG has been explained analytically in [10] (using their notation, Κ₂/Κ₁≈0.651780083√3).

6. Optimal FHG

An analysis similar to that in Sections 3–5 can be given for FHG with phase matching diagram Fig. 1(b). The slowly varying amplitude equations are,

 

Fig. 4. (1.9MB) THG amplitudes near ψ_crit≈41.42°. The red, green and blue curves are Re{A ₁}, Im{A ₂} and Re{A ₃} respectively (the other components are zero for standard THG solutions). (a) ψ=41.41°. (b) ψ=41.43°.

Download Full Size | PDF

Unlike standard THG solutions, those for FHG have no critical dependence on ψ. Instead, the intensity of the 4th harmonic monotonically increases with distance into the crystal and a “best” value for ψ depends on one’s criteria. For example, using the criteria of shortest crystal for 99% conversion, one finds that for lithium niobate at 140° C, λ=1.536 µm, and 30MW/cm² incident intensity, the minimum crystal length is l=0.698 cm at ψ=30.7°. For 95% conversion, l=0.465 cm at ψ=36.6°.

 

Fig. 5. (a) (1.8MB) FHG energy fluxes U_q, q=1, 2, 4. (b) (1.8MB) FHG amplitudes Re{A ₁}, Im{A ₂}, and Im{A ₄}. Red, green and blue curves are for 1st, 2nd and 4th harmonics respectively.

Download Full Size | PDF

7. Conclusion

We have determined optimal poling patterns for THG and FHG assuming perfectly phase matched infinite plane wave harmonics. We have not considered the effects of tuning (phase mismatch) or beam width (walk-off). We find that there exists a critical THG solution for which 100% conversion (at x=∞) is predicted. This solution is also optimal for finite length crystals, and there corresponds a unique poling pattern (the Σ-pattern in Fig. 2(b)) for which the THG process takes place most rapidly. The solutions for FHG and THG are qualitatively different. For FHG we suggest using the optimality criteria that 99% (or 95%) conversion is attained in the shortest length of crystal.

We note there exist infinitely many THG poling patterns for which 100% conversion (at x=∞) is predicted, however, the conversion process is most rapid for the Σ-pattern. For all these patterns the Fourier coefficients have the critical ratio $\frac{{\sigma }_2}{{\sigma }_1}\approx 0.6518\sqrt{n_3\cos \frac{{\theta }_3}{n_1}}$ . Poling patterns such as these will almost certainly be more practical than Σ-patterns for use as photonic crystal designs because the poling process normally results in poled regions with boundaries coinciding with crystal planes of the material, unlike the curved boundaries in Σ-patterns. Nevertheless, the optimal THG and FHG Σ-patterns provide a theoretical baseline for evaluating the efficiency of a photonic crystal design and might also be used as initial design guesses. The support of the Australian Research Council for the Centre of excellence for Ultrahigh bandwidth Devices for Optical Systems (CUDOS) is acknowledged.