Keywords: Recollision, Three-Step Model, Strong-Field Ionization, Tunnel Ionization, Long and Short Trajectories, laser Induced Electron Diffraction, High Harmonic Generation, Double Ionization. Contents 1.1 Soft-Coulomb Potential, Single Active Electron Model, Split Operator Method, Photoelectron Spectra, Volkov States. Contents 2.1 Electron Diffraction, Interference, Recollision, Direct ionization, LIED spectrum. Contents 3.1 Strong Field Approximation, Laser Induced Electron Diffraction, Volkov States, Transition Amplitudes, Direct Ionization Probability, Recollision Amplitude, LCAO, Orbital Imaging, Gaussian-type Orbital, Slater-type Orbital Collective Effects, Lorentz-Lorenz Shift, Cooperative Lamb Shift, Dicke Narrowing, Superradiance, Superradiant Laser Maxwell's Equations, Optical Bloch Equations, Non-Hermitian Wave packets, Two and Multilevel Systems, Nano-Layer, Susceptibility, Resonances, Finite-Difference Time Domain, Perfectly Matched Layers Contents 2.1 Two Level Systems, Resonances, Lorentz-Lorenz Shift, Collective Effects, Reflection, Transmission, Extinction, Susceptibility. Dipole Induced Electromagnetic Transparency Three Level Systems, Overlapping Resonances, Dipole Induced Electromagnetic Transparency, Rubidium Atoms, Strong and Weak Interaction Regime, Lorentz-Lorenz Shift, Collective Effects, Reflection, Transmission, Extinction, Susceptibility Contents 4.1 DIET, Slow Light, Refractive Index, Group Velocity, Group Index, Pulse Shaping, Transmitted Pulse, Reflected Pulse

succinct en Français

Les interférences quantiques apparaissant lors de la superposition cohérente d'états quantiques de la matière sont à l'origine de la compréhension et du contrôle de nombreux processus élémentaires au niveau microscopique.

Dans cette thèse, deux problèmes distincts, qui ont pour origine de tels effets, sont discutés avec leurs applications potentielles :

1. La diffraction électronique induite par Laser (LIED) et l'imagerie des orbitales moléculaires que l'on peut réaliser grâce à ce processus, 2. Les effets collectifs dans des vapeurs atomiques ou moléculaires denses et un effet de transparence électromagnétique induite par interaction dipôle-dipôle (DIET) qui apparaît dans ces systèmes denses.

Le manuscrit est donc constitué de deux parties distinctes qui peuvent lues indépendamment l'une de l'autre. La première partie de cette thèse traite du mécanisme de recollision dans des molécules linéaires simples lorsque le système est exposé à un champ laser infrarouge de forte intensité. Cette interaction provoque une ionisation tunnel du système moléculaire, conduisant à la création d'un paquet d'ondes électronique dans le continuum. Ce paquet d'ondes suit une trajectoire oscillante, dirigée par le champ laser. Cela provoque une collision avec l'ion parent qui lui a donné naissance. Ce processus de diffraction peut être de nature inélastique, engendrant la génération d'harmoniques d'ordre élevé (HHG) ou l'ionisation double non-séquentielle, ou de nature élastique, processus que l'on appelle généralement " diffraction électronique induite par laser". La LIED porte des informations sur la molécule et sur l'état initial à partir duquel les électrons sont arrachés sous forme de motifs de diffraction formés en raison de l'interférence entre différentes voies de diffraction. Dans ce projet, une méthode est développée pour l'imagerie des orbitales moléculaires, reposant sur des spectres de photoélectrons obtenus par LIED. Cette méthode est basée sur le fait que la fonction d'ondes du continuum conserve la mémoire de l'objet à partir duquel elle a été diffractée. Un modèle analytique basé sur l'approximation de champ fort (SFA) est développé pour des molécules simples linéaires et appliqué aux orbitales moléculaires HOMO et HOMO-1 du dioxyde de carbone. L'interprétation et l'extraction des informations orbitalaires imprimées dans les spectres de photoélectrons sont présentées en détail. Par ailleurs, nous estimons que ce type d'approche pourrait être étendu à l'imagerie de la dynamique électronucléaire de tels systèmes.

La deuxième partie de cette thèse traite des effets collectifs dans des vapeurs atomiques ou moléculaires denses. L'action de la lumière sur ces gaz crée des dipôles induits xvi localisés qui oscillent temporellement et produisent des ondes électromagnétiques secondaires. Lorsque les particules constitutives de ce gaz sont assez proches, ces ondes secondaires peuvent coupler les dipôles induits entre-eux, et lorsque cette corrélation devient prépondérante la réponse du gaz devient une réponse collective : les atomes ou molécules sont alors fortement couplées. Ceci conduit à des effets uniques pour de tels systèmes, comme l'effet Dicke, la super-radiance, ou les décalages spectraux de Lorentz-Lorenz ou de Lamb. A cette liste d'effets collectifs, nous avons ajouté un effet de transparence induite dans l'échantillon. Cet effet collectif a été appelé "transparence électromagnétique induite par interaction dipôle-dipôle", ou DIET. La nature collective de l'excitation du gaz dense réduit la vitesse de groupe de la lumière transmise à quelques dizaines de mètre par seconde, créant ainsi une lumière dite "lente". Ces effets, qui ont été prédits théoriquement dans le cas d'un système modèle, sont démontrés ensuite pour les transitions de type D1 du Rubidium-85, et d'autres applications potentielles sont également discutées, en particulier un effet de façonnage spectral d'impulsions laser ainsi que les applications liées au ralentissement très important de la vitesse de propagation de la lumière dans le milieu.

The first part deals with the recollision mechanism in molecules when the system is exposed to high intensity infrared laser fields. The interaction with the intense field will tunnel ionize the system, creating an electron wave packet in the continuum. This wave packet follows an oscillatory trajectory directed by the laser field. This results in a collision with the parent ion from which the wave packet was formed. This scattering process can end up in different channels including either inelastic scattering resulting in high harmonic generation (HHG) and non-sequential double ionization, or elastic scattering often called laser induced electron diffraction. LIED carries information about the molecule and about the initial state from which the electron was born as diffraction patterns formed due to the interference between different diffraction pathways. In this project, a method is developed for imaging molecular orbitals relying on scattered photoelectron spectra obtained via LIED. It is based on the fact that the scattering wave function keeps the memory of the object from which it has been scattered. An analytical model based on the strong field approximation (SFA) was developed for linear molecules and applied to the HOMO and HOMO-1 molecular orbitals of carbon dioxide. Extraction of orbital information imprinted in the photoelectron spectra is presented in detail. It is anticipated that it could be extended to image the electro-nuclear dynamics of such systems.

The second part of the thesis deals with collective effects in dense atomic or molecular vapors. The action of light on the vapor samples creates dipoles which oscillate and produce secondary electro-magnetic waves. When the constituent particles are close enough and exposed to a common exciting field, the induced dipoles can affect one another, setting up a correlation which forbids them from responding independently towards the external field. The result is a cooperative response leading to effects unique to such systems which include Dicke narrowing, superradiance, Lorentz-Lorenz and Lamb shifts. To this list of collective effects, one more candidate has been added, which was revealed during this study: an induced transparency in the sample. This transparency, induced by dipole-dipole interactions, is named "dipole-induced electromagnetic transparency". The collective nature of the dense vapor excitation reduces the group velocity of the transmitted light to a few tens of meter per second resulting in 'slow' light. These effects are demonstrated for the D1 transitions of 85Rb and other potential applications are also discussed.

A Brief History of Light

The interaction of matter with light is one of the fundamental topic of research since the beginning of scientific culture. Understanding the fundamentals of the dynamics of physical systems such as atoms and molecules due to the interaction with light provides a path towards engineering and manipulating the systems. Their electromagnetic response led to many breakthroughs in science and technology and extended the area of our knowledge by providing new methods and tools to explore more in all science fields including chemistry, cosmology, biology etc... Nowadays, light and related technologies became an integral part of our daily life. The UN General Assembly on the 20 th of December 2013 proclaimed 2015 as the "International Year of Light and Light-based Technologies" by "... recognizing the importance of light and light-based technologies in the lives of the citizens of the world and for the future development of global society on many levels .... This year, in fact, ... coincides with the anniversaries of a series of important milestones in the history of the science of light, including the works on optics by Ibn Al-Haytham in 1015, the notion of light as a wave proposed by Fresnel in 1815, the electromagnetic theory of light propagation proposed by Maxwell in 1865, Einstein's theory of the photoelectric effect in 1905 and of the embedding of light in cosmology through general relativity in 1915, the discovery of the cosmic microwave background by Penzias and Wilson and Kao's achievements concerning the transmission of light in fibres for optical communication, both in 1965, ..." (a) The most important conceptual milestones that changed our 'vision about light' are recalled here.

From Corpuscles to Wave

Light propagation has always been an interesting topic of scientific studies. Sir Issac Newton (1642-1726) developed the corpuscular theory of light in his book Opticks (1704) to understand the classical behavior of light, including the straight line propagation of light through a medium by assuming it as a collection of perfectly elastic particles. In (a) Quoted from the Resolution adopted by the UN General Assembly on 20 December 2013:

A/68/440/Add.2; Globalization and interdependence: science and technology for development -Report of the Second Committee 1678 Christian Huygens (1629-1695) communicated to the Académie des sciences his alternative longitudinal wave model of light and published it in 1690 in his Traité de la lumière which did not get a lot of attention at that time. As explained by Huygens, light was then considered as a disturbance in a luminiferous aether. Almost a century after Newton's Opticks, in 1801, the experiments on interference of light by Thomas Young (1773-1829) were explained by the concepts of wave theory. Later in 1818, Augustin-Jean Fresnel (1788-1827) used the idea of a wave character of light to explain the phenomenon of diffraction. He also showed mathematically that the polarization of light can be explained only if light is a transverse wave. Those studies confirmed the wave character of light and the corpuscular theory was rejected completely.

The Electromagnetic Character of Light

In 1845, Michael Faraday (1791-1867) demonstrated in the presence of a dielectric, that linearly polarized light can be rotated by the application of a magnetic field. This was an indication of the relation of light with electricity and magnetism. James Clark Maxwell (1831-1879) unified the electric and magnetic fields via his famous four equations first appeared in 1865 in [1] describing the self-propagating electromagnetic field. The theory suggested that electromagnetic (EM) waves travel with the speed of light though vacuum and Maxwell himself proposed that light is a disturbance of the same hypothetical medium, aether. In 1887 Albert A. Michelson (1852Michelson ( -1931) ) and Edward W. Morley (1838Morley ( -1923) ) disproved in their famous experiment the existence of this unrealistic hypothetical aether. The existence of electromagnetic waves was then confirmed by Heinrich Rudolf Hertz (1857-1894) in 1888 [2].

The dual nature of Light

In 1900 the theoretical model of light was again revised by Max Planck (1858-1947) while explaining the black body spectrum. He proposed that light waves could loose or gain energy only as an integral multiple of a fundamental unit of energy, called quantum, which is proportional to the frequency of the light wave. This quantization of energy exchange explained the black body spectrum. This theoretical proposal of the dual nature of light was used by Albert Einstein (1879-1955) for giving a convincing explanation for the photo-electric effect in 1905, which in turn proved the existence of photons. Later, in 1927, the wave-particle duality of matter, as suggested by Louis de Broglie in his PhD thesis (1924), was confirmed experimentally by Clinton Davisson (1881Davisson ( -1958) ) and Lester Germer and independently by Sir George Paget Thomson (1892Thomson ( -1975)). These discoveries triggered the development of the quantum mechanical description of matter.

The complete and exact description of Light-Matter interactions requires quantized theories of light and matter, an exercise which is often extremely difficult in practice. But in many cases, and for example for the propagation of weak EM fields through a medium, the quantum model of light can be neglected. In those problems light will be described by simple oscillatory functions which obey Maxwell's equations. This semi-classical treatment of the light-matter interaction is used in this thesis.

Maxwell's Equations

Maxwell's Equations are very influential because of their wide applicability in different domains of physics. The most interesting and important remark we obtained from Maxwell's equations is the insight regarding the electromagnetic character of light: "Changing electric fields produce magnetic fields, and changing magnetic fields produce electric fields. Thus the fields can animate one another in turn, giving birth to self-reproducing disturbances that travel at the speed of light. Ever since Maxwell, we understand that these disturbances are what light is" -Frank Wilczek. Those equations are used in this thesis because of the classical behavior of light considered throughout the investigation.

Maxwell's equations in the presence of a dielectric medium can be written as [3] × E = -∂B ∂t (1.1a)

× H = J + ∂D ∂t (1.1b) • D = ρ (1.1c) • B = 0 (1.1d)
where E is the electric field, ρ is the free charge density, J is the free charge current density, H is the magnetic field, D is the displacement current and B is the magnetic flux.

Chapter 1. General Introduction

In a linear dielectric, those electrodynamic quantities are related via

D = εE (1.2a) H = B µ (1.2b)
where ε = ε 0 (1 + χ e ) is the electric permittivity and µ = µ 0 (1 + χ m ) is the magnetic permeability of the medium. ε 0 is the absolute permittivity and µ 0 is the absolute permeability of free space. χ e and χ m are the electric and magnetic susceptibilities that tell us how the system responds to an applied electromagnetic field.

The electric field E and magnetic field B are related to the scalar and vector potentials as

E = -V - ∂A ∂t (1.3a) B = × A (1.3b)
where A is the vector potential and V is the scalar (electrostatic) potential.

The electromagnetic force acting on a charge q moving under the influence of an electromagnetic field with a velocity v is given by the Lorentz force

F = q(E + v × B) .
(1.4)

Maxwell's equations together with the Lorentz force explain classical electrodynamics completely, provided that the boundary conditions related to the electromagnetic fields are given at the interfaces. The interaction of a system with the electromagnetic fields can be studied using these equations (under the semiclassical approximations) by modeling the system under consideration properly so that the quantum behavior of the observables of the system are taken into account.

Structure of the Thesis

In this thesis, a non-relativistic semiclassical approach is used to study the interaction of light with dilute and dense samples of simple quantum emitters, such as atoms or molecules.

My thesis is divided into two major parts: 1) Interaction of strong electromagnetic field with single atoms and molecules, and 2) Interaction of weak electromagnetic field with dense layer of atoms and molecules.

The first part is devoted to the electron dynamics in atoms and molecules induced by the action of intense infrared (IR) laser fields. Due to the high energy applied to the system, the electrons are ionized and are driven back and forth by the oscillating external electric field. This leads to multiple collisions of ionized electrons with the residual ion system, a phenomenon also known as "recollision". The action of strong laser pulses and the induced process of recollision lead to many phenomena, including Above Threshold Ionization (ATI), High Harmonic Generations (HHG), and Laser Induced Electron Diffraction (LIED).

Laser induced electron diffraction can be used to analyse the initial electron density of the system. In the first part of the thesis, the possibility for extracting informations regarding the initial orbital (or electron density) from the LIED signal is studied. This tomography technique is applied to simple linear molecules in two dimensions. Simple approximate expressions are derived to calculate the ionization signal, and are compared with "exact" numerical calculations, so that the essential information for the reconstruction of the initial wave function can be extracted.

In the second part of the thesis, the response of a dilute or dense collection of quantum emitters (atoms or molecules) towards an applied laser field is studied. The influence of the inter-particle interactions on the collective response of the system to the applied field is investigated in detail. A numerical approach based on the Liouville-von Neumann equation is developed for the calculation of the response of the system, and a simple analytical model is derived to analyse the results obtained from the numerical simulations.

Introduction

The dynamics of atomic and molecular systems was studied widely using relatively weak and long laser pulses since the 1960s. Until the end of the 1980s, many interesting phenomena were out of the scope due to the short time scales of the nuclear and electronic motions. The nuclear dynamics is occurring in the tens of femtosecond (fs) time scale while the electronic dynamics is taking place at shorter time scale, of a few hundred attoseconds (as). The introduction of Q-switching in 1961 [1] and mode-locking in 1965 [2] were two great tools that helped for decreasing the pulse duration of such laser systems. To achieve much shorter pulses, say in sub-femtosecond or attosecond regime, one has to go for a combination of high harmonics which can be achieved only via electron dynamics. To turn on nonlinear electron dynamics in atoms and molecules that can eventually produce harmonic generation demand intense laser fields.

The electric field experienced by the ground state electron of the hydrogen atom can be estimated as about 5.1 × 10 11 V/m. Depending on the size of the orbitals in which the electrons are located, the Coulomb attraction experienced by each electron will be different from one another. The field intensity corresponding to the electric field experienced by the ground state electron of the Hydrogen atom is about 3.5 × 10 16 W/cm 2 . Thus having an intense laser field is essential in order to excite an atomic system to such high energies that lead finally to the production of high harmonics. Getting a higher intensity, say 10 13 W/cm 2 at least, became an issue because of the saturation of the system. Amplifying electromagnetic fields to values comparable with the field at atomic sites will ionize and damage laser systems permanently. Chirped Pulse Amplification (CPA) [3] revolutionized the race for high intensity by reducing the cost and the size of laser systems. With the help of the CPA technique, the field can be amplified up to the ultra-relativistic regime.

Basics of Strong Field Physics

Increasing the intensity and reducing the pulse duration triggered the development of Strong field physics and attosecond science [4]. In 1979, Agostini et. al published their results on the multiphoton ionization in xenon atoms at a pressure of 5 × 10 -5 Torr using a linearly polarized 15 ns pulse with an intensity of the order of 10 13 W/cm 2 [5]. For the second harmonic of the laser applied on the system, they observed a signature of above threshold ionization for the first time as an additional peak appeared exactly at the separation of a photon energy from the six-photon ionization peak. Since then, the studies on quantum systems in the presence of strong electromagnetic fields became a hot topic of research. Later, Paul Corkum introduced a semi-classical description for the understanding of the electron dynamics in the presence of strong electromagnetic fields [6].

Electron in intense laser field

The states |Φ v k of electrons in a linearly polarized plane-wave electric field were analyzed within the framework of quantum mechanics by D. M. Volkov as early as 1935 [7]. They are solutions of the time-dependent Schrödinger equation

i ∂ t |Φ v k = Ĥv |Φ v k , ( 1.1)
where |Φ v k are the Volkov states and Ĥv is the Volkov Hamiltonian given in Eq. (4.11a). Volkov states, the quantum states of charged particles in the presence of an oscillating electric field, are widely used in studies including Compton scattering, photoionization or bremsstrahlung effects. Here, the atomic or molecular system is subjected to intense laser fields and ionized. For a general view, it is easier to consider the ionized electron as a classical point charge in the laser field with certain initial velocity and direction, determined by the phase and strength of the field at the time of ionization. Now, let us consider a free electron in the presence of a linearly polarized sinusoidal electric field of amplitude E 0 and frequency ω L . The force acting on the electron within the frame of classical mechanics can be written as

F (t) = -eE 0 sin(ω L t).
(1.

2)

The classical trajectory followed by the electron under the action of the applied field is given by

r(t) = eE 0 mω 2 L sin(ω L t) + V d t + r 0 , (1.3)
where m denotes the mass of the electron, r 0 its initial position and V d the drift velocity. Depending on the ionization time, the number of solutions will vary. They give the times at which the electron in the laser field returns to its initial position. The return of electron to its initial position is referred as a recollision [6]. If the electron is ionized from an atom referred as the core located at r = r 0 , the field will drive back the electron to the core (parent atom or ion) that results in this recollision process. In such a recollision event, the field contributed by the core will be important compared to the case discussed here. In such an oscillating field, the average electron energy is given by the ponderomotive potential

U p = e 2 E 2 0 4mω 2 L .
(1.4)

In the presence of strong fields, an electron wave packet can therefore be formed in the continuum, where the motion of the wave packet is mainly determined by the electric field.

Strong field ionization

In the presence of intense fields, nonlinear effects will show up in the laser-matter interaction. They were initially treated using perturbation theory [8,9]. But the experimental results of [5] where not fitting with the well established lowest order perturbation theory, indicating the breakdown of the approximation and the need for introducing higher order terms which became equally important as the lower order terms. It was because of the dynamic shiftknown as AC-Stark shift -in the atomic energy levels.

The problem became more complex as the laser intensity increased to values at which the first ionization peak starts to diminish in favor of higher order ionization events, disagreeing with the power laws of perturbation theory [8,10,11]. This disappearance or suppression of first order peak is due to the AC-Stark shift in the continuum. If the intensity is large enough, the shift in the ionization potential can almost be equal to the ponderomotive energy U p . That means, in addition to the number of photons required to overcome the ionization energy n I ≈ I p /( ω L ), the system has to absorb additional photons n S ≈ U p /( ω L ) to overcome the AC-Stark shift. As a result, the first appearance of ATI peaks will be at the energy (n S + n P ) ω L -(I p + U p ).

In the infrared regime, strong field ionization can take place due to the formation of a potential barrier in the quasi-static picture of the system potential coupled to the intense laser field. If the field is varying slowly, i.e. if ω L I p / , as the field gets stronger the total potential will be distorted forming a barrier through which the bound electrons can tunnel out. It was first proposed by Keldysh in 1965 [12,13] and later on, it was realized that it was the key process behind strong field ionization in low frequency fields. Fig. 1.2 (a) shows the quasi-static picture of an atomic potential coupled to such an intense field. Since the barrier formation is an intensity dependent phenomenon, above a certain intensity, the barrier formed can be smaller and lower resulting in the exposure of the bound state to the continuum, which eventually will end in the ionization of the system. This type of strong field ionization is known as above the barrier ionization [14]. It is illustrated in Fig. 1.2 (b). A review of the Keldysh theory of strong field ionization can be found in [15].

The three-step model of recollision physics

The three-step model introduced by Paul Corkum [6] is based on three major events in the dynamics of the system while it is exposed to an intense infrared laser field. The action of the electromagnetic field alters the tails of the Coulomb potential and when the oscillating field becomes high enough, a barrier appears, depending on the intensity of the applied field. Fig. 1.3 at t = π/(2ω L ) shows such a situation for the peak value of the oscillating electric field. Depending on the height of the barrier formed, the bound state electron wave function can be ionized either through tunnel ionization or via above the barrier ionization.

The ionized wave packet will be accelerated in the applied field until the field cancels out at t = π/ω L . As the driving field changes sign, the wave packet starts to return to the parent ion from which it has been ionized initially. Since the barrier formation is a slow process compared to the time scales of the electronic dynamics, the bound electrons will be tunneling out when the barrier is thin enough and the ionization rate will increase as the barrier gets thinner over a half cycle and will progressively be stopped when the field changes sign. This behavior of tunnel ionization will give birth to electron wave packets in the continuum at different times t i . During tunnel ionization, the electron wave packets are formed in the continuum, outside of the Coulomb dominant potential well, so that the effect of the Coulomb core can be neglected in the presence of the applied field. This is the main assumption of the Strong Field Approximation (SFA).

Fig. 1.4 shows the classical trajectories followed by electrons born at different moments of an optical cycle. Electrons born early are shown in blue and the ones born later are in red. Curves with colors in between are the trajectories of electrons born in between. Each electron born at time t i will follow the electric field and at certain time t r it will return to its initial position. Early born electrons are traveling very far from the core and hence will return to the core very late. These trajectories, which go very far and return late to the core, are called long trajectories. On the other hand, those born late during the first half cycle can not go far from the core. Relatively, a short interval of time after the birth of these electrons the field will change its sign and will direct the electron back to the core. The kinetic energy distribution of these different electron trajectories upon recollision can be estimated from the classical equations of motion. Being transcendental equations, exact analytical solutions for the kinetic energy of electrons at the moment of recollision are not possible. Fig. 1.5 shows a typical recollision energy distribution obtained from numerical analysis of long and short trajectories. The recolliding electrons are distributed over an energy range between zero and a maximum of 3.17 U p . This peak value of the recolliding energy is possessed by those electrons born at about t i 0.3 (2π/ω L ). This value of energy is the upper limit for the energy accumulated in recolliding process. Thus, as U p increases, electrons will acquire more energy from the field and will posses higher energy at the time of recollision with the core. This rescattering process leads to different outcomes described hereafter: High harmonic generation, non-sequential double ionization and laser induced electron diffraction.

High Harmonic Generation (HHG)

The kinetic energy accumulated during the excursion in the continuum can be emitted in a recombination process. This emitted energy will be of frequencies of an odd multiple of the frequency of the driving field. This process is known as high harmonic generation. In a quantized picture, the process can be seen as the multiple excitation of the electron due to the presence of an intense low-frequency field. The electron absorbs many photons, a process which can lead to above threshold ionization (ATI) [16]. On the course of recollision, the energy absorbed from the laser field can also be emitted as a single photon, constituting the HHG signal.

The generation of harmonics was first observed in an experiment using a ruby laser of intensity 10 3 W/cm 2 [17] and a quartz crystal as nonlinear medium. This system being in the pertubative (weak intensity) regime, the amplitude of the harmonics was decreasing with the energy. In 1980, harmonic orders up to 27 were reported for much higher intensities (about 5×10 14 W/cm 2 ) [18,19]. Later, in 1987, an experiment on rare gases exposed to an intense (3×10 13 W/cm 2 ) picosecond Nd:YAG laser radiation at 1064 nm showed production of high harmonics up to the 21th order in the case of xenon [20] with considerable intensity. Contrary to relatively weak field experiments, this experiment showed the presence of a plateau of constant amplitude up to the 17th order. After the 17th harmonic, the amplitude of the HHG signal drops very quickly. This effect is known as the cut off [21] of the HHG signal. Accurate theoretical studies and understanding of the underlying physics came out later [22][23][24][25], from the three-step mechanism discussed previously. Studies on HHG led to many discoveries and became an important tool for studying the electronic or nuclear dynamics in intense laser fields. HHG is used for the development of coherent sources of light in the XUV regime [26][27][28][29][30] and serves as a key process for developing attosecond pulses and attosecond physics [31][32][33][34][35][36][37]. HHG is also widely used in molecular orbital tomography [38][39][40][41].

Non-Sequential Double Ionization (NSDI)

Double ionization is a process in which two electrons leave the system simultaneously as the liberation of the second electron is assisted by the first one via electron correlations. It was first observed in alkaline earth elements in 1975 [42] and then observed in noble gases in the presence of low frequency fields of high intensity [43,44]. In such a field, the energy of the recolliding electron can be shared with another bound electron of the system, leading to non-sequential double ionization. Multiple ionization was also observed in rare gas atoms in strong low frequency laser fields [45,46]. Studying double ionization yields gives insights to the electron correlations in the system, which can be used for probing the inner atomic electron dynamics that is occurring in the attosecond time scale. Some studies showed that NSDI can occur through two different channels. Either the returning electron has acquired enough energy so that both electrons end up in the continuum after the recollision. In this case NSDI occurs instantaneously and the two electrons leave the system in the same direction. If the energy acquired by the recolliding electron is not sufficient for the ionization of the second electron, the former can still excite the ion so that the second electron can be ionized by the field. This process is delayed and the two electrons leave in opposite directions due to the electron-electron repulsion [47].

Laser Induced Electron Diffraction (LIED)

In conventional electron diffraction, a few kV electron beam from an external source is diffracted by a molecular system via elastic scattering. Likewise, the coherent electron wave packet recolliding with the core can be scattered elastically. This leads to electron diffraction known as laser induced electron diffraction. As in the case of conventional electron diffraction, LIED can be used for the estimation of the structural details of the molecules and imaging the molecular dynamics [48,[START_REF] Safavi-Naeini | Electromagnetically induced transparency and slow light with optomechanics[END_REF]. It can also be used for the reconstruction of molecular orbitals. This is the subject of this part of the thesis.

Orbital Imaging

Understanding the dynamics of molecular systems is an important problem. Properties of atoms and molecules are depending on how the electrons are distributed around nuclei. Hence, understanding changes in the electron densities are of importance especially to study the reaction dynamics of the molecular systems. Due to the ultrashort time scales of the electron dynamics, it is extremely hard to investigate changes in the electron density during a reaction dynamics. A way to image the molecular orbitals is via orbital tomography.

Tomography is a general technique for reconstructing an object from a set of sectional images. The idea of the orbital tomography technique also lies in the same principle of extracting informations about the electron density. It was shown experimentally that such reconstructions of molecular orbitals are possible in the case of aligned molecules in the gas phase using the HHG spectrum obtained from selective ionization [38]. The spectrum of the emitted harmonics is recorded from spatially aligned molecules excited with strong near infrared laser fields. The characteristics of the oriented molecular orbital will be imprinted in this HHG spectrum since the spectrum is a projection of the recolliding electron wave packet on the ground state wave function. Hence, by repeating the procedure at different relative orientations (molecule with respect to the applied field), one can get more details about the orbitals in space. These spectra provide sectional images and then, following the tomographic technique [START_REF] Hau | Light speed reduction to 17 metres per second in an ultracold atomic gas[END_REF], the molecular orbitals can be reconstructed.

The experimental realization of orbital tomography enhanced research in the related fields and there were many modifications done on improving the reconstruction procedure [41,[51][52][53][54][55]. The molecular orbital imaging technique shows a possibility for imaging the chemical reactions and reaction dynamics occurring at the molecular level in femtosecond timescales [48,[START_REF] Safavi-Naeini | Electromagnetically induced transparency and slow light with optomechanics[END_REF]56,57].

In this part of the thesis, molecular orbital imaging using LIED is discussed with the help of a theoretical model based on the strong field approximation (SFA). The model is applied for reconstructing of the highest occupied molecular orbitals HOMO and HOMO-1 of carbon dioxide. The effect of the degree of alignment is also discussed.

The system, numerical method and relevant details for the calculation of photoelectron spectra are discussed in the forthcoming chapter. An analytical model based on SFA is then developed and the model is applied for extracting the information from the calculated photoelectron signals that are relevant for the reconstruction of molecular orbitals.

[49] M. Peters, T. T. Nguyen-Dang, E.

Introduction

During the interaction of strong fields with dilute samples of atoms and/or molecules in which the applied field can no more be treated as a simple perturbation to the system, the ionized electron wave packet on the course of recollision with its parent ion will end up in different output channels as explained in the previous chapter. The process considered here is the elastic scattering of the returning wave packet from the ionic core. On scattering, the wave packet will be diffracted from the core that gives rise to the diffraction pattern. This spectrum contains many information about the system and the ionization processes. The goal of this project is to extract those informations that are relevant for the reconstruction of the initial orbital from which the electron was ionized.

In this chapter, the numerical model is discussed by which the system is treated in order to get the relevant quantities.

System and interaction potential

In general the systems considered are multi-particle systems like atoms or molecules.

Treating exactly such systems in intense laser field is not practical because of the complexity in the potentials and dynamics. But there are some approximations that can be used for transforming the problem tractable both numerically and experimentally. The main ingredients chosen for modeling the problem of laser-matter interaction in the intense nearinfrared are the following: Single active electron approximation, soft-Coulomb potentials, frozen nuclei, dipole approximation.

Single active electron approximation

An atom (molecule) is a multi-particle system composed of a heavy nucleus (nuclei) surrounded by lighter electrons. Depending on the complexity of the system, the number of particles involved in the electronic dynamics can vary from one to a few tens. For such

Chapter 2. System & Numerical Implementation

an atom the potential of the system can be written as

V c = N ı=1 - 1 |r ı | + N =ı+1 1 |r ı | , (2.1)
where r ı is the position of the ı th electron and r ı = r ı -r  are the relative positions of the electrons in the system. An exact quantum treatment of such a system beyond the simplest cases of hydrogen or helium by considering all interaction terms between particles makes it intractable.

In a single-active electron approximation (SAE), we will assume that the laser field will interact only with the most weakly-bound electron in the system. One can then reduce the multi-particle problem to a two body problem, provided the central bound potential is modified such that it accounts for the influence of other interactions empirically. In this approach all other electrons are considered as a frozen [1][2][3].

By choosing the SAE model, the actual system is approximated to a single electron atom similar to the hydrogen atom but with a modified core potential. The potential can thus be reduced to

V c (r ) = - V r , ( 2.2)
where r is a function of the spatial coordinate of the active electron and the coefficient V is chosen empirically such that the potential behaves well (see below).

The model assumes that the response of the system is entirely due to the active electron in this effective potential given by Eq. (2.2). Similar to atoms, SAE models are widely used for understanding the dynamics of simple molecules in intense fields [4][5][6][7].

It is counter-intuitive when one considers the real physical picture in which there are many equivalent electrons that can interact with the electromagnetic field in a similar way. But in practice there are many cases where all interaction terms in the system are not important in determining the results of the laser-matter interaction and in explaining experimental results. The reason for the wide applicability of SAE-models for explaining experimental results lies in the fact that multiply excited states are usually well separated from singly excited states [6,[8][9][10]. In the case of closely spaced levels, the approximation collapses [11]. In the case of processes like non-sequential double ionization, the correlation between the electrons is important and has to be discussed for a real physical picture [12][13][14]. In this specific case the SAE-model is obviously inappropriate. There are other studies showing evidences of the breakdown of the SAE approximation [5,[15][16][17][18] and many attempts for including multi-electron effects were proposed [19][20][21][22].

Soft-Coulomb potential

There is yet another major problem to be fixed in order to be ready with the model system for numerical calculations: The form of the Coulomb potential causes numerical difficulties due to the divergence at the origin. A simple and elegant way to skip this numerical challenge is the choice of a soft-core potential that can mimics the main features of the real potential [23]. In general, soft-core potentials are defined by choosing the following form for r in Eq. (2.2)

r = |r| l + a l 1/l , (2.3)
where r is the coordinate of the active electron and a is a parameter chosen specifically to remove the singularity of the Coulomb potential at the origin. It can also be seen as a parameter that empirically represents the screening effects of the other electrons in the system. l is a power that can be chosen depending on the problem. Here we have chosen the usual value of l = 2 [24, 25] (a) .

The exact value of a is chosen such that the energy of the lowest eigenvalue of the potential matches the ground state energy of the system. In addition, the non-zero value chosen for a converts the infinite Coulomb well to a finite well which falls like -1/|r| asymptotically in order to represent the long range part of the potential accurately. Thus the effective soft-Coulomb potential for an atom under the SAE approximation can be written as

V c (r) = -V |r| 2 + a 2 . (2.4)
In general, for molecules composed of N atoms, the soft-Coulomb potential is given by

V c (r) = N =1 -V  |r -R  | 2 + a 2  , (2.5)
where R  is the position vector of  th nucleus. The functions V  and the parameters a  are chosen depending on the system.

In this part of the thesis, the strong field ionization and the associated electron scattering from a symmetric CO 2 molecule are discussed in detail. The system is prepared by aligning it along the so-called y-axis. The C-O bond length is denoted by R and a linearly polarized laser pulse is applied perpendicular to the alignment direction.

The choice of molecule fixes the potential and is chosen from [26], in which the same system was studied. The potential is given by

V c (r) = 3 =1 -V  (r) |r -R  | 2 + a 2  , ( 2.6)
where

V  (r) = V ∞  + V 0  -V ∞  exp - |r -R  | 2 σ 2  , (2.7)
where V ∞  is the he effective nuclear charge of  th nucleus as seen by the electron at infinite distance, V 0  is the bare charge of the  th nucleus and the decrease of the effective charge of the nucleus with distance due to the distance-dependent electron-electron screening effects is taken into account via the parameter σ  . The pulse is chosen such that its duration is smaller than the time scale of the nuclear dynamics of the molecule. We therefore further simplify the problem by freezing the nuclear coordinates. If the polarization of the applied field is taken as the x-axis, as shown in Fig. 2.1, the three dimensional problem is reduced to two dimensions with all relevant dynamics of the system limited in the xy-plane.

Interaction potential

Being in the long wavelength regime (IR domain), the interaction of the system with the applied field is taken into account under the dipole approximation. Thus the total potential energy of the system can be written as

V (r, t) = 3 =1 -V  (r) |r -R  | 2 + a 2  -µ.E(t) , (2.8)
where µ denotes the electric dipole. The Hamiltonian of the system can then be written a

Ĥ(t) = -∇ 2 /2 + V (r, t). (2.9)
It is important to note that the interaction potential and the Hamiltonian are written here in the length gauge [27]. Being gauge invariant, the results of exact quantum mechanical calculations obtained in different gauges should be the same. But the calculations based on approximate theoretical models may be gauge-dependent [28][29][30]. It is chosen that the length gauge will be employed for all calculations in the present thesis.

Numerical Method

As introduced previously, the system considered here is the symmetric CO 2 molecule. The HOMO is of 1π g symmetry and the HOMO-1 is 1π u . Fig. 2.2 shows these two orbitals for the equilibrium geometry. For the HOMO, which has an ionization energy of 13.8 eV at equilibrium separation, the electron density is concentrated around the two oxygen atoms and is characterized by two nodal planes. For the HOMO-1 of ionization energy 17.6 eV, the electron density is distributed in all atoms [31].

Being a polyatomic molecule, CO 2 possesses additional degrees of freedoms including vibrational degrees. Since the molecule in the following studies is assumed to have frozen nuclei, one has to be careful with the choice of pulse duration. The CO 2 molecule has three vibrational modes. Out of them, asymmetric bending has a time scale of about 15 fs. Thus the constraint on the nuclear degrees of freedom restricts us from choosing a pulse that lasts for more than this time. It is also important to keep in mind that choosing a pulse longer than the time scale of asymmetric vibration will essentially invalidate the assumption of a symmetric molecule. The following numerical model is prepared for treating this particular problem which can be adapted to other systems very easily.

Solving the time-dependent Schrödinger equation (TDSE) using numerical techniques is the key for studying the dynamics of such a quantum system in the presence of the time dependent potential V (r, t). In general, solving the TDSE numerically is based on three steps:

1. The representation of the initial wave function and operators related to the problem on a finite spatial grid.

2. The propagation of the wave function in time.

3. The analysis of the wave function at a specific time as required by the problem for calculating the desired quantities.

These steps are discussed in the following subsections.

Spatial grid

The implementation of the numerical simulation relies on the definition of the grid which will be used for describing the system and the operators. The grid, in general, needs not be an equidistant grid. There are many methods using non-uniform spatial grids for solving the TDSE like [32] for instance. The advantage of having non-uniform grids is that it gives a freedom to choose a very fine mesh in regions where the wave function is of more interest or where it varies quickly, and to choose a relatively wide mesh in regions of least interest or where the wave function varies smoothly. But the efficiency of the FFT algorithm associated with uniform grids is a strong argument, in terms of computation times, in favor of uniform grids without compromising on the accuracy of the computation.

The spatial grid is defined as follows:

x ı = x min + (ı -1) ∆x, ı = 1, 2, • • • , N x , (2.10)
where ı is the index associated with a grid point and N x is the total number of grid points in the x direction. The grid step is given by

∆x = (x max -x min ) /(N x -1), (2.11)
where x max and x min are the upper and lower limits of the grid in x. A similar grid is used in the y direction, giving finally a two dimensional spatial grid. In principle, the space in x that one has to take into account for the exact solution of the TDSE is continuous and infinite.

Discretization in space converts all continuous integrals into finite sums. Indeed, any x integral of a function g(x) can be expressed approximately as

∞ -∞ g(x)dx ı g(x ı ) ∆x . (2.
12)

The methods also requires to discretize all operators and wave functions on the defined grid for the calculations. Out of many ways to solve the TDSE, the system under the SAE is simulated here using the split operator method.

The code is prepared as described in [33,34], and detailed below.

Split operator method

The split operator method was used first in the context of the paraxial wave equation in optics which is quite similar to the TDSE in the mathematical point of view [35]. The method was then adapted for solving the TDSE to make use of its easiness in applying fast Fourier transforms (FFT). As the name indicates this method relies on the effective splitting of the operator used for propagating the wave function.

A formal solution for the TDSE at any time t can be written as

Ψ(t) = Û (t, t 0 )Ψ(t 0 ) (2.13)
where t 0 is the initial time and t is the arbitrary time. Û (t, t 0 ) denotes the propagation operator. Solving Eq. (2.13) is as difficult as solving the TDSE. But for time dependent systems, if the Hamiltonian varies only a little over a chosen time step ∆t, Eq. (2.13) can be approximated to

Û (t, t 0 ) ≈ Nt =1 exp -i Ĥ ∆t (2.14)
where Ĥ denotes the Hamiltonian at time t  = t 0 + ( -1/2) ∆t and ∆t = (t -t 0 )/N t is the time step used while propagating the wave function in time. Thus it is clear that the value of ∆t has to be very small for an accurate description of the dynamics. In principle the Hamiltonian of the system can be partitioned in different ways depending on the circumstances, most probably to find approximate solutions for the problem. Here in this case of solving the TDSE numerically, the Hamiltonian is split in two parts: The first part contains all explicit functions of the momentum p and the other part contains all explicit functions of the coordinate r, i.e

Ĥ(t) = T + V (r, t) (2.15)
Using the partition used above, the components of Eq. (2.14) can be written as

e -i Ĥ∆t = e -i 2 V ∆t e -i T ∆t e -i i 2 V ∆t + O(∆t 3 ) (2.

16)

The advantage of separating the spatial from the momentum dependence of the Hamiltonian is the fact that non-local quantum mechanical operators in the spatial representation are local in the momentum representation and vice versa. Indeed, V (r, t) is diagonal in the coordinate representation and the operator T is diagonal in the momentum representation. The effect of a localized operator on the wave function represented in the appropriate space is then a simple multiplication. Hence calculating the action of the first exponent e -i 2 V ∆t on the initial wave function Ψ(r, t = 0) is easy when the wave function is properly described on the spatial grid. The resultant wave function is then transformed to the momentum space via Fourier transform (numerically using a fast Fourier transform or FFT) for the application of e -i T ∆t . The wave function obtained after this is Fourier transformed back to the spatial representation and the final operator e -i 2 V ∆t is applied. This procedure is repeated for all successive time steps (N t ) to get the final wave function.

Note that the ordering V /2, T , V /2 can also be changed for T /2, V , T /2 for the application of above procedure. Both methods are accurate to the second order in ∆t and the error is proportional to [[ T , V ], Ĥ]∆t 3 . The choice between these two approaches is merely a question of the implementation of the problem, the two approaches being as efficient and accurate.

Absorbing boundaries

A common difficulty with propagating wave functions in a finite grid is the unwanted reflection of outgoing wave functions at the boundary of the simulation region (end of the grid). Here, the FFT algorithm used for Fourier transforms will give rise to unphysical reflections which have to be avoided for the proper description of the dynamics. This can be done by patching the boundaries of the grid with an absorbing term C 2D (r), or mask function, which will absorb the outgoing wave function.

The procedure for including such boundaries is rather complicated because of the fact that any kind of action on the wave function will add additional momentum components to the wave packet that ultimately affect the accuracy of the calculation. This error on the calculation cannot be completely avoided but can be reduced by optimizing various parameters of the mask function. C 2D (r) can be any normalized piece-wise function which has value 1 in the regions where the dynamics is taking place (simulation region) and goes to zero smoothly in the regions close to the boundaries. Multiplying such a function with the wave function on the course of the evolution at each time step will solve the problem of unwanted reflections.

The mask function C 2D (r) chosen in the following calculations is

C 2D (r) = C(x) × C(y) (2.17)
where r = √ x 2 + y 2 . The one dimensional components of the mask function are given by

C(x) =            0 |x| 2 x s 1 + cos π x-xs xs /2 x s < |x| < 2 x s 1 |x| x s (2.18)
where x s is a pre-fixed position in space defining the range of the simulation region. It should be adjusted with the specific problem under consideration. This value depends on the laser parameters used and on the type of dynamics considered. x s also defines the edges of the grid since in our case we can obviously use x max = -x min 2 x s because the wave function is completely absorbed beyond this limit. The function defined above is one out of many possible forms that are used widely [1]. For example in [1], the mask function is a similar cosine function but with a power of 1/8. For the problems which will be considering here we have verified numerically that using different forms of C(x) does not have a dramatic effect on the results obtained. Figure 2.3 shows the spatial dependence of the mask function in one dimension. The piece-wise parts are clearly marked using shaded regions. The smooth variation of the mask function over two wide ranges of the grid that are marked as regions 2 will make sure that the fast momenta added to the wave packet due to the action of the mask function are minimized so that the dynamics of the system is not affected with artifacts. While trying to smooth absorbing boundaries one should keep in mind that, because of the broad range of energy components in the wave packets, it is not possible to design a perfect absorbing boundary mask function. The simulation region which is marked as 3 is the part of the grid where the dynamics of the system is taking place. The dynamics can lead to the elongation of wave packets due to the action of strong laser fields and one has to make sure that the simulation region covers this dynamics effectively. Such optimization of the grid will help in fixing x s . The grid should be designed wide enough so that the boundaries are beyond the slowly varying part of the mask function. An alternative way to implement an absorbing boundary condition is by adding an artificial purely imaginary component to the potential part of the Hamiltonian, similar to an imaginary optical potential in optics. In split operator methods this technique helps to dampen the wave function where the imaginary terms are defined. It has been shown that both methods, adding an imaginary part to the Hamiltonian or using a well defined mask function, give similar results.

Asymptotic analysis of the wave function

The wave packet is propagated as described above and the procedure for extracting the photo-electron spectrum is based on analyzing the wave function near the grid boundaries, where the Coulomb potential is negligible.At each time step, the total wave function is split in two parts: The inner part Ψ I (r, t) and the asymptotic part Ψ A (r, t). It can be done by multiplying the wave function Ψ(r, t) with C 2D (r). The inner part of the wave function, which is the part of wave packet that interacts with the attractive Coulomb potential, can be separated out as follows

Ψ I (r, t) = C 2D (r)Ψ(r, t).
(

2.19)

The asymptotic part is then

Ψ A (r, t) = [1 -C 2D (r)] Ψ(r, t). (2.20)
This asymptotic part of the wave function can be considered as evolving under a simplified Hamiltonian where the Coulomb potential is neglected. Hence it can be seen as a free electron wave packet propagating in an electric field.

The asymptotic part of the wave function can thus be expanded in Volkov states in order to propagate it analytically. For this purpose, the numerical procedure we are following is described in details in Refs. [33,34]. It allows to calculate accurately and efficiently the accumulated ionized wave packet at any time t after the end of the pulse in the momentum representation. The photo-electron spectrum is directly obtained from the square modulus of this propagated momentum wave function.

In this approach, the calculated ionization amplitude is approximate. It is in part because of the fact that there are very slow electrons which will take infinite time to reach the asymptotic boundary. The wave function is therefore propagated for a long period of time after the end of the pulse to collect all electrons, whatever their energy. This is tested numerically until convergence is achieved.

Preliminary Results

The CO 2 molecule is assumed to be aligned along the y axis. The applied field is acting perpendicular to the molecular axis. This scheme sets a natural orthonormal coordinate system for the problem. The linearly polarized applied field ionizes the system both in the positive and negative directions of this axis and drives back and forth the ionized electron via the interaction potential. Hence the electron motion is taking place in the xy-plane.

The numerical code was prepared with a spatial grid in x and y such that -250 a.u x 250 a.u and -250 a.u y 250 a.u with a time step ∆t 0.45 a.u. The potential parameters for CO 2 are taken from [26].

Imaginary time propagation: Calculation of the initial state

The initial state Ψ(r, 0) of the system can be any of the eigenstates Ψ n (r) of the system Hamiltonian. These states are calculated using the split operator method in imaginary time t = -iτ . The time evolution of an arbitrary wave packet in imaginary time can be written as

Ψ(r, τ ) = n a n e -Enτ Ψ n (r) . (2.21)
It is clear from the above equation that as the propagation goes on, the wave packet will concentrate to the lowest energy eigenstate. This method is very efficient in calculating the ground state of the system. But here it is also aimed at getting orbitals other than the ground state. It can be obtained by recognizing the possible solutions with their symmetry. The solutions should be correctly symmetrized and orthonormalized at each time step so that the right solutions are picked from different numerical solutions. This can be done via the Gram-Schmidt orthogonalization procedure. More details of the imaginary time propagation (ITP) technique are given in [36]. While calculating the initial state, it is important to look at the symmetry of the solution obtained via ITP. It is equally important to verify the convergence of the calculation to make sure that the obtained solution is accurate within the model used. Here, the convergence is checked by looking at the energy of the calculated state.

Let Ψ i (r) be the initial guess function which is propagated using ITP. The energy associated with this state can be calculated via

E 0 = Ψ i | Ĥ|Ψ i (2.22)
The function is then propagated in imaginary time to find the initial state. Let us thus denote Ψ(r, t n ) as the state obtained at time t n from this initial guess function. Its associated energy is given by

E n = Ψ(t n )| Ĥ|Ψ(t n ) (2.

23)

The calculation of the initial state is said to be converged if the following condition is satisfied

E n+1 -E n E n = ∆ ε , (2.24)
where ε 1 is a convergence parameter. If the time step used for the ITP is too large, the propagation can ultimately lead to an inaccuracy of the calculated state or even to an increase in the energy of the propagated wave function. For better convergence, the time step should thus be regularly decreased during the computation. This change of time step will be repeated until the condition (2.24) is fully achieved.

A graph of this optimization problem is given in Fig. 2.4. Panel (a) shows the relative error (2.24) as a function of the iteration step in log scale and panel (b) is the calculated energy (2.23) as a function of the iteration step. The relative error is calculated at every step and it is decreased as the wave function is propagated longer in imaginary time: When convergence is first achieved the time step is decreased and another round of optimization starts. This procedure will create a spike-like structure in the graph shown Panel (a) and a step-like structure in the graph of the energy shown Panel (b). These regular decreases of the time step will be repeated till the calculation is fully converged.

Ionization of the HOMO and HOMO-1 of CO 2

The HOMO and HOMO-1 of CO 2 are calculated using ITP. Electron density plots corresponding to these states are given in Fig. 2.5. The equilibrium bond length for the molecule is 1.2 Å, but for the sake of clarify in the present figure an elongated molecule is considered with C-O bonds of length 2.5 Å. These orbitals and their dynamics will be discussed in detail in the forthcoming chapters.

The system is excited with a strong laser pulse E(t) = -∂ t A x (t), where A x (t) is the potential vector given by

A x (t) = E 0 ω L sin 2 πt T cos(ω L t + φ) (2.25)
and hence the associated electric field is given by

E(t) = E 0 sin 2 πt T sin(ω L t + φ) - π ω L T sin 2πt T cos(ω L t + φ) , (2.26)
where E 0 ∝ √ I 0 is the field amplitude, T is the pulse duration, ω L = 2πc/λ L is the laser frequency and φ is the career envelope phase. In the present study we will choose the IR wavelength λ L = 2.1 µm and intensities around I 0 = 1 × 10 14 W/cm 2 . A typical electric pulse obtained from Eq.(2.26) is shown Fig. 2.6. The choice made on the form of the electric field given by Eq. (2.26) ensures that the vector potential goes to zero when the field vanishes at time t = T , in agreement with Maxwell's equations. This is important since the asymptotic analysis is build upon the value of the vector potential A x (T ) at the end of the pulse [33,34]. Non-zero values of the vector potential at t = T give rise to unphysical DC components, which can be unphysically present in the LIED spectrum as a global shift.

The pre-aligned target molecule is exposed to an intense laser pulse lasing at a wavelength of 2.1 µm with an intensity of the order of 10 14 W/cm 2 . A schematic representation of the ionization process in strong IR fields is given in Fig. 2.7. Fig. 2.8 shows typical ionization probabilities obtained for the HOMO and HOMO-1, calculated from the norm of the electronic wave function, together with the associated electric field. Both the HOMO and HOMO-1 ionization probabilities show a quite similar time dependence. There are here two jumps in the ionization signal followed by a plateau. This is because ionization of the system in intense infrared laser fields is occurring only over time windows during which the intensity reaches values large enough to tunnel ionize the system. The ionization bursts do not show up when the field is maximum, this is because it takes some time for the ionized electrons to reach the grid boundary where they will be counted as ionized.

The two main peaks of the exciting pulse trigger two ionization events which are imprinted in the ionization signal as a function of time. At the end of the pulse, the remaining continuum wave packet is left with a distribution of kinetic energies. This wave packet will take a certain time to reach the detector. Thus, the asymptotic analysis should be done for certain time even after the pulse. Once all high energy components reach the detector, the signal starts to form the second plateau which will keep increasing with a small slope due to the arrival of low kinetic energy electrons. Due to the presence of very small kinetic energy components in the wave packet, a negligible increment in the signal W/cm 2 . The electric field is given as a blue solid line and the corresponding ionization probabilities are shown as a red dotted line for the HOMO and as a red solid line for the HOMO-1.

will sustain forever. The convergence of the calculation can, hence, be checked by looking at the ionization probability versus time. Once the increment in the signal is small enough, i.e, smaller than a fixed threshold, the propagation and the asymptotic analysis can be stopped and the LIED (photo-electron) spectrum can be computed.

Conclusion and Outlook

In this chapter, the system and the numerical model that will be used for studying the dynamics of CO 2 in strong IR laser pulses are discussed in detail. The system considered here is an aligned carbon dioxide molecule. The strong field ionization and the associated elastic scattering of the ionized wave packet will be studied using a numerical solution of the time-dependent Schrödinger equation obtained from the split operator method. The states from which electrons are ionized are the HOMO and HOMO-1. These initial states are calculated using the imaginary time propagation technique. The code is optimized for the system and tested in terms of convergence.

The numerical results are discussed in the next chapter.

Introduction

In the previous chapter, all essential details were given for calculating the SAE dynamics of of a model CO 2 molecule. In the current chapter, the elastic scattering of the ionized wave packet is discussed.

In brief, the photoionization of CO 2 is calculated by propagating the initial wave function using the split operator method. The asymptotic part of the wave function is analyzed and the photoelectron spectrum is calculated. Photoelectrons obtained via strong infrared field excitation can be classified in two main categories: Directly ionized and rescattered electrons. These photoelectrons are carrying essential information about the system.

High energy electrons and high frequency electromagnetic waves are used widely for studying the structure of matter [1,2]. After a scattering event, the diffraction patterns drawn by the electrons or photons can be analyzed, and some structural information on the system can be extracted. Laser induced electron diffraction can be seen in the same perspective [3,4]. Unlike traditional scattering processes, in LIED the scattering beam of electrons is taken from the molecule itself. After ionization, the resulting outgoing waves carry information about the scattering centers. Just like in the case of photons, the electrons which are scattered from the ionic core can be considered as an image of the system in the reciprocal space. If it is possible to transfer this information back to the spatial coordinate space, one should gain an image of the system taken by the rescattering electrons.

In the case of a single polyatomic molecule, the recolliding process can simply be related to the well-known multiple-slit experiment. Here, the slits are replaced by the attraction due to the different nuclei in the molecule. In the case of CO 2 , there are a maximum of three potential slits corresponding to the three atoms. But the spectra obtained via LIED is complicated due to many physical processes involved throughout the formation of the spectra. These are discussed in the following sections.

Interference Patterns of the LIED Spectrum

LIED spectra will be also referred as photoelectron spectra. They give a picture of the momentum distribution of the ionized electron. A typical photoelectron spectrum obtained from the solution of the TDSE for the HOMO of CO 2 is given in Fig. 3.1 for a square pulse consisting of a single optical cycle. Panel (a) shows the spectrum in a linear scale while panel (b) shows the same spectrum using a logarithmic scale.

The spectrum is elongated in the k x direction along which the electric field is acting. Two successive ionization events create a continuum wave packet which is driven very far from the molecule. These two events happen along the direction of the field, giving photoelectrons distributed with momenta as shown in the figure. Two ionization events and a single recollision event occur over a single optical cycle and give rise to interferences between the photoelectrons in momentum space via two major types of ionization pathways: Directly ionized electron, and electrons ionized after a recollision event. The spectrum will be more complicated as the number of cycles increases. This is due to the repeated ionization and recollision events which can be seen as an overlay of many single optical cycle spectra. Thus it is interesting to understand the basic structures that can be found in the simplest LIED spectrum that we obtain with a single optical cycle. It is well-known from optical physics that the diffracted waves from an object will carry the image of the scattering centers in the reciprocal space. These images, when transformed back to the real space, give an image of the scattering object. This property of diffraction is transposable to matter waves and has been used in a wide variety of experiments and studies proving the correspondence of matter waves with optical waves. It has also been used for practical purposes [2,5,6]. These interference patterns carry information about the structure of the molecule [7][8][9]. Apart from multiple-slit like diffraction patterns, there are other interference patterns appearing in the spectra. Indeed, additional ring-like structures can be seen, which also come from the interference between the rescattered electron wave packets. These ring-like structures appearing in the photoelectron spectra (marked with faint dotted lines in Fig. 3.2) are due to the interference between the long and short trajectories of the recolliding electrons.

As shown in Fig. 1.4 and Fig. 1.5, for a specific energy we can find two trajectories: A long and a short trajectory which will necessarily interfere in the photoelectron spectrum. This results in the ring-like structures seen in the LIED spectrum [10].

Another interesting interference in the LIED spectrum is due to the superposition of the direct and recollision pathways of the strong field process. This holographic interference of the electron wave occurs only over a small window of k y , due to the limited spread of directly ionized electrons in the transverse direction [11][12][13]. To see distinct structures due to this particular type of interference, it is worth to use, in the numerical simulation, a hypothetical electric field which consists of a half-cycle pulse with very high intensity and an associated counter part with very low amplitude, sustaining for a longer time [14]. Such a pulse does not lead to recollision. Mathematically, it can be designed by taking an appropriate time-dependent envelop. The long half-cycle with very small amplitude will impose a zero potential vector at the end of the pulse. The strong field ionization will be triggered only when the intense half-cycle interacts with the system. Being very weak, the second half of the optical cycle will be unable to brings back the ionized electron to the core. The spectrum obtained will be a spectrum build up of directly ionized electrons only. These interference patterns are widely used for extracting information about molecules. LIED was suggested as a tool for extracting the structural information by T. Zuo, A. D. Bandrauk and P. B. Corkum in [3]. Later on, it was demonstrated both theoretically and experimentally that LIED is a reliable tool for retrieving structural information [4, 8, 9, 11-13, 15, 16].

An important comment on the LIED spectrum is that the final energies are not equal to the recollision energy. The laser field keeps accelerating the electron after the recollision. This additional velocity gained by the photoelectrons is playing a role in the interference between different electron waves recolliding at different times with different energies. For the problem considered here which is aimed at retrieving the information about the molecule from the LIED spectrum, out of different interferences in this spectrum, the multiple-slit type interference is a rather good candidate, because it is at the origin of the robust interference structure seen in the k y momentum distribution which is mainly due to the rescattering with the ionic core. The most interesting information about the system, including the position of the nuclei and the orbital of origin of the photoelectron, is imprinted in the k y momentum distribution, along the molecular axis.

Thus, to get a simpler spectrum for analyzing these details, the least interesting part of the spectrum, that is also affected by many different interference processes and distortions induced by the laser field itself, is averaged out so that a simplified 1D photoelectron spectrum S(k y ) is obtained. Averaging out the k x component of the 2D spectrum S(k x , k y ) yields the averaged 1D spectrum

S(k y ) = ∞ -∞ S(k x , k y ) dk x . (3.1)
One has to keep in mind that there are some cases where the spectrum is not resolved enough for extracting information. This is the case when U p is too small, leading to small recolliding energies. Low energy electrons are characterized by large de Broglie wavelengths which are inappropriate for resolving the short distances in the diffraction process. This problem can be avoided to a certain extend by using low frequency fields and high intensities so that the recolliding energy increases. The dependence of the LIED spectrum on different laser parameters is analyzed in the next section.

Dependence of the LIED on various Laser Parameters

LIED, as its name is indicating, has a strong dependence on the laser parameters including intensity, wavelength, carrier envelop phase (CEP), pulse duration and state of polarization (ellipticity). For instance, there are many works discussing strong field ionization in elliptically or circularly polarized fields [17,18]. For our study, linearly polarized fields of intensities around 10 14 W/cm -2 are used. In the following subsections we explore the dependence of the LIED spectrum on the wavelength, pulse duration and CEP.

Dependence of the LIED Spectrum on the Wavelength

The quasi-static picture of laser induced ionization holds well for infrared frequencies. But the dependence of the recollision energy on the wavelength of the applied field has a strong influence on the outcome of the process. Being proportional to the square of the wavelength, the recollision energy increases quickly with longer wavelengths. In the case of LIED, gaining more energy while in the continuum will increase the resolution of the photoelectron spectrum because the recolliding electron is then characterized by a smaller de Broglie wavelength. Interference patterns will be more pronounced and will be extended to higher kinetic energies. The interference structures are well pronounced at 1.3µm and 2.1µm. At 800 nm, the interference patterns are more difficult to read. Averaged 1D spectra S(k y ) obtained from Eq. (3.1) are plotted in Fig. 3.5. Being primarily due to the recollision, the distribution of the photoelectrons in k y shows a dependence on the wavelength. The blue solid line is the averaged spectrum at 800 nm. The red and green solid curves show the same for 1.3 and 2.1 µm respectively.

- Here, λ = 2.1 µm will be chosen in order to have a better resolution and in order to provide a broader range of energies where the LIED spectrum can be analyzed. This is of course important for a reliable analysis.

2 -1 0 1 2 -4 -3 -2 -1 0 k y (a.u.) S(k y ) (log. scale)

Dependence of the LIED Spectrum on the Pulse Duration

For few optical cycle pulses, the LIED spectrum has a strong dependence on the pulse duration. An important and evident change in the LIED spectrum is the shape of the spectrum itself. With a single optical cycle pulse, the spectrum has a shape which indicates the directionality of the process of ionization and associated recollision. Fig. 3.6 shows the LIED spectrum S(k x , k y ) obtained with a single optical cycle pulse of λ = 2.1 µm from the HOMO at the internuclear distance R = 3 Å.

In a single optical cycle pulse, there are two ionization events and a recollision event. Ionized electrons from the first half cycle will be recolliding with the core on the phase change of the field. While the recollision is occurring, a second ionization event is also launched. This creates an asymmetry in the spectrum. In addition, the maximum kinetic energy that recolliding electrons can achieve is 3.17U p . Thus in the momentum space, the photoelectron spectrum will drop quickly to zero beyond this energy limit. The momentum distribution should therefore be mainly contained inside a circle of radius 3.17 U p , as seen in Fig. 3.6. This circular shape associated with the asymmetry discussed above are the main reasons behind the peculiar brain-shape of this photoelectron spectrum. This shape will of course be modified by the action of a multi-cycle field. The asymmetry in the spectrum, due to unbalanced ionization and recollision events, gradually disappears in the case of multi-cycle pulses. The ionization process thus becomes more and symmetric with large pulse durations.

An important assumption made on modeling the system was on the nuclear degrees of freedom of the molecule. Nuclear coordinates participate in both rotational and vibrational motions of the molecule. Rotational degrees of freedom will be out of the time scales considered in the present case and they are neglected. The vibrational degrees of freedom of an isolated CO 2 molecule are characterized by three vibrational modes: A symmetric, an anti-symmetric and a bending mode. Out of these, the anti-symmetric stretch has the shortest period, around 15 fs. For the pulses used in this calculation, the 4-optical cycle pulse, which lasts for 10 fs (FWHM), is already a bit long but one can also note that the anti-symmetric stretch can only be excited if the symmetry of the system is broken, but this is not the case here.

Dependence of the LIED Spectrum on the CEP

In general, the phase of any oscillating function can be written as

θ(t) = ω L t + φ , (3.2)
where ω L is the frequency of the oscillation and φ is a constant phase. In the case of pulsed oscillations of long durations, i.e. oscillations with a slowly varying envelop, this additional constant φ does not affect the oscillatory characteristics, but shifts the positions of the maxima and minima in the pulse. This phase φ is known as the carrier envelop phase (CEP). In other words, in the case of a sinus pulse, φ is translating the oscillation smoothly from a sine to a cosine wave. In addition, since strong field processes are highly nonlinear, the CEP becomes an important parameter while using ultra-short pulses. Fig. 3.8 shows typical single optical cycle pulses for the CEPs φ = 0 and φ = π/2. We clearly see that the blue curve corresponding to the CEP φ = 0 shows two identical maxima pointing in opposite directions. On the other hand, the red curve corresponding to the CEP φ = π/2 shows the presence of a main negative peak. One can thus expect that the associated LIED spectra will be very different.

Fig. 3.9 depicts these LIED spectra. The spectrum obtained for φ = 0 with two identical maxima pointing in opposite directions has more high kinetic energy electrons than what is obtained when the recolliding mechanism is absent as for φ = π/2 which shows the presence of a main negative peak in the associated electric field.

-4 -2 Panel (c) of Fig. 3.9 shows the corresponding 1D spectra for these two cases. The blue solid curve corresponds panel (a) and the red solid curve to panel (b). The k y distribution of the photoelectrons dies of very quickly for the cosine pulse. Again, this is due to the weak recollision probability. In the case of a sine wave, the recollision process enhances the signal for large values of k y .

0 2 4 -4 -2 0 2 4 k x (a.u.) k y (a.u.) (a) -4 -2 0 2 4 k x (a.u.) (b) -2 -1 0 1 2 -4 -3 -2 -1 0 k y (a.u.) S(k y ) (log scale) (c)
The dependence of the LIED spectrum with the CEP is much less dramatic with more optical cycles because in this case a larger number of ionization-followed-by-recollision events takes place. To illustrate this effect, Fig. 3.10 shows the average 1D spectra S(k y ) for the CEPs φ = 0 and φ = π/2 with a 4-optical cycle pulse. As it can be seen, the two spectra coincide rather well when compared with Fig. 3.9 (c).

In the forthcoming discussions, for such longer pulses, we will chose a CEP fixed at φ = 0.

LIED Spectra of the HOMO and HOMO-1

All the discussions done in the previous sections are for the HOMO of the molecule. It is interesting to see the modifications in the LIED spectra due to the initial state of the photoelectrons. Indeed, the key ideas for an attempt to extract the orbital information from the LIED spectra must rely in this modification.

In the LCAO approximation, the HOMO of CO 2 , an antisymmetric orbital, is mainly formed by mixing the 2p x orbitals of the oxygen atoms. The zero contribution of the carbon atom orbitals arises from the presence of a nodal plane perpendicular to the molecular axis (see Figs. 2.2 and 2.5 of Chapter 2). In the case of an aligned molecule as considered here, this nodal plane is parallel to the polarization direction of the field. This anti-symmetry implies that the electron necessarily leaves the molecule with a non-zero k y component of its momentum. Thus it indicates the presence of a nodal plane parallel to the field polarization in the reciprocal space of (k x , k y ). A signature of this symmetry information can be seen in the 2D spectra S(k x , k y ) and in the 1D spectra S(k y ) of the HOMO (see Fig. 3.9 for instance). The antisymmetric nature of the initial state also gives a vanishing recollision probability at the origin, where the carbon atom is located. The mono-periodic interference structure seen for the HOMO along k y with a continuously decreasing amplitude (see Fig. 3.10 for instance) is therefore an indication of a double-slit type interference effect reminiscent of a standard Young's type experiment.

A different interference pattern will be obtained for the HOMO-1. This is because this orbital is a symmetric mixture of the 2p x orbitals of both oxygens and carbon. As explained in details in the next chapter, this symmetry leads to a doubly-periodic structure in the k y dependence of the LIED spectrum, as shown in Fig. 3.11. This doubly-periodic structure is reminiscent of a triple-slit Young's type experiment. The doubly-periodic structure is clearly seen in the successive periodic appearance of peaks of higher and lower amplitudes in the red curve of panel (c), associated with a 4-optical cycle pulsed ionization. The period between the higher peaks is the double of that between the higher and lower peaks. The later coincides with the fringe separation between the peaks obtained for the HOMO, indicating that this smaller inter-fringe separation corresponds to the distance between the two oxygen atoms. The larger separation is doubled compared to the previous one, and it therefore corresponds to the C-O bond length.

A noticeable difference is therefore seen in both the 2D and 1D photoelectron spectra of the HOMO-1 in comparison with the ones of the HOMO. This is due to the different symmetries of these orbitals and therefore to the additional contribution of the 2p x orbital of the carbon atom in the HOMO-1. With this qualitative understanding of the photoelectron LIED spectra, it is clear that these spectra should now be analyzed in more details, quantitatively, using an analytical model, hoping that the structural information about the molecular orbitals could be retrieved with a reasonable accuracy using such an analytical approach.

This will be done in the next chapter. A choice of pulse duration much smaller than the vibrational timescale will help to analyze the system within the frozen nuclei approximation. Such an approximate analysis of the dependence of the LIED spectrum on the C-O bond length can be done with independent calculations at different internuclear distances. Because of the change in the position of the scattering centers, the spectrum will necessarily be affected.

Fig. 3.12 shows the 2D LIED spectra obtained for internuclear distances of 2, 3, 4 and 5 Å. The left column shows the photoelectron spectra of the HOMO and the right column shows the same for the HOMO-1. For both the HOMO and HOMO-1, the LIED spectra show a modification in the periodicity along k y which is connected to the arrangement of the nuclei in the molecule. As the bond length increases, the fringes become thinner and are getting closer. One should note that it has already been shown that for the case considered here, the bond length can directly be extracted from the interference patterns [7,19].

Introduction

The numerical calculations and some preliminary results obtained for the HOMO and HOMO-1 of aligned CO 2 molecules were discussed in Chapters 2 and 3.

The 2D LIED spectrum S(k x , k y ) calculated by solving the time-dependent Schrödinger equation contains some information about the molecule and therefore about the initial electronic wave function.

We will try to extract this information from the momentum distribution in the direction y parallel to the molecular axis since the momentum distribution along the direction x parallel to the field is dominated by the interaction with the strong linearly polarized field.

It is therefore of a lower interest for the imaging purpose we intend to develop.

Averaging S(k x , k y ) along the field direction provides a 1D spectrum S(k y ) from which the information will be extracted. But to extract this information, we need to derive an analytical model that will be compared, at some stage, with the calculated 1D spectra.

The model we develop is a very simplified approximate analytical model for LIED based on the strong field approximation (SFA). It is first developed in the general case of a linear molecule aligned in a direction normal to the polarization vector. It is then applied to the HOMO and HOMO-1 of the CO 2 molecule.

This model will finally be used to image the molecular orbitals at different internuclear distances to estimate if it is of interest for imaging a dynamical process such as the symmetric dissociation of the CO 2 molecule.

Building Blocks for the SFA Model

The CO 2 molecule is aligned perpendicular to the field. Ionization and the associated dynamics take predominantly place in the plane defined by the orthogonal system of coordinates consisting of the molecular axis and the polarization of the applied field. The field direction is taken as the x-axis and the molecular axis as the y-axis. The system is exposed to an intense IR electric field of the form

E = f (t) sin(ω L t + φ) x , ( 4.1)
where f (t) is the envelop of the pulse. As a first simplification, during the derivation of the SFA model we will consider that f (t) = E 0 is constant over an optical cycle. The vector potential associated to a single optical cycle is then given by

A(t) = E 0 ω L cos(ω L t) -1 x . (4.2)
The system evolves according to the TDSE and the state of the system |Ψ t at time t can be written formally as

|Ψ t = Û (t ← t i ) |Ψ t i , ( 4.3)
where |Ψ t i is the initial state and Û (t ← t i ) is the evolution operator obeying the TDSE

i ∂ t Û (t ← t i ) = Ĥ(t) Û (t ← t i ) , ( 4.4)
where Ĥ(t) is the system Hamiltonian given in Eq. (2.9). It contains a kinetic energy part and a time dependent potential part. Depending on the situation, one of them could be more influential than the other and could decide for the outcome of the dynamical process [1]. To take such an idea into account, we partition the Hamiltonian Ĥ(t) in two parts as

Ĥ(t) = Ĥ1 (t) + Ĥ2 (t) , ( 4.5)
without specifying, at the moment, which parts of the Hamiltonian are included in Ĥ1 (t) and which parts are included in Ĥ2 (t). It is important to note that neither Ĥ1 nor Ĥ2 (t) stand necessarily for the field free Hamiltonian Ĥ0 (t). They can stand for any part of the total Hamiltonian Ĥ(t). In practice, one of these two Hamiltonians, say Ĥ2 (t) for instance, is usually chosen as the part of Ĥ(t) which is responsible for the particular physical process under study.

With such an arbitrary partition of the Hamiltonian, the total evolution operator Û (t ← t i ) can be related to the evolution operator Û1 (t ← t i ) associated with Ĥ1 (t) and to Ĥ2 (t)

via the Dyson equation [2] Û (t ← t i ) = Û1 (t ← t i ) -i t t i Û (t ← t ) Ĥ2 (t ) Û1 (t ← t i ) dt . (4.6)
It is easy to prove that this equation is exact simply by verifying that this formal expression is a solution Eq. (4.4). Here, Û1 (t ← t i ) is the evolution operator obeying

i ∂ t Û1 (t ← t i ) = Ĥ1 (t) Û1 (t ← t i ) . (4.7)

Formal Exact Solution of the TDSE

Eq. (4.6) gives a formal exact solution of the TDSE. Following Eq. (2.9), the system Hamiltonian can be written as

Ĥ(t) = -∇2 /2 + Vc (r) -μ.E(t) , ( 4.8)
where ∇ = i k is the gradient operator and Vc (r) is the Coulomb or soft-Coulomb SAE potential of the molecule given in Eq. (2.5). The last term denotes the length gauge coupling with the field written within the dipole approximation.

The simplest realistic picture of strong field ionization including the essential ingredients of tunnel ionization followed by recollision, requires to consider at least a complete optical cycle. For the derivation of the model we limit ourself to a single optical cycle.

We now specify the partition we use to describe the tunnel ionization process as [3,4] Ĥ1

≡ Ĥ0 = k2 2 + Vc (r) (4.9a) Ĥ2 (t) ≡ -μ.E(t) (4.9b)
Under such a partition the Dyson equation (4.6) can be written as

Û (t ← t i ) = Û0 (t ← t i ) -i t t i Û (t ← t ) -μ.E(t ) Û0 (t ← t i ) dt . (4.10)
The physical interpretation of this equation is simple: Starting from the right, the system evolves field free and accumulates a phase from the initial time t i until the excitation by the field at time t . The resulting ionization process forms an electron wave packet propagating in the field under the total Hamiltonian until an arbitrary final time t. The integral over t simply means that the ionization process can take place at any time between t i and t.

The Dyson equation (4.10) involves the total evolution operator (also called propagator in the present study) Û (t ← t ) between the time of ionization t and the final time t.

During this time interval a recollision event may take place. In order to emphasize this recollision event in the expression of Û (t ← t ), we now adopt another partition of the total Hamiltonian as

Ĥ1 (t) ≡ Ĥv (t) = k2 2 -μ.E(t) (4.11a) Ĥ2 ≡ Vc (r) (4.11b)
In this case Ĥ1 (t) ≡ Ĥv (t) is the Volkov Hamiltonian that we have already introduced in Eq. (1.1). Using this new partition, the evolution operator Û (t ← t ) can be written as

Û (t ← t ) = Ûv (t ← t ) -i t t Û (t ← t ) Vc (r) Ûv (t ← t ) dt , (4.12)
where Ûv (t ← t ) is the propagator associated with Ĥv (t). It verifies

i ∂ t Ûv (t ← t ) = Ĥv (t) Ûv (t ← t ) . ( 4.13)
Substituting Eq.(4.12) in Eq.(4.10) we get

Û (t ← t i ) = Û0 (t ← t i ) + Ûd (t ← t i ) + Ûr (t ← t i ) (4.14)
where

Û0 (t ← t i ) = exp -i Ĥ0 (t -t i ) (4.15a) Ûd (t ← t i ) = -i t t i dt Ûv (t ← t ) -μ.E(t ) Û0 (t ← t i ) (4.15b) Ûr (t ← t i ) = - t t i dt t t dt Û (t ← t ) Vc (r) Ûv (t ← t ) -μ.E(t ) Û0 (t ← t i )(4.15c)
In this expression Û0 (t ← t i ) is the field free evolution of the system, Ûd (t

←

Exact Transition Amplitude

Once the formal solution of the problem is established, the transition amplitude between the state at the end of the pulse and the continuum can be calculated by projecting the formal solution given by Eqs. (4.16) and (4.17) at the end of the pulse (t = t f ) on the field-free ionized states. Let |Ψ + k be the elastically scattered outgoing waves in the direction of the electron wave vector k = (k, θ k ) for a prescribed asymptotic kinetic energy k 2 /2. The transition amplitude can be calculated as

a(k) = Ψ + k |Ψ t f (4.18a) = Ψ + k |Ψ 0 t + Ψ + k |Ψ d t + Ψ + k |Ψ r t (4.18b)
The dynamics is starting at the initial time t i from an eigenstate |Φ 0 of the field-free Hamiltonian Ĥ0 , of energy -I p where I p is the ionization potential. Therefore

|Ψ 0 t f = exp i I p (t f -t i ) |Φ 0 . (4.19)
Since the ionized states |Ψ + k are orthogonal to |Φ 0 , the transition amplitude reduces to a sum of two terms

a(k) = a d (k) + a r (k) (4.20)
where

a d (k) = Ψ + k |Ψ d t f (4.21a) a r (k) = Ψ + k |Ψ r t f (4.21b) a d (k)
is the transition amplitude associated to direct ionization and a r (k) includes the recollision process. Hence, the energy-resolved and angle-resolved ionization probability, or in other words the 2D LIED spectrum, can be written as

S(k x , k y ) = a d (k) 2 + 2 Re a * d (k)a r (k) + a r (k) 2 (4.22)
where

a d (k) = -i t f t i Ψ + k | Ûv (t f ← t ) -μ.E(t ) Û0 (t ← t i ) |Φ 0 dt , ( 4.23)
and

a r (k) = - t f t i dt t f t dt Ψ + k | Û (t f ← t ) Vc Ûv (t ← t ) × -μ.E(t ) Û0 (t ← t i ) |Φ 0 . (4.24)
These formal expressions where the total evolution operator Û (t f ← t ) is still involved, are exact for the calculation of the 2D LIED spectrum S(k x , k y ).

Approximate Transition Amplitude

To take into account recollision processes one needs to evaluate the amplitude a r (k) given in Eq. (4.24). However, because of the appearance of Û (t f ← t ) in this expression, it appears that evaluating a r (k) is as difficult as solving the TDSE. To render this evaluation simpler, an approximation should be made. For this purpose, we use the Strong Field Approximation (SFA) [4,5]. At this point, the evolution operator Û (t ← t ) in Eq. (4.24) is replaced by the Volkov evolution operator Ûv (t ← t ). Thus a r (k) -

t f t i dt t f t dt Ψ + k | Ûv (t f ← t ) Vc Ûv (t ← t ) × -μ.E(t ) Û0 (t ← t i ) |Φ 0 . (4.25)
Replacing Û (t ← t ) by Ûv (t ← t ) in Eq. (4.24) means that after the first recollision event we neglect the Coulomb interaction compared to the interacting with the strong IR field. It means that we are not taking into account any recollision process beyond the first one. If a second recollision takes place the approximation should break down. These events are however relatively rare and the approximation seems relatively reasonable. In addition, the wave packets are analyzed in the asymptotic region where SFA should hold within a reasonable accuracy in comparison with the exact treatment of the problem.

The second important approximation we are going to make is that the outgoing waves |Ψ + k are approximated by simple plane waves |Φ pw k , whose coordinate representation is

Φ pw k (r) = r|Φ pw k = e ik.r 2π . (4.26)
This approximation is obviously valid at large distances where the Coulomb interaction is weak compared to the kinetic energy of the wave packet but it is not very accurate in the vicinity of the ionic core. The reason behind this approximation is that the action of the Volkov evolution operator Ûv (t ← t ) on a plane wave can be evaluated analytically using an expansion in a basis of Volkov states {|Φ v k (t) } (see Eq. (1.1)). Within an arbitrary time interval [t 1 , t 2 ], the Volkov propagator can be written as

Ûv (t 2 ← t 1 ) = dk |Φ v k (t 2 ) Φ v k (t 1 )| (4.27)
where

r|Φ v k (t) = Φ v k (r, t) = e -i S(k,t) e i [k+
Φ pw k | Ûv (t 2 ← t 1 ) = e -i S(k,t 1 ,t 2 ) Φ pw k | (4.30)
where

k = k -A(t 2 ) + A(t 1 ) (4.31) and S(k, t 1 , t 2 ) = 1 2 t 2 t 1 [k -A(t 2 ) + A(τ )] 2 dτ . ( 4
a d (k) -i t f t i e -i S(k,t ,t f )-Ip(t -t i ) Φ pw k | -μ.E(t )|Φ 0 dt (4.33)
where

k = k -A(t f ) + A(t ) (4.34)
and to a r (k) -

t f t i dt t f t dt e -i S(k,t ,t f )-Ip(t -t i ) × Φ pw k | Vc Ûv (t ← t ) [-μ.E(t )] |Φ 0 (4.35)
where

k = k -A(t f ) + A(t ) (4.36)
In a realistic case, the potential vector A(t f ) will be zero at the end of the pulse, thus simplifying Eqs. (4.34) and (4.36). In addition, the electrons with a large recollision energy are ionized around t 0.3 (2π/ω L ), near the maximum of the field, as explained in section 1.2.3. Thus the value of the potential vector A(t ) at the time t of ionization will be relatively small for these high kinetic energy electrons and we are going to neglect this quantity in Eq. (4.34). The same high kinetic energy electron recollides with the ionic core at a time t close to a minimum of the field, corresponding to a maximum of the vector potential A(t ) ±E 0 /ω L x [6,7]. We thus finally end up with

k k (4.37a) k k ± E 0 /ω L x (4.37b)
The potential vector A(t ) therefore induces a strong shift on the parallel component k x of the electron momentum. Since our aim is to derive an approximate expression for the 1D spectrum along k y we will safely assume that this shift is of no real importance since it will be averaged out in the calculation of S(ky) and we will therefore not take it into account in the following.

To avoid large discrepancies between the approximate SFA spectrum and the spectrum obtained from the solution of the TDSE we should keep in mind that we should restrict the analysis to the highest electron momentum components. Moreover, it is clear from the numerical calculations that the multiple-slit interference patterns are not well resolved in the LIED spectra for slow electrons. These approximations allow us to separate the temporal from the spatial integral in Eq. (4.33) as

a d (k) A d Φ pw k | x |Φ 0 (4.38)
where

A d = -i t f t i E(t ) e -i S(k,t ,t f )-Ip(t -t i ) dt (4.39)
Similarly, the recollision amplitude becomes a r (k) -

t f t i dt t f t dt E(t ) e -i S(k,t ,t f )-Ip(t -t i ) Φ pw k | Vc Ûv (t ← t ) x |Φ 0 (4.40)
The completion of the general approximate formula for the calculation of a r (k) relies on the evaluation of the integrals shown in (4.40). This calculation requires the knowledge of the wave function at the time of recollision

|Φ rec (t ) = Ûv (t ← t ) x |Φ 0 (4.41)
This wave function can be evaluated by making use of the closure relation of the complete basis set of plane waves

{|Φ pw k } |Φ rec (t ) = dk Φ pw k | Ûv (t ← t ) x |Φ 0 |Φ pw k (4.42a) dk e -i k 2 (t -t )/2 Φ pw k | x |Φ 0 |Φ pw k (a) (4.42b)
where we have again dropped the momentum shift along the parallel component k x since it will be averaged out in the calculation of S(ky). Finally, as seen previously, high kinetic energy electrons have an excursion time in the continuum between ionization and recollision which lasts for about t -t = ∆t 0.7 (2π/ω L ) [6,7]. We will therefore approximate the recolliding wave function as

|Φ rec dk e -i k 2 ∆t/2 Φ pw k | x |Φ 0 |Φ pw k (b) (4.43)
an expression which does not depend anymore on t .

The temporal and spatial integrals can then be separated in the expression (4.40) of the recollision amplitude as

a r (k) A r Φ pw k | Vc |Φ rec (4.44)
where (a) This integral can also be evaluated in the coordinate representation, yielding the convolution

Φ rec (r, t ) ∝ dr e i (r-r ) 2 2(t -t ) x Φ 0 (r ) .
(b) Note that this evolution is reminiscent from a free wave packet evolution during the time interval ∆t. This is because the momentum shift in the k x direction was dropped.

A r = -

t f t i dt t f t dt E(t ) e -i S(k,t ,t f )-Ip(t -t i ) (4.45)
Thus the approximate, but still general expression for total transition amplitude is given by

a(k) = A d Φ pw k | x |Φ 0 + A r Φ pw k | Vc |Φ rec (4.46)
and the approximate 2D photoelectron spectrum can be calculated as

S(k x , k y ) = A d Φ pw k | x |Φ 0 + A r Φ pw k | Vc |Φ rec 2 (4.47)
It is possible to calculate the coefficients A d and A r via the saddle point method. We will not perform this calculation because we would rather prefer to use the coefficients as two adjustable parameters that we will be able to tune in order to match the approximate formula (4.47) with the spectrum calculated from the TDSE.

Application of the SFA model on CO 2

The general expression obtained for the transition amplitude should now be applied to the specific case of CO 2 discussed so far. With the adapted expression of the LIED spectra for CO 2 we will be able to compare the approximate spectra with the ones obtained from the numerical solution of the TDSE in order to verify the validity of the model.

A choice of initial state

We now need to specify more precisely the initial state |Φ 0 of the system. In quantum chemistry, there are many basis set ansatz for representing localized wave functions like Gaussian-or Slater -type orbitals, which can facilitate the calculation of bound states. Here we will choose a similar but simplified approach. The HOMO and HOMO-1 of CO 2 were presented in Chapter 2. Here, appropriate functional forms will be chosen that can represent these wave functions with a reasonable precision.

In order to obtain a simple analytical form of the 1D ionization probability S(k y ) in the case of CO 2 , the initial state can be assumed, in a crude approximation, to be a simple combination of 2p x Slater type orbitals (STO), whose coordinate representation is

Φ s 2px (r) = r|Φ s 2px = N s r cos(θ r ) e -ζ r (4.48)
in two dimensions, with the normalization factor

N s = 8 3π ζ 2 , (4.49)
ζ being the so-called Slater orbital exponent. This analytical form, which constitutes a reasonable candidate for representing the orbitals, will be used for the present imaging and reconstruction purpose.

In quantum chemistry, STOs are sometimes used as basis set functions for developing molecular orbitals due to their similarity with the eigenfunctions of the hydrogen atom.

A linear combination of these atomic orbitals with the right symmetry, following the LCAO theory, can therefore provide a realistic approximate molecular orbital. One of the disadvantage of STOs lies however in the difficulty of the evaluation of multi-center integrals with such analytical forms. This was an important reason, historically, for the introduction of Gaussian-type orbitals (GTO). If we express CO 2 molecular orbitals as linear combinations of atomic orbitals centered on different nuclei, such integrals will show up in the evaluation of Eq. (4.44). We will therefore also use 2p x Gaussian type orbitals, whose coordinate representation is at short range, as show Fig. 4.2. We will use these two types of basis functions in the following in order to derive a simple analytical formula for S(k y ).

Φ g 2px (r) = r|Φ g 2px = N g r cos(θ r ) e -α

Ionization Amplitude

Here there are three atoms and two of them are identical. Thus, in general, the HOMO and HOMO-1 molecular orbitals can be represented using STOs or GTOs as the following linear combination

|Φ 0 = ξ o 1 |Φ - 2p x + ξ c |Φ 2p x + ξ o 2 |Φ + 2p x , ( 4.52)
where ξ o 1 and ξ o 2 are distinct weighting factors applied on the two oxygen atoms and ξ c is the weighting factor of carbon. These parameters are assumed to be real. They can be either positive or negative. |Φ 2p x denotes the 2p x orbital centered on carbon. |Φ - 2p x is the 2p x orbital centered on the first oxygen atom at r = -R and |Φ + 2p x is the 2p x orbital centered on the second oxygen atom at r = +R. We thus have

|Φ ± 2p x = T±R |Φ 2p x (4.53)
where T±R denotes the translation operator by the displacement vector ±R.

The expression for the transition amplitude a(k) given in Eq. (4.46) can be expanded on this initial guess.

a(k) = A d ξ o 1 Φ pw k | x |Φ - 2p x + ξ c Φ pw k | x |Φ 2p x + ξ o 2 Φ pw k | x |Φ + 2p x + A r ξ o 1 Φ pw k | Vc |Φ - rec + ξ c Φ pw k | Vc |Φ 0 rec + ξ o 2 Φ pw k | Vc |Φ + rec (4.54)
where, following Eq. (4.43), |Φ 0 rec and |Φ ± rec are given by

|Φ 0 rec = dk e -i k 2 ∆t/2 Φ pw k | x |Φ 2p x |Φ pw k (4.55a) |Φ ± rec = dk e -i k 2 ∆t/2 Φ pw k | x |Φ ± 2p x |Φ pw k (4.55b)
The above general expression for a(k), Eq. (4.54), can be further simplified for the aligned CO 2 molecule. Looking at the properties of the integrals in Eq. (4.54), the first three integrals represent direct ionization from displaced orbitals and the last three are the ionization amplitudes after a recollision event. From this expression, it appears that the three direct ionization amplitudes are nothing but Fourier transforms of products of the dipole with displaced 2p x orbitals. In momentum space, this spatial translation becomes a simple phase shift of the form e ±ikyR . Taking this simplification into account, Eq. (4.54) can be reduced to

a(k) = A d ξ o 1 e i kyR + ξ c + ξ o 2 e -i kyR Φ pw k | x |Φ 2p x + A r ξ o 1 Φ pw k | Vc |Φ - rec + ξ c Φ pw k | Vc |Φ 0 rec + ξ o 2 Φ pw k | Vc |Φ + rec (4.56)
Compared to Eq. (4.54), the number of integrals that we have to evaluate is reduced to four. Out of them, the last three are very similar. The evaluation of the first integral using a GTO orbital yields

Φ pw k | x |Φ g 2px ∝ (k 2 x -2α) exp - k 2 x + k 2 y 4α , ( 4.57)
and the direct ionization amplitude using a GTO becomes

a g d (k) = A d ξ o 1 e i kyR + ξ c + ξ o 2 e -i kyR (k 2 x -2α) exp - k 2 x + k 2 y 4α , (4.58)
where the proportionality factor of Eq. (4.57) has been included in A d .

The calculation of the recollision amplitude a r (k) is more involved since it needs the knowledge of the functional form of the recolliding wave functions |Φ - rec , |Φ 0 rec and |Φ + rec . Using a GTO initial state we obtain

Φ 0 rec (r) = r|Φ 0 rec ∝ α -iβ -2β 2 x 2 (α -iβ) 3 exp i αβ α -iβ r 2 , ( 4.59)
where β -1 = 2∆t in atomic units. From Eq. (4.55a) we see that the initial 2p x orbital used for the calculation of the wave function Φ 0 rec (r) is centered at the origin where the carbon atom is located. Eq. (4.59) therefore gives the recolliding wave function for electrons emitted from the 2p x orbital of the carbon atom. The two other wave functions Φ ± rec (r) are the recolliding wave functions for electrons emitted from the 2p x orbitals of the two oxygen atoms. In a similar maner as what we have obtained for the direct ionization amplitude, these two additional wave functions are identical to the one of carbon apart from an additional phase shift. The total recolliding wave function is therefore

Φ rec (r) = Φ - rec (r) + Φ 0 rec (r) + Φ + rec (r) = ξ o 1 e i αβ α-iβ (R 2 +2yR) + ξ c + ξ o 2 e i αβ α-iβ (R 2 -2yR) Φ 0 rec (r) , (4.60)
Φ 0 rec (r) being given by Eq. (4.59). It is also interesting to note that in the near IR (800 nm to 2.5 µm) the parameter β is in the range 10 -2 to 10 -3 a.u. In comparison, the Gaussian orbital exponent α is usually of the order of 1 a.u. These orders of magnitudes will be used in following in order to simplify the expression of the ionization amplitude.

To evaluate precisely the amplitude a r (k) we have 9 integrals to calculate since ionization may originate from any of the three atoms and recollision may take place on any of the three atoms. Fortunately, these 9 integrals are similar. To understand the origin of these 9 terms, we can consider a simple picture, as shown in Fig. 4.3.

In general, all three atoms are contributing to the ionization signal. On recollision, the contribution from the first oxygen atom O 1 will scatter from the parent atom O 1 itself (shown as O 1 -O 1 in Fig. 4.3) as well as from the two neighboring atoms: from the carbon atom C (shown as O 1 -C) and from the second oxygen atom O 2 (shown as O 1 -O 2 ). This part of the rescattering amplitude shown in the figure can be written as

O 1 -O 1 O 1 -C O 1 -O 2 y x O 2 C O 1
a O 1 r (k) = Φ pw k | Vc |Φ - rec . (4.61)
where the potential operator Vc , supposed to be of a Coulomb form, is given in the coordinate representation by

V c (r) = - 1 |r + R| - 1 |r| - 1 |r -R| (4.62)
On the course of recollision one is interested by the scattering taking place in the vicinity of the three Coulombic cores. At first order near the singularities of the potential wells, i.e. for x → 0 and y → { -R, 0, R }, and taking into account Eqs. (4.59) and (4.60), the above integral can be reduced to

a O 1 r (k) ∝ - ξ o 1 |k y | e ikyR + e iβR 2 + e -ikyR e iβ4R 2 . (4.63)
Similarly, we obtain for the second oxygen atom

a O 2 r (k) = Φ pw k | Vc |Φ + rec ∝ - ξ o 2 |k y | e -ikyR + e iβR 2 + e ikyR e iβ4R 2 , (4.64)
and for the carbon atom

a C r (k) = Φ pw k | Vc |Φ 0 rec ∝ - ξ c |k y | 2 e iβR 2 cos(k y R) + 1 . (4.65)
Finally, the total recollision amplitude is

a r (k) = A r a O 1 r (k) + a C r (k) + a O 2 r (k) , (4.66)
where the proportionality factors have been included in A r . Note that within this approximation the recolliding amplitude a r (k) depends only on k y and not on k x . In this simplified case there is no need for integrating the recolliding amplitude over k x to obtain S(ky). It significantly simplifies our quest for a 1D LIED spectrum. We will now apply this model to the specific cases of the HOMO and HOMO-1 of CO 2 .

Approximate LIED Spectrum of the HOMO

The HOMO of CO 2 can be seen as a combination of two 2p x orbitals of opposite sign localized on the oxygen atoms of the molecule. This structure of the HOMO can be imitated by choosing the weighting parameters of the model as follows:

ξ c = 0, ξ o 1 = 1 and ξ o 2 = -1.
This choice will reduce the general expressions obtained for the direct and rescattering amplitudes significantly. The direct ionization amplitude becomes

a HOMO d (k) = A d (k 2 x -2α) e -k 2
x +k 2 y 4α sin(k y R) .

(4.67)

The recollision amplitude is reduced to

a HOMO r (k) = A r 1 -e i4βR 2 |k y | sin(k y R) , (4.68)
and the 1D spectrum can finally be written as The HOMO-1 of CO 2 has contributions from all atoms of the molecule. It can be seen as a combination of three 2p x orbitals of same sign localized on the oxygens and on the carbon. Hence the weighting parameters of the model are chosen as follows: ξ c = 1 and

S HOMO (ky) = |A d | 2 e -k 2 y 2α + |A r | 2 k 2 y sin 2 (k y R) . ( 4
ξ o 1 = ξ o 2 = ξ o = 0.
With this choice, the direct ionization probability is given by

a HOMO-1 d (k) = A d (k 2 x -2α) e -k 2
x +k 2 y 4α 2ξ o cos(k y R) + 1 .

(4.70)

The recolliding amplitude can be written as

a HOMO-1 r (k) = A r |k y | cos(k y R) e iβR 2 + ξ o e i4βR 2 + ξ o + 1 2 + ξ o e iβR 2 , (4.71)
and the 1D spectrum can finally be written as As expected, both spectra are dropping to zero as they go to higher kinetic energies. Globally, since they show the right behavior, the functions S HOMO (ky) and S HOMO-1 (ky) are promising candidates in terms of possibilities to find a set of parameters that could adequately reproduce the LIED spectra calculated from the TDSE.

S HOMO-1 (ky) = A 2 d e -k 2 y 2α 2ξ o cos(k y R) + 1 2 + A 2 r k 2 y cos(k y R) e iβR 2 +

Reconstruction of the Molecular Orbitals

Within the ansatz we have chosen for the initial state of the HOMO, there are two key parameters that are essential if we want to reconstruct the orbital from the 1D spectra: first the Gaussian orbital exponent α, and second the internuclear separation R. A third important parameters appears for the HOMO-1: the weighting parameter ξ o .

As we can figure out from the sin 2 (k y R) function of the approximate model for the HOMO, the internuclear distance can be extracted directly from the LIED spectrum by measuring the position of the fringes. The minima of the spectrum correspond to a destructive interference and they can hence can be equated to the momenta

k n = n π R , (4.73)
where n is the order of diffraction. Once the bond length R is obtained from the minima of the spectrum, it can be set in the model as a known parameter and we can then perform a fitting procedure in order to search for other parameters such as α or ξ o .

One should remember that the fringes will be profoundly imprinted in the LIED spectra only if the de Broglie wavelength of the ionized electrons is small enough to resolve the different scattering centers. One of the problems we will be facing is that the criterion set by the de Broglie wavelength of the electrons will become an important issue if the bond length is small. For example, in the case of CO 2 , the equilibrium value for R is about 1.16 Å. For this distance the fringe separation is ∆k = π/R 2.7 a.u. One has therefore to have recolliding electrons with momenta larger than 5.4 a.u. in order to see a minimum of two diffraction peaks. This corresponds to a recolliding kinetic energy of 14.6 a.u. If the tunnel ionization process is induced by a 2 µm radiation, such energies are obtained for laser intensities larger than 2.5 × 10 14 W/cm 2 . This is therefore relatively demanding in terms of laser performance if we take into account the fact that the reconstruction procedure also requires ultra-short pulses. In such cases of small internuclear separations, the extraction of R from the LIED spectrum is difficult, as we will see in the following sections.

In a first step, the reconstruction procedure will be shown considering an elongated CO 2 molecule. We will then check wether the method could be used as an imaging process for the dissociation of CO 2 .

The HOMO of an elongated molecule

Let us first consider an elongated CO 2 molecule, with R = 5 Å. This large value of R helps to obtain a well contrasted LIED spectrum. In order to quantify the quality of the reconstruction, we also compute the overlaps of the reconstructed HOMO orbitals with the initial state used in the TDSE calculation. They are defined as

O G = Φ g 0 |Φ T DSE 0 (4.74a) O S = Φ s 0 |Φ T DSE 0 (4.74b)
The values we have obtained here are O G 91.1 % and O S 94.7 %. Obviously, the Slater HOMO orbital has a better overlap with the TDSE initial state. This is not too surprising since this state shows a correct asymptotic behavior while the Gaussian HOMO orbital does not. As expected from the overlaps, the Gaussian orbital has more discrepancy compared to the Slater orbital. One can also note that, apart from the inaccuracy in the value of α, the small mismatch in the bond length R will shift the peaks of the reconstructed orbitals from the expected position by 0.12 Å. This creates a small distortion that one can see in panels With the satisfactory results obtained here for an elongated CO 2 molecule, the model will now be applied to more complicated situations. Here, in this demonstration of orbital reconstruction, the spectrum is well resolved with many maxima and minima that helped to extract information with a very good accuracy (see the TDSE spectrum shown in blue in Fig. 4.5). The position of the maxima and minima are the key for fitting accurately the analytical model with the numerical result. Once several maxima and minima are clearly seen, fitting all parameters is relatively easy.

Imaging the dissociation dynamics of CO 2

In the present section, the molecule is assumed to be following a dissociation dynamics. We suppose that the experiment is realized with ultrashort pulses such that the frozen nuclei approximation can still be used to obtain snapshots of the dynamics at different times, and therefore at different internuclear distances. The HOMO and HOMO-1 are leaving different signatures in the LIED spectrum and the model developed will be fitted to the spectrum, as demonstrated in Section 4.7.1. The fitting procedure will be performed for both the HOMO and HOMO-1 at internuclear distances R equal to 2.0, 3.0, 4.0 and 5.0 Å. In table 4.1, the internuclear distance obtained by the fitting procedure is denoted by R , while R is the exact value of the internuclear distance used in the TDSE calculation. The relative error achieved in the retrieval of the internuclear distance is also given. Fig. 4.9 finally shows the reconstructed GTO and STO orbitals along with the exact initial state. The color map used for these plots is the same as the one given in Fig. 4.6. The orbitals in the left column are the initial states used in the TDSE. Those in the middle column are the reconstructed orbitals using the Gaussian approximated initial states, and the right column contains the approximate Slater HOMO orbitals.

Imaging the HOMO orbital dynamics

The internuclear distances R are retrieved with very good accuracy, even from spectra with relatively small internuclear distances. The accuracy of R increases as the spectrum gets more and more resolved. The overlaps of the reconstructed orbitals are good and the Slater type orbitals, as expected, show better overlap with the TDSE initial states.

Imaging the HOMO-1 orbital dynamics

The HOMO-1 of CO 2 is an interesting case since, unlike the HOMO, it has contributions from all atoms of the molecule. The symmetry of this molecular orbital with the contribution from the carbon atom create a double period structure in the 1D LIED spectra (see Fig. 3.11). The photoelectron spectrum obtained from the model also possesses this double period structure (see the red curve in Fig. 4.4). This analytical model will be fitted to the numerical spectra at different internuclear distances for retrieving the set of parameters R , α and ξ o .

The additional parameter ξ o is a crucial quantity that determines how much is the contribution from the 2p x of the oxygen atoms compared to the 2p x of the carbon atom in the HOMO-1. Since it is a non-bonding orbital, the structure of the HOMO remains unchanged when the internuclear distance is varied, as shown in Fig. 4.9. This is no more the case for the HOMO-1, which is a bonding orbital whose dependence on the internuclear distance is strong. 4.2. This table gives the numerical values of all relevant quantities obtained from the fitting procedure. The reconstructed orbitals are finally given in Fig. 4.11. The color map used for these plots is the same as the one given in Fig. 4.6. The orbitals in the left column are the initial states used in the TDSE. Those in the middle column are the reconstructed orbitals using the Gaussian approximated initial states, and the right column contains the approximate Slater HOMO orbitals. As expected, the HOMO-1 is mainly made up of the atomic orbital of carbon and the contributions from the oxygens are relatively small. This contribution drops quickly as the bond length increases. Beyond 3 Å, the HOMO-1 is simply the 2p x orbital located on carbon. But the potential wells of the oxygens still play an important role during recollision, and leave their signature in the LIED spectra. If there was not such contributions from the oxygens, the 1D LIED spectra would not posses any periodicity. The analytical model we have derived is able to extract the correct information from the LIED spectrum, as shown in the reconstructed images of Fig. 4.11.

R α ζ ξ o R |R -R | R (%) O G (%) O S (%)

Conclusion

In this chapter, a simple analytical model for extracting the molecular information for imaging the orbitals was developed. The model deployed is based on the strong field approximation. It separates the direct and recolliding parts of the wave function in the evaluation of the ionization spectrum. The general model is applied to the HOMO and HOMO-1 orbitals of CO 2 . The orbitals of CO 2 are assumed as linear combinations of Gaussian-type orbitals for the evaluation of the integrals associated with the model. The Gaussian-type orbitals are then fitted with more realistic Slater-type orbitals to obtain an approximate form that can mimic the HOMO and HOMO-1 with a good accuracy.

The model is used to fit the results obtained from the TDSE calculations. The Gaussian and the Slater orbital exponents are extracted along with the bond length and other parameters.

The overlaps between the initial states used in the TDSE and the reconstructed orbitals are in the range 90 -95%, depending on the value of the internuclear distance. The bond lengths are extracted with a maximum error of 6%.

A slow dissociation of the CO 2 molecule can therefore be imaged using this model, which showed its ability to extract the orbital information over a wide range of internuclear distances (2 Å and beyond).

Conclusion and Outlook-Part I

This chapter concludes the Part I of this thesis.

Conclusion

In this part of the thesis, Laser Induced Electron Diffraction (LIED) from a CO 2 molecule initiated by an intense ultra-short linearly polarized infrared laser field was studied in detail.

The molecule was assumed to be initially pre-aligned in a direction perpendicular to the field polarization. Spectra of rescattered photoionized electrons were calculated by solving a 2D time dependent Schrödinger equation (TDSE) using the split operator method.

The diffraction patterns obtained via LIED imprint information about the molecular structure and orbitals which are of interest for imaging the molecular dynamics in real time.

The calculations were performed for different bond lengths and orientations of the system. Extraction of the bond length seems relatively accessible with current laser facilities for distances larger than 2 Å, from clear diffraction patterns observed in the LIED spectra.

In order to extract the information about the orbital, an analytical model was developed based on the Strong Field Approximation (SFA) using a Linear Combination of Atomic Orbitals (LCAO) as a guess initial molecular wave function.

The 2D LIED spectra obtained by solving the TDSE were averaged along the polarization direction of the applied field to obtain 1D photoelectron spectra.

The SFA model we have developed includes a limited number of adjusting parameters. These parameters were adjusted to match the TDSE 1D spectra in order to retrieve the geometrical and orbital parameters characterizing the initial molecular orbital.

The parameters obtained in this fitting procedure are in good agreement with the exact initial molecular orbital, proving that the approximation we have derived has the potential to be used for imaging purposes based on LIED.

This analytical model allowed us to extract R, the bond length, α, the Gaussian orbital exponent, and ξ o the relative contribution of the 2p x atomic orbitals of the oxygen atoms in comparison with the contribution from the 2p x atomic orbital of the carbon atom in the molecular orbitals we have analyzed: the HOMO and HOMO-1.

R could be extracted with an error that never exceeds 6%. The Gaussian orbital exponents that we obtained can be used to design Slater-type orbitals in order to obtain an accurate reconstructed molecular orbital. The reconstructed molecular orbitals are in good agreement with the exact initial states, with overlaps in the range of 90 to 95%.

Outlook

Imaging the dissociation dynamics of a molecular system is of wide interest in the community.

The main difficulty that is faced on this matter is the short time scale of typical molecular fragmentation processes.

The recollision process that occurs within a very few femtoseconds, hence, has the potential to reveal these dynamics which is hidden from traditional imaging techniques. Aside of the popular HHG based imaging techniques, LIED can also be an alternative candidate in the race for imaging fast reaction dynamics.

LIED spectra have a strong dependence on the geometry of the scattering centers and hence the nuclear motions can be related to modifications seen in the LIED spectra. Thus the mapping of time resolved LIED spectra to the nuclear dynamics will give the possibility for imaging the reaction dynamics with very good spatial and temporal resolutions.

However, the LIED reconstructing method still suffers from important limitations.

For a reliable reconstruction of the molecular orbitals, high energy photoelectrons are required. With high kinetic energy electrons, the diffraction patterns show clear signatures of the internuclear distance which are the key information to analyze the reaction dynamics.

In contrast to this inevitable need, typical LIED spectra are mainly built with directly ionized slow electrons and this reduces the quality of the spectrum significantly.

Another disadvantage of LIED based imaging techniques is the need for pre-aligned molecules. We are currently working on this topic, and we hope that in the near future the LIED reconstructing method could be extended to molecular systems which are not perfectly aligned. Introduction 111

Collective Effects 111 1.3 Towards New Directions 114

1.4

Introduction

In a clear contrast with the first part of this thesis, in which interactions of isolated atoms and molecules with intense laser fields and following dynamics were discussed, this second part deals with dense samples of atoms or molecules, that we will denote as quantum emitters, in the weak field regime.

The aim of this part is to understand the effect of dipole-dipole interactions in dense samples of quantum emitters coupled to a weak electromagnetic field. The action of the laser field will induce dipoles which will in turn act as local sources of electromagnetic fields oscillating at a frequency determined by the applied field and the collective nature of the medium. Dipole fields produced by the excited particles in the system will keep them coupled to each other. If the couplings are strong enough, the entire system will response cooperatively such that the field will not see the system as a collection of independent particles. The response of such systems towards the light will be modified dramatically via cooperative effects, which opens a broad regime with a lot of interesting practical applications, including slow light, storage of light, elements for light based memory devices and quantum information technologies, radiation transport in light-harvesting systems, sources of coherent radiation and even cooperative effect based lasers [1][2][3][4][5][6][7][8][9][10]. Investigations for designing attosecond sources based on high harmonic generations are nowadays extended to dense atomic and molecular systems and even to solids [11][12][13][14]. These diversities in interest of research and concomitant technical developments demand an intensive study on how quantum emitters interact with each other as they transit from weakly dense systems, where they behave almost independently, to dense systems, where they behave cooperatively. Research on the cooperative behavior of dense systems is hence a hot area owing to the fact that its complexity leads to the observation of new dynamics and better engineering of materials for future technologies. Before proceeding further, a brief account of collective effects is given below.

Collective Effects

The interaction of a laser field with an ensemble of collectively responding quantum emitters is a complex physical problem. These complexities are arising from dipole-dipole interactions in the system which lead to many well-known collective effects. These effects are profound when the average inter-particle separation is smaller than the transition wavelengths of the constituent particles. At this condition, the cooperativeness between quantum emitters will induce shifts in the resonance frequencies and other associated effects.

Lorentz-Lorenz (or LL) Shift

A detailed analysis of how a dense system respond to an electric field in the weak field regime dates back to 1888. In the phenomenal works of H. A. Lorentz and L. Lorenz, they introduced a simple and elegant way to correct the driving field of an arbitrary quantum emitter in a dense system which is now known as the Lorentz-Lorenz correction. It is a mean field correction used to account for electric fields from the induced sources in the system. It will induce, as it is known for coupled oscillators, a shift in the resonance frequency which is known as the Lorentz-Lorenz (or LL) shift

∆ LL = N µ 2 3 ε 0 , ( 1.1)
where N is the number density of the system and µ is the electric dipole. This shift in the resonance frequency of the system is clearly a signature that the quantum emitters in the system are well coupled via dipole-dipole interactions.

This predicted LL shift was observed experimentally by different research teams of the domain, validating the applicability of this mean field approach while determining the optical properties of such systems in normal conditions [15,16]. But in the case of cold atomic samples, the shift in the resonance frequency is still a matter of debate since there are some experimental results defying the models based on the mean field approach. The interesting fact here is that the expected mean field behavior is retrieved in such systems when some incoherence broadenings are added [17,18].

Cooperative Lamb Shift

An additional shift of resonance frequency in dense systems is the cooperative Lamb shift.

The standard Lamb shift is the shift introduced to the resonance frequency of an atom due to its interaction with the vacuum, through which the atom absorbs and emits virtual photons [19].

Similar to the interaction of an atom with the vacuum, quantum emitters in a dense sample can absorb and emit virtual photons among themselves, keeping the constituent emitters coupled, and resulting in a shift of the resonance frequency known as the collective or cooperative Lamb shift (CLS). Because of the exchange of virtual photons among the atoms or molecules in the system, the CLS is a geometry dependent shift which can be given as a function of the LL shift. To give an example, for a system of atoms confined in a rectangular slab the CLS is given by [20]

∆ CLS = 3 4 ∆ LL 1 - sin(2k ) 2k . (1.
2)

The existence of the CLS in dense systems was observed in various systems [21][22][23].

Dicke Narrowing

Robert H. Dicke predicted in 1952 that, if a system is prepared with a large number of quantum emitters, frequent coherence-preserving collisions could decrease the Doppler width of a spectral line of the system [24]. This narrowing of the spectral width, which is now known as the Dicke narrowing, occurs when the mean free path of an atom gets smaller than the resonance wavelength of the system. In a system that preserves coherence during collision processes, it leads to an averaging over different Doppler states which results in a narrower line width than the Doppler width. The spectrum S D of the system can then be written as

S D = Γ + ηΓ D ∆ 2 + (Γ + ηΓ D ) 2 , (1.3)
where Γ D is the Doppler width of the spectrum, ∆ is the detuning and η = 2πΛ/λ is the Dicke parameter, where Λ is the mean free path and λ is the wavelength. The understanding of Dicke narrowing is very important in spectroscopy and for analyzing the optical response of a many particle systems for processes like electromagnetically induced transparency-EIT [25][26][27][28][29] for instance.

Superradiance

In a dense sample of atoms where the average separation between the constituent particles is smaller than the radiative transition wavelength of the system, spontaneous emission can be enhanced via cooperative effects.

If an electric field excites a system prepared in the cooperative regime, the excited emitters will be coupled via dipole-dipole interactions. This interaction between particles will force them to emit photons collectively via spontaneous emission which will scale as the square of the number of total excited atoms in the system. This enhancement in spontaneous emission is know as superradiance [30][31][32]. Superradiance was observed in many different systems, see [33][34][35][36] for instance, and it continues to be a topic of investigation.

Towards New Directions

There are many works on these topics of collective effects giving insight to the laser-matter interaction in the cooperative regime [37,38]. Collective effects start to be used for developing new technologies which can revolutionize current technology. Cooperative effects in coherent and incoherent systems are studied extensively by many research teams [39][40][41]. Out of them, one of the recent developments was based on superradiance for designing a new laser system as suggested in [42]. In traditional laser systems using a resonance cavity to create a standing wave, the coherence is maintained by the photons trapped in this cavity. Because of the vibrations of the cavity mirrors, the monochromaticity of emitted photons is affected and will give a spread in the frequency around the resonance. But if one uses the superradiance in which the coherence is stored in the coupled excited atoms, the problem related to the vibrating mirrors is solved. In 2012, a group of scientists from the University of Colorado designed a laser based on superradiance [6,42] and there are many propositions and studies on this topic, hoping that superradiance can be efficiently used for lasing purposes [43,44]. The superradiant laser offers a highly stable and monochromatic source of light that could influence all fields of science, from defining time to detecting gravitational waves.

Other interesting developments in the field deal with the use of cooperative effects for designing essential elements for light storage, for light based memory devices or for making use of the ability of photons to transfer quantum information effectively [1,[45][46][47][48]. These are based on the slowing down of light inside a medium [START_REF] Safavi-Naeini | Electromagnetically induced transparency and slow light with optomechanics[END_REF]. Slow light, which was demonstrated both in cold and thermal systems, is caused by a dramatic increase in the group index of refraction of a dispersive system [START_REF] Hau | Light speed reduction to 17 metres per second in an ultracold atomic gas[END_REF]51] due to the nonlinearity. Besides slow light, there are experimental realizations on stopping and storing of light up to one minute and then retrieving the information stored, thus indicating the possibility for developing photon-based memory devices [52][53][54][55][56][57]. These developments, which are based on electromagnetic induced transparency, have however small fidelity due to the distortion of the photon [START_REF] Shapiro | Single-photon Kerr nonlinearities do not help quantum computation[END_REF][START_REF] Gea-Banacloche | Impossibility of large phase shifts via the giant Kerr effect with single-photon wave packets[END_REF][START_REF] Viswanathan | Multimode analysis of a conditional phase gate based on second-order nonlinearity[END_REF]. On this ground, the interest is now turned on improving the mechanism for storing photons for useful purposes [START_REF] Zheng | Waveguide-QED-Based Photonic Quantum Computation[END_REF][START_REF] Paredes-Barato | All-Optical Quantum Information Processing Using Rydberg Gates[END_REF].

Nonlinearity can thus be achieved by preparing the system in the cooperative regime. The purpose of this part of the PhD work was to investigate collective scattering of light from a dense sample of quantum emitters for understanding the response of such system towards the electromagnetic field in different interaction regimes. The collective response can be observed in any dense medium composed of homogeneous as well as heterogeneous quantum emitters. But the mixing of different types of quantum emitters with different physical properties and energy level structures can make the system extremely complicated. Thus, in this part, the system is first assumed to be composed of only one type of quantum emitters. These studies can be extended in a straightforward way to a mixture of more than one type of quantum emitters.

Introduction

This part of the thesis deals with a collection of millions of quantum emitters confined in a cell of nanometric size exposed to weak electromagnetic fields. In principle, a full description of such a system needs a complete quantum approach based on quantum electrodynamics. In this quantum picture, the electromagnetic field is a collection of discretized wave modes propagating through space: the photons. Even if the system we consider contains only a few number of atoms or molecules and a small number of photons, it is already difficult to treat the system by preserving the quantum description of light.

In practice, for a laser power of 1 mW and a photon energy of 1.5 eV, the approximate number of photons encountering the surface of the nano-layer per second can be estimated as about 10 11 . Treating this large number of photons individually becomes a herculean task. At this point, it is important to make an approximation in order to exclude the difficulty for treating the system exactly. Therefore, as we discussed in chapter 1, we use a semi-classical approach in which the material system is described quantum mechanically (in order to include the quantized behavior of the observables of interest), while explaining the light-mater interaction by keeping the electromagnetic field under the classical picture governed by Maxwell's equations.

There is another complementary approximation to the system by treating it as a statistical ensemble: quantities like position, velocity, or electric dipole are averaged to produce some continuous distribution functions of those quantities. This approach is justified when the particle density is not too weak. It is the case in our theoretical study since the dynamics we are interested in concerns the collective response of the quantum system which shows up for relatively large densities.

Maxwell's Equations in Source Free Isotropic Medium

The system under consideration is assumed to be non-magnetic and to have no free electric charges, so that χ m = 0, and ρ = 0 (see Chapter 1).

In that case

J = ∂P ∂t (2.1)
where P is the macroscopic polarization of the system induced by the field. In weak fields the polarization varies linearly with E P = ε 0 χ e E .

(2.2)

Maxwell's equations become

× E = -µ 0 ∂H ∂t (2.3a) × H = ∂P ∂t + ε 0 ∂E ∂t (2.3b
)

• E = 0 (2.3c) • H = 0 (2.3d)
Since the system is approximated as a continuous medium, the field obtained by solving Maxwell's equations is a macroscopic field. The wave equation for the macroscopic electric field can be written as

2 -ε 0 µ 0 ∂ 2 ∂t 2 E = µ 0 ∂ 2 P ∂t 2 .
(2.4)

Substituting for the polarization P and taking into account Eq.(2.4), in the frequency domain we obtain 2 Ẽ -ε 0 µ 0 ω 2 1 + χ e (ω) Ẽ = 0 (2.5)

where Ẽ(ω) denotes the Fourier transform of E(t). The factor ε 0 µ 0 [1 + χ e (ω)] can be equated to 1/v 2 , where v is the velocity of light through the medium. It can also be shown that in the time domain the complex functions

E ± (r, t) = E 0 e ∓i(k•r-ωt) ê , ( 2.6)
where ω 2 /v 2 = |k| 2 = k 2 and ê is a unit vector, are solutions of the wave equation whose Fourier transform is given in Eq.(2.5), k being the wave vector. The magnetic field can finaly be calculated from the electric field using Eq.(2.3a) as

H ± (r, t) = H 0 e ∓i(k•r-ωt) ĥ (2.7)
where ĥ = (k/k) × ê and H 0 = k E 0 /(µ 0 ω). The positive and negative signs on the fields represent the positive and negative frequency solutions respectively.

Note that the vectors ê, h and k form an orthonormal system. The modes of the electromagnetic field described here, which contain the transverse electric and magnetic fields (TEM modes), will be used for studying the dense system of quantum emitters. Note that the electric field corresponding to the wave equation can also be written as a linear combination of the positive and negative parts

E(r, t) = 1 2 E + (r, t) + E -(r, t) , (2.8a) H(r, t) = 1 2 H + (r, t) + H -(r, t) . (2.8b)
Note finaly that a general temporal shape of the propagating wave can be obtained by superimposing different spectral components together, i.e,

E + (r, t) = 1 2π ∞ 0 Ẽ+ (r, ω) e i(ωt-k•r) dω , ( 2.9a
)

H + (r, t) = 1 2π ∞ 0 H+ (r, ω) e i(ωt-k•r) dω .
(2.9b)

Correction to the Electric Field

The equations detailed here are sufficient to treat the propagation of the macroscopic electromagnetic fields. While doing the assumption of classical propagation of light through the medium, the system is considered as a continuous ensemble of quantum emitters. This smoothing of a system comprised of discretized constituent particles (atoms and/or molecules), approximates the system to a continuous medium as in the "Polarium" model of R. J. Glauber [1]. The averaging on the position of the particles also suggests that in the presence of an applied electromagnetic field, the polarization of the system can be written as a smooth (continuous) function that depends on the position vector r.

While solving the dynamical equations at any point inside the system, one has to account for the electric fields at the position of calculation due to the neighboring polarization elements. Since the averaging process removes the structural information of the system, the solutions of Maxwell's equations give only the macroscopic field, which is different from the microscopic field experienced by individual particles. The separation between the particles at the microscopic level modifies the field depending on the inter-particle distance and therefore on the number of particles per unit volume [2]. This mutual influence of the constituent particles can be included by doing a correction to the electric field locally. An elegant way to include this correction is via the Lorentz-Lorenz correction term that can be easily derived for a dense system as a function of the macroscopic polarization.

The electric field experienced by a constituent particle, the so-called local field E local (r, t), can be written as the sum of the macroscopic field E(r, t) and the internal field produced by the constituent dipoles in the system E int (r, t). It can be shown that [3,4] E int (r, t) = P(r, t) 3ε 0 .

(2.10)

Thus, inside the system

E local (r, t) = E(r, t) + P(r, t) 3ε 0 . ( 2

.11)

This correction to the macroscopic field is essential for explaining the collective macroscopic response of a dense system, including light scattering, reflection, absorption and transmission of incident electromagnetic waves. This correction and related formula were discovered by H. A. Lorentz [5] and independently by L. Lorenz [6] at the end of the 19 th century. Thus, the correction made on the local field is known as the Lorentz-Lorenz correction, or LL correction.

Energy flow and reflection or transmission spectra

Our ultimate aim is to study the response of a system of interest towards an applied electromagnetic field. From a fundamental point of view, the action of the electromagnetic field will turn on the energy flow into and out of the system. A part of this energy will be dissipated in the system due to decay effects and incoherent processes which can be estimated classically from the well-known Beer-Lambert's law [7]. In addition, energy transfer in the system directly relates to the excitation dynamics and hence by analyzing the energy flow, the dynamics of the system can be studied.

The electromagnetic energy flow can be represented by the Poynting vector, first proposed by J. H. Poynting [8], which gives the rate of energy flux per unit area of an electromagnetic field. The law of energy conservation is written

∂u ∂t = -• S (2.12)
where u(r, t) is the energy density of the electromagnetic field and S(r, t) is the Poynting vector, given by S(r, t) = E(r, t) × H(r, t) .

(2.13)

Being a vector, S(r, t) gives the energy flux at a particular point and in a particular direction. Thus, knowing the electric and magnetic fields, the reflection and transmission spectra characterizing the response of the system can be defined as the ratio

F(ω) = |S(r, ω)| |S 0 (r, ω)| (2.14)
where |S(r, ω)| is the magnitude of energy flux measured at a particular position and in a particular direction, and |S 0 (ω)| is the incident energy flux at the same position and direction.

Maxwell-Bloch Equations

The electromagnetic field equations given in section 2.2 are used to describe the field propagation through the system. Describing the system itself by retaining its quantum behavior alongside the classically described fields including the particle-field interaction forms the basic platform for performing the calculations to study the dynamics and response of the system. Quantum mechanically, the dynamics of a system under the influence of any potential is described by the time dependent Schrödinger equation (TDSE) or by the Liouville-von Neumann equation. From the fundamental point of view, a weak field acting on a quantum emitter can excite the system to a higher energy state which is resonant with the frequency of the field. First, it is better to understand the dynamics by considering the simplest case of a two-level system that can encapsulate the basic physics that we will explore later on with more complicated systems.

The two-level system

In a real physical world, there are many and sometimes infinite energy levels in a quantum emitter. In some cases a multilevel system can be approximated by a two level system that provides a simplified physical insight to its basics quantum behavior. An ideal monochromatic field acting on a system having an energy level resonant with the field can usually be treated as a two level system (see Fig. 2.1). A system that contains a large number of such atoms is of interest to explain the basic physics of dense samples of ideal two-level systems. Relatively large number of photons are also necessary to keep the electromagnetic field classical in practice. Since the system is assumed to have a very large number of interacting quantum emitters, it is advisable to approach the system via the density matrix formalism in order to include decay and decoherence processes associated, for instance, with atomic or molecular collisions.

Let ω 2 -ω 1 = ω 12 be the transition frequency of the two level system under consideration, where the frequency ω 1 corresponds to state |1 and ω 2 to state |2 (see Fig. 2.1). The field free Hamiltonian of the system can de written as

H 0 = ω 1 |1 1| + ω 2 |2 2| (2.15)
Assuming for simplicity that the applied electric field is monochromatic and linearly polarized along ŷ(a) , the interaction of the system with the applied field can be written using the electric dipole approximation as

H I = -µ• E = -|1 2| + |2 1| 2|µ ŷ|1 E 0 cos(ωt) (2.16)
where µ is the electric dipole and E 0 cos(ωt) is the ŷ component of the electric field given by Eq.(2.8a). µ ŷ is the ŷ component of the electric dipole and 2|µ ŷ|1 is the associated dipole matrix element. Note that we have assumed here that, for symmetry reasons, 1|µ ŷ|1 = 2|µ ŷ|2 = 0.

(a) i.e. ê = ŷ in Eq.(2.6).

The Rabi frequency, which can be used to define the frequency of the oscillations of the populations of states |1 and |2 due to the coupling with the field, can be written as

Ω 12 = -2|µ ŷ|1 E 0 (2.17)
and substituting it in Eq. (2.16) we obtain

H I = Ω 12 cos(ωt) |1 2| + |2 1| . (2.18)
Note that this expression simplifies within the rotating wave approximation (RWA) [9], where we only keep the slowly varying components in the interaction Hamiltonian

H RW A I = Ω 12 2 |1 2| e iωt + |2 1| e -iωt . (2.

19)

In order to introduce the density matrix formulation of the same two-state problem, let us now assume that the system is initially in a pure state with

|Ψ = C 1 |1 + C 2 |2 , ( 2
ρ 11 = C * 1 C 1 (2.22a) ρ 12 = C * 1 C 2 = ρ * 21 (2.22b) ρ 21 = C * 2 C 1 = ρ * 12
(2.22c)

ρ 22 = C * 2 C 2 (2.22d)
The evolution of the system can be evaluated using the von-Neumann equation, written in atomic units as

ρ = -i H, ρ (2.23)
where the dot symbol denotes the first order time derivative and where

H = H 0 + H I . (2.

24)

Writing explicitly the evolution equations for each element in the density matrix we get

ρ11 = i Ω 12 cos(ωt) (ρ 12 -ρ 21 ) (2.25a) ρ12 = i Ω 12 cos(ωt) (ρ 11 -ρ 22 ) + i ω 12 ρ 12 (2.25b) ρ21 = i Ω 12 cos(ωt) (ρ 22 -ρ 11 ) -i ω 12 ρ 21 (2.25c) ρ22 = i Ω 12 cos(ωt) (ρ 21 -ρ 12 ) (2.25d)
The equations (2.25a) and (2.25d) describe the evolution of the populations and the equations (2.25b) and (2.25c) show the evolution of the coherences. The two-level system described by these equations was assumed to have no decay. But in practice it is almost impossible to find such a system without any interaction with its surroundings (or reservoir).

In the cases where the coupling of the system with the reservoir is weak, so that any change in the reservoir due to the interaction with the system is negligible, the master equation which governs the evolution of the system can be written in the interaction picture as

ρ = -i H, ρ + Λ ρ (2.26)
where ρ now denotes the quantum emitter reduced density matrix which is calculated by taking the partial trace over the reservoir degrees of freedom of the density matrix associated with the total system { quantum emitter + reservoir }. The approximation performed here is known as the Born approximation.

The operator Λ contains the self correlations in the reservoir which carry the information regarding the interaction of the reservoir with the system, i.e, the memory regarding the interaction with the system. If the characteristic time scales over which the reservoir correlation function decays are much smaller than the characteristic time scales of the system, memory effects in the dynamics of reduced density matrix can be neglected. This approximation is known as the Markov approximation [10,11]. Now, the positivity of the reduced density matrix ensures the Lindblad form of the master equation, first introduced by V. Gorini, A. Kossakowski, E. Sudarshan [12,13] and independently by G. Lindblad [14].

Finally, the operator Λ contains all decay effects in the system due to the interaction with the reservoir. The non-diagonal elements of Λ include a pure dephasing rate γ and the diagonal elements account for the radiationless decay rate Γ of the excited state |2 . This leads to the following set of evolution equations of the system, including decay and dephasing [15,16]

ρ11 = i Ω 12 cos(ωt) (ρ 12 -ρ 21 ) + Γ ρ 22 (2.27a) ρ12 = i Ω 12 cos(ωt) (ρ 11 -ρ 22 ) + [iω 12 -(γ + Γ/2)] ρ 12 (2.27b) ρ21 = i Ω 12 cos(ωt) (ρ 22 -ρ 11 ) -[iω 12 + (γ + Γ/2)] ρ 21 (2.27c) ρ22 = i Ω 12 cos(ωt) (ρ 21 -ρ 12 ) -Γ ρ 22 (2.27d)
The set of equations given in (2.27) is known as the Optical Bloch Equations [17].

Multi-level system

In reality, the incident electromagnetic field consists in a laser pulse which is characterized by a certain spread in frequency and the physical system may have many energy levels close enough so that the action of such a laser pulse can drive the system to a coherent superposition of all those states (see Fig. 2.2). Thus, the model should be extended to multilevel systems. However, the physical insight gained from two-level system simulations can still be used to analyse multi-level results. The evolution of a multi-level system can be calculated using Eq.(2.26), provided that the total Hamiltonian H is given by

H = H 0 + H I (2.28)
where

H 0 =  ω  | | (2.29)
and

H I (r, t) = ,ı Ω ı (r, t) | ı| . (2.30)
Note that, compared to Eq.(2.17), we have used here a generalized definition of an instantaneous Rabi frequency associated with the value of the electric field at position r and time t

Ω ı (r, t) = -| µ |ı • E(r, t) (2.

31)

The density matrix ρ for multilevel systems can be defined from its elements

ρ ı = |ρ|ı . (2.32)
It is important to note that for this study, the electromagnetic field that will be used to turn on the dynamics of the system will be weak enough so that stimulated emission and multiple excitations can be neglected, such that

ρ 00  ρ  (2.33)
Following the same weak field assumption, we will also neglect coherences between two excited states. This assumption reduces the number of equations to be solved for analyzing the dynamics of the system. Following now the same development as the one detailed for two-level systems, one can easily show that the N -level density matrix elements obey the following optical Bloch equations

ρ11 = >1 i Ω 1 (r, t) [ρ 1 -ρ 1 ] + Γ ρ  (2.34a) ρ1 = i Ω 1 (r, t) [ρ 11 -ρ  ] + [iω 1 -(γ + Γ/2)] ρ 1 (2.34b) ρ1 = i Ω 1 (r, t) [ρ  -ρ 11 ] -[iω 1 + (γ + Γ/2)] ρ 1 (2.34c) ρ = i Ω 1 (r, t) [ρ 1 -ρ 1 ] -Γ ρ  (2.34d)
A detailed report on the density-matrix approach to the dynamics of multilevel atoms in electromagnetic fields can be found for instance in [18]. These equations describe the evolution of the system. They allow the evaluation of the density matrix elements provided the decay rates of the system and the field acting on the system are known. This system of equations is necessarily coupled to Maxwell's equation, describing the propagation of the electromagnetic field in the medium.

Connection with the electromagnetic field

The action of the electromagnetic field produces a macroscopic polarization in the medium which can be calculated, in the linear regime, via

P(r, t) = N μ , (2.35)
where N is the number of quantum emitters per unit volume and μ is the expectation value of the transition dipole calculated from the density matrix at position r and time t using the following definition

μ = Tr[ρ(r, t) µ] . (2.36)
Solving Eqs.(2.34) requires to know the electromagnetic field since the Rabi frequencies demand the knowledge of the electric field as a function of t and r. In turn, solving Maxwell's equations (2.3) at any position and time requires to know the density matrix at the same position and time in order to calculate the macroscopic polarization. Thus, the equations (2.34) and (2.3) have to be solved self-consistently. These set of equations are called Maxwell-Bloch equations.

Once the electromagnetic field is calculated, the energy flow into and out of the system can be estimated by calculating the Poynting vector from which the reflection and transmission spectra are obtained. The difference between the incident and the transmitted plus reflected energy flows gives finally access to the energy loss in the medium. This procedure can therefore be used to calculate the extinction spectrum associated with the system.

Non-Hermitian Wave Packet Propagation Technique

An alternative approach based on non-Hermitian Hamiltonians was suggested recently [19] for describing the dynamics of the system under the influence a of weak field. In principle, the energies of the ground state |1 and of the excited state |2 are real, and correspond to the frequencies ω 1 and ω 2 . The idea presented in [19] is to introduce some arbitrary imaginary parts to these energies so that by fixing them in comparison with the optical Bloch equations describing the two-level system given by Eq. ( 2.27), one can evaluate the dissipative dynamics in a simpler way.

Let us therefore assume that the ground state energy includes an imaginary part (+ γ 1 /2) and that the excited state energy includes (γ 2 /2) in its imaginary part. Writing down the TDSE including the new terms yields in atomic units

i Ċ1 = (ω 1 + iγ 1 /2) C 1 (t) + Ω 12 (t) C 2 (t) (2.38a) i Ċ2 = Ω 12 (t) C 1 (t) + (ω 2 -iγ 2 /2) C 2 (t) (2.38b)
It is also possible to derive a density matrix form of evolution from Eq. (2.38) where

C * ı C  = ρ N H ı .
The superscript N H stands for the non-Hermitian form of the Hamiltonian.

For the populations, the evolution of the density matrix elements takes the form

ρNH 11 = i Ω 12 (t) (ρ N H 12 -ρ N H 21 ) + γ 1 ρ N H

11

(2.39a)

ρNH 22 = i Ω 12 (t) (ρ N H 21 -ρ N H 12 ) -γ 2 ρ N H 22 (2.39b)
where ρNH 11 + ρNH 22 = 0 due to the conservation of the norm. It therefore comes that, necessarily,

γ 1 ρ N H 11 (t) = γ 2 ρ N H 22 (t) (2.40)
and it appears that the ratio (γ 1 /γ 2 ) has therefore to be time-dependent.

We now write the complete set of evolution equations for the density matrix elements of the system:

ρNH 11 = i Ω 12 (t) (ρ N H 12 -ρ N H 21 ) + γ 2 ρ N H

22

(2.41a)

ρNH 12 = i Ω 12 (t) (ρ N H 11 -ρ N H 22 ) + [i ω 12 -(γ 2 -γ 1 )/2] ρ N H 12 (2.41b) ρNH 21 = i Ω 12 (t) (ρ N H 22 -ρ N H 11 ) -[i ω 12 + (γ 2 -γ 1 )/2] ρ N H

21

(2.41c)

ρNH 22 = i Ω 12 (t) (ρ N H 21 -ρ N H 12 ) -γ 2 ρ N H 22 (2.41d)
An immediate comparison with the optical Bloch equations (2.34) gives two options for choosing the additional imaginary parts γ 1 and γ 2 introduced in the energies:

-Choosing γ 2 = Γ allows to describe correctly the population dynamics at the cost of a bad description of the coherences ρ 12 (t) and ρ 21 (t).

-Choosing γ 2 -γ 1 = 2γ + Γ allows to describe correctly the coherence dynamics while loosing the accuracy in the description of the populations ρ 11 (t) and ρ 22 (t).

At this point it is interesting to note that, unlike in the optical Bloch equations, the decay terms which appear in Eq.(2.41) are due to the imaginary terms introduced in the state energies. This is in contrast with the optical Bloch equations where they are defined by the super-operator Λ (2.26) which contains all decay effects in the system.

Note also that the two options mentioned previously are not compatible if one wants to maintain the conservation of the norm, as defined in Eq. (2.40). A choice has therefore to be made between these two options. Now assuming that the system is initially in its ground state, and looking at Eq.( 2.41) we can see that the evolution of the populations [Eqs. (2.41a) and (2.41d)] involves terms of second order in Ω 12 while the evolution of the coherences [Eqs. (2.41b) and (2.41c)] involves terms of first order in Ω 12 . We hence adopt the second option

γ 2 -γ 1 = 2γ + Γ. (2.42)
In the linear regime, this choice should be able to describe the dynamics of the system adequately.

From Eqs.(2.40) and (2.42), we obtain

γ 1 (t) = (2γ + Γ) |C 2 (t)| 2 |C 1 (t)| 2 -|C 2 (t)| 2 (2.43a) γ 2 (t) = (2γ + Γ) |C 1 (t)| 2 |C 1 (t)| 2 -|C 2 (t)| 2 (2.43b)
Another difference between the decay terms in the optical Bloch equations and in the non-Hermitian model is that in the later they are time dependent. The time dependence is coming into the picture described by Eq.(2.43) via the populations

|C 1 (t)| 2 and |C 2 (t)| 2 .
As mentioned before, the reason for the intrinsic wrong description of the populations in the chosen model is due to the fact that they are formulated as quantities that help to adapt the equations for describing correctly the coherences, which are essential for analyzing the Chapter 2. Atom-Field Interactions in Dense Media dynamics of the system. This non-Hermitian model is supposed to hold if

|C 1 (t)| 2 |C 2 (t)| 2 (2.44)
at any time t.

With this choice of γ 1 (t) and γ 2 (t), the evolution of the coefficients C 1 (t) and C 2 (t) can finally be written in atomic units as

i Ċ1 = ω 1 + i (γ + Γ/2) |C 2 (t)| 2 |C 1 (t)| 2 -|C 2 (t)| 2 C 1 (t) + Ω 12 (t) C 2 (t) (2.45a) i Ċ2 = Ω 12 (t) C 1 (t) + ω 2 -i (γ + Γ/2) |C 1 (t)| 2 |C 1 (t)| 2 -|C 2 (t)| 2 C 2 (t) (2.45b)
It is possible to generalize the non-Hermitian model to multi-level systems by following the same path, i.e, by adding imaginary parts: (+ γ 1 /2) to the ground state and (γ  /2) to the  th excited state, where  = 2, 3, ..., N , and solving the equations in comparison with the optical Bloch equations given in Eq. (2.34). One arrives at the following expressions for the time-dependent decay rates

γ 1 (t) = (2γ + Γ) k 2 |C k (t)| 2 |C 1 (t)| 2 -k 2 |C k (t)| 2 (2.46) γ  (t) = (2γ + Γ) |C 1 (t)| 2 |C 1 (t)| 2 -k 2 |C k (t)| 2 (2.47)
and the evolution of the coherences can then be written in atomic units as

i Ċ1 = ω 1 + iγ 1 (t)/2 C 1 (t) +  2 Ω 1 (t) C  (t) (2.48a) i Ċ = Ω 1 (t) C 1 (t) + ω  -iγ  (t)/2 C  (t) (2.

48b)

The advantage of this model is clear while implementing the numerical scheme. Optical Bloch equations demand the evaluation of a full density matrix which is in practice a set of first order differential equations. The density matrix elements have to be calculated for each couple of states, i.e, three equations for each transition, which is time consuming if the quantum emitter has many levels. On the contrary, the non-Hermitian model requires only two equations (corresponding to the coefficients C 1 (t) and C  (t)) to be solved for each transition. There is therefore a clear practical advantage for the non-Hermitian wave packet propagation if there are many allowed transitions in the system. It is however limited to weak fields in order to fulfill the condition (2.44).

System

Before going into the details of the numerical evaluation of the coupled Maxwell-Bloch equations described above, I will describe here the type of system and its geometry.

The system is a nanometric thick infinite layer of quantum emitters having two flat surfaces (see Fig. 2.3). The quantum emitter density N and the thickness of the layer are fixed so that the effects arisen in the system can only be due to the fundamental dynamics of the quantum emitters themselves. These quantum emitters could be atoms, molecules or nano-particles. The potential transitions involved are, for example, transitions between electronic, vibrational or rotational levels in a molecule. Another candidate could be electronic transitions in atoms, including if necessary the fine and hyperfine structures of the atom, such as in the D 1 and D 2 lines of rubidium for instance. To fulfill the condition (2.33) the molecules/atoms have to be initially in their ground state.

The prepared system is exposed to a ŷ-polarized electric field acting normal to the nano-layer (see Fig. 2.3). The electric field drives the quantum emitters to excited states. This field driven dynamics produces a macroscopic polarization in the layer which is the cause of back-actions to the field which can be studied by analyzing the reflected and transmitted fields. The associated transmission and reflection spectra can be used to extract some information regarding the dynamics turned on by the external field in the nano-layer.

In practice, the electromagnetic field is applied to the system through a laser source. By focusing the laser radiation to a spot whose size is of a few µm, a layer of quantum emitters with a transverse surface of about 1 cm 2 , and with a finite nano-thickness in the direction of the propagation of the field, can be easily approximated to an infinite layer of quantum emitters. Such a system reduces the dimension of the problem to the one of a 1D system. Indeed, in this case Maxwell's equations (2.3a) and (2.3b) reduce to

µ 0 ∂ t H z (x, t) = -∂ x E y (x, t) (2.49a) ε 0 ∂ t E y (x, t) = -∂ x H z (x, t) -∂ t P y (x, t) (2.

49b)

Note also that if, numerically, Maxwell's equations are treated in 1D, an additional factor 1/ √ 3 should be introduced in the polarization P y (x, t) (or equivalently in the transition dipole) in order to take into account the averaging performed over the other directions x and ẑ [4].

Numerical implementation

The system is described quantum mechanically either via optical Bloch equations or via non-Hermitian wave packet propagations. These descriptions of the system's dynamics are nothing but some sets of first order differential equations which can be solved using the fourth order Runge-Kutta method. This method is discussed in Appendix B. Solving Maxwell's equations (2.3a) and (2.3b) in 3D or (2.49a) and (2.49b) in 1D can be done through a discretization of space and time. This can be done by using a finite difference time domain (FDTD) method [20] within which both the electric and magnetic fields, E and H, are propagated in time and space on a spatial and temporal grid. This method ensures the explicit description of the electromagnetic field and it is numerically very stable. The boundary conditions are automatically taken into account through a staggered grid which was developed initially by K. S. Yee [21,22] and later named as the Yee cell. Maxwell's equations are discretized on the Yee cell using the central finite difference scheme.

n x -3/2 n x -1 n x -1/2 n x n x + 1/2 n x + 1
Spatial Grid ( in units of ∆x)

n t -1/2 n t n t + 1/2
Temporal Grid ( in units of ∆t) ). The grid in red is for the electric field and the one in black is for the magnetic field. The grids for the electric and magnetic fields are shifted by ∆x/2 in space and ∆t/2 in time.

E y H z E y
The electric field is calculated at spatial points (n x -1)∆x, n x ∆x, (n x +1)∆x, etc... for the times (n t -3/2)∆t, (n t -1/2)∆t, (n t +1/2)∆t, etc... The magnetic field is calculated at spatial points (n x -3/2)∆x, (n x -1/2)∆x, (n x + 1/2)∆x, etc... for the times (n t -1)∆t, n t ∆t, (n t + 1)∆t, etc...

Since this discretization necessarily induces some error in the evaluation of the derivatives, the spatial and temporal grid steps have to be chosen carefully. Let ∆t be the time taken by the electromagnetic wave to travel by ∆x in the vacuum, then ∆t = ∆x/c where c is the speed of light in vacuum. If ∆x is the spatial grid step then the temporal grid step ∆t should be less than ∆t , i.e, ∆t ≤ ∆x c .

(2.50)

Eq.(2.50) is known as the Courant-Friedrichs-Lewy condition (CFL), which appeared in 1928 in the seminal paper of R. Courant, K. Friedrichs and H. Lewy explaining the conditional stability of explicit time-dependent finite difference schemes [23,20]. In our simulations, we use the temporal step ∆t = ∆t /2, i.e. for a spatial grid step ∆x = 1 nm the temporal grid step is ∆t = 1.67 as.

Yee's scheme consists in considering the electric and magnetic components shifted spatially by half the cell step (∆x/2) and by half the temporal step (∆t/2) when considering a central difference approximation of the derivatives. Figure 2.4 shows the Yee cell in one dimension used for discretizing the electromagnetic field equations in space and time.

The algorithm can be summarized in three steps:

1. Replace all the derivatives in Eq.(2.49) with finite differences. Discretize space and time so that the electric and magnetic fields are staggered in both dimensions, i.e (b) ,

E nt+1/2,nx+1 y -E nt+1/2,nx y ∆x = -µ 0 H nt+1,nx+1/2 z -H nt,nx+1/2 z ∆t (2.51a) H nt,nx+1/2 z -H nt,nx-1/2 z ∆x = -ε 0 E nt+1/2,nx y -E nt-1/2,nx y ∆t -∂ t P nt-1/2,nx y (2.51b)
2. Solve (2.51) to obtain "updated" electromagnetic fields from past (known) fields using

E nt+1/2,nx y = E nt-1/2,nx y + ∆t ε 0 H nt,nx-1/2 z -H nt,nx+1/2 z ∆x -∂ t P nt-1/2,nx y (2.52a) H nt+1,nx+1/2 z = H nt,nx+1/2 z + ∆t µ 0 ∆x E nt+1/2,nx y -E nt+1/2,nx+1 y (2.52b)
3. Update the electric and magnetic field for each time step by using Eq. (2.52) until the fields have been obtained over the desired duration.

Figure 2.5 shows a schematic view of the simulation region for better understanding. The nano-layer of length is placed at the center of the spatial grid (in x), within which Maxwell's equations are discretized along with a temporal grid.

In order to avoid unphysical numerical reflections of outgoing electromagnetic fields back to the simulation region, the boundaries are covered with artificial non-physical absorbing layers which can efficiently absorb any incident electromagnetic field. It can be done by using the most powerful perfectly matched layers (PML) [24,25]. The numerical approach used for this purpose is the most economical variant of this method known as convolutional perfectly matched layers (CPML) [26].

(b) Note that ∂ t P y can be calculated directly from the quantum evolution equations of the system. : Schematic view of the numerical implementation: a nano-layer of thickness , exposed normal to a ŷ-polarized electric field, is placed in a spatio-temporal grid. The incident, reflected and transmitted fields are calculated by solving Maxwell's equations using the finite-difference time domain method. Polarization inside the system is calculated using one of the models described in the previous sections. Perfectly matched layers are used at the boundaries of the simulation region in order to avoid unphysical reflections of the electromagnetic fields back to the simulation region.

The ŷ-polarized electromagnetic field is acting normal to the nano-layer. A part of the incident field is reflected from the interface of the nano-layer. The rest of the field is transmitted through the layer. A part of the field is however lost in the layer due to the dissipative forces in the system. The calculation scheme of the cross-coupled equations, either Maxwell-Bloch equations or Maxwell-non-Hermitian wave packet equations, is given below [27] (see Fig. The magnetic field H is calculated in accordance with Eq.(2.49a) and then by using Eq.(2.49b), the electric field E is updated. The knowledge of the macroscopic polarization P calculated via the density matrix elements ρ ı (or via the wave packet coefficients C  ) from the previous time step is utilized in this step of the calculation. The density matrix elements (or the wave packet coefficients) are then updated at all spatial grid points using the electric field calculated previously. Finally, from the calculated electromagnetic field components and the density matrix elements (or the wave packet coefficients), the macroscopic polarization P is updated at each grid point.

Comparison of the two quantum models

The self consistent model described above was implemented for both Maxwell-Bloch equations and Maxwell-non-Hermitian wave packet equations as described in [19].

As a preliminary step we compare the two models by calculating coherences and populations as a function of time. The calculation is performed for a number densities N = 2.5 × 10 21 cm -3 in a nano-layer of thickness = 400 nm consisting of two level atoms having a radiationless decay rate Γ = 0.5 THz and a pure dephasing rate γ = 10 THz for a transition energy of 2 eV with a transition dipole of 2 D. The system is excited with a weak delta kick (in practice a 180 as Gaussian pulse of central frequency tuned near the transition frequency). Fig. 2.7 shows the results obtained with these two models. The evolution of the population of the excited state is shown in panel (a). The solid blue line shows the excited state population calculated via a non-Hermitian wave packet propagation and the red dashed line shows the same calculated by using optical Bloch equations. As expected, for the non-Hermitian model, the population shows a clear mismatch with the optical Bloch equations, which are, by definition, describing them correctly. This is because of the arbitrary choice made on γ  (t) for  = 1, 2 [Eq. (2.42)]. The evolution of the ground state gain factor γ 1 (t) for this system is shown in Fig. 2.8. One can easily notice the similar time-dependence of γ 1 (t) [Fig. 2.8] and of the excited state population |C 2 (t)| 2 [Fig. 2.7, panel (a), blue line]. This is not surprising in view of the equation (2.43a). For a system prepared according to Eq.(2.44) and probed by a weak field, the ground state gain factor in the non-Hermitian wave packet model γ 1 (t) is nothing but the excited state population scaled by (2γ + Γ). Meanwhile, the wrong description of the populations ensures the accuracy in the description of the coherences as it can be noticed from Fig. 2.7, panel (b). Thus, for the analysis of coherence dynamics, especially for the calculation of the macroscopic polarization P y (x, t) that we feed into Eq.(2.52), one can choose any of these two models in weak fields.

For a multilevel system there is however a clear advantage for choosing the non-Hermitian wave packet model in terms of computation time. Fig. 2.9 shows a comparison of the CPU time used by a non-Hermitian wave packet simulation compared to solving the optical Bloch equations, as a function of the number of quantum states involved.

The graph shows the CPU time taken by both models for a 1D spatial grid of total size 2.56 µm with a spatial step ∆x = 1 nm. The coupled equations were propagated for a nano-layer of width = 400 nm for a total propagation time of 1.67 ps with a time step of 1.67 as. The excitation pulse is a delta kick similar to the one used for Figs. 2.7 and 2.8. These estimations were done on a Intel Xeon E5-1650 processor. The blue line shows the time taken by the non-Hermitian wave packet propagation and the red line is the time taken to solve the optical Bloch equations.

There are some situations in which the multilevel system is characterized by low damping rates. In this case, the equations have to be propagated for longer times, and it is advisable to choose the non-Hermitian wave packet propagation scheme instead of the "exact" optical Bloch equation scheme.

Conclusion

The basic theory describing the dynamics of quantum emitters in the presence of weak electromagnetic fields was discussed in this chapter in order to introduce the basic concepts and the fundamental assumptions made for the present research.

The coupled equations known as Maxwell-Bloch equations were introduced along with the numerical implementation of these dynamical equations in a spatial and temporal grid.

The system studied was introduced. It is a nano-layer composed of quantum emitters (atoms, molecules, nano particles, ...) driven by an electromagnetic field acting normal to the layer surface. The finite difference time domain (FDTD) method was described with relevant details for such geometry.

An alternative model for the description of the coherences in the system, the non-Hermitian wave packet model, was also discussed. This approximate model can reduce the computation time when dealing with multilevel systems. Preliminary calculations comparing optical Bloch equations with the non-Hermitian approach were done.

The theory discussed here will be used in the forthcoming chapters for analyzing both two-state model systems and more complicated realistic quantum multilevel systems.

Introduction

The simplest case for studying the interaction of light with matter as described in the previous chapter is a two-level system where only one of the many allowed transitions (of frequency ω 12 ) of a quantum emitter is in resonance with the applied field. The dynamics of such system can be studied using Eq.(2.34). Reflected and transmitted electromagnetic energies are calculated from Maxwell's equations and by normalizing them with respect to the incident energy, the reflection and transmission probabilities can be calculated.

Since the system under consideration is a collection of two-level emitters, there are many other effects that should be considered for understanding the spectra calculated from the numerical simulations which includes broadening mechanisms and other density effects. Temperature dependent broadening mechanisms like Doppler and pressure (collisional) broadening can be reduced to a certain extend in a dense system.

Spectral Broadening, Lineshapes and Shifts

Let us assume that there are N particles (two-level quantum emitters) per unit volume. In weak fields, the response of the system can be described in terms of the electric susceptibility χ e (ω) which is a complex quantity that does not depend on the geometry of the system.

For the present discussion, we are only concerned with the time-dependence of a real monochromatic optical wave of linear polarization that we write as

E(t) = E 0 cos(ωt) ŷ = E 0 2 e iωt + e -iωt ŷ (3.1)
Being in the linear regime, the polarization can be obtained from Eq.(2.2) as

P = 1 2 ε 0 E 0 χ e (ω) e iωt + χ e (-ω) e -iωt ŷ (3.2)
since the field contains two conjugate frequencies ±ω.

The polarization can also be expressed using Eq.(2.35) as written below

P = N μ . ( 3.3)
Knowing from Eq.(2.16) that

μ = |1 2| + |2 1| µ 12 ŷ (3.4)
where µ 12 denotes the dipole matrix element 2| µ ŷ |1 , we finally obtain

P = N ρ 12 + ρ 21 µ 12 ŷ . ( 3.5)
Let us now introduce the slowly-varying density matrix elements

ρ12 = ρ 12 (t) e iωt (3.6a) ρ21 = ρ 21 (t) e -iωt (3.6b)
in Eq.(3.5) to obtain

P = N ρ21 e iωt + ρ12 e -iωt µ 12 ŷ . (3.7)
The comparison of Eqs.(

yields

χ e (ω) = 2µ 12 ε 0 E 0 N ρ21 = 2µ 12 ε 0 E 0 N ρ 21 (t) e -iωt . (3.8)
It is now possible to solve Eq.(2.34c) (with  = 2) analytically if we assume that ρ21 varies slowly with time (i.e. ρ21 ≈ 0) and that we are in the linear regime (i.e. ρ 11 -ρ 22 ≈ 1 if we start initially from the ground state). In this case, within the RWA [see Eq.(2.19)] we obtain

ρ21 = ρ 21 (t) e -iωt ≈ -iΩ 12 /2 i(ω 12 -ω) + (γ + Γ/2) . ( 3.9)
From these equations and from the definition (2.17) of the Rabi frequency it comes that

χ e (ω) ≈ N µ 2 12 ε 0 1 (ω 12 -ω) -i (γ + Γ/2) . ( 3.10)
The expression (3.10) gives the RWA approximation of the susceptibility, thus neglecting the (very small) counter-rotating contribution. Note also that in this expression, as mentioned at the end of Section 2.5, the factor 1/3 that one may find in the literature (due to an average over the orientation of the dipole) is here already included in µ 2 12 .

Separating the real from the imaginary part, we have Since in weak fields the transition probability between states |1 and |2 is proportional to |χ e (ω)| 2 , it appears from Eq.(3.10) that this transition is characterized by a Lorentzian shape. It can be (and in practice will be) modified by other physical effects like inhomogeneous broadening [1] (which assumes that different atoms see different electromagnetic fields), homogeneous broadening [2], or local field effects [3] for instance.

χ e (ω) ≈ N µ 2 12 ε 0 ω 12 -ω (ω 12 -ω) 2 + (γ + Γ/2) 2 + i γ + Γ/2 (ω 12 -ω) 2 + (γ + Γ/2) 2 . (3.
In general the particles are distributed within the nano-layer with different velocities. Assuming a Maxwell-Boltzmann distribution of the particles of mass m [4], the number density N (v ) associated with the speed v along the direction k of the propagation of the electromagnetic field can be written as

N (v ) = N m 2πk B T 1/2 exp - v 2 v 2 p , (3.12)
where N is the total number density for the uniform distribution of particles and

v p = 2k B T m (3.13)
is the most probable velocity at temperature T . This velocity distribution which is due to the distribution of kinetic energies possessed by the particles has to be considered when dealing with a realistic sample. At room temperature, say for noble gases around 25 • C, the most probable velocity varies between 250 to 1500 m/s depending on the mass of the atom. This thermal motion induces a shift in the resonance frequency which is nothing but the Doppler shift which can be written as

∆ Doppler = -k v (3.14)
where k is the propagation constant of the applied field. Eq.(3.12) being the particle velocity distribution in the direction of the propagation of the applied field, this thermal distribution will induce an inhomogeneous broadening in the transition spectrum that should, in principle, be taken into account. The spectral shape induced by the Gaussian distribution of particles will also be Gaussian.

Combining the natural broadening with the inhomogeneous Doppler broadening, the susceptibility will become a function of velocity and will have a temperature dependence

χ e (ω, v ) = N µ 2 12 ε 0 1 (ω 12 -ω -k v ) -i (γ + Γ/2) m 2πk B T 1/2 exp - v 2 v 2 p (3.15)
to be compared with Eq. (3.10).

Averaging over all velocities in the nano-layer

χ e (ω) = +∞ -∞ χ e (ω, v ) dv (3.16)
is equivalent to the calculation of the convolution of a Lorentzian shape with a Gaussian shape. It is well known that the lineshape offered by Eq.(3.16) is a Voigt profile [5]. The Gaussian lineshape associated with the thermal (Doppler) broadening contributes more around the resonance frequency and dies of quickly in the wings in comparison with the Lorentzian shape due to the natural decay. Thus the Voigt profile has a Gaussian core and a Lorentzian tail.

There are many other broadening mechanisms like pressure broadening or collisional broadening in the system which can be important when one deals with systems with large number of particles. The source of the collisional broadening is the close encounters between the quantum emitters carrying an electric dipole due to the polarization. Calculation of precise expressions for these type of broadening mechanisms is complicated since it involves many particles in different quantum states. But the collisional broadening is often treated as a Lorentzian-type broadening [2,6] which makes the calculation of the susceptibility not so hard in comparison with the inhomogeneous effects. Even though the timescales of collisional broadening and natural broadening are different it is proven that they usually have comparable values [2,6].

The dependence of the pressure broadening can be estimated from classical considerations. Let τ be the relaxation time for a mean free path l = 1/(N σ) where σ is the collision cross section, then

τ = l v p , ∝ 1 N √ T . ( 3.17)
Thus the shift in the frequency is proportional to the number density N and to the square root of the temperature. A detailed review on the broadening mechanisms and associated frequency shifts in two level systems can be found in [7].

From this point we will assume that we deal with relatively cold atoms or molecules (temperatures typically below 1 K) and the constituent particles are assumed to be frozen. The Doppler broadening will therefore be neglected in order to investigate the collective effects due to the local fields. This also helps for keeping the profile simple and easily tractable. On the contrary, the contributions from the collisional broadening are included in our model using the Lindblad formalism of the operator Λ (see chapter 2, Eqs.(2.26) and (2.27)) [8,9].

Lorentz Model for Two Level Systems

An alternative classical way of deriving an expression for the susceptibility is by modeling the system as a collection of a large number of oscillating dipoles driven by an external field. This is the Lorentz oscillator model, and it takes all relevant physical processes into account classically. This approach is usually valid in the linear regime where the susceptibility does not dependent on the applied field.

The basic assumption of the Lorentz model is that an incident light wave propagating through the medium of quantum emitters excites the system. The quantum emitters are considered as identical damped harmonic oscillators. The evolution of the induced oscillating dipoles is due to a displacement of the positive and negative charge densities in Thus the plasma frequency given by Eq.(3.23), which is due to the interaction with the neighboring quantum emitters, can be written as [3]

12 = 2N ω 12 ε 0 µ 12 , (3.30)
and Eq.(3.25) becomes

ω 2 12 = ω 2 12 - 2N ω 12 µ 2 12 3 ε 0 . ( 3.31)
Let us now define ∆ 12 as the dipole shift induced in the resonance frequency as

ω 12 = ω 12 -∆ 12 , (3.32)
it appears from Eq.(3.31) that

∆ 12 = ω 12   1 -1 - 2N µ 2 12 3 ω 12 ε 0   (3.33)
If ∆ 12 ω 12 , the first order approximation on Eq.(3.33) will reduce ∆ 12 to

∆ 12 ≈ N µ 2 12 3 ε 0 , ( 3.34)
which is nothing but the well-known Lorentz-Lorenz (LL) shift [10,11] introduced in Eq.(1.1). This LL shift and its contribution towards collective effects will be discussed later in this chapter and the forthcoming chapter.

In the Lorentz model, the frequency 12 is known as the plasma frequency that explains the oscillation of charge centers in conducting media such as metals and plasmas. In the present case, it represents the natural collective oscillation frequency of the "sea" of dipoles in the system, not of individual dipoles. This collective oscillation of the coupled dipoles matters for the response of the system when the dipole-dipole interaction is not negligible. Indeed, these collective or cooperative behavior of the dipoles (which arise from the fact that the electromagnetic field experienced by an arbitrary quantum emitter in the system is modified by its neighboring dipoles) can modify the response of the system towards the electromagnetic field drastically [12][13][14].

Detuning and Scaling

There are some quantities that help to generalize the equations and variables to an adimensional form. The collective effects are represented by the plasma (or dipole coupling) frequency 12 . Eq.(3.30) can be rewritten in terms of the LL shift,

2 12 ≈ 6∆ 12 ω 12 . ( 3.35)
The Lorentz-Lorenz shift contains the essential part in the coupling frequency that describes how strongly the system is interacting. It will be shown in the forthcoming section that the quantity N µ 2 12 is the key parameter for seeing the collective effects in the nano-layer of the two-level systems under consideration. Hence a dimensionless quantity called interaction strength S 12 is defined as

S 12 = ∆ 12 γ + Γ/2 (3.36)
that helps to quantify the interaction regime. Note that the relation between the classical damping rate γ m and the decay and decoherence rates is γ + Γ/2 = γ m /2. For large values of S 12 the dipole coupling frequency dominates over the damping rate, and the constituent particles interact strongly.

In the forthcoming section a reduced detuning is also used to analyze the optical response of the dense two-level system. This reduced detuning is defined as

δ = ω -ω 12 γ + Γ/2 . ( 3.37)

Optical Response of a Two-Level System

When the dipoles are close enough, the radiation from a dipole couples it with its neighbors. This effect can be explored by calculating probabilities for reflection and transmission at different interaction regimes characterized by the value of S 12 [15].

The reflection probability R can be calculated from Eq.(2.14) by measuring the reflected energy, i.e,

R(ω) = |S(ω)| R |S 0 (ω)| (3.38)
where |S(ω)| R is the reflected energy.

The transmission probability T can be calculated in the same manner,

T (ω) = |S(ω)| T |S 0 (ω)| . ( 3.39)
where |S(ω)| T is the transmitted energy. Here the indexes R and T denote calculations of the energy flux performed in front of (x = x R < 0) and beyond (x = x T > ) the nano-layer.

A part of the energy given to the system will be lost (extinction) while propagating through the layer. It is not recovered in the reflection or in the transmission. The probability of extinction E can thus be calculated from the conservation of energy as

E(ω) = 1 -R(ω) + T (ω) . (3.40)
This loss of energy is related to the imaginary part of the propagation constant k inside the nano-layer. Reflection from a layer of weakly interacting identical dipoles is feeble.

Those are therefore multiplied by some constant factors given alongside with the same color codes for a better comparison.

Fig. 3.2 shows some reflection spectra from a 600 nm thick layer composed of identical two-level emitters for different values of S 12 [15]. The system is assumed to have a transition energy of 2 eV with a damping rate γ + Γ/2 = γ m /2 = 10.5 THz due to the decay and decoherence processes.

The incident electromagnetic fields are reflected selectively from the layer around the transition frequency [16,17]. The reflection increases with the increase in the interaction strength. For S 12 0.01, the system shows very weak reflection around the transition frequency and almost all energy given to the sample is transmitted or extinct. The response of the system towards the electromagnetic field is very close to the independent particle response characterized by a Lorentzian profile in the frequency domain. It is because in weak densities, the effect of the polarization of the medium barely modifies the local electric field experienced by the constituent dipoles. Hence they respond to the field as if they were independent. The reflection gets stronger as the interaction between the induced dipoles increases. It starts to show some additional features that will be explained later as S 12 approaches 1.

For S 12 10, the constituent particles feel very strongly the presence of neighboring dipoles that modify the system response. As the interaction strength increases, the layer becomes opaque, reflects most of the incident energy and absorbs the remaining. The energy absorbed by the system can be dissipated due to many physical phenomena represented by the decay constants γ and Γ. In the strong interaction regime (S 12 > 1), the presence of other dipoles becomes really significant. The system loses its transparency and behaves like a mirror around the transition frequency.

The transmission through the layer then falls to zero. Fig. 3.3 shows the extinction [panel (a)] and transmission [panel (b)] corresponding to weakly (S 12 = 10 -3 ) and strongly (S 12 = 12) interacting dipoles.

For S 12 = 12, the system behaves like a mirror by reflecting the major portion (≈ 80%) of the incident energy, the remaining portion being absorbed and dissipated. We note that the window of reflection within which the incident field is mainly reflected measures 3∆ 12 (see Fig. 3.2). This measure will be explained in the next section. The interaction between the induced dipoles enhances the reflection probability so that in this regime, the system (nano-layer of two-level quantum emitters) behaves like a metal. Note that a similar problem was studied analytically by R. J. Glauber and S. Prasad in 2000 [13,18] as an eigenvalue problem for a dense layer of quantum emitters.

Lorentz-Lorenz shift and reflection window

As demonstrated in Section 3.2.1, the LL correction on the local field experienced by an arbitrary quantum emitter induces a shift in the resonance frequency. We have shown using the Lorentz model that the frequency shift is ∆ 12 and that it originates from the presence of neighboring interacting dipoles. Indeed, as shown schematically Fig. 3.4, two neighboring oscillating dipoles interact through the electric field they create.

The induced LL-shift has been known for a long time and the initial theoretical prediction of Lorentz and Lorenz has been confirmed by several experiments [19][20][21][22]. It was for instance measured while studying the optical response of dense atomic vapors (for N of the order of 10 17 cm -3 ) in the vicinity of the fundamental atomic resonance frequency [19,23]. The dependence of the shift on the degree of excitation was also predicted [7] and experimental evidences showed the presence of such a dependence [24,25]. The validity of the LL shift has thus been shown by many different experimental groups. These measurement therefore validate the theories based on the LL correction introduced in the local field experienced by the quantum emitters [14,26].

In the calculation executed here for studying the linear optical response of dense two-level systems, the LL red-shift is obtained in the associated susceptibility χ e (ω), as shown in Eqs. (3.24), (3.32) and (3.34).

Eq.(3.24) is able to reproduce exactly the susceptibility extracted from the full quantum calculation [3]. Fig. 3.5 shows a comparison between the susceptibilities obtained from the quantum calculation and the semi-classical analytical model developed in this chapter [Eq. (3.24)] for the parameters used to obtain S 12 = 12 in Fig. 3.3, as a function of the reduced detuning δ. The blue solid line is the susceptibility obtained from the semi-classical model and the red dashed line shows the susceptibility obtained from the numerical integration of the optical Bloch equations. They mach in great accuracy, showing the following essential feature: the resonance frequency is red-shifted by the LL shift ∆ 12 . Analyzing the susceptibility, one can also estimate the width of the reflection window obtained for large densities (see Fig. 3.2). Apart from the energy dissipated in the layer, this system reflects a large part of the incident radiation from the surface which receives it, where the reflectivity

R = n(ω) -1 n(ω) + 1 2 (3.41)
is close to unity.

In this expression, n(ω) = Re[ 1 + χ e (ω)] is the real part of the refractive index of the medium. At maximum reflection R(ω) = 1, and the real part of the electric susceptibility goes to -1. Equating the real part of Eq.(3.24) with -1 gives the range of frequencies at which the reflection is maximum:

Re χ e (ω) = -1 ⇐⇒ ω 12 -∆ 12 < ω < ω 12 + 2∆ Even though the real part of the refractive index is zero, there is a non-vanishing imaginary part associated that explains the absorption and dissipation in the system and the associated zero transmission through the system. At R = 1, the nano-layer therefore behaves like a dissipative mirror. Zero transmission through the layer within the window of reflection means that the electromagnetic field cannot penetrate through the entire layer of thickness = 600 nm. The field can however propagate though the layer to a certain extend which is determined by the system parameters. The field inside the nano-layer is exponentially decaying and the distance from the interface to the point inside the layer at which the exponentially decaying evanescent field falls to 1/e is known as the skin depth, ζ(ω). By definition this skin depth is given by

ζ(ω) = 1 κ = c ω 1 Im 1 + χ e (ω) (3.43)
where κ is the imaginary part of the propagation constant. (3.32). The skin depth at ω 12 is very small (≈ 15 nm) due to the collective response of the two-level system. Fig. 3.7 shows the skin depth for the strongly interacting system (S 12 = 12) as a function of the reduced detuning δ around the shifted resonance frequency ω 12 . The skin depth reaches a minimum of about 15 nm at the shifted resonance frequency ω = ω 12 (corresponding to δ ≈ -12.4 in Fig. 3.7) i.e, the electromagnetic field of this frequency enters the layer and dies off very quickly, within a distance of 15 nm. In this situation the calculation of polarization, electromagnetic fields, populations and coherences should be done within this skin depth limit where the field is not yet zero.

Fabry-Pérot Modes

The sidebands seen in Fig. 3.6 in the tails of the reflection spectrum arise from the interference between multiply reflected and transmitted electromagnetic fields and can be interpreted as Fabry-Pérot modes [26,27]. Indeed, systems with more than one interface can produce interferences between multiply reflected and transmitted fields. The system studied here can be seen as a Fabry-Pérot étalon consisting of two parallel partially reflecting planes. Fig. 3.8 illustrates these multiple reflections and transmissions in the layer. The reflected electromagnetic fields are calculated at a point x R in region I (x R < 0) and the transmitted fields are calculated at x T in region III (x T > ). Inside the layer (region II), a part of the incident electromagnetic field is reflected back and forth, and is therefore associated with partial transmissions. These partially transmitted electromagnetic fields contribute to the reflected and transmitted fields at positions x R and x T . Due to the path differences between different contributions induced by the additional geometrical paths traveled inside the layer, the interference between them gives maxima and minima. These interference fringes appear as sidebands in the reflection spectrum (see Fig. Approaching the problem formally, the (complex) geometrical path difference between any successive fields can be written as

D = 2k (3.44)
where k = k r + iκ is the wave vector of the electromagnetic field inside the layer. Let t 1 be the transmissivity of the interface I → II and t 2 be the same for the interface II → III, then the primary field noted 1 st in Fig. 3.8 in the transmission region can be written as

E 1 st = t 1 t 2 e -κ e ikx (3.45)
where k is the wave vector in the vacuum. Similarly, the n th order electromagnetic field transmitted to region III due to the everlasting multiple reflections inside the layer can be written as

E n th = r 2 2 e -2κ e iD n-1 E 1 st (3.46)
where

r 2 2 = 1 -t 2 2 .
Eq. (3.46) suggests that the transmitted field contributions behave like a geometric progression. Thus the total transmitted field can be calculated by summing up all the contributions. Taking the ratio between the incident amplitude and the transmitted amplitude of the electric field, the collective transmissivity can then be extracted as

t(ω) = 4ne -κ (n + 1) 2 -(n -1) 2 e -i2nk
(3.47)

since we know the classical relations between the refractive index n and t 1 , t 2 and r 2

t 1 (ω) = 2 1 + n(ω) (3.48a) t 2 (ω) = 2n(ω) 1 + n(ω) (3.48b) r 2 (ω) = n(ω) -1 1 + n(ω) . ( 3

.48c)

Following the same path, the global coefficient of collective reflection can be derived as Fig. 3.9 shows a comparison between the calculated reflection spectrum obtained from the numerical solution of the optical Bloch equations and the analytical form of the reflection spectrum obtained via the semi-classical Lorentz model including the Fabry-Pérot formula treatment. The blue dotted line is the calculated spectrum and the red solid line is the analytical form of the reflection spectrum given by Eq. (3.50). The reflection window and the sidebands are represented quite accurately by the analytical model, giving a very accurate analytical tool for analyzing the calculated dynamics in more complex systems.

r(ω) = -2i(n 2 -1)e -ink sin(nk ) (n + 1) 2 -(n -1) 2 e -i2nk (

Conclusion and Outlook

The collective response of a nanometric layer composed of two-level quantum emitters was calculated by integrating the coupled Maxwell-Bloch equations. It was also analyzed using a semi-classical model based on the Lorentz oscillator which can describe the linear optical response of dielectrics.

A clear collective behavior of the system is obtained for large densities in which the predominant selective reflection associated with zero transmission is described by the collective oscillation frequency of the induced dipoles in the system.

The analytical model developed in this thesis is accurate enough to explain the collective effects observed in a dense layer of quantum emitters exposed to weak electromagnetic fields. Fig. 3.10 shows for instance the relative error in the susceptibility given by Eq. (3.24) in comparison with the electric susceptibility calculated with the optical Bloch equations. The blue solid curve is the relative error for weakly coupled dipoles and the red solid line is the same for strongly interacting dipoles. The error is maximum around the transition frequency, where the absorption of electromagnetic energy is more important and hence where the susceptibility varies abruptly. The shift in the curve for strongly interacting dipoles is due to the LL shift [15]. The relative errors in both cases are very small (of the order of a few 10 -4 ). Thus the analytical model developed here is accurate enough to explain the physics of both dilute and dense systems. All features including the broad reflection window for strongly interacting dipoles and the Frabry-Pérot sidebands were explained from classical optics considerations and with the help of the mean-field Lorentz-Lorenz correction. The model is therefore convincing for treating more complicated systems and this will be subject of the forthcoming chapters.

Introduction

The collective response of a nano-layer of two level quantum emitters was discussed in the previous chapter. It was shown that, for a dense layer of identical coupled dipoles, reflection becomes dominant within a window around the transition frequency along with a redshift [1][2][3].

For a nano-layer comprised of three-level quantum emitters one can expect a similar behavior with two predominant reflection windows associated with the two transitions with their own LL shifts whose amplitude will depend on the corresponding transition dipoles.

Depending on the energy splitting between the excited states, the system could behave differently. In this chapter, the optical response of such systems is discussed.

Multilevel Systems

The light sources and their associated bandwidths usually allow energetically more than one transition in a quantum emitter. Limiting our study to two-level systems only would therefore be very restrictive.

In a more realistic and general picture, one has to treat more than one excited states when dealing with light-matter interactions. In the previous chapter, a simple analytical model was developed for two level systems in order to analyze their optical response. This model will now be extended to the case of multilevel transitions within the same framework.

In a multilevel system, we consider (N -1) transition dipoles between the ground state |1 and the excited states | , with  = 2, ..., N . When the system is subjected to an external electromagnetic field, these dipoles will get excited, thus creating a macroscopic polarization in the system. Here the system is a bit complicated due to the presence of different types of macroscopic polarizations that we note P 1 (t) due to the transition dipoles µ 1 , where  = 2, ..., N .

As shown in the previous chapter, within the limits of linear optics, and thus for weak excitation where ρ 11

 2
ρ  , the equation of evolution of the macroscopic polarization associated with the |1 → | transition in the system can be written as

P1 + γ m Ṗ1 + ω 2 1 P 1 = ε 0 2 1 E local (t), (4.1)
where 1 is the coupling (or plasma) frequency associated with the  th transition, defined as in Eqs. (3.23) and (3.30). Note that we are still considering here an incident electromagnetic field of linear polarization.

The driving field E local (t) includes the fields radiated by all dipoles. Thus the constituent particles with many levels can be imagined as a collection of dipoles of different kinds immersed in a collective field of these dipoles along with the applied field.

Thus in a system with more than one type of dipoles, additional couplings arise between dipoles of different kinds which will be taken into account via the total polarization included in the local field

E local (t) = E(t) + 1 3ε 0 N =2 P 1 (t) . (4.2)
By taking into account the modification in the local field given by Eq.(4.2) and subsequently looking for the steady state solution, one can take the Fourier decomposition of Eq.(4.1)

W 1 (ω) P1 (ω) = ε 0 2 1 Ẽ(ω) + 2 1 3 N k=2 P1k (ω) , ( 4.3)
where

W 1 (ω) = ω 2 1 -ω 2 + iγ m ω (4.4) and 1 = 2N ω 1 ε 0 µ 1 . ( 4.5)
In Eq.( 4.3) we recognize the total polarization

P (ω) = N k=2 P1k (ω) (4.6)
and therefore W 1 (ω) +

W 1 (ω) P1 (ω) = ε 0 2 1 Ẽ(ω) +

1

3W 1 (ω) χ(ω) , (4.8) where χ(ω) is the total susceptibility of the system. Summing over all excited states we finally obtain the total susceptibility of the system as

χ(ω) = 3 X(ω) 3 -X(ω) , ( 4.9)
where

X(ω) = N k=2 2 1k
W 1k (ω) (4.10)

Eq.(4.9) gives the total susceptibility associated with a system having N levels in which all (N -1) excited states are coupled to the ground state only. There exists some systems which can potentially be candidates for such studies, like vibrational states of molecules [4] or electronic states of atomic species [2] for instance. The later will be discussed in detail in this chapter.

Dipole Induced Electromagnetic Transparency (DIET)

Eq.(4.9) is a general expression for an N -level system. Now looking at the specific case of a three-level system, this equation can be further simplified to get a clear picture. Indeed, the susceptibility for a three-level system can be simplified to

χ(ω) = 2 13 W 12 + 2 12 W 13 W 12 W 13 -1 3 ( 2 13 W 12 + 2 12 W 13 ) . ( 4

.11)

The resonance frequencies of the three-level system are shifted due to the coupling between dipoles of same kind as well as between dipoles of different kinds, i.e, due to the transitions |1 → |2 and |1 → |3 . Thus the resonance frequencies are expected to be shifted in a complex way. These resonance frequencies, the eigenfrequencies of the coupled oscillators at which the collective dipoles absorb the incident radiation, can be found when the real part of the denominator of Eq. (4.11) is zero. This is a simple second order equation in ω that one can easily solve for a specific system. It therefore appears that the contributions from different dipoles to the dynamics are not independent. The system responds to the field collectively with two eigenfrequencies, that we will denote by ω -and ω + .

When the splitting between the two excited states is large enough, the two dipoles will respond to the applied field independently, which means they produce their own reflection windows with certain width determined by the coupling terms. But the response is more complex for closely spaced transitions. Indeed, when the two dipoles are oscillating at similar resonance frequencies, the widths associated the large density could lead to an overlapping of these two transitions that potentially enables quantum interference effects between these two transitions. One such case is discussed below. The calculated reflection from a three-level system of transition frequencies ω 12 = 484 THz and ω 13 = 532 THz having transition dipoles µ 12 = µ 13 = 2 D is given in Fig. 4.1 as a function of the reduced detuning. The reduced detuning is defined with respect to the first excited state. These two transitions broaden due to dipole-dipole interactions and result in the overlapping of the two reflection windows.

Approximating the system to a collection of two types of dipoles, it is possible to calculate the expected reflection probability in the absence of coupling between dipoles of different kinds. For this we perform two separate calculations with a single dipole and we multiply the calculated reflection probabilities. The result, shown as a red solid line in Fig. 4.1, peaks around the resonance frequency ω + where both dipoles contribute towards the collective reflection. But in reality, these dipoles are coupled via dipole fields that change the response of the system dramatically with two well separated windows of reflection.

The blue curve shows the reflection probability from a system of coupled dipoles. Due to the couplings the reflection windows are broadened and overlap. In the region where they overlap, the reflection probability is maximum for an uncoupled system. The coupling between the dipoles leads however to an interference effect that produces a spectral hole, as seen Fig. 4.1. Thus, instead of a maximum in the reflection, the dipole coupling produces a minimum.

It can be seen as the cancellation of the induced dipoles in the medium [1]. When the dipoles in the system are aligned and out of phase, the resulting macroscopic polarization has a minimum value. Depending on the coupling, this minimum value could be zero when the system has no effective total dipole. The susceptibility associated with each dipole possesses independently a Lorentzian shape similar to the case discussed in the previous chapter. But here, in the presence of dipole-dipole couplings, the line shape is modified. Fig. 4.2 shows the susceptibility in the case of the three-level system discussed here.

The susceptibility shows a double peak structure that corresponds to the two transitions in the system. The first transition corresponds to the shifted resonance frequency ω -. It has a Lorentzian profile and the second peak shows a Fano profile characterizing the quantum interference between the two indistinguishable excitation pathways [5]. The minimum in the Fano profile corresponds to the cancellation of the two dipoles resulting in the minimum of the reflection window.

The frequency ω h at which the reflection has a minimum is

ω h = ω 2 12 2
13 + ω 2 (4.12)

It corresponds to the frequency at which the real part of the numerator of Eq. (4.11) is zero. If the second dipole is stronger than the first one (µ 12 < µ 13 ), ω h shifts towards the transition frequency ω 12 of the first excited state. In the opposite case (µ 12 > µ 13 ), ω h shifts towards ω 13 . Thus the position of the minimum reflection in the case of a three-level system is somewhere in between the transition frequencies of the system [1], depending on the coupling strengths. The most interesting information from the analysis done here is that a field of frequency ω h incident on such a system, will pass and will propagate though the layer without any collective reflective effects. The medium is lossy because of the propagation losses associated with the imaginary part of the refractive index.Thus the layer behaves like a lossy and transparent medium. Transparent in the sense that the real part of the refractive index is one. Thus for a layer of low decay effects, the field corresponding to the spectral hole ω h in the reflection spectrum can be observed in the transmitted signal. This transparency is induced in the system due to the destructive interference of the overlapping excitation pathways and hence named as Dipole Induced Electromagnetic Transparency [1,2]. Similar situation has previously predicted in the case of overlapping transitions in isolated atom or molecule [6,7]. Fig. 4.3 shows the transmission spectrum of the three-level system considered here. The peak corresponding to DIET is exactly at the position predicted by Eq. (4.12). This figure shows the transmission spectra for two different thicknesses. The blue curve is the transmission spectrum for a layer of thickness = 600 nm. Even though the field corresponding to the DIET frequency enters the system, feeling no reflection at the interface, a part of it is lost in the medium while propagating from one interface to the other via the decay mechanisms of the system characterized by γ m = 10.5 THz. This propagation loss depends on the thickness of the medium via the Beer-Lambert's law [8]. Thus a reduction in the layer thickness can increase the transmitted signal. The red curve in Fig. 4.3 is the transmission from a = 300 nm layer of the same three-level quantum emitters. Apart from the thickness, choosing a system with low damping rates, like with colder atoms for instance, can improve the signal remarkably. This is shown in Fig. 4.4 for a layer of thickness = 300 nm. The decay rate γ m is reduced to 5.5 THz in comparison with the previous calculation. The peak corresponding to DIET is increased compared to that of γ m = 10.5 THz shown in Fig. 4.3. The modification in the scale of reduced detuning is due to the difference in γ m .

Since the transmittance depends on the decay effects and hence on the imaginary part of the susceptibility of the system, the DIET signal could be absent in the transmitted signal of a system with high damping effects or large thicknesses even though the destructive interference effect is present. The presence of this destructive interference can be identified by looking at the reflected signal. The presence of a spectral hole in the reflection spectrum is a clear signature of the destructive interference between the radiation emitted by the coupled dipoles and hence it constitutes a signature of DIET [1].

Rubidium Atoms

These studies can be further extended to a realistic multilevel system having more than two excited states. In the strong dipole interaction (i.e. dense) regime, the Lorentz-Lorenz correction to the local field felt by the quantum emitters leads to redshifts in all transitions, as shown in the previous section. The most interesting case, discussed above, is when two nearby transitions are significantly overlapping. This leads to interferences between these two transitions which modify the optical response of the system. It leads to DIET, which is due to the interference between two indistinguishable excitation pathways [1]. A series of such narrow windows can thus surely be obtained in the case of a multilevel system.

The DIET effect in a realistic multilevel system can be studied in a nano-layer of dense rubidium atoms for instance. Rb is widely used both in thermal gases experiments [9] and in cold atomic gases experiments [10]. In addition, the Rb D line transitions are characterized by large S to P transition dipoles [11]. Having large transition dipoles is of special interest since the dipole-dipole interaction depends on the square of the transition dipoles. The same scales linearly with the number density. Thus Rb allows to have relatively low density with strong interactions. Among most common isotopes of Rb, 85 Rb is considered here. It is possible to choose either the D 1 or D 2 transitions in 85 Rb. Here the D 1 transitions are considered, i.e, the system under consideration is a nano-layer of dense 85 Rb atoms exposed to an electromagnetic field tuned near the 5 2 S 1/2 to 5 2 P 1/2 transitions (see Fig. 4.5) [11].

5 2 S 1/2 is the ground state of the system and 5 2 P 1/2 is the excited state. These levels are split into two sub-levels, as shown in Fig. 4.5, due to the hyperfine coupling [12] that effectively turn the system into a multilevel system with four dipole-allowed transitions. Transitions between the hyperfine levels of same electronic state are obviously forbidden. Thus the D 1 transitions of the 85 Rb isotope can be considered as a mixture of two three-level systems having same excited states but different ground states. From Eqs. (4.9) and (4.10), the susceptibility of such a system can be written as

χ(ω) =  2 1 /W 1 +  2 1  /W 1  1 -1 3  2 1 /W 1 +  2 1  /W 1  , ( 4.13)
where the indices 1 correspond to the transitions from the sub-level F = 2 of 5 2 S 1/2 and 1  is for the transitions from F = 3. Now the possible overlapping that can be achieved in such a system for observing the DIET effect are of two different kinds [2].

When the dipole interactions are strong enough the excited states can broaden, leading to overlapping resonances that will induce DIET. There are two such possible DIET signals in the system since transitions are possible to two different excited states from two different ground states. If the interaction between the dipoles extends further, it may exceed the ground state splitting, giving rise to a new DIET effect due to the overlap of the two different ground states.

The LL shifts for the D 1 transitions of Rb can be written as

∆ F F = ∆ 0 S 2 F F (4.14)
where ∆ 0 = N µ 2 SP /(3 ε 0 ) is the LL shift defined in the absence of hyperfine splitting and S F F is a measure of the relative strength of the transitions [11] that should be multiplied with the transition dipole µ SP in the absence of hyperfine splitting in order to obtain the transition dipoles to the hyperfine-splitted excited states. The interaction strength for the D 1 transitions is defined as

S SP = ∆ 0 γ m (4.15)
The electric susceptibility and the associated reflection spectrum from a = 600 nm thick 85 Rb atomic layer for S SP = 3.8 × 10 -3 are shown in Fig. 4.6. The sample is assumed to have negligible Doppler and collisional shifts. The reduced detuning δ is defined with respect to the D 1 transition frequency in the absence of hyperfine splitting [11]. presents the reflection spectrum and panel (c) shows the transmission. The expected overlapping of the two transitions originating from F = 3 (or from F = 2 also) result in two Fano profiles for the susceptibility. This explains the two sharp minima in the reflection spectrum. Concurrently, two narrow peaks appear in the transmission spectrum. Those are characteristic from a DIET effect [2]. We have here a first theoretical demonstration of the existence of DIET in a realistic multi-level system.

An additional increase of the dipole-dipole interaction can lead to the mixing of the F = 2 and F = 3 ground states into a (approximate) single excitation state which is coupled to the excited states. This leads to two broadened excitations from the two different ground states to an approximate single excited state which can overlap to result in a destructive interference. Fig. 4.8 shows this effect in very strongly interacting samples of 85 Rb atoms, with S SP = 200.

Panel (a) in Fig. 4.8 is the modulus of the susceptibility. Panel (b) is the reflection spectrum for a layer of thicknesses = 600 nm. Panel (c) is the corresponding transmission spectrum. The susceptibility of the system, in comparison with the previous case, is modified significantly. The system has an enhanced Lorentzian profile and one dominant and two faint Fano-type profiles. The reflection and transmission spectra of the system have a background where DIET is observed. This is because of the Fabry-Pérot etalon effect discussed previously. Due to the mismatching between the layer thickness = 600 nm and the transition wavelength λ 0 795 nm, long tails are created near the reflection window. This can be improved by using a thicker slab.

Conclusion

In this chapter, the semi-classical model of the Lorentz oscillator was extended to a multilevel system. This extended model can explain the calculated response of a realistic multi-level system towards weak electromagnetic fields.

A new collective effect named dipole induced electromagnetic transparency (DIET) was discussed. It is an induced electromagnetic transparency in a system of strongly interacting quantum emitters with two (or more) closely spaced excitations. The mechanism lies on the destructive interference between two overlapping resonances. The frequency ω h of the induced transmission window can be calculated using Eq. (4.12). An electromagnetic field of this frequency, incident on the system, will be transmitted through the layer with an associated extinction obeying Beer-Lambert's law.

DIET was discussed in the case of the D 1 transitions of a gas of 85 Rb atoms at low temperature. Different DIET effects were calculated due to the mixing of the different energy levels of the considered transition.

The signature of DIET is imprinted very clearly in the reflection spectrum, which shows a spectral hole at the transmission frequency ω h . Several applications can be suggested in which DIET could play an important role. Some of those will be discussed in the next chapter.

Introduction

An interference between two indistinguishable pathways leads to an induced electromagnetic transparency called dipole-induced electromagnetic transparency in a dense sample of multilevel particules. This effect was discussed in detail in the previous chapter. It provides some interesting applications. Two of the potential applications of DIET are pulse shaping and slow light. These applications are discussed in this chapter.

Pulse Shaping

Shaping a laser pulse in a desired way is an interesting application [1][2][3]. Pulse shaping or waveform synthesis is the process of changing the pulse shape (envelope) to a form that is optimized for a desired application. There are many ways for achieving this goal.

In some of the techniques the pulse is modified by using a beam shaper [2], spatial light modulators [4], deformable mirrors [5], etc... Pulse shaping techniques showed their strong impact as an essential experimental control tool in ultrafast spectroscopy, nonlinear fiber optics, and strong-field physics. There are many research fields where pulse shaping pulse shaping is used, like superconducting qubits [6], laser filament experiments [7], holographic and nonlinear pulse processing [4], chemical reactions, etc... The goal can be, for instance, to find an optimal pulse shape in order to minimize unwanted effects, for example to minimize the space-charge induced emittance growth in RF photoinjector cathodes [8], or to maximize the desired processes in experiments.

Input Output

The cancellation of the dipoles due to the destructive interference of the overlapping resonances seen in DIET can also be used for pulse shaping. If the layer of interacting three-level quantum emitters is characterized by a reflection window which is wider than the FWHM of the incident pulse ∆ω, the collective response of the quantum emitters results in a selective reflection and transmission of the incident pulse. Fig. 5.2 shows the reflected pulse from a layer of width = 3.5 µm of three-level systems with the same parameters as the one used in section 4.3, for an incident pulse of ∆ω = 40π THz. The blue curve is the reflected pulse shape from the layer which is compared with the incident pulse envelope shown by the green curve. The spectral hole is present at the position where DIET occurs. The associated transmission is shown in the inset of the figure.

The transmission is very small due to the high damping factor used in the calculations. All these plots are normalized with the peak intensity of the pulse. The reflected pulse has an important modification compared to the incident field. The pulse shape now looks like the superposition of two distinct pulses. If the position of the spectral hole that separates these two pulses can be controlled to a certain extend, it can be of great interest to eliminate a part of the frequencies from the pulse that could be of interest for many purposes.

Rewriting the equation for the frequency of the spectral hole

ω h = ω 2 12 2
13 + ω 2 it is clear that the position of the spectral hole depends on the coupling between the dipoles and on the energies of the excited states. Thus once the excited states and other parameters of the system are chosen, the control of the position of the hole over the reflection spectrum is limited. Since the number density N can be controlled by modifying the pressure of the gas [9], it opens a possibility to control the coupling between the dipoles. For example, if the system is prepared by mixing two different types of atoms having closely spaced energy levels, the coupling between the dipoles can be modified by changing the number density of each type of emitters. The ratio of these couplings is where N and N are the number densities of the first and the second types of atoms.

Controlling the dipole-dipole interaction in a dense medium by adjusting material parameters offers an additional way that helps to engineer the collective effects in such systems. Fig. 5.3 shows the result for a mixture of two two-level quantum emitters of very close transition energies. A simple change of the relative densities of the two types of quantum emitters changes the position of the spectral hole as expected. The blue solid curve shows the reflected pulse envelope for 2 13 = 10 2 12 . For this case the spectral hole appears near the transition frequency of the first quantum emitter. The red line is for 2 12 = 2 13 and the spectral hole appears exactly in the middle. For 2 12 = 10 2 13 which is represented with the green line, the reflected pulse is imprinted with the hole near the second quantum emitter transition. Thus by changing the dipole-dipole coupling, the frequency of DIET and hence the spectral hole in the reflected signal can be controlled between the two transitions. The upper and lower limits of this window within which DIET can be controlled are fixed by the choice of quantum emitters. In the case of a multilevel system, DIET will leave a succession of spectral holes in the reflected pulse. Fig. 5.4 shows the reflected pulse shape from such a system with ten transitions regularly spaced by 0.01 eV. The detuning is defined with respect to the first transition energy of the system, which is at 2 eV.

Slow Light

The speed of light is a fascinating quantity. For an electromagnetic field, the dependence of the field propagation on the material parameters is a consequence of Maxwell's equations [10]. The refractive index of a system, which gives access to the speed of the light in a medium, depends on the resonances in the system under consideration. In addition, near the resonances a strong dispersion is observed. Precisely, the refractive index is defined as a ratio of phase velocities. For a pulse propagating in a dense gas, the group velocity v g can be written as

v g (ω) = c n g (ω) , ( 5.3)
where n g is called the group index of the medium and is defined as

n g (ω) = n(ω) + ω dn dω . ( 5

.4)

In general, the dispersion in a medium dn/dω is small. Large values of the dispersion are usually associated with resonances and hence with absorption in the medium. The interesting phenomenon of slow light, that is associated with large values of the group index, is hence not easy to achieve in practice. The slow light effect in physical systems [11,12] was observed in weak pulse propagation through amplifying media in the late 1960's [13,14] and early 1970's [15]. Later, researchers managed to decrease the speed of light to a few tens of m/s and then stopped it [16][17][18]. There are many applications for these effects in different fields including communications, quantum information [19,20], optical memories [21,22], etc... The physics behind slow light is based on electromagnetic induced transparency (EIT). Similar to EIT, DIET can also induce the slowing down of light at frequencies near the overlapping resonances [23,24]. To see this effect using DIET, we consider a 600 nm layer of strongly interacting three-level systems. Keeping all other parameters as in the previous section, the group velocity of light can be estimated. Fig. 5.5 shows the group velocity and the transmitted field for the system considered. The group velocity of the transmitted pulse is here reduced by a factor of 210, i.e, the group index of the model system is 210.

The blue solid curve is the transmitted pulse due to the induced transparency using DIET. The transmission is not very efficient due to the inevitable losses in a dense medium. The red curve shows the group velocity as a function of the reduced detuning. It is clear that in the vicinity of DIET the group velocity is small compared to the speed of light in the vacuum. Even though the model used here is not so efficient as already achieved experiments using EIT, for example [18], it shows that the transmitted part of the incident field at DIET frequency will be traveling with a reduced speed.

It is also interesting to examine the effect in a realistic sample. For that, we consider the D 1 transitions of 85 Rb atoms. We suppose that the pressure is such that the system is in the strong interaction regime, with S SF = 200, so that the interaction between all dipoles in the layer are strong enough to induce associated transparencies via destructive interferences. A similar kind of reduction in the velocity of the transmitted pulse is obtained. Fig. 5.6 shows the relevant results for a layer of thickness 2µm.

Panel (a) shows the group index of the system. The group index is very high, a few 10 6 , similar to recent experimental studies on slow light in plasmonic waveguides [12]. The local minimum in the group index, corresponding to n g ≈ 4.10 6 around δ = 15, corresponds to the DIET transmission window. Within this window, the imaginary part of the refractive index n(ω), shown with the red curve in the panel (b), shows a window of minimum, making the layer less dissipative. The real part of the refractive index, represented by the blue solid curve in panel (b), shows a point of inflation at the frequency of DIET.

The red curve in panel (c) is the transmitted pulse through the layer. It lies exactly in the window of minimum dissipation. Even though the imaginary part of the refractive index and hence the absorption coefficient is small, due to the relatively long distance ( = 2 µm) traveled in the layer, the transmittance is reduced to 45%. The group velocity, represented by the blue curve in panel (c), is very small, a few tens of m/s. A zoom of the plot showing this very small group velocity is given in Fig. 5.7. At the maximum of the group index, the group velocity drops below 20 m/s. At the peak of the transmission, the group velocity reaches 80 m/s. The transit time for the pulse through the layer which is defined in general as T = /v g is also an interesting quantity to look at. Being the reciprocal of the group velocity, the transit time behaves in the similar way as the group index as a function of the reduced detuning and the value of the transit time for a particular value of the reduced detuning is given by the group index multiplied by the factor /c. For the pulse considered here, the transit time in the DIET regime lies around 30 ns, to be compared with a transit time of about 6 fs in the vacuum.

Conclusion

In this chapter, two potential applications of DIET were discussed.

Pulse shaping was discussed in the case of a three-level system and of a multilevel system in the strong interaction regime corresponding to a dense gas. A simple way for controlling the position of the spectral hole in the reflected pulse was also discussed by varying the relative pressure of two different atomic gases. By changing the position of the spectral hole that corresponds to DIET, it is possible to exclude the desired frequency from the deformed reflected pulse. This may be of interest for different purposes.

The slowing down of the speed of light inside the layer can be seen in systems showing DIET. Essentially, the pulse transmitted due to DIET is experiencing a large dispersion and hence a large group index that can slow down the pulse down to a few tens of m/s. Even

Conclusion and Outlook-Part II

This chapter concludes the work described in the second part of the thesis. It also describes some recent results obtained with stronger incident pulses.

Keywords:

Populations and Coherences, Non-Hermitian Wave Packets

Conclusion

In this part of the thesis, the optical properties of multilevel systems confined in a nanometric cell (layer) was investigated in detail within the weak field limit.

The optical response of a dense sample of two-level systems was calculated in both the weak (dilute) and strong (dense) interaction regimes of the constituent particles. In the weak interaction regime, the system is characterized by a weak Lorentzian lineshape that confirms the absence of collective behavior of the constituent particles. A simple Lorentz model describing the system is used for interpreting the physics. Many interesting collective effects are observed in the strong interaction regime, including the enhancement and broadening of selective reflection due to the coupling between the dipoles.

In the presence of additional closely spaced quantum states, the response of the system changes dramatically. Broadened interactions due to dipole-dipole couplings induce an overlap of different resonances, leading to the cancellation of the dipoles in the system at a specific wavelength. The field then propagates through the medium, with minimum refractive effects. This transparency induced due to the dipole coupling in the system was named "dipole-induced electromagnetic transparency" (DIET).

These effects were explained using model systems and realistic calculations were also designed on the basis of the D 1 transitions of 85 Rb atoms at low temperature. An anticipated transparency is obtained in the case of Rb atoms, giving a possibility for verifying the existence of DIET experimentally. Since there are many groups working on the field using Rb samples, this may be achieved in the near future.

DIET resembles the well-known electromagnetically induced transparency process (EIT). Cancellation of the dipoles happens at a frequency between the two nearest transitions in the system. The frequency at which the transparency occurs can be controlled by changing the relative number density of different types of atoms. Just like in EIT, the light passing through the system is delayed due to a high dispersion effect. This slow light effect was calculated in the case of Rb atoms.

DIET modifies the reflected spectrum, which shows one or several spectral holes. This effect can be used for excluding desired frequencies from a laser pulse. This type of pulse shaping can be of experimental interests and it was discussed in the case of multilevel model systems.

All these effects were discussed exclusively within the limits of a weak probe. Moving out of this limit will affect the system significantly by populating excited states with a significant probability and by turning on different kinds of responses to the applied fields. These effects are briefly introduced in this last chapter of the thesis.

Non-Hermitian Model Revisited

For treating the system quantum mechanically, either Bloch equations or non-Hermitian wave packet equations [1] can be used. The later was of great interest while dealing with multilevel systems since the method used speeds up the calculations. But the equations for the non-Hermitian wave packet propagation will be diverging in strong fields due to the large number of excitations in the system. The equations describing the coherences of a multilevel system in the non-Hermitian wave packet model were given in Eq. (2.48). We rewrite them here: The problem for extending this model beyond the weak probe limit is clear from Eq. (6.2).

i ∂C 1 ∂t = ω 1 + i γ 1 (t) 2 C 1 +  2 Ω 1 (t) C  (6.1a) i ∂C  ∂t = Ω 1 (t) C 1 + ω  -i γ  (
In the presence of strong fields, the excitation probabilities are not negligible. Since the definition of γ ı (t) depends on the population difference, these quantities will diverge as the population of the excited states grows beyond 50%. This divergence is inevitable in this model. It can however be avoided if one relaxes the constraint of having a normalized wave function. In this case, both the ground and excited state populations can vary accordingly so that the pole is eliminated from the calculation. Relaxing the condition on the conservation of the norm leads to another non-Hermitian model that, as we will see later, is valid in the strong field limit. With this approach, and in contrast with Eq. (6.2), it is clear that the rates γ 1 (t) and γ 2 (t) should not diverge anymore.

6.2.1

Comparison of the Different Non-Hermitian Models

We now compare the different models for different field strengths. A comparison between the optical Bloch equations and the non-Hermitian model obtained for norm-conserved wave packets was already detailed in Chapter 2. Here the calculations are performed for a 600 nm nano-layer composed of two levels systems with a transition energy of 2 eV and with S 12 = 1.26. But the case will be different as the incident field intensity increases so that the system is disturbed significantly. In this case the condition |C 1 (t)| 2 |C 2 (t)| 2 is no more valid. Fig. 6.3 shows the population and the coherence and Fig. 6.4 shows the optical response calculated from the three models describing the dynamics of the two-level system exposed to a strong field. A 10 fs Gaussian pulse of amplitude E 0 = 5 × 10 9 V/m ≈ 0.01 au is used. In the presence of such a strong field, the optical Bloch equations are used as a reference for comparing the results obtained from the non-Hermitian wave packet models.

The thick blue curves in Figs. 6.3 and 6.4 are the results obtained from the optical Bloch equations. The green thin curves are those obtained from the non-Hermitian model-1 (with conservation of the norm) and the red dashed curves are those from the non-Hermitian model-2 (revised model). The panel (a) in Fig. 6.3 shows the population of the excited state calculated from the three models. None of the two wave packet methods are able to describe accurately the evolution of this population even if the second method is doing a slightly better job. An interesting feature for the population calculated from model-2 is that its variation is similar to the one of the coherence of the system which is given in panel (b). The coherence calculated from model-1 disagrees clearly with the optical Bloch equations while the coherence calculated from model-2 is very accurate. A similar conclusion can be reached from the real part of the coherence.

The well-behaved coherence obtained from model-2 in strong fields ensures a better accuracy for the calculation of the response of the system compared to model-1. Fig 6.4 shows the comparison between the responses obtained from the different models. Panel (a) shows the reflection, panel (b) shows transmission and panel (c) shows the extinction. The color codes used in this plot are the same as in Fig. 6.3. The non-Hermitian model-1 shows a clear disagreement with the optical Bloch equations while the model-2 reproduces accurately the "exact" results. This method could therefore be used to predict the response of the system in strong fields.

Outlook

DIET is of experimental interest. It can be used in many ways for practical purposes. The system considered in this thesis is a uniform distribution of quantum emitters. The pulse transmitted through DIET systems constitutes slow light. In the case of 85 Rb, we have shown that the transmitted pulse can travel with velocities of a few tens of m/s. Thus the addition (or the presence) of a mechanical wave (oscillation), like a sound wave for instance, could affect the quality of DIET. For example, the action of a longitudinal wave will set the constituent particles in oscillation in the direction of propagation. If the acoustic wave and the applied electromagnetic fields are propagating in the same direction, the emission of dipolar fields will happen while the dipoles are oscillating. Since the group velocity of the electromagnetic field is very small, there may be an important influence on DIET either in terms of frequency or in terms of efficiency. It is a problem that could be of interest for the future.

The development of the non-Hermitian code for a multilevel system in the strong field regime is also interesting. Currently the model, named non-Hermitian model-2, works well for two-level systems. The problem lies with the choice of the empirical decay and decoherence rates introduced in the model. An appropriate definition of these empirical quantities could lead to an accurate multilevel model. It could then be used for studying strong field phenomena in dense gases. Some groups reported high harmonic generation in dense systems. It may be of interest to look at similar systems to explore more about the strong field effects in a high density sample of quantum emitters. The advantage of having a wave packet model is that it is computationally much more efficient than a time dependent density matrix approach when many quantum states are populated.

III

Appendices

A

Atomic Unit System

It is important to use an appropriate unit system that ease the calculations we want to do. The most commonly used unit system, the Système International (SI) is not always practical. For example, to describe physical quantities in relativistic calculations, string theory, etc.. the Lorentz-Heaviside system of units is more useful. When we deal with quantum gravity, the most practical unit system is based on Planck units. Atomic units (a.u.) is a system of units based on universal atomic physical constants related to the ground state of the hydrogen atom. They simplify the equations and calculations of physical quantities in atomic, molecular and optical (AMO) physics.

The standards of atomic units are described below.

A.1 Concept of Atomic Units

From the fundamental concept of charge in the low energy regime, the charge that can be carried by any particle is an integer multiple of the electron charge e which is about 1.60217657 × 10 -19 C. Or in other words, in atomic and molecular physics, charges can be specified as units of the fundamental charge e. Thus it is convenient to set e as the unit of charge when dealing with atoms and molecules. Therefore, e = 1.60217657 × 10 -19 C = 1 a.u. of charge.

(A.1)

In the SI system, the time-dependent Schrödinger equation (TDSE) for a one electron atom with infinite nuclear mass can be written as

i ∂ t ψ = - 2 2m e ∆ψ + V ψ (A.2)
where m e is the electron mass. If we rescale the coordinates associated with the equation in atomic units: c ≈ 137 a.u. of speed.

The table below summarizes the atomic unit system.

B

Runge-Kutta Method

The optical Bloch equations and the non-Hermitian wave packet approach used for describing the system quantum mechanically are a set of first-order differential equations of the form dy/dx = f (x, y). First order differential equations of this form can be solved numerically by the well-known fourth-order Runge-Kutta method (the RK-4 method) provided the initial conditions are known (i.e, y(x 0 ) = y 0 ). This method was developed by the mathematician C. Runge and extended by another mathematician, M. W. Kutta.

The differential equation that should be solved is

dy dx = f (x, y(x)) (B.1)
The value of y at a point x + h can be written as In those cases where the integration in Eq.(B.2) cannot be done analytically, the only way to do it is by looking at possible ways to calculate them numerically. For this numerical purpose, the integration can be replaced by the quadrature of order N y(x + h) ≈ y(x) + h

N i=1 C i f (x + α i h, y(x + α i h)) (B.3)
where C i and α i are some numbers that will be determined later. In principle, the summation is an infinite sum that is truncated at a certain finite value N for numerical reasons. The accuracy on y(x + h) strictly depends on the order of this quadrature.

Solving the quadrature is again difficult since it depends on the values of y at other positions (x + α i h). The unknown quantities α i and C i can be estimated using a Taylor expansion of Eq.(B.3). For the first term we get

h C 1 f (x + α 1 h, y(x + α 1 h)) = C 1 K 1 (B.

Fig. 1 .

 1 Fig.1.1 shows as red dots the graphical solutions for the equation |r(t) -r 0 | = 0. Depending on the ionization time, the number of solutions will vary. They give the times at which the electron in the laser field returns to its initial position. The return of electron to its initial position is referred as a recollision[6]. If the electron is ionized from an atom

Figure 1 . 1 .:

 11 Figure 1.1.: Typical classical trajectory of an electron in the presence of an oscillating field.

Figure 1 . 2 .

 12 Figure 1.2.: Quasi-static picture of an atomic potential coupled to an intense infrared laser field. Panel (a): Tunneling regime. Panel (b): Above the barrier ionization.

0Figure 1 . 3 .:

 13 Figure 1.3.: Three step model. At t = 0: The bound state electron. At t = π/(2ω L ): Formation of the potential barrier and ionization of an electron. At t = π/ω L : Excursion of the ionized electron in the laser field. At t = 3π/(2ω L ): The ionized electron is driven back to the parent ion. At t = 2π/ω L : Recollision with the parent ion.

Chapter 1 .Figure 1 . 4 .:

 114 Figure 1.4.: Classical trajectories of electrons ionized at different intervals of time. Electron born at earlier time is shown in blue and that born late in the cycle is shown in red. Other electrons born in between are shown in order with colors from blue to red

Figure 1 . 5 .:

 15 Figure 1.5.: Kinetic energy at the moment of recollision obtained from the classical equations of motion for ionized electrons in the laser field as a function of the ionization time.

Figure 2 . 1 .

 21 Figure 2.1.: CO 2 molecule aligned along the y axis. A linearly polarized (in the x direction) electric field is applied to the molecule.

Figure 2 . 2 .

 22 Figure 2.2.: HOMO & HOMO-1 of the CO 2 molecule.

1 ∆Figure 2 . 4 .:

 124 Figure 2.4.: Convergence scheme chosen in the ITP. Panel (a): Optimization of the relative error (2.24). A zoom of a portion of the graph is given in the inset. Panel (b): Convergence of the state energy (2.23). A zoom of a portion of the graph is given in the inset.

Figure 2 . 5 .:

 25 Figure 2.5.: Square modulus of the HOMO and HOMO-1 of a symmetric CO 2 molecule with R = 2.5 Å calculated by imaginary time propagation.

Figure 2 . 6 .:

 26 Figure 2.6.: Typical normalized electric field E(t)/E 0 from Eq.(2.26) used for inducing the dynamics in the system.

2 Figure 2 . 7 .:

 227 Figure 2.7.: Schematic drawing of LIED from a CO 2 molecule.

Figure 2 . 8 .:

 28 Figure 2.8.: Time dependence of the HOMO and HOMO-1 ionization signals for a pulse of carrier wavelength 2.1 µm and intensity 10 14W/cm 2 . The electric field is given as a blue solid line and the corresponding ionization probabilities are shown as a red dotted line for the HOMO and as a red solid line for the HOMO-1.

Figure 3 . 1 .

 31 Figure 3.1.: Typical Photoelectron spectrum S(k x , k y ) in linear (a) and log (b) scale for a square pulse consisting of a single optical cycle with a wavelength of 2.1 µm and an intensity of 10 14 W/cm 2 .

Fig. 3 .

 3 Fig.3.1 (b) contains bright and dark regions, which are due to the interference between the ionized wave packets. The horizontal diffraction patterns seen in the spectra are due to the interference between the different components of the electron wave packet that are scattered from different atoms in the molecule.

Fig. 3 .

 3 Fig. 3.2 gives an emphasized zoom of the LIED spectrum shown in Fig. 3.1 (b). The black parallel dashed lines show the multiple slit-like interference due to the superposition of scattered wave packets from different atomic centers.

Figure 3 . 2 .:

 32 Figure 3.2.: Zoom of the photoelectron spectrum S(k x , k y ) given in Fig. 3.1 (b) showing interferences between electron wave packets in momentum space. The dashed horizontal lines mark some multiple-slit like interference patterns due to the interference between wave packets rescattered by different atomic centers. The faint circular dotted curves mark interferences between long and short trajectories.

Figure 3 . 3 .:

 33 Figure 3.3.: Holographic structure in the LIED spectrum S(k x , k y ).Recolliding wave packets interfere with the directly ionized wave packet creating rich interference patterns around the region where the directly ionized electrons are distributed.

Fig. 3 .

 3 Fig.3.3 (a) shows the log scale LIED spectrum obtained for such a half-cycle calculation. In the case of a full cycle excitation, holographic interferences take place in this energy range. They can be easily identified by comparing the photoelectron spectra of Panel (a) and Panel (b), which shows the log scale photoelectron spectrum with a full cycle where the holographic interference is taking place. The branchings appearing in the LIED spectrum given in Fig.3.3 (b) are the holographic structures. It should be possible to get information about the molecules from these structures but they are very hard to measure because they are shifted for molecules seeing different field intensities because located at different positions within the focal volume of the laser beam.

Fig. 3 .

 3 Fig. 3.4 shows for instance the 2D LIED spectrum S(k x , k y ) (in log scale) obtained for a single optical cycle pulse of λ = 800 nm, 1300 nm and 2100 nm in panels (a), (b) and (c) respectively. All spectra are plotted with the same color map, i.e., between 10 -8 and 1.

Figure 3 . 4 .

 34 Figure 3.4.: LIED spectrum S(k x , k y ) obtained for a single optical cycle pulse and an internuclear distance R = 3 Å. The field intensity is 10 14 W/cm -2 . Panel (a) is for 800 nm, (b) is for 1300 nm and panel (c) is for 2100 nm.

Figure 3 . 5 .:

 35 Figure 3.5.: Normalized average photoelectron spectrum S(k y ) for the same laser parameters as in Fig. 3.4. The blue, red and green solid curves are the spectra for 800 nm, 1.3 µm and 2.1 µm.

Figure 3 . 6 .

 36 Figure 3.6.: LIED spectrum obtained for a single optical cycle pulse. The blue circle represents the maximum recolliding energy of 3.17U p .

Fig. 3

 3 .7 depicts the LIED spectra S(k x , k y ) for multi-cycle pulses. The laser parameters used are the same as in Fig. 3.6 except for the pulse duration. Panel (a) is the LIED spectrum for a 2-optical cycle pulse while (b) and (c) are for 3-and 4-optical cycle pulses.

Figure 3 . 7 .

 37 Figure 3.7.: LIED spectra S(k x , k y ) obtained for (a) 2 optical cycles (b) 3 optical cycles and (c) 4 optical cycles.

Figure 3 . 8 .:

 38 Figure 3.8.: Normalized electric fields E(t)/E 0 with CEP φ = 0 (blue curve) and CEP φ = π/2 (red curve) in the case of a single optical cycle laser pulse.

Figure 3 . 9 .

 39 Figure 3.9.: LIED spectra S(k x , k y ) obtained for a single optical cycle for the CEPs φ = 0 (a) and φ = π/2 (b). Panel (c) shows the averaged 1D spectra S(k y ): The blue curve is for φ = 0 and the red curve is for φ = π/2.

Figure 3 .

 3 Figure 3.10.: 1D spectra S(k y ) obtained for a 4-optical cycle pulse with CEP φ = 0 (blue curve) and CEP φ = π/2 (red curve).

Figure 3 .

 3 Figure 3.11.: LIED spectra from the HOMO-1 of CO 2 . Panel (a) is the 2D spectrum S(k x , k y ) obtained for a single optical cycle pulse and panel (b) shows the same for a 4-optical cycle pulse. Panel (c) shows the corresponding averaged 1D LIED spectra S(k y ). The blue curve is for one optical cycle and the red curve is for 4 optical cycles. The wavelength is 2.1 µm and the intensity 10 14 W/cm 2 .

Figure 3 .

 3 Figure 3.12.: Photoelectron spectra for different internuclear distances: The left column is for the HOMO and the right column is for the HOMO-1 of CO 2 . The wavelength is 2.1 µm, the intensity is 10 14 W/cm 2 with a 4 optical cycle pulse.

Figure 4 . 2 .

 42 Figure 4.2.: Panel (a): 2p x Slater type orbital (STO). Panel (b): 2p x Gaussian orbital (GTO) for α 0.2134 ζ 2 .

Figure 4 . 3 .:

 43 Figure 4.3.: Schematic of the recollision process in CO 2 .

.69) 4 . 6 . 2 . 2

 4622 Approximate LIED Spectrum of the HOMO-1

Fig 4 . 4

 44 Fig 4.4 shows the 1D spectra for the HOMO and for the HOMO-1, S HOMO (ky) and S HOMO-1 (ky), as a function of k y , for an arbitrary set of parameters. The HOMO spectrum, shown in blue, has a regular periodicity similar to that of the TDSE results (see for instance Fig.3.10 for the HOMO). This periodicity is related to the sin 2 (k y R) function in the approximate model of Eq. (4.69). The HOMO-1 spectrum, shown in red, has a double periodicity consisting of relatively high peaks separated by small peaks. These are the characteristic features of the LIED spectra from the HOMO-1, as observed in the TDSE calculations (see for instance Fig.3.11 for the HOMO-1).

Figure 4 . 4 .

 44 Figure 4.4.: 1D Approximate spectra for the HOMO and HOMO-1 of CO 2 . Both signals are estimated for an arbitrary set of parameters.

Fig. 4 .

 4 Fig. 4.5 shows the 1D spectrum S HOMO (ky) and the fitted model for obtaining R and α.The result of the TDSE is given in blue and the fit obtained with our approximate formula (4.69) is shown in red. The value of R that we have extracted from the fitting procedure is 5.12 Å. The error obtained in the extraction of R is therefore of 2.4 % only.

Figure 4 . 5 .:

 45 Figure 4.5.: Fit obtained for the HOMO in the case of R = 5 Å. The blue solid line is the TDSE result and the red dashed line is the model fitted to the calculated spectrum. The wavelength used is 2.1 µm for a total pulse duration of 4 optical cycles and an intensity of 10 14 W/cm 2 .

  (a) and (b) of Fig. 4.7. This can be verified by calculating the same differences for the fitted values α = 0.92 a.u and ζ 2.08 a.u and for the exact internuclear distance R = 5 Å. Panels (c) and (d) of Fig. 4.7 show the result of this additional estimation. Compared to panels (a) and (b) of Fig. 4.7 where both errors in R and ζ or α are considered, panels (c) and (d) of the same figure explicitly show the discrepancy in the reconstructed orbital due to the error in the estimated value of ζ or α. The distortion induced due to the error in R can be seen by comparing both figures.

Figure 4 . 6 .

 46 Figure 4.6.: Reconstructed HOMO for R = 5 Å. Panel (a) is the initial state used in the TDSE calculation of the LIED spectrum shown in Fig. 4.5. Panel (b) is the reconstructed Slater HOMO orbital and panel (c) is the reconstructed Gaussian HOMO orbital.

Figure 4 . 7 .:

 47 Figure 4.7.: Difference between the TDSE initial state and the reconstructed orbitals: Panel (a) shows (Φ TDSE 0 -Φ g 0 ) and panel (b) shows (Φ TDSE 0 -Φ s 0 ). Panels (c) and (d) show the same quantity but for the exact internuclear distance R = 5 Å.

Fig. 4 .

 4 Fig. 4.8 shows the 1D LIED spectra for the HOMO at different internuclear distances. The blue solid curves are the numerical TDSE results and the red dashed lines are the fit of our analytical model. Panel (a) shows the results for R = 2 Å, (b) for R = 3 Å, (c) for R = 4 Å and (d) for R = 5 Å.

Figure 4 . 8 .:Figure 4 . 9 .

 4849 Figure 4.8.: Calculated and fitted 1D LIED spectra for the HOMO orbital at different internuclear distances. Panel (a) is for R = 2.0 Å, (b) shows R = 3.0 Å, (c) shows R = 4.0 Å, and (d) is at R = 5.0 Å. The blue solid lines are the numerical results of the TDSE and the red dashed lines are the results of the analytical model. The wavelength used is 2.1 µm for a total pulse duration of 4 optical cycles and an intensity of 10 14 W/cm 2 .

Fig. 4 .

 4 Fig. 4.10 shows the different fits obtained. The parameters obtained from the plots are used for the reconstruction of the HOMO-1 and for the calculation of the overlaps. The retrieved values of bond lengths are in very good agreement with the values used in the TDSE calculations. The orbitals reconstructed are also of good quality, as shown by the values of the overlaps in table4.2. This table gives the numerical values of all relevant quantities obtained from the fitting procedure.

Figure 4 .

 4 Figure 4.10.: Calculated and fitted LIED spectra for the HOMO-1 at different internuclear distances. Panel (a) is for R = 2.0 Å, (b) shows R = 3.0 Å, (c) shows R = 4.0 Å, and (d) is at R = 5.0 Å. The blue solid lines are the TDSE results and the red dashed lines are the analytical model.

Figure 4 .

 4 Figure 4.11.: Reconstructed HOMO-1 orbitals at different internuclear distances. The left column contains the initial states used in the TDSE. The middle column contains the reconstructed Gaussian orbitals and the right column is for the reconstructed Slater type orbitals.

Chapter 5 .

 5 Conclusion & Outlook-Part I

. 20 )

 20 so that the associated density matrix isρ = |Ψ Ψ| =

12 Figure 2 . 2 .

 1222 Figure 2.2.: N -level system of transition frequencies ω 1 coupled to an electromagnetic field of limited duration. The action of the electromagnetic pulse can populate more than one excited state.

Figure 2 . 3 .:

 23 Figure 2.3.: Schematic view of the system: a nano-layer of thickness is exposed to a normal incident ŷ-polarized electric field represented in blue in this figure. The excitation dynamics of the quantum emitters is solved and the calculated reflected (in red) and transmitted (in dark green) fields are used to extract the reflection and transmission spectra associated to this system.

Figure 2 . 4 .:

 24 Figure 2.4.: Yee cell used in the simulations of the 1D electromagnetic field equations(2.49). The grid in red is for the electric field and the one in black is for the magnetic field. The grids for the electric and magnetic fields are shifted by ∆x/2 in space and ∆t/2 in time. The electric field is calculated at spatial points (n x -1)∆x, n x ∆x, (n x +1)∆x, etc... for the times (n t -3/2)∆t, (n t -1/2)∆t, (n t +1/2)∆t, etc... The magnetic field is calculated at spatial points (n x -3/2)∆x, (n x -1/2)∆x, (n x + 1/2)∆x, etc... for the times (n t -1)∆t, n t ∆t, (n t + 1)∆t, etc...

Figure 2 . 5 .

 25 Figure 2.5.: Schematic view of the numerical implementation: a nano-layer of thickness , exposed normal to a ŷ-polarized electric field, is placed in a spatio-temporal grid. The incident, reflected and transmitted fields are calculated by solving Maxwell's equations using the finite-difference time domain method. Polarization inside the system is calculated using one of the models described in the previous sections. Perfectly matched layers are used at the boundaries of the simulation region in order to avoid unphysical reflections of the electromagnetic fields back to the simulation region.

Figure 2 . 6 .:

 26 Figure 2.6.: Scheme used for the numerical implementation of the nano-layer and electromagnetic field dynamics.

Figure 2 . 7 .

 27 Figure 2.7.: Panel (a) shows the populations and panel (b) shows the coherence as a function of time for a 400 nm thick layer of two level atoms of density N = 2.5 × 10 21 cm -3 . The solid blue lines are obtained from non-Hermitian wave packet propagations and the dashed red line are from optical Bloch equations (see text for details).

Figure 2 . 8 .:

 28 Figure 2.8.: Dependence of the ground state gain factor γ 1 (t) of the two-level system described in the non-Hermitian wave packet model. All parameters are the same as in Fig. 2.7.

Figure 2 . 9 .:

 29 Figure 2.9.: Comparison between the CPU times taken by the non-Hermitian wave packet model and by optical Bloch equations as a function of the number of states involved. The blue line shows the time taken by the non-Hermitian wave packet propagation and the red line is the time taken by optical Bloch equations.

11 )

 11 Typical variations of the real and imaginary parts of χ e (ω) are shown figure 3.1.

Figure 3 . 1 .:

 31 Figure 3.1.: Real and imaginary parts of the susceptibility in a two-state system as a function of the adimensional detuning.

Figure 3 . 2 .:

 32 Figure 3.2.:Reflection spectra from a layer of identical two-level quantum emitters of thickness = 600 nm as a function of the reduced detuning δ and interaction strength S 12 . γ is 10 THz and Γ = 1 THz. Reflection from a layer of weakly interacting identical dipoles is feeble. Those are therefore multiplied by some constant factors given alongside with the same color codes for a better comparison.

2 EFigure 3 . 3 .

 233 Figure 3.3.: Extinction (panel a) and transmission (panel b) probabilities for the weak (blue solid line) and strong (brown solid line) interaction regime specified in Fig. 3.2. The y-scale on the left is for weakly (S 12 = 10 -3 ) interacting dipoles and the one on the right is for strongly (S 12 = 12) interacting dipoles (see text for details).

p 1 p 2 Figure 3 . 4 .:

 234 Figure 3.4.: Dipole fields produced by two arbitrary oscillating dipoles p 1 and p 2 . The presence of the two oscillating fields keep them coupled.

Figure 3 . 5 .:Figure 3 . 6 .:

 3536 Figure 3.5.: Susceptibility of a dense two-level system: the blue solid line is the susceptibility given by Eq.(3.24) and the red dashed line is the same calculated by integrating the optical Bloch equations numerically.

Figure 3 . 7 .:

 37 Figure 3.7.: Skin depth ζ(ω) as a function of the reduced detuning δ for S 12 = 12. ζ(ω) is minimum at the shifted resonance frequency ω = ω 12 of Eq.(3.32). The skin depth at ω 12 is very small (≈ 15 nm) due to the collective response of the two-level system.

Figure 3 . 8 .:

 38 Figure3.8.: Schematic view of multiple reflections and transmissions in the system: Reflection is measured at a point x R in region I and transmission is measured at a point x T in region III. Region II covers the nano-layer of thickness . Linearly polarized light is incident on the interface between I and II that results in primary reflection and transmission. A part of the field is reflected back and forth inside the layer in region II and contributes to the reflection and transmission measurements.

  3.6) and in the transmission spectrum (see Fig.3.3, panel b).

3 . 49 )Figure 3 . 9 .:

 34939 Figure 3.9.: Comparison between the calculated reflection spectrum obtained from the numerical solution of the optical Bloch equations and the reflection probability obtained from Fabry-Pérot interference scheme (Eq. 3.50) as a function of the reduced detuning δ for S 12 = 12.

Figure 3 .

 3 Figure 3.10.: Relative error in the electric susceptibilities |χ e (ω)| obtained from the analytical model when compared to the numerical integration of the optical Bloch equations as a function of the reduced detuning δ. The blue curve is for S 12 = 10 -3 and the red curve for S 12 = 12.

7 )

 7 If we divide this equation by W 1 (ω)ε 0 Ẽ(ω) we can finally extract the partial susceptibilityχ 1 (ω) associated with the |1 → | transition χ 1 (ω) = 2 1

Figure 4 . 1 .:

 41 Figure 4.1.: Reflection spectrum for S 12 = S 13 ≈ 12 as a function of the reduced detuning defined with respect to the first excited state.The second excited state is shifted by ω 13 -ω 12 = 48 THz. The red line shows the expected reflection probability in the case of independent dipoles and the blue line shows the same for coupled dipoles. The black dashed vertical lines show the positions of the shifted eigen-frequencies ω -and ω + .

  |χ(ω)| (Log scale)

Figure 4 . 2 .:

 42 Figure 4.2.: Susceptibility of the three-level system for S 12 = S 13 ≈ 12 as a function of the reduced detuning defined with respect to the first excited state. The two peaks correspond to the positions of the shifted eigenfrequencies ω -and ω + . The first one has a Lorentzian shape and the second one is a Fano profile which is due to the presence of overlapping transitions.

Figure 4 . 3 .:

 43 Figure 4.3.: Transmission spectrum for S 12 = S 13 ≈ 12 as a function of the reduced detuning. The blue line is for a layer of thickness = 600 nm and the red line is for = 300 nm. The total decay rate is γ m = 10.5 THz. The black dashed vertical line shows ω h calculated from Eq. (4.12). The peaks in the middle correspond to DIET.

Figure 4 . 4 .:

 44 Figure 4.4.: Transmission spectrum for γ m = 5.5 THz as a function of the reduced detuning for a layer of thickness = 300 nm.

5 2 S 1 / 2 5 2 P 1

 21221

Figure 4 . 5 .

 45 Figure 4.5.: 85 Rb-D1 Transitions

Figure 4 . 6 .Figure 4 . 7 .

 4647 Figure 4.6.: Panel (a): Absolute value of the susceptibility in log scale. Panel (b): Reflection spectrum from a 600 nm thick layer of 85 Rb atoms with negligible Doppler and collisional broadening with S SP = 3.8×10 -3 . The transitions are marked with red dashed vertical lines and with the corresponding (F, F ) values

Figure 4 . 8 .

 48 Figure 4.8.: Panel (a): Absolute value of the susceptibility (log scale) of 85 Rb atoms at S SP = 200. Panel (b): Reflection spectrum. Panel (c): Transmission spectrum. The width of the sample is 600 nm.

Figure 5 . 1 .:

 51 Figure 5.1.: Schematic view: Pulse shaping.

3 Figure 5 . 2 .:

 352 Figure 5.2.:The green line is the incident pulse shape normalized with respect to the peak intensity and the red line is the normalized reflected pulse envelope. The corresponding transmission through the layer is shown in the inset.

13 Figure 5 . 3 .: 12 = 2 13 and the green line is for 2 12 < 2 13 .

 13531213 Figure 5.3.:Reflected pulse shapes from a mixture of two two-level systems. The blue line is the reflection for2 12 > 2 13 , the red line is for2 12 = 2 13 and the green line is for 2 12 < 2 13 . Two vertical lines show the transition frequencies ω 12 and ω 13 .

Figure 5 . 4 .:

 54 Figure 5.4.: The blue line is the incident pulse shape normalized with respect to the peak intensity and the red line is the normalized reflected pulse envelope.

Figure 5 . 5 .:

 55 Figure 5.5.: The blue line is the transmitted pulse shape normalized with respect to the peak intensity and the red line is the group velocity as a function of the reduced detuning.

Figure 5 . 6 .

 56 Figure 5.6.: Panel (a): Group index n g (ω). Panel (b): Real and imaginary parts of the refractive index n(ω) of the system in the vicinity of DIET. Panel (c): Group velocity (left y-axis) and normalized transmitted pulse (right y-axis) for a layer of 2µm of 85 Rb atoms.

Figure 5 . 7 .:

 57 Figure 5.7.: Group velocity as a function of the reduced detuning δ for the D 1 transitions of 85 Rb atoms.

Fig 6 . 1 Figure 6 . 1 .:Figure 6 . 3 .:

 616163 Figure 6.1.: Response of a 600 nm layer of two-level systems calculated using optical Bloch equations (thick blue curve), the non-Hermitian wavepacket model-1 (thin red curve) and the non-Hermitian wavepacket mode-2 (dashed green curve) as a function of the reduced detuning for S 12 = 1.26. Panel (a): Reflection. Panel (b): Transmission. Panel (c): Extinction, with a weak 10 fs pulse.

Figure 6 . 4 .:

 64 Figure 6.4.: Optical response calculated for a 600 nm layer of twolevel systems at S 12 = 1.26 excited by a strong 10 fs gaussian pulse. Panel (a): Reflection. Panel (b): Transmission. Panel (c): Extinction. The blue thick lines are from the optical Bloch equations, the green thin lines are from model-1 and the red dashed curves are from model-2.

  y(x + h) = y(x) + x+h x f (x , y(x ))dx (B.2)

4 )

 4 Interférences quantiques dans la dynamique d'atomes et molécules dans un champ électromagnétiqueMots-clés: Recollision et Champs Intenses, Effets collectifs, Imagerie orbitalaire, Diffraction Electronique Induite par Laser, Transparence Electromagnétique Induite par Interaction Dipôle-Dipôle, Lumière Lente Dans cette thèse, nous discutons deux cas spécifiques d'interaction laser-matière où des effets d'interférences quantiques jouent un rôle clé. Dans la première partie nous montrons qu'on peut obtenir une image d'une molécule par photo-ionisation par champ laser intense infrarouge. On crée ainsi un paquet d'ondes électronique oscillant dans le continuum. Ces électrons subissent une collision élastique avec l'ion moléculaire parent, qui peut être utilisée pour imager la molécule et ses orbitales. Le processus est démontré dans le cas de la HOMO et de la HOMO-1 de CO 2 . La deuxième partie traîte des effets collectifs dans des vapeurs atomiques denses confinées dans des cellules de taille nanométrique. Les atomes répondent alors collectivement à un champ laser incident, ce qui conduit à un décalage et à un élargissement des transitions. Lorsque deux résonances se recouvrent, une interférence quantique destructive se produit à une fréquence caractéristique. Ceci induit un phénomène de transmission électromagnétique induite par interaction dipolaire, ou DIET. Quantum Interferences in the Dynamics of Atoms and Molecules in Electromagnetic Fields Keywords: Recollision and Strong Fields, Collective Effects, Orbital Imaging, Laser Induced Electron Diffraction, Dipole Induced Electromagnetic Transparency, Slow Light In this thesis, two specific cases of lasermatter interaction where quantum interference plays a key role are discussed. The first part deals with imaging a molecule using its own electrons. An intense infrared laser field acting on the molecule can create laser-driven electron wave packets in the continuum. These electrons undergo elastic collision with the parent molecule that can be used for imaging the molecule and its orbitals. The process is demonstrated in the case of both the HOMO and HOMO-1 of symmetric CO 2 . The second part discusses collective effects and overlapping resonances in dense atomic vapors confined in nano-cells. Closely spaced atoms respond collectively to a common exciting field that leads to shifts and broadening of resonances. Once two or more transitions are broad enough to overlap, quantum interference will take place and it leads to destructive interference corresponding to a characteristic frequency determined by the system parameters. It will lead to a transmission peak called "Dipole-Induced Electromagnetic Transparency" or DIET. Université Paris-Saclay Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France

  .32)

	The substitution of the exact outgoing waves |Ψ + k with the plane waves |Φ pw k	and the
	insertion of Eqs. (4.19) and (4.30) in Eqs. (4.23) and (4.25) lead to

  ξ o e i4βR 2 + ξ o + 1 2 + ξ o e iβR 2 2

. (4.72) 4.6.2.3 Presentation of the Approximate 1D LIED Spectra

Table 4 .

 4 1 summarizes the retrieved values of the parameters obtained by the fit, especially α and ζ, for different internuclear distances. It also shows the corresponding overlaps between the TDSE wave function and the GTO and STO reconstructed initial states.

	R	α	ζ	R	|R -R | R	(%) O G (%) O S (%)
	2.0 1.06 2.23 2.13	6.4%	89.1%	92.7%
	3.0 0.98 2.15 3.11	3.6%	90.5%	94.2%
	4.0 0.94 2.10 4.12	2.9%	90.8%	94.4%
	5.0 0.92 2.08 5.12	2.3%	91.7%	94.7%

Table 4 .

 4

1.: Retrieved parameters for the HOMO orbital at different internuclear distances.

Table 4 .

 4

2.: Retrieved parameters for the HOMO-1 orbital at different internuclear distances. Compared to Table 4.1, a column has been added with the value of ξ o .

  Two-level system of transition frequency ω 12 = ω 2 -ω 1 coupled to a monochromatic electromagnetic field E(t) at frequency ω = ω 12 . The ground state |1 is associated to the energy ω 1 and the excited state |2 to ω 2 .

	|2	hω 2
	E(t)	hω 12
	|1	hω 1
	Figure 2.1.:

  The time-dependent Schrödinger equation (TDSE) for a two-level system described by the state vector |Ψ(t) given in Eq. (2.20) with time-dependent coefficients C 1 (t) and C 2 (t)

	can be written as	i	∂ ∂t	|Ψ(t) = Ĥ|Ψ(t) .	(2.37)

  From Eq.(3.42) it is clear that the frequency window having maximum reflection is span over 3∆ 12 which gives a reflection window of width 3∆ 12 (see Fig.3.6).

		10 4
								Semi-classical model
		10 3						Optical-Bloch equations
	|χ e (ω)| 2 log scale	10 0 10 1 10 2				∆ 12
		10 -1
		-100 10 -2	-80	-60	-40	-20	0	20	40	60	80	100
						Reduced Detuning δ
											12 .	(3.42)

  |C k (t)| 2 |C 1 (t)| 2k 2 |C k (t)| 2 , (6.2a) γ  (t) = (2γ + Γ)|C 1 (t)| 2 |C 1 (t)| 2k 2 |C k (t)| 2 .(6.2b)These equations were derived from the time dependent Schrödinger equation for a non-Hermitian Hamiltonian satisfied by a normalized wave function, with coefficients C ı where ı = 1, 2, 3..N .

		t) 2	C 	(6.1b)
	where
	γ 1 (t) =	(2γ + Γ) k 2

  Improved Non-Hermitian Model for a Two-Level System From Eq. (2.27) it is clear that the coherences depend on the difference in populations between the ground and excited states. Thus it is important to describe the quantity|C 1 (t)| 2 -|C 2 (t)| 2accurately in a two-level system. The aim is therefore to find a set of dynamical equations for C 1 (t) and C 2 (t) that describes this population difference correctly, and not necessarily the norm|C 1 (t)| 2 + |C 2 (t)| 2 .As seen from the previous approach, the condition on the empirical gain and decay factors γ 1 (t) and γ 2 (t) should beγ 2 (t) -γ 1 (t) = 2γ + Γ.(6.3)In addition, an accurate description of the evolution of the population difference requires2Γ |C 2 (t)| 2 = γ 1 (t) |C 1 (t)| 2 + γ 2 (t) |C 2 (t)| 2 (6.4) |C 2 (t)| 2 |C 1 (t)| 2 + |C 2 (t)| 2(6.5a) γ 2 (t) = (2γ + Γ) γ 1 (t) + 2Γ |C 2 (t)| 2 |C 1 (t)| 2 + |C 2 (t)| 2 (6.5b)

	From Eqs. (6.3) and (6.4), we obtain
	γ 1 (t) =	(Γ -2γ)

Table A .

 A 1.: Conversion Factors for Atomic Units ) 2 3 /m e e 4 2.41888432.10 -17 s Energy m e e 4 /(4πε 0 ) 2 4.35974417.10 -18 J Electric Field m 2 e e 5 /(4πε 0 ) 3 4 5.14220652.10 11 V/m Intensity cε 0 m 4 e e 10 /(4πε 0 ) 6 8 3.50944758.10 16 W.cm -2 A short history of the SI units in electricity by Ludwik Kowalski in The Physics Teacher 24, 97 (1986). -Scaling Mount Planck I: A view from the bottom by Frank Wilczek in Physics Today 54, 6 (2001).

	Physical Quantity	Expression in the SI system	Conversion Factor
	Mass	m e (electron rest mass)	9.10938291.10 -31 kg
	Electric Charge	e (electron charge)	1.60217657.10 -19 C
	Angular Momentum			1.05457173.10 -34 J.s
	Coulomb Constant	1/4πε 0	8.98755178.10 9 m/F
	Length	4πε 0	2 /m e e 2	5.29177210.10 -11 m
	Time (4πε 0 Fine Structure Constant e 2 /4πε 0 c	≈ 1/137
	Electric Dipole Moment	4πε 0	2 /m e e	8.47835326.10 -30 C.m
	Magnetic Dipole Moment	e /2m e	9.27400968.10 -24 J.T -1
	For More Details

-

References

(a) Reference[25] gives a nice account on soft-Coulomb potentials

We have shown during this study that the cooperativeness of systems with at least two excited states coupled to a ground state via a laser field can posses an interesting phenomenon which is named as dipole-induced electromagnetic transparency. In the coming chapters of the thesis, this topic will be discussed in detail with a theoretical demonstration by using a realistic system of rubidium atoms.

the particle, compared to their equilibrium positions. This results in a polarization of the system, which in turn modifies the external electric field.

The differential equation for the displacement y(t) of the charged particles from their equilibrium positions due to the action of the local electric field E local (t) can be written as m ÿ = -mω 2 12 y -mγ m ẏ + q E local (t) (3.18) where m and q are the mass and electric charge, respectively, of the particles considered. ω 12 is the frequency of the harmonic oscillations and γ m is the damping constant of these oscillations. In this expression -mω 2 12 y(t) denotes the binding force, -mγ m ẏ(t) is the friction force and q E local (t) is Coulomb's force.

Let us now remember that the polarization is modeled according to Eq.(2.35) as

and that one can relate the dipole operator μ with the displacement operator ŷ μ = q ŷ . (

Assuming now that the expectation value of the displacement operator can be assimilated with its classical counterpart y(t) we obtain a classical equation of motion for the polarization component P y (t) [3] m Py = -mω (ω 2 12 -ω 2 ) + i ωγ m (3.24) where

We see here that the resonance frequency is shifted by a term which depends on the density N and on the ratio q 2 /m. Reminding the argument for introducing the local field correction in Section 2.2.1, one can safely interpret that this shift is induced by the dipole-dipole coupling between the constituent quantum emitters in the system.

The difficulty for estimating the coupling frequency 12 in the system is because of the presence of the parameter q 2 /m which can be estimated by quantizing the dipole oscillations. For this purpose let us consider an oscillating transition dipole µ 12 (t) characterized by the charge q, the amplitude of oscillation x 0 and the frequency ω 12 . We assume for the moment that there is no damping. The length of the dipole is then x(t) = x 0 cos(ω 12 t) and µ 12 (t) = q x(t) = q x 0 cos(ω 12 t). The RMS value of the dipole is therefore given by

If the oscillation energy is quantized to the zero-point energy ω 12 /2 of the harmonic oscillator, it is easy to show that the amplitude of oscillation x 0 is nothing else than the characteristic length of the harmonic oscillator response calculated using the three models match very well in the weak probe regime.

The choice made on the description of the quantity |C 1 (t)| 2 -|C 2 (t)| 2 for excluding the numerical difficulty as the system approaches the population inversion does not bring any noticeable discrepancy. Thus the non-Hermitian wave packet model-2 can be used for studying dense two-level systems within the weak field limit. Chapter A. Atomic Unit System (A.2) using the Bohr radius a 0 , we get ∆ ≡ a 2 0 ∆ and the TDSE becomes

The LHS of equation (A.3) reveals the appearance of a characteristic time scale t 0 = m e a 2 0 / while the RHS of the same equation reveals the appearance of a characteristic energy scale E 0 = 2 /(m e a 2 0 ). If we define t = t/t 0 and V = V /E 0 Eq.(A.3) can be reorganized as the following dimensionless equation

We have the intrinsic scales of length, energy and time, a 0 , E 0 and t 0 that simplify the TDSE. In the atomic system of units, the Bohr radius is defined as 1 a.u. of length and the quantities E 0 and t 0 define the atomic units of energy and time. Note that the ionization energy of the hydrogen atom is 0.5 a.u. of energy and that the Bohr orbital period is 2π a.u. of time.

A.2 Other Units

Comparing equations (A.2) and (A.4), we can see that the definition of atomic units are compatible with the choices = 1 a.u. for the unit of angular momentum and m e = 1 a.u. for the unit of mass.

Defining now the atomic unit of electric field F 0 as the field strength due to the hydrogen nucleus (a proton) at a distance a 0 , we get

The intensity of a linearly polarized electric field of amplitude F 0 can now be written as

One can finaly note that since dimensionless physical constants keep their value unchanged in any unit system, the fine structure constant α = e 2 4πε 0 c ≈ 1 137 is preserved in the atomic unit system. The reciprocal of the fine structure constant therefore gives the speed of light

Chapter B. Runge-Kutta Method

where K 1 = h f (x, y(x)). The second term can be written as follows

where K 2 = h f (x + α 2 , y(x) + α 2 K 1 ). This procedure can be extended to any order i up to i = N . If one uses only two terms from this quadrature, we obtain the second order RK-2 method. Here, for reasons related to a balance between accuracy and numerical efficiency, we use the RK-4 method, and four terms are included in the quadrature. In the case of RK-4, the four expressions of the K i 's are given by