diff --git a/contents/quantum_systems/code/julia/energy.jl b/contents/quantum_systems/code/julia/energy.jl new file mode 100644 index 000000000..d045d073c --- /dev/null +++ b/contents/quantum_systems/code/julia/energy.jl @@ -0,0 +1,18 @@ +# We are calculating the energy to check +function calculate_energy(wfc, H_k, H_r, dx) + # Creating momentum and conjugate wavefunctions + wfc_k = fft(wfc_r) + wfc_c = conj(wfc_r) + + # Finding the momentum and real-space energy terms + energy_k = wfc_c.*ifft((H_k) .* wfc_k) + energy_r = wfc_c.* H_r .* wfc_r + + # Integrating over all space + energy_final = 0 + for i = 1:length(energy_k) + energy_final += real(energy_k[i] + energy_r[i]) + end + + return energy_final*dx +end diff --git a/contents/quantum_systems/quantum_systems.md b/contents/quantum_systems/quantum_systems.md index 661db87b0..a53f45c26 100644 --- a/contents/quantum_systems/quantum_systems.md +++ b/contents/quantum_systems/quantum_systems.md @@ -2,7 +2,9 @@ As I am sure you have heard, the quantum world is weird. As you deal with progressively smaller and smaller systems, at some point, it becomes less accurate to describe objects as particles. -Instead, it is better to describe objects as probability waves. +Instead, it is better to describe objects as probability densities. +These densities are easiest to understand in terms of _wavefunctions_, which are complex functions characterizing a quantum system's behaviour. + Again, this is pretty common knowledge; however, there is a distinct lack of readable literature on how to simulate quantum systems, even though there are numerous methods for exactly that! This section will deal with the computation of quantum states with classical machines. Now, I know what you are thinking, "Wait. Why are we simulating quantum systems on classical computers? Why not simulate it with some sort of experiment or with quantum computers?" @@ -13,20 +15,21 @@ A _quantum computer_ is the quantum analog to a classical computer, replacing bi Quantum computers are usually thought of as a way to use quantum mechanics to eventually solve real-world problems with new quantum algorithms. Both Grover's and Shor's algorithms are good examples of cases where quantum computation could greatly change the landscape of modern computation as we know it! -_Quantum Simulators_ on the other hand are quantum systems used to better un -derstand quantum mechanics -Because supercomputers are not great at performing quantum computations, quantum simulators exist as a building block for quantum computation; however, their purpose is not explicity for quantum information theory. -Often times a _universal quantum simulator_ is often called a quantum computer. +_Quantum simulators_ on the other hand are quantum systems used to better understand quantum mechanics. +These will often come in the form of experimental quantum systems that express quantum behaviour and allow us to better understand other areas of quantum systems. +In other words, quantum simulators are general techniques to study quantum systems on quantum hardware; however, quantum computers are quantum hardware used for the explicit purpose of quantum computation with qubits. +Because supercomputers are not great at performing quantum computations, certain quantum simulators exist as a building block for quantum computation. +A _universal quantum simulator_ is often called a quantum computer for this reason. -The truth is that until we have real quantum simulators, simulating quantum systems on classical hardware is as good as we can do. -This section is devoted to all the different methods currently used to solve complex quantum systems, so let's start with the Schrodinger Equation, which has many different fomulations. +The truth is that quantum simulators are hard to make in laboratories, so simulating quantum systems on classical hardware is as good as we can do in most cases. +This section is devoted to all the different methods currently used to solve complex quantum systems, so let's start with the Schrödinger Equation, which has many different formulations. Here is the easiest one to explain: $$ i \hbar \frac{\partial \Psi(\mathbf{r},t)}{\partial t} = \left[-\frac{\hbar^2}{2m} \nabla^2 + V(\mathbf{r},t) \right] \Psi(\mathbf{r},t) $$ -Where $$\Psi(\mathbf{r},t)$$ is your quantum _wavefunction_, $$V(\mathbf{r},t)$$ is a _trapping potential_, $$\nabla^2$$ is a _laplacian_, $$\mathbf{r}$$ is some sort of spatial component, and $$t$$ is time. +Where $$\Psi(\mathbf{r},t)$$ is a quantum wavefunction, $$V(\mathbf{r},t)$$ is a _trapping potential_, $$\nabla^2$$ is a _laplacian_, $$\mathbf{r}$$ is some sort of spatial component, and $$t$$ is time. There is a lot to take in here; however, it's ultimately just some time derivative on the left-hand side and a spatial derivative (with some extra steps) on the right-hand side. In this way, it isn't too different from the diffusion (heat) equation: @@ -35,11 +38,14 @@ $$ $$ where $$D$$ is some positive definite matrix and $$\phi(\mathbf{r},t)$$ is the density (or temperature) of the system. -In fact, this is why one of the most common types of quantum simulation is via _diffusion monte carlo_. -There really isn't that much of a difference between the two systems in terms of classical simulation. +In fact, this is why one of the most common types of quantum simulation is sometimes called _diffusion monte carlo_. +There really isn't that much of a difference between the two systems in terms of how they are simulated on classical hardware... but we are getting ahead of ourselves. +For now, let's talk about how quantum mechanics differs from classical mechanics and how we can use this to our advantage. + +## Probability Density -As a note: quantum mechanics works fundamentally differently than classical mechanics in physics. -The wavefunction is essentially a set of all possible states for an object to be in, where there is some probability for the particle to be found in each state. +Quantum mechanics works fundamentally differently than classical mechanics in physics. +The wavefunction can be thought of as a set of all possible states for an object to be in, where there is some probability for the particle to be found in each state. This means that it is not possible to say that a particle is at a particular location, and instead we often say that it could be at any location with probability, as shown in the _probability density_: $$ @@ -48,8 +54,8 @@ $$ Here, there are 2 things to note: -1. The absolute value squared of a complex parameter $$\Psi(\mathbf{r},t)$$ is a dot product (inner product) between a complext function and it's Hermitian conjugate -2. As you have probably heard, once a wavefunction is observed it collapses onto a single state +1. The absolute value squared of a complex parameter $$\Psi(\mathbf{r},t)$$ is a dot product (inner product) between a complex function and it's Hermitian conjugate. This means the value will always be real, while the wavefunction, itself, might not be. +2. As you have probably heard, once a wavefunction is observed it collapses onto a single state. This can be simply interpreted as absolute knowledge of the particle's location. A probability density doesn't make sense if we know where the particle is! Now, to be clear: the probabilities must all sum to 1, or (more formally): @@ -57,42 +63,187 @@ $$ \int_{-\infty}^{+\infty}|\Psi(\mathbf{r},t)|^2 d\mathbf{r} = 1 $$ -As another note: Just like position space can be parameterized by a position vector $$\textbf{x}$$, wavefunctions can be parameterized by a _wave_vector $$\textbf{k}$$ in frequency space. +This simply means that the probability of finding our quantum particle *somewhere in real space* is 1. +In other words, our particle must exist somewhere in the known universe. + +As another note: Just like position space can be parameterized by a position vector $$\textbf{x}$$, wavefunctions can also be parameterized by a _wave_ vector $$\textbf{k}$$ in frequency space. +Any wavevector $$\textbf{k}$$ has the same units as reciprocal space and is thus analogous to angular frequency $$\omega$$. Often times, the wavevector space is called _momentum_ space, which makes sense when considering the de Broglie formula: $$ -p = \frac{h}{\lambda} +p = \frac{h}{\lambda} = \frac{2 \pi h}{2 \pi \lambda} = \hbar k $$ where $$h$$ is Planck's constant and $$\lambda$$ is the wavelength. -This means that we can ultimately move between real and momentum space by using [Fourier Transforms](../cooley_tukey/cooley_tukey.md), which is incredibly useful in a number of cases! +This means that we can ultimately move between position and momentum space by using [Fourier Transforms](../../algorithms/cooley_tukey/cooley_tukey.md), which is incredibly useful in a number of cases! + +Even though the relation between position and momentum space is an essential cornerstone of understanding modern quantum mechanics, it is difficult to understand at a fundamental level. +Position space and momentum space are related by a Fourier transform; however, the rather hand-wavey argument above might not have been convincing enough and it does not offer any intuitive description of _why_ the Fourier transform comes into this discussion at all. +The easiest way to understand this might be to look at the _Heisenberg uncertainty principle_, which is a fundamental relation between position and momentum space. + +## Heisenberg Uncertainty Principle + +Simply put, the Heisenberg uncertainty principle states that we cannot definitely know both the position and momentum of a quantum particle. +In particular, it says: + +$$ +\sigma_x \sigma_p \geq \frac{\hbar}{2} +$$ + +where $$\hbar$$ is Planck's constant and $$\sigma_q = \sqrt{\frac{1}{N}\sum_{i=1}^{N}(q_i-\mu)^2}$$. +In this case, $$\sigma$$ is the standard deviation, $$\mu$$ is the statistical mean of your distribution, $$N$$ is the number of points sampled, $$q_i$$ is the value for each point $$i$$, and $$q$$ stands for $$r$$ or $$p$$. +Ultimately, this means that if we have a higher precision in position space, we will have a lower precision in momentum space. +The converse is also true: a higher precision in momentum space will lead to a lower precision in position space. + +This makes the most sense if we imagine having a gaussian-like probability density ($$|\Psi(x)|^2$$) in position space, which will provide a gaussian-like density when in momentum space. +Here, we see that if we have a broader distribution in one space, we must have a thinner distribution in the opposite space, as shown here: + +

+ +

+ + +Because the density can be interpreted as "the probability of finding a quantum particle at any provided location in position ($$x_i$$) or momentum ($$k_i$$) space, the interpretation is clear: the more we understand about a particle's position, the less we understand about it's momentum. +This is a powerful statement and should be given some thought. + +To me, the most interesting part of this description is not the physical interpretation, but the fact that this act of transforming between larger and smaller gaussians is precisely what Fourier transforms do! +This further strengthens our argument from before. +Position and momentum space are related by the Fourier transform! + +This is the heart of several algorithms for simulating quantum systems, including the [Split-operator method](../../algorithms/split-operator_method/split-operator_method). + +At least for me, I found this description to be intuitive, but not complete. +There is still something missing from the picture that should be described in more depth, and to do that, we need to dive deeper into the heart of quantum mechanics and into _Hamiltonians_. + +## Hamiltonian + +Here is the Schrödinger equation again: + +$$ +i \hbar \frac{\partial \Psi(\mathbf{r},t)}{\partial t} = \left[-\frac{\hbar^2}{2m} \nabla^2 + V(\mathbf{r},t) \right] \Psi(\mathbf{r},t) +$$ + +We described it in the initial section of this chapter. +For the most part, when we are trying to solve this equation the left-hand side does not change. +It's always $$i \hbar \frac{\partial \Psi(\mathbf{r},t)}{\partial t}$$. +On the other hand, the right-hand side can change a lot depending on the situation. +That is to say that we can easily simplify the Schrödinger equation by using a mathematical formalism known as the _Hamiltonian_. + +To be clear, Hamiltonian mechanics is not a quantum-specific idea. +It's everywhere in statistical physics and is often taught in classical physics courses as an analogue to another notational form known as Lagrangian mechanics. +For simplicity, we will restrict our discussion here to interpreting Hamiltonians physically. +We can basically say that the Hamiltonian is a measure of the energy of our quantum system. +More specifically, we can say that the Hamiltonian is a set of energy _operators_ that act on our wavefunction. + +In the case of a 1D particle in a harmonic trap, we might use the following definitions: + +$$ +\begin{align} +\hat H &= \hat T + \hat V \\ +\hat T &= \frac{p^2}{2m} \\ +\hat V &= \frac{1}{2}\omega x^2 +\end{align} +$$ + +where $$p = -i\hbar \nabla$$ is the _momentum operator_ and $$\omega$$ is the _trapping frequency_ indicating how confined our quantum system will be. +In this case, $$\hat T$$ is an operator that works on our wavefunction in momentum space, while $$\hat V$$ acts in position space. +Both of these are operators. +That is to say that they _operate_ on our quantum system by transforming it in some way. +Ultimately, this means that the operators are not meant to be interpreted on their own without acting on some other object, in this case, the wavefunction $$\Psi(x)$$. -Unfortunately, the interpretation of quantum simulation is rather tricky and is ultimately easier to understand with slightly different notation. +In the end, we can update our Schrödinger equation to be + +$$ +i \hbar \frac{\partial \Psi(\mathbf{r},t)}{\partial t} = \hat H \Psi(\mathbf{r},t) +$$ + +Which is a lot cleaner and more general. +Now, the Schrödinger equation can solve any quantum system so long as it can be written in terms of Hamiltonian mechanics! + +When looking at the operators, it is clear that the $$\hat V$$ operator is in position space. +We can clearly see that it operates as a function of $$x$$. +That said, it is not immediately obvious why the $$\hat T$$ is in momentum space. +This is not an easy question to answer, but it is definitely important and will be covered in more depth when we discuss spectral methods. + +For now, we will blanketly say + +$$ +\frac{\partial f}{\partial x} = \mathcal{F}^{-1}\left( 2\pi i k \mathcal{F}\left( f \right)\right) +$$ + +In other words, we can derive a function by performing a Fourier transform on the function, multiplying by some momentum-space grid, and then inverse-transforming it back. +Because this operation inherently involves a transform into momentum space before transformation, it is a momentum-space operator. + +This is the most intuitive reasoning I can find; however, I am sure there are more intuitive explanations of why the derivatives are always momentum-space operations. +This section will be updated further when we discuss spectral methods, but if you have better descriptions, please let me know! + +## Bra Ket Notation + +Unfortunately, the interpretation of quantum simulation is rather tricky and is sometimes easier to understand with slightly different notation. This notation is called _braket_ notation, where a _ket_ looks like this: $$ \lvert A \rangle $$ -And basically describes $$A$$ as a column vector. +and basically describes $$A$$ as a column vector. The _bra_ represents the Hermitian conjucate of the ket and looks like this: $$ \langle B \rvert $$ -It is often represented as a row vector for $$B$$. Because of this, $$ \langle B \rvert A \rangle $$ represents the inner product of the two vectors and $$ \lvert A \rangle \langle B \rvert $$ represents the outer product. +The ket is often represented as a row vector for $$B$$. +Because of this, $$ \langle B \rvert A \rangle $$ represents the inner product of the two vectors and $$ \lvert A \rangle \langle B \rvert $$ represents the outer product. Now, to this point, the braket notation does not have any particularly quantum-like features; however, it becomes useful when describing actual quantum phenomenon. For example, if we want to indicate the probability of a wavefunction $$\psi$$ collapsing onto state $$\phi$$, we might write: $$\langle \phi \rvert \psi \rangle$$, which is precisely the same as the probability density defined above. +Now that we have a basic understanding of the notation, we should go through several other important quantum mechanical ideas and properties. + +## Eigenstates +As mentioned, the wavefunction $$\Psi(x)$$ is complex and has both real and imaginary parts; however, there are certain states that are eclusively real. +These states are _eigenstates_ of the system, and are often described as the constituent states that make up all other possible wavefunctions. +In other words, + +$$ +\lvert \Psi(x)\rangle = \sum_i c_i \lvert \Psi_i \rangle +$$ + +Where $$c_i$$ is some constant describing _how much_ of a given eigenstate $$i$$ is in the full wavefunction. +As you might expect, all of the $$c_i$$'s should sum to 1. + +## Energy Calculations + +When it comes to quantum systems, there is no quantity more important than energy. +Basically, every eigenstate of the system has a different energy associated with it, and you can find this energy with a simple calculation: + +$$ +E = \langle \Psi \lvert \hat H \lvert \Psi \rangle +$$ + +Which can be done rather trivially in code by finding the conjugate of the wavefunction and multiplying it with another wavefunction after operation in position and momentum space. +This ultimately looks like this: + +{% method %} +{% sample lang="jl" %} +[import, lang:"julia"](code/julia/energy.jl) +{% endmethod %} + +This calculation will be used in many different simulations of quantum systems to check our results. +In the end, many quantum simulations are focused on the _ground_ state, which is the lowest energy state ($$\Psi_0$$); however, sometimes higher energy states are desired. + +## The Future + As we proceed to add new algorithms to simulate quantum systems, I will add more and more notation to this section; however, there are already huge textbooks out there related to understanding and studying quantum systems. We don't want to re-invent the wheel here. -Instead, we want to focus on an area that is often not considered with too much detail -- the algorithms and methods researchers use to ascertain new knowedge about quantum mechanics, like the split-operator method, DMRG, quantum monte carlo, exact diagonalization, and many more. +Instead, we want to focus on an area that is often not considered with too much detail: algorithms and methods researchers use to ascertain new knowedge about quantum mechanics, like the split-operator method, DMRG, quantum monte carlo, exact diagonalization, and many more. -Quantum mechanics is one of those areas of physics that really does push the boundary of human knowledge in a number of different areas and computing is no different. +Quantum mechanics is one of those areas of physics that really does push the boundary of human knowledge in a number of different areas and computing is one of those areas. In fact, [quantum information theory](../quantum_information/quantum_information.md) is currently set to be the next innovation to radically change the landscape of modern computation as we know it! Of course, because of the large-scale effects that this will likely have on the industry, it deserved it's own section. +As always, if there is something that you feel is missing from this section, please feel free to contact me or create an issue on GitHub and we'll get to it as soon as we can! + diff --git a/contents/quantum_systems/res/gaussian.gif b/contents/quantum_systems/res/gaussian.gif new file mode 100644 index 000000000..d6faa20f0 Binary files /dev/null and b/contents/quantum_systems/res/gaussian.gif differ diff --git a/contents/split-operator_method/code/julia/split_op.jl b/contents/split-operator_method/code/julia/split_op.jl index b77bae5f6..7fd771601 100644 --- a/contents/split-operator_method/code/julia/split_op.jl +++ b/contents/split-operator_method/code/julia/split_op.jl @@ -1,7 +1,11 @@ -using Plots -pyplot() +#------------split_op.jl-------------------------------------------------------# +# +# Plotting: to plot individual timesteps, use gnuplot like so: +# p "output00000.dat" u 1:2 w l +# rep "output00000.dat" u 1:3 w l +# +#------------------------------------------------------------------------------# -# struct to hold all parameters for simulation struct Param xmax::Float64 res::Int64 @@ -11,34 +15,48 @@ struct Param x::Vector{Float64} dk::Float64 k::Vector{Float64} + im_time::Bool Param() = new(10.0, 512, 0.05, 1000, 2 * 10.0/512, Vector{Float64}(-10.0 + 10.0/512 : 20.0/512 : 10.0), - pi / 10.0, - Vector{Float64}(vcat(0:512/2 - 1, -512/2 : -1) * pi/10.0)) - Param(xmax::Float64, res::Int64, dt::Float64, timesteps::Int64) = new( + pi / 10.0, + Vector{Float64}(vcat(0:512/2 - 1, -512/2 : -1) * pi/10.0), + false) + Param(xmax::Float64, res::Int64, dt::Float64, timesteps::Int64, + im_val::Bool) = new( xmax, res, dt, timesteps, 2*xmax/res, Vector{Float64}(-xmax+xmax/res:2*xmax/res:xmax), - pi/xmax, Vector{Float64}(vcat(0:res/2-1, -res/2:-1)*pi/xmax) + pi/xmax, Vector{Float64}(vcat(0:res/2-1, -res/2:-1)*pi/(xmax)), + im_val ) end -# struct to hold all operators mutable struct Operators V::Vector{Complex{Float64}} PE::Vector{Complex{Float64}} KE::Vector{Complex{Float64}} wfc::Vector{Complex{Float64}} + + Operators(res) = new(Vector{Complex{Float64}}(res), + Vector{Complex{Float64}}(res), + Vector{Complex{Float64}}(res), + Vector{Complex{Float64}}(res)) end # Function to initialize the wfc and potential function init(par::Param, voffset::Float64, wfcoffset::Float64) - V = 0.5 * (par.x - voffset).^2 - wfc = 3 * exp.(-(par.x - wfcoffset).^2/2) - PE = exp.(-0.5*im*V*par.dt) - KE = exp.(-0.5*im*par.k.^2*par.dt) + opr = Operators(length(par.x)) + opr.V = 0.5 * (par.x - voffset).^2 + opr.wfc = exp.(-(par.x - wfcoffset).^2/2) + if (par.im_time) + opr. = exp.(-0.5*par.k.^2*par.dt) + opr.PE = exp.(-0.5*opr.V*par.dt) + else + opr. = exp.(-im*0.5*par.k.^2*par.dt) + opr.PE = exp.(-im*0.5*opr.V*par.dt) + end - opr = Operators(V, PE, KE, wfc) + return opr end # Function for the split-operator loop @@ -48,11 +66,11 @@ function split_op(par::Param, opr::Operators) # Half-step in real space opr.wfc = opr.wfc .* opr.PE - # fft to phase space + # fft to momentum space opr.wfc = fft(opr.wfc) - # Full step in phase space - opr.wfc = opr.wfc .* opr.KE + # Full step in momentum space + opr.wfc = opr.wfc .* opr. # ifft back opr.wfc = ifft(opr.wfc) @@ -60,20 +78,66 @@ function split_op(par::Param, opr::Operators) # final half-step in real space opr.wfc = opr.wfc .* opr.PE - # plotting density and potential + # density for plotting and potential density = abs2.(opr.wfc) - plot([density, real(opr.V)]) - savefig("density" * string(lpad(i, 5, 0)) * ".png") - println(i) + # renormalizing for imaginary time + if (par.im_time) + renorm_factor = sum(density) * par.dx + + for j = 1:length(opr.wfc) + opr.wfc[j] /= sqrt(renorm_factor) + end + end + + # Outputting data to file. Plotting can also be done in a similar way + # This is set to output exactly 100 files, no matter how many timesteps + if ((i-1) % div(par.timesteps, 100) == 0) + outfile = open("output" *string(lpad(i-1, 5, 0))* ".dat","w") + + # Outputting for gnuplot. Any plotter will do. + for j = 1:length(density) + write(outfile, string(par.x[j]) * "\t" + * string(density[j]) * "\t" + * string(real(opr.V[j])) * "\n") + end + + close(outfile) + println("Outputting step: ", i) + end end end +# We are calculating the energy to check +function calculate_energy(par, opr) + # Creating real, momentum, and conjugate wavefunctions + wfc_r = opr.wfc + wfc_k = fft(wfc_r) + wfc_c = conj(wfc_r) + + # Finding the momentum and real-space energy terms + energy_k = 0.5*wfc_c.*ifft((par.k.^2) .* wfc_k) + energy_r = wfc_c.*opr.V .* wfc_r + + # Integrating over all space + energy_final = 0 + for i = 1:length(energy_k) + energy_final += real(energy_k[i] + energy_r[i]) + end + + return energy_final*par.dx +end + # main function function main() - par = Param(10.0, 512, 0.05, 1000) - opr = init(par, 0.0, 1.0) + par = Param(5.0, 256, 0.05, 100, true) + + # Starting wavefunction slightly offset so we can see it change + opr = init(par, 0.0, -1.00) split_op(par, opr) + + energy = calculate_energy(par, opr) + println("Energy is: ", energy) end main() diff --git a/contents/split-operator_method/res/Psi_t.png b/contents/split-operator_method/res/Psi_t.png new file mode 100644 index 000000000..e83d14c8a Binary files /dev/null and b/contents/split-operator_method/res/Psi_t.png differ diff --git a/contents/split-operator_method/res/Psi_tdt.png b/contents/split-operator_method/res/Psi_tdt.png new file mode 100644 index 000000000..05a83a893 Binary files /dev/null and b/contents/split-operator_method/res/Psi_tdt.png differ diff --git a/contents/split-operator_method/res/imaginary_time.gif b/contents/split-operator_method/res/imaginary_time.gif new file mode 100644 index 000000000..27cf83fb5 Binary files /dev/null and b/contents/split-operator_method/res/imaginary_time.gif differ diff --git a/contents/split-operator_method/res/imaginary_time.png b/contents/split-operator_method/res/imaginary_time.png new file mode 100644 index 000000000..a01704cfe Binary files /dev/null and b/contents/split-operator_method/res/imaginary_time.png differ diff --git a/contents/split-operator_method/res/imaginary_time_0.png b/contents/split-operator_method/res/imaginary_time_0.png new file mode 100644 index 000000000..813bea737 Binary files /dev/null and b/contents/split-operator_method/res/imaginary_time_0.png differ diff --git a/contents/split-operator_method/res/imaginary_time_75.png b/contents/split-operator_method/res/imaginary_time_75.png new file mode 100644 index 000000000..5550932ac Binary files /dev/null and b/contents/split-operator_method/res/imaginary_time_75.png differ diff --git a/contents/split-operator_method/res/real_time.gif b/contents/split-operator_method/res/real_time.gif new file mode 100644 index 000000000..8a4c68b50 Binary files /dev/null and b/contents/split-operator_method/res/real_time.gif differ diff --git a/contents/split-operator_method/res/real_time.png b/contents/split-operator_method/res/real_time.png new file mode 100644 index 000000000..586594ab7 Binary files /dev/null and b/contents/split-operator_method/res/real_time.png differ diff --git a/contents/split-operator_method/res/real_time_0.png b/contents/split-operator_method/res/real_time_0.png new file mode 100644 index 000000000..87ffa67a5 Binary files /dev/null and b/contents/split-operator_method/res/real_time_0.png differ diff --git a/contents/split-operator_method/res/real_time_75.png b/contents/split-operator_method/res/real_time_75.png new file mode 100644 index 000000000..2824069f0 Binary files /dev/null and b/contents/split-operator_method/res/real_time_75.png differ diff --git a/contents/split-operator_method/res/split_op_method.svg b/contents/split-operator_method/res/split_op_method.svg new file mode 100644 index 000000000..0626eea04 --- /dev/null +++ b/contents/split-operator_method/res/split_op_method.svg @@ -0,0 +1,2119 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + + iFFT + Momentum + Position + + + + FFT + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Position + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/contents/split-operator_method/res/times_Uk.png b/contents/split-operator_method/res/times_Uk.png new file mode 100644 index 000000000..bc901cda3 Binary files /dev/null and b/contents/split-operator_method/res/times_Uk.png differ diff --git a/contents/split-operator_method/res/times_Ur.png b/contents/split-operator_method/res/times_Ur.png new file mode 100644 index 000000000..534712acd Binary files /dev/null and b/contents/split-operator_method/res/times_Ur.png differ diff --git a/contents/split-operator_method/split-operator_method.md b/contents/split-operator_method/split-operator_method.md index 524cd6d18..429218866 100644 --- a/contents/split-operator_method/split-operator_method.md +++ b/contents/split-operator_method/split-operator_method.md @@ -1,87 +1,145 @@ # The Split-Operator Method -The Split-Operator Method (also called the Split-Step Method), was actually the primary method I used to solve the Schrodinger equation during my PhD. -It is one of the simplest and fastest methods for this purpose and is widely used throughout modern quantum research in the area -- in particular when dealing with the Non-linear Schrodinger Equation: +The Split-Operator Method (also called the Split-Step Method), was actually the primary method I used to solve the Schrödinger equation during my PhD. +It is one of the simplest and fastest methods for this purpose and is widely used throughout modern quantum research in the area, in particular when dealing with the Non-linear Schrödinger Equation (NLSE): $$ -i \hbar \frac{\partial \Psi(\mathbf{r},t)}{\partial t} = \left[-\frac{\hbar^2}{2m}\nabla^2 + g|\Psi(\mathbf{r},t)|^2 \right] \Psi(\mathbf{r},t), -$$ - -which follows from the notation provided in the [quantum systems](../quantum_systems/quantum_systems.md) chapter: $$\Psi(\mathbf{r},t)$$ is a quantum wave-function with spatial ($$\mathbf{r}$$) and time ($$t$$) dependence and $$\nabla^2$$ is a laplacian; however, in this case, we also add an interaction tern $$g$$ next to a nonlinear $$|\Psi(\mathbf{r},t)|^2$$ term. -By adding in the $$V(\mathbf{r})$$ term, we get an equation used to study superfluid Bose--Einstein Condensate (BEC) systems: - -$$ -i \hbar \frac{\partial \Psi(\mathbf{r},t)}{\partial t} = \left[-\frac{\hbar^2}{2m}\nabla^2 + V(\mathbf{r}) + g|\Psi(\mathbf{r},t)|^2 \right] \Psi(\mathbf{r},t). +i \hbar \frac{\partial \Psi(\mathbf{r},t)}{\partial t} = \left[-\frac{\hbar^2}{2m}\nabla^2 + V(\mathbf{r}) + g|\Psi(\mathbf{r},t)|^2 \right] \Psi(\mathbf{r},t), $$ +which follows from the notation provided in the [quantum systems](../quantum_systems/quantum_systems.md) chapter: $$\Psi(\mathbf{r},t)$$ is a quantum wave-function with spatial ($$\mathbf{r}$$) and time ($$t$$) dependence, $$\nabla^2$$ is a laplacian, and $$V(\mathbf{r})$$ is a potential of some sort (like $$\omega x^2$$ or something). +In this case, we also add an interaction term $$g$$ next to a nonlinear $$|\Psi(\mathbf{r},t)|^2$$ term. This is the system I studied for most of my PhD (granted, we played a few tricks with parallelization and such, so it was _slightly_ more complicated). -At it's heart, the split-op method is nothing more than a pseudo-spectral differential equation solver... That is to say, it solves the Schrodinger equation with [FFT's](../cooley_tukey/cooley_tukey.md). +At its heart, the split-op method is nothing more than a pseudo-spectral differential equation solver... That is to say, it solves the Schrödinger equation with [FFT's](../cooley_tukey/cooley_tukey.md). In fact, there is a large class of spectral and pseudo-spectral methods used to solve a number of different physical systems, and we'll definitely be covering those in the future. -As mentioned in the [quantum systems](../quantum_systems/quantum_systems.md) section, we can represent a a quantum wavefunction in momentum space, which is parameterized with the wavevector $$k$$. -In the hamiltonian shown above, we can split our system into real-space components, $$\hat{H}_R = \left[V(\mathbf{r}) + g|\Psi(\mathbf{r},t)|^2 \right] \Psi(\mathbf{r},t)$$, and momentum space components, $$\hat{H}_M = \left[-\frac{\hbar^2}{2m}\nabla^2 \right]\Psi(\mathbf{r},t)$$. +As mentioned in the [quantum systems](../quantum_systems/quantum_systems.md) section, we can represent a quantum wavefunction in momentum space, which is parameterized with the wavevector $$k$$. +In the Hamiltonian shown above, we can split our system into position space components, $$\hat{H}_r = \left[V(\mathbf{r}) + g|\Psi(\mathbf{r},t)|^2 \right] \Psi(\mathbf{r},t)$$, and momentum space components, $$\hat{H}_k = \left[-\frac{\hbar^2}{2m}\nabla^2 \right]\Psi(\mathbf{r},t)$$. +I'll be honest, I didn't know what notation to use for $$\hat H_r$$ because $$p$$ is used to describe momentum. +I settled on $$r$$ for _real space_, but that is somewhat notationally ambiguous. +In addition, $$k$$ will indicate momentum space because it is a sum of all wavevectors, typically notated as $$k$$. +Bad notation aside, let's continue. + If we assume a somewhat general solution to our quantum system: $$ -\Psi(\mathbf{r},t + dt) = \left[e^{-\frac{i\hat{H}dt}{\hbar}}\right]\Psi(\mathbf{r},t) = \left[e^{-\frac{i(\hat{H}_R + \hat{H}_M)dt}{\hbar}}\right]\Psi(\mathbf{r},t) +\Psi(\mathbf{r},t + dt) = \left[e^{-\frac{i\hat{H}dt}{\hbar}}\right]\Psi(\mathbf{r},t) = \left[e^{-\frac{i(\hat{H}_r + \hat{H}_k)dt}{\hbar}}\right]\Psi(\mathbf{r},t) $$ -and assume we are simulating our system by a series of small tiemsteps ($$dt$$), we can perform similar splitting by using the Baker-Campbell-Housdorff formula: +and assume we are simulating our system by a series of small timesteps ($$dt$$), we can perform similar splitting by using the Baker-Campbell-Housdorff formula: $$ -\Psi(\mathbf{r},t+dt) = \left[e^{-\frac{i\hat{H}_Rdt}{\hbar}}e^{-\frac{i\hat{H}_Mdt}{\hbar}}e^{-\frac{[i\hat{H}_R, i\hat{H}_M]dt^2}{2}}\right]\Psi(\mathbf{r},t) +\Psi(\mathbf{r},t+dt) = \left[e^{-\frac{i\hat{H}_rdt}{\hbar}}e^{-\frac{i\hat{H}_kdt}{\hbar}}e^{-\frac{[i\hat{H}_r, i\hat{H}_k]dt^2}{2}}\right]\Psi(\mathbf{r},t) $$ -This accrues a small amount of error ($$dt^2$$) related to the commutation of the real and momentum-space components of the Hamiltonian. That's not okay. -In order to change the $$dt^2$$ error to $$dt^3$$, we can split the system by performing a half-step in real space before doing a full-step in momentum space, through a process called _Strang Splitting_ like so: +This accrues a small amount of error ($$dt^2$$) related to the commutation of the real and momentum-space components of the Hamiltonian. +This is a relatively large error and that's not okay. +In order to change the $$dt^2$$ error to $$dt^3$$, we can split the system by performing a half-step in position space before doing a full-step in momentum space, through a process called _Strang Splitting_ like so: $$ -\Psi(\mathbf{r},t+dt) = \left[e^{-\frac{i\hat{H}_Rdt}{2\hbar}}e^{-\frac{i\hat{H}_Mdt}{\hbar}}e^{-\frac{i\hat{H}_Rdt}{2\hbar}} \right]\Psi(\mathbf{r},t) + \mathcal{O}(dt^3) +\Psi(\mathbf{r},t+dt) = \left[e^{-\frac{i\hat{H}_rdt}{2\hbar}}e^{-\frac{i\hat{H}_kdt}{\hbar}}e^{-\frac{i\hat{H}_rdt}{2\hbar}} \right]\Psi(\mathbf{r},t) + \mathcal{O}(dt^3) $$ -We can then address each part of this solution in chunks, first in real space, then in momentum space, then in real space again by using [Fourier Transforms](../cooley_tukey/cooley_tukey.md). +We can then address each part of this solution in chunks, first in position space, then in momentum space, then in position space again by using [Fourier Transforms](../cooley_tukey/cooley_tukey.md). Which looks something like this: $$ -\Psi(\mathcal{r}, t+dt) = \left[\hat{U}_R(\frac{dt}{2})\mathcal{F}\left[\hat{U}_M(dt) \mathcal{F} \left[\hat{U}_R(\frac{dt}{2}) \Psi(\mathbf{r},t) \right] \right] \right] + \mathcal{O}(dt^3) +\Psi(\mathcal{r}, t+dt) = \left[\hat{U}_r\left(\frac{dt}{2}\right)\mathcal{F}^{-1}\left[\hat{U}_k(dt) \mathcal{F} \left[\hat{U}_r\left(\frac{dt}{2}\right) \Psi(\mathbf{r},t) \right] \right] \right] + \mathcal{O}(dt^3) $$ -where $$\hat{U}_R = e^{-\frac{i\hat{H}_Rdt}{\hbar}}$$, $$\hat{U}_M = e^{-\frac{i\hat{H}_Mdt}{\hbar}}$$, and $$\mathcal{F}$$ and $$\mathcal{F}^{-1}$$ indicate forward and inverse Fourier Transforms. +where $$\hat{U}_r = e^{-\frac{i\hat{H}_rdt}{\hbar}}$$, $$\hat{U}_k = e^{-\frac{i\hat{H}_kdt}{\hbar}}$$, and $$\mathcal{F}$$ and $$\mathcal{F}^{-1}$$ indicate forward and inverse Fourier Transforms. +Here's a flowchart of what we are looking for every timestep: + +

+ +

-As a small concession here, using this method enforces periodic boundary conditions, where the wavefunction will simply slide from one side of your simulation box to the other, but that's fine for most cases. + +For the most part, that's it: +1. Multiply the wavefunction in real space with the real-space operator. +2. Flip to momentum space with a Fourier transform. +3. Multiply the momentum-space wavefuntion by the momentum-space operator. +4. Flip to position space with an inverse Fourier transform. +5. Repeat 1-4 until satisfied. + +If we guess that our initial wavefunction is gaussian-like and is slightly offset from the center or the trap, this should allow us to see our wavefunction "sloshing" back and forth in our trap, like so: + +

+ +

+ +As a small concession, using this method enforces periodic boundary conditions, where the wavefunction will simply slide from one side of your simulation box to the other, but that's fine for most cases. In fact, for many cases (such as large-scale turbulence models) it's ideal. +That said, there is more to the story. +As we mentioned in the [quantum systems](../quantum_systems/quantum_systems.md) section, many simulations of quantum systems desire to find the ground state of our system. +The split-operator method can be used for that too! +If we run this simulation in _imaginary time_, by simply setting $$\tau = it$$ and stepping through $$\tau$$ instead of $$t$$, we will no longer see an "real-world" example of how the atoms should behave, but will instead see an exponential decay of higher-energy states. +If we run the simulation for long enough with a small enough timestep, all higher energy states will vanish. +This means that we can find the ground state of our system by running the simulation in imaginary time, which is an incredibly useful feature! +If we run the same simulation as above in imaginary time, we should see our wavefunction smoothly move to the center of our trap (the lowest energy position), like so: + +

+ +

+ + +## The Algorithm + Luckily, the code in this case is pretty straightforward. -Frist, we need to set all the initial parameters, including the initial grids in real and momentum space: +As a note before starting, we will be using normalized units in this simulation where $$\hbar = c = 1$$. +These units are often called _natural_ units. +Many of you (*cough* experimentalists *cough*) will probably think that these units are completely unphysical, and they are; however, they allow us to output fractions and whole numbers. +For example, if we are trying to find the energy of the ground state of atoms in a simple harmonic oscillator, we know it should be $$\frac{1}{2}\hbar \omega$$, where $$\omega$$ is the coefficient in front of the $$x^2$$ term known as the _frequency_ of the trap. +If we were to calculate the energy in real units, our simulation would output $$5.272859 \times 10^{-35}$$, which is hard to interpret. +By instead using natural units, we get precisely $$\frac{1}{2}$$ and we know that those are in units of $$\hbar\omega$$. +There is no doubt that it makes the simulation easier to understand (albeit a little misleading in the end). + +Regardless, we first need to set all the initial parameters, including the initial grids in real and momentum space: {% method %} {% sample lang="jl" %} -[import:4-32, lang:"julia"](code/julia/split_op.jl) +[import:9-32, lang:"julia"](code/julia/split_op.jl) {% endmethod %} As a note, when we generate our grid in momentum space `k`, we need to split the grid into two lines, one that is going from `0` to `-kmax` and is then discontinuous and goes from `kmax` to `0`. This is simply because the FFT will naturally assume that the `0` in our grid is at the left side of the simulation, so we shift k-space to match this expectation. +Also, for this code we will be using the notation from above: `opr.R` will be the real space operators and `opr.M` will be the momentum space operators. +There is another boolean value here called `im_time`, which is for imaginary time evolution. + Afterwards, we turn them into operators: {% method %} {% sample lang="jl" %} -[import:34-42, lang:"julia"](code/julia/split_op.jl) +[import:34-60, lang:"julia"](code/julia/split_op.jl) {% endmethod %} Here, we use a standard harmonic potential for the atoms to sit in and a gaussian distribution for an initial guess for the probability distribution. -As a note, if we run this simulation in _imaginary time_, by simply setting $$\tau = it$$ and stepping through $$\tau$$, we will no longer see an "real-world" example of how the atoms should behave, but will instead see an exponential decay of higher-energy states. -This means that we can find the ground state of our system by running the simulation in imaginary time, which is an incredibly useful feature! +If we give either the trap or the atoms a slight offset (so the gaussian distribution of atoms does not *quite* rest at the bottom of the $$x^2$$ potential, we can see the atoms moving back and forth in the potential as we move the simulation forward in time. +This means that we can easily see the dynamics of our quantum system! +If we run the simulation in imaginary time, we will see the gaussian distribution of atoms move towards the center of the potential, which is the location with the lowest energy. +Both of these have been shown in the figures above. -And finally go step-by-step through the simulation: +The final step is to do the iteration, itself. {% method %} {% sample lang="jl" %} -[import:44-70, lang:"julia"](code/julia/split_op.jl) +[import:63-109, lang:"julia"](code/julia/split_op.jl) {% endmethod %} And that's it. + +There is something a bit odd about the simulation in imaginary time, though. +Basically, in imaginary time, we see an exponential decay of all the higher energy states, which means we are technically losing a large amount of our wavefunction density every timestep! +To solve this issue, we _renormalize_ by enforcing that $$\int_{-\infty}^{+\infty}\Psi^\ast\Psi dx = 1$$. +As you can see from the code, this involves summing the density, multiplying that sum by `dx`, and then dividing each element in the wavefunction by the `sqrt()` of that value. + The Split-Operator method is one of the most commonly used quantum simulation algorithms because of how straightforward it is to code and how quickly you can start really digging into the physics of the simulation results! ## Example Code +This example code is a simulation of a gaussian distribution of atoms slightly offset in a harmonic trap in imaginary time. +So long as the code is written appropriately, this means that the atoms should move towards the center of the trap and the energy should decay to $$\frac{1}{2}\hbar\omega$$, which will be simply $$\frac{1}{2}$$ in this simulation. +Checking to make sure your code can output the correct energy for a harmonic trap is a good test to make sure it is all working under-the-hood before simulating systems with more complicated Hamiltonians. + {% method %} {% sample lang="jl" %} [import, lang:"julia"](code/julia/split_op.jl)