Excited states dynamics in time-dependent density functional theory: high-field molecular dissociation and harmonic generation.
Alberto Castro$^{1,2,3}$
M. A. L. Marques$^2$
Julio A. Alonso$^1$
George F. Bertsch$^3$ and Angel Rubio$^2$
We present a theoretical description of femtosecond laser induced dynamics of the hydrogen molecule and of singly ionised sodium dimers, based on a real-space, real-time, implementation of time-dependent density functional theory (TDDFT). High harmonic generation, Coulomb explosion and laser induced photo-dissociation are observed. The scheme also describes non-adiabatic effects, such as the appearance of even harmonics for homopolar but isotopically asymmetric dimers, even if the ions were treated classically. This TDDFT-based method is reliable, scalable, and extensible to other phenomena such as photoisomerization, molecular transport and chemical reactivity.
dimers.tex
\documentclass[pra,twocolumn,showpacs,footinbib]{revtex4}
\usepackage{graphicx,rotating}

%%%%%%%%%%%%%%%% Useful Macros
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\bea}{\begin{eqnarray}}
\newcommand{\eea}{\end{eqnarray}}
\newcommand{\ba}{\begin{array}}
\newcommand{\ea}{\end{array}}
\newcommand{\mint}[1]{\int\! d^{3} #1 \, }
\newcommand{\mdint}[2]{\mint{#1}\!\!\mint{#2}}
\newcommand{\mkint}[1]{\int\!\frac{d^3 #1}{(2\pi)^3}}
\newcommand{\mdkint}[2]{\mkint{#1}\!\!\mkint{#2}}
\newcommand{\sgn}{{\rm sgn}}
\newcommand{\Tr}{{\rm Tr}}
\newcommand{\da}{{\dagger}}
\newcommand{\cD}{\hat {\cal D}}
\newcommand{\unitmatrix}{{\mathbbm{1}}}
\newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\fd}[2]{\frac{\delta #1}{\delta #2}}
\newcommand{\T}{^{\rm T}}
% The bold faces
\newcommand{\bA}{{\bf A}}
\newcommand{\bB}{{\bf B}}
\newcommand{\bk}{{\bf k}}
\newcommand{\bp}{{\bf p}}
\newcommand{\bq}{{\bf q}}
\newcommand{\bQ}{{\bf Q}}
\newcommand{\br}{{\bf r}}
\newcommand{\bR}{{\bf R}}
\newcommand{\bs}{{\bf s}}
\newcommand{\bu}{{\bf u}}
\newcommand{\bv}{{\bf v}}
\newcommand{\bx}{{\bf x}}

\newcommand{\bpi}{{\bf \pi}}
\newcommand{\bg}{{\bf \hat \gamma}}
%%%%%%%%%%%% End Macros

\begin{document}
\title{Excited states dynamics in time-dependent density functional theory:
      high-field molecular dissociation and harmonic generation.}

\author{Alberto Castro$^{1,2,3}$, M. A. L. Marques$^2$, Julio A. Alonso$^1$, George F. Bertsch$^3$ 
        and Angel Rubio$^2$}
\affiliation{
$^1$Departamento de F\'{\i}sica Te\'{o}rica, Universidad de Valladolid, 47011 Valladolid (Spain)\\
$^2$Departamento de F\'{\i}sica de Materiales, Facultad de Qu\'{\i}micas, Universidad del Pa\'{\i}s Vasco;\\
    Centro Mixto CSIC-UPV/EHU and Donostia International Physics Centre (DIPC), 20080 San Sebasti\'{a}n (Spain)\\
$^3$Physics Department and Institute for Nuclear Theory, University of Washington, Seattle WA 98195 (USA)}

%\date{}

\def\be{\begin{equation}}
\def\ee{\end{equation}}

\begin{abstract}
We present a theoretical description of femtosecond laser induced dynamics
of the hydrogen molecule and of singly ionised sodium dimers, 
based on a real-space, real-time,
implementation of time-dependent density functional theory (TDDFT). 
High harmonic generation, Coulomb explosion and 
laser induced photo-dissociation are observed. 
The scheme also describes non-adiabatic effects, such as the appearance of even harmonics
for homopolar but isotopically asymmetric dimers,
even if the ions were treated classically.
This TDDFT-based method is reliable, scalable, and extensible
to other phenomena such as photoisomerization,
molecular transport and chemical reactivity.
\end{abstract}
\pacs{33.80.Gj, 33.80.Wz, 42.50.Hz}

\maketitle

It is now possible to study electron and molecular dynamics in 
real time 
using various experimental techniques employing
intense ultra-short laser sources \cite{pBrabec00_yotros}.
%Still needed is a practical 
%computational framework to
%describe the excited-state electron-ion dynamics theoretically. 
Some examples of such investigations
%of processes
%we would like to describe 
include X-ray photoelectron spectroscopy of molecules \cite{pLeone00_yotros},
pump-probe ionisation measurements \cite{faraday_discussions}, production of high harmonics
as a source of soft X-rays \cite{x-ray-citations}, the measurement of 
electron-phonon interactions
in thin films \cite{pProbst97}, and the estimation of the onset of Coulomb screening 
\cite{pHuber01}.
%A feasible computational framework, with solid theoretical foundations, is 
%nowadays clearly necessary to help understanding the wealth of new phenonema
%that have emerged in the last years since the advent of intense, ultra-short
%laser sources \cite{pBrabec00_yotros}.
%With the new femto- and sub-femtosecond lasers equipping
%modern laboratories, it is now common, for example, to probe the electron dynamics
%in many-body systems, or to study non-linear optical processes in
%condensed matter.
%Applications include X-ray photoelectron spectroscopy in molecular 
%systems \cite{pLeone00_yotros}, 
%above threshold ionisation, production of high harmonics as a
%source of ultra-violet or soft X-ray lasers, the study of the electron-phonon
%interactions on the surface of organic thin films \cite{pProbst97} 
%and the estimation of the onset of Coulomb screening \cite{pHuber01}.
%Also there
%have been extensive studies of above threshold ionisation and
%alignment of molecules in laser fields (see, e.g. Ref. \cite{pMacholm01}, and papers cited therein).
A technologically important and very active field of research is the application of ultra-short laser
pulses to induce, control and monitor
chemical reactions \cite{pDion96_yotros,pGerber01,pYamanouchi02}. 
Whenever the intensity of the laser field is comparable
to the molecular electronic fields, perturbative 
expansions break down and new processes appear, 
which are not fully understood from a microscopical point of view \cite{new_processes}.
A practical and accurate computational framework to descibe excited-state electron-ion dynamics
is therefore still needed.


Not surprisingly, the smallest systems have attracted particular attention 
from both experimentalists and theoreticians, as a bench-horse to improve
our understanding of electron dynamics at the femtosecond scale \cite{pKnight97, atoms_and_dimers_studies}.
%Due to their computational
%simplicity, atoms \cite{pKnight97} and dimers have been extensively studied
%\cite{atoms_and_dimers_studies}.
However, the methods used in these calculations
can not be easily extended to larger and more realistic systems. The exact quantum
mechanical solution of a 3D system of more than three particles
is certainly not feasible with state-of-the-art computers. 
1D models are much easier to handle, but they can not really be used as predictive
tools for problems involving the interaction of lasers with large clusters
or solid-state systems of technological relevance.

To tackle such a problem, time-dependent density functional 
theory %\cite{pRungeGross84,bRunge96,pAngel02}
(TDDFT) \cite{tddft} appears as a valuable tool.
% making a compromise between 
%computational ease and accuracy.
Even with the simplest approximation to the exchange-correlation potential, the
adiabatic local density approximation (ALDA), one obtains a very good compromise
between computational ease and accuracy \cite{ALDA_problems}.
TDDFT can certainly be applied to large systems
in non perturbative regimes, while providing a consistent treatment
of electron correlation. It has been well tested in the study
of electron excitations, like the optical absorption spectra in
the linear regime \cite{pYabana96, review}.
%In its simplest form, the adiabatic local density approximation (ALDA),
%it has some well known deficiencies, which can nevertheless
%be corrected by the swarm of exchange-correlation functionals now available \cite{marques}.
Although almost all applications of TDDFT in the field of laser physics
have only involved electronic dynamics, recent attempts 
have also been made at describing the coupled nuclear and electronic 
motion in laser fields \cite{pCalvayrac00_yotros}, accounting for the nuclear motion classically.
A full quantum mechanical treatment of the system
could in principle be done within
a multi-component TDDFT,
although it has not been tried for
more than three particles \cite{pKreibich01-2}.
However, since many vibrational quanta are coherently excited, there is good motivation
for the classical treatment of the nuclei.
The purpose of this work is to illustrate a general method
to study many-electron systems subject to strong laser fields.
It is based on the quantum mechanical propagation of the electronic
wave packet -- described within TDDFT --
combined with classical motion of the nuclei.
As the laser field populates the excited Born-Oppenheimer
surfaces, this scheme includes diabatic effects, while maintaining
a good scaling with the size of the system.
%This will allow us to study complex processes such
%as photoisomerization or photochemical reactivity, in systems
%of tens of atoms.
As an illustration
we focused on one and two electron dimers, namely
Na$_2^+$ and the hydrogen molecule.
%The former is a good
%test, since it has been exhaustively studied using
%a diverseness of approaches \cite{pMagnier99,pMagnier96,pMacholm96,pGerber01}.
%In particular, a recent experiment \cite{pGerber01} focused on
%the photofragmentation of Na$_2^+$ in intense femtosecond laser fields.
%Using a pump-probe technique, the authors discovered that Na$_2^+$
%dissociated in four different channels, ranging from simple field ionisation
%followed by Coulomb explosion, to photodissociation on light-induced 
%potentials. For the H$_2$ molecule, Kreibich and coworkers
%recently reported 1D exact calculations beyond the Born-Oppenheimer 
%approximation \cite{pKreibich01}.
%Changing the isotopic composition to HD, they observed the generation of the
%parity forbidden even peaks in the harmonic spectrum, concluding that
%non-adiabatic effects can alter crucially the emission spectrum.

The equations of motion may be derived from the Lagrangian:
\def\a{\alpha}
\begin{eqnarray}
  \nonumber
  {\cal L} = 
  \sum_\a \left[ {1\over 2} m_\a  \dot{\vec{R}}_\a^2
  + Z_\a \;\vec R_\a \cdot \vec {\cal E} (t) \right]
  - \sum_{\a<{\a'}} {Z_\a Z_{\a'} \over | \vec{R}_\a - \vec{R}_{\a'}|}
\\
  \sum_i \langle \phi_i |
  i\;\hbar{\partial \over \partial t}
  -e \; \vec x \cdot \vec {\cal E} (t)
  |\phi_i\rangle
  - E_{DFT} (\{\phi\}, \{\vec R\})
  \,,
\end{eqnarray}
where $E_{DFT}$ is the usual Kohn-Sham density functional, depending on the
electron orbitals $\{\phi\}$ and the nuclear coordinates $\{\vec R\}$, and $\vec{\cal E}(t)$
is the time-dependent electric field from the laser pulse. Variation of the
Lagrangian then yields Newton's equations for the nuclear coordinates,
\begin{eqnarray}
  \nonumber
  m_\a {d^2 \vec{R}_\a \over d t^2} & = & - \vec{\nabla}_{R_\a} E_{DFT}
  \\ \label{eq:forces}
  & & -\sum_{\a'}{(\vec {R}_\a
  -\vec R_{\a'})Z_\a Z_{\a'}\over |\vec{R}_\a -\vec{R}_{\a'}|^3} + Z_\a \; \vec {\cal E} (t),  
\end{eqnarray}
and the usual TDDFT equations for the orbital variables \cite{tddft}.
%Changing the isotopic composition to HD, they observed the generation of the
%parity forbidden even peaks in the harmonic spectrum, concluding that
%non-adiabatic effects can alter crucially the emission spectrum. 
We solved these equations in real time, following the method of
Yabana and Bertsch \cite{pYabana96},  using a real-space grid
representation of the orbitals \cite{pFlocard78,pChelikowsky94_yotros}. This scheme has the advantage
that the Kohn-Sham Hamiltonian is a very sparse matrix.  The forces in Eq. (\ref{eq:forces})
are calculated with the help of a generalised Hellmann-Feynman theorem,
\be
 - \vec{\nabla}_{\vec{R}_\a} E_{DFT}= -\sum_i \langle \phi_i |\vec{\nabla}_\a H_{KS}|\phi_i\rangle.
\ee
where $H_{KS} = \delta E_{DFT} / \delta n({\bf r},t)$ is the Kohn-Sham
Hamiltonian.
For numerical reasons we represent the electron-ion interaction by
norm-conserving non-local Troullier-Martins pseudopotentials 
\cite{pTroullierMartins91}. 
% ATTENTION:
%A possible vector potential could be included
%in the model, although the non-local parts of the pseudopotentials must
%be handled with care to ensure gauge invariance \cite{pBertsch00}. 
%A useful numerical check on the computation is
%energy conservation, which requires that the excitation energy of the
%molecular at the end of the time evolutions equals the work done by
%the external field. {\bf Can we say this?}




%Technically, we have implemented our TDDFT scheme
%within a real-space, real-time framework. This
%technic has been extensively applied to the calculation of optical
%response in clusters \cite{pChelikowsky94,pYabana96,pCastro02} and 
%solids \cite{pBertsch02}. The Kohn-Sham states (or the exact single-particle
%states in the case of Na$_2^+$) were propagated with a Krank-Nicholson
%like scheme \cite{pFlocard78}, whereas
%the time-dependent forces on the classical nuclei were calculated from Ehrenfest theorem:
%\be
%  \vec{F}_\alpha(t) = \left<\Psi(t)\right|\vec{\nabla}_{R_\alpha}\hat H(t)
% \left|\Psi(t)\right>\,.
%\ee
%$\left|\Psi(t)\right>$ denotes the Kohn-Sham Slater determinant
%at time $t$, $\hat H$ is the Kohn-Sham Hamiltonian of the system
%-- including the time-dependent external field --
%and $\vec{R}_\alpha$ represents the position of ion $\alpha$. 
%The ion-electron
%interaction was modelled by a norm-conserving non-local 
%Troullier-Martins pseudopotential \cite{pTroullierMartins91} in the case
%of sodium, and by a soft Coulomb in the case of hydrogen. 
%Our tipical setups had grid spacings of 0.3\,\AA\ for Na$_2^+$
%and 0.14\,\AA\ for the hydrogen molecules, and the mesh points were
%contained inside a sphere of at least 10\AA\ radius. With these
%parameters the eigenvalues were converged to better than $0.1$\,eV.
%The time-steps for the time propagation were 0.005\,$\hbar\,eV^{-1}$ 
%for Na$_2^+$ and 0.002\,$\hbar\,eV^{-1}$ for hydrogen.
%In all
%simulations we used as an envelope for the laser electric field a
%half-cycle of a cosine. 


We have studied two different classes of time-dependent problems:
photofragmentation and high harmonic generation.
We now discuss the first case, the Na$_2^+$ dimer in 
a femtosecond laser field.  
It is a good test, since it has been exhaustively studied using
a diverseness of approaches \cite{pMagnier99,pMagnier96_yotras,pGerber01}.
In particular, a recent experiment \cite{pGerber01} focused on
the photofragmentation of Na$_2^+$ in intense femtosecond laser fields.
Using a pump-probe technique, the authors discovered that Na$_2^+$
dissociated in four different channels, ranging from simple field ionisation
followed by Coulomb explosion, to photodissociation on light-induced 
potentials.
For these calculations, we used an uniform grid
spacing of 0.3~\AA, and the system was confined to a sphere of
radius 10~\AA. 
In these one-electron calculations, we omitted the core-valence exchange-correlation.
%Since there is only a single electron, we omit the
%exchange-correlation term in the energy functional $E_{DFT}$. 
In this case the total electronic energy is given by the Kohn-Sham eigenvalues,
and they can be used to compute the adiabatic potential energy surfaces
shown in Fig. \ref{fig:pes}. 
The two lowest single-photon transitions from the $1^2\Sigma_g^+$ ground state
are at 2.5\,eV (to the $1^2\Sigma_u^+$ state) and at 3.2\,eV (to the $1^2\Pi_u$ state).
The latter is achieved by a laser polarized perpendicularly to the internuclear axis.
These energies accord well with the observed
single-photon transitions \cite{bonacic}. For the time-dependent calculations,
we start with the dimer in its ground state, which is then propagated with
a modified Krank-Nicholson scheme \cite{pFlocard78}.
The time step for the time integration was 
0.005 $\hbar$ eV$^{-1} \approx 0.003$ fs.  A simple check on
the implementation of the time evolution operator consists of calculating
the linear photo-absorption spectrum, using a weak 
$\delta$-function external field, as in Ref. \cite{pYabana96}.  Almost all the
spectral weight is concentrated in two peaks (see inset in Fig.~\ref{fig:pes}),
which are at energies
corresponding exactly to the vertical transitions between energy
surfaces.  We note that this exact correspondance is only obtained
for one-electron systems: in general the TDDFT spectra will have
shifts from the energy surfaces determined by the Kohn-Sham eigenvalues.
\begin{figure}
\includegraphics[scale=0.33]{pes.ps}
\caption{\label{fig:pes}
Adiabatic energy surfaces of the Na$_2^{+}$ dimer,
as obtained by our three-dimensional real-space code.
Similar results can be found in Ref. \cite{pMagnier_Oise96}.
In the inset, the computed photoabsorption cross-section of the
same molecule.}
\end{figure}
%\begin{figure}
%\includegraphics[scale=0.33]{plot1.eps}
%\caption{
%Adiabatic energy surfaces of the Na$_2^{+}$ dimer,
%as obtained by our three-dimensional real-space code.
%Similar results can be found in Ref. \cite{pMagnier_Oise96}.
%In the inset, photoabsorption cross-section of the
%same molecule.}
%\end{figure}
\begin{figure*}[t]
\includegraphics[scale=0.43]{dis.eps}
\caption{\label{fig:resonant_dissociation}
Evolution of internuclear distance (top panel) and
electronic charge in simulation region (bottom panel) for
the Na$_2^+$ molecule. The dimers are excited with
laser pulses of 2.5, 3.2 and
1.57~eV in columns (a), (b) and (c) respectively.
}
\end{figure*}
\begin{figure}[b]
\includegraphics[scale=0.32]{even.eps}
\caption{\label{fig:harmonics.ps}
Harmonic spectra of HD (left panels) and H$_2$ (right panels). 
The nuclear masses used in the calculation are $m'=\alpha m$, being
$m$ the real mass. In this way, top
plots were made using for the nuclear masses their real values
whereas bottom plots were made using a hundredth of their real values.
%The laser pulse frequency was 770~nm, and its total length 30 cycles.
%The shape was a linear ramp until the 10th cycle, and constant
%thereafter.
}
\end{figure}

Next we examine the evolution of the dimer under high-field excitation.
We consider external fields of the form:
\begin{equation}
 \vec {\cal E}(t) = \left( \frac{8\pi}{c} I_0\right)^{(1/2)} 
 \sin \left( \pi \; \frac{t}{\tau} \right) \sin(\omega t) {\hat e},
 \;\;\;
 0 < t <\tau,  
\end{equation}
where $I_0$ is the maximum intensity of the pulse,
${\hat e}$ is the polarization vector,
and $\tau$ is the 
pulse duration, taken as $\tau = 80 $ fs.  As a first case, we
examine the effect of excitation at the lower resonant frequency,
$\omega = 2.5$ eV.  In Fig. 2(a), we present a series of
runs at different intensities, ranging from weak ($10^{10}$ W/cm$^2$)
to moderate ($2.1\times 10^{12}$ W/cm$^2$).  Since the $1^2\Sigma_u^+$
surface is anti-bonding, excitation at this resonant frequency should 
lead to dissociation, even at moderate intensities.  This is indeed confirmed
by our calculations.
The upper panel depicts the internuclear separation 
of the dimer, which exhibits an acceleration during the laser pulse and
a nearly constant velocity expansion thereafter.  Clearly, the dimer
dissociates at all field levels that we applied.  To examine the ionisation
of the dimer,
we assumed that any density reaching the edges of the simulation box corresponds
to unbound electrons. By absorbing this density at the boundaries, we can thus
define the ionisation probability as $I(t) = 1 - N(t)$, where $N(t)$ is the charge
that remains inside the simulation box at time $t$. The lower panel of Fig. 2(a) shows $N$ as a 
function of $t$.
%we put an absorptive boundary condition at the outer
%boundary of the spherical grid at 10 \AA.  
%Assuming that any electron that reaches the boundaries is unbound, the total
%absorbed density corresponds to the ionisation probability.
%The lower panel shows the probability for the electron remaining inside
%the simulation box as a function of time.  
We see that there is practically no ionisation
for the lower fields, and only a 20\% ionization probability for the 
$2.1 \times 10^{12}$ W/cm$^2$ field.  Thus, in this range of intensities, the
laser dissociates the dimer without ionising it.

%In the first case we modeled the
%femtosecond-laser induced
%photodissociation
%of the Na$_2^+$ dimer, both at resonant and off-resonant frequencies.
%As within pseudopotential theory Na$_2^+$ is an one-electron system,
%its excitation energies can be obtained through differences of eigenvalues.
%As a cross-check of our computer implementation, we also calculated the optical excitation spectrum,
%which can easily be done by using this
%real-space technic \cite{pChelikowsky94,pYabana96,pCastro02}
%(shown as an inset in Fig.~\ref{fig:pes}). 
%Almost all spectral weight is concentrated in two peaks
%at about 2.5~eV and 3.2~eV, corresponding to
%excitations parallel and perpendicular to the the dimer axis, 
%respectively.
%These results 
%compare favourably to experimental data \cite{pEllert95}.
%%%% Check this reference (JCP 104 1427 (1996)!!!
%\begin{figure}[t]
%\includegraphics[scale=0.3,angle=-90]{graphs/linear_response.ps}
%\caption{\label{fig:optical_response}
%Photoabsorption cross-section of the Na$_2^+$ dimer.}
%\end{figure}
%Absorption at 2.5~eV corresponds to the 
%$1^2\Sigma_{g}^{+} \to \; 1^2\Sigma_{u}^{+}$
%transition (see Fig. \ref{fig:pes}).
%\begin{figure}
%\includegraphics[scale=0.33]{graphs/pes.ps}
%\caption{\label{fig:pes}
%Adiabatic energy surfaces of the Na$_2^{+}$ dimer,
%as obtained by our three-dimensional real-space code.
%Similar results can be found in Ref. \cite{pMagnier_Oise96}.
%In the inset, photoabsorption cross-section of the
%same molecule.}
%\end{figure}
%The $1^2\Sigma_{u}^{+}$ surface is anti-bonding, so
%exciting the ground state using
%a 2.5~eV laser pulse, polarized in the direction of
%the dimer axis, should 
%populate this excited state, and lead to
%dissociation, even at moderate intensities. 
%This is indeed the result we obtained. 
%In Fig.~\ref{fig:resonant_dissociation}~(a),
%we present a series of simulations using a 2.5\,eV laser,
%at intensities ranging from the weak $10^{10}{\rm\,W/cm^2}$
%to the moderate $2.1\times10^{12}{\rm\,W/cm^2}$.
%The total duration of the pulse was
%80\,fs (which corresponds to a spectral 
%width of about 0.1\,eV), and the envelope was a semicycle cosine.
%In the upper panel, we plot the internuclear separation of the dimer, and
%in the lower panel, the total electronic charge that remains
%inside the simulation box -- we assume that the electron
%density that exits this box corresponds to ionisation.
%Even at low intensity when no ionisation occurs,
%we obtain dissociation of the dimer. This is a clear
%signature of photodissociation on light-induced potentials.
%Above $10^{12}{\rm\,W/cm^2}$, as
%charge starts to leave the system, field ionisation  
%appears as a competing process.
%It should be noted that, due to the classical treatment
%of the nuclear dynamics, we can not distinguish all
%possible dissociation channels induced by a given laser. 
%First attempts to overcome this limitation, 
%by turning to a full quantum mechanical description, are currently 
%underway.


We next consider the excitation at the upper resonance frequency, 
$\omega = $ 3.2\,eV, corresponding
to an electric field perpendicular to the dimer axis. 
%Tuning now the frequency to the second resonance,
%at 3.2\,eV (corresponding to the $1^2\Sigma_{g}^{+} \to \; 1^2\Pi_{u}$
%transition), and polarizing the laser along a 
%direction perpendicular to the dimer axis,
%qualitatively different results are obtained  [see Fig.~\ref{fig:resonant_dissociation}~(b)].
Since the $1^2\Pi_{u}$ surface is bonding, 
[see Fig.~\ref{fig:resonant_dissociation}~(b)]
no dissociation
is expected unless the Coulomb explosion channel is opened through
ionisation. We see that the dimer remains bound over the entire range
of intensities that produced dissociation at the lower resonant frequency.

Finally, we also performed simulations at the
non-resonant frequency $\omega =$ 1.57\,eV,
the one used in Ref.
\cite{pGerber01}. Fig.~\ref{fig:resonant_dissociation}~(c)
shows how dissociation now occurs only at much higher intensities,
and it is mainly due to ionisation (almost absent in the resonant
calculations for the range of intensities used): 
since Na$_2^{2+}$ has no bonding states, ionisation is followed by 
Coulomb explosion.

Another process in which the nuclear motion may play an important role is high harmonic
generation.  Even harmonics may be created by irradiating HD with an intense laser pulse,
but not by irradiating H$_2$:  even harmonic generation
is forbidden for a centrosymmetric molecule.
In an adiabatic
treatment of the nuclear coordinates, the nuclear masses play no role and
the even harmonics can not appear.
This is no longer the case if non-adiabatic effects are taken into account, for the 
different masses of H and D break the symmetry.
Kreibich et al.  \cite{pKreibich01} studied this  process in 
a 1-D model with a full quantum mechanical treatment of the nuclear
motion, finding strong even harmonics at high
harmonic number.  To discern whether the classical treatment of nuclear
motion also produces these harmonics, 
we studied the same 1D problem within our framework.
%we have applied Eq.(1) to
%the same system.  
As in Ref.~\cite{pKreibich01},
we took the laser field to have a frequency 
%corresponding to photon energy 
of 1.6 eV, and an intensity that rises linearly to $10^{14}$ W/cm$^2$
over an interval of 10 optical cycles, and is held constant thereafter.
We then calculated the spectral intensity of the generated harmonics, $H(\omega)$:
\be \label{eq:harmonic_spectrum}
  H(\omega) \sim \left|\int\!\!d t\: e^{i\omega\: t}\:
  \frac{d^2}{d t^2} \left<\Psi(t)\right|\hat{e} \cdot \vec D \left|\Psi(t)\right>
  \right|^2\,.
\ee
We find that the
classical treatment does indeed produce even harmonics, but
much smaller than the quantum treatment.  The results
are shown in Fig. \ref{fig:harmonics.ps}.  The top left panel depicts
the harmonic spectrum for HD, and only odd harmonics are apparent.
However, it may be proved that the HD Hamiltonian already violates centrosymmetry
within our classical treatment, through a term of the form:
$$-\frac{1}{2}\left(\frac{1}{M_H}-\frac{1}{M_D}\right)P(t)\left(\hat{p}_1 + 
\hat{p}_2\right),
$$
where $P(t)=\frac{1}{2}\left(P_H(t)-P_D(t)\right)$ is the relative time-dependent
nuclear momentum
and $\hat{p}_i$ are the 
electronic momentum operators \cite{ref_frame}.  Its effect can
be enhanced by decreasing the nuclear masses.  In the bottom left
panel, the H and D masses have been decreased by a factor 100, and
then the second- and fourth-order harmonics become visible.  As a
qualitative check of the numerics, we also show the same graphs for
H$_2$, in which no even harmonics can occur.

Thus we see that on qualitative level the non-adiabatic dynamics generating
even harmonics are obtained with the classical treatment of the
nuclear coordinates. However, the quantum treatment may
be needed for a quantitative result.  
%In the full quantum theory, 
By describing the nuclei quantum mechanically,
the ground 
state violates centrosymmetry and the even harmonics can be generated
even if the nuclear motion is frozen.  In contrast, in the classical
treatment the ground state is symmetric and the symmetry violation only builds
up as the nuclei move.

%In order to further assess to which extent diabatic effects are included
%in this framework, we studied the harmonic spectra of the H$_2$ and HD
%molecules -- addressing the same problem studied by Kreibich {\em et al.}
%within a 1D, numerically exact quantum-mechanical scheme \cite{pKreibich01}:
%Changing the isotopic composition to HD, they observed the generation of the
%parity forbidden even peaks in the harmonic spectrum \cite{adiabatic_forbidden}, 
%concluding that non-adiabatic effects can alter crucially the emission spectrum.
%As in this work, we employed a 1.6\,eV, 10$^{14}{\rm\,W/cm^2}$ laser, switched on linearly
%during the first 10 optical cycles, and held constant afterwards.
%We also made use of the 1D approximation the sake of comparison
%\cite{3Dcalculations}.
%The harmonic spectrum was calculated from
%the dipole acceleration expectation value:
%\be \label{eq:harmonic_spectrum}
%  S(\omega) \sim \left|\int\!\!d t\: e^{i\omega\: t}\:
% \frac{d^2}{d t^2} \left<\Psi(t)\right|\hat{e} \cdot \vec D \left|\Psi(t)\right>
%  \right|^2\,.
%\ee
%Here $\vec D$ denotes is the total dipole operator of the system, including both
%electrons and nuclei.
%The results are shown in Fig. \ref{fig:harmonics.ps}. 
%No even harmonics seem to be present in the HD case
%(see top-left panel). This is not due to the absence
%of non-adiabatic couplings in our equations, as it is easy to proove that our
%field-free Hamiltonian is not symmetric under parity: Changing the reference
%frame to nuclear geometrical center, the Hamiltonian is broken into an even term, and
%a term of the form:
%$-\frac{1}{2}\left(\frac{1}{M_H}-\frac{1}{M_D}\right)P(t)\left(\hat{p}_1 + \hat{p}_2\right)$,
%where $P(t)$ is the classical nuclei momentum, and $\hat{p}_i$ are the electronic momenta.
%This latter term is odd, in contrast to the rest of the hamiltonian, and should
%lead to the presence of even harmonics. This is made clear in the botton panels, where
%the masses have been reduced by two orders of magnitude -- so that the magnitude
%of the odd term is enhanced. 
%Even harmonics are now clearly present in the HD case, as opposed to the fully 
%symmetric H$_2$ (on the right panels of the figure, corresponding calculations for H$_2$
%are shown).
%Therefore, even harmonics are present
%yet very weak in the present model, in contrast to the  full quantum mechanical
%description. This weakness may be due to (i) the fact that HD has a permanent dipole
%moment beyond our classical description of the ions, and (ii) a different treatment
%of the non-adiabaticity during the time evolution. 

In summary, we have examined the computational feasibility of
including nuclear dynamics in time-dependent density functional
theory 
%or time-dependent self-consistent field theory 
using
a pseudopotential code to study the femtosecond laser
induced dynamics of sodium dimers. 
Using this approach for treating the Na$_2^+$ dimer, 
we were able to distinguish different phodissociation regimes, ranging
from dissociation on light induced potentials, to field ionisation 
followed by Coulomb explosion. 
Electronic and ionic degrees of freedom are thus coupled, so that one can
observe the electron-phonon transfer of energy.
We also found, with another example, that non-adiabatic effects are 
present in the general treatment based on Eq. (1).
One of the major attractivenesses 
of this method resides in its reasonable
scaling behaviour when applied to larger systems.
We thus expect to be able to tackle problems like
photoisomerization or even photochemical reactivity in
systems of dozens of atoms in the near future.

This work was supported by the RTN program
of the European Union NANOPHASE (contract HPRN-CT-2000-00167),
Basque Country University,
Iberdrola S.A. and DGESIC (PB98-0345).
GB acknowledges support by the US Department of Energy under
Contract Nr. E-FG-06-90ER-411132.
Computer time was kindly provided by the CEPBA.
We thank E.~K.~U. Gross for enlightenment discussions.
AC thanks the University of Washington and the DIPC for kind hospitality.


\bibliographystyle{prsty}
\bibliography{abbrevs,biblio}

\end{document}

dimers.bbl
\begin{thebibliography}{10}

\bibitem{pBrabec00_yotros}
For a review, see T. Brabec and F. Krausz, Rev.~Mod.~Phys. {\bf 72},  545  (2000).
For recent developments, see
M. Drescher {\it et~al.}, Science {\bf 291},  1923  (2001) and
M. Hentschel {\it et~al.}, Nature {\bf 414},  509  (2001).

%\bibitem{pDreschner01}
%M. Dreschner {\it et~al.}, Science {\bf 291},  1923  (2001).
%
%\bibitem{pHentschel01}
%M. Hentschel {\it et~al.}, Nature {\bf 414},  509  (2001).

\bibitem{pLeone00_yotros}
L. N.-Glandorf {\it et~al.}, Phys.~Rev.~A {\bf 62},  023812  (2000);
Phys.~Rev.~Lett. {\bf 87},  193002  (2001).

\bibitem{faraday_discussions}
See, for example, articles in Faraday Discussions {\bf 115} (2000).

\bibitem{x-ray-citations}
Z. Chang {\it et~al}, Phys. Rev. Lett. {\bf 79}, 2967 (1997); C. Spiedelmann {\it et al},
Science {\bf 278}, 671 (1997).

\bibitem{pProbst97}
M. Probst and R. Haight, Appl.~Phys.~Lett. {\bf 71},  202  (1997).

\bibitem{pHuber01}
R. Huber {\it et~al.}, Nature {\bf 414},  286  (2001).

%\bibitem{pMacholm01}
%M. Machholm and N.~E. Henriksen, Phys.~Rev.~Lett. {\bf 87},  193001  (2001).

\bibitem{pDion96_yotros}
C. Dion {\it et~al.}, J.~Chem.~Phys. {\bf 105},  9083  (1996);
R.~W. Shoenlein {\it et~al.}, %L.~A. Peteanu, R.~A. Mathies, and C.~V. Shank, 
Science {\bf 254},  412  (1991).

\bibitem{pGerber01}
A. Assion, T. Baumert, U. Weichmann, and G. Gerber, Phys.~Rev.~Lett. {\bf 86},
  5695  (2001).

\bibitem{pYamanouchi02}
K.~Yamanouchi, Science {\bf 295}, 1659 (2002).

\bibitem{new_processes}
Examples commonly invoked are bond softening, vibrational population trapping,
molecular alignment and above threshold dissociation.

\bibitem{pKnight97}
M. Protopapas, C.~H. Keitel, and P.~L. Knight, Rep. Prog. Phys. {\bf 60},  389
  (1997).

\bibitem{atoms_and_dimers_studies}
Some of the works include exact 3D calculations for one and two electron 
atoms \cite{pParker96}, for one electron dimers \cite{pMagnier99}, and 1D models
for two electron systems \cite{pGrobe93_yotros,pLein01,pKreibich01}.

\bibitem{pParker96}
J. Parker {\it et~al.},% K.~T. Taylor, C.~W. Clark, and S. Blodgett-Ford, 
J.~Phys.~B {\bf 29},  L33  (1996).

\bibitem{pMagnier99}
S. Magnier, M. Persico, and N. Rahman, J.~Phys.~Chem.~A {\bf 103},  10691
  (1999).

\bibitem{pGrobe93_yotros}
R. Grobe and J.~H. Eberly, Phys.~Rev.~A {\bf 48}, 4664 (1993); D.~G. Lappas
{\it et al}, J.~Phys.~B {\bf 29}, L619 (1996).

\bibitem{pLein01}
M. Lein, E.~K.~U. Gross, and V. Engel, Phys.~Rev.~A {\bf 64},  023406  (2001).

\bibitem{pKreibich01}
T. Kreibich {\it et~al.}, %M. Lein, V. Engel, and E.~K.~U. Gross, 
Phys.~Rev.~Lett. {\bf 87}, 103901  (2001).

\bibitem{tddft}
E. Runge and E.~K.~U. Gross, Phys.~Rev.~Lett. {\bf 52},  997  (1984);
E.~K.~U. Gross, J.~F. Dobson, and M. Petersilka, in
{\it Density Functional Theory II}, edited by R.F. Nalewajski,
``Topics in Current Chemistry", Vol 181 (Springer, Berlin, 1996).

\bibitem{ALDA_problems}
The ALDA has nevertheless some well known problems, but most of them can be
solved by the swarm of exchange-correlation potential now available -- 
see, for example,  
M.~A.~L. Marques, A. Castro, and A. Rubio, J.~Chem.~Phys. {\bf 115}, 3006 (2001),
and references therein.

%\bibitem{pRungeGross84}
%E. Runge and E.~K.~U. Gross, Phys.~Rev.~Lett. {\bf 52},  997  (1984).

%\bibitem{bRunge96}
%E.~K.~U. Gross, J.~F. Dobson, and M. Petersilka, in
%{\it Density Functional Theory II}, edited by R.F. Nalewajski,
%``Topics in Current Chemistry", Vol 181 (Springer, Berlin, 1996).

%\bibitem{pAngel02}
%G. Onida, L. Reining, and A. Rubio, Rev.~Mod.~Phys. (2002).

\bibitem{pYabana96}
K. Yabana and G.~F. Bertsch, Phys.~Rev.~B {\bf 54},  4484  (1996);
Int.~J.~Quantum~Chem. {\bf 75}, 55 (1999).


%\bibitem{bertsch_yotros}
%A.~Rubio {\it et al.},
%Phys. Rev. Lett. {\bf 77} 247 (1996);
%K.~Yabana, and G.~F.~Bertsch,
%Z.~Phys.~D (Atom. Mol. Cl.) {\bf 42}, 219 (1997);
%Phys.~Rev.~A {\bf 60}, 1271 (1999);
%{\bf 58}, 2604 (1999);
%J.~I.~Iwata., K.~Yabana, and G.~F. Bertsch,
%J.~Chem.~Phys. {\bf 115}, 8773 (2001);
%G.~F.~Bertsch, A.~Schnell, and K.~Yabana,
%J.~Chem.~Phys. {\bf 115}, (2001);
%M.~A.~L. Marques, A. Castro, and A. Rubio, J.~Chem.~Phys. {\bf 115},  3006
%  (2001);
%A.~Castro {\it et al}, J. Chem. Phys. {\bf 116}, 1930 (2002).

\bibitem{review}
G. Onida, L. Reining, and A. Rubio, Rev.~Mod.~Phys. (2002) and references therein.

%\bibitem{marques}
%M.~A.~L. Marques, A. Castro, and A. Rubio, J.~Chem.~Phys. {\bf 115}, 3006 (2001).

\bibitem{pCalvayrac00_yotros}
F. Calvayrac {\it et~al.}, %P.-G. Reinhard, E. Suraud, and C.~A. Ullrich, 
Phys. Rep. {\bf 337},  493  (2000);
E. Suraud and P.~G. Reinhard, Phys.~Rev.~Lett. {\bf 85},  2296  (2000).

\bibitem{pKreibich01-2}
T. Kreibich and E.~K.~U. Gross, Phys.~Rev.~Lett. {\bf 86},  2984  (2001).

\bibitem{pFlocard78}
H. Flocard, S.~E. Koonin, and M.~S. Weiss, Phys.~Rev.~C {\bf 17},  1682 (1978).

\bibitem{pChelikowsky94_yotros}
  J.~R. Chelikowsky, N. Trouiller, and Y. Saad, Phys.~Rev.~Lett. {\bf 72},  1240
  (1994); A.~Rubio {\it et al.}, Phys. Rev. Lett. {\bf 77} 247 (1996); 
  T.~L.~Beck, Rev.~Mod.~Phys. {\bf 72}, 1041 (2000).

%\bibitem{pFlocard78}
%H. Flocard, S.~E. Koonin, and M.~S. Weiss, Phys.~Rev.~C {\bf 17},  1682
%  (1978).

\bibitem{pTroullierMartins91}
N. Troullier and J.~L. Martins, Phys.~Rev.~B {\bf 43},  1993  (1991).

%\bibitem{pBertsch00}
%G.~F.~Bertsch, J.-I.~Iwata, Angel Rubio and K.~Yabana, Phys.~Rev.~B {\bf 62},
%7998 (2000).

\bibitem{pMagnier96_yotras}
S.~Magnier, M.~Persico, and N.~Rahman, Chem.~Phys.~Lett. {\bf 262},  747
  (1996);
M.~Machholm and A.~Suzor-Weiner, J.~Chem.~Phys. {\bf 105},  971  (1996).

%\bibitem{pEllert95}
%Ellert, Phys.~Rev.~Lett. {\bf 75},  1731  (1995).
\bibitem{bonacic}
V. Bona\v{c}i\'c-Kouteck\'{y} {\it et~al.}, %, J. Pittner, C. Fuchs, P. Fantucci, M. F. Guest, and J. Kouteck\'{y},
J. Chem. Phys. {\bf 104}, 1427 (1996).



\bibitem{pMagnier_Oise96}
S. Magnier, F. Oise, and M. Seeuws, Mol.~Phys. {\bf 89},  711  (1996).

%\bibitem{adiabatic_forbidden}
%Even harmonics are forbidden if an adiabatic approximation is used.

%\bibitem{3Dcalculations}
%3D calculations are not
%computationally forbidden in our model and should yield qualitatively similar results.

\bibitem{ref_frame}
A canonical transformation to the nuclear geometrical center has been performed.

\end{thebibliography}
