Tailored pump-probe transient spectroscopy with time-dependent density-functional theory: controlling absorption spectra Jessica Walkenhorst,1 Umberto De Giovannini,1,∗ Alberto Castro,2,† and Angel Rubio1,3,‡ 1Nano-Bio Spectroscopy Group and ETSF Scientific Development Center, Departamento de Quimica, Universidad del Pa´ıs Vasco UPV/EHU, Avenida de Tolosa 72, E-20018, San Sebastia´n, Spain§ 2ARAID Foundation - Institute for Biocomputation and Physics of Complex Systems, University of Zaragoza Mariano Esquillor Gomez s/n, 50018 Zaragoza, (Spain) 3Max Planck Institute for the Structure and Dynamics of Matter Hamburg, Germany (Dated: January 19, 2016) 6 Recent advances inlasertechnology allowusto followelectronicmotion atitsnaturaltime-scale 1 withultra-fasttimeresolution,leadingthewaytowardsattosecondphysicsexperimentsofextreme 0 precision. In this work, we assess the use of tailored pumps in order to enhance (or reduce) some 2 given features of the probe absorption (for example, absorption in the visible range of otherwise transparentsamples). Thistypeofmanipulationofthesystemresponsecouldbehelpfulforitsfull n characterization, since it would allow to visualize transitions that are dark when using unshaped a J pulses. In order to investigate these possibilities, we perform first a theoretical analysis of the non-equilibrium response function in this context, aided by one simple numerical model of the 8 Hydrogen atom. Then, we proceed to investigate the feasibility of using time-dependent density- 1 functionaltheoryasameanstoimplement,theoretically,thisabsorption-optimizationidea,formore complex atoms or molecules. We conclude that the proposed idea could in principle be brought to ] r the laboratory: tailored pump pulses can excite systems into light-absorbing states. However, we e also highlight the severe numerical and theoretical difficulties posed by the problem: large-scale h non-equilibrium quantum dynamics are cumbersome, even with TDDFT, and the shortcomings of t o state-of-the-art TDDFT functionals may still be serious for these out-of-equilibrium situations. . t a m I. INTRODUCTION trary, one wants to study the electronic dynamics only, - disentangling it from the vibronic degrees of freedom, d Time-resolved pump-probe experiments are power- then one must perform TAS with attosecond pulses [7], n ful techniques to study the dynamics of atoms and a possibility recently demonstrated [8, 9]. o c molecules: the pump pulse triggers the dynamics, which The theoretical description of these processes, which [ is then monitored by measuring the time-dependent re- involve non-linear light-matter interaction, and the en- 1 sponseoftheexcitedsystemtoaprobepulse. Thetime- suing non-equilibrium electron dynamics, is challenging. v resolution of this technique has increased over the years, Time-dependent density functional theory (TDDFT) 4 and nowadays, it canbe usedto observe theelectron dy- [10–12]isawell-establishedtooltocomputetheresponse 4 namics in real time, giving rise to the field of attosecond of a many-electron system to arbitrary perturbations. 5 physics [1, 2]. Traditionally, the vast majority of TDDFT applications 4 A suitable setup to observe charge-neutral excitations have addressed the first-order response of the ground- 0 isthetime-resolvedphotoabsorptionortransientabsorp- state system to weak electric fields – which can provide . 1 tion spectroscopy (TAS), where the time-dependent op- theabsorptionspectrum,theoptically-allowedexcitation 0 tical absorption of the probe is measured. TAS can energies and oscillator strengths, etc. Nevertheless, the 6 of course be used to look at longer time resolutions: extension of TDDFT to the description of excited state 1 if we look at molecular reaction on the scale of tens spectral properties and its ability to simulate transient : v or hundreds of femtoseconds, the atomic structure will absorptionspectroscopy(TAS)hasrecentlybeendemon- Xi have time to re-arrange. These techniques are thus strated [13, 14]. mainly employed in femtochemistry [3, 4] to observe and r In this work, we are not only interested in simulating a control modification, creation, or destruction of bonds. attosecond TAS of atoms and molecules, but in study- TAS has been successfully employed, for example, to ing the possibility of tailoring the pump to control the watch the first photo-synthetic events in cholorophylls spectra. In fact, the measurement and control of ultra- and carotenoids [5] (a review describing the essentials of fast processes are inherently intertwined: quantum op- this technique can be found in Ref. [6]). If, on the con- timal control theory (QOCT) [15, 16] can be viewed as the inverse of theoretical spectroscopy: rather than at- tempting to predict the reaction of a quantum system to ∗ [email protected] aperturbation, it attemptstofind the perturbation that † acastro@bifi.es induces a given reaction in a given quantum system. It ‡ [email protected] is the quantum version of the more general control the- § [email protected] ory [17–21], which was needed given the fast advances in 2 experimental quantum control [22–31]. where Hˆ is the static Hamiltonian, that describes the systemitselfandF(t)Dˆ isthecouplingtoaprobepulse The possibility of combining QOCT with TDDFT has µ via the dipole operator been established recently [32]. Furthermore it has been shown,thatitcanbeusedtooptimizestrong-fieldioniza- N tion[33],photo-induceddissociation[34]andiscompati- Dˆ =−(cid:88)rˆ(i) (2) ble with Ehrenfest dynamics [35]. Very recently, Krieger µ µ i=1 et al. used TDDFT to study the influence of laser inten- sity,frequencyanddurationonthelaser-induceddemag- inwhichN isthenumberofelectronsinthesystem,and netization process in bulk materials, which takes place µ = x,y,z determines the polarization direction. Note, on time scales of < 20 fs [36]. that the implementation of the coupling via the dipole Here, we take the first steps towards the use of this operator is an approximation and could be removed in combination of TDDFT with QOCT to control excited practical implementations since the theory handles non- statespectraoffinitesystems. Thisideaisverymuchre- dipolar fields. Also note, that we work in the length latedtotheconceptofelectromagneticallyinducedtrans- gauge all over the paper. parency [37, 38]. Control of the absorption spectra may If,atthetimet=T,thesystemhasbeendrivenbythe mean its elimination or reduction, or its increase. Our previous pump E to the state |Ψ[E](T)(cid:105), the complete gedanken setup throughout this paper is the following: dipole-dipole response function for the perturbation at for a certain time interval [0,T] a quantum system is later times is given by: driven by a “classical” pump pulse E(t) whose precise shape can be manipulated. After the pump has ended, χ [E](t,t(cid:48))= (3) the (linear) response of the system to some later per- Dˆµ,Dˆν (cid:104) (cid:105) turbation is calculated. Our goal is to design the shape −iθ(t−t(cid:48))(cid:104)Ψ[E](T)| Dˆ (t),Dˆ (t(cid:48)) |Ψ[E](T)(cid:105), µI νI of the pump pulse in such a way, that the response to somelaterperturbationisoptimal insomegivenway. In where the operators are expressed in the interaction pic- particular,wedemonstratehowthetailoredpumppulses ture. maybeusedtotransformatransparentatomormolecule ThedifferencebetweenEq.(3)andtheequilibriumre- into an excited one that absorbs in the visible. sponse function [39] or the response function of a system This paper is structured as follows. In Section II, we in a many-body eigenstate is that the two times t and analyze the optical linear response of a system in an ex- t(cid:48) cannot be reduced to only one by making use of the cited state, looking at the location and the shape of the time-translational invariance. Since |Ψ[E](T)(cid:105) depends resulting spectral peaks and their time-dependence. In on the pump, the first-order response of the system does Section III we present the theory of quantum optimal explicitly depend on both the pump E and the probe F. control and how it can be applied to optimize spectral It is given by: properties of systems in excited states. We then bring these concepts into application in Section IV. In Sec- (cid:90) t tion IVA we illustrate the conclusions gained in Sec- Dµ(1ν)[E,F](t)= dt(cid:48)F(t(cid:48))χDˆµ,Dˆν [E](t,t(cid:48)). (4) tion II using the analyticly solvable Hydrogen atom. We T then proceed to demonstrate control of the excited state An intuitive physical meaning can be gained from this propertiesinthissystem. InSectionIVBwefinallycom- equation for the response function: it is the first order of bineourmethodologywithtime-dependentdensityfunc- thesystemresponse,ifweconsiderasuddenperturbation tionaltheory, firstfortheHeliumatom, andthenforthe att(cid:48) =T+τ :F(t)=δ(t−(T+τ)),whereτ denotesthe Methane molecule. We conclude our work in Section V. delaybetweentheendofthepumpandtheperturbation: χ [E](t,T +τ)=D(1)[E,δ ](t). (5) Dˆµ,Dˆν µν T+τ This object contains all the necessary information about II. SHORT REVIEW OF OUT-OF-EQUILIBRIUM (PUMPED) the interacting system to compute the absorption of any ABSORPTION SPECTRA givenprobe,aslongasitisweakenoughfortheresponse to be linear. In order to analyze it, it is useful to take the Fourier transform with respect to the variable t, and In pump-probe spectroscopy, the probe may arrive af- expandEq.(3)inaneigenbasisofthestaticHamiltonian ter, during, or even before the pump; in this work, we Hˆ consider a non-overlapping regime in which the probe arrives after the pump has vanished. The time evolution (cid:90)∞ ofasystemaftertheendofthepumpisdescribedbythe |Ψ[E](T)(cid:105)=(cid:88)γ |Φ (cid:105), Hˆ|Φ (cid:105)=ε |Φ (cid:105), (6) j j j j j Hamiltonian (atomic units will be used hereafter): j=1 Hˆ(t)=Hˆ+F(t)Dˆ , (1) obtaining a Lehmann representation for the time- µ 3 dependent non-equilibrium response function: ready pointed out in Ref. [40], in this non-overlapping regime the dependence of the spectrum on both the χ [E](ω,T +τ)= (7) Dˆµ,Dˆν pump-pulseandthedelaytimeentersthroughthemodifi- (cid:88)(cid:90) dµ dν (cid:40) γjγk∗eiωjkτ − γj∗γke−iωjkτ (cid:41) cpaetaikonpoofsitthioenpseaarkeaimntprliintsuidcepsraonpdersthieaspeosfetxhcelumsiavneyly-b:otdhye jm mk ω+ωjm+iΓ/2 ω−ωjm+iΓ/2 system. Second,fortheanalysisoftheeffectofthepump jkm pulse and of the pump-probe delay on the spectrum, we intermsoftheexactenergydifferencesω =ε −ε ,the can distinguish between three cases: (i) the system is in jk k j dipole-matrix elements dµ = (cid:104)Φ |Dˆ |Φ (cid:105) ∈ R and the its ground state, (ii) the system is in an excited eigen- jm j µ m pump-probedelayτ. Asimilarrepresentationinthecase state,(iii)thesystemisinanon-stationarystate,i.e.ina of non-equilibrium spectroscopy has been recently pre- linear combination of non-degenerate eigenstates. In all sented in a similar context [40]. By writing γj =|γj|eiϕj cases we focus on the positive part of the energy range, Eq. (7) turns into which we denote by ω+. The following shape analysis in terms of Lorentzian and Rayleigh contributions is for χ [E](ω,T +τ)= (8) discretepeaksonly. Wedenotethisbyreplacing(cid:80)(cid:82) by(cid:80) Dˆµ,Dˆν (cid:40) (cid:41) (we will comment on the continuum part later on in this (cid:88)(cid:90) dµ dν |γ γ | eiΘkj(τ) − e−iΘkj(τ) section). jm mk j k ω+ωjm+iΓ/2 ω−ωjm+iΓ/2 jkm In case (i) γ = δ and (cid:61)χ [E](ω+) reduces to i i0 Dˆ,Dˆ the usual Lehmann representation for the ground state with spectrum: Θkj(τ)=ϕj −ϕk−ωkjτ. (9) (cid:61)χ [E](ω+)=(cid:88)|d |2L(ω+−ω ). (15) Dˆ,Dˆ 0m 0m Since the absorption depends on the imaginary part of m the response function, we get from Eq. (8): All peaks are positive and have Lorentzian shape. In (cid:61)χ [E](ω,T +τ)= (cid:88)(cid:90) dµ dν |γ γ | case (ii), where γi = δiξ, the system is in an excited Dˆµ,Dˆν jm mk j k eigenstate Φξ, and therefore in the positive energy part jkm of the spectrum we can find both positive and negative (cid:8) peaks of Lorentzian shape: cos(Θ (τ))L(ω−ω )+sin(Θ (τ))R(ω−ω ) kj jm kj jm −cos(Θ (τ))L(ω+ω )+sin(Θ (τ))R(ω+ω )(cid:9) (cid:61)χ [E](ω+)= (16) kj jm kj jm Dˆ,Dˆ (10) (cid:88) |d |2L(ω+−ω )− (cid:88) |d |2L(ω++ω ). ξm ξm ξm ξm ξ<m ξ>m with the Rayleigh peaks Note,thatinbothcases(i)and(ii),thespectrumistime- ω¯ R(ω¯)= (11) independentandhasonlyLorentziancontributions. This ω¯2+Γ2/4 isthemaindifferencetothecase(iii), fornon-stationary states. Inthiscase,thespectrumcanbedividedintotwo and the Lorentzian peaks parts,onetime-independentandoneoscillatorypartdue Γ/2 to interferences between the involved states: L(ω¯)= . (12) ω¯2+Γ2/4 (cid:61)χ [E](ω+,T +τ)= (17) Dˆ,Dˆ The absorption cross section tensor is: (cid:61)χ0 [E](ω+)+(cid:61)χINT [E](ω+,T +τ) Dˆ,Dˆ Dˆ,Dˆ 4πω The equilibrium term consists of the sum over the sta- σ (ω)= (cid:61)χ (ω), (13) µ,ν c DˆµDˆν tionary state spectra of the eigenstates involved, scaled by their occupations: and for a random sample, the absorption will be its ori- entational average, i.e. the absorption coefficient: (cid:61)χ0 [E](ω+)=(cid:88)|γ |2× (18) Dˆ,Dˆ j 1 j σ¯(ω)= Trσ(ω). (14) 3 (cid:88) (cid:88) |d |2L(ω−ω )− |d |2L(ω+ω ) . jm jm jm jm Since we are only interested in the trace, in the follow- j<m j>m ing, we concentrate on the diagonal terms, we will omit hereafter the orientation indexes in order to ease the no- ItspeaksarealwaysofLorentzianshapeanddependnei- tation. ther on time nor on the initial phase difference ϕ −ϕ . j k It is influenced by the pump laser only through the oc- We now analyze Eq. (10) in more detail. First, as al- cupations |γ |2. The phase- and time-dependency of the j 4 spectrum enters through the interference term where x is the full set of quantum coordinates, and E is the control field, an external potential applied to the (cid:88) (cid:61)χIDˆN,TDˆ [E](ω,T +τ)= djmdmk|γjγk|{ system(inourcase,thepumppulse). Inordertoperform j(cid:54)=k;m optimizations the field must be discretized, for example cosΘ (τ)L(ω−ω )+sinΘ (τ)R(ω−ω ) with the help of a sine Fourier basis. In our numerical kj jm kj jm (cid:9) simulations: − cosΘ (τ)L(ω+ω )+sinΘ (τ)R(ω+ω ) kj jm kj jm (19) M (cid:88) E (t)= c sin(ω t) (21) c n n which in turn is governed by the phase differences Θ , kj n=1 which have contributions from both the phase difference ϕ −ϕ at the end of the laser and from its time evolu- where M is the dimension of the optimization search j k tion ω τ. Θ mixes real and imaginary part of the re- space, and c is the set of all the parameters that de- kj kj sponse function and converts Lorentzian line shapes into termine the field: c = c1,...,cM. The frequencies, and Rayleighlineshapesandviceversa. Thisconversionhap- theirmaximumvalueorcut-offfrequency,maybechosen pens periodically with the frequency given by the energy at will. differencesωkj betweentheoccupiedstatesinvolved. The The specification of E, together with an initial value time-dependence of a spectrum is therefore a clear sign condition, Ψ(0) = Ψ determines the full evolution of 0 of a non-stationary state. Experimentally, this periodic thesystem,Ψ[E],viathepropagationoftheSchr¨odinger beatingpatternwasrecentlyobservedbyGoulielmakiset equation. The behaviour of the system must then be al [8]. This demonstrates, how using a pump to imprint measured by defining a “target functional” F, whose aninternalphasedifferenceϕj−ϕk ontoastateandcon- value is high if the system evolves according to our goal, trolling the delay time τ between pump and probe laser and small otherwise. In many cases, it is split into two canbeusedtochangeaspectrum,convertingabsorption parts, F[Ψ,E] = J [Ψ]+J [E], so that J only depends 1 2 1 into emission peaks (and vice versa) as well as chang- on the state of the system, and J , called the “penalty”, 2 ing the overall shape of the lines. In Section IVA we dependsexplicitlyonthecontrolE. RegardingJ ,itmay 1 demonstrate these line shape changes using the example dependonthefullevolutionofthesystemduringthetime of an exactly solvable Hydrogen atom. Furthermore we interval[0,T],oronlyonthesystemstateattimeT,asit demonstrate,howtousealasertocontrolthesefeatures. is the case in this work. Often, the functional is defined Note, that in the discussion above, the lineshape anal- through the expectation value of an observable Oˆ: ysis is valid for isolated peaks without contribution from continuum states. If coupling to continuum states is in- JT[Ψ(T)]=(cid:104)Ψ(T)|Oˆ|Ψ(T)(cid:105). (22) 1 volved,anadditionalshapingcomesfromthedependence ofthematrixelementsdjm ontheenergy. Thisise.g.the The mathematical problem is then reduced to the case for Fano line-shapes which may acquire a complex problem of maximizing a real-valued function G: Fano q factor [41]. G[c]=F[Ψ[E ],E ]. (23) c c The absorption of light is related to the average absorp- III. QUANTUM OPTIMAL CONTROL OF EXCITED STATE SPECTRA tion coefficient σ¯[Ec](E) [Eq. (14)]. The larger the ab- sorption coefficient at a certain energy, the more light is absorbed at this energy. In order to find a laser pulse to In this work we employ Quantum Optimal Control make a system, that is transparent in its ground state, Theory (QOCT) to optimize the response of a system in absorb as much light as possible, we therefore optimize the situation described in the previous section. QOCT the absorption coefficient in the visible by taking the in- is concerned with studying the optimal Hamiltonian (in tegral of σ¯[E ](E) over the respective energy range. We practice, a portion of the Hamiltonian, such as the tem- c employed two different control targets: poral profile of the coupling of an atom or molecule to a laser pulse) that induces a target system behaviour. (cid:90) Emax In the following, we present its specific application to GAτ [c]= dE σ¯τ[Ec](E), (24a) the problem of optimizing response functions of excited Emin sretadtuecse.dWtoeawsimllaalllsmoosdheolw, ithcoawn,bifetshoelvepdroabnleamlyticcaanllyb.e GBτ [c]=(cid:90) EmadxE σ¯τ[Ec](E)e(cid:16)−γN0−NN0T[E](cid:17), (24b) Emin Let us consider a quantum mechanical system gov- erned by the Schr¨odinger equation during the time in- whereσ¯τ[Ec](E)[inthefollowingwewillcallitjustσ¯(E)] terval [0,T]: is the average absorption coefficient of the system at a given time delay τ after the pump pulse E(t), and E min i∂Ψ(x,t)=Hˆ[E,t]Ψ(x,t), (20a) and Emax define the optimization region - the energy ∂t range, where the absorption is optimized. In the second Ψ(x,0)=Ψ (x), (20b) targetfunctionwehaveintroducedanexponentialfactor 0 5 that depends on N and N , the number of electrons in Our goal is to find a laser pulse that drives the system 0 T the system at the beginning and the end of the pump from state |Ψ(t=0)(cid:105)=|Φ (cid:105) into a target state |Ψ¯(cid:105) a pulse, respectively. The reason to introduce this factor is to avoid ionization, i.e. we wish to lead the system to |Ψ¯(cid:105)=α|Φa(cid:105)+β|Φb(cid:105)+γ|Φc(cid:105) (27) a state with the desired absorption properties, but keep- ing the ionization probability low. Keeping the ioniza- in a given time T – α,β, and γ are complex coefficients. tionlowisparticularlyimportantintheTDDFTcalcula- Sincethespectralpropertiesofthisstatecanbetheneas- tions, if performed with adiabatic functionals, since with ilyobtainedusingEqs.(17), (14), theproblemoffinding current state-of-the-art adiabatic functionals, ionization a pulse giving the desired optical properties translates of the system leads to unphysical shifts in the position to the one of maximizing the overlap |(cid:104)Ψ(T)|Ψ¯(cid:105)|2 while of the absorption peaks. The term exp(cid:16)−γN0−NT[E](cid:17) keepingthefunctionalformofthelaserfixed–i.e.chang- N0 ing only ω1,2, ϕ1,2, and ε1,2. If we choose ω1,2 resonant therefore inflicts a penalty, whose strength can be mod- with the transition frequencies ω ,ω and assume they ab ac ulated by γ, to pump pulses that produce strong ioniza- are sufficiently separated in energy we can apply the ro- tion. We implement ionization using absorbing bound- tating wave approximation and obtain the laser parame- ariesandthusthetotalnumberofelectronsis,ingeneral, tersasfunctionofα,β,γ as(seeAppendixBfordetails): notconservedduringtime. Inpractice,onecanalsocom- bine the two target functions: one may start optimiza- tfrioonmstuhseinpgreGvAτio,usanodptliamteurmc.ontinue with GBτ , restarting ε1 = T2 sina(racrccocso(|sα(||α)|)d|βab| (28a) 2 arccos(|α|) |γ| Once the target is defined one is left with the prob- ε = (28b) lem of choosing an optimization algorithm to find the 2 T sin(arccos(|α|)dac maximum (or maxima) of G. Two broad families can be and distinguished: gradient-free procedures, which only re- quire the computation of G given a control input E, and ϕ =ϕ −π+ω T (28c) β 1 ba gradient-based procedures, that also require the compu- ϕ =ϕ −π+ω T . (28d) tationofthegradientofGwithrespecttoE. QOCTpro- γ 2 ca vides an expression for the gradient that can be adapted We will come back to this example below in Sec. IVA. forthiscase(seeappendixAfordetails). Thisapproach, however is numerically unfeasible for the target covered inthispaper. Forthisreason,inoursimulations,weem- ployed the gradient-free Simplex-Downhill algorithm by IV. APPLICATIONS Nelder and Mead [42]. In principle, if the system can be reduced to a few- Any QOCT formulation is constructed on top of a level model, the optimal fields can be found analytically. given model for the physics of the process under study. To illustrate this approach, below we briefly illustrate a In this paper we study and optimize the absorption simple example of controlling the absorption properties spectra of atoms and molecules using either analyti- ofasingleHydrogenatom. Insteadofdirectlyoptimizing callysolvablemodelHamiltoniansorobtainingthespec- GA,B weherederivealaserthatdrivesthesystemintoa τ tra by using time-dependent density functional theory state with the wanted optical properties. Let us suppose (TDDFT) [11, 12] – the time-dependent counterpart of that the situation can be approximated by a three-level DFT [43]. HamiltonianHˆ witheigenstates|Φ (cid:105),|Φ (cid:105)and|Φ (cid:105)and a b c Based on the Runge-Gross theorem [10] TDDFT es- the corresponding eigenenergies ε , ε and ε . We define a b c tablishes a one-to-one correspondence between the time- the transition energies ω = ε −ε , ω = ε −ε and ab b a bc c b dependent density and the time-dependent external po- ω =ε −ε . The dipole coupling between the states is ac c a tential of a many-electron system. Together with the givenbyd andd (bothassumedtoberealnumbers), ab ac Kohn-Sham (KS) scheme [44] it allows us to recast the and we further consider the case where the coupling be- many-body time-dependent problem into a simpler one tween the states |Φ (cid:105) and |Φ (cid:105) is dipole forbidden. b c where the interacting electrons are replaced by a ficti- Thesystemispumpedbyalaserfieldcomposedoftwo tioussetofnon-interactingelectronswiththesametime- carrier frequencies ω1,2 of the form: dependent density. This system of non-interacting elec- trons can then be represented with a single Slater deter- E(t)=ε˜ (t)cos(ω t+ϕ )+ε˜ (t)cos(ω t+ϕ ), (25) 1 1 1 2 2 2 minant formed by a set of KS orbitals leading to great computational simplifications. with phases ϕ , amplitudes ε and envelope ε˜ (t) 1,2 1,2 1,2 In the following we will work with spin-compensated defined by systems of N electrons doubly occupying N/2 spatial (cid:18) (cid:19) t orbitals. The time evolution of these orbitals ϕ (i = ε˜ (t)=2ε sin2 π . (26) i 1,2 1,2 T 1,N/2), is governed by the time-dependent Kohn-Sham 6 equations ∂ 1 i ϕ (r,t)=− ∇2ϕ (r,t)+v [n](r,t)ϕ (r,t),(29) ∂t i 2 i KS i N/2 (cid:88) n(r,t)=2 |ϕ (r,t)|2, (30) i i=1 where v [n](r,t) is the KS potential. It is, in general, KS a functional of the density and is defined as v [n](r,t)=v (r)+v(r,t)+v [n](r,t)+v [n](r,t), KS 0 H xc (31) where v (r) represents the static (ionic) external poten- 0 tial, v(r,t) = E(t)·r is the coupling to the time depen- dent electric field E(t) in the dipole approximation (in thelengthgauge),v [n](r,t)=(cid:82)d3r(cid:48)n(r,t)/|r(cid:48)−r|isthe H Figure 1. Absorption coefficient σ¯(ω) of the state defined classicalelectrostaticHartreepotential,andv [n](r,t)is xc in Eq. (32) with ϕ = 0, (1/2)π, π and (3/2)π. The total theexchangeandcorrelationpotentialaccountingforthe spectrum (black line) is the sum of two phase-independent manyelectroneffects[11,45]. Inoursimulationstheions terms0.4σ¯ (redshaded)and0.6σ¯ (blueshaded)coming areclampedtotheirequilibriumpositions. Allnumerical fromtheex2cpiztedstatespectraofther3epszpectivestates,plusthe calculations were performed using the octopus code [46]. phase-dependent interference term σ¯INT(ω,ϕ) (green dashed line),whichisresponsibleforthechangeofthespectrumwith the delay time. A. One Electron Systems: the Hydrogen Atom ther purely Lorentzian (ϕ = 0, π) or purely Rayleigh In Section II we have discussed how the amplitudes (ϕ = 1/2π, 3/2π). The shaded areas indicate the and shapes of excited state absorption spectra depend weightedequilibriumcontributions,thedottedlineshows on the relative phases ϕi of the expansion coefficients γi. the interference terms, and the solid line the complete Here, we illustrate this effect in a Hydrogen atom, which spectrum. The energy range shown includes the transi- is initially pumped into the state tions from n = 2 to all higher states and from n = 3 to √ √ all higher states and to n = 2. Transitions to the ground |Ψ¯(cid:105)= 0.4|2pz(cid:105)+ 0.6eiϕ|3pz(cid:105). (32) state lie outside of this region. As can be learnt from Eq. (34), the interference terms Let us first vary ϕ to show the effect of the phase on the require the existence of states which are dipole-coupled finalspectrum. Forourpumpedstatethestationarypart toboth2p and3p . Thisisthecasefors-andd-orbitals. of the spectrum is composed of the weighted stationary z z This means, that e.g. for Hydrogen in a linear combina- state spectra coming from |2p (cid:105) and |3p (cid:105): z z tion of the states 2s and 4f, all the interference terms σ¯0(ω)=0.4σ¯ (ω)+0.6σ¯ (ω), (33) wouldvanishandthespectrumwouldbepurelythesum 2pz 3pz of the weighted equilibrium contributions. whereas the phase-dependent interference term is Let us take a closer look at the structure of the in- terference terms. We start with the interference term at σ¯INT(ω,ϕ)=0.4·0.6· 4πω (cid:104)(cid:88)(cid:90) d d ω23 = 0.069 Ha having contributions from terms with 3c 2p,m m,3p m=2s, m=3s and m=3d. All contributions have dif- m ferent prefactors with the ones coming from the 2s-state (cid:8) cosϕL(ω−ω3p,m)+sinϕR(ω−ω3p,m) having the opposite sign compared to the ones coming −cosϕL(ω−ω )+sinϕR(ω−ω )(cid:9) from 3s and 3d states. For the other peaks, the inter- 3p,m 3p,m (cid:90) ference terms are much smaller at the energies ω than (cid:88) 3n + d2p,mdm,3p their counterparts at ω2n (compare the purely blue to m the purely red peaks in Fig. 1). From Eq. (34), it is ap- (cid:8)cosϕL(ω−ω )−sinϕR(ω−ω ) parent that the amplitude of each interference term is 2p,m 2p,m the same for ω and ω with the same n. The differ- (cid:9)(cid:105) 2n 3n −cosϕL(ω−ω2p,m)−sinϕR(ω−ω2p,m) (34) ence comes purely from the factor 4πω – note that the 3c sign of the Rayleigh contributions is opposite in these Note the change in the sign of the Rayleigh terms in pairs of peaks. This variation of amplitude has the fol- both sums. Fig. 1 shows the different contributions and lowing consequences for the change of the overall spec- the complete spectrum for ϕ=0, (1/2)π, π and (3/2)π, trum: At ω the spectrum has positive contributions 3n which are the cases, where the interference term is ei- from σ¯ and contributions from the interference terms, 3pz 7 but since the interference terms are much smaller than sign of a pump pulse driving the system into a state σ¯ ,thespectrumchangesonlyslightlyfordifferentϕ’s. with specific optical properties. For this problem, we 3pz This is different for the peaks at energies ω . Here, the will use the three-levels model, and the analytical equa- 2n spectrum has positive, phase-independent contributions tions of control presented in Section III. The target from σ¯ (ω), but the contributions from the interfer- state will again be the one defined in Eq. (32) |Ψ¯(cid:105) = 2pz √ √ ence terms are much larger and dominate the spectrum 0.4|2p (cid:105)+ 0.6eiϕ|3p (cid:105) with a relative phase of ϕ=0: √z √ z leading to a strong dependence of the spectrum in this |Ψ¯(cid:105)= 0.4|2p (cid:105)+ 0.6|3p (cid:105). Thethreeactivestatesare z z energy range on the phase ϕ. For ϕ = 0 and ϕ = π, then|1s(cid:105),|2p (cid:105)and|3p (cid:105). Notethatsince|2p (cid:105)and|3p (cid:105) z z z z σ¯INT(ω,ϕ) only contains Lorentzian peaks and conse- havethesamesymmetry,theyaredecoupledinthedipole quently the whole spectrum only contains Lorentzians. approximation, and in consequence the system fits into Nevertheless,σ¯INT(ω,ϕ)changessignbetweenϕ=0and theframeworkdescribedinSectionIII.Wemaytherefore ϕ = π, switching the sign of all peaks at ω2n. This is a writedowntheshapeofacontrolpulse,assumingatotal demonstration of, how the manipulation of the internal pulse time of T =3200 a.u: phase ϕ can lead to a switch from gain (negative peaks) (cid:32) √ to loss (positive peaks) regime and vice versa. Finally, 2π 0.4 E(t)= cos(ω (t−T)+π) (35) for ϕ = (1/2)π and ϕ = (3/2)π, the interference spec- T d 1s→2p 1s→2p trum contains purely Rayleigh peaks. Together with the small contributions from the stationary-state contribu- √0.6 (cid:33) (cid:18) t (cid:19) + cos(ω (t−T)+π) sin2 π . tions, the final spectrum consists of slightly asymmetric d 1s→3p T 1s→3p Rayleighpeaks,againwithdifferentsignsforϕ=(1/2)π and ϕ = (3/2)π. One can therefore not only change We numerically solved the TSDE in order to check the peaksfromemissiontoabsorptionpeaks,butalsomanip- validity of the three-level approximation. To this end ulate their shape. The phases ϕ therefore play a critical wediscretizedtheequationsonasphericalgridofradius role in the spectral weights and the peaks of the photo- R=60 a.u., spacing of ∆x=0.435 a.u. and with 20 a.u. absorption spectrum. wide absorbing boundaries placed at the edges. The re- We now look at the variation of the spectrum with time, assuming that an initial, yet unknown pump laser created the state of Eq. (32) with ϕ = ϕ = 0 at t = 32 T and we probe the system at different delay times τ. Figure2showsthecorrespondingtime-resolvedspectrum s σ¯(ω,τ) of |Ψ¯(cid:105). Since the eigenenergies of |2p (cid:105) and |3p (cid:105) n0.8 z z o are different, the phase Θ (τ) in Eq. (9) evolves with i 32 t the frequency ω32. At τ = 0, τ = 2ωπ32, τ = ωπ32 and ula0.4 τ = 3π , the spectra of Fig. 1 are reproduced. One sees p the s2tωro32ng changes of σ¯ in the energy range of the peaks Po 0 ω , whilethepeaksω remainalmostunchanged. The 2n 3n spectrum is periodic with T = 2π ≈ 91 a.u.. ω32 0 1000 2000 3000 t [a.u.] 50 30 Figure3. Time-evolutionofthepopulationsofthe1s−,2p − 10 z and 3p − state. Dashed lines show the analytic model, solid z -10 lines the numerical results. The total pump-laser (upper panel,black)hastwocarrier-frequencies,oneresonanttothe -30 0 transition |Ψ (cid:105) → |Ψ (cid:105) (green), the other resonant to the 1s 2p 20 0.14 transition |Ψ (cid:105) → |Ψ (cid:105) (blue). The lower panel shows the 40 0.08 1s 3p 60 0.06 phase difference ϕ3p−ϕ2p. 80 0.02 100 sults are collected in Figure 3 where we show the time- evolution of the populations |a(t)|2, |b(t)|2 and |c(t)|2 of F√igure 2. √Time-resolved spectrum of the initial state thestates|1s(cid:105), |2pz(cid:105)and|3pz(cid:105)respectively. Thenumeri- 0.4|2p (cid:105)+ 0.6|3p (cid:105) of Hydrogen. Because the phases of calvalues(solidlines)followcloselytheonescorrespond- z z the |2p(cid:105) and |3p(cid:105) states evolve with different velocities, the ing to the model (B10) (dashed lines) except for a small spectral weights of each of the peaks changes with time, superimposed oscillatory behavior. A frequency analysis leading to a time-dependent spectrum with a periodicity of of the additional oscillations shows, that they are due to T = 2π ≈91a.u.. ω32 the components neglected in the rotating wave approx- imation. The small deviation in the final populations We now move on to the control problem – i.e. the de- from the analytic prediction comes from the excitation 8 into the 3d-states (not shown). The coupling to these lengths λ = 800 nm and λ = 1450 nm and their first orbitals was neglected in the three-level approximation. nine odd harmonics as shown in Figure 4(a). Optimiza- This population transfer to the 3d-states nonetheless is tion is achieved transferring population from the ground less than 4% , and we achieve a transfer into the target state (at (cid:15) = −2.238 a.u.) into the first excited state 0 wave function |Ψ¯(cid:105) of 96%. Furthermore the transfer is (at (cid:15) = −1.705) with the help of the 9th harmonic of 1 obtained precisely with the desired relative phase ϕ = 0 λ = 1450 nm at ω = 0.534 a.u.. Due to this popula- P13 as reported in the bottom panel of Figure 3. tion transfer, the peak at ω = 0.533 a.u. turns from 0→1 positive to negative and peaks coming from the first ex- cited state (ω = 0.076 a.u. and ω = 0.159 a.u.) 1→2 1→4 B. More Than one Electron: Results Based on arise in the excited-state spectrum, where the peak at TDDFT ω1→2 is located in the visible part of the energy range. At the same time, population is transfered into the sec- ond excited state ((cid:15) = −1.629 a.u.), leading to e.g. the We here turn to systems with more than one electron, 2 peaks at ω =0.062 a.u. and ω =0.103 a.u.. This and investigate the possibility to drive the absorption of 2→3 2→5 interpretation is confirmed by the population analysis in atoms and molecules into the visible using a laser pulse Figure 4(d) where |(cid:104)Ψ(t)|Φ (cid:105)|2 is plotted over time. In optimized with the gradient-free optimization algorithm i particular it is apparent that at the end of the pump presented in Section III in combination with TDDFT. only ≈ 8% of the electrons remain in the ground state, whereas the rest has been transfered into higher lying states thus explaining the appearance of the new peaks 1. Helium in the spectrum. To complete the picture in Figure 4(c) we show the evolution of the control function GA with As a first example we study the one-dimensional soft- the number of iterations. As can be seen, GA shows a Coulomb Helium atom. This model is defined by the steady increase during the optimization. Hamiltonian: H(t)=T +V (t)+V , (36) ext ee where 2 2 V (t)=− − +E(t)(x +x ) (37) ext (cid:112) (cid:112) 1 2 1+x2 1+x2 1 2 is the external potential with softened Coulomb inter- action and the dipolar coupling to the external time- dependent field E(t) The electron-electron interaction is also described by a soft-Coulomb function V = ee √ 1 . For the optimization we solve the equa- 1+(x1−x2)2 tionsdiscretizedonaregulargridofsizeL=100a.u.and spacing ∆x=0.2 a.u. with 20 a.u. absorbing boundaries Figure 4. Optimization of the absorption of one-dimensional at the borders of the simulation box. Results obtained Helium: TDSEvs.TDEXX.(a)Powerspectrumoftheinitial withtheoptimizedpulsewerefurtherconvergedinabox and optimized laser pulses: The green line shows the initial of size L = 200 a.u. with 70 a.u. absorbing boundaries. laser used for the TDSE; the green shaded area shows the Time was discretized with a time step of ∆t=0.025 a.u. initial laser used for the TDEXX optimization (these two for a maximum propagation time of 1250 a.u. during op- only differ by the indicated peaks at ω = 0.534 a.u. and timization and 2250 a.u. for convergence. The duration ω=0.549a.u.);theblueandredlinesaretheoptimalpulses of the pump pulse was chosen to be TP = 800 a.u. and obtained when using TDSE and TDEXX, respectively. (b) thedelaybetweenpumpandprobewassettoτ =50a.u. Ground state (dashed line) and excited state (shaded) spec- for all the calculations. Finally the target region for op- tra of optimized one-dimensional Helium, in blue and red for timization was chosen between 0.06 a.u. and 0.23 a.u. theTDSEandTDEXXcases,respectively. Theexcitedstate (≈ 200 nm and 800 nm). We carried out optimizations transitions of the exact calculations are indicated. (c) The controlfunctionGA asafunctionofthenumberofiterations, attwodifferenttheorylevels: exact(TDSE)andTDDFT alsoinblueandredfortheTDSEandTDEXX,respectively. with the adiabatic EXX functional (TDEXX) [47]. (d) The populations |(cid:104)Ψ(t)|Ψ (cid:105)|2 of the exact time propaga- Let us first focus on the optimization obtained by i tion under the influence of the optimized pump pulse for the solving the exact TDSE as illustrated in Figure 4(b) (red) ground state, (green) first excited state, (blue) second where the ground state spectra are compared to the excitedstate,(pink)thirdexcitedstateand(turquoise)fourth spectra of the systems excited by the optimized pump- excited state. pulses E (t) obtained after 100 iterations. In the exact 100 case, the search space was constructed from two wave For the TDEXX case we adapted the search space 9 laser pulse resonant with the excitation energy from the ground state into the first excited state. After the pulse thesystemsareinasuperpositionofthesetwostatesand the spectra should contain time-dependent interference terms, which oscillate with the period time T = 2π , (cid:15)1−(cid:15)0 which is T = 11.76 a.u. for TDSE and T = 11.26 a.u. for TDEXX. However, on the scale of Figure 5(a) the TDSE spectrum hardly presents any oscillation. There- fore, in Figure 5(c) we report cuts at ω = 0.2, 0.4, 0.6 n and 0.8 a.u. From the figure it is apparent that, albeit with different phases, each cut presents oscillations with the expected period of T =11.76 a.u.. The TDEXX cal- culations, in Figure 5(b) and (d), present a different pic- ture. Firstofall,theamplitudeoftheoscillationsismuch largerthanintheexactcaseand,second,theoscillations are two times faster than expected. We conclude, that Figure 5. (Top) Transient Absorption Spectrum of Helium aftertheexcitationwitha45cyclesin2laserpulseofintensity theTDEXXdescriptionseemstohaveasimilarstructure I =5.26·1011Wcm−2withacarrierfrequencyresonanttothe totheexactcase, inthesense, thattheenergydifference excitation energy from the ground to the first excited state oftheinvolvedstatesisreflectedintheperiodicityofthe for exact (a) and adiabatic EXX (b) with ωexact =0.534 a.u. oscillations of the spectrum. Nonetheless, there are ma- and ωEXX = 0.549 a.u.. Time-evolution of the absorption jor differences in the behaviour, which is reflected in the cross-section at selected energies ωn = 0.2 (red), 0.4 (blue), factor of two in the periodicity. 0.6 (purple) 0.8 (turquoise) a.u. for exact (c) and adiabatic EXX (d). In the exact case the curve at 0.6 a.u. is offset by -0.1forclarity. Inallcases,thetimeintervalT =2π/((cid:15) −(cid:15) ) 1 0 is shown. 2. Methane Dication Finally, we apply our scheme to a poly-atomic by replacing the laser component at the carrier fre- molecule: doubly-ionizedMethane,CH+2. Thegoalhere 4 quency ω = 0.534 a.u. (in resonance with the first ex- is to design a laser capable to turn this molecule, trans- citation in the TDSE case) by a laser component with parent in nature, visible. To this end we used the same ω = 0.549 a.u., which is its TDEXX equivalent. In our strategy as we did before for Helium, namely we opti- experience failing to meet this requirement resulted in mize the laser on a small simulation box and then con- poor optimizations. The resulting optimization is shown verge the results with the optimized laser on a larger in Figure 4(b). The results follow a trend similar to box. Duringtheoptimizationroutinethesimulationbox the exact case. However the TDEXX optimization is has a radius of R = 15 a.u., including 5 a.u. absorbing smaller and the 1 → 2 peak present in TDSE seems to boundaries while the results are converged in a box of be missing. The difference between TDEXX and TDSE R = 30 a.u. with 15 a.u. absorbing boundaries. We dis- can be tracked down to a known problem of the adia- cretizetheTDDFTequationsonathree-dimensionalgrid batic approximation in TDDFT. In particular, the lack with a spacing of ∆x=0.3 a.u.. The reason for this box of memory in the adiabatic approximation, causing a choice is the fact that the computational costs of three- spurious time dependence of the exchange potential, is dimensional calculations scale with the third power of responsible for the poor population transfer and the ex- the simulation box radius. The maximum propagation cess of asymmetric peaks in the spectrum [14, 48]. This time is 850 a.u. during the optimization and 1600 a.u. problem is further amplified by the ionization of the sys- for convergence. In all cases, the pump duration was tem, which results in an unphysical shift of the peaks 600 a.u., the time step was ∆t=0.04 a.u. and the delay tohigherenergies(comparethegroundstateandtheex- was τ =0 a.u.. citedstatespectruminFigure4). Theseeffects,however, Theoptimizationregionwaschosenastheintervalbe- strongly depend on the fraction of the total density that tween 0.057 a.u. and 0.139 a.u. (328 and 750 nm) and in getsdrivenoutofequilibriumandthereforebecomemore ordertodiscouragethealgorithmfromexcitingtoomany dominant with decreasing size of the system – with He- electrons into the continuum, we used the target func- lium being the worst case. In large molecules with many tionalGB (24b),whichincludesanexponential“penalty” τ electrons we expect the error to be greatly reduced (as for ionization. has been empirically shown in studies of light induced To obtain a good description of states close to the charge transfer in organic photovoltaic blends [49, 50]). ionization threshold, we employed the average density A different perspective on the same problem can be self-interaction corrected (ADSIC) LDA functional [51], obtained by comparing the time evolution of an excited which is asymptotically correct. state spectrum in TDSE and TDEXX as shown in Fig- The inclusion of resonant frequencies in the search ure 5. The systems were excited by a 45 cycle sin2 space is a good practice that facilitates transitions be- 10 tween eigenstates and enables the optimization algo- cited states ((cid:15) to (cid:15) ). The idea is that the system could 3 8 rithm to populate excited eigenstates. These molecular be excited into the ionization threshold, and then relax excitation frequencies can easily be obtained from the into one of the target states. The laser frequencies, that ground-state spectrum reported in Figure 6. By popu- wereincludedinthesearchspacesandthecorresponding lating the correct eigenstates, the system might absorb resonances are summarized in Appendix C. in the visible region: consider two eigenstates with en- ergies (cid:15) and (cid:15) , that differ by an energy in the visible: h T 0.057 a.u. ≤ (cid:15) −(cid:15) ≤ 0.139 a.u.. By exciting the sys- h T tem into the lower “target” state (cid:15) one might obtain T transition peaks in the visible, due to the transition to thehigherone. Note,however,thatthisfactisnotguar- anteed since the transition might be dipole forbidden. We cannot rule out this possibility since our ground- state linear-response TDDFT calculation does not pro- videthisinformation. Thegroundstatespectrumshows, Figure 7. Ground state (black line) and excited state (shaded) spectra of doubly-ionized Methane (CH2+) for two 4 different pump pulses EI (red) and EII (blue). The grey shadedboxmarkstheoptimizationarea,theblueshadedarea the visible region of the spectrum. TheoptimizedspectraareshowninFigure7. Itcanbe seen that both search spaces include optimal lasers that cause the molecule to loose its transparency and absorb in the visible. The achieved opacity can be quantified in terms of the control function GA (24b) being the inte- Figure 6. Ground state spectrum of doubly-ionized Methane gral over the absorption spectrum in the visible range of CH2+. Peaks are numbered for later reference (as discussed the spectrum. Comparing the opacity achieved in search 4 in the text). The shaded grey area marks the optimization space I (GA = 0.017) with the one achieved in search I range. space II (GA = 0.020), we conclude that search space II II is better suited for the pursued optimization. Thus, that the first possible target state is (cid:15) . The energy dif- including energy levels at the ionization threshold in the 3 ference between (cid:15) = 0.690 a.u. and (cid:15) = 0.816 a.u. is search space might be a useful strategy in further opti- 3 4 ω = 0.126 a.u. and lies – with 362 nm – at the red mizations. 3→4 end of the visible spectrum. Also (cid:15) provides a tran- 4 sition in the visible range – into (cid:15) = 0.938 a.u. with 5 ω = 0.122 a.u. = 373 nm. Starting from (cid:15) , the V. CONCLUSIONS 4→5 5 states have even more than one transition in the visi- ble. We must therefore choose a frequency search space, In this work, we assessed the possibility of using tai- thatallowstheconstructionofapumppulse,thatexcites lored pumps in order to enhance some given features of electrons from the ground state into (cid:15)3 and higher lying the probe absorption – for example, the absorption in states either directly or by successive excitations. the visible range of otherwise transparent samples. We Here we present results for two possible search spaces. firstdetailedatheoreticalanalysisofthenon-equilibrium Thefirstsearchspaceincludesfrequenciesthatareeither response function in this context, aided by one simple resonanttothegroundstateexcitationenergies(cid:15) , orto numerical model of the Hydrogen atom. Then, we inves- n excitedstateexcitations(cid:15) −(cid:15) . Toavoidionization,all tigatedthefeasibilityofusingTDDFTtheoryasameans m n carrier frequencies are smaller than (cid:15) = 1.0 a.u.. The toimplement,theoretically,thisabsorption-optimization 7 secondfrequencysearchspacewasdesignedusingtheion- idea, for more complex atoms or molecules. ization potential I of the system (which is equal to mi- Thetheoreticalanalysisoftheresponsefunctioncanbe P nustheenergy(cid:15) ofthehighestoccupiedKSorbitalϕ ) done by writing it in a generalized form of the Lehmann H H and the energy differences I −(cid:15) . One frequency is I representation, validforsystemsthathavebeenpumped P n P itself,whileallothersareresonantwithrelaxationsbring- out of equilibrium by a first pulse, and whose response ingdownstatesatthationizationthresholdtoboundex- to a probe pulse (in our case, assumed non-overlapping