Second-Order Twirl & Quadratic Symmetries
###########################################

This page details the computation of the second-order Pauli twirl,
a key operation in quantum information theory for characterizing noise
and benchmarking quantum devices. The calculation relies on the concept
of "quadratic symmetries" of the underlying gate set.

Theoretical Background
======================

The second-order twirl of a quantum channel or operator :math:`M` over a group :math:`G`
is given by its projection onto the commutant of the twofold tensor product representation
of the group, :math:`U^{\otimes 2}`.

Supplemental Theorem 1 of `arXiv:2310.11505 <https://arxiv.org/abs/2310.11505>`_
provides a powerful way to compute this by constructing an orthogonal basis for this space,
called the **quadratic symmetries**, :math:`Q_{kj}`.

**Definitions:**

*   **Linear Symmetries (:math:`L_j`):** The set of Pauli strings that commute with all generators of the system :math:`G`. This is found using ``system.get_commutants()``.
*   **Commutator Graph Components (:math:`C_k`):** The connected components of the graph of all Pauli strings, where edges are defined by commutation with the generators of :math:`G`. This is found using the internal ``system.get_commutator_graph_components()`` method.

The quadratic symmetries are then constructed as:

.. math::

   Q_{kj} = \sum_{S \in C_k} S \otimes L_j S

The final second-moment twirl of an operator :math:`M` is its projection onto this basis:

.. math::

   T^{(2)}(M) = \sum_{j,k} \frac{\text{Tr}(Q_{kj} M)}{\text{Tr}(Q_{kj}^2)} Q_{kj}

``PauLie`` provides the ``second_moment`` application function to perform this entire calculation.

Usage Example
=============

The main user-facing function is ``second_moment``, which takes a system definition and an operator to twirl. Let's use a simple system generated by the Pauli-Z operator and see how it affects the operator :math:`X \otimes I`.

.. code-block:: python

   from paulie.common.pauli_string_factory import get_pauli_string as p
   from paulie.application.second_moment import second_moment

   # Define the system generators (a group generated by Z on 1 qubit)
   system = p(["Z"])
   
   # Define an operator M to twirl.
   # Note: M must be a 2n-qubit operator. Since our system is n=1, M is 2-qubit.
   M = p([(1.0, "XI")])

   # Compute the second-order twirl. The function handles the construction
   # of the quadratic symmetry basis internally.
   twirled_M = second_moment(M, system)

   print(f"Original operator M: {M}")
   print(f"Twirled operator T^(2)(M): {twirled_M}")

**Output:**

.. code-block:: text

   Original operator M: XI
   Twirled operator T^(2)(M): 0.0*II

In this case, the result is the zero operator because the operator :math:`X \otimes I` has no overlap with any of the quadratic symmetries generated by the Z-group. The ``second_moment`` function correctly projects it to zero.