Closes #5895 Add matrix support for fermionic objects (#5920)

**Context:**
The [fermi](https://docs.pennylane.ai/en/stable/code/qml_fermi.html)
module includes functions and classes for creating and manipulating
fermionic operators: `FermiWord` and `FermiSentence` class. It would be
beneficial to add support for obtaining the matrix representation of
these operators, similar to the `to_mat` instance method in the
`PauliWord` and `PauliSentence` classes.

**Description of the Change:**

- `to_mat` instance method is added to `FermiWord` and `FermiSentence`
classes; the implementation for two methods is nearly identical -- using
Jordan-Wigner transform to convert an instance of `FermiWord` or
`FermiSentence` to a `PauliWord` or `PauliSentence` instance, then using
the existing `to_mat` method in the two classes to compute the matrix
representation.
- `to_mat` takes `n_orbitals` as an optional input, which determines the
matrix size $2^{n_\mathrm{orbitals}} \times 2^{n_\mathrm{orbitals}}$

**Benefits:**
The new method `to_mat` in `FermiWord` and `FermiSentence` class allows
computing the matrix representation of a given Fermi operator.

**Possible Drawbacks:**
- Currently, `to_mat` in `FermiWord` and `FermiSentence` class only
supports dense matrices. The exponential growth of matrix dimensions
(a.k.a Fock space) with respect to `n_orbitals` (or the largest orbital
index in the input fermionic operator) may limit its applicability.

**Related GitHub Issues:**
[#5895](#5895)

### Why Jordan-Wigner transform?
The choice of Jordan-Wigner mapping was primarily motivated by its
direct and straightforward representation, where **each Slater
determinant has a one-to-one mapping with the computational basis states
of qubits**. This feature simplifies the interpretation of matrix
elements significantly, as **each matrix element directly corresponds to
the expectation values of the fermionic operator between these
determinants**. Such a direct correspondence facilitates easier
understanding and analysis of quantum states, making it intuitive and
useful in the context of quantum chemical applications, as well as in
educational and debugging scenarios.

In contrast, while Parity and Bravyi-Kitaev mappings offer computational
efficiencies for larger systems and specific quantum hardware
considerations, they encode fermionic information into qubit states in
ways that are more complex and less intuitive. Parity mapping encodes
the parity of the number of occupied orbitals up to each position, and
Bravyi-Kitaev uses a hierarchical structure.
- _Parity_: Matrix elements in the parity-mapped Hamiltonian correspond
to interactions that conserve or flip these parity properties.
Physically, each element can be thought of as involving transitions that
maintain or alter the parity characteristics of the fermionic system,
rather than direct particle occupations.
- _Bravyi-Kitaev_: Matrix elements in a Bravyi-Kitaev transformed
Hamiltonian reflect operations on a hierarchical structure of fermionic
occupations. Each qubit state relates to a specific combination of
fermion occupation and non-occupation across different levels of this
hierarchy, with matrix elements representing more complex interactions
across these levels.

The physical interpretation of each matrix element in the two mappings
requires a deeper understanding of the encoded structures, making them
less transparent for straightforward interpretations (e.g. A
transformation matrix may be necessary to establish the correspondence
between qubit computational basis states and Slater determinants).

Including the option to select different mappings could potentially
offer greater flexibility and optimization for users dealing with
various types of quantum hardware or larger systems. However, it would
also introduce additional complexity into our codebase and could
complicate the user experience for those who are not as familiar with
the subtleties of each mapping method. Therefore, I chose to focus
exclusively on the Jordan-Wigner mapping for its clarity and ease of
interpretation.


#### Notes on the design (OUTDATED, `_to_mat` function has been removed)
The primary motivation for centralizing the main logic in the `_to_mat`
function, with both `FermiWord.to_mat()` and `FermiSentence.to_mat()`
methods calling this function, is to reduce the code duplication and
enhance maintainability. The implementation for these methods is nearly
identical; hence, modularizing the code into a separate, reusable
function that accepts either a `FermiWord` or `FermiSentence` instance
seemed like the most efficient approach. Here are several advantages of
the current implementation:

1. **Modularity and Reduced Code Duplication**: As mentioned, this
approach helps in avoiding repetition of similar code across multiple
classes, adhering to the DRY (Don't Repeat Yourself) principle.

2. **Flexibility in Usage**: By potentially making this function public
(renaming it to `to_mat` and importing it in the `__init__.py`), users
gain flexibility. They can either use the method associated with an
instance:
   ```python
   word = FermiWord({(0, 0): "+", (1, 1): "-"})
   word.to_mat()
   ```
   Or directly call the function:
   ```python
   from pennylane.fermi.fermionic import to_mat
   to_mat(word)
   ```
Extending this idea, we could further develop it into a versatile
`to_mat` function that accepts both Fermi and Pauli operators. This
function would have the signature:
```python
def to_mat(operator: Union[FermiWord, FermiSentence, PauliWord, PauliSentence], **kwargs) -> ndarray
```
Such a design aligns with Python best practices of maintaining
lightweight classes, as functions are generally easier to test than
instance methods because they are less dependent on instance state. This
approach not only simplifies testing but also enhances code usability
and readability.

---------

Co-authored-by: RenkeHuang <renke_huang@outlook.com>
Co-authored-by: Thomas R. Bromley <49409390+trbromley@users.noreply.github.com>
Co-authored-by: soranjh <soran.jahangiri@gmail.com>

doc/releases/changelog-dev.md

@@ -21,6 +21,10 @@
 
 * `SProd.terms` now flattens out the terms if the base is a multi-term observable.
   [(#5885)](https://github.com/PennyLaneAI/pennylane/pull/5885)
+
+* A new method `to_mat` has been added to the `FermiWord` and `FermiSentence` classes, which allows
+  computing the matrix representation of these Fermi operators.
+  [(#5920)](https://github.com/PennyLaneAI/pennylane/pull/5920)
 
 <h3>Improvements 🛠</h3>
 
@@ -146,6 +150,7 @@ Ahmed Darwish,
 Lillian M. A. Frederiksen,
 Pietropaolo Frisoni,
 Emiliano Godinez,
+Renke Huang,
 Christina Lee,
 Austin Huang,
 Christina Lee,

pennylane/fermi/fermionic.py

@@ -273,6 +273,36 @@ def __pow__(self, value):
 
         return operator
 
+    def to_mat(self, n_orbitals=None):
+        r"""Return the matrix representation.
+
+        Args:
+            n_orbitals (int or None): Number of orbitals. If not provided, it will be inferred from
+                the largest orbital index in the Fermi operator.
+
+        Returns:
+            NumpyArray: Matrix representation of the :class:`~.FermiWord`.
+
+        **Example**
+
+        >>> w = FermiWord({(0, 0): '+', (1, 1): '-'})
+        >>> w.to_mat()
+        array([0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
+              [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
+              [0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j],
+              [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j])
+        """
+        largest_orb_id = max(key[1] for key in self.keys()) + 1
+        if n_orbitals and n_orbitals < largest_orb_id:
+            raise ValueError(
+                f"n_orbitals cannot be smaller than {largest_orb_id}, got: {n_orbitals}."
+            )
+
+        largest_order = n_orbitals or largest_orb_id
+        mat = qml.jordan_wigner(self, ps=True).to_mat(wire_order=list(range(largest_order)))
+
+        return mat
+
 
 # pylint: disable=useless-super-delegation
 class FermiSentence(dict):
@@ -463,6 +493,36 @@ def simplify(self, tol=1e-8):
             if abs(coeff) <= tol:
                 del self[fw]
 
+    def to_mat(self, n_orbitals=None):
+        r"""Return the matrix representation.
+
+        Args:
+            n_orbitals (int or None): Number of orbitals. If not provided, it will be inferred from
+                the largest orbital index in the Fermi operator
+
+        Returns:
+            NumpyArray: Matrix representation of the :class:`~.FermiSentence`.
+
+        **Example**
+
+        >>> fs = FermiSentence({FermiWord({(0, 0): "+", (1, 1): "-"}): 1.2, FermiWord({(0, 0): "+", (1, 0): "-"}): 3.1})
+        >>> fs.to_mat()
+        array([0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j],
+              [0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j],
+              [0.0 + 0.0j, 1.2 + 0.0j, 3.1 + 0.0j, 0.0 + 0.0j],
+              [0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, 3.1 + 0.0j])
+        """
+        largest_orb_id = max(key[1] for fermi_word in self.keys() for key in fermi_word.keys()) + 1
+        if n_orbitals and n_orbitals < largest_orb_id:
+            raise ValueError(
+                f"n_orbitals cannot be smaller than {largest_orb_id}, got: {n_orbitals}."
+            )
+
+        largest_order = n_orbitals or largest_orb_id
+        mat = qml.jordan_wigner(self, ps=True).to_mat(wire_order=list(range(largest_order)))
+
+        return mat
+
 
 def from_string(fermi_string):
     r"""Return a fermionic operator object from its string representation.

tests/fermi/test_fermionic.py

@@ -125,6 +125,22 @@ def test_init_error(self, operator):
         with pytest.raises(ValueError, match="The operator indices must belong to the set"):
             FermiWord(operator)
 
+    def test_to_mat(self):
+        """Test that the matrix representation of FermiWord is correct."""
+
+        expected_mat = np.zeros((4, 4), dtype=complex)
+        expected_mat[2, 1] = 1.0
+
+        mat = fw1.to_mat()
+        assert np.allclose(mat, expected_mat)
+
+    def test_to_mat_error(self):
+        """Test that an error is raised if the requested matrix dimension is smaller than the
+        dimension inferred from the largest orbital index.
+        """
+        with pytest.raises(ValueError, match="n_orbitals cannot be smaller than 2"):
+            fw1.to_mat(n_orbitals=1)
+
 
 class TestFermiWordArithmetic:
     WORDS_MUL = (
@@ -442,6 +458,14 @@ def test_array_must_not_exceed_length_1(self, method_name):
 fs3 = FermiSentence({fw3: -0.5, fw4: 1})
 fs4 = FermiSentence({fw4: 1})
 fs5 = FermiSentence({})
+fs6 = FermiSentence({fw1: 1.2, fw2: 3.1})
+fs7 = FermiSentence(
+    {
+        FermiWord({(0, 0): "+", (1, 1): "-"}): 1.23,  # a+(0) a(1)
+        FermiWord({(0, 0): "+", (1, 0): "-"}): 4.0j,  # a+(0) a(0) = n(0) (number operator)
+        FermiWord({(0, 0): "+", (1, 2): "-", (2, 1): "+"}): -0.5,  # a+(0) a(2) a+(1)
+    }
+)
 
 fs1_x_fs2 = FermiSentence(  # fs1 * fs1, computed by hand
     {
@@ -600,6 +624,25 @@ def test_pickling(self):
         new_fs = pickle.loads(serialization)
         assert fs == new_fs
 
+    def test_to_mat(self):
+        """Test that the matrix representation of FermiSentence is correct."""
+        expected_mat = np.zeros((8, 8), dtype=complex)
+        expected_mat[4, 2] = 1.23 + 0j
+        expected_mat[5, 3] = 1.23 + 0j
+        for i in [4, 5, 6, 7]:
+            expected_mat[i, i] = 4.0j
+        expected_mat[6, 1] = 0.5 + 0j
+
+        mat = fs7.to_mat()
+        assert np.allclose(mat, expected_mat)
+
+    def test_to_mat_error(self):
+        """Test that an error is raised if the requested matrix dimension is smaller than the
+        dimension inferred from the largest orbital index.
+        """
+        with pytest.raises(ValueError, match="n_orbitals cannot be smaller than 3"):
+            fs7.to_mat(n_orbitals=2)
+
 
 class TestFermiSentenceArithmetic:
     tup_fs_mult = (  # computed by hand