diff --git a/exercises/cv_submissions/cv_tutorial_luis.ipynb b/exercises/cv_submissions/cv_tutorial_luis.ipynb new file mode 100644 index 0000000..bd11680 --- /dev/null +++ b/exercises/cv_submissions/cv_tutorial_luis.ipynb @@ -0,0 +1,736 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "7858dc18", + "metadata": {}, + "source": [ + "# Continuous Variable Quantum Information (One Mode)\n", + "In continuous variable quantum information the information carriers are the modes of a quantum harmonic oscillator. In general one mode can be described with the Hamiltonian\n", + "$$\n", + "\\hat{H} = \\hbar \\omega \\left( \\hat{n} +\\frac{1}{2} \\right),\n", + "$$\n", + "where $\\omega$ is the frequency of the oscillation and $\\hat{n}$ is the number operator.\n", + "One natural option for states to work with is the number or Fock basis $\\ket{n}$, it consists of the eigenstates of the Hamiltonian. Any one mode state can be decomposed in this basis as\n", + "$$\n", + "\\ket{\\psi} = \\sum_{n\\geq0} c_n \\ket{n}, \\qquad \\hat{H}\\ket{\\psi} = \\sum_{n\\geq0} n \\bullet c_n \\ket{n}\n", + "$$\n", + "### Phasespace\n", + "The ladder operators $(\\hat{a}, \\hat{a}^\\dagger)$ fulfill the relations\n", + "$$\n", + "\\hat{a}\\ket{n} = \\sqrt{n}\\ket{n-1}, \\qquad \\hat{a}^\\dagger\\ket{n} = \\sqrt{n+1}\\ket{n+1}.\n", + "$$\n", + "Most importantly they fulfill the commutation relation $[\\hat{a},\\hat{a}^\\dagger]=1$. With these one can define the quadrature operators $(\\hbar=1)$\n", + "$$\n", + "\\hat{q}(\\theta) = \\frac{1}{2} \\left( e^{-i\\theta}\\hat{a} + e^{i\\theta}\\hat{a}^\\dagger \\right).\n", + "$$\n", + "The quadratures $\\hat{x}\\equiv\\hat{q}(0)$ and $\\hat{q}\\equiv\\hat{q}(\\frac{\\pi}{2})$ are called position/momentum or amplitude/phase quadratures. They fulfill the commutation relation $[\\hat{x},\\hat{p}]=\\frac{1}{2}$, similar to canonical conjugate variables in classical mechanics.\n", + "## State representation\n", + "CV quantum states can be represented in different ways. In the following some are discussed.\n", + "### Density Matrix\n", + "Any state can be represented in the form of a density matrix operator\n", + "$$\n", + "\\hat{\\rho} = \\ket{\\psi}\\bra{\\psi} = \\sum_{n,m\\geq0} c_n \\bar{c}_m \\ket{n}\\bra{m},\n", + "$$\n", + "thus the state can be represented as a infinitely long matrix. Usually to represent the state one truncates the matrix with a cutoff $N$\n", + "$$\n", + "\\hat{\\rho} \\approx \\sum_{n,m=0}^{N} c_n \\bar{c}_m \\ket{n}\\bra{m}.\n", + "$$\n", + "### Wigner Function\n", + "The Wigner function is a different but equivalent way to describe CV states. It is a quasiprobability function in phasespace $W(x,p)$ and can be written in terms of the density matrix\n", + "$$\n", + "W(x,p)=\\frac{1}{(2\\pi)^2} \\int_{-\\infty}^{\\infty}\\int_{-\\infty}^{\\infty} e^{-i(xa-pb)} Tr \\left[ \\hat{\\rho}e^{i(a\\hat{x}-b\\hat{p})} \\right] da db.\n", + "$$\n", + "### Wavefunction\n", + "The wavefunction of a state is defined as\n", + "$$\n", + "\\psi(x) = \\braket{x|\\psi},\n", + "$$\n", + "where $\\ket{x}$ is the eigenstate of the quadrature operator $\\hat{x}$. It can also be defined as the integral\n", + "$$\n", + "\\psi(x) = \\int_{-\\infty}^{\\infty} W(x,p) dp.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "45898f24", + "metadata": {}, + "source": [ + "### Import necesary libraries" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "b4b108df", + "metadata": {}, + "outputs": [], + "source": [ + "# Import libraries\n", + "import numpy as np\n", + "import math\n", + "\n", + "DEFAULT_CUTOFF = 10\n", + "\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib.colors import TwoSlopeNorm\n", + "\n", + "# Laguerre Polynomials L_k^a\n", + "from scipy.special import genlaguerre\n", + "\n", + "# Interactive widgets\n", + "from ipywidgets import (\n", + " IntSlider, FloatRangeSlider, IntText, FloatText,\n", + " Button, VBox, HBox, Output, GridspecLayout, Tab\n", + ")\n", + "from IPython.display import display, clear_output" + ] + }, + { + "cell_type": "markdown", + "id": "a1f4fb37", + "metadata": {}, + "source": [ + "## States from Density matrix\n", + "Here it is used that the Wigner function of the single $\\ket{n}\\bra{m}$ can be analitically derived to be:\n", + "$$\n", + "W_{nm}(x,p) = \\frac{1}{\\pi} e^{-x^2-p^2} (-1)^n (x - ip)^{m-n} \\sqrt{2^{m-n} \\frac{n!}{m!}} L^{m-n}_n (2x^2 + 2p^2).\n", + "$$\n", + "This expression can be efficiently calculated for any pair of (n,m). This can be used to find the Wigner function of any state that is given by the density matrix in the fock basis up to some cutoff $N$\n", + "$$\n", + "W_{\\hat{\\rho}} \\approx \\sum_{n,m=0}^{N} \\rho_{nm}W_{nm}\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "6b4cda34-aa71-4cf2-9164-df6f5a176331", + "metadata": {}, + "outputs": [], + "source": [ + "# Obtain the Wigner function of a Fock state W_|n> tol:\n", + "\n", + " W_nm = Wigner_Fock(n,m)\n", + "\n", + " if n == m:\n", + " terms.append((rho[n,m], W_nm))\n", + " else:\n", + " # use W_mn = conj(W_nm)\n", + " terms.append((2*rho[n,m], W_nm))\n", + "\n", + " def W(x,p):\n", + " result = 0\n", + " for coeff, W_nm in terms:\n", + " result += np.real(coeff * W_nm(x,p))\n", + " return result\n", + "\n", + " return W" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "da428858", + "metadata": {}, + "outputs": [], + "source": [ + "# Evaluate the Wigner function on a grid of points\n", + "def evaluate_Wigner(W, X_points: int = 100, X_range: list = [-5.0, 5.0], \n", + " P_points: int = None, P_range: list = None):\n", + " X_axis = np.linspace(X_range[0], X_range[1], X_points)\n", + " P_axis = np.linspace(P_range[0], P_range[1], P_points) if P_range is not None else X_axis\n", + " W_values = np.array([[W(x, p).real for x in X_axis] for p in P_axis])\n", + " return X_axis, P_axis, W_values\n", + "\n", + "# Plot the Wigner function on a grid of points\n", + "def plot_Wigner(W, \n", + " X_min: float, X_max: float, X_points: int, \n", + " P_min: float = None, P_max: float = None, P_points: int = None, \n", + " title='Wigner function'):\n", + " X_axis, P_axis, W_values = evaluate_Wigner(W, X_points=X_points, X_range=[X_min, X_max], \n", + " P_points=P_points, \n", + " P_range=[P_min, P_max] if P_min is not None and P_max is not None else None)\n", + " # Set colormap with white for 0 and positive to red (negative to blue)\n", + " m = np.max(np.abs(W_values))\n", + " norm = TwoSlopeNorm(vmin=-m, vcenter=0.0, vmax=m)\n", + "\n", + " plt.figure(figsize=(8, 6))\n", + " plt.contourf(X_axis, P_axis, W_values, levels=100, cmap='RdBu_r', norm=norm)\n", + " plt.colorbar(label=r'$W(x, p)$')\n", + " plt.xlabel(r'$x$')\n", + " plt.ylabel(r'$p$')\n", + " plt.title(title)\n", + " plt.show()\n", + "\n", + "# Plot the density matrix of a state in the Fock basis\n", + "def plot_density_matrix(rho, title=r'Density matrix \\hat{\\rho} in Fock basis'):\n", + " plt.figure(figsize=(6, 5))\n", + " plt.imshow(np.abs(rho), cmap='viridis', interpolation='nearest')\n", + " plt.colorbar(label=r'$|\\rho_{nm}|$')\n", + " plt.xlabel(r'$n$')\n", + " plt.ylabel(r'$m$')\n", + " plt.title(title)\n", + " plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "8d987fbe", + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "f4a3ef2f511c449dac397c62f52e4f6d", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "VBox(children=(IntSlider(value=2, description='Size N:', max=10, min=1), FloatRangeSlider(value=(-5.0, 5.0), d…" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Interactive plotting cell\n", + "################################################\n", + "# Size of the density matrix / state vector\n", + "size_slider = IntSlider(\n", + " value=2,\n", + " min=1,\n", + " max=10,\n", + " step=1,\n", + " description='Size N:'\n", + ")\n", + "\n", + "################################################\n", + "# Plotting parameters\n", + "X_range_slider = FloatRangeSlider(\n", + " value=[-5.0, 5.0],\n", + " min=-10.0,\n", + " max=10.0,\n", + " step=0.1,\n", + " description='X range:'\n", + ")\n", + "\n", + "P_range_slider = FloatRangeSlider(\n", + " value=[-5.0, 5.0],\n", + " min=-10.0,\n", + " max=10.0,\n", + " step=0.1,\n", + " description='P range:'\n", + ")\n", + "\n", + "num_points_slider = IntText(\n", + " value=100,\n", + " description='Num points:'\n", + ")\n", + "\n", + "################################################\n", + "# Containers\n", + "\n", + "matrix_container = VBox()\n", + "sv_container = VBox()\n", + "\n", + "output_DM = Output()\n", + "output_SV = Output()\n", + "\n", + "\n", + "################################################\n", + "# Density matrix widgets\n", + "\n", + "def create_density_matrix_widgets(N):\n", + "\n", + " grid = GridspecLayout(N, N)\n", + "\n", + " widgets = []\n", + "\n", + " for i in range(N):\n", + "\n", + " row = []\n", + "\n", + " for j in range(N):\n", + "\n", + " widget = FloatText(\n", + " value=0.0,\n", + " description=f'ρ[{i},{j}]',\n", + " layout={'width': '140px'}\n", + " )\n", + "\n", + " grid[i, j] = widget\n", + " row.append(widget)\n", + "\n", + " widgets.append(row)\n", + "\n", + " return widgets, grid\n", + "\n", + "\n", + "################################################\n", + "# State vector widgets\n", + "\n", + "def create_state_vector_widgets(N):\n", + "\n", + " # Two columns:\n", + " # Re[x_j], Im[x_j]\n", + "\n", + " grid = GridspecLayout(N + 1, 2)\n", + "\n", + " widgets = []\n", + "\n", + " # Column headers\n", + " grid[0, 0] = HBox([])\n", + " grid[0, 1] = HBox([])\n", + "\n", + " for j in range(N):\n", + "\n", + " re_widget = FloatText(\n", + " value=0.0,\n", + " description=f'Re[x{j}]',\n", + " layout={'width': '160px'}\n", + " )\n", + "\n", + " im_widget = FloatText(\n", + " value=0.0,\n", + " description=f'Im[x{j}]',\n", + " layout={'width': '160px'}\n", + " )\n", + "\n", + " grid[j + 1, 0] = re_widget\n", + " grid[j + 1, 1] = im_widget\n", + "\n", + " widgets.append(\n", + " (re_widget, im_widget)\n", + " )\n", + "\n", + " return widgets, grid\n", + "\n", + "\n", + "################################################\n", + "# Update density matrix\n", + "\n", + "def update_matrix(change=None):\n", + "\n", + " global rho_widgets\n", + "\n", + " N = size_slider.value\n", + "\n", + " rho_widgets, grid = create_density_matrix_widgets(N)\n", + "\n", + " matrix_container.children = [grid]\n", + "\n", + "################################################\n", + "# Update state vector\n", + "\n", + "def update_state_vector(change=None):\n", + "\n", + " global sv_widgets\n", + "\n", + " N = size_slider.value\n", + "\n", + " sv_widgets, grid = create_state_vector_widgets(N)\n", + "\n", + " sv_container.children = [grid]\n", + "\n", + "\n", + "################################################\n", + "# Update both grids whenever N changes\n", + "\n", + "size_slider.observe(\n", + " update_matrix,\n", + " names='value'\n", + ")\n", + "\n", + "size_slider.observe(\n", + " update_state_vector,\n", + " names='value'\n", + ")\n", + "\n", + "\n", + "# Create initial widgets\n", + "update_matrix()\n", + "update_state_vector()\n", + "\n", + "\n", + "################################################\n", + "# Density Matrix: Run button\n", + "\n", + "run_button_DM = Button(\n", + " description='Run Wigner calculation',\n", + " button_style='success'\n", + ")\n", + "\n", + "\n", + "def on_run_button_clicked_DM(b):\n", + "\n", + " N = size_slider.value\n", + "\n", + " # Construct density matrix\n", + "\n", + " rho = np.array([[rho_widgets[i][j].value\n", + " for j in range(N)] for i in range(N)],dtype=complex)\n", + "\n", + "\n", + " # Plotting parameters\n", + "\n", + " X_range = X_range_slider.value\n", + " P_range = P_range_slider.value\n", + "\n", + " X_points = num_points_slider.value\n", + " P_points = num_points_slider.value\n", + "\n", + " # Calculate and plot\n", + "\n", + " with output_DM:\n", + "\n", + " clear_output(wait=True)\n", + "\n", + " print(\"Density matrix:\")\n", + " print(rho)\n", + "\n", + " # Create Wigner representation\n", + " W = Wigner_DensityMatrix(rho)\n", + "\n", + " # Evaluate Wigner function\n", + " X, P, W_values = evaluate_Wigner(\n", + " W,\n", + " X_points=X_points,\n", + " X_range=X_range,\n", + " P_points=P_points,\n", + " P_range=P_range\n", + " )\n", + "\n", + " # Plot\n", + " plot_Wigner(\n", + " W,\n", + " X_min=X_range[0],\n", + " X_max=X_range[1],\n", + " P_min=P_range[0],\n", + " P_max=P_range[1],\n", + " X_points=X_points,\n", + " P_points=P_points\n", + " )\n", + "\n", + "\n", + "run_button_DM.on_click(\n", + " on_run_button_clicked_DM\n", + ")\n", + "\n", + "\n", + "################################################\n", + "# State Vector: Run button\n", + "\n", + "run_button_SV = Button(\n", + " description='Run Wigner calculation',\n", + " button_style='success'\n", + ")\n", + "\n", + "\n", + "def on_run_button_clicked_SV(b):\n", + "\n", + " N = size_slider.value\n", + " # Construct state vector\n", + "\n", + " sv = np.array(\n", + " [\n", + " sv_widgets[j][0].value\n", + " + 1j * sv_widgets[j][1].value\n", + " for j in range(N)\n", + " ],\n", + " dtype=complex\n", + " )\n", + "\n", + " # Check that the state vector is not zero\n", + "\n", + " norm = np.linalg.norm(sv)\n", + "\n", + " with output_SV:\n", + "\n", + " clear_output(wait=True)\n", + "\n", + " if norm == 0:\n", + "\n", + " print(\n", + " \"Error: the state vector cannot be the zero vector.\"\n", + " )\n", + "\n", + " return\n", + "\n", + " # Normalize state vector\n", + "\n", + " sv = sv / norm\n", + " \n", + " # Construct density matrix rho = |psi>log for further details." + ] + } + ], + "source": [ + "# Gaussian states\n", + "\n", + "##################################################################\n", + "# Gaussian operations: Symplectic matrices\n", + "def One_Mode_Squeeze(r,theta = 0):\n", + " \"\"\"Single mode squeezing symplectic transformation\"\"\"\n", + " S = np.array([[np.cosh(r)-np.sinh(r)*np.cos(theta), -np.sinh(r)*np.sin(theta)],\n", + " [-np.sinh(r)*np.sin(theta), np.cosh(r)+np.sinh(r)*np.cos(theta)]])\n", + " return S\n", + "\n", + "def Phase_rotation(theta):\n", + " \"\"\"Single mode phase rotation symplectic transformation\"\"\"\n", + " R = np.array([[np.cos(theta), -np.sin(theta)],\n", + " [np.sin(theta), np.cos(theta)]])\n", + " return R\n", + "\n", + "def One_Mode_Symplectic(theta, r, phi):\n", + " \"\"\"General single mode symplectic transformation\"\"\"\n", + " return np.einsum('ij,jk,kl->il', Phase_rotation(theta), One_Mode_Squeeze(r), Phase_rotation(phi))\n", + "\n", + "def Beam_splitter(theta):\n", + " \"\"\"Two mode beam splitter symplectic transformation\"\"\"\n", + " tau = np.cos(theta)**2\n", + " BS = np.array([[np.sqrt(tau), 0, np.sqrt(1-tau), 0],\n", + " [0, np.sqrt(tau), 0, np.sqrt(1-tau)],\n", + " [-np.sqrt(1-tau), 0, np.sqrt(tau), 0],\n", + " [0, -np.sqrt(1-tau), 0, np.sqrt(tau)]])\n", + " return BS\n", + "\n", + "def Two_Mode_Squeeze(r, theta = 0):\n", + " \"\"\"Two mode squeezing symplectic transformation\"\"\"\n", + " S = np.array([[np.cosh(r), 0, np.sinh(r)*np.cos(theta), np.sinh(r)*np.sin(theta)],\n", + " [0, np.cosh(r), np.sinh(r)*np.sin(theta), -np.sinh(r)*np.cos(theta)],\n", + " [np.sinh(r)*np.cos(theta), np.sinh(r)*np.sin(theta), np.cosh(r), 0],\n", + " [np.sinh(r)*np.sin(theta), -np.sinh(r)*np.cos(theta), 0, np.cosh(r)]])\n", + " return S\n", + "\n", + "def Controlled_Z(phi):\n", + " \"\"\"Two mode controlled-Z gate symplectic transformation\"\"\"\n", + " CZ = np.array([[1, 0, 0, 0],\n", + " [0, 1, phi, 0],\n", + " [0, 0, 1, 0],\n", + " [phi, 0, 0, 1]])\n", + " return CZ\n", + "\n", + "##################################################################\n", + "\n", + "# Obtain the Wigner function of a one mode Gaussian state W_(V, d) \n", + "def Wigner_Gaussian(V,d=None):\n", + " N = V.shape[0] // 2\n", + " if d is None:\n", + " d = np.zeros(V.shape[0])\n", + " detV = np.linalg.det(V)\n", + " invV = np.linalg.inv(V)\n", + " W_G = lambda xi: (1 / (np.pi**N * np.sqrt(detV))) \\\n", + " * np.exp(-0.5 * (xi - d).T @ invV @ (xi - d))\n", + " def W_G_mode(m, x, p):\n", + " xi = np.zeros(2 * N)\n", + " xi[2 * m] = x\n", + " xi[2 * m + 1] = p\n", + " return W_G(xi)\n", + " return W_G_mode\n", + "\n", + "\n", + "##################################################################\n", + "# Obtain the density matrix of a Gaussian up to cutoff\n", + "def DensityMatrix_Gaussian(V, d=None, mode=0, cutoff=DEFAULT_CUTOFF):\n", + " \"\"\"\n", + " Obtain the density matrix of the selected Gaussian mode in the Fock basis.\n", + " \"\"\"\n", + " # Todo" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "66d34474", + "metadata": {}, + "outputs": [], + "source": [ + "# Interactive plotting cell for Gaussian states\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "01080458", + "metadata": {}, + "source": [ + "### Functions to make plots" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bd1f16c9", + "metadata": {}, + "outputs": [ + { + "ename": "", + "evalue": "", + "output_type": "error", + "traceback": [ + "\u001b[1;31mnotebook controller is DISPOSED. \n", + "\u001b[1;31mView Jupyter log for further details." + ] + } + ], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python (CV)", + "language": "python", + "name": "sf-manual" + }, + "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.14.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}