Collisions of excited Na atoms with N2, C2H2 and N2O molecules. Surface-hopping calculation of the Na(3p) fine structure population

The results of optical collision experiments: fine structure population ratios and polarization alignment angles obtained from scattering of Na atoms which are optically 3s → 3p excited during the collision with N2, C2H2 and N2O molecules, have been simulated numerically and analysed using quantum chemically computed potential energy surfaces, classical trajectory calculations and the semiclassical surface hopping method combined with a specially developed diabatization procedure.


Optical collisions are scattering processes where the projectile is optically excited during its interaction with the target.

This paper is concerned with the optical collision experiments of Na atoms with a series of molecules which have been carried out in the past few years by Grosser and coworkers1–3 following similar experiments with atomic targets.1,4,5

The experiments can briefly be explained as follows.

The collisions take place in a crossed beam arrangement.

The collision complexes are excited with polarized laser light and a detection laser is used to separate the excited Na atoms from the background.

The excitation laser is detuned by a few 100 cm−1 from the resonant 3s → 3p transition of Na, as a rule, but not always to higher frequency.

Therefore optical transitions occur only at those geometries, the so-called Condon configurations, where the energy difference in the collision complex corresponds just to the detuned laser frequency.

Excitation probability as a function of detuning and polarization angle is recorded for a few selected scattering angles.

From flight time measurement the velocity of the Na atoms after collision is obtained; this allows to reconstruct the kinematics of the observed scattering processes.

The goal of the experiments is to test the interaction potentials of excited Na with the target molecules and to obtain “snapshots” of the collision complexes, in particular of the direction of the 3p orbital after the process of optical excitation.

The results are presented as polar diagrams of the probability as a function of an orientation angle; they are analysed with the aid of numerical simulations using quantum chemically calculated potential surfaces and classical trajectories.2,3,6

A related experiment consists in measuring the ratio of excited Na atoms leaving the complex in the lower 3p1/2 fine structure state compared to the sum 3p1/2 + 3p3/2, the two levels being separated by spin–orbit interaction.3,5,7

Under the conditions the experiments have been carried out, the non-zero results for this ratio can be explained only by assuming non-adiabatic coupling effects between different excited electronic states.

Therefore precise knowledge of the interaction between the potentials is needed for a detailed understanding of the underlying process.

It is the aim of the present work to explain the observed fine structure ratios as functions of the kinetic energy after the collision numerically by combining classical trajectory methods with a semiclassical surface hopping procedure.

Optical collisions

Upon collision of an excited Na(3p) atom with a target molecule the degeneracy of the atomic 3p levels is removed to an extent depending on the separation and arrangement of the collision complex.

In the case of a linear closed-shell target molecule as treated here three excited potential energy surfaces can be distinguished, two of 2A′ (electronic) symmetry and one 2A″.

The 3p charge distribution is no longer spherically symmetric as in the unperturbed atom.

In each of the excited components the 3p orbital is described by a real wave function and accordingly the excited Na atom has a quadrupole moment (of about −16.8 au).

The interaction of Na(3p) with the target is determined by the interaction of this quadrupole with the molecular multipole moments as well as by the dispersion effect.

Quantum chemically the latter is an effect of electron correlation whereas the multipole interactions can be treated already on the Hartree–Fock (HF) level.

It is important to have the correct orientation of the 3p atomic orbitals already in low order of the theoretical treatment.

In our calculations we therefore choose a slightly distorted (fixed) geometry for the target molecules such that the HF multipoles correctly reproduce the experimental moments.

Alternatively the molecular moments obtained from a correlated wave function could be used; we hope that the experimental moments also account for the average over the configurations of the rotating and vibrating target molecule.

In a previous paper6 we have explained how a qualitative understanding of the excited potential surfaces can be obtained if the most important contributions to the multipole interactions are discussed.

For every geometrical arrangement of the collision complex each one of the three 3p surfaces is characterized by the spatial orientation of the corresponding 3p orbital.

Whereas in the 2A″ state the 3p atomic orbital is perpendicular to the plane defined by the linear target and the Na atom, the two 2A′ states have in-plane 3p orbitals.

The expectation value of the 3p orbital direction after excitation leads to the so-called alignment angle, defined with respect to the center-of-mass (CM) trajectory of the incoming particle.

With quantum chemistry the 3p direction angles can easily be obtained by calculating the 3s → 3p transition moment vectors; except for very close separations they are mutually orthogonal.

For the descriptions of the experiments we also need the potential energy surface of the 3s ground state.

In this state the Na atom is approximately spherical; the surface is nearly structureless down to rather short distances where the Pauli repulsion sets in.

Furthermore on the three excited surfaces the spin–orbit interaction of the Na(3p) electron must be taken into account which contributes small corrections to the energies determined by the interaction of Na(3p) and the target molecule, in the order of the spin–orbit coupling constant of the Na(3p) state (11.5 cm−1).

Thus numerical simulation of the scattering experiments is carried out by running classical trajectories on the surfaces calculated with ab-initio quantum chemistry including electron correlation and the spin–orbit corrections.

The polarization of the laser radiation is taken into account by projecting onto the 3s → 3p transition dipole vector.

For this projection only the orientation of the transition vector is needed along the trajectories, not its value.

Therefore our method to obtain the correct multipole interaction on the HF level is very useful and saving computer time.

The three potential surfaces derived from the Na(3p) components may penetrate each other at certain geometrical arrangements, i.e. change their energetic ordering.

This situation favours non-adiabatic couplings.

Penetration of the two surfaces of 2A′ symmetry leads to conical intersections.

The 2A′ surfaces of the Na(3p) + HCCH system have such intersection points on the molecular axis, about 8 a0 off the position of the H atoms.

In the case of the N2, N2O (and CO2) target molecules, whose quadrupole moments are negative in contrast to HCCH, the penetration occurs on rings around the molecular axis with radii between 10 a0 (CO2) and 16.5 a0 (N2).

Evidently the probability that a Na trajectory passes close to an intersection point such that non-adiabatic processes can be expected is much higher for a ring-shaped intersection line than in the case of only two axial intersection points.


In this section we shall investigate the dynamics of the collision process taking place on the quantum chemically calculated potential energy surfaces.

The dynamics of nuclear motion is treated according to the concept of classical trajectories.

The relative velocity of the colliding particles at large distance and the rotational temperature of the molecule are specified by the experimental setup.

The impact parameter b and the initial rotational state of the molecule are determined by a random number generator for every trajectory.

The bond distances and angles of the target molecule are kept fixed at linear geometry, thus the surfaces become two dimensional and only five coordinates remain describing the dynamics of the system in the center of mass system.

At a distance depending on the detuning of the laser frequency the collision complex is excited to one of the states of the 3p-manifold.

Quantum chemical calculations

The potential energy and transition dipole moment surfaces of the collision complexes Na + N2 and Na + C2H2 have been presented previously.2,6

Graphical representations of their transition vectors to the uppermost surfaces and of the Condon configurations for several detuning frequencies are shown below (Figs. 1 and 2).

The value of the transition moments is practically constant for all transitions and every configuration.

The direction of these vectors characterizes the surfaces and is important for the analysis.

The sudden change of orientation in the figures indicates the position of the conical intersection where two surfaces penetrate each other.

The quadrupole moments of the molecules were computed in HF approximation leading to the following results: 5.44 au for C2H2 and −0.93 au for N2.

The quantum chemical calculations for the system Na + N2O have been performed in a similar way.

The surfaces are functions of only two geometrical parameters, the center of mass distance R between the atom and the target molecule and the angle ϑ between R and the molecular axis.

The fixed bond lengths are 2.18 a0 for NN and 2.17 a0 for NO.

This corresponds essentially to shifting the central atom by 0.06 a0 towards the oxygen atom, as compared to the experimental geometry8 and yields a dipole moment of 0.07 au (0.17 D, O being the negative pole, exp.9

0.167 D) and a quadrupole moment of −2.85 au.

For the calculations a 7s7p4d3f ANO (atomic natural orbitals) basis set from Widmark et al. has been used10 for the target molecule and a selfmade 8s6p4d3f basis for the Na atom.

The dipole vector surfaces have been calculated in HF approximation; the uppermost dipole vector surface is shown in Fig. 3.

For the interaction energies we used the MCCEPA (multiconfiguration coupled electron pair approximation) program of Fink and Staemmler11 and single determinant (HF) reference functions.

In order to justify the approximation for the orientation of the transition dipoles as introduced in section 2 we have carried out test calculations for the Na + N2 and Na + N2O systems, evaluating the transition dipole vectors with correlated wave functions.

It turned out that although the values of the moments may change by up to 10 or 20% (mainly due to the reduced weight of the reference configurations in the correlated function) the orientation angles remain constant within 1° even for geometries close to a conical intersection where the electronic states are considerably mixed.

The potential and transition dipole surfaces were interpolated by two-dimensional spline functions.

The derivatives of the potential resulted from local expansions for every grid point and subsequent interpolation by spline functions.

We found the achieved smoothness sufficient for our Monte Carlo based simulation method.

The surfaces can be obtained from the authors on request.

In order to simulate the population of the fine structure states the adiabatic potential surfaces have to be modified, so that they include the spin–orbit interaction.

This can be achieved approximately by diagonalizing the matrix for every point of the 3p-manifold.

In the asymptotic region the three adiabatic energies e1, e2 and e3 are degenerate and the correct splitting a with the spin–orbit coupling constant a = 11.5 cm−1 results for Na(3p).

Usually with positive detuning the optical excitation of the collision complex leads exclusively to the population of the highest of the 3p states.

Therefore under condition of strictly adiabatic dynamics all sodium atoms would be in the 2P3/2 state after the collision.

The involved surfaces however show conical intersections which lead to non-adiabatic transitions and to the population of the lower 2P1/2 state.

In the asymptotic limit of infinite relative velocity of the colliding particles the 2P1/2/2P3/2 ratio is given by projection of the molecular state after excitation onto the final atomlike states 2PJ.5

Thus for high collision energy the relative 2P1/2 population is expected to approach 33%.

The semiclassical trajectory surface-hopping method

We use the semiclassical trajectory surface-hopping method developed by Preston and Tully12–15 which will be reviewed shortly in the following.

For every trajectory the probability of being on the excited surface number j is given by the squared absolute value of a time dependent coefficient cj(t).

The equation for the time derivative ċj(t) can be obtained by inserting the expansion for the electronic wave function Ψel(r,R(t)) into the time-dependent Schrödinger equation.14,16

The nuclear coordinates symbolized by R(t) are assumed to be functions of time t obeying classical equations of motion; r stands for the coordinates of the electrons.

In general the resulting expression for ċj(t) contains matrix elements 〈ψj|Hel|ψk〉 of the electronic Hamiltonian as well as matrix elements of the nuclear coordinate derivatives 〈ψj|∂/∂R|ψk〉.

Conventional quantum chemical programs yield adiabatic wave functions ψa that diagonalize the electronic Hamiltonian within the space of the electronic wave functions.

The calculation of matrix elements of ∂/∂R with these wave functions (neglected in the Born–Oppenheimer approximation) is very time consuming.

On the other hand these matrix elements become large in the region of a conical intersection or avoided crossing and are responsible for the radiationless transitions between different (adiabatic) electronic states.

A unitary transformation leads to a (quasi)-diabatic basis set ψd of electronic functions,17–20 which change as smoothly as possible with the nuclear coordinates and where the matrix elements of ∂/∂R are often assumed to be negligibly small.

In the case of only two involved states ψd1 and ψd2 this transformation is determined by the definition of an angle ϕ describing the mixing of the adiabatic states and depending on the nuclear coordinates.21

In the diabatic picture the matrix of the electronic Hamiltonian is no longer diagonal, but the matrix elements H12 = 〈ψd1|Hel|ψd2〉 can easily be calculated14,21 as soon as the mixing angle is known.

The expression for the off-diagonal element reads H12 = (E1 − E2)sin ϕ cos ϕ, where E1 and E2 are the adiabatic energies; the corresponding diagonal expressions are H11 = E1cos2ϕ + E2sin2ϕ and H22 = E1sin2ϕ + E2cos2ϕ.

With aj and bj we denote the complex time-dependent coefficients of eqn. (2) in the diabatic and adiabatic basis respectively.

Using the mixing angle ϕ the coefficients of the two different pictures can immediately be transformed into each other.

Including the partial derivatives with respect to the nuclear coordinates the following system of differential equations is obtained16,21 where the matrix elements Hjk(R(τ)) and 〈∂/∂Rjk(R(τ)) are implicit functions of the runtime τ on the trajectory.

These semiclassical equations and the classical equations of motion are solved simultaneously by numerical integration.

After transformation to (quasi)diabatic electronic wave functions the matrix elements 〈∂/∂Rjk may be neglected for intermediate relative velocities .

However in our examples of optical collisions the correct 2P1/2 population for the asymptotic case of high velocity can only be obtained by maintaining the full operator.

A simple approximation for 〈∂/∂Rjk will be introduced in section 4.1.

The transition to the uppermost 3p surface takes place when a trajectory crosses the line (Condon configuration) corresponding to the appropriate energy difference defined by the detuning.

At this moment the time-dependent coefficients in the adiabatic picture are set to b1 = 1 and b2 = 0.

Every trajectory is computed twice, corresponding to the two possibilities of excitation on the incoming and on the outgoing part of the path.

We denote these two cases with type 1 and type 2 respectively.

The (complex valued) diabatic coefficients aj are computed by numerical integration using eqn. (3) and transformed to b1 = a1cos ϕ + a2sin ϕ and b2 = −a1sin ϕ + a2cos ϕ after every time step.

When the probability |bi|2 falls below 0.5 a random decision is made whether to hop to the other surface or not, and then the coefficients are reset accordingly.

In the case of a hop the momentum of nuclear motion has to be corrected in order to ensure energy conservation.

Following this procedure several hops between the two in plane states including the spin–orbit correction (1) are possible per trajectory.

As is to be expected most hops occur not too far from the conical intersection.

The number of trajectories ending in the two states 2P1/2 and 2P3/2 have to be weighted with the values of the impact parameter before being interpreted as probability.


In the neighbourhood of a conical intersection the adiabatic wave functions are rapidly changing.

This is illustrated e.g. by the behaviour of the orientation of the transition moment vector to the higher 2A′ state (shown in Figs. 1–3) which in fact maps the orientation of the Na(3p) orbital in the final state.

Due to the orthogonality constraint the 3p orbital in the lower excited 2A′ state (not shown) has an orientation rotated by 90°.

Diabatic wave functions are obtained by a transformation of the mixed adiabatic states (in our case a 2 × 2 rotation) such that a suitable molecular property behaves as smoothly as possible.

For that we compare the orientation of the 3p orbitals with the orientation of a classical 3p quadrupole Q in the field of the electric multipoles of the target molecule.

The classical quadrupole–dipole term of the interaction potential e.g. reads where μ is the dipole moment of the molecule.

The position of the sodium atom relative to the molecule is given in polar coordinates (R,ϑ) and α is the direction angle of the Na quadrupole (with α = 0 defined as the direction of the vector connecting the molecule and the sodium atom).

Analogous expressions hold for the quadrupole-quadrupole interaction and for higher interaction terms; in general the potential V(R,ϑ,α) is a sum of such terms.

For energy calculations this approach would not be sufficient because the polarization effects are missing.

But as we need only orientation angles and since polarizabilities are essentially isotropic and do not affect the 3p orientation our approximation seems to be justified.

In the classical picture we cannot distinguish two interacting electronic states.

Hence there is no conical intersection, but there are two mutually orthogonal angles α1 and α2 corresponding to minimal and maximal interaction energy.

The angles αi can be determined by the condition The direction α1 as a function of R and ϑ is also shown in Figs. 1 to 3.

In the asymptotic region far away from the conical intersection the angles αi coincide with the two directions of the transition dipole moments to the ground state, but they behave smoothly also near the intersection points.

Therefore we force the 3p orbitals of the diabatic wave functions to coincide with these classical angles everywhere.

In other words the deviation of the transition dipole vector direction to the highest 3p surface from the angle α1 can serve as a “diabatization angle” ϕ defining the transformation to the (quasi)-diabatic states.

Note that the classical orientation angles are only weakly dependent on the exact values of the molecular multipoles.

In all interaction terms (e.g. eqn. (4)) the Na quadrupole moment is only a factor which does not affect the α values of stationary energy.

Furthermore it is rather the ratio of the molecular multipole moments than their absolute value which determines the orientation angles.

The function α(R, ϑ) has been approximated using the values of the molecular dipole and quadrupole moments given in section 3.1.

In addition in order to improve the agreement of the two angles in the areas where there is no mixing of the adiabatic states we have adjusted in a semiempirical way the hexadecupole terms for all three investigated molecules.

The Figs. 1 to 3 show the directions of α1 and of the transition dipole vectors for the three molecules.

The mixing angle ϕ has been determined as difference between these two directions.

Thus in principle the transformation can directly be read from the figures.

In our examples with only two involved states a necessary condition for a considerable change of probability in (3) is that the absolute values of H12 and thus of sin ϕ and cos ϕ are not too small, and this is the case whenever the two adiabatic states are considerably mixed.

It can easily be seen in Figs. 1 to 3 that the necessary mixing angles around 45° are found only in rather restricted regions close to the conical intersections.

Results and discussion

Population of the Na(3p) fine structure states

Fig. 4 shows the resulting fine structure populations for optical collisions of Na with the molecules N2 and N2O for the rotational temperature 100 K and at a detuning value of 120 cm−1.

To facilitate comparison of our trajectory simulations with the experimental results3 (and Fig. 5), measured at a laboratory scattering angle of 18°, the abscissa is chosen to be the relative kinetic energy after collision.

Our 2P1/2 fractions represented by solid lines were computed neglecting the matrix elements of ∂/∂R in the diabatic picture.

Instead in order to ensure the correct asymptotic probability (33%) for high relative velocity a simple correction was applied: the matrix element 〈∂/∂R〉(R,ϑ) in (3) was assumed to be proportional to the radial derivative of the mixing angle ϕ.

The proportionality constant f in 〈∂/∂R12 ≈ fṘϕ(R,ϑ)/∂R was adjusted individually for N2 and N2O at high velocities.

Maintaining the same constant in the full range of relative motion leads to the dashed curves in Fig. 4.

The error bars in the figures were determined by separate statistical test runs and do not represent any systematic error of our calculation method.

A qualitative discussion of the results has to consider the position of the conical intersection as well as the Na-molecule distance at which the optical excitation occurs.

For Na + N2 the location of the conical intersection is definitely the most favorable for a high 2P1/2 fraction.

Our calculated percentage, which is slightly decreasing in the investigated range of collision energy, is in excellent agreement with measurement.3

The relative contribution of type 2 trajectories on the 2P1/2 surface decreases from about 50% to about 25% with increasing relative velocity.

For the Na + N2O system the 120 cm−1 Condon line is situated at a radial distance which lies just inside the conical intersection.

Thus less than 10% of the 2P1/2 trajectories are of type 2 except for very low collisional energy.

The experimental population of the lower spin–orbit state is reproduced satisfactorily by our computations for this collision complex too.

Only excitations to the uppermost 3p surface were taken into account in the calculations for Na + N2 and Na + N2O. In the case of Na + C2H2 corresponding simulations lead to a population of the 2P1/2 state which is considerably much lower than according to the experiment.3

The conical intersection is located in collinear geometry for C2H2, what means that much fewer trajectories pass through this region of strong mixing than for N2 and N2O, where the non-adiabatic coupling is strong at T-shaped geometry.

Further the Condon configuration for the detuning by 120 cm−1 is situated clearly outside of the conical intersection region on the C2H2 surface.

For that reason no non-adiabatic transitions to the 2P1/2 state can result from type 2 trajectories for C2H2.

The calculated 2P1/2 populations for Na + C2H2 at a detuning value of 120 cm−1 and for the rotational temperature 100 K are shown in Fig. 6.

The low probability of non-adiabatic transitions in the case of Na + C2H2 represented by the dashed curves is not able to explain the rather high measured 2P1/2 fraction.3

But in contrast to the other two investigated molecules the energy differences between the 3p surfaces are much smaller, so that excitations of the lower two 3p states are possible.

At a detuning of 120 cm−1 the number of excitations to the uppermost surface is about three times higher than to each one of the lower surfaces.

After excitation to the lower one of the two mixed A′ surfaces most trajectories will contribute to the 2P1/2 population because of the low non-adiabatic transition likelihood.

The third state (A″) with the 2P3/2 asymptotic limit was assumed not to mix with the other two states.

The results represented by solid lines in Fig. 6 were obtained including excitations to all three surfaces and are in good agreement with experiment.

A test calculation for Na + N2O showed that less than 15% of the excitations take place to the lower A′ state and many of these trajectories jump back to the uppermost surface.

So it seems reasonable to consider only excitations to the highest 3p surface except for the system Na + C2H2.

Polarization and alignment

As mentioned in section 2 the intensity of the Na(3s → 3p) transition depends on the polarization direction angle β of the excitation laser.

Table 1 shows classical trajectory results for the scattering angle ΘCM = 40° and several values of the detuning frequency.

A qualitative discussion of this method can be found in .ref. 2

β is measured relative to the incidence direction of the Na atom.

For example β = ΘCM denotes the situation where polarization is parallel to the Na momentum vector after collision.

For all three examined molecules the weight of the type 1 trajectories is dominant, leading to an optimal polarization direction in backward direction (β > 90°).

Figs. 1, 2 and 3 are useful for a qualitative understanding: For large detuning frequencies the transition dipole vector points more or less in radial direction towards the center of mass of the molecule for all possible Condon points.

The distance between the condon points for type 1 and type 2 trajectories decreases with increasing detuning, resulting in alignment angles near 90° and large alignment ratios.

For 120 cm−1 detuning frequency the direction of the transition dipole vector shows more variation for Na + N2O and Na + C2H2, explaining the rather small alignment ratios.

The calculated alignment results may be compared with measurements for Na + N2 and Na + C2H22,3 (and Table 1).

While the experimental data are reproduced satisfactorily for the N2 system the rotation of the alignment angle and the decrease of the alignment ratio at small detuning frequency are clearly underestimated for C2H2 by our calculations.

In the case of C2H2 excitations to all three surfaces of the 3p-manifold were included.

Summary and conclusion

The presented results for the population of the Na(3p) fine structure states after an optical collision with a linear molecule are in good agreement with the experiment.3

Our simulations are based on quantum chemically calculated potential energy surfaces and dipole moment functions.

Nuclear dynamics was treated in the framework of the semiclassical surface hopping method.

In order to test the effect of molecular vibration we introduced a harmonic potential for N2 allowing the exchange of vibrational energy with the other degrees of freedom.

Neither for the population of the Na(3p) fine structure states nor for the alignment direction and ratio a significant change within the statistical uncertainty could be found.

We therefore conclude that our assumption of fixed molecular geometry is a reasonable approximation.

The analysis of the results for the 2P1/2 fraction after collision shows a sensitive dependence on the details of the energy surfaces, the non-adiabatic mixing of states and especially on the location of the conical intersections.

The experimentally observed dependence of the Na(3s → 3p) transition intensity on the polarization direction angle of the excitation laser may be understood qualitatively from the properties of the potential and dipole moment surfaces.

The deviations between experiment and simulation are largest for small detuning frequencies.

The success of our simulations of the optical collision experiments is a strong indication for the correct understanding of the underlying physical processes.