Quantum Systems

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 densities. These densities are easiest to understand in terms of wavefunctions, which are complex functions characterizing a quantum system's behavior.

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?" Well, here's where the notation gets really sticky.

There is a clear difference between quantum computers and quantum simulators. A quantum computer is the quantum analog to a classical computer, replacing bits with qubits by using quantum information theory. 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 understand quantum mechanics. These will often come in the form of experimental quantum systems that express quantum behavior 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 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 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:

$$\frac{\partial\phi(\mathbf{r},t)}{\partial t} = D \nabla^2 \phi(\mathbf{r},t)$$

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 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

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:

$$P(\mathbf{r}, t) = |\Psi(\mathbf{r},t)|^2 = \Psi(\mathbf{r},t)^{*}\Psi(\mathbf{r},t)$$

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 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):

$$\int_{-\infty}^{+\infty}|\Psi(\mathbf{r},t)|^2 d\mathbf{r} = 1$$

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} = \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 position and momentum space by using Fourier Transforms, 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.

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)$.

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. The bra represents the Hermitian conjugate of the ket and looks like this:

$$\langle B \rvert$$

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 exclusively 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:

# We are calculating the energy to check <Psi|H|Psi>
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

calculateEnergy :: Parameters -> Operators -> Double
calculateEnergy param ops = (* dx param) . sum . map realPart $ elems totalE
  where
    totalE = potentialE .+ kineticE
    potentialE = wfcConj .* v ops .* wfc ops
    kineticOp = liftArray ((/ 2) . (^ 2)) (ks param)
    kineticE = wfcConj .* idft (kineticOp .* dft (wfc ops))
    wfcConj = liftArray conjugate $ wfc ops

double calculate_energy(struct params par, struct operators opr) {
    double complex wfc_r[opr.size];
    double complex wfc_k[opr.size];
    double complex wfc_c[opr.size];
    memcpy(wfc_r, opr.wfc, sizeof(wfc_r));

    memcpy(wfc_k, opr.wfc, sizeof(wfc_k));
    fft(wfc_k, opr.size, false);

    for (size_t i = 0; i < opr.size; ++i) {
        wfc_c[i] = conj(wfc_r[i]);
    }

    double complex energy_k[opr.size];
    double complex energy_r[opr.size];

    for (size_t i = 0; i < opr.size; ++i) {
        energy_k[i] = wfc_k[i] * cpow(par.k[i] + 0.0*I, 2);
    }

    fft(energy_k, opr.size, true);

    for (size_t i = 0; i < opr.size; ++i) {
        energy_k[i] *= 0.5 * wfc_c[i];
        energy_r[i] = wfc_c[i] * opr.v[i] * wfc_r[i];
    }

    double energy_final = 0;

    for (size_t i = 0; i < opr.size; ++i) {
        energy_final += creal(energy_k[i] + energy_r[i]);
    }

    return energy_final * par.dx;
}

double calculate_energy(Params &par, Operators &opr) {
    vector_complex wfc_r(opr.wfc);
    vector_complex wfc_k(opr.wfc);
    vector_complex wfc_c(opr.size);
    fft(wfc_k, false);

    for (size_t i = 0; i < opr.size; ++i) {
        wfc_c[i] = conj(wfc_r[i]);
    }

    vector_complex energy_k(opr.size);
    vector_complex energy_r(opr.size);

    for (size_t i = 0; i < opr.size; ++i) {
        energy_k[i] = wfc_k[i] * pow(complex(par.k[i], 0.0), 2);
    }

    fft(energy_k, true);

    for (size_t i = 0; i < opr.size; ++i) {
        energy_k[i] *= 0.5 * wfc_c[i];
        energy_r[i] = wfc_c[i] * opr.v[i] * wfc_r[i];
    }

    double energy_final = 0;

    for (size_t i = 0; i < opr.size; ++i) {
        energy_final += real(energy_k[i] + energy_r[i]);
    }

    return energy_final * par.dx;
}

def calculate_energy(par: Param, opr: Operators) -> float:
    """Calculate the energy <Psi|H|Psi>."""
    # Creating real, momentum, and conjugate wavefunctions.
    wfc_r = opr.wfc
    wfc_k = np.fft.fft(wfc_r)
    wfc_c = np.conj(wfc_r)

    # Finding the momentum and real-space energy terms
    energy_k = 0.5 * wfc_c * np.fft.ifft((par.k ** 2) * wfc_k)
    energy_r = wfc_c * opr.V * wfc_r

    # Integrating over all space
    energy_final = sum(energy_k + energy_r).real

    return energy_final * par.dx

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: algorithms and methods researchers use to ascertain new knowledge 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 one of those areas. In fact, quantum information theory 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!

License

Code Examples

The code examples are licensed under the MIT license (found in LICENSE.md).

Text

The text of this chapter was written by James Schloss and is licensed under the Creative Commons Attribution-ShareAlike 4.0 International License.

Images/Graphics

• The animation "FTgaussian" was created by James Schloss and is licensed under the Creative Commons Attribution-ShareAlike 4.0 International License.

Pull Requests

After initial licensing (#560), the following pull requests have modified the text or graphics of this chapter:

• none