API documentation

Migration from the fermion-only API

Import PauliDecomposition from second_quantization.fermions instead of second_quantization.pauli_strings.
hilbert_space.basis_operators returns (fermion_dict, boson_dict). Pass this pair unchanged as operator_dict when reusing it in to_matrix.
Matrix conversion, symbolic_basis, and partial traces place fermions before bosons, preserving the input order within each type. The rightmost mode varies fastest.
Bosonic truncation t includes occupations 0 through t, giving t + 1 states. An integer applies to every boson; a list follows the bosons' input order. Fermion-only calls need no truncation.
Parity is 0 for even and 1 for odd total occupation for fermions, bosons, and mixtures. This preserves the fermion-only convention.

Symbolic parameters can be kept as dictionary keys for all particle types. Pauli decomposition also accepts matrices with symbolic entries. Bosonic inversion operates on numerical matrices and drops coefficients at or below 1e-10 by default; it reconstructs operators within the chosen finite cutoff.

Hilbert Space

`basis_operators(operators, sparse, truncation=None)`

Split fermionic and bosonic operators and build their matrix bases.

Fermionic operators are mapped through fermion_basis; bosonic operators are mapped through boson_basis with the provided truncation. The return value is a pair of dictionaries whose keys are 1 and the operator symbols/powers that are explicitly constructed.

Raises

ValueError If non-fermion/boson symbols are provided or bosonic truncation is omitted when bosons are present.

Source code in second_quantization/hilbert_space.py

def basis_operators(operators, sparse, truncation=None):
    """Split fermionic and bosonic operators and build their matrix bases.

    Fermionic operators are mapped through ``fermion_basis``; bosonic operators
    are mapped through ``boson_basis`` with the provided truncation. The return
    value is a pair of dictionaries whose keys are ``1`` and the operator
    symbols/powers that are explicitly constructed.

    Raises
    ------
    ValueError
        If non-fermion/boson symbols are provided or bosonic truncation is
        omitted when bosons are present.
    """
    fermion_ops, boson_ops = _split_operators(operators)

    fermion_dict = (
        fermions.fermion_basis(fermion_ops, sparse)
        if fermion_ops
        else {1: np.array([1])}
    )
    boson_dict = (
        bosons.boson_basis(boson_ops, truncation, sparse)
        if boson_ops
        else {1: np.array([1])}
    )
    return fermion_dict, boson_dict

`make_dict_callable(hamiltonian_dict)`

Create a callable function from a dictionary of SymPy expressions and NumPy arrays.

This function takes a dictionary where keys are symbolic expressions (containing free symbols) and values are NumPy arrays, and creates a callable function that evaluates the weighted sum of the arrays based on the symbolic expressions.

Parameters:

Name	Type	Description	Default
`hamiltonian_dict`	`dict[Expr, ndarray]`	A dictionary where keys are SymPy expressions (which may contain free symbols) and values are NumPy arrays of the same shape.	required

Returns:

Type	Description
`callable`	A callable function that takes values for all free symbols found in the
`callable`	dictionary keys and returns the weighted sum: Σ(expr_value * array) where
`callable`	expr_value is the numerical evaluation of each symbolic expression.

Raises:

Type	Description
`ValueError`	If any symbol names would be converted to Dummy variables by SymPy's lambdify function. This typically happens with special characters or reserved names.

Example

import sympy as sp
import numpy as np

# Define symbolic parameters
x, y = sp.symbols('x y')

# Create dictionary with symbolic expressions and matrices
ham_dict = {
    x: np.array([[1, 0], [0, 0]]),
    y: np.array([[0, 1], [1, 0]])
}

# Create callable function
func = make_dict_callable(ham_dict)

# Evaluate at specific parameter values
result = func(x=2.0, y=1.5)  # Returns 2.0 * first_matrix + 1.5 * second_matrix

Source code in second_quantization/hilbert_space.py

def make_dict_callable(
    hamiltonian_dict: dict[sympy.Expr, np.ndarray],
) -> callable:
    """Create a callable function from a dictionary of SymPy expressions and NumPy arrays.

    This function takes a dictionary where keys are symbolic expressions (containing
    free symbols) and values are NumPy arrays, and creates a callable function that
    evaluates the weighted sum of the arrays based on the symbolic expressions.

    Args:
        hamiltonian_dict: A dictionary where keys are SymPy expressions (which may
            contain free symbols) and values are NumPy arrays of the same shape.

    Returns:
        A callable function that takes values for all free symbols found in the
        dictionary keys and returns the weighted sum: Σ(expr_value * array) where
        expr_value is the numerical evaluation of each symbolic expression.

    Raises:
        ValueError: If any symbol names would be converted to Dummy variables by
            SymPy's lambdify function. This typically happens with special characters
            or reserved names.

    Example:
        ```python
        import sympy as sp
        import numpy as np

        # Define symbolic parameters
        x, y = sp.symbols('x y')

        # Create dictionary with symbolic expressions and matrices
        ham_dict = {
            x: np.array([[1, 0], [0, 0]]),
            y: np.array([[0, 1], [1, 0]])
        }

        # Create callable function
        func = make_dict_callable(ham_dict)

        # Evaluate at specific parameter values
        result = func(x=2.0, y=1.5)  # Returns 2.0 * first_matrix + 1.5 * second_matrix
        ```
    """
    all_symbols = list(
        sympy.ordered(sum([k for k in hamiltonian_dict.keys()]).free_symbols)
    )

    array_of_expr = sympy.Array(list(hamiltonian_dict.keys()))
    callable_expr = sympy.lambdify(all_symbols, array_of_expr, modules="numpy")
    tensor = np.array([v for v in list(hamiltonian_dict.values())])
    sum_rule = [None for _ in range(len(tensor.shape) - 1)]

    # check if the arguments of the callable contain dummys
    arg_names = list(inspect.signature(callable_expr).parameters.keys())
    if np.any(["Dummy" in name for name in arg_names]):
        raise ValueError(
            "Variable names must be chosen such that they are not lambdified to Dummy variables. Avoid special characters"
        )

    def func(*args, **kwargs):
        prefac = callable_expr(*args, **kwargs)[:, *sum_rule]
        return np.sum(prefac * tensor, axis=0)

    # adjust signature
    func.__signature__ = inspect.signature(callable_expr)

    return func

`parity_operator(operators, sparse, truncation=None)`

Construct total occupation parity: 0 for even and 1 for odd.

Fermion and boson occupations contribute modulo two. Tensor factors follow the order fermions then bosons, as in to_matrix.

Source code in second_quantization/hilbert_space.py

def parity_operator(
    operators: list[sympy.Expr], sparse: bool, truncation: int | list[int] | None = None
):
    """Construct total occupation parity: 0 for even and 1 for odd.

    Fermion and boson occupations contribute modulo two. Tensor factors follow
    the order fermions then bosons, as in ``to_matrix``.
    """

    fermion_ops, boson_ops = _split_operators(operators)

    f_parity = (
        fermions.fermion_parity(fermion_ops, sparse)
        if fermion_ops
        else np.zeros((1, 1), dtype=int)
    )
    b_parity = (
        bosons.boson_parity(boson_ops, truncation, sparse)
        if boson_ops
        else np.zeros((1, 1), dtype=int)
    )

    parity_diag = np.bitwise_xor.outer(f_parity.diagonal(), b_parity.diagonal()).ravel()
    if sparse:
        return scipy.sparse.diags(parity_diag, format="csr", dtype=parity_diag.dtype)
    else:
        return np.diag(parity_diag)

`partial_trace_generators(subset, all_operators, sparse, operator_dict=None, truncation=None)`

Projectors for tracing out fermionic and/or bosonic modes.

Constructs vectors that project onto each configuration of the complement of subset while fixing subset in vacuum, enabling partial traces via the recipe documented in the return value description.

Parameters:

Name	Type	Description	Default
`subset`	`list[Expr]`	Modes to remove via tracing.	required
`all_operators`	`list[Expr]`	All modes in the system, ordered as in `to_matrix`.	required
`sparse`	`bool`	Whether to keep intermediate matrices sparse.	required
`operator_dict`	`dict`	Optional cached operator dictionaries.	`None`
`truncation`	`int \| list[int] \| None`	Maximum bosonic occupations, or inferred from `operator_dict`.	`None`

Returns:

Type	Description
`ndarray`	Array of projectors `p` such that `sum(p[:, i, :] @ M @ p[:, i, :].conj().T)`
`ndarray`	traces out `subset` from matrix `M`. Remaining modes follow the
`ndarray`	fermions-first order of `symbolic_basis`.

Source code in second_quantization/hilbert_space.py

def partial_trace_generators(
    subset: list[sympy.Expr],
    all_operators: list[sympy.Expr],
    sparse: bool,
    operator_dict: dict = None,
    truncation: int | list[int] | None = None,
) -> np.ndarray:
    """Projectors for tracing out fermionic and/or bosonic modes.

    Constructs vectors that project onto each configuration of the complement of
    ``subset`` while fixing ``subset`` in vacuum, enabling partial traces via the
    recipe documented in the return value description.

    Args:
        subset: Modes to remove via tracing.
        all_operators: All modes in the system, ordered as in ``to_matrix``.
        sparse: Whether to keep intermediate matrices sparse.
        operator_dict: Optional cached operator dictionaries.
        truncation: Maximum bosonic occupations, or inferred from ``operator_dict``.

    Returns:
        Array of projectors ``p`` such that ``sum(p[:, i, :] @ M @ p[:, i, :].conj().T)``
        traces out ``subset`` from matrix ``M``. Remaining modes follow the
        fermions-first order of ``symbolic_basis``.
    """
    fermion_ops, boson_ops = _split_operators(all_operators)
    all_operators = fermion_ops + boson_ops
    set_subset = set([subset] if not isinstance(subset, list) else subset)
    complement = [element for element in all_operators if element not in set_subset]

    subset_operators, subset_powers = _get_all_combinations(
        subset=subset,
        all_operators=all_operators,
        sparse=sparse,
        operator_dict=operator_dict,
        truncation=truncation,
        return_powers=True,
    )
    complement_operators, complement_powers = _get_all_combinations(
        subset=complement,
        all_operators=all_operators,
        sparse=sparse,
        operator_dict=operator_dict,
        truncation=truncation,
        return_powers=True,
    )

    size = list(subset_operators.values())[0].shape[0]
    vacuum = np.zeros(size)
    vacuum[0] = 1

    vecs = []

    def _boson_norm(ops, powers):
        norm = 1.0
        for op, power in zip(ops, powers):
            if isinstance(op, BosonOp):
                norm *= math.sqrt(math.factorial(int(power)))
        return norm

    subset_items = list(subset_operators.values())
    subset_scaled = [
        mat / _boson_norm(subset, powers)
        for mat, powers in zip(subset_items, subset_powers)
    ]

    complement_items = list(complement_operators.values())
    complement_scaled = [
        mat / _boson_norm(complement, powers)
        for mat, powers in zip(complement_items, complement_powers)
    ]

    for element in complement_scaled:
        aux = []
        for g in subset_scaled:
            generator = g @ element
            aux.append(vacuum @ generator)
        vecs.append(aux)

    return np.array(vecs)

`symbolic_basis(operators, truncation=None)`

Generate creation monomials in matrix basis order.

Fermions precede bosons, preserving the order within each particle type, as in to_matrix. The rightmost mode varies fastest.

Parameters:

Name	Type	Description	Default
`operators`	`list[FermionOp \| BosonOp]`	Fermionic and/or bosonic modes.	required
`truncation`	`int \| list[int] \| None`	Maximum bosonic occupation, either shared or per boson. Required only when bosons are present.	`None`

Returns:

Type	Description
`list[Expr]`	Creation monomials acting on vacuum, starting with `sympy.S.One`.
`list[Expr]`	Bosonic powers are unnormalized: occupation `n` contributes a factor
`list[Expr]`	of `sqrt(n!)` when applied to vacuum.

Example

For 2 fermions [c, d], returns: - [1, d†, c†, c†*d†] representing vacuum, single occupations, and double occupation.

Source code in second_quantization/hilbert_space.py

def symbolic_basis(
    operators: list[FermionOp | BosonOp], truncation: int | list[int] | None = None
) -> list[sympy.Expr]:
    """Generate creation monomials in matrix basis order.

    Fermions precede bosons, preserving the order within each particle type,
    as in ``to_matrix``. The rightmost mode varies fastest.

    Args:
        operators: Fermionic and/or bosonic modes.
        truncation: Maximum bosonic occupation, either shared or per boson.
            Required only when bosons are present.

    Returns:
        Creation monomials acting on vacuum, starting with ``sympy.S.One``.
        Bosonic powers are unnormalized: occupation ``n`` contributes a factor
        of ``sqrt(n!)`` when applied to vacuum.

    Example:
        For 2 fermions `[c, d]`, returns:
        - `[1, d†, c†, c†*d†]` representing vacuum, single occupations, and double occupation.
    """
    fermion_ops, boson_ops = _split_operators(operators)
    operators = fermion_ops + boson_ops
    truncation = bosons._truncation_list(truncation, len(boson_ops))
    symbols = []
    strings = _generate_all_powers(operators, truncation)

    for string in strings:
        sub_symbol = sympy.S.One
        for idx, power in enumerate(string):
            if power == 1:
                sub_symbol *= Dagger(operators[idx])
            elif power > 1:
                sub_symbol *= Dagger(operators[idx]) ** power
        symbols.append(sub_symbol)
    return symbols

`to_matrix(expression, operators, sparse, operator_dict=None, truncation=None)`

Convert a symbolic operator expression to matrices.

Expands expression into normal products of the provided fermionic and/or bosonic operators, builds matrix representations via tensor products, and groups terms by their purely symbolic prefactor. Creation operators are inferred via Dagger and implemented as conjugate-transposed annihilation matrices.

Parameters:

Name	Type	Description	Default
`expression`	`Expr`	SymPy expression containing the operators to expand.	required
`operators`	`list[FermionOp \| BosonOp]`	Modes whose stable fermion partition, followed by their stable boson partition, sets the tensor-product ordering.	required
`sparse`	`bool`	If True, use sparse matrices for all operations.	required
`operator_dict`	`dict`	Optional cached `(fermion_dict, boson_dict)` from `basis_operators` to avoid recomputation.	`None`
`truncation`	`int \| list[int]`	Truncation passed to bosonic basis construction when needed.	`None`

Returns:

Type	Description
`dict[Expr, ndarray \| csr_array]`	Dict mapping symbolic prefactors to the summed matrix representation of
`dict[Expr, ndarray \| csr_array]`	all terms that share that prefactor.

Source code in second_quantization/hilbert_space.py

def to_matrix(
    expression: sympy.Expr,
    operators: list[
        sympy.physics.quantum.fermion.FermionOp | sympy.physics.quantum.boson.BosonOp
    ],
    sparse: bool,
    operator_dict: dict = None,
    truncation: int | list[int] = None,
) -> dict[sympy.Expr, np.ndarray | scipy.sparse.csr_array]:
    """Convert a symbolic operator expression to matrices.

    Expands ``expression`` into normal products of the provided fermionic and/or
    bosonic operators, builds matrix representations via tensor products, and
    groups terms by their purely symbolic prefactor. Creation operators are
    inferred via ``Dagger`` and implemented as conjugate-transposed annihilation
    matrices.

    Args:
        expression: SymPy expression containing the operators to expand.
        operators: Modes whose stable fermion partition, followed by their stable
            boson partition, sets the tensor-product ordering.
        sparse: If True, use sparse matrices for all operations.
        operator_dict: Optional cached ``(fermion_dict, boson_dict)`` from
            ``basis_operators`` to avoid recomputation.
        truncation: Truncation passed to bosonic basis construction when needed.

    Returns:
        Dict mapping symbolic prefactors to the summed matrix representation of
        all terms that share that prefactor.
    """

    if operator_dict is None:
        fermion_dict, boson_dict = basis_operators(operators, sparse, truncation)
    else:
        fermion_dict, boson_dict = operator_dict

    if sparse:

        def kron(a, b):
            return scipy.sparse.kron(a, b, format="csr")
    else:
        kron = np.kron

    dict_matrices = {}
    # Expand the expression to ensure all terms are separated
    expression = expression.expand()
    # Iterate over each term in the expression
    for term, coeff in expression.as_coefficients_dict().items():
        term_mat = float(coeff)
        symbol = S.One

        fermion_factors = fermion_dict[1].copy()
        boson_factors = boson_dict[1].copy()

        for factor in term.as_ordered_factors():
            base, power = factor.as_base_exp()
            if factor in fermion_dict:
                fermion_factors = fermion_factors @ fermion_dict[factor]
            elif Dagger(factor) in fermion_dict:
                fermion_factors = (
                    fermion_factors @ fermion_dict[Dagger(factor)].conj().T
                )
            elif factor in boson_dict:
                boson_factors = boson_factors @ boson_dict[factor]
            elif Dagger(factor) in boson_dict:
                boson_factors = boson_factors @ boson_dict[Dagger(factor)].conj().T
            elif (
                isinstance(base, BosonOp)
                and (base in operators or Dagger(base) in operators)
                and power.is_Integer
                and power > 0
            ):
                # Powers beyond the truncated bosonic basis vanish.
                boson_factors = 0 * boson_factors
            else:
                symbol *= factor

        term_mat = term_mat * kron(fermion_factors, boson_factors)

        if symbol in dict_matrices.keys():
            dict_matrices[symbol] += term_mat
        else:
            dict_matrices[symbol] = term_mat
    return dict_matrices

`to_operators(matrix, basis, truncation=None)`

Recover a normal-ordered operator expression from a matrix.

Dispatches to the appropriate inversion backend based on the operator types present in basis:

Fermionic only – Pauli decomposition + Jordan–Wigner inversion (original behaviour, requires dim = 2^n).
Bosonic only – diagonal falling-factorial decomposition with finite-difference inversion (requires truncation).
Mixed – Pauli decomposition on the fermionic tensor structure followed by bosonic inversion on each coefficient block (requires truncation).

If a dictionary of matrices is supplied, each entry is converted and scaled by its symbolic key before summing.

Parameters

matrix : ndarray, csr_array, or dict thereof Square matrix (or dict of matrices) in the Hilbert space defined by basis and the tensor-product ordering of :func:to_matrix. basis : list of FermionOp / BosonOp Operators defining the Hilbert space. They are stably partitioned into fermions followed by bosons, matching :func:to_matrix. truncation : int, list of int, or None Required when basis contains bosonic operators; passed to the bosonic inversion backend.

Returns

sympy.Expr Normal-ordered symbolic expression.

Source code in second_quantization/hilbert_space.py

def to_operators(
    matrix: dict[sympy.Expr, np.ndarray | scipy.sparse.csr_array]
    | np.ndarray
    | scipy.sparse.csr_array,
    basis: list[sympy.Expr],
    truncation: int | list[int] | None = None,
) -> sympy.Expr:
    """Recover a normal-ordered operator expression from a matrix.

    Dispatches to the appropriate inversion backend based on the operator
    types present in ``basis``:

    * **Fermionic only** – Pauli decomposition + Jordan–Wigner inversion
      (original behaviour, requires dim = 2^n).
    * **Bosonic only** – diagonal falling-factorial decomposition with
      finite-difference inversion (requires ``truncation``).
    * **Mixed** – Pauli decomposition on the fermionic tensor structure
      followed by bosonic inversion on each coefficient block (requires
      ``truncation``).

    If a dictionary of matrices is supplied, each entry is converted and
    scaled by its symbolic key before summing.

    Parameters
    ----------
    matrix : ndarray, csr_array, or dict thereof
        Square matrix (or dict of matrices) in the Hilbert space defined by
        ``basis`` and the tensor-product ordering of :func:`to_matrix`.
    basis : list of FermionOp / BosonOp
        Operators defining the Hilbert space. They are stably partitioned into
        fermions followed by bosons, matching :func:`to_matrix`.
    truncation : int, list of int, or None
        Required when ``basis`` contains bosonic operators; passed to the
        bosonic inversion backend.

    Returns
    -------
    sympy.Expr
        Normal-ordered symbolic expression.
    """
    fermion_ops, boson_ops = _split_operators(basis)
    truncation = bosons._truncation_list(truncation, len(boson_ops))
    dimension = 2 ** len(fermion_ops) * math.prod(t + 1 for t in truncation)
    matrices = matrix.items() if isinstance(matrix, dict) else ((S.One, matrix),)
    fermion_cache = {"": S.One}
    terms = {}

    for prefactor, value in matrices:
        if value.shape != (dimension, dimension):
            raise ValueError("Matrix dimension is incompatible with basis.")
        for pauli_string, block in fermions._pauli_coefficients(
            value, len(fermion_ops)
        ).items():
            if pauli_string not in fermion_cache:
                fermion_cache[pauli_string] = fermions.string_to_fermion_operators(
                    pauli_string, fermion_ops
                )
            coefficients = bosons._boson_coefficients(
                block, truncation, 1e-10 if boson_ops else 0.0
            )
            for powers, coefficient in coefficients.items():
                coefficient = bosons._symbolic_coefficient(coefficient)
                terms[powers] = terms.get(powers, S.Zero) + (
                    prefactor * coefficient * fermion_cache[pauli_string]
                )

    result = S.Zero
    for powers, fermion_expr in terms.items():
        if fermion_ops:
            fermion_expr = normal_ordered_form(
                sympy.expand(fermion_expr), independent=True
            )
        result += fermion_expr * bosons._boson_monomial(powers, boson_ops)
    return sympy.expand(result)

Fermions

`PauliDecomposition(matrix, threshold=0.0)`

Decompose a matrix into Pauli-string coefficients.

Parameters

matrix : np.ndarray | scipy.sparse.spmatrix | sympy.Matrix Square matrix with power-of-two dimension. threshold : float Omit coefficients whose absolute value is not larger than this value. Retain symbolic coefficients whose magnitude cannot be determined.

Returns

dict[str, complex | float | sympy.Basic] Pauli strings mapped to scalar coefficients.

Raises

ValueError If matrix is not square or its dimension is not a power of two.

Source code in second_quantization/fermions.py

def PauliDecomposition(
    matrix: np.ndarray | scipy.sparse.spmatrix | sympy.Matrix,
    threshold: float = 0.0,
) -> dict[str, complex | float | sympy.Basic]:
    """Decompose a matrix into Pauli-string coefficients.

    Parameters
    ----------
    matrix : np.ndarray | scipy.sparse.spmatrix | sympy.Matrix
        Square matrix with power-of-two dimension.
    threshold : float
        Omit coefficients whose absolute value is not larger than this value.
        Retain symbolic coefficients whose magnitude cannot be determined.

    Returns
    -------
    dict[str, complex | float | sympy.Basic]
        Pauli strings mapped to scalar coefficients.

    Raises
    ------
    ValueError
        If ``matrix`` is not square or its dimension is not a power of two.
    """
    if matrix.shape[0] != matrix.shape[1]:
        raise ValueError("Matrix is not square.")
    size = matrix.shape[0]
    if size == 0 or size & (size - 1):
        raise ValueError("Matrix dimension is not a power of 2.")
    blocks = _pauli_coefficients(matrix, size.bit_length() - 1, threshold)
    return {key: block[0, 0] for key, block in blocks.items()}

`fermion_basis(operators, sparse)`

Build Jordan–Wigner annihilation matrices.

Parameters

operators : list[sympy.Expr] Fermionic modes in tensor-product order. sparse : bool Return CSR arrays instead of dense arrays.

Returns

dict[sympy.Expr, np.ndarray | scipy.sparse.csr_array] Full-space identity and one annihilation matrix per mode.

Source code in second_quantization/fermions.py

def fermion_basis(
    operators: list[sympy.Expr], sparse: bool
) -> dict[sympy.Expr, np.ndarray | scipy.sparse.csr_array :]:
    """Build Jordan–Wigner annihilation matrices.

    Parameters
    ----------
    operators : list[sympy.Expr]
        Fermionic modes in tensor-product order.
    sparse : bool
        Return CSR arrays instead of dense arrays.

    Returns
    -------
    dict[sympy.Expr, np.ndarray | scipy.sparse.csr_array]
        Full-space identity and one annihilation matrix per mode.
    """

    if sparse:

        def kron(a, b):
            return scipy.sparse.kron(a, b, format="csr")

        def eye(N):
            return scipy.sparse.eye(N, format="csr")
    else:
        kron = np.kron
        eye = np.eye

    size = len(operators)
    basis_operators = {1: eye(2**size)}
    for i, operator in enumerate(operators):
        basis_operators[operator] = reduce(
            kron, [sigma_z] * i + [sigma_minus] + [eye(2 ** (size - i - 1))]
        )

    return basis_operators

`fermion_parity(operators, sparse)`

Build the fermionic occupation-parity matrix.

Parameters

operators : list[sympy.Expr] Fermionic modes defining the register. sparse : bool Return a CSR array instead of a dense array.

Returns

np.ndarray | scipy.sparse.csr_array Diagonal matrix with 0 for even and 1 for odd occupation.

Source code in second_quantization/fermions.py

def fermion_parity(
    operators: list[sympy.Expr], sparse: bool
) -> np.ndarray | scipy.sparse.csr_array:
    """Build the fermionic occupation-parity matrix.

    Parameters
    ----------
    operators : list[sympy.Expr]
        Fermionic modes defining the register.
    sparse : bool
        Return a CSR array instead of a dense array.

    Returns
    -------
    np.ndarray | scipy.sparse.csr_array
        Diagonal matrix with 0 for even and 1 for odd occupation.
    """
    size = 2 ** len(operators)
    parity_diag = np.bitwise_count(np.arange(size)) % 2
    if sparse:
        return scipy.sparse.diags(parity_diag, format="csr", dtype=parity_diag.dtype)
    else:
        return np.diag(parity_diag)

`string_to_fermion_operators(pauli_str, fermion_ops)`

Convert a Pauli string through the inverse Jordan–Wigner map.

Parameters

pauli_str : str String containing I, X, Y, and Z labels. fermion_ops : list[sympy.Expr] Fermionic modes corresponding to the string positions.

Returns

sympy.Expr Fermionic operator expression before normal ordering.

Source code in second_quantization/fermions.py

def string_to_fermion_operators(pauli_str: str, fermion_ops: list) -> sympy.Expr:
    """Convert a Pauli string through the inverse Jordan–Wigner map.

    Parameters
    ----------
    pauli_str : str
        String containing ``I``, ``X``, ``Y``, and ``Z`` labels.
    fermion_ops : list[sympy.Expr]
        Fermionic modes corresponding to the string positions.

    Returns
    -------
    sympy.Expr
        Fermionic operator expression before normal ordering.
    """
    expr = S.One
    parity = S.One

    for key, op in zip(pauli_str, fermion_ops):
        if key == "X":
            expr = expr * parity * (Dagger(op) + op)
        elif key == "Y":
            expr = expr * parity * (1j * (Dagger(op) - op))
        elif key == "Z":
            expr = expr * (-(2 * Dagger(op) * op - 1))
        parity = parity * (-(2 * Dagger(op) * op - 1))

    return expr

Bosons

`BosonDecomposition(matrix, truncation, threshold=1e-10)`

Decompose a matrix into normal-ordered bosonic coefficients.

Parameters

matrix : np.ndarray | scipy.sparse.csr_array Square matrix with dimension prod(t + 1 for t in truncation). truncation : list[int] | int Per-mode cutoffs; a scalar describes one mode. threshold : float Omit coefficients whose absolute value is not larger than this value.

Returns

dict[tuple[tuple[int, int], ...], complex | float] Per-mode (creation, annihilation) powers mapped to coefficients.

Raises

ValueError If the matrix shape and truncation are incompatible.

Source code in second_quantization/bosons.py

def BosonDecomposition(
    matrix: np.ndarray | scipy.sparse.csr_array,
    truncation: list[int] | int,
    threshold: float = 1e-10,
) -> dict[tuple[tuple[int, int], ...], complex | float]:
    """Decompose a matrix into normal-ordered bosonic coefficients.

    Parameters
    ----------
    matrix : np.ndarray | scipy.sparse.csr_array
        Square matrix with dimension ``prod(t + 1 for t in truncation)``.
    truncation : list[int] | int
        Per-mode cutoffs; a scalar describes one mode.
    threshold : float
        Omit coefficients whose absolute value is not larger than this value.

    Returns
    -------
    dict[tuple[tuple[int, int], ...], complex | float]
        Per-mode ``(creation, annihilation)`` powers mapped to coefficients.

    Raises
    ------
    ValueError
        If the matrix shape and truncation are incompatible.
    """

    truncation = [truncation] if isinstance(truncation, int) else truncation
    if matrix.shape[0] != matrix.shape[1]:
        raise ValueError("Matrix is not square.")
    if matrix.shape[0] != prod(t + 1 for t in truncation):
        raise ValueError("Matrix dimension is incompatible with truncation.")
    return _boson_coefficients(matrix, truncation, threshold)

`boson_basis(operators, truncation, sparse)`

Build truncated bosonic annihilation-power matrices.

Parameters

operators : list[sympy.Expr] Bosonic modes in tensor-product order. truncation : list[int] | int Per-mode cutoffs; a scalar applies to every mode. sparse : bool Return CSR arrays instead of dense arrays.

Returns

dict[sympy.Expr, np.ndarray | scipy.sparse.csr_array] Full-space matrices for powers zero through each mode's cutoff.

Raises

ValueError If the number of cutoffs does not match the number of modes.

Source code in second_quantization/bosons.py

def boson_basis(
    operators: list[sympy.Expr], truncation: list[int] | int, sparse: bool
) -> dict[sympy.Expr, np.ndarray | scipy.sparse.csr_array]:
    """Build truncated bosonic annihilation-power matrices.

    Parameters
    ----------
    operators : list[sympy.Expr]
        Bosonic modes in tensor-product order.
    truncation : list[int] | int
        Per-mode cutoffs; a scalar applies to every mode.
    sparse : bool
        Return CSR arrays instead of dense arrays.

    Returns
    -------
    dict[sympy.Expr, np.ndarray | scipy.sparse.csr_array]
        Full-space matrices for powers zero through each mode's cutoff.

    Raises
    ------
    ValueError
        If the number of cutoffs does not match the number of modes.
    """
    truncation = _truncation_list(truncation, len(operators))

    if sparse:
        kron = partial(scipy.sparse.kron, format="csr")
        eye = partial(scipy.sparse.eye, format="csr")
    else:
        kron = np.kron
        eye = np.eye

    def annihilation_powers(dim: int, max_power: int):
        identity = eye(dim)
        powers = [identity]
        if max_power == 0:
            return powers

        factors = np.ones(dim, dtype=float)
        for p in range(1, max_power + 1):
            factors[p:] *= np.sqrt(np.arange(1, dim - p + 1, dtype=float))
            data = factors[p:].copy()
            if sparse:
                mat = scipy.sparse.diags([data], [p], shape=(dim, dim), format="csr")
            else:
                mat = np.diag(data, k=p)
            powers.append(mat)
        return powers

    boson_basis = {}

    identity_matrices = [eye(truncation[n] + 1) for n in range(len(operators))]

    for i, operator in enumerate(operators):
        base_operators = identity_matrices.copy()

        powers = annihilation_powers(truncation[i] + 1, truncation[i])
        for power in range(truncation[i] + 1):
            base_operators[i] = powers[power]
            boson_basis[operator**power] = reduce(kron, base_operators)
            base_operators[i] = identity_matrices[i]

    return boson_basis

`boson_parity(operators, truncation, sparse)`

Build the bosonic occupation-parity matrix.

Parameters

operators : list[sympy.Expr] Bosonic modes in tensor-product order. truncation : list[int] | int Per-mode cutoffs; a scalar applies to every mode. sparse : bool Return a CSR array instead of a dense array.

Returns

np.ndarray | scipy.sparse.csr_array Diagonal matrix with 0 for even and 1 for odd total occupation.

Raises

ValueError If the number of cutoffs does not match the number of modes.

Source code in second_quantization/bosons.py

def boson_parity(
    operators: list[sympy.Expr], truncation: list[int] | int, sparse: bool
) -> np.ndarray | scipy.sparse.csr_array:
    """Build the bosonic occupation-parity matrix.

    Parameters
    ----------
    operators : list[sympy.Expr]
        Bosonic modes in tensor-product order.
    truncation : list[int] | int
        Per-mode cutoffs; a scalar applies to every mode.
    sparse : bool
        Return a CSR array instead of a dense array.

    Returns
    -------
    np.ndarray | scipy.sparse.csr_array
        Diagonal matrix with 0 for even and 1 for odd total occupation.

    Raises
    ------
    ValueError
        If the number of cutoffs does not match the number of modes.
    """
    truncation = _truncation_list(truncation, len(operators))

    dims = np.array([t + 1 for t in truncation])
    total_dim = np.prod(dims)

    indices = np.arange(total_dim)

    total_occupation = np.zeros(total_dim, dtype=int)
    temp_indices = indices.copy()

    # The rightmost Kronecker factor is the fastest-varying index.
    for dim in dims[::-1]:
        occupation = temp_indices % dim
        total_occupation += occupation
        temp_indices //= dim
    parity_diag = total_occupation % 2

    if sparse:
        return scipy.sparse.diags(parity_diag, format="csr", dtype=parity_diag.dtype)
    else:
        return np.diag(parity_diag)

`boson_to_operators(matrix, operators, truncation, threshold=1e-10)`

Recover a normal-ordered bosonic expression from a matrix.

Parameters

matrix : np.ndarray | scipy.sparse.csr_array Square matrix in the truncated bosonic basis. operators : list[sympy.Expr] Bosonic modes in tensor-product order. truncation : list[int] | int Per-mode cutoffs; a scalar applies to every mode. threshold : float Omit coefficients whose absolute value is not larger than this value.

Returns

sympy.Expr Normal-ordered expression in operators.

Source code in second_quantization/bosons.py

def boson_to_operators(
    matrix: np.ndarray | scipy.sparse.csr_array,
    operators: list[sympy.Expr],
    truncation: list[int] | int,
    threshold: float = 1e-10,
) -> sympy.Expr:
    """Recover a normal-ordered bosonic expression from a matrix.

    Parameters
    ----------
    matrix : np.ndarray | scipy.sparse.csr_array
        Square matrix in the truncated bosonic basis.
    operators : list[sympy.Expr]
        Bosonic modes in tensor-product order.
    truncation : list[int] | int
        Per-mode cutoffs; a scalar applies to every mode.
    threshold : float
        Omit coefficients whose absolute value is not larger than this value.

    Returns
    -------
    sympy.Expr
        Normal-ordered expression in ``operators``.
    """

    truncation = _truncation_list(truncation, len(operators))
    return sympy.expand(
        sum(
            _symbolic_coefficient(coefficient) * _boson_monomial(powers, operators)
            for powers, coefficient in BosonDecomposition(
                matrix, truncation, threshold
            ).items()
        )
    )