Summer Research Fellowship Programme of India's Science Academies

Theoretical simulation of High Harmonic Generation by time-dependent Schrödinger equation

Souvick Chakraborty

Department of Physics, Indian Institute of Technology Bhubaneswar, Argul, Khordha 752050

Dr. Atanu Bhattacharya

Department of Inorganic and Physical Chemistry, Indian Institute of Science, Bangalore, CV Raman Rd, Bengaluru, Karnataka 560012


A dielectric medium is an electrical insulator that can be polarised by an applied electric field. If the external field is weak (i.e., weaker than the field binding the electrons to the nuclei), the input radiation acts as a small perturbation to the medium and this produces a polarisation that is directly proportional to the applied electric field. This is called linear optical effect. However, if the external Electric field strength is comparable to that of the field binding the electrons then the relationship between the polarisation and the electric field is no longer linear and is expressed as a Taylor series expansion, involving only the magnitude. Harmonic generation in dielectric solids is well understood and extensively used in modern laser physics where only low order harmonics can be generated. However, in some cases, when a gas or plasma is illuminated by an intense laser pulse, the sample emits high order harmonics of the generation beam. High Harmonic Generation is an extreme nonlinear process in which 10th or even 100th harmonics are generated. The intensity of the laser beam is in the order 1014 W/m2 which can be obtained by focusing a high-power femtosecond laser beam. The generation of the attosecond pulse using HHG can be explained using a three-step semiclassical model.

Keywords: High Harmonic Generation, attosecond spectroscopy, molecular electrostatic potential, time-dependent, Schrodinger equation


The abbreviations that will be used in this project are as follows:

 HHGHigh Harmonic Generation 
 VUVVacuum Ultraviolet 
 IP Ionization Potential 
 ESPElectrostatic Potential 



HHG was discovered by scientists who were interested in the response of atoms to intense laser fields. In nonlinear optics, the general (and quite reasonable) understanding was that high-order processes are less probable than low-order processes. Consequently, on a not so successful route toward vacuum ultraviolet light, researchers used ultraviolet radiation in order to start with VUV photons. At the same time, a lot of progress in the laser technology allowed scientists to focus radiation at very high intensities of the order 1013-1015 W/cm2. Atoms when exposed to such intense radiation, get ionized by absorbing several photons. As a result, above-threshold-ionization and multiple-ionization were observed as unexpected results of such experiments. First results for HHG from VUV radiation were obtained by Rhodes and his co-workers, using a 248-nm excimer laser with high laser intensity [1]. Strong fluorescence was seen from Ar, Kr and Xe ions, together with high harmonics up to the 17th harmonic in Ne. Shortly afterward, high-order harmonics were observed in Xe, Kr and Ar using a 30-ps Nd:YAG laser (1064 nm). The HHG spectra exhibited a characteristic behavior, with a rapid intensity decrease for the first orders, a plateau from the seventh harmonic up to a very high order (e.g., the 29th in Ar) and an abrupt cutoff. A typical spectrum for the High harmonic generation is represented in ​​Fig 1​​.

    Schematic Representation of the HHG Spectrum

    The generation of the attosecond pulse using HHG can be explained using a three-step semiclassical model. The model was proposed by Corkum and Kulander independently in the year 1993 ​[2]​​[3]​. It consists of three steps - ionization, acceleration, and recombination ​[4] ​​[5]​.

    Statement of the Problem

    The purpose of this project is to simulate the HHG spectrum for a number of atoms and molecules and compare it with the experiments. From the HHG spectrum, we can find out the cutoff frequency for each atom and hence the maximum kinetic energy the electron can have. The cutoff frequency determines the maximum harmonic order upto which the HHG can be achieved. The cutoff frequency depends on the laser parameters used in the experiment, and the type of atoms or molecules (i.e. the potential experienced by the electron).

    First, the electrostatic potential of the atom or molecule is to be determined. Then the effects of the high-intensity Laser pulse is studied by simulation using the time-dependent Schrödinger equation and the spectrum is obtained. By matching the theoretically obtained spectrum with the experimental one we can verify if we guessed the correct potential of the atom. Elements or compounds with different ionization energies have different cutoff frequencies. Higher the ionization energy, higher will be the cutoff. The three-step model can be used to determine the electron trajectories for different harmonic orders.

    Objectives of the Research

    The objectives of this project are:

    • To determine the potential of the atoms and molecules.
    • To simulate the HHG for different atoms and molecules.


    High Harmonic Generation is widely used to produce attosecond pulses which is a prerequisite for attosecond physics. It can be used to produce a wide range of harmonics. The harmonic cut-off varies linearly with increasing laser intensity up until the saturation intensity Isat where harmonic generation stops. The cutoff can be increased by changing the atomic species to lighter noble gases but these have a lower conversion efficiency, so there is a balance to be found out. High harmonics are emitted co-linearly with the driving laser and can have a very tight angular confinement, sometimes with less divergence than that of the fundamental field and near Gaussian beam profiles ​[6]​. It can also be used as an alternative x-ray light source.



    The first experiment on HHG was done by Rhodes and his colleagues showed strong fluorescence for noble gas ions like Ar, Kr, and Xe along with high order harmonics for Ne using a 248-nm excimer laser ​[1]​. Shortly afterward, higher order harmonics were observed in Xe, Kr and Ar too using a 30-ps Nd:YAG laser (1064nm) ​[7]​. From the various experiments it was revealed that the HHG spectrum is characterized by a rapid intensity decrease for the first orders, a plateau from the seventh harmonic up to a very high order (e.g., the 29th in Ar) and an abrupt cutoff. Theorists tried to reproduce the experimental results by numerically solving the time-dependent Schrödinger equation for a single active electron, considering only the single atom response and succeeded. However, the early single atom results did not show a locked HHG phase which is a necessary condition for generation of attosecond pulse. In 1993, a simple physical picture of the HHG process was proposed by Corkum and Kulander independently [2][3]. According to this model, the electron tunnels through the Coulomb energy barrier, which is suppressed by the presence of the linearly polarised laser electric field. The electron then undergoes classical oscillations in the laser field. The influence of the Coulomb force from the nucleus is practically negligible at this time. If the electron comes back to the vicinity of the nucleus, it may recombine back to the ground state, thus producing a photon with energy IP (the ionization potential) plus the kinetic energy acquired during the oscillatory motion. The trajectories are obtained by solving the Newton's equation:

    mdvdt=eE\displaystyle m\frac{d\boldsymbol v}{dt}=-e\boldsymbol E

    assuming that the electrons are born with zero velocity at t=0. The energy of emitted HHG photons is then simply

    Eph=12mv2+IP\displaystyle E_{ph}=\frac12mv^2+IP

    where IP is the ionization energy. This simple classical calculation allows us to understand the cutoff behavior at

    Ecutoff=IP+3.17UP\displaystyle E_{cutoff}=IP+3.17U_P

    where UP=e2E24mω2U_P=\frac{e^2E^2}{4m\omega^2}is the ponderomotive energy.

    In 1999, the quantum theory of the three-step model was presented by M. Yu. Kuchiev and V. N. Ostrovsky ​[8]​ and was accepted. Later in 2007, Xue-Shen Liu and Na-na Li theoretically studied the interaction of a He+ ion with fundamental laser pulses and different combined laser pulses [9]. They found out that for low laser intensity the cut-off energy is extended; however, for higher laser intensity the cut-off of the enhanced harmonics is extended under the combination of two fundamental laser pulses. These are some of the works which are relevant to this project.



    To study the interaction of atoms and molecules with high-intensity laser pulses, we have to solve the one-dimensional time-dependent Schrödinger equation, which can be expressed as (atomic units are used otherwise stated):

    itψ(x,t)=[122x2+V(x)ε(t)x]ψ(x,t)\displaystyle i\frac\partial{\partial t}\psi(x,t)=\lbrack-\frac12\frac{\partial^2}{\partial x^2}+V(x)-\varepsilon(t)x\rbrack\psi(x,t)

    where x is the electronic coordinate, ε(t) is the electric field of the laser pulses which is linearly polarised. Here we have used two types of potential functions-one with four parameters and another with seven parameters.

    V1(x)=b(xa)2+cd(xe)2+fy0\displaystyle V_1(x)=-\frac b{\sqrt{(x-a)^2+c}}-\frac d{\sqrt{(x-e)^2+f}}-y_0
    V2(x)=b(xa)2+cy0\displaystyle V_2(x)=-\frac b{\sqrt{(x-a)^2+c}}-y_0

    These parameters can be adjusted in order to get the same ground state binding energy of the respective atoms, ions or molecules. The electric field of the simulated HHG is given by

    ε(t)=ε0f(t)sin(ω0t)\displaystyle \varepsilon(t)=\varepsilon_0f(t)\sin(\omega_0t)

    where f(t) is the pulse envelope function for the fundamental Laser pulse given as

    f(t)=sin2(ω0t2n)\displaystyle f(t)=\sin^2(\frac{\omega_0t}{2n})

    where n is the number of cycles in the envelope. The intensity of the Electric field used is 5.0×1014 W/cm2. The number of pulses in the enevelope is 46. However, for atoms or molecules with low IP, n is lowered to 7. The value of ω0 is taken to be 0.056 (for a wavelength of 750nm).​

    d(t)=<ψ(x,t)dV(x)dx+ε(t)ψ(x,t)>\displaystyle d(t)=<\left.\psi(x,t)\right|-\frac{dV(x)}{dx}+\varepsilon(t)\left|\psi(x,t)>\right.

    The harmonic-generation spectrum is proportional to the modulus squared of the Fourier transform of d (t). The harmonic order of the cutoff is given by

    Ncutoff=(IP+3.17UP)/ω0\displaystyle N_{cutoff}=(I_P+3.17U_P)/\omega_0


    The time-dependent Schrödinger equation is solved and the harmonic spectra is obtained theoretically using a python program for the noble atoms (He, Ne, Ar, Kr), Helium ion, different types of molecules( H2O, and N2). At first, the potential at various points along the required direction is determined using the Gaussian ​[10]​ and plotted as a function of distance from the center of the nuclei. Then the points are fitted using the two potential functions and the different parameters are obtained. Then we put the values of the parameters in the Python program. The python program returns the first four energy eigenvalues and also plots the harmonic spectra as a function of the harmonic order. The calculated value of the cutoff is then compared to that obtained from the simulation.


    • Helium Ion (He+): Xue-shen Liu and Na-na Li in their paper ​[9]​ simulated the HHG of Helium ion using a 'soft Coulomb potential', V(x)=Zx2+aV(x)=-\frac Z{\sqrt{x^2+a}}. But the ESP obtained from Gaussian differs from this 'soft Coulomb potential'. Our aim is to simulate the spectrum using the two potential functions in Equation 5a & 5b​​.
      Plot of the ESP from gaussian and the 'soft Coulombic potential' vs distance from the nucleus for Helium ion.
        The HHG spectrum obtained using the soft Coulombic potential V(x)=Zx2+aV(x)=-\frac Z{\sqrt{x^2+a}}​with Z=2 and a=0.5

        The two potential functions were used to fit the obtained ESP and then the HHG spectrum was obtained.

        ESP fitted using the potential function in ​Equation 5
          ESP fitted using the potential function in ​Equation 5
            ESP obtained from Gaussian fitted using the two potential functions in​ Equation 5a & 5b

            The values of different parameters obtained after curve-fitting were used to simulate the HHG spectrum. Here, we present a comparison of the spectra obtained by using the three types of potential functions.

              Comparison of the HHG spectrum of Helium ion obtained using three different kinds of potential functions. The two fitted and the one fitted functions are respectively the potentials given by​ Equation 5a & 5b​.

              The IPs obtained from the Python program using the three different potentials are given below:

              IPs obtained from different potential functions
              Type of Potential Function Ionization Potential Values(in eV)
               Soft Coulomb Potential 54.4
              One Fitted Potential  54.98
              Two Fitted Potential 54.28 

              We also used the two potential functions to simulate the HHG spectra for He, Ne, Ar and Kr.

              HHG spectrum of He atom using the one fitted function. 
                HHG spectrum of He atom using the two fitted function
                  HHG spectrum of He atom along with the parameters used for fitting. The HHG spectrum for Helium ion using the corresponding potential functions is also plotted to show the difference in cutoff of the two. 

                  The IP of Helium atom obtained from the one-fitted function is 25.88 eV and that from the two-fitted function is 24.97 eV. The experimental value of the IP is 24.6 eV.

                  HHG spectra generated using the one fitted potential
                    HHG spectra generated using the two fitted potential
                      HHG spectra of He, Ne, Ar and Kr atoms using the two types of potential functions.
                      IPs of He, Ne, Ar, Kr using the two types of potential functions.
                      AtomExperimental IP(eV)Calculated IP using One-Fitted Function (eV)Calculated IP using Two-Fitted Function (eV)
                      He 24.6  25.8824.97 
                       Ne 21.6 20.05 21.63
                       Ar 15.7 16.72 15.68
                       Kr14.115.79 14.03 
                      The values of different parameters obtained after curve fitting using the two potential functions.
                      AtomOne-Fitted Function Two-Fitted Function
                       He y=0.46(x0.13)2+0.030.01y=-\frac{0.46}{\sqrt{(x-0.13)^2+0.03}}-0.01y=0.046(x0.41)2+0.010.48(x(0.03))2+0.050.0001y=-\frac{0.046}{\sqrt{(x-0.41)^2+0.01}}-\frac{0.48}{\sqrt{(x-(-0.03))^2+0.05}}-0.0001
                       Ney=0.34(x0.34)2+0.020.03y=-\frac{0.34}{\sqrt{(x-0.34)^2+0.02}}-0.03 y=0.033(x0.53)2+0.0010.41(x0.13)2+0.050.01y=-\frac{0.033}{\sqrt{(x-0.53)^2+0.001}}-\frac{0.41}{\sqrt{(x-0.13)^2+0.05}}-0.01
                       Ar y=0.24(x0.5)2+0.0050.03y=-\frac{0.24}{\sqrt{(x-0.5)^2+0.005}}-0.03 y=0.06(x0.74)2+0.0020.37(x0.11)2+0.170.01y=-\frac{0.06}{\sqrt{(x-0.74)^2+0.002}}-\frac{0.37}{\sqrt{(x-0.11)^2+0.17}}-0.01
                       Kr y=0.22(x0.53)2+0.0020.01y=-\frac{0.22}{\sqrt{(x-0.53)^2+0.002}}-0.01 y=0.02(x2.39)2+0.20.24(x0.67))2+0.010.03y=-\frac{0.02}{\sqrt{(x-2.39)^2+0.2}}-\frac{0.24}{\sqrt{(x-0.67))^2+0.01}}-0.03
                      • Nitrogen Molecule: The two potential functions were used to generate the HHG spectrum of the Nitrogen molecule by using ESP determined along three directions- along the bond, perpendicular to the N-atom and perpendicular to the bond.
                        Geometry of the Nitrogen molecule
                        ESP obtained along three different directions of Nitrogen molecule
                          HHG spectra of Nitrogen generated from the obtained ESP using the one fitted potential function
                            HHG spectra of Nitrogen molecule obtained for ionization of the electron a)along the bond, b)perpendicular to one of the N atoms, and c)perpendicular to the bond

                            The experimental IP of the Nitrogen molecule is 15.58 eV. The calculated IPs and the different parameters obtained from curve-fitting is given below.

                            Calculated IP and parameters of the one-fitted function along different directions of the nitrogen molecule
                            DirectionCalculated IP(eV)Parameters of the One-Fitted Function 
                            Along the N-N bond 15.36 =0.3(x0.43)2+0.0270.02y=-\frac{0.3}{\sqrt{(x-0.43)^2+0.027}}-0.02
                            Perpendicular to the N atom  15.64 y=0.3(x0.44)2+0.0270.03y=-\frac{0.3}{\sqrt{(x-0.44)^2+0.027}}-0.03
                            Perpendicular to the N-N bond  15.71 y=0.34(x0.27)2+0.040.03y=-\frac{0.34}{\sqrt{(x-0.27)^2+0.04}}-0.03
                            • Water (H2O): The ionization potential of water is near to the IP of nitrogen molecule. So, we can expect their cutoffs to be near as well. We generated the HHG spectrum of water. Here we present a comparative study of the HHG spectrum of water and nitrogen molecule.
                              HHG spectra of water and Nitrogen molecule using the one-fitted function.
                              Comparison of IPs of water and nitrogen molecules
                               MoleculeExperimental IP (eV) Calculated IP (eV)
                               H2O 12.6212.66 
                              N2 15.58 15.49 

                              The function used for fitting the ESP in case of water is y=0.22(x0.39)2+0.010.02y=-\frac{0.22}{\sqrt{(x-0.39)^2+0.01}}-0.02​.

                              CONCLUSION AND RECOMMENDATIONS

                              To understand the HHG process from the atomic point of view, it is essential to know how it generates and how the spectra are obtained. Here we studied the theoretical simulation of the HHG process. We have used electrostatic potential (ESP) for the study of atoms and molecules. This helps to get the HHG spectrum for any kind of molecule; and obtaining HHG spectrum for molecule is one of the leading scientific problems to the HHG spectroscopy community. After investigation we found following key features:

                              • The ESP can be used to construct a model potential for atom or for molecule.
                              • The HHG spectrum generated by using ESP can reproduce the structure and the cut off of the HHG spectrum.
                              • Cut off comparison has been made with He, Ne, Ar and Kr atom where IP gradually decreases and as a result the cut off also decreases.
                              • For molecules, N2 and H2O have been compared with respect to their IPs and this reflects in its HHG spectrum as well.

                              Though it has some limitations and we should rectify it by introducing a few terms during the solution of time dependent Schrödinger equation. For example, various people have observed minima on Argon HHG spectrum, which we do not see in our work. To understand more about specific features and introduce the corresponding term in the TDSE solution, detailed quantum mechanical simulation is required.


                              To start with, I would like to express my sincere gratitude to my guide, Dr. Atanu Bhattacharya for his excellent guidance, caring, patience and providing me with an excellent atmosphere for doing research. I also extend my gratitude to all the faculty members of the Inorganic and Physical Chemistry department for providing me a wonderful opportunity to get to know the higher end research here in IISc.

                              Also, I would like to thank my mentor, Mr. Sankhabrata Chandra for guiding me through every step in this project and helping me figure out the various intricacies involved. Without his constant help, this work could not have been completed.

                              I am extremely grateful to the Indian Academy of Sciences for providing me with the fellowship and such a great platform to do a fruitful summer project.

                              I would like to thank my friends and teachers at IIT Bhubaneswar who have continuously helped me and supported me since the last year.

                              Most importantly, I thank my family who have always supported and helped me to pursue my dream. My greatest debt of all to my Mom, the best educators of life lessons, is acknowledged in the dedication.


                              • Marc Vrakking, 2014, Attosecond and XUV Physics: Ultrafast Dynamics and Spectroscopy, Attosecond and XUV Physics, pp. 321-335

                              • Zenghu Chang, 2016, Fundamentals of Attosecond Optics, pp. 165-172

                              • McPherson, A., Gibson, G., Jara, H., Johann, U., Luk, T.S., McIntyre, I., Boyer, K., and Rhodes, C.K.(1987) Studies of multiphoton production of vacuumultraviolet radiation in the rare gases. J. Opt. Soc. Am. B, 4, 595 ​

                              • Corkum, P.B. (1993) Plasma perspective on strong-field multiphoton ionization. Phys. Rev. Lett., 71, 1994

                              • Schafer, K.J., Yang, B., DiMauro, L.F., Kulander, K.C. (1993) Above threshold ionization beyond the high harmonic cutoff, Phys. Rev. Lett. 70, 1599;see also Krause, J.L., Schafer K.J., and Kulander, K.C. (1992), High-order harmonic generation from atoms and ions in the high intensity regime. Phys. Rev. Lett., 68, 3535.

                              • L'Huillier, A.; Schafer, K. J.; Kulander, K. C. (1991). "Theoretical aspects of intense field harmonic generation". Journal of Physics B: Atomic, Molecular and Optical Physics. doi:10.1088/0953-4075/24/15/004

                              • Ferray, M., L’Huillier, A., Li, X.F., Lompré, L.A., Mainfray, G., and Manus, C. (1988): Multiple-harmonic conversion of 1064 nm radiation in rare gases. J. Phys. B, 21, 31.

                              • Kuchiev, M. Y., and Ostrovsky, V. N. (1999) Quantum theory of high harmonic generation as a three-step process. Phys. Rev. A ,60, 3111

                              • Liu, X., and Li, N.(2007) Efficient extension and enhancement of high-order harmonics of He+ by combined laser pulses. J. Phys. B, 41,015602

                              • Gaussian 09, Revision A.02, M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, T. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman, and D. J. Fox, Gaussian, Inc., Wallingford CT, 2016.


                              • Fig 3: Liu, X., and Li, N.(2007) Efficient extension and enhancement of high-order harmonics of He+ by combined laser pulses. J. Phys. B, 41,015602
                              Written, reviewed, revised, proofed and published with