Fully explicit nonlinear optics model in a particle-in-cell framework
Journal of Computational Physics 250 (2013) 388–402
Contents lists available at SciVerse ScienceDirect
Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp
Fully explicit nonlinear optics model in a particle-in-cell framework D.F. Gordon ⇑, M.H. Helle, J.R. Peñano Plasma Physics Division, Naval Research Laboratory, Washington, DC 20375, United States
a r t i c l e
i n f o
Article history: Received 10 August 2012 Received in revised form 25 January 2013 Accepted 5 May 2013 Available online 24 May 2013 Keywords: Particle-in-cell Nonlinear optics Finite difference time domain
a b s t r a c t A numerical technique which incorporates the nonlinear optics of anisotropic crystals into a particle-in-cell framework is described. The model is useful for simulating interactions between crystals, ultra-short laser pulses, intense relativistic electron bunches, plasmas, or any combination thereof. The frequency content of the incident and scattered radiation is limited only by the resolution of the spatial and temporal grid. A numerical stability analysis indicates that the Courant condition is more stringent than in the vacuum case. Numerical experiments are carried out illustrating the electro-optic effect, soliton propagation, and the generation of fields in a crystal by a relativistic electron bunch. Published by Elsevier Inc.
1. Introduction Nonlinear optics is often described in the frequency domain so that dispersive effects can be treated in a simple way. The convolution integrals associated with nonlinear effects are simplified by assuming that the radiation spectrum takes the form of a discrete set of narrow peaks. In the case of few cycle pulses, the convolution integrals cannot be reduced, and a time domain model becomes attractive. The finite-difference-time-domain (FDTD) technique is a well established method for solving the exact Maxwell equations in a bounded or periodic domain [1]. It may be used in connection with the particle-in-cell (PIC) technique to self-consistently solve for the motions of a large number of charged particles in an electromagnetic field [2]. This approach is often used to model fully nonlinear laser-plasma interactions and beam-plasma interactions. In this work, it is extended to account for nonlinear laser-crystal and beam-crystal interactions. This is accomplished by incorporating a model for bound particles into the PIC code turboWAVE [3]. The model is fully parallelized, and runs in up to three dimensions. The PIC aspect of turboWAVE has much in common with a number of other codes designed to model laser-plasma interactions [4–10]. The nonlinear optics aspect is related to a number of other codes designed to model ultrashort pulse propagation in a nonlinear medium [11–16]. Both types of codes can be divided into those that describe the radiation by means of a complex envelope, and those that are fully explicit, i.e., those that resolve the optical time scale. TurboWAVE supports both models, but in this work only the fully explicit model is used. As a result, there is no assumption about the frequency content of the radiation, and given the right model for the dielectric response, all nonlinear and dispersive effects are accounted for. Our model for the dielectric response generalizes the auxiliary equation technique described in Refs. [13,15]. In particular, bound charges in the dielectric are represented by effective particles, whose contribution to the four-current is computed using PIC techniques. The effective particles are subjected to forces arising from a superposition of macroscopic and microscopic fields. The macroscopic field is the usual electromagnetic field that is computed in any FDTD code, while the ⇑ Corresponding author. Tel.: +1 202 767 5036. E-mail address: [email protected] (D.F. Gordon). 0021-9991/$ - see front matter Published by Elsevier Inc. http://dx.doi.org/10.1016/j.jcp.2013.05.014
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402
389
microscopic field is a three dimensional electrostatic potential representing, e.g., an atomic binding potential. The resulting equation of motion is expanded as an anharmonic oscillator equation. By using multiple species of effective particles, each satisfying a different oscillator equation, the material response can be tailored to match that of real materials over a broad range of frequencies. Free charges are self-consistently incorporated into the model by linear superposition of the free and bound sources in the FDTD field solver.1 The resulting model is capable of modeling interactions among laser pulses, particle beams, plasmas, and dielectric crystals, in a fully dispersive, nonlinear, anisotropic, and kinetic way. To our knowledge, no other model has been implemented which simultaneously accounts for all these effects. One can argue that the auxiliary equation technique could be similarly extended, since, after all, one can recover any physics with the right ‘‘auxiliary equations.’’ The effective particle model described here has the advantage that all of the required physics falls out of a single unifying physical principle. In Ref. [17], the turboWAVE extensions described in this paper were introduced for the first time. In Ref. [18], the model was applied to a novel electro-optic diagnostic technique. In the present paper, the model is described in more detail, and a formal analysis of numerical stability is given. The model is extended to include third order nonlinearities, and the implementation is demonstrated via three numerical experiments. 2. Description of the model The numerical model described here solves the exact Maxwell equations:
@D @t @B rE¼ @t rD¼q
rH¼Jþ
ð1aÞ ð1bÞ ð1cÞ
rB¼0
ð1dÞ
Here, B ¼ l0 H þ M and D ¼ 0 E þ P, where P is the polarization and M is the magnetization. In the present work, we take M ¼ 0. The polarization can be expressed in terms of a spatially smoothed four-current due to bound charges, denoted by ðhgi; hjiÞ. This results in [19]
r B ¼ l0 ðJ þ hjiÞ þ
1 @E c2 @t
@B rE¼ @t r E ¼ ðq þ hgiÞ=0 rB¼0
ð2aÞ ð2bÞ ð2cÞ ð2dÞ
With this formulation, any field solver designed to account for free charges can be easily adapted to account for bound charges by making the substitution
ðq; JÞ ! ðq þ hgi; J þ hjiÞ
ð3Þ
Our model for ðhgi; hjiÞ is based on calculating the motions of effective particles that respond to the superposition of a microscopic binding potential and a macroscopic electric field E. Denoting the displacement of an effective particle from its equilibrium position by r, and expanding the microscopic potential in a Taylor series, results in the effective particle equation of motion 2
d ri dt
2
þ
X dr j q 2Cij þ ðX2 Þij r j þ aijk rj r k bri r j r j ¼ Ei m dt jk
ð4Þ
where the subscripts vary over Cartesian coordinates, q=m is the charge to mass ratio, and the v B force is neglected. The anisotropy of the medium is represented by the tensors C; X2 , and a. These are, respectively, the damping rate, the square of the resonant frequency, and the first anharmonic coefficient. At present, the second anharmonic coefficient, b, is assumed scalar. By means of a coordinate transformation, X2 can always be made diagonal. The basis vectors which induce this transformation will be denoted by ðeu ; ev ; ew Þ. The basis vectors used for the general calculation will be denoted ðex ; ey ; ez Þ. The crystallographic basis vectors, which generally coincide with ðeu ; ev ; ew Þ, are denoted (h1 0 0i, h0 1 0i, h0 0 1i). In general, the parameters in Eq. (4) depend not only on tensor indices, but also on an index p that identifies a particular particle. In other words, each particle may move in a unique potential well, and may therefore satisfy a unique equation of motion. The system of the particle together with its equation of motion is called an oscillator. In practical implementations, it is convenient to group all particles that satisfy the same equation of motion into an oscillator species with index s. Then the 1 The coupling of free and bound charges is through the macroscopic field only, i.e., collisions are neglected. As a result, the model is most useful when the free charges do not impinge on the dielectric.
390
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402
material parameters, q; m; X2 ; C; a, and b, depend on s, while only r is explicitly dependent on p. As discussed below, the macroscopic sources ðhgi; hjiÞ can be derived from the set rp using a variation on source deposition techniques developed for PIC codes. 3. Analytical dispersion relation In order to analyze the model described above, consider the locally averaged displacement hrs i associated with each oscillator species. In the case of a uniform medium, this is related to the spatially smoothed four-current by
ðhgi; hjiÞ ¼
X s
@ qs Ns r fhrs i; fhrs i @t
ð5Þ
where N s is the density of oscillators, and the tensor f is an oscillator strength. Assuming local phase coherence, hrs i satisfies the same equation as rp , which in the ðeu ; ev ; ew Þ basis, after linearizing, is
! @2 @ qs @Ai @/ 2 : þ 2 C þ X i ¼ þ hr iis is iis @t ms @t @ni @t2
ð6Þ
where i is the coordinate index, ni is the ith coordinate, A is the vector potential, and / is the scalar potential. Using the Coulomb gauge, the radiation is described by
c2 r2
! X q Ns @hris i @ 2 / @2 Ai ¼ fiis s þ 2 0 @t @t@ni @t s
ð7Þ
Consider any mode with a wave-vector that is collinear with one of the principal axes, ðeu ; ev ; ew Þ. Then, it can be shown that / ¼ 0, and the transverse components of A satisfy the dispersion relation 2
x
! X fiis x2ps 2 c2 k ¼ 0 1þ Dis ðxÞ s
ð8Þ
where x2ps ¼ N s q2s =0 ms , and
Dis ðxÞ ¼ X2iis 2ixCiis x2
ð9Þ
In the case of a medium comprised of a single, lossless, isotropic species with unit oscillator strength, the dispersion relation is 2
ðX2 x2 Þðx2 c2 k Þ þ x2 x2p ¼ 0
ð10Þ
The primary feature of this dispersion relation is a stop-band between the resonant frequency, X, and the cutoff frequency, qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X2 þ x2p . 4. Oscillator parameters Measured dielectric properties are usually given as susceptibilities in the frequency domain. Hence, in order to apply the time dependent model described above to real materials, it is necessary to connect the constants X2 ; C; f; a, and b, with the susceptibilities vðnÞ ðx1 ; x2 ; . . .Þ. In the ðeu ; ev ; ew Þ basis, the linear susceptibility has only the diagonal elements
vii ðxÞ ¼
X fiis x2ps Dis ðxÞ s
ð11Þ
Following Ref. [20], the first anharmonic coefficient is related to
Xq fs x2ps s ðaijk Þs ms Ds ðx1 ÞDs ðx2 ÞDs ðx3 Þ s
ð2Þ vijk ð x1 ; x2 ; x3 Þ ¼
vð2Þ by ð12Þ
where for simplicity, a cubic crystal is assumed. The nonlinear susceptibility is related to the electro-optic coefficient rijk via ð2Þ vijk ðx þ dx; x; dxÞ ¼
1X il ðxÞrlmk ðx; dxÞmj ðxÞ 2 lm
ð13Þ
Here, the factor of 1=2 accounts for the degeneracy of sum and difference generation when dx ! 0, and is the relative permittivity. Finally, the nonlinear refractive index is given by
n2 ðxÞ ¼
x2ps 3 g0 X q2s bs 4 n0 ðxÞ2 s m2s Ds ðxÞDs ðxÞ3
ð14Þ
391
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402 Table 1 Oscillator parameters for gallium phosphide. Parameter
Lattice oscillator
Electronic oscillator
Oscillator strength, f Resonance frequency, X
1
1
6:90 1013 rad/s
Damping frequency, c
6:25 1010 rad/s
6:38 1015 rad=s 0
Plasma frequency, xp
9:27 1013 rad/s
1:78 1016 rad/s
Anharmonic coefficient, a123
3:28 1037 m1 s2
4:1 1041 m1 s2
Fig. 1. Nonlinear Lorentz model for GaP. Panel (a) displays RðÞ, where the solid curve is the two-oscillator (2O) model and the points are from the ð2Þ expressions in Ref. [21]. Panel (b) is similar to (a) except IðÞ is displayed. Panel (c) displays v123 ðxL dx; xL ; dxÞ, where xL is a probe laser frequency, given by 2pc=xL ¼ 0:6328 lm. The solid curve is the 2O model, and the dashed curve is the expression in Ref. [22]. The triangles indicate the infra-red frequencies where the data points in [22] are located. Panel (d) displays r 123 ðx; dxÞ ¼ r 41 ðx; dxÞ, where dx ! 0. The curve is the 2O model, and the points are the data from Ref. [23]. The 2O model cannot be made consistent with both Refs. [22,23].
where n0 is the linear refractive index and g0 is the impedance of free space (see Appendix). Note that the charge to mass ratio qs =ms amounts to an extraneous free parameter that could have been absorbed into qs N s ; as , and bs . In this work, the electronic charge to mass ratio is always assumed. As an example, consider the cubic crystal gallium phosphide (GaP), which is useful for electro-optic sensing applications [24–28]. Expressions fitting the measured linear dispersion have already been given in Ref. [21]. There, two separate expressions were used for the optical and THz regimes. A two-oscillator model based on the parameters given in Table 1, agrees with both expressions, as illustrated in Fig. 1(a) and (b). Once the linear dispersion is determined, the only remaining free parameters are the anharmonic coefficients associated with each oscillator. In order to choose these parameters, an attempt was made to simultaneously fit data from Refs. [22,23]. ð2Þ In Ref. [22], vijk ðxL dx; xL ; dxÞ is reported, where xL is the frequency of a probing helium–neon laser, and dx is one of several infra-red frequencies produced by a multi-line H2þO2 laser. In Ref. [23], rijk ðx; dxÞ is measured, where dx is in the radio frequency (RF) range, and x takes several values in the optical range. It was found that the two-oscillator model cannot be made to agree with both experiments. The parameters of Table 1 are chosen to obtain agreement with the measured rijk ðx; dxÞ. The resulting dispersion is displayed in Fig. 1(c) and (d). In order to make the two-oscillator model agree with vð2Þ ðxL dx; xL ; dxÞ, a much larger absolute value of the lattice anharmonic coefficient must be used. This leads, roughly, to the oscillator parameters given in Ref. [17], although for the simulations presented there and in Ref. [18], the lattice oscillator was taken to be linear.
5. Numerical technique 5.1. Oscillator equation and source deposition 1
1
Most electromagnetic PIC codes use the leap-frog technique. More specifically, Enþ1 is computed using En ; Bnþ2 and Jnþ2 , 3 1 where n is the time level. Then, Bnþ2 is computed using Bnþ2 and Enþ1 . In some cases a Poisson solver is used to refine
392
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402
Fig. 2. Charge Deposition. (a) Geometry of a grid cell in relation to the points on the Yee mesh. The three dimensional arrows point to the mesh points where the corresponding electric field components are known. (b) Quadratic Weighting. The bound charges are represented by charge clouds (solid curve) whose equilibrium position (dashed curve) is the centroid of the cell volume. For an electron, the charge in a cell is increased (decreased) by the area of the intersection of the blue (red) area with the cell. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Enþ1 using qnþ1 . As discussed above, any such field solver can be used in the presence of an oscillator population by making the substitution ðq; JÞ ! ðq þ hgi; J þ hjiÞ. This implies that hgi is known at integral time levels, and hji is known at half-integral time levels. 1 To compute hginþ1 and hjinþ2 , the displacements rnp and rnþ1 are needed. Let T represent the matrix of transformation of a p vector from the ðeu ; ev ; ew Þ basis to the ðex ; ey ; ez Þ basis. If all quantities characterizing the oscillators are stored in the ðeu ; ev ; ew Þ basis, the displacements can be updated using2
rnþ1 2rni þ rn1 r nþ1 r ni F ni 2 n i i i þ 2 C þ ð X Þ r ¼ ii i ii Dt m Dt 2
ð15Þ
where Dt is the time step and
Fn ¼ qT1 En marn rn þ mbðrn rn Þrn
ð16Þ rnp
The macroscopic four-current ðhgi; hjiÞ associated with the microscopic orbits can be computed using PIC techniques. These techniques regard each particle as a ‘‘cloud’’ with dimensions on the order of a cell size. Knowing the orbit of each cloud’s centroid allows one to compute the amount of charge in each cell as a function of time, as well as the current passing through each cell wall. Various schemes have been developed to facilitate this calculation [2,29,30]. In the case of bound charges, these schemes have to be applied with care in order to minimize round-off errors associated with very small displacements. To deposit the sources due to bound charges, place one oscillator, per species, in each grid cell.3 The cell geometry is shown in Fig. 2(a). The deposition of charge, in one dimension, is illustrated in Fig. 2(b). The oscillator is represented by two oppositely charged density distributions. The equilibrium distribution, shown as a dashed line, is immobile. The displaced distribution, shown as a solid line, executes the orbit described by Eq. (15). The net charge in a cell is the integral over the cell volume of the difference between the distributions. To minimize round-off error, this integral should not be carried out by re-using a PIC routine that deposits each distribution separately. Instead, expressions should be written for the net charge in a cell, and expanded to lowest order in the displacement. Then, for the oscillator in cell ði; j; kÞ, the net charge density is distributed in the surrounding cells as follows:
N s qs 2 Dx n N s qs ¼ ey Tfrp 2 Dy N s qs ¼ ez Tfrnp 2 Dz
hgini1;j;k ¼ ex Tfrnp
ð17aÞ
hgini;j1;k
ð17bÞ
hgini;j;k1
ð17cÞ
Here, the cell dimension is Dx Dy Dz. The current density through the surrounding cell walls follows directly from charge conservation:
2 3
Forward differencing of the damping term improves numerical stability, see Eq. (23). In the present work, we assume local phase coherence. By using many effective particles per cell, this assumption could be removed.
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402
Dx nþ1 n ex hjii12;j;k ¼ hginþ1 i1;j;k hgii1;j;k Dt 2 Dy nþ12 nþ1 n ey hjii;j1;k ¼ hgii;j1;k hgii;j1;k Dt 2 Dz nþ12 nþ1 n ez hjii;j;k1 ¼ hgii;j;k1 hgii;j;k1 Dt 2
393
ð18aÞ ð18bÞ ð18cÞ
Note that the above expressions are particular to the choice of quadratic weighting. This choice provides adequate smoothing, while also allowing for abrupt vacuum-dielectric transitions. A useful fact is that the current density due to an effective particle, through any cell wall enclosing that particle, is the normal component of 12 N s qs v p , where
vp
rnþ1 rnp p ¼ Tf Dt
! ð19Þ
is the velocity of the effective particle in the ðex ; ey ; ez Þ basis. 5.2. Linear stability and accuracy In an ordinary PIC code, the primary stability criterion is the Courant condition, which in one dimension reads cDt < Dz, where Dt is the time step and Dz is the grid spacing. It is clear that in the presence of a dispersionless dielectric, this would be modified to read cDt=n < Dz, where n is the refractive index. In the case of a fully dispersive model such as the one described above, the modification of the Courant condition is more subtle. Expressing Eqs. (6) and (7) in finite difference form, assuming propagation along a principal axis, and inserting exponential forms for Ai and hris i, gives the numerical dispersion relation
cos xDt ¼ VðkÞ þ P i ðxÞ
ð20Þ
where
VðkÞ ¼ 1 2c2 Dt 2
2 1 k Dx 2 x Dx2
sin
þ
2 1 k Dy 2 y Dy 2
sin
þ
2 1 k Dz 2 z Dz 2
sin
! ð21Þ
and
P i ðxÞ ¼
X s
x2ps Dt2 sin2 12 xDt 2 1 ðX Þiis Dt ð1 þ Ciis DtÞ sin 12 xDt 2i Ciis Dt sin xDt 4 1 f 2 iis
2
2
ð22Þ
Note that when the density of oscillators vanishes, P i ! 0 and the well known vacuum equation is recovered. The effect of P i can be estimated by assuming X2 Dt 2 1 and kx Dx ky Dy kz Dz xDt p. These assumptions are motivated by the fact that the resonance frequency should be well resolved, and by the expectation that numerical instabilities are most severe at the Nyquist frequency. Under these assumptions, the stability criterion becomes
c2 c2 c2 Dt K þ 2þ 2 2 Dx Dy Dz
1=2
2 1 X fiis x2ps Dt 1 4 s 1 þ Ciis Dt
!1=2 ð23Þ
This is more stringent than the usual Courant condition. Numerical experiments show that Eq. (23) is accurate in practice. In order to evaluate the accuracy and stability of the differencing scheme more precisely, the numerical dispersion relation is evaluated numerically for the case of a single, lossless oscillator species. Since the dispersion relation scales, we normalize all frequencies to the resonant frequency, X, taking xp ¼ 2X; Dz ¼ 0:5c=X, and Dx ¼ Dy ! 1. The results are shown in Fig. 3 for three choices of the time step. The case Dt ¼ 0:2=X is shown in panel (a), the case Dt ¼ 0:4=X is shown in panel (b), and the case Dt ¼ 0:48=X is shown in panel (c). The analytical dispersion relation (10) is overlayed for comparison. As usual, the best accuracy is obtained for small values of x and k. As cDt=Dz increases, the accuracy improves, until numerical instability sets in. This behavior is similar to the vacuum case. However, unlike the vacuum case, instability occurs for cDt=Dz < 1, as shown in panel (c). Better accuracy may be obtained by decreasing both Dz and Dt proportionately. 5.3. Nonlinear stability In the case of the fully nonlinear system, a rigorous stability analysis is difficult. However, if the oscillator equation is viewed as describing an effective particle in a potential well, the stability criterion may be viewed as the requirement that the particle never enter any region of runaway acceleration. Such a region exists if b > 0, for in this case there is a divergent repulsive force as r ! 1. Instability is also possible if b ¼ 0 and a – 0. Assuming there are no metastable regions, the stability criterion can be written
FNL r < mX2 r r
ð24Þ
394
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402
Fig. 3. Numerical dispersion for (a) cDt=Dz ¼ 0:4, (b) cDt=Dz ¼ 0:8, and (c) cDt=Dz ¼ 0:96. Solid curves are the analytical dispersion relation (10), solid circles are Rðx=XÞ, and open circles connected by a thin line are Iðx=XÞ. In all cases, xp =X ¼ 2, and Dz ¼ 0:5c=X. Panel (c) shows that numerical instability may occur for cDt=Dz < 1.
where
FNL ¼ mbðr rÞr marr
ð25Þ
If one is interested primarily in second order effects, the system can be stabilized by taking b < 0. The magnitude of b is chosen such that the order of magnitude of the first, second, and third order forces are the same at some critical point. This leads to
b
a2
X2
ð26Þ
where a and X are typical values of the elements of a and X. If third order effects are important physically, a stabilizing fifth order nonlinearity can be easily introduced. In this case,
FNL ¼ mdðr rÞ2 r þ mbðr rÞr marr
ð27Þ
where d is the fifth order anharmonic coefficient. Demanding that the first, third, and fifth order forces be equal at some critical point gives 2
d
b
X2
ð28Þ
The ultimate nonlinear stabilization technique would utilize physical effects, such as ionization in the case of the electronic oscillator, or material damage in the case of the lattice oscillator. However, accounting for these processes is beyond the scope of the present work.
5.4. False resonance technique for accelerating computation The addition of one effective particle per grid cell to a vacuum FDTD field solver increases the computational workload by only about 25%. However, it sometimes happens that the accuracy and stability criteria, X2 Dt 2 1 and x2p Dt 2 1, force one to use a very large number of time steps to simulate a given process. This limitation can be mitigated if the frequency range of interest is well below X. In such cases, it is possible to artificially shift X to a lower frequency while preserving the dispersion relation in the region of interest. In this way, the ‘‘true resonance’’ is replaced by a ‘‘false resonance’’ which is easier to resolve, but leads to the same physics. As an example of a false resonance model, consider again GaP, which has the true resonance parameters displayed in Table 1. When modeling an electron bunch with a characteristic time of, say, 100 fs, the frequency range of interest is well below the electronic resonance, which has 2p=X 1 fs. One may then use a false resonance model such as the one described by Table 2. The linear and nonlinear dispersion curves associated with this model are shown in Fig. 4. In the THz range, the dispersion relation remains accurate. Note that the dispersion in this range could be fine tuned by adding additional false resonances. In Ref. [17], computation was accelerated by sub-cycling the oscillator equation with respect to a super-cycle time step used for all other equations. In order to make this stable, one must choose the super-cycle time step to be an integral number of periods of the under-resolved frequency. Apart from placing needless constraints on the simulation parameters, this scheme was found to give results that do not satisfactorily conserve energy, and which eventually diverge. In the case of the false resonance approach, one can prove linear stability, and as is shown below, energy is accurately conserved.
395
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402 Table 2 False resonance model for gallium phosphide. Parameter
Lattice oscillator
False resonance
Oscillator strength, f Resonance frequency, X
1
1
6:90 1013 rad/s
Damping frequency, C
6:25 1010 rad/s
8:91 1014 rad/s 0
Plasma frequency, xp
9:27 1013 rad/s
2:46 1015 rad/s
Anharmonic coefficient, a123
3:28 1037 m1 s2
1:23 1040 m1 s2
Fig. 4. False resonance model for GaP. Panel (a) displays RðÞ, where the solid curve is the false resonance model and the points are from the expressions in ð2Þ Ref. [21]. Panel (b) displays v123 ð2dx; dx; dxÞ, where the curve is the false resonance model and the points are the true resonance model (Table 1).
5.5. Boundary conditions and moving window In simulations of laser-plasma interactions, the moving window technique [31] is often used. In this technique, the mesh points are shifted by one cell in between time levels such that the computational region, or ‘‘window,’’ moves at the speed of light. Boundary conditions are greatly simplified by the requirements of causality. When the plasma density is low, the laser pulse group velocity is nearly the speed of light, and the pulse stays in the window for a long time. Hence, the window only has to be as long as the laser pulse, even if the interaction region is much longer. When the group velocity is significantly smaller than the speed of light, the pulse quickly falls behind a light frame window, and one might as well work in the lab frame. If the window speed is the group velocity, a fully dispersive model might allow precursor signals to reach the front of the window. Moreover, signals reflected from the back of the window might work their way toward the region of interest. One solution is to use perfectly matched layers (PML) [32] to absorb waves near the boundaries of the window. In order for this to be effective, it is necessary to use a large number of layers due to the fact that the moving window involves shifting the field data by one cell between time levels. In order that this shift should weakly perturb the solution, the PML conductivity has to change very gradually from one cell to the next. In the case of lab frame simulations, another option is to use Lindman’s absorbing boundary condition [33] in the longitudinal direction, and periodic boundary conditions transversely. If a Poisson equation has to be solved, the solution in the box can be matched to a decaying solution outside the box [30]. TurboWAVE supports all the options mentioned above. A simulation that uses Lindman’s boundary condition in the lab frame is discussed in Section 6.1. A simulation that uses PML absorbers with a moving window that moves at the group velocity is discussed in Section 6.2. A simulation that uses a moving window that moves at the speed of light is discussed in Section 6.3. 6. Numerical experiments The model described above has been implemented as a module within the turboWAVE [3] framework. In this section, three numerical experiments are presented. First, the electro-optic effect in GaP is demonstrated by comparing the numerical change in polarization with that predicted by the low frequency theory. Second, soliton propagation in fused silica is demonstrated. Finally, to take full advantage of the PIC framework, a relativistic electron bunch is passed near a GaP crystal, and the induced fields are examined.
396
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402 Table 3 Parameters for simulation of electro-optic effect. The fields are measured inside the dielectric. The THz field is a half-cycle pulse with base-to-base width sH . Parameter
Symbol
Value
Time step Space step Cells
cDt Dz Nz
0.012 lm 0.017 lm
Steps
Nt
THz field THz half wave Laser field Laser FWHM Laser wavelength Crystal length Crystal orientation
Ey
217 68 kV/cm 2.2 ps 18 kV/cm 220 fs 0.81 lm 120 lm Eq. 30
sH jE x j
sL k0 L T
213
6.1. Electro-optic effect In this numerical experiment, an x-polarized laser pulse is copropagated with a y-polarized THz half-wave in a GaP crystal, and the change in the polarization state of the laser pulse is measured. The length of the THz half-wave is chosen to be long enough so that the standard picture of the electro-optic effect applies. In this picture, the change in the impermeability4 tensor due to a DC electric field is
Dgij ¼ r ijk Ek
ð29Þ
For any Zinc–Blende type crystal, r ijk ¼ r 123 ¼ r41 if all three indices are distinct, and r ijk ¼ 0 otherwise. The crystal orientation is represented by the operator T, which can be viewed as a sequence of rotations that re-orients the crystal starting with the ðeu ; ev ; ew Þ and ðex ; ey ; ez Þ bases aligned. It turns out that the electro-optic effect is maximized if [34]
T ¼ Rz ðp=2ÞRx ðp=2ÞRz ðp=4Þ
ð30Þ
with the laser field polarized in the x-direction and the DC field polarized in the y-direction. Here, Ri is an active right-handed rotation about the ith coordinate axis. As an example, if E and Dg are given in the ðex ; ey ; ez Þ basis, and r is given in the ðeu ; ev ; ew Þ basis, then Dg ¼ TrT1 ET1 . To calculate the change in polarization due to the electro-optic effect, a principal axis transformation is usually made [34]. Alternatively, one may make direct use of the solution of the slowly varying envelope equation,
EðzÞ ¼ exp
ipDz Eðz ¼ 0Þ nk0
ð31Þ
where E is the complex electric field envelope, k0 is the free space laser wavelength, n is the refractive index of the unperturbed crystal, and D is the change in permittivity induced by the DC field. This equation remains valid even when D is a tensor and E is a vector, and is easily evaluated in a software environment that supports the exponential of a matrix. We now compare the prediction of Eq. (31) with the results of a turboWAVE simulation. In the simulation, an x-polarized laser pulse and a y-polarized THz half-wave are incident on a GaP crystal oriented as in Eq. (30). The simulation parameters are given in Tables 1–3. The lattice nonlinearity was stabilized with b ¼ a2 =X2 . The results from the simulation are shown in Fig. 5. The data are evaluated 60 lm into the crystal to avoid interference from reflections generated at either vacuumcrystal interface. Since the turboWAVE model uses a universal electric field, the optical and THz frequency components have to be extracted by means of a Fourier filter. The envelope of the optical wave is constructed by applying a hard-edged bandpass filter centered at the optical frequency, þx0 . The result is then shifted in frequency by x0 , multiplied by 2 (to account for the energy in the negative frequencies) and the inverse Fourier transform is applied. Based on Fig. 5, the polarization ratio at z ¼ 60 lm, and at times prior to the generation of internal reflections, is jE y =E x j2 0:0018. Applying the analytical formula with the inputs Ey ¼ 68 kV/cm, n ¼ 3:16, and r123 ¼ 0:89 pm/V gives jE y =E x j2 0:0019. Thus, the simulation model closely agrees with the theoretical prediction. An interesting feature of the simulation data is that the second y-polarized optical pulse is larger than the first. This appears to happen because the reflected x-polarized optical pulse continues to be rotated by the reflected y-polarized THz pulse. This leads to a coherent reinforcement of the reflected y-polarized optical pulse. 6.2. Soliton propagation In this numerical experiment, a 20 fs, 2.1 lm pulse is propagated in fused silica in both the linear and soliton regimes. The Sellmeier formula for fused silica can be matched exactly to a three oscillator model, the parameters of which are displayed 4
Impermeability is defined by g ¼ 1 , where
is the relative permittivity.
397
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402
Fig. 5. Simulated electric field vs. time at z ¼ 60 lm from the entrance plane of the crystal. (a) y-polarized THz field (low-pass Fourier filter) (b) x-polarized optical field envelope (band-pass Fourier filter with shift) (c) y-polarized optical field envelope (band-pass Fourier filter with shift). The incident optical field is a single, purely x-polarized pulse, and the incident THz field is a single, purely y-polarized half-wave. The additional optical pulses, and the late-time distortion of the THz field, are due to reflections.
Table 4 Oscillator parameters for fused silica. Parameter
Osc. 1
Osc. 2
Osc. 3
Oscillator strength, f Resonance frequency, X
1
1
1
2:35 1014 0
1:63 1016 0
2:78 1016 0
1:79 1014 0
1:06 1016 0
2:30 1016
rad/s rad/s
4:72 1054
m2 s2
Damping frequency, C Plasma frequency, xp Anharmonic coefficient, b
Unit rad/s
in Table 4. The second anharmonic coefficient is chosen to give a nonlinear refractive index of n2 ¼ 3 1016 cm2 =W. The choice to take b – 0 only in the third oscillator has no effect other than on the frequency dependence of n2 . The amplitude envelope of the fundamental soliton has the form sechðsÞ, where s ¼ ðt z=v g Þ=T 0 ; v g is the group velocity, and T 0 is the pulse width parameter. In simulations, it is convenient to have a function whose support is strictly in a finite region. Furthermore, in a model of the type used here, the initial conditions depend on z, and the pulse starts in vacuum. Hence, the appropriate initial condition is
Ex ðt ¼ 0; zÞ ¼ 1=4 jE 0 jsechðf=cT 0 Þ sinðx0 f=cÞCðf=cT 0 Þ
ð32Þ
where Ex is the field that is initialized in the vacuum region, jE 0 j is the desired peak field after propagating into the dielectric, f ¼ z z0 ; z0 fixes the spatial origin, and CðxÞ is a suitable clipping function. The clipping function used here is
CðxÞ ¼
8 > < > :
jxj 6 4; 4 < jxj < 6;
1 cos2 ½ðjxj 4Þp=4
ð33Þ
jxj P 6; 0
It should be noted that when the pulse enters the dielectric there is a small reflection that is not accounted for in Eq. (32). The parameters of the pulse, and the numerical parameters associated with the simulation, are displayed in Table 5. The intensity and electric field are related by I0 ¼ 1=2 jE 0 j2 =2g0 . The free space wavelength is k0 ¼ 2pc=x0 . The group velocity dispersion (GVD) coefficient is defined by cb2 ¼ @ 2 ð1=2 xÞ=@t 2 jx¼x0 . The dispersion length, LD ¼ T 20 =jb2 j, measures the propagation length needed to observe group velocity dispersion (GVD) effects. The nonlinear length, LNL ¼ c=x0 n2 I0 , measures the propagation length needed to observe nonlinear pulse distortions. The conditions b2 < 0 and LD ¼ LNL result in stable propagation of the pulse without distortion [20,35]. This is illustrated in Fig. 6. Comparison of panels (a) and (b) shows that when I0 ¼ 0:4 TW/cm2, the pulse maintains its original form after propagating nearly 3LD . Panel (c) shows that when I0 is reduced by a factor of 100, the pulse broadens significantly over the same propagation distance. 6.3. Electron bunch passing near a crystal The electro-optic effect can be used to measure the temporal structure of relativistic electron bunches [18,26–28]. The general principle is that the bunch self-fields scatter a probing optical pulse into a new polarization and frequency [18,21,24]. In order to make use of such a principle, the form of the bunch fields in the crystal has to be understood. In the simplest interpretation, the self-fields of the bunch are treated as the equivalent of a planar terahertz half-wave incident on the crystal. In reality, the picture may be complicated by effects such as Cherenkov-like wakes [17,36] and diffraction radiation [37,38]. In order to fully characterize these effects, we carry out simulations of the propagation of an electron
398
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402 Table 5 Parameters for simulation of soliton propagation. Parameter
Symbol
Value
Time step Space step Cells
cDt Dz Nz
0.014 lm 0.034 lm
Steps
Nt
Window speed PML thickness Intensity Pulse width Free space wavelength Dispersion length Nonlinear length GVD coefficient
vg
4 105 0:68c 8.7 lm 0:4 TW=cm2 20 fs 2.1 lm 2.8 mm 2.8 mm
– I0 T0 k0 LD LNL b2
212
1:42 1025 s2/m
Fig. 6. Simulation of 1D soliton propagation in fused silica. Plots show 1=2 E2x =g0 vs. v g t z, where v g t is held fixed in each plot, and v g is the velocity of the moving window. Panel (a) shows the initial pulse, (b) shows the pulse after 8 mm propagation, and (c) shows the result from (b) in the case where the initial intensity is I0 =100.
bunch near a GaP crystal. The parameters are given in Table 6, and the geometry is shown in Fig. 7. The simulations differ from the one carried out in Ref. [17] in that the bunch charge is varied until significant nonlinear distortions are observed in the bunch fields, and in that the false resonance technique is used in place of sub-cycling the oscillator equation. The false resonance oscillator parameters are the same as those given in Table 2, except that a stabilizing third order nonlinearity is added to each oscillator. The electron bunch is initialized in the presence of an immobile neutralizing charge, and is abruptly accelerated at the beginning of the simulation. This generates spurious transient fields that have to be attenuated. For this reason, the bunch is allowed to propagate for t0 ¼ 15 picoseconds before encountering the GaP crystal, whose plane of incidence is located at z ¼ 4200 lm. Convergence to the proper steady state is facilitated by placing PML absorbers at the transverse boundaries of the simulation region. The spurious self-force exerted on the bunch during the transient phase is small in the present context. Fig. 8 displays the electric field after the bunch has moved about 400 lm beyond the plane of incidence of the GaP crystal. More specifically, Fig. 8(a) and (b) show Ey ðt t0 ; x ¼ 0; y; zÞ for the cases Q ¼ 0:28 pC and Q ¼ 2:8 pC, where Q is the bunch
Table 6 Parameters for simulation of fields induced by an electron bunch. Parameter
Symbol
Value
Time step Space step Cells
cDt Dx ¼ Dy ¼ Dz Nx Ny Nz
0.168 lm 0.672 lm
Steps
Nt
Window speed Bunch diameter Bunch length Electron energy Crystal size Crystal orientation
v rb
3 104 c 8 lm 80 fs 250 MeV 605 360 1 lm3 Eq. (30)
sb ðc 1Þmc2 Lx Ly Lz T
210 210 210
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402
399
Fig. 7. Geometry of the simulation of an electron bunch passing near an electro-optic crystal. The electron bunch is loaded into a Gaussian weighted ellipsoid (green) and propagated in the z-direction adjacent to a crystal whose internal structure is oriented as in Eq. (30) (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this article.) .
charge. Fig. 8(c) and (d) show the corresponding line-outs at y ¼ 0. The primary features in Fig. 8(a) and (b) are the vertical phase fronts located in the range 400 < z ct < 300 lm, and the diagonal phase front that appears at all z ct. As discussed previously [17], these can be roughly associated with coherent transition radiation (CTR) and Cherenkov radiation, respectively. In electro-optic sensing applications, the CTR-like feature is supposed to retain the form of the vacuum bunch fields. However, as is well known, dispersive effects lead to distortions of the type shown in the figure. When the bunch charge is high enough to excite large nonlinear polarization currents in the material, the fields induced in the crystal may be further distorted. Such a case is shown in Fig. 8(b) and (d). For the chosen crystal orientation, the nonlinear polarization is ð2Þ
ð2Þ
Pð2Þ ¼ v123 E2y ex þ 2v123 Ex Ey ey
ð34Þ
Fig. 8. Simulation of fields induced by an electron bunch with (a,c) Q ¼ 0:28 pC and (b,d) Q ¼ 2:8 pC. The low charge case is essentially the linear result, while the high charge case exhibits nonlinear distortions. The three dimensional data is evaluated at x ¼ 0 in order to present a two dimensional plot. The green oval represents the location of the electron bunch, which propagates from left to right. The boundary between the crystal and vacuum is located at y ¼ 200 lm. The plane of incidence is located at z ct ¼ 435 lm. The line-outs in (c) and (d) are evaluated at y ¼ 0.
400
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402
Fig. 9. Conservation of energy during simulated bunch-crystal interaction.
Hence, a y-polarized THz field induces an x-polarized field with second harmonic components. This field then mixes with the original y-polarized field to produce a distorted y-polarized field. This distortion can be clearly seen by comparing panels (c) and (d) of Fig. 8. Qualitatively, the effect is to produce narrower features in the waveform. In an electro-optic sensing application, this could adversely affect the fidelity of the diagnostic. Finally, consider the flow of energy in the crystal during the simulation. For the specific form of the dielectric response considered in this paper, Poynting’s theorem can be written as
Z
dn S þ
0
@X
d 2 dt
Z
3
d rðE2 þ c2 B2 Þ þ
X
XdU p þ Lp ¼ 0 dt p2X
ð35Þ
where the first term is the energy radiated through a surface @ X, the second is the rate of change of field energy in the volume X, and the third is the rate of change of the energy of all the effective particles in X. Here, the effective particle energy is
Up ¼
2 1 dr p mp þ qp /p ðr p Þ 2 dt
ð36Þ
where mp is the mass, rp is the displacement, qp is the charge, and /p ðr p Þ is the microscopic electrostatic potential of the displaced particle. The rate of frictional damping for an effective particle is
Lp ¼ 2Cp mp
2 dr p dt
ð37Þ
Fig. 9 displays the three terms in Poynting’s theorem, integrated over time, along with the sum. The ‘‘radiated’’ curve is negative, corresponding to the fact that the bunch fields propagate into the crystal. The steep initial slope corresponds to the generation of the CTR-like feature, while the more gradual slope corresponds to the generation of the Cherenkov-like feature. The ‘‘field’’ and ‘‘material’’ curves show that half the incident energy goes to the fields, while the other half goes to the effective particles. The ‘‘sum’’ curve shows that energy is conserved, to the extent that it vanishes. One notices a slight discrepancy at the moment the bunch fields first enter the crystal. This may be related to the numerical process of generating a reflected wave.
7. Conclusions The physics of nonlinear crystals excited by laser pulses, plasmas, particle beams, or any combination of these, can be numerically simulated using an effective particle model based on particle-in-cell techniques. In this model, bound charges respond to a superposition of microscopic and macroscopic fields, where the former are chosen to satisfy the known dispersion characteristics of the material, and the latter are self-consistently computed using the usual FDTD technique. The resulting numerical model is extremely flexible and makes very few assumptions. It is, however, computationally expensive. Nevertheless, it is now possible to carry out three dimensional simulations of electro-optic sensing, where the bunch fields induced in an electro-optic crystal are computed with unprecedented detail. Other applications of this model might include modeling the nonlinear optics of nano-composites, or frequency mixing with ultra-short pulses.
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402
401
Acknowledgments This work was supported by NRL Base Funds and the Department of Energy. We acknowledge useful conversations with B. Hafizi and A. Ting. Appendix A In this appendix, the relationship between the nonlinear refractive index, n2 , and the second anharmonic coefficient, b, is derived. For simplicity, all quantities are reduced to scalars. Consider a perturbation expansion of the anharmonic oscillator equation. Suppressing the species index, and expressing all quantities in frequency space, the first order equation is
DðxÞr ð1Þ ¼
q EðxÞ m
ð38Þ
and the third order equation is
DðxÞr ð3Þ ¼ b r ð1Þ r ð1Þ r ð1Þ
ð39Þ
Here, the asterisk denotes convolution over frequency, and the electric field is in the form
EðxÞ ¼
E0 ½dðx x0 Þ þ dðx þ x0 Þ
2
ð40Þ
Inserting r ð1Þ into the third order equation gives
3 qE0 FðxÞ rð3Þ ¼ b 2m DðxÞ
ð41Þ
where
FðxÞ ¼
dðx 3x0 Þ 3dðx x0 Þ 3dðx þ x0 Þ dðx þ 3x0 Þ þ 2 þ þ D 3 ð x0 Þ D ðx0 ÞDðx0 Þ D2 ðx0 ÞDðx0 Þ D3 ðx0 Þ
ð42Þ
Restoring species indices, the nonlinear polarization is
X Ns qs rð3Þs
Pð3Þ ¼
ð43Þ
s
This quantity must now be connected to the refractive index. The refractive index is defined by n ¼ 1=2 , with the relative permittivity. Using D ¼ 0 E, this can be written as
n¼
0 E þ P 0 E
1=2 ð44Þ
For a weakly nonlinear interaction, P ð3Þ Pð1Þ , and
n n0 þ
Pð3Þ 2n0 0 E
ð45Þ
where n0 is the linear index, given by n20 ¼ 1 þ P ð1Þ =0 E. The nonlinear index, n2 , is conventionally defined by
n ¼ n0 þ n2 I0
ð46Þ
where I0 is the cycle averaged intensity. Choosing the Fourier transform convention such that E0 is the peak value of the electric field, one has the relationship I0 ¼ n0 E20 =2g0 . Finally, equating Eqs. (45) and (46), and evaluating at x0 , gives
n2 ðx0 Þ ¼
x2ps 3 g0 X q2s bs 2 2 4 n0 ðx0 Þ s ms Ds ðx0 ÞDs ðx0 Þ3
ð47Þ
References [1] K. Yee, Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media, IEEE Trans. Antennas Propag. 14 (1966) 302. [2] J.P. Verboncoeur, particle simulation of plasmas: review and advances, Plasma Phys. Control. Fusion 47 (2005) A231–A260. [3] D. Gordon, Improved ponderomotive guiding center algorithm, IEEE Trans. Plasma Sci. 35 (2007) 1486–1488. [4] D.W. Forslund, J.M. Kindel, E.L. Lindman, Plasma simulation studies of stimulated scattering processes in laser-irradiated plasmas, Phys. Fluids 18 (1975) 1017–1030. [5] C.D. Decker, W.B. Mori, T. Katsouleas, Particle-in-cell simulations of raman forward scattering from short-pulse high-intensity lasers, Phys. Rev. E 50 (1994) 3338. [6] P. Mora, T. Antonsen, Kinetic modeling of intense, short laser pulses propagating in tenuous plasmas, Phys. Plasmas 4 (1997) 217.
402
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402
[7] A. Pukhov, J. Meyer-ter-Vehn, Relativistic laser-plasma interaction by multi-dimensional particle-in-cell simulations, Phys. Plasmas 5 (1998) 1880– 1886. [8] R. Fonseca, L. Silva, F. Tsung, V. Decyk, W. Lu, C. Ren, W. Mori, S. Deng, S. Lee, T. Katsouleas, J. Adam, OSIRIS: a three dimensional, fully relativistic particle in cell code for modeling plasma based accelerators, Lect. Notes Comput. Sci. 2331 (2002) 342–351. [9] C. Nieter, J. Cary, VORPAL: a versatile plasma simulation code, J. Comput. Phys. 196 (2004) 448–473. [10] C. Huang, V. Decyk, C. Ren, M. Zhou, W. Lu, W. Mori, J. Cooley, T. Antonsen, T. Katsouleas, QUICKPIC: a highly efficient particle-in-cell code for modeling wakefield acceleration in plasmas, J. Comput. Phys. 217 (2006) 658–679. [11] M. Mlejnek, E. Wright, J. Moloney, Dynamic spatial replenishment of femtosecond pulses propagating in air, Opt. Lett. 23 (1998) 382–384. [12] J. Peñano, P. Sprangle, B. Hafizi, A. Ting, D. Gordon, C. Kapetanakos, Propagation of ultra-short, intense laser pulses in air, Phys. Plasmas 11 (2004) 2865–2874. [13] M. Fujii, M. Tahara, I. Sakagami, W. Freude, P. Russer, High-order FDTD and auxiliary differential equation formulation of optical pulse propagation in two-dimensional Kerr and Raman nonlinear dispersive media, IEEE J. Quantum Electron. 40 (2004) 175–182. [14] J. Peñano, P. Sprangle, B. Hafizi, W. Manheimer, A. Zigler, Transmission of intense femtosecond laser pulses into dielectrics, Phys. Rev. E 72 (2005) 36412-1–36412-7. [15] J. Greene, A. Taflove, General vector auxiliary differential equation finite-difference time-domain method for nonlinear optics, Opt. Exp. 14 (2006) 8305–8310. [16] S. Champeaux, L. Berge, D. Gordon, A. Ting, J. Peñano, P. Sprangle, (3 + 1)-Dimensional numerical simulations of femtosecond laser filaments in air: toward a quantitative agreement with experiments, Phys. Rev. E 77 (2008). [17] D. Gordon, M. Helle, D. Kaganovich, A. Ting, Electro-optic and terahertz diagnostics, in: S. Gold, G. Nusinovich (Eds.), AIP Conf. Proc. no. 1299, November 4, 2010, pp. 67–75. [18] M. Helle, D. Gordon, D. Kaganovich, A. Ting, Extending electro-optic detection to ultrashort electron beams, Phys. Rev. ST/AB 15 (2012) 052801-1– 052801-11. [19] J. Jackson, Classical Electrodynamics, John Wiley and Sons, New York, 1975. [20] R. Boyd, Nonlinear Optics, second ed., Academic Press, San Diego, 2003. [21] S. Casalbuoni, H. Schlarb, B. Schmidt, P. Schmuser, B. Steffen, A. Winter, Numerical studies on the electro-optic detection of femtosecond electron bunches, Phys. Rev. ST/AB 11 (2008) 072802-1–072802-18. [22] W. Faust, C. Henry, R. Eick, Dispersion in the nonlinear susceptibility of GaP near the reststrahl band, Phys. Rev. 173 (1968) 781–786. [23] Y. Berozashvili, S. Machavariani, A. Natsvlishvili, A. Chirakadze, Dispersion of the linear electro-optic coefficients and the non-linear susceptibility in GaP, J. Phys. D: Appl. Phys. 22 (1988) 682–686. [24] G. Gallot, D. Grischkowsky, Electro-optic detection of terahertz radiation, J. Opt. Soc. Am. B 16 (1999) 1204–1212. [25] S. Jamison, J. Shen, A. MacLeod, W. Gillespie, D. Jaroszynski, High-temporal-resolution, single-shot characterization of terahertz pulses, Opt. Lett. 28 (2003) 1710–1712. [26] G. Berden, S. Jamison, A. MacLeod, W. Gillespie, B. Redlich, A. van der Meer, Electro-optic technique with improved time resolution for real-time, nondestructive, single-shot measurements of femtosecond electron bunch profiles, Phys. Rev. Lett. 93 (2004) 114802-1–114802-4. [27] G. Berden, W. Gillespie, S. Jamison, E.-A. Knabbe, A. MacLeod, A. van der Meer, P. Phillips, H. Schlarb, B. Schmidt, P. Schmuser, B. Steffen, Benchmarking of electro-optic monitors for femtosecond electron bunches, Phys. Rev. Lett. 99 (2007) 164801-1–164801-4. [28] A. Debus, M. Bussmann, U. Schramm, R. Sauerbrey, C. Murphy, Z. Major, R. Horlein, L. Veisz, K. Schmid, J. Schreiber, K. Witte, S. Jamison, J. Gallacher, D. Jaroszynski, M. Kaluza, B. Hidding, S. Kiselev, R. Heathcote, P. Foster, D. Neely, E. Divall, C. Hooker, J. Smith, K. Ertel, A. Langley, P. Norreys, J. Collier, S. Karsch, Electron bunch length measurements from laser-accelerated electrons using single-shot THz time-domain interferometry, Phys. Rev. Lett. 104 (2010) 084802-1–084802-4. [29] T. Esirkepov, Exact charge conservation scheme for particle-in-cell simulation with an arbitrary form-factor, Comput. Phys. Commun. 135 (2001) 144– 153. [30] C.K. Birdsall, A.B. Langdon, Plasma Physics via Computer Simulation, Institute of Physics Publishing, Bristol and Philadelphia, 1991. [31] C. Decker, W. Mori, Group velocity of large amplitude electromagnetic waves in a plasma, Phys. Rev. Lett. 72 (1994) 490–493. [32] J.-P. Berenger, Three-dimensional perfectly matched layer for the absorption of electromagnetic waves, J. Comput. Phys. 127 (1996) 363–379. [33] E. Lindman, Free space boundary conditions for the time dependent wave equation, J. Comput. Phys. 18 (1975) 66–78. [34] A. Yariv, P. Yeh, Optical Waves in Crystals, John Wiley and Sons, 2003. [35] G. Agrawal, Nonlinear Fiber Optics. Optics and Photonics, second ed., Academic Press, San Diego, 1995. [36] C. Scoby, P. Musumeci, J. Moody, M. Gutierrez, Electro-optic sampling at 90° interaction geometry for time-of-arrival stamping of ultrafast relativistic electron diffraction, Phys. Rev. ST/AB 13 (2010) 022801-1–022801-7. [37] R. Fiorito, A. Shkvarunets, T. Watanabe, V. Yakimenko, Interference of diffraction and transition radiation and its application as a beam divergence diagnostic, Phys. Rev. ST/AB 9 (2006) 052802-1–052802-17. [38] A. Shkvarunets, R. Fiorito, Vector electromagnetic theory of transition and diffraction radiation with application to the measurement of longitudinal bunch size, Phys. Rev. ST/AB 11 (2008) 012801-1–012801-14.
Contents lists available at SciVerse ScienceDirect
Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp
Fully explicit nonlinear optics model in a particle-in-cell framework D.F. Gordon ⇑, M.H. Helle, J.R. Peñano Plasma Physics Division, Naval Research Laboratory, Washington, DC 20375, United States
a r t i c l e
i n f o
Article history: Received 10 August 2012 Received in revised form 25 January 2013 Accepted 5 May 2013 Available online 24 May 2013 Keywords: Particle-in-cell Nonlinear optics Finite difference time domain
a b s t r a c t A numerical technique which incorporates the nonlinear optics of anisotropic crystals into a particle-in-cell framework is described. The model is useful for simulating interactions between crystals, ultra-short laser pulses, intense relativistic electron bunches, plasmas, or any combination thereof. The frequency content of the incident and scattered radiation is limited only by the resolution of the spatial and temporal grid. A numerical stability analysis indicates that the Courant condition is more stringent than in the vacuum case. Numerical experiments are carried out illustrating the electro-optic effect, soliton propagation, and the generation of fields in a crystal by a relativistic electron bunch. Published by Elsevier Inc.
1. Introduction Nonlinear optics is often described in the frequency domain so that dispersive effects can be treated in a simple way. The convolution integrals associated with nonlinear effects are simplified by assuming that the radiation spectrum takes the form of a discrete set of narrow peaks. In the case of few cycle pulses, the convolution integrals cannot be reduced, and a time domain model becomes attractive. The finite-difference-time-domain (FDTD) technique is a well established method for solving the exact Maxwell equations in a bounded or periodic domain [1]. It may be used in connection with the particle-in-cell (PIC) technique to self-consistently solve for the motions of a large number of charged particles in an electromagnetic field [2]. This approach is often used to model fully nonlinear laser-plasma interactions and beam-plasma interactions. In this work, it is extended to account for nonlinear laser-crystal and beam-crystal interactions. This is accomplished by incorporating a model for bound particles into the PIC code turboWAVE [3]. The model is fully parallelized, and runs in up to three dimensions. The PIC aspect of turboWAVE has much in common with a number of other codes designed to model laser-plasma interactions [4–10]. The nonlinear optics aspect is related to a number of other codes designed to model ultrashort pulse propagation in a nonlinear medium [11–16]. Both types of codes can be divided into those that describe the radiation by means of a complex envelope, and those that are fully explicit, i.e., those that resolve the optical time scale. TurboWAVE supports both models, but in this work only the fully explicit model is used. As a result, there is no assumption about the frequency content of the radiation, and given the right model for the dielectric response, all nonlinear and dispersive effects are accounted for. Our model for the dielectric response generalizes the auxiliary equation technique described in Refs. [13,15]. In particular, bound charges in the dielectric are represented by effective particles, whose contribution to the four-current is computed using PIC techniques. The effective particles are subjected to forces arising from a superposition of macroscopic and microscopic fields. The macroscopic field is the usual electromagnetic field that is computed in any FDTD code, while the ⇑ Corresponding author. Tel.: +1 202 767 5036. E-mail address: [email protected] (D.F. Gordon). 0021-9991/$ - see front matter Published by Elsevier Inc. http://dx.doi.org/10.1016/j.jcp.2013.05.014
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402
389
microscopic field is a three dimensional electrostatic potential representing, e.g., an atomic binding potential. The resulting equation of motion is expanded as an anharmonic oscillator equation. By using multiple species of effective particles, each satisfying a different oscillator equation, the material response can be tailored to match that of real materials over a broad range of frequencies. Free charges are self-consistently incorporated into the model by linear superposition of the free and bound sources in the FDTD field solver.1 The resulting model is capable of modeling interactions among laser pulses, particle beams, plasmas, and dielectric crystals, in a fully dispersive, nonlinear, anisotropic, and kinetic way. To our knowledge, no other model has been implemented which simultaneously accounts for all these effects. One can argue that the auxiliary equation technique could be similarly extended, since, after all, one can recover any physics with the right ‘‘auxiliary equations.’’ The effective particle model described here has the advantage that all of the required physics falls out of a single unifying physical principle. In Ref. [17], the turboWAVE extensions described in this paper were introduced for the first time. In Ref. [18], the model was applied to a novel electro-optic diagnostic technique. In the present paper, the model is described in more detail, and a formal analysis of numerical stability is given. The model is extended to include third order nonlinearities, and the implementation is demonstrated via three numerical experiments. 2. Description of the model The numerical model described here solves the exact Maxwell equations:
@D @t @B rE¼ @t rD¼q
rH¼Jþ
ð1aÞ ð1bÞ ð1cÞ
rB¼0
ð1dÞ
Here, B ¼ l0 H þ M and D ¼ 0 E þ P, where P is the polarization and M is the magnetization. In the present work, we take M ¼ 0. The polarization can be expressed in terms of a spatially smoothed four-current due to bound charges, denoted by ðhgi; hjiÞ. This results in [19]
r B ¼ l0 ðJ þ hjiÞ þ
1 @E c2 @t
@B rE¼ @t r E ¼ ðq þ hgiÞ=0 rB¼0
ð2aÞ ð2bÞ ð2cÞ ð2dÞ
With this formulation, any field solver designed to account for free charges can be easily adapted to account for bound charges by making the substitution
ðq; JÞ ! ðq þ hgi; J þ hjiÞ
ð3Þ
Our model for ðhgi; hjiÞ is based on calculating the motions of effective particles that respond to the superposition of a microscopic binding potential and a macroscopic electric field E. Denoting the displacement of an effective particle from its equilibrium position by r, and expanding the microscopic potential in a Taylor series, results in the effective particle equation of motion 2
d ri dt
2
þ
X dr j q 2Cij þ ðX2 Þij r j þ aijk rj r k bri r j r j ¼ Ei m dt jk
ð4Þ
where the subscripts vary over Cartesian coordinates, q=m is the charge to mass ratio, and the v B force is neglected. The anisotropy of the medium is represented by the tensors C; X2 , and a. These are, respectively, the damping rate, the square of the resonant frequency, and the first anharmonic coefficient. At present, the second anharmonic coefficient, b, is assumed scalar. By means of a coordinate transformation, X2 can always be made diagonal. The basis vectors which induce this transformation will be denoted by ðeu ; ev ; ew Þ. The basis vectors used for the general calculation will be denoted ðex ; ey ; ez Þ. The crystallographic basis vectors, which generally coincide with ðeu ; ev ; ew Þ, are denoted (h1 0 0i, h0 1 0i, h0 0 1i). In general, the parameters in Eq. (4) depend not only on tensor indices, but also on an index p that identifies a particular particle. In other words, each particle may move in a unique potential well, and may therefore satisfy a unique equation of motion. The system of the particle together with its equation of motion is called an oscillator. In practical implementations, it is convenient to group all particles that satisfy the same equation of motion into an oscillator species with index s. Then the 1 The coupling of free and bound charges is through the macroscopic field only, i.e., collisions are neglected. As a result, the model is most useful when the free charges do not impinge on the dielectric.
390
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402
material parameters, q; m; X2 ; C; a, and b, depend on s, while only r is explicitly dependent on p. As discussed below, the macroscopic sources ðhgi; hjiÞ can be derived from the set rp using a variation on source deposition techniques developed for PIC codes. 3. Analytical dispersion relation In order to analyze the model described above, consider the locally averaged displacement hrs i associated with each oscillator species. In the case of a uniform medium, this is related to the spatially smoothed four-current by
ðhgi; hjiÞ ¼
X s
@ qs Ns r fhrs i; fhrs i @t
ð5Þ
where N s is the density of oscillators, and the tensor f is an oscillator strength. Assuming local phase coherence, hrs i satisfies the same equation as rp , which in the ðeu ; ev ; ew Þ basis, after linearizing, is
! @2 @ qs @Ai @/ 2 : þ 2 C þ X i ¼ þ hr iis is iis @t ms @t @ni @t2
ð6Þ
where i is the coordinate index, ni is the ith coordinate, A is the vector potential, and / is the scalar potential. Using the Coulomb gauge, the radiation is described by
c2 r2
! X q Ns @hris i @ 2 / @2 Ai ¼ fiis s þ 2 0 @t @t@ni @t s
ð7Þ
Consider any mode with a wave-vector that is collinear with one of the principal axes, ðeu ; ev ; ew Þ. Then, it can be shown that / ¼ 0, and the transverse components of A satisfy the dispersion relation 2
x
! X fiis x2ps 2 c2 k ¼ 0 1þ Dis ðxÞ s
ð8Þ
where x2ps ¼ N s q2s =0 ms , and
Dis ðxÞ ¼ X2iis 2ixCiis x2
ð9Þ
In the case of a medium comprised of a single, lossless, isotropic species with unit oscillator strength, the dispersion relation is 2
ðX2 x2 Þðx2 c2 k Þ þ x2 x2p ¼ 0
ð10Þ
The primary feature of this dispersion relation is a stop-band between the resonant frequency, X, and the cutoff frequency, qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X2 þ x2p . 4. Oscillator parameters Measured dielectric properties are usually given as susceptibilities in the frequency domain. Hence, in order to apply the time dependent model described above to real materials, it is necessary to connect the constants X2 ; C; f; a, and b, with the susceptibilities vðnÞ ðx1 ; x2 ; . . .Þ. In the ðeu ; ev ; ew Þ basis, the linear susceptibility has only the diagonal elements
vii ðxÞ ¼
X fiis x2ps Dis ðxÞ s
ð11Þ
Following Ref. [20], the first anharmonic coefficient is related to
Xq fs x2ps s ðaijk Þs ms Ds ðx1 ÞDs ðx2 ÞDs ðx3 Þ s
ð2Þ vijk ð x1 ; x2 ; x3 Þ ¼
vð2Þ by ð12Þ
where for simplicity, a cubic crystal is assumed. The nonlinear susceptibility is related to the electro-optic coefficient rijk via ð2Þ vijk ðx þ dx; x; dxÞ ¼
1X il ðxÞrlmk ðx; dxÞmj ðxÞ 2 lm
ð13Þ
Here, the factor of 1=2 accounts for the degeneracy of sum and difference generation when dx ! 0, and is the relative permittivity. Finally, the nonlinear refractive index is given by
n2 ðxÞ ¼
x2ps 3 g0 X q2s bs 4 n0 ðxÞ2 s m2s Ds ðxÞDs ðxÞ3
ð14Þ
391
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402 Table 1 Oscillator parameters for gallium phosphide. Parameter
Lattice oscillator
Electronic oscillator
Oscillator strength, f Resonance frequency, X
1
1
6:90 1013 rad/s
Damping frequency, c
6:25 1010 rad/s
6:38 1015 rad=s 0
Plasma frequency, xp
9:27 1013 rad/s
1:78 1016 rad/s
Anharmonic coefficient, a123
3:28 1037 m1 s2
4:1 1041 m1 s2
Fig. 1. Nonlinear Lorentz model for GaP. Panel (a) displays RðÞ, where the solid curve is the two-oscillator (2O) model and the points are from the ð2Þ expressions in Ref. [21]. Panel (b) is similar to (a) except IðÞ is displayed. Panel (c) displays v123 ðxL dx; xL ; dxÞ, where xL is a probe laser frequency, given by 2pc=xL ¼ 0:6328 lm. The solid curve is the 2O model, and the dashed curve is the expression in Ref. [22]. The triangles indicate the infra-red frequencies where the data points in [22] are located. Panel (d) displays r 123 ðx; dxÞ ¼ r 41 ðx; dxÞ, where dx ! 0. The curve is the 2O model, and the points are the data from Ref. [23]. The 2O model cannot be made consistent with both Refs. [22,23].
where n0 is the linear refractive index and g0 is the impedance of free space (see Appendix). Note that the charge to mass ratio qs =ms amounts to an extraneous free parameter that could have been absorbed into qs N s ; as , and bs . In this work, the electronic charge to mass ratio is always assumed. As an example, consider the cubic crystal gallium phosphide (GaP), which is useful for electro-optic sensing applications [24–28]. Expressions fitting the measured linear dispersion have already been given in Ref. [21]. There, two separate expressions were used for the optical and THz regimes. A two-oscillator model based on the parameters given in Table 1, agrees with both expressions, as illustrated in Fig. 1(a) and (b). Once the linear dispersion is determined, the only remaining free parameters are the anharmonic coefficients associated with each oscillator. In order to choose these parameters, an attempt was made to simultaneously fit data from Refs. [22,23]. ð2Þ In Ref. [22], vijk ðxL dx; xL ; dxÞ is reported, where xL is the frequency of a probing helium–neon laser, and dx is one of several infra-red frequencies produced by a multi-line H2þO2 laser. In Ref. [23], rijk ðx; dxÞ is measured, where dx is in the radio frequency (RF) range, and x takes several values in the optical range. It was found that the two-oscillator model cannot be made to agree with both experiments. The parameters of Table 1 are chosen to obtain agreement with the measured rijk ðx; dxÞ. The resulting dispersion is displayed in Fig. 1(c) and (d). In order to make the two-oscillator model agree with vð2Þ ðxL dx; xL ; dxÞ, a much larger absolute value of the lattice anharmonic coefficient must be used. This leads, roughly, to the oscillator parameters given in Ref. [17], although for the simulations presented there and in Ref. [18], the lattice oscillator was taken to be linear.
5. Numerical technique 5.1. Oscillator equation and source deposition 1
1
Most electromagnetic PIC codes use the leap-frog technique. More specifically, Enþ1 is computed using En ; Bnþ2 and Jnþ2 , 3 1 where n is the time level. Then, Bnþ2 is computed using Bnþ2 and Enþ1 . In some cases a Poisson solver is used to refine
392
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402
Fig. 2. Charge Deposition. (a) Geometry of a grid cell in relation to the points on the Yee mesh. The three dimensional arrows point to the mesh points where the corresponding electric field components are known. (b) Quadratic Weighting. The bound charges are represented by charge clouds (solid curve) whose equilibrium position (dashed curve) is the centroid of the cell volume. For an electron, the charge in a cell is increased (decreased) by the area of the intersection of the blue (red) area with the cell. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Enþ1 using qnþ1 . As discussed above, any such field solver can be used in the presence of an oscillator population by making the substitution ðq; JÞ ! ðq þ hgi; J þ hjiÞ. This implies that hgi is known at integral time levels, and hji is known at half-integral time levels. 1 To compute hginþ1 and hjinþ2 , the displacements rnp and rnþ1 are needed. Let T represent the matrix of transformation of a p vector from the ðeu ; ev ; ew Þ basis to the ðex ; ey ; ez Þ basis. If all quantities characterizing the oscillators are stored in the ðeu ; ev ; ew Þ basis, the displacements can be updated using2
rnþ1 2rni þ rn1 r nþ1 r ni F ni 2 n i i i þ 2 C þ ð X Þ r ¼ ii i ii Dt m Dt 2
ð15Þ
where Dt is the time step and
Fn ¼ qT1 En marn rn þ mbðrn rn Þrn
ð16Þ rnp
The macroscopic four-current ðhgi; hjiÞ associated with the microscopic orbits can be computed using PIC techniques. These techniques regard each particle as a ‘‘cloud’’ with dimensions on the order of a cell size. Knowing the orbit of each cloud’s centroid allows one to compute the amount of charge in each cell as a function of time, as well as the current passing through each cell wall. Various schemes have been developed to facilitate this calculation [2,29,30]. In the case of bound charges, these schemes have to be applied with care in order to minimize round-off errors associated with very small displacements. To deposit the sources due to bound charges, place one oscillator, per species, in each grid cell.3 The cell geometry is shown in Fig. 2(a). The deposition of charge, in one dimension, is illustrated in Fig. 2(b). The oscillator is represented by two oppositely charged density distributions. The equilibrium distribution, shown as a dashed line, is immobile. The displaced distribution, shown as a solid line, executes the orbit described by Eq. (15). The net charge in a cell is the integral over the cell volume of the difference between the distributions. To minimize round-off error, this integral should not be carried out by re-using a PIC routine that deposits each distribution separately. Instead, expressions should be written for the net charge in a cell, and expanded to lowest order in the displacement. Then, for the oscillator in cell ði; j; kÞ, the net charge density is distributed in the surrounding cells as follows:
N s qs 2 Dx n N s qs ¼ ey Tfrp 2 Dy N s qs ¼ ez Tfrnp 2 Dz
hgini1;j;k ¼ ex Tfrnp
ð17aÞ
hgini;j1;k
ð17bÞ
hgini;j;k1
ð17cÞ
Here, the cell dimension is Dx Dy Dz. The current density through the surrounding cell walls follows directly from charge conservation:
2 3
Forward differencing of the damping term improves numerical stability, see Eq. (23). In the present work, we assume local phase coherence. By using many effective particles per cell, this assumption could be removed.
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402
Dx nþ1 n ex hjii12;j;k ¼ hginþ1 i1;j;k hgii1;j;k Dt 2 Dy nþ12 nþ1 n ey hjii;j1;k ¼ hgii;j1;k hgii;j1;k Dt 2 Dz nþ12 nþ1 n ez hjii;j;k1 ¼ hgii;j;k1 hgii;j;k1 Dt 2
393
ð18aÞ ð18bÞ ð18cÞ
Note that the above expressions are particular to the choice of quadratic weighting. This choice provides adequate smoothing, while also allowing for abrupt vacuum-dielectric transitions. A useful fact is that the current density due to an effective particle, through any cell wall enclosing that particle, is the normal component of 12 N s qs v p , where
vp
rnþ1 rnp p ¼ Tf Dt
! ð19Þ
is the velocity of the effective particle in the ðex ; ey ; ez Þ basis. 5.2. Linear stability and accuracy In an ordinary PIC code, the primary stability criterion is the Courant condition, which in one dimension reads cDt < Dz, where Dt is the time step and Dz is the grid spacing. It is clear that in the presence of a dispersionless dielectric, this would be modified to read cDt=n < Dz, where n is the refractive index. In the case of a fully dispersive model such as the one described above, the modification of the Courant condition is more subtle. Expressing Eqs. (6) and (7) in finite difference form, assuming propagation along a principal axis, and inserting exponential forms for Ai and hris i, gives the numerical dispersion relation
cos xDt ¼ VðkÞ þ P i ðxÞ
ð20Þ
where
VðkÞ ¼ 1 2c2 Dt 2
2 1 k Dx 2 x Dx2
sin
þ
2 1 k Dy 2 y Dy 2
sin
þ
2 1 k Dz 2 z Dz 2
sin
! ð21Þ
and
P i ðxÞ ¼
X s
x2ps Dt2 sin2 12 xDt 2 1 ðX Þiis Dt ð1 þ Ciis DtÞ sin 12 xDt 2i Ciis Dt sin xDt 4 1 f 2 iis
2
2
ð22Þ
Note that when the density of oscillators vanishes, P i ! 0 and the well known vacuum equation is recovered. The effect of P i can be estimated by assuming X2 Dt 2 1 and kx Dx ky Dy kz Dz xDt p. These assumptions are motivated by the fact that the resonance frequency should be well resolved, and by the expectation that numerical instabilities are most severe at the Nyquist frequency. Under these assumptions, the stability criterion becomes
c2 c2 c2 Dt K þ 2þ 2 2 Dx Dy Dz
1=2
2 1 X fiis x2ps Dt 1 4 s 1 þ Ciis Dt
!1=2 ð23Þ
This is more stringent than the usual Courant condition. Numerical experiments show that Eq. (23) is accurate in practice. In order to evaluate the accuracy and stability of the differencing scheme more precisely, the numerical dispersion relation is evaluated numerically for the case of a single, lossless oscillator species. Since the dispersion relation scales, we normalize all frequencies to the resonant frequency, X, taking xp ¼ 2X; Dz ¼ 0:5c=X, and Dx ¼ Dy ! 1. The results are shown in Fig. 3 for three choices of the time step. The case Dt ¼ 0:2=X is shown in panel (a), the case Dt ¼ 0:4=X is shown in panel (b), and the case Dt ¼ 0:48=X is shown in panel (c). The analytical dispersion relation (10) is overlayed for comparison. As usual, the best accuracy is obtained for small values of x and k. As cDt=Dz increases, the accuracy improves, until numerical instability sets in. This behavior is similar to the vacuum case. However, unlike the vacuum case, instability occurs for cDt=Dz < 1, as shown in panel (c). Better accuracy may be obtained by decreasing both Dz and Dt proportionately. 5.3. Nonlinear stability In the case of the fully nonlinear system, a rigorous stability analysis is difficult. However, if the oscillator equation is viewed as describing an effective particle in a potential well, the stability criterion may be viewed as the requirement that the particle never enter any region of runaway acceleration. Such a region exists if b > 0, for in this case there is a divergent repulsive force as r ! 1. Instability is also possible if b ¼ 0 and a – 0. Assuming there are no metastable regions, the stability criterion can be written
FNL r < mX2 r r
ð24Þ
394
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402
Fig. 3. Numerical dispersion for (a) cDt=Dz ¼ 0:4, (b) cDt=Dz ¼ 0:8, and (c) cDt=Dz ¼ 0:96. Solid curves are the analytical dispersion relation (10), solid circles are Rðx=XÞ, and open circles connected by a thin line are Iðx=XÞ. In all cases, xp =X ¼ 2, and Dz ¼ 0:5c=X. Panel (c) shows that numerical instability may occur for cDt=Dz < 1.
where
FNL ¼ mbðr rÞr marr
ð25Þ
If one is interested primarily in second order effects, the system can be stabilized by taking b < 0. The magnitude of b is chosen such that the order of magnitude of the first, second, and third order forces are the same at some critical point. This leads to
b
a2
X2
ð26Þ
where a and X are typical values of the elements of a and X. If third order effects are important physically, a stabilizing fifth order nonlinearity can be easily introduced. In this case,
FNL ¼ mdðr rÞ2 r þ mbðr rÞr marr
ð27Þ
where d is the fifth order anharmonic coefficient. Demanding that the first, third, and fifth order forces be equal at some critical point gives 2
d
b
X2
ð28Þ
The ultimate nonlinear stabilization technique would utilize physical effects, such as ionization in the case of the electronic oscillator, or material damage in the case of the lattice oscillator. However, accounting for these processes is beyond the scope of the present work.
5.4. False resonance technique for accelerating computation The addition of one effective particle per grid cell to a vacuum FDTD field solver increases the computational workload by only about 25%. However, it sometimes happens that the accuracy and stability criteria, X2 Dt 2 1 and x2p Dt 2 1, force one to use a very large number of time steps to simulate a given process. This limitation can be mitigated if the frequency range of interest is well below X. In such cases, it is possible to artificially shift X to a lower frequency while preserving the dispersion relation in the region of interest. In this way, the ‘‘true resonance’’ is replaced by a ‘‘false resonance’’ which is easier to resolve, but leads to the same physics. As an example of a false resonance model, consider again GaP, which has the true resonance parameters displayed in Table 1. When modeling an electron bunch with a characteristic time of, say, 100 fs, the frequency range of interest is well below the electronic resonance, which has 2p=X 1 fs. One may then use a false resonance model such as the one described by Table 2. The linear and nonlinear dispersion curves associated with this model are shown in Fig. 4. In the THz range, the dispersion relation remains accurate. Note that the dispersion in this range could be fine tuned by adding additional false resonances. In Ref. [17], computation was accelerated by sub-cycling the oscillator equation with respect to a super-cycle time step used for all other equations. In order to make this stable, one must choose the super-cycle time step to be an integral number of periods of the under-resolved frequency. Apart from placing needless constraints on the simulation parameters, this scheme was found to give results that do not satisfactorily conserve energy, and which eventually diverge. In the case of the false resonance approach, one can prove linear stability, and as is shown below, energy is accurately conserved.
395
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402 Table 2 False resonance model for gallium phosphide. Parameter
Lattice oscillator
False resonance
Oscillator strength, f Resonance frequency, X
1
1
6:90 1013 rad/s
Damping frequency, C
6:25 1010 rad/s
8:91 1014 rad/s 0
Plasma frequency, xp
9:27 1013 rad/s
2:46 1015 rad/s
Anharmonic coefficient, a123
3:28 1037 m1 s2
1:23 1040 m1 s2
Fig. 4. False resonance model for GaP. Panel (a) displays RðÞ, where the solid curve is the false resonance model and the points are from the expressions in ð2Þ Ref. [21]. Panel (b) displays v123 ð2dx; dx; dxÞ, where the curve is the false resonance model and the points are the true resonance model (Table 1).
5.5. Boundary conditions and moving window In simulations of laser-plasma interactions, the moving window technique [31] is often used. In this technique, the mesh points are shifted by one cell in between time levels such that the computational region, or ‘‘window,’’ moves at the speed of light. Boundary conditions are greatly simplified by the requirements of causality. When the plasma density is low, the laser pulse group velocity is nearly the speed of light, and the pulse stays in the window for a long time. Hence, the window only has to be as long as the laser pulse, even if the interaction region is much longer. When the group velocity is significantly smaller than the speed of light, the pulse quickly falls behind a light frame window, and one might as well work in the lab frame. If the window speed is the group velocity, a fully dispersive model might allow precursor signals to reach the front of the window. Moreover, signals reflected from the back of the window might work their way toward the region of interest. One solution is to use perfectly matched layers (PML) [32] to absorb waves near the boundaries of the window. In order for this to be effective, it is necessary to use a large number of layers due to the fact that the moving window involves shifting the field data by one cell between time levels. In order that this shift should weakly perturb the solution, the PML conductivity has to change very gradually from one cell to the next. In the case of lab frame simulations, another option is to use Lindman’s absorbing boundary condition [33] in the longitudinal direction, and periodic boundary conditions transversely. If a Poisson equation has to be solved, the solution in the box can be matched to a decaying solution outside the box [30]. TurboWAVE supports all the options mentioned above. A simulation that uses Lindman’s boundary condition in the lab frame is discussed in Section 6.1. A simulation that uses PML absorbers with a moving window that moves at the group velocity is discussed in Section 6.2. A simulation that uses a moving window that moves at the speed of light is discussed in Section 6.3. 6. Numerical experiments The model described above has been implemented as a module within the turboWAVE [3] framework. In this section, three numerical experiments are presented. First, the electro-optic effect in GaP is demonstrated by comparing the numerical change in polarization with that predicted by the low frequency theory. Second, soliton propagation in fused silica is demonstrated. Finally, to take full advantage of the PIC framework, a relativistic electron bunch is passed near a GaP crystal, and the induced fields are examined.
396
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402 Table 3 Parameters for simulation of electro-optic effect. The fields are measured inside the dielectric. The THz field is a half-cycle pulse with base-to-base width sH . Parameter
Symbol
Value
Time step Space step Cells
cDt Dz Nz
0.012 lm 0.017 lm
Steps
Nt
THz field THz half wave Laser field Laser FWHM Laser wavelength Crystal length Crystal orientation
Ey
217 68 kV/cm 2.2 ps 18 kV/cm 220 fs 0.81 lm 120 lm Eq. 30
sH jE x j
sL k0 L T
213
6.1. Electro-optic effect In this numerical experiment, an x-polarized laser pulse is copropagated with a y-polarized THz half-wave in a GaP crystal, and the change in the polarization state of the laser pulse is measured. The length of the THz half-wave is chosen to be long enough so that the standard picture of the electro-optic effect applies. In this picture, the change in the impermeability4 tensor due to a DC electric field is
Dgij ¼ r ijk Ek
ð29Þ
For any Zinc–Blende type crystal, r ijk ¼ r 123 ¼ r41 if all three indices are distinct, and r ijk ¼ 0 otherwise. The crystal orientation is represented by the operator T, which can be viewed as a sequence of rotations that re-orients the crystal starting with the ðeu ; ev ; ew Þ and ðex ; ey ; ez Þ bases aligned. It turns out that the electro-optic effect is maximized if [34]
T ¼ Rz ðp=2ÞRx ðp=2ÞRz ðp=4Þ
ð30Þ
with the laser field polarized in the x-direction and the DC field polarized in the y-direction. Here, Ri is an active right-handed rotation about the ith coordinate axis. As an example, if E and Dg are given in the ðex ; ey ; ez Þ basis, and r is given in the ðeu ; ev ; ew Þ basis, then Dg ¼ TrT1 ET1 . To calculate the change in polarization due to the electro-optic effect, a principal axis transformation is usually made [34]. Alternatively, one may make direct use of the solution of the slowly varying envelope equation,
EðzÞ ¼ exp
ipDz Eðz ¼ 0Þ nk0
ð31Þ
where E is the complex electric field envelope, k0 is the free space laser wavelength, n is the refractive index of the unperturbed crystal, and D is the change in permittivity induced by the DC field. This equation remains valid even when D is a tensor and E is a vector, and is easily evaluated in a software environment that supports the exponential of a matrix. We now compare the prediction of Eq. (31) with the results of a turboWAVE simulation. In the simulation, an x-polarized laser pulse and a y-polarized THz half-wave are incident on a GaP crystal oriented as in Eq. (30). The simulation parameters are given in Tables 1–3. The lattice nonlinearity was stabilized with b ¼ a2 =X2 . The results from the simulation are shown in Fig. 5. The data are evaluated 60 lm into the crystal to avoid interference from reflections generated at either vacuumcrystal interface. Since the turboWAVE model uses a universal electric field, the optical and THz frequency components have to be extracted by means of a Fourier filter. The envelope of the optical wave is constructed by applying a hard-edged bandpass filter centered at the optical frequency, þx0 . The result is then shifted in frequency by x0 , multiplied by 2 (to account for the energy in the negative frequencies) and the inverse Fourier transform is applied. Based on Fig. 5, the polarization ratio at z ¼ 60 lm, and at times prior to the generation of internal reflections, is jE y =E x j2 0:0018. Applying the analytical formula with the inputs Ey ¼ 68 kV/cm, n ¼ 3:16, and r123 ¼ 0:89 pm/V gives jE y =E x j2 0:0019. Thus, the simulation model closely agrees with the theoretical prediction. An interesting feature of the simulation data is that the second y-polarized optical pulse is larger than the first. This appears to happen because the reflected x-polarized optical pulse continues to be rotated by the reflected y-polarized THz pulse. This leads to a coherent reinforcement of the reflected y-polarized optical pulse. 6.2. Soliton propagation In this numerical experiment, a 20 fs, 2.1 lm pulse is propagated in fused silica in both the linear and soliton regimes. The Sellmeier formula for fused silica can be matched exactly to a three oscillator model, the parameters of which are displayed 4
Impermeability is defined by g ¼ 1 , where
is the relative permittivity.
397
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402
Fig. 5. Simulated electric field vs. time at z ¼ 60 lm from the entrance plane of the crystal. (a) y-polarized THz field (low-pass Fourier filter) (b) x-polarized optical field envelope (band-pass Fourier filter with shift) (c) y-polarized optical field envelope (band-pass Fourier filter with shift). The incident optical field is a single, purely x-polarized pulse, and the incident THz field is a single, purely y-polarized half-wave. The additional optical pulses, and the late-time distortion of the THz field, are due to reflections.
Table 4 Oscillator parameters for fused silica. Parameter
Osc. 1
Osc. 2
Osc. 3
Oscillator strength, f Resonance frequency, X
1
1
1
2:35 1014 0
1:63 1016 0
2:78 1016 0
1:79 1014 0
1:06 1016 0
2:30 1016
rad/s rad/s
4:72 1054
m2 s2
Damping frequency, C Plasma frequency, xp Anharmonic coefficient, b
Unit rad/s
in Table 4. The second anharmonic coefficient is chosen to give a nonlinear refractive index of n2 ¼ 3 1016 cm2 =W. The choice to take b – 0 only in the third oscillator has no effect other than on the frequency dependence of n2 . The amplitude envelope of the fundamental soliton has the form sechðsÞ, where s ¼ ðt z=v g Þ=T 0 ; v g is the group velocity, and T 0 is the pulse width parameter. In simulations, it is convenient to have a function whose support is strictly in a finite region. Furthermore, in a model of the type used here, the initial conditions depend on z, and the pulse starts in vacuum. Hence, the appropriate initial condition is
Ex ðt ¼ 0; zÞ ¼ 1=4 jE 0 jsechðf=cT 0 Þ sinðx0 f=cÞCðf=cT 0 Þ
ð32Þ
where Ex is the field that is initialized in the vacuum region, jE 0 j is the desired peak field after propagating into the dielectric, f ¼ z z0 ; z0 fixes the spatial origin, and CðxÞ is a suitable clipping function. The clipping function used here is
CðxÞ ¼
8 > < > :
jxj 6 4; 4 < jxj < 6;
1 cos2 ½ðjxj 4Þp=4
ð33Þ
jxj P 6; 0
It should be noted that when the pulse enters the dielectric there is a small reflection that is not accounted for in Eq. (32). The parameters of the pulse, and the numerical parameters associated with the simulation, are displayed in Table 5. The intensity and electric field are related by I0 ¼ 1=2 jE 0 j2 =2g0 . The free space wavelength is k0 ¼ 2pc=x0 . The group velocity dispersion (GVD) coefficient is defined by cb2 ¼ @ 2 ð1=2 xÞ=@t 2 jx¼x0 . The dispersion length, LD ¼ T 20 =jb2 j, measures the propagation length needed to observe group velocity dispersion (GVD) effects. The nonlinear length, LNL ¼ c=x0 n2 I0 , measures the propagation length needed to observe nonlinear pulse distortions. The conditions b2 < 0 and LD ¼ LNL result in stable propagation of the pulse without distortion [20,35]. This is illustrated in Fig. 6. Comparison of panels (a) and (b) shows that when I0 ¼ 0:4 TW/cm2, the pulse maintains its original form after propagating nearly 3LD . Panel (c) shows that when I0 is reduced by a factor of 100, the pulse broadens significantly over the same propagation distance. 6.3. Electron bunch passing near a crystal The electro-optic effect can be used to measure the temporal structure of relativistic electron bunches [18,26–28]. The general principle is that the bunch self-fields scatter a probing optical pulse into a new polarization and frequency [18,21,24]. In order to make use of such a principle, the form of the bunch fields in the crystal has to be understood. In the simplest interpretation, the self-fields of the bunch are treated as the equivalent of a planar terahertz half-wave incident on the crystal. In reality, the picture may be complicated by effects such as Cherenkov-like wakes [17,36] and diffraction radiation [37,38]. In order to fully characterize these effects, we carry out simulations of the propagation of an electron
398
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402 Table 5 Parameters for simulation of soliton propagation. Parameter
Symbol
Value
Time step Space step Cells
cDt Dz Nz
0.014 lm 0.034 lm
Steps
Nt
Window speed PML thickness Intensity Pulse width Free space wavelength Dispersion length Nonlinear length GVD coefficient
vg
4 105 0:68c 8.7 lm 0:4 TW=cm2 20 fs 2.1 lm 2.8 mm 2.8 mm
– I0 T0 k0 LD LNL b2
212
1:42 1025 s2/m
Fig. 6. Simulation of 1D soliton propagation in fused silica. Plots show 1=2 E2x =g0 vs. v g t z, where v g t is held fixed in each plot, and v g is the velocity of the moving window. Panel (a) shows the initial pulse, (b) shows the pulse after 8 mm propagation, and (c) shows the result from (b) in the case where the initial intensity is I0 =100.
bunch near a GaP crystal. The parameters are given in Table 6, and the geometry is shown in Fig. 7. The simulations differ from the one carried out in Ref. [17] in that the bunch charge is varied until significant nonlinear distortions are observed in the bunch fields, and in that the false resonance technique is used in place of sub-cycling the oscillator equation. The false resonance oscillator parameters are the same as those given in Table 2, except that a stabilizing third order nonlinearity is added to each oscillator. The electron bunch is initialized in the presence of an immobile neutralizing charge, and is abruptly accelerated at the beginning of the simulation. This generates spurious transient fields that have to be attenuated. For this reason, the bunch is allowed to propagate for t0 ¼ 15 picoseconds before encountering the GaP crystal, whose plane of incidence is located at z ¼ 4200 lm. Convergence to the proper steady state is facilitated by placing PML absorbers at the transverse boundaries of the simulation region. The spurious self-force exerted on the bunch during the transient phase is small in the present context. Fig. 8 displays the electric field after the bunch has moved about 400 lm beyond the plane of incidence of the GaP crystal. More specifically, Fig. 8(a) and (b) show Ey ðt t0 ; x ¼ 0; y; zÞ for the cases Q ¼ 0:28 pC and Q ¼ 2:8 pC, where Q is the bunch
Table 6 Parameters for simulation of fields induced by an electron bunch. Parameter
Symbol
Value
Time step Space step Cells
cDt Dx ¼ Dy ¼ Dz Nx Ny Nz
0.168 lm 0.672 lm
Steps
Nt
Window speed Bunch diameter Bunch length Electron energy Crystal size Crystal orientation
v rb
3 104 c 8 lm 80 fs 250 MeV 605 360 1 lm3 Eq. (30)
sb ðc 1Þmc2 Lx Ly Lz T
210 210 210
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402
399
Fig. 7. Geometry of the simulation of an electron bunch passing near an electro-optic crystal. The electron bunch is loaded into a Gaussian weighted ellipsoid (green) and propagated in the z-direction adjacent to a crystal whose internal structure is oriented as in Eq. (30) (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this article.) .
charge. Fig. 8(c) and (d) show the corresponding line-outs at y ¼ 0. The primary features in Fig. 8(a) and (b) are the vertical phase fronts located in the range 400 < z ct < 300 lm, and the diagonal phase front that appears at all z ct. As discussed previously [17], these can be roughly associated with coherent transition radiation (CTR) and Cherenkov radiation, respectively. In electro-optic sensing applications, the CTR-like feature is supposed to retain the form of the vacuum bunch fields. However, as is well known, dispersive effects lead to distortions of the type shown in the figure. When the bunch charge is high enough to excite large nonlinear polarization currents in the material, the fields induced in the crystal may be further distorted. Such a case is shown in Fig. 8(b) and (d). For the chosen crystal orientation, the nonlinear polarization is ð2Þ
ð2Þ
Pð2Þ ¼ v123 E2y ex þ 2v123 Ex Ey ey
ð34Þ
Fig. 8. Simulation of fields induced by an electron bunch with (a,c) Q ¼ 0:28 pC and (b,d) Q ¼ 2:8 pC. The low charge case is essentially the linear result, while the high charge case exhibits nonlinear distortions. The three dimensional data is evaluated at x ¼ 0 in order to present a two dimensional plot. The green oval represents the location of the electron bunch, which propagates from left to right. The boundary between the crystal and vacuum is located at y ¼ 200 lm. The plane of incidence is located at z ct ¼ 435 lm. The line-outs in (c) and (d) are evaluated at y ¼ 0.
400
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402
Fig. 9. Conservation of energy during simulated bunch-crystal interaction.
Hence, a y-polarized THz field induces an x-polarized field with second harmonic components. This field then mixes with the original y-polarized field to produce a distorted y-polarized field. This distortion can be clearly seen by comparing panels (c) and (d) of Fig. 8. Qualitatively, the effect is to produce narrower features in the waveform. In an electro-optic sensing application, this could adversely affect the fidelity of the diagnostic. Finally, consider the flow of energy in the crystal during the simulation. For the specific form of the dielectric response considered in this paper, Poynting’s theorem can be written as
Z
dn S þ
0
@X
d 2 dt
Z
3
d rðE2 þ c2 B2 Þ þ
X
XdU p þ Lp ¼ 0 dt p2X
ð35Þ
where the first term is the energy radiated through a surface @ X, the second is the rate of change of field energy in the volume X, and the third is the rate of change of the energy of all the effective particles in X. Here, the effective particle energy is
Up ¼
2 1 dr p mp þ qp /p ðr p Þ 2 dt
ð36Þ
where mp is the mass, rp is the displacement, qp is the charge, and /p ðr p Þ is the microscopic electrostatic potential of the displaced particle. The rate of frictional damping for an effective particle is
Lp ¼ 2Cp mp
2 dr p dt
ð37Þ
Fig. 9 displays the three terms in Poynting’s theorem, integrated over time, along with the sum. The ‘‘radiated’’ curve is negative, corresponding to the fact that the bunch fields propagate into the crystal. The steep initial slope corresponds to the generation of the CTR-like feature, while the more gradual slope corresponds to the generation of the Cherenkov-like feature. The ‘‘field’’ and ‘‘material’’ curves show that half the incident energy goes to the fields, while the other half goes to the effective particles. The ‘‘sum’’ curve shows that energy is conserved, to the extent that it vanishes. One notices a slight discrepancy at the moment the bunch fields first enter the crystal. This may be related to the numerical process of generating a reflected wave.
7. Conclusions The physics of nonlinear crystals excited by laser pulses, plasmas, particle beams, or any combination of these, can be numerically simulated using an effective particle model based on particle-in-cell techniques. In this model, bound charges respond to a superposition of microscopic and macroscopic fields, where the former are chosen to satisfy the known dispersion characteristics of the material, and the latter are self-consistently computed using the usual FDTD technique. The resulting numerical model is extremely flexible and makes very few assumptions. It is, however, computationally expensive. Nevertheless, it is now possible to carry out three dimensional simulations of electro-optic sensing, where the bunch fields induced in an electro-optic crystal are computed with unprecedented detail. Other applications of this model might include modeling the nonlinear optics of nano-composites, or frequency mixing with ultra-short pulses.
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402
401
Acknowledgments This work was supported by NRL Base Funds and the Department of Energy. We acknowledge useful conversations with B. Hafizi and A. Ting. Appendix A In this appendix, the relationship between the nonlinear refractive index, n2 , and the second anharmonic coefficient, b, is derived. For simplicity, all quantities are reduced to scalars. Consider a perturbation expansion of the anharmonic oscillator equation. Suppressing the species index, and expressing all quantities in frequency space, the first order equation is
DðxÞr ð1Þ ¼
q EðxÞ m
ð38Þ
and the third order equation is
DðxÞr ð3Þ ¼ b r ð1Þ r ð1Þ r ð1Þ
ð39Þ
Here, the asterisk denotes convolution over frequency, and the electric field is in the form
EðxÞ ¼
E0 ½dðx x0 Þ þ dðx þ x0 Þ
2
ð40Þ
Inserting r ð1Þ into the third order equation gives
3 qE0 FðxÞ rð3Þ ¼ b 2m DðxÞ
ð41Þ
where
FðxÞ ¼
dðx 3x0 Þ 3dðx x0 Þ 3dðx þ x0 Þ dðx þ 3x0 Þ þ 2 þ þ D 3 ð x0 Þ D ðx0 ÞDðx0 Þ D2 ðx0 ÞDðx0 Þ D3 ðx0 Þ
ð42Þ
Restoring species indices, the nonlinear polarization is
X Ns qs rð3Þs
Pð3Þ ¼
ð43Þ
s
This quantity must now be connected to the refractive index. The refractive index is defined by n ¼ 1=2 , with the relative permittivity. Using D ¼ 0 E, this can be written as
n¼
0 E þ P 0 E
1=2 ð44Þ
For a weakly nonlinear interaction, P ð3Þ Pð1Þ , and
n n0 þ
Pð3Þ 2n0 0 E
ð45Þ
where n0 is the linear index, given by n20 ¼ 1 þ P ð1Þ =0 E. The nonlinear index, n2 , is conventionally defined by
n ¼ n0 þ n2 I0
ð46Þ
where I0 is the cycle averaged intensity. Choosing the Fourier transform convention such that E0 is the peak value of the electric field, one has the relationship I0 ¼ n0 E20 =2g0 . Finally, equating Eqs. (45) and (46), and evaluating at x0 , gives
n2 ðx0 Þ ¼
x2ps 3 g0 X q2s bs 2 2 4 n0 ðx0 Þ s ms Ds ðx0 ÞDs ðx0 Þ3
ð47Þ
References [1] K. Yee, Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media, IEEE Trans. Antennas Propag. 14 (1966) 302. [2] J.P. Verboncoeur, particle simulation of plasmas: review and advances, Plasma Phys. Control. Fusion 47 (2005) A231–A260. [3] D. Gordon, Improved ponderomotive guiding center algorithm, IEEE Trans. Plasma Sci. 35 (2007) 1486–1488. [4] D.W. Forslund, J.M. Kindel, E.L. Lindman, Plasma simulation studies of stimulated scattering processes in laser-irradiated plasmas, Phys. Fluids 18 (1975) 1017–1030. [5] C.D. Decker, W.B. Mori, T. Katsouleas, Particle-in-cell simulations of raman forward scattering from short-pulse high-intensity lasers, Phys. Rev. E 50 (1994) 3338. [6] P. Mora, T. Antonsen, Kinetic modeling of intense, short laser pulses propagating in tenuous plasmas, Phys. Plasmas 4 (1997) 217.
402
D.F. Gordon et al. / Journal of Computational Physics 250 (2013) 388–402
[7] A. Pukhov, J. Meyer-ter-Vehn, Relativistic laser-plasma interaction by multi-dimensional particle-in-cell simulations, Phys. Plasmas 5 (1998) 1880– 1886. [8] R. Fonseca, L. Silva, F. Tsung, V. Decyk, W. Lu, C. Ren, W. Mori, S. Deng, S. Lee, T. Katsouleas, J. Adam, OSIRIS: a three dimensional, fully relativistic particle in cell code for modeling plasma based accelerators, Lect. Notes Comput. Sci. 2331 (2002) 342–351. [9] C. Nieter, J. Cary, VORPAL: a versatile plasma simulation code, J. Comput. Phys. 196 (2004) 448–473. [10] C. Huang, V. Decyk, C. Ren, M. Zhou, W. Lu, W. Mori, J. Cooley, T. Antonsen, T. Katsouleas, QUICKPIC: a highly efficient particle-in-cell code for modeling wakefield acceleration in plasmas, J. Comput. Phys. 217 (2006) 658–679. [11] M. Mlejnek, E. Wright, J. Moloney, Dynamic spatial replenishment of femtosecond pulses propagating in air, Opt. Lett. 23 (1998) 382–384. [12] J. Peñano, P. Sprangle, B. Hafizi, A. Ting, D. Gordon, C. Kapetanakos, Propagation of ultra-short, intense laser pulses in air, Phys. Plasmas 11 (2004) 2865–2874. [13] M. Fujii, M. Tahara, I. Sakagami, W. Freude, P. Russer, High-order FDTD and auxiliary differential equation formulation of optical pulse propagation in two-dimensional Kerr and Raman nonlinear dispersive media, IEEE J. Quantum Electron. 40 (2004) 175–182. [14] J. Peñano, P. Sprangle, B. Hafizi, W. Manheimer, A. Zigler, Transmission of intense femtosecond laser pulses into dielectrics, Phys. Rev. E 72 (2005) 36412-1–36412-7. [15] J. Greene, A. Taflove, General vector auxiliary differential equation finite-difference time-domain method for nonlinear optics, Opt. Exp. 14 (2006) 8305–8310. [16] S. Champeaux, L. Berge, D. Gordon, A. Ting, J. Peñano, P. Sprangle, (3 + 1)-Dimensional numerical simulations of femtosecond laser filaments in air: toward a quantitative agreement with experiments, Phys. Rev. E 77 (2008). [17] D. Gordon, M. Helle, D. Kaganovich, A. Ting, Electro-optic and terahertz diagnostics, in: S. Gold, G. Nusinovich (Eds.), AIP Conf. Proc. no. 1299, November 4, 2010, pp. 67–75. [18] M. Helle, D. Gordon, D. Kaganovich, A. Ting, Extending electro-optic detection to ultrashort electron beams, Phys. Rev. ST/AB 15 (2012) 052801-1– 052801-11. [19] J. Jackson, Classical Electrodynamics, John Wiley and Sons, New York, 1975. [20] R. Boyd, Nonlinear Optics, second ed., Academic Press, San Diego, 2003. [21] S. Casalbuoni, H. Schlarb, B. Schmidt, P. Schmuser, B. Steffen, A. Winter, Numerical studies on the electro-optic detection of femtosecond electron bunches, Phys. Rev. ST/AB 11 (2008) 072802-1–072802-18. [22] W. Faust, C. Henry, R. Eick, Dispersion in the nonlinear susceptibility of GaP near the reststrahl band, Phys. Rev. 173 (1968) 781–786. [23] Y. Berozashvili, S. Machavariani, A. Natsvlishvili, A. Chirakadze, Dispersion of the linear electro-optic coefficients and the non-linear susceptibility in GaP, J. Phys. D: Appl. Phys. 22 (1988) 682–686. [24] G. Gallot, D. Grischkowsky, Electro-optic detection of terahertz radiation, J. Opt. Soc. Am. B 16 (1999) 1204–1212. [25] S. Jamison, J. Shen, A. MacLeod, W. Gillespie, D. Jaroszynski, High-temporal-resolution, single-shot characterization of terahertz pulses, Opt. Lett. 28 (2003) 1710–1712. [26] G. Berden, S. Jamison, A. MacLeod, W. Gillespie, B. Redlich, A. van der Meer, Electro-optic technique with improved time resolution for real-time, nondestructive, single-shot measurements of femtosecond electron bunch profiles, Phys. Rev. Lett. 93 (2004) 114802-1–114802-4. [27] G. Berden, W. Gillespie, S. Jamison, E.-A. Knabbe, A. MacLeod, A. van der Meer, P. Phillips, H. Schlarb, B. Schmidt, P. Schmuser, B. Steffen, Benchmarking of electro-optic monitors for femtosecond electron bunches, Phys. Rev. Lett. 99 (2007) 164801-1–164801-4. [28] A. Debus, M. Bussmann, U. Schramm, R. Sauerbrey, C. Murphy, Z. Major, R. Horlein, L. Veisz, K. Schmid, J. Schreiber, K. Witte, S. Jamison, J. Gallacher, D. Jaroszynski, M. Kaluza, B. Hidding, S. Kiselev, R. Heathcote, P. Foster, D. Neely, E. Divall, C. Hooker, J. Smith, K. Ertel, A. Langley, P. Norreys, J. Collier, S. Karsch, Electron bunch length measurements from laser-accelerated electrons using single-shot THz time-domain interferometry, Phys. Rev. Lett. 104 (2010) 084802-1–084802-4. [29] T. Esirkepov, Exact charge conservation scheme for particle-in-cell simulation with an arbitrary form-factor, Comput. Phys. Commun. 135 (2001) 144– 153. [30] C.K. Birdsall, A.B. Langdon, Plasma Physics via Computer Simulation, Institute of Physics Publishing, Bristol and Philadelphia, 1991. [31] C. Decker, W. Mori, Group velocity of large amplitude electromagnetic waves in a plasma, Phys. Rev. Lett. 72 (1994) 490–493. [32] J.-P. Berenger, Three-dimensional perfectly matched layer for the absorption of electromagnetic waves, J. Comput. Phys. 127 (1996) 363–379. [33] E. Lindman, Free space boundary conditions for the time dependent wave equation, J. Comput. Phys. 18 (1975) 66–78. [34] A. Yariv, P. Yeh, Optical Waves in Crystals, John Wiley and Sons, 2003. [35] G. Agrawal, Nonlinear Fiber Optics. Optics and Photonics, second ed., Academic Press, San Diego, 1995. [36] C. Scoby, P. Musumeci, J. Moody, M. Gutierrez, Electro-optic sampling at 90° interaction geometry for time-of-arrival stamping of ultrafast relativistic electron diffraction, Phys. Rev. ST/AB 13 (2010) 022801-1–022801-7. [37] R. Fiorito, A. Shkvarunets, T. Watanabe, V. Yakimenko, Interference of diffraction and transition radiation and its application as a beam divergence diagnostic, Phys. Rev. ST/AB 9 (2006) 052802-1–052802-17. [38] A. Shkvarunets, R. Fiorito, Vector electromagnetic theory of transition and diffraction radiation with application to the measurement of longitudinal bunch size, Phys. Rev. ST/AB 11 (2008) 012801-1–012801-14.
