The strong-field control of IBr photodissociation re-visited

Abstract

The photodissociation of IBr is a paradigm for a process that can be controlled by a strong, non-resonant electric field, known as the non-resonant dynamic stark effect (NRDSE). As shown by B. J. Sussman et al., Science, 2006, 314, 278, a carefully timed intense infra-red pulse can enhance or reduce the flux into the different dissociation channels. This was supported by quantum dynamics simulations using a 3-state model of IBr, but these were unable to reproduce the experimental time-scales. In this paper, we revisit this pump-control scenario using quantum dynamics simulations including all 36-states of IBr in a coupled manifold, with potentials and couplings depending on the applied field strength, i.e. including the light-molecule interaction to all orders. The results reproduce the features of the experimental control, with a better fit to the time-scale than previous simulations. The mechanism by which the control operates is then found be a combination of excited-state excitation and modulation of the avoided crossing on the dissociation pathway.

---

1. Introduction

The photo-dissociation of IBr provides a classic example of a reaction governed by non-adiabatic interactions. The excited electronic states accessed after photo-excitation are usually described by a dominant pair of states with potential energy curves that exhibit an avoided crossing due to the strong spin–orbit coupling in bromine.^1–6 This scheme can successfully describe the branching reaction which results in bromine being produced in either its ground- or first-excited electronic state.³

In a seminal experiment, Sussman at al.⁵ used a strong non-resonant laser pulse to control the photodissociation branching ratio. Using velocity map imaging (VMI), they detected iodine atoms appearing with two distinct velocities and the relative amounts of the two could be changed by changing the time delay between the resonant excitation and non-resonant control pulses. Conservation of momentum arguments were then used to convert these measurements into the branching ratio for bromine being produced in either the ground electronic state Br or first excited electronic state Br* . For short time delays between the excitation and control pulse an increase in the production of Br was seen, while for longer delays Br* was increased.

Simulations using a 3-state model (the ground state and dominant excited-state pair) supported this assignment.^4,7 These gave an explanation for the effect by a shifting of the avoided crossing by the non-resonant pulse which leads to control of the outgoing wavepacket into the dissociation channels. The simulations, however, were unable to reproduce the timescales of the experiment. In the experiments, the control pulse was active out to delays of 500 fs and the crossing point between increased Br and increased Br* was around 200 fs. In contrast, the simulations showed that the control was only active out to 250 fs with the branching ratio crossing point around 100 fs. A key point is that in the simulations the outgoing wavepacket hits the avoided crossing 80 fs after excitation, which is not easily reconciled with the experimental result if the control is due purely to manipulating the dynamics at the intersection.

In this paper, we revisit the non-resonant control of IBr with the aim of explaining this discrepancy. When all the angular momentum states of iodine and bromine are taken into consideration, the manifold of valence orbitals actually contains 36 states, coupled by strong spin–orbit coupling. This raises the question of how realistic is the simple 3-state model? The full set of potential curves was first calculated by Patchkovskii⁸ and used to extract spectroscopic data. These curves were recalculated by Sanz–Sanz and Worth⁹ as a function of applied electric field with a mind to simulating the interaction between IBr and a light field. Initial simulations with a simplified light-molecule interaction term, showed that this 36-states model do indeed show the same basic dissociation dynamics of the 3-states model despite many more states being involved.

In the following, we simulate the control experiments including the full light-molecule interaction and coupling between all 36 states. As in the experiments, it is found that the outgoing wavepacket splits into 2 components: fast and slow. The change in ratio of fast and slow components as a function of control pulse delay matches the timescales of the experiment much better than the previous calculations. The resulting mechanism for the control is found to be related not only to modulating the avoided crossing between the exit channels, but also to transitions between the manifold of states, either in the bound region (before the avoided crossing) or in the asymptotic region (after the avoided crossing).

2. Theoretical details

2.1. Electronic structure and fitting of ab initio points

The potential energy curves for the full spin–orbit coupled valence manifold of IBr were obtained as a function of field strength and angle. The description of the calculations and the fitting procedure is detailed in our previous work.⁹ Here only a brief description is included of the ab initio calculations.

The potential energy points were calculated using the Molpro¹⁰ package. The calculations were made in two steps: multi-configurational self consistent field (MCSCF) followed by a multi-reference configuration interaction (MRCI) calculation. The active space consists of 32 A′ and 13 A′′ orbitals considering C_s symmetry. The total number of states dissociating to the ground electronic states of I and Br are 18, all included in the calculations. When the spin–orbit coupling is considered the two irreducible representations double the number of states, obtaining a 36 × 36 spin–orbit matrix.

The basis set used was the diffuse augmented correlation consistent polarized Valence Quadruple-ζ basis set¹¹ for Br and the relativistic all-electron Quadruple-ζ correlation-consistent basis set^12,13 for I. The Gaussian basis sets were contracted as (22s 17p 13d 3f 2g)/[8s 7p 5d 3f 2g] for the Br atom and (28s 23p 17d 5f 3g 1h)/[12s 10p 8d 5f 3g 1h] for the I atom. This differs from the basis set used in the previous calculations⁹ and provides a better description of the system, with the two lowest spin–orbit dissociation channels at 14 325 cm⁻¹ and 17 382 cm⁻¹ above the minimum of the ground electronic state, and deeper minima for the bound excited-states (2685 cm⁻¹). It does not, however, change the ordering of the states or the response of the potential surfaces to a applied electric field.

The potential surfaces for the 36 field-free states are shown in Fig. 1. The avoided crossing responsible for the branching between the lower two dissociation channels is highlighted by picking out the ³Π_0⁺ and states involved. This pair of states are often referred to in the spectroscopic literature¹⁴ as B̃(³Π_0⁺) and Y(0⁺) and this avoided crossing is the main factor in the well-studied predissociation of IBr.^2,3,15

Fig. 1 Spin–orbit coupled adiabatic curves of the 36 states dissociating to the four spin–orbit states of I and Br. The curves in bold are the 3 electronic state model considered for comparison. The excitation pulse central frequency in the simulations is at 2.58 eV = 20 810 cm⁻¹.

The full matrix was calculated under the effect of a static electric field with intensities running from 0.000 to 0.025 a.u. and orientations between 0 and 180 degrees. The effect on the potential energy curves by the application of the electric field is shown in Fig. 2 of our previous work,⁹ as a function of field strength and orientation. The potential energy curves are strongly affected by the interaction with the field modifying the position of crossing points, depth of the minima and the relative separation within the states. This is the effect that is attributed to the control of the channel dissociation, affecting the crossings among the states via a control pulse that is applied at different time delays with respect to the excitation pump pulse.

Fig. 2 Time-slice of the momentum density after 350 fs of propagation (top panel) and momentum density (lower panel) as a function of the propagation time without control pulse for an excitation pulse with field strength s = 15 and excitation energy of ω = 2.58 eV.

The fitting of the ab initio points was a challenge for several reasons: the number of elements of the spin–orbit matrix, the sudden changes appearing in the spin–orbit elements due to the phase change of the wavefunction along the electronic structure calculations, and the swapping of states at certain configurations. These problems are exacerbated by the effect of the static electric field on the electronic properties. The procedure used to overcome all these problems was the implementation of a code that ensures continuity of the non-diagonal elements of the spin–orbit matrix using an optimisation method. The procedure is based on the idea that the transformation between spin-free and spin–orbit states is local and is not affected by the changes of the relative phase of the wavefunction. At each value of bond length, field strength and field orientation angle, the pure-spin potential matrix, W, is related to the diagonal spin–orbit potential matrix, V, by a unitary transformation

W(R, E, θ) = U^†V(R, E, θ)U (1)

with

W(R, E, θ) = W⁽⁰⁾(R) + D(R, E, θ) (2)

where D(R, E, θ) is the correction due to the applied field, which is explicitly calculated including the field in the electronic Hamiltonian at each bond-length and angle.

The diagonal elements of W(R, E, θ) and V(R, E, θ) are those calculated using Molpro.¹⁰ The off-diagonal elements of W(R, E, θ) are optimised for each electric field strength and orientation angle along the bond length distance. The optimisation method takes the calculated values of the off-diagonal matrix elements at the dissociation limit, and propagates these by using for the optimisation guess at each point the optimised off-diagonal elements of the closest bond length point. All the details of this procedure can be found in our previous work.⁹

2.2. Quantum dynamics calculations

The system dynamics calculations are obtained by solving the time-dependent Schrödinger equation using wavepacket propagation as implemented in Quantics package of programs.^16,17

As mentioned above, the potential energy surfaces of the system are calculated as a function of the IBr internuclear distance, field strength and orientation angle between the applied electric field and the internuclear molecular axis. In our preliminary results,⁹ the propagation was run in the full 36 electronic states, The potentials, however, were not dependent on the strength and/or orientation and the light-molecule interaction was modelled using dipole-field and polarizability-field interaction terms: the calculated field-free dipole moments where used along with a constant polarisability, obtained from the field-dependent potentials as the derivative of the dipole moment with respect to the field strength in the vicinity of the avoided crossing. The results did not provide the control time-scales seen in the experiments.

In the new propagations presented here, the potential energy surfaces directly include the dependence on the field strength. The orientation angle, however, was not included, because one complete propagation including field strength and orientation took several months to run, and considering that several propagation must be done to observe the effect of the control pulse at different time delays, the computational effort was impracticable. The angle was thus fixed so that the electric field lies along the z-axis in the energetically favourable direction with respect to the IBr dipole moment.

The Hamiltonian used for the time-dependent propagation can thus be written as follows

Ĥ_jk(R, t) = T̂δ_jk + W⁽⁰⁾_jk(R) + V_jk(R, E(t)) (3)

where the pair of subscripts denote that the quantities are matrix elements. T̂ is the kinetic energy operator, which is taken to be diagonal in the pure-spin basis, W⁽⁰⁾_jk(R) is the field-free potential matrix with the pure-spin adiabatic potential energies in the diagonal elements and the field-free spin–orbit couplings in the off-diagonal elements. The time-dependent potential function, V_jk(R, E(t)), is the field-driven correction to the field-free potentials, the matrix D in eqn (2) with the time-dependent field given as a sum of two pulses,

E(t) = E₁(t) + E₂(t) (4)

a pump pulse (E₁(t)) and a control pulse (E₂(t)). Both pulses are described by a time-dependent envelope (t), a central frequency ω, and a central time, t₀

E_i(t) = s_i[scr E, script letter E]_i(t)cos(ω_i(t − t_i0)) (5)

The parameters for the pulses were taken from Stolow and Sussman.⁵ In the experiments, the frequency of the pump pulse lies in the visible (520 nm = 2.38 eV) and the control pulse in the IR (1700 nm = 0.73 eV). For the simulations, however, the excitation pump pulse was given a frequency of 2.58 eV as this is where the calculated ³Π_0⁺ state lies at the Franck–Condon point. The parameter s_i controls the strength of the pulse. For each pulse, the envelope was a normalised Gaussian function with width τ_i defined by the full width half maximum (FWHM) of the pulse, i.e.The width of the pump pulse was τ₁ = 100 fs and the width of the control pulse τ₂ = 150 fs. Strengths of s₁ = 15 and s₂ = 75 were taken for the excitation and control pulse respectively, which correspond to maximum field strengths of 0.00340 a.u. (1.20 × 10¹⁴ W cm⁻²) and 0.0114 a.u. (4.00 × 10¹⁴ W cm⁻²). These are stronger than the experimental pulses, but this strength was needed as the excitation is very weak, and lower strengths result in the signal being lost in the numerical noise.

The light-molecule coupling, V_jk(R, E(t)), is thus defined differently from our previous work.⁹ The pulses interact directly with the potentials and the non-resonant interaction within the electronic states comes from the value of the off-diagonal field-dependent couplings.

As the potentials are a function of bond-length and field-strength, the propagations are 2-dimensional in the 36 spin–orbit coupled electronic states. The potentials and couplings for the full manifold were provided by spline fits to the calculated points, and the applied light-field treated as a time-dependent function that moves along the coordinate for the electric-field. Effectively, the one-dimensional potentials due to the IBr bond are smoothly changed as the field changes. This is an extension to coupled states of the electric-nuclear Born–Oppenheimer approach of Balint–Kurti and co-workers^18,19 and includes the light-molecule interaction to all orders.

The grids used for all the 2-dimensional propagations are: for the internuclear distance a fast Fourier transform (FFT) collocation, with a total of 512 points running between 3.80 and 18.00 Bohr and for the field strength a sine-DVR with 101 points running between −0.025 and 0.025 a.u. The propagation used a Runge–Kutta integrator of order 5. To avoid the use of even long grids, the wavepacket is absorbed before it reaches the end of the grid using a complex absorbing potential (CAP). The CAP has the form −iW(Q) = −iη(Q − Q_c)^bθ(Q − Q_c), where θ(x) denotes the heaviside's step function and Q_c is the grid point where the CAP is switched on. The optimised CAP used in this study has an exponent set to b = 3, a strength η = 0.05793 Hartree and is placed at R = 17.00 Bohr.

The initial wavepacket is formed in the same way as in our previous work using energy relaxation, propagating an initial guess in imaginary time until convergence is reached. Propagations including the pulses were run from −500 fs to +600 fs. The large number of coupled states, long grid for the dissociative coordinate and long propagation time means that these calculations are quite expensive, requiring around 31 hours of CPU time on a Dell PowerEdge R450 Server with Intel Xeon Silver 4316 CPUs.

3. Results

3.1. Propagation with no control pulse

In our previous work⁹ the natural dynamics of the 36-state was investigated and it was found that, in line with experiments, IBr excited around 500 nm dissociates into the two lower channels. In the following, channel 1 refers to the lowest dissociation channel leading to I + Br, and channel 2 to the first excited dissociation channel, which leads to I + Br*. It was also found that the strong avoided crossing between the ³Π_0⁺ and states, where crossing from states in the lower to upper channel takes place, is reached approximately 80 fs after excitation after which time the wavepacket enters the dissociation channel.

In the NRDSE experiments of Sussman et al.⁵ the observable measured is the velocity of the iodine atoms after dissociation which relates directly to the momentum of the outgoing wavepacket. In Fig. 2 the density of the IBr wavepacket in the momentum representation is shown as a function of time after excitation. At t = 0 the wavepacket is formed in the excited-state manifold and immediately picks up a momentum of 3.4 a.u. In the first 100 fs this splits into two streams with different momenta, centered around 3.0 and 4.2 a.u. Using the usual conservation of energy arguments the fast channel is assigned to dissociation in the lower channel 1 (I + Br), and the slow to dissociation in the upper channel 2 (I + Br*). The density is attenuated due to absorption of the wavepacket by the CAP. This happens for the slow channel by 500 fs, and the fast channel 400 fs.

A node is seen in the density for the slower stream as it splits from the density of the initial wavepacket and moves into the dissociation channel 2. This is due to having to cross the avoided crossing to reach channel 2. The long-time momentum distribution, taken at 350 fs, is shown in the top panel of Fig. 2. As in experimental distribution, the slow channel is larger than the fast channel. The difference is, however, not as pronounced as in the experiments, which indicates that the coupling in the model is somewhat too strong, overestimating the adiabatic dissociation via the channel 1.

3.1.1. Propagation with a control pulse. The simulation of the IBr photodissociation was repeated with the same 100 fs excitation pulse with energy 2.58 eV and strength s = 15 and a control pulse with energy 0.73 eV of 150 fs width and strength s = 75 at various time delays, from Δt = −80 to 200 fs in 10 fs intervals, from 200 to 280 fs in 20 fs intervals and then at 320 and 400 fs to cover the long time delays. The momentum density from some representative simulations is shown in Fig. 3. A delay of Δt = 30 fs has the control pulse arriving before the wavepacket gets to the avoided crossing. Δt = 80 fs is when it is at the avoided crossing. Δt = 120, 180 fs are when the wavepacket is moving away from the crossing, and Δt = 240, 320 fs is when the wavepacket is in the dissociation channels and away from the interaction region.

Fig. 3 Momentum density as a function of the propagation time for different control pulse time delays.

The Fig. 3 shows a different behavior for short, medium, and long time delays of the control pulse. At short time delays, Δt < 120 fs, the control pulse further excites the wavepacket to other states in the manifold, leading to higher momentum. Compared to the no-control density in Fig. 2, the slow momentum channel at 3 a.u. is clearly depleted more than the high momentum channel at 4.2 a.u. With medium control delays between Δt = 80 fs and Δt = 180 fs it can be seen that the channel from the initial wavepacket to the slow momentum channel opens and the node in the stream disappears. Finally, at long delay times, Δt > 240 fs, the secondary excitations do not appear and the momentum density returns to the initial two channel picture, but the slow channel is enhanced, presumably due to excitation by the pulse between states within the dissociation channels.

The short and long-time effect of the control pulse is thus due to excitations between states, either in the interaction or dissociation regions of the potentials. The medium-time behaviour, where the slow channel opens up, is the result of the non-resonant pulse tuning the gap of the avoided crossing. This is actually a counter-intuitive behaviour. As shown in the previous paper,⁹ the gap widens as the electric field increases and this should lead to less population crossing to the slow upper channel and the electric field must provide the coupling between states to overcome this.

To compare the simulations directly with the experimental results the momentum density at long times as a function of control pulse delay is required. The long-time was taken at 350 fs when the outgoing wavepacket is in the dissociation channels but not yet reached the CAP. At long times the momenta will not change and this will relate directly to the measured velocities. It should be noted that at 350 fs the effect of the control pulse at delay times greater than 300 fs will stay play a role, but the cost of the simulations prohibits longer times being taken.

The plot of the long-time momentum density. Shown in Fig. 4, has the same features of the experimental plot of velocities as a function of control pulse time delay (Fig. 2 in ref. 5). When the control pulse strongly overlaps the excitation pulse, around t = 0 fs, the density in both channels is severely depleted. Once the centre of the control pulse is after the excitation pulse, t = 100 fs, the density in the fast channel is seen to increase more than the slow channel. At later times, t = 150 fs, the slow channel density is enhanced. Thus the simulations capture the effect of the non-resonant control pulse.

Fig. 4 Long-time (350 fs) momentum density as a function of the control time delay out to 300 fs.

This behaviour can be related to the momentum densities plotted at specific time delays in Fig. 3. The depletion of signal at short time delays is due to the further excitation of the wavepacket due to the overlap for the control pulse with the pump pulse. That the fast channel density rises faster than the slow is then due to the fact this excitation depletes the slow channel more. The opening up of the slow channel then leads to the rise of the slow channel density. At medium time delays this is due to modulating the avoided crossing, while at longer delays this is due to excitation between the states in the dissociation channel.

Finally, to summarize the effect of the control pulse we integrate the momentum density under the two peaks at long time. This provides the amount of system, P_i, going into the different channels, here the fast and slow components. The ratio of these amounts is the branching ratio.

The change in this branching ratio as a function of control time, defined as
where P₀ is the branching ratio with no control pulse. This is shown in Fig. 5.

Fig. 5 Integral of slow and fast components of the momentum grid population as function of the control time delay. Four types of excitation pulses are included, two field strengths s = 10, 15 and two excitation energies (frequencies) ω = 2.38, 2.58 eV (521, 481 nm).

This plot again has a similar form to the experimental curve in Fig. 4 of ref. 5 and follows the changes in momentum density with control pulse delay times described above. At short times the branching ratio changes in a negative direction, indicating an enhancement of the fast channel. This, like the experimental plot, peaks around 80 fs before the slow channel opens and starts to redress the balance. After 150 fs, the positive deviation shows that the slow channel is being enhanced, as in the experiments. That this does not rise as much, or for as long as in the experiments, may be due to the fact the simulations become unreliable at longer times due to the timescale of the propagations.

Unlike the experiments, there is also a positive peak starting at negative times. This indicates that the overlap of the control pulse starting before the pump pulse is having a significant influence on the excitation, meaning the pump pulse is too strong to be perturbative. A relatively strong pump pulse, however, was required due to the weak nature of the excitation.

4. Conclusions

A Hamiltonian has been set up for the study of IBr interacting with an electric field, including the light-molecule interaction directly in the potential functions and couplings. All 36 valence states are included. This Hamiltonian has then been used for the simulation of IBr photodissociation controlled after photo-excitation by a timed non-resonant IR pulse.

After photo-excitation above the second dissociation limit, the observed behaviour mirrors that seen experimentally. Despite the complexity of the many coupled states, with transitions occurring between all of them, the dissociation with just a pump pulse is found to take place with two distinct streams of different momenta. Simulations including the control pulse show an enhancement of the fast momentum channel at short time delays, followed by an enhancement of the slow channel at longer time delays.

Examination of the simulations show that the effect of the control pulse can be divided into 3 phases. In the short time, the control pulse depletes the channels, affecting the slow channel more resulting in an enhancement of the fast channel. When the delay time is similar to the time it takes for the out-bound wavepacket to reach the avoided crossing that controls the branching ratio, the control pulse modulates the avoided crossing, opening the slow channel to change the enhancement. At longer time delays, after the crossing, the slow channel is enhanced by transitions in the exit channel from the fast stream.

In summary, we have set up a complete ab initio model for the interaction of IBr with a light field and found it able to describe and explain the strong-field control of this molecule.

Author contributions

CSS ran all the calculations and the analysis. GAW provided the code required for the propagations. Both authors discussed the results and wrote the paper.

Conflicts of interest

There are no conflicts to declare.

Data availability

The Quantics code used for the calculations is open source and available on request to the authors (https://www.chem.ucl.ac.uk/quantics). The input and potential files for the calculations are available from the UCL data repository (https://rdr.ucl.ac.uk) at the DOI: https://doi.org/10.5522/04/29205701.

Acknowledgements

This work was initially funded by a grant from the EPSRC EP/G014124/1. CSS would like to acknowledge the Ministry of Science and Innovation (MICINN, Spain) for the funded grant PID2021-122549NB-C22. GAW and CSS both acknowledge the EPSRC for funding under the COSMOS programme grant EP/X026973/1.

Notes and references

1. W. G. Brown, Phys. Rev., 1932, 42, 355–363.
2. M. S. de Vries, N. J. A. van Veen, M. Hutchinson and A. E. de Vries, Chem. Phys., 1980, 51, 159–168.
3. H. Guo, J. Chem. Phys., 1993, 99, 1685–1692.
4. M. Y. I. Benjamin, J. Sussman and A. Stolow, Phys. Rev. A, 2005, 71, 051401.
5. B. J. Sussman, D. Townsend, M. Y. Ivanov and A. Stolow, Science, 2006, 314, 278–281.
6. Y. Kobayashi, K. F. Chang, T. Zeng, D. M. Neumark and S. R. Leone, Science, 2019, 365, 79–83.
7. C. Sanz-Sanz, G. W. Richings and G. A. Worth, Faraday Discuss., 2011, 153, 275–291.
8. S. Patchkovskii, Phys. Chem. Chem. Phys., 2006, 8, 926–940.
9. C. Sanz-Sanz and G. A. Worth, Phys. Chem. Chem. Phys., 2019, 21, 14429–14439.
10. H.-J. Werner and P. J. Knowles and G. Knizia and F. R. Manby and M. Schütz and others, MOLPRO, version 2009.1, a package of ab initio programs, 2009, 2009, https://www.molpro.net.
11. A. Wilson, D. Woon, K. Peterson and J. T. H. Dunning, J. Chem. Phys., 1999, 110, 7667–7676.
12. G. Schoendorff and J. A. Boatz, J. Phys. Chem. A, 2022, 126, 4848–4861.
13. D. Bross and K. Peterson, Theor. Chem. Acc., 2014, 133, 1434.
14. M. S. Child and R. B. Bernstein, J. Chem. Phys., 1973, 59, 5916–5925.
15. A. B. Petersen and I. W. M. Smith, Chem. Phys., 1978, 30, 407–413.
16. G. A. Worth, Comp. Phys. Commun., 2020, 248, 107040.
17. G. A. Worth, K. Giri, G. W. Richings, M. H. Beck, A. Jäckle and H.-D. Meyer, Quantics package, Version 1.1, 2015.
18. G. G. Balint-Kurti, F. R. Manby, Q. Ren, M. Artamonov, T.-S. Ho and H. Rabitz, J. Chem. Phys., 2005, 122, 84110.
19. S. Zou, C. Sanz and G. G. Balint-Kurti, J. Chem. Phys., 2008, 129, 124307.

---