{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "### Calculations to accompany Section II.B of:\n", "## A computational introduction to quantum statistics using harmonically trapped particle\n", "\n", "#### Martin Ligare, Department of Physics & Astronomy, Bucknell Univeristy\n", "\n", " Abstract\n", " In a 1997 paper Moore and Schroeder argued that the development of student understanding\n", "of thermal physics could be enhanced by computational exercises that highlight the link between\n", "the statistical definition of entropy and the second law of thermodynamics [Am. J. Phys. 65, 26\n", "(1997)]. I introduce examples of similar computational exercises for systems in which the quantum\n", "statistics of identical particles plays an important role. I treat isolated systems of small numbers of\n", "particles confined in a common harmonic potential, and use a computer to enumerate all possible\n", "occupation-number configurations and multiplicities. The examples illustrate the effect of quantum\n", "statistics on the sharing of energy between weakly interacting subsystems, as well as the distribution\n", "of energy within subsystems. The examples also highlight the onset of Bose-Einstein condensation\n", "in small systems.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "Section II.B first considers $N_B = 20$ bosons interacting with $N_D=20$ distinguishable particles in a 1-D harmonic trap, with $E=40\\epsilon$ units of energy, where \n", "$\\epsilon = \\hbar\\omega$. The calculations are analogous to those done in Moore and Schroeder's 1997 paper (and the Bucknell PHYS 211 entropy lab), except that the mutiplicities of the macrostates for the bosons is different from that for distinguishable particles. Similar calculations are then done for $N_B=20$ bosons interacting with $N_F=20$ fermions. The energy is not divided equally between the \n", "systems.\n", "\n", "As discussed in the paper, the multipicities for distinguishable particles is\n", "\n", "$$\n", "\\Omega^{\\textrm{D})}(N,q) = \\binom{q + N - 1}{q} = \\frac{(q+N-1)!}{q!(N-1)!},\n", "$$\n", "\n", "while that for bosons is\n", "\n", "$$\n", "\\Omega^{(\\textrm{B})}(N,q) = p( q\\, |\\, \\mbox{number of parts $\\leq N$}),\n", "$$\n", "\n", "where $p$ is the function returning the number of integer partions of $q$ with the number\n", "of parts less than $N$, and that for fermions is\n", "\n", "\\begin{eqnarray*}\n", "\\Omega^{(\\textrm{F})}(N,q) &=& \\Omega^{(\\textrm{B})}(N,q-q^{(F)}_{T=0})\n", " \\nonumber \\\\\n", " &=& p( q-N(N-1)/2\\, |\\, \\mbox{number of parts $\\leq N$})\n", "\\end{eqnarray*}\n", "\n", "This notebook uses the `partitions` generator from the `utilities.iterables` submodule\n", "of `sympy`. There are some \"home-grown\" functions that are faster. \n", "See posts from Chris Smith at \n", "https://code.activestate.com/recipes/218332-generator-for-integer-partitions/. This is faster thanthe function from `sympy.utilities.iterables`.\n", "\n", "(Another approach is given at\n", "https://math.stackexchange.com/questions/18659/algorithm-for-generating-integer-partitions; See comment from `user8548`. This returns partitions in form of Mathematica, cf. dictionary returned by `sympy` or by\n", "`partitions()` function defined below.)\n", "\n", "Mathematica is also faster at this; I haven't explored Sage.\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from scipy import special\n", "from sympy.utilities.iterables import partitions\n", "\n", "import matplotlib as mpl\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Following is an Ipython magic command that puts figures in notebook.\n", "%matplotlib notebook\n", " \n", "# M.L. modifications of matplotlib defaults\n", "# Changes can also be put in matplotlibrc file, \n", "# or effected using mpl.rcParams[]\n", "mpl.style.use('classic') \n", "plt.rc('figure', figsize = (6, 4.5)) # Reduces overall size of figures\n", "plt.rc('axes', labelsize=16, titlesize=14)\n", "plt.rc('figure', autolayout = True) # Adjusts supblot parameters for new size" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Functions returning mulitplicities of $n$ paticles with $q$ units of energy in 1-D harmonic trap" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def omega_D(q,n):\n", " '''Multiplicity of n distinguishable particles with q units of \n", " energy; Eq. (3) of manuscript.'''\n", " return special.binom(q+n-1,q)\n", "\n", "def omega_B(q,n):\n", " '''Multiplicity of n bosons with q units of energy; \n", " Eq. (6) of manuscript'''\n", " count = 0\n", " for p in partitions(q,n):\n", " count += 1\n", " return count\n", "\n", "def omega_F(q,n):\n", " '''Mulitplicity of n fermions with q units of energy;\n", " Eq. (7) of manuscript'''\n", " eF = n*(n-1)/2 # Fermi energy\n", " return omega_B(int(q-eF), n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Calculate multiplicities and entropies" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 4min 41s, sys: 216 ms, total: 4min 41s\n", "Wall time: 4min 41s\n" ] } ], "source": [ "%%time\n", "n_particles = 20\n", "energy_total = 80\n", "eF = int(n_particles*(n_particles -1)/2) # Fermi energy\n", "\n", "omega_D_vals = np.zeros(energy_total+1)\n", "omega_B_vals = np.zeros(energy_total+1)\n", "omega_F_vals = np.zeros(energy_total+1)\n", "\n", "for q in range(energy_total+1):\n", " omega_D_vals[q] = omega_D(q, n_particles)\n", " omega_B_vals[q] = omega_B(q, n_particles)\n", " omega_F_vals[q] = omega_F(q + eF, n_particles)\n", " #print(i) # progress monitor\n", " \n", "s_D = np.log(omega_D_vals)\n", "s_B = np.log(omega_B_vals)\n", "s_F = np.log(omega_F_vals)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "scrolled": false }, "outputs": [ { "data": { "application/javascript": [ "/* Put everything inside the global mpl namespace */\n", "window.mpl = {};\n", "\n", "\n", "mpl.get_websocket_type = function() {\n", " if (typeof(WebSocket) !== 'undefined') {\n", " return WebSocket;\n", " } else if (typeof(MozWebSocket) !== 'undefined') {\n", " return MozWebSocket;\n", " } else {\n", " alert('Your browser does not have WebSocket support. ' +\n", " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", " 'Firefox 4 and 5 are also supported but you ' +\n", " 'have to enable WebSockets in about:config.');\n", " };\n", "}\n", "\n", "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", " this.id = figure_id;\n", "\n", " this.ws = websocket;\n", "\n", " this.supports_binary = (this.ws.binaryType != undefined);\n", "\n", " if (!this.supports_binary) {\n", " var warnings = document.getElementById(\"mpl-warnings\");\n", " if (warnings) {\n", " warnings.style.display = 'block';\n", " warnings.textContent = (\n", " \"This browser does not support binary websocket messages. \" +\n", " \"Performance may be slow.\");\n", " }\n", " }\n", "\n", " this.imageObj = new Image();\n", "\n", " this.context = undefined;\n", " this.message = undefined;\n", " this.canvas = undefined;\n", " this.rubberband_canvas = undefined;\n", " this.rubberband_context = undefined;\n", " this.format_dropdown = undefined;\n", "\n", " this.image_mode = 'full';\n", "\n", " this.root = $('
');\n", " this._root_extra_style(this.root)\n", " this.root.attr('style', 'display: inline-block');\n", "\n", " $(parent_element).append(this.root);\n", "\n", " this._init_header(this);\n", " this._init_canvas(this);\n", " this._init_toolbar(this);\n", "\n", " var fig = this;\n", "\n", " this.waiting = false;\n", "\n", " this.ws.onopen = function () {\n", " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", " fig.send_message(\"send_image_mode\", {});\n", " if (mpl.ratio != 1) {\n", " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", " }\n", " fig.send_message(\"refresh\", {});\n", " }\n", "\n", " this.imageObj.onload = function() {\n", " if (fig.image_mode == 'full') {\n", " // Full images could contain transparency (where diff images\n", " // almost always do), so we need to clear the canvas so that\n", " // there is no ghosting.\n", " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", " }\n", " fig.context.drawImage(fig.imageObj, 0, 0);\n", " };\n", "\n", " this.imageObj.onunload = function() {\n", " fig.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", "\n", " this.ondownload = ondownload;\n", "}\n", "\n", "mpl.figure.prototype._init_header = function() {\n", " var titlebar = $(\n", " '');\n", " var titletext = $(\n", " '');\n", " titlebar.append(titletext)\n", " this.root.append(titlebar);\n", " this.header = titletext[0];\n", "}\n", "\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "\n", "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "mpl.figure.prototype._init_canvas = function() {\n", " var fig = this;\n", "\n", " var canvas_div = $('');\n", "\n", " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", "\n", " function canvas_keyboard_event(event) {\n", " return fig.key_event(event, event['data']);\n", " }\n", "\n", " canvas_div.keydown('key_press', canvas_keyboard_event);\n", " canvas_div.keyup('key_release', canvas_keyboard_event);\n", " this.canvas_div = canvas_div\n", " this._canvas_extra_style(canvas_div)\n", " this.root.append(canvas_div);\n", "\n", " var canvas = $('');\n", " canvas.addClass('mpl-canvas');\n", " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", "\n", " this.canvas = canvas[0];\n", " this.context = canvas[0].getContext(\"2d\");\n", "\n", " var backingStore = this.context.backingStorePixelRatio ||\n", "\tthis.context.webkitBackingStorePixelRatio ||\n", "\tthis.context.mozBackingStorePixelRatio ||\n", "\tthis.context.msBackingStorePixelRatio ||\n", "\tthis.context.oBackingStorePixelRatio ||\n", "\tthis.context.backingStorePixelRatio || 1;\n", "\n", " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", "\n", " var rubberband = $('');\n", " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", "\n", " var pass_mouse_events = true;\n", "\n", " canvas_div.resizable({\n", " start: function(event, ui) {\n", " pass_mouse_events = false;\n", " },\n", " resize: function(event, ui) {\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " stop: function(event, ui) {\n", " pass_mouse_events = true;\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " });\n", "\n", " function mouse_event_fn(event) {\n", " if (pass_mouse_events)\n", " return fig.mouse_event(event, event['data']);\n", " }\n", "\n", " rubberband.mousedown('button_press', mouse_event_fn);\n", " rubberband.mouseup('button_release', mouse_event_fn);\n", " // Throttle sequential mouse events to 1 every 20ms.\n", " rubberband.mousemove('motion_notify', mouse_event_fn);\n", "\n", " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", "\n", " canvas_div.on(\"wheel\", function (event) {\n", " event = event.originalEvent;\n", " event['data'] = 'scroll'\n", " if (event.deltaY < 0) {\n", " event.step = 1;\n", " } else {\n", " event.step = -1;\n", " }\n", " mouse_event_fn(event);\n", " });\n", "\n", " canvas_div.append(canvas);\n", " canvas_div.append(rubberband);\n", "\n", " this.rubberband = rubberband;\n", " this.rubberband_canvas = rubberband[0];\n", " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", " this.rubberband_context.strokeStyle = \"#000000\";\n", "\n", " this._resize_canvas = function(width, height) {\n", " // Keep the size of the canvas, canvas container, and rubber band\n", " // canvas in synch.\n", " canvas_div.css('width', width)\n", " canvas_div.css('height', height)\n", "\n", " canvas.attr('width', width * mpl.ratio);\n", " canvas.attr('height', height * mpl.ratio);\n", " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", "\n", " rubberband.attr('width', width);\n", " rubberband.attr('height', height);\n", " }\n", "\n", " // Set the figure to an initial 600x600px, this will subsequently be updated\n", " // upon first draw.\n", " this._resize_canvas(600, 600);\n", "\n", " // Disable right mouse context menu.\n", " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", " return false;\n", " });\n", "\n", " function set_focus () {\n", " canvas.focus();\n", " canvas_div.focus();\n", " }\n", "\n", " window.setTimeout(set_focus, 100);\n", "}\n", "\n", "mpl.figure.prototype._init_toolbar = function() {\n", " var fig = this;\n", "\n", " var nav_element = $('');\n", " nav_element.attr('style', 'width: 100%');\n", " this.root.append(nav_element);\n", "\n", " // Define a callback function for later on.\n", " function toolbar_event(event) {\n", " return fig.toolbar_button_onclick(event['data']);\n", " }\n", " function toolbar_mouse_event(event) {\n", " return fig.toolbar_button_onmouseover(event['data']);\n", " }\n", "\n", " for(var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " // put a spacer in here.\n", " continue;\n", " }\n", " var button = $('');\n", " button.addClass('ui-button ui-widget ui-state-default ui-corner-all ' +\n", " 'ui-button-icon-only');\n", " button.attr('role', 'button');\n", " button.attr('aria-disabled', 'false');\n", " button.click(method_name, toolbar_event);\n", " button.mouseover(tooltip, toolbar_mouse_event);\n", "\n", " var icon_img = $('');\n", " icon_img.addClass('ui-button-icon-primary ui-icon');\n", " icon_img.addClass(image);\n", " icon_img.addClass('ui-corner-all');\n", "\n", " var tooltip_span = $('');\n", " tooltip_span.addClass('ui-button-text');\n", " tooltip_span.html(tooltip);\n", "\n", " button.append(icon_img);\n", " button.append(tooltip_span);\n", "\n", " nav_element.append(button);\n", " }\n", "\n", " var fmt_picker_span = $('');\n", "\n", " var fmt_picker = $('');\n", " fmt_picker.addClass('mpl-toolbar-option ui-widget ui-widget-content');\n", " fmt_picker_span.append(fmt_picker);\n", " nav_element.append(fmt_picker_span);\n", " this.format_dropdown = fmt_picker[0];\n", "\n", " for (var ind in mpl.extensions) {\n", " var fmt = mpl.extensions[ind];\n", " var option = $(\n", " '', {selected: fmt === mpl.default_extension}).html(fmt);\n", " fmt_picker.append(option);\n", " }\n", "\n", " // Add hover states to the ui-buttons\n", " $( \".ui-button\" ).hover(\n", " function() { $(this).addClass(\"ui-state-hover\");},\n", " function() { $(this).removeClass(\"ui-state-hover\");}\n", " );\n", "\n", " var status_bar = $('