qten.topology.z2

Module reference for qten.topology.z2.

z2

Two- and three-dimensional \(\mathbb{Z}_2\) invariants of time-reversal invariant insulators.

A gapped 2-D occupied bundle with even filling carries one Kane--Mele index \(\nu\). A gapped 3-D bundle carries four Fu--Kane indices \((\nu_0; \nu_1\nu_2\nu_3)\). This module evaluates them from a rank-3 Tensor with dims (MomentumSpace, HilbertSpace, HilbertSpace) (or a MomentumBlockSpace of diagonal \((k,k)\) blocks whose momenta form a complete 2-D or 3-D grid).

The input mesh is Fourier-interpolated to a tight-binding hopping tensor so that time-reversal invariant momenta (TRIM) and Wilson-loop strings can be sampled independently of whether those points sit on the original grid.

Core API

• z2_indices Fu--Kane inversion parities at the \(2^d\) TRIM, hybrid-Wannier Wilson loops, or both.
• Z2ParityResult, Z2WilsonResult, Z2CombinedResult Result mappings returned by method="parity", "wilson", and "both".

Mathematical convention

TRIM are the points \(\Gamma_i=n/2\) with \(n\in\{0,1\}^d\). At each TRIM the occupied inversion eigenvalues \(\xi_n(\Gamma_i)=\pm 1\) come in Kramers pairs. Their pair product is

\[ \delta(\Gamma_i)=\prod_{m=1}^{N_{\mathrm{occ}}/2}\xi_{2m}(\Gamma_i) =(-1)^{N_-(\Gamma_i)/2}, \]

where \(N_-\) is the number of occupied negative parities. In two dimensions the returned index satisfies \((-1)^\nu=\prod_{i=1}^{4}\delta(\Gamma_i)\). In three dimensions

\[ (-1)^{\nu_0}=\prod_{i=1}^{8}\delta(\Gamma_i), \qquad (-1)^{\nu_j}=\prod_{\Gamma_i:\,k_j=\pi}\delta(\Gamma_i). \]

Without inversion, hybrid Wannier charge centers give the same invariants. Along a closed string at fixed \(k_\perp\),

\[ W(k_\perp)=\prod_\ell \operatorname{polar}\!\left[ U^\dagger(k_\ell)\,e^{-2\pi i\,\tau\cdot\Delta k}\,U(k_{\ell+1}) \right], \qquad \bar{x}_n(k_\perp)=\frac{\operatorname{Arg}\lambda_n(W)}{2\pi}\bmod 1, \]

where \(U(k)\) holds occupied eigenvectors and \(\tau\) is the orbital fractional offset. The plane invariant is the Soluyanov--Vanderbilt largest-gap crossing count of those centers. In 2-D that single plane is \(\nu\); the two loop orientations should agree. In 3-D each TRIM plane \(k_j=0\) or \(k_j=\pi\) carries a 2-D invariant, and

\[ \nu_0=\nu(k_j=0)+\nu(k_j=\pi)\pmod{2}, \qquad \nu_j=\nu(k_j=\pi). \]

If the three 3-D axes disagree on \(\nu_0\), the majority vote is returned.

Fourier interpolation

Sampled Bloch matrices are placed on the rectangular reciprocal quotient and inverted with an FFT. Evaluation at fractional \(k\) is the trigonometric polynomial

\[ H(k)=\sum_R t(R)\,e^{-2\pi i\,k\cdot R}, \qquad t=\mathcal{F}^{-1}[H_{\mathrm{mesh}}]. \]

The same interpolant is used for an explicit inversion tensor. When inversion is assembled from orbital offsets about a center \(c\),

\[ I_{\alpha\beta}(k)=(I_0)_{\alpha\beta} \exp\bigl(-i[(r_\alpha-c)+(r_\beta-c)]\cdot k_{\mathrm{cart}}\bigr). \]

The periodic cell must be diagonal in the primitive basis.

Numerical methods

• method="parity" uses Fu--Kane products of inversion eigenvalues at the TRIM. It requires an inversion operator: either an explicit rank-3 tensor whose matrices are paired to \(H(k)\) by momentum labels, or orbital Offset labels from which spatial inversion about inversion_center is assembled.
• method="wilson" tracks hybrid Wannier charge centers. It does not need inversion symmetry.
• method="auto" tries parity first and falls back to Wilson loops if inversion cannot be resolved.
• method="both" runs both constructions. Parity must succeed; the returned indices are the parity values.

Z2ParityTrimDiagnostics

Bases: TypedDict

Inversion-parity diagnostics at one time-reversal invariant momentum.

This mapping is one value in Z2ParityResult "diagnostics". TRIM are labeled by bits \(n\in\{0,1\}^d\) with \(k=n/2\).

Attributes:

Name	Type	Description
delta	int	Fu--Kane pair-parity product \(\delta(\Gamma)=\pm 1\). Equal to \((-1)^{N_-/2}\), where \(N_-\) is the number of occupied negative inversion eigenvalues.
parity_eigenvalues	Tensor	Occupied inversion eigenvalues \(\xi_n(\Gamma)\) as a labeled Tensor with dims (IndexSpace(n_occupied),) and shape (n_occupied,).
commutator_error	float	Relative residual \(\|HI-IH\|/\|H\|\) at this TRIM. Large values mean the supplied or assembled inversion does not commute with \(H(\Gamma)\).
direct_gap	float	Occupied-to-empty direct gap at this TRIM. nan if the occupied count leaves no empty band.

delta instance-attribute

delta: int

parity_eigenvalues instance-attribute

parity_eigenvalues: Tensor

commutator_error instance-attribute

commutator_error: float

direct_gap instance-attribute

direct_gap: float

Z2ParityResult

Bases: TypedDict

Result of z2_indices(..., method="parity").

Fu--Kane indices from inversion eigenvalues at the \(2^d\) TRIM. The runtime object is a plain dict; keys below are required.

Attributes:

Name	Type	Description
indices	tuple[int, ...]	\(\mathbb{Z}_2\) indices as integers in \(\{0,1\}\). Length 1 in two dimensions, \((\nu,)\). Length 4 in three dimensions, \((\nu_0, \nu_1, \nu_2, \nu_3)\).
method	{'parity'}	Construction that produced indices.
parity_products	dict[tuple[int, ...], int]	TRIM bit-tuple \(n\) to \(\delta(\Gamma_n)=\pm 1\). Each key has one 0/1 entry per spatial axis, with \(k_j=n_j/2\).
diagnostics	dict[tuple[int, ...], Z2ParityTrimDiagnostics]	Per-TRIM Z2ParityTrimDiagnostics with the same keys as parity_products.
direct_gap	float	Minimum finite occupied-to-empty gap over the TRIM. nan if none of those gaps are finite.

See Also

z2_indices Public constructor of this mapping.

indices instance-attribute

indices: tuple[int, ...]

method instance-attribute

method: Literal['parity']

parity_products instance-attribute

parity_products: dict[tuple[int, ...], int]

diagnostics instance-attribute

diagnostics: dict[tuple[int, ...], Z2ParityTrimDiagnostics]

direct_gap instance-attribute

direct_gap: float

Z2WilsonPlaneResult

Bases: TypedDict

Hybrid-Wannier data on one Wilson-loop plane.

In two dimensions this is one loop orientation over the Brillouin zone. In three dimensions it is one TRIM plane \(k_{\mathrm{normal}}=0\) or \(1/2\).

Attributes:

Name	Type	Description
z2	int	Plane \(\mathbb{Z}_2\) invariant in \(\{0,1\}\), from the Soluyanov--Vanderbilt largest-gap crossing count of the Wannier centers.
wcc	Tensor	Hybrid Wannier charge centers \(\bar x_n(k_\perp)\in[0,1)\) as a labeled Tensor with dims (IndexSpace(n_perp), IndexSpace(n_occupied)) and shape (n_perp, n_occupied). The first axis follows sweep.
gap_pos	Tensor	Largest-gap position on the Wannier circle at each sweep sample, as a labeled tensor with dims (IndexSpace(n_perp),).
sweep	Tensor	Fractional \(k_\perp\) samples from \(0\) to \(1/2\), labeled by the same IndexSpace(n_perp) as wcc and gap_pos.
min_gap	float	Minimum occupied-to-empty direct gap along the Wilson strings on this plane. nan if no finite gap is available.
kramers_resolved	bool	Whether Wannier centers at the TRIM-plane endpoints pair into Kramers partners within kramers_tolerance.

z2 instance-attribute

z2: int

wcc instance-attribute

wcc: Tensor

gap_pos instance-attribute

gap_pos: Tensor

sweep instance-attribute

sweep: Tensor

min_gap instance-attribute

min_gap: float

kramers_resolved instance-attribute

kramers_resolved: bool

Z2WilsonResult

Bases: TypedDict

Result of z2_indices(..., method="wilson").

Hybrid-Wannier \(\mathbb{Z}_2\) indices. The runtime object is a plain dict; keys below are required.

Attributes:

Name	Type	Description
indices	tuple[int, ...]	Same layout as Z2ParityResult.indices: \((\nu,)\) in 2-D or \((\nu_0, \nu_1, \nu_2, \nu_3)\) in 3-D.
method	{'wilson'}	Construction that produced indices.
planes	dict[tuple[int, float], Z2WilsonPlaneResult]	Plane-resolved hybrid-Wannier data. In 2-D the key is (loop_axis, 0.0) for each loop orientation. In 3-D the key is (normal, trim) with trim in {0.0, 0.5}.
axis_z2	tuple[tuple[int, ...], ...]	Per-axis plane invariants. In 2-D each entry is (ν,) for one loop orientation. In 3-D each entry is (ν(k_j=0), ν(k_j=π)).
min_gap	float	Minimum min_gap over planes. nan if none are finite.

See Also

Z2WilsonPlaneResult Value type stored in planes. z2_indices Public constructor of this mapping.

indices instance-attribute

indices: tuple[int, ...]

method instance-attribute

method: Literal['wilson']

planes instance-attribute

planes: dict[tuple[int, float], Z2WilsonPlaneResult]

axis_z2 instance-attribute

axis_z2: tuple[tuple[int, ...], ...]

min_gap instance-attribute

min_gap: float

Z2CombinedResult

Bases: TypedDict

Result of z2_indices(..., method="both").

Both constructions are run. indices follows the Fu--Kane parity values; a mismatch with Wilson emits a RuntimeWarning.

Attributes:

Name	Type	Description
indices	tuple[int, ...]	Copy of parity["indices"].
method	{'both'}	Construction tag for this combined mapping.
parity	Z2ParityResult	Full Fu--Kane Z2ParityResult.
wilson	Z2WilsonResult	Full hybrid-Wannier Z2WilsonResult.

See Also

z2_indices Public constructor of this mapping.

indices instance-attribute

indices: tuple[int, ...]

method instance-attribute

method: Literal['both']

parity instance-attribute

parity: Z2ParityResult

wilson instance-attribute

wilson: Z2WilsonResult

z2_indices

z2_indices(
    bloch_hamiltonian: Tensor,
    n_occupied: int | None = None,
    *,
    method: Literal["auto"] = "auto",
    inversion: Tensor | None = None,
    inversion_center: Offset
    | Sequence[float]
    | None = None,
    n_loop: int = 32,
    n_perp: int = 17,
    parity_tolerance: float = 1e-05,
    kramers_tolerance: float = 0.08,
    gap_tolerance: float = 1e-08,
) -> Z2ParityResult | Z2WilsonResult

z2_indices(
    bloch_hamiltonian: Tensor,
    n_occupied: int | None = None,
    *,
    method: Literal["parity"],
    inversion: Tensor | None = None,
    inversion_center: Offset
    | Sequence[float]
    | None = None,
    n_loop: int = 32,
    n_perp: int = 17,
    parity_tolerance: float = 1e-05,
    kramers_tolerance: float = 0.08,
    gap_tolerance: float = 1e-08,
) -> Z2ParityResult

z2_indices(
    bloch_hamiltonian: Tensor,
    n_occupied: int | None = None,
    *,
    method: Literal["wilson"],
    inversion: Tensor | None = None,
    inversion_center: Offset
    | Sequence[float]
    | None = None,
    n_loop: int = 32,
    n_perp: int = 17,
    parity_tolerance: float = 1e-05,
    kramers_tolerance: float = 0.08,
    gap_tolerance: float = 1e-08,
) -> Z2WilsonResult

z2_indices(
    bloch_hamiltonian: Tensor,
    n_occupied: int | None = None,
    *,
    method: Literal["both"],
    inversion: Tensor | None = None,
    inversion_center: Offset
    | Sequence[float]
    | None = None,
    n_loop: int = 32,
    n_perp: int = 17,
    parity_tolerance: float = 1e-05,
    kramers_tolerance: float = 0.08,
    gap_tolerance: float = 1e-08,
) -> Z2CombinedResult

Compute the 2-D or 3-D \(\mathbb{Z}_2\) indices of an occupied band subspace.

The n_occupied lowest-energy eigenstates, which must form an even number of Kramers pairs, define an occupied bundle over a complete two- or three-dimensional periodic momentum grid. Two numerical methods are available:

• method="parity" evaluates Fu--Kane inversion eigenvalues at the \(2^d\) TRIM. Two dimensions return \((\nu,)\); three dimensions return \((\nu_0; \nu_1\nu_2\nu_3)\).
• method="wilson" computes hybrid Wannier charge centers. In 2-D the two loop orientations should agree on \(\nu\). In 3-D the strong index is \(\nu_0=\nu(k_i=0)+\nu(k_i=\pi)\bmod 2\) and the weak indices are the three \(k_i=\pi\) plane invariants. If the three axes disagree on \(\nu_0\), the majority vote is returned.

The Hamiltonian is Fourier-interpolated from the supplied mesh, so TRIM and Wilson strings need not coincide with sampled \(k\)-points. This is the construction used for odd meshes such as \(27^3\) or \(9^2\). The periodic cell must be diagonal in the primitive basis.

Parameters:

Name	Type	Description	Default
bloch_hamiltonian	Tensor	Rank-3 Hermitian Tensor with dims (MomentumSpace, HilbertSpace, HilbertSpace), or with a MomentumBlockSpace of diagonal \((k,k)\) blocks whose momenta form a complete 2-D or 3-D reciprocal quotient. The last two axes are square Bloch-Hamiltonian matrices and are aligned onto a common Hilbert space.	required
n_occupied	int | None	Number of lowest-energy bands defining the occupied subspace. It must be even and lie strictly between zero and the total band count. Defaults to half the bands using integer division.	None
method	('auto', 'parity', 'wilson', 'both')	Numerical construction. "auto" tries Fu--Kane parity and falls back to Wilson loops if inversion cannot be resolved (RuntimeError). "both" requires parity to succeed. Defaults to "auto".	"auto"
inversion	Tensor | None	Optional rank-3 inversion operator with dims (MomentumSpace, HilbertSpace, HilbertSpace), or with a diagonal MomentumBlockSpace of \((k,k)\) blocks. Band axes are aligned onto the Hamiltonian Hilbert space. Each momentum is paired with \(H(k)\) by its label, not by data-axis order, and the labeled points must form the same complete reciprocal quotient. If omitted, spatial inversion is assembled from orbital Offset labels about inversion_center.	None
inversion_center	Offset | Sequence[float] | None	Fixed point of spatial inversion, as an Offset or a \(d\)-vector in the Hamiltonian's direct-lattice coordinates. Defaults to the centroid of the unique orbital offsets.	None
n_loop	int	Number of Wilson-loop samples around each closed \(k\)-string. Must be at least 8 when Wilson loops are evaluated. Defaults to 32.	32
n_perp	int	Number of hybrid-Wannier samples from a TRIM plane's \(k_\perp=0\) edge to \(k_\perp=\pi\). Must be at least 5 when Wilson loops are evaluated. Defaults to 17.	17
parity_tolerance	float	Maximum relative \([H,I]\) commutator and inversion-eigenvalue deviation accepted at a TRIM. Defaults to 1e-5.	1e-05
kramers_tolerance	float	Maximum Wannier-center separation allowed when pairing Kramers partners on TRIM-plane endpoints. Defaults to 0.08.	0.08
gap_tolerance	float	Warning threshold for the minimum sampled occupied-to-empty direct gap. Defaults to 1e-8.	1e-08

Returns:

Type

Description

Z2ParityResult or Z2WilsonResult or Z2CombinedResult

Result mapping. Every method returns: "indices": \((\nu,)\) in 2-D or \((\nu_0, \nu_1, \nu_2, \nu_3)\) in 3-D, as integers in \(\{0,1\}\). "method": the construction that produced those indices. For method="parity" the mapping also contains Fu--Kane "parity_products" at each TRIM, per-TRIM "diagnostics" (including a labeled "parity_eigenvalues" Tensor of shape (n_occupied,)), and "direct_gap". For method="wilson" it contains hybrid-Wannier "planes" whose "wcc" (shape (n_perp, n_occupied)), "gap_pos", and "sweep" values are labeled Tensor objects, per-axis plane invariants "axis_z2", and "min_gap". For method="both" it contains both "parity" and "wilson" sub-results; "indices" follows the parity values.

Raises:

Type	Description
TypeError	If the Hamiltonian or inversion first dimension is not a MomentumSpace or MomentumBlockSpace, or either Hamiltonian matrix dimension is not a HilbertSpace.
ValueError	If method is unsupported; the input is not a rank-3 square Bloch Hamiltonian; n_occupied is invalid; the momentum space is not two- or three-dimensional and periodic with a diagonal cell; a MomentumBlockSpace contains off-diagonal \((k,k')\) blocks; Hamiltonian or inversion momenta do not form a unique complete reciprocal quotient; inversion band axes do not span the Hamiltonian Hilbert space; or n_loop / n_perp are below the Wilson-loop minima.
RuntimeError	For method="parity" or method="both", if inversion cannot be constructed or is not resolved at a TRIM.

Warns:

Type

Description

RuntimeWarning

If the sampled minimum direct gap is no larger than gap_tolerance; if method="auto" falls back from parity to Wilson loops; if the two 2-D Wilson orientations disagree; if the three 3-D Wilson axes disagree on \(\nu_0\) (majority vote is used); if Kramers pairing of Wannier centers is unresolved; if parity and Wilson indices disagree; or if the sampled spectrum is not time-reversal / Kramers consistent.

Notes

In 2-D, Fu--Kane gives \((-1)^\nu=\prod_i\delta(\Gamma_i)\). In 3-D, \((-1)^{\nu_0}=\prod_i\delta(\Gamma_i)\) and \((-1)^{\nu_j}=\prod_{k_j=\pi}\delta(\Gamma_i)\). Wilson indices use the hybrid-Wannier plane invariants described in the module docstring. Both constructions evaluate the Fourier interpolant of the input mesh rather than requiring TRIM or Wilson strings to sit on sampled \(k\)-points.

Examples:

Use Fu--Kane parities when an inversion tensor is available:

result = z2_indices(hamiltonian, n_occupied=2, inversion=inversion, method="parity")
indices = result["indices"]

Fall back to Wilson loops on a system without inversion:

wilson = z2_indices(hamiltonian, n_occupied=2, method="wilson")

See Also

chern_number First Chern number of a 2-D occupied bundle.

Source code in src/qten/topology/z2.py

def z2_indices(
    bloch_hamiltonian: Tensor,
    n_occupied: int | None = None,
    *,
    method: Literal["auto", "parity", "wilson", "both"] = "auto",
    inversion: Tensor | None = None,
    inversion_center: Offset | Sequence[float] | None = None,
    n_loop: int = 32,
    n_perp: int = 17,
    parity_tolerance: float = 1e-5,
    kramers_tolerance: float = 0.08,
    gap_tolerance: float = 1e-8,
) -> Z2ParityResult | Z2WilsonResult | Z2CombinedResult:
    r"""Compute the 2-D or 3-D \(\mathbb{Z}_2\) indices of an occupied band
    subspace.

    The ``n_occupied`` lowest-energy eigenstates, which must form an even
    number of Kramers pairs, define an occupied bundle over a complete
    two- or three-dimensional periodic momentum grid. Two numerical methods
    are available:

    - ``method="parity"`` evaluates Fu--Kane inversion eigenvalues at the
      \(2^d\) TRIM. Two dimensions return \((\nu,)\); three dimensions return
      \((\nu_0; \nu_1\nu_2\nu_3)\).
    - ``method="wilson"`` computes hybrid Wannier charge centers. In 2-D the
      two loop orientations should agree on \(\nu\). In 3-D the strong index
      is \(\nu_0=\nu(k_i=0)+\nu(k_i=\pi)\bmod 2\) and the weak indices are the
      three \(k_i=\pi\) plane invariants. If the three axes disagree on
      \(\nu_0\), the majority vote is returned.

    The Hamiltonian is Fourier-interpolated from the supplied mesh, so TRIM
    and Wilson strings need not coincide with sampled \(k\)-points. This is
    the construction used for odd meshes such as \(27^3\) or \(9^2\). The
    periodic cell must be diagonal in the primitive basis.

    Parameters
    ----------
    bloch_hamiltonian : Tensor
        Rank-3 Hermitian [`Tensor`][qten.linalg.tensors.Tensor] with dims
        ``(MomentumSpace, HilbertSpace, HilbertSpace)``, or with a
        [`MomentumBlockSpace`][qten.symbolics.state_space.MomentumBlockSpace]
        of diagonal \((k,k)\) blocks whose momenta form a complete 2-D or 3-D
        reciprocal quotient. The last two axes are square Bloch-Hamiltonian
        matrices and are aligned onto a common Hilbert space.
    n_occupied : int | None, optional
        Number of lowest-energy bands defining the occupied subspace. It must
        be even and lie strictly between zero and the total band count.
        Defaults to half the bands using integer division.
    method : {"auto", "parity", "wilson", "both"}, optional
        Numerical construction. ``"auto"`` tries Fu--Kane parity and falls
        back to Wilson loops if inversion cannot be resolved
        (`RuntimeError`). ``"both"`` requires parity to succeed. Defaults to
        ``"auto"``.
    inversion : Tensor | None, optional
        Optional rank-3 inversion operator with dims
        ``(MomentumSpace, HilbertSpace, HilbertSpace)``, or with a diagonal
        `MomentumBlockSpace` of \((k,k)\) blocks. Band axes are aligned onto
        the Hamiltonian Hilbert space. Each momentum is paired with
        \(H(k)\) by its label, not by data-axis order, and the labeled
        points must form the same complete reciprocal quotient. If omitted,
        spatial inversion is assembled from orbital
        [`Offset`][qten.geometries.spatials.Offset] labels about
        ``inversion_center``.
    inversion_center : Offset | Sequence[float] | None, optional
        Fixed point of spatial inversion, as an `Offset` or a \(d\)-vector in
        the Hamiltonian's direct-lattice coordinates. Defaults to the centroid
        of the unique orbital offsets.
    n_loop : int, optional
        Number of Wilson-loop samples around each closed \(k\)-string.
        Must be at least 8 when Wilson loops are evaluated. Defaults to 32.
    n_perp : int, optional
        Number of hybrid-Wannier samples from a TRIM plane's \(k_\perp=0\)
        edge to \(k_\perp=\pi\). Must be at least 5 when Wilson loops are
        evaluated. Defaults to 17.
    parity_tolerance : float, optional
        Maximum relative \([H,I]\) commutator and inversion-eigenvalue
        deviation accepted at a TRIM. Defaults to ``1e-5``.
    kramers_tolerance : float, optional
        Maximum Wannier-center separation allowed when pairing Kramers
        partners on TRIM-plane endpoints. Defaults to ``0.08``.
    gap_tolerance : float, optional
        Warning threshold for the minimum sampled occupied-to-empty direct
        gap. Defaults to ``1e-8``.

    Returns
    -------
    Z2ParityResult or Z2WilsonResult or Z2CombinedResult
        Result mapping. Every method returns:

        - ``"indices"``: \((\nu,)\) in 2-D or \((\nu_0, \nu_1, \nu_2, \nu_3)\)
          in 3-D, as integers in \(\{0,1\}\).
        - ``"method"``: the construction that produced those indices.

        For ``method="parity"`` the mapping also contains Fu--Kane
        ``"parity_products"`` at each TRIM, per-TRIM ``"diagnostics"``
        (including a labeled ``"parity_eigenvalues"``
        [`Tensor`][qten.linalg.tensors.Tensor] of shape
        ``(n_occupied,)``), and ``"direct_gap"``.

        For ``method="wilson"`` it contains hybrid-Wannier ``"planes"`` whose
        ``"wcc"`` (shape ``(n_perp, n_occupied)``), ``"gap_pos"``, and
        ``"sweep"`` values are labeled
        [`Tensor`][qten.linalg.tensors.Tensor] objects, per-axis plane
        invariants ``"axis_z2"``, and ``"min_gap"``.

        For ``method="both"`` it contains both ``"parity"`` and ``"wilson"``
        sub-results; ``"indices"`` follows the parity values.

    Raises
    ------
    TypeError
        If the Hamiltonian or inversion first dimension is not a
        `MomentumSpace` or `MomentumBlockSpace`, or either Hamiltonian
        matrix dimension is not a `HilbertSpace`.
    ValueError
        If ``method`` is unsupported; the input is not a rank-3 square Bloch
        Hamiltonian; ``n_occupied`` is invalid; the momentum space is not
        two- or three-dimensional and periodic with a diagonal cell; a
        `MomentumBlockSpace` contains off-diagonal \((k,k')\) blocks;
        Hamiltonian or inversion momenta do not form a unique complete
        reciprocal quotient; inversion band axes do not span the Hamiltonian
        Hilbert space; or ``n_loop`` / ``n_perp`` are below the Wilson-loop
        minima.
    RuntimeError
        For ``method="parity"`` or ``method="both"``, if inversion cannot be
        constructed or is not resolved at a TRIM.

    Warns
    -----
    RuntimeWarning
        If the sampled minimum direct gap is no larger than ``gap_tolerance``;
        if ``method="auto"`` falls back from parity to Wilson loops; if the
        two 2-D Wilson orientations disagree; if the three 3-D Wilson axes
        disagree on \(\nu_0\) (majority vote is used); if Kramers pairing of
        Wannier centers is unresolved; if parity and Wilson indices disagree;
        or if the sampled spectrum is not time-reversal / Kramers consistent.

    Notes
    -----
    In 2-D, Fu--Kane gives \((-1)^\nu=\prod_i\delta(\Gamma_i)\). In 3-D,
    \((-1)^{\nu_0}=\prod_i\delta(\Gamma_i)\) and
    \((-1)^{\nu_j}=\prod_{k_j=\pi}\delta(\Gamma_i)\). Wilson indices use the
    hybrid-Wannier plane invariants described in the module docstring. Both
    constructions evaluate the Fourier interpolant of the input mesh rather
    than requiring TRIM or Wilson strings to sit on sampled \(k\)-points.

    Examples
    --------
    Use Fu--Kane parities when an inversion tensor is available:

    ```python
    result = z2_indices(hamiltonian, n_occupied=2, inversion=inversion, method="parity")
    indices = result["indices"]
    ```

    Fall back to Wilson loops on a system without inversion:

    ```python
    wilson = z2_indices(hamiltonian, n_occupied=2, method="wilson")
    ```

    See Also
    --------
    [`chern_number`][qten.topology.chern_number]
        First Chern number of a 2-D occupied bundle.
    """
    method_name = str(method).lower()
    if method_name not in {"auto", "parity", "wilson", "both"}:
        raise ValueError("method must be 'auto', 'parity', 'wilson', or 'both'.")

    engine = _build_engine(bloch_hamiltonian, n_occupied, inversion, inversion_center)
    engine.check_time_reversal(stacklevel=2)
    parity_result: Z2ParityResult | None = None
    wilson_result: Z2WilsonResult | None = None

    if method_name in {"parity", "both"}:
        parity_result = engine.run_parity(parity_tolerance)
    elif method_name == "auto":
        try:
            parity_result = engine.run_parity(parity_tolerance)
        except RuntimeError as exc:
            warnings.warn(
                f"Parity method unavailable ({exc}). Falling back to Wilson loops.",
                RuntimeWarning,
                stacklevel=2,
            )

    if (
        method_name == "wilson"
        or method_name == "both"
        or (method_name == "auto" and parity_result is None)
    ):
        wilson_result = engine.run_wilson(int(n_loop), int(n_perp), kramers_tolerance)

    if method_name == "both":
        if parity_result is None or wilson_result is None:
            raise RuntimeError("method='both' requires parity and Wilson results.")
        if parity_result["indices"] != wilson_result["indices"]:
            warnings.warn(
                f"Parity {parity_result['indices']} and Wilson "
                f"{wilson_result['indices']} disagree.",
                RuntimeWarning,
                stacklevel=2,
            )
        chosen: Z2ParityResult | Z2WilsonResult | Z2CombinedResult = {
            "indices": parity_result["indices"],
            "method": "both",
            "parity": parity_result,
            "wilson": wilson_result,
        }
    elif parity_result is not None and wilson_result is None:
        chosen = parity_result
    elif wilson_result is not None:
        chosen = wilson_result
    else:
        raise RuntimeError("Z2 calculation produced no result.")

    if chosen["method"] == "parity":
        min_gap = chosen["direct_gap"]
    elif chosen["method"] == "wilson":
        min_gap = chosen["min_gap"]
    else:
        sampled = [chosen["parity"]["direct_gap"], chosen["wilson"]["min_gap"]]
        finite = [gap for gap in sampled if math.isfinite(gap)]
        min_gap = min(finite) if finite else float("nan")
    if math.isfinite(min_gap) and min_gap <= gap_tolerance:
        warnings.warn(
            f"Minimum sampled direct gap is {min_gap:.6e}; the occupied "
            "bundle is not isolated, so its Z2 indices are not well-defined.",
            RuntimeWarning,
            stacklevel=2,
        )
    return chosen