 Ja olen decided not to give a tutorial because I somehow felt that, well, I didn't know what to do. So, but I will show some snippets of pseudo code and if we have time at least after today's lectures, I will show some programs that I put online and maybe show you a little bit how to run them if you want to try yourselves. Anyway, so today I will talk about this stochastic series expansion method and again I will focus on quantum spins. This method can be used for other systems as well, but in these lectures I will basically just talk about quantum spins. The second lecture will be on a different kind of method, so this is a finite temperature method, although you can go to basically the ground state as well if you go to low enough temperature. But there's another way to more directly go to the ground state by projection and that can actually be done for some quantum spin systems very conveniently if you go to the valence bond basis of the spins. So I will talk about that. And then I don't want to just talk about the basic techniques to generate numerical data on our finite lattices. I also want to talk about how to analyze those results to draw conclusions about the thermodynamic limit, which is normally what we want to do. So I want to talk about several things, but in particular a way which is quite convenient for analyzing critical points. And then the last lecture I want to talk about a slightly different topic, but it's still Monte Carlo and still spins, but it has to do with out of equilibrium physics. So in particular things like quenches, ramps and annealing and quantum annealing. And there's the connections to quantum computing there as well, which I want to talk about. Okay, but today, Stochastic series expansion, I can mention that I have this pretty long review article which discusses quantum spins, methods and also some physics. 207 pages, something like that, and it's published in some conference series, but you can easily find it on archive. Okay, so today's, or let's say this first lecture today. I want to talk about the models a little bit first before jumping into the methods just so that we have an idea of what systems we are actually going to study. So just briefly about the Heisenberg model, which is the basic model that everything starts from here. And then although I will talk about the Stochastic series expansion method, I first want to give a brief introduction to path integrals on the lattice. You have already seen path integrals during this course, but maybe not on the lattice. And there's a quite neat relationship between the Stochastic series and the path integrals. Convey that. And then I get into the actual Stochastic series expansion, which is again the alternative to doing a path integral. And I want to show in a little bit of detail how to implement that for the spin one half Heisenberg model in 2D in this case, but it can be done easily in other dimensions as well. And if you have time, I will look at a program that I put online and I will tell you where you can find that program as well. Okay, so quantum anti ferromagnet. So many materials are well described by the Heisenberg model, which looks like this. So this is a spin rotationally invariant interaction. And of course the spin on each side can be anything, but today I will focus on spin one half. And in the simplest case, this lattice where these spins sit on is something like a square lattice and you have only nearest neighbor interactions. But of course you can have other lattices as well. And the lattices in general you can classify into two kinds or I should say the interactions combined with the lattice. So in a bipartite system, these interacting neighbors always sit on different sub lattices. So it's somehow possible to decompose the lattice into two groups of sites. Here I denote them by the white and red dots. So they are forming this checkerboard pattern. And you see that if I have nearest neighbor interactions here, this I and J here is always on different sub lattices. Okay, so that's called bipartite. And that's, if we have an anti ferromagnet, that's compatible with nail order. So this anti ferromagnetic pattern. So here I put these arrows to indicate the spins. And you can see that of course on the white sites I have one orientation of the spins. And on the red sites the orientation is the opposite. So always if you have a bipartite situation, you can have a kind of anti ferromagnet on that system. Okay, well I mean you can have anti ferromagnets in other cases too. But it's easy to see that if you have a classical system, for example, it goes into that state. But you can also have other states. So even if you have some interaction like that, and okay maybe we later add some other terms. Even if the system can have nail order, in the quantum case it doesn't necessarily have. So I want to say a few more words about that later. But let's look at the bipartite case first. So in that case we cannot satisfy the conditions mentioned here. So the simplest case of that is just the triangle. Or you can think of an extended triangular lattice but I just show one triangle here. So then we have what's also called frustration. So if you think about wanting to form an anti ferromagnetic state it's just not possible. So here if I have just icing spins and I can orient them as these colored dots here. Say I have down and up. You see that there's one bond which is not anti ferromagnetic oriented. And it's not possible to have all bonds satisfied, happy with the anti ferromagnetic interaction. So then you can get some other kind of order. So in fact in the triangular lattice you get something like 120 degree anti ferromagnet. Where they are not pointing opposite to each other. They form this angle. And in general if you have some more complicated lattices you can have all kinds of states including spin liquids which have no order at all. Okay. So today I will focus on bipartite lattices and that's because those are the ones we can actually do quantum Monte Carlo simulations. As I will mention a little bit further along for these non bipartite cases we have the so called sine problem which makes simulations very hard. For these systems and not only for the Heisenberg model one can put other kinds of interactions here too which somehow have a bipartite property. We can do very nice simulations. Okay so that means that we are normally studying anti ferromagnetism and we may be studying the destruction of anti ferromagnetism when we get to some other state. The order parameter we are dealing with is what we call the sub lattice magnetization or sometimes we call it the stagger magnetization. So it's just the operator is just summing up all the spins but putting this phase factor in which corresponds to rotating half of the spins. So if you have a classical system and you go to zero temperature. For example if you just look at this case here then this would be the ground state configuration so then this expectation value. And if you take the absolute value of it because it could be you know oriented any in any direction in spin space the spins of course have three components. So if I take the just the magnitude of the order parameter it will always be as classically at zero temperature. But if you look at a finite system or if you really look at the finite system and then really think of the partition function as including all possible ground state configurations. Of course there's an infinite number of them so if you break the symmetry we get into one of them but in principle in a finite lattice we should include all of them. So then the magnetization expectation value itself is formally zero but we can always avoid that subtlety by looking at the squared order parameter which doesn't care about how the spins are overall oriented. What what angle they form. Okay okay so again I'm still talking about classical systems here if we raise the temperature then of course this order parameter is is reduced okay. And at some point it will vanish in a phase transition right okay but if you look at quantum systems and we stay at zero temperatures we look at the ground state. Then you can see that this kind of state is not even an eigen state of the Hamiltonian if you I didn't repeat the Hamiltonian here but you can just try and see this is not an eigen state of the Hamiltonian. So that means that if there is anti ferromagnetic order in the ground state there must be some fluctuations on top of it so that means that this order parameter is reduced from its maximal classical value even at zero temperature. And in fact the order can completely go away at zero temperature as well. If there are enough quantum fluctuations and we will talk about in particular in the third lecture I will you know talk about such things more the first two lectures will be mainly focused on on methods but I will maybe mention some. Some of these issues further along as well okay but there's an important theorem which applies to. To breaking of a continuous symmetry so in this Heisenberg model the Hamiltonian has the full rotational invariance right so if we. Break the symmetry we break a continuous symmetry and the Mermin Wagner theorem tells us that. Let me not talk about where this comes from but I'm sure many of you know it already. That if you take the Heisenberg quantum model on a one D chain there can be no nail order that order is destroyed by quantum fluctuations in one dimension. And if we go to two dimensions we can have order only at zero temperature the order is destroyed by thermal fluctuations so that's true you know both classically and quantum mechanically. Okay and in 3D we can have order at zero temperature as well as a finite temperature so. So this is an important thing to keep in mind when we look at systems in different dimensions. And we have this full rotational invariance of the Hamiltonian no spin and isotropy. Okay so now let's come to quantum Monte Carlo. So what we are doing in Monte Carlo is we are taking this quantum mechanical you know thermal expectation value. Which is defined with you know the operator H and some operator A that we want to calculate the expectation value of we have to somehow rewrite this operator expression. In a form which is what we use in important sampling in Monte Carlo simulations as Werner talked about yesterday. So we have to have some space of configurations. So C denotes some configuration and now in the quantum case it will the configurations can look more complicated than you know the lattice you start from and the spins that you put there. It can be something else. There's just something that gives us some configuration space with some weight function and then this observable gets some value in each of those configurations. So then we can sample this with the weight W as you would do in a classical simulation. Okay right now in the quantum case there is the possibility that actually we may not know how to do this in a way which. Gives all positive weights and in Monte Carlo simulations of course you use the weights as relative probabilities. So if they are not positive definite you have a problem and that's called the sign problem and it's it's really bad. It's it's the you know it's it's plaguing you know the whole field. On the other hand you can say that well of course things cannot be so easy so you know this is just what we have to deal with. All the time people find you know more systems where we know how to do this in a positive definite way. Okay so that that's what we have to do and there are many different ways to do it. You can apply all kinds of transformations or or tricks here to to get something that works. So I just mentioned a couple of ways here so word lines that's basically you know the imaginary time path integral which I will talk about a little bit. You can do the stochastic series expansion and for fermions where you know these methods in general have a sign problems. There are at least in some cases some other kinds of methods that that can deal with these systems but in general fermions are still of course difficult. Okay so this is if we want to do you know finite temperature calculations. So if we want to just go to the ground state there is an alternative or several alternatives I should say. So you can do projection so that that's basically starting from a trial state which by the way is a bad name. I will say a few more things about that later but this is what people normally call it in this field. So you take a trial state which basically could be anything or almost anything and then you apply some operator with. When some parameter is large in this case this power M or in this case this number beta. When those are large one can show that you just get the ground state filtered out. I will say more about that in the next lecture. Okay so in the end it turns out that these two are quite similar. So I will talk about this kind of finite temperature methods today. This first lecture and this in the second lecture and in the second lecture I will also show you explicitly that they are almost the same in terms of technical implementations. It's basically just a question of a boundary effect. Here you have a kind of periodic time boundaries and here you have an open time boundaries. But okay what that means will be more clear in a moment. Okay so these things can be done particularly efficiently for spin one half models. You know these Heisenberg models and extended kinds of Heisenberg models. Yes please. Okay so I could answer that now but since that's the topic of next lecture I will leave it the next lecture. But I can just say that okay there's under some conditions. If the eigen value of the lowest state is the largest in magnitude you get the ground state. Otherwise you could also in principle get the highest state. But you can always fix that with some constant. But I will return to that question in the next lecture. Okay so actually often we write the Hamiltonian in this way. You will see that later today. Then you can recognize this as a singlet projector. So the Heisenberg model is basically a sum of singlet projectors. And those have some very nice properties that make it easy. Particularly easy to study these spin one half models. And that's nice because there are also many many interesting open exciting problems for spin one half models. Even for bipartite lattices. There was some time some years ago everybody said everything that has been done with. Everything that can be done without the sign problem has already been done. And that's still not true. It wasn't true 20 years ago when people told me that and it's not true today. So there's actually a huge amount of stuff to do. And we can do it really well because we can go to big lattices. Okay so I hope I can convey some of that as well. Okay now let's talk a little bit about path integrals then. So again we want to compute this thermal expectation value. Where beta is the inverse temperature and we also may also in this formulation. We may want to go to zero temperature. So the main difficulty here of course is how to deal with this exponential operator. Because it's not so easy in practice to take the exponential of a huge operator. If you have a big system the size of this Hilbert space is 2 to the n. Where n is the number of spins in this spin one half case. So how do we do that? So in the path integral what you do first you slice this partition function or this exponential operator. So first you just write it as a product of smaller operators. Meaning that there's a delta tau here which is beta divided by L. And you just do a product of those operators. And that's exact and it seems like we have done nothing. But now we can start to work with this because now if L is large. This delta tau is small and we can start to do some approximations. Okay and later we will get rid of the approximations. But just to get started we do some approximations. Okay so what we do first is we take a complete set of states. And then insert between all these exponential operators. So now you start to see where Monte Carlo can come into play here. Because we have a lot of states here and there's a lot of summations. And of course that's exactly what Monte Carlo can do very well. To do sums over a huge configuration spaces. So eventually there will be some of these states done stochastically. Anyway this is still completely exact here. And we can choose whatever basis which is practical for us to work with here. So for spin systems of the kind I'm discussing here. You can just think of the up and down z component basis. Okay so but now we still cannot easily work with this. Because these are still big matrices. And now one can proceed in different ways. Many of you probably have seen something like this before using the Trotter approximation. I want to use something which is easier. So imagine that this delta tau is very small. Then we can just do a Taylor expansion of these exponentials. And we just keep the leading order, the linear order in tau contributions. So then we get something like that. Okay so that of course has a big error. The error in each term here is of the order of delta tau. Because it's delta tau squared in each factor. But we have L factors. And if you look at that that cancels one of those. So it's order delta tau. But it turns out that we can actually very easily. And this was realized only maybe 20 years ago. So that one can actually very easily or pretty easily let's say. Take the delta tau goes to zero limit in the program itself. So that one can eventually get rid of the error. But for now let's just keep a finite delta tau and see what we can do with this. Okay and let's do an example then it's normally clearer. So let's look at a system of hardcore bosons. Meaning that you have a lattice with bosons. But on each site you can have only zero or one boson. So this is just hopping here, boson hopping. And this is actually equivalent to the XY model. I'm sure many of you have seen this as well. So you can just do a mapping where the empty site corresponds to say minus one half. And the occupied site corresponds to plus one half. And then these hopping operators just correspond to these spin raising and lowering operators. There's just a factor of two coming in front. And you can identify this with the X, X, Y, Y part of the interaction. So this is the equivalent of the spin one half XY model. Okay, so let's see what we can do with that. So this was the expression we had. And now if the Hamiltonian is this Hamiltonian, which consists of all these raising and lowering operators, which have the effect, well let's look at it for bosons. It's a little bit easier just in terms of pictures. So we have these bosons that are jumping between their neighbors. So now it's useful to do things graphically. So what I do is this state here is represented down here. So this is the state zero indicated there. I have three bosons here on 12 sites or whatever it is. And then the next matrix element here corresponds to the following, or the next state, the following line and so on. Okay, so you see we have the same state here and here. So we have the same state here and here. Okay, and what can happen here? Well, this Hamiltonian can jump these bosons. So given picture here corresponding to a given term here, only one thing can happen, so to speak, at the time if I really break down this into all possible paths. So at each step here, one term of the Hamiltonian is active. So here, for example, one boson is jumping like there. Occasionally none of them is active because we also have the one here. So you see there are some cases where nothing happens between two slices. It just stays as it is. Okay, so this is the graphical representation of one term of this thing when we replace H by the whole sum of terms. Okay, so let's focus on the left one for now. So you have a starting state and then the bosons fluctuate and then eventually they come back again to where they were. And now if you look at the weight here. So this product of matrix elements, that's the weight of a given configuration. And you can see if nothing happens, the weight is one on each slice. And if something happens, it's just delta tau. Well, it's minus delta tau times the matrix element of the particular term being active. That matrix element is always one because it's just hopping one boson and there's a minus sign which cancels that minus sign. So the weight is just the number of these jumps that you have in the configuration. That's the path weight. I'm running out of batteries here so I'm going to try to change the batteries while I continue talking. Now if you look at the right picture there, that's an interesting concept that that picture is illustrating, namely the winding number. So if you have periodic boundary conditions, then one boson can go out in one end of the system and it can enter at the other boundary. So we have a ring, the system is forming a ring. And then we can have cyclic permutations of the system. So what happens here is that here... Oops, hit the wrong button here. Sorry about that. So here we still have the same configuration here and here because these bosons are indistinguishable. If I had different colors of these bosons, for example, then I would have permuted my colors here so then this would not be allowed because I have to have exactly the same state here and here. But because these are indistinguishable, this is still exactly the same state here as here. Although the bosons have permuted during their travel in imaginary time. So we call this by the way the imaginary time direction because this e to the minus beta h is formally the Schrodinger time evolution operator in imaginary time. So each e to the minus delta h propagates imaginary time delta. Okay, so these configurations which wrap around the boundaries like that are said to have winding numbers. So this has winding number zero and this has winding number one. And this winding number is actually an important concept that will also appear in the stochastic series matter and it allows us to calculate, for example, in the case of bosons, the superfluid stiffness of the system. So the ability of the system to do a lot of winding is related to its superfluid stiffness. And in the spin system that corresponds to what's called the spin stiffness, which is basically corresponding to an elastic modulus of the system, the energy cost of twisting the spin direction in an ordered system. Okay, I will not say much more about winding numbers in terms of their physical significance, but basically this is an important thing to sample because in the end we want to sample these paths stochastically and we would like to make sure, at least if we want to measure the superfluid density, that we can also generate fluctuations in this winding number. So the kind of moves that one could imagine using for sampling here would be some local deformations of these world lines, so these are what's called world lines. So you can introduce and remove kinks in various ways and that changes the weight because this Nk, which is the number of kinks, is decreasing or increasing and then that you can use in a metropolis, accept, reject, probability. But these local updates, you see, they cannot really change this winding number. They cannot change the winding number because they only act locally, so they cannot make these global permutations. So to do changes in the winding numbers, one needs some global updates, so that will be an important issue. Okay, so how about expectation values? So now we were just, for now, sampling basically the normalization of a system or the partition function, I should say, in this case. So how about some expectation value? Well, if you have an expectation value, we have a very similar thing here, so now it's one over z. Z, we just discussed. And now there's a similar thing here, but we have one change in one of these matrix elements here where the operator A that we are interested in is located. And again, we want to write the whole thing in a suitable form on the color sampling. So we have already talked about the weight and now we should see what is this estimator for some physical quantities. Now it's very easy if you have a quantity which is diagonal in the occupation number basis because then this just becomes a number when you act on the state. And of course you can put the operator in any of these matrix elements here. This is completely cyclically invariant, as you saw from the pictures as well, or at least can imagine from the pictures. So we can just average that over all these time slices of the system if we like. So that's very easy. But sometimes of course we're interested in off diagonal operators as well and then it's more tricky and I will not discuss the most general case. I will just talk about one case which is the simplest one. Namely if the operator that you are interested in is some part of the Hamiltonian. So for example if we are interested in the kinetic energy which in the case of the bosons we are looking at is actually the full energy because we haven't introduced any interactions yet. Then what do you do? Well so now you have this kinetic energy here and it turns out if you just even neglect the exponential operator if delta is small and so you just consider having k here that doesn't change the order of the error of the whole thing so it's perfectly fine to do. So then what we have here is everything exactly the same as in the way except the last matrix element where we just have this k and let's even just consider one particular hopping term between two sides. So then we have one matrix element like that and the rest is the same as before. So then we can multiply and divide that with the weight that was present in the partition function. So then everything cancels out except these two. And here I have again you put in the expression to order delta tau. Now you can see that this has a very easy evaluates to something very easy. So either if you look at some where the boson jumps exactly on the side that you are considering so then this matrix element is non-zero and it's just one and in that case you also get a contribution from that one from exactly the same term in that one so that just becomes one of a delta but if there's nothing there, if there's no kink there then this is zero so you just get zero. So that means that to measure the kinetic energy what you have to do is just look at your picture here or of course you have it in the computer and count how many of these jumps do you have in the configuration. That's related to the kinetic energy you put in all factors as you should and so on then the total kinetic energy is just going to be minus the number of those jumps divided by beta. So what that then implies is since we know that the kinetic energy is proportional to the system size and it means that the number of these jumps should always be proportional to beta times n. So that means actually the density of these jumps if you think of the density, the probability of having it on one bond you get by just dividing by beta and n so the density is some constant, the density of jumps. So that's an important concept. So for example what that means again if I go to the picture is that if I make delta smaller and smaller I get more and more slices but at some point the number of these kinks will not grow anymore. The total number of kinks will be constant and on most of the slices nothing will happen. So you will just have a relatively speaking small number of events and a lot of time nothing happens. So I hope that's clear. Is that clear? Yeah please. Okay so this was one over delta and then that should be cancelled by something. Let's see now. Yeah it's okay. Let's see now where is that cancelled out. Yeah so you're right it should be, oh it's because that's where the beta comes from because this was on, if I put it on one of these so I'm again doing the trick that I can average it over all the places. So in the end I count all of them but then I have to divide by how many I have and the number I have is beta divided by delta tau. So that will be cancelled out. Yeah it's good point. It may look confusing at first but that's exactly because eventually the probability of having a kink will be very small. So the probability of having a kink will also be proportional to one over delta. I think you can be proportional to delta. So when you have a kink you will get a contribution one over delta but it's cancelled out by the fact that the probability of having one is very small on the order of delta. Thanks. Okay so now okay so I will just quickly say what happens if you include interactions here. So now if the Hamiltonian has two part, kinetic and potential part, then I can actually make an approximation so the potential part is V. So this would be exact for normal numbers but since these are operators that don't commute you make an error here but that's fine because it goes away in the continuum limit and then the potential part is diagonal in this basis so it comes outside, it's just a number and then if we do this at all time slices we get exactly what we had before for the weight times an exponential of all these interactions on the slices. So then when we try to do something to the word lines we have to take into account how the potential energy changes as well. Okay so how exactly do that let me not get to because I just want to basically introduce the formal aspects of the path integral so that we can make some comparison with the stochastic series method where I will talk about things in more detail but I just want to point out still that you can take the continuum limit so then the configurations start to look like this there's no slices anymore there's just a continuum of lines which have straight segments and horizontal segments where a boson suddenly jumps to its neighbor and the density of these jumps will be proportional to the kinetic energy. So that also means that the storage requirements of this method may not be as bad as you may have imagined initially so if we look here it looks like okay we should store information similar to an icing model and the lattice size is the number of slices and the number of lattice sites and of course if we make the number of slices go to infinity that looks like it should be not possible but in fact it is possible because when we go to the continuum limit we can actually just store the location of these events where something happens so we just need to keep a list of events. Okay and then we of course need to do some updates let me just not talk much about how we update these word lines but basically again we would introduce kinks and remove kinks in a local way but actually nobody does that because people have developed very smart ways to make very big changes in these things and those are called loop updates we will talk about loop updates in the context of the SSC method in a while so I will not say what they are but basically loop updates cannot just change things locally but they can basically change big pieces of the configuration at the same time, yes. Actually you will sample over all of those so all possible paths are allowed including what you say I mean in some site maybe there is no kink at all it could just go straight and in other ones there is a lot and when you do sampling of them you will get the important ones automatically so then of course on average it will be the same but there will be some fluctuations I don't know if that's exactly answering your question but you do it all at the same time so it's not like you sample one site at the time you basically, it's like the icing model if you simulate the icing model you make changes in many different places and in the end the whole configuration changes and you do that here you can do some series of small updates considered everywhere but you can do some more drastic changes as well but it may become more clear when we discuss the SSC method where I think many things are actually a bit easier okay so let's switch to the series expansion representation then so what we do then is we start from actually a Taylor expansion of this exponential operator that we have to deal with and so what I've done here again I look at the partition function Taylor expand the exponential and I also express the trace as a sum of a diagonal matrix element in whatever basis is convenient for us to use and now I want to deal with these powers of the Hamiltonian and the Hamiltonian of course consists of many terms and I'm going to break up the power of H into sequences of all possible combinations of these operators so formally to denote those I introduce an index sequence so this is a sequence of indices so you see I have a power n so I want to have a n operators and these a's here they refer to these indices here which normally would denote where on the lattice the operator sits but it can be further decomposed as well but there are some terms and there are some indices referring to those so this index again goes between 1 to m where m is just the number of different operators you have okay so now if we break H to the n into all possible strings of these individual terms I can write the thing as this so now I just have another sum here which is summing over all these index sequences which is the same as then taking the power H to the n so now you see what's done here is you have some state and then there are operators acting on that state and eventually well of course when the operators act on the state the state can change but eventually the state has to be coming back to where it started so again this is a kind of periodic structure here now when we decompose the Hamiltonian into operators here what we require in this context and this is basically automatic for the spin one-half models I'm talking about but it's not always automatic you may have to think a little bit but it's always easy that when you act with one term on one of these basis states you should get another basis state it could be the same as well but the point is that there should not be a sum of many states on this side just one that's how these operators are defined so you can always split up the operations of H into such terms then there will be actually something very similar as in the path integral here you have a state which is evolving in time so you can think of this dimension here as time as well and there's a formal relationship between them but I will not talk about that but if we again look at the hardcore bosons then when these operators act if they are allowed to act so actually not all strings are allowed because some of the strings could try to jump a boson on top of another boson that's not legal and also we have to satisfy this periodicity here so not all the paths or not all the operator strings are allowed but for those that are allowed in the case of the hardcore boson system the matrix element is always one so the weight of the path is just beta to the n divided by n factorial which came from the Taylor expansion okay now we can really make the relationships to the path integral clear by doing a bit more manipulation which you don't really have to do if you make a program but it's useful to see define a partially propagated states I take this state and just propagate it with the first p operators and then I call that alpha p so that's some state that is in general different from that one and then what I do is I just insert those states between the Hamiltonians and you know that doesn't do anything it's just a notational thing because H1 on A0 that is H1 so if we stick that oh sorry there's a typo that should be H1 as well and this should be a typo here but I should correct that before we post it online anyway so then it looks even more like like the path integral so we can actually draw identical pictures so I will not repeat it here I could draw exactly identical pictures to what we do for word lines for path integrals although it looks like it's a bit different because the starting point was different and there's some different factors and so on but in the end there's a clear relationship between these two and they are equivalent in some mathematical sense but in practice it's often easier to deal with it in this way let me derive some estimator here let's derive the estimator for the energy because that will tell us something interesting and important again so if I want to calculate I haven't told you how to do sampling yet but I can just formally look at things if I want to calculate the expectation value of the Hamiltonian I have something exactly as I had in the partition function but now I have another H there so I had H to the N from the Taylor expansion and another H just from my expectation value of H and then I can just do a simple relabeling so I say this is H to the N plus 1 so let's just call N plus 1 N and just relabel these sums then I have to correct in this factor because if I put N plus 1 there I have to divide by N and multiply by N plus 1 sorry, divide by beta and multiply by N plus 1 up here but when I relabel that to N it's just N divided by beta and there should be a minus sign as well because there's a minus beta here so then we can easily see that okay, now this is exactly the same as the weight was in the partition function except that there's this N divided by beta there so then this whole expression is exactly of the Monte Carlo form and we can evaluate the energy just as the expectation value of N divided by beta so if we imagine that in the simulation this N is fluctuating we are sampling also overall these powers we should actually measure the expectation value of that power and it's also actually very easy to just take the temperature derivative of this expression and calculate the specific heat and it turns out the specific heat is a kind of fluctuation of N so that would be the variance of N here but it's minus N so what follows from that if you think that the specific heat at least if you go to zero temperature the specific heat is normally zero so if this is zero it means that the variance equals basically the average average times beta so from this follows that if we think of okay when we do this Taylor expansion how high orders do we have to go to the order we have to go to for this to converge is actually because the energy is proportional to N so then that means that the expectation value will be proportional to beta times N and from this if we set the temperature to zero at least it means that the width of the distribution is square root of beta N okay so that may now look crazy that we are doing Taylor expansions to the order of beta N because beta could be easily let's say 10,000 and N could be 10,000 so then we are doing Taylor expansions to order 10 million and it's actually not too hard to do that stochastically and you may worry about some convergence properties and so on here but trust me there's nothing to worry about because we are working on a lattice the spectrum is bounded that means this Taylor expansion will always converge for finite beta, finite lattice the Taylor expansion always converges and this is the order of N for which it converges actually this let me show something simpler here which you can relate to so if I just take a number e to the x and I write the Taylor expansion of that just N x to the N N factorial and I just call this I just call this W of N so it's the weight of the N term in the Taylor expansion of this what does this look like well let's say that x is relatively large then it's easier to draw it actually will look something like that it's actually a binomial distribution of course and this average N is x and the width of the distribution is square root of x so actually this is like a Poisson distribution it's exactly the Poisson distribution so this has the probabilities of that so what we are doing up here is just doing this but for non-commuting operators so we don't have just numbers we have to deal with these strings of operators so this property somehow still goes through in some average sense that this average length of the sequence is related to the average of our operator so that's quite interesting so again, see is the specific heat well I think you know the physical meaning of that so it was derived here just by in this expression if you just take the temperature derivative this will actually come out automatically you can try that so normally unless we have some strange system the specific heat should go to zero when the temperature goes to zero so that's why I'm saying this here it's not completely true if the temperature if c is not zero but it's almost true oh that's the system size yeah I should have said that n is the total number of spins or lattice sites in my system yep sorry I should have I think I said it before I should write it out somewhere okay let me manipulate this a little bit more before we do something with it and actually I'm going much slower than I had hoped so let's see how far we can actually get okay so if we do sampling of this thing n is going to fluctuate so we have to have some Monte Carlo updates which increase or decrease n and as you can imagine it's a little bit annoying to work with a configuration space where the configuration itself is fluctuating although it actually can be done but it's a bit easier if we do a scheme which has a fixed string length and that can be actually done very easily so what we can do is we can truncate the expansion at some point and you know you can see it from here since it's very similar from this if we okay this is of course really discrete I just drew it as a continuous line this eventually decays exponentially so if I do you know some truncation here at L you know nobody is going to carry I'm never going to see you know the difference in practice because it's exponentially small so we can actually truncate the expansion at some also in the quantum case at some n equals L and then we can actually do a trick so if I have let's say here n is 10 and I have these operators okay here I apologize I call it M I have somehow mixed pictures from different presentations and not noticed but okay this L should be this M here let's say I decide the max is going to be 14 then what I can do is I can actually put in some unit operators in this string just to make it look like it has length 14 so I can think of the unit operator I or 1 or whatever we call it as in a similar way as I think of the Hamiltonian terms so I can just define H0 equals I and then I can formally regard it as part of the sampling of the terms of the Hamiltonian but then I do an over counting because this should just be one term it should really be just that one and if I distribute in all possible ways these unit operators then I have to compensate for that over counting and over the counting is just of course the combinatorial factor L over N L choose N so I just have to compensate for that so then this becomes my partition function another type of here this should be N not M right so this came from the compensating factors so the N factorial is gone and there's an L factorial instead which is just a constant and then the summation over N formally went away it's just that now N here has a different meaning N is just the number of these operators which is not zero so N is still there but it's not formally a sum here it's something which is implicit you just have to know each instance how many actual operators are here if you don't count these unit operators so this makes the configuration space a little bit easier and it turns out and I will show that this cut off here if you think of it here you can actually choose it automatically you don't have to worry about the program can detect what it should be to be big enough okay so now we are ready to actually do an implementation which I think will make the whole thing a bit clearer so let's look at the 2D Heisenberg model and I'm drawing pictures for the 2D Heisenberg model but actually I'm formulating it in a way which is completely independent of the lattice so I write my Hamiltonian as a sum over bonds so in the 2D case which I use for illustration I have numbered my sites in some way in some natural scheme and I have numbered my bonds the bond consists of two interacting spins so bond one connects spins one and two and so on so I can write my Hamiltonian as a sum of these bonds instead of the sum of spins if I introduce some mapping BJB to correspond to these so for example bond 23 I 23 would be 7 and J would be 11 so it refers to those sites that are connected by the bonds so if you have defined your lattice in the computer by a list of sites connected by bonds then you can define the Hamiltonian like this and the method I will talk to you about can be used for any model essentially without changes because it only needs that list of sites connected by bonds but let's have the 2D case in mind because it's more concrete that way okay so you remember before I told you that we can write the Hamiltonian as in terms of these operators one quarter minus si.sj so I have a minus sign here and a minus sign here so this is what I'm going to do now as well but I only wrote it down for the diagonal and off diagonal terms explicitly so the Hamiltonian as a diagonal term and now this is this one actually yeah so if I define the Hamiltonian with a minus j then you see that this is the diagonal part of it up to a constant which we don't care about of course and then the off diagonal part the si interaction is again these racing and lowering operators and now that comes with a minus sign because I had pulled a minus sign out in front so for the Heisenberg model I define for each bond two types of operators one is the diagonal and the two is the off diagonal part okay so this NB corresponds to the number of bonds and I have for each bond I have two types of operators okay the reason I'm adding this or writing the interaction in that way which has the effect of that one quarter there is that then the matrix elements are all the same all the non-cero matrix elements are the same and these are all the non-cero matrix elements so there are only non-cero matrix elements of these staggered configurations when the spins are up down or down up on in the bra and the cat and they are all one half and that's something which will make the scheme much much easier than if it's not the case if it's not the case one can still do what I'm telling you about here but it's a little bit more complicated but still relatively easy but I'm going to just talk about the easiest case okay so then formally I write the partition function again like this but this is the product now of this string of operators but I refer to the operators with two indices because of these here so if this A is one it's diagonal if A is two it's off diagonal and then I pulled out the minus sign here so the minus sign cancel except for the cancel except for the sign in here so for each off diagonal term if I count how many off diagonal terms I have here and I call that N2 then I have a minus one to the N2 here so that's what eventually can cause a sign problem because we would like everything to be positive and everything is positive here the only thing that could be negative is that sign there but we will see later on that for bipartite lattices it's nothing to worry about okay and again we do the fixed length scheme we have the propagated states and now I can draw a picture which hopefully makes things clearer so here I have an eight side system and this is now some representation that you can use in the computer as well and there's a graph going with it so I have eight spins I denote them by sigma and they take values plus or minus one although they are plus minus one half but I can just call them plus minus one and then this is the starting state I call it just alpha here and it's some up and down spins, red and white and then I have some operators acting on the state and now I denote the operators by open boxes for the diagonal operators and these black solid boxes for off diagonal operators and you see here here's the diagonal operator acting nothing happens, the state is the same an off diagonal acts and then these two spins are flipped here there is nothing so that corresponds to one of these unit operators that are just sitting there doing nothing and then some more things happen and eventually there are four spin flipping events here and the state is propagated back exactly to itself so this is a valid configuration of the SSC for this case and the way it's represented in the computer is that you have to store this state here which is the start or end state it's the same and then you have to store this A and B which tell you what type of operator and where it's located and you can actually do that in a more compact way which not only saves memory but it's just well it saves a bit time as well I think so you can just pack this into one single entry but before I talk about that let me again mention this sign problem here so this sign if you have a bipartite lattice it goes away and this is an example of that that you have to have an even number of spin flips to be able to propagate the state back to itself if you have an odd number you cannot satisfy the periodicity condition so that's not true if you have a triangular lattice so for example on a single triangle I show a case here a spin configuration and then here I just draw it in this direction I act on these two spins, I flip them then I flip those two, those two and now I'm back to the original one and I did it with three operations so that means minus one to three that's minus one so that configuration comes with a minus sign so we cannot do frustrated systems with this method because it causes sign problems formally we can deal with it but in practice it doesn't work so let's stick to bipartite lattices where if you think a little bit it's true that we always need an even number of operators there okay so this index which is used in the program which is available online that corresponds to two times B plus A minus one so if that index is even it corresponds to a diagonal operator and then you just divide it by two to get the location if it's odd it's an off diagonal and then you I guess subtract one and divide by two to get the location okay so this is how the configurations can be represented just by a spin state and this string of numbers referring to operators and we don't need to store all these states in between we will later see that we occasionally need to store a bit more information but basically if you have this information you can recreate all the spin states when you need them okay so here you can see that this is a complete discreet representation of the quantum statistical mechanics problem but it represents the imaginary time continuum because there's no approximation here this is completely exact well we did the truncation but I would argue that that's not an approximation because it can be arbitrarily good and never visible in practice and again in principle one can relate this dimension here to imaginary time so we call it also here the time dimension so this is an integer representation of the time continuum so it's good to have the time continuum because everything is exact so it's good to have a discreet representation instead of working with continuous times because it's just more practical and faster in the computer alright so another thing to note here is that these are like the events that we talked about in the path integral and between events nothing really happens so in some part of the program we actually want to store these events and how they are related to each other and we will do that by links so we will call these operators and the spins that they are connected to we will call those vertices so these are the possible vertices that we have in the Heisenberg model again they correspond to those non-seromatrix elements so you see it's up and down on either side of the operator and we call these the legs of the vertex and I number the legs in this way for convenience so now a more compact storage of that whole spin configuration I showed you before is this one I have the starting state which is the same as the end state and then I have these vertices and I have erased the spins that are not changed between events and in the computer these will be represented by links because it turns out that for some of the updates that I want to do on this system I will want to move among these vertices in the way they are connected to each other so you can think of them really as connected to each other by these lines and these connections can be stored in the computer I just want to really go into details because it takes a long time and it's just confusing but basically there's some table that will allow you to jump for example from here to here and then you may want to move there and jump and go across the boundaries and all those things can be done in some simple tables so those will be used for what we call loop updates of the system okay so how do we actually do Monte Carlo sampling here well I indicated already that we will move around in this configuration but we actually have two types of updates in this case and the first one is well let me first say in general so this was the configuration weight okay and this is just a picture of a configuration and as you remember now from Werner's lecture yesterday you have an accept probability which includes the weight so this is the metropolis basically which includes the weight ratios if I change the configuration to something with primes on it I have the weight ratio but also what I think Werner called it a priori probabilities the a priori probabilities should also come in there so this is the a priori probability of going from this configuration to that one or from that one to that one if you do something like an icing model simulation where the a priori step corresponds to just choosing a spin at random for example the one to flip then that probability is always the same so those will just cancel out here and in most Monte Carlo meters those will just cancel out but here as you will see it's not always true so we actually need to keep that okay that will be clear when we consider what we actually do so let's talk about what we call the diagonal update so you see if we want to change something in this configuration what can we do well these diagonal operators they don't do anything to the spin they are not associated with that many constraints the only constraint is that diagonal operator can only appear between unequal spins so I could not put it here but I could put it there for example and I can always remove a diagonal operator that still leaves the configuration and a loud one so I can remove diagonal operators and I can insert for example one's diagonal operators where there is nothing formally introduced unit operators so I could introduce something there so the first step of the simulation is just to insert and remove diagonal operators so we can go from zero zero which corresponds to the nothing here to some diagonal operator at bond B okay so to do that we need to know what state we have because otherwise we don't know if it's allowed to put it in there or not again we cannot put it between parallel spins and so to have those states we can start here and then we can start to do that here and then we can go up and whenever we encounter these off diagonal operators we just flip the spin so we just one by one generate these states and whenever we can we do that kind of update okay but now there's an important thing to notice here if you insert an operator in principle you can insert it at N different places okay some of them will eventually not be legal but a priori there are N different places to put them but if you remove one there's only one way to remove it so you don't have balance there you have N ways of putting it in but only one way of taking it out and that's why to get balance detailed balance you have to compensate by also including those select probabilities so here the select probability is 1 over N because there are N places to put it going this way the select probability is just one because there's just one way to do it okay so that is said here and the weight ratios are easy to compute because N here just increases or decreases by one so in the end we have some acceptance probabilities of inserting an operator or removing an operator and those are very easy to do again sequentially just go up and try it wherever possible and if it's not possible it means there's one of those guys there's then you just flip the corresponding spins okay so I have a little piece of pseudo code how to do it too I think I should not well we have a little time so maybe I can say just a little bit so you loop over so it's basically what I just said you loop over all those locations in the space, in the time space and then you look at this S which encodes the type and location of the operator so if it's zero that means there's nothing there so then you can generate a random bond and try to insert the operator at that bond but if the spins are equal there it's illegal so then you just cycle that means that you go to the next point of the loop so you did nothing but if that's not the case then you can insert it by generating a random number and if that is successful then you update this holder, the index corresponding to the operator and you increase N okay if there was something there already you try to take it out with some probability and you decrease N by one and if none of those hold that means there's an off diagonal operator there and the bond you extract from this S and then you flip the spins at that bond so you have made these lists of sites connected by bonds and you work with those so that's very easy to do the diagonal updates okay but that's just part of the story but that's already a big part because what that will do it will fluctuate N, the actual number of operators and it will put in these diagonal operators at different places so it will sort of sample locations of operators, operators and their locations so the next update we are going to consider will actually not change the number or location of the operators it will just change the types of operators meaning diagonal ones can become off diagonal ones and vice versa so those we call off diagonal operator updates the easiest way to do that would be just to consider pairs of operators where that's allowed so for example if you look at these two off diagonal operators I can just change them to diagonal operators and in the process the spins between them will change but this is still a completely legal configuration but you can only do that when there is nothing between here, so for example in this case if I try to do that you see that there's something in between which would no longer be legal if I change the spins along those lines so you would have to go and check where that's allowed and you could do it actually in this case we can still do it if we consider the periodic boundary conditions and make the changes under and above instead then it's allowed in this case but those are local operators and they are not local updates and they are not very efficient and in particular these winding numbers that we talked about cannot be changed these correspond more or less to the small kinks that we removed and introduced in the path integral case okay but here comes the nice thing one can do something called loop updates or in this case operator loop updates and that's as follows so you start somewhere let's say you start here and then you always go to the next spin the next leg on the same vertex and then you jump to the one where it's linked to in this representation where you think of these as physical links you jump there and then again you go to the neighboring spin you move along the links and you see that eventually that will trace out a closed loop in the system it will always no matter where you start you can do this and it will make a closed loop for you and then what you can do you can flip all the spins along the loop so you should think of it as flipping all these spins that you hit including you flip all the spins here in between although you don't see them because we don't need to represent them but they are still there if you consider this representation with all the states but in the process of changing the spins also the operators change so whenever you hit an operator you change two spins the type of the operator changes as well because this was diagonal but now when you flip the spins it's still okay if the operator is off diagonal actually it's a bit redundant to even draw these operators here I don't even have to do it because just the four spins on the vertex completely specify what it is it's just a little bit sort of added clarity hopefully what it's doing and you should also know that some operators maybe hit twice and then the type doesn't change because it changes twice so first it changes here to off diagonal but then it changed back to diagonal when it was hit again so this is really very efficient these loop updates because the loops can become very long and they can also change these winding numbers and so on so this is what makes the scheme really efficient okay so I also have some pseudo code for these loop updates but since I only have six minutes left I don't want to do it in really in detail but the point is that one can make all the loops so this is just showing one loop but if I keep track of where I have been in this space if I start somewhere else and make another loop and so on eventually I will have covered the whole system by loops and when I construct them I also flip them with probability one half so if you have heard about the Swensen Wang cluster algorithm which you know it's very common for simulating the icing model this is very similar to the Swensen Wang cluster update but it's a loop update but the principle is exactly the same and Werner will talk about that tomorrow I guess he will talk about Swensen Wang so you can keep in mind that you have already heard about a cluster type algorithm for a quantum spin system okay so here this is talking about this pseudo code about how you cover the whole system with loops and you put in some information in this so I have made this list that links these vertices and that's the list that you use to move in the space and basically you destroy the list as you go along it in it and put in some markers to indicate if the loop has been flipped or not and then you later will use that information to possibly change the spin state as well and okay this looks a bit messy and I don't really have time to do it in detail but it turns out if you want to do this really efficiently and also quite simply you can actually work with bit operations because when you have represented the operators in the way I do with even and odd numbers corresponding to diagonal and off diagonal operators you can just change from diagonal and vice versa by flipping the zero bit of that thing so that's a very fast operation and also some movements in this list can be done with bit flips but I think I will not talk about that in detail it just gets boring but there are papers of course discussing these things and I will also be here all weeks if somebody is interested I can talk more about it and also there are these programs online which I will point you to in a moment okay so after the loop updates have been done this spin state that we store the starting and end state which are the same that also can possibly be changed because the loops can actually go across that state and change the spin so that has to be done explicitly because that state we have to store and it has to be up to date and okay I also discuss it how to construct the linked vertex list but this is more for people who want to look at it later and try to understand what it's doing it's not hard but it's you know you can figure it out yourself if you have this kind of storage what you would actually do but it's available for you okay let me just now show some illustrations of how this is working I mentioned that this cutoff L can be done by the program automatically so we have the equilibration part of the simulation the burn in somebody called it yesterday so we can just start with basically an empty string and some random state and some small n and then we keep track of this n you know this order of the Taylor expansion and if okay so we have all chosen some small L as well so if n comes close to L that means L is not enough so then we just increase L so we can say that L should always be something like n times one and a third of let's say the maximum n that has happened so far so here's an example of that in an actual simulation of a 16 by 16 system here I started with L equals 20 I think so L is this solid black line and the actual n after each step is the black dots and the maximum n I reached so far is the red line and you see this increases very quickly in the beginning and then everything sort of flattens out and fluctuates shows the probability distribution of n in a long simulation after I had fixed L to whatever it was here and you see here that it fluctuates around 4800 well you see it here too and you know you never reach the maximum which is 6500 the maximum is you know out somewhere here clearly that doesn't introduce any error at all okay but you also have to make sure everything that these things are working as Werner mentioned yesterday as well things can in principle go wrong and of course you can make some coding errors you may have some long time scales that you don't know about but actually in this case the time scales are normally not that bad but still you have to check things carefully so one way to check is to check with exact diagonalization so here I show something where we did a 4x4 system exactly and in the 1D Heisenberg chain you can also compare with the beta ansats there's an exact solution so let me just flash you quickly the results so here the solid curve is the susceptibility versus temperature of this 4x4 lattice and the red dots are the simulation results and the error bars are so small that you really have to just take the difference of this and try it a lot to see the difference so the difference is on the order of 10 to the minus 6 and you can see that to that accuracy there's no difference just consistent with statistical errors that you have to compute of course and this was a really long simulation it took well several days or maybe weeks here are some results for some long 1D chains which illustrate so that you can get the ground state because we can get the ground state energy from the beta ansats these lines for two pretty long systems 1000, 4000 sites and as I increase beta and this is on a logarithmic scale you see that the energy becomes very consistent with those exact values so this definitely can give you unbiased results where there are only statistical errors no systematical errors okay just quickly this will be posted online so you can get the URL from here but basically I have a couple of programs with a little bit of instructions here which do Heisenberg 1D and 2D simulations I can just quickly show it just maybe flash it just to show you that actually it's not a very long program that's the main point I want to make so here's the actual quantum Monte Carlo program there's some instructions here some comments okay let me not even say what the different things are but okay there's the diagonal update that we discussed there's making the linked vertex list so you see that's maybe I don't know 40 lines there's the loop update which actually is very short this is just the loop update but okay it calls to smaller subroutines but those are very small so it's not much longer than that pseudo code that I showed you and then it's doing some measurements and writing some results and so on but this is the whole program it's maybe I don't know 300 lines and this is basically research quality code so you can actually even do some research probably with this code if you like okay I think I'm a little bit over time so 3 minutes over time so or just one I guess because we started 2 minutes late okay so thank you