Quantum computing for atomic and molecular resonances

Resonances are intermediate or quasi-stationary states that exist during unique atomic processes such as when an excited atom autoionizes, an excited molecule disassociates unimolecularly, or a molecule attracts an electron and then the ion disassociates into stable ionic and neutral subsystems.¹ The characteristics of resonances, such as energy and lifetime, can be revealed by experiments or predicted by theory. One theoretical method to compute properties associated with such resonances is called the complex-scaling method, developed in Refs. 2–7. This method is based on the Balslev–Combes theorem, which is valid for dilation-analytic potentials and can be extended for non-dilation-analytic potential energies.⁸ Additionally, several variants have been developed to study problems such as Stark resonances^9–11 induced by an external electric field. The real space extension of this method uses standard quantum chemistry packages and stabilization graphs.¹² Its main applications are to study the decay of metastable states existing above the ionization threshold of the Li center in open-shell systems such as LiHe,¹³ in the computation of transition amplitudes among metastable states,¹⁴ and in explaining Autler–Townes splitting of spectral lines.¹⁵ 

The complex-scaling method usually requires a large basis set to predict resonances with good accuracy. For example, the helium ¹S resonance uses 32 Hylleraas type functions for basis construction,¹⁶ and the $H2−$ ²$Σu+(σg2σu)$ resonance takes a total of 38 constructed Gaussian atomic bases.⁸ The computational overhead will become overwhelming if more basis functions need to be considered, such as when simulating larger molecular systems or requiring higher accuracy. Moreover, dimensional scaling and large-order dimensional perturbation theory have been applied for complex eigenvalues using the complex-scaling method.^17,18 As for bound states,^19–23 quantum computing algorithms can overcome the above computational limitation problem for resonances. However, most algorithms cannot be directly adapted to resonance calculation with the complex-scaling method because the complex-rotated Hamiltonian is non-Hermitian. For example, the propagator $e−iH(reiθ)t$ in the conventional phase estimation algorithm (PEA) with trotterization²⁴ will be non-unitary, and it cannot be implemented in a quantum circuit directly. In this way, a quantum algorithm for resonance calculation that can work with non-Hermitian Hamiltonians is needed. Daskin et al.²⁵ proposed a circuit design that can solve complex eigenvalues of a general non-unitary matrix. The method applies the matrix rows to an input state one by one and estimates complex eigenvalues via an iterative PEA process. However, for molecular Hamiltonians, the gate complexity of this general design is exponential in system size. In our previous publication,²¹ we briefly mentioned that our direct measurement method can solve complex eigenvalues of non-Hermitian Hamiltonians with polynomial gates. This study extends the direct measurement method and applies it to simple molecular systems as benchmark tests to obtain resonance properties. In particular, we will use IBM’s Qiskit²⁶ simulators and their quantum computers to calculate these resonances.

In Secs. II–V, we first show how to obtain the complex-scaled Hamiltonian for molecular systems and transform it into the Pauli operator form. Then, we introduce the direct measurement method that can derive the Hamiltonian’s complex eigenvalues. Finally, we apply this method to do resonance calculation for a simple model system and a benchmark test system $H2−$ using simulators and IBM quantum computers.

This section presents the steps needed to convert the complex-rotated Hamiltonian to a suitable form that can be simulated on a quantum computer. In the Born–Oppenheimer approximation, the electronic Hamiltonian of a molecular system can be written as a sum of electronic kinetic energy and potential energy of the form

where Z_σ is the σ_th nucleus’s charge, R_σ is the σ_th nucleu’s position, and r_i and r_j represents the i_th and j_th electrons’ position. The complex-scaling method is applied to the study of molecular resonances within the framework of Born–Oppenheimer approximation. Following Moiseyev et al.,²⁷ the electronic coordinates are dilated independently of the nuclear coordinates. Given such a Hamiltonian H(r) in Eq. (1), where r represents electrons’ coordinates, the complex-scaling method rotates r into the complex plane by θ, r→re^iθ. Thus, the Hamiltonian becomes H(re^iθ). After a complex rotation by θ, each electron’s position r becomes r/η, where η = e^−iθ, and thus, the new Hamiltonian from Eq. (1) becomes

It is shown that the system’s resonance state’s energy E and width $Γ=1τ$⁠, where τ is the life time, are related to the corresponding complex eigenvalue of H(re^iθ),^3,28

When doing exact calculations in an infinite basis limit, E_θ in Eq. (5) is not a function of θ. However, there would be dependence in reality because only a truncated basis set is always used in practice. The best resonance estimate is when the complex energy E_θ pauses or slows down in its trajectory^28,29 in the (E_θ, θ) plane or $dEθdθ=0$⁠. In this way, E and Γ can be obtained by solving the new Hamiltonian’s eigenvalues for θ trajectories and looking for the pause. A scaling parameter α is commonly used in the complex rotation process to locate better resonances, which makes η = αe^−iθ. We refer the readers to the book on non-Hermitian quantum mechanics by Moiseyev for more details and method applications.²⁷ 

After choosing a proper orthogonal basis set {ψ_i(r)}, the Hamiltonian can be converted into a second-quantization form,

In the equation, $ai†$ and a_i are fermionic creation and annihilation operators. The coefficients h_ij and h_ijkl can be calculated by

With the Jordan–Wigner transformation,³⁰ 

in which X, Y, and Z are the Pauli operators and

and the Hamiltonian in Eq. (6) will be further transformed into Pauli operators as

In the summation, c_i represents a complex coefficient and P_i represents a k-local tensor product of Pauli operators, where k ≤ n and n is the size of the basis set. Alternatively, the Bravyi–Kitaev transformation³⁰ or parity transformation can also be used in the final step for obtaining the Hamiltonian in the qubit space.

The above process is the same as the conventional Hamiltonian derivation in quantum computing for electronic structure calculations of bound states.^19,31–34 Here, for resonance calculations, to make the Hamiltonian more compatible with the direct measurement method, we rewrite Eq. (10) as

where n_a = ⌈log₂ L⌉. The coefficient β_i and the operator V_i are determined in the following ways:

The direct measurement method is inspired by the direct application of the phase estimation algorithm³⁵ as briefly discussed in our previous publication.²¹ Here, the basic idea is to apply the complex-rotated Hamiltonian to the state of the molecular system and obtain the complex energy information from the output state. Since the original non-Hermitian Hamiltonian cannot be directly implemented in a quantum circuit, this direct measurement method embeds it into a larger dimensional unitary operator.

Assuming n spin orbitals need to be considered for the system, the direct measurement method requires n_s = n qubits to prepare the state of the model system $ϕrs$ and an extra n_a ancilla qubits to enlarge the non-Hermitian Hamiltonian to be a unitary operator. The quantum circuit is shown in Fig. 1.

The B and V gates in the circuit are designed to have the following properties:

which means B transforms the initial ancilla qubits’ state to a vector of coefficients and V applies all V_i on system qubits based on ancilla qubits’ states. One construction choice for B could be implementing the unitary operator

As for V, a series of multi-controlled V_i gates will do the work. If $ϕrs$ is chosen as an eigenstate and we apply the whole circuit of B, V, and B^†,

on it, the output state will be

where Ee^iφ (E ≥ 0) is the corresponding eigenvalue and |Φ^⊥⟩ is a state whose ancilla qubits’ state is perpendicular to $0a$⁠. Then, we can derive E by measuring the output state. To obtain the phase φ, we apply a similar circuit for $Hθ′=xI⊗n+Hθ$⁠, where x is a selected real number, and perform the measurements. The calculation details are found in  Appendix C.

In this section, we calculate the resonance properties of a model system using the direct measurement method. This system is the following one-dimensional potential:²⁸ 

Parameters are chosen as λ = 0.1 and J = 0.8. The potential plot is in Fig. 2. This potential is used to model some resonance phenomena in diatomic molecules. We only consider one electron under this potential. The original Hamiltonian and the complex-rotated Hamiltonian can be written as

To make the setting consistent with the original literature, η is chosen to be e^−iθ and the scaling parameter α is embedded in n Gaussian basis functions

The {χ_k(α)} basis set is not orthogonal, so we apply the Gram–Schmidt process and iteratively construct an orthogonal basis set {ψ_i} as follows:

Since there is only one electron, we do not consider spin interactions. This {ψ_i} basis set is used in the second-quantization step to get the final Hamiltonian in the Pauli matrix form. The resonance eigenvalue found in Ref. 28 with n = 10 basis functions is E_θ = 2.124 − 0.019i hartree. We will try to get the same resonance by applying the direct measurement method using the Qiskit package. The Qiskit package supports different backends, including a statevector simulator that executes ideal circuits, a QASM simulator that provides noisy gate simulation, and various quantum computers. In what follows, we show the results when the basis function number is n = 5 and n = 2. In particular, the former n = 5 case shows how θ trajectories locate the best resonance estimate, and the latter n = 2 cases show how to further simplify the quantum circuit for the direct measurement method and run it on IBM quantum computers.

C1 in Table I is our primary example where we follow the above steps in Secs. II and III for n = 5. An example of the complex-rotated Hamiltonian is shown in  Appendix A. Figure 3 shows a sweep of scaling parameters α for statevector simulations of θ trajectories. Most trajectories pause around the point, E_θ = 2.1265 − 0.0203i hartree, when α = 0.65 and θ = 0.160. Based on Eq. (5), this indicates that the resonance energy and width are E = 2.1265 hartree and Γ = 0.0406 hartree, respectively, close to the resonance energy from Ref. 28 The IBM quantum computer cannot perform the method due to a large number of standard gates in the circuit. Instead, we used the QASM simulator for 4 * 10⁴ shots and obtained the system’s resonance energy at α = 0.65, θ = 0.160, and E_θ = 2.1005 − 0.3862i hartree. This result has an error of around 0.3 hartree but can be augmented by more sample measurements.

When taking n = 2 for the basis function, we are not able to locate the best resonance estimate (see Fig. 3) based on direct diagonalization. Hence, we only use the direct measurement method to calculate the complex eigenenergy when α = 0.65 and θ = 0.160, where the best location is at n = 5. We run the direct measurement method using simulators first and then try to reduce the number of ancilla qubits to make the resulting circuit short enough to be executed in the IBM quantum computers.

C2 in Table I is the case when we follow the steps for n = 2 in Secs. II and III. The Hamiltonian H_θ and how to calculate its complex eigenvalue are shown in Appendix D 1 [Eq. (D1)]. Figure 4 gives the quantum circuit for H_θ. This circuit can be executed in simulators with the results listed in Table II.

However, it is too complicated to be successfully run in IBM quantum computers. For C3 in Table I, we simplify the quantum circuit by calculating the complex eigenvalue for the Hamiltonian H_θ in Appendix D 2 [Eq. (D3)]. Because there are only four terms left, two ancilla qubits are enough for the method. The simplified quantum circuit is then shown in Fig. 5. To avoid introducing more ancilla qubits, instead of $Hθ′=Hθ+xII$⁠, we can run a similar four-qubit circuit for $Hθ′=Hθ+Hθ3$⁠, which has the same terms of tensor products as H_θ with different coefficients. This circuit can be executed successfully in the simulators and the IBM quantum computers. However, it costs around 200 gates in the IBM quantum computers, leading to a significant error. The resulting resonance eigenenergies and errors can be seen in Table III.

For the Hamiltonian in Eq. (D3), a simpler circuit can be constructed if we try to calculate the complex eigenvalue of its square [Eq. (D6) in Appendix D 3]. This is C4 in Table I. The quantum circuit for this $Hθ2$ is shown in Fig. 6.

We can also run a similar three-qubit circuit for $(Hθ2)′$ = $Hθ2+Hθ4$⁠. The implementation of the circuit costs nine gates in the IBM quantum computers after circuit optimization. The resulting eigenenergies are in Table IV.

This section presents a proof of concept that by using our quantum algorithm, the direct measurement method, one can calculate molecular resonances on a quantum computer. We focus on the resonances of a simple diatomic molecule, $H2−$ ²$Σu+(σg2σu)$⁠. Moiseyev and Corcoran⁸ showed how to obtain this molecule’s resonance using a variational method based on the (5s, 3p, 1d/3s, 2p, 1d) contracted Gaussian atomic basis, which contains a total of 76 atomic orbitals for $H2−$⁠. They picked around 45 configurations of natural orbitals as a final basis for the resonance calculation. Here, however, we are not going to use their contracted Gaussian atomic basis that needs 76 system qubits with additional ancilla qubits, which is too large to be simulated by classical computers. The number of gates would also be overwhelming. One may try an iterative diagonalization approach to get a few eigenvalues without constructing matrices or vectors. Another possible solution could be using tensor network simulators. Recent studies by Ellerbrock and Martinez show that tensor network simulators are able to efficiently and accurately simulate over 100-qubit circuits with moderate entanglement.³⁶ In another study, Zhou et al. showed that even strongly entangled systems (as those generated by 2D random circuits) can be simulated by matrix product states comparably accurate to modern quantum devices.³⁷ However, building those simulators for our system is beyond the scope of this paper. In this way, we picked small basis sets, 6-31g and cc-pVDZ, for our simulations. We used the Born–Oppenheimer approximation followed by complex rotation, as shown in Sec. II, “COMPLEX-SCALED HAMILTONIAN” and mapped the Hamiltonian to the qubit space, as shown in  Appendix B. We then apply the direct measurement method to the Hamiltonian to obtain complex eigenvalues. An example quantum circuit to run the direct measurement method can be found in  Appendix E.

Figure 7 shows one complex eigenvalue’s θ trajectories at α = 1.00 under different basis sets after running the algorithm. Figure 7(a) is simulated using the 6-31g basis set. Eight spin orbitals are considered in our self-defined simulator, and 16 qubits are needed to run the algorithm. In this case, if we fix η = αe^−iθ at the lowest point in the figure, which has α = 1.00, θ = 0.18, the resonance energy obtained by the direct measurement method is E_θ = −0.995 102 − 0.046 236i hartree. This complex energy is close to that obtained in Ref. 8, E_θ = −1.0995 − 0.0432i hartree, especially the imaginary part. Figure 7(b) is simulated using the cc-pVDZ basis set. We only considered the s and p_z basis functions for H atoms for easier simulation. 12 spin orbitals are considered in our self-defined simulator, and a total of 23 qubits are needed to run the algorithm in quantum computers. The results show that the resonance energy at the lowest point in the figure, which has α = 1.00, θ = 0.22, is E_θ = −1.045 083 − 0.044 513i hartree. This is even closer to that obtained in Ref. 8. However, we want to note that the lowest points in Figs. 7(a) and 7(b) are not pause points. In addition, they do not reveal real resonance properties. Even after shifting different α in simulations, we cannot find a consistent pause point in θ trajectories to locate the best resonance estimation. The reason may be related to our selected basis. Compared with the literature,⁸ our basis set is much smaller and is not optimized for the resonance state. Still, this application gives a proof of concept and shows that one can calculate molecular resonances on a quantum computer. In the future, if more qubits are available in quantum computers, a large basis can be used, and we may be able to show finer structures in trajectories that can locate the best resonance point. In addition, a larger basis set should lead to a more accurate resonance calculation.

In this paper, we construct and show a proof of concept for a quantum algorithm that calculates atomic and molecular resonances. We first presented the complex-scaling method to calculate molecular resonances. Then, we introduced the direct measurement method, which embeds a molecular system’s complex-rotated Hamiltonian into the quantum circuit and calculates the resonance energy and lifetime from the measurement results. These results represent the first applications of the complex-scaling Hamiltonian to molecular resonances on a quantum computer. The method is proven to be accurate when applied to a simple one-dimensional quantum system that exhibits shape resonances. We tested our algorithm on quantum simulators and IBM quantum computers. Furthermore, when compared to the exponential time complexity in traditional matrix-vector multiplication calculations, this method only requires O(n⁵) standard gates, where n is the size of the basis set. These findings show this method’s potential to be used in a more complicated molecular system and for better accuracy in the future when more and better qubit machines are available.

We would like to thank Rongxin Xia, Zixuan Hu, and Manas Sajjan for useful discussions. We also like to acknowledge financial support from the National Science Foundation under Award No. 1955907.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Table V shows an example of the model system's complex-rotated Hamiltonian.

Table VI shows an example of the $H2−$ system's complex-rotated Hamiltonian.

If the output state equation (17) is measured many times, the possibility of obtaining the $0a$ state, p, is related to E by the following equation:

which reveals $|E|=pA$⁠. To obtain the phase, one way is that we apply a similar circuit for $Hθ′=xI⊗n+Hθ$⁠, where x is a selected real number. Then, the updated U_r′ leads us to

By applying $|E|=pA$ to Eq. (C2), we can solve the phase φ and finally the complex eigenvalue as

If we expand the exponential term in Eq. (C3), it becomes

Since the measurement errors for p and p′, i.e., Δ(p) and Δp′, are $O(1N)$⁠, based on Eq. (C4), the error for the complex eigenvalue Ee^iφ is

The larger the sampling size, the more accurate the obtained complex eigenvalues are.

There are also other choices to obtain the phase. For example, instead of adding the I^⊗n part, we can try building U_r′ based on $Hθ+Hθ2$ or $Hθ+Hθ3$ to get an equation such as Eq. (C2) containing phase information. This equation together with Eq. (C1) will reveal the complex eigenvalue for the input eigenstate with another expression.

The complex-rotated Hamiltonian of the model system is

By running the circuit shown in Fig. 4 for H_θ and a similar circuit for $Hθ′=xII+Hθ$⁠, the complex eigenvalue can be derived by

where A and A′ can be obtained from the absolute value of coefficients in H_θ and $Hθ′$ and p and p′ can be obtained from the measurement results.

The complex-rotated Hamiltonian of the model system without the II term is

If we choose $Hθ′=Hθ+Hθ3$⁠, which has the same terms of tensor products as H_θ with different coefficients, by running Fig. 5, the complex eigenvalue for the original Hamiltonian can be represented by

or

where A and A′ can be obtained from the absolute value of coefficients in H_θ and $Hθ′$ and p and p′ can be obtained from the measurement results.

The square of the Hamiltonian in Eq. (D3) is

If we choose $(Hθ2)′=Hθ2+Hθ4$⁠, by running Fig. 6, the complex eigenvalue for the original Hamiltonian is

or

where A and A′ can be obtained from the absolute value of coefficients in $Hθ2$ and $Hθ2+Hθ4$ and p and p′ can be obtained from their measurement results.

The complex-scaled Hamiltonian of $H2−$ at θ = 0.18, α = 1.00 in  Appendix B can be written as

We would like to mention that the terms explicitly shown in Eq. (E1) are following the order in  Appendix B. It is a coincident that their phases are similar. For example, one term we did not show in the Hamiltonian is 0.021 284 * e^1.542696iIIIIYYII, which has a different phase.

To construct the quantum circuit for the direct measurement method, we need to create the B gate and V gate. The B gate can be prepared by the coefficients from the Hamiltonian in Eq. (E1),

as shown in Eq. (15). The V gate can be constructed by a series of controlled-V_i gates, where V_i are

The whole circuit is shown in Fig. 8. The encoding of control qubits is based on the binary form of V_i’s index i. For example, V₃ is applied to $ψs$ if the ancilla qubit state is $3a=00000011a$⁠.