2.5 Split Operator Method
4.2.1 High-Order Harmonic Generation (HHG)
sin2( πt
2T0) sin(ω0t) t≤ T0
sin(ω0t) T0 ≤ t
(4.2.2)
where T0 is the period from unstable output to stable output.
The intense field often refers to the laser. The laser wavelength λ and intense I were 775nm and 3× 1013W cm−2. We converted the λ and I to ω0 and E0 in atomic unit.
ω0 = 2π× 3 × 108
λ× 4.1341 × 1016 = 5.88× 10−2 (4.2.3) E0 =−
√ I
3.51× 1016 =−2.92 × 10−2 (4.2.4)
4.2.1 High-Order Harmonic Generation (HHG)
To investigate HHG, the most straightforward strategy is the direct solution of TDSE.
The specification for the calculation of TDSE is stated as follow. We used the mapping function (2.1.15) to discretize the finite interval [0, rmax] on the GLL grids. There were 201 GLL grids without the end points, and we did not use the end points, x =−1 and x = 1, in the present calculations. The maximum radius rmax was 200 a.u., and the length scale L was 100. The L in the present calculation is much larger than the L in the time-independent calculation in Section 4.1 because the wave spreads in the intense field. The wave function is expanded in Legendre polynomials, and the maximum l is 19.
A mask function was applied to the present calculations to avoid the reflection of the spreading wave. The mask function fm(r) is the cosine form,
fm(r) =
1 r≤ r0
cos14( π(r− r0)
2(rmax− r0)) r0 ≤ r. (4.2.5) The fm(r) makes the wave function gradually decay to 0 from r0 to rmax. We used r0 = 100 a.u. for the present calculations.
For the time propagation, we turned on the laser field. First 10 optical cycles were for turning on the laser, T0 = 10. The wave function was propagated for 30 optical cycles under the intense field of the constant peak. There were total 40 optical cycles. A optical cycle was divided into 2000 time steps. The procedure of time propagation is identical to the procedure in Subsection 2.5.4.
Figure 4.2.1 shows the remaining probability of a hydrogen atom in the intense field E. The remaining probability Pr(t) is calculated by
Pr(t) =⟨
Ψ(r, t)|Ψ(r, t)⟩
. (4.2.6)
0.9982 0.9984 0.9986 0.9988 0.9990 0.9992 0.9994 0.9996 0.9998 1.0000
0 5 10 15 20 25 30 35 40
Pr(t)=⟨ Ψ(t)|Ψ(t)⟩
time t/Tlaser
Figure 4.2.1: Remaining probability Pr(t) of a hydrogen atom in the laser field of the wavelength 775 nm and intensity 3× 1013W/cm2 duration 40 optical cycles.
The Pr does not change for the first 10 optical cycles. After 10 optical cycles, the Pr decays slowly and linearly. When the time propagation end, t = 40T , Pr still keeps over 99 percent, and few electrons escape to free space. It means that we do not need to consider the electronic emission from the hydrogen atom.
Figure 4.2.2 displays the dipole moment induced in a hydrogen atom by the laser field E. The dipole moments of length form and acceleration from can be calculated from the wave function ψ(r, t) at each time step as
dL(t) =⟨
ψ(r, t) z ψ(r, t)⟩
=
∫
ψ(r, t)r cos θψ∗(r, t)d3r, (4.2.7)
dA(t) =⟨
ψ(r, t) − z
r3 + E0fs(t) ψ(r, t)⟩
=
∫
ψ(r, t)[−z
r3 + E0fs(t)ψ∗(r, t)]d3r. (4.2.8)
−0.15
−0.10
−0.05 0.00 0.05 0.10 0.15
0 5 10 15 20 25 30 35 40
dL(t)
time t/Tlaser
Figure 4.2.2: Length-form dipole dL(t) induced in a hydrogen atom by a laser field with the wavelength 775 nm and intensity 3× 1013W/cm2.
The dipole moment in length form oscillates with the variation of fs(t). The ampli-tude of oscillation is increasing during the first 10 optical cycles and keeps the constant amplitude during the last 30 optical cycles.
Figure 4.2.3 shows the power (harmonic) spectrum of HHG of a hydrogen atom by the laser field. The power spectra, PL(ω) and PA(ω), are square of Fourier transform of
the dipole momentums, dL(t) and dA(t).
Figure 4.2.3: Harmonic spectrum obtained from the Fourier transform that takes 15 cycles from the 20th to the 35th cycle in figure 4.2.2.
We took 15 cycles from the 20th to the 35th cycle to perform the Fourier transform and calculate the power spectra.
First we make sure the physics process of ionization. For the time-dependent system, the ponderomotive energy Up and ionization energy EI are Up ≈ 0.062 and EI = 0.5 for a hydrogen atom. By the Keldysh parameter,
γ =
the ionization is a multiphoton process and may lead to HHG.
In physics, the spectra (see Figure 4.2.3) provide some features about HHG. The difference between peaks is twice as large as the photon ¯hω0. Only odd harmonic peaks occur in the power spectra. Both are the features of HHG and unchanged even if we change the parameters of simulation over a adapted range. The cut-off law can determine the maximum order by Up and EI,
Ecut = 3.17Up+ EI. (4.2.12)
For the spectra in Figure 4.2.3,
nω0 = 3.17× 0.062 + 0.5 = 0.696 (4.2.13) n is about 12. It matches the fact that the values of peaks decay quickly from n = 13.
Tong and Chu (1997b) list the the peak values in their article. It can be compared with us to make sure the accuracy. The peak values in Figure 4.2.3 and the results of Tong and Chu (1997b) are listed in Table 4.2.1. The peak values, PL and PA,converge to the same values. The peak values (n = 3, 5, 9, 11) on the plateau are in good agreement with the work (Tong and Chu, 1997b).
Table 4.2.1: Peak values in the HHG spectrum in Figure 4.2.3
n PL(nω0) PA(nω0) PL(nω0)a PA(nω0)a 3 2.77[−06] 2.78[−06] 2.77[−06] 2.75[−06]
5 1.48[−06] 1.48[−06] 1.04[−06] 1.04[−06]
7 2.71[−07] 2.71[−07] 6.46[−08] 6.54[−08]
9 1.51[−07] 1.51[−07] 1.65[−07] 1.65[−07]
11 1.54[−07] 1.54[−07] 1.51[−07] 1.51[−07]
13 1.25[−08] 1.27[−08] 1.34[−08] 1.34[−08]
15 2.59[−10] 2.51[−10] 2.37[−10] 2.43[−10]
17 2.74[−12] 2.51[−12] 1.98[−12] 1.83[−12]
a Work of Tong and Chu (1997b)
The extra topic is the computational speed. For the present calculation, there were 80000 time steps in a time propagation. The time propagation only took about 100 seconds by Graphic Processing Unit (GPU) accelerator.
Bibliography
Abramowitz, M. and I. A. Stegun (1972). Handbook of Mathematical Functions With Formulas, Graphs, and Mathematical Tables. Applied Mathematics Series. Dover Pub-lications.
Bayliss, A., A. Class, and B. J. Matkowsky (1995). Roundoff error in computing derivatives using the chebyshev differentiation matrix. Journal of Computational Physics 116 (2), 380 – 383.
Canuto, C., M. Y. Hussaini, A. Quarteroni, and T. A. Zang (1988). Spectral Methods in Fluid Dynamics. Springer Series in Computational Physics. Springer Berlin Heidelberg.
Canuto, C., A. Quarteroni, M. Y. Hussaini, and T. A. Z. Jr. (2007). Spectral methods : evolution to complex geometries and applications to fluid dynamics. Scientific Compu-tation. Springer Berlin Heidelberg.
Ceperley, D. (1978, Oct). Ground state of the fermion one-component plasma: A monte carlo study in two and three dimensions. Phys. Rev. B 18, 3126–3138.
Ceperley, D. M. and B. J. Alder (1980, Aug). Ground state of the electron gas by a stochastic method. Phys. Rev. Lett. 45, 566–569.
Chen, J., J. B. Krieger, Y. Li, and G. J. Iafrate (1996, Nov). Kohn-sham calculations with self-interaction-corrected local-spin-density exchange-correlation energy functional for atomic systems. Phys. Rev. A 54, 3939–3947.
Davidson, E. R., S. A. Hagstrom, S. J. Chakravorty, V. M. Umar, and C. F. Fischer (1991, Dec). Ground-state correlation energies for two- to ten-electron atomic ions. Phys. Rev.
A 44, 7071–7083.
Hohenberg, P. and W. Kohn (1964, Nov). Inhomogeneous electron gas. Phys. Rev. 136, B864–B871.
Jiang, T.-F. and S.-I. Chu (1992, Dec). High-order harmonic generation in atomic hy-drogen at 248 nm: Dipole-moment versus acceleration spectrum. Phys. Rev. A 46, 7322–7324.
Jiang, T. F., X.-M. Tong, and S.-I. Chu (2001, Jan). Self-interaction-free density-functional theoretical study of the electronic structure of spherical and vertical quantum dots. Phys. Rev. B 63, 045317.
Kohn, W. and L. J. Sham (1965, Nov). Self-consistent equations including exchange and correlation effects. Phys. Rev. 140, A1133–A1138.
Kramida, A., Yu. Ralchenko, J. Reader, and and NIST ASD Team (2014). NIST Atomic Spectra Database (ver. 5.2), [Online]. Available: http://physics.nist.gov/asd [2014, November 26]. National Institute of Standards and Technology, Gaithersburg, MD.
Krause, J. L., K. J. Schafer, and K. C. Kulander (1992, Jun). High-order harmonic generation from atoms and ions in the high intensity regime. Phys. Rev. Lett. 68, 3535–3538.
Krieger, J., Y. Li, and G. Iafrate (1990). Derivation and application of an accurate kohn-sham potential with integer discontinuity. Physics Letters A 146 (5), 256 – 260.
Krieger, J. B., J. Chen, G. J. Iafrate, and A. Savin (1999). Construction of an accurate self-interaction-corrected correlation energy functional based on an electron gas with a gap. In A. Gonis, N. Kioussis, and M. Ciftan (Eds.), Electron Correlations and Materials Properties. Springer US.
Krieger, J. B., Y. Li, and G. J. Iafrate (1992a, Jan). Construction and application of an accurate local spin-polarized kohn-sham potential with integer discontinuity: Exchange-only theory. Phys. Rev. A 45, 101–126.
Krieger, J. B., Y. Li, and G. J. Iafrate (1992b, Nov). Systematic approximations to the optimized effective potential: Application to orbital-density-functional theory. Phys.
Rev. A 46, 5453–5458.
Kryachko, E. S. and E. V. Ludena (1990). Energy Density Functional Theory of Many-Electron Systems. Understanding Chemical Reactivity. Springer.
Kurth, S., J. P. Perdew, and P. Blaha (1999). Molecular and solid-state tests of den-sity functional approximations: Lsd, ggas, and meta-ggas. International Journal of Quantum Chemistry 75 (4-5), 889–909.
Li, Y., J. B. Krieger, and G. J. Iafrate (1993, Jan). Self-consistent calculations of atomic properties using self-interaction-free exchange-only kohn-sham potentials. Phys. Rev.
A 47, 165–181.
Martin, R. M. (2004). Electronic Structure: Basic Theory and Practical Methods. Cam-bridge University Press.
Midorikawa, K. (2011). High-order harmonic generation and attosecond science. Japanese Journal of Applied Physics 50 (9R), 090001.
Oliver, G. L. and J. P. Perdew (1979, Aug). Spin-density gradient expansion for the kinetic energy. Phys. Rev. A 20, 397–403.
Parr, R. G. and W. Yang (1989). Density-Functional Theory of Atoms and Molecules.
International Series of Monographs on Chemistry. Oxford University Press.
Perdew, J. P., K. Burke, and M. Ernzerhof (1996, Oct). Generalized gradient approxima-tion made simple. Phys. Rev. Lett. 77, 3865–3868.
Perdew, J. P., K. Burke, and M. Ernzerhof (1997, Feb). Generalized gradient approxima-tion made simple [phys. rev. lett. 77, 3865 (1996)]. Phys. Rev. Lett. 78, 1396–1396.
Perdew, J. P. and A. Zunger (1981, May). Self-interaction correction to density-functional approximations for many-electron systems. Phys. Rev. B 23, 5048–5079.
Peyret, R. (2002). Spectral Methods for Incompressible Viscous Flow. Applied Mathemat-ical Sciences. Springer.
Slater, J. C., T. M. Wilson, and J. H. Wood (1969, Mar). Comparison of several exchange potentials for electrons in the cu+ ion. Phys. Rev. 179, 28–38.
Strang, G. (1968). On the construction and comparison of difference schemes. SIAM Journal on Numerical Analysis 5 (3), 506–517.
Telnov, D. A. and S.-I. Chu (1999, Apr). Multiphoton detachment of H− near the one-photon threshold: Exterior complex-scaling—generalized pseudospectral method for complex quasienergy resonances. Phys. Rev. A 59, 2864–2874.
Tong, X.-M. and S.-I. Chu (1997a, May). Density-functional theory with optimized ef-fective potential and self-interaction correction for ground states and autoionizing res-onances. Phys. Rev. A 55, 3406–3416.
Tong, X.-M. and S.-I. Chu (1997b). Theoretical study of multiple high-order harmonic generation by intense ultrashort pulsed laser fields: A new generalized pseudospectral time-dependent method. Chemical Physics 217 (23), 119 – 130. Dynamics of Driven Quantum Systems.
Trotter, H. F. (1959). On the product of semi-groups of operators. Proc. Amer. Math.
Soc. 10, 545–551.
von Barth, U. and L. Hedin (1972). A local exchange-correlation potential for the spin polarized case. i. Journal of Physics C: Solid State Physics 5 (13), 1629.
Wang, J., S.-I. Chu, and C. Laughlin (1994, Oct). Multiphoton detachment of H−. ii.
intensity-dependent photodetachment rates and threshold behavior—complex-scaling generalized pseudospectral method. Phys. Rev. A 50, 3208–3215.
Wang, T.-J. and Z.-W. Yan (2012, Feb). A note to legendre-gauss-lobatto nodes. Journal of Henan University of Science and Technology: Natural Science 33 (1).
Appendix A
Legendre Polynomials
Let
x = cos θ
Properties
PN(±1) = (±1)N (A.0.1)
PN′ (±1) = (±1)N +1N (N + 1)
2 (A.0.2)
Recurrence Relation Following formula in term of Pl−1(x) and Pl−2(x) can be used to find the Legendre polynomial with any l through iteration.
Pl(x) = 1
l[(2l− 1)xPl−1(x)− (l − 1)Pl−2(x)] (A.0.3) Following formula can be used to find the 1st derivative of Legendre polynomial with any l by given Pl(x) and Pl−1(x).
Pl′(x) = 1
1− x2[−lxPl(x) + lPl−1(x)] (A.0.4) Following formula can be used to find the 2nd derivative of Legendre polynomial with any l by given Pl′(x) and Pl(x).
Pl′′(x) = 1
1− x2[2xPl′(x)− N(N + 1)Pl(x)] (A.0.5) Following relation is used in the proof.
d
dx[(1− x2)Pl′(x)] = [(1− x2)Pl′(x)]′ =−N(N + 1)Pl(x) (A.0.6)
Orthogonality ∫ 1
−1
Pl(x)Pl′(x)dx = 2
2l + 1δll′ (A.0.7)
Special Case Expand the differential term on left side of (A.0.6).
(1− 2x)Pl′(x) + (1− x2)Pl′′(x) =−N(N + 1)Pl(x) (A.0.8) On Gauss-Legendre-Lobatto grids, x = xi ̸= x0 ̸= xN, the first term is zero because of Pl′(xi) = 0.
(1− x2i)Pl′′(xi) =−N(N + 1)Pl(xi) (A.0.9) Rearrange the equation.
Pl′′(xi) = N (N + 1)Pl(xi)
(x2− 1) (A.0.10)
Appendix B
Numerical Calculations of Gradient and Laplacian
In Kohn-Sham calculations, we need to numerically calculate the gradient of density ρ(r) or Laplacian of orbital ψ(r) for the parameters of generalized gradient approximation (GGA)
In spherical coordinates, the gradient and Laplacian generally read
∇ = ∂
∂rˆer+1 r
∂
∂θˆeθ+ 1 r sin θ
∂
∂ϕeˆϕ, (B.0.1)
and
∇2 = 1 r2
∂
∂r(r2 ∂
∂r) + 1 r2sin θ
∂
∂θ(sin θ ∂
∂θ) + ∂2
∂ϕ2. (B.0.2)
For the KS calculations in the present research, the orbitals and densities only have two variables, r and θ. The numerical method to find the derivative respect to radius r is provided in Subsection 2.2.4. This appendix provides the numerical differential with respect to polar angle θ.
First, we expand the density ρ(r) in spherical harmonic with m = 0.
ρ(r) = ρ(r, θ) =∑
l
ρl(r)Yl,0(θ, ϕ) (B.0.3)
where the coefficient ρl(r) is
ρl(r) =
∫
Yl,0(θ, ϕ)ρ(r, θ) sin θdΩ. (B.0.4) The spherical harmonics with m = 0 can be in term of the Legendre polynomials Pl(x)
Yl0 =
√2l + 1
4π Pl(x) (B.0.5)
where x = cos θ. The expansion becomes ρ(r, θ) =∑
l
2l + 1
4π ρl(r)Pl(x). (B.0.6)
In the expansion, the radial and angular parts are separated. The differential operators with respect to θ only work on the angular part. Apply the chain rule
dPl(x)
The 2nd-order differential operator works on Pl(x).
d2Pl(x)
In (B.0.2), there is a 2nd-order differential operator with respect to θ 1
Let the 2nd-order differential operator works on Pl(x). With (B.0.7) and (B.0.8), 1 The numerical differential with respect to θ only need to perform once.