Non Perturbative Time-Dependent Density Functional Theory, TDDFT: Study of Ionization and Harmonic Generation in Linear Di-(N2) and Tri-(CO2, OCS, CS2) Atomic Molecules with Ultrashort Intense Laser Pulses-Orientational Effects

High-order harmonic generation, HOHG, resulting from the interaction of intense ultrashort laser pulses with atoms has been extensively studied in recent years offering potential application as a source of coherent ultrashort radiation in the extreme ultraviolet and soft Xray regions.1-3 For molecules, molecular high order harmonic generation, MHOHG3 offers the possibility of synthesizing attosecond (1 attosecond=10-18 second) pulses for creating and controlling coherent electron wavepackets, CEWP.4 An intuitive and efficient theoretical picture based on a classical three step electron recombination trajectory5-6 model in the laser field has helped to elucidate this process. According to this model, an electron tunnels out from the atom or molecule and may recombine with the parent5,7 or neighbouring8 ion emitting a high energy photon, after undergoing laser–driven motion in the continuum. If the tunnelling electron under the laser field does not return to the parent ion position, it is completely ionized and this is referred to as molecular ionization. The electron can also recombine with a neighbouring ion, an example of laser induced electron transfer, LIET.9 Theoretical and experimental work on HOHG has been mostly devoted to atoms. The study of MHOHG in molecules is at the early stages. In contrast to atoms, for molecules, the returning wave packet in the semiclassical picture of HOHG encounters a core comprising two or more nuclei, which are presumed to behave as pointlike source potentials leading to interference in the MHOHG spectrum. Additionally, it has been suggested that the electron recollision crosssections responsible for MHOHG relates to the projection of the valence orbitals with respect to the direction of the propagation of the recolliding electron, making possible a new molecular orbital tomography10. It is therefore, natural to ask if there is any orientation


Introduction
High-order harmonic generation, HOHG, resulting from the interaction of intense ultrashort laser pulses with atoms has been extensively studied in recent years offering potential application as a source of coherent ultrashort radiation in the extreme ultraviolet and soft Xray regions. [1][2][3] For molecules, molecular high order harmonic generation, MHOHG 3 offers the possibility of synthesizing attosecond (1 attosecond=10 -18 second) pulses for creating and controlling coherent electron wavepackets, CEWP. 4 An intuitive and efficient theoretical picture based on a classical three step electron recombination trajectory 5-6 model in the laser field has helped to elucidate this process. According to this model, an electron tunnels out from the atom or molecule and may recombine with the parent 5,7 or neighbouring 8 ion emitting a high energy photon, after undergoing laser-driven motion in the continuum. If the tunnelling electron under the laser field does not return to the parent ion position, it is completely ionized and this is referred to as molecular ionization. The electron can also recombine with a neighbouring ion, an example of laser induced electron transfer, LIET. 9 Theoretical and experimental work on HOHG has been mostly devoted to atoms. The study of MHOHG in molecules is at the early stages. In contrast to atoms, for molecules, the returning wave packet in the semiclassical picture of HOHG encounters a core comprising two or more nuclei, which are presumed to behave as pointlike source potentials leading to interference in the MHOHG spectrum. Additionally, it has been suggested that the electron recollision crosssections responsible for MHOHG relates to the projection of the valence orbitals with respect to the direction of the propagation of the recolliding electron, making possible a new molecular orbital tomography 10 . It is therefore, natural to ask if there is any orientation 494 dependence of the closely related process of MHOHG spectra. Also, due to the shapes and the symmetries of different molecular orbitals, we are interested to know how molecular orbitals can come into play into the MHOHG and molecular ionization yield processes when different laser intensities and molecular angle versus laser polarisation are taken into account. While theoretical approaches based on the accurate 3D time dependent Schrödinger equation, TDSE, are limited to the one electron sytem 9, 11-12 and 2D models for two electrons system 13 because of limitation of memory and computational time, the major difficulties in the theoretical study of molecular systems in a laser field resides on the multi-center nature of the molecule and also the treatment of the multi electron-electron interactions. Here, we focus on time-dependent density functional theory, TDDFT, methods beyond linear response theory 14 as a tool for studying the nonperturbative response of molecules to intense ultrashort laser pulses. 15 The time independent Density Functional Theory, DFT, has become an ubiquitous method 16 of solving ground state electronic problems in atoms and molecules. Extension to the nonperturbative regime through time-dependent methods has emerged by the existence of a rigorous theorem relating the exact time dependent density to external time-dependent potentials 17 and has been extended to linear response via a linear perturbative version of TDDFT. 18 We focus here on a non-linear non perturbative TDDFT. 15 For the visual understanding of electron motion in our molecules during the time propagation under strong laser pulses, we will use the time dependent approach of the electron localization function, TDELF. [19][20][21][22] Conceptually, the physical idea of the electron localization function, ELF, is based on the fact that a highly localized electron repels other electrons with same spin very strongly due to the Pauli exclusion. The ELF uses the probability to find a second spin-like electron in the vicinity of the first one as a measure for electronic localization. Numerical values of the probability are conveniently mapped on the interval ]0,1] facilitating analysis and interpretation. So, in areas where two electrons of the same spin have a high probability to be found (in reference to the homogeneous gas) the function should tend to 0 by construction. It follows that, in areas where the probability of finding two electrons of opposite spin is high (in reference to the homogeneous gas) the function tends to 1. Consequently, a region of the space with a high value of the ELF corresponds to a region of chemical bond or lone pair of the Lewis theory. Thus, TDELF will allow us to make a clear separation between the HOMO of our molecules and the rest of the valence molecular orbitals 23 , and therefore will enable us to visualise their behaviour in the presence of a laser pulse, such as an example laser induced electron transfer, LIET, which has recently been shown to be controllable through the pulse shape and phase. 9,24 In the present work, non perturbative TDDFT methods 15 as opposed to linear response TDDFT 18 are used to study the orientational dependence of ionization and molecular high order harmonic generation, MHOHG, in tri-atoms CO 2 , OCS, CS 2 and the di-atomic N 2 molecules as a function of laser intensity, I 0 ≥ 3.54 x 10 14 W/cm 2 =10 -2 au (1au; Io=3.54x10 16 W/cm 2 , E o =5x10 9 V/cm) for few cycle 800 nm laser pulses. This work is organized as follow: First we briefly recall the TDDFT formalism, followed by the procedure for calculating ionization rates and the time dependent electron localization function, TDELF. The next sections are devoted to results and discussions. The paper is ended by some concluding remarks for future research.

495
introduced by Hohenberg and Kohn 25 , and extended by Kohn and Sham 26 (KS) is based on the existence of an exact mapping between one-particle density and external potentials. This leads to the density of the interacting system being obtained from the density of an auxiliary system of non-interacting particles moving in an effective local single-particle potential. A time dependent generalisation of DFT was provided by Runge and Gross 17 , showing that there is a one-to-one correspondence between the external (time dependent) potential, (,) ext rt υ , and the time dependent one-electron density,n(r,t), for many-body systems evolving from a fixed initial state. The time dependent electronic density is written as: 27 2 ,, where, N σ is the number of the occupied KS orbital in the spin state  and (,) i rt σ ψ is the orbital obtained through the time dependent KS equations (in a.u.): The first term is the external potential, due to the interaction of the electron with an external laser field and the nuclei, while the second term accounts for the classical Hartree electrostatic interaction between electrons, and the third term is the exchange correlation (xc) potential includes all non-trivial many body effects, and has an extremely complex (and essentially unknown) functional dependence on the density. When an intense electric with a low frequency field is applied to a molecule, electrons tend to be displaced out by ionization; so, the higher ionization potential, IP, or barriers, the more this electronic displacement will be hindered. Therefore, the choice of an approach which accurately reproduces the experimental IP is useful for our analysis. Local density approximation (LDA) has been widely used in strong fields 15 due to its simplicity and applicability to various systems with relatively lower computational cost. However, LDA is constructed from the homogenous gas system and suffers from a wrong asymptotic behaviour originating from the incomplete cancellation of the self-interactions, SI. For these reasons, we have used the Van Leeuwen and Baerends 28 potential, LB94, which introduces a gradient correction to the LDA, exchange correlation so as to reproduce the Coulomb asymptotic behaviour of the potential. (3-c) www.intechopen.com

Coherence and Ultrashort Pulse Laser Emission 496
A convenience of the LB94 potential lies in its explicit dependency on the local density and its gradient. This potential yields excellent estimates of the IP for atoms and small molecules, especially in the framework of the linear response TDDFT excitation energy. 29 This model potential significantly corrects the higher moments of the density, notably for diffuse outer orbitals involved in excitations and in our case in ionization. The LB94 also includes partial and considerable self interaction correction, SIC, resulting from an interaction of an electron with itself, of importance to electronic transport. 30 Such corrections are known to produce spurious oscillations in photoionization calculations, 31 but little is known about the corrections for strong fields in the nonperturbative regime. Our results show that the LB94 potential offers accurate IP's for all orbitals; an essential starting point for accurate ionization. Furthermore, at high intensities, excited states have very short lifetime thus, reducing considerably non adiabatic effects. The neglect of nonadiabatic effects, which are a major difficulty of linear 32 TDDFT are implicit in the strong field approximation, SFA 6 , a current simple nonpertubative model.

a. Ionization -KS orbitals
In the present work, the TDDFT KS equation, Eq.2, was discretized in space using finitedifference (FD) grid techniques. For numerical efficiency we represent the electron-ion interaction by the norm conserving non-local Troullier-Martins pseudopotential 33 , which allows to explicitly treat many fewer electrons, and to avoid the fast oscillations produced by the core electrons. The pseudopotentials are generated such that the pseudowavefunctions imitate the all-electron atomic wavefunctions. The bond length has been set to its experimental 1, 34-36 value, for OCS bonds length are CS=0.156 nm, CO =0.1157 nm, CS= 0.155 nm for CS 2 , CO=0.1162 nm for CO 2 and NN=0.1041nm for N2. The molecule is placed in a large three-dimensional (3D) rectangular grid cell (atomic units au, are used) of dimension |z max | = 90 au (1a.u. = 0.0529nm) and |x max | =|y max | =25 au if the laser is oriented parallel to the molecule. The uniform FD grid spacing, Δz = Δx= Δy= 0.28 au was used and the convergence of the calculations was checked against results which use a smaller grid spacing. The size of the grid is determined by the maximum radius α = E/mω 2 of an electron in a time dependent electric field of maximum amplitude E. Thus at E = 0.1 a.u. corresponding to an intensity I= 3.50×10 14 W/cm 2 (1a.u: Io= 3.50×10 16 W/cm 2 , Eo= 5×10 9 V/cm), α = 30 au at λ = 800 nm, ω = 0.057 au. Thus since α is the size of the recolliding electron trajectory, z max = 90 a.u. >> α ensures that all recolliding electrons responsible for MHOHG are collected in the numerical procedure. The laser intensity is related to the field strength by for pulses satisfying Maxwell's equations. 4,37 The external potential created by an intense laser field is assumed to arise from an oscillating electric field with a cosine envelope. When the laser is linearly polarized parallel to the molecular axis z, (,) ext rt υ is given by T determines the laser pulse duration. The time-dependent equations, Eq.(2) are solved using the Crank-Nicholson scheme 34,36 with Δt=0.02 a.u (1a.u=24x10 -18 s) as the time step. All calculations have been done using eight optical cycles pulse duration on a multiprocessors parallel machine. To absorb the ionized charge density and avoid reflections at the boundary of the box due to ionization during the propagation, we employ a mask function. By using this technique, the KS wave-function is multiplied at each time-step at grid boundaries by a function M(z) which is 1 in the inner simulation region and gradually goes to 0 at the boundaries. The masking function in one dimension is defined as: |z max |=90 au is the extension of the box in the z direction, a=10 au is the width of the mask function. It is assumed that the electrons in the max zz domain are completely ionized and cannot return to the nuclei. A different approach to remove outgoing flux density consists in adding an imaginary potential 38 which causes an exponential decay of the wavefunction or in using the exterior complex scaling coordinate transformation as absorbing boundary. All methods produce identical results for sufficiently large grids. The time-dependent ionization probability, P i,σ (t), of an individual spin-orbital is calculated as: 15 and the normalised remaining electron population as: where IP is the ionization potential, N m is the order of the harmonics and U p = I /4mω 2 is the ponderomotive energy. 5,8 Calculations of MHOHG with ultrashort (few cycles) intense, I ≥ 10 14 W/cm 2 requires special attention to grid sizes and gauges as mentioned in recent work. 39 In the present work, the MHOHG spectrum is calculated during the dynamical propagation as follows. 15 The power spectrum of the dipole acceleration ( ) z dt or equivalently a timedependent Hellman-Feyman theorem in a given direction yields the predicted MHOHG spectra: 7, 39 Recent work has shown that the most accurate method for generating the MHOHG spectra when dealing with molecular systems. 4,12 is to calculate the dipole acceleration, ( ) z dt , from the exact time dependent function according to Ehrenfest's theorem. 12

c. Time-Dependent Electron Localization Function (TDELF)
The ELF is a method for the mapping of electron pair probability in multi-electron systems. It was first introduced by Becke and Edgecombe 22 as a descriptor of electronic localization, using arguments based on the conditional same spin pair probability function. Here, we briefly recall the definitions of the extension of the ELF for time dependent processes. [20][21] The time-dependent ELF or TDELF, η(r,t), is given as: 2 j σ is the squared modulus of the current density and is the time-dependent correction to the original ELF. 22 C(r,t) measures the probability of finding one electron in the near vicinity of a reference like-spin electron at position r and time t. Thus, the Pauli repulsion between two like-spin electrons, described by the smallness of C(r,t), is taken as a measure of the electron localization. C h (r,t) is the corresponding factor found for a uniform electron gas, i.e., a homogeneous electron gas whose density is the same as that of the real system. The TDELF will be used for the time-resolved observation of the molecular ionization, and can thus provide a visual understanding of the electron density transfer dynamics (most notably lone pair and bond electrons) under strong laser pulses. www.intechopen.com

Electronic analysis and ionization
One of the relevant information regarding the ground state geometry of these molecules is that they are linear. The electronic structure from the DFT-KS calculation is depicted in Fig.1, where the three highest occupied KS MO's HOMO, HOMO-1, HOMO-2, with their phases and energies are presented. For the tri-atomic molecules, one finds that the highest occupied molecular orbital, HOMO, is an anti-bonding type, 1π g , formed of π lobes lying perpendicular to the inter-nuclear axis and located on the sulphur or oxygen atoms for respectively CS 2 and CO2, while the anti-bonding orbitals forming the HOMO, 2π , of the OCS is located on the sulphur and the oxygen atom. The HOMO of the diatomic N 2 molecule is a bonding σ g formed mostly by the overlap between the 2p z orbitals of the nitrogen atoms. The computed values of the IP from various ab-initio theory levels are presented in Table 1.  A survey of Table 1 reveals that for the tri-atomic molecules, the IP calculated with the HF approach are higher (around 7%) than the experimental data. This deviation is most related to the absence of the correlation in the HF approach. Using DFT, the serious underestimation of the ionization potential, IP at the LDA level is clearly noticeable. Computed values are approximately 43% lower than the experimental ones. The reason of this large discrepancy lies as we mentioned in the previous theoretical section 2), is the wrong asymptotic behaviour of the LDA exchange-correlation potential which originates from the SI, as opposed to the HF approach (SI free eliminated by the exact exchange). Results using partially corrected SI through the hybrid functional B3LYP moderately improve the calculated value of the IP, the relative error drops from 43% with the LDA to 25%. For the N 2 molecule, the HOMO and HOMO-1 are nearly degenerate in energy, as consequence, the HOMO and HOMO-1 order is inverted in the HF results and we pointed out that CASSCF do not remove this error in N 2 . 55 Again, the LDA and B3LYP functional underestimate by at least 20% the experimental value of the IP. These observations reflect the increasing correlation effect when molecular orbitals are nearly degenerate. Regarding the results with the LB94 potential, one sees that it performs much better, the global agreement with the experiment being satisfied within 3% of error. A convenience of the LB94 potential lies in its explicit dependence on the local density and a gradient correction to the LDA exchange potential so as to reproduce the correct asymptotic behaviour. We note that for the OCS molecule the HOMO-1 and HOMO-2 eigenvalues are nearly degenerate with an inverse molecular orbital order with respect to the experimental data, whereas the symmetric CS 2 (as in N 2 ) respects the experimental order with negative orbital energies nearly equal to the corresponding IP's, thus suggesting that these orbitals should be close to Dyson orbitals. 40 a. Tri-atomic molecules (CO 2 , CS 2 , OCS,) To gain a better understanding of the molecular ionization, we have investigated the time evolution of some relevant KS orbital using Eq.6. Three angular orientations are considered =0 o , 45 o , 90 o . Results using three laser intensities range from lower 3.50×10 14 , to higher 1.40×10 15 and highest 2.99×10 15 W/cm 2 are displayed in Fig.2a, Fig.2b and Fig.2c for respectively the CO 2, CS 2 and the OCS molecule. Evidence of a strong dependence of the KS remaining orbital population, N i, (t), with the laser intensity is observed. It emerges as expected that in general the KS MO ionization is most important for CS 2 due to its lower IP as compared to OCS. The smaller is the ionization potential, IP, of the electronic shell, the easier it can be ionized. When =0 o , for symmetry reason the two components of each 2π or 1π g orbital have the same behaviour as they are degenerate. For the lower laser intensity Io= 3.50×10 14 W/cm 2 , one sees that the HOMO, shows as expected the dominant response to the laser field. It is followed respectively by the inner shell orbitals 3σ/3σ u (HOMO-1 for OCS and CO 2 and HOMO-2 for CS 2 ). As the laser intensity increases, it emerges that the HOMO does not remain the most affected by the external laser field. Instead, our TDDFT/LB94 calculation shows that the ionization of the inner 3σ/3σ u molecular orbital dominates so that it exceeds that of the HOMO and thus presents the dominant response to the field. This behaviour of the HOMO under high laser intensity is related to their symmetry which has a nodal plane containing the molecular axis. This gives rise to a low ionization yield when aligned parallel to the laser electric field, while the 3σ u /3σ has the right symmetry (MO shape) to interact with the laser field as its density is maximum parallel to the field. One notes also that the ionization of the inner bonding 1π u /1π  (HOMO-2 for OCS and CO 2 and HOMO-1 for CS 2 ) is increasing with the laser intensity. At angles =45 o and =90 o , the symmetry behaviour of the different components of the π orbitals is broken (degeneracy is removed) and we have renamed them as follows: components of the 1π g of CS 2 and CO 2 are 1π g x, 1π g y and those of the 1π of the OCS are 1πx, 1πy. When =45 o , independently of the laser intensity applied on these molecules, evidence of the ionization enhancement (compare to =0 o case) of the HOMO is clearly visible. The high ionization of the HOMO is followed by that of the inner 1π u /1π and 3σ u /3σ MO orbital. The enhanced ionization of the HOMO is the result of the combined effects of its lower IP and its alignment along the laser polarization axis at =45 o . The total orbital ionization yield follows the trend CS 2 >OCS > CO 2 inversely proportional to their IP. Finally, when =90 o , once again, one finds that the KS MO's of CS 2 ionize more than the others. One finds also that independently of the laser intensity, the HOMO remains the most affected by the field. For higher laser peak intensities, 1.40×10 15 and 3.0×10 15 W/cm 2 , the ionization of the 3σ u /3σ HOMO-2 of CS 2 or HOMO-1 of OCS and CS 2 is so important that, it tends to exceed that of the HOMO. More remarkably, for CS 2 and CO 2 , our TDDFT calculations predict that KS orbital ionization is lower for =90 o than for =0 o as the laser intensity is increasing mainly due to the contribution of the bonding 1π u/ 1π orbital. Although this orbital has a larger IP than the HOMO 1π g /2π, its density dominates at 90 o whereas the HOMO, 1π g , has a node (zero density) at 90 o in the center of the molecule. We conclude that orbitals which have maximum density in the direction of laser polarization, such orbitals will influence or dominate in ionization and in the MHOHG processes. However, for the OCS molecule, one sees that as the laser is increasing, the two components of the HOMO (2πy, 2πx) present the dominant response to the laser field. This is because in comparison with the CS 2 and CO 2 MOs', the HOMO and the HOMO-2 of the OCS have almost the same shape (density), as consequence, its HOMO ionizes faster as the laser intensity is increasing due to its lower IP.
In conclusion, we find that the KS HOMO ionization increases from =0 o and reaches a maximum when the molecule is aligned at an angle around 45 o with respect to the direction of the laser polarization, and then decreases at larger angles towards 90 o . To see how the molecular orientation as well as each KS orbital contributes to the total molecular ionization, we display in Table II the total molecular ionization yield (see Eq.9) as a function of the laser intensity and the angle molecular orientation, .
I=3.5x10 14  In general the ionization yield follows the trend CS 2 >OCS>CO 2 due to their IP difference. The ionization is significantly enhanced at higher laser intensity in the case =0 o than the case of =90 o . This is attributed to the KS MO geometry of the linear molecule. In fact, at =0 o the molecular ionization yield is high due to the fact that the two components of the HOMO behave in the same manner (same ionization rate) and due to the additive contribution played by the inner 3σ u /3σ MO whose densities are aligned along the laser polarization axis whereas for =90 o , only one component of the HOMO as well as HOMO-2 is very active in the presence of the laser field due to degeneracy removal.

b. N 2 molecule
We add for comparison the two center N 2 molecule for which the LB94 potential MO IP's are in excellent agreement with experiment (Table 1). Fig.2d illustrates ionization calculations done at I=3.50×10 14 , 1.40×10 15 , and 3.0×10 15 W/cm 2 laser peak intensity for =0 o , 45 o and 90 o . At 3.50×10 14 W/cm 2 and =0 o , the bonding HOMO, 3σ g , shows as expected the dominant response to the field followed by the inner, 2σ u , HOMO-2, anti-bonding, while the ionization of the 1π u HOMO-1 is not relevant. At angles =45 o and =90 o , the symmetry behaviour of the different components of the 1π u orbital is broken (degeneracy is removed) and we have renamed them as 1π u x, 1π u y. For =45 o , again the HOMO, 3σ g , is again the most affected by the laser field followed by the inner 2σ u , HOMO-2 which is closely also followed by the 1π u y component of the HOMO-1. At =90 o , while the HOMO, 3σ g , is still the most influenced by the electric field, it is now followed by the 1π u y component of the HOMO-1 and the contribution of the others are irrelevant. At I=1.40×10 15 W/cm 2 , for =0 o , both the HOMO, 3σ g , and 2σ u , HOMO-2 present the dominant response to the field during the first four optical cycle and finally the ionization of the 3σ g exceeds that of the inner 2σ u during the rest of the propagation. At the highest, I=3.0×10 15 W/cm 2 laser peak intensity, the ionization of both the 3σ g , and 2σ u , MOs ionization remain nearly the same despite the 3eV energy gap between them. It turns out that when aligned parallel to the laser field, the 3σ g , and 2σ u , whose electron density is maximal along the laser polarization are the most affected by the field. For the lower I= 3.50×10 14 W/cm 2 , the ionization of the HOMO, 3σ g, is most important due to its lower IP, however, at high laser intensities I= 1.40×10 15 and 3,0×10 15 W/cm 2 , the 2σ u HOMO-2 presents the dominant response to the field due to its antibonding nature (ionized fast) and the symmetry coupling with the HOMO 3σ g . At =90 o and higher laser intensities, I=1.40×10 15 and 2.99×10 15 W/cm 2 , the 1π u y component of the HOMO-1 show the dominant response to the field because its electron density is maximal along the laser polarization axis, similar to the CO 2 results.

TDELF analysis for tri-atomic molecules
The ground state ELF images are shown in Fig.3a for = 0.85 (i) corresponding to highly localized pair electrons such as lone pairs and the two-dimensional representation of the ELF isosurfaces in the plane containing the three atoms of the molecule (ii). At = 0.85 (i), the image exhibits only the lone pair (noted a) of the oxygen atom while the CO bond localization is absent. This is because the lone pair electron on the oxygen is the most localized in the CO 2 molecule. By representing the ELF in the plane containing the three atoms (ii), one can more easily identify the CO chemical bond region (note b) Fig.3b exhibits for comparison in a two-dimensional representation the ELF of CO 2 , OCS and CS 2 molecule. The highly localized bond region (CO, SO, CS) and the lone pair region are clearly identifiable. From the electronic structure analysis (Fig.1), it is evident that the HOMO, 1πg/2π is essentially located in a diffuse lone pair region (a) while the HOMO-1 and the HOMO-2 are essentially located in localized bonding region (b). In Fig.3b, the lobe representing the lone pair on the oxygen (one lobe on the main axis of the three atoms) and the sulphur (two lobes symmetrically distributed along the main axis of the three atoms) and the laser time evolution are shown in detail. Such clear separation between HOMO and relevant inner orbitals enables us to get a deep insight on the main factors (shape and IP) of the KS ionization process as described in the KS ionization section-4. The TDELF images evolve in real time and space for =0 o and =90 o at 3.50×10 14 W/cm 2 laser peak intensity respectively for the molecules to follow. a. CO 2 and CS 2 When =0 o , at the initial time step i.e. 0 optical cycle (o.c.), the ELF picture is that of the ground state. As the time propagation evolves, the first major perturbation occurs on the lone pair region, i.e. the oxygen for CO 2 or the sulphur for CS 2 after 1.75 optical cycles. The electric field being negative LIET is thus expected to occur from left to right. The major perturbation consists of a lengthening of the lone pair region followed by the appearance of a hole, reflecting the depletion/ionization of the electron density. The ionization is found as expected higher in CS 2 (large hole) than in CO 2 . The highest perturbation in the CS 2 is due to its lower IP. At 2.25 optical cycles, (the electric field is positive), a perturbation as exactly described below arises now on the other atom of the oxygen or sulphur. The same asymmetric perturbation behaviour is observed respectively at 3.75 and 4.25 optical cycles with the perturbation getting larger. The asymmetric distribution of the electron density on the two oxygens (for CO 2 ) or sulphur (CS 2 ) atoms shows that the two atoms do not experience simultaneously or instantaneously the same laser field strength, and that, depending on the sign of the electric field, the perturbation affects successively each of the oxygen/sulfur atom (the picture shows that at 1.75 o.c, for negative electric field the atom at previous KS ionization analysis, the picture shows that the major part of the perturbation is located on the region corresponding to the lone pair region (1π g ). As the time propagation evolves the intensity of the perturbation on the lone pair region increases as well, reaching a maximum for the extrema (minimum i.e 3.75 o.c (E<0) or maximum i.e 4.25 o.c (E>0)) of the external electric field and seeming to disappear around zero values of the external field (e.g at 3.0 o.c). The weak deformation of the CO bond domain reflects the weak contribution of the inner orbital to the molecular ionization process in agreement with the previous KS ionization finding. For the perpendicular configuration =90 o , the first noticeable www.intechopen.com Non Perturbative Time-Dependent Density Functional Theory, TDDFT, Study of Ionization and Harmonic Generation in Linear Di-(N 2 ) and Tri-(CO 2 , OCS, CS 2 )… 507 perturbation occurs after 1.75 o.c and corresponds to the lengthening or increase in diffuseness of the lone pair region on the oxygen atoms of CO 2 or the sulphur atoms for CS 2 . The perturbation is symmetrically distributed on the system; this effect reflects the fact that the two oxygen or sulphur atoms feel simultaneously the same effect of the laser field. As the time propagation evolves (3.75 o.c and 4.25 o.c), one also sees that the dominant response comes from the lone pair region. Most notably, one observes that during the propagation, the electronic perturbation is mainly characterized by an interaction between the outer lone pair and the inner bonding CO region. This may be regarded as the manifestation of the enhancement of the ionization of the bonding 1π u molecular orbital localized on the CO bond region. Evidence of the laser field intensity influence is noticeable, the stronger is the intensity peak value, and higher is the perturbation. Once again, one observes that the ionization is found as expected higher on CS 2 than CO 2 , and the perturbation from the bonding region is also higher. This suggests that for the case of CS 2 , both the HOMO and the inner orbitals are significantly affected at I=3.5×10 14 W/cm 2 . We next illustrate results at higher I= 3.0×10 15 W/cm 2 laser peak intensity in Fig.3e for CO 2 . For =0 o , one finds that the first major perturbation starts closer to the beginning of the time propagation and it is mostly localized on one of the lone pair regions of the molecule.
Thereafter, the asymmetric perturbation appears instantaneously on the CO bond and the lone pair domain (1.75o.c), its effect increasing with the duration of the propagation. Therefore, it is difficult to identify which part of the system presents the dominant response to the field; this reflects the contribution of both the HOMO KS orbital (due to its lower IP) and the inner 3σ u MO (due to its shape or its electron density which is maximal along the polarization of the electric field). These results corroborate our finding in the previous section concerning the KS ionization. When =90 o , at the beginning of the time propagation, the electron depletion occurs on the lone pair region and is characterised by diffuseness of this region. The image shows that at 1.75 optical cycles for which E<0, the interaction between the lone pair density and that of the CO bond is visible and an amount of electron depletion is noticeable in the lone pair region. The symmetric distribution of the electron density perturbation on the two oxygen atoms is clearly evident. As the time propagation evolves, the bi-partition of the perturbation domain of the lone pair orbital of each oxygen atom is observable and tends to enhance the overlap between the lone pair and the CO bonding region. This can be seen as the strong perturbation (ionization) of the inner 1π u KS orbital. These observations are consistent with our results in Fig.2a and Fig.2b which show that the IP and the shape of the KS orbital are the main factors which define its ionization. At lower laser peak intensity, the IP is the main factor which governed the ionization process while at higher laser peak intensity; the ionization of the KS orbital is influenced by its shape.

b. OCS
The most interesting feature in the study of the non symmetric OCS molecule with the TDELF function is the fact that the HOMO is located on two different atoms; the oxygen and www.intechopen.com Non Perturbative Time-Dependent Density Functional Theory, TDDFT, Study of Ionization and Harmonic Generation in Linear Di-(N 2 ) and Tri-(CO 2 , OCS, CS 2 )… 509 the sulphur. As found in the previous section, the two atoms ionize differently due to their IP difference. The oxygen is easily recognized with one lobe while the sulphur bears two. For =0 o , at 1.75 (E<0) optical cycles, the perturbation of the electron density appears on the lone pair position of the oxygen atom, while its homolog lone pair on the sulphur atom seems less sensitive to the negative electric field. The ELF after 2.25 (E>0) optical cycle shows the perturbation appears now on the lone pair region of the sulphur atom while the corresponding lone pair region on the oxygen atom remains nearly unchanged. As the propagation is evolving, the intensity of the perturbation on the lone pair region is increasing as well; reaching a maximum at the electric field maxima or minimal. For the orientation =90 o , an overview of the TDELF image at 1.75 and 2.25 optical cycle shows an amount of electron density leaving the system. While the two atoms are feeling exactly the same laser intensity at the same time, the ionized electron density is much higher from the sulphur moiety than the oxygen one, due to the electronegativity difference of the sulphur and the oxygen atom. As expected, the perturbation on the different lone pair region is getting higher as the laser intensity is increasing. Fig.3f shows that the OCS behaves in the presence of the laser as superposition of the moieties, one from CS 2 and the other from CO 2 , giving rise to a non-symmetric perturbation on the system when the laser is oriented parallel or perpendicular to the laser field.
In conclusion, we find that at the lower I=3.50×10 14 W/cm 2 (10 -2 a.u.) laser peak intensity, the lone pair region (HOMO) of linear molecules is the most affected by the field, in good agreement with our finding in the KS-MO ionization analysis in section-4. We show next that harmonics emitted by the molecule at that particular laser peak intensity is also mostly due to the HOMO of each molecule, whereas at higher intensities, inner shell ionization must be considered. 41

Molecular High Order Harmonic Generation-MHOHG
In Fig.4a, Fig.4b, we present the MHOHG power spectra obtained by, Eq.11, emitted respectively by CS 2 and CO 2 at =0˚, 45˚ and 90˚. Calculations are done at I=3.50×10 14 W/cm 2 laser intensity and the ionization yield is given in Table II to see the connection between the molecular ionization and MHOHG spectrum signal. One finds in general that the overall shape of the power spectra resemble each other. For all orientations, features of these spectra strongly resemble those of harmonic spectra from atoms, i.e., a sharp decrease of the first few harmonics from bound orbital electrons, followed by a "plateau" and ending with the 'cut-off' due to recollision of an ionized electron that usually determines the highest harmonic order achievable 5,[8][9]42 and is given by the classical law, N m ω=IP+3.17U p . The computed cut-off number, N m , which is in agreement with the classical law is found approximately at the 55 th harmonic for CS 2 and 57 th for CO 2 ; the difference lies in their IP. The figures show that the harmonic peaks in the plateau region are in general smaller in the case of = 0˚ (parallel polarization) than other orientations. The difference can be as high as one order of magnitude, which is the case for the harmonics between the 35 th and 45 th for CO 2 . The most interesting feature of these spectra appears when one has a close look at the plateau region. There is a broad region with suppressed harmonic emission rate (indicated by the arrow on the figures). This has been interpreted in terms of two center interference between the contributions from two atomic centers. [11][12][43][44] In fact, a CO 2 or CS 2 molecule can be regarded from their HOMO as an elongated diatomic molecule where the point emitters are located on the two sulphur nuclei for CS 2 or the two oxygen nuclei for Fig. 4a. (Color online) TDDFT/LB94 harmonic spectra of CO 2 for some angles, between the molecular axis and the laser polarization direction. Laser pulse duration, eight optical cycles, peak intensity I= 3.50×10 14 W/cm 2 , and wavelength, 800nm. Color online) TDDFT/LB94 harmonic spectra of CS 2 for some angles, between the molecular axis and the laser polarization direction. Laser pulse duration, eight optical cycles, peak intensity I= 3.50×10 14 W/cm 2 , and wavelength, 800nm.
Non Perturbative Time-Dependent Density Functional Theory, TDDFT, Study of Ionization and Harmonic Generation in Linear Di-(N 2 ) and Tri-(CO 2 , OCS, CS 2 )… 511 CO 2 . The suppression comes from destructive interference between the contributions from the two atomic centers discussed below. a. CO 2 The total molecular ionization yield (Table II) is respectively 2.0, 2.3 and 1.2 % when =0ºo, 45º, 90º. From Fig.4a, one sees that when =0°, the spectrum shows the presence of a low signal in agreement with the experimental observation [45][46][47] , but with slight different laser peak intensities centered around the 29th harmonic. When the molecular axis is 45˚ with respect to the laser polarisation, two plots are done, showing the projection of the time dependent dipole acceleration along the z and y axis. The MHOHG spectrum is split into contributions from the z-axis (noted by 45_z on the figure) and the corresponding perpendicular x-axis (noted by 45_x). One sees that the minimum (the destructive interference between the contributions from the two atomic centers) is shifted to the harmonic order N around 35 for both. Furthermore the MHOHG signal increases significantly compared to what is obtained at =0˚ due to the strong HOMO 1π g orbital ionization. Evidence of the elliptical polarization of the MHOHG spectra is visible due to the different signal amplitude from the two component of the harmonic spectrum as shown in the Fig.4a. The interference minimum of the z component of the harmonic spectrum is located near the harmonic order N=31 while it shifts to N= 37 for the x component. It also appears that the harmonic spectrum signal from the z component is higher than that of the corresponding x component, due to the different contribution of the inner 3σ u and 1π u molecular orbitals. When =90˚, no obvious minimum is observable, in comparison to the case =0˚. Instead, one observes an enhancement of the signal of harmonics above N=25. This is mainly the result of the constructive interference from the two oxygen atoms. In fact, they have the same potential and are feeling exactly at the same time, the same electric field. b. CS 2 For CS 2 , the total molecular ionization yield (Table II) is respectively 13.3, 12.9 and 9.1% when =0 o , 45 o , 90 o . The harmonic spectrum at 3.5x10 14 W/cm 2 is given in Fig.4b. A close view reveals that the destructive interference (weak MHOHG signal spectrum amplitude) takes place near N=25 when =0˚, it moves to harmonic order 31 when =45˚. Evidence of the elliptical polarization of the MHOHG spectra is again visible due to the different signal amplitude from the z and x component of the harmonic spectrum near the cut-off. The data shows for harmonics order N<31, the signal from the x component is higher than that of the z component, while the opposite trend is found for N near the cut-off. When =90˚, a shallow minimum is observed near the cut-off, in contrary to the case =0˚, so the evidence of the constructive interference is not clearly observable on the spectrum while the two sulphur atom have the same potential and are experiencing the same laser intensity. This comes from the lower IP of the KS MOs of CS 2 , which gives rise to high ionization involving inner MOs at I=3.50×10 14 W/cm 2 laser peak intensity as shown in table II. To check this, we have further investigated the influence of the laser intensity or the inner orbitals on the interference effect, by repeating the calculation on CS 2 at the low laser intensity Io=1.72 x 10 14 W/cm 2 for =90 o and 0 o in Fig.4c. The computed total ionization yield, Eq.9 is found around 6.5% for the parallel case and less than 4.4% for the perpendicular one. We assume that for this lower laser intensity, the contribution of the inner orbital is unimportant as shown on the orbital population ionization analysis in Fig.4c. The corresponding MHOHG spectrum shows a cut-off around the 35 th harmonic (Fig.4d). As in Fig.4a, the harmonic signals in the plateau region are higher for the perpendicular orientation than the parallel one. For =0 o , a pronounced minimum is visible around the 27 th  harmonic while when =90 o the harmonic spectrum signal is significantly higher as expected resulting from a constructive interference as found in the CO 2 and the H 2 + case. 48 c. OCS Considering the nonsymmetric OCS case, one finds that the amount of ionization, computed using Eq.9 is nearly the same at the end of the eight cycle laser pulse (around 7%) independent of the laser molecule angle . However, the difference between theses orientations on the harmonic spectra (Fig.4e) is remarkable. For the parallel orientation, =0 o , the spectra shows in the plateau region a striking pronounced minimum or destructive interference as for the CS 2 and CO 2 case around N= 39. For =90 o , unlike what we find for the case of CS 2 , a minimum is found around N=25. This interference minimum is rather broad and as deep as is found for =0 o . Although the oxygen and the sulfur atom experience exactly the same external electric field, the minimum observed is due to the electronegativity and the IP difference between the two atoms. To explain these results, we invoke the semiclassical picture of MHOHG. When the driving laser field is polarized perpendicular to the molecular axis, if the two nuclei are identical (CS 2 and CO 2 ), they experience the same tunnelling ionization and the returning electron wavepacket can interact simultaneously with all the ionic cores. Then, both nuclei at the same time emit radiation of the same intensity, frequency and phase 12,49 , thus leading to enhanced harmonic yield. However, if the two nuclei are different (different atomic electronegativity or atomic ionization potential, IP as for the OCS case), although they experience the same laser intensity, both nuclei emit different radiation of different intensity, giving rise to the presence of both minima and maxima in the spectrum. In contrast, the electron wavepacket if they are driven by a laser field parallel to the molecular axis can miss the molecular core. In fact the nuclei screen each other from the returning electron wavepacket. Thus the interference of harmonics emitted from each nucleus occurs over a wide spread of harmonic order thus giving rise to a smaller efficiency for HHG in the case of parallel orientation. Fig. 4e. (Color online) TDDFT/LB94 harmonic spectra of OCS for some angles, between the molecular axis and the laser polarization direction. Laser pulse duration, eight optical cycles, peak intensity I= 3.50×10 14 W/cm 2 , and wavelength, 800nm.

d. N 2
In Fig.4f we present the MHOHG power spectra obtained by, Eq.11, emitted respectively by N 2 at =0˚and 90˚. Calculations are done at I=3.50×10 14 W/cm 2 laser intensity and the ionization yield is respectively 3% at =0˚ and 5% at =90˚. The 'cut-off' due to recollision of an ionized electron that usually determines the highest harmonic order achievable, N m , is found approximately at the 55 th harmonic. When the molecule is aligned parallel to the laser polarization ( =0˚), the spectrum displays two shallow minima, the first at N=25 in agreement with the experimental measurements 45-47 at I=2.3 x 10 14 W/cm 2 and the second one near the cut-off at N=43. This latter is absent in the experimental data most probably due to the different laser intensity used or due to the fact that experimentally the noise level is much higher near the cut-off. The difference with the corresponding harmonic spectrum amplitude for =90˚ is remarkable. An interesting future is the complete absence of a minimum interference and the harmonic signal is strongest near the cut-off region of the spectrum. This reflects the constructive interference in the MHOHG recombination step that results in a maximum in the high harmonic signal originating from the HOMO-1 of symmetry π u (although it is bound by an additional 1.6 eV) due to its shape and also its large recombination dipole as shown in Fig.1 and Fig.3d. Fig. 4f. ((Color online) TDDFT/LB94 harmonic spectra of N 2 for some angles, between the molecular axis and the laser polarization direction. Laser pulse duration, eight optical cycles, peak intensity I= 3.50×10 14 W/cm 2 , and wavelength, 800nm.

e. Comparative MHOHG
In figures 4a-f, we have illustrated MHOHG spectra for CO 2 , CS 2 , OCS and N 2 , i.e. 3-center and 2-center molecules. The TDELF images, Fig 3a-f show that OCS behaves like a superposition of two 2-center moieties, CO and CS. The dominant features of the MHOHG spectra are intensity minima at certain laser-molecule angles, . Thus in Fig.4a, CO 2 exhibits a minimum around N=37 at =45 o . CS 2 z and y harmonic components have minima around N=30, Fig.4c and similarly for OCS Fig.4d. For comparison the 2-center N 2 shows a minimum around N=25 and 41 at =0 o which disappears for perpendicular ionization and recombination at =90 o , Fig.4f. CO 2 (Fig.4a), CS 2 (Fig 4b) and N 2 (Fig.4f) show no significant minima in their MHOHG spectrum at =90 o whereas OCS (Fig.4d) manifests at =90 o a broad minimum centered around N=30. The symmetric molecules N 2 , CO 2 , CS 2 can be treated as 2-center ionization and recombination systems with emitter interference patterns for the multiphoton ionization and photon recombination. 4,12 Thus for molecular bonding orbitals sums of atomic orbital, σ g and π u , both ionization and recombination amplitude will have 2-center interference modulation proportional to cos(p e •R/2) whereas antibonding Non Perturbative Time-Dependent Density Functional Theory, TDDFT, Study of Ionization and Harmonic Generation in Linear Di-(N 2 ) and Tri-(CO 2 , OCS, CS 2 )… 515 antisymmetric combinations σ u and π g produce a sin(p e •R/2) interference 4 for electron momentum p e and internuclear distance R. Thus for perpendicular laser-molecule orientation, =90 o , cos(p e •R/2) =1 and sin(p e •R/2)=0 for electrons ionized or recombining along the laser polarization. CO 2 , OCS and CS 2 ionization probabilities are in fact minimal at =90 o , whereas for N 2 at higher intensities these are equal (Table 2). This result is confirmed by the HOMO antibonding combination of atomic orbitals of the triatomics π g as compared to the bonding 3σ g HOMO of N 2 . Similarly, lower higher IP bonding π u orbitals in the triatomics have dominant ionization at =90 o and higher intensities where cos (p e •R/2) = 1, Fig.2. In the case of harmonics, one observes in general little structure or interference at =90 o for symmetric molecule due to both cos(p e •R/2) and sin(p e •R/2) modulating factors. Clear MHOG intensity minima are found for CO 2 around N=37 at =45 o Fig.4a, for CS 2 and N=30 and =45 o Fig.4b, and N 2 at N=25 for =0 o , Fig.4f. Harmonics at orders N=25 are produced by electrons with momentum p e =1.69, since p e 2 /2=Nω, with ω=0.057 au (800nm). The corresponding electron wavelength is λ e = 2π/p e =0.197 nm, or equivalently 2R for N 2 (R=0.104 nm). Since recombination occurs to the bonding 3σ g , HOMO, cos(p e •R/2) =cos(πR/λ e )=cos(π/2)=0, thus explaining the minimum. The same exercise at =45 o for CO 2 gives λ e =R(CO 2 )=0.2 nm and sin(p e •R/2)= sin(πR/λ e )=0. The same result is found for CS2, since in both cases recombination occurs in an antibonding π g HOMO. In these three symmetric molecular cases, two center interferences, cos(p e •R/2) for bonding and sin(p e •R/2) for antibonding HOMO's regulate MHOHG spectral intensities. The nonsymmetric OCS MHOHG spectra present anomalous intensities as seen in Fig.4d. A clear minimum appears at order N=35 at =0 o and an unexpected minimum around N=30 at =90 o . The charge asymmetry in this molecule is responsible for more complex recombination and MHOHG spectra.

Conclusion
The nonlinear nonperturbative TDDFT calculations presented in this paper for symmetric CO 2 , CS 2 and non-symmetric OCS tri-atomics and the di-atomic N 2 address the TDELF analysis and the MHOHG process occurring from MO ionization at high laser intensities where linear TDDFT is not applicable. As major results, we find that at equilibrium distance R and at intensities I 0 > 3.5x10 14 W/cm 2 , lower inner highest occupied molecular orbitals contribute significantly to ionization and to the MHOHG process. This is due to the symmetry of these orbitals. Even though such lower inner shell orbitals have higher ionization potentials and MHOHG processes occur when orbital densities are maximal with laser polarization direction. Our simulations also reveal that the direction of the laser polarization, , with respect to the molecular axis of the linear molecule can have a significant effect on MHOHG. At some angles, ionization probabilities of different orbitals cross in magnitude. 50 So, while maxima (constructive) and minima (destructive) due to intramolecular interference are found in the dependence of the harmonic intensities on the orientation of the asymmetric OCS molecule, only maxima interference are found for symmetric CS 2 and CO 2 molecules when the laser polarization is perpendicular to the molecular axis. The relative position of the minima interference increases to high harmonic order, N, when the laser-molecule angle, , increases. This is explained in section (6-e) by a simple general formula based on recombination of a continuum electron with a HOMO. 4 These findings are confirmed with the time dependent electron localization function, TDELF, representation through the analysis in term of density perturbations appearing on the TDELF images of each molecule. 23 For < 90 o and at lower laser intensity (Io =10 14 W/cm 2 ), one sees that the HOMO is the most affected by the laser field and a large asymmetry density is found, i.e., we clearly see that during each half cycle, the perturbation occurs alternatively from one nucleus to another (favouring minimal interference) while for =90 o , both nuclei simultaneously feel the same perturbation from the laser field (favouring maximal interference). 51 The present TDDFT simulations show also that the local LB94 potential gives reliable ionization energies 52 , IP, of individual orbitals, thus suggesting that the KS orbitals are close to Dyson orbitals 40 which are the exact many body orbitals for ionization. The relation between these Dyson and the LB94 KS orbitals requires further exploration. Furthermore, the MHOHG spectra show strong dependence on the laser intensity and the laser-molecular angle . This suggests that MHOHG can be a potentially powerful way of studying molecular structure such as orbital tomography. 10 This is made possible by observing and analyzing the shape of the MHOHG spectrum in the plateau region whereas electron recombination is predicted to occur near zero electric field. 49 We also show that elliptical polarization of the MHOHG spectra are influenced by inner shell ionization. In general, for intensities above I=3.5x10 14 W/cm 2 , inner shell orbitals, i.e. lower highest occupied molecular orbitals HOMO-1 and HOMO-2 with larger ionization potentials, IP, than the highest occupied orbital, HOMO, can contribute considerably to total ionization. The main reason is that these inner-valence orbitals have fewer nodes than the HOMO and therefore for certain laser polarization, the density of these lower (but higher IP) orbitals at the laser-molecule angle is much larger than the HOMO. At angle =45 o , the degeneracy of π orbitals is removed by the laser, resulting in different MHOHG polarization components. The phase dependence of these different component harmonics and their relation to electron dynamics has not yet been explored. 49 As a general rule, we note that the ionization is generally preferred as the molecule is aligned along the major axis of the electron distribution in the active molecular orbital. It should be emphasized that while the local LB94 potential yields good molecular orbital ionization potentials IP, further study of MHOHG including the nucleus separation distance effects such as recently reported 53 for H 2 and using more accurate exchange potentials based on long range-short range corrected model 54 offer scope for future accurate characterization of MHOHG processes. The effects of the atomic position (bond length) on the molecular ionization and MHOHG signal are left for future study in order to treat new nonlinear phenomena in molecules such as Charge Resonance Enhanced Ionization, CREI. 51