Calculation of DC Stark Resonances for the Ammonia Molecule

A model potential previously developed for the ammonia molecule is treated in a single-center partial-wave approximation in analogy with a self-consistent field method developed by Moccia. The latter was used in a number of collision studies. The model potential is used to calculate DC Stark resonance parameters, i.e., resonance positions and shifts using the exterior complex scaling method for the radial coordinate. Three molecular valence orbitals are investigated for fields along the three Cartesian coordinates, i.e., along the molecular axis and in two perpendicular directions. The work extends previous work on the planar-geometry water molecule for which non-monotonic shifts were observed. We find such non-monotonic shifts for fields along the molecular axis. For perpendicular fields, we report the splitting of the 1e orbitals into a fast- and a slow-ionizing orbital.


I. INTRODUCTION
Research into ionization of ammonia molecules (NH 3 ) has been an ongoing topic of interest with recently proposed new ideas about multiple ionization (or rather the lack thereof) in the context of experimental fragmentation studies in proton-ammonia collisions at intermediate and high energies [1].The fragmentation study puts (perhaps) in doubt the validity of an independent electron model approach which was used in various forms to explain total (net) ionization cross sections, as well as doubly differential cross sections which represent the emission properties of ionized electrons, i.e., their distribution over polar angle and energy.Such studies of net ionization differential cross sections were reported in the Born approximation [2], and a continuum distorted wave method [3], which described the earlier measurements [4] quite well.The fact that these experimental data were consistent with net ionization was demonstrated in yet another study by showing the contributions from particular molecular orbitals (MOs) [5].Most of these studies employed the single-centre Slater-type orbital based Hartree-Fock calculations of Moccia [6].
From a theoretical perspective the role of multiple ionization in proton collisions with the 'isoelectronic' molecules water (H 2 O) [7][8][9], methane (CH 4 ) [10], and ammonia [1] was analyzed in the framework of the independent-atom model [11].An accurate representation of proton-water molecule differential cross sections at an intermediate energy (250 keV) was obtained with a classical trajectory Monte-Carlo method [12] which was based on a threecenter model potential.Collision calculations using this model potential have been carried out recently for the ammonia molecule [13].These works, and the problem with interpreting fragmentation cross sections [1] serve as a motivation to extend our previous studies of Stark resonance parameters for the water molecule within a model potential approach to the case of the ammonia molecule.The main idea of the model potential approach is to avoid the technical difficulties of a self-consistent effective potential.
We are studying the molecule for fixed orientation, i.e., the rotational (and vibrational) degrees of freedom are ignored.In collision problems at intermediate and high energies this approach is justified by the time scale of the collision process, and orientation averaging is applied when computing probabilities or cross sections, such as, e.g., in Refs.[11,12].For the dc Stark problem the fixed orientation with respect to the external field does represent a more serious issue, since over longer time scales the field would act to orient a molecule with permanent dipole moment [14,15].An exception would be if the molecule was found in a matrix isolation environment, e.g., by trapping in a cold rare gas matrix.This subject has received recently renewed attention in the context of proposals to measure the electron electric dipole moment using diatomic molecules [16,17].Given that strong electric fields are potentially going to be applied (cf.Ref. [18]) the problem of Stark ionization should be researched in this context, as well.
The Stark resonance problem is addressed in this paper by following the exterior complex scaling (ECS) method which implements a derivative discontinuity at the radial distance where complex scaling sets in [19].The ECS methodology was developed over the years and has been compared to the complex absorbing potential (CAP) method [20,21].Following Moiseyev [22,23] one can argue that the smooth ECS and CAP methods are equivalent.
They share the features that starting at some critical radial distance r s either a gradual continuation of the real r-axis into the complex plane is carried out, or a complex absorber is implemented for r > r s .Both methods show some dependence on either the scaling angle θ s which extends the path into the complex plane, or on the strength parameter of the CAP.
Perturbative corrections can be employed in the case of the CAP.
The hard ECS method was developed further by Scrinzi [24] as an effective absorber for time-dependent problems.A derivative discontinuity in the wave function at the scaling radius r s needs to be implemented, and it allows for the choice of scaling angles close to the critical value of θ s = π/2.This, in turn, allows to use a reduced region r s < r < r max to compute the tails of the resonance states.At the outer boundary r max the Dirichlet condition of vanishing wave function is applied.We used his methodology previously for the planar water molecule [25].The extension from a planar geometry does not pose additional problems for the present case of the ammonia molecule: a partial wave expansion of the orbitals is implemented, and as before, the radial functions are solved using a finite element method.
The choices for scaling radius r s = 16.2 a.u., and r max = 24.3a.u.were made in this work combined with a scaling angle θ s = 0.9 π/2.Atomic units (ℏ = m e = e = 4πϵ 0 = 1) are used throughout this work.
The geometry of the NH 3 molecule is shown in Fig. 1 together with the three directions along which electric fields are applied.The arrows indicate the force directions that are applied individually, i.e., one Cartesian direction at a time.The force directions F i are opposite to the electric field directions E i due to the negative charge of the electron.

II. MODEL
The model potential for the ammonia molecule is a straightforward extension of previous modelling of the water molecule (Refs.[26][27][28][29][30]).The model combines three spherically symmetric potentials for the atomic constituents.Each part contains a screening contribution, and the parameters are adjusted such that the overall potential falls of as −1/r at large distances, as is required to avoid contributions from electronic self-interaction.We keep the model parameters for the hydrogen atoms, and model the central nitrogen atom in analogy to the oxygen atom in H 2 O.
Thus, the potential is defined as follows: The scalar variables r j (with j = 1, 2, 3) represent the electron distances from the protons.
The hydrogenic parameters α H = 0.6170 and N H = 0.9075 are taken from previous works for the water molecule.The latter choice then fixes the potential parameter N N = 6.2775 to yield the appropriate asymptotic effective potential at large r, as 3 (1 is the long-range effective charge contribution to the potential from the hydrogen atoms.
The nitrogen atom screening parameter was chosen as α N = 1.525 in order for the model to yield orbital energies that follow closely values obtained in the Hartree-Fock approximation.
The geometry of NH 3 is adopted from the work of Moccia [6], with a N-H bond length of 1.928 a.u., polar angle θ p = 108.9degrees, and azimuthal angles ϕ j for the three hydrogenic protons spaced apart by 120 degrees.In particular, the azimuthal angles of the hydrogen atoms are chosen to be 90, 210, 330 degrees, which singles out the y − z plane as containing one of the protons.
The Schrödinger equation for the MOs and an electric field in the ẑ direction can be written as with r 0 ≡ r, while the Z i (r i ) are screening functions for the constituent atomic centers, i.e., Z 0 (r) = rV N (r) for the nitrogen atom, and Z i (r) = rV H (r) with i = 1, 2, 3 for the hydrogens, as defined in eq. ( 2).Note that the electric field component is E z = −F z , i.e., our notation F z refers to the force experienced by a free electron.
The MO wavefunctions ψ ≡ ψ ν are expanded in complex-valued spherical harmonics, and the radial functions are expanded using a finite-element method (FEM).The functions f in are local basis functions on interval i of the radial interval 0 < r < r max .The index n labels the polynomial basis functions [19].The Schrödinger equation is solved as outlined in Ref. [25,31] and leads to a matrix eigenvalue problem with the c inℓm being elements of the eigenvectors.With this discretization technique we are solving the three-dimensional problem, which is defined in eq. ( 3) for a field along the ẑ direction.The force direction due to the external dc field as experienced by the electron is controlled by the sign of F z as explained in Fig. 1.
The FEM approach from Ref. [25,31], and outlined in Ref. [19,32] was used with the partial wave expansion of ψ truncated at ℓ max = 3 to test how the MO eigenvalues respond to changes in the one free screening parameter contained in the model potential.The partial wave expansion allowed the spherical components of the matrix elements to be calculated using a Wigner 3j coefficient package [33] (which can be found on the author(s) homepage1 ) as before in our work with water [25,31,32].The hydrogenic potentials are expanded in spherical harmonics which allows for the use of Wigner 3j coefficients rather than evaluating three-dimensional integrals numerically.
Comparison with the SCF eigenvalues of Moccia shows that the three outermost MOs can be reproduced well with the simple model potential.The 2a 1 MO is too weakly bound at the level of 10%, which is deemed acceptable, since it is expected to contribute less to the overall molecular ionization rate.The comparison of eigenvalues obtained for ℓ max = 3 and ℓ max = 5 is provided, since the resonance parameter calculations are performed with ℓ max = 3 only.The table also contains results from a localized Hartree-Fock method as implemented in Turbomole [34][35][36][37].
Ref. [ MO eigenvalues for the model potential as compared to the SCF eigenvalues of Moccia (Ref.[6]).The E 1e energy appears twice, i.e. for the MOs 1e 1 and 1e 2 .The fourth row shows the localized HF method [34,35] eigenvalues based on the optimized geometry in HF approximations as calculated in Turbomole [37] using the def2-QZVPPD basis set, while the fifth row gives the eigenvalues from the optimized effective potential method [36] using the d-aug-cc-pVTZ-oep basis.

A. Resonance parameters for fields along the vertical ±ẑ-direction
We begin the discussion of resonance parameters with fields along ±ẑ for which the two degenerate MOs 1e 1 and 1e 2 should yield identical results.The dominant contribution to dc field ionization is expected from the weakest bound orbital (3a 1 ).For this orbital the combined molecular and external electric field leads to over-the-barrier ionization at the strongest fields calculated (|F z | of order 0.1 a.u.).
In Fig. 2 the top row shows the resonance position (left column) and resonance half-width (right column) for this orbital as a function of field strength F z .Positive values F z > 0 correspond to the field direction pushing electrons out in the direction from the hydrogen atom plane past the nitrogen atom, while negative values F z < 0 are for electrons ejected from the hydrogen plane (at negative z) away from the nitrogen atom which is located at The change of the resonance position with field strength can be described as monotonically stronger binding for F z > 0, since electron density is transferred from the hydrogen atoms in the direction of the central nitrogen atom.For the opposite field direction (F z < 0) we observe non-monotonic behavior.First one expects marginally weaker binding when transferring electron density from a nitrogen to a hydrogen atom.
In the molecule the shift in electron density will be towards regions around the partially shielded protons where the electron binding is weaker.This feature becomes apparent at strong fields (over-the-barrier regime), but there is an intermediate range of field strengths (−0.08 < F z < −0.04 a.u.)where there is a non-monotonic variation of the resonance position with field strength.
The resonance widths are obviously small in the tunneling regime.They change by orders of magnitude as the field is increased, and the ionization rate for emission from the hydrogen plane (F z < 0) is stronger than in the opposite direction, by more than a factor of two.
For field strengths of the order of 0.1 a.u.saturation in the ionization rate sets in, which is associated with the over-the-barrier regime.At these field strengths one may reach the limitations of the exponential decay model, and, thus, results for stronger fields are not reported.
-0.We note that the behavior of the resonance position is consistent with the change in ionization rate (or resonance width) as a function of field direction.In the strong field regime (at about 0.1 a.u. and beyond) the binding energies are quite different and the ionization rates change by a factor of two when the field is reversed.An interesting observation for F z > 0 is the rise in the ionization rate even though the resonance position indicates stronger binding.This phenomenon is associated with density being driven by the field towards the barrier region.
In the middle panel the results are shown for the two degenerate 1e MOs.The dependence of the resonance position on field strength F z is monotonic in this case, and varies only at the 5% level in the given field strength range.The corresponding decay rates are weaker by orders of magnitude as compared to the 3a 1 MO, and remain in the tunneling regime.This conclusion will be supported further below by probability density plots for |F z | = 0.1 a.u..
The bottom panel shows results for the more deeply bound 2a 1 MO.Here the variation in the resonance position is only at the level of 3%, and the ionization rate is suppressed by two to three additional orders of magnitude.The shape of the dc Stark shift (left panel) as a function of field strength and orientation is similar to what is observed for the 1e pair of MOs.A small asymmetry can be observed in the decay rates, with a small enhancement for In Fig. 3 we illustrate the situation with probability density contour plots of the MOs.The field-free case is shown in the middle row.The outermost MO (3a 1 ) shown on the left has an asymmetric probability density with respect to z = 0 with higher probability values on the nitrogen side.When the dc field is pushing electrons out on this side the nitrogen potential provides attraction, and causes some concentration of probability in this distribution, as shown in the top left panel (strong red drop-like shape at z > 0, and also at z < 0).The interpretation of the density plots is that they describe steady-state decay.
The bottom left panel shows the case of a strong field pushing in the direction past the hydrogen atoms.The probability distribution is more diffuse, showing that the outflow on the side of the hydrogen plane is hindered less.This observation is consistent with the decay rate results shown in the top right panel of Fig.In Fig. 4 results are presented for fields in a perpendicular direction relative to the axis connecting the N atom with the hydrogen atom plane.The arrangement of the three hydrogen atoms is such that one resides on the y-axis, i.e., field emission occurs along the direction of the H-H bond perpendicular to this axis.The degeneracy between the 1e 1 and 1e 2 MO energies is expected to be broken when a dc field is applied in this x (or the perpendicular ŷ) direction.
The top row for the outermost MO 3a 1 shows a symmetric behaviour in the dc Stark shift (left panel) and likewise a symmetric ionization rate with respect to reversal of the field direction.The change in the rise of the ionization rate at strong fields indicates that one is only approaching the over-the-barrier regime, i.e., saturation has not set in yet at . The increase in binding is at the 10% level for the strongest fields.The ionization rates for these fields are smaller than for ionization along the ẑ axis by a substantial factor (about three or six, depending on the field direction ±ẑ) The middle row shows the different behaviors for the 1e 1 and 1e 2 MOs.We classify the two orbitals as fast-vs slow-ionizing under x oriented fields (orange vs blue markers).The behavior is symmetric with respect to field orientation.The 1e 1 MO (orange dots) shifts towards less binding for both field orientations, and the 1e 2 (blue dots) is bound more deeply as the field is increased in either direction.
The ionization rates (right panel) are also symmetrical with respect to field reversal.The 1e 1 ionizes more readily by almost a factor of two for this field orientation.It is remarkable that these MOs ionize easily at strong fields with the 1e 1 MO displaying a rate which is moving towards that of the more weakly bound 3a 1 MO.Comparing the ionization rates for the 1e 1 and 1e 2 orbitals for fields along x versus ẑ we notice an order-of-magnitude increase.
The reverse trend is true for the 3a 1 MO.
The bottom row shows that the results for the 2a 1 MO are symmetric with respect to field orientation (as for 3a1).The decay rates are somewhat larger than in the case of z-oriented fields, even though the MO is bound more deeply with increasing field strength.
We support our resonance parameter values again with selective plots of probability densities from the ECS approach in Fig. 5.The middle row is identical with that in Fig. 5, but is included for direct comparison of the cases with electric field.It is immediately apparent that three MOs (3a 1 , 1e 1 , and also 1e 2 which is not shown), contribute strongly to ionization of the molecule for this field orientation.
For the 3a 1 MO we observe outflow in the form of two jets directed above and below the hydrogen atom plane.For the 1e 1 MO we find that the apparent asymmetry in the x − z plane for the field-free case has no repercussions for the outflow in the case with fields of either direction: both cases show very symmetric probabilities under these conditions which is consistent with the findings for the resonance parameters.
For the 2a 1 MO shown in the right column symmetry with respect to field orientation is expected.This is evident from the density plots by comparing the top and bottom panels.
We note the relatively strong effect the field has on this relatively deeply bound orbital.
C. Resonance parameters for fields along the ±ŷ-direction -0.In Fig. 6 results are given for field orientation along ŷ.Given the triangular nature of the hydrogen atom plane these results differ strongly from those in the previous section.By choice of azimuthal angles ϕ i = 90, 210, 330 degrees, a proton is located on the positive y axis, and asymmetry is expected when reversing the field direction.
The top row shows that this asymmetry plays a very small role for the outermost MO 3a 1 at weak fields, and is barely noticeable.The tabulated data (cf.Appendix C) show that the shifts differ by less than a percent.
The 1e 1 and 1e 2 MOs (blue and orange dots, middle row) show markedly different behavior when compared to ±x oriented fields.The shifts follow monotonic curves, as there is no longer symmetry under field orientation reversal.The 1e 2 MOs is ionizing more rapidly as compared to 1e 1 by about two orders of magnitude.The variation in resonance position is clearly at odds with this result, i.e., it apparently does not play a role here.As discussed before, the actual value of the MO binding energy is not the deciding factor, but rather how the electron density is driven towards the potential barrier by the external field.
The 1e 1 MO ionizes very weakly, and its rate is comparable to those obtained for ±ẑ oriented fields.Thus, one may conclude that the 1e 2 MO is affected by this field orientation dramatically.
In the third row we give results for the deeply bound 2a 1 MO which remains deeply in the tunneling regime for the given field strengths.It shows a small amount of asymmetry in the dc Stark shifts and in the decay rates.
The probability density plots shown in Fig. 7 again help to understand the finding for the parameter values.The results for MO 3a 1 are very close to the corresponding plots in Fig. 5 and are not shown.Instead we show over the x − y plane probability densities for the strongly ionizing 1e 2 MO (left column), the much more weakly ionizing 1e 1 MO (middle column) and 2a 1 in the right column.
For the field-free case the density plots show that the MOs 1e 2 and 1e 1 are actually not perfectly aligned with our x-and y-axes.This is related to the fact that they share the same probability density shapes, except that they are rotated by 90 degrees, and this is incompatible with the three-fold symmetry.Once the strong field is turned on along either the x-or the y-axis, however, the densities respect the symmetry of the external field, and one can identify 1e 1 = 1e x and 1e 2 = 1e y .
The results for the 1e 2 MO (left column) show an asymmetry in the outflow for the case F y > 0 vs F y < 0. It is interesting to observe that while the shape of the outflow is very different for both cases, the ionization rates are actually differ only at the level of up to 50 % In the middle column for MO 1e 1 we observe a symmetry in the outflows despite the fact that the arrangement of hydrogen atoms is asymmetric.For the case of MO 2a 1 we find that the central parts of the density are very different for the two field directions, but the parts showing the outflow are again very similar, and this is similar to what one observes for MO 1e 1 .

IV. CONCLUSIONS
We have extended our previous model calculation for H 2 O to the case of NH 3 .A simple model potential was applied to approximate what an SCF model (e.g., local density functional theory) might obtain for dc field ionization.In fact, the MO eigenvalues were shown to be comparable to the LHF results, as shown in Table I.
The results demonstrate that depending on the orientation of the electric field the three outer MOs of ammonia, i.e., 3a 1 , 1e 1 , and 1e 2 can have appreciable ionization rates.This supports the case for multiple ionization being important whether in ion-molecule collisions or in strong-field laser-molecule interactions.
Interesting details emerge from our model calculations: (i) non-monotonic dc shifts for the 3a1 MO for fields along ẑ; (ii) for perpendicular fields to the molecular axis we observe that the 1e orbitals separate into fast-and slow-ionizing ones; the (iii) fast-ionizing 1e orbital can acquire a comparable ionization rate ti the outermost 3a1 orbital.
It would be of interest to test these prediction with more sophisticated models, such as Hartree-Fock theory for MO ionization rates, density functional theory (DFT) with electron correlation, or even coupled-cluster theory of net ionization which is possible within a recently developed quantum chemistry code [38][39][40].Work is in progress to replace the model potential by exchange-correlation potentials obtained from DFT.

FIG. 1 .
FIG. 1.The geometry for the NH 3 molecule as implemented in this work showing the nitrogen atom (in green) and the three hydrogen atoms (in blue) schematically.Electric fields are applied pushing electrons out along the x-direction (black arrow), the y-direction (brown arrow), and the z-direction (red arrow), and will be denoted by positive values of F x , F y , F z respectively.Negative values of F x , F y , F z correspond to fields pushing in the opposite directions.The coordinates are given in atomic units.

FFIG. 2 .
FIG.2.Resonance positions (left panels) and half-widths (right panels) in atomic units for the outer MOs for electric fields along the axis connecting the nitrogen atom with the atomic hydrogen plane.F z > 0 values correspond to fields pushing electrons out on the nitrogen side, while F z < 0 corresponds to emission from the side of the hydrogen atom plane.Top row: 3a 1 , middle row: the doubly degenerate 1e MO (identical results for 1e 1 and 1e 2 ), bottom row for 2a 1 .

2 .FIG. 3 .
FIG. 3. Probability density contour plots for the MOs 3a 1 , 1e 1 , 2a 1 (left to right) in the y = 0 plane, i.e., as a function of x and z.Middle row: field-free case; top row: F z = 0.1 a.u.; bottom row: F z = −0.1 a.u..The contours are at values: [0.00, 0.005, 0.01, 0.02, 0.04, 0.06, 0.08, 0.1, 0.12, 0.14].The positions of the atomic nuclei are indicated by white dots with the N atom at z = 0 and the three proton locations projected onto the x − z plane.The central dot corresponds to the proton residing below the positive y-axis.

FFIG. 4 .
FIG. 4. Same as in Fig. 2 but for electric dc fields along the x−axis.In the middle panel results are shown for the MOs 1e 1 (orange dots) and 1e 2 (blue dots).The 1e 1 MO has the larger ionization rate for this field orientation.

FFIG. 6 .
FIG.6.Same as in Fig.4but for electric dc fields along the y−axis.In the middle panel results are shown for the MOs 1e 1 (blue dots) and 1e 2 (orange dots) which are clasified as the slow-vs fast-ionizing 1e orbital respectively.Note that the y − z plane contains a hydrogen atom at y > 0 and this causes an asymmetry in the resonance parameters with respect to the sign of F y .