Table Of ContentTailored 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
