Distributed source model for the full-wave electromagnetic simulation of nonlinear terahertz generation

doi:10.1364/OE.20.018397

Distributed source model for the full-wave electromagnetic simulation of nonlinear terahertz generation

Christophe Fumeaux, Hungyen Lin, Kazunori Serita, Withawat Withayachumnankul, Thomas Kaufmann, Masayoshi Tonouchi, and Derek Abbott »View Author Affiliations

Optics Express, Vol. 20, Issue 16, pp. 18397-18414 (2012)
http://dx.doi.org/10.1364/OE.20.018397

View Full Text Article

Acrobat PDF (1726 KB)

Journal Search
Article Lookup

Article Tools

Citations

Alert me when this article is cited
Export Citation/Save

Full Text
Article Info
References (43)
Cited By
Figures (10)
Multimedia (1)
Metrics

Abstract

The process of terahertz generation through optical rectification in a nonlinear crystal is modeled using discretized equivalent current sources. The equivalent terahertz sources are distributed in the active volume and computed based on a separately modeled near-infrared pump beam. This approach can be used to define an appropriate excitation for full-wave electromagnetic numerical simulations of the generated terahertz radiation. This enables predictive modeling of the near-field interactions of the terahertz beam with micro-structured samples, e.g. in a near-field time-resolved microscopy system. The distributed source model is described in detail, and an implementation in a particular full-wave simulation tool is presented. The numerical results are then validated through a series of measurements on square apertures. The general principle can be applied to other nonlinear processes with possible implementation in any full-wave numerical electromagnetic solver.

1. Introduction

The interaction of optical beams with nonlinear crystals can be exploited for the generation of terahertz radiation, for example through the processes of optical rectification [1

A. Rice, Y. Jin, X. F. Ma, X.-C. Zhang, D. Bliss, J. Larkin, and M. Alexander, “Terahertz optical rectification from <110> zinc-blende crystals,” Appl. Phys. Lett. 64(11), 1324–1326 (1994). [CrossRef]

–3

J. Z. Xu and X.-C. Zhang, “Optical rectification in an area with a diameter comparable to or smaller than the center wavelength of terahertz radiation,” Opt. Lett. 27(12), 1067–1069 (2002). [CrossRef] [PubMed]

] or difference frequency generation (DFG) [4

J. B. Khurgin, M. W. Pruessner, T. H. Stievater, and W. S. Rabinovich, “Suspended AlGaAs waveguides for tunable difference frequency generation in mid-infrared,” Opt. Lett. 33(24), 2904–2906 (2008). [CrossRef] [PubMed]

, 5

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]

]. For both those nonlinear processes, the wavelength of the optical pump beams is commonly more than two orders of magnitude shorter than the wavelength in the band of the generated terahertz beam. In the case of optical rectification, which is the main focus of this paper and can be seen as a particular case of DFG, the short wavelength of the pump effectively allows the creation of a source with sub-wavelength dimensions for near-field terahertz microscopy [6

Q. Chen, Z. Jiang, G. X. Xu, and X.-C. Zhang, “Near-field terahertz imaging with a dynamic aperture,” Opt. Lett. 25(15), 1122–1124 (2000). [CrossRef] [PubMed]

–9

R. Lecaque, S. Grésillon, and C. Boccara, “THz emission microscopy with sub-wavelength broadband source,” Opt. Express 16(7), 4731–4738 (2008). [CrossRef] [PubMed]

Realistic models for these types of electro-optical sources are necessary to provide physical insight into the electromagnetic field distribution of the generated terahertz beam. Analytical methods describing nonlinear rectification [10

Y. S. Lee, Principles of Terahertz Science and Technology (Springer, 2008).

] provide an invaluable basis for the physical understanding of the terahertz generation process. They are however only applicable for canonical shapes as the mathematical intricacies rapidly limit the geometrical complexity that can be handled theoretically. Therefore, numerical techniques which solve Maxwell’s equations in discrete form need to be considered for general-purpose solutions in the near field, e.g. for the prediction of near-field interactions between a probing terahertz beam generated through optical rectification and micro-structured samples. In this case, the accuracy of simulated results is crucially dependent on the suitability of the terahertz source modeling.

Among the advanced numerical techniques in widespread application for full-wave electromagnetic simulations, the most prominent are the Finite-Integration Technique (FIT) [11

T. Weiland, “A discretization method for the solution of Maxwell’s equations for six-component fields,” AEU, Int. J. Electron. Commun. 31(3), 116–120 (1977).

], the Finite-Difference Time-Domain (FDTD) method [12

A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed. (Artech House 2005).

], the Finite-Element Method (FEM) [13

J.-M. Jin, The Finite Element Method in Electromagnetics, 2nd ed. (Wiley 2002).

] or the solution of Integral Equations (IE) with the Method of Moments (MoM) [14

W. C. Gibson, The Method of Moments in Electromagnetics (Chapman & Hall/CRC, 2008).

]. All those methods are implemented as the foundation to powerful commercially available electromagnetic simulation tools. Despite inherent difficulties in implementation and validation of nonlinear material models, some of those numerical techniques – applied in the time domain – are in principle capable of handling nonlinear effects [15

C. M. Dissanayake, M. Premaratne, I. D. Rukhlenko, and G. P. Agrawal, “FDTD modeling of anisotropic nonlinear optical phenomena in silicon waveguides,” Opt. Express 18(20), 21427–21448 (2010). [CrossRef] [PubMed]

, 16

I. Ahmed, E. H. Khoo, O. Kurniawan, and E. P. Li, “Modeling and simulation of active plasmonics with the FDTD method by using solid state and Lorentz–Drude dispersive model,” J. Opt. Soc. Am. B 28(3), 352–359 (2011). [CrossRef]

Most full-wave simulations require a discretizing grid (or mesh) with a cell resolution typically below a tenth of a wavelength (< λ/10) – which must remain applicable across the whole frequency spectrum of interest. For optical rectification and DFG however, a difficulty arises because the pump beam and the generated wave have wavelength scales differing by orders of magnitudes. As a consequence, the discretization of large three-dimensional (3D) problems with a mesh suitable for both the pump beam and the generated beam is not feasible with reasonable computer resources. To illustrate this point, assuming a pump beam frequency 100 times larger than the generated beam frequency, a simulation of the terahertz beam including a full-wave description of the pump beam would increase the memory requirement by a factor of 100³, i.e. by a factor of one million, as compared to a simulation of the terahertz beam alone. An efficient alternative to full-wave simulations including all involved wavelength scales is provided by co-simulations, i.e. coupled treatment of separate specific models for pump beam and terahertz radiation. This is the approach chosen in this paper where an analytically modeled infrared (IR) pump beam provides the sources for full-wave terahertz wave simulations. Terahertz discretized sources are distributed inside the crystal and introduced onto a series of closely-spaced source planes, which re-create the 3D non-linear process inside the nonlinear medium. This scheme, which is based on the assumption of a weak influence of the terahertz beam onto the pump beam, provides an accurate model of the generated terahertz radiation, while effectively bypassing the full-wave nonlinear treatment in the generation process. Similar distributed source models have been previously employed for the analysis of terahertz traveling-wave photonic mixers using transmission line equations [17

M. Neshat, D. Saeedkia, and S. Safavi-Naeini, “Semi-analytical calculation of terahertz signal generated from photocurrent radiation in traveling-wave photonic mixers,” Int. J. Infrared Millim. Waves 29(9), 809–822 (2008). [CrossRef]

In the second part of this paper, a specific implementation of the proposed spatially distributed excitation model is described for an electromagnetic simulation code based on a particular time-domain numerical technique, namely the Finite-Volume Time-Domain (FVTD) method [18

P. Bonnet, X. Ferrieres, B. L. Michielsen, P. Klotz, and J. L. Roumiguières, “Finite-volume time domain method,” in Time Domain Electromagnetics (S. M. Rao, Ed. San Diego, CA: Academic Press, 1999).

, 19

C. Fumeaux, D. Baumann, P. Leuchtmann, and R. Vahldieck, “A generalized local time-step scheme for efficient FVTD simulations in strongly inhomogeneous meshes,” IEEE Trans. Microw. Theory Tech. 52(3), 1067–1076 (2004). [CrossRef]

]. However, the underlying approach is general enough to be applied in conjunction with any other full-wave simulation technique.

As a practical validation example, the distributed source model is applied to the modeling of terahertz generation in a GaAs <100> crystal with sub-wavelength rectangular metallic apertures. The achieved agreement between simulations and experiments demonstrates the possibility of predictive numerical modeling for the generated terahertz beam, while taking into account the pump beam configuration, the velocity mismatch between the pump and the generated beam, as well as diffraction and refraction effects. Therefore, the scheme can be applied to support the optimization of near-field terahertz systems based on optical rectification or DFG.

2. Distributed Source Model

2.1 General source concept

The proposed distributed source modeling avoids the challenging and computationally expensive treatment of the entire nonlinear problem in a full-wave simulation. Instead, specific modeling approaches are used for the different wavelength scales involved in the nonlinear process: the propagation of the near-IR pump beam is efficiently computed using an analytical model or an asymptotic (optical) method, whereas the terahertz beam is modeled in full-wave electromagnetic simulations. Specifically in the present case, a pulsed Gaussian beam model obtained from paraxial approximation provides the near-IR pump transient transverse field distribution along the propagation axis, which is used to compute an equivalent source distribution for the full-wave terahertz model. To mimic the distributed nonlinear process inside the crystal, distributed sources are introduced in discrete form on closely spaced transverse planes. Each of these planes represents nonlinear generation in a thin crystal slice.

Equivalent sources in one of these transverse planes are defined similarly as introduced for the Gaussian aperture source model previously applied to optical rectification configurations, both in the far field [20

H. Lin, C. Fumeaux, B. M. Fischer, and D. Abbott, “Modelling of sub-wavelength THz sources as Gaussian apertures,” Opt. Express 18(17), 17672–17683 (2010). [CrossRef] [PubMed]

] and the near field [21

H. Lin, C. Fumeaux, B. Seam Yu Ung, and D. Abbott, “Comprehensive modeling of THz microscope with a sub-wavelength source,” Opt. Express 19(6), 5327–5338 (2011). [CrossRef] [PubMed]

]. In those two papers, the terahertz wave generation was approximated as originating from a single aperture with a Gaussian profile located inside the crystal. This was demonstrated to qualitatively describe the behavior of the generated terahertz beam. However, that approach required finding an appropriate effective location for the Gaussian aperture in the crystal. The extension to distributed sources, as proposed in the present paper, is based on the (phased) superposition of terahertz radiated field contributions originating from a large number of Gaussian apertures distributed in the pump beam path. In contrast to a single Gaussian aperture, the present model is able to describe the generation process a priori, i.e. without requiring adjustment of any free parameter. Therefore, predictive modeling of a terahertz generation system in the near and far field can be accomplished provided that both the material parameters and geometrical configuration are accurately known. The remainder of this section describes the main features of the co-simulation approach, taking as example the terahertz generation through optical rectification in a 500 μm thick GaAs crystal with <100> orientation, pumped with a pulsed near-IR beam.

2.2 Pump Beam Model

The pump beam is produced by a pulsed near-IR femtosecond fiber laser (Toptica), with a carrier free-space wavelength λ_0,pump = 1.55 μm. The major features of the pump beam, which is used as input for the terahertz distributed sources, are described in the following and schematically depicted in Fig. 1 .

Fig. 1 Illustration of the paraxial modeling of the pump beam in the GaAs crystal, including first reflection. (a) Schematic of the spatial beam profile w inside the crystal, showing an exaggerated divergence with the incident beam profile in solid blue line and the reflected beam profile in red dashed line. (b) Schematic of the transient pulse propagation showing envelope attenuation and modeled normalized pump beam power profile at three different points in time.

Pump beam profile along the propagation axis. The near-IR pulsed beam can be accurately approximated spatially by a Gaussian beam field distribution, as described by the paraxial approximation [22

B. E. A. Saleh and M. C. Teich, Fundamentals of Photonics (John Wiley & Sons, 1991).

]. In the considered system, the beam is focused on the GaAs crystal input surface with a waist (radius) of approximately w₀ = 25 μm. The propagation model of this pulsed beam in the GaAs crystal needs then to take into account the refractive index of the medium, which according to [23

E. D. Palik, ed., Handbook of Optical Constants of Solids (Academic Press, 1998)

] is equal to n_IR = 3.3737 at the frequency of operation, while neglecting dispersion of the pulse. The main effect from a denser propagation medium on the Gaussian beam is an increase in the Rayleigh range, and a decrease in the beam divergence. This reduction of the off-axis scattering in a dielectric medium can be intuitively understood by considering the shortening of the wavelength in the dielectric crystal

λ_{ε, I R} = λ_{0} / n_{I R}

. Consequently, a given beam waist, expressed in terms of the wavelength, becomes larger by a factor of

n_{I R}

. In the GaAs crystal with a thickness D = 500 μm, for the considered beam waist of 25 μm, the increase of beam radius can be estimated (using standard paraxial formulas for Gaussian beams) to be around 0.2 μm for a single pass through the crystal, and around 1.5 μm for a three-way path, 3D = 1500 μm, including two reflections. Considering that this beam radius variation is well within the waist measurement uncertainty, a constant Gaussian beam profile is assumed throughout the crystal length. It would be however straightforward to include any more significant variations of the beam radius into the model for cases where the beam divergence is more pronounced.

Pulse envelope attenuation. The near-IR pump beam is approximated as a single ultrashort impulse with a duration below 100 fs. Considering that the near-IR pulse shape is not critical to the model, it can be described as a delta-function impulse with power

P_{IR}

. Along the beam path, the pulse is attenuated by losses and nonlinear conversion in the GaAs crystal. With an extinction coefficient κ_IR = 1.33∙10⁻⁴ [23

E. D. Palik, ed., Handbook of Optical Constants of Solids (Academic Press, 1998)

] for 1.55 μm radiation in GaAs, an attenuation coefficient α_IR = 537 Np/m is calculated. This corresponds to a damping of the power by a factor of 0.765 along the 500-μm-thick GaAs crystal. The magnitude of this attenuation is confirmed experimentally through direct power measurement before and after the crystal. This independent measurement indicates an approximate power reduction by a similar factor of 0.72, after compensation of reflection losses. Considering the uncertainty associated with the experimental procedure, the attenuation coefficient derived from the material data in [23

E. D. Palik, ed., Handbook of Optical Constants of Solids (Academic Press, 1998)

] is assessed as describing reliably the pulse decay.

Pump pulse propagation velocity. The near-IR beam impulse propagates in the GaAs crystal with a velocity determined by the refractive index of the medium, i.e. c_IR = c₀/n_IR, where c₀ is the free-space light velocity. In the source model described later in Sec. 2.3, this velocity of propagation defines the propagation of the numerical source terms. This is essential to accurately model the velocity mismatch between pump beam and generated terahertz beam.

Multiple Reflections. When the pump beam pulse reaches the GaAs crystal output surface, it is partly reflected back into the crystal. The first path through the crystal builds the main contribution to terahertz generation and results in a strong initial pulse, whereas reflections of the pump pulse create delayed secondary terahertz pulses. The initial pulse can be isolated from subsequent reflections in a time-gated system. Thin crystals however require considering multiple reflections. To this end, the pump pulse model can also be extended to include propagation after reflection from the output interface to air. In the test case considered, the calculation of the Fresnel reflection coefficient at normal incidence reveals that around R = 29% of the incident power is reflected from each GaAs-Air interface. The damped and reflected pulse in turn also contributes to optical rectification in its direction of propagation and creates delayed secondary pulses. Further contributions from multiple near-IR beam reflections can be similarly taken into account in the propagation model. The number of reflections is denoted in the following as m_refl. Considering that the power is reduced by a substantial amount at each reflection, only two relevant scenarios are considered: firstly, a single pass of the IR beam through the crystal (m_refl = 0), which creates the primary generated terahertz pulse that can be isolated in measurement through time-gating. The second considered scenario involves two inner reflections (m_refl = 2), i.e. three passes through the crystal, resulting in the emission of the primary terahertz pulse and one secondary terahertz pulse from the crystal output surface.

2.3 Distributed terahertz sources

The generation of terahertz radiation through optical rectification in GaAs crystal is a distributed process discretized in the volume as a number of layered sources. The on-axis and transversal formulation is presented in the following and depicted in Fig. 2 .

Fig. 2 Schematic principle of the distributed source modeling for nonlinear rectification. (a) On-axis discretization of the nonlinear crystal in slices modeled as radiating apertures, with transverse discretization of sources schematically depicted as dots. The dots are evenly distributed for representation, but can also be arranged on an unstructured grid. (b) Equivalence of spatial source amplitude distribution to instantaneous IR power density. (c) Time-domain linear transfer function from the IR pulse to the generated terahertz pulse.

On-axis source discretization. The volume of a nonlinear crystal with a thickness of D = 500 μm is divided into infinitesimally thin slices along the beam propagation axis, where each of the slices is a source of terahertz radiation. Each individual slice can then be approximated as a two-dimensional transverse aperture, as depicted in Fig. 2(a). Each radiating aperture can be handled as shown in [21

H. Lin, C. Fumeaux, B. Seam Yu Ung, and D. Abbott, “Comprehensive modeling of THz microscope with a sub-wavelength source,” Opt. Express 19(6), 5327–5338 (2011). [CrossRef] [PubMed]

], and the contributions from all the apertures, in time domain, can be superposed as described next.

Transverse source discretization. The generated terahertz field in each aperture can be expressed according to the equivalence principle [24

C. A. Balanis, Antenna Theory, 3rd Edition (Wiley, 2005).

] as distributed electric and magnetic current sources

{\vec{J}}_{THz}

and

{\vec{M}}_{THz}

. Those equivalent sources at time t and location

\vec{r}

are related to the generated transverse components of the terahertz electric and magnetic fields (denoted as

{\vec{E}}_{THz}

and

{\vec{H}}_{THz}

) according to the well-known equations

\begin{matrix} {\vec{J}}_{THz} (t, \vec{r}) = \hat{n} \times {\vec{H}}_{THz} (t, \vec{r}) \\ {\vec{M}}_{THz} (t, \vec{r}) = - \hat{n} \times {\vec{E}}_{THz} (t, \vec{r}), \end{matrix}

(1)

where

\hat{n}

represents the unit vector in the direction of beam propagation. This notation represents a vector form of the Huygens’ principle, which in the case of a Gaussian aperture takes the form of infinitesimal crossed-dipole sources made out of superposed equivalent electric and magnetic current elements, as described in [20

H. Lin, C. Fumeaux, B. M. Fischer, and D. Abbott, “Modelling of sub-wavelength THz sources as Gaussian apertures,” Opt. Express 18(17), 17672–17683 (2010). [CrossRef] [PubMed]

] and depicted in Fig. 2(b). To mimic the square law nonlinear rectification process [25

D. Li and G. Ma, “Pump-wavelength dependence of terahertz radiation via optical rectification in (110)-oriented ZnTe crystal,” J. Appl. Phys. 103(12), 123101 (2008). [CrossRef]

], the time-dependent magnitude of the generated terahertz fields is obtained through convolution of the instantaneous pump beam power density

P_{IR}

with a linear transfer function

T (t)

, see Fig. 2(c), i.e.

\begin{matrix} E_{THz} (t, \vec{r}) = P_{IR} (t, \vec{r}) * T (t) \\ H_{THz} (t, \vec{r}) = P_{IR} (t, \vec{r}) * T (t) / Z_{ε_{r}, THz}, \end{matrix}

(2)

where

Z_{ε_{r}, THz}

is the intrinsic impedance of the crystal at terahertz frequency, * denotes the convolution, and as mentioned,

P_{IR} (t, \vec{r})

is approximated as a propagating delta function. This impulse response can be also interpreted as analogy to a moving particle with a dipole moment [26

D. H. Auston, “Subpicosecond electro-optic shock waves,” Appl. Phys. Lett. 43(8), 713–715 (1983). [CrossRef]

, 27

D. A. Kleinman and D. H. Auston, “Theory of electrooptic shock radiation in nonlinear optical media,” IEEE J. Quantum Electron. 20(8), 964–970 (1984). [CrossRef]

]. The linear relationship of Eq. (2) between pump power density and generated field amplitudes provides the crucial coupling link between the near-IR model and the full-wave terahertz simulation. The equivalent crossed-dipole source distribution Eq. (1) is then constructed on the basis of the propagating pump pulse as described in the following.

Volumetric source. The volumetric nature of the source is re-created by superposing all contributions from the apertures, each 2D aperture representing a thin 3D crystal slice. In a time-domain solver, this is simply accomplished by introducing time-dependent electric and magnetic equivalent current sources at each discrete source location. The current source amplitudes are determined by the instantaneous pump pulse power density according to Eq. (1) and Eq. (2). According to Fig. 1, we consider a crystal with a thickness D, with the input surface located at the origin. The ultrashort near-IR pump pulse with Gaussian transverse profile propagates in the positive z direction, and its waist is located at the crystal input surface. Terahertz sources are then defined by considering M_refl internal reflections of the pump beam within the crystal, which translates into M_refl delayed source contributions denoted by index m_refl = 1,…, M_refl. Therefore, the magnitudes of the electric and magnetic current sources at time t and at location

\vec{r} = (x, y, z)

can be written as a sum

\begin{matrix} {\vec{J}}_{THz} (t, x, y, z) = \sum_{m_{refl} = 0}^{M_{refl}} \underset{(a)}{\underset{︸}{η_{optrect} P_{IR,max}}} \cdot \underset{(b)}{\underset{︸}{e^{- {(t_{D} / 2 σ_{G})}^{2}} \cdot \sin (2 π f_{c} t_{D})}} \\ \cdot \underset{(c)}{\underset{︸}{e^{- α_{IR} z_{p}}}} \cdot \underset{(d)}{\underset{︸}{R^{m_{refl}}}} \cdot \underset{(e)}{\underset{︸}{e^{- 2 (x^{2} + y^{2}) / w^{2} (z_{p})}}} \cdot \underset{(f)}{\underset{︸}{(\begin{matrix} ρ_{x} \\ ρ_{y} \\ 0 \end{matrix})}} \end{matrix}

(3)

{\vec{M}}_{THz} (t, x, y, z) = \sum_{m_{refl} = 0}^{M_{refl}} \underset{(g)}{\underset{︸}{{(| {\vec{J}}_{THz} (t, \vec{r}) | \cdot Z_{ε_{r}, THz})}_{m_{refl}}}} \cdot \underset{(h)}{\underset{︸}{{(- 1)}^{m_{refl}} \cdot (\begin{matrix} - ρ_{y} \\ ρ_{x} \\ 0 \end{matrix})}},

(4)

where

z_{p}

is the accumulated (physical) path inside the crystal, as measured from the input surface and considering the number of reflections

m_{refl}

to this point

z_{p} = {\begin{matrix} z + m_{refl} \cdot D for m_{refl} even \\ (D - z) + m_{refl} \cdot D for m_{refl} odd \end{matrix},

(5)

and where

t_{D}

is the retardation time

t_{D} = t - 3 σ_{G} - \frac{z_{p}}{c_{0} / n_{IR}} .

(6)

The different terms and variables in the Eqs. (3)-(5) are explicitly described in the following:

(a) The first term defines the magnitude of the sources. It includes a proportionality factor $η_{optrect}$ that describes the efficiency of the nonlinear process, i.e. the direct link between the power of the near-IR pump beam and the terahertz source magnitude. In the present notation, the maximum power density in the centre of the beam at the waist is denoted as $P_{IR,max}$ . Under practical measurement settings, the terahertz signal after probing a sample is normalized using a baseline measurement, i.e. a reference signal without sample. This normalization approach allows canceling of this term (a), and is also used in the simulations presented in this paper.
(b) The second factor models the transient behavior of the broadband terahertz pulse, which basically describes the transfer function $T (t)$ from the ultrashort IR pulse to the generated terahertz pulse. In the presented case, this time-dependence is chosen as a sine-modulated Gaussian pulse with bandwidth $1 / π σ_{G}$ , where $σ_{G}$ is the temporal standard deviation of the generated pulse. The spectrum of the pulse is centered at the frequency $f_{c}$ and its standard deviation $σ_{G}$ is selected to cover the 0.1–2.5 THz frequency band. It should be emphasized that the transient shape of this pulse may, but does not have to, correspond to the actual shape of the generated terahertz pulse. The fundamental requirement is that the simulated pulse bandwidth needs to cover the terahertz spectrum of interest, and so provides sufficient signal to (numerical) noise ratio for Fourier analysis.
(c) The third factor describes the attenuation of the beam power density along the beam path (Sec. 2.1), normalized with respect to the maximum input power density $P_{IR,max}$ . It is emphasized that the attenuation term has to take into account the accumulated path $z_{p}$ inside the crystal, i.e. considering also the path after reflections.
(d) The reflections at the crystal inner surfaces after the first path of the pump pulse are taken into account using the Fresnel power reflection coefficient R, raised to a power equal to the number of reflections m_refl. For the initial path of the incident pulse through the crystal, m_refl equals to 0.
(e) The Gaussian transverse field distribution is described by the fifth factor, which is obtained based on the pump beam power profile with radius $w (z_{p})$ . Without a loss of generality, this radius is approximated as a constant $w (z_{p}) = w_{0}$ in our example (Sec. 2.1). It is noted that the square law relationship introduces the factor of 2 in the exponent.
(f) The last term of Eq. (3) is the polarization vector $\vec{ρ} = {(ρ_{x}, ρ_{y}, 0)}^{T}$ of the generated terahertz pulse. For propagation in the + z direction, this unit vector indicates a direction of the generated terahertz electric field in the xy-plane. It is clear from Eq. (1) that the polarization is along the same direction as the equivalent electric currents. In the example shown later, a polarization $\vec{ρ} = {(1, 0, 0)}^{T}$ is selected.
(g) In Eq. (4), the first term describes the proportionality between the amplitudes of the magnetic and electric current sources. In the crossed-dipole model that excites the Gaussian aperture [28
C. Fumeaux, D. Baumann, S. Atakaramians, and E. Li, “Considerations on paraxial Gaussian beam source conditions for time-domain full-wave simulations,” Annual Rev. of Progress in Appl. Computational Electromagnetics, 401–406 (2009).
], this proportionality is determined by the intrinsic impedance $Z_{ε_{r}, THz}$ of the medium (GaAs) at terahertz frequencies. Therefore, the summation terms for the magnetic current sources are the same terms (a)-(e) as in Eq. (4) multiplied by $Z_{ε_{r}, THz}$ .
(h) The vector term in Eq. (4) defines the direction of the generated magnetic field in the transverse plane, and is therefore orthogonal to the polarization vector $\vec{ρ}$ (f). A + sign in front the vector indicates sources launching a wave propagating along the beam axis (towards + z, i.e. for the initial pulse and for the third pass), whereas the – sign results in a generated wave launched in the opposite direction (towards –z, e.g. after first inner reflection).

The retardation time

t_{D}

given in Eq. (6) is composed of the reference time t and two additional terms contributing to retardation. The first retardation term ensures that the sine-modulated Gaussian pulse is launched with a negligible initial value in the whole computational domain at

t = 0

. A value of this retardation term corresponding to 3 to 4 times

σ_{G}

is generally deemed as acceptable. The second retardation time is the spatial retardation, which describes the time that the pump near-IR pulse needs to travel from the crystal input surface at

z = 0

to the source location after an accumulated path

z_{p}

. This retardation time naturally introduces the propagation velocity of the pump beam within the crystal into the discretized terahertz sources.

2.4 Full-wave time-domain modeling of the terahertz beam

The distributed discrete sources introduced in the previous paragraph can be straightforwardly implemented in any time-domain full-wave simulation methods. This can be for example realized as excitation currents according to Eq. (3) and Eq. (4). According to Eq. (1), an equivalent and alternative excitation method consists of adding generated fields to the locally computed field, as soft sources. The particularity of the numerical simulation technique will decide on the specific implementation method.

The particular implementation presented in this paper is realized with an in-house code based on the FVTD method [18

], which is briefly described here for the sake of completeness. FVTD has been introduced in the early 1990’s [29

N. K. Madsen and R. W. Ziolkowski, “A three-dimensional modified finite volume technique for Maxwell’s equations,” Electromagnetics 10(1-2), 147–161 (1990). [CrossRef]

, 30

V. Shankar, A. H. Mohammadian, and W. F. Hall, “A time-domain, finite-volume treatment for the Maxwell equations,” Electromagnetics 10(1-2), 127–145 (1990). [CrossRef]

] as a full-wave time-domain method implemented in unstructured meshes, e.g. in a tetrahedral mesh. This feature makes the method particularly well suited for conformal and multi-scale problems, considering that a tetrahedral discretization allows adapting locally the mesh density to the spatial resolution required by the geometry. To improve the efficiency of the simulation, the spatial mesh inhomogeneity can be combined with an inhomogeneity in temporal discretization realized as a local time stepping scheme [19

A schematic of the discretized problem (realized using Altair HyperMesh) is shown in Fig. 3 . The longitudinal distribution of sources along the beam path is introduced in the FVTD mesh as auxiliary source surfaces, each discretized with a triangular mesh. On each of these apertures, sources terms are implemented in the FVTD algorithm as extra fluxes, added to the standard FVTD fluxes originating from the solution of Maxwell’s equations. The aperture radius (200 μm) is much larger than the pump beam radius. For the distribution of sources along the 500 μm thick crystal, 50 apertures with inter-plane distance of 10 μm are selected (Fig. 3(b)). Simulation results are not critically dependent on this spacing, as long as it is chosen as a fraction of the shortest wavelength of interest, in the present case

\sim λ_{ε, \min} / 3

. Reduction of the computational domain volume is possible provided there are symmetries in the geometry of the pump beam and probed structures. In the present case, as shown in Fig. 3, both electric and magnetic wall symmetries are used to reduce the problem size by a factor of four. A spherical outer surface with absorbing boundary condition (ABC) realized as quasi-perfectly-matched absorber is used for nearly reflection-free domain truncation [31

C. Fumeaux, K. Sankaran, and R. Vahldieck, “Spherical perfectly matched absorber for finite-volume 3-D domain truncation,” IEEE Trans. Microw. Theory Tech. 55(12), 2773–2781 (2007). [CrossRef]

Fig. 3 Full-wave simulation set-up. (a) General view of interfaces in the FVTD model, including spherical absorbing boundary conditions (ABC), as well as perfect electric conductor (PEC) and perfect magnetic conductor (PMC) symmetries. The materials are indicated in bold fonts. (b) Partial view of one discretized source aperture inside the crystal, the outline of other source planes being shown as thick lines. The tetrahedral volume discretization is not shown for the sake of clarity. (c) Schematic of one triangular element on a source aperture, with indication of the typical dimension, as well as the orientation of normal vector and current source vectors.

The relative permittivity of GaAs at terahertz frequencies is taken as

ε_{r, THz} =

12.96 [23

E. D. Palik, ed., Handbook of Optical Constants of Solids (Academic Press, 1998)

]. Dispersion is neglected in the present case because of the stability of the values in the range of interest, but dispersive materials could be handled as described in [32

D. Baumann, C. Fumeaux, C. Hafner, and E. P. Li, “A modular implementation of dispersive materials for time-domain simulations with application to gold nanospheres at optical frequencies,” Opt. Express 17(17), 15186–15200 (2009). [CrossRef] [PubMed]

]. The relatively high value of permittivity requires the discretization inside the crystal to be very fine, as cell dimensions have to be a small fraction of the shortest dielectric wavelength of interest

λ_{ε, \min}

in the terahertz range. Therefore, a discretization of about 2 μm (

λ_{ε, \min} / 16

) is chosen in the present case in the crystal around the axis of the generated beam, as illustrated in Fig. 3(c). Considering that the dielectric length of the crystal alone represents 15

λ_{ε, \min}

, the problem becomes computationally expensive. To reduce the computational cost, the discretization is progressively increased to around 10 μm for the free-space cells outside the crystal. Under such configuration, around 10 million tetrahedral cells that translate into 7 GB of computer memory (RAM) are required.

To retrieve the spectral behavior of the terahertz pulse, the transformation of time-domain data to the frequency domain is obtained through a Discrete Fourier Transform performed on the fly during the FVTD time iteration [12

A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed. (Artech House 2005).

]. The frequency-domain data is used for the modeling of the detection process, as described in the following.

2.5 Far-field detection

In the considered optical rectification system, the transmitted terahertz radiation after near-field interaction is detected in the far-field region with a photoconductive antenna. Only a part of the diffracted radiation is collected by the first parabolic mirror of a standard far-field terahertz time-domain spectroscopy (TDS) detection setup and reaches the detector after being guided along an optical train.

In the present simulation, the whole far-field detection process is modeled as depicted in Fig. 4 . It starts by performing a near-field to far-field transformation [12

A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed. (Artech House 2005).

], using an infinite Huygens’ plane located 100 μm after the output crystal surface. This provides far-field radiation patterns describing the angular dependence of the radiated field amplitude and phase, for all frequencies of interest. Integrating the patterns coherently over the acceptance solid angle covered by the first parabolic mirror (or lens) yields the detected signal.

Fig. 4 Schematic of far-field detection process, (a) far-field radiation pattern and acceptance angle by the first parabolic mirror or lens, (b) example of far-field amplitude pattern at 2 THz for the case considered, and illustration of the acceptance angle inside which the detected signal is integrated.

2.6 Present model limitations and possible extensions

The distributed source approach provides a simple modeling technique for the simulation of nonlinear optical rectification. The main simplifying assumption in this implementation is that the effect of the generated terahertz beam on the pump beam is negligible. In other words, the generated beam model does not influence the pump beam model. The typical conversion efficiencies (< 10⁻⁴) in this type of systems [33

M. C. Hoffmann, K.-L. Yeh, J. Hebling, and K. A. Nelson, “Efficient terahertz generation by optical rectification at 1035 nm,” Opt. Express 15(18), 11706–11713 (2007). [CrossRef] [PubMed]

] validates this assumption. In the following, some other present limitations are listed, with possible extensions aiming at increasing the sophistication of the physical description and refining accuracy of the approach.

Available data. The accuracy of a computer-generated solution is in principle fundamentally dependent on the numerical method accuracy. The numerical methods available for electromagnetic simulations are very mature and are thoroughly validated. So using best practices, the confidence in the accuracy of the solution for a standard problem is usually very high. However, another major aspect that influences the agreement of a numerical solution with experimental results is the reliability of the available geometry data and material characteristics. For problems in the terahertz range, the fabrication tolerances are often not negligible compared to the wavelength of operation and also the materials electromagnetic properties are not known precisely. Examples of uncertainties include the thickness of the crystal, which is only specified to a few microns accuracy, and the electromagnetic properties of the crystal. In this sense, some assumptions that are made in the presented simulations can be refined through better knowledge of the actual problem. As an example, the dispersion of GaAs properties is neglected in the simulations presented below, due to a greater material properties uncertainty compared to dispersive effects. An independent spectral measurement of these properties on the actual crystal sample could assist in refining the model.

Pump beam model. A near-IR ultrashort pulse model with Gaussian spatial distribution is assumed in the pump beam model, as it describes the beam accurately. There are, however, no inherent limitation to the generality of the source model that would preclude extension to other field distributions. For example, a propagating IR beam computed using a semi-analytical model [34

E. K. Rahani and T. Kundu, “Electromagnetic THz radiation modeling by DPSM,” Int. J. Infrared Millim. Waves 33(3), 376–390 (2012). [CrossRef]

] could be employed to define the distributed sources.

Conversion efficiency. In the present simulations, the relative field distribution of the generated beam is obtained a priori. Prior knowledge in the conversion efficiency term

η_{optrect}

is bypassed by normalizing a sample with a reference measurement. A physical conversion efficiency, obtained through an analytical model [35

D. Côté, J. E. Sipe, and H. M. van Driel, “Simple method for calculating the propagation of terahertz radiation in experimental geometries,” J. Opt. Soc. Am. B 20(6), 1374–1385 (2003). [CrossRef]

], could be included in the model to extend the simulation towards absolute results.

Application in frequency-domain full-wave methods. Equations presented in Sec. 2.3 are tailored for the source definition in time-domain full-wave methods, as these are conceptually close to the pulse-based experimental procedure. It is, however, straightforward to adapt these equations for implementation in frequency-domain full-wave methods, by translating the time delay arising from the propagation of the pump pulse into a phase shift.

Application in commercial codes. Application of the presented source technique is applicable in commercial codes, provided the possibility of inputting custom-defined sources on multiple planes. As non-exhaustive examples, the implementation in CST Microwave Studio and SEMCAD X appears feasible within the framework of available source definitions. Specific implementations within those codes will be subject of further communications.

Pre-computation of source. For a particular system, the generated beam can be in principle pre-computed to yield a field distribution at the output surface (either in the time domain or in the frequency domain). This pre-computation of the source can be realized partly analytically or numerically and can provide aggregated information about the distributed process. This can potentially lower the computational cost significantly, especially for simulations of different samples with a given system.

3. Experimental Validation

3.1 Test setup with sub-wavelength apertures

Characterizing the generated terahertz beam in the near-field region is very challenging in an experimental setup, as a near-field probe can significantly affect the field distribution [36

K. Wang, D. M. Mittleman, N. C. J. van der Valk, and P. C. M. Planken, “Antenna effects in terahertz apertureless near-field optical microscopy,” Appl. Phys. Lett. 85(14), 2715–2717 (2004). [CrossRef]

]. Therefore, an indirect method is chosen to probe the beam through a standard far-field detection procedure. For this purpose, the GaAs crystal output surface is metallized with a 100 nm gold layer. This thickness is sufficient to block all terahertz radiation, but can still be considered infinitesimal in simulation. Small square apertures are then etched into this metal layer [37

A. J. L. Adam, J. M. Brok, M. A. Seo, K. J. Ahn, D. S. Kim, J. H. Kang, Q. H. Park, M. Nagel, and P. C. Planken, “Advanced terahertz electric near-field measurements at sub-wavelength diameter metallic apertures,” Opt. Express 16(10), 7407–7417 (2008). [CrossRef] [PubMed]

, 38

M. A. Seo, A. J. L. Adam, J. H. Kang, J. W. Lee, K. J. Ahn, Q. H. Park, P. C. M. Planken, and D. S. Kim, “Near field imaging of terahertz focusing onto rectangular apertures,” Opt. Express 16(25), 20484–20489 (2008). [CrossRef] [PubMed]

], with three different sizes as shown in Fig. 5(a) . Such sub-wavelength aperture structures are currently within the scope of many investigations, especially in conjunction with extraordinary transmission phenomena [39

A. Bitzer and M. Walther, “Terahertz near-field imaging of metallic subwavelength holes and hole arrays,” Appl. Phys. Lett. 92(23), 231101 (2008). [CrossRef]

]. In the simulation, a single discretized crystal model includes the geometry of the three squares on the crystal output surface, with the possibility to change the material property from metallic (approximated here as perfectly conducting layer) to a transparent GaAs-to-air material boundary, Fig. 5(b). The three apertures are then evaluated in separate simulations. A rounding of the square aperture corners is included in the geometry, with a curvature radius of 7 μm estimated from inspection of the fabricated samples.

Fig. 5 Probing apertures on a metallized GaAs crystal. (a) Photographs of manufactured apertures with three different sizes. (b) Triangular surface discretization of the FVTD model including the three possible apertures. The size of the aperture can be selected by selecting the surface properties either as metal or transparent GaAs-Air interface.

In the measurements, the pump beam is focused onto the input crystal surface with a beam waist of 25 μm and the generated terahertz field is measured with a photoconductive antenna in a far-field setup. System details can be found in [40

K. Serita, S. Mizuno, H. Murakami, I. Kawayama, Y. Takahashi, M. Yoshimura, Y. Mori, J. Darmo, and M. Tonouchi, “Scanning laser terahertz near-field imaging system,” Opt. Express 20(12), 12959–12965 (2012). [CrossRef] [PubMed]

]. As mentioned, only the radiated field collected by the first parabolic mirror (or lens) reaches the detector, Fig. 4(a). The frequency dependence of the terahertz wave form is obtained through Fourier transform of the measured pulse, using standard TDS techniques with appropriate time-gating to either include or suppress the secondary (reflected) pulses. The measurements performed on the GaAs with square apertures are normalized to a reference measurement obtained from a bare GaAs crystal (i.e. without metallization). As described in Sec. 2, simulation of the system takes into account all relevant experimental aspects, from the pump beam to the far-field detection and normalization process. It should be emphasized that simulations using plane waves, as commonly available in commercial codes and applied for this type of problems, would not provide the information needed for normalization.

3.2 Time-domain results

A representation of the simulated time evolution of the generated terahertz pulse is given in Fig. 6 , for the near-IR pulse with Gaussian distribution and waist w₀ = 25 μm. The first three snapshots starting from the left-hand side in Fig. 6(a) show the terahertz pulse propagating in the crystal, being amplified by the distributed sources along the pump beam path, while being subject to diffraction effects arising from the sub-wavelength transverse dimensions of the source distribution. The snapshot at 5.88 ps shows the terahertz pulse reaching the crystal output surface. The last snapshot shows the primary radiated pulse, as well as the part of the pulse reflected in the crystal, travelling in opposite directions.

Fig. 6 (a) Instantaneous images of the E-field magnitude at different times, showing the terahertz pulse amplification and propagation through the nonlinear crystal. Only one half of the symmetric field distribution (for x > 0) is shown, with a PEC symmetry at x = 0. The excitation pulse is a sine-modulated Gaussian pulse as described in Sec. 2.3. (Media 1). (b) Zoom on the (rescaled) simulated THz pulse in the GaAs illustrating the Čherenkov cone with an angle close to 20.42 degrees, arising from the velocity mismatch between the IR and terahertz beam (c) Simulated pulse detected at a sensor located 400 μm in front of the aperture.

From the field distributions in Fig. 6(a) and the online movie, the Čherenkov cone of the radiating terahertz wave can clearly be observed. The effect is ascribed to the difference in the phase velocity of terahertz radiation and the group velocity of the optical pulse. This phenomenon, illustrated by the proposed discrete model, is qualitatively agreeable with the theoretical prediction [26

D. H. Auston, “Subpicosecond electro-optic shock waves,” Appl. Phys. Lett. 43(8), 713–715 (1983). [CrossRef]

, 27

D. A. Kleinman and D. H. Auston, “Theory of electrooptic shock radiation in nonlinear optical media,” IEEE J. Quantum Electron. 20(8), 964–970 (1984). [CrossRef]

] and experimental observation [41

D. H. Auston, K. P. Cheung, J. A. Valdmanis, and D. A. Kleinman, “Čherenkov radiation from sfemtosecond optical pulses in electro-optic media,” Phys. Rev. Lett. 53(16), 1555–1558 (1984). [CrossRef]

]. Further investigation, as shown in the snapshot of the pulse in Fig. 6(b), shows that the angle

α_{\overset{⌣}{C}}

of the Čherenkov cone apparent in the numerical model is very close to the theoretical value of 20.42 degrees determined from the permittivities

ε_{r, IR} =

11.38 and

ε_{r, THz} =

12.96 according to

α_{\overset{⌣}{C}} = arc \cos (\frac{\sqrt{ε_{r, IR}}}{\sqrt{ε_{r, THz}}}) .

(7)

The simulated time-domain shape of the generated terahertz radiation at a distance of 400 μm in front of the crystal is shown in Fig. 6(c). It is noted that a direct comparison of the time-domain results for the simulated and measured pulses is not possible, as simulated time-domain data are only available in the near field, whereas measured time-domain data are detected in the far field. A time-domain near-field to far-field transformation, as described in the Sec. 8.6 of [12

A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed. (Artech House 2005).

], could be in principle implemented in the simulations to resolve the radiated pulse shape for any selected direction in the far field. However, as mentioned in Sec. 2.3, term (b), it should be noted that, for frequency-domain analysis, the pulse shape and spectrum do not need to be identical in simulation and measurements. The shape of the curve in Fig. 6(c) illustrates the generated initial terahertz pulse as well as the first secondary delayed pulse.

3.3 Frequency-domain results

The frequency-domain results are obtained after Fourier transformation of the measured pulses. They describe the combined spectral behavior arising from the source and from the transmission through the aperture. These spectral results are therefore representative of the typical operation of a terahertz TDS system in transmission mode. The experimental system for far-field detection is modeled in the simulation as described in Sec. 2.5. In both experiment and simulation, the signal from a terahertz pulse generated from a bare GaAs crystal is first measured to provide a reference spectrum. Subsequently, the outputs from the metallized GaAs crystal with various square apertures can be normalized to this reference spectrum. The normalized results for the three apertures with edge lengths of 25, 50 and 75 μm are compared in Fig. 7 .

Fig. 7 Comparison of measured and simulated normalized transmission E-field spectrum for the three apertures considered (25, 50 and 75 μm) when the pump beam is focused on the crystal input surface with a waist of 25 μm. (a) Spectrum including the initial terahertz pulse only. (b) Spectrum including the first secondary pulse.

In Fig. 7(a), the Fourier-transformed results are shown for a time gating of ~12 ps, see Fig. 6(c), so that only the initial pulse is included in the Fourier transform. A satisfying agreement is observed between the measured data and the simulated results. Discrepancies can be observed at low and high frequencies:

• At low frequencies, increase in measured transmission below 0.1-0.3 THz is explained by the low-frequency breakdown, which causes the measured signal to decrease below the system noise dynamic range [42
P. U. Jepsen and B. M. Fischer, “Dynamic range in terahertz time-domain transmission and reflection spectroscopy,” Opt. Lett. 30(1), 29–31 (2005). [CrossRef] [PubMed]
]. This breakdown happens at a higher frequency for lower signals, i.e. for the smaller apertures. Eventually all measured curves converge to a value of 1 at dc, indicating that both the stronger reference and the aperture signals have reached the noise floor.
• At high frequencies, the measured curves exhibit a lower relative transmission compared to the simulation. There are two reasons for this discrepancy: firstly in simulation, the results above 2 GHz become sensitive to the value of the GaAs permittivity (taken here from [23
E. D. Palik, ed., Handbook of Optical Constants of Solids (Academic Press, 1998)
]). Secondly, the measurements are very sensitive to alignment, as discussed later in this section.

An interesting feature of the simulated curves is that the normalized curves can become larger than unity. This can be understood physically by the fact that only a part of the angular radiation pattern is integrated for detection (Fig. 4(b)). Especially at higher frequencies, the aperture can block out-of-phase contributions which would destructively interfere with broadside radiation. This can increase the normalized received field magnitude past 1.0.

The comparison between the measurements and simulations is shown in Fig. 7(b) for the case where the time-gating is extended to include the secondary pulse (~24 ps). The oscillations of the signal arising from interferences of the primary and secondary pulses are reproduced qualitatively in the simulations. There is a discrepancy in the amplitude of these oscillations, but our experience indicates that these oscillations are very sensitive to alignment during measurement. Figure 8 compares the simulated results for the 75 μm aperture to those obtained in two independent measurements. Figure 8(a) shows the same results as the black curve in Fig. 7(b), whereas Fig. 8(b) shows the same results obtained with a different alignment. Figure 8(b) shows an almost perfect match qualitatively and quantitatively, both in terms of oscillations and general behavior of the curve. It is noticed that this confirms experimentally the possibility of measuring values larger than unity above 2 THz. The synchronicity of the oscillations is not conserved between 0.5 and 1.0 THz, which is attributed to the fact that the simulation neglects dispersion of the permittivity of the GaAs in the terahertz range. For consistency, the set of measurements corresponding to Fig. 8(a) is chosen as primary set for this paper, as it is complete.

Fig. 8 Normalized transmission E-field spectrum in two independent measurements, compared to the simulation. The curves (a) correspond to the top curves of Fig. 7(b).

As a further validation, the focus point of the pump beam is displaced from the input surface. To this end, the crystal is moved away from the source by a distance calculated to increase the beam radius on the input surface to 50 μm. The comparison of simulations and measurements in this case is shown in Fig. 9 , and shows a satisfying agreement for both cases of time gating considered.

Fig. 9 Comparison of measured and simulated normalized transmission E-field spectrum for the three apertures considered (25, 50 and 75 μm) for the pump beam defocused, with a beam radius of around 50 μm at the input surface. (a) Spectrum including the initial terahertz pulse only. (b) Spectrum including the first secondary pulse.

4. Application to optical rectification system design

The knowledge gained from the full-wave simulation of optical generation system can be very helpful for supporting the design and analysis of an optical rectification system, by providing information on the field behavior within the crystal and in the immediate near-field region. It is not possible to attain such information from conventional simulation models where only a point source or plane wave excitation is available. Figure 10(a) illustrates the growth of the terahertz field in the crystal for two selected frequencies. It can be immediately inferred from the field distribution that at some distance in the crystal, the field reaches saturation, suggesting the upper bound of the crystal thickness. The time-resolved terahertz field pattern inside the crystal might suggest optimal approaches to coupling the beam into free space.

Fig. 10 Illustrative results from the source simulations. (a) Amplitude growth in the crystal. In those graphs, the velocity mismatch between pump and terahertz radiation is the same as in the previous validation example. The standing wave pattern is due to reflections. (b) Field amplitude at 2 THz for different distances (0.5, 25 and 50 μm) from the output surface. There is a localized field enhancement in the small aperture, and therefore, the scale for the top right plot is capped to about 1/3 of the achieved value.

The proposed model also provides data about the frequency-dependent divergence of the beam away from the surface of the crystal. As an example, the field amplitude at different distances from the crystal output surface is shown in Fig. 10(b) at a selected frequency. Left-hand side plots illustrate the rapid widening of the beam with increasing distance due to diffraction. This particular information on field distribution assists in effective manipulation of the generated beam in the near field. It can help in designing optical components to collect and direct the terahertz power to the target. The right-hand side pictures show for comparison the same beam after transmission through the 25 μm aperture. A field enhancement localized to the immediate vicinity of the aperture is clearly observed. It is possible to model other sophisticated structures, such as metamaterials or plasmonic hole arrays, and predict their responses due to the non-uniform Gaussian beam excitation. In fact, the technique has already been successfully applied to the prediction of frequency-dependent beam size in the design of a metamaterial-based thin-film sensor excited by a near-field source [43

W. Withayachumnankul, H. Lin, K. Serita, C. M. Shah, S. Sriram, M. Bhaskaran, M. Tonouchi, C. Fumeaux, and D. Abbott, “Sub-diffraction thin-film sensing with planar terahertz metamaterials,” Opt. Express 20(3), 3345–3352 (2012). [CrossRef] [PubMed]

5. Conclusion

This paper has demonstrated a three-dimensional source technique for full-wave simulation of terahertz radiation generated through optical rectification in a nonlinear crystal. The approach has been described in detail and demonstrated in a particular implementation. This work is based on the Hugyens’ principle and the superposition principle, and therefore it is general enough for straightforward implementation in other numerical techniques. Experimental validation through TDS terahertz measurement of sub-wavelength apertures with various sizes has demonstrated the suitability of the technique. An important feature of this model is that it allows predictive modeling of the frequency-dependent field distribution for the generated terahertz radiation, based on knowledge of the geometry and materials of the system. Beyond the particular case of optical rectification, the same principle can be refined and also extended to other nonlinear processes such as DFG. It is anticipated that similar approaches, integrated into commercial full-wave simulation codes, can become very helpful in the design of near-field terahertz system based on nonlinear crystal generation.

Acknowledgment

CF acknowledges the ARC Future Fellowship funding scheme under FT100100585. HL acknowledges the Prime Minister’s Australia Asia Endeavour Award for research collaboration with Japan. WW acknowledges Australian Post-Doctoral Fellowship from the Australian Research Council (ARC) through Discovery Project DP1095151. MT acknowledges funding by Grant-in-Aid A22246043, JSPS, and Industry-Academia Collaborative R&D, JST. DA acknowledges funding via the Australian Research Council (ARC) Discovery Projects DP0881445 and DP1097281.

References and links

1.	A. Rice, Y. Jin, X. F. Ma, X.-C. Zhang, D. Bliss, J. Larkin, and M. Alexander, “Terahertz optical rectification from <110> zinc-blende crystals,” Appl. Phys. Lett. 64(11), 1324–1326 (1994). [CrossRef]
2.	G. Dakovski, B. Kubera, and J. Shan, “Localized terahertz generation via optical rectification in ZnTe,” J. Opt. Soc. Am. B 22(8), 1667–1670 (2005). [CrossRef]
3.	J. Z. Xu and X.-C. Zhang, “Optical rectification in an area with a diameter comparable to or smaller than the center wavelength of terahertz radiation,” Opt. Lett. 27(12), 1067–1069 (2002). [CrossRef] [PubMed]
4.	J. B. Khurgin, M. W. Pruessner, T. H. Stievater, and W. S. Rabinovich, “Suspended AlGaAs waveguides for tunable difference frequency generation in mid-infrared,” Opt. Lett. 33(24), 2904–2906 (2008). [CrossRef] [PubMed]
5.	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]
6.	Q. Chen, Z. Jiang, G. X. Xu, and X.-C. Zhang, “Near-field terahertz imaging with a dynamic aperture,” Opt. Lett. 25(15), 1122–1124 (2000). [CrossRef] [PubMed]
7.	T. Kiwa, M. Tonouchi, M. Yamashita, and K. Kawase, “Laser terahertz-emission microscope for inspecting electrical faults in integrated circuits,” Opt. Lett. 28(21), 2058–2060 (2003). [CrossRef] [PubMed]
8.	W. Withayachumnankul, G. M. Png, X. Yin, S. Atakaramians, I. Jones, H. Lin, B. S. Y. Ung, J. Balakrishnan, B. W.-H. Ng, B. Ferguson, S. P. Mickan, B. M. Fischer, and D. Abbott, “T-ray sensing and imaging,” Proc. IEEE 95(8), 1528–1558 (2007). [CrossRef]
9.	R. Lecaque, S. Grésillon, and C. Boccara, “THz emission microscopy with sub-wavelength broadband source,” Opt. Express 16(7), 4731–4738 (2008). [CrossRef] [PubMed]
10.	Y. S. Lee, Principles of Terahertz Science and Technology (Springer, 2008).
11.	T. Weiland, “A discretization method for the solution of Maxwell’s equations for six-component fields,” AEU, Int. J. Electron. Commun. 31(3), 116–120 (1977).
12.	A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed. (Artech House 2005).
13.	J.-M. Jin, The Finite Element Method in Electromagnetics, 2nd ed. (Wiley 2002).
14.	W. C. Gibson, The Method of Moments in Electromagnetics (Chapman & Hall/CRC, 2008).
15.	C. M. Dissanayake, M. Premaratne, I. D. Rukhlenko, and G. P. Agrawal, “FDTD modeling of anisotropic nonlinear optical phenomena in silicon waveguides,” Opt. Express 18(20), 21427–21448 (2010). [CrossRef] [PubMed]
16.	I. Ahmed, E. H. Khoo, O. Kurniawan, and E. P. Li, “Modeling and simulation of active plasmonics with the FDTD method by using solid state and Lorentz–Drude dispersive model,” J. Opt. Soc. Am. B 28(3), 352–359 (2011). [CrossRef]
17.	M. Neshat, D. Saeedkia, and S. Safavi-Naeini, “Semi-analytical calculation of terahertz signal generated from photocurrent radiation in traveling-wave photonic mixers,” Int. J. Infrared Millim. Waves 29(9), 809–822 (2008). [CrossRef]
18.	P. Bonnet, X. Ferrieres, B. L. Michielsen, P. Klotz, and J. L. Roumiguières, “Finite-volume time domain method,” in Time Domain Electromagnetics (S. M. Rao, Ed. San Diego, CA: Academic Press, 1999).
19.	C. Fumeaux, D. Baumann, P. Leuchtmann, and R. Vahldieck, “A generalized local time-step scheme for efficient FVTD simulations in strongly inhomogeneous meshes,” IEEE Trans. Microw. Theory Tech. 52(3), 1067–1076 (2004). [CrossRef]
20.	H. Lin, C. Fumeaux, B. M. Fischer, and D. Abbott, “Modelling of sub-wavelength THz sources as Gaussian apertures,” Opt. Express 18(17), 17672–17683 (2010). [CrossRef] [PubMed]
21.	H. Lin, C. Fumeaux, B. Seam Yu Ung, and D. Abbott, “Comprehensive modeling of THz microscope with a sub-wavelength source,” Opt. Express 19(6), 5327–5338 (2011). [CrossRef] [PubMed]
22.	B. E. A. Saleh and M. C. Teich, Fundamentals of Photonics (John Wiley & Sons, 1991).
23.	E. D. Palik, ed., Handbook of Optical Constants of Solids (Academic Press, 1998)
24.	C. A. Balanis, Antenna Theory, 3rd Edition (Wiley, 2005).
25.	D. Li and G. Ma, “Pump-wavelength dependence of terahertz radiation via optical rectification in (110)-oriented ZnTe crystal,” J. Appl. Phys. 103(12), 123101 (2008). [CrossRef]
26.	D. H. Auston, “Subpicosecond electro-optic shock waves,” Appl. Phys. Lett. 43(8), 713–715 (1983). [CrossRef]
27.	D. A. Kleinman and D. H. Auston, “Theory of electrooptic shock radiation in nonlinear optical media,” IEEE J. Quantum Electron. 20(8), 964–970 (1984). [CrossRef]
28.	C. Fumeaux, D. Baumann, S. Atakaramians, and E. Li, “Considerations on paraxial Gaussian beam source conditions for time-domain full-wave simulations,” Annual Rev. of Progress in Appl. Computational Electromagnetics, 401–406 (2009).
29.	N. K. Madsen and R. W. Ziolkowski, “A three-dimensional modified finite volume technique for Maxwell’s equations,” Electromagnetics 10(1-2), 147–161 (1990). [CrossRef]
30.	V. Shankar, A. H. Mohammadian, and W. F. Hall, “A time-domain, finite-volume treatment for the Maxwell equations,” Electromagnetics 10(1-2), 127–145 (1990). [CrossRef]
31.	C. Fumeaux, K. Sankaran, and R. Vahldieck, “Spherical perfectly matched absorber for finite-volume 3-D domain truncation,” IEEE Trans. Microw. Theory Tech. 55(12), 2773–2781 (2007). [CrossRef]
32.	D. Baumann, C. Fumeaux, C. Hafner, and E. P. Li, “A modular implementation of dispersive materials for time-domain simulations with application to gold nanospheres at optical frequencies,” Opt. Express 17(17), 15186–15200 (2009). [CrossRef] [PubMed]
33.	M. C. Hoffmann, K.-L. Yeh, J. Hebling, and K. A. Nelson, “Efficient terahertz generation by optical rectification at 1035 nm,” Opt. Express 15(18), 11706–11713 (2007). [CrossRef] [PubMed]
34.	E. K. Rahani and T. Kundu, “Electromagnetic THz radiation modeling by DPSM,” Int. J. Infrared Millim. Waves 33(3), 376–390 (2012). [CrossRef]
35.	D. Côté, J. E. Sipe, and H. M. van Driel, “Simple method for calculating the propagation of terahertz radiation in experimental geometries,” J. Opt. Soc. Am. B 20(6), 1374–1385 (2003). [CrossRef]
36.	K. Wang, D. M. Mittleman, N. C. J. van der Valk, and P. C. M. Planken, “Antenna effects in terahertz apertureless near-field optical microscopy,” Appl. Phys. Lett. 85(14), 2715–2717 (2004). [CrossRef]
37.	A. J. L. Adam, J. M. Brok, M. A. Seo, K. J. Ahn, D. S. Kim, J. H. Kang, Q. H. Park, M. Nagel, and P. C. Planken, “Advanced terahertz electric near-field measurements at sub-wavelength diameter metallic apertures,” Opt. Express 16(10), 7407–7417 (2008). [CrossRef] [PubMed]
38.	M. A. Seo, A. J. L. Adam, J. H. Kang, J. W. Lee, K. J. Ahn, Q. H. Park, P. C. M. Planken, and D. S. Kim, “Near field imaging of terahertz focusing onto rectangular apertures,” Opt. Express 16(25), 20484–20489 (2008). [CrossRef] [PubMed]
39.	A. Bitzer and M. Walther, “Terahertz near-field imaging of metallic subwavelength holes and hole arrays,” Appl. Phys. Lett. 92(23), 231101 (2008). [CrossRef]
40.	K. Serita, S. Mizuno, H. Murakami, I. Kawayama, Y. Takahashi, M. Yoshimura, Y. Mori, J. Darmo, and M. Tonouchi, “Scanning laser terahertz near-field imaging system,” Opt. Express 20(12), 12959–12965 (2012). [CrossRef] [PubMed]
41.	D. H. Auston, K. P. Cheung, J. A. Valdmanis, and D. A. Kleinman, “Čherenkov radiation from sfemtosecond optical pulses in electro-optic media,” Phys. Rev. Lett. 53(16), 1555–1558 (1984). [CrossRef]
42.	P. U. Jepsen and B. M. Fischer, “Dynamic range in terahertz time-domain transmission and reflection spectroscopy,” Opt. Lett. 30(1), 29–31 (2005). [CrossRef] [PubMed]
43.	W. Withayachumnankul, H. Lin, K. Serita, C. M. Shah, S. Sriram, M. Bhaskaran, M. Tonouchi, C. Fumeaux, and D. Abbott, “Sub-diffraction thin-film sensing with planar terahertz metamaterials,” Opt. Express 20(3), 3345–3352 (2012). [CrossRef] [PubMed]

OCIS Codes
(050.1220) Diffraction and gratings : Apertures
(180.4243) Microscopy : Near-field microscopy
(300.6495) Spectroscopy : Spectroscopy, teraherz
(110.6795) Imaging systems : Terahertz imaging

ToC Category:
Physical Optics

History
Original Manuscript: May 31, 2012
Revised Manuscript: July 20, 2012
Manuscript Accepted: July 21, 2012
Published: July 26, 2012

Virtual Issues
Vol. 7, Iss. 9 Virtual Journal for Biomedical Optics

Citation
Christophe Fumeaux, Hungyen Lin, Kazunori Serita, Withawat Withayachumnankul, Thomas Kaufmann, Masayoshi Tonouchi, and Derek Abbott, "Distributed source model for the full-wave electromagnetic simulation of nonlinear terahertz generation," Opt. Express 20, 18397-18414 (2012)
http://www.opticsinfobase.org/vjbo/abstract.cfm?URI=oe-20-16-18397

Sort: Author | Year | Journal | Reset

References

A. Rice, Y. Jin, X. F. Ma, X.-C. Zhang, D. Bliss, J. Larkin, and M. Alexander, “Terahertz optical rectification from <110> zinc-blende crystals,” Appl. Phys. Lett.64(11), 1324–1326 (1994). [CrossRef]
G. Dakovski, B. Kubera, and J. Shan, “Localized terahertz generation via optical rectification in ZnTe,” J. Opt. Soc. Am. B22(8), 1667–1670 (2005). [CrossRef]
J. Z. Xu and X.-C. Zhang, “Optical rectification in an area with a diameter comparable to or smaller than the center wavelength of terahertz radiation,” Opt. Lett.27(12), 1067–1069 (2002). [CrossRef] [PubMed]
J. B. Khurgin, M. W. Pruessner, T. H. Stievater, and W. S. Rabinovich, “Suspended AlGaAs waveguides for tunable difference frequency generation in mid-infrared,” Opt. Lett.33(24), 2904–2906 (2008). [CrossRef] [PubMed]
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. Express17(16), 13502–13515 (2009). [CrossRef] [PubMed]
Q. Chen, Z. Jiang, G. X. Xu, and X.-C. Zhang, “Near-field terahertz imaging with a dynamic aperture,” Opt. Lett.25(15), 1122–1124 (2000). [CrossRef] [PubMed]
T. Kiwa, M. Tonouchi, M. Yamashita, and K. Kawase, “Laser terahertz-emission microscope for inspecting electrical faults in integrated circuits,” Opt. Lett.28(21), 2058–2060 (2003). [CrossRef] [PubMed]
W. Withayachumnankul, G. M. Png, X. Yin, S. Atakaramians, I. Jones, H. Lin, B. S. Y. Ung, J. Balakrishnan, B. W.-H. Ng, B. Ferguson, S. P. Mickan, B. M. Fischer, and D. Abbott, “T-ray sensing and imaging,” Proc. IEEE95(8), 1528–1558 (2007). [CrossRef]
R. Lecaque, S. Grésillon, and C. Boccara, “THz emission microscopy with sub-wavelength broadband source,” Opt. Express16(7), 4731–4738 (2008). [CrossRef] [PubMed]
Y. S. Lee, Principles of Terahertz Science and Technology (Springer, 2008).
T. Weiland, “A discretization method for the solution of Maxwell’s equations for six-component fields,” AEU, Int. J. Electron. Commun.31(3), 116–120 (1977).
A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed. (Artech House 2005).
J.-M. Jin, The Finite Element Method in Electromagnetics, 2nd ed. (Wiley 2002).
W. C. Gibson, The Method of Moments in Electromagnetics (Chapman & Hall/CRC, 2008).
C. M. Dissanayake, M. Premaratne, I. D. Rukhlenko, and G. P. Agrawal, “FDTD modeling of anisotropic nonlinear optical phenomena in silicon waveguides,” Opt. Express18(20), 21427–21448 (2010). [CrossRef] [PubMed]
I. Ahmed, E. H. Khoo, O. Kurniawan, and E. P. Li, “Modeling and simulation of active plasmonics with the FDTD method by using solid state and Lorentz–Drude dispersive model,” J. Opt. Soc. Am. B28(3), 352–359 (2011). [CrossRef]
M. Neshat, D. Saeedkia, and S. Safavi-Naeini, “Semi-analytical calculation of terahertz signal generated from photocurrent radiation in traveling-wave photonic mixers,” Int. J. Infrared Millim. Waves29(9), 809–822 (2008). [CrossRef]
P. Bonnet, X. Ferrieres, B. L. Michielsen, P. Klotz, and J. L. Roumiguières, “Finite-volume time domain method,” in Time Domain Electromagnetics (S. M. Rao, Ed. San Diego, CA: Academic Press, 1999).
C. Fumeaux, D. Baumann, P. Leuchtmann, and R. Vahldieck, “A generalized local time-step scheme for efficient FVTD simulations in strongly inhomogeneous meshes,” IEEE Trans. Microw. Theory Tech.52(3), 1067–1076 (2004). [CrossRef]
H. Lin, C. Fumeaux, B. M. Fischer, and D. Abbott, “Modelling of sub-wavelength THz sources as Gaussian apertures,” Opt. Express18(17), 17672–17683 (2010). [CrossRef] [PubMed]
H. Lin, C. Fumeaux, B. Seam Yu Ung, and D. Abbott, “Comprehensive modeling of THz microscope with a sub-wavelength source,” Opt. Express19(6), 5327–5338 (2011). [CrossRef] [PubMed]
B. E. A. Saleh and M. C. Teich, Fundamentals of Photonics (John Wiley & Sons, 1991).
E. D. Palik, ed., Handbook of Optical Constants of Solids (Academic Press, 1998)
C. A. Balanis, Antenna Theory, 3rd Edition (Wiley, 2005).
D. Li and G. Ma, “Pump-wavelength dependence of terahertz radiation via optical rectification in (110)-oriented ZnTe crystal,” J. Appl. Phys.103(12), 123101 (2008). [CrossRef]
D. H. Auston, “Subpicosecond electro-optic shock waves,” Appl. Phys. Lett.43(8), 713–715 (1983). [CrossRef]
D. A. Kleinman and D. H. Auston, “Theory of electrooptic shock radiation in nonlinear optical media,” IEEE J. Quantum Electron.20(8), 964–970 (1984). [CrossRef]
C. Fumeaux, D. Baumann, S. Atakaramians, and E. Li, “Considerations on paraxial Gaussian beam source conditions for time-domain full-wave simulations,” Annual Rev. of Progress in Appl. Computational Electromagnetics, 401–406 (2009).
N. K. Madsen and R. W. Ziolkowski, “A three-dimensional modified finite volume technique for Maxwell’s equations,” Electromagnetics10(1-2), 147–161 (1990). [CrossRef]
V. Shankar, A. H. Mohammadian, and W. F. Hall, “A time-domain, finite-volume treatment for the Maxwell equations,” Electromagnetics10(1-2), 127–145 (1990). [CrossRef]
C. Fumeaux, K. Sankaran, and R. Vahldieck, “Spherical perfectly matched absorber for finite-volume 3-D domain truncation,” IEEE Trans. Microw. Theory Tech.55(12), 2773–2781 (2007). [CrossRef]
D. Baumann, C. Fumeaux, C. Hafner, and E. P. Li, “A modular implementation of dispersive materials for time-domain simulations with application to gold nanospheres at optical frequencies,” Opt. Express17(17), 15186–15200 (2009). [CrossRef] [PubMed]
M. C. Hoffmann, K.-L. Yeh, J. Hebling, and K. A. Nelson, “Efficient terahertz generation by optical rectification at 1035 nm,” Opt. Express15(18), 11706–11713 (2007). [CrossRef] [PubMed]
E. K. Rahani and T. Kundu, “Electromagnetic THz radiation modeling by DPSM,” Int. J. Infrared Millim. Waves33(3), 376–390 (2012). [CrossRef]
D. Côté, J. E. Sipe, and H. M. van Driel, “Simple method for calculating the propagation of terahertz radiation in experimental geometries,” J. Opt. Soc. Am. B20(6), 1374–1385 (2003). [CrossRef]
K. Wang, D. M. Mittleman, N. C. J. van der Valk, and P. C. M. Planken, “Antenna effects in terahertz apertureless near-field optical microscopy,” Appl. Phys. Lett.85(14), 2715–2717 (2004). [CrossRef]
A. J. L. Adam, J. M. Brok, M. A. Seo, K. J. Ahn, D. S. Kim, J. H. Kang, Q. H. Park, M. Nagel, and P. C. Planken, “Advanced terahertz electric near-field measurements at sub-wavelength diameter metallic apertures,” Opt. Express16(10), 7407–7417 (2008). [CrossRef] [PubMed]
M. A. Seo, A. J. L. Adam, J. H. Kang, J. W. Lee, K. J. Ahn, Q. H. Park, P. C. M. Planken, and D. S. Kim, “Near field imaging of terahertz focusing onto rectangular apertures,” Opt. Express16(25), 20484–20489 (2008). [CrossRef] [PubMed]
A. Bitzer and M. Walther, “Terahertz near-field imaging of metallic subwavelength holes and hole arrays,” Appl. Phys. Lett.92(23), 231101 (2008). [CrossRef]
K. Serita, S. Mizuno, H. Murakami, I. Kawayama, Y. Takahashi, M. Yoshimura, Y. Mori, J. Darmo, and M. Tonouchi, “Scanning laser terahertz near-field imaging system,” Opt. Express20(12), 12959–12965 (2012). [CrossRef] [PubMed]
D. H. Auston, K. P. Cheung, J. A. Valdmanis, and D. A. Kleinman, “Čherenkov radiation from sfemtosecond optical pulses in electro-optic media,” Phys. Rev. Lett.53(16), 1555–1558 (1984). [CrossRef]
P. U. Jepsen and B. M. Fischer, “Dynamic range in terahertz time-domain transmission and reflection spectroscopy,” Opt. Lett.30(1), 29–31 (2005). [CrossRef] [PubMed]
W. Withayachumnankul, H. Lin, K. Serita, C. M. Shah, S. Sriram, M. Bhaskaran, M. Tonouchi, C. Fumeaux, and D. Abbott, “Sub-diffraction thin-film sensing with planar terahertz metamaterials,” Opt. Express20(3), 3345–3352 (2012). [CrossRef] [PubMed]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Click here to see a list of articles that cite this paper