Repository: qiskit-community/ffsim (python)

Issue:
Incorrect `norb` computation in `qiskit.jordan_wigner` when `n_qubits` is specified
### Issue
The function `jordan_wigner` computes `norb` using the expression in [line 46](https://github.com/qiskit-community/ffsim/blob/main/python/ffsim/qiskit/jordan_wigner.py#L46):

`norb = 1 + max(orb for term in op for _, _, orb in term)`

This determines the number of spatial orbitals based on the maximum orbital index found in `op`. However, when the user explicitly specifies `n_qubits`, the function should instead derive `norb` from `n_qubits // 2`.

Currently, this leads to incorrect indexing when applying `_qubit_action` in [line 62](https://github.com/qiskit-community/ffsim/blob/main/python/ffsim/qiskit/jordan_wigner.py#L62)

`qubit_op @= _qubit_action(action, orb + spin * norb, n_qubits)`

### Proposed solution
If `n_qubits` is provided, `norb` should be computed as:

```
if n_qubits is not None:
    norb = n_qubits // 2
```

### Minimal example
```
ffsim.qiskit.jordan_wigner(
    op={
        # ((True, False, 5), (True, False, 2), (False, False, 8), (False, False, 3)): (
        #     0.00022631523064525592 + 0j
        # ),
        ((True, False, 4), (True, True, 3), (False, True, 5), (False, False, 5)): (
            0.0015093838884580508 + 0j
        ),
    },
    n_qubits=2 * norb,
).to_sparse_list()
```

gives

```
[('YXYZX', [4, 5, 9, 10, 11], np.complex128(-9.433649302862818e-05+0j)),
 ('YYYZX', [4, 5, 9, 10, 11], np.complex128(-9.433649302862818e-05j)),
 ('YXYZY', [4, 5, 9, 10, 11], np.complex128(-9.433649302862818e-05j)),
 ('YYYZY', [4, 5, 9, 10, 11], np.complex128(9.433649302862818e-05+0j)),
 ('YXXZX', [4, 5, 9, 10, 11], np.complex128(-9.433649302862818e-05j)),
 ('YYXZX', [4, 5, 9, 10, 11], np.complex128(9.433649302862818e-05+0j)),
 ('YXXZY', [4, 5, 9, 10, 11], np.complex128(9.433649302862818e-05+0j)),
 ('YYXZY', [4, 5, 9, 10, 11], np.complex128(9.433649302862818e-05j)),
 ('XXYZX', [4, 5, 9, 10, 11], np.complex128(-9.433649302862818e-05j)),
 ('XYYZX', [4, 5, 9, 10, 11], np.complex128(9.433649302862818e-05+0j)),
 ('XXYZY', [4, 5, 9, 10, 11], np.complex128(9.433649302862818e-05+0j)),
 ('XYYZY', [4, 5, 9, 10, 11], np.complex128(9.433649302862818e-05j)),
 ('XXXZX', [4, 5, 9, 10, 11], np.complex128(9.433649302862818e-05+0j)),
 ('XYXZX', [4, 5, 9, 10, 11], np.complex128(9.433649302862818e-05j)),
 ('XXXZY', [4, 5, 9, 10, 11], np.complex128(9.433649302862818e-05j)),
 ('XYXZY', [4, 5, 9, 10, 11], np.complex128(-9.433649302862818e-05+0j))]
```

but un-commenting the first term gives (note the different support)

```
 [('YXYZX', [4, 5, 12, 13, 14], np.complex128(-9.433649302862818e-05+0j)),
 ('YYYZX', [4, 5, 12, 13, 14], np.complex128(-9.433649302862818e-05j)),
 ('YXYZY', [4, 5, 12, 13, 14], np.complex128(-9.433649302862818e-05j)),
 ('YYYZY', [4, 5, 12, 13, 14], np.complex128(9.433649302862818e-05+0j)),
 ('YXXZX', [4, 5, 12, 13, 14], np.complex128(-9.433649302862818e-05j)),
 ('YYXZX', [4, 5, 12, 13, 14], np.complex128(9.433649302862818e-05+0j)),
 ('YXXZY', [4, 5, 12, 13, 14], np.complex128(9.433649302862818e-05+0j)),
 ('YYXZY', [4, 5, 12, 13, 14], np.complex128(9.433649302862818e-05j)),
 ('XXYZX', [4, 5, 12, 13, 14], np.complex128(-9.433649302862818e-05j)),
 ('XYYZX', [4, 5, 12, 13, 14], np.complex128(9.433649302862818e-05+0j)),
 ('XXYZY', [4, 5, 12, 13, 14], np.complex128(9.433649302862818e-05+0j)),
 ('XYYZY', [4, 5, 12, 13, 14], np.complex128(9.433649302862818e-05j)),
 ('XXXZX', [4, 5, 12, 13, 14], np.complex128(9.433649302862818e-05+0j)),
 ('XYXZX', [4, 5, 12, 13, 14], np.complex128(9.433649302862818e-05j)),
 ('XXXZY', [4, 5, 12, 13, 14], np.complex128(9.433649302862818e-05j)),
 ('XYXZY', [4, 5, 12, 13, 14], np.complex128(-9.433649302862818e-05+0j))]
```

which shows how the `norb` for second term is derived by the max. of the whole term.

### Environment
MacOS, `ffsim 0.0.51`

Pull request:
Replace n_qubits argument in ffsim.qiskit.jordan_wigner with norb
Fixes #373 .

Originally, I envisioned `n_qubits` allowing the user to pad the operator with extra qubits. I don't see an actual use case for this though, and #373 shows that it's more important to be able to specify the actual number of orbitals. In the new code, the number of qubits is always twice the number of spatial orbitals. If necessary, the user can pad with extra qubits in post-processing.

Patch:
diff --git a/python/ffsim/qiskit/jordan_wigner.py b/python/ffsim/qiskit/jordan_wigner.py
index 3b1f1316..6239c9db 100644
--- a/python/ffsim/qiskit/jordan_wigner.py
+++ b/python/ffsim/qiskit/jordan_wigner.py
@@ -19,59 +19,75 @@ from qiskit.quantum_info import SparsePauliOp
 from ffsim.operators import FermionOperator
 
 
-def jordan_wigner(op: FermionOperator, n_qubits: int | None = None) -> SparsePauliOp:
-    """Jordan-Wigner transformation.
+def jordan_wigner(op: FermionOperator, norb: int | None = None) -> SparsePauliOp:
+    r"""Jordan-Wigner transformation.
 
     Transform a fermion operator to a qubit operator using the Jordan-Wigner
-    transformation.
+    transformation. The Jordan-Wigner transformation maps fermionic annihilation
+    operators to qubits as follows:
+
+    .. math::
+
+        a_p \mapsto \frac12 (X_p + iY_p)Z_1 \cdots Z_{p-1}
+
+    In the transformed operator, the first ``norb`` qubits represent spin-up (alpha)
+    orbitals, and the latter ``norb`` qubits represent spin-down (beta) orbitals. As a
+    result of this convention, the qubit index that an orbital is mapped to depends on
+    the total number of spatial orbitals. By default, the total number of spatial
+    orbitals is automatically determined by the largest-index orbital present in the
+    operator, but you can manually specify the number using the `norb` argument.
 
     Args:
         op: The fermion operator to transform.
-        n_qubits: The number of qubits to include in the output qubit operator. If not
-            specified, the minimum number of qubits needed to accommodate the fermion
-            operator will be used. Must be non-negative.
+        norb: The total number of spatial orbitals. If not specified, it is determined
+            by the largest-index orbital present in the operator.
 
     Returns:
         The qubit operator as a Qiskit SparsePauliOp.
 
     Raises:
-        ValueError: Number of qubits was negative.
-        ValueError: Number of qubits was not enough to accommodate the fermion operator.
+        ValueError: Number of spatial orbitals was negative.
+        ValueError: Number of spatial orbitals was fewer than the number detected in the
+            operator.
     """
-    if n_qubits and n_qubits < 0:
-        raise ValueError(f"Number of qubits must be non-negative. Got {n_qubits}.")
+    if norb and norb < 0:
+        raise ValueError(
+            f"Number of spatial orbitals must be non-negative. Got {norb}."
+        )
     if not op:
-        return SparsePauliOp.from_sparse_list([("", [], 0.0)], num_qubits=n_qubits or 0)
+        return SparsePauliOp.from_sparse_list(
+            [("", [], 0.0)], num_qubits=2 * (norb or 0)
+        )
 
-    norb = 1 + max(orb for term in op for _, _, orb in term)
-    if n_qubits is None:
-        n_qubits = 2 * norb
-    if n_qubits < 2 * norb:
+    norb_in_op = 1 + max(orb for term in op for _, _, orb in term)
+    if norb is None:
+        norb = norb_in_op
+    if norb < norb_in_op:
         raise ValueError(
-            "Number of qubits is not enough to accommodate the fermion operator. "
-            f"The fermion operator has {norb} spatial orbitals, so at least {2 * norb} "
-            f"qubits is needed, but got {n_qubits}."
+            "Number of spatial orbitals specified is fewer than the number detected in "
+            f"the operator. The operator has {norb_in_op} spatial orbitals, but "
+            f"only {norb} were specified."
         )
 
-    qubit_terms = [SparsePauliOp.from_sparse_list([("", [], 0.0)], num_qubits=n_qubits)]
+    qubit_terms = [SparsePauliOp.from_sparse_list([("", [], 0.0)], num_qubits=2 * norb)]
     for term, coeff in op.items():
         qubit_op = SparsePauliOp.from_sparse_list(
-            [("", [], coeff)], num_qubits=n_qubits
+            [("", [], coeff)], num_qubits=2 * norb
         )
         for action, spin, orb in term:
-            qubit_op @= _qubit_action(action, orb + spin * norb, n_qubits)
+            qubit_op @= _qubit_action(action, orb + spin * norb, norb)
         qubit_terms.append(qubit_op)
 
     return SparsePauliOp.sum(qubit_terms).simplify()
 
 
 @functools.cache
-def _qubit_action(action: bool, qubit: int, n_qubits: int):
+def _qubit_action(action: bool, qubit: int, norb: int):
     qubits = list(range(qubit + 1))
     return SparsePauliOp.from_sparse_list(
         [
             ("Z" * qubit + "X", qubits, 0.5),
             ("Z" * qubit + "Y", qubits, -0.5j if action else 0.5j),
         ],
-        num_qubits=n_qubits,
+        num_qubits=2 * norb,
     )


Test patch:
diff --git a/tests/python/qiskit/jordan_wigner_test.py b/tests/python/qiskit/jordan_wigner_test.py
index b79979f0..95830d0a 100644
--- a/tests/python/qiskit/jordan_wigner_test.py
+++ b/tests/python/qiskit/jordan_wigner_test.py
@@ -24,16 +24,16 @@ def test_random(norb: int, nelec: tuple[int, int]):
     op = ffsim.random.random_fermion_hamiltonian(norb, seed=rng)
     linop = ffsim.linear_operator(op, norb=norb, nelec=nelec)
     vec = ffsim.random.random_state_vector(ffsim.dim(norb, nelec), seed=rng)
+    expected_result = ffsim.qiskit.ffsim_vec_to_qiskit_vec(linop @ vec, norb, nelec)
 
     qubit_op = ffsim.qiskit.jordan_wigner(op)
     qubit_op_sparse = qubit_op.to_matrix(sparse=True)
     actual_result = qubit_op_sparse @ ffsim.qiskit.ffsim_vec_to_qiskit_vec(
         vec, norb, nelec
     )
-    expected_result = ffsim.qiskit.ffsim_vec_to_qiskit_vec(linop @ vec, norb, nelec)
     np.testing.assert_allclose(actual_result, expected_result, atol=1e-12)
 
-    qubit_op = ffsim.qiskit.jordan_wigner(op, n_qubits=2 * norb)
+    qubit_op = ffsim.qiskit.jordan_wigner(op, norb=norb)
     qubit_op_sparse = qubit_op.to_matrix(sparse=True)
     actual_result = qubit_op_sparse @ ffsim.qiskit.ffsim_vec_to_qiskit_vec(
         vec, norb, nelec
@@ -41,13 +41,36 @@ def test_random(norb: int, nelec: tuple[int, int]):
     np.testing.assert_allclose(actual_result, expected_result, atol=1e-12)
 
 
-def test_bad_qubit_number():
-    """Test passing bad number of qubits raises errors."""
+def test_more_norb():
+    """Test specifying more spatial orbitals than are present in the operator."""
+    op = ffsim.FermionOperator(
+        {
+            (ffsim.cre_a(1), ffsim.cre_b(1)): 1.0,
+            (ffsim.cre_a(2), ffsim.cre_b(2)): 1.0,
+        }
+    )
+    qubit_op = ffsim.qiskit.jordan_wigner(op, norb=5)
+    assert sorted(qubit_op.to_sparse_list()) == sorted(
+        [
+            ("YZZZZX", [1, 2, 3, 4, 5, 6], np.complex128(-0.25j)),
+            ("YZZZZY", [1, 2, 3, 4, 5, 6], np.complex128(-0.25 + 0j)),
+            ("XZZZZX", [1, 2, 3, 4, 5, 6], np.complex128(0.25 + 0j)),
+            ("XZZZZY", [1, 2, 3, 4, 5, 6], np.complex128(-0.25j)),
+            ("YZZZZX", [2, 3, 4, 5, 6, 7], np.complex128(-0.25j)),
+            ("YZZZZY", [2, 3, 4, 5, 6, 7], np.complex128(-0.25 + 0j)),
+            ("XZZZZX", [2, 3, 4, 5, 6, 7], np.complex128(0.25 + 0j)),
+            ("XZZZZY", [2, 3, 4, 5, 6, 7], np.complex128(-0.25j)),
+        ]
+    )
+
+
+def test_bad_norb():
+    """Test passing bad number of spatial orbitals raises errors."""
     op = ffsim.FermionOperator({(ffsim.cre_a(3),): 1.0})
     with pytest.raises(ValueError, match="non-negative"):
-        _ = ffsim.qiskit.jordan_wigner(op, n_qubits=-1)
-    with pytest.raises(ValueError, match="enough"):
-        _ = ffsim.qiskit.jordan_wigner(op, n_qubits=7)
+        _ = ffsim.qiskit.jordan_wigner(op, norb=-1)
+    with pytest.raises(ValueError, match="fewer"):
+        _ = ffsim.qiskit.jordan_wigner(op, norb=3)
 
 
 def test_hubbard():