Torsional effects in strong-ﬁeld ionization of molecules

We investigate how the ionization rate of molecules consisting of a pair of phenyl rings depends on the torsional dihedral angle between these two rings. We consider different orientations of the molecules with respect to the external ﬁeld direction and ﬁnd differences in the rates in biphenyl and substituted biphenyl. The evaluation of the rates is made possible by application of the integral representation of the weak-ﬁeld asymptotic theory for tunneling ionization. We report a strong variation of the ionization rates with dihedral angle with contributions from several orbitals and a large difference between biphenyl and substituted biphenyl related to the permanent dipole of the latter system.


I. INTRODUCTION
Ionization by a tunnelinglike process is key in much of strong-field and attosecond physics [1]. This is seen, for example, by its prominent role as the first step in the three-step model of high-order harmonic generation (HHG) [2], where the electron first tunnels through the effective barrier formed by the external field and the atomic or molecular potential, then propagates in the field, and finally recombines with the parent ion under the emission of high-frequency coherent radiation. Therefore, in connection with HHG, tunneling is a fundamental ingredient in characterizing the process itself, including its efficiency, but also central in the application of HHG for investigations of molecular properties such as molecular tomography [3] or HHG spectroscopy (see, e.g., Refs. [4,5]). Tunneling ionization is also the first step in laser-induced electron diffraction, where the electron that is released by a tunneling ionization process diffracts on the parent ion in a rescattering event and thereby records target structure information and evolution of the system typically in the vibrational degree of freedom from the time of ionization and until rescattering [6][7][8][9][10][11][12]. Another example is strongfield holography, where interpretation of the photoelectron momentum distribution in terms of an interference between a signal and reference electron wave, whose creation is influenced by tunneling, allows extraction of ultrafast dynamics [13][14][15][16], as well as in tunneling-ionization imaging of molecular photoexcitation, where the products of dissociation are recorded to detect changes occurring in the excited orbital [17]. Accordingly, in strong-field ionization by intense (near) infrared laser pulses there has been and still is focus on the sensitivity of the tunneling ionization dynamics on different molecular properties such as internuclear distance (see, e.g. Refs. [18][19][20][21]), molecular alignment and orientation (see, e.g, Refs. [22][23][24][25][26]), isotope effect [27,28], multielectron effects (see, e.g., Refs. [29][30][31][32]), and contribution from multiple orbitals (see, e.g., Ref. [33][34][35]). As will be discussed below, developments in the formulation and implementation of tunneling theory allow the precise prediction of tunneling rates for systems larger than have been possible to treat with accurate methods before. This progress means that new directions in the study of tunneling ionization of molecules open up and, equally important, it allows us to address cases where tunneling ionization is known to play a role but where its influence has been difficult to access due to the lack of accurate tunneling rates. Torsional effects in strong-field ionization is one such example and this is the topic of the present work.
We will study strong-field ionization of biphenyl and a substituted biphenyl molecule 3,5-difluoro-3 , 5dibromobiphenyl (DFDBrBPh). These systems consist of two phenyl rings connected by a C-C bond. Figures 1  and 6 show the geometries of biphenyl and DFDBrBPh, respectively. The two phenyl rings of the molecules are allowed to rotate around the central C-C bond axis, thereby forming different conformers. The internal rotation between the rings is described by the dihedral angle. The torsional motion has attracted theoretical and experimental interest. In particular, experimental works have shown, by alignment and Coulomb explosion techniques, that nonresonant femtosecond laser pulses can induce torsional motion via the interaction between the induced dipole moment and the laser field [36,37]. The use of this effective interaction to manipulate and coherently control vibration and torsion is termed dynamic Stark control [38]. It was discussed that laser-induced torsional motion can result in a transition from one enantiomer of a chiral molecule, such as DFDBrBPh FIG. 1. Orientation of the biphenyl molecule with respect to the electric field F . The axes x and z determine the MF of biphenyl, while the primed axis refers to the LF, in which the electric field F is directed along the z axis. The rotation from the LF to the MF is specified by the Euler angles (0, β, γ ). The MPA of the biphenyl molecule coincides with the z axis. The SMPA is illustrated in Fig. 2. [36,37] or 3,5-difluoro-3 , 5 -dibromo-4 -cyanobiphenyl [38,39], to its mirror form by light-induced deracemization, thereby producing enantiomeric excess. More recently, the possibility for laser-induced enantiomeric conversion of DFDBrBPh based on dynamic Stark control was investigated theoretically [40,41], including considerations of possibilities associated with invoking of artificial neural networks [42]. In the latter theoretical studies, the intensities considered were up to tens of TW/cm 2 , but the role of ionization in the dynamics was not considered. For the formulation of reliable control schemes for torsional motion, e.g., involving deracemization protocols, it is desirable to be able to include the effects of ionization.
In this work we investigate the effects of torsion on fieldinduced tunneling ionization of biphenyl and DFDBrBPh. We will evaluate the ionization rates in the presence of an electric field F by the weak-field asymptotic theory (WFAT) [43], which expresses the rate as an asymptotic expansion in the weak-field limit F → 0. The theory was originally developed in its tail representation relying on the asymptotic form of wave functions. However, implementation of the tail representation of the WFAT, or any other tunneling theory based on wave functions at large distance, for polyatomic molecules is technically extremely demanding. To resolve this issue, the integral representation of the WFAT was developed [44] (see also Ref. [45]). That approach allows evaluation of tunneling rates for polyatomic molecules and its implementation was discussed and demonstrated by treating benzene and naphthalene [46], and we will use that methodology in the following.
Adiabatic theory [47] justifies the use of static field ionization rates in modeling ionization dynamics induced by an intense laser pulse in the regime where the photon energy of the driving laser is much smaller than the binding energy of the ionizing electron. This regime is easily reached for the molecules considered in intense near-infrared or midinfrared laser pulses. Indeed, much of strong-field physics is conducted in this regime and static field ionization rates obtained by the WFAT have already been successfully applied in the interpretation of experimental data [4,5,17,28]. The possibility of determining accurate ionization rates for biphenyl and substituted biphenyl as a function of dihedral angle is therefore relevant for experimental studies resolving laserinduced torsional motion by strong-field Coulomb explosion [36][37][38][39]. To capture the effect of ionization on the torsional motion, one possible avenue is to describe the ionization as a loss process described by the tunneling rate of ionization. This rate can then be included in density matrix approaches or Monte Carlo wave packet approaches for the nuclear motion, as done, e.g., in modeling dissociative double ionization of hydrogen molecules [48].
The paper is organized as follows. The WFAT of tunneling ionization is briefly summarized in Sec. II. The results are discussed in Sec. III. Section IV summarizes the paper.

II. WEAK-FIELD ASYMPTOTIC THEORY
We consider tunneling ionization from a particular selected molecular orbital described within the Hartree-Fock (HF) approximation. The orbital ψ 0 (r) and its energy E 0 are defined by an effective one-electron Schrödinger equation (atomic units are used throughout) with a generally nonlocal potentialV (r) describing the interaction of an active electron with the parent molecular ion. The potential is an operator of the formV (r)ψ 0 (r) = V (r, r )ψ 0 (r )dr and includes exchange effects. The potential satisfiesV where Z is the total charge of the parent ion. Equation (1) is written in a molecular frame (MF). The ionizing electric field of strength F is assumed to point in the +z direction of the laboratory frame (LF). The orientation of the molecule with respect to the field is specified by the three Euler angles (α, β, γ ) defining a passive rotationR from the LF to the MF [49]. The rotation operatorR is a matrix such that v MF = Rv LF , where v LF and v MF are columns representing the same vector in the LF and MF, respectively. The ionization rate does not depend on α, so we set α = 0. Within the leading-order approximation of the WFAT, the ionization rate is given by [43] ( where is the structure factor defining the dependence of the rate on the orientation of the molecule with respect to the ionizing field and is the field factor defining its dependence on the field strength.
All quantities needed to implement these formulas are determined by properties of the unperturbed orbital, namely, is the same vector in the MF, g 00 (β, γ ) is a coefficient appearing in the asymptotic tail of ψ 0 (r) at r → ∞, and the subscript 00 specifies parabolic quantum numbers of the dominant ionization channel [43]. While the orbital energy E 0 and dipole μ can be calculated routinely by means of standard quantum chemistry packages, the calculation of the asymptotic coefficient g 00 (β, γ ) is a nontrivial part of the implementation. Within the integral-equation approach to the WFAT [44], this coefficient was expressed in terms of the orbital in an integral form [45] g 00 (β, γ ) = * HereV is a core potential [see Eq.
(2)] and 00 (r) is an analytically known function given by where M(a, b, x) is a confluent hypergeometric function [50]. We calculate g 00 (β, γ ) by evaluating the integral in Eq. (7) using a technique developed in Ref. [46]. The results reported below are obtained by an implementation of the above equations in the GAMESS program [51]. Before we consider the results, let us discuss the region of applicability of the leading-order WFAT. The leadingorder WFAT and the consideration of only the 00 channel has been shown to be accurate when the applied field F is smaller than a critical field F c , which is a boundary between the tunneling and the over-the-barrier regimes. The explicit expression for F c is given in Eq. (62) of Ref. [52] and in the present case F c reads F c κ 4 /(8|2Z − κ|). For the present systems, the values of κ 0.8 give a value of F c ∼ 0.04, which is much higher than the maximal field considered in the calculations, F = 0.02. Hence we are working in a regime where the application of the leading-order theory is accurate. First-order corrections to the leading-order WFAT include correction to the dominant 00 channel and contributions from the next-to-dominant channels 01. Such consideration includes an account for the field distortion of the ionizing orbitals, which is beyond the present scope (see Ref. [53] for an application to H 2 + molecules). The role of the next-todominant channels 01 close to nodal surfaces was discussed in Ref. [54] and shown to be down with a factor of F/2κ 2 , corresponding to ∼0.015 for the maximal field considered here. In the present mean-field HF approach, the total rate is the incoherent sum of contributions from several orbitals [30,55]. For biphenyl and substituted biphenyl, we will see that several orbitals contribute to the rate, and close to nodal surfaces of a particular orbital, the dominant 00 contributions from other orbitals not having a nodal surface will dominate over a possible 01 contribution from an orbital with a nodal surface. The only exception is in the results at a dihedral angle ϕ d = 0 • , where all the orbitals considered give a zero value for the structure factor. This case should be treated as explained in Ref. [54] [see Eq. (14) and Fig. 7 therein] by considering the 01 channel, which is beyond the leading-order WFAT. In this work, however, we restrict the treatment to the leading-order WFAT and consider only the 00 channel. Hence, the results in a very narrow window close to ϕ d = 0 • will be modified by contributions beyond the leading order. Such modifications are left for a future study. Finally, we note that many-electron WFAT has been formulated [30]. It is therefore possible to go beyond the present HF level of theory if correlated many-electron calculations can be performed. As an example we mention that small molecules like H 2 and LiH were considered in Ref. [56] using wave functions obtained from a configuration-interaction method. For biphenyl and substituted biphenyl the application of correlated methods poses a challenge for the future.

III. RESULTS AND DISCUSSION
Experimental work has shown that it is possible to threedimensionally adiabatically align biphenyl and substituted biphenyl molecules and to induce torsional motion [36][37][38][39]. It was shown in Ref. [37] that the lowest vibrational mode corresponds primarily to the torsional motion of the phenyl rings, i.e., dynamics in the dihedral angle between the rings. It was demonstrated that it is possible for a kick pulse to induce dynamics in the dihedral angle. This was achieved by Coulomb explosion with a strong laser pulse delayed with respect to the kick pulse. We now consider the ionization rate as a function of this dihedral angle for different orientations of biphenyl and DFDBrBPh with respect to the field. Inspired by the experiments, we focus on two directions of the external electric field with respect to the field: (i) along the most polarizable axis (MPA), which is along the C-C bond connecting the two rings (see Figs. 1 and 6), and (ii) along the second most polarizable axis (SMPA), which is perpendicular to the MPA (see the insets in Figs. 2 and 7). Since perfect alignment can be difficult to achieve experimentally, we will also consider the case of a finite small angle between the MPA and the direction of the field. As representative examples of field strengths in the tunneling ionization process, we choose the field strengths of F = 1.6 × 10 −2 and 2 × 10 −2 corresponding to intensities of ∼10 13 W/cm 2 . These field strengths provide sufficiently low rates of ionization to justify the use of the WFAT. See also the discussion about critical field strength in Sec. II. We note that the corresponding intensities are of the same order of magnitude as those considered in very recent theoretical work on laser-induced deracemization, neglecting effects of ionization [40][41][42].
The torsional frequency (calculated value 66 cm −1 ) and the typical frequencies of the C-C stretch (1300 cm −1 ) and C-H stretch (3000 cm −1 ) correspond to timescales of these motions of 500, 25, and 11 fs, respectively, so the nuclei can be assumed to be stationary within, say, 1 fs. For typical pulse durations of the ionizing pulse (10-25 fs for an 800-nm pulse), there will hence be time for the C-H bond to vibrate. The effects of vibration can be considered in tunneling theory, and in the absence of dissociation such vibrations do not significantly change the rates [57]. In the following calculations, we therefore only consider the variation in the dihedral angle as it would be induced by external laser pulses and otherwise fix the nuclear positions.
The reported results are obtained within the HF approximation using the uncontracted polarization-consistent upc-2 basis set [58] for both biphenyl and DFDBrBPh. In Ref. [46] the convergence of the structure factors, and hence rates, with the basis set was studied in detail. For He, Ne, Ar, and Kr, where exact HF results are available, the upc-2 results were found to be accurate within 0.1-3 %. Moreover, for CO 2 , where calculations with upc-3 and upc-4 basis sets were possible, the upc-2 basis set results fully captured the variation in the structure factor with alignment angle. We stress that the good convergence of the structure factor with the basis set is a very attractive property of the weak-field asymptotic theory formulated in the present integral representation approach.

A. Torsional and alignment effects in strong-field ionization of biphenyl
Biphenyl is a symmetric molecule with a zero total dipole moment. The dihedral angle between the planes of the two rings is denoted by ϕ d (0 • ϕ d 90 • ). The C-C bond axis passes through the planes of both phenyl rings and presents the MPA of the molecule. Figure 1 depicts a biphenyl molecule for vanishing dihedral angle (ϕ d = 0 • ) placed in the electric field such that the z axis of the MF is directed along the MPA and the x axis of the MF is parallel to the plane of one of the rings, while the electric field F points along the z axis of the LF specified above.
The difference in the total energy of various conformations of biphenyl with respect to its ϕ d = 0 conformer forms the dihedral potential of the molecule. The dihedral potential calculated within the HF approximation using the polarizationconsistent upc-2 basis set [58] as a function of ϕ d is shown in Fig. 2. The minimum of the dihedral potential is attained at a dihedral angle of ϕ d = 48 • . This differs from a more accurate theoretical prediction of ϕ d = 39 • obtained with density functional theory and the coupled cluster method [59]. At present, the theory is only implemented at the HF level of theory. Nevertheless, the leading-order WFAT should provide qualitatively correct predictions of the torsional effects in the strong-field ionization of molecules (see the discussion at the end of Sec. II).
The SMPA is located exactly between the two phenyl rings as indicated in Fig. 2. We consider tunneling ionization from three molecular orbitals which present the highest occupied molecular orbitals (HOMOs) of the ϕ d = 0 • angle conformer. We consider these three orbitals since their binding energies are of similar magnitude and since the order of magnitude of the ionization rate within the WFAT is determined mainly by a binding energy E 0 of the orbital via the field factor W 00 [Eq. (5)]. Although the energy ordering of the orbitals changes as the dihedral angle varies from 0 • and 90 • , we nevertheless denote them according to their ordering at ϕ d = 0 • : HOMO, HOMO-1, and HOMO-2, respectively (see Fig. 3). As an example of the interchange of energy position, we see from Fig. 3 that the orbital denoted HOMO-2 becomes the highest molecular orbital at ϕ d = 90 • . We note that the experimentally determined ionization energies vary from 0.293 to 0.323 [60]. This variation is of similar magnitude to the HOMO energies in Fig. 3.
The structure factors G 00 of Eq. (4) and ionization rates of Eq. (3) for the two field strengths considered are shown in Fig. 4  of relative orientation between the molecule and the direction of the external field, we consider the cases of the field directed along the MPA (β = 0 • ) and along the SMPA (β = 90 • ). In the case ϕ d = 0 • , the three orbitals considered have nodal surfaces in the plane of the molecule. Since the field when directed along both the MPA and the SMPA for ϕ d = 0 • lies in the plane of the rings, the corresponding structure factors G 00 and partial rates for the dominant ionization channel 00 [43] vanish for these three HOMOs. For the ϕ d = 90 • conformation, the nodal surfaces of the three orbitals considered contain the MPA of the conformer, so the corresponding structure factor G 00 turns to zero, when the field is directed along the MPA. Finally, the MPA of the HOMO and the SMPA of the HOMO-2 are completely contained in the corresponding nodal surfaces for any conformer of biphenyl, so the corresponding structure factors and rates equal zero for fields along the MPA and SMPA and are not shown in Fig. 4. All other conformations have generally nonzero structure factors and rates when the field is directed along the MPA or SMPA. Figures 4(b) and 4(c) illustrate the sensitivity to field strengths with typical magnitudes of ∼10 −7 a.u. for F = 0.016 and HOMO-2 (see Fig. 3) decreases with ϕ d , which means that field factor increases with dihedral angle and since the rate is the product of the field factor and the structure factor this results in a shifting of the maximum in G 00 to a larger ϕ d value in the rate. Similarly, the shift to smaller dihedral angles from the structure factor to the rates for the HOMO-1 with the field along the MPA can be explained by the decrease in the orbital energy with dihedral angle and the corresponding decrease in the field factor (see Fig. 3).
We now turn to consideration of the case when the MPA is not perfectly aligned with the field direction, i.e., the case of a finite β in Fig. 1. We consider the value β = 10 • (and γ = 0 • for definiteness) and calculate the corresponding structure factors G 00 (β = 10 • , γ = 0 • ) and rates (β = 10 • , γ = 0 • ) for the three HOMOs considered. The results for the two values of the field strength are shown in Fig. 5. It can be seen that the orders of magnitude of the ionization rates from the HOMO-1 and HOMO-2 corresponding to different field strengths are similar to those in Fig. 4. For β = 10 • the electric field does not lie in the nodal surfaces of the HOMO for conformations with nonzero ϕ d values. Hence, unlike the β = 0 • case, the structure factors and rates of the HOMO are generally nonzero (black curves in Fig. 5). For the HOMO, the combination of an increasing structure factor and a decreasing energy (Fig. 3) results in an ionization rate attaining its maxima at around 48 • for the field F = 0.016 a.u. and around 57 • for the field F = 0.02 a.u. We also see shifts in the peak positions between the structure factors and the rates for the HOMO-1 and HOMO-2 and the reason is also in this case the behavior of the field factor with binding energy.

B. Torsional and orientational effects in strong-field ionization of DFDBrBPh
In this section we consider the case of the substituted biphenyl molecule DFDBrBPh. The DFDBrBPh molecule  is an axially chiral biphenyl derivative with two hydrogens replaced with two bromine atoms at one phenyl ring and two hydrogens replaced with two fluorine atoms at another ring. Like for biphenyl, the C-C bond axis connecting the two phenyl rings presents the MPA. Figure 6 depicts the molecule in the MF and the electric field F points along the z axis of the LF. The z axis of the MF is parallel to the MPA of DFDBrBPh, while the x axis is parallel to one of the phenyl rings. The two bromine atoms are located in the positive direction of the z axis and the two fluorine atoms are located in its negative direction. The dihedral potential, the dihedral angle, and the SMPA are shown in Fig. 7. The value of the dihedral angle ϕ d = ϕ Br − ϕ F determines the conformation of the molecule. For the conformer of ϕ d = 39 • the SMPA is perpendicular to the MPA and passes between the two phenyl rings 11 • away from the Br-phenyl ring and 28 • away from the F-phenyl ring [36]. As for biphenyl, we consider three HOMOs due to the mixing of their relative energy positions as a function of ϕ d . The dependence of the orbital energies is shown in Fig. 8. As was the case for biphenyl, we see that the ordering of the orbital energies is a function of the dihedral angle.
The structure factors and ionization rates for the cases with the external field along the MPA and SMPA are shown in Fig. 9. For the ϕ d = 0 • conformation, the plane of the rings coincides with the nodal surfaces of the three orbitals considered. Accordingly, when the field is along the MPA or the SMPA, the structure factors, and therefore the rates, vanish. For other conformations, the different orbitals give quite different contributions to ionization. For example, the contribution of the HOMO-1 to ionization along the MPA exceeds that of the HOMO-2 by approximately six orders of magnitude, which cannot be explained alone based on differences in orbital energies and corresponding field factors. Rather, the large magnitude of the structure factors and ionization rates (β = 0 • , γ = 0 • ) of the HOMO-1 is explained by the unusually large dipole moment of about 4.8 a.u. of this orbital. The dipole of this orbital points from the bromine to the fluorine end and leads, in the presence of the field, to an exponential increase of the structure factor [see Eq. (4)], associated with a field-induced upward linear Stark shift of the orbital energy (see Ref. [43] for a thorough discussion of dipole effects in tunneling ionization of molecules). Compared with the HOMO-1 results, the structure factor and ionization rates of the HOMO-2 turn out to be relatively small. As was the case for biphenyl, we see changes in the behavior of the curves going from structure factors to rates. These changes can again be understood by the dependence of the field factor on the orbital energy.
As for biphenyl, we also consider for DFDBrBPh a case where the field is along neither the MPA nor the SMPA. To this end, Fig. 10 shows structure factors and ionization rates when the electric field points in the direction corresponding to the Euler angles β = 10 • and γ = 0 • . In this case the contribution of the HOMO-1 again prevails because of the large dipole of the orbital. For conformations with ϕ d around 90 • , however, the dipole of the HOMO sharply increases from 0.3 to 3.7, which results in large values of the structure factors (black line in Fig. 10) and a corresponding increase in the rates from the HOMO the neighborhood of a point ϕ d = 90 • . We emphasize that DFDBrBPh shows a different behavior of ionization than biphenyl due to large nonzero dipole moments 023018-7 of the three HOMOs investigated, e.g., for the field strength F = 0.02 the ionization rate of DFDBrBPh reaches 10 −3 a.u., while for biphenyl it is just about 10 −5 a.u. even though the field-free orbital energies are similar in the two systems.

IV. CONCLUSION AND OUTLOOK
With this work we have initiated a study of torsional effects in strong-field ionization of molecules. We have investigated biphenyl and substituted biphenyl as illustrative examples motivated by the fact that torsional motion has been elucidated in alignment resolved experiments in such systems [36][37][38][39] and that these molecules have been explored theoretically in connection with laser-based protocols for laser-induced deracemization [40][41][42]. The determination of the molecular property related to the strong-field tunneling ionization process, the structure factor of the WFAT [43], was made possible by application of the newly implemented methodology [46] for the integral representation of WFAT [44] (see also Ref. [45]). Using this approach, we obtained accurate tunneling rates and investigated the dependence of these on the dihedral angle between the two phenyl rings, as well as on the orientation of the molecules with respect to the ionizing field. The results showed almost up to an order of magnitude variation in the ionization rates with the dihedral angle. Furthermore, we found that biphenyl and the substituted biphenyl considered, DFDBrBPh, have three highlying occupied orbitals with similar energies. These three orbitals all contribute significantly to the ionization process. The relative contributions from the orbitals vary with the dihedral angle in a way determined by the variation in binding energy, field factor, and structure factor. In this connection, it was found that a permanent dipole, such as that present in DFDBrBPh, can have significant effects on the ionization dynamics. In addition, this dipole effect was found to be strongly dependent on the dihedral angle.
The accurate determination of strong-field tunneling rates for biphenyl and axial chiral systems like DFDBrBPh raises some interesting possibilities for the future. To name a few, such rates should make absolute comparisons with measurements exploring torsional motion by Coulomb explosion feasible, i.e., the theoretical calculation of the variation in the rate with dihedral angle should allow one to capture the variations in measured absolute Coulomb explosion yields with torsional angle. Moreover, these rates will allow one to consider effects of electron ionization dynamics in efforts related to laser control of torsional motion and laser manipulation of axial chiral systems. Finally, these rates will be valuable input in efforts aiming at resolving nuclear motion with, e.g., high-order harmonic generation or laser-induced electron diffraction approaches in larger systems than have been considered so far.