{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "is_executing": true, "tags": [ "remove-cell" ] }, "outputs": [], "source": [ "from qutip import *\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import matplotlib.transforms as transforms\n", "import matplotlib.animation as animation\n", "from IPython.display import HTML" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "is_executing": true, "tags": [ "remove-cell" ] }, "outputs": [], "source": [ "plt.rcParams['figure.dpi'] = 150" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Wigner functions of the Quantum Harmonic Oscillator\n", "\n", "In this notebook, you will explore the Wigner functions (Wigner quasiprobability distributions) of various states of the quantum Harmonic oscillator. \n", "\n", "As discussed in the first half of the lecture, the Wigner function is defined as:\n", "\n", "$$\n", "W(x,p) = \\frac{1}{\\pi \\hbar}\n", "\\int_{-\\infty}^{\\infty}\n", "\\psi^*(x+y)\\psi(x-y)\n", "e^{2i p y / \\hbar}\n", "dy\n", "$$\n", "\n", "Here we will calculate the Wigner function of quantum states built from the eigenstates of the Harmonic oscillator. The stationary states of the harmonic oscillator potential defined by the eigenstates of the Hamiltonian:\n", "\n", "$$\n", "\\hat H = \\frac{1}{2} m \\omega^2 \\hat x^2 + \\frac{1}{2m} \\hat p^2\n", "$$\n", "\n", "They are labelled by a number $n \\ge 0$, with wavefunctions $\\psi_n(x)$, and are typically written as $|n\\rangle$, with $n$ an integer, in Dirac notation. \n", "\n", "$$\n", "\\psi_n(x) = \\frac{1}{\\sqrt{2^n\\,n!}} \\cdot \\left(\\frac{m\\omega}{\\pi \\hbar}\\right)^{1/4} \\cdot e^{\n", "- \\frac{m\\omega x^2}{2 \\hbar}} \\cdot H_n\\left(\\sqrt{\\frac{m\\omega}{\\hbar}} x \\right), \\qquad n = 0,1,2,\\ldots. \n", "$$\n", "\n", "where $H_n(x)$ are the Hermite polynomials. \n", "\n", "Another name for a Harmonic oscillator eigenstate is a Fock state: $|1\\rangle$ is a $n=1$ photon Fock state, $|10\\rangle$ is an $n=10$ photon Fock state. \n", "\n", "For this notebook, we will use the python software package QuTiP (QUantum Toolbox In Python) to plot the Wigner functions of different states of the quantum harmonic oscillator. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "hide_input": false, "is_executing": true, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "def plot_wigner_psi_phi(psi,alpha_max = 7.5):\n", " fig = plt.figure(figsize=(9,9))\n", "\n", " widths = [6,3]\n", " heights = [6,3]\n", " spec = fig.add_gridspec(ncols=2, nrows=2, width_ratios=widths,\n", " height_ratios=heights)\n", "\n", " x = np.linspace(-alpha_max,alpha_max,200)\n", " wig = wigner(psi,x,x)\n", " psi_x = np.sum(wig,axis=0)\n", " psi_p = np.sum(wig,axis=1)\n", "\n", "\n", " ax = fig.add_subplot(spec[0,0])\n", " plot_wigner(psi,x,x,fig=fig,ax=ax)\n", " ax = fig.add_subplot(spec[0,1])\n", " base = plt.gca().transData\n", " rot = transforms.Affine2D().rotate_deg(90)\n", " ax.plot(x,-psi_p, transform = rot+base)\n", " ax.set_xticks([])\n", " ax.set_ylim(-alpha_max,alpha_max)\n", " \n", " ax = fig.add_subplot(spec[1,0])\n", " ax.plot(x,psi_x)\n", " ax.set_yticks([]);\n", " ax.set_xlim(-alpha_max,alpha_max)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The Harmonic Oscillator Ground State\n", "\n", "Here below, we will create a \"quantum state\" (actually, something in QuTiP called a \"quantum object\") and make a plot of the Wigner function.\n", "\n", "QuTiP works by storing numerical representations of quantum states in a discrete (countable) Hilbert space of finite size: quantum states are represented by vectors. \n", "\n", "In the code below, the variable `N` is used to tell QuTiP how big the Hilbert space should be. In general, picking the size of the Hilbert space is one of the trickiest things when performing numerical simulations: if you pick one that is too small, your result will be physically incorrect, while if you pick a Hilbert space that is too big, then your computation will become exponentially slower. Fortunately, QuTiP makes good use of efficient sparse matrix routines from pre-compiled libraries.\n", "\n", "In this case, since we are working only with the ground state, we will pick `N = 1` to minimize the computation time.\n", "\n", "Qutip contains many handy routines for making quantum states: one example is the function `fock()` that returns a Fock state. Of course, Fock states are not so hard (in a vector, they are just all zeros and one `1`), but it's handy anyway. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-02-12T14:24:54.565144900Z", "start_time": "2024-02-12T14:24:53.994144900Z" }, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "N = 2\n", "psi = fock(N,0)\n", "plot_wigner_psi_phi(psi)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A coherent state\n", "\n", "Qutip can also automatically create coherent states with a given $\\alpha$ using the function `coherent()`, which is handy since the formula for a coherent is quite messy:\n", "\n", "$$\n", "|\\alpha\\rangle = e^{-\\frac{1}{2} |\\alpha|^2} \\sum_{n=0}^\\infty\\frac{\\alpha^n}{\\sqrt{n!}} |n\\rangle\n", "$$\n", "\n", "Technically, the sum extends to infinity: to create any true coherent state (aside from $\\alpha = 0$), we would need an infinite dimensional Hilbert space. However, as we know, the square of the coefficients follows a Poisson distribution, so for a small enough coherent state, one can \"truncate\" the Hilbert space without resulting in too much error. \n", "\n", "Here, we will work with $\\alpha = 3$, corresponding to an average photon number of $\\langle n \\rangle = |\\alpha|^2 = 9$. Truncating the Hilbert space at $N = 30$ gives a reasonable looking coherent state.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-02-12T14:24:55.539190900Z", "start_time": "2024-02-12T14:24:54.560147400Z" }, "lines_to_next_cell": 2, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "N = 30\n", "psi = coherent(N,3)\n", "plot_wigner_psi_phi(psi,alpha_max=7)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*(You may asky: why does specifying a coherent state with $\\alpha = 3$ seems to produce a plot of a state with $\\alpha = 4.5$? Possibly a bug in the software? Yes, sort of: it turns out that qutip scales $\\alpha$ by a factor of $\\sqrt{2}$ for an unknown reason when calculating Wigner functions...I've submitted a report to the issue tracker system...)*" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "source": [ "## Change phase of $\\alpha$\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-02-12T14:25:31.833865600Z", "start_time": "2024-02-12T14:24:55.568189800Z" }, "collapsed": false, "jupyter": { "outputs_hidden": false }, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "N = 30\n", "plt.rcParams[\"animation.html\"] = \"jshtml\"\n", "plt.rcParams['animation.embed_limit'] = 100\n", "\n", "fig, ax = plt.subplots(1, 1, figsize=(4, 4))\n", "\n", "def update(phi):\n", " ax.cla()\n", " psi = coherent(N, 3 * np.exp(1j * phi))\n", " plot_wigner(psi,fig=fig,ax=ax)\n", " \n", "# Create the animation\n", "# FuncAnimation iterates over these frames, passing each frame value to the provided update function as an argument\n", "phi_values = np.linspace(0, 2*np.pi, 100)\n", "anim = animation.FuncAnimation(fig, update, frames=phi_values)\n", "plt.close(fig) # Close the figure to prevent it from displaying inline\n", "\n", "HTML(anim.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What if I add them (cat states)?\n", "\n", "This one is particularly interesting!\n", "\n", "We know the Shroedinger equation is linear: the sum of two states is equal to the sum of two wavefunction.\n", "\n", "What about the Wigner function? If I add two states, what does the Wigner function of the *sum* of the two state look like? \n", "\n", "Let's try it with these two:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-02-12T14:25:32.696863600Z", "start_time": "2024-02-12T14:25:31.774852400Z" }, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "N = 30\n", "fig,ax = plt.subplots(1,2,figsize=(12,6))\n", "plot_wigner(coherent(N,3j), fig=fig, ax=ax[0])\n", "plot_wigner(coherent(N,-3j), fig=fig, ax=ax[1]);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Add $\\psi$: do the Wigner functions just add? \n", "\n", "*(Note convenient feature of adding quantum objects in QuTiP: it automatically normalizes them for you!)*" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-02-12T14:25:33.506856400Z", "start_time": "2024-02-12T14:25:32.698860700Z" }, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "N = 30\n", "psi = coherent(N,-3j) + coherent(N,3j)\n", "plot_wigner_psi_phi(psi,alpha_max=7)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nearly! And the difference is super-interesting, and has a particularly nice physical interpretation:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "start_time": "2024-02-12T14:25:33.530859700Z" }, "collapsed": false, "is_executing": true, "jupyter": { "outputs_hidden": false }, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "N = 30\n", "fig, ax = plt.subplots(1, 1, figsize=(4, 4))\n", "\n", "def update(phi):\n", " ax.cla()\n", " psi = coherent(N,3*np.exp(1j*phi)) + coherent(N,-3*np.exp(1j*phi))\n", " plot_wigner(psi,fig=fig,ax=ax)\n", "\n", "phi_values = np.linspace(0, 2*np.pi, 100)\n", "anim = animation.FuncAnimation(fig, update, frames=phi_values)\n", "plt.close(fig) # Close the figure to prevent it from displaying inline\n", "\n", "HTML(anim.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*(Foreshadowing to mixed states in future lectures: you can also make Wigner functions for density matrices, and the mixed state corresponding to two coherent state density matrices misses the cool stuff above, and the Wigner function is just the sum of the two Wigner functions...)*\n", "\n", "## Simpler sum: $|0\\rangle + |1\\rangle$\n", "We can also look at a simpler example: just adding $|0\\rangle$ and $|1\\rangle$:\n", "\n", "$$\n", "|\\psi\\rangle = \\frac{1}{\\sqrt{2}}\n", "\\left(\n", "|0\\rangle + e^{i \\phi}|1\\rangle\n", "\\right)\n", "$$\n", "\n", "This is an example you have likely done yourself in a previous quantum course, and the answer is that the particle now has a non-zero $\\langle x \\rangle$ and / or $\\langle p \\rangle$, depending on the phase. \n", "\n", "In fact, if you look at the wavefunction $|\\psi(t=0)\\rangle = \\frac{1}{\\sqrt{2}}(|0\\rangle + |1\\rangle)$, the wavefunction at a later time $t$ will be given by:\n", "\n", "$$\n", "|\\psi(t)\\rangle = \n", "\\frac{e^{-i\\omega t/2}}{\\sqrt{2}} \n", "\\left(\n", "|0\\rangle + e^{-i \\omega t}|1\\rangle \n", "\\right)\n", "$$\n", "\n", "The relative phase $\\phi$ between the two states evolves linearly in time: $\\phi = -\\omega t$. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "is_executing": true, "jupyter": { "outputs_hidden": false }, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1, figsize=(4, 4))\n", "\n", "def update(phi):\n", " ax.cla()\n", " N = 2\n", " psi0 = fock(N,0)\n", " psi1 = fock(N,1)\n", " psi = psi0 + np.exp(1j*phi)*psi1\n", " plot_wigner(psi,fig=fig,ax=ax)\n", "\n", "phi_values = np.linspace(0, 2*np.pi, 100)\n", "anim = animation.FuncAnimation(fig, update, frames=phi_values)\n", "plt.close(fig) # Close the figure to prevent it from displaying inline\n", "\n", "HTML(anim.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you look at $|\\Psi(x,t)|^2$, it should look familiar from what you expect!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Other \"Fock\" states ( = photon number states = stationary states)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Zero photons" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "is_executing": true, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "N = 30\n", "psi = fock(N,0)\n", "plot_wigner_psi_phi(psi)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One photon" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "is_executing": true, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "N = 30\n", "psi = fock(N,1)\n", "plot_wigner_psi_phi(psi)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Take a look at states up to $n=30$:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "is_executing": true, "jupyter": { "outputs_hidden": false }, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "N = 30\n", "fig, ax = plt.subplots(1, 1, figsize=(4, 4))\n", "\n", "def update(n=0):\n", " ax.cla()\n", " psi = fock(N,n)\n", " plot_wigner(psi,fig=fig,ax=ax)\n", "\n", "anim = animation.FuncAnimation(fig, update, frames=range(30))\n", "plt.close(fig) # Close the figure to prevent it from displaying inline\n", "\n", "HTML(anim.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## \"Photon shot noise\"\n", "\n", "If we go back to the coherent states above, what is particularly strange is that the quantum noise of the coherent state in the Wigner plot is always the same: it is always a Heisenberg-limited ball, independent of it's amplitude. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "is_executing": true, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "N = 30\n", "psi = coherent(N,3)\n", "plot_wigner_psi_phi(psi)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "But wait a minute: didn't we learn that coherent states have \"shot noise\" that increases in amplitude when $\\alpha$ get larger. Why don't I see any of this increasing noise in the Wigner plot above? \n", "\n", "The answer is that the \"shot noise\" is not the \"quantum fluctuation\" noise we see above. The quantum fluctuations of the momentum and position (the \"quadratures\") of the coherent state is indeed independent of its amplitude, and in fact always Heisenberg limited. \n", "\n", "The important difference is that if we measure \"photon number\", then we project onto a Fock state. And, as we know, we need a Poissionian distribution of Fock states to construct a coherent state. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "is_executing": true, "jupyter": { "outputs_hidden": false }, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "N = 30\n", "fig, ax = plt.subplots(1,2,figsize=(9, 3.75))\n", "\n", "def update(n):\n", " psi = coherent(N,3)\n", " psi.get_data()[n:] = 0\n", " psi.unit()\n", "\n", " plot_fock_distribution(psi,fig=fig,ax=ax[0],unit_y_range=False)\n", " plot_wigner(psi, fig=fig, ax=ax[1])\n", "\n", "n_values = list(range(1, 29))\n", "anim = animation.FuncAnimation(fig, update, frames=n_values)\n", "\n", "plt.close(fig) # Close the figure to prevent it from displaying inline\n", "\n", "HTML(anim.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Squeezed states\n", "\n", "In the end, if you are interested in sensing a quadrature of an oscillator (for example, if you are working for LIGO trying to detect gravitational waves), you might want to be able to detect signals that are smaller than $x_{zpf}$ in position. \n", "\n", "At first, it might seem that this is impossible: quantum mechanics doesn't allow it! The ground state and coherent state are Heisenberg limited! \n", "\n", "And you would be correct if you restricted yourself to the ground state or a coherent state. (And don't even bother with any of the higher Fock states, the have only bigger $\\sigma_x$ and are not even Heisenberg limited!)\n", "\n", "Is everything hopeless? \n", "\n", "It turns out not: Heisenberg only limits the product $\\sigma_x \\sigma_p > \\hbar / 2$. I can actually get lower uncertainty in my position measurement as long as I don't mind sacrificing some uncertainty in my momentum: in this case, the \"spherical quantum fluctuation ball\" of my coherent state gets \"squeezed\" in the vertical direction, becoming a cigar shape. These will actually be implemented in the next generation of LIGO to enhance the sensitivity of the interferometer for detecting gravitational waves! \n", "\n", "You can find lots of information about squeezed states on wikipedia, and you can also easily generate them in QuTiP using the `squeeze` operator, which squeezes a state along a particular quadrature direction. Squeezed light has many applications, including the sensing one mentioned above, but also for reducing noise in photon counting experiments, and also in quantum communication. \n", "\n", "Squeezing is so importanat it even has two wikipedia pages :\n", "\n", "[Squeezed coherent state](https://en.wikipedia.org/wiki/Squeezed_coherent_state)\n", "\n", "[Squeezed states of light](https://en.wikipedia.org/wiki/Squeezed_states_of_light)\n", "\n", "In the code below, We will first create a coherent state in a slightly different way by \"displacing\" the ground state:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "is_executing": true, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "N = 30\n", "psi = displace(N,3) * basis(N, 0) \n", "plot_wigner_psi_phi(psi)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And we can also add a \"squeezing\" operation to produce a squeezed state:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "is_executing": true, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "N = 50\n", "psi = displace(N,3) * squeeze(N, 1) * fock(N,0)\n", "plot_wigner_psi_phi(psi)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The second parameter in the function is the squeezing parameter: the larger it is, the more squeezed the state becomes. You can play with it to see how it works. \n", "\n", "(Be careful: for high squeezing you may need a bigger Hilbert space...)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "is_executing": true, "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "N = 50\n", "a = destroy(N)\n", "x = a+a.dag()\n", "p = -1j*(a-a.dag())\n", "psi_coh = fock(N,0)\n", "psi_sq = squeeze(N, 1) * fock(N,0)\n", "print(expect(x**2,psi_coh)-expect(x,psi_coh)**2)\n", "print(expect(p**2,psi_coh)-expect(p,psi_coh)**2)\n", "print(expect(x**2,psi_sq)-expect(x,psi_sq)**2)\n", "print(expect(p**2,psi_sq)-expect(p,psi_sq)**2)\n", "print((expect(x**2,psi_sq)-expect(x,psi_sq)**2)*(expect(p**2,psi_sq)-expect(p,psi_sq)**2))" ] } ], "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.12.5" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }