Orbital rotations

The orbital rotation is a fundamental operation in the simulation of fermions. This page discusses orbital rotations and their properties.

Definitions of orbital rotation

There are two related ways to specify an orbital rotation, by providing either a unitary matrix or an antihermitian matrix. Generally, ffsim represents orbital rotations using unitary matrices.

Unitary matrix definition

An orbital rotation is described by an $N\times N$ unitary matrix $\mathbf{U}$ (here $N$ is the number spatial orbitals), and we denote the corresponding operator as $\mathcal{U}$. This operator has the following action on the fermionic creation operators $\{a_{\sigma ,i}^{†}\}$:

$$\begin{array}{r}\mathcal{U}{a}_{\sigma ,i}^{†}{\mathcal{U}}^{†}=\sum _j{\mathbf{U}}_{ji}{a}_{\sigma ,j}^{†}.\end{array}$$

That is, $a_{\sigma ,i}^{†}$ is mapped to a new operator $b_{\sigma ,i}^{†}$ where $b_{\sigma ,i}^{†}$ is a linear combination of the operators $\{a_{\sigma ,i}^{†}\}$ with coefficients given by the $i$-th column of $\mathbf{U}$. The $\{b_{\sigma ,i}^{†}\}$ also satisfy the fermionic anticommutation relations, so they are creation operators in a rotated basis.

The mapping $\mathbf{U}\mapsto \mathcal{U}$ satisfies the important homomorphism property:

$$\begin{array}{r}\begin{array}{rl}{\mathbf{U}}^{†}& \mapsto {\mathcal{U}}^{†},\\ {\mathbf{U}}_1{\mathbf{U}}_2& \mapsto {\mathcal{U}}_1{\mathcal{U}}_2\end{array}\end{array}$$

for any unitary matrices $\mathbf{U}$, ${\mathbf{U}}_1$, and ${\mathbf{U}}_2$.

Antihermitian matrix definition

Instead of the unitary matrix $\mathbf{U}$, an alternative description of an orbital rotation can be given by providing an antihermitian matrix $\mathbf{K}$. In this case, the orbital rotation is defined as the unitary

$$\begin{array}{r}\mathcal{U}=\prod _{\sigma }\exp \left(\sum _{ij}{\mathbf{K}}_{ij}{a}_{\sigma ,i}^{†}{a}_{\sigma ,j}\right)\end{array}$$

It is a theorem that if $\mathbf{K}=\log (\mathbf{U})$, then the two definitions coincide.

Independent orbital rotations for spin-up and spin-down

The definitions given above assume that the same orbital rotation applies to both the spin-up and spin-down orbitals. However, it is possible for each spin sector to have an independent orbital rotation, in which case a pair of unitary matrices $({\mathbf{U}}^{\alpha },{\mathbf{U}}^{\beta })$ needs to be specified. Letting ${\mathbf{K}}^{\sigma }=\log ({\mathbf{U}}^{\sigma })$, the definitions above then become

$$\begin{array}{r}\begin{array}{rl}\mathcal{U}{a}_{\sigma ,i}^{†}{\mathcal{U}}^{†}& =\sum _j{\mathbf{U}}_{ji}^{\sigma }{a}_{\sigma ,j}^{†},\\ \mathcal{U}& =\prod _{\sigma }\exp \left(\sum _{ij}{\mathbf{K}}_{ij}^{\sigma }{a}_{\sigma ,i}^{†}{a}_{\sigma ,j}\right).\end{array}\end{array}$$

Application to time evolution by a quadratic Hamiltonian

A quadratic Hamiltonian is an operator of the form (here we restrict the definition to Hamiltonians with particle number and spin Z symmetry)

$$\mathcal{M}=\sum _{\sigma ,ij}{\mathbf{M}}_{ij}{a}_{\sigma ,i}^{†}{a}_{\sigma ,j}$$

where $\mathbf{M}$ is a Hermitian matrix. Time evolution by this Hamiltonian is given by the operator $e^{-i\mathcal{M}t}$, where $t$ is the evolution time. Orbital rotations can be used to implement time evolution by a quadratic Hamiltonian using two different methods.

Method 1: Diagonalize the quadratic Hamiltonian

A quadratic Hamiltonian can always be rewritten as

$$\mathcal{M}=\mathcal{U}\left(\sum _{\sigma ,i}{\lambda }_i{n}_{\sigma ,i}\right){\mathcal{U}}^{†}$$

where the $\{{\lambda }_i\}$ are real numbers called orbital energies, $\mathcal{U}$ is an orbital rotation, and $n_{\sigma ,i}=a_{\sigma ,i}^{†}{a}_{\sigma ,i}$ is the occupation number operator. The $\{{\lambda }_i\}$ and the unitary matrix $\mathbf{U}$ describing the orbital rotation are obtained from an eigendecomposition of $\mathbf{M}$:

$${\mathbf{M}}_{ij}=\sum _k{\lambda }_k{\mathbf{U}}_{ik}{\mathbf{U}}_{jk}^{\ast }.$$

Time evolution by $\mathcal{M}$ can be implemented with the following steps:

• Compute the orbital energies $\{{\lambda }_i\}$ and the orbital rotation matrix $\mathbf{U}$ by performing an eigendecomposition of $\mathbf{M}$.

• Perform the orbital rotation ${\mathcal{U}}^{†}$, which corresponds to the matrix ${\mathbf{U}}^{†}$.

• Perform time evolution by the operator $\sum _{\sigma ,i}{\lambda }_i{n}_{\sigma ,i}$.

• Perform the orbital rotation ${\mathcal{U}}^{†}$, which corresponds to the matrix $\mathbf{U}$.

Method 2: Implement directly as an orbital rotation

The antihermitian matrix definition of an orbital rotation shows that the time evolution operator $e^{-i\mathcal{M}t}$ is itself an orbital rotation defined by $\mathbf{K}=-i\mathbf{M}t$, or $\mathbf{U}=e^{-i\mathbf{M}t}$. Thus, the time evolution operator can be implemented directly as an orbital rotation.

The following code cell demonstrates the equivalence of the two methods by applying them to a randomly generated state vector and checking that the result is the same in both cases.

import numpy as np
import scipy.linalg

import ffsim

rng = np.random.default_rng(12345)
norb = 6
nelec = (3, 3)

# Generate a random Hermitian matrix M to represent the quadratic Hamiltonian
mat = ffsim.random.random_hermitian(norb, seed=rng)

# Generate a random state vector and set the evolution time
vec = ffsim.random.random_state_vector(ffsim.dim(norb, nelec), seed=rng)
time = 1.0

# Apply time evolution via diagonalization
eigs, orbital_rotation = scipy.linalg.eigh(mat)
result1 = ffsim.apply_num_op_sum_evolution(
    vec, eigs, time=time, norb=norb, nelec=nelec, orbital_rotation=orbital_rotation
)

# Apply time evolution directly as an orbital rotation
evolution_mat = scipy.linalg.expm(-1j * time * mat)
result2 = ffsim.apply_orbital_rotation(vec, evolution_mat, norb=norb, nelec=nelec)

# Check that the results match
np.testing.assert_allclose(result1, result2)

Quantum circuit implementation

An orbital rotation can be implemented as a quantum circuit by decomposing it into Givens rotations. Under the Jordan-Wigner transformation, the circuit can be implemented as a dense brickwork pattern of two-qubit rotations, followed by a single layer of single-qubit phase gates. The following code cell constructs and displays a Qiskit circuit for implementing an orbital rotation. Notice that the circuit consists of two independent sub-circuits acting on each spin sector.

from qiskit.circuit import QuantumCircuit, QuantumRegister

norb = 6
nelec = (3, 3)

# Generate a random orbital rotation
orbital_rotation = ffsim.random.random_unitary(norb, seed=rng)

# Create a quantum circuit that implements this orbital rotation
qubits = QuantumRegister(2 * norb, name="q")
circuit = QuantumCircuit(qubits)
circuit.append(ffsim.qiskit.OrbitalRotationJW(norb, orbital_rotation), qubits)

# Decompose the orbital rotation into basic gates before drawing
circuit.decompose().draw("mpl", scale=0.6)

When an orbital rotation is applied to a single electronic configuration (in the qubit picture, a computational basis state), the quantum circuit implementation can be optimized to use fewer gates, resulting in a “diamond” pattern of gates rather than a dense “brickwork” pattern. For additional information, see this paper.

The following code cell constructs and displays an example of the optimized circuit.

# Create a quantum circuit that applies an orbital rotation to the Hartree-Fock state
circuit = QuantumCircuit(qubits)
circuit.append(ffsim.qiskit.PrepareHartreeFockJW(norb, nelec), qubits)
circuit.append(ffsim.qiskit.OrbitalRotationJW(norb, orbital_rotation), qubits)

# Optimize the orbital rotation circuit
circuit = ffsim.qiskit.PRE_INIT.run(circuit)

# Decompose the orbital rotation into basic gates before drawing
circuit.decompose().draw("mpl", scale=0.6)