Is the Jellium model Hamiltonian incorrect?

[petergthatsme]

I have a quick question related to the correctness of the Jellium (plane wave basis) Hamiltonian - in particular the potential term.

A minimal code snippet to generate the potential term of a spinless Jellium Hamiltonian (with the plane wave basis states living on a 2D, 3x3 grid) reads:

from openfermion.hamiltonians.jellium import plane_wave_potential
from openfermion import Grid
import math
grid = Grid(dimensions=2, length=3, scale=2*math.pi)
potential = plane_wave_potential(grid, spinless=True, e_cutoff = None)

There are two potential issues that I see (and I'm curious if I'm missing something maybe?):

1. openfermion tries to implement the second line of Eq 7 in [1]. Looking at plane_wave_potential, the outer loop over \nu (or omega_indices in the code defined on line ~146 of OpenFermion/src/openfermion/hamiltonians/jellium.py), for, say the example above, will iterate over momenta values: [-1. -1.], [-1. 0.], [-1. 1.], [ 0. -1.], [0,0], [0. 1.], [ 1. -1.], [1. 0.], [1. 1.] This seems potentially incomplete to me. I would think there should be other momenta values here like [-2, 0], for example. This is because (see Eq 7 in [1]), orbitals that correspond to index, (for example): \nu + p = [-2, 0] + [1,0] = [-1, 0] are still on the grid of allowable basis states. In other words, the \nu index (in code, effectively related to omega_indices) should include cases corresponding to higher momenta values than what the base 3x3 grid includes.

2. The authors of [1] write eq 7 for a 3d case. The code, however, seems to be written in a general way, and allows other dimensions (say dimension==2). However, plane_wave_potential scales each term of the potential energy by a 1/momenta_squared (see line ~162 of OpenFermion/src/openfermion/hamiltonians/jellium.py) factor, where momenta_squared is calculated from the outer index (i.e index \nu). This seems incorrect for dimensions other than 3. For example, book [2] (see equations 1.19 and 1.68) shows that this scaling should go like 1/momenta (i.e., there is no square) in 2 dimensions.

I'd love to know if I'm missing something or if these are issues with the current Jellium implementation in the code. thanks

[1] Babbush et al, PRX 8, 011044 (2018) [2] Giuliani, G. & Vignale, G. Quantum Theory of the Electron Liquid

[petergthatsme]

Regarding point (1) above: Going through the paper ([1] from above post) in more detail, I am guessing the "missing terms" may have to do with the introduced cutoff on \nu discussed just below eq 7?

I think point (2) still holds though.

[fdmalone]

This sounds like a bug or at least the code should not be this generic. The fourier transform of the Coulomb potential should go like 1/q like you say, not 1/q^2