Using Machine Learning to Understand the Causes of Quantum Decoherence in Solution-Phase Bond-Breaking Reactions

Decoherence is a fundamental phenomenon that occurs when an entangled quantum state interacts with its environment, leading to collapse of the wave function. The inevitability of decoherence provides one of the most intrinsic limits of quantum computing. However, there has been little study of the precise chemical motions from the environment that cause decoherence. Here, we use quantum molecular dynamics simulations to explore the photodissociation of Na2+ in liquid Ar, in which solvent fluctuations induce decoherence and thus determine the products of chemical bond breaking. We use machine learning to characterize the solute–solvent environment as a high-dimensional feature space that allows us to predict when and onto which photofragment the bonding electron will localize. We find that reaching a requisite photofragment separation and experiencing out-of-phase solvent collisions underlie decoherence during chemical bond breaking. Our work highlights the utility of machine learning for interpreting complex solution-phase chemical processes as well as identifies the molecular underpinnings of decoherence.

an appropriate solvent density at the simulation temperatures (1.26 g/mL at 120 ± 2 K) so that the system was well in the center of the liquid region of the L-J phase diagram.
In all simulations, periodic boundary conditions were implemented with minimum image convention 1 and all interactions were tapered smoothly to zero at 16 Å over a 2 Å range with a center of mass-based switching function according to Steinhauser. 2 All simulations were performed in the microcanonical ensemble.
Because all interactions were taken to be pair-wise additive, the full Hamiltonian of the system is Ĥ = H cl + Ĥqm .In atomic units, the classical portion of the Hamiltonian is given by where v χ i is the velocity of the i th sodium (χ = Na + ) or solvent atom (χ = solv) at position R χ i and mass m χ .U χ−γ is a classical potential between atom types χ and γ (χ, γ = Na + or solv).
The quantum Hamiltonian is given by where pi and r i are the momentum and position operators for electron i, respectively, and V χ is the pseudopotential representing the interaction between an electron and atom type χ.
The classical interaction between the two Na + cations was modeled through a pointcharge Coulomb potential, U Na + −Na + (R) = 1/R, since the short-range repulsion between the cores is negligible around the internuclear separation of Na + 2 .
Classical interactions were modeled with Lennard-Jones potentials: where r ij is the distance between the i th and j th solvent/Na + site, q i is the charge on the i th site, ϵ ij is the potential well depth, and σ ij is the finite distance at which the inter-particle potential is zero.The Lennard-Jones parameters 3 used in this study are listed in Table S1.
Phillips-Kleinman (PK) pseudopotentials were used to account for the interactions be-tween the classical particles and the quantum mechanical electron. 46][7] For the e − -Ar interaction, V e−Ar , we used a modified version of the pseudopotential developed by Gervais et al., 8 described in detail in our previous work. 9 For the e − -Na + interaction, we used rigorously-derived psedopotentials previously developed by our group, the details of which can be found in Refs.6 and 7, respectively.The final pseudopotential fits are presented here in Tables S2 and S3.

Simulation Setup
The eigenstates of the quantum mechanical valence bonding electron were expanded on a three dimensional grid.We used a grid that contained 32 × 32 × 32 grid points.These dimensions were chosen to keep the basis set as small as possible for each system while still capturing the spatial extent of the electronic wave function.We centered the grid in the middle of the simulation cell and shifted all classical particles relative to the grid every 500 fs to avoid leakage of the wave function off the edges of the grid.In this way, the wave function was always located roughly in the center of the simulation cell.The classical particles were shifted an integer number of grid spaces to avoid discontinuities in the quantum energy that would prevent total energy of the simulation from being conserved. 10We used the velocity Verlet algorithm 1 to propagate the classical degrees of freedom (v Na + i , v solv i , R Na + i , and R solv i ) of the Hamiltonian in Eqs. 1 and 2 in the microcanonical (N, V, E) ensemble.We determined the forces from the sum of the classical-classical and classical-quantum interactions described above.We used the implicit restart Lanczos method to iteratively solve the TISE for the ground state wavefunction at every 4 fs time step. 11The quantum forces on the classical particles were then found using the Hellman-Feynman theorem: where, F Q i is the quantum force on classical particle i at position R i .Because the wave function is expanded in a basis that does not functionally depend on the position of the classical particles, Eq. 4 is formally exact (in other words, there are no issues with Pulay forces from the basis functions changing with time). 12r this paper, we collected data from 210 non-equilibrium dissociation trajectories of Na + 2 in liquid argon.Our initial configurations were generated from uncorrelated groundstate configurations and placed onto the electronic excited state at time zero of our nonequilibrium trajectories.The dynamics were propagated with nonadiabatic surface hopping, using Tully's fewest switches surface hopping (FSSH) algorithm. 13For our analysis in the main text, all data is collected prior to any instances of surface hopping onto the electronic ground state and the decoherence of interest occurred adiabatically on the excited state before any surface hops took place in any of the trajectories.

Difference in Local Solvent Potential
The difference in local solvent potential measures the difference in integrated potential felt by the electron within a 2.6 radius around each Na + .This is calculated using the value of the pseudopotential on our grid basis.We take all the grid points within a 2.6 radius of each Na + and sum the value of the potential at these points.The orange curve plotted in Figure 1 shows the absolute difference between the integrated pseudopotentials around each Na + in the 180-fs window prior to localization.Figure 1 shows that around the ∼ 60 fs prior to localization, the value of the absolute difference in the integrated pseudopotential between each Na + sharply increases.The trend suggests that the local environments around each Na + become significantly different as the system approaches localization.The units of |∆V e − −Ar | are expressed in the reduced units of the pseudopotential, where A = 0.28 Hartrees.

Local Solvent Density
The local density around each Na + is calculated using a density field method from the work of Willard and Chandler, 14 with a coarse-graining length of 1.7 .The grid used in the calculation is the electron grid basis.The curves shown in Figure 2

Collective Solvent Velocity
To visualize the solute-solvent collisions coinciding with decoherence we calculated the components of the collective solvent velocities moving away or towards each Na + .We weight the distance of each solvent atom to each grid point using the same weighting function as that in the density field.To account for the motions of the Na + we subtract off its velocity vector from all argon velocity vectors to view the velocity field in a locally stationary frame.
We then sum the distance weighted instantaneous relative velocity vectors to measure the collective argon velocities around each grid point.For instance, two argons at the same distance from a grid point with exactly opposite instantaneous velocities will return a collective velocity of zero at that grid point.
In order to measure net direction of collective solvent motion (ϕ), the grid is generated as a spherical shell around each Na + with a radius of 2. In Figure 3 panel a, the net collective solvent velocity is negative in the first 100 fs after excitation corresponding to the inward motion of solvent velocity as the fragments dissociate.
After the caging event which occurs at ∼100 fs after excitation for this system, we observe a positive value for ϕ, indicating the net motion of solvent away from each Na + .In this time regime, the two Na + atoms do not experience significant differences in their local solvent motions and suggest that the caging collisions with the two fragments occur simultaneously, on average.When viewed as the time to localization, ϕ shows significant differences between the localized and unlocalized Na + .Most notably, ϕ is negative for the 60 fs leading up to localization for the Na + that receives the electron, showing the net inward motion of solvent in this time regime.As expected, the opposing Na + has a positive value for ϕ, suggesting a lack of collisions is needed for localization.

Machine Learning Details
Balanced Random Forest Classification

Classifier Implementation
Our Balanced Random Forest (BRF) classifier was implemented in Python version 3.9.4 using imblearn version 0.7.0 along with scikit-learn version 0.24.2.Table 3 shows the features used to train the model, along with a feature description.Hyper-parameter optimization was done using a grid-search method to tune n_estimators and the max_depth of the random forest model.The chosen hyper-parameters are shown in Table 4.
Features that were not used in the final model but were included in the initial feature set included atom-centered-symmetry functions (ACSFs), spherical pseudopotential values, and effective volume values for each sodium.Spherical pseudopotential values were identical to the integrated pseudopotential, except that the values were integrated on the surface of a sphere around each sodium rather than the entire volume of the sphere.ACSFs were implemented following the formalism of Behler 15 and hyperparameters were tuned using a random search method.

Feature Set & Training/Testing
Since every 9-tuple of our data set was taken from a single MD trajectory, we split our data into train and test splits by allocating certain trajectories for training and testing.
This eliminates the possibility of data leakage where correlated examples from the same trajectory end up in both the training and testing data, potentially leading to inflated model performance.We also trained models using a random test/train split methodology and nearly identical results were obtained.
The learning curve for the BRF classifier is shown in Figure 5. Table 5 shows the balanced accuracy scores for the various ways we evaluated model performance.We also tested model

Principal Component Analysis (PCA)
We additionally did exploratory principal component analysis on the classifier data set.This is shown in Figures 11-14.

Regression Implementation
Our Gaussian Process (GP) regression model was implemented in Python version 3.9.4 using scikit-learn version 0.24.2.A radial basis function (squared-exponential) kernel was used, with an initial length scale setting of 1.0, and a minimum and maximum length scale of 1E-2 and 1E2 respectively.

Feature Set & Training/Testing
Our feature space for regression included several physically motivated encodings of the solute local environment.Atom-centered symmetry functions (ACSFs) corresponding to equations 5 and 6 with the cutoff function given by equation 7 were computed using an in-house Mathematica code implemented in Mathematica 13.1.For all ACSFs, R ij and R ik denote sodium-argon distances, R ik denotes argon-argon distances, and θ ijk corresponds to argonsodium-argon angles.The hyper-parameters η, S , λ, and ζ were optimized using a random search method.
In addition to the ACSFs for each sodium atom, we included dimer bond distance, spherical sodium-argon pseudopotential, and effective volumes around each sodium, to our feature

8 .
The relative solvent velocity field is calculated on this grid and we take the component of the collective velocity along a normal vector to the grid surface.A visualization of the spherical grid with the normal components of the collective velocity vectors are shown in Figure3panel c(localized) and panel d(unlocalized) for a single timestep during localization.The net direction of the collective solvent velocities is then the sum of all normal vectors on the shown spherical surface.The plots in Figure 3 panel a and panel b are the ensemble averaged net collective solvent velocities in the specified time regimes.

Figures 8 -
Figures 6 and 7 shows the standardized distributions of each feature for the training and test sets, respectively.

Figure 6 :
Figure 6: Distributions of the standardized features of the training set for the balanced random forest model.

Figure 7 :
Figure 7: Distributions of the standardized features of the test set for the balanced random forest model.

Table 1 :
Lennard-Jones and Coulomb Potential Parameters for the Solute-Solvent Systems Studied in This Work.

Table 2 :
Parameters used in the construction of the sodium-electron pseudopotential, ϕ(r) = 3 i=1 c i e −α i r 2 .Further details of the sodium-electron pseupotential can be found in Ref 6.

Table 4 :
Chosen hyper-parameters for the Balanced Random Forest classifier model, evaluated using a grid-search method.

Table 5 :
Balanced accuracy scores over the various model validation methodologies for the BRF classifier.A single train/test split, cross-validation, and replicate train/test splits produced very similar balanced accuracy scores for the BRF model.This indicates that the reported balanced accuracy score is robust to sampling biases.Note that the error shown in standard error of the mean.

Table 6 :
Machine learning features for the GP regression model along with feature descriptions Na-Ar pseudopotential integrated over the surface of a 2.8 angstrom sphere centered on the Na + (2) core effective volume Na + (1) Free volume around the Na + (1) core effective volume Na + (2) Free volume around the Na + (2) core G radial,N a + (1) G radial evaluated from the Na + (1) core on the nearest 50 argon atoms G radial,N a + (2) G radial evaluated from the Na + (2) core on the nearest 50 argon atoms G angular,N a + (1) G angular evaluated from the Na + (1) core on the nearest 50 argon atoms G angular,N a + (2)