Four-center matrix elements

Projecting an interacting continuum problem onto localized orbitals produces a four-index interaction tensor. In Hubbard-like and multiorbital models, its entries become the direct repulsion, exchange, pair-hopping, correlated-hopping, and longer-range interaction parameters that control the many-body Hamiltonian. Computing that tensor requires the spatial orbital products, not only the orbital densities.

QUARTIC2D evaluates

\[ U_{1234}=\iint d^2\mathbf r\,d^2\mathbf r'\, \phi_1^*(\mathbf r)\phi_2^*(\mathbf r') U(|\mathbf r-\mathbf r'|) \phi_3(\mathbf r)\phi_4(\mathbf r'). \]

The orbitals \(\phi_i\) are localized scalar functions in two dimensions. They may be real or complex. The kernel \(U\) is assumed to be translationally invariant and radial in the in-plane coordinates.

Define two transition fields,

\[ \rho_{13}(\mathbf r)=\phi_1^*(\mathbf r)\phi_3(\mathbf r), \qquad \rho_{42}(\mathbf r)=\phi_4^*(\mathbf r)\phi_2(\mathbf r). \]

The matrix element becomes

\[ U_{1234}=\iint d^2\mathbf r\,d^2\mathbf r'\, \rho_{13}(\mathbf r) U(|\mathbf r-\mathbf r'|) \rho_{42}^*(\mathbf r'). \]

This is the input structure used by Interaction. Neither transition field is required to be positive, real, or equal to the other. The many-body label attached to a matrix element follows from the orbital indices, not from a separate numerical mode in QUARTIC2D.

Common orbital-index choices

Physical term	Four-index element	First transition field	Second transition field
direct / density-density	\(U_{ijij}\)	\(\phi_i^*\phi_i\)	\(\phi_j^*\phi_j\)
exchange	\(U_{ijji}\)	\(\phi_i^*\phi_j\)	\(\phi_i^*\phi_j\)
pair hopping	\(U_{iijj}\)	\(\phi_i^*\phi_j\)	\(\phi_j^*\phi_i\)
correlated hopping	e.g. \(U_{iiij}\)	\(\phi_i^*\phi_i\)	\(\phi_j^*\phi_i\)

The table lists common many-body names. QUARTIC2D receives only the two transition fields, the radial kernel, and the relative displacement; no separate interaction-type switch is required.

Displaced localized states

The two transition fields may be centered at different sites. QUARTIC2D represents their relative in-plane displacement by

\[ \boldsymbol\delta=(\delta_x,\delta_y) =\delta(\cos\phi_\delta,\sin\phi_\delta). \]

Interaction accepts an array of Cartesian vectors with shape (D, 2):

deltas = np.array([
    [0.0, 0.0],
    [1.0, 0.0],
    [0.0, 1.0],
])
print(deltas.shape)

Output:

(3, 2)

Each row is evaluated independently. Angular dependence is retained through the harmonic-pair phase factors, so two vectors with the same magnitude can give different values for anisotropic transition fields.

Kernel convention

The package expects the momentum-space radial kernel \(U(q)\), with

\[ U(\mathbf R)=\int\frac{d^2\mathbf q}{(2\pi)^2} U(q)e^{i\mathbf q\cdot\mathbf R}. \]

Examples include

\[ U(q)\propto \frac{1}{q} \]

for a bare 2D Coulomb form, screened kernels such as Rytova-Keldysh [Keldysh, 1979, Rytova, 1967], gate-screened interactions, interlayer kernels, Yukawa forms, and other scalar radial functions.

The callable must accept NumPy arrays:

kappa = 0.4


def U_q(q):
    return 2.0 * np.pi / np.sqrt(q**2 + kappa**2)

print(U_q(np.array([0.0, 1.0])))

Output:

[15.70796327  5.83379110]

QUARTIC2D does not impose a unit system. The real-space coordinates, momentum coordinates, kernel, and orbital normalization must use mutually consistent units. If real-space coordinates are measured in nanometers, for example, the momentum variable is in inverse nanometers and the kernel must be written in the corresponding convention.

See Fourier and coordinate conventions for the complete Fourier normalization used by the package.

The benchmark suite contains the explicit kernel families and parameter values in (7) and (8). Use those definitions when comparing a production kernel with the numerical regimes discussed in Choosing numerical methods.