I. INTRODUCTION
Density functional theory (DFT)Reference Hohenberg and Kohn1 opened the door for computational material design through the solution of the Kohn–Sham equationReference Kohn and Sham2 of electron wave functions in actual materials. The properties of several phases of materials can be determined from DFT calculations without any empirical parameters. The total energy and force formalismReference Ihm, Zunger and Cohen3 within DFT is a standard tool for examining the stability and thermal properties of condensed matter in the electronic ground state. On the other hand, experimental developments using femtosecond lasers for monitoring ultrafast dynamics,Reference Niikura, Légaré, Hasbani, Bandrauk, Ivanov, Vilieneuve and Corkum4 controlling chemical reactions,Reference Zare5 and fabricating materialsReference Shimotsuma, Hirao, Kazansky and Qiu6 have necessitated the elucidation of atomic-scale dynamics under electronic excitation. One of the powerful tools for simulating dynamics under electronic excitation is the application of time-dependent density functional theory (TDDFT)Reference Runge and Gross7to molecular dynamics (MD) within an Ehrenfest approximation.Reference Ehrenfest8
In this invited paper, recent computational approaches for simulating laser-induced dynamics of molecules encapsulated in carbon nanotubes (CNTs) are discussed. Since the chiral structure of CNTsReference Iijima9, Reference Iijima and Ichihashi10 and their capillary effectReference Ajayan and Iijima11 were first identified, the nature of the nanoscale hollow space inside CNTs has attracted considerable attention. CNTs have been capable of encapsulating not only metallic atoms but also fullerene,Reference Smith, Monthioux and Luzzi12 metallofullerene,Reference Kitaura and Shinohara13 water,Reference Maniwa, Kataura, Abe, Udaka, Suzuki, Achiba, Kira, Matsuda, Kadowaki and Okabe14 and coronene molecules.Reference Okazaki, Iizumi, Okubo, Kataura, Liu, Suenaga, Tahara, Yudasaka, Okada and Iijima15 The next point of interest was then whether these encapsulated molecules could interact with light irradiated from outside the CNT, and a method to measure the optical absorption of photoharvesting molecules encapsulated in CNTs became available in the mid 2000s.Reference Yanagi, Iakoubovskii, Kazoui, Minami, Maniwa, Miyata and Kataura16These experimental developments have motivated us to study the photoinduced microscopic dynamics of molecules encapsulated in CNTs, which is difficult to observe directly through experiments.
II. THEORY AND COMPUTATIONAL SCHEME
Since the purpose of the present simulation was to examine MD induced as a result of optical excitation by short strong laser pulses, the situation must be different from ordinary photoexcitation, which can be described within the framework of perturbation theory by using the Fermi golden rule. In our case, the intensity of ultrashort pulses is strong, and the time constant of the optical interaction is short (in the order of several femtoseconds). In such short interactions, the excited states are transient. Furthermore, the positions of atoms show significant displacement within short time intervals, often in the order of less than 1 ps. Based on these facts, the assumption of steady states for the electron wave functions is no longer valid.
An approach alternative to traditional perturbation theory is the direct calculation of real-time propagation of electron wave functions as well as atomic motion during irradiation with light. Thus, we consider numerical integration of the time-dependent Kohn–Sham equationReference Runge and Gross7in the presence of a laser field:
$$i\hbar {{d\psi _n^{{\rm{KS}}} (r,t)} \over {dt}} \,= (H_{{\rm{KS}}} (r,t) + V_{{\rm{ext}}} (r,t))\psi _n^{{\rm{KS}}} (r,t)\quad.$$ Here, we use time-dependent scalar potential 
$V_{{\rm{ext}}} (r,t)$ to express the optical field. 
$H_{{\rm{KS}}} (r,t)$ and 
$\psi _n^{{\rm{KS}}} (r,t)$ are the Kohn–Sham Hamiltonian and orbitals, respectively. Expressing optical fields by vector potentials is also useful for investigating electron dynamics.Reference Otobe, Yamagiwa, Iwata, Yabana and Nakatsukasa17At the same time, for electron–ion dynamics, the use of a scalar potential enables us to simultaneously account for interactions of light with both electrons and ions.Reference Castro, Marques, Alonso, Bertsch and Rubio18
The combination of MD of ions and real-time propagation of electron wave functions requires a numerically robust scheme implemented in plane wave code,Reference Sugino and Miyamoto19 which uses a split operator schemeReference Suzuki20 to compute the time evolution operator of the Kohn–Sham Hamiltonian including nonlocal pseudopotentials. The robustness of the simulation was tested using the rule of total energy conservation throughout the MD simulation. Here, we restrict ourselves to an adiabatic exchange-correlation approximation of TDDFT. For this purpose, we must compute the work performed by the laser field. If we introduce a time-dependent fictitious charge 
${\rm{\rho }}_{{\rm{ext}}} (r,t)$ to express the scalar potential of Eq. (1)as:
$$V_{{\rm{ext}}} (r,t) = \int {{{{\rm{\rho }}_{{\rm{ext}}} (r\prime,t)} \over {tnq#x007C;r\prime - rtnq#x007C;}}dr\prime\quad,} $$ we can compute the time derivative of the Coulomb interaction between the fictitious charge and ions as well as the interaction between the fictitious charge and valence electrons, under the condition that 
${\rm{\rho }}_{{\rm{ext}}} (r,t)$ is a product of both space- and time-dependent functions [i.e., as 
${\rm{\rho }}_{{\rm{ext}}}^0 (r)f(t)$]. The time derivative of the work can be expressed as:
$${{df(t)} \over {dt}}\int \!\!{{\rm{\rho }}_{{\rm{ext}}}^0 } (r) \left( \int {{{\rm{\rho }}(r \prime, t) + {\rm{\rho }}_{{\rm{ext}}}^0 (r\prime)f(t)} \over {tnq#x007C;r - r\primetnq#x007C;}}dr \prime\right. \left+ \mathop\sum\limits_I {{{Z_{\rm{I}} } \over {tnq#x007C; {r - R_{\rm{I}} (t)} tnq#x007C;} } } \right )dr\quad,$$ where 
${\rm{\rho }}(r,t)$ is the valence electron density and R I(t) and Z I are the position and valence charge of ion I. By numerically integrating Eq. (3) with respect to time, we can obtain the work performed by the optical field. In this case, the sum of the DFT total energyReference Ihm, Zunger and Cohen3 with the energy of all ions added and the work subtracted must be constant with respect to time.Reference Miyamoto and Zhang21Note that this formula contains the self-energy of the fictitious charge, which is not, however, an issue for the purposes of checking the numerical stability based on the conservation rule. The amount of energy increase due to the optical field can be evaluated only at the time when the optical field vanishes completely.Reference Miyamoto and Zhang22
The practical model is explained below. Since current methodsReference Sugino and Miyamoto19, Reference Miyamoto and Zhang21 use a plane wave basis set to express valence wave functions, three-dimensional periodic boundary conditions must be used. A schematic representation of our unit cell is shown in Fig. 1, where the location of a CNT is shown in a box of vacuum.
Schematic representation of the unit cell used in the present numerical simulation.
The x-axis in Fig. 1 corresponds to the direction of polarization of light, in which the polarity of the scalar potential is rapidly inverted to a direction away from the CNT region in the vacuum region. This sawtooth-type potential varies with time. For practical settings, we set the corresponding vacuum distance (the distance from the wall of one CNT to that in the next neighboring unit cell, i.e., the interwall distance) to 20 Å. At the same time, on the y-axis, the interwall distance between CNTs was set to 10 Å to avoid artificial intertube interactions. The z-axis in Fig. 1 is parallel to the CNT axis. We considered four or two k points of the Bloch wave function in the direction of the z-axis when the period of the CNT placement was set to be the same as or double the primitive period along with the z-axis, respectively. To express the interaction between valence electrons and ions, norm-conserving nonlocal pseudopotentialsReference Troullier and Martins23 were used. Regarding the exchange-correlation potential of DFT, the local density approximationReference Perdew and Zunger24 fitted with numerical data for homogeneous gasReference Ceperley and Aldar25 was used. The cutoff energy of the plane wave basis was set to 60 Ry. MD simulations were performed in accordance with the Hellmann–Feynman forcesReference Ihm, Zunger and Cohen3 within the Ehrenfest approximation,Reference Ehrenfest8 which is equivalent to following the MD trajectory of the average of several adiabatic potential energy surfaces.Reference Tully26 The results of validation and the limitations of this approximation are discussed in relevant literature introduced in Ref. Reference Marques, Maitra, Nogueira, Gross and Rubio27.
III. DISCUSSION
In this section, the first subsection presents the penetration of the optical field into a semiconducting CNT, which shows the feasibility of optical excitation of molecules encapsulated in CNTs. The second and the last subsections present the dynamics induced by ultrashort laser pulses to hydrochloride molecules and acetylene molecules located inside CNTs.
A. Electric field modulation inside semiconducting CNTs
Before presenting the interaction between light and molecules encapsulated in CNTs, it is valuable to discuss the penetration of optical fields into the hollow space of a CNT. Depending on the chirality of the CNT, its electronic structure can be semiconducting, semimetallic, or metallic.Reference Hamada, Sawada and Oshiyama28, Reference Saito, Fujita, Dresselhaus and Dresselhaus29 A strong screening effect emerges in metallic CNTsReference Ajiki and Ando30 when the polarization of the optical field is perpendicular to the CNT axis (cross polarization). In this regard, optical excitation with cross polarization was discovered in semiconducting CNTs,Reference Lefebvre and Finnie31raising questions about the ability of optical fields to penetrate into semiconducting CNTs.Reference Zhang and Miyamoto32
Figure 2(a) shows the schematics of the penetration of an optical field into a (14,0) semiconducting CNT. The polarization of light is set to be perpendicular to the CNT axis, matching the condition shown in Fig. 1. The frequency of the light corresponds to a wave length of 822 nm, in resonance with the excitation energy of (14,0) CNTs. The induced field originating from electron dynamics on the CNT wall emerges in a direction opposite to the applied filed. When the optical frequency is not in resonance with the CNT excitation energy, the amplitude of the induced field is smaller than that of the applied field and its direction is opposite to that of the applied field. Therefore, the optical field becomes weaker inside the hollow space of the CNT. However, as apparent from Fig. 2(b), when the optical frequency is in resonance with the CNT excitation energy, the amplitude of the induced field exceeds the amplitude of the applied field, such that the optical field in the hollow space of the CNT is stronger than that of the applied field due to an antiphase response to the applied field.
(a) Schematic representation of an optical field in the vicinity of a semiconducting CNT. Blue arrow denotes the applied field; red arrow denotes the combined applied and induced fields in the hollow region of the CNT. (b) Time evolution of the applied field and the combined applied and induced fields for a (14,0) semiconducting CNT.
The amplification of optical fields in the hollow space of a (14,0) CNT can also be understood as the accumulation of potential energy upon irradiation of light. Figure 3 shows the time evolution of the total energy in the presence of work performed by the optical field and the conserved amount in accordance with the energy conservation rule.Reference Miyamoto and Zhang21As mentioned above, the self-energy term of a fictitious charge representing the optical field is included in the potential charge that causes a variation in potential following the optical frequency. At the same time, the drift of the center of mass of the potential indicates an actual increase in potential energy due to the resonance.
Increase in potential energy of a (14,0) CNT upon irradiation with a pulse with an optical wave length of 822 nm. The conserved quantity (combined potential energy and kinetic energy of ions, with the work performed by the optical field subtracted) is also denoted with a dotted line.
Thus, penetration of optical fields into semiconducting CNTs has been demonstrated, showing optimistic perspectives for laser-induced photochemical reactions inside CNTs. In Sec. III. B and III. C, the dynamics induced by compressed laser beams (laser pulses) is presented.
B. Case of hydrochloride molecules under ultrashort pulses
Ultrashort laser pulses with full-width at half-maximum (FWHM) of less than 5 fsReference Adachi, Yamane, Morita and Yamashita33 can destroy chemical bonds inside molecules. According to the Fourier transform of the pulse along the time axis, higher energy components are also present in the pulse, despite the fact that the typical wave length is around 800 nm. A hydrochloride (HCl) molecule was tested within the present computational scheme, and spontaneous breaking of the bond between hydrogen (H) and chlorine (Cl) was observed with a pulse shape as displayed in Fig. 4(a) at a wave length of 800 nm and a FWHM of 2 fs. The maximum intensity of the pulse was 12 V/Å. Such disintegration was also observed when a HCl molecule was encapsulated in a (8,0) CNT, as shown in Fig. 4(b).Reference Miyamoto, Zhang and Rubio34
(a) Pulse shape of an ultrashort laser pulse. The dotted line is a Gaussian envelope with a FWHM of 2 fs. (b) Dynamics of hydrochloride molecule inside a semiconducting (8,0) CNT between 0 and 40 fs. Open circles denote carbon atoms; hatched large and small circles denote chlorine and hydrogen atoms, respectively. Red arrow shows the polarization of the laser field; small arrows in the right panel denote the direction of significant atomic displacements.
The HCl molecule can take one of several (meta)stable orientations inside a (8,0) CNT, and Fig. 4(b) shows the case when the orientation of the HCl molecule is perpendicular to the CNT axis. The simulation shows that after the pulse, a H atom is spontaneously ejected, while the Cl atom remains stationary. The arrows in the right panel of Fig. 4(b) denote the trajectory of the current MD, showing significant displacement of carbon (C) atoms due to collisions between ejected H atoms and the CNT wall. Upon continuing the simulation, the formation of a monovacancy initiated by a collision was observed. At the same time, when the HCl molecule was assumed to be parallel to the CNT axis, the H–Cl bond showed stretching but no disintegration with the same pulse as that in Fig. 4(a). If the H–Cl axis of the H–Cl molecule was tilted from the orientation of that shown in Fig. 4(b), at the moment of laser shot, the H–Cl dissociation also occurred but the ejected H atom did not strongly collide with CNT wall,Reference Miyamoto, Zhang and Rubio34 resulting in no damage to the CNT. The dynamics of the HCl molecule was thus concluded to depend on its orientation with respect to the CNT axis. We also tested a HCl molecule encapsulated in a boron nitride (8,0) tube and found that the BN tube disturbed the dynamics in the case of the same pulse.Reference Miyamoto, Zhang and Rubio34Although the energy gap of the BN tube is wider than that of the CNT, the BN tube was not completely transparent to the pulse. This was attributed to accidental absorption peaks in higher energy regions of the BN tube, which may match the Fourier components of the pulse in the frequency domain.
C. Case of acetylene molecules under ultrashort pulses
From the computationally obtained knowledge of ejection of H atoms from HCl molecules inside CNTs, high-yield production of hydrogen can be expected from any other H-containing molecules encapsulated inside CNTs at high density. Acetylene (C2H2) molecules were selected as our targetReference Miyamoto, Zhang and Rubio35 after HCl. C2H2 molecules can be exothermically encapsulated in (8,0) and (14,0) CNTs. Especially in the case of (14,0) CNTs, C2H2 molecules are stable, forming dimers as illustrated in the upper panel of Fig. 5. Under irradiation of short laser pulses with a maximum intensity of 12 V/Å, as illustrated in Fig. 4(a), C2H2 molecules undergo rotation around their C–C axes in a synchronized manner, as shown in the lower panel of Fig. 5. At 90 fs after pulse irradiation, one of the C–H bonds of the C2H2dimer is completely broken, causing the ejection of H atoms, as in the case of HCl molecules encapsulated inside (8,0) CNT. In contrast, ejection of hydrogen atoms was not observed during the MD simulation when the maximum intensity of the laser pulse was attenuated to 10 V/Å.
(Upper) Two acetylene molecules encapsulated in a (14,0) CNT in one of several (meta)stable configurations. (Lower) Dynamics of these two molecules upon irradiation with the laser pulse shown in Fig. 4(a), with a maximum intensity of the field of 12 V/Å. In the simulation, the frame of the (14,0) CNT is not shown to make the acetylene molecules visible.
The synchronized motion of C2H2 molecules was also numerically analyzed. Figure 6(a) shows a plot of the four C–H bond lengths of C2H2 dimers in the first 50 fs of the simulation. In this panel, no correlation can be observed between the four C–H lengths beyond 10 fs. At the same time, Fig. 6(b) shows the correlated motion of the C–C polar angles of these two molecules, which are illustrated in the inset of Fig. 6(b), for a longer time interval. The synchronized dynamics lasted up to 100 fs into the simulation and broke with the ejection of the H atom.
(a) Time evolution of the four C–H bonds of C2H2 dimer inside a (14,0) CNT. (b) Time evolution of C–C polar angles (hatched areas in the inset) of the C2H2molecules.
It should also be noted that this type of synchronized motion is unique to the case of molecules encapsulated inside CNTs. Free C2H2 dimers were also examined, and the detailed dynamics were different from those inside a CNT,Reference Miyamoto, Zhang and Rubio35clearly indicating the modulation of the optical field by the CNT.Reference Zhang and Miyamoto32
IV. CONCLUSION
MD simulations coupled with TDDFT were applied for simulating dynamics induced by strong laser pulses. Optical fields can penetrate semiconducting CNTs, reaching the hollow region, which can be used to turn semiconducting CNTs into nanoscale test tubes for photochemistry. Two examples of MD simulations for molecules encapsulated in CNT were presented. In both cases, the CNT itself also changed its shape depending on the strength of the laser field. In the case of a HCl molecule inside an (8,0) CNT, the damage on the CNT was sensitive to the molecular orientation of the HCl molecule with respect to the axis of the tube. The larger diameter (14,0) CNT can result in a more complex reaction with the encapsulated molecule, which is highly unlikely to occur in gas phase due to the almost instantaneous dimerization of free atoms into diatomic molecules. Further computations are expected to result in the development of exotic chemical reactions occurring in such confined spaces, paving the way to new horizons in photochemistry using CNTs as test tubes. If simulation times can be extended to the order of picoseconds, we might be able to answer the fundamental question of whether TDDFT shows the semiequilibrium energy distribution of excited electrons under finite electronic temperature, as experimentally observed.Reference Hertel and Moos36 The reason for a potential negative answer to this question might be the limitation of single particle representation of the Kohn–Sham orbitals and the adiabatic assumption on the exchange-correlation potential. In that case, TDDFT-MD simulations must be connected to semiequilibrium simulationsReference Jeschke, Garcia and Bennemann37assuming the Fermi energy distribution in electron state occupation at a given electron temperature.
ACKNOWLEDGMENTS
The present work was conducted under collaboration with H. Zhang and A. Rubio. All calculations presented here were performed by using the Earth Simulator. The author is also grateful for the support provided by the Research Organization of Information Science and Technology at Tokyo office.
