{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ff8465b6",
   "metadata": {},
   "source": [
    "| [**Lattice QCD I**](https://homepages.uni-regensburg.de/~lec17310/teaching/wise2122/lqft.html) | - | [**Christoph Lehner**](https://go.ur.de/lehner) |\n",
    "|:-|--|-:|\n",
    "\n",
    "# Chapter 11: symmetries of the Wilson action\n",
    "\n",
    "In this chapter we study in more detail the symmetries of the lattice action of QCD.  In order to do so, we first remind ourselves of the full lattice action with gauge links $U_\\mu$ and a single fermion field $\\psi$.\n",
    "We have\n",
    "$$\n",
    "\\begin{align}\n",
    "S = \\sum_x \\bar\\psi(x) \\left( \\sum_\\mu \\gamma_\\mu D_\\mu + m - \\frac12 \\sum_\\mu \\overleftarrow{D}_\\mu\\overrightarrow{D}_\\mu\\right) \\psi(x) + \\beta \\sum_x \\sum_{\\mu<\\nu} ( 1 - P_{\\mu\\nu}(x)) \n",
    "\\end{align}\n",
    "$$\n",
    "with\n",
    "$$\n",
    "\\begin{align}\n",
    "D_\\mu &= \\frac12 \\left(\\overleftarrow{D}_\\mu + \\overrightarrow{D}_\\mu\\right)\n",
    "= \\frac12 \\left(C_\\mu^+ - C_\\mu^-\\right) \\,, \\\\\n",
    "\\overleftarrow{D}_\\mu \\overrightarrow{D}_\\mu &= C_\\mu^+ + C_\\mu^- - 2 \\,, \\\\\n",
    "C^+_\\mu \\psi(x) &= U_\\mu(x) \\psi(x + a \\hat{\\mu}) \\,, \\\\\n",
    "C^-_\\mu \\psi(x) &= U^\\dagger_\\mu(x - a\\hat{\\mu}) \\psi(x - a \\hat{\\mu}) \\,, \\\\\n",
    "P_{\\mu\\nu}(x) &= \\frac1{N_c} \\textrm{Re}\\: \\textrm{Tr}\\: U_{\\mu\\nu}(x) \\,, \\\\\n",
    "U_{\\mu\\nu}(x) &= U_\\mu(x) U_\\nu(x + a \\hat{\\mu})U_\\mu(x+a \\hat{\\nu})^\\dagger U_\\nu(x)^\\dagger  \\,.\n",
    "\\end{align}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e1e68af",
   "metadata": {},
   "source": [
    "## Particle and antiparticle solutions in Minkowski space\n",
    "\n",
    "In chapter 5 we learned that the wavefunction of a free Dirac particle in Minkowski space\n",
    "satisfies\n",
    "$$\n",
    "\\begin{align}\n",
    "(i \\gamma^0 \\partial_0 - m)\\psi(t) = 0\n",
    "\\end{align}\n",
    "$$\n",
    "if we assume that it is invariant in space, i.e., $\\partial_i \\psi=0$.\n",
    "This equation has solutions\n",
    "$$\n",
    "\\begin{align}\n",
    "\\psi^+(t) = (1-\\gamma^0)e^{imt} \\psi^+_0\n",
    "\\end{align}\n",
    "$$\n",
    "and\n",
    "$$\n",
    "\\begin{align}\n",
    "\\psi^-(t) = (1+\\gamma^0)e^{-imt} \\psi^-_0\n",
    "\\end{align}\n",
    "$$\n",
    "since $\\gamma^0(1-\\gamma^0)=-(1-\\gamma^0)$ and $\\gamma^0(1+\\gamma^0)=1+\\gamma^0$.\n",
    "Here $\\psi^\\pm_0$ are vectors selecting a spin configuration.\n",
    "We can think of $\\psi^+$ as the **particle** and $\\psi^-$ as the **antiparticle** solution.\n",
    "\n",
    "We can now consider the **charge conjugation** transformation $C$ that maps particles to antiparticles and vice versa.\n",
    "Since $(\\gamma^0)^* = \\gamma^0$ and\n",
    "$$\n",
    "i \\gamma^2(1 \\pm \\gamma^0)i\\gamma^2 =  1 \\mp \\gamma^0\n",
    "$$\n",
    "we find\n",
    "$$\n",
    "\\begin{align}\n",
    "C\\psi^\\pm(t) &\\equiv -i\\gamma^2(\\psi^\\pm(t))^* = -i\\gamma^2 (1\\mp \\gamma^0)e^{\\mp imt} (\\psi^\\pm_0)^* \\\\\n",
    "&=i\\gamma^2 (1\\mp \\gamma^0)i\\gamma^2 e^{\\mp imt} i\\gamma^2 (\\psi^\\pm_0)^* \\\\\n",
    "&= (1 \\pm \\gamma^0) e^{\\mp imt} i\\gamma^2 (\\psi^\\pm_0)^* \\\\\n",
    "&=\\psi^\\mp(t)\n",
    "\\end{align}\n",
    "$$\n",
    "if we identify\n",
    "$$\n",
    "\\begin{align}\n",
    "i\\gamma^2 (\\psi^\\pm_0)^* = \\psi^\\mp_0 \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "The charge conjugation operation defined in this way therefore relates particle solutions of a given spin to antiparticle solutions of a different spin configuration.\n",
    "\n",
    "We can generalize this discussion beyond free-field solutions studying the charge conjugation operation\n",
    "$C$ that replaces\n",
    "$$\n",
    "\\begin{align}\n",
    "\\psi \\to \\psi^C &= -i\\gamma^2 \\psi^* = -i\\gamma^2 (\\psi^\\dagger)^T = -i\\gamma^2 (\\bar{\\psi}\\gamma^0)^T \\\\\n",
    "&= -i\\gamma^2 \\gamma^0 \\bar{\\psi}^T \\,.\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b11094c5",
   "metadata": {},
   "source": [
    "## Charge-conjugation symmetry\n",
    "\n",
    "In Euclidean space this transformation can be written as\n",
    "$$\n",
    "\\begin{align}\n",
    "\\psi \\to \\psi^C = C^{-1} \\bar\\psi^T \\,,\n",
    "\\end{align}\n",
    "$$\n",
    "where $C^{-1}=\\gamma^2 \\gamma^4$ with Euclidean gamma matrices.  Note that it is somewhat more\n",
    "convenient to define instead\n",
    "$$\n",
    "\\begin{align}\n",
    "C^{-1} = i\\gamma^2 \\gamma^4\n",
    "\\end{align}\n",
    "$$\n",
    "by including a global phase factor such that $C^{-1}=C=C^{\\dagger}=-C^T$.  It is straightforward\n",
    "to show that\n",
    "$$\n",
    "\\begin{align}\n",
    "C \\gamma^\\mu C = -(\\gamma^\\mu)^T\n",
    "\\end{align}\n",
    "$$\n",
    "which for the free theory again yields a symmetry of the action under\n",
    "$$\n",
    "\\begin{align}\n",
    "\\psi &\\to C \\bar\\psi^T \n",
    "\\end{align}\n",
    "$$\n",
    "such that\n",
    "$$\n",
    "\\begin{align}\n",
    "\\bar\\psi &\\to (C \\bar\\psi^T)^\\dagger \\gamma^4 = (\\bar\\psi^T)^\\dagger C \\gamma^4 \\\\\n",
    "&=-(\\bar\\psi^T)^\\dagger (\\gamma^4)^T C =-(\\gamma^4\\bar\\psi^\\dagger)^T C \\\\\n",
    "&=-\\psi^T C \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "The free action therefore transforms as\n",
    "$$\n",
    "\\begin{align}\n",
    "S &= \\sum_x \\bar\\psi(x) \\left( \\sum_\\mu \\gamma_\\mu \\partial_\\mu + m - \\frac12 \\sum_\\mu \\overleftarrow{\\partial}_\\mu\\overrightarrow{\\partial}_\\mu\\right) \\psi(x) \\\\\n",
    "\\to S^C &\\equiv -\\sum_x \\psi(x)^T C \\left( \\sum_\\mu \\gamma_\\mu \\partial_\\mu + m - \\frac12 \\sum_\\mu \\overleftarrow{\\partial}_\\mu\\overrightarrow{\\partial}_\\mu\\right) C \\bar\\psi(x)^T \\\\\n",
    "&= -\\sum_x \\psi(x)^T \\left( \\sum_\\mu (-\\gamma_\\mu)^T \\partial_\\mu + m - \\frac12 \\sum_\\mu \\overleftarrow{\\partial}_\\mu\\overrightarrow{\\partial}_\\mu\\right) \\bar\\psi(x)^T \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Since $\\psi$ and $\\bar\\psi$ are Grassmannian variables, they anticommute.  In addition\n",
    "for commuting variables\n",
    "$$\n",
    "\\begin{align}\n",
    "\\sum_x f(x) g(x+\\hat{\\mu}) = \\sum_x g(x)f(x-\\hat{\\mu})\n",
    "\\end{align}\n",
    "$$\n",
    "such that the first derivative term gets an additional minus sign.  Finally $v^T X^T w^T=w X v$.  In total, we find\n",
    "$$\n",
    "\\begin{align}\n",
    " S^C &= \\sum_x \\bar\\psi(x) \\left( \\sum_\\mu \\gamma_\\mu \\partial_\\mu + m - \\frac12 \\sum_\\mu \\overleftarrow{\\partial}_\\mu\\overrightarrow{\\partial}_\\mu\\right) \\psi(x) = S \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "The free Wilson fermion action is therefore invariant under charge conjugation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c08075e",
   "metadata": {},
   "source": [
    "In order to preserve a symmetry in the interacting case, we need \n",
    "the covariant shifts $C^\\pm$ to transform as\n",
    "$$\n",
    "\\begin{align}\n",
    "\\sum_x f(x) C^\\pm g(x) &\\to \\sum_x g(x)^T C^\\mp f(x)^T \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "Concretely the transformed gauge fields \n",
    "$$\n",
    "U_\\mu \\to U_\\mu^C\n",
    "$$ \n",
    "in $C^\\pm$ need to obey\n",
    "$$\n",
    "\\begin{align}\n",
    " \\sum_x f(x) U^C_\\mu(x) g(x + a \\hat{\\mu})\n",
    " =\\sum_x g(x)^T U^\\dagger_\\mu(x - a\\hat{\\mu}) f(x - a \\hat{\\mu})^T\n",
    " =\\sum_x g(x+a\\hat{\\mu})^T U^\\dagger_\\mu(x) f(x)^T\n",
    " =\\sum_x f(x) U^*_\\mu(x) g(x+a\\hat{\\mu}) \\,.\n",
    " \\end{align}\n",
    "$$\n",
    "This is satisfied if we demand that the gauge fields transform under\n",
    "charge conjugation symmetry as\n",
    "$$\n",
    "\\begin{align}\n",
    "U_\\mu &\\to U_\\mu^C = U_\\mu^* \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Finally, we need to check if the gauge-field part of the action also preserves this transformation.\n",
    "We note that under $U_\\mu \\to U_\\mu^* = (U_\\mu^T)^\\dagger$\n",
    "we have\n",
    "$$\n",
    "\\begin{align}\n",
    "U_{\\mu\\nu}(x) &\\to (U^\\dagger_\\mu(x))^T (U^\\dagger_\\nu(x + a \\hat{\\mu}))^T U_\\mu(x+a \\hat{\\nu})^T U_\\nu(x)^T \\\\\n",
    "&= (U_\\nu(x) U_\\mu(x+a \\hat{\\nu}) (U^\\dagger_\\nu(x + a \\hat{\\mu}))   (U^\\dagger_\\mu(x)))^T \\\\\n",
    "&=U_{\\nu\\mu}(x)^T\n",
    "\\end{align}\n",
    "$$\n",
    "such that\n",
    "$$\n",
    "\\begin{align}\n",
    "\\textrm{Tr}\\: U_{\\mu\\nu}(x)\n",
    "&\\to \\textrm{Tr}\\: U_{\\nu\\mu}(x) \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "Finally, we note that $U_{\\mu\\nu}(x)^\\dagger = U_{\\nu\\mu}(x)$ such that due to the real part\n",
    "$P_{\\mu\\nu}$ and therefore also the full Wilson action is invariant under charge conjugation symmetry."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8123849",
   "metadata": {},
   "source": [
    "## $\\gamma_5$-Hermiticity\n",
    "\n",
    "Next, we prove the $\\gamma_5$-Hermiticity (see chapter 10) of the Wilson Dirac operator.  The Dirac operator is given by\n",
    "$$\n",
    "\\begin{align}\n",
    "D(U) = \\sum_\\mu \\gamma_\\mu D_\\mu + m - \\frac12 \\sum_\\mu \\overleftarrow{D}_\\mu \\overrightarrow{D}_\\mu \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "We first note that\n",
    "$$\n",
    "\\begin{align}\n",
    "(C_\\mu^+)^\\dagger = (U_\\mu S_\\mu^+)^\\dagger = (S_\\mu^+)^\\dagger U_\\mu^\\dagger = S_\\mu^- U_\\mu^\\dagger = C_\\mu^-\n",
    "\\end{align}\n",
    "$$\n",
    "since for any $f$ and $g$\n",
    "$$\n",
    "\\begin{align}\n",
    "\\sum_x f(x)^\\dagger (S_\\mu^+)^\\dagger g(x) &=\n",
    "\\sum_x (g(x)^\\dagger S_\\mu^+ f(x))^\\dagger = \\sum_x (g(x)^\\dagger f(x+a\\hat{\\mu}))^\\dagger \n",
    "= \\sum_x f(x+a\\hat{\\mu})^\\dagger g(x) \\\\\n",
    "&= \\sum_x f(x)^\\dagger g(x-a\\hat{\\mu}) = \\sum f(x)^\\dagger S_\\mu^- g(x)\n",
    "\\end{align}\n",
    "$$\n",
    "and therefore $(S_\\mu^+)^\\dagger = S_\\mu^-$.\n",
    "\n",
    "Equipped with these identities, we find that\n",
    "$$\n",
    "\\begin{align}\n",
    "\\gamma_5 D(U)^\\dagger \\gamma_5 &= \\gamma_5\\left(\\sum_\\mu \\gamma_\\mu (D_\\mu)^\\dagger + m - \\frac12 \\sum_\\mu (\\overleftarrow{D}_\\mu \\overrightarrow{D}_\\mu)^\\dagger \\right)\\gamma_5 \\\\\n",
    "&=\\left(-\\sum_\\mu \\gamma_\\mu (D_\\mu)^\\dagger + m - \\frac12 \\sum_\\mu (\\overleftarrow{D}_\\mu \\overrightarrow{D}_\\mu)^\\dagger \\right) \\\\\n",
    "&=\\left(\\sum_\\mu \\gamma_\\mu D_\\mu + m - \\frac12 \\sum_\\mu \\overleftarrow{D}_\\mu \\overrightarrow{D}_\\mu \\right) \\\\\n",
    "&=D(U) \\,.\n",
    "\\end{align}\n",
    "$$\n",
    "This is equivalent to\n",
    "$$\n",
    "\\begin{align}\n",
    "(\\gamma_5 D(U))^\\dagger = \\gamma_5 D(U)\n",
    "\\end{align}\n",
    "$$\n",
    "and\n",
    "$$\n",
    "\\begin{align}\n",
    "(D(U)\\gamma_5 )^\\dagger =  D(U)\\gamma_5 \\,.\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc1f4922",
   "metadata": {},
   "source": [
    "## Space-time symmetries\n",
    "\n",
    "Let us consider an $n$-dimensional lattice $\\Lambda=\\{ \\vec{n} \\in \\mathbb{Z}^n \\}$\n",
    "and operations $g:\\Lambda \\to \\Lambda$ that can be expressed as a sequence of axis exchange operations $A_{\\mu\\nu}$\n",
    "and axis reflection operations $R_\\mu$\n",
    "with\n",
    "$$\n",
    "\\begin{align}\n",
    "(A_{\\mu\\nu} \\vec{n})_\\rho &= \\begin{cases}\n",
    " \\vec{n}_\\rho & \\text{if } \\rho \\neq \\mu \\land \\rho \\neq \\nu \\,, \\\\\n",
    " \\vec{n}_\\mu & \\text{if } \\rho = \\nu \\,, \\\\\n",
    " \\vec{n}_\\nu & \\text{if } \\rho = \\mu \\,,\n",
    "\\end{cases} \\,, \\\\\n",
    "(R_{\\mu} \\vec{n})_\\rho &=\\begin{cases}\n",
    " \\vec{n}_\\rho & \\text{if } \\rho \\neq \\mu\\,, \\\\\n",
    " -\\vec{n}_\\rho & \\text{if } \\rho = \\mu \\,,\n",
    "\\end{cases} \\,, \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "where it is sufficient to consider $A_{\\mu\\nu}$ with $\\mu < \\nu$.\n",
    "\n",
    "Rotations, e.g., can be written as a combination of both.  Take the rotation along the $z$ axis of a three-dimensional vector\n",
    "$$\n",
    "\\begin{align}\n",
    "\\begin{pmatrix}\n",
    "x \\\\ \n",
    "y \\\\\n",
    "z\n",
    "\\end{pmatrix}\n",
    "\\to\n",
    "\\begin{pmatrix}\n",
    "y \\\\ \n",
    "-x \\\\\n",
    "z\n",
    "\\end{pmatrix}\n",
    "\\end{align}\n",
    "$$\n",
    "which can be written as a $R_1 \\circ A_{01}$.\n",
    "\n",
    "In $n$ dimensions we therefore have\n",
    "$$\n",
    "\\begin{align}\n",
    "n! 2^n\n",
    "\\end{align}\n",
    "$$\n",
    "unique non-translation space-time symmetries on a lattice.  If we consider\n",
    "the symmetries of the Hamiltonian, $n=3$ is relevant and we find $48$ unique operations.\n",
    "The corresponding group is called the full cubic group $O_h$ and the study of its\n",
    "irreducible representations yields good quantum numbers of the Hamiltonian.  We will\n",
    "return to this later.\n",
    "\n",
    "For the action $n=4$ is relevant and we have 384 unique operations.  Since all of them can\n",
    "be written as combinations of axis-exchange and reflection operations it is sufficient to\n",
    "study the behavior of the action under this restricted set of operations.  The corresponding\n",
    "group is called hypercubic group.\n",
    "\n",
    "Under exchange of axis, we re-label the Lorentz index of the\n",
    "gluon fields $U_\\mu$ accordingly and apply $A_{\\mu\\nu}$ to the coordinates $x$.\n",
    "The symmetry of the gluon Wilson action under axis exchange symmetry \n",
    "is trivial to see.  Finally, we need to consider\n",
    "transformations of the fermion fields.  We find axis exchange symmetry of the entire action if\n",
    "we transform $\\psi$ and $\\bar\\psi$ such that\n",
    "$$\n",
    "\\begin{align}\n",
    "\\bar\\psi \\gamma_\\rho \\psi = \\bar\\psi \\gamma_{\\rho'} \\psi\n",
    "\\end{align}\n",
    "$$\n",
    "where $\\rho'=\\rho$ if $\\mu \\neq \\rho \\land \\nu \\neq \\rho$ and $\\rho'=\\mu$ if $\\rho=\\nu$ and $\\rho'=\\nu$ if $\\rho=\\mu$.\n",
    "It can be shown that using a change of basis in spinor space for $\\psi$ and $\\bar\\psi$\n",
    "such transformations can always be found.\n",
    "\n",
    "The above symmetries also apply to a lattice with finite number of spatial sites $N_\\mu$ in each direction.\n",
    "For the reflection symmetry it is convenient to consider lattice sites to run from $-N_\\mu/2+1$ to $N_\\mu/2$.\n",
    "In this convention, the action also exhibits a reflection symmetry $R_\\mu$.\n",
    "It is, however, more convenient to study instead the set of reflections of all orthogonal directions to $\\mu$ that we write as\n",
    "$P_\\mu$.  As an example in four dimensions $P_0 = R_1 R_2 R_3$ which can be inverted to $R_0 = P_1 P_2 P_3$.\n",
    "The action can be shown to be invariant under all $P_\\mu$ if\n",
    "$$\n",
    "\\begin{align}\n",
    "\\psi(x) & \\to \\gamma_\\mu \\psi(P_\\mu x) \\,, \\\\\n",
    "\\bar\\psi(x) & \\to \\bar\\psi(P_\\mu x)\\gamma_\\mu  \\,, \\\\\n",
    "U_\\mu(x) &\\to U_\\mu(P_\\mu x) \\,,\\\\\n",
    "U_\\nu(x) &\\to U_\\nu(P_\\mu x - a \\hat{\\nu})^\\dagger\n",
    "\\end{align}\n",
    "$$\n",
    "for $\\nu \\neq \\mu$.  The corresponding calculation is straightforward.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "facd86ba",
   "metadata": {},
   "source": [
    "## Study of the full cubic group $O_h$ - implications for the eigenstates of the Hamiltonian"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "210208fa",
   "metadata": {},
   "source": [
    "In the following, we will study the finite group $O_h$ in more detail in an attempt to understand\n",
    "the resulting good quantum numbers.  We will introduce group-theoretic concepts as needed.\n",
    "\n",
    "As a reminder:  A **group** is the tuple $(G,\\circ)$ with set $G$ and $a\\circ b:G\\times G\\to G$ such that\n",
    "- ($\\forall a,b,c \\in G$)($(a \\circ b) \\circ c = a \\circ (b \\circ c)$)\n",
    "- ($\\exists e \\in G$)($\\forall a \\in G$)($e \\circ a = a \\circ e = a$); $e$ is unique\n",
    "- ($\\forall a \\in G$)($\\exists a^{-1} \\in G$)($a \\circ a^{-1} = a^{-1} \\circ a = e$)\n",
    "\n",
    "A group is called **finite** if $\\vert G \\vert \\in \\mathbb{N}$.\n",
    "\n",
    "For finite groups these studies can be performed numerically in a straightforward manner.  To do so, we first\n",
    "need to construct a set of all group elements that are generated by axis interchange and reflection.  We do so by considering the action on a complete basis of $\\Lambda$ for $n=3$ and demonstrate that the set is closed under the generating operations.\n",
    "\n",
    "Our first task is to associate for each $g_i\\in G$ an index $i\\in [0,\\ldots,\\vert G \\vert - 1]$\n",
    "and to create a **multiplication table** $M \\in [0,\\ldots,\\vert G \\vert - 1]^{\\vert G \\vert \\times \\vert G \\vert}$\n",
    "with elements\n",
    "$$\n",
    " M_{ij} = l   \\quad \\leftrightarrow \\quad g_j \\circ g_i = g_l \\,.\n",
    "$$\n",
    "\n",
    "To simplify the discussion, we first study the group $O$ that is made up of all rotations before studying the full cubic group $O_h$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e31371b0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "24\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "basis = [[1,0,0],[0,1,0],[0,0,1]]\n",
    "\n",
    "def interchange(x,a,b):\n",
    "    xp=[xi for xi in x]\n",
    "    xp[a]=x[b]\n",
    "    xp[b]=x[a]\n",
    "    return xp\n",
    "\n",
    "def reflection(x,a):\n",
    "    xp=[xi for xi in x]\n",
    "    xp[a]=-x[a]\n",
    "    return xp\n",
    "\n",
    "def rotateZ90(x):\n",
    "    return [-x[1],x[0],x[2]]\n",
    "\n",
    "def rotateX90(x):\n",
    "    return [x[0],-x[2],x[1]]\n",
    "\n",
    "def find_closure(actions):\n",
    "    elements = [basis]\n",
    "    sequence = { str(basis) : [] }\n",
    "    n = 0\n",
    "    while n != len(elements):\n",
    "        n = len(elements)\n",
    "        for a in actions:\n",
    "            for e in elements:\n",
    "                ep = a(e)\n",
    "                if ep not in elements:\n",
    "                    elements.append(ep)\n",
    "                    sequence[str(ep)] = sequence[str(e)] + [a]\n",
    "\n",
    "    return elements, sequence\n",
    "    \n",
    "# First only study group of rotations, i.e., cubic group O\n",
    "elements = find_closure([\n",
    "    #lambda e: [interchange(b,0,1) for b in e],\n",
    "    #lambda e: [interchange(b,0,2) for b in e],\n",
    "    #lambda e: [interchange(b,1,2) for b in e],\n",
    "    #lambda e: [reflection(b,0) for b in e],\n",
    "    #lambda e: [reflection(b,1) for b in e],\n",
    "    #lambda e: [reflection(b,2) for b in e],\n",
    "    lambda e: [rotateZ90(b) for b in e],\n",
    "    lambda e: [rotateX90(b) for b in e]\n",
    "])\n",
    "\n",
    "print(len(elements[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f7cbe0b2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23]\n",
      " [ 1  2  3  0  5  6  7  4  9 10 11  8 13 14 15 12 17 18 19 16 21 22 23 20]\n",
      " [ 2  3  0  1  6  7  4  5 10 11  8  9 14 15 12 13 18 19 16 17 22 23 20 21]\n",
      " [ 3  0  1  2  7  4  5  6 11  8  9 10 15 12 13 14 19 16 17 18 23 20 21 22]\n",
      " [ 4 16 14 22  8 17  2 21 12 18  6 20  0 19 10 23 11 15  3  7  1 13  9  5]\n",
      " [ 5 17 15 23  9 18  3 22 13 19  7 21  1 16 11 20  8 12  0  4  2 14 10  6]\n",
      " [ 6 18 12 20 10 19  0 23 14 16  4 22  2 17  8 21  9 13  1  5  3 15 11  7]\n",
      " [ 7 19 13 21 11 16  1 20 15 17  5 23  3 18  9 22 10 14  2  6  0 12  8  4]\n",
      " [ 8 11 10  9 12 15 14 13  0  3  2  1  4  7  6  5 20 23 22 21 16 19 18 17]\n",
      " [ 9  8 11 10 13 12 15 14  1  0  3  2  5  4  7  6 21 20 23 22 17 16 19 18]\n",
      " [10  9  8 11 14 13 12 15  2  1  0  3  6  5  4  7 22 21 20 23 18 17 16 19]\n",
      " [11 10  9  8 15 14 13 12  3  2  1  0  7  6  5  4 23 22 21 20 19 18 17 16]\n",
      " [12 20  6 18  0 23 10 19  4 22 14 16  8 21  2 17  1  5  9 13 11  7  3 15]\n",
      " [13 21  7 19  1 20 11 16  5 23 15 17  9 22  3 18  2  6 10 14  8  4  0 12]\n",
      " [14 22  4 16  2 21  8 17  6 20 12 18 10 23  0 19  3  7 11 15  9  5  1 13]\n",
      " [15 23  5 17  3 22  9 18  7 21 13 19 11 20  1 16  0  4  8 12 10  6  2 14]\n",
      " [16 14 22  4 17  2 21  8 18  6 20 12 19 10 23  0 15  3  7 11 13  9  5  1]\n",
      " [17 15 23  5 18  3 22  9 19  7 21 13 16 11 20  1 12  0  4  8 14 10  6  2]\n",
      " [18 12 20  6 19  0 23 10 16  4 22 14 17  8 21  2 13  1  5  9 15 11  7  3]\n",
      " [19 13 21  7 16  1 20 11 17  5 23 15 18  9 22  3 14  2  6 10 12  8  4  0]\n",
      " [20  6 18 12 23 10 19  0 22 14 16  4 21  2 17  8  5  9 13  1  7  3 15 11]\n",
      " [21  7 19 13 20 11 16  1 23 15 17  5 22  3 18  9  6 10 14  2  4  0 12  8]\n",
      " [22  4 16 14 21  8 17  2 20 12 18  6 23  0 19 10  7 11 15  3  5  1 13  9]\n",
      " [23  5 17 15 22  9 18  3 21 13 19  7 20  1 16 11  4  8 12  0  6  2 14 10]]\n"
     ]
    }
   ],
   "source": [
    "def apply(el, e0, idx):\n",
    "    for acc in el[1][str(el[0][idx])]:\n",
    "        e0 = acc(e0)\n",
    "    return e0\n",
    "\n",
    "def multiplication_table(el, ba):\n",
    "    n=len(el[0])\n",
    "    M=[[-1 for i in range(n)] for j in range(n)]\n",
    "    for i in range(n):\n",
    "        ei = apply(el, ba, i)\n",
    "        for j in range(n):\n",
    "            eji = apply(el, ei, j)\n",
    "            M[i][j] = el[0].index(eji)\n",
    "    return M\n",
    "\n",
    "    \n",
    "mt = multiplication_table(elements, basis)\n",
    "\n",
    "print(np.array(mt))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "4e37aee4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "[0, 3, 2, 1, 12, 18, 6, 20, 8, 9, 10, 11, 4, 22, 14, 16, 15, 17, 5, 23, 7, 21, 13, 19]\n"
     ]
    }
   ],
   "source": [
    "def find_identity(el, ba):\n",
    "    return el[0].index(ba)\n",
    "\n",
    "ident = find_identity(elements, basis)\n",
    "assert ident == 0 # for convenience always have first element be the identity\n",
    "\n",
    "print(ident)\n",
    "\n",
    "def inverse_list(mt, e):\n",
    "    return [x.index(e) for x in mt]\n",
    "\n",
    "inv = inverse_list(mt, ident)\n",
    "\n",
    "print(inv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "3468aa48",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]\n",
      "[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]\n"
     ]
    }
   ],
   "source": [
    "# test inverses\n",
    "print([mt[i][inv[i]] for i in range(len(inv))])\n",
    "print([mt[inv[i]][i] for i in range(len(inv))])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d15092be",
   "metadata": {},
   "source": [
    "It is convenient to identify the index 0 with the identity element $e$.\n",
    "\n",
    "Next, we consider $\\Gamma_r : G \\to \\text{GL}(V)$ that maps group\n",
    "elements to invertible matrices acting on a finite-dimensional vector space $V$.\n",
    "We call $\\Gamma_r$ a **representation** of the group if it preserves its structure, i.e.,\n",
    "$$\n",
    "(\\forall a,b \\in G)(\\Gamma_r(a) \\Gamma_r(b) = \\Gamma_r(a \\circ b)) \\,.\n",
    "$$\n",
    "In other words, it is a group homomorphism between $G$ and the general linear group over vector space $V$.\n",
    "\n",
    "If there is a lower-dimensional subspace $V'$ of $V$ for which ($\\forall g \\in G, v \\in V'$)($\\Gamma_r(g) v \\in V'$)\n",
    "we call the representation $\\Gamma_r$ **reducible** and the projection of $\\Gamma_r$ to $V'$ a **sub-representation**.\n",
    "In other words $\\Gamma_r$ has a block-matrix form in the appropriate basis.  If no such sub-representation exists,\n",
    "we call $\\Gamma_r$ **irreducible**.\n",
    "\n",
    "If $G$ describes a symmetry of a Hamiltonian $H$, the Hamiltonian has a corresponding block structure and\n",
    "the spectrum of $H$ can be classified following the irreducible representations.  This will be one of the main\n",
    "uses of representation theory in the remainder of this lecture.\n",
    "\n",
    "So the next question is, given our group $G$, how many irreducible representations are there?  For a finite group\n",
    "there will always be a finite number and there is a straightforward method to find out how many there are.\n",
    "\n",
    "To this end, we first define a **conjugacy class** $C$ as follows:  For all $a,b \\in G$\n",
    "we say $a$ and $b$ are **conjugate** ($a \\sim b$) if there is\n",
    "a group element $g\\in G$ such that $g a g^{-1} = b$.  This defines an equivalence relation, i.e.,\n",
    "if $a \\sim b$ and $b\\sim c$ then $a\\sim c$.  Therefore it splits the group $G$ in separate\n",
    "conjugacy classes, where a member of each class is conjguate to all others but not conjugate to members\n",
    "of a different class.\n",
    "\n",
    "It is straightforward to construct the conjugacy classes of a finite group, as we will do below, and one can show that the number of conjugacy classes is exactly the number of irreducible representations (irreps) of the group!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5db47160",
   "metadata": {},
   "outputs": [],
   "source": [
    "# find conjugation classes / invariant subspaces\n",
    "def find_conjugation_classes(mt, inv, ident):\n",
    "    n=len(inv)\n",
    "    tostudy=list(range(n))\n",
    "    classes = []\n",
    "    while len(tostudy) > 0:\n",
    "        x = tostudy[0]\n",
    "        conjugate_to = []\n",
    "        for i in range(n):\n",
    "            conjugate_to.append(mt[mt[i][x]][inv[i]])\n",
    "        conjugate_to = set(conjugate_to)\n",
    "        for y in conjugate_to:\n",
    "            tostudy.remove(y)\n",
    "        classes.append(conjugate_to)\n",
    "    return classes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f6c8941a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[{0}, {1, 3, 4, 12, 19, 23}, {8, 2, 10}, {5, 7, 13, 15, 16, 18, 20, 22}, {6, 9, 11, 14, 17, 21}]\n"
     ]
    }
   ],
   "source": [
    "classes = find_conjugation_classes(mt, inv, ident)\n",
    "print(classes)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97961e63",
   "metadata": {},
   "source": [
    "So there will be five irreps for $O$.  Next, we define a few obvious representations and test them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "298918b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# find a matrix representation of individual elements\n",
    "def matrix_vector(l):\n",
    "    # transform a lattice vector\n",
    "    return np.matrix(elements[0][l]).T\n",
    "\n",
    "def matrix_trivial(l):\n",
    "    return np.matrix([[1]])\n",
    "\n",
    "def matrix_determinant(l):\n",
    "    return np.matrix([[np.linalg.det(matrix_vector(l))]], dtype=np.int32)\n",
    "\n",
    "def matrix_regular(l):\n",
    "    M = np.zeros(shape=(len(mt),len(mt)), dtype=np.int32)\n",
    "    for i in range(len(mt)):\n",
    "        M[i,mt[l][i]] = 1\n",
    "    return np.matrix(M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "4f9627c4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "matrix_vector\n",
      "matrix_regular\n",
      "matrix_trivial\n"
     ]
    }
   ],
   "source": [
    "# test if matrix is a group representation\n",
    "def test_representation(matrix):\n",
    "    print(matrix.__name__)\n",
    "    for i in range(len(mt)):\n",
    "        for j in range(len(mt)):\n",
    "            eps = np.linalg.norm(matrix(j) * matrix(i) - matrix(mt[i][j]))\n",
    "            assert eps < 1e-13\n",
    "            \n",
    "test_representation(matrix_vector)\n",
    "test_representation(matrix_regular)\n",
    "test_representation(matrix_trivial)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9dc35f9",
   "metadata": {},
   "source": [
    "Next, we need a test if a representation is irreducible.  To this end,\n",
    "we first define the **character** of a representation $r$\n",
    "$$\n",
    "\\chi_r(g) = \\text{Tr}\\:\\Gamma_r(g) \\,.\n",
    "$$\n",
    "It is straightforward to see that the character of two elements $a,b\\in G$\n",
    "is identical if they belong to the same conjugacy class.  We call such functions\n",
    "**class functions**.\n",
    "\n",
    "Let $C$ be the set of all conjugacy classes of group $G$, then the matrix\n",
    "$$\n",
    "(\\chi_r)_c\n",
    "$$\n",
    "for $c \\in C$ and irreps $r$ is the so-called **character table** of the group.\n",
    "One writes $c$ in the columns and $r$ in the rows.\n",
    "\n",
    "Next, we\n",
    "define the inner product of two class functions $\\chi: G \\to \\mathbb{R}$ and $\\phi: G\\to \\mathbb{R}$\n",
    "as\n",
    "$$\n",
    "\\langle \\chi, \\phi \\rangle = \\frac1{\\vert G \\vert} \\sum_{g \\in G} \\chi(g) \\phi(g) \\,.\n",
    "$$\n",
    "\n",
    "Finally, one can show that $\\langle \\chi_r , \\chi_r \\rangle = 1$ if and only if $r$ is irreducible.\n",
    "\n",
    "We therefore now know how many irreducible representations there are and how to check if a representation is irreducible.  We next check the representations that we already constructed above and then figure out a way to construct the remaining irreps.\n",
    "\n",
    "Before doing so, we note that from Schur's orthogonality theorem it follows that two rows in the\n",
    "character table are orthonormal with respect to the above inner product.  We also note that\n",
    "$\\chi_r(e)$ is the dimensionality of the irreducible representation.  By convention one writes the conjugacy\n",
    "class containing $e$ in the first column of the character table and the trivial representation that\n",
    "maps everything to $1$ to the first row."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "4a46a2ef",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n",
      "1.0\n",
      "1.0\n",
      "|C|   1,   6,   3,   8,   6\n",
      "------------------------------\n",
      "A1    1,   1,   1,   1,   1\n",
      "T1    3,   1,  -1,   0,  -1\n"
     ]
    }
   ],
   "source": [
    "def fmt(row):\n",
    "    return \", \".join([ \"%3.0f\" % i for i in row ])\n",
    "\n",
    "def character_table_row(irrep):\n",
    "    return [np.trace(irrep(list(c)[0])) for c in classes]\n",
    "\n",
    "def inner_product(a, b):\n",
    "    return sum([ai*bi*len(ci) for ai, bi, ci in zip(a, b, classes)]) / len(mt)\n",
    "\n",
    "characters_trivial = character_table_row(matrix_trivial)\n",
    "characters_vector = character_table_row(matrix_vector)\n",
    "\n",
    "print(inner_product(characters_trivial, characters_vector))\n",
    "print(inner_product(characters_trivial, characters_trivial)) # == 1 -> irreducible\n",
    "print(inner_product(characters_vector, characters_vector)) # == 1 -> irreducible\n",
    "\n",
    "print(\"|C|\",fmt([len(c) for c in classes]))\n",
    "print(\"-\" * 30)\n",
    "print(\"A1 \",fmt(characters_trivial))\n",
    "print(\"T1 \",fmt(characters_vector))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "9a4abca6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "24.0\n"
     ]
    }
   ],
   "source": [
    "characters_regular = character_table_row(matrix_regular)\n",
    "print(inner_product(characters_regular, characters_regular)) # not irreducible"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22476fe1",
   "metadata": {},
   "source": [
    "The names A1 for the trivial representation and T1 for the representation transforming as a vector\n",
    "are convention but the rows in the character table are the unique fingerprint that defines them.\n",
    "\n",
    "In order to find all irreps, we can study tensor products of existing representations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ae43acb7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "matrix_vector_X_matrix_vector\n",
      "4.0\n"
     ]
    }
   ],
   "source": [
    "# product representation\n",
    "def product_matrix(A, B):\n",
    "    def _mat(i):\n",
    "        return np.kron(A(i),B(i))\n",
    "    _mat.__name__ = A.__name__ + \"_X_\" + B.__name__\n",
    "    return _mat\n",
    "\n",
    "matrix_vector_vector = product_matrix(matrix_vector, matrix_vector)\n",
    "test_representation(matrix_vector_vector)\n",
    "characters_vector_vector = character_table_row(matrix_vector_vector)\n",
    "print(inner_product(characters_vector_vector, characters_vector_vector)) # not irreducible\n",
    "\n",
    "#tA=[[1,2],[3,4]]\n",
    "#tB=[[1.4,2.2],[1.2,3.9]]\n",
    "#print(np.kron(np.matrix(tA),np.matrix(tB))-np.matrix([[tA[i1][j1]*tB[i2][j2] for j1 in range(2) for j2 in range(2)] for i1 in range(2) for i2 in range(2)]))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "041cde12",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n",
      "1.0\n"
     ]
    }
   ],
   "source": [
    "print(inner_product(characters_vector_vector, characters_vector)) # contains both T1 and A1\n",
    "print(inner_product(characters_vector_vector, characters_trivial))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8620add9",
   "metadata": {},
   "source": [
    "In order to find irreducible representations again, we need to study projections of the product space down to invariant subspaces.  One can show that given a representation $\\Gamma_r:G \\to \\text{GL}(V)$\n",
    "and a character $\\chi_s$ of an irrep $s$,\n",
    "$$\n",
    "\\begin{align}\n",
    "P(\\chi_s) = \\frac{\\text{dim}(V)}{\\vert G \\vert} \\sum_{g \\in G} \\Gamma_r(g) \\chi_s(g)\n",
    "\\end{align}\n",
    "$$\n",
    "projects to the corresponding invariant subspace $V'$ of the irreducible representation.\n",
    "Therefore\n",
    "$$\n",
    "\\begin{align}\n",
    "P(\\chi_s)\\Gamma_r P(\\chi_s)\n",
    "\\end{align}\n",
    "$$\n",
    "is either zero or the irreducible representation $s$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "a30afd5a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 0. 0.]\n",
      " [0. 1. 0.]\n",
      " [0. 0. 1.]]\n",
      "[[0. 0. 0.]\n",
      " [0. 0. 0.]\n",
      " [0. 0. 0.]]\n",
      "[[0.]]\n",
      "[[1.]]\n"
     ]
    }
   ],
   "source": [
    "def projection(matrix, character_row):\n",
    "    ret = 0.0 * matrix(ident)\n",
    "    for chi, c in zip(character_row, classes):\n",
    "        for ci in c:\n",
    "            ret += matrix(ci)*chi\n",
    "    P = np.matrix(ret * character_row[0] / len(mt))\n",
    "    assert np.linalg.norm(P*P - P) < 1e-14\n",
    "    return P\n",
    "        \n",
    "print(projection(matrix_vector, characters_vector))\n",
    "print(projection(matrix_vector, characters_trivial))\n",
    "\n",
    "print(projection(matrix_trivial, characters_vector))\n",
    "print(projection(matrix_trivial, characters_trivial))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "55502422",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.   0.   0.   0.   0.   0.   0.   0.   0. ]\n",
      " [ 0.   0.5  0.  -0.5  0.   0.   0.   0.   0. ]\n",
      " [ 0.   0.   0.5  0.   0.   0.  -0.5  0.   0. ]\n",
      " [ 0.  -0.5  0.   0.5  0.   0.   0.   0.   0. ]\n",
      " [ 0.   0.   0.   0.   0.   0.   0.   0.   0. ]\n",
      " [ 0.   0.   0.   0.   0.   0.5  0.  -0.5  0. ]\n",
      " [ 0.   0.  -0.5  0.   0.   0.   0.5  0.   0. ]\n",
      " [ 0.   0.   0.   0.   0.  -0.5  0.   0.5  0. ]\n",
      " [ 0.   0.   0.   0.   0.   0.   0.   0.   0. ]]\n",
      "[[0.33333333 0.         0.         0.         0.33333333 0.\n",
      "  0.         0.         0.33333333]\n",
      " [0.         0.         0.         0.         0.         0.\n",
      "  0.         0.         0.        ]\n",
      " [0.         0.         0.         0.         0.         0.\n",
      "  0.         0.         0.        ]\n",
      " [0.         0.         0.         0.         0.         0.\n",
      "  0.         0.         0.        ]\n",
      " [0.33333333 0.         0.         0.         0.33333333 0.\n",
      "  0.         0.         0.33333333]\n",
      " [0.         0.         0.         0.         0.         0.\n",
      "  0.         0.         0.        ]\n",
      " [0.         0.         0.         0.         0.         0.\n",
      "  0.         0.         0.        ]\n",
      " [0.         0.         0.         0.         0.         0.\n",
      "  0.         0.         0.        ]\n",
      " [0.33333333 0.         0.         0.         0.33333333 0.\n",
      "  0.         0.         0.33333333]]\n"
     ]
    }
   ],
   "source": [
    "print(projection(matrix_vector_vector, characters_vector))\n",
    "print(projection(matrix_vector_vector, characters_trivial))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ada20739",
   "metadata": {},
   "source": [
    "It is straightforward to show that properties such as symmetrization/antisymmetrization of basis vectors in the product space create invariant subspaces.  These subspaces may or may not correspond to irreducible representations, however, they already reduce the dimensionality of the relevant space and therefore are powerful tools to identify the remaining irreps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "812ca49a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def antisymmetric_projector(n):\n",
    "    ei = [np.zeros(shape=(n,)) for i in range(n)]\n",
    "    for i in range(n):\n",
    "        ei[i][i] = 1\n",
    "        \n",
    "    # 00 01 10 11\n",
    "    # (|01> - |10>)(<01| - <10|) / 2\n",
    "    # v_i v_j -> v_i v_j - v_j v_i\n",
    "    P = sum([\n",
    "        np.matrix(np.outer(np.kron(ei[i],ei[j]) - np.kron(ei[j],ei[i]),\n",
    "                           np.kron(ei[i],ei[j]) - np.kron(ei[j],ei[i])))/2\n",
    "     for i in range(n) for j in range(i)])\n",
    "        \n",
    "    assert np.linalg.norm(P*P - P) < 1e-14\n",
    "    return P\n",
    "\n",
    "def symmetric_projector(n):\n",
    "    ei = [np.zeros(shape=(n,)) for i in range(n)]\n",
    "    for i in range(n):\n",
    "        ei[i][i] = 1\n",
    "        \n",
    "    # 00 01 10 11\n",
    "    # (|01> + |10>)(<01| + <10|) / 2 + |00><00| + |11><11|\n",
    "    # v_i v_j -> v_i v_j + v_j v_i\n",
    "    P = sum([\n",
    "        np.matrix(np.outer(np.kron(ei[i],ei[j]) + np.kron(ei[j],ei[i]),\n",
    "                           np.kron(ei[i],ei[j]) + np.kron(ei[j],ei[i])))/4.\n",
    "     for i in range(n) for j in range(n)])\n",
    "        \n",
    "    assert np.linalg.norm(P*P - P) < 1e-14\n",
    "    return P\n",
    "\n",
    "def diagonal_projector(n):\n",
    "    ei = [np.zeros(shape=(n,)) for i in range(n)]\n",
    "    for i in range(n):\n",
    "        ei[i][i] = 1\n",
    "        \n",
    "    # |00><00| + |11><11|\n",
    "    # v_i v_j -> v_i v_j \\delta_{ij}\n",
    "    P = sum([\n",
    "        np.matrix(np.outer(np.kron(ei[i],ei[i]),\n",
    "                           np.kron(ei[i],ei[i])))\n",
    "     for i in range(n)])\n",
    "        \n",
    "    assert np.linalg.norm(P*P - P) < 1e-14\n",
    "    return P\n",
    "\n",
    "def invert_projector(projector):\n",
    "    def _proj(n):\n",
    "        P = projector(n)\n",
    "        P = np.matrix(np.eye(n*n)) - P\n",
    "        assert np.linalg.norm(P*P - P) < 1e-14\n",
    "        return P\n",
    "    return _proj\n",
    "        \n",
    "def traceless_projector(n):\n",
    "    ei = [np.zeros(shape=(n,)) for i in range(n)]\n",
    "    for i in range(n):\n",
    "        ei[i][i] = 1\n",
    "        \n",
    "    # (|ii>-\\sum_l |ll> / ndim)(<ii|-\\sum_l <ll| / ndim) + sum_{i!=j} |ij><ij|\n",
    "    # v_i v_j -> v_i v_j - v_l v_l \\delta_{ij} / ndim\n",
    "    P = sum([\n",
    "        np.matrix(np.outer(np.kron(ei[i],ei[j]),np.kron(ei[i],ei[j])))\n",
    "     for i in range(n) for j in range(n) if i != j])\n",
    "    \n",
    "    sll = sum([np.kron(ei[l],ei[l]) for l in range(n)]) / n\n",
    "    \n",
    "    P += sum([\n",
    "        np.matrix(np.outer(np.kron(ei[i],ei[i]) - sll,np.kron(ei[i],ei[i]) - sll))\n",
    "     for i in range(n)])\n",
    "    \n",
    "    assert np.linalg.norm(P*P - P) < 1e-14\n",
    "    return P\n",
    "\n",
    "# product representation\n",
    "def projected_product_matrix(A, B, projectors):\n",
    "    \n",
    "    n1 = len(A(0))\n",
    "    n2 = len(B(0))\n",
    "    assert n1 == n2\n",
    "    \n",
    "    P = projectors[0](n1)\n",
    "    for p in projectors[1:]:\n",
    "        P = P * p(n1)\n",
    "             \n",
    "    def _mat(i):\n",
    "        return P*np.kron(A(i),B(i))*P\n",
    "    _mat.__name__ = A.__name__ + \"_projected_\" + B.__name__\n",
    "    return _mat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "b2843edc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "matrix_vector_projected_matrix_vector\n",
      "1.0\n",
      "|C|   1,   6,   3,   8,   6\n",
      "------------------------------\n",
      "A1    1,   1,   1,   1,   1\n",
      "T1    3,   1,  -1,   0,  -1\n",
      "AVV   3,   1,  -1,   0,  -1\n"
     ]
    }
   ],
   "source": [
    "matrix_a_vector_vector = projected_product_matrix(matrix_vector, matrix_vector, \n",
    "                                                         [antisymmetric_projector])\n",
    "test_representation(matrix_a_vector_vector)\n",
    "characters_a_vector_vector = character_table_row(matrix_a_vector_vector)\n",
    "print(inner_product(characters_a_vector_vector, characters_a_vector_vector)) # is irreducible\n",
    "\n",
    "\n",
    "print(\"|C|\",fmt([len(c) for c in classes]))\n",
    "print(\"-\" * 30)\n",
    "print(\"A1 \",fmt(characters_trivial))\n",
    "print(\"T1 \",fmt(characters_vector))\n",
    "print(\"AVV\",fmt(characters_a_vector_vector)) # but unfortunately already known, it is the T1 representation!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "bcedd29b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "matrix_vector_projected_matrix_vector\n",
      "3.0\n"
     ]
    }
   ],
   "source": [
    "matrix_s_vector_vector = projected_product_matrix(matrix_vector, matrix_vector, \n",
    "                                                         [symmetric_projector])\n",
    "test_representation(matrix_s_vector_vector)\n",
    "characters_s_vector_vector = character_table_row(matrix_s_vector_vector)\n",
    "print(inner_product(characters_s_vector_vector, characters_s_vector_vector)) # reducible"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "45de16c2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "matrix_vector_projected_matrix_vector\n",
      "1.0\n",
      "|C|   1,   6,   3,   8,   6\n",
      "------------------------------\n",
      "A1    1,   1,   1,   1,   1\n",
      "T1    3,   1,  -1,   0,  -1\n",
      "T2    3,  -1,  -1,   0,   1\n"
     ]
    }
   ],
   "source": [
    "# only look at off-diagonal part\n",
    "matrix_s_offdiag_vector_vector = projected_product_matrix(matrix_vector, matrix_vector, \n",
    "                                                         [invert_projector(diagonal_projector),symmetric_projector])\n",
    "test_representation(matrix_s_offdiag_vector_vector)\n",
    "characters_s_offdiag_vector_vector = character_table_row(matrix_s_offdiag_vector_vector)\n",
    "print(inner_product(characters_s_offdiag_vector_vector, characters_s_offdiag_vector_vector)) # irreducible\n",
    "\n",
    "print(\"|C|\",fmt([len(c) for c in classes]))\n",
    "print(\"-\" * 30)\n",
    "print(\"A1 \",fmt(characters_trivial))\n",
    "print(\"T1 \",fmt(characters_vector))\n",
    "print(\"T2 \",fmt(characters_s_offdiag_vector_vector)) # new irrep!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "d279d894",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "matrix_vector_projected_matrix_vector\n",
      "2.0\n"
     ]
    }
   ],
   "source": [
    "# now look at diagonal part\n",
    "matrix_diag_vector_vector = projected_product_matrix(matrix_vector, matrix_vector, \n",
    "                                                         [diagonal_projector])\n",
    "test_representation(matrix_diag_vector_vector)\n",
    "characters_diag_vector_vector = character_table_row(matrix_diag_vector_vector)\n",
    "print(inner_product(characters_diag_vector_vector, characters_diag_vector_vector)) # reducible"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "afe156c9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "matrix_vector_projected_matrix_vector\n",
      "1.0000000000000004\n",
      "|C|   1,   6,   3,   8,   6\n",
      "------------------------------\n",
      "A1    1,   1,   1,   1,   1\n",
      "T1    3,   1,  -1,   0,  -1\n",
      "T2    3,  -1,  -1,   0,   1\n",
      "E     2,   0,   2,  -1,   0\n"
     ]
    }
   ],
   "source": [
    "# only the traceless diagonal part\n",
    "matrix_diag_traceless_vector_vector = projected_product_matrix(matrix_vector, matrix_vector, \n",
    "                                                         [diagonal_projector,traceless_projector])\n",
    "test_representation(matrix_diag_traceless_vector_vector)\n",
    "characters_diag_traceless_vector_vector = character_table_row(matrix_diag_traceless_vector_vector)\n",
    "print(inner_product(characters_diag_traceless_vector_vector, characters_diag_traceless_vector_vector)) # irreducible\n",
    "\n",
    "print(\"|C|\",fmt([len(c) for c in classes]))\n",
    "print(\"-\" * 30)\n",
    "print(\"A1 \",fmt(characters_trivial))\n",
    "print(\"T1 \",fmt(characters_vector))\n",
    "print(\"T2 \",fmt(characters_s_offdiag_vector_vector))\n",
    "print(\"E  \",fmt(characters_diag_traceless_vector_vector)) # new irrep!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "799077b1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "matrix_vector_projected_matrix_vector\n",
      "1.0\n",
      "|C|   1,   6,   3,   8,   6\n",
      "------------------------------\n",
      "A1    1,   1,   1,   1,   1\n",
      "T1    3,   1,  -1,   0,  -1\n",
      "T2    3,  -1,  -1,   0,   1\n",
      "E     2,   0,   2,  -1,   0\n",
      "??    1,   1,   1,   1,   1\n"
     ]
    }
   ],
   "source": [
    "# what was the trace diagonal part?  sounds like trivial representaton\n",
    "matrix_diag_trace_vector_vector = projected_product_matrix(matrix_vector, matrix_vector, \n",
    "                                                         [diagonal_projector,invert_projector(traceless_projector)])\n",
    "test_representation(matrix_diag_trace_vector_vector)\n",
    "characters_diag_trace_vector_vector = character_table_row(matrix_diag_trace_vector_vector)\n",
    "print(inner_product(characters_diag_trace_vector_vector, characters_diag_trace_vector_vector)) # irreducible\n",
    "\n",
    "print(\"|C|\",fmt([len(c) for c in classes]))\n",
    "print(\"-\" * 30)\n",
    "print(\"A1 \",fmt(characters_trivial))\n",
    "print(\"T1 \",fmt(characters_vector))\n",
    "print(\"T2 \",fmt(characters_s_offdiag_vector_vector))\n",
    "print(\"E  \",fmt(characters_diag_traceless_vector_vector))\n",
    "print(\"?? \",fmt(characters_diag_trace_vector_vector)) # A1 again"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "5865b48b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "matrix_vector_X_matrix_vector_X_matrix_vector_projected\n",
      "1.0000000000000004\n",
      "|C|   1,   6,   3,   8,   6\n",
      "------------------------------\n",
      "A1    1,   1,   1,   1,   1\n",
      "T1    3,   1,  -1,   0,  -1\n",
      "T2    3,  -1,  -1,   0,   1\n",
      "E     2,   0,   2,  -1,   0\n",
      "A2    1,  -1,   1,   1,  -1\n"
     ]
    }
   ],
   "source": [
    "# we still need one more, take direct product of two irreps and subtract all known irreps\n",
    "# should not matter unless no new irrep is in product (T2 \\otimes T2 does, e.g., not work)\n",
    "def orthogonal_matrix(A, character_rows):\n",
    "    P = A(ident) * 0.0\n",
    "    for cr in character_rows:\n",
    "        P += np.matrix(projection(A, cr))\n",
    "        \n",
    "    P = np.matrix(np.eye(len(P))) - P\n",
    "    assert np.linalg.norm(P*P - P) < 1e-14\n",
    "    \n",
    "    def _mat(i):\n",
    "        return P*A(i)*P\n",
    "    _mat.__name__ = A.__name__ + \"_projected\"\n",
    "    return _mat\n",
    "\n",
    "matrix_last_irrep = orthogonal_matrix(\n",
    "    #product_matrix(matrix_diag_traceless_vector_vector,matrix_diag_traceless_vector_vector),\n",
    "    #product_matrix(matrix_s_offdiag_vector_vector, matrix_vector),\n",
    "    product_matrix(matrix_vector,product_matrix(matrix_vector,matrix_vector)),\n",
    "[\n",
    "    characters_trivial,\n",
    "    characters_vector,\n",
    "    characters_s_offdiag_vector_vector,\n",
    "    characters_diag_traceless_vector_vector\n",
    "])\n",
    "\n",
    "test_representation(matrix_last_irrep)\n",
    "characters_last_irrep = character_table_row(matrix_last_irrep)\n",
    "print(inner_product(characters_last_irrep, characters_last_irrep)) # irreducible!\n",
    "\n",
    "print(\"|C|\",fmt([len(c) for c in classes]))\n",
    "print(\"-\" * 30)\n",
    "print(\"A1 \",fmt(characters_trivial))\n",
    "print(\"T1 \",fmt(characters_vector))\n",
    "print(\"T2 \",fmt(characters_s_offdiag_vector_vector))\n",
    "print(\"E  \",fmt(characters_diag_traceless_vector_vector))\n",
    "print(\"A2 \",fmt(characters_last_irrep)) # looks like an irrep!  we are done"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73b56291",
   "metadata": {},
   "source": [
    "Now we have the complete character table of $O$!  This allows us to project the irreps out\n",
    "of any representation of $G$.  This will be useful, when we construct operators to study the\n",
    "spectrum of $H$.\n",
    "\n",
    "We have shown above, e.g., that $V_i = \\bar{\\psi}\\gamma_i \\psi$\n",
    "transforms as a vector under $O$, i.e., they transform in the T1 representation.\n",
    "\n",
    "We also have already seen in the last chapter that a bilinear $V_i$ creates\n",
    "a meson-type state (made out of one quark and one anti-quark).  If we now\n",
    "want to systematically study the spectrum of two suchs mesons, the above discussion\n",
    "would tell us that the operators\n",
    "$$\n",
    "\\begin{align}\n",
    "V^{(2),\\text{T2}}_{ij} = \\frac12 (V_i V_j + V_j V_i) (1 - \\delta_{ij})\n",
    "\\end{align}\n",
    "$$\n",
    "transform in T2.  Note that $(i,j)=(x,y),(x,z),(y,z)$ are all independent components\n",
    "of the three-dimensional T2 representation (check first column of the character table).\n",
    "Similarly we can create operators that transform in all\n",
    "other irreps and therefore study the complete spectrum of two-meson states."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07a9477e",
   "metadata": {},
   "source": [
    " Finally, we now include the reflections and complete this chapter by studying the full cubic group $O_h$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "2425f878",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "48\n",
      "[{0}, {32, 1, 2, 4, 40, 19}, {33, 3, 35, 5, 39, 41, 21, 23}, {24, 12, 6}, {7, 8, 13, 16, 26, 28}, {9, 11, 45, 15, 47, 17, 27, 29}, {10, 43, 44, 46, 14, 25}, {18, 36, 30}, {34, 37, 38, 20, 22, 31}, {42}]\n"
     ]
    }
   ],
   "source": [
    "elements = find_closure([\n",
    "    lambda e: [interchange(b,0,1) for b in e],\n",
    "    lambda e: [interchange(b,0,2) for b in e],\n",
    "    lambda e: [interchange(b,1,2) for b in e],\n",
    "    lambda e: [reflection(b,0) for b in e],\n",
    "    lambda e: [reflection(b,1) for b in e],\n",
    "    lambda e: [reflection(b,2) for b in e]\n",
    "])\n",
    "\n",
    "mt = multiplication_table(elements, basis)\n",
    "ident = find_identity(elements, basis)\n",
    "assert ident == 0 # for convenience always have first element be the identity\n",
    "inv = inverse_list(mt, ident)\n",
    "\n",
    "print(len(elements[0]))\n",
    "classes = find_conjugation_classes(mt, inv, ident)\n",
    "print(classes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "633863fe",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n",
      "1.0\n",
      "1.0\n",
      "|C|   1,   6,   8,   3,   6,   8,   6,   3,   6,   1\n",
      "------------------------------------------------------------\n",
      "A1+   1,   1,   1,   1,   1,   1,   1,   1,   1,   1\n",
      "A1-   1,  -1,   1,  -1,   1,  -1,   1,   1,  -1,  -1\n",
      "T1-   3,   1,   0,   1,   1,   0,  -1,  -1,  -1,  -3\n"
     ]
    }
   ],
   "source": [
    "characters_trivial = character_table_row(matrix_trivial)\n",
    "characters_determinant = character_table_row(matrix_determinant)\n",
    "characters_vector = character_table_row(matrix_vector)\n",
    "\n",
    "print(inner_product(characters_trivial, characters_trivial)) # == 1 -> irreducible\n",
    "print(inner_product(characters_determinant, characters_determinant)) # == 1 -> irreducible\n",
    "print(inner_product(characters_vector, characters_vector)) # == 1 -> irreducible\n",
    "\n",
    "print(\"|C|\",fmt([len(c) for c in classes]))\n",
    "print(\"-\" * 60)\n",
    "print(\"A1+\",fmt(characters_trivial))\n",
    "print(\"A1-\",fmt(characters_determinant))\n",
    "print(\"T1-\",fmt(characters_vector)) # plus/minus for relative sign of first and last character of 1d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "7fa68d19",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "matrix_vector_projected_matrix_vector\n",
      "1.0\n",
      "|C|   1,   6,   8,   3,   6,   8,   6,   3,   6,   1\n",
      "------------------------------------------------------------\n",
      "A1+   1,   1,   1,   1,   1,   1,   1,   1,   1,   1\n",
      "A1-   1,  -1,   1,  -1,   1,  -1,   1,   1,  -1,  -1\n",
      "T1-   3,   1,   0,   1,   1,   0,  -1,  -1,  -1,  -3\n",
      "T1+   3,  -1,   0,  -1,   1,   0,  -1,  -1,   1,   3\n"
     ]
    }
   ],
   "source": [
    "matrix_a_vector_vector = projected_product_matrix(matrix_vector, matrix_vector, \n",
    "                                                         [antisymmetric_projector])\n",
    "test_representation(matrix_a_vector_vector)\n",
    "characters_a_vector_vector = character_table_row(matrix_a_vector_vector)\n",
    "print(inner_product(characters_a_vector_vector, characters_a_vector_vector)) # is irreducible\n",
    "\n",
    "print(\"|C|\",fmt([len(c) for c in classes]))\n",
    "print(\"-\" * 60)\n",
    "print(\"A1+\",fmt(characters_trivial))\n",
    "print(\"A1-\",fmt(characters_determinant))\n",
    "print(\"T1-\",fmt(characters_vector)) \n",
    "print(\"T1+\",fmt(characters_a_vector_vector))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "daa718e8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "matrix_vector_projected_matrix_vector\n",
      "1.0\n",
      "|C|   1,   6,   8,   3,   6,   8,   6,   3,   6,   1\n",
      "------------------------------------------------------------\n",
      "A1+   1,   1,   1,   1,   1,   1,   1,   1,   1,   1\n",
      "A1-   1,  -1,   1,  -1,   1,  -1,   1,   1,  -1,  -1\n",
      "T1-   3,   1,   0,   1,   1,   0,  -1,  -1,  -1,  -3\n",
      "T1+   3,  -1,   0,  -1,   1,   0,  -1,  -1,   1,   3\n",
      "T2+   3,   1,   0,  -1,  -1,   0,   1,  -1,  -1,   3\n"
     ]
    }
   ],
   "source": [
    "matrix_s_offdiag_vector_vector = projected_product_matrix(matrix_vector, matrix_vector, \n",
    "                                                         [symmetric_projector,invert_projector(diagonal_projector)])\n",
    "test_representation(matrix_s_offdiag_vector_vector)\n",
    "characters_s_offdiag_vector_vector = character_table_row(matrix_s_offdiag_vector_vector)\n",
    "print(inner_product(characters_s_offdiag_vector_vector, characters_s_offdiag_vector_vector)) # is irreducible\n",
    "\n",
    "print(\"|C|\",fmt([len(c) for c in classes]))\n",
    "print(\"-\" * 60)\n",
    "print(\"A1+\",fmt(characters_trivial))\n",
    "print(\"A1-\",fmt(characters_determinant))\n",
    "print(\"T1-\",fmt(characters_vector)) \n",
    "print(\"T1+\",fmt(characters_a_vector_vector))\n",
    "print(\"T2+\",fmt(characters_s_offdiag_vector_vector))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "be4eaf96",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "matrix_vector_projected_matrix_vector_X_matrix_determinant\n",
      "1.0\n",
      "|C|   1,   6,   8,   3,   6,   8,   6,   3,   6,   1\n",
      "------------------------------------------------------------\n",
      "A1+   1,   1,   1,   1,   1,   1,   1,   1,   1,   1\n",
      "A1-   1,  -1,   1,  -1,   1,  -1,   1,   1,  -1,  -1\n",
      "T1-   3,   1,   0,   1,   1,   0,  -1,  -1,  -1,  -3\n",
      "T1+   3,  -1,   0,  -1,   1,   0,  -1,  -1,   1,   3\n",
      "T2+   3,   1,   0,  -1,  -1,   0,   1,  -1,  -1,   3\n",
      "T2-   3,  -1,   0,   1,  -1,   0,   1,  -1,   1,  -3\n"
     ]
    }
   ],
   "source": [
    "# can find a partner irrep by multiplying any irrep with a 1d irrep\n",
    "matrix_det_s_offdiag_vector_vector = product_matrix(matrix_s_offdiag_vector_vector, matrix_determinant)\n",
    "test_representation(matrix_det_s_offdiag_vector_vector)\n",
    "characters_det_s_offdiag_vector_vector = character_table_row(matrix_det_s_offdiag_vector_vector)\n",
    "print(inner_product(characters_det_s_offdiag_vector_vector, characters_det_s_offdiag_vector_vector)) # is irreducible\n",
    "\n",
    "print(\"|C|\",fmt([len(c) for c in classes]))\n",
    "print(\"-\" * 60)\n",
    "print(\"A1+\",fmt(characters_trivial))\n",
    "print(\"A1-\",fmt(characters_determinant))\n",
    "print(\"T1-\",fmt(characters_vector)) \n",
    "print(\"T1+\",fmt(characters_a_vector_vector))\n",
    "print(\"T2+\",fmt(characters_s_offdiag_vector_vector))\n",
    "print(\"T2-\",fmt(characters_det_s_offdiag_vector_vector))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "b73367f1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "matrix_vector_projected_matrix_vector\n",
      "1.0000000000000004\n",
      "matrix_vector_projected_matrix_vector_X_matrix_determinant\n",
      "1.0000000000000004\n",
      "|C|   1,   6,   8,   3,   6,   8,   6,   3,   6,   1\n",
      "------------------------------------------------------------\n",
      "A1+   1,   1,   1,   1,   1,   1,   1,   1,   1,   1\n",
      "A1-   1,  -1,   1,  -1,   1,  -1,   1,   1,  -1,  -1\n",
      "T1-   3,   1,   0,   1,   1,   0,  -1,  -1,  -1,  -3\n",
      "T1+   3,  -1,   0,  -1,   1,   0,  -1,  -1,   1,   3\n",
      "T2+   3,   1,   0,  -1,  -1,   0,   1,  -1,  -1,   3\n",
      "T2-   3,  -1,   0,   1,  -1,   0,   1,  -1,   1,  -3\n",
      "E+    2,   0,  -1,   2,   0,  -1,   0,   2,   0,   2\n",
      "E-    2,   0,  -1,  -2,   0,   1,   0,   2,   0,  -2\n"
     ]
    }
   ],
   "source": [
    "matrix_diag_traceless_vector_vector = projected_product_matrix(matrix_vector, matrix_vector, \n",
    "                                                         [traceless_projector,diagonal_projector])\n",
    "test_representation(matrix_diag_traceless_vector_vector)\n",
    "characters_diag_traceless_vector_vector = character_table_row(matrix_diag_traceless_vector_vector)\n",
    "print(inner_product(characters_diag_traceless_vector_vector, characters_diag_traceless_vector_vector)) # is irreducible\n",
    "\n",
    "matrix_det_diag_traceless_vector_vector = product_matrix(matrix_diag_traceless_vector_vector, matrix_determinant)\n",
    "test_representation(matrix_det_diag_traceless_vector_vector)\n",
    "characters_det_diag_traceless_vector_vector = character_table_row(matrix_det_diag_traceless_vector_vector)\n",
    "print(inner_product(characters_det_diag_traceless_vector_vector, characters_det_diag_traceless_vector_vector)) # is irreducible\n",
    "\n",
    "print(\"|C|\",fmt([len(c) for c in classes]))\n",
    "print(\"-\" * 60)\n",
    "print(\"A1+\",fmt(characters_trivial))\n",
    "print(\"A1-\",fmt(characters_determinant))\n",
    "print(\"T1-\",fmt(characters_vector)) \n",
    "print(\"T1+\",fmt(characters_a_vector_vector))\n",
    "print(\"T2+\",fmt(characters_s_offdiag_vector_vector))\n",
    "print(\"T2-\",fmt(characters_det_s_offdiag_vector_vector))\n",
    "print(\"E+ \",fmt(characters_diag_traceless_vector_vector))\n",
    "print(\"E- \",fmt(characters_det_diag_traceless_vector_vector))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "fd65239e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "matrix_vector_X_matrix_vector_X_matrix_vector_projected\n",
      "1.0000000000000007\n",
      "matrix_vector_X_matrix_vector_X_matrix_vector_projected_X_matrix_determinant\n",
      "1.0000000000000007\n",
      "|C|   1,   6,   8,   3,   6,   8,   6,   3,   6,   1\n",
      "------------------------------------------------------------\n",
      "A1+   1,   1,   1,   1,   1,   1,   1,   1,   1,   1\n",
      "A1-   1,  -1,   1,  -1,   1,  -1,   1,   1,  -1,  -1\n",
      "T1-   3,   1,   0,   1,   1,   0,  -1,  -1,  -1,  -3\n",
      "T1+   3,  -1,   0,  -1,   1,   0,  -1,  -1,   1,   3\n",
      "T2+   3,   1,   0,  -1,  -1,   0,   1,  -1,  -1,   3\n",
      "T2-   3,  -1,   0,   1,  -1,   0,   1,  -1,   1,  -3\n",
      "E+    2,   0,  -1,   2,   0,  -1,   0,   2,   0,   2\n",
      "E-    2,   0,  -1,  -2,   0,   1,   0,   2,   0,  -2\n",
      "A2-   1,   1,   1,  -1,  -1,  -1,  -1,   1,   1,  -1\n",
      "A2+   1,  -1,   1,   1,  -1,   1,  -1,   1,  -1,   1\n"
     ]
    }
   ],
   "source": [
    "matrix_last_irrep_minus = orthogonal_matrix(\n",
    "    product_matrix(product_matrix(matrix_vector, matrix_vector), matrix_vector),\n",
    "[\n",
    "    characters_trivial,\n",
    "    characters_determinant,\n",
    "    characters_vector,\n",
    "    characters_a_vector_vector,\n",
    "    characters_s_offdiag_vector_vector,\n",
    "    characters_diag_traceless_vector_vector,\n",
    "    characters_det_s_offdiag_vector_vector,\n",
    "    characters_diag_traceless_vector_vector,\n",
    "    characters_det_diag_traceless_vector_vector\n",
    "])\n",
    "\n",
    "test_representation(matrix_last_irrep_minus)\n",
    "characters_last_irrep_minus = character_table_row(matrix_last_irrep_minus)\n",
    "print(inner_product(characters_last_irrep_minus, characters_last_irrep_minus)) # irreducible!\n",
    "\n",
    "matrix_last_irrep_plus = product_matrix(matrix_last_irrep_minus, matrix_determinant)\n",
    "test_representation(matrix_last_irrep_plus)\n",
    "characters_last_irrep_plus = character_table_row(matrix_last_irrep_plus)\n",
    "print(inner_product(characters_last_irrep_plus, characters_last_irrep_plus)) # is irreducible\n",
    "\n",
    "\n",
    "print(\"|C|\",fmt([len(c) for c in classes]))\n",
    "print(\"-\" * 60)\n",
    "print(\"A1+\",fmt(characters_trivial))\n",
    "print(\"A1-\",fmt(characters_determinant))\n",
    "print(\"T1-\",fmt(characters_vector)) \n",
    "print(\"T1+\",fmt(characters_a_vector_vector))\n",
    "print(\"T2+\",fmt(characters_s_offdiag_vector_vector))\n",
    "print(\"T2-\",fmt(characters_det_s_offdiag_vector_vector))\n",
    "print(\"E+ \",fmt(characters_diag_traceless_vector_vector))\n",
    "print(\"E- \",fmt(characters_det_diag_traceless_vector_vector))\n",
    "print(\"A2-\",fmt(characters_last_irrep_minus))\n",
    "print(\"A2+\",fmt(characters_last_irrep_plus))\n",
    "\n",
    "# And we have the complete character table with all irreps of the full cubic group!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "d2b38f46",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We end this chapter by demonstrating the orthonormality of the irreps\n",
    "characters = [\n",
    "    characters_trivial,\n",
    "    characters_determinant,\n",
    "    characters_vector,\n",
    "    characters_a_vector_vector,\n",
    "    characters_s_offdiag_vector_vector,\n",
    "    characters_det_s_offdiag_vector_vector,\n",
    "    characters_diag_traceless_vector_vector,\n",
    "    characters_det_diag_traceless_vector_vector,\n",
    "    characters_last_irrep_minus,\n",
    "    characters_last_irrep_plus\n",
    "]\n",
    "\n",
    "for i in range(len(characters)):\n",
    "    for j in range(len(characters)):\n",
    "        eps = abs(inner_product(characters[i],characters[j]) - (1 if i == j else 0))\n",
    "        assert eps < 1e-13\n",
    "        \n",
    "# We can now use this character table to project out any irrep of any given product representation!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9cbdbca",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
