qlat.hlbl_contract — Hadronic Light-by-Light Contractions

Source: qlat/qlat/hlbl_contract.pyx

Note: Update this document when updating the source file.

Outline

  1. Overview

  2. Muon-Line Field

  3. Pair Label Generators

  4. Local Current from Propagators

  5. Point-Selected Probability Fields

  6. CurrentMoments Class

  7. Four-Pair Contraction

  8. Two-Plus-Two Pair Contraction

  9. Examples

Overview

hlbl_contract implements the hadronic light-by-light (HLbL) scattering contraction routines needed for the muon anomalous magnetic moment (g-2) calculation on the lattice. The module provides:

  • Muon-line fields (mk_m_z_field_tag) — the electromagnetic vertex kernel at a point z, used to connect the muon line to the hadronic vacuum polarization.

  • Local currents (mk_local_current_from_props) — the electromagnetic current J_mu(x) constructed from quark propagator pairs.

  • CurrentMoments — moments of the current distribution used in the multipole expansion of the HLbL tensor.

  • Four-pair contractions (contract_four_pair_no_glb_sum) — the dominant HLbL diagram where two photon propagators connect to a single quark loop.

  • Two-plus-two pair contractions (contract_two_plus_two_pair_no_glb_sum) — the subdominant HLbL diagram with two separate quark loops connected by two photon propagators.

All contraction functions return results as NumPy arrays with shape (n_labels, s_limit, l_limit) where s and l are angular momentum quantum numbers in the multipole expansion. Global sums over MPI ranks are not performed inside these functions; the caller must sum the returned arrays.

import qlat as q

q.begin_with_mpi([[1, 1, 1, 1]])
# ... (see examples below)
q.end_with_mpi()

Muon-Line Field

mk_m_z_field_tag(psel_d, xg_x, xg_y, a, tag) -> SelectedPointsRealD

Compute the muon-line field at vertex z for sampled source point x and sink point y.

Parameter

Type

Description

psel_d

PointsSelection

Point selection for the d (detector) sites

xg_x

Coordinate

Global coordinate of source point x

xg_y

Coordinate

Global coordinate of sink point y

a

float

Lattice spacing (in muon-mass units; use muon_mass * a in lattice units)

tag

int

0 = subtraction (preferred), 1 = no subtraction

Returns a SelectedPointsRealD containing the muon-line field values at the selected points.

Pair Label Generators

contract_four_pair_labels(tags: list[str]) -> list

Generate the contraction labels for the four-pair HLbL diagram. Each tag in tags specifies a reference point configuration:

  • "ref-far" — reference point far from the interaction region

  • "ref-center" — reference point at the center

  • "ref-close" — reference point close to the interaction region

contract_two_plus_two_pair_labels() -> list

Generate the contraction labels for the two-plus-two pair HLbL diagram.

Local Current from Propagators

mk_local_current_from_props(sprop1, sprop2) -> SelectedPointsWilsonMatrix

Construct the local electromagnetic current from a pair of point-selected propagators:

J(x) = gamma5 * sprop2^+ * gamma5 * gamma_mu * sprop1

Parameter

Type

Description

sprop1

PselProp

First propagator (forward)

sprop2

PselProp

Second propagator (conjugated with gamma5)

Returns a SelectedPointsWilsonMatrix containing the current at each selected point.

Point-Selected Probability Fields

mk_psel_d_prob_xy(psel_prob, psel_d_prob, idx_xg_x, idx_xg_y) -> tuple

Compute the probability-weighted field for a pair of source/sink points.

Parameter

Type

Description

psel_prob

SelectedPointsRealD

Probability weights for the primary point selection

psel_d_prob

SelectedPointsRealD

Probability weights for the detector point selection

idx_xg_x

int

Index of source point x in psel

idx_xg_y

int

Index of sink point y in psel

Returns (prob_pair, psel_d_prob_xy) where prob_pair is the combined probability for the (x, y) pair and psel_d_prob_xy is the per-detector-point probability field.

CurrentMoments Class

CurrentMoments stores the multipole moments of the electromagnetic current distribution, used in the HLbL tensor decomposition.

Constructor

CurrentMoments()

Create an uninitialized instance.

CurrentMoments(current: SelectedPointsWilsonMatrix, psel_d_prob_xy: SelectedPointsRealD)

Create and compute moments from a local current field weighted by detector-point probabilities.

Copy and Assignment

copy(is_copying_data: bool = True) -> CurrentMoments

Return a copy. If is_copying_data is False, return an uninitialized instance.

__imatmul__(v1: CurrentMoments)

In-place assignment: self @= v1.

Compute Moments

set_from_current(current: SelectedPointsWilsonMatrix, psel_d_prob_xy: SelectedPointsRealD)

Compute moments from the given current and probability fields.

Global Sum

glb_sum() -> CurrentMoments

Return a new CurrentMoments with all entries summed across MPI ranks. Does not modify self.

Four-Pair Contraction

contract_four_pair_no_glb_sum(coef, psel_prob, psel_d_prob, idx_xg_x, idx_xg_y, smf_d, sc_xy, sc_yx, cm_xy, cm_yx, inv_type, tags, r_sq_limit, muon_mass, z_v) -> numpy.ndarray

Perform the four-pair HLbL contraction for a given (x, y) source/sink pair. Returns a 3-D NumPy array of shape (n_labels, s_limit, l_limit).

Parameter

Type

Description

coef

complex

Overall coefficient (default 1.0)

psel_prob

SelectedPointsRealD

Probability weights for primary point selection

psel_d_prob

SelectedPointsRealD

Probability weights for detector point selection

idx_xg_x

int

Index of source point x in psel

idx_xg_y

int

Index of sink point y in psel

smf_d

SelectedPointsRealD

Muon-line field at vertex z

sc_xy

SelectedPointsWilsonMatrix

Local current J(x) for x->y direction

sc_yx

SelectedPointsWilsonMatrix

Local current J(x) for y->x direction

cm_xy

CurrentMoments

Current moments for x->y direction

cm_yx

CurrentMoments

Current moments for y->x direction

inv_type

int

0 = light quark, 1 = strange quark

tags

list[str]

Reference point tags (e.g., ["ref-far", "ref-center"])

r_sq_limit

int

Maximum r^2 cutoff

muon_mass

float

Muon mass

z_v

float

Vector coupling constant

Note: The returned array has not been globally summed; the caller must sum across MPI ranks.

Two-Plus-Two Pair Contraction

contract_two_plus_two_pair_no_glb_sum(coef, psel_prob, psel_lps_prob, idx_xg_x, lps_hvp_x, edl_list_c, r_sq_limit, muon_mass, z_v) -> tuple

Perform the two-plus-two pair HLbL contraction. This diagram involves two quark loops connected by two photon propagators.

Returns (n_points_selected, n_points_computed, lsl_arr) where lsl_arr has shape (n_labels, s_limit, l_limit).

Parameter

Type

Description

coef

complex

Overall coefficient

psel_prob

SelectedPointsRealD

Probability weights for primary point selection

psel_lps_prob

SelectedPointsRealD

Probability weights for the loop point selection

idx_xg_x

int

Index of source point x in psel

lps_hvp_x

SelectedPointsComplexD

HVP tensor at point x

edl_list_c

SelectedPointsComplexD

External double-line list

r_sq_limit

int

Maximum r^2 cutoff

muon_mass

float

Muon mass

z_v

float

Vector coupling constant

Note: The returned array has not been globally summed; the caller must sum across MPI ranks.

Examples

Compute Muon-Line Field

import qlat as q
import os
import tempfile

size_node_list = [[1, 1, 1, 1]]

q.begin_with_mpi(size_node_list)

total_site = q.Coordinate([4, 4, 4, 8])
geo = q.Geometry(total_site)

psel = q.PointsSelection(total_site, [[0, 0, 0, 0]])
psel_d = q.PointsSelection(total_site, [[1, 0, 0, 0],
                                         [0, 1, 0, 0]])

xg_x = q.Coordinate([0, 0, 0, 0])
xg_y = q.Coordinate([2, 0, 0, 0])
a = 0.1  # lattice spacing in muon-mass units
tag = 0  # subtraction (preferred)

# Compute muon-line interpolation data (required by mk_m_z_field_tag).
# In production, use precomputed tables loaded via
# load_multiple_muonline_interpolations. Here we compute a small test
# sample for demonstration purposes.
# Note: compute_save_muonline_interpolation requires at least 2 MPI ranks.
interp_path = os.path.join(tempfile.gettempdir(), "muon-line-interp")
q.compute_save_muonline_interpolation(
    f"{interp_path}/0000000000",
    [3, 2, 2, 2, 2],
    (1e-8, 1e-3, 1024 * 16, 1024 * 1024 * 16),
)

# Set extrapolation weights to use only the single interpolation table
# at index 0. The default weights reference multiple tables at different
# grid resolutions that are not available in this standalone example.
q.set_muon_line_m_extra_weights([[1.0], [1.0]])

smf_d = q.hlbl_contract.mk_m_z_field_tag(psel_d, xg_x, xg_y, a, tag)
print(f"muon-line field shape: {smf_d.n_points}")

q.end_with_mpi()

Local Current from Propagators

import qlat as q

size_node_list = [[1, 1, 1, 1]]

q.begin_with_mpi(size_node_list)

total_site = q.Coordinate([4, 4, 4, 8])
geo = q.Geometry(total_site)

psel = q.PointsSelection(total_site, [[0, 0, 0, 0]])

sprop1 = q.PselProp(psel)
sprop2 = q.PselProp(psel)
sprop1.set_zero()
sprop2.set_zero()

scf = q.hlbl_contract.mk_local_current_from_props(sprop1, sprop2)
print(f"local current computed at {scf.n_points} points")

q.end_with_mpi()

Four-Pair Contraction Workflow

import qlat as q
import numpy as np

size_node_list = [[1, 1, 1, 1]]

q.begin_with_mpi(size_node_list)

total_site = q.Coordinate([4, 4, 4, 8])
geo = q.Geometry(total_site)

# Create minimal selections
psel = q.PointsSelection(total_site, [[0, 0, 0, 0],
                                [1, 0, 0, 0]])
psel_d = q.PointsSelection(total_site, [[0, 0, 0, 0]])

# Probability weights (uniform for demonstration)
psel_prob = q.SelectedPointsRealD(psel, 1)
psel_d_prob = q.SelectedPointsRealD(psel_d, 1)
psel_prob.set_zero()
psel_d_prob.set_zero()

# Generate pair labels
tags = ["ref-far", "ref-center"]
labels = q.hlbl_contract.contract_four_pair_labels(tags)
print(f"number of labels: {len(labels)}")

# Generate two-plus-two labels
labels_22 = q.hlbl_contract.contract_two_plus_two_pair_labels()
print(f"number of 2+2 labels: {len(labels_22)}")

q.end_with_mpi()

CurrentMoments

import qlat as q

size_node_list = [[1, 1, 1, 1]]

q.begin_with_mpi(size_node_list)

total_site = q.Coordinate([4, 4, 4, 8])
geo = q.Geometry(total_site)

psel = q.PointsSelection(total_site, [[0, 0, 0, 0]])
psel_d = q.PointsSelection(total_site, [[0, 0, 0, 0]])

# Build a current and probability field
sprop1 = q.PselProp(psel)
sprop2 = q.PselProp(psel)
sprop1.set_zero()
sprop2.set_zero()

scf = q.hlbl_contract.mk_local_current_from_props(sprop1, sprop2)
psel_d_prob = q.SelectedPointsRealD(psel_d, 1)
psel_d_prob.set_zero()

cm = q.hlbl_contract.CurrentMoments(scf, psel_d_prob)
cm_glb = cm.glb_sum()
print("CurrentMoments computed and globally summed")

q.end_with_mpi()