Orbital rotations and quadratic Hamiltonians¶

This page discusses orbital rotations and how they can be used to implement time evolution by a quadratic Hamiltonian.

Orbital rotations¶

The orbital rotation is a fundamental operation in the simulation of a system of fermionic modes. An orbital rotation is described by an \(N \times N\) unitary matrix \(\mathbf{W}\) (here \(N\) is the number spatial orbitals), and we denote the corresponding operator as \(\mathcal{W}\). This operator has the following action on the fermionic creation operators \(\set{a^\dagger_{\sigma, i}}\):

\[\mathcal{W} a^\dagger_{\sigma, i} \mathcal{W}^\dagger = \sum_j \mathbf{W}_{ji} a^\dagger_{\sigma, j}.\]

That is, \(a^\dagger_{\sigma, i}\) is mapped to a new operator \(b^\dagger_{\sigma, i}\) where \(b^\dagger_{\sigma, i}\) is a linear combination of the operators \(\set{a^\dagger_{\sigma, i}}\) with coefficients given by the \(i\)-th column of \(\mathbf{W}\). The \(\set{b^\dagger_{\sigma, i}}\) also satisfy the fermionic anticommutation relations, so they are creation operators in a rotated basis. The mapping \(\mathbf{W} \mapsto \mathcal{W}\) satisfies the properties

\[\begin{split}\begin{align*} \mathbf{W}^\dagger &\mapsto \mathcal{W}^\dagger, \\ \mathbf{W}_1 \mathbf{W}_2 &\mapsto \mathcal{W}_1 \mathcal{W}_2 \end{align*}\end{split}\]

for any unitary matrices \(\mathbf{W}\), \(\mathbf{W}_1\), and \(\mathbf{W}_2\).

Time evolution by a quadratic Hamiltonian¶

Orbital rotations can be used to implement time evolution by a quadratic Hamiltonian. A quadratic Hamiltonian is an operator of the form (here we only consider Hamiltonians with particle number and spin Z symmetry)

\[\mathcal{M} = \sum_{\sigma, ij} \mathbf{M}_{ij} a^\dagger_{\sigma, i} a_{\sigma, j}\]

where \(\mathbf{M}\) is a Hermitian matrix. A quadratic Hamiltonian can always be rewritten as

\[\mathcal{M} = \mathcal{W} \left(\sum_{\sigma, i} \lambda_i n_{\sigma, i}\right)\mathcal{W}^\dagger\]

where the \(\set{\lambda_i}\) are real numbers which we’ll refer to as orbital energies, \(\mathcal{W}\) is an orbital rotation, and we have introduced the occupation number operator \(n_{\sigma, i} = a^\dagger_{\sigma, i} a_{\sigma, i}\). The \(\set{\lambda_i}\) and the unitary matrix \(\mathbf{W}\) describing the orbital rotation are obtained from an eigendecomposition of \(\mathbf{M}\):

\[\mathbf{M}_{ij} = \sum_k \lambda_k \mathbf{W}_{ik} \mathbf{W}_{jk}^*.\]

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

Compute the orbital energies \(\set{\lambda_i}\) and the orbital rotation matrix \(\mathbf{W}\) by performing an eigendecomposition of \(\mathbf{M}\).
Perform the orbital rotation \(\mathcal{W}^\dagger\), which corresponds to the matrix \(\mathbf{W}^\dagger\).
Perform time evolution by the operator \(\sum_{\sigma, i} \lambda_i n_{\sigma, i}\).
Perform the orbital rotation \(\mathcal{W}^\dagger\), which corresponds to the matrix \(\mathbf{W}\).

This logic can be implemented in ffsim as follows:

[1]:

import numpy as np

import ffsim


def apply_quad_ham_evolution(
    vec: np.ndarray, mat: np.ndarray, time: float, norb: int, nelec: tuple[int, int]
) -> np.ndarray:
    """Apply time evolution by a quadratic Hamiltonian to a state vector."""
    energies, orbital_rotation = np.linalg.eigh(mat)
    vec = ffsim.apply_orbital_rotation(
        vec, orbital_rotation.T.conj(), norb=norb, nelec=nelec
    )
    vec = ffsim.apply_num_op_sum_evolution(vec, energies, time, norb=norb, nelec=nelec)
    vec = ffsim.apply_orbital_rotation(vec, orbital_rotation, norb=norb, nelec=nelec)
    return vec

ffsim already includes a function called apply_num_op_sum_evolution for performing this operation, but it accepts the orbital energies and rotation as arguments rather than the matrix \(M\). The reason for this design is that many applications involve applying time evolution by the same Hamiltonian repeatedly as part of a subroutine. Since the orbital energies and rotation are the same every time, they should be computed only once at the beginning, and then passed to the function that applies the time evolution.