Long-Lasting Orientation of Symmetric-top Molecules Excited by Two-Color Femtosecond Pulses

Impulsive orientation of symmetric-top molecules excited by two-color femtosecond pulses is considered. In addition to the well-known transient orientation appearing immediately after the pulse and then reemerging periodically due to quantum revivals, we report the phenomenon of field-free long-lasting orientation. Long-lasting means that the time averaged orientation remains non-zero until destroyed by other physical effects, e.g. intermolecular collisions. The effect is caused by the combined action of the field-polarizability and field-hyperpolarizability interactions. The dependence of degree of long-lasting orientation on temperature and pulse's parameters is considered. The effect can be measured by means of second (or higher-order) harmonic generation, and may be used to control the deflection of molecules traveling through inhomogeneous electrostatic fields.


I. INTRODUCTION
Over the years, diverse optical methods have been developed to align and orient molecules of varying complexity and many applications related to studies of molecular and photon-induced processes are based on the ability to control the absolute orientation of the molecules. For reviews, see [1][2][3][4][5][6].
Here, we investigate the orientation dynamics of symmetric-top molecules excited by single two-color femtosecond laser pulses. In addition to the wellknown transient orientation appearing immediately as a response to the laser excitation, we predict the existence of long-lasting orientation. Long-lasting means that the time-averaged orientation remains non-zero, within the model, indefinitely or until destroyed by additional physical effects, e.g. by collisions. The long-lasting orientation induced by a two-color pulse has an intricate dependence on both the molecular polarizability and hyperpolarizability. Related effects have been recently observed in chiral molecules excited by one-color laser pulses with twisted polarization [30,32] and investigated in non-linear molecules excited by THz pulses [23,33].
The paper is organized as follows. In the next section, we describe our numerical approaches for simulating the laser-driven molecular rotational dynamics. In Sec. III, we present the long-lasting orientation, which is the main result of this work. Section IV is devoted to a qualitative analysis of the effect, and a derivation of the approximate classical formula for the degree of longlasting orientation. Additional results are presented in Sec. V.

II. NUMERICAL METHODS
In this work, the rotational dynamics of symmetric-top molecules is treated within the rigid rotor approximation. We performed both classical and quantum mechanical simulations of molecular rotation driven by two-color laser fields.
This section outlines the theoretical approaches used in both cases.

A. Classical simulation
In the classical limit, the rotational dynamics of a single rigid top is described by Euler's equations [50] IΩ = (IΩ) × Ω + T, where I = diag(I a , I b , I c ) is the moment of inertia tensor, Ω = (Ω a , Ω b , Ω c ) is the angular velocity, and T = (T a , T b , T c ) is the external torque resulting from the interaction between field-induced dipole moment and the electric field. All the quantities in Eq.
(1) are expressed in the rotating molecular frame of reference, equipped with a basis set including the three principal axes of inertia, a, b, and c.
In the laboratory frame of reference, the electric field of a two-color laser pulse is defined by where the two terms correspond to the FW and its SH, respectively. ω is the carrier frequency of the FW field, ϕ is the relative phase of the second harmonic, ε n (t) = ε n,0 exp[−2 ln 2 (t/σ n ) 2 ], n = 1, 2, is the field's envelope with ε n,0 as the peak amplitude, and σ n is the full width at half maximum (FWHM) of the laser pulse intensity profile. The polarization direction of the SH field is given by e SH = cos(φ SH )e Z + sin(φ SH )e X , where φ SH is its angle with respect to the Z axis, e Z and e X are the unit vectors along laboratory Z and X axes, respectively. The electric field in the molecular frame of reference can be expressed as where e 1 = Qe Z and e 2 = Qe SH are the unit vectors expressed in the molecular frame of reference. Q is a 3 × 3 time-dependent orthogonal matrix relating the laboratory and the molecular frames of reference. It is parametrized by a quaternion, q which has an equation of motionq = qΩ/2, with Ω = (0, Ω) being a pure quaternion [51,52]. Considering laser-polarizability and laser-hyperpolarizabiltiy interactions, the torque induced by a two-color field has two contributions T = T α + T β , where [53] T α = (αE × E) = ε 2 1 2 (αe 1 ) × e 1 + ε 2 2 2 (αe 2 ) × e 2 , (4) Here, the overline (· · · ) represents averaging over the optical cycle, α and β are the polarizability and hyperpolarizability tensors, respectively. ǫ ijk is the Levi-Civita symbol, β mnj is a component of the hyperpolarizability tensor, e 1m and e 2m are the components of the FW and SH fields, respectively.
To simulate the behavior of an ensemble of noninteracting molecules, we use the Monte Carlo approach. For each molecule, the Euler's equations [Eq. (1)] with the torques in Eqs. (4) and (5) are solved numerically using the standard fourth order Runge-Kutta algorithm. In the simulations we used ensembles consisting of N ≫ 1 molecules. The initial uniform random quaternions, representing isotropically distributed molecules, were generated using the recipe from [54]. Initial angular velocities are distributed according to the Boltzmann distribution, where i = a, b, c, T is the temperature and k B is the Boltzmann constant.

B. Quantum simulation
The Hamiltonian describing the rotational degrees of freedom of a molecule and the molecular polarizability and hyperpolarizability couplings to external timedependent fields can be written as where H r is the field-free Hamiltonian [55], and H int (t) is the molecule-field interaction potential, with two Here E i , α ij , and β ijk are the components of the field vector, polarizability tensor α, and hyperpolarizability tensor β, respectively. Since the optical carrier frequency of the laser fields, ω [see Eq.
(2)], is several orders of magnitude larger than a typical rotational frequency of small molecules, the energy contribution due to the interaction with the molecular permanent dipole, µ, −µ · E(t) is negligible. We use the eigenstates of H r , |JKM , describing the field-free motion of quantum symmetric-top [55], as the basis set in our numerical simulations. The three quantum numbers are J, K and M , where J is the total angular momentum, while K and M are its projections on the molecular a axis and the laboratory-fixed Z axis, respectively. The time-dependent Schrödinger equation i ∂ t |Ψ(t) = H(t)|Ψ(t) is solved by numerical exponentiation of the Hamiltonian matrix (see Expokit [57]) with the initial state being one of the fieldfree eigenstates, |Ψ(t = 0) = |JKM . The degree of molecular orientation is derived by calculating the induced polarization, the expectation value of the dipole projection. The polarization along each of the axes in the laboratory-fixed frame of reference is given by where e i represents one of the unit vectors e X , e Y , e Z . Thermal effects are accounted for by computing the incoherent average of the time-dependent polarizations obtained for the various initial states |JKM . The relative weight of each of the projections µ (J,K,M) i (t) is defined by the Boltzmann distribution, where Z = J,K,M ǫ K exp (−ε J,K,M /k B T ) is the partiton function, and ε J,K,M is the energy/eigenvalue corresponding to |JKM state. For molecules with two or more identical atoms, an additional statistical factor ǫ K must be included in the distribution [58]. For the case of methyl fluoride (CH 3 F) molecule considered in this work, ǫ K is given by where I H = 1/2. In our simulations, the basis set included all the states with J ≤ 30. For our sample molecule, CH 3 F, at initial temperature of T = 5 K, this means that initial states with J ≤ 8 were included. Additional details about the numerical simulations, including the matrix elements of the interaction Hamiltonian H int , can be found in Appendix A.

III. LONG-LASTING ORIENTATION
We continue to consider the methyl fluoride (CH 3 F), as an example for a symmetric-top molecule. The molecule is excited by a two-color pulse in which the polarizations of the FW and SH are parallel and along Z axis [φ SH = 0, see Eq. (2)]. In addition, here we set the relative phase between them to be zero (ϕ = 0). Later on we discuss what changes when this phase changes. Table I summarizes the molecular properties of CH 3 F. Moments of inertia, dipole moment, and polarizability tensor components are taken from NIST, where they were computed within the density functional theory (DFT, method CAM-B3LYP/aug-cc-pVTZ) [59]. The hyperpolarizability values are literature values taken from [60]. Figure 1 shows the projection of the dipole moment along the laboratory Z axis, µ Z , calculated classically and quantum mechanically (see Methods Section II). In the classical case, the angle brackets · · · denote ensemble average, that is the average of the dipole projections of N = 10 8 molecules, initially isotropically distributed in space and having random angular velocities [see Eq. (6)].
In the quantum case, · · · denotes incoherent average of initially populated rotational states [see Eq. (9)]. Note that the averages µ X and µ Y are zero. Here, the initial temperature is T = 5 K, the peak intensities of the FW and SH fields are I FW = 8 × 10 13 W/cm 2 and I SH = 3 × 10 13 W/cm 2 , respectively, and the duration (FWHM) of the pulses are σ 1 = σ 2 = 120 fs [see Eq.
In the case of symmetric-top molecules considered here, we observe long-lasting (persistent) orientation, a previously unreported phenomenon in two-color orientation schemes. The inset in Fig. 1 demonstrates that after the initial oscillations are washed out, the classical polarization/degree of orientation attains a constant, nonzero value. In the quantum case too, despite its being partially masked by the revivals, the sliding time average of the signal is approximately constant and it persists indefinitely within the adopted model. This long-lasting orientation is one of the main results of this work.
Several comments are in order. Additional physical effects can distort the long-term field-free picture of identical periodically appearing revivals seen in Fig. 1. These include the centrifugal distortion and the radiation Z-projection of the dipole moment, µZ and the orientation factor, cos(θµZ) ≡ µZ /µ as a function of time for CH3F molecule at initial rotational temperature T = 5 K. Here µ is the magnitude of the dipole moment and θµZ denotes the angle between the dipole moment and laboratory Z axis. The solid blue and dotted red lines represent the results of quantum and classical simulations, respectively. The solid green line is the time average defined by µZ (t) = (∆t) −1 t+∆t/2 t−∆t/2 dt ′ µZ (t ′ ), where ∆t = 19.6 ps. The inset shows a magnified portion of the signals. emission due to rapidly rotating molecular permanent dipole moment. Dephasing of the rotational states caused by the centrifugal distortion leads to the eventual decay of the revivals' peaks [21,22]. Nevertheless, the average dipole remains almost unchanged (see [23]). The radiative emission results in the gradual decrease of the rotational energy [21,22].
However, for a rarefied molecular gas, the estimated relative energy loss during a single revival is very small. The proper description of the behavior on an even longer timescale (nanoseconds), requires the inclusion of collisions and fine structure effects [66,67], which is beyond the scope of the current work. Furthermore, it should be noted that higher laser pulse intensities lead to higher degree of orientation, but when the intensity is high enough for molecular ionization, another effect kicks in, namely orientation mechanism due to selective molecular ionization of molecules with specific orientation [43,44]. Considering this kind of orientation is also beyond the scope of this work.

IV. LONG-LASTING ORIENTATION -A QUALITATIVE DESCRIPTION
An explicit form of the interaction potential [Eq. (7)] can be obtained by expressing the electric field vector in the rotating molecular frame of reference. For the sake of the current discussion, this can be done conveniently by using an orthogonal rotation matrix parameterized by the three Euler angles, R(φ, θ, χ). We use the definition convention adopted in [55], according to which, φ and θ are the standard azimuth and polar angles defining the orientation of the molecular frame z axis, and χ is the additional rotation angle about z axis. The basis set in the rotating molecular frame of reference consists of the three principal axes of inertia, a, b, c. For molecules belonging to the C 3v symmetry group (such as CH 3 F), there are three non-zero polarizability components (two of them are equal), and 11 non-zero hyperpolarizability components (three of which are independent) [68]. For definiteness, we associate the axis of the three-fold rotational symmetry with the most polarizable molecular principal axis a (z axis in the rotating frame), having the smallest moment of inertia, I a . In this case, the non-zero polarizability elements are α aa > α bb = α cc , and the independent hyperpolarizability elements are β aaa , β abb = β acc , β bbb = −β bcc . The other non-zero hyperpolarizability elements are obtained by permuting the indices of the independent elements [68]. The parameters of the CH 3 F molecule are listed in Table I.
We consider the case of a two-color pulse in which both the FW and SH are polarized along Z axis. The interaction potential is obtained by carrying out the summation in Eq. (7) (where all the quantities are expressed in the basis of principal axes of inertia) and we average over the optical cycle. The resulting potential Angle dependence of the potential V ⋆ β (θ), cos(ϕ)[3b sin 2 (θ) cos(θ) + cos 3 (θ)], see Eq. (13), for different relative phases, ϕ. Here, b = β abb /βaaa = 2/3. has two contributions, To facilitate the qualitative discussion in this section, we let β bbb = β bcc = 0. These elements of the hyperpolarizability tensor are the smallest (see Table  I), and their omission does not affect the qualitative features of the discussed phenomena. Thus, the hyperpolarizability interaction becomes And V = V α + V ⋆ β is a function of a single variable θ-the polar angle between the symmetry axis of the molecule (a axis) and the laboratory Z axis (axis of laser polarization).

A. Approximate classical formula
In the case of weak excitation, we can derive an approximate classical formula for the degree of longlasting orientation. Since the long-lasting orientation manifests itself under field-free conditions, we begin by considering the free motion of a single classical symmetric top. The free motion of the unit vector a, pointing along the rotational symmetry axis of the molecule, is given by a simple vectorial differential equationȧ = (L/I) × a.
Here, L is the conserved angular momentum vector, I is the moment of inertia along the orthogonal axes b and c (I a < I b = I c ≡ I). The solution of this equation is given by where L is the magnitude of angular momentum and a(0) is vector a at t = 0. The above equation describes precession of a around L at a rate L/I (see Fig. 3). In the special case of a linear molecule, L a = L · a(0) = 0, so that Eq. (14) reduces to Equation (15) describes a uniform rotation of a in a plane perpendicular to the angular momentum vector L.
The degree of long-lasting orientation (see the inset of Fig. 1) can be obtained by considering the ensemble average projection of the molecular axis a on the laboratory Z axis, a Z = e Z · a, and then evaluating its time average.
where t = 0 defines the end of the two-color pulse (when the free motion begins). Note that for CH 3 F the molecular dipole, µ points along −a (see Table I).
Next, we exchange the order of the ensemble and time averaging. The time average of a Z is obtained from Eq. (14) and it reads where L Z = e Z · L, L a = a · L, and subindex f denotes that all the quantities are taken after the pulse. With the potential in Eqs. (11) and (13), both φ and χ are cyclic coordinates. Therefore, the canonically conjugate angular momenta L Z and L a are conserved. As a consequence, Eq. (17) becomes where L Z and L a are taken before the pulse. At this stage, we can conclude that the long-lasting orientation is strictly zero when the initial temperature is zero and/or in the limit of a linear rotor. In the first case, L Z = L a = 0, while in the second case L a = 0, because I a = 0 for linear molecules.
For the ensemble averaging, it is advantageous to express all the quantities in the basis of principal axes of inertia. The magnitude of the angular momentum after the pulse, L 2 f , is given by where L a , L b , L c are the values before the pulse, while δL b and δL c are the changes in angular momentum components due to laser excitation. Explicit expressions for δL b and δL c can be obtained using the impulsive approximation. In this approximation, we assume that the duration of the two-color pulse is much shorter than the typical period of molecular rotation, such that the molecular orientation remains unchanged during the pulse. Using this approximation and the Euler-Lagrange equations, we derive the explicit expressions for δL b and δL c (the details are summarized in Appendix B), where f (θ) = P 1 sin(2θ) +P 2 sin(θ)[(3 cos(2θ) + 1)β abb − 2 cos 2 (θ)β aaa ], (22) and Here σ = σ 1 = σ 2 , ε 1,0 and ε 2,0 are the peak amplitudes of the FW and SH, respectively. L Z is expressed in terms of the molecular frame components L a,b,c using the rotation matrix R(φ, θ, χ), such that L Z = − sin(θ) cos(χ)L b +sin(θ) sin(χ)L c +L a cos(θ). (25) Finally, we carry out the ensemble average where w = I/I a , andĨ L (w) is a monotonic function of w > 1. In the limit of a linear molecule (w → ∞), I L (w) → 0 [see Eq. (C10) and Fig. 7]. For the polarizability interaction alone, f (θ) = P 1 sin(2θ), and f 2n (θ) sin(2θ) dθ = 0. For the hyperpolarizability interaction alone, f (θ) = P 2 sin(θ)[(3 cos(2θ) + 1)β abb − 2 cos 2 (θ)β aaa ] which is a symmetric function (about θ = π/2), and therefore f n (θ) sin(2θ) dθ = 0 for all n. Only when both polarizability and hyperpolarizability interactions are included, a Z = 0. In this case, where P 1 and P 2 are given by Eqs. (23) and (24), and I L (w) is given by Eq. (C10). The details of the derivation of Eq. (28) are summarized in Appendix C. According to Eq. (28), the degree of long-lasting orientation scales as σ 2 /T . In Fig. 4, we compare  Fig. 1. Cases of fully time-dependent field (solid blue) and using impulsive approximation (dashed red) are compared. The dotted green line is obtained using the approximate formula in Eq. (28) with µZ = −µ aZ (dipole moment points againts a axis, see Table I).
the temperature (panel A) and pulse duration (panel B) dependencies of the long-lasting orientation obtained using the approximate formula in Eq. (28) with the numerical results obtained by evaluating the formula in Eq.
(18) using the Monte Carlo approach as described in the Methods Section II (using the impulsive approximation, see Appendix B).
There is a good agreement between the numerical results and the results obtained using the approximate formula, especially at higher temperatures (higher initial angular momenta), where the assumption |δL b /L b |, |δL c /L c | ≪ 1 is well satisfied. The pulse duration dependence shows the connection to the energy gained by the molecule from the laser pulse. In the limit of weak excitation (low pulse intensity and/or high temperature), the approximate formula also reveals the more involved dependence on the fields' amplitudes, according to Eqs. (23) and (24).
The hyperpolarizability part of the interaction potential, V ⋆ β (θ) [see Eq. (13)], is an asymmetric function of θ (about θ = π/2, see Fig. 2), similar to the orienting potential, which is proportional to − cos(θ), due to a single THz pulse interacting with the molecular dipole, µ. As we show here and as it was shown in [23], excitation by such orienting potentials results in transient orientation followed by residual long-lasting orientation.
Despite the similarity, the mechanisms behind the long-lasting orientation induced by a femtosecond twocolor and a picosecond THz pulse are not the same. It was shown in [23] that in the limit of vanishing THz pulse duration, the induced long-lasting orientation tends to zero. In other words, a δ-kick by a purely orienting potential doesn't lead to long-lasting orientation. In contrast, here we show that a δ-kick by a combined, aligning and orienting, potentials results in a long-lasting orientation [see the discussion under Eq. (27), also see   Figure 5 depicts the long-lasting orientation of the dipole moment as a function of temperature for the case of collinearly polarized two-color pulse. Due to the short pulse duration (σ = 120 fs), the results of the fully time-dependent simulation (solid blue line) are well reproduced using the impulsive approximation (dashed red line). The impulsive approximation is described in Appendix B. As mentioned in Sec. IV, the long-lasting orientation vanishes at T = 0 K, [see Eq. (18)]. At high temperatures, the long-lasting orientation decreases as ∝ T −1 . Therefore, there should be an optimal temperature for which the long-lasting orientation is maximal. As is shown in Fig. 5, for the field parameters used here, the optimal temperature is T ≈ 20 K. In the derivation leading to the approximate formula in Eq. (28), we assumed β bbb = β bcc = 0 . Nevertheless, the results obtained using Eq. (28) qualitatively agree with the numerical results at higher temperatures as well (dotted green line).

V. TEMPERATURE AND POLARIZATION DEPENDENCE OF THE LONG-LASTING ORIENTATION
As an additional example, we consider the case of a cross-polarized two-color pulse in which the polarizations of the FW and SH are along Z and X axes, respectively. Figure 6 shows the dipole signal along the laboratory X axis. Note that the Y and Z-projections of the dipole moment, µ Y,Z are exactly zero. In this case, a similar transient dipole response along the polarization direction of SH can be seen. On the long time scale, it is followed by the long-lasting orientation. Notice that for the field parameters used here, the achieved degree of both transient and long-lasting orientation is higher in the case of cross-polarized FW and SH compared with Fig. 1 (also see [46]).

VI. CONCLUSIONS
We have theoretically demonstrated a new phenomenon of long-lasting (persistent) orientation of symmetric-top molecules excited by a single two-color femtosecond pulse. The residual orientation was shown to last indefinitely (within the adopted model), or until destroyed by other physical effects, e.g. intermolecular collisions.
We derived an approximate classical expression revealing several qualitative features of the phenomenon, including the scaling with temperature, pulse duration, and other field and molecular parameters. The predictions of the formula are in full agreement with the results of numerical simulations in the limit of weak excitation. A quick check for different polarizations showed that in the case of cross-polarized FW and SH, the achieved degree of both transient and long-lasting orientation may be higher than in the case of parallel configuration. Further and careful optimization of the parameters of the two-color pulse may give rise to even higher degrees of long-lasting orientation. Moreover, it has been demonstrated before that the degree of molecular alignment/orientation can be enhanced when a sequence of several laser pulses is used instead of a single pulse [17,[69][70][71][72][73][74], and a similar approach may be beneficial for increasing the degree of the long-lasting orientation. The orientation may be measured with the help of second (or higher-order) harmonic generation [42] and may be utilized in deflection experiments using inhomogeneous electrostatic fields [75,76].

Appendix A: Quantum-mechanical simulation
We use free symmetric-top wave functions |JKM , the eigenfunctions of kinetic energy Hamiltonian, H r [55], as the basis set for numerical calculations. In this basis, the rotational kinetic energy of a prolate top is given by [55] JKM |H r |JKM = CJ(J + 1) + (A − C)K 2 , (A1) where the rotational constants are C = 2 /2I b = 2 /2I c , A = 2 /2I a with A > C. Using a spherical basis, the transformation between tensors in laboratory-and molecule-fixed reference frames is given by [55] T (k) where D k pq * (R) represents the complex conjugate of the Wigner D-matrix, and R is the set of three Euler angles. T (k) p and T (k) q denote the spherical tensors of rank k in the laboratory-and molecule-fixed reference frames, respectively. Using the transformation in Eq. (A2), V α [see Eq. (7)] (averaged over the optical cycle) can be expressed as [31,77] where γ is the angle between the instantaneous field polarization and the Z axis (according to Eq. (2), γ 1 = 0 and γ 2 = φ SH ). The repeated index n = −2, · · · , 2 is summed over. The spherical tensor elements of the molecular polarizability are given by In addition, V β in Eq. (7) (averaged over the optical cycle) can be recast into lab,m and β (j,j12) lab,m are the spherical tensor elements of the electric field and the molecular hyperpolarizability in the laboratory reference frame, respectively.
The spherical tensor A (j3,j12) lab can be constructed by performing dyadic product of A (1) twice [78,79], namely where 1m 1 1m 2 |j 12 m 12 and j 12 m 12 1m 3 |jm represent the Clebsch-Gordan coefficients. A (1) is the electric field's spherical tensor of rank one and its components are given by A (1) For the two-color field defined in Eq. (2), the nonzero components of A while the hyperpolarizability in the laboratory reference frame is given by β The matrix elements of the molecule-field interaction potential, JKM |H int |J ′ K ′ M ′ can be obtained with the help of the relation [55] JKM |D s where e (1) The components of spherical tensor representing the dipole moment in the molecular reference frame are defined by by µ Appendix B: Derivation of formulas in Eqs. (20,21) The Lagrangian of the classical symmetric-top molecule excited by a two-color pulse reads where the set of canonical variables includes the three Euler angles and their time derivatives, (θ, φ, χ;θ,φ,χ).
The interaction potential V (θ, t) is given by Eqs. (11) and (13). In terms of the canonical variables, the components of angular momentum (in the basis of molecular principal axes of inertia) are given by [55] L b = I −φ sin(θ) cos(χ) +θ sin(χ) , L c = I φ sin(θ) sin(χ) +θ cos(χ) , such that the Lagrangian in Eq. (B1) becomes Assuming the angles remain constant during the impulsive excitation, the changes of angular momentum components are given by δL b = I −δφ sin(θ) cos(χ) + δθ sin(χ) , δL c = I δφ sin(θ) sin(χ) + δθ cos(χ) , δL a = I a δφ cos(θ) + δχ , where δg = g(t f ) − g(t i ) denotes the difference between a quantity after (t f ) and before (t i ) the pulse. once and solving the resulting system of equations for δθ, δφ, δχ. During the integration, the angles are assumed to be constant, so that integrals such as such that δL a = 0.

(B7)
Appendix C: Derivation of formula in Eq. (28) In this section, we carry out the intermediate steps in the derivation of Eqs. (27) and (28). After dφ integration, Eq. (26) becomes Next, we expand 1/L 2 f in powers of f (θ) and carry out dχ integration. Only terms proportional to even powers of f (θ) survive the integration. The first non-vanishing term is proportional to f 2 (θ), such that where The temperature dependence is isolated in I L and can be exposed by change of variables. For this, we substitute where w = I/I a . Changing to spherical coordinates, and integrating dφ L produces Next, we introduce where I Z (w) = l 2 exp −l 2 [1 + cos 2 (θ L )(w − 1)] Thus, wherẽ 4(w − 1) 5/2 , (C10) see Fig. 7. When both polarizability and hyperpolarizability interactions are included, f (θ) is given by Eq. (22), and f 2 (θ) sin 2θ dθ = 64 105 P 1 P 2 (2β abb − 3β aaa ) .

CONFLICT OF INTEREST STATEMENT
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
All the authors participated in formulating the problem and initiating this study.
L.X. and I.T. equally contributed to the calculations and numerical simulations. All the authors participated in analyzing the results and writing the manuscript. Y.P. and I.A. supervised and guided the work.