2c:null 2d:["$","title",null,{"children":"02 - Quantum particle tunnelling through a barrier | Matthew Lilley | Observable"}] 32:[] 30:{"id":"a609eaf88f13e77e","slug":"quantum-particle-tunnelling-through-a-barrier","title":"02 - Quantum particle tunnelling through a barrier","thumbnail":"4f2145a9deed49b327d1f900a96ce5d28d15052058be9f2a8df39961e171a8bb","access_level":"public","create_time":"$D2020-10-27T15:59:09.966Z","update_time":"$D2024-11-13T15:20:37.823Z","view_time":null,"first_public_time":"$D2020-12-02T22:43:44.182Z","fork":false,"stars":3,"starred":false,"owner":{"type":"team","id":"4a28ea9ed826d0a0","login":"mklilley","name":"Matthew Lilley","avatar_url":"https://avatars-next.observableusercontent.com/5b842c5b0f775dd6c0697d886b371e2b0ef06d5a7fdd2064fec88b6200598484","tier":"starter_2024"},"version":1602,"publish_version":null,"role":null,"authored":false,"authors":[{"type":"individual","id":"dcd27f799cb1f0fc","login":"mklilley","name":"Matthew Lilley","avatar_url":"https://avatars-next.observableusercontent.com/5b842c5b0f775dd6c0697d886b371e2b0ef06d5a7fdd2064fec88b6200598484","tier":null}],"members":[{"type":"individual","id":"dcd27f799cb1f0fc","login":"mklilley","name":"Matthew Lilley","avatar_url":"https://avatars-next.observableusercontent.com/5b842c5b0f775dd6c0697d886b371e2b0ef06d5a7fdd2064fec88b6200598484","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":true,"trashed":false,"subscribed":null,"forks":2,"collection":{"id":"61c3b2fecc25e94b","slug":"quantum-physics","title":"Quantum physics - Particle in a box","owner":{"id":"4a28ea9ed826d0a0","login":"mklilley","name":"Matthew Lilley"}},"type":"notebook","body":{"title":"02 - Quantum particle tunnelling through a barrier","cells":[{"id":0,"value":"md`# 02 - Quantum particle tunnelling through a barrier`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":802,"value":"md`\n## What?\nIn this second notebook of the series on quantum particles inside \"boxes\", we are once again solving the time-dependent [Schrödinger equation](http://hyperphysics.phy-astr.gsu.edu/hbase/quantum/schrcn.html#c1) in 1-D, i.e.\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":804,"value":"md`\n${tex.block`i\\hbar\\frac{\\partial\\psi}{\\partial t} = -\\frac{\\hbar^2}{2m}\\frac{\\partial^2\\psi}{\\partial x^2} + V\\psi`}\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":809,"value":"md` where we'll be setting ${tex`\\hbar=1`} and ${tex`m=1`} to visually simplify things. Everything is the same as the [first notebook](https://observablehq.com/@mklilley/quantum-particle-moving-in-a-harmonic-potential?collection=@mklilley/quantum-physics) apart from the shape of the potential which is no longer harmonic but instead is a [square potential well](https://en.wikipedia.org/wiki/Finite_potential_well).\n\nThe particle we are going to simulate can be thought of as e.g. an alpha particle trapped inside a nuclues by the nuclear force.\n\nWe will initially consider the particle to have zero average momentum\n\n\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1224,"value":"md`\n\n**p0** = ${viewof p0.bind(html``)}\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1217,"value":"md`\nThe particle should then be trapped by the potential barrier, but as you run the simulation you'll see that the particle can actually find its way out of the potential in a process known as [quantum tunnelling](https://phys.libretexts.org/Bookshelves/University_Physics/Book%3A_University_Physics_%28OpenStax%29/Map%3A_University_Physics_III_-_Optics_and_Modern_Physics_%28OpenStax%29/07%3A_Quantum_Mechanics/7.07%3A_Quantum_Tunneling_of_Particles_through_Potential_Barriers).\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":36,"value":"md`Try it out for yourself.\n\nRun simulation ${viewof simulate.bind(html``)}\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":40,"value":"graph = plot([x,psi0_squared,V_plot],[ \"x\", \"|Psi|^2\",\"V_norm\"]) // Initialise the simulation plot","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1249,"value":"{\n let button\n \n if(simulate){\n button = html``;\n }\n else{\n button = html``;\n }\n\n button.onclick = () => {\n reset()\n }\n\n return button;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":963,"value":"md`Below, we keep track of the probability of finding the particle inside the well and can see that (at least by eye) it decays exponentially in time.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":311,"value":"graph_P_well = plot([nj.array([[0,0]])],[ \"time\", \"P_well\"]) // Initialise the simulation plot","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1531,"value":"md`## Why?\nThe exponential behaviour we see above is key to understanding the [half lives](http://hyperphysics.phy-astr.gsu.edu/hbase/Nuclear/halfli.html#c1) of radioactive nuclei and was first used by [Gamow](http://web.ihep.su/dbserv/compas/src/gamow28/eng.pdf) to understand [alpha emission](http://hyperphysics.phy-astr.gsu.edu/hbase/Nuclear/radact.html#c2). It's also intimately conneted the process of nuclear fusion in which two positive nuclei must tunnel through a [Coulomb barrier](http://hyperphysics.phy-astr.gsu.edu/hbase/NucEne/coubar.html) in order to reach proximities at which the nuclear force can bind them together.\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1527,"value":"md`## How?\nTo understand the nuts and bolts of the algorithm being used here, head back over to the [first notebook](https://observablehq.com/@mklilley/quantum-particle-moving-in-a-harmonic-potential). Everything below is more or less the same, the main differences are the potential, the initial conditions and an option to plot at specific time intervals.\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":901,"value":"md`**Simulation parameters**`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1084,"value":"md`Be aware that this simulation box is much larger than last time. As you change parameters you might need to wait a few seconds for the simulation grid to catch-up...alas JavaScript has its limits.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":50,"value":"md`Time step:\n\n**dt** = ${viewof dt.bind(html``)}`\n\n","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1449,"value":"md`Time steps between each plotting frame:\n\n**n** = ${viewof n.bind(html``)}`\n\n","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":52,"value":"md`\nSize of simulation box: \n\n**xmax** = ${viewof xmax.bind(html``)}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":54,"value":"md`Simulation grid spacing:\n\n**dx** = ${viewof dx.bind(html``)}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":56,"value":"numx = {\n if(dx){\n if(dx !=0){\n return Math.round(2*xmax/dx +1)\n }\n else { // This is in case the input is 0, we don't want the whole notebook to show errors and make the screen jump all over the place\n return 0\n }\n }\n else{ // This is in case the input is empty, we don't want the whole notebook to show errors and make the screen jump all over the place\n return 0\n }\n}\n","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":866,"value":"md`Absorbing boundary strength:\n\n**${tex`\\lambda`}** = ${viewof lambda.bind(html``)}\n
\n\nAbsorbing boundary width:\n\n**${tex`w/x_{max}`}** = ${viewof wOverXmax.bind(html``)}\n`\n\n","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":875,"value":"w = wOverXmax*xmax","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":898,"value":"md`**Build simulation grid**\n\nFirst we need to make a column of ${tex`x`} values`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":68,"value":"x = nj.arange(numx).multiply(dx).subtract(xmax).reshape(numx,1)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1080,"value":"md`From these ${tex`x`} values we can calculate the corresponding ${tex`k`} values for Fourier space as we did last time.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":70,"value":"k = fftfreq(x)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":904,"value":"md`**Potential**\n\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1138,"value":"md`We construct the potential by \"gluing\" together several constant pieces to create a [piecewise function](https://en.wikipedia.org/wiki/Piecewise)`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1096,"value":"md`\n**barrier_width** = ${viewof barrier_width.bind(html``)}`\n\n","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1115,"value":"md`\n**barrier_height** = ${viewof barrier_height.bind(html``)}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1097,"value":"md`\n\n**well_width** = ${viewof well_width.bind(html``)}`\n\n","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1046,"value":"V = {\nlet outside_right = nj.zeros(Math.round((numx-1)/2)) // this should make the potential free region start at 0\nlet barrier = nj.ones(Math.round(barrier_width/dx)).multiply(barrier_height)\nlet well = nj.zeros(Math.round(well_width/dx))\nlet outside_left = nj.ones(numx - Math.round(well_width/dx) - Math.round(barrier_width/dx) - Math.round((numx-1)/2) ).multiply(barrier_height)\n\n\nlet V = nj.concatenate(outside_left,well,barrier,outside_right)\nreturn V.reshape(numx,1)\n\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1515,"value":"md`Later on, we'll need to integrate the probability density inside the well, so it will be convenient to know the index values where it begins and ends.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1255,"value":"well_start = numx - Math.round(well_width/dx) - Math.round(barrier_width/dx) - Math.round((numx-1)/2)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1260,"value":"well_end = well_start + Math.round(well_width/dx)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1213,"value":"md`Let's also create a normalised version of the potential so that we can always plot it alongside the probability density - it will always have a height of 1 on the graph.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1193,"value":"V_plot = V.divide(barrier_height)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":895,"value":"md`**Exponentials**\n\nWe can now create the complex exponentials that allow us to advance the solution as in the \\`\\`psi_dt\\`\\` function - just as we did last time.\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":72,"value":"expK = cexp(k.pow(2).multiply(-dt*(2*Math.PI)**2.0))","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1134,"value":"expV = cexp(V.multiply(-dt/2))","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":419,"value":"expABC = {\n let absorb = nj.tanh(x.subtract(xmax).divide(w)).add(2).subtract(nj.tanh(x.add(xmax).divide(w))).multiply(-lambda*dt)\n \n let expABC = nj.exp(absorb)\n \n expABC = nj.concatenate(expABC, nj.zeros([numx,1]) ) // need to add 2nd column of zeros for the zero imaginary part\n \n return expABC\n\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":922,"value":"md`**Initial conditions**`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1147,"value":"md`We can represent a particle moving with average momentum momentum ${tex`p_0`} localised in space around ${tex`x_0`} by the following [wave packet](http://farside.ph.utexas.edu/teaching/qmech/Quantum/node25.html):\n\n${tex.block`\\psi_0 = \\left(\\frac{2}{\\pi \\sigma^2}\\right)^{1/4} e^{-(x-x_0)^2/\\sigma^2 + ip_0x}`}\n\n\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1332,"value":"md`We'll place the particle in the centre of the potential well so that`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1330,"value":"x0 = - barrier_width - well_width/2","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1336,"value":"psi0 = {\n let sigma2 = 2\n let norm_constant = Math.pow(2 / (Math.PI * sigma2), 0.25);\n let psi0 = nj.concatenate(nj.exp(x.subtract(x0).pow(2).multiply(-1).divide(sigma2)), nj.zeros([numx,1]) )\n psi0 = cmultiply(psi0,cexp(x.multiply(p0)))\n\n return psi0.multiply(norm_constant);\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1380,"value":"md` Now we can create ${tex`|\\psi_0|^2`} for initialising the first plot ` ","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1369,"value":"psi0_squared = cnorm2(psi0)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1398,"value":"md`and ${tex`P_{\\text{well}} = \\int _{\\text{well}} |\\psi_0|^2 dx`} for the second plot` ","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1415,"value":"P_well0 = psi0_squared.slice([well_start,well_end]).sum()*dx","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":925,"value":"md`**Simulation**\n\nFirst, let's set up the variables that we'll need.\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1344,"value":"mutable psi=psi0 // This sets the initial condition for the simulation of psi\n\n// Note, for more info on \"mutable\" see https://observablehq.com/@mbostock/mutable-values","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1374,"value":"md`We create an additional variable psi_plot so that we have the flexibility visualise the simulation every **n** timesteps.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1349,"value":"mutable psi_plot = psi0","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1436,"value":"md`We need to create two arrays to hold evolving values of P_well and also the time for the second plot`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1407,"value":"mutable P_well = nj.array([P_well0]) // this creates an array to store the evolving values of P_well","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1424,"value":"mutable time = nj.array([0]) // this creates an array to store the evolving time for plotting along with P_well","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1433,"value":"mutable i = 0 // this keeps track of the number of the elapsed timesteps","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1463,"value":"md` Now let's bring all the initialisations together in a function that can be called using the reset button.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1462,"value":"reset = function(){\n mutable psi = psi0;\n mutable psi_plot = psi0;\n mutable time = nj.array([0])\n mutable P_well = nj.array([P_well0])\n mutable i =0\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1473,"value":"md`And now for the main simulation loop.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":80,"value":"{ \n if(simulate===true) { // Only simulate when simulation checkbox has been checked\n mutable psi = psi_dt(psi, expV, expK, expABC)\n \n // Update the timestep\n mutable i+=1\n \n // Only update the psi_plot for plotting every n time-steps\n if(i%n==0) {\n \n mutable psi_plot = psi // Note this line will trigger psi_squared for the plots and also for the next step below\n \n // Add new data to the P_well and time arrays by concatenating old array with new value\n mutable P_well = nj.concatenate(mutable P_well,nj.array([psi_squared.slice([well_start,well_end]).sum()*dx])) \n mutable time = nj.concatenate(mutable time, nj.array([i*dt]))\n \n }\n\n\n }\n return \"Indefinite simulation loop\"\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":85,"value":"psi_squared = cnorm2(psi_plot)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1477,"value":"md`The last step is to update the graphs at the top of the page with new values.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":941,"value":"graph.update([x,psi_squared,V_plot]) // Update the simulation plot when new values of psi_plot are available","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1007,"value":"graph_P_well.update([nj.stack([time, P_well], -1)]) // Update the simulation plot when new values of P_well are available","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":5,"value":"md`---`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":4,"value":"md`## Appendix`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1545,"value":"md`### numjs\nImport the numjs library\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":10,"value":"nj = require('https://cdn.jsdelivr.net/gh/nicolaspanel/numjs@0.16/dist/numjs.min.js')","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1562,"value":"md`### Tutorial specific functions\n\nWe'll need to import several functions that we created in the [first tutorial](https://observablehq.com/@mklilley/quantum-particle-moving-in-a-harmonic-potential).\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":679,"value":"import {fftfreq, plot, cmultiply, cexp, cnorm2, psi_dt} from \"@mklilley/quantum-particle-moving-in-a-harmonic-potential\"","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1551,"value":"md`### Observable \nHere is some observable specific code that makes it conveient to update parameters. See the [first tutorial](https://observablehq.com/@mklilley/quantum-particle-moving-in-a-harmonic-potential) for more detail.\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":28,"value":"import {View} from \"@mbostock/synchronized-views\"","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":44,"value":"viewof simulate = new View(false)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":58,"value":"viewof dx = new View(0.1)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":60,"value":"viewof xmax = new View(30)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":62,"value":"viewof dt = new View(0.01)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":869,"value":"viewof wOverXmax = new View(0.05)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":872,"value":"viewof lambda = new View(20)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1104,"value":"viewof well_width = new View(2)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1106,"value":"viewof barrier_width = new View(1)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1119,"value":"viewof barrier_height = new View(2)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1238,"value":"viewof p0 = new View(0)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":1455,"value":"viewof n = new View(5)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"}],"files":"$W32","stdlib":"1","titling":"auto"}} 31:"$30:body"