Computational Prediction of the Structure of Coulomb Exploded Molecules

Julia Bellamy, Hobart and William Smith Colleges/Dartmouth College, Physics & Engineering Major
Mentored by Dr. Loren Greenman

Coulomb Explosion Imaging (CEI) is a method in which femtosecond X-ray free-electron (XFEL) and strong IR lasers can charge up molecules through sequential multiphoton ionization that then explode in charged fragments [1]. The final momenta of these charged fragments are recorded, and contain information about their initial geometry that are obscured in other imaging methods. There are many methods to determine the geometry of molecules such as studying reaction pathways, analysis through molecular orbital theory, performing spectroscopy, solving the molecular Schrodinger equation, and so on. These all have their advantages, but a common limitation is that no current method can directly determine the geometry from CEI experimental data.

Fig 1

Fig. 1: The second order differential equation relating the mass and the second derivative of the position to the sum of the Coulomb force that can be used to predict the initial geometries of CEI molecules.

The goal of this research project is to predict the geometry utilizing Newtonian dynamics. Relating Newton’s Second Law to the Coulomb force, the second order differential equation (Figure 1) can be solved using standard integration methods. To solve a second order differential equation there is a need for two initial conditions: the position and the momenta. This poses a problem since the initial position is what we are trying to predict, so the development of a code utilizing an inversion method is necessary.

The code developed for this research has three main classes. A Molecule class that holds information about the charges and masses of each charged fragment, a Testing class that simulates CEI final momenta data (that eventually will be replaced with a Data class to utilize experimental data from the JRM Lab), and an Inverter class that solves for the initial starting geometry.

Fig 2

Fig. 2: A depiction of the the inversion method, where the final momenta from a three body simulated CEI is partitioned into its initial potential and kinetic energy. The inversion of the potential energy and the kinetic energy into the initial position and velocity provide the initial conditions to solve the second order differential equation.

 

The Inverter class is aptly named for running the Inverter method that is the baseline of this code. Given a three body simulated CEI explosion, the final momenta can be converted into its final energy, which is all kinetic. A partitioning scheme is then applied to the final kinetic energy, where the final energy is divided into the initial energy conditions. The total initial energy for a three body charge system has three potential energy terms and six kinetic energy terms. If the final energy were partitioned 15 ways, then 15 would be divided in 9 ways (or 3 ways if neglecting initial kinetic energy). The partitioning method takes into account that the potential energy between fragments with different charges and masses are unique, so there is also a permutation of all possible partitions generated. After assigning all the final energy into its fractional initial energy conditions, the potential energy can be rearranged (a.k.a. inverted) to solve for the initial positions and the kinetic energy can be inverted for the initial velocities. These then given the initial conditions to the second order differential equation and a method to determine the initial geometry. This inversion method is depicted graphically in Figure 2.

Fig 3

Fig. 3: Simulation of OCS CEI for the O+C+S+ charge channel where the color indicates the time of the explosion in attoseconds. Left: x and y components of the position of the three fragments over time, where the purple indicates the final position. Right: x and y components of the momenta of the three fragments over time, where the purple indicates the final momenta.

The inversion method was tested by simulating OCS CEI explosion of the O+C+S+ charge channel. The mass of O, C, and S were used and a randomly assigned initial geometry was generated. The final momenta after the simulated propagation of CEI was used in the inverter method described above (see Figure 3). Partitions of 15, 30, and 60 were run under the assumption that there was no initial kinetic energy. It can be seen in Figure 4 that as the partitions increased the number of solutions also increased (91, 406, and 1711 respectively). There seemed to have been close initial guesses to get final momenta similar to the simulated final momenta for one charged fragment but not the other two.

Fig 4

Fig. 4: The x and y components of the final momenta of interest. The open blue circles indicates the x and y final momenta of the simulated OCS explosion. The smaller, colorful dots are all x and y components final momenta solutions given 15, 30, and 60 partitions respectively.

From each solution point, it has its associated initial geometry that was compared to the simulated starting geometry of OCS. It can be seen that the initial geometries for 60 partitions in Figure 5, that the lowest non-normalized RMS error is 1.23 (red line) and differs in position and shape to the simulated geometry. Comparing to 15 partitions (not pictured), the error only decreased from 1.24 to 1.23. This is an order of 10 increase in the preciseness of the energy partitions but only a 0.01 change in the error.

Fig 5

Fig. 5: The x and y positions of the initial position of the simulated OCS geometry (blue line) and the inverter solutions of the initial geometry, depicting the lowest non-normalized RMS error (dotted lines) rotated to the same COM frame as the simulated geometry. The smallest error is 1.23 given by the dotted red line.

To determine if this method has validity in trying to predict the geometry of CEI molecules, 1000 simulations of OCS were run under the conditions described above, all with different randomly assigned initial geometries. A comparison of the non-normalized RMS geometry error to the momenta error indicated a high density of points at low momentum error and similar geometry error (see Figure 6). This suggests that there is a consistency to the solutions and incorporating further physics (initial kinetic energy, charge-up model, etc.), computational power (more precise energy partitions, N-fragment framework, 3D space, etc.), and comparison to experimental data would bring the inversion method up to predicting the initial geometries of CEI molecules.

Figure 6

Fig. 6: Comparison of the error between the simulated initial position and the solution initial position to the error between the simulated final momenta and the solution final momenta for 1000 simulations of O+C+S+. All errors are non-normalized RMS errors.

References

[1] X. Li et. al., Phys. Rev. Research 4 (2022)

 

Acknowledgments

Firstly, I would like to thank Dr. Loren Greenman for being a great advisor for this research project and helping to develop the code. I would also like to thank Dr. Daniel Rolles for sharing his experimental CEI data and providing a link between this project and Lana Chaleunrath-Pham's REU project. I am also grateful to Jeremy Kamman for his camaraderie during this research. I would also like to extended my appreciation to the rest of the Physics REU 2022 cohort for the kindness and shared learning experience that we got to experience this summer. Additionally, thank you to Bret Flanders, Loren Greenman, and Kim Coy for organizing and managing this program, as well as Kansas State University and the National Science Foundation for funding of the program.

Final Presentation