Double-factorized representation of the molecular Hamiltonian

This page discusses the double-factorized representation of the molecular Hamiltonian.

Double-factorized representation

The molecular Hamiltonian is commonly written in the form

$$H=\sum _{\sigma ,pq}{h}_{pq}{a}_{\sigma ,p}^{†}{a}_{\sigma ,q}+\frac{1}{2}\sum _{\sigma \tau ,pqrs}{h}_{pqrs}{a}_{\sigma ,p}^{†}{a}_{\tau ,r}^{†}{a}_{\tau ,s}{a}_{\sigma ,q}+\text{constant}.$$

This representation of the Hamiltonian is daunting for quantum simulations because the number of terms in the two-body part scales as $N^4$ where $N$ is the number of spatial orbitals. An alternative representation can be obtained by performing a “double-factorization” of the two-body tensor $h_{pqrs}$:

$$H=\sum _{\sigma ,pq}{h}_{pq}^{\prime }{a}_{\sigma ,p}^{†}{a}_{\sigma ,q}+\sum _{k=1}^L{\mathcal{U}}_k{\mathcal{J}}_k{\mathcal{U}}_k^{†}+{\text{constant}}^{\prime }.$$

Here each ${\mathcal{U}}_k$ is an orbital rotation and each ${\mathcal{J}}_k$ is a so-called diagonal Coulomb operator of the form

$$\mathcal{J}=\frac{1}{2}\sum _{\sigma \tau ,ij}{\mathbf{J}}_{ij}{n}_{\sigma ,i}{n}_{\tau ,j},$$

where $n_{\sigma ,i}=a_{\sigma ,i}^{†}{a}_{\sigma ,i}$ is the occupation number operator and ${\mathbf{J}}_{ij}$ is a real symmetric matrix. The one-body tensor and the constant have been updated to accomodate extra terms that arise from reordering fermionic ladder operators.

For an exact factorization, $L$ is proportional to $N^2$. However, the two-body tensor often has low-rank structure, and the decomposition can be truncated while still maintaining an accurate representation. Thus, in practice, $L$ often scales only linearly with $N$, resulting in a much more compact representation that gives rise to more efficient schemes for simulating and measuring the Hamiltonian, which utilize efficient quantum circuits for orbital rotations and time evolution by a diagonal Coulomb operator.

Application to time evolution via Trotter-Suzuki formulas

In this section, we discuss how the double-factorized Hamiltonian can be simulated using Trotter-Suzuki formulas.

Brief background on Trotter-Suzuki formulas

Trotter-Suzuki formulas are used to approximate time evolution by a Hamiltonian $H$ which is decomposed as a sum of terms:

$$H=\sum _k{H}_k.$$

Time evolution by time $t$ is given by the unitary operator

$$e^{iHt}.$$

To approximate this operator, the total evolution time is first divided into a number of smaller time steps, called “Trotter steps”:

$$e^{iHt}=(e^{iHt/r}{)}^r.$$

The time evolution for a single Trotter step is then approximated using a product formula, which approximates the exponential of a sum of terms by a product of exponentials of the individual terms. The formulas are approximate because the terms do not in general commute. A first-order asymmetric product formula has the form

$$e^{iH\tau }\approx \prod _k{e}^{i{H}_k\tau }.$$

Higher-order formulas can be derived which yield better approximations.

Application to the double-factorized Hamiltonian

Using the double-factorized representation, the Hamiltonian is decomposed into $L+1$ terms (ignoring the constant term),

$$H=\sum _{k=0}^L{H}_k,$$

where

$$H_0=\sum _{\sigma ,pq}{h}_{pq}^{\prime }{a}_{\sigma ,p}^{†}{a}_{\sigma ,q}$$

and for $k=1,\dots ,L$,

$$H_k={\mathcal{U}}_k{\mathcal{J}}_k{\mathcal{U}}_k^{†}.$$

We have that

• $H_0$ is a quadratic Hamiltonian and can be simulated as described in Orbital rotations and quadratic Hamiltonians, and

• $H_k$ is a diagonal Coulomb operator (up to an orbital rotation) for $k=1,\dots ,L$, and ffsim includes the function apply_diag_coulomb_evolution for simulating time evolution by such an operator.

Given the ability to simulate each of the individual terms, we can implement Trotter-Suzuki formulas for time evolution by the entire Hamiltonian. This is implemented in ffsim by the function simulate_trotter_double_factorized, which implements higher-order Trotter-Suzuki formulas in addition to the first-order asymmetric formula mentioned previously. The first-order asymmetric formula corresponds to setting the argument order=0, which is the default. order=1 corresponds to the first-order symmetric (commonly known as the second-order) formula, order=2 corresponds to the second-order symmetric (fourth-order) formula, and so on.