Dynamics of He2∗ triplet state excimer bubbles in superfluid 4He

Chemical Physics 332 (2007) 304–312 www.elsevier.com/locate/chemphys

Dynamics of He2 triplet state excimer bubbles in superfluid 4He J. Eloranta

*

Department of Chemistry and Biochemistry, California State University at Northridge, 18111 Nordhoff Street, Northridge, CA 91330, United States Received 28 September 2006; accepted 14 December 2006 Available online 21 December 2006

Abstract Time-dependent density functional theory calculations for bulk superfluid 4He were carried out to model dynamics around He2 excimers after optical excitation from the 3a to 3d state. The liquid dynamics occurring after a sudden change in the helium–liquid interaction results in interfacial dynamics, which can be divided into three different modes: (1) non-linear processes yielding shock and solitonic progressions, (2) fast interfacial dynamics related to thinning of the liquid–gas interface that occurs within few picoseonds and (3) slow ˚ ) in spherical breathing motion of the liquid–gas interface with recursion times up to 110 ps. The long-range repulsive tail (R > 12 A the He–He2 interaction is found to play an important role in determining the recursion time of the solvent cavity breathing mode. As energy differences of just few wavenumbers in this region are sufficient to produce large changes in the recursion time, none of the pair potentials derived from the first principles could reproduce the experimental data [V.A. Benderskii, J. Eloranta, R. Zadoyan, V.A. Apkarian, J. Chem. Phys. 117 (2002) 1201]. Therefore it is concluded that the pump–probe experiments measure energy differences that are not possible to calculate using the current electronic structure methods. The results obtained from the density functional theory calculations are consistent with the proposed experimental scheme.  2006 Elsevier B.V. All rights reserved. Keywords: He2 excimer; Density functional theory; Superfluid helium; Dynamics

1. Introduction Bulk liquid helium has been subject to many theoretical studies after the discovery of its superfluidity by Kapitza [1]. Shortly after this, Landau developed the well-known phenomenological two-fluid model to describe superfluid helium [2,3]. This theory assumes that a fluid dynamic description can be applied for two separate fluid fractions. The superfluid fraction has no viscosity associated with it and is described by the Euler equation. The liquid viscosity is completely due to the normal fluid fraction, which can be modeled by the Navier–Stokes equation. At low temperatures the superfluid fraction dominates, but upon increasing temperature the contribution of the normal fluid fraction increases and finally reaches unity at the k-point. The superfluid fraction constitutes the ground state of a

*

Tel.: +1 818 677 2677; fax: +1 818 677 4068. E-mail address: [email protected]

0301-0104/$ - see front matter  2006 Elsevier B.V. All rights reserved. doi:10.1016/j.chemphys.2006.12.010

high-dimensional quantum system, which is very difficult to model theoretically. Approximate analytical treatments have been especially valuable in developing understanding of the basic properties of the quantum liquid [4]. To obtain more accurate and generic description, numerical methods must usually be applied. Even in this case the high dimensionality of the problem poses a serious challenge. Quantum Monte Carlo (QMC) based methods overcome the problem by efficient sampling of the high-dimensional space and avoid explicitly storing the many-body wavefunction. If, however, explicit description of extended bulk superfluid is required, the number of particles quickly becomes too large for any QMC based method. The most successful approaches have been variational Monte Carlo (VMC), diffusion Monte Carlo (DMC), and path integral Monte Carlo (PIMC) methods [5,6]. The basic forms of these methods can be used for finding the ground state solution but extensions have also been formulated for finding the excited states. The PIMC method can also consider finite temperature effects and thus can provide a detailed

J. Eloranta / Chemical Physics 332 (2007) 304–312

microscopic understanding of superfluidity. A general deficiency of these methods is that they can only provide stationary solutions or at most dynamics in imaginary time. Another common approach for dealing with the dimensionality problem is to find a dimension reduction mapping such that the problem becomes effectively a single particle Schro¨dinger equation with a non-linear potential term. This is the general approach taken in bosonic density functional theory for 4He (DFT) [6,7]. It can be shown through the Madelung transformation that time-dependent DFT is very closely related to the fluid dynamic treatment proposed by Landau [8]. In the present study, time-dependent DFT is applied to describe dynamics in superfluid 4He. Among the first experiments providing information about the bulk liquid helium on molecular scales were the absorption and emission experiments on solvated electrons, He* and He2 excimers [9–13]. By comparing the spectral line shifts and shapes with the gas phase data, it was possible to obtain information about the solvation cavity structures surrounding these ‘‘impurities’’. The shapes of such solvation cavities are essentially dictated by the electronic structure of the impurities [14–16]. The optical pump–probe experiments by Benderskii et al. have demonstrated that it is possible to interrogate superfluid helium bath dynamics by using molecular probes (e.g. He2 excimers) [17]. By changing the sample cell temperature, it was further shown that the recurrence time of the pump– probe signal tracked the super/normal ratio known from the macroscopic scale experiments. Hence it was shown that the macroscopic concept of superfluidity, in terms of the two-fluid model, also applies on the microscopic scale. The proposed experimental scheme considered that the spherical 3a and 3d Rydberg states of He2 excimers were responsible for creating the observed liquid bath motion. Previous time-independent DFT calculations have shown that the 3a and 3d excimers occupy spherical cavities (‘‘bub˚ , respectively bles’’) in the liquid with radii of 7 and 13 A [18]. An optical transition from the 3a to 3d state would therefore generate an oscillating spherical bubble as the surrounding liquid attempts to find its new equilibrium geometry around the expanded 3d state excimer. The pump–probe signal was proposed to originate from the time-dependent multi-photon ionization process (e.g. the fluorescence depletion scheme). When the bubble is in its contracted shape, the ionization pathway is suppressed and the fluorescence can be detected whereas the extended cavity shape allows ionization of the excimers and consequently suppression of the fluorescence [17]. Thus the pump–probe experiment determines the recursion time for the bubble breathing motion, which was observed to be in the order of 150 ps depending on the temperature. The motion was damped after one cycle indicating very efficient energy dissipation via sound emission as well as viscous drag on the bubble edge. The theoretical interpretation of the experiment was carried out using the two-fluid model with an effective fluid dynamic potential describing the excimer–bulk helium system. This potential

305

was obtained by fitting the theoretical recursion times to the experimental results at various temperatures and pressures. However, the effective fluid dynamic potential is not straightforward to interpret. For example, it contains the convoluted impurity – liquid potential, surface dynamic components, surface tension terms, and partially missing quantum effects (e.g. quantum pressure) [16]. For this reason, the present study attempts to model the observed bubble dynamics at 0 K by applying time-dependent DFT and verify the details of the proposed experimental scheme. 2. Methods The experiments reported in Ref. [17] used short laser pulses to prepare the 3d state of He2 . In the following, the mechanism of robust optical preparation of the 3d state from a metastable 3a He2 excimer population in liquid 4He is considered. If the exciting laser pulse (e.g. pump) is sufficiently short (80 fs) then it can be assumed that the surrounding liquid is stationary during the pulse. The liquid does, however, modify the energetics of the He2 excimer Rydberg progression significantly [18]. The interaction between the Rydberg electron and the liquid is strongly repulsive because of the Pauli exclusion principle. Thus the high-lying Rydberg states have larger spatial electron extents and they experience increasing solvent induced spectral blue shifts. This is shown in Table 1 where the energies of the relevant Rydberg states have been calculated by the time-independent DFT method. The diatomic constants describing the He2 potential energy surfaces as well as the transition dipole moments connecting the electronic states were assumed to have their gas phase values [19,20]. The present numerical calculation employs the standard semi-classical field–matter interaction Hamiltonian: H ¼ H el ðr; RÞ þ H nucl ðRÞ þ EðtÞ  r

ð1Þ

where r refers to the electronic and R to the nuclear degrees of freedom. The first term represents the electronic Hamiltonian (i.e. a static n-level system), the second term is the nuclear Hamiltonian corresponding to Eq. (5), and the last term accounts for the time-dependent oscillating electric field produced by a short laser pulse (transform-limited). The wavefunction is taken to be a product of the electronic and the nuclear parts, which implies the Born–Oppenheimer approximation. A total of six electronic states (e.g. Table 1 Energies (cm1) and transition dipole moments (Debye) for the He2 electronic Rydberg states State

Gas phase

Liquid phase

Transitions

3

0 4768 10,890 20,392 21,507 21,551

0 4778 11,040 23,270 28,900 62,300

b(9.5), c(7.6), e(3.7) a(9.5), d(3.1), f(4.5) a(7.6), d(6.4), f(8.5) b(3.1), c(6.4), e(22.1) a(3.7), d(22.1) b(4.5), c(8.5)

a b 3 c 3 d 3 e 3 f (R) 3

306

J. Eloranta / Chemical Physics 332 (2007) 304–312

3

a, 3b, 3c, 3d, 3e and 3f), each with 20 vibrational levels, were included in the calculations. Time propagation of the resulting discrete-level system was carried out by using a similar predictor–corrector method that was applied in the following DFT calculations [16]. Finally, it should be noted that this calculation does not involve the rotating wave approximation (RWA). Dynamics of the helium bath was modeled by timedependent DFT. The energy functional of Dupont-Roc et al. (DR) for describing superfluid 4He is defined as (He2 impurity center included; atomic units assumed throughout,  h = 1) [21]: Z 1 1 2 E½W ¼ jrWðrÞj d3 r þ 2M He 2 Z Z c  qðrÞU l ðjr  r0 jÞqðr0 Þd3 r0 d3 r þ 2 Z Z  qðrÞðqðrÞÞ1þc d3 r þ qðrÞU He2 –He ðrÞd3 r ð2Þ

with q(r) = jW(r)j2 and with the following auxiliary functions defined: 8 h  i < 4e r 12  r 6 when r P h r r U l ðrÞ ¼   4 : U l ðhÞ hr when r < h Z ð3Þ qðrÞ ¼ qðr0 ÞPh ðjr  r0 jÞd3 r0 where  3 when r 6 h 3 Ph ðrÞ ¼ 4ph 0 when r > h ˚ and the model where the screening distance h is 2.377 A ˚ 3+3c and 2.8, respecparameters c and c are 10,455,400 K A tively [21]. The Ul term assumes a He–He Lennard-Jones interaction with screening at short distances (e = 10.22 K ˚ ). The second auxiliary function in Eq. and r = 2.556 A (3) provides spherical averaging for a given density. U He2 –He denotes the He2 –He pair potential and assumes pair-wise additivity. An extension of Eq. (2) has been proposed by Stringari et al. (SD), which reproduces the known superfluid 4He dispersion relation more exactly than DR [22,23]. The backflow functional term in the SD model was found earlier to have negligible effect on bubble edge dynamics and was therefore ignored in the present calculations [16]. Stationary solutions (e.g. dE/dW* = 0) to Eq. (2) can be obtained via self-consistent treatment of the time-independent DFT equation (formally non-linear time-independent Schro¨dinger equation): 

1 DWðrÞ þ U ðq; rÞWðrÞ ¼ lWðrÞ 2M He

approach was used in obtaining the static liquid structure around the initial state He2 (3a) excimer. In the present 1D calculations, Eq. (4) is subject to homogenous Dirichlet ˚ ) and Neuboundary condition close to the origin (3.0 A ˚ mann condition at the outer boundary (40 A). After the static solution W(r, 0) around the 3a state excimer was determined, time evolution of the liquid surrounding the optically prepared 3d state was obtained by solving the time-dependent DFT equation: i

oWðr; tÞ 1 ¼ DWðr; tÞ þ U ðqðr; tÞ; rÞWðr; tÞ ot 2M He  lWðr; tÞ

ð5Þ

This equation results from the least action principle, which requires that the action is stationary with respect to variations dW*(r, t). For a numerical solution of Eq. (5), the previously outlined predictor–corrector method was used in spherical coordinates [16,24,25]. In the present case, the interaction between the excited state He2 (3d) and He atoms is isotropic and therefore restriction to 1-D treatment involves no approximations. A time step of 1 fs and a regular spatial ˚ with 0.2 A ˚ grid step were applied grid extending over 500 A ˚ of the spatial grid was in the calculations. The last 50 A allocated for absorbing boundary layer as described earlier [16]. It should be noted that solution of the present problem in 3-D is computationally prohibitive. The slow bubble edge breathing motion will generate long-wavelength phonons and therefore a large spatial sphere radius is required for describing the dynamics properly. All calculations were carried out on dual processor Linux systems (AMD Athlon). For the 3a state the He2 –L–He potential can be obtained from the ab initio pair potential because the many-body effects are small [20]. However, for the 3d state many-body effects appear to be important and therefore various different He2 ð3 dÞ–He pair potentials were tested: He2 (3d)–He ab initio potential [20], experimentally motivated effective fluid dynamic potential [17], adjusted excimer–liquid pair potential, and the electron–helium pseudopotential method [26]. The adjusted He2 (3d)–He pair potential was parameterized as follows: U He2 –He ðrÞ ¼ c1 ec2 r þ c3 ec4 r

2

ð6Þ

where c1, c2, c3, and c4 are potential parameters listed in Table 2. In the pseudopotential approach the electronic part is described by an atomic radial Schro¨dinger equation for a single electron with a core charge of +1 and an effec-

ð4Þ

where W(r) is a function related to one particle density by q = jWj2, U(q, r) is the effective many-body potential (or mean-field) corresponding to either DR or SD, and l is the chemical potential (e.g. average energy per boson). For the ground state solution W(r) can be taken as a non-negative real function. This time-independent

Table 2 Parameters for potentials p1, p2, and p3 (Eq. (6)) Potential

c1

c2

c3

c4

p1 p2 p3

16 18 20

0.1 0.1 0.1

1000 1000 1000

0.08 0.08 0.08

˚. Energy units are in cm1 and distance in A

J. Eloranta / Chemical Physics 332 (2007) 304–312

tive core potential (ECP) for describing the inner shell electrons. The He2 –He interaction was calculated by convoluting the electron–helium scattering pseudopotential with the He2 Rydberg electron density [26]. For use with the density functional method of Eq. (5), the liquid–excimer potential is defined as Z U He2 –He ðrÞ ¼ qel ðr0 ÞU el–He ðjr  r0 jÞd3 r0 ð7Þ where qel is the normalized Rydberg electron density and Uel–He is the electron–helium pseudopotential [26]. Alternatively, a local pseudopotential approach could be also used in Eq. (7) [27]. Polarization of the ground state He atoms by the Heþ 2 core is ignored in the calculations, as the liquid tends to reside far away from the charged core. At each time step of Eq. (5), a time-independent radial equation for the Rydberg electron is solved: Z 1  Dr Wel ðrÞ þ qðr0 ÞU el–He ðjr  r0 jÞd3 r0 Wel ðrÞ 2M el þ U ECP ðrÞWel ðrÞ ¼ kWel ðrÞ

ð8Þ

where Wel(r) is the Rydberg electron wavefunction and k is the corresponding eigenvalue. This equation assumes that the electron follows the liquid density q adiabatically during the simulation. Here UECP denotes the ECP due to the remaining inner shell electrons of the excimer and it is represented as U ECP ðrÞ ¼ A expðBrc Þ

ð9Þ

where A, B, and C are the ECP model parameters. Their values were determined by least squares fitting the resulting one-electron densities to the ab initio derived Rydberg electron natural orbital densities [20]. In Eq. (8) the 3a state corresponds to the ‘‘1s’’ ground state and 3d the ‘‘2s’’ state (second lowest eigenstate). Details of the numerical treatment of Eq. (8) have been given elsewhere [16]. At the origin homogenous Neumann boundary condition and at the outer boundary homogenous Dirichlet boundary condition were applied. The spatial electron grid was chosen to correspond to the liquid 4He grid. For analyzing the time-dependent spherical liquid density profiles around the excimer, it is convenient to define the bubble radius as Z Rb Z 1 qðrÞd3 r ¼ ðq0  qðrÞÞd3 r ð10Þ 0

307

file by using Rb and a. Previously, it has been shown that this function is able to approximate the liquid profile well with an exception of possible density oscillations at the liquid–vacuum interface [18]. By least squares fitting Eq. (11) to the instantaneous density obtained from time propagation of Eq. (5), the dynamics can be rationalized in terms of Rb and a. This allows separation of the liquid interface motion into two different coordinates and furthermore defines an approximate trajectory (Rb(t), a(t)) for the bubble interface motion. 3. Results The present calculations aim at explaining the timedependent dynamics of He2 excimers in superfluid helium as observed in the experiments of Benderskii et al. [17]. For this reason the liquid configuration around the 3a state was chosen as the initial state for calculating the liquid dynamics. The static structure for this state has been calculated earlier and the solvation bubble radius was found to ˚ [18]. The final state is taken to mimic the 3d be Rb  7 A state, which is responsible for initiating the radial liquid bubble dynamics. Optical preparation (e.g. the pump) of the 3d state from a metastable 3a population must involve a non-linear absorption process. This was modeled by solving the time-dependent behavior of a 6 (electronic) · 20 (vibrational)-level system with parameters given in Table 1 (see Eq. (1)). This system was subject to an oscillating electric field originating from a short laser pulse. By setting the laser intensity to 1013 W/cm2, the laser pulse time envelope to Gaussian with full width at half maximum (FWHM) 80 fs and wavelength to 790 nm, the time-dependent populations for the 3a, 3c, and 3d states were obtained. The result from this calculation is plotted in Fig. 1. The final population of the 3d state was observed not to depend strongly on the laser pulse characteristics as its contribution in the resulting superposition state was always in the range of 10–20% over the intensity range 1011–1014 W/cm2

Rb

where Rb denotes the bubble radius. A simplified description for the liquid density profile is given by [26]:  0 when r 6 R0 qðr; R0 ; aÞ ¼ aðrR0 Þ q0 ð1  ½1 þ aðr  R0 Þe Þ when r > R0 ð11Þ where R0 specifies the position of the density profile and a the thickness of the profile. To a good approximation Rb  R0 + 2/a holds, which enables us to describe the pro-

Fig. 1. Time-dependency of the 3a, 3c, and 3d electronic states during the exciting laser pulse is shown. Laser pulse width is 80 fs (Gaussian), intensity 1013 W/cm2, and average energy 12,658 cm1.

308

J. Eloranta / Chemical Physics 332 (2007) 304–312

and pulse widths 80–150 fs. These calculations demonstrate the robustness of the 3d state preparation via intense short laser pulses, which were used in the earlier pump–probe experiments [17]. For the electron–helium pseudopotential method, the applied ECP in Eq. (9) was first calibrated. The ECP parameters were obtained as A = 10.6, B = 0.39 and C = 3.6 (a.u.). Comparison of the resulting Rydberg electron densities for both 3a and 3d states is shown in Fig. 2a and b. The results are very accurate for the 3a state and also qualitatively correct for the 3d state. This ECP model calculation predicts the ‘‘1s’’–‘‘2s’’ energy difference as 2.6 eV which is close to the known He2 3 a–3 d energy gap (2.5 eV) [19,20]. Calculation of the liquid structure around the He2 excimer using the electron pseudopotential approach leads to comparable bubble radii that were obtained from the pair-wise approximation using the ab initio pair potentials [20]. A general trend is, however, that this method ˚ for 3a state and yields slightly larger bubble radii, 0.8 A 3 3 ˚ for d state. For the a state the Heþ core polariza0.5 A 2 tion may partially explain this difference whereas for the 3d state this effect is negligible. The pseudopotential model should account properly for the compressibility of the

Fig. 2. Comparison of the Rydberg electron densities obtained from the earlier MRCI calculations (ab initio) and the radial one-electron treatment of Eqs. (8) and (9) (ECP) [20]. The upper panel (a) shows results for the 3a state and lower (b) for the 3d state.

Fig. 3. Calculated liquid densities around He2 excimer by using the pair potential approximation (ab initio) and the one-electron treatment of Eqs. (4) and (8) (ECP). Results for both 3a and 3d states are shown. For definitions of the models, see Table 2 and the text.

Rydberg electron orbital in a spherically symmetric environment whereas the He2 –He ab initio based potential cannot describe such many-body effects. Results from both the pseudopotential and ab initio based calculations are shown in Fig. 3. A summary of the results from the time-dependent DFT calculations (e.g. Eq. (5)) are given in Table 3 with the He2 – He potential parameters defined in Table 2. An overview of the applied pair potentials is shown in Fig. 4. As an example of a bubble edge trajectory, results from SD-p3 and SDab initio calculations are plotted in Fig. 5. Note the rapid damping of the bubble oscillatory motion for the steep ab initio potential. The effect of pressure on the bubble edge recursion time was calculated by using the experimentally known liquid densities at the given pressure and 2 K temperature in the simulations [28]. The results are given in Table 4 and a comparison between the calculated and the experimental data is plotted in Fig. 6 [17]. The pseudopotential method yields a similar recursion time as was obtained by using the ab initio pair potential (see Table 3). Most importantly, the bubble edge recursion times in both cases remain well below 100 ps. Results from the simulation employing the pseudopotential method are shown in Figs. 7 and 8 and the one-electron eigenvalues as a function of time in Fig. 9. At 0 K emission of sound (240 m/s) is the only source for energy dissipation because the superfluid fraction has no internal friction (e.g. viscosity). Prior to sound emission, generation of shock waves were observed for all potential models. These shocks had essentially the same properties as described previously [16]. When the liquid experienced a very repulsive potential, a soliton-type progression similar to Ref. [16] was also observed. Characteristic to solitons, no dispersion occurs during its time propagation. Finally, it must be remembered that in the present case these non-linear effects involve rather strong and sharp variations in the liquid density and therefore it is not clear how well the underlying DFT method can describe these effects.

J. Eloranta / Chemical Physics 332 (2007) 304–312

309

Table 3 Summary of results for the time-dependent DFT calculations ˚) ˚) Max. Rb (A Method Eq. Rb (A

t for max. Rb (ps)

˚) Min. Rb (A

t for min. Rb (ps)

Dsa (ps)

DR-fluid dynamics DR-ab initio SD-fluid dynamics SD-ab initio SD-pseudopotential SD-p1 SD-p2 SD-p3

36 17 34 15 20 45 47 49

10.0 13.2 9.9 11.5 12.5 10.8 11.6 12.2

84 43 79 40 49 100 106 111

2.8 1.6 2.7 1.5 1.6 1.0 1.1 1.2

12.3 12.4 11.9 12.2 12.9 12.8 13.7 14.5

16.6 17.3 15.9 15.3 15.4 16.8 18.2 19.5

The columns have following meanings: ‘‘Eq. Rb’’ is the bubble equilibrium radius, ‘‘Max. Rb’’ and ‘‘t for max. Rb’’ indicate the maximum bubble radius that was reached during the simulation and the time how long it took to reach this geometry, respectively; ‘‘Min. Rb’’ and ‘‘t for min. Rb’’ give the corresponding minimum values give for the bubble, which has returned to the inner turning point; Dsa is approximate time scale for interfacial (a) motion. Parameters for potentials p1, p2, and p3 are defined in Table 2.

Table 4 Pressure dependency of the bubble edge dynamics obtained from SD-p3 calculation Pressure (psi)

Max. ˚) Rb (A

t for max. Rb (ps)

Min. ˚) Rb (A

t for min. Rb (ps)

Dsa (ps)

Eq. Rb ˚) (A

0 20 40

20.0 17.6 15.4

52 41 33

12.4 11.5 10.5

116 93 74

0.8 0.8 0.7

13.8 12.8 11.5

For explanation of the columns see Table 3.

Fig. 4. An overview of the pair potentials used in the time-dependent DFT simulations. The ab initio potential was taken from Ref. [20] and the effective fluid dynamic potential from Ref. [17]. For definitions of the potentials, see Table 2 and the text.

Fig. 6. Experimental vs. calculated (SD-p3) pressure dependency for the bubble edge recursion time. The experimental data was taken from previously published results [17].

4. Discussion

Fig. 5. Bubble edge trajectories from SD-p3 and DR-ab initio calculations are shown. Note the difference in the bubble edge recursion time.

A steady state of triplet ground state He2 excimers (3a) in liquid 4He can be produced via the strong field ionization method [29]. Due to the favorable energy-level structure and large transition dipole moments in the triplet electronic manifold, excitation with intense laser pulses leads to efficient multi-photon processes. In the gas phase this is likely to proceed to complete ionization with high efficiency. However, in the liquid phase the excited

310

J. Eloranta / Chemical Physics 332 (2007) 304–312

Fig. 7. Bubble edge trajectory from a simulation using both electronic and liquid degrees of freedom (Eqs. (5) and (8); ECP). The liquid dynamics occurred around the ‘‘2s’’ (e.g. 3d) electronic state. The data shown corresponds to SD-pseudopotential.

Fig. 9. Time-dependency of the one-electron energies for different states is shown. These values are the eigenvalues of Eq. (8). The data shown corresponds to Figs. 7 and 8. The dynamics was run on the ‘‘2s’’ (e.g. 3d) state. 3

Fig. 8. Progression of the Rydberg electron density (‘‘2s’’ or 3d) during the simulation, where the time separation between each line is 1 ps, is shown. The data shown corresponds to the liquid trajectory of Fig. 7. The arrows indicate the direction of propagation for the Rydberg electron density and the inset shows magnification of the electron density at large distances.

Rydberg states experience large blue shifts with respect to the gas phase. Even though the number density of the Rydberg states increases towards the ionization limit, their energetics is greatly modified by the surrounding liquid. In this case, it becomes more difficult to find suitable resonances for efficiently climbing up the electronic manifold towards the ionization limit. The final state (e.g. 3d) can be considered to represent an ionization bottleneck state when the surrounding liquid is in a geometry dictated by the 3a initial state. This means that a resonant multi-photon ionization process would be strongly suppressed when the solvation bubble is its initial geometry. Identification of the 3d state as the bottleneck state was based on inspecting the spatial extents of the Rydberg electron densities and calculating the energies for the states when the liquid remains in the 3a state configuration (see Table 1). When the bubble proceeds closer to the new equilibrium geometry around the

d state, the higher excited Rydberg states are solvated even more efficiently as shown in Fig. 9. If a second laser pulse (e.g. the probe) is present at this point, the energy-level structure may now have reached a point, where ionization can proceed efficiently via the resonant multi-photon process. The present calculations demonstrate that a multi-photon process can populate the 3d state efficiently (see Fig. 1). Even though the levels are not exactly in resonance with the laser (k = 790 nm), the intense electric field and large transition dipole moments between the states allow efficient population of the excited states. The outcome of the single pulse excitation process is a superposition of many states, which after a random bath perturbation, collapse into a single state with the ensemble probability given by the corresponding basis function amplitude. From the experimental point of view, it is important that the preparation method is not sensitive to the laser pulse characteristics. This, of course, will ensure proper repeatability of the experiment. By varying the laser pulse intensity and the pulse width it was observed that the final 3d state population exhibits only small variations with respect to these parameters. This shows that the two-photon excitation scheme is a stable method for preparing the 3d state from a given 3a population. Once the 3d state is populated, the liquid begins to relax towards its new equilibrium geometry. Based on the simulations, the resulting bubble dynamics can be classified into three distinct processes: (1) highly non-linear shock wave and soliton generation, (2) bubble edge thinning (e.g. changes in a in Eq. (11)) and (3) a slow bubble edge breathing motion (e.g. Rb in Eq. (11)). The contribution of the initial non-linear dynamics increases with the repulsiveness of the excited state potential. This process provides a fast initial energy dissipation pathway. The solitonic waves, which were observed with the most repulsive potentials, had velocities 260 m/s, which is slightly higher than the speed of sound in liquid helium. They typically exhibit large devi-

J. Eloranta / Chemical Physics 332 (2007) 304–312

ations from the bulk density (up to 15%) and show essentially no dispersion during the propagation. Thus the repulsive head of the bubble edge runs away from the rest of the interface and this part can be mostly ignored in the simplified bubble edge dynamics. However, it must be kept in mind that this initial non-linear contribution is included in the effective fluid dynamic potential and complicates analysis of such data. The second process, namely the changes in the bubble edge thickness, occurs within the couple of first ps after which a remains fairly constant and only small changes in a occur at the turning points of the bubble edge trajectory (Rb(t)). Actually, the first and second processes cannot be completely separated as the initial bubble edge thinning also generates some nonlinear liquid excitations. The time scale for the bubble edge breathing motion (the third process) is essentially dictated by the long tail of the He2 –He pair potential and the bubble radius (see Fig. 4). In order to reach long bubble edge recursion times, the potential must be very shallow and ˚ radius (for the bubble Rb > 12 A ˚ ). This extend over 10 A will yield enough momentum for the bubble edge to reach recursion times greater than 100 ps as demonstrated in Fig. 5. It is also interesting to note that similar recursion times have been calculated for single-electron bubbles, which were proceeding towards the ground state equilibrium [16]. The present method does not allow an easy way to fit pair potentials to yield desired properties for the bubble trajectories and therefore a simple parameterization of Eq. (6) (‘‘adjusted potential’’) was used to match the experiment and to obtain a qualitative understanding of the system. In general, too steep potential leads to too high acceleration of the liquid edge and consequently strong energy dissipation via emission of sound. This effect is demonstrated in Fig. 5. Based on this consideration, the ab initio pair potential appears to be too repulsive in order to yield long time-scale (>100 ps) bubble breathing dynamics [20]. Simulations using the effective fluid dynamic potential yield recursion times slightly less than 100 ps whereas the previous fluid dynamic calculations gave 150 ps [17]. It should be noted that extraction of the pair potential from the effective fluid dynamic potential could only be carried out approximately. This assumes, for example, a constant liquid profile with a Heaviside function shape and the classical liquid surface tension term 4pR2b c ˚ 2). The difference in the recursion times (c  0.18 cm1 A originates from the approximations involved in this procedure. The response of the liquid around a 3d Rydberg electron was also obtained by solving both the one-electron (Eq. (8)) and time-dependent DFT (Eq. (5)) equations simultaneously. Before the simulations, the effective core potential (Eq. (9); see also Figs. 2 and 3) was calibrated. The oneelectron model was found to be only in qualitative agreement with the ab initio data. As shown in Fig. 2b, it was not possible to obtain an exact fit for the ab initio and the pseudopotential electron densities for the 3d state. This discrepancy is further reflected in the calculated liquid den-

311

sity profiles as shown in Fig. 3. The many-body effects should also result in small differences in the density profiles, but considering the limited accuracy of the pseudopotential model, it is not possible to extract this information from the calculations. The main motivation in considering the electronic degrees of freedom in the present calculations is related to the idea of ‘‘electron spring’’, which might be able to provide a smooth push of the surrounding liquid over extended range of bubble radii. It is obvious from the potentials shown in Fig. 4 and the recursion time data in Table 3 that this type of very smooth interaction would be required to obtain long time-scale bubble edge dynamics. However, as shown in Fig. 7, the recursion times remain close to the ones obtained using the ab initio pair potential. Thus, at least in the present calculation, the ‘‘electron spring’’ did not work as intended. This may be related to the fact that the basis set used in the multi-reference configuration interaction (MRCI) ab initio calculations provided a poor description of the Rydberg electron density at large distances. When designing Gaussian basis sets for chemical purposes, a good description of the electron density far away from the nuclei is not very important, as its contribution to the total energy is small. In the present case the situation is exactly the opposite. Since the ECP was calibrated against this data, it would naturally mean that the one-electron model suffers from the same weakness. An accuracy of few wavenumbers for the pair poten˚ ), which is very tial is required in this region (R > 12 A difficult to achieve with any theoretical approach. Behavior of the electron density as a function of time is shown in Fig. 8, where it can be seen that the 3d Rydberg electron orbital begins to expand as the liquid is moving out. During the bubble expansion, the electronic states exhibit varying degree of solvation. This can be seen in Fig. 9 where the energies of the ‘‘1s’’, ‘‘2s’’, and ‘‘3s’’ states are shown as a function of time. It can be clearly seen that the excited states exhibit more dramatic changes in energy. Based on these results, the Rydberg electron resonant ionization efficiency should be modulated as a function of the bubble radius Rb(t). Increase in the external pressure is expected to affect the bubble dynamics since this essentially introduces an opposing force that acts on the bubble edge. Therefore the bubble edge recursion time should decrease as the external pressure increases. By using the adjusted potential (p3 in Table 2) as an example, it can be seen that the bubble recursion time decreases as the external pressure is increased. This is demonstrated in Fig. 6 (see also Table 4). Furthermore, the pressure dependency of the recursion time tracks qualitatively the experimentally observed bubble recursion times [17]. It should be noted that the experimentally observed recursion times should be longer than the current 0 K calculations predict because the normal fraction imposes a viscous drag on the bubble edge. The present calculations demonstrate that it is possible to obtain the experimentally observed long recursion times for objects, which have sizes that are consistent with He2 excimers.

312

J. Eloranta / Chemical Physics 332 (2007) 304–312

By comparing the bubble edge recursion times (see Table 3), it can be seen that the SD functional predicts slightly smaller recursion times than DR. The difference is only 3–5 ps. To model the bubble edge breathing dynamics, it is sufficient to use methods that can describe only the initial linear part of the bulk helium dispersion relation properly (e.g. the Gross-Pitaevskii theory). The classical fluid dynamic treatment can essentially do this and, as such, it provides a very efficient way of simulating the bubble edge dynamics [17]. This model will not, however, be able to describe changes in the bubble edge thickness, or the initial highly dissipative non-linear dynamics. Even the current DFT methods may have problems in describing the non-linear dynamics since it involves liquid densities that vary strongly from the bulk density. 5. Conclusions This work provided theoretical verification for the experimental scheme of Benderskii et al. [17]. In that study, the dynamics was assigned to originate from the triplet He2 excimers, which were optically excited from the 3a to 3d state. The surrounding liquid was then allowed to relax and the 3d ! 3b fluorescence intensity was observed as a function of the delay time between the two laser pulses (e.g. pump and probe). The present calculations provide support for the robust optical preparation of the 3d state He2 excimers using short laser pulses, but fail to explain the observed long (>100 ps) bubble edge recursion times properly by using the pair potentials obtained from the first principles. The simulations yielded the correct time scales only when adjusted pair potentials, which had a longer spatial extent than the calculated ones, were used. If the accuracy problems with the ab initio potential were ignored, the electron bubble dynamics would fit in the experimental picture better [16,17]. However, thermalized electrons have conduction band (1 eV) well below the exciting laser energies and therefore it is unclear how solvated electrons could be probed in the pump–probe experiment. Theory at the current level can only provide a conclusion that the size of the object responsible for the observed dynamics is between a 3d state excimer and a free electron [17]. In fact, this is outcome is not be surprising as the accuracy required for the pair potentials is in the order of few wavenumbers ˚ . This clearly highlights that at distances greater than 10 A

the experiments are measuring interactions that cannot be accurately described by the present theoretical methods. References [1] [2] [3] [4] [5] [6]

[7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29]

P.L. Kapitza, Nature (London) 141 (1938) 74. L.D. Landau, J. Phys. (Moscow) 5 (1941) 185. L.D. Landau, J. Phys. (Moscow) 11 (1947) 91. R.P. Feynmann, Phys. Rev. 90 (1953) 1116. D.M. Ceperley, Rev. Mod. Phys. 67 (1995) 279. E. Krotscheck, J. Navarro (Eds.), Microscopic Approaches to Quantum Liquids in Confined Geometries, World Scientific, New Jersey, 2002. M. Barranco, R. Guardiola, S. Hernandez, R. Mayol, J. Navarro, M. Pi, J. Low Temp. Phys. 142 (2006) 1. E. Madelung, Z. Phys. 40 (1927) 322. M.A. Woof, Phys. Rev. Lett. 15 (1965) 235. W.S. Dennis, E. Durbin Jr., W.A. Fitzsimmons, O. Heybey, G.K. Walters, Phys. Rev. Lett. 23 (1969) 1083. B. Sethumadhavan, W. Yao, Y.H. Huang, R.E. Lanou, H.J. Maris, G.M. Seidel, Phys. Rev. Lett. 97 (2006) 015301. A. Ghosh, H.J. Maris, Phys. Rev. Lett. 95 (2005) 265301. S.G. Kafanov, Ya. Parshin, Ia.A. Todoshchenko, J. Exp. Theor. Phys. 91 (2000) 991. A.P. Hickman, N.F. Lane, Phys. Rev. Lett. 26 (1971) 1216. A.P. Hickman, W. Steets, N.F. Lane, Phys. Rev. B 12 (1975) 3705. J. Eloranta, V.A. Apkarian, J. Chem. Phys. 117 (2002) 10139. A.V. Benderskii, J. Eloranta, R. Zadoyan, V.A. Apkarian, J. Chem. Phys. 117 (2002) 1201. J. Eloranta, N. Schwentner, V.A. Apkarian, J. Chem. Phys. 116 (2002) 4039. G. Herzberg, Spectra of Diatomic Molecules, van Nostrand, New York, 1950. J. Eloranta, V.A. Apkarian, J. Chem. Phys. 115 (2001) 752. J. Dupont-Roc, M. Himbert, N. Pavloff, J. Treiner, J. Low Temp. Phys. 81 (1990) 31. F. Dalfovo, A. Lastri, L. Pricaupenko, S. Stringari, J. Treiner, Phys. Rev. B 52 (1995) 1193. M. Casas, F. Dalfovo, A. Lastri, L. Serra, S. Stringari, Z. Phys. 35 (1995) 67. W. Press, W. Vetterling, S. Teukalosky, B. Flannery, Numerical Recipes in C, 2nd ed., Cambridge University Press, New York, 1999. F. Scheid, Numerical Analysis, McGraw-Hill Book Company, New York, 1968. J. Jortner, N.R. Kestner, S.A. Rice, M.H. Cohen, J. Chem. Phys. 43 (1965) 2614. V. Grau, M. Barranco, R. Mayol, M. Pi, Phys. Rev. B 73 (2006) 064502. J. Wilks, The Properties of Liquid and Solid Helium, Clarendon Press, Oxford, 1967. A.V. Benderskii, R. Zadoyan, N. Schwentner, V.A. Apkarian, J. Chem. Phys. 110 (1999) 1542.