.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/rmt/paper_replications/syk_demo.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_api_gallery_rmt_paper_replications_syk_demo.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_api_gallery_rmt_paper_replications_syk_demo.py:


Kitaev's SYK Model: Majorana Fermions, Wigner-Dyson Statistics
=================================================================

Unlike every classical ensemble in this package, the Sachdev-Ye-Kitaev
(SYK) model is not built from i.i.d. matrix entries: it is an explicit
many-body Hamiltonian acting on :math:`N` Majorana fermions
:math:`\gamma_1, \dots, \gamma_N` (Hermitian operators satisfying the
Clifford algebra :math:`\{\gamma_a, \gamma_b\} = 2\delta_{ab}`),

.. math::

    H = i^{q/2} \sum_{i_1 < \cdots < i_q} J_{i_1 \cdots i_q}\,
    \gamma_{i_1} \cdots \gamma_{i_q},

with independent random couplings :math:`J_{i_1\cdots i_q} \sim
\mathcal{N}\!\left(0,\, \tfrac{(q-1)!\, J^2}{N^{q-1}}\right)` for every
:math:`q`-index tuple (:math:`q=4` here, the standard maximally chaotic
case). This example builds the Sachdev-Ye-Kitaev Hamiltonian directly
from ``physicskit.rmt.ensembles.syk.majorana_operators`` and verifies
the Clifford algebra :math:`\{\gamma_a, \gamma_b\} = 2\delta_{ab}` it
must satisfy to machine precision. Then it shows the qualitative
signature the SYK model is famous for: even though :math:`H` is not
built from i.i.d. matrix entries at all, its many-body level statistics
-- once restricted to a single fermion-parity sector, the two of which
are otherwise uncorrelated and would wash out any repulsion signal if
pooled together -- are nonetheless Wigner-Dyson, not Poisson.

References:
A. Kitaev, unpublished KITP talks (2015).
Y.-Z. You, A. W. W. Ludwig, C. Xu, Phys. Rev. B 95 (2017) 115150 --
the N mod 8 symmetry-class periodicity (N=16 here falls in the
GOE class).

Run:
    python examples/paper_replications/syk_demo.py

.. GENERATED FROM PYTHON SOURCE LINES 38-138



.. image-sg:: /api/gallery/rmt/paper_replications/images/sphx_glr_syk_demo_001.png
   :alt: Kitaev's SYK model: a genuine many-body Hamiltonian with Wigner-Dyson statistics, Many-body spectrum, single parity sector (N=16, 20 realizations), Many-body spectrum: Wigner-Dyson, not Poisson
   :srcset: /api/gallery/rmt/paper_replications/images/sphx_glr_syk_demo_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Clifford algebra {gamma_a, gamma_b} = 2*delta_ab verified for N=8: max deviation 0.00e+00 (machine precision)
    Saved syk_replication.png






|

.. code-block:: Python


    import matplotlib.pyplot as plt
    import numpy as np
    from scipy.stats import kstest

    import physicskit.rmt as rmt
    from physicskit.rmt.ensembles.syk import _sample_syk_hamiltonian, majorana_operators

    SEED = 2026

    # --- The Clifford algebra, verified directly (not plotted: the deviation
    # is machine precision everywhere, so a heatmap of it is just one flat
    # color and shows nothing) ---
    N_CLIFFORD = 8
    gammas = majorana_operators(N_CLIFFORD)
    anticommutator_error = 0.0
    for a in range(N_CLIFFORD):
        for b in range(N_CLIFFORD):
            target = 2.0 * np.eye(gammas[0].shape[0]) if a == b else 0.0
            anticommutator_error = max(
                anticommutator_error,
                np.abs(gammas[a] @ gammas[b] + gammas[b] @ gammas[a] - target).max(),
            )
    assert anticommutator_error < 1e-10
    print(f"Clifford algebra {{gamma_a, gamma_b}} = 2*delta_ab verified for N={N_CLIFFORD}: max deviation {anticommutator_error:.2e} (machine precision)")

    fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))

    # --- Both panels below share one set of many-body spectra: the density
    # of states (panel 1) and the level-spacing ratios drawn from it
    # (panel 2), each computed once per trial in a single parity sector ---
    N = 16
    Q = 4
    N_TRIALS = 20

    gammas = majorana_operators(N)
    dim = gammas[0].shape[0]
    # Fermion parity operator: the Hermitian product of all N Majoranas
    # (same i^(N/2)-phase convention the module uses to Hermitize odd
    # products of Majoranas); H (built from even q=4 products) commutes
    # with it, so its +1 eigenspace is an invariant, physically meaningful
    # sector -- pooling both parity sectors instead would superpose two
    # mutually uncorrelated spectra and wash out level repulsion.
    parity = np.eye(dim, dtype=complex)
    for gamma in gammas:
        parity = parity @ gamma
    parity = parity * (1j ** (N // 2))
    parity_eigvals, parity_eigvecs = np.linalg.eigh(parity)
    sector = parity_eigvecs[:, np.round(parity_eigvals.real) == 1]

    rng = np.random.default_rng(SEED)
    ratios = []
    pooled_eigs = []
    for _ in range(N_TRIALS):
        h = _sample_syk_hamiltonian(N, Q, 1.0, rng, gammas)
        h_sector = sector.conj().T @ h @ sector
        h_sector = (h_sector + h_sector.conj().T) / 2.0
        eigs = np.linalg.eigvalsh(h_sector)
        spacings = np.diff(eigs)
        ratios.append(np.minimum(spacings[:-1], spacings[1:]) / np.maximum(spacings[:-1], spacings[1:]))
        # Standardize each disorder realization (zero mean, unit std) before
        # pooling, so trial-to-trial fluctuations in overall spectral shift
        # and width don't smear out the pooled density of states.
        pooled_eigs.append((eigs - eigs.mean()) / eigs.std())
    ratios = np.concatenate(ratios)
    pooled_eigs = np.concatenate(pooled_eigs)

    goe_surmise = rmt.stats.RatioSurmise(beta=1)

    ax = axes[0]
    ax.hist(pooled_eigs, bins=60, density=True, color="steelblue", alpha=0.5)
    ax.set_xlabel("E (standardized per disorder realization)")
    ax.set_ylabel("density of states")
    ax.set_title(f"Many-body spectrum, single parity sector (N={N}, {N_TRIALS} realizations)")


    def poisson_ratio_cdf(r: np.ndarray) -> np.ndarray:
        return 2.0 * r / (1.0 + r)


    ks_goe = kstest(ratios, goe_surmise.cdf).statistic
    ks_poisson = kstest(ratios, poisson_ratio_cdf).statistic

    r_grid = np.linspace(0, 1, 400)
    ax = axes[1]
    ax.hist(ratios, bins=40, density=True, alpha=0.5, color="steelblue", label=f"SYK (N={N}, single parity sector)")
    ax.plot(r_grid, goe_surmise.pdf(r_grid), "k-", lw=2, label=f"GOE surmise (KS={ks_goe:.3f})")
    ax.plot(r_grid, 2.0 / (1.0 + r_grid) ** 2, "k--", lw=2, label=f"Poisson (KS={ks_poisson:.3f})")
    ax.set_xlabel("r (consecutive spacing ratio)")
    ax.set_ylabel("density P(r)")
    ax.set_title("Many-body spectrum: Wigner-Dyson, not Poisson")
    ax.legend(fontsize=8)

    fig.suptitle(
        "Kitaev's SYK model: a genuine many-body Hamiltonian with Wigner-Dyson statistics",
    )
    fig.tight_layout()
    out_path = "syk_replication.png"
    fig.savefig(out_path, dpi=150)
    print(f"Saved {out_path}")


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 48.857 seconds)


.. _sphx_glr_download_api_gallery_rmt_paper_replications_syk_demo.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: syk_demo.ipynb <syk_demo.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: syk_demo.py <syk_demo.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: syk_demo.zip <syk_demo.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_