API documentation

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 = [], []
    for operator in operators:
        if isinstance(operator, sympy.physics.quantum.fermion.FermionOp):
            fermion_ops.append(operator)
        elif isinstance(operator, sympy.physics.quantum.boson.BosonOp):
            boson_ops.append(operator)
        else:
            raise ValueError("Can't submit non operator types")
    if not truncation and len(boson_ops) > 0:
        raise ValueError("Please provide a truncation for your bosons")

    fermion_dict = (
        fermions.fermion_basis(fermion_ops, sparse)
        if len(fermion_ops) != 0
        else {1: np.array([1])}
    )
    boson_dict = (
        bosons.boson_basis(boson_ops, truncation, sparse)
        if len(boson_ops) != 0
        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 the combined parity operator for fermionic/bosonic mixtures.

Fermionic modes contribute their occupation-parity bit; bosonic modes contribute (-1) raised to their occupation number (requiring a truncation). The individual parities are tensored in the order fermions then bosons.

Source code in second_quantization/hilbert_space.py

def parity_operator(
    operators: list[sympy.Expr], sparse: bool, truncation: int | list[int] | None = None
):
    """Construct the combined parity operator for fermionic/bosonic mixtures.

    Fermionic modes contribute their occupation-parity bit; bosonic modes
    contribute ``(-1)`` raised to their occupation number (requiring a
    truncation). The individual parities are tensored in the order fermions
    then bosons.
    """

    fermion_ops = [
        op
        for op in operators
        if isinstance(op, sympy.physics.quantum.fermion.FermionOp)
    ]
    boson_ops = [
        op for op in operators if isinstance(op, sympy.physics.quantum.boson.BosonOp)
    ]

    if boson_ops and truncation is None:
        raise ValueError("Please provide a truncation for bosonic parity operators")

    f_parity = (
        fermions.fermion_parity(fermion_ops, sparse)
        if fermion_ops
        else (scipy.sparse.eye(1, format="csr") if sparse else np.eye(1))
    )
    b_parity = (
        bosons.boson_parity(boson_ops, truncation, sparse)
        if boson_ops
        else (scipy.sparse.eye(1, format="csr") if sparse else np.eye(1))
    )

    if sparse:
        result = scipy.sparse.kron(f_parity, b_parity, format="csr")
    else:
        result = np.kron(f_parity, b_parity)
    return result

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

Projectors for tracing out a fermionic subset.

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]`	Fermionic modes to remove via tracing.	required
`all_operators`	`list[Expr]`	Ordered list of all fermionic modes in the system.	required
`sparse`	`bool`	Whether to keep intermediate matrices sparse.	required
`operator_dict`	`dict`	Optional cached operator dictionaries.	`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`.

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 a fermionic subset.

    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: Fermionic modes to remove via tracing.
        all_operators: Ordered list of all fermionic modes in the system.
        sparse: Whether to keep intermediate matrices sparse.
        operator_dict: Optional cached operator dictionaries.

    Returns:
        Array of projectors ``p`` such that ``sum(p[:, i, :] @ M @ p[:, i, :].conj().T)``
        traces out ``subset`` from matrix ``M``.
    """
    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)`

Generate a list of symbolic operators corresponding to their occupation.

This function creates the symbolic representation of all possible occupation states in the fermionic Fock space for the given set of fermionic operators.

Parameters:

Name	Type	Description	Default
`fermions`		List of fermionic operators.	required

Returns:

Type	Description
`list[Expr]`	List of symbolic fermionic operators representing all possible occupation states.
`list[Expr]`	The first element is always the vacuum state (represented as `sympy.S.One`),
`list[Expr]`	followed by all single-particle states, two-particle states, etc.

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]
) -> list[sympy.Expr]:
    """Generate a list of symbolic operators corresponding to their occupation.

    This function creates the symbolic representation of all possible occupation states
    in the fermionic Fock space for the given set of fermionic operators.

    Args:
        fermions: List of fermionic operators.

    Returns:
        List of symbolic fermionic operators representing all possible occupation states.
        The first element is always the vacuum state (represented as `sympy.S.One`),
        followed by all single-particle states, two-particle states, etc.

    Example:
        For 2 fermions `[c, d]`, returns:
        - `[1, d†, c†, c†*d†]` representing vacuum, single occupations, and double occupation.
    """
    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]`	Ordered basis of modes (fermions, bosons, or both) that 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: Ordered basis of modes (fermions, bosons, or both) that 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():
            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
            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 Ordered operators defining the Hilbert space. Fermionic modes must come before bosonic modes (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
        Ordered operators defining the Hilbert space.  Fermionic modes must
        come before bosonic modes (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.
    """
    if isinstance(matrix, dict):
        return sum([k * to_operators(v, basis, truncation) for k, v in matrix.items()])

    fermion_ops = [op for op in basis if isinstance(op, FermionOp)]
    boson_ops = [op for op in basis if isinstance(op, BosonOp)]

    if boson_ops and truncation is None:
        raise ValueError("truncation is required when bosonic operators are present")

    trunc_list = (
        []
        if not boson_ops
        else (
            [truncation] * len(boson_ops)
            if isinstance(truncation, int)
            else list(truncation)
        )
    )

    if fermion_ops:
        if boson_ops:
            blocks = _fermion_blocks_from_mixed(matrix, len(fermion_ops))
            pauli_blocks = dict(
                PauliDecomposition(blocks, recursion_depth=len(fermion_ops))[0]
            )
        else:
            pauli_blocks = dict(PauliDecomposition(matrix)[0])
    else:
        pauli_blocks = {"": matrix}

    term = S.Zero
    for pauli_str, coeff in pauli_blocks.items():
        fermion_expr = (
            normal_ordered_form(
                string_to_fermion_operators(pauli_str, fermion_ops), independent=True
            )
            if fermion_ops
            else S.One
        )
        if boson_ops:
            boson_expr = bosons.boson_to_operators(coeff, boson_ops, trunc_list)
        else:
            boson_expr = sympy.nsimplify(coeff, rational=False)
        term += fermion_expr * boson_expr

    term = sympy.nsimplify(term).expand()

    # SymPy's normal_ordered_form does not preserve mixed products such as
    # c * Dagger(a): it applies the same commutation rule across species. The
    # individual fermionic and bosonic factors above are already normal ordered.
    if fermion_ops and boson_ops:
        return term
    return normal_ordered_form(term, independent=True)

Pauli strings

`PauliDecomposition(matrix, PauliStringInit='', recursion_depth=None, threshold=0.0)`

Decompose a square matrix into Pauli strings.

Recursively slices the matrix into 2x2 blocks, decomposes each block, and accumulates tensor-product labels composed of I, X, Y, and Z. Works with dense numpy arrays, scipy.sparse matrices, and sympy.Matrix so long as the dimension is a power of two.

Parameters:

Name	Type	Description	Default
`matrix`		Square matrix to decompose.	required
`PauliStringInit`		Prefix used internally for recursion.	`''`

Returns: defaultdict mapping Pauli strings (e.g., "IXZ") to their complex or symbolic coefficients.

Source code in second_quantization/pauli_strings.py

def PauliDecomposition(matrix, PauliStringInit="", recursion_depth=None, threshold=0.0):
    """Decompose a square matrix into Pauli strings.

    Recursively slices the matrix into 2x2 blocks, decomposes each block, and
    accumulates tensor-product labels composed of ``I``, ``X``, ``Y``, and
    ``Z``. Works with dense numpy arrays, ``scipy.sparse`` matrices, and
    ``sympy.Matrix`` so long as the dimension is a power of two.

    Args:
        matrix: Square matrix to decompose.
        PauliStringInit: Prefix used internally for recursion.
    Returns:
        ``defaultdict`` mapping Pauli strings (e.g., ``"IXZ"``) to their
        complex or symbolic coefficients.
    """
    # Dimension check
    if matrix.shape[0] != matrix.shape[1]:
        raise ValueError("Matrix is not square.")
    if (qBitDim := np.log2(matrix.shape[0])) % 1 != 0:
        raise ValueError("Matrix dimension is not a power of 2.")
    qBitDim = int(qBitDim)
    if recursion_depth is None:
        recursion_depth = qBitDim
    if recursion_depth < 0:
        raise ValueError("recursion_depth must be non-negative.")
    I = sympy.I if isinstance(matrix, sympy.Matrix) else 1j  # noqa: E741

    decomposition = defaultdict(lambda: 0)

    # Output for dimension 1
    if qBitDim == 0:
        if _max_abs(matrix[0, 0]) > threshold:
            decomposition[PauliStringInit] = matrix[0, 0]
        remainder = matrix * 0
        return decomposition, remainder

    if recursion_depth == 0:
        return decomposition, matrix

    # Calculates the tensor product coefficients via the sliced submatrices.
    halfDim = int(2 ** (qBitDim - 1))

    coeff1 = (matrix[0:halfDim, 0:halfDim] + matrix[halfDim:, halfDim:]) / 2
    coeffX = (matrix[halfDim:, 0:halfDim] + matrix[0:halfDim, halfDim:]) / 2
    coeffY = -I * (matrix[halfDim:, 0:halfDim] - matrix[0:halfDim, halfDim:]) / 2
    coeffZ = (matrix[0:halfDim, 0:halfDim] - matrix[halfDim:, halfDim:]) / 2

    coefficients = {"I": coeff1, "X": coeffX, "Y": coeffY, "Z": coeffZ}
    remainder_coeffs = {}

    # Recursion for the Submatrices
    for c, submatrix in coefficients.items():
        subDec, subRem = PauliDecomposition(
            submatrix,
            f"{PauliStringInit}{c}",
            recursion_depth=recursion_depth - 1,
            threshold=threshold,
        )
        for key, value in subDec.items():
            decomposition[key] += value
        remainder_coeffs[c] = subRem

    remI = remainder_coeffs["I"]
    remX = remainder_coeffs["X"]
    remY = remainder_coeffs["Y"]
    remZ = remainder_coeffs["Z"]

    TL = remI + remZ
    BR = remI - remZ
    TR = remX + I * remY
    BL = remX - I * remY

    if isinstance(matrix, sympy.Matrix):
        remainder = TL.row_join(TR).col_join(BL.row_join(BR))
    elif sparse.issparse(matrix):
        remainder = sparse.bmat([[TL, TR], [BL, BR]], format=matrix.format)
    else:
        remainder = np.block([[TL, TR], [BL, BR]])

    return decomposition, remainder

`string_to_fermion_operators(pauli_str, fermion_ops)`

Convert a Pauli string to its Jordan-Wigner fermionic expression.

Implements the standard JW inverse with accumulated parity strings

X_i -> (Z_0 x ... x Z_{i-1}) (c^dagger_i + c_i) Y_i -> (Z_0 x ... x Z_{i-1}) i(c^dagger_i - c_i) Z_i -> -(2 c^dagger_i c_i - 1) I_i -> identity (but still contributes to the parity string)

Parameters

pauli_str : str Pauli string of length len(fermion_ops), e.g. "IXZ". fermion_ops : list of FermionOp Ordered fermionic modes.

Returns

sympy.Expr Fermionic operator expression (not yet normal-ordered).

Source code in second_quantization/pauli_strings.py

def string_to_fermion_operators(pauli_str: str, fermion_ops: list) -> sympy.Expr:
    """Convert a Pauli string to its Jordan-Wigner fermionic expression.

    Implements the standard JW inverse with accumulated parity strings:
        X_i  ->  (Z_0 x ... x Z_{i-1}) (c^dagger_i + c_i)
        Y_i  ->  (Z_0 x ... x Z_{i-1}) i(c^dagger_i - c_i)
        Z_i  ->  -(2 c^dagger_i c_i - 1)
        I_i  ->  identity (but still contributes to the parity string)

    Parameters
    ----------
    pauli_str : str
        Pauli string of length ``len(fermion_ops)``, e.g. ``"IXZ"``.
    fermion_ops : list of FermionOp
        Ordered fermionic modes.

    Returns
    -------
    sympy.Expr
        Fermionic operator expression (not yet normal-ordered).
    """
    expr = S.One
    parity = S.One  # accumulated Z-string (JW parity prefactor)

    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))
        # accumulate parity string for the next site
        parity = parity * (-(2 * Dagger(op) * op - 1))

    return expr