Deformation of an inner valence molecular orbital in ethanol by an intense laser field

Experimental and analytical details

A linearly polarized Ti:Sapphire laser pulse (~60 fs, ~800 nm, and ~140 μJ) was converted into a circularly polarized pulse by passing through an achromatic quarter wave plate and was focused on an effusive ethanol beam with an off-axis parabolic mirror (f = 200 mm). The effusive beam of deuterated ethanol CH₃CD₂OH vapor was continuously supplied into a vacuum chamber through a microsyringe (70-μm inner diameter) and a skimmer (0.2-mm orifice diameter; Beam Dynamics model 2). The base pressure of the chamber without the sample was below 1 × 10⁻⁸ Pa.

The three-dimensional momentum vectors of an electron and an ion from an identical molecule were measured in coincidence (18, 21). Ions and electrons created in the focal region were accelerated by an electrostatic lens (28) toward two microchannel plate detectors with delay line position encoding (RoentDek HEX80) on opposite ends of the vacuum chamber. Two-dimensional positions and time of flights were recorded using time-to-digital converters with a resolution of 25 ps (RoentDek TDC8HP). The laser repetition rate was 1 kHz, and the detection count rate was set to be less than 0.3 counts per laser shot.

The measured ion and electron momentum vectors (p→ion and p→ele) were used to obtain the relative angle φ_rel between the projected vectors onto the polarization plane (p→ionpol and p→elepol shown in fig. S1C). Laboratory frame electron momentum distributions have a torus shape, indicating that electrons are emitted mainly along the polarization plane (fig. S1, A and B). This electron motion suggests that the ionization proceeds in the tunnel ionization regime (25). RFPAD of the tunnel ionization was derived with φ_rel by taking account of the electron drift by the circularly polarized laser field (p→elepol⊥p→tunnel). Defining φ_rel to be positive in going in the same direction as the E-field rotation, we obtained the RFPAD as a function of φ_RFPAD = φ_rel − 90° (fig. S1C). The laser intensity was determined by a least-squares fit to a theoretical expression of the electron momentum distribution in a circularly polarized laser field (29).

The laboratory frame momentum distributions of the CD₂OH⁺ and CH₃CD₂⁺ fragment ions are almost isotropic and have a peak at the center (p_ion = 0) (fig. S2, A and B). The orientation of the parent molecule was determined from the direction of the fragment recoil based on the axial recoil approximation. In the present analysis, we selectively analyzed the coincidence events for the fragment ions satisfying the following two conditions. The first condition was that the out-of-plane angle of the ion emission with respect to the polarization plane was smaller than 10° (aionout‐of‐plane < 10°) (see also fig. S2D). This condition extracted the events producing the fragment ion recoiling along the polarization plane. The other condition concerned the fragment ion velocity υ_ion, which is the sum of the recoil velocity and the initial velocity of the parent molecule in the thermal distribution. The initial velocity spread of the parent molecule blurred fragment recoil direction from which the orientation of the parent molecule was determined. The velocity distribution of the parent molecule was approximated by that of the parent ion (fig. S2, A and B), because the recoil momentum [<1 atomic unit (au) as shown in fig. S1A] given by the electron tunneling and the circularly polarized laser field is smaller than the initial momentum spread of the parent ion [ΔpCH3CD2OH+jet= 14.3 au and ΔpCH3CD2OH+TOF= 2.4 au in full width at half maximum (FWHM)]. The velocity distribution indicated that most of the parent ions were slower than 2 × 10⁻⁴ au in the present experiment. Therefore, we imposed the condition of υ_ion ≥ 2 × 10⁻⁴ au (fig. S2, C and D).

Simulation method of tunnel ionization probability in laser electric field

The simulation procedure of the tunnel ionization probability [W_10a′(Φ, Θ)] is similar to the previous one (4). Briefly, we used the Kohn-Sham formalism in the DFT to calculate the tunnel ionization probability. Here, we considered a case in which the external field changes very slowly in time. More precisely, we assumed that the tunneling occurs much faster than one cycle of the external field. We also assumed that tunnel or above barrier processes are dominant and that the multiphoton process is negligible. Under these circumstances, the time-dependent Kohn-Sham equation (30) yields the static Kohn-Sham equation with an external dipole field, V_ext(r→) = −eE→⋅r→, where E is the strength of the electric field at a certain time, as{h[n(r→)]+Vext(r→)}φi(r→)=ϵiφi(r→)(3)

Here, the electrons in the molecule are continuously emitted so that the static Kohn-Sham orbitals φ_i(r→) must satisfy the outgoing boundary condition without any incident waves. This is the so-called Gamow state (31). Because of the outgoing boundary condition, the orbital eigenvalues, ϵ_i, are complex numbersϵi=ϵiR+iΓi(4)where ϵiR and Γ_i are the real and the imaginary parts of ϵ_i.

The imaginary part of the eigenvalue, Γ_i, is related to the ionization probability. To see it, we multiplied φi* to Eq. 3 and subtracted its complex conjugate−ℏ22m(φi*∇2φi−φi∇2φi*)=2iΓi∣φi∣2(5)where m is the electron mass. We defined the current density of the ith orbital J→i as usualJ→i=−iℏ2m(φi*∇→φi−φi∇→φi*)(6)

Then, we found∇→⋅J→i=−2ℏΓi∣φi∣2(7)

Integrating both sides over the volume V, which includes the molecule inside and using the Gauss theorem, we had∫Sn→⋅J→idS=−2ℏΓi∫V∣φi∣2dr→(8)where S is the surface of the volume V and n→ is a normal vector to the surface S. The ionization probability of the orbital i, w_i, is defined and is related to Γ_i aswi=∫Sn→⋅J→idS∫V∣φi∣2dr→=−2ℏΓi(9)

The electrons emitted to the continuum, in principle, should contribute to the self-consistent potential. However, if the ionization probability is very small, then the contribution of emitted electrons to the potential is negligible. Under the condition that the ionization probability is sufficiently small, we calculated ionization probability through the following two steps: First, we solved the static Kohn-Sham equation under the static external field, eEz. In this step, the tunnel ionization was forced to be prohibited by placing infinite wall potential outside the barrier. The problem becomes a usual static Kohn-Sham problem except the appearance of the external dipole field and infinite wall potential. In the second step, we calculated the Gamow state solution (31) for each Kohn-Sham orbital using the Hamiltonian obtained at the first step with the outgoing boundary condition.

The ionization probability is sensitive to the asymptotic behavior of the potential and the binding energies of the occupied orbitals. To appropriately incorporate the asymptotic behavior, we used the exchange-correlation potential, which takes account of the self-interaction correction. We adopted an approximate construction of the optimized effective potential including the self-interaction correction, which was proposed by Krieger, Li, and Iafrate (KLI) (32, 33). In this treatment, it has been shown that the ionization potentials of atoms and molecules approximately coincide with the energies of the highest occupied orbitals. The potential in this model also has a correct asymptotic behavior, −e²/r, for neutral molecules, where r is the radial distance from the center of the molecule.

To calculate the Gamow states, we needed to solve the static Kohn-Sham equation (Eq. 3) with the outgoing boundary condition. For systems without spherical symmetry, the treatment of the outgoing boundary condition is not simple. Instead of imposing the outgoing boundary condition explicitly, we used the absorbing boundary condition (ABC). The absorbing potential was placed in the spatial region outside a certain radius r = R with a thickness of ΔR. Outside the region of the absorbing potential, r > R + ΔR, the wave functions were set to vanish. If the absorbing potential works ideally, then there exist only the outgoing waves just inside the absorbing potential.

In the present calculation, we used the following spherical absorbing potential with a linear radial dependence−iW(r)={0(r<R)−iW0r−RR(R<r<R+ΔR)(10)where R was set beyond the barrier region. As listed in table S1, we determined the height W₀ (>0) and the thickness ΔR for which the electrons coming into the region r > R are absorbed completely. For the electrons with kinetic energy E_ele, the parameters satisfied the following condition (34)20Eele1/2ΔR8m<W0<110Eele3/28mΔR(11)

Using the absorbing potential, the Gamow state was obtained by solving the following Kohn-Sham equation{h[n(r→)]+Vext(r→)−iW(r→)}φi(r→)=(ϵiR+iΓi)φi(r→)(12)with the vanishing boundary condition, φ_i(r→) = 0 for |r→| > R + ΔR.

We also derived an expression for the ionization probability with the absorbing potential, W(r→). We started with Eq. 9. We took a sphere of radius R for the volume V and applied the Gauss theorem. Then, we obtained the following expression for Γ_iΓi=−∫r>Rdr→∣φi∣2W(r→)∫r<Rdr→∣φi∣2(13)

We assumed that the denominator of Eq. 13 is equal to a normalization constant, unity, so long as the ionization probability is sufficiently small. Then, the ionization probability of the orbital i is expressed as the following integral form using the absorbing potential W(r→) and the electron density niele(r→) = |φ_i(r→)|²wi=2ℏ∫r>Rdr→niele(r→)W(r→)(14)

We checked whether the ionization probabilities calculated with Eq. 14 coincide with those obtained from the imaginary part of the eigenvalue within the numerical error.

To express the orbital wave functions, we used the real-space grid method. This is a convenient representation in the Kohn-Sham theory, since the potential is almost local in the coordinate representation. To impose the ABC, one must treat large spatial region far outside the molecule. The number of grid points becomes substantially large. To save the computational effort, we reduced the number of grid points using the adaptive grid.

The adaptive grid was generated as follows. We introduced the following coordinate transformation from (x, y, z) to (u, v, w) for each Cartesian coordinatex=ku1+(k−1)(uasinh(u/a))l(15)y=kv1+(k−1)(vasinh(v/a))l(16)z=kw1+(k−1)(wasinh(w/a))l(17)where a, k, and l determine the property of the transformation. By this transformation, there holds x ~ u for small x (x << a) and x ~ ku for large x (x >> a). We then discretized uniformly the variables (u, v, w) with a constant interval of h. This produced a uniform grid for small x (x << a) and a coarse grid for large x (x >> a). The transformation parameters used in the calculations are also summarized in table S1. We examined carefully that the results are not sensitive to the choice of the parameters.

Simulation method of tunnel ionization probability for field-free MO

The numerical method described in the previous section cannot be applied to the ionization probability from the undeformed MOs in the field-free condition. Instead, we evaluated the ionization probability from the number of electrons outside the potential barrier formed by the external field. We defined the potential barrier based on the Kohn-Sham potential, V_KS = V_H + V_xc, where V_H is the Hartree potential and V_xc is the exchange-correlation potential, together with the potential by the electric field, V_ext = −eE→⋅r→. The position r_s was defined in the outer region of the barrier as the point where the total potential, V_KS + V_ext + V_ion, where V_ion is the pseudopotential of ions (35, 36), is equal to the binding energy of each MO. We estimated the ionization probability from the summation of the electron density outside the plane orthogonal to the electric field at r = r_s.

The present DFT calculation of the field-free MO was based on the same procedure used for the calculation of the molecule in the E-field. We used the static Kohn-Sham formalism in the DFT without an external dipole field, which is expressed in Eq. 3 with V_ext(r→) = 0. We solved the static Kohn-Sham equation with the optimized effective potential including the KLI self-interaction correction (32, 33) to obtain each Kohn-Sham orbital. Each orbital wave function was expressed using the same adaptive grid represented in Eqs. 15 to 17.

Derivation of theoretical RFPAD and integrated probability from angular-dependent ionization probability

We considered three vectors of the laser electric field E→, the momenta of the tunneled electron p→tunnel, and recoil ion p→ion. As an initial coordinate, we took a laboratory frame (X_L, Y_L, Z_L), where the Z_L axis is parallel to the propagation direction of the laser field and the E-field rotates in the X_LY_L plane with the azimuth angle Φ_L (fig. S4A). Assuming that the electron tunnels in the direction opposite to the E-field, the tunneling direction is defined only by the azimuth angle φ_L, which fulfills the condition of φ_L = Φ_L + 180° (fig. S4A). The orientation of the ethanol molecule is defined by the Euler angles (α, β, γ) relating the laboratory frame (X_L, Y_L, Z_L) shown in fig. S4A and the molecular-recoil frame (X_MR, Y_MR, Z_MR) shown in fig. S4C, where the Z_MR axis is parallel to p→ion and the C─C─O plane is on the X_MRZ_MR plane. Here, Φ_LR and φ_RFPAD are defined as the angles of the E-field and electron tunneling directions, respectively, from X_LR axis in the laboratory recoil frame (X_LR, Y_LR, Z_LR) (fig. S4B), which is obtained by the coordinate rotation represented with the direct cosine matrix A(α, 0°, 0°) (37). In the laboratory recoil frame, β and γ determine the three-dimensional molecular orientation without any ambiguity. Thus, the probability of the electron tunneling as a function of φ_RFPAD (=Φ_LR + 180°) is described asΩ10a′(φRFPAD)=∫02πdγ∫π2−Δβπ2+Δβdβ sin βPE(ΦLR)Pmol(β,γ)σ10a′(ΦLR,β,γ)(18)where P_E(Φ_LR) is the probability of the E-field direction, for which we substituted 1/(2π) because the laser pulse duration is long enough for the E-field to be equally distributed in the polarization plane; P_mol(β, γ) is the probability of the molecular orientation direction, for which we substituted 1/(4π) because the ethanol molecules are randomly oriented with respect to β and γ. The integral with respect to β is limited by the small out-of-plane angle Δβ (=10°) (fig. S2D).

The ionization probability σ_10a^′(Φ_LR, β, γ) was obtained from the ionization probability W_10a^′(Φ, Θ) in the molecular frame (X, Y, Z) (fig. S4D) by two successive rotational transformations. We prepared three-dimensional grid points (Φ_LR, β, γ) with equal intervals for the respective coordinates. The laboratory recoil frame was transformed to the molecular frame (fig. S4D) by the coordinate rotation represented with the two direct cosine matrices A(0°, β, γ) (fig. S4, B and C) and A(−90°, −90°, −η_CD2OH⁺) (fig. S4, C′ and D), in which η_CD2OH⁺ (=158.6°) is the angle between the recoil direction of CD₂OH⁺ and the C─C (X) axis. The ionization probability σ_10a′(Φ_LR, β, γ) at each grid point was equal to the probability W_10a′(Φ, Θ) at the transformed point in the molecular frame. We interpolated discrete data points of W_10a′(Φ, Θ) calculated by the DFT to obtain the probability at the transformed point. Thus, the theoretical RFPAD Ω_10a′(φ_RFPAD) was obtained by calculating the integral of Eq. 18.

To compare the theoretical RFPAD with the experimental one, the solid angle of the photoelectron emission was taken into account. We approximated the solid angle of photoelectron emission on the polarization plane to be the same as the out-of-plane distribution (fig. S3). The calculated RFPAD was convoluted with the Gaussian curve with its width of 31° (FWHM).

When we calculated the excitation fraction f3a″exc expressed in Eq. 2, we used the ion yields obtained from the time-of-flight mass spectrum. Here, we could not select the ionizing molecules in the limited range of the out-of-plane angle β. The probability P_10a′ was calculated by the following integralP10a′=18π2∫02πdΦLR∫02πdγ∫0πdβ sin βσ10a'(ΦLR,β,γ)(19)

Note that the analytical procedure mentioned in this section is valid not only for the 10a′ but also for all other orbitals.