How to perform Physical Simulation - RPIQuantumComputing/QuantumCircuits GitHub Wiki

Quantum Computing + Physical Simulation

"Nature isn't classical, dammit, and, if you want to make a simulation of nature, you'd better make it quantum mechanical, and, by golly, it's a wonderful problem, because it doesn't look so easy." - Richard Feynman

Quantum Computing has an obvious application-simulating quantum mechanics, but the implementation of such simulations is not as obvious. And, in general, the problem of general quantum mechanics simulations, as alluted to above by Feynman, isn't easy.

First, I have to introduce an illusive concept in physics: The Hamiltonian! You know how your high school physics teacher told you that energy is such a useful concept that it is used all throughout physics. Well, I am here to tell you-they were not lying. In fact, you can throw out almost everything you know about forces beyond knowing how to integrate individual ones across a given path and get pretty far in physics.

Let us consider a hard problem to solve in classical mechanics that, when examined in a cursory fashion, looks too easy:

The image is from Wikipedia. Let us compute the energy of this system. What force is acting on the system: Gravity! What is the energy of gravity? $U_{g} = mgh$, where $m$ is the object's mass, g is the acceleration of gravity (9.81 m/s^2), and $h$ is one's height from the ground.

Now, we know the length of each line segment, the mass of the ball, and a natural parameter for describing the system would be the angle relative to the pivot for the first mass and then the angle from the x-axis defined at the current position of the other ball for the second mass.

$x_{1} = l_{1}\sin{(\theta_{1})}$

$y_{1} = -l_{1}\cos{(\theta_{1})}$

$x_{2} = l_{1}\sin{(\theta_{1})} + l_{2}\sin{(\theta_{2})}$

$y_{2} = -l_{1}\cos{(\theta_{1})} - l_{2}\cos{(\theta_{2})}$

Above, I simply expressed the positions in terms of these "generalized" variables, called generalized as I can define any weird coordinate system for which they would range. Now, there is also something called kinetic energy, the energy associated with the current movement of the masses. Its energy equation is $U_{KE} = \frac{mv^{2}}{2}$. So, we need to also find the $v$s, which definitionally are the derivatives of our positional coordinates.

$x_{1}' = l_{1}\cos{(\theta_{1})}$

$y_{1}' = l_{1}\theta_{1}'\sin{(\theta_{1})}$

$x_{2}' = l_{1}\cos{(\theta_{1})} + l_{2}\cos{(\theta_{2})}$

$y_{2}' = l_{1}\theta_{1}'\sin{(\theta_{1})} + l_{2}\theta_{2}'\sin{(\theta_{2})}$

And, well, from an energy perspective, we are done! Now, let us combine everything, as the description of the system is defined by the energy interactions. T will denote kinetic energy while U will denote potential energy.

$T = \frac{m_{1}v_{1}^{2}}{2} + \frac{m_{2}v_{2}^{2}}{2}$

$T = \frac{m_{1}}{2}(x_{1}'^{2} + y_{1}^{2}) + \frac{m_{2}}{2}(x_{2}'^{2} + y_{2}^{2})$

$T = \frac{m_{1}l_{1}^{2}\theta_{1}'^{2}}{2} + \frac{m_{2}}{2}(l_{1}^{2}\theta_{1}'^{2} + l_{2}^{2}\theta_{2}'^{2} + 2l_{1}l_{2}\theta_{1}'\theta_{2}'\cos{(\theta_{1} - \theta_{2})})$

The last step used $\cos{(\theta_{1})}\cos{(\theta_{2})} + \sin{(\theta_{1})}\sin{(\theta_{2})} = \cos{(\theta_{1} - \theta_{2})}$.

$V = m_{1}gy_{1} + m_{2}gy_{2}$

$V = -m_{1}gl_{1}\cos{(\theta_{1})} - m_{2}g(l_{1}\cos{(\theta_{1})} + l_{2}\cos{(\theta_{2})}$

$V = -(m_{1} + m_{2})gl_{1}\cos{(\theta_{1})} - m_{2}gl_{2}\cos{(\theta_{2})}$

Then, we can define something called the Lagrangian: $L = T - U$.

To convert to the Hamiltonian, ones does the following:

$H = \sum_{i=1}^{2} \theta_{i}'p_{theta_{i}} - L$

$p_{theta_{i}}= \frac{\partial L}{\partial \theta_{i}'}$

But, I will spare you the tedious calculus exercise and just show you the result, like any good professor. Now, this is fun and all, but how does it help us simulate the system?

Well, nature is lazy and, so, it tends to ensure $L = T - U$ is as minimal as possible, so it becomes an exercise in finding what $\theta_{1}(t), \theta_{2}(t)$ implies the smallest change to $L$. As much as I hate to do so, I shalt not derive this result from scratch and thus leave it to someone much better at explaining advanced calculus than I am-some guy on YouTube: Video.

$\frac{d}{dt}(\frac{\partial L}{\partial \theta_{i}'}) = \frac{\partial L}{\partial \theta_{i}}$

I will say, however, that one can numerically solve the above equation, so what the formulation really does is give us something we can tell a computer to do. Now, the Hamiltonian is just $L = T + U$ for most problems, so it isn't necessarily a harder thing to derive in most cases. But, it gives some more physical understanding. One, it allows one to define momentum in the usual way in our coordinate system which implies a system of first order that is typically easier than the second order system given by $L$. It also allows us to see which momentas are time-invariant, classify the symmetries of our problem such as being invariant under a coordinate translation, and get a gague for how one may construct the force equations using it, if one so desired.

As complicated as our Hamiltonian looks, it should amaze your that simply given the inital velocities and positions for the two balls and plugging in our complicate expression to a computer, we can compute the following complex path image: [Wikipedia](https://en.wikipedia.org/wiki/Double_pendulum :

I digress. I simply wanted to introduce the Hamiltonian in a non-quantum context, as it really is a non-quantum thing. It is just that physists are sorta lazy and don't like to reinvent the wheel when possible. Hence, the formalism follows with us into the weird, weird world of quantum mechanics.

Suppose we have a lattice of atoms defining a material. If the grid repeats, one can simplify the problem by supposing our sample is so large that is is basically infinite. Now, that may seem weird-why consider a more difficult problem? Well, in some sense, it is easier as, now, one can say whatever interaction occurs for one small segment at the top must also follow the one at the bottom, due to being completely symmetrical and infinite in grid size. This is called a periodic boundry condition.

I am about to destory your dreams of simulating the nucleus-it is kind of hard. But, for most problems, it is just kind of there and really only adds mass and enacts the necessary force of pulling electrons close to it into orbitals defining the atoms. So, treating its center as one object with those properties as an approximation for it being decomposable into neutrons and protons is both a fair and necessary assumption as otherwise the problem becomes just too difficult too quickly.

Most of the time, one will use something like Density Functional Theory or Density Renormalization Group Theory or some other fancy sounding simulation method. But, those are just more approximations-we want to actually simulate the material with the fullest generality.

Now, what is a material if not just a repeating grid of buckets, being the protons, with electrons being the objects within those buckets? Laugh as you may, but this is a logically equivelent way of conceptualizing a material. If we define the process of electron exchange (or tunneling) as the removal from one bin to the addition of another using something called creation and anihillation operations one associated with each electron, we can define interactions as both a removal and addition plus an addition and removal from the two interacting sites (weird, I know, but essentially it allows one to interact from within one's own bin to another's bin), the chemical potential (intrinsic energy of bonds) as the energy required to add and then remove each site, and the spin energy (electrons are little bar magnets, so they have something called spin up and spin down) as the product of the spin up and spin down chemical potential terms which is the same as the Coulomb potential (the electromagnetic interaction between the spins) up to some units. Obviously, one needs to imply other constraints such as electron number conservation, but this fully describes the electrical aspect of the material. But, from a chemistry perspective, all that matters is the electrons baby!

$H = - t \sum_{i=0}^{N-1} \sum_{\sigma=\uparrow, \downarrow} (a^\dagger_{i, \sigma} a_{i+1, \sigma} + a^\dagger_{i+1, \sigma} a_{i, \sigma}) - \mu \sum_{i=0}^{N-1} \sum_{\sigma=\uparrow, \downarrow} a^\dagger_{i, \sigma} a_{i, \sigma} + U \sum_{i=0}^{N-1} a^\dagger_{i, \uparrow} a_{i, \uparrow} a^\dagger_{i, \downarrow} a_{i, \downarrow}$

Now, this is just one particular ansatz and there may be more complication when considering a model. For example, maybe an external magnetic field is present, all the atoms are different, or some exotic interaction happens like superconductivity that requires a more advanced model, but, luckily for us, developing the Hamiltonian is usually the easy/route part. The challenge is actually in solving it!

Why is this good? It means we can solve one problem which, in effect, solves them all! How powerful of an idea is that?

Now, I said before that one may have to consider some exotic interaction. What if we expose our system to a weird enviornment which changes the electromagnetic potential around the material? Now, it might be hard to compute the energy of this new contribution directly, but what if we examined how it effects the creation and annihilation of electrons at each site through experiement indirectly by examining the distributional shift in positions with its incremental addition?

As we should still have the same number of electrons, all this field can do shift the rate of creation and annihilation at sites meaning the transformation is unitary, norm preserving, as each element is the number of electrons at position i, implying that one is not changing electron number. Thus, we are dealing with a linear transformation of our system and the experiment, in effect, simply quantifies a transformation matrix.

$c_{2i}^{\dagger} = a_{i, \text{up}}^{\dagger}$

$c_{2i + 1}^{\dagger} = a_{i, \text{down}}^{\dagger}$

$\bar{c_{i}^{\dagger}} = \sum_{k = 0}^{2N - 1} u_{ik}c_{k}^{\dagger}$

$\bar{c_{i}} = \sum_{k = 0}^{2N - 1} u_{ik}^{*} c_{k}$

Now, the model given is...well...not used that commonly in practice. The reason being is it forces a geometry, assumes a very regular and certain exchange mechanism, assume periodicity, and kind of completely neglects the differences our wells, atoms, might have. But, with this change, it becomes similarly difficult.

The specific Hamiltonian used here was originally found by the article author by QuanaSys's Quantum Challenge.

Now, let us introduce a way of solving this Hamtilonian:

In the above diagram, made by the good folks at PennyLane, one can see that this method involves a quantum computer and a classical computer. Now, before I can really get into this method, I have to motivate the approach.

When one wants to compute $f(x) = \sin{x}$, one could just find a way of computing $\sin{x}$ directly or one could try to find a similar function to it that is easier to work exactly with. We will be doing the latter, as converting arbitrary physical systems to other arbitrary physical systems by direct analogy is pretty hard! Before you hurl bricks at me arguing that DFT was mocked by me as approximations even though it tends to also do a similar process of projecting the system empirically onto another system, I have to defend the mocking by pointing out a bigger issue than this mapping: DFT restricts the Hamtiltonian's complexity and usually doesn't simulate all possible configurations. Instead, it simulates a number of them, adds new terms to simulate other types of more complicated exchanges, and relies on its restrictions. Not to say that doesn't work well, but, regardless, it limits the size of the system and even moderately sized simulations are known to take on the order of hours to days to run.

Hopefully, I have avoided the stoning, but that will remain to be seen. How do we do this conversion process? Simple, we guess a function: $g(\theta, f(x))$ where $x$ is our system's inital configuration (i.e. a random guess obeying the electron number restriction), $f$ is just some function to do something to imbed the configuration into physically valid situation relative to the system $g$, $g$ is a representation for our solution, and $\theta$ is our learned parameters. We could try making $g(\theta, f(x)) \approx H$, but that is kind of indirect as, what we really want, is the solution to the set of differential equations implied by $H$. $g(\theta, f(x))$ in effect will be the solution function itself!

That seems a little odd as...well...where are our dynamic parameters? Why describe our system as a function? This is where my "we do the same thing in quantum mechanics" is going to bite me a bit. We actually consider things called wavefunctions, functions who can be acted upon with specific operators to get the desired value. Thus, one applies a function to our function which indirectly encodes the situation being enacted on our wavefunction to get the answer. This is because a definite value for "position" and "velocity" begins to lose its meaning in the quantum context and instead it is easier to say we get out a probability distribution for these as our output. So, if we want a function as output, we better encode our dynamic parameters as arguments to create our operators which convert our wavefunction into the probability distribution we desire. The reason we can do this without caring specifically about how complicate our applied operation will be is simply due to the fact that physics is just complexity on simple ideas most of the time and, so far, it has worked without incidient.

I encourage one just to accept that observations can be decomposed as linear operators, as proving so is beyond my ability and I am not quite sure how believable and complete an explaination one can even give to such a valid, thought-provoking question.

But, if one does that, our approach simply needs a way for us to find $\theta$. Our inital state of the system is going to decomposed into the inital state of its units (our inital wells), we have some function we parameterize which takes it to give us our wavefunction, and now we wish to "judge" how well it behaves. But, if you recall, I said $L = T - V$ and $H = T + V$, but $H$ is more than just $H = 2V + L$! It is total energy, so wouldn't it just be swell if we were to minimize our energy by applying $H$ to the produced product state of our inital wells? Now, now, I can really hear the chior singing a horrific tone at that. I should explain that knowing the ground state allows us to find the excited states, by simply applying the operations used to add energy to the system.

To do this, we need to apply an arbitrary basis to measure our qubits-easy to do with some premade tools (reduces to converting our operators into Pauli operators and possibly our inputs units into multiple qubits if more than 2 basis states are required for each). This will measure our input relative to $H$, then we use a classical optimizer to generate new guesses for $\theta$.

Okay, sure, this variational approach, called Quantum Variational Eigensolving (as finding the ground state corresponds to solving our smallest eigenvalue if our Hamiltonian is a matrix, mathematically), solves our problem. Why is it better? As alluded to before, it is operating in a space of possibilities equal to that of the true system not an arbitrarily truncated one, so all interactions should be mathematically representable and thus enactable. But, there is the added benefit of the fact that measuring an arbitrary $H$ takes time proportional to the terms in $H$, polynomial for physical systems, state preparation is usually a fairly trivial operation taking a polynomial amount of time and usually a square or less amount of time for practical instances, and the ansatz can be whatever we wish it to be. But, who is to say our ansatz can be enacted in poly time? Well, its job is to approximate our system, so its growth should be similar to that of our system and our system grows in poly time. But, we would never choose an ansatz that couldn't be done in poly time to simulate. Luckily for us, any sufficiently complicated and symmetry respecting ansatz can sufficiently approximate any physically possible quantum system in polynomial time.

The reason for the final line has to do with something called the area growth law of Hilbert spaces in quantum mechanics, but I will leave that to other works to tell you more about. Thus, the hard work is over-you now should have a good conceptual understanding for how to perform physical simulation on a quantum computer from the ground up! What will you do with this power?

In the above figure, one can see an example of running VQE on a real life quantum computer! There is a few things to note: 1) the ideal solution here is actually 3x lower than either of these, 2) the difference between theoretical and experimental tends to be within a constant factor of each other (due to VQE requiring some of its parameterization to serve as noise mitigation), and 3) smoothing and checking for convergence is necessary for VQE as noise can easily overwealm-especially gradient-based approaches-the signal making the results meaningless. The code for VQE is starkly similar to the VQE code demo quantum circuits has and is an adviseable starting point for those exploring VQE. On thing to note, the poor performance of our ansatz here is the reason why our results are less than idea and typically one either needs to add more layers of an ansatz, better allign the ansatz with the hardware, or reconsider the ansatz used in one's approach. This is mainly an trial-and-error type of thing, but, as with other optimization tasks, that is part of the fun!