Collisional S-Matrix for the Vibrational Dynamics of H+H2 by Quantum Computing

An algorithm and a system of quantum circuits is developed and applied to compute accurately the S matrix for the transitions between vibrational states of H2 for collisions with H. The algorithm was applied to 100 eV laboratory collision energy at a quantum circuit simulator. The effects of the discretized dissociative continuum to the transition cross sections are carefully studied and accuracy and convergence of the results with the chosen parameters of the algorithm and the collision system are verified by comparison with a solution of the time-dependent Schrodinger equation using the classical algorithm as well as comparison with a few results from the literature.


INTRODUCTION
A Hydrogen atom colliding with a hydrogen molecule is the most fundamental neutral atommolecule collision systems. Detailed knowledge of the inelastic and elastic processes is of high importance in the star formation processes from molecular clouds, but also in energy and particle redistribution of fusion divertor plasma, especially when hydrogen molecules are vibrationally excited.
Vibrationally excited molecules have significantly enhanced dissociation. Vibrationally excited states of H2 are formed by associative desorption, collisional neutralization, and excitation on the plasmafacing material surfaces, and possibly, by the three-body diatomic association in H+H+H collision. In addition, the elastic scattering at excited vibrational states of abundant hydrogen molecules can play important role in the transport and dissipation od of the divertor plasma momentum [1][2][3][4][5][6].
A typical collision event evolves through dynamically coupled electronic, vibrational, and rotational degrees of freedom. Historically, these systems have been studied thoroughly only for the processes from the ground states (electronic, vibrational, rotational) as motivated by applications, and as well as limited by experimental and theoretical capabilities of the time [3]. Also, the corresponding calculations were usually done within a manifold of bound vibrational states, thus neglecting the possible importance of inelastic processes through dissociative channels, as well as the dynamic change of the dissociative continuum edge with the relative position of the projectile and target. Importance of the effects of dissociative continuum to the transitions between the bound vibrational and electronic states were documented in [5,6].
Paul Benioff and Richard Feynman [7,8] have proposed the idea of quantum computing in the 1980s, predicting the advantages of quantum algorithms in a quantum computer for simulating the complex quantum mechanical systems. Classical computing faces difficulties in, for example, eigenvalue problems of many-body quantum systems, and transition dynamics among excited states.
Despite all the challenges and difficulties to build a satisfying fault-free quantum computer in the near future, an increasing amount of effort is invested to develop quantum algorithms with this new perspective. However, the quantum computing algorithm approach for calculation of a full S matrix for the transition dynamics between bound and continuum states is still mussing.to solve the vibrational dynamics of a system in collision with the inclusion of dissociative continuum is still missing.
In this work, we aimed to build develop an algorithm and develop the system of quantum circuits for a quantum computer that would be able to calculate the evolution operator of the system, from which we obtain the -matrix [5,6,8] and relevant cross-sections for the vibrational transition between dissociating vibrational states at the ground electronic potential curve of H2 ( 1 + 1 ) in collision with H(1s) . be able to calculate the collisional S-matrix [5,6,9] for transitions between dissociating vibrational states at the ground electronic surface of in collision with H(1s). A peculiarity of the present calculation is the inclusion of the dissociative continuum in the calculation. The finite set of the bound vibrational states can be considered as a closed quantum system in presence of the interactions with the dissociative "environment", in which case the calculated matrix will be for an open quantum system [16]. Thus, the current calculation is a possible test ground (not pursued here) for quantum computing methods for the open systems approximations in the Lindblad type of equations.
In Section 2 we define our method and prepare the algorithm to run at a quantum computer and quantum circuits for calculation of the matrix for the vibrational dynamics of H2 in collision with H, including the dissociative continuum. In Section 3 we test the algorithm by simulation using quantum circuits to calculate show the S matrix and results for the relevant transition cross sections of transitions with bound vibrational states, as well as for the dissociation. These results are obtained from simulation using quantum circuits and verified by comparison with benchmark. The benchmark is obtained by solving time dependent Schrodinger equations as system of coupled ordinary differential equations (ODE), using the expansion of the H+H2 wavefunction in the chosen vibrational basis of pseudostates, Finally, our conclusions are presented in Section 4.
Atomic system of units, ℏ = = = 1 , is used throughout the paper unless specified differently.

METHOD
The H2 ( 1 + 1 ) together with H(1s), for various distances of the H atoms, constitutes the ground electronic surface of H3 molecule. For the H3 ground surface, we adopt the analytical fit of Boothroyd [10], augmented with HF-CI results of Krstic et al [1,2] at distances between any two atoms smaller than 1 a.u. The evolution operator of a quantum dynamic system is commonly approximated by trotterization method [17,18], as explained in detail in Section 2.3. For a basis for expansion of the exponential operators we use eigenstates of the unperturbed vibrational Hamiltonian 0 of the where is the vibrational coordinate and  is the reduced mass of H2 (~0.5 amu). This eigenvalue problem was solved with finite quantization "volume" boundary conditions, pseudo-states. Although the relevant vibronic continua are discretized, with this choice of max the density of continuum states stays high enough even for several eV above the continuum edge. The discretization of the dissociative continuum was used earlier by Krstic et al [5,6] to calculate dissociation in the H3 + collision system at collision energies below 10 eV.
The excitation energy of the first vibrational state of H2 is ~0.5 , corresponding to a characteristic vibration time exceeding 50 a.u., which is comparable to the collision time of the order of ( 2 ) 1 2 ⁄ for energies < 100 . The collision time is still short compared to the time scale of molecular rotations (excitation energy 0.01 , i.e., the characteristic times >1000 a.u.), thus enabling one to consider the direction of the diatomic internuclear axis of H2 as fixed. The consequence is that the angle  , defined as the angle between ⃗⃗ (defined from the center of mass of the molecule to the projectile nucleus) and diatomic internuclear axis ⃗, stays constant during a collision, i.e., enters the theory as a parameter. In effect, the equations of motion are completely decupled as far as the angular variables are concerned. For the purposes of testing of our algorithm, we adopt in this work  = 90°, as a typical representant of the atom-molecule collisional configuration, sufficient for the algorithm testing purposes. We note that for accurate calculations of an atom-diatomic molecule collision in a frozen rotation approximation, one would need to average the cross sections for various  in range [0, ] rad.
We use Jacobi coordinates ( ⃗⃗ , ⃗, ) for the collision geometry in this work.
The fully quantal approach, that would besides the electronic motion also quantize the motion of the hydrogen projectile is not practical at higher energies (few tens eV) because this would involve too large number of partial waves. An alternative is to replace the partial wave with the impact To enable the quantum approach for solving this type of problems, Eq. 3 is firstly solved on a numerical mesh, employing the split-operator technique [19][20] in the energy representation. Evolution of the state vector in time from some large initial "− ": which can be iterated with a small timestep , starting from 0 = (− ), and following the ordering of time, as where = − + . Since is moving classically, = ( ), and = √ 2 + 2 2 = 20 is chosen large enough so that is not interacting at = ± with 2 . Then (− ) = 0 is an unperturbed initial vibrational state, which could be any bound ( ) or a linear combination of ( )'s.
In this paper we do not consider association of two atoms of H in a presence of the H projectile and dissociative state of H2 is not used as an initial state in the calculation of the S matrix. τ is chosen small enough ( = 0.01) to have desired transition probabilities |⟨ | ( )⟩| 2 from initial (− ) deviate from the benchmark ODE solution less than a given error , which will be discussed in detail in Sec. 3.

Preparing qubit Hamiltonian
With the discrete -state vibrational basis {| ⟩} , defined by Eq. 2, one can write the Hamiltonian operator as: where the matrix elements of are defined by One can use the second quantization of the Hamiltonian, which replaces the excitation operator | ⟩⟨ | from state to state with the corresponding operator of creation and annihilation of a population of the vibrational states, i.e., with + . Such procedure would require N qubits, which would require significant quantum (and also classical) computing resources when N=256. In this work, we utilize the Qubit Efficient Encoding (QEE) [11] method to map the vibrational Hamiltonian to qubit Hamiltonian. With the QEE mapping, a -state system can be encoded using = log 2 = 8 qubits, making full use of all the qubit computational basis, and establishing obvious quantum advantage of use of the computational resources. With the chosen vibrational basis set {| ⃗⃗⟩}, the ℎ state of the system | ⟩ is encoded to the ℎ qubit computational basis of qubits | ⟩, which can be expressed in individual qubits | 〉 = | n−1 , … , 0 〉 with ∈ {0,1}. Hence one can write the qubit Hamiltonian as: | ⟩⟨ | is now the qubit excitation operator which enables the transition from state | ⟩ to state| 〉: (| ⟩⟨ |)| 〉 = | 〉 . It can be further factorized to individual qubits with the expression: ⊗ , where is the index of qubits, which can be replaced by Pauli operations { , , , } to obtain a full qubit Hamiltonian: After simplification, the qubit Hamiltonian can be written as a linear combination of complex coefficients with Pauli quantum operations: where is the linear combination of ℎ and ℎ coefficients and ∈ { , , , } ⨂ represent the exponential correlation with the number of qubits used to encode the system: = 4 , which is the number of total combinations of -fold Pauli operators.

Calculation of the matrix elements of Hamiltonian
The basic building elements of our algorithm are matrix elements of the total Hamiltonian , defined in Eq. 2, using extended vibrational basis, defined by Eq. 1 where ( = 0, … , − 1) are eigenenergies of the vibrational Hamiltonian obtained as described   In order to obtain ℎ ( ( )) at all time steps, the obtained ℎ ( ) were fitted to the desired times in the time evolution using fourth order splines. We note that the evolution time of the system was varied from = − to = , where = ( ) = √ max 2 − 2 ⁄ and max = 20. Thus, for = 0.01, we find that ( max ) < 10 −6 . The fitting to the desired time instances was done for each impact parameter, used in calculation of the cross sections.

Trotterization of the evolution operator
To simulate the time dependent evolution on a quantum computer, the evolution operator at time , ( , − ) = − ( ( )) , needs to be applied for each timestep , strictly following the time ordering. For simplicity, ( , − ) will be written as ( ) from now on. Then, following Eq. 5b, one can combine ( ) at each discrete timestep together, and the overall time evolution operator is obtained as where − is the initial time, is final time and 's are the discrete time points.
We note that ( 1 ) and ( 2 ) do not commute if 1 ≠ 2 . However, in order to calculate ( ) one needs to factor it into individual Pauli operators. In general, the Pauli terms in a qubit Hamiltonian do not commute with each other. Hence the evolution operator − . A proper approximation of the evolution operator is utilized here by Trotterization, i.e., by Suzuki-Trotter approximation [21]. This approximation splits the exponential operator in a time-loaded product of simpler exponentials operators, which can be calculated for example using their eigenstates. The accuracy of the approximation depends on the small timestep τ in the evolution procedure, as well as on the number of factors used to simplify the exponential operator. Thus, for qubit Hamiltonian with terms, the second order Suzuki-Trotter gives [21]: Each term in Eq. 14, − 2 , can be applied by quantum circuits, discussed in detail in Subsection 2.4.
Applying trotterization of a higher order can improve the accuracy, but at the cost of more terms leading to much longer computational time [12]. The choice of trotterization order should be made with careful tradeoff.

Quantum circuit for simulating time evolution
The operator − 2 in Eq. 14 with ∈ { , , , } ⨂ as the -fold tensor product of Pauli operators and can be applied using quantum circuits. However, the exponentiation of the -fold is a complex quantum operation which requires a further decomposition to the elementary qubit operations such as and the single qubit gates, which can be executed directly by quantum computers. The decomposition is done in two steps.
where and are corresponded Pauli matrices. If , = or , = or and no single qubit operation will be added on ℎ qubit. The expression of the rotational exponential operators, like in Eq.
15, in terms of a finite matrix is an advantageous consequence of the qubit encoding of the Hamiltonian.
For example, in a 4-qubit case the operator − 2 ( ⊗ ⊗ ⊗ ) can be decomposed to a circuit corresponding to the operator − 2 ( ⊗ ⊗ ⊗ ) and several single qubit gates, shown in Fig. 3.
This can be easily achieved by applying one ( ) and connecting neighboring ℎ and ℎ qubits by gate at both sides of ( ) if = = . If = , it will be omitted by gates. An example of the circuit for decomposing the operator − 2 ( ⊗ ⊗ ⊗ ) is shown in Fig. 4. Fig. 4. An example of decomposing − 2 ( ⊗ ⊗ ⊗ ) to a circuit of gate and gates.
After two steps of decomposition, the operator − 2 is a quantum circuit ready to be applied on a quantum computer. The evolution operator ( ) = − ( ( )) at time , which is after trotterization a product of many − ( ( )) 2 terms, can be expressed by the quantum circuit.
Furthermore, by grouping all evolution operators for each time step as a product following Eq. 13, the quantum circuit of overall evolution operator is obtained and shown in Fig. 5. With initial state defined, the probability distribution of the final states after time evolution can be obtained upon measurements, which are used in calculations in Sec. 3.

Calculation of the -matrix
The -matrix is a unitary matrix that relates the initial and final states of a physical system which undergoes a scattering process. It is important to stress that initial and final states of the system related by the S-matrix are defined in absence of the interaction which causes the transitions. Each column of the -matrix contains transition amplitudes from one initial state to all possible final states.
If the initial state of the system is a mixture of populations of various eigenstates of the system, the final state will be described by the coherent combination of various -matrix elements, implying the interference of the complex transition amplitudes. In this work, the S matrix element for transition from the state | ⟩ to the state | ⟩ is: where = − + , is the number of time steps and is the unitary operator of the overall time evolution.
The S-matrix operator is expressed as: which in quantum circuit representation is equivalent to the overall time evolution operator introduced in Subsection 2.4. Following qubit efficient encoding, the -matrix operator can be written as: where | ⟩ and | ⟩ are encoded qubit computational basis states and each element of the -matrix is In case that initial state is a pure vibrational state, the method, shown in Section 2.4 can result to only | | 2 i.e., to probabilities from a state | ⟩ to all states | ⟩. In this section, we propose an original idea of a quantum module to calculate the full -matrix on a quantum computer, including both phases and amplitudes, using a quantum computer, which can be considered as a variation of the Hadamard test [22]. The quantum circuit to calculate the real part of expectation value ⟨ | | ⟩ is shown in Fig.   6. It requires two quantum registers: one register with qubits (e.g., for 256-state system, = 8) enough to load the quantum states | ⟩ and | ⟩, and another register with an ancilla qubit. At the end, only ancilla qubit is measured to give the value of [⟨ | | ⟩]. At the first step, the quantum circuit is initialized in the state: The controlled-operation is applied from ancilla to other qubits, resulting in the state 1 √2 (|0⟩ ⊗ | ⟩ + |1⟩ ⊗ | ⟩) (22) Note that is the overall evolution operator.
After a Hadamard gate applied to the ancilla qubit in Fig. 6 the quantum state becomes: The probability of measuring ancilla qubit in state |0⟩ is:

RESULTS
The test simulations of the algorithm in Section 2 are carried out following the quantum algorithm shown in Fig. 6 and the results are obtained by using gate-based simulator Pennylane [15].
We calculated the transition probabilities for specific cases of initial and finals states, | ⟩ and | ⟩, for various values of impact parameter using quantum circuits shown in Fig. 5 We also calculate the cross sections for dissociation of a bound vibrational state | ⟩ by summation of the probabilities of transition from | ⟩ to all included states of dissociative continuum | ⟩: It is important to investigate the convergence of the bound-bound cross sections with number of included states of discrete dissociative continuum. By increasing stepwise the number of continuum states from 2, we find that this number is max = 113 . The − ( ) for two typical cases of the initial state | ⟩, = 0 and = 7 are shown in Fig. 7a and 7b. The max states span positive energy from 0 to ~5 eV. In Fig. 7c we show ( ) for a few typical bound-bound transitions, calculated with max continuum states. We note that unitarity of the matrix has been preserved in all calculations (to the 6 digits).      One more factor of uncertainty deserves attention: the deviations of the results obtained using quantum circuit simulation and the benchmark results obtained by high accuracy ODE solver, as explained in Subsection 2.6. Fig. 12 shows maximal relative deviation (%, in absolute value) of the final state probabilities after evolution between these two families at = 0.01. Each symbol in Fig. 12 is obtained by taking the maximal value of the relative deviations taking into account of the probabilities of all states. The datapoints are also collected from different sizes of the problem ( = 15 to 128).
Interestingly, we find that the main contributor to the deviation (maximal relative deviation) is the ground vibrational state, which contributes more than 5% to the uncertainty at = 128. We also provided the maximal relative deviation (%) of the final state probabilities excluding the ground state data, shown in Fig. 12 as circle symbol. All other states contribute to case of calculation with 128 states at around 1.5%.

CONCLUSIONS
The algorithm for a quantum computer and corresponding quantum circuits were developed to

S2. S-matrix calculation by quantum algorithm
In this section, we show a full -matrix, with both its real and imaginary parts, for bound-bound transitions, obtained for the H+H2 transitions in a set of 15 bound vibrational states. The simulation is carried out following the quantum algorithm shown in Fig. 6 and the results are obtained from gatebased simulator, Pennylane [15].