2.5 Scrambling and Winding.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to Make a Wormhole\n",
    "\n",
    "## Part 5: Scrambling and Winding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider a bunch of qubits evolving under some random unitary operator. If we wanted to run this on a quantum computer, for example, we'd want to compile down that unitary into a series of, say, one and two qubit operators that we apply in a sequence since these are simpler to implement. Consider a single qubit in our group, which starts off separable from the rest. As time goes on, maybe it interacts with another qubit: then, generically, the two qubits will get entangled. Each of those qubits then goes off and gets entangled with more qubits, and now it's possible all four are entangled. And so on. In other words, entanglement spreads between quantum systems like an epidemic.\n",
    "\n",
    "Here's an easy way to get a handle on this. Below is a python implementation of above: a bunch of qubits start of separable, and we evolve them under some randomly selected time evolution. We want to track how one qubit gets entangled with the rest. One easy way to do this is to instead start off with a maximally entangled pair of qubits. One stays outside the group, the other one evolves with the rest of the group. We can then calculate the entanglement between the one who stays outside and all the ones in the group in turn, and track it over time. Think of it like we have two qubits connected by Ariadne's thread; we toss one in and it gets scrambled into the rest, but we can keep track of it by following the thread home.\n",
    "\n",
    "What's a good measure of entanglement though? A classical measure is the quantum mutual information:\n",
    "\n",
    "$$ QMI(A, B) = S(A) + S(B) - S(A, B)$$\n",
    "$$ S = - tr( \\rho \\ ln \\rho ) = - \\sum_{\\lambda} \\lambda  \\ ln \\lambda $$\n",
    "\n",
    "\\\\(S\\\\) is called the Von Neumann entropy. In the latter part the sum is over eigenvalues \\\\( \\lambda \\\\) of \\\\( \\rho \\\\)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import qutip as qt\n",
    "\n",
    "n = 4\n",
    "dt = 0.01\n",
    "t_max = 1000\n",
    "\n",
    "ref, rest = qt.bell_state(\"00\"), qt.tensor(*[qt.rand_ket(2) for i in range(n)])\n",
    "state = qt.tensor(ref, rest)\n",
    "state.dims = [[2]*(n+2), [1]*(n+2)]\n",
    "\n",
    "E = qt.rand_herm(2**(n+1))\n",
    "U = qt.tensor(qt.identity(2), (-1j*dt*E).expm())\n",
    "U.dims = [[2]*(n+2), [2]*(n+2)]\n",
    "\n",
    "entropies = [qt.entropy_mutual((state*state.dag()).ptrace((0, i)), 0, 1) for i in range(1, n+2)]\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.set_xlabel('qubit')\n",
    "ax.set_ylabel('entanglement')\n",
    "\n",
    "bar = ax.bar(list(range(1, n+2)), entropies)\n",
    "\n",
    "for i in range(t_max):\n",
    "    state = U*state\n",
    "    entropies = [qt.entropy_mutual((state*state.dag()).ptrace((0, i)), 0, 1) for i in range(1, n+2)]\n",
    "    for b, e in zip(bar, entropies):\n",
    "        b.set_height(e)\n",
    "    fig.canvas.draw()\n",
    "    plt.pause(0.0001)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just like we can track a qubit as it gets scrambled among the rest, we could also keep track of an operator, watch how it starts off localized and then gets spread out over the whole Hilbert space.\n",
    "\n",
    "So we're interested in tracking in some sense the size/complexity of the operator \\\\( O(t) = e^{iEt} O e^{-iEt} \\\\) as it evolves in time.\n",
    "\n",
    "Now just as quantum states form a vector space, operators also form a vector space. For example, I can write any 2x2 Hermitian matrix in the following 4-vector form:\n",
    "\n",
    "$$ H = tI + xX + yY + zZ $$\n",
    "\n",
    "So given an operator H I can expand it in terms of a real (t, x, y, z) vector. One way to do this is to use the generalization of the inner product for matrices: \\\\( (tr(H), tr(X^{\\dagger}H), tr(Y^{\\dagger}H), tr(Z^{\\dagger}H)) \\\\). But we know another way: We could multiply I, X, Y, Z and H all by the cup state and then take the normal inner product: \\\\( ( \\langle I\\cup \\mid H \\cup \\rangle , \\langle X \\cup \\mid H \\cup \\rangle , \\langle Y \\cup \\mid H \\cup \\rangle , \\langle Z \\cup \\mid H \\cup \\rangle ) \\\\).\n",
    "\n",
    "I could decompose a 2x2 unitary in the same way: the components would then be generally complex. Cool, but does the same work for the nxn case? Yes, there's always some basis for Hermitian matrices. But we want to consider spaces of qubits and it turns out that the following 16 operators forms a basis for two qubit operators:\n",
    "\n",
    "$$ \\begin{matrix}II & IX & IY & IZ \\\\\n",
    "                 XI & XX & XY & XZ \\\\\n",
    "                 YI & YX & YY & YZ \\\\\n",
    "                 ZI & ZX & ZY & ZZ\\end{matrix} $$\n",
    "\n",
    "And then same goes for the n-qubit case: we take all arrangements of the Pauli operators, and we can expand any operator in the so-called Pauli basis. The point is that in the Pauli basis we can define the notion of the size of an operator. For the 1 qubit case, \\\\(I\\\\) has size 0, while \\\\(X, Y, Z\\\\) each have size 1. In the 2 qubit case, we have 1 operator of size 0, 6 operators of size 1, and 9 operators of size 2. We'll write \\\\( |P| \\\\) for the size of a Pauli operator P. \n",
    "\n",
    "So we have the idea that generically operators grow in size over time. Let's watch it happen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "import matplotlib.pyplot as plt\n",
    "import qutip as qt\n",
    "import numpy as np\n",
    "from itertools import product\n",
    "\n",
    "def pauli_basis(n):\n",
    "    IXYZ = {\"I\": qt.identity(2), \"X\": qt.sigmax(), \"Y\": qt.sigmay(), \"Z\": qt.sigmaz()}\n",
    "    names, ops = [], []\n",
    "    for P in product(IXYZ, repeat=n):\n",
    "        names.append(\"\".join(P))\n",
    "        ops.append(qt.tensor(*[IXYZ[p] for p in P]))\n",
    "    return names, ops\n",
    "\n",
    "def to_pauli(op, Pops):\n",
    "    return np.array([(o.dag()*op).tr() for o in Pops])/np.sqrt(len(Pops))\n",
    "\n",
    "n = 2\n",
    "dt = 0.01\n",
    "Pnames, Pops = pauli_basis(n)\n",
    "\n",
    "H = qt.rand_herm(2**n)\n",
    "H.dims = [[2]*n, [2]*n]\n",
    "U = (-1j*dt*H).expm()\n",
    "O = qt.tensor(*[qt.sigmax() if i == 0 else qt.identity(2) for i in range(n)])\n",
    "\n",
    "probabilities = [(a*a.conj()).real for a in to_pauli(O, Pops)]\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.set_xlabel('paulis')\n",
    "ax.set_ylabel('probabilities')\n",
    "plt.xticks(np.arange(len(Pnames)), Pnames)\n",
    "bar = ax.bar(list(range(len(Pnames))), probabilities)\n",
    "\n",
    "t_max = 5000\n",
    "for i in range(t_max):\n",
    "    O = U.dag()*O*U\n",
    "    probabilities = [(a*a.conj()).real for a in to_pauli(O, Pops)]\n",
    "    for b, p in zip(bar, probabilities):\n",
    "        b.set_height(p)\n",
    "    fig.canvas.draw()\n",
    "    plt.pause(0.0001)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In fact, we could think about the operator as a quantum particle undergoing its own Schrodinger evolution but defined on a graph where nodes are Pauli matrices and node are connected by edges when they are one operator different from another.\n",
    "\n",
    "<hr>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So now we can put some interesting facts together. We know that:\n",
    "\n",
    "$$ O_{L} (-t) \\mid TFD \\rangle \\to \\rho_{\\beta, R}^{\\frac{1}{2}} O_{R}^{T} (t) \\mid \\cup \\rangle $$\n",
    "$$ O_{R}^{T} (t) \\mid TFD \\rangle \\to O_{R}^{T} (t) \\rho_{\\beta, R}^{\\frac{1}{2}} \\mid \\cup \\rangle $$\n",
    "\n",
    "And our general setup is:\n",
    "\n",
    "![](img/freed_operator_and_state.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The question is what properties does our \\\\(E)\\\\ have to have for this to work? In the state transfer picture, we throw a qubit in on the left in the past, which then should get completely scrambled over the whole left side. To the extend that the qubit is scrambled over the whole thing, we could say that the qubit is \"deeper inside\" the black hole. Then our coupling should somehow transfer our deep qubit from the left to the right, and then the forward time evolution on the right side should unscramble it and Bob's your uncle.\n",
    "\n",
    "Following [arXiv:1911.06314](https://arxiv.org/abs/1911.06314), it turns out that energy operators that have the right scrambling properties exhibit something called \"perfect size winding.\" And even if the size winding isn't perfect, you can still often get stuff through.\n",
    "\n",
    "The authors of the aforementioned paper provide the following ansatz:\n",
    "\n",
    "$$ \\big( \\rho_{\\beta, R}^{\\frac{1}{2}} O_{R}^{T} (t) \\big)(t) = \\sum_{Paulis} e^{i \\alpha \\frac{|P|}{n}}r_{p}P $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here the sum is over all the n-qubit Paulis where n is the number of qubits on the right side! Since we can see that this is an operator that only acts on the right. \\\\(r_{p}\\\\) is some real coefficient, \\\\(P\\\\) is the given Pauli, and \\\\( \\alpha \\\\) is some overall phase. (Incidentally, in the limit of large size for a random P, \\\\( \\alpha = \\frac{4}{3}g\\\\)--we'll come back to the \\\\(g\\\\).) We can see that in the Pauli basis, each Pauli basis vector has a phase proportional to its size. In other words, for perfect sizing winding Hamiltonians, in the Pauli basis, the phases wind around at a speed proportional to its size. \n",
    "\n",
    "Considering that \\\\(  \\big( \\rho_{\\beta, R}^{\\frac{1}{2}} O_{R}^{T} (t) \\big)^{\\dagger} =  O_{R}^{T} (t) \\rho_{\\beta, R}^{\\frac{1}{2}}  \\\\), we can say:\n",
    "\n",
    "$$ O_{L} (-t) \\mid TFD \\rangle = \\sum_{Paulis} e^{i \\alpha \\frac{|P|}{n}}r_{p}\\mid P \\rangle$$\n",
    "$$ O_{R}^{T} (t) \\mid TFD \\rangle = \\sum_{Paulis} e^{-i \\alpha \\frac{|P|}{n}}r_{p}\\mid P \\rangle$$\n",
    "\n",
    "Which is of course:\n",
    "\n",
    "$$ O_{L} (-t) \\mid TFD \\rangle = \\sum_{Paulis} e^{i \\alpha \\frac{|P|}{n}}r_{p}P_{L}\\mid \\cup \\rangle$$\n",
    "$$ O_{R}^{T} (t) \\mid TFD \\rangle = \\sum_{Paulis} e^{-i \\alpha \\frac{|P|}{n}}r_{p}P_{L}\\mid \\cup \\rangle$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The only difference between \\\\( O_{L} (-t)  \\\\) and \\\\( O_{R}^{T} (t) \\\\) is the direction of the winding, whether clockwise or counterclockwise.\n",
    "\n",
    "Recall our choice of V operator from before.\n",
    "\n",
    "$$ V = \\frac{1}{b-a} \\sum_{i}^{b} Z_{i, L} Z_{b-i, R} $$\n",
    "\n",
    "Its role is simply to unwind the phases. It turns out that the ZZ coupling is sensitive enough to size to do the job. \n",
    "\n",
    "Finally, when \\\\( \\beta = 0 \\\\), in other words, at infinite temperature when the TFD reduces down to the \\\\(\\cup \\\\), \\\\( \\rho_{\\beta, R}^{\\frac{1}{2}} O_{R}^{T} \\\\) is Hermitian and has real expansion coefficients in the Pauli basis. And so no unwinding is necessary!\n",
    "\n",
    "<hr />"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "VPython",
   "language": "python",
   "name": "vpython"
  },
  "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}