  umbre                                                                                                                                                                                                                              ි වහ්අවියයක මා඼ක් අරය දැණර කිරයක මඩාරපරයක මදයක වොලුගු හැකකරයකතයකතයකරියක. වපිජිනකයක මෙ඾ාලුක අනමමීරයක. Here on this side we had models with independent entries and here we had models with rotational invariance and we argued that the intersection contains the Gaussian ensembles so the ensembles characterized by a joint distribution of the entries of this form well in general. So these models do have independent entries but also the property of rotational invariance and by virtue of the Rosenzweig and Porter theorem these are the only ensembles that enjoy both properties. Okay so remember that these ensembles can come in three varieties real symmetric complex Hermitian and quaternion self dual characterized by a Dyson index beta equal to one two or four so since we only have one ensemble in this intersection it makes sense to give some special names to those. These are names that you will find everywhere in the literature on the subject. So we have the G O E this stands for Gaussian orthogonal ensemble beta equals to one so a G O E is an ensemble of real symmetric matrices so beta equal to one with Gaussian probability distribution of the entries. Okay so the name is again slightly misleading so Gaussian orthogonal ensemble does not contain orthogonal matrices it contains real symmetric matrices that are diagonalized by orthogonal matrices. Okay it's called Gaussian orthogonal ensemble but it does not contain orthogonal matrices similarly we have Gaussian unitary ensemble G U E beta equals to two this ensemble contains complex Hermitian matrices which are diagonalized by unitary matrices and with independent Gaussian entries. And then of course you can guess we have a Gaussian simple symplectic ensemble for beta equals to four. Okay these are the only ensembles that lie in the intersection of these two categories models with independent entries and models with rotational invariance. Okay this classification scheme is very useful so every time you come across a random matrix with real spectrum you should be able immediately to locate it on this chart. Okay and you know from this chart that if your model has independent entries well then almost surely it does not have the property of rotational invariance. And if your model has a property of rotational invariance then almost surely it does not have independent entries. Okay now can I erase everything here? Good so let's go back to our 2 by 2 example because it gives us some good mileage on the general n by n case. So if you take a 2 by 2 G O E matrix then that's exactly the example we started off from. So this is a matrix of this type with the joint distribution of the three variables of the form exponential of minus one of trace of x square. Now the question you might have is we have the joint probability density of the three independent variables these are real numbers. So what is the joint probability density function of the two eigenvalues? The two eigenvalues are random variables they will be described by a joint probability density. Can we compute it starting from this input? Remember this morning we computed starting from here we computed the probability density function of the spacing but not of the individual eigenvalues separately. Now if you wanted really to compute this object in the long way let's say what you should do is compute this integral. You take the joint distribution of the entries and then you stick in two delta functions one that says that lambda one is its expression in terms of the entries. And then the other one that says that the second eigenvalue and here you put the expression of the two eigenvalues in terms of the of the entries. And then you perform this calculation this will be quite long and boring actually normally you don't do it this way because there are quicker ways to do this calculation. But if you wanted in principle this is the way to go. So if you carry out this calculation what you will find is that the joint probability density of the two eigenvalues has this form. It will be proportional to exponential minus one half lambda one square plus lambda two square. And then you will have an absolute value here downstairs of lambda two minus lambda one. Now this formula basically contains all we need to know about eigenvalues of random matrices. It has in itself all the information and all the important points that we will need in the in the following you see here you have Gaussian Gaussian weight. So if this term here was not there then the two eigenvalues will be just simply independent Gaussian Gaussian variables. Which means that basically the eigenvalues if this term was not there but even when when it is there it means that the eigenvalues don't like to stay very far away from the origin. So configurations in which you have two eigenvalues that are very large in absolute value are suppressed they have a very low probability. But on top of this of this object you have a term here that correlates the two eigenvalues. It is a term that is proportional to the distance between the two eigenvalues and now we are understanding why the eigenvalues repel each other. Because every configuration of the two eigenvalues where they get too close to each other have a vanishing probability. So here again so this information explains in some in some ways why we have level repulsion in a in a G O E. There is this term here lambda 2 minus lambda 1 in absolute value. Actually you can now from this joint distribution try to recompute the level spacing distribution it will be much quicker from this from this route. And you see here this object here will give rise to the factor S in the level spacing that is exactly responsible for the level repulsion. S times exponential of minus S square exponential minus S square comes from here but S is already in here spacing. Okay so how does this generalize to the case of n by n matrix. So here we have a 2 by 2 G O E how is the joint distribution of the n eigenvalues for an n by n G O E. Well I will not prove it but starting from this formula you can sort of imagine how to generalize it to the n by n case. So I will give you this result without proof. So we have that for Gaussian matrices. So let's say G O E, G U E and G S E n by n. So the corresponding result is the following the joint pdf of the n eigenvalues is proportional to exponential of minus 1 half summation I 1 to n. Lambda I square. So you have the Gaussian term that we we expect. It's the same is the analog of this object for 2 by 2. Times you have a product J less than K absolute value of lambda J minus lambda K to the power beta. So I recall that beta is equal to 1, 2, 4 for G O E, G U E and G S E. So this product if you write it explicitly is for n equal to 2 this will be lambda 2 minus lambda 1 or lambda 1 minus lambda 2 in absolute value. For n equals to 3 this would be lambda 1 minus lambda 2, lambda 1 minus lambda 3, lambda 2 minus lambda 3. So you see that this term here is quite interesting because it is formed by product of all pairs of eigenvalues. It means that if you have n eigenvalues here on this line this eigenvalue here the smallest one knows about all the others including the one that is the farthest from it. So this system of random variable is the opposite of independent random variables. It is a strongly correlated system where every random variable knows about all the others. This makes random matrix theory interesting in itself because it is a theory of strongly correlated random variables. So all the classical tools, low of large numbers, central limit theorem, everything you know from the study of IAD random variables stop working here. We need to devise new tools and clearly if you take beta equals to 1 and capital N equal to 2 then you recover this expression here. So I asked myself the question when did this formula first appear. So when was this actually discovered and there is some interesting historical observation. So this formula is a corollary of theorem 2 in a paper by HSU, don't know how to pronounce it, which was published in the unknowns of eugenics in 1939. That's a pretty scary title and you can find it in the handouts on page 8, 9, 8, 9 and 10 and on page 9 there is basically the disclaimer. So the editors of the unknowns of eugenics are basically saying well you know publication of this material online is for scholarly research purposes and it's not an endorsement or promotion of the views expressed in any of these articles or eugenics in general. This is the earliest appearance of this joint distribution. It's good to know for historical reasons. Okay, so now that we have something interesting here to work with, so we know at least that for the Gaussian ensemble we have the joint distribution of the eigenvalues. Actually I wanted us to make another observation which you can understand at the level of a 2 by 2 matrix. You see here you have 3 independent variables in the matrix. These are 3 independent entries and you only have 2 eigenvalues. In general this happens also in higher dimensions, so you have an order n square entries and an order n eigenvalues. So you see basically this object here that correlates the 2 eigenvalues can be seen as the price you have to pay in reducing the complexity from an order n square in random variables may be independent to an order n eigenvalues. So you decrease your complexity. You are basically throwing away a lot of random variables. But the price to pay is that now you are correlating them all. So even if you start with independent entries the eigenvalues will always be correlated. So most of the first part of this course will be about how to deal with ensembles of strongly correlated random variables for which the standard tools don't work. We need an entirely new set of tools. Good. Can I erase everything? Excellent. So now we ask what kind of quantities we can compute, what kind of quantities are interesting to compute in random matrix theory. So of course we start with the simplest. So quantities interest. Well suppose that we want to do like a numerical experiment which is also by the way in one of the handouts. So what I do is I generate an n by n real symmetric matrix. For example a GOE matrix but we can do it with any kind of probability distribution for the entries. So this matrix is real symmetric. We produce it, we produce one sample by sampling all the entries in the upper triangle independently. The diagonal with variance one and the off diagonal with variance one half. Remember you need to do that. If you take them all with the same variance this does not belong to the GOE. So you produce this sample. This is one line in Matlab. And then you compute the n eigenvalues. This will be n real numbers corresponding to the first sample. So I'm just doing it step by step like what you really have to do on a computer. First sample n eigenvalues then you do the same thing here. Second sample so the entries will be in general different and the eigenvalues will be in general different. n real numbers corresponding to the second sample and you do it m times. So this n by n and this is lambda one m lambda n m. You produced all all these numbers. So you produced n times m numbers real numbers at this operation. Now all you have to do is you put all these numbers together in a in a mega vector of size n times m. So you produce a mega vector lambda one one lambda n one up to lambda one m lambda n m. You have a mega vector of real numbers and now you produce a histogram a normalized histogram of of these numbers. So your histogram will be peaked around the position around which you find most of the time an eigenvalue. So if you do that for example for the GOE you will find a situation like this. You have zero and then you will have some sort of shape of your histogram like this. Actually you can find a more precise version of the sketch in the in the numerical hands up handouts. It will be on page five. This is basically the sketch that you will that you will get with a lot of wiggles if you try to do this experiment in practice. So there is no cheating here no fuzzy maths. You just in math lab it's probably six or seven lines. You find this lot of wiggles. So what kind of features do we do we have for this histogram. So this histogram will be like symmetric around zero. So it means that on average you are as likely to find a positive eigenvalue as you are finding a negative eigenvalue in all this in all this mega mega vector. The other the other feature is that here basically the your histogram drops to zero at the position around square root of 2n. Which means that if you take n like 100 it will be very unlikely to sample to have any entry in this in this vector that is larger than root of 200 which I have absolutely no idea how much is it. Okay. So this is also interesting because this edge point grows with n. So if you take a matrix 10 by 10 and a matrix 20 by 20 there are two differences. The first one is that in a matrix 20 by 20 you have 20 eigenvalues and in a matrix 10 by 10 you have 10 eigenvalues. This is trivial. Less trivial is the fact that in the 20 by 20 matrix the eigenvalues are bigger. So you got many many more but they are also bigger. They cover a wider fraction of the real axis. If you think about it this is less much less trivial to prove but it is true. You got more eigenvalues but they are also bigger. Okay. Now the shape of this histogram for m the number of samples very large is called basically the average spectral density. It has this object has a lot of names can be like average spectral density spectral density level density that it has a lot of names. It is it is the same the same object. It is the shape of this normalized histogram when you're not the number of samples is sufficiently large. Okay. And we call it row n of lambda. So row n of lambda which is the shape of this of this histogram in in words is the probability to find. An eigenvalue around the point lambda. So this is the the probability to find an eigenvalue around the point lambda an eigenvalue in this mega mega vector. So it is a probability regardless of which sample it came from. So we don't want we don't want the probability that the third sample that the second eigenvalue of the third sample was around the point lambda. We want the probability of finding an eigenvalue anywhere in this in the sequence and anywhere in this sorted order around point lambda. Okay. And this this is given by this histogram because we are putting all the eigenvalues together. So we forget where they where they come from. The average spectral density satisfies it is a it is a probability density. So it satisfies the normalization to one. So the area under this Instagram must be equal to one. Okay. So this is one object that we we know how to compute numerically. Now let's let's ask the question. Can we compute it analytically? Can I can I remove? Okay. So let's let's leave it there for the moment. So first remark. So the spectral density row or the average specs spectral density row and of lambda is a highly known universal quantity. And I explain what what this means. So this the spectral density this this the shape of this histogram is a highly known universal quantity. This means that if you change the probability distribution of the entries in general this shape will change and often it will change a lot. Okay. So every ensemble every joint distribution of the entries will have its own distribution of the of the eigenvalues. Now the question is can we compute can we compute it. So the problem that we have is given the joint PDF of the entries. Can we compute row n of lambda. I give you this object which defines the ensemble and I asked the question can I compute from this object the average spectral density. In principle if I do everything correctly I should be able to compare this result with numerical diagonalization and find a good agreement. If I did everything correctly because I can sample random matrices according to this distribution diagonalize them put everything together form the histogram and the result that I have from Matlab should match this function that I have computed. Okay. So in principle we have a tool to check if we are doing something good or bad. Okay. Now the situation becomes unfortunately a bit less less nice. So remember our original scheme matrices with independent entries matrices with rotational invariance and Gaussian ensembles. Now it turns out that random matrices work in the opposite way as you might expect from your experience on random variables like independent random variables in the case of independent random variables is the easiest among all random variables. If you have independence then that's an amazing property that allows you to do a lot of calculations. The case of random matrices is completely different. So matrices with independent entries are super hard. So we don't know nearly as much as we would like to about matrices with independent entries. So in general rho n of lambda the average spectral density cannot be computed. So the spectral density although you can simulate these matrices easily and produce this histogram easily for n equals 17 or 23 or 34 you don't know how to compute you don't have an expression we don't even have tools to compute this this object. So this is something for you to bring you know to bring this field forward we will need some some tools that maybe some of you can can create to compute this this object. What we can do is sometimes we can compute rho n to infinity of lambda sometimes so sometimes for very large random matrices. We have a few tools not so many that in some specific cases allow us to compute but I put it in inverted commas and I will explain why the density of eigenvalues the spectral density for n to infinity. This will be the topic of my second second week so in during my second week I will try to deal with this type of problems. The main the main tool that we we use in this case is replicas so tools that are borrowed from the physics of spin glasses. But the field of independent matrices with independent entries is open and in the sense that we don't we don't know much about about them. The case of matrices with rotational invariance is much nicer so in the case of matrices with rotational invariance the spectral density the other spectral density can be computed for finite n in closed form. And also by taking the limit n to infinity also for large large matrices clearly so this the main tool here will be the method of orthogonal polynomials which I'm going to describe for you tomorrow. If you use the method of orthogonal polynomials you will find let me see there should be here yes on page 20 so on page 20 of the numerical handout. I wrote like a small matlab code so you have a GUE so beta equals to Gaussian ensemble and here you have the spectral density the dots are the histogram there so numerical diagonalization of matrices and the continuous the continuous curve is the theoretical the theoretical result so we have an exact formula for the spectral density and the two match very. Very well so the topic of tomorrow's lecture would be to how to compute this isn't it nice I don't see enthusiastic faces it's good right so you just diagonalize matrices and then you have a function and it works come on be happy good. Okay so this is very super exciting. So before. So okay this was just a recap on the on the scheme. I just raised everything so tomorrow I will introduce the orthogonal polynomial technique. Next week I mean two weeks we will deal with matrices with independent entries using replicas. And now I wanted to introduce another technique that can be used to compute the spectral density sometimes in the limit and to infinity. Okay so you see n equal to two and n equal to infinity are often the two cases that we can we can deal with easily or more more easily and equals to because we can do calculations explicitly and maybe try to infer from what happens. For n equals to in the cases of larger and and n equals to infinity because sometimes there are simplifications that the take take place so now I will try to introduce a technique that is used in this in this situation. So first of all let me define the average spectral density I defined it in in words as the probability of finding and again value around the position lambda but how can we define it in mathematical terms well I will give you the definition and then. So you have a sum of delta functions where lambda must be equal to the position of each eigenvalue. And this the sum is average over the ensemble so in this is in theory let's see in practice how it works so suppose that you have the first sample of of your random matrices so this object here. Will produce a set of spikes at the location of each eigenvalue lambda one of the first sample lambda two of the first sample lambda n of the first sample okay. Just just this object here you just have to to sum all the spikes at the position of each eigenvalue. Now you do the same thing for the second sample and you will have spikes somewhere else in general you will have the same number of spikes one two three four five six seven three four five six and seven and down to the. And sample yeah I mean the the the point is that when you when you take this this average this function which is highly singular because it is basically a set of discrete. So sample by sample for finite and you can identify exactly the value at which one eigenvalue is the problem is that when you when you take the average. So taking the average means that basically you are putting all these spikes together and the density profile will be basically higher when you have when you have several spikes spikes around a certain position. So this object becomes a smooth function yeah so this this object becomes basically a smooth function of of lambda and so it makes sense to interpret it as a probability density function. So the probability that you will find one eigenvalue around so between lambda and lambda plus the lambda okay. Note that this this doesn't have it has nothing to do with the limit and to infinity this is still valid for for finite and what is what the problem is the is the average. So when you take when you take this in this direction a lot of samples in the limit and to infinity not the limit and to infinity okay. Good so the question is can we compute this this object for for a given matrix model in the limit and to infinity. Well we can sometimes using a technique that is called resolvent resolvent method which I'm trying to describe now so the resolvent method is quite interesting. So first I give you a definition so the definition is we define a certain function g n of Z which is defined as follows one over N the trace of one over Z times the identity minus the matrix H. Or if you want in terms of the eigenvalues is one over N summation I one to N one over Z minus lambda I okay. So we define this object X doesn't need to be a random matrix it can be a you know a fixed matrix you can define it for any fixed matrix. It is the trace of basically the inverse matrix of Z the identity minus X which in terms of the eigenvalues is written in this in this form. So so far there is nothing random the important thing is that Z is taken to be a complex variable. So Z belongs to the complex number minus so to which we have excluded the set of the eigenvalues because otherwise this object is not is not defined. So it is a function that is defined on the complex plane except the set of corresponding to the eigenvalues and this is defined for any fixed matrix. Now what we can do if X is a random random matrix what we can do is we can take the average of this object it becomes because it becomes a random function. So we can take the average over the ensemble and then we can take the limit and to infinity. And this object here we define a mathematical G of Z and this object is called resolvent. It has several names it can be called greens function it has several names but these two are are enough for our purposes. So we have a random function random complex function we take the average over the ensemble and then we take the limit and to infinity. So the claim is that this object is equal to the integral of the spectral density in the limit and to infinity times Z minus lambda. So I will show why this claim makes sense but this formula establishes a connection between the greens function or the resolvent and the average spectral density in the limit and to infinity. So if we are able to compute this object for our ensemble and if we are able to invert this relation which is an integral relation. So it is not obvious that we are able to invert it then we have access to this object. So it turns out that in many problems it is easier to compute this object than to compute directly the spectral density. So we use this indirect method to access the spectral density. Of course we have two problems so first to prove that this is true. And second how do I invert this relation to find the spectral density from the resolvent. You see this direction is easy. If I have the spectral density I can always compute an integral and find the resolvent. But this direction is hard because if I know this function how do I reconstruct the integrand. But it can be done. Absolutely yes. There is a very specific connection. This I think we will discuss a bit in two lectures. Because we will see this in a completely different electrostatic context. So this will become clear why the name is like if you can wait just a couple of lectures this will become clear. So let me give you a proof of this. Okay so first we start from a finite n version of this of the right hand side. So we take the finite n spectral density we divide by z minus lambda and we integrate over lambda. Then we inject here the definition of the finite n spectral density. Definition is integral d lambda divided by z minus lambda and then we have the average of 1 over n summation over i delta of lambda minus lambda i. I am not doing anything I am just inserting here the definition of the spectral density for finite n. Now the average and the integral would commute because that's a linear operation. So you have that this is equal to 1 over n summation over i integral d lambda delta of lambda minus lambda i divided by z minus lambda. This delta would kill this integral and reconstruct a z minus lambda i. So this object here is just the average of g n of z and then you take the limit n to infinity to the left and to the right. So this is a two minutes proof that mathematicians would not like it but it works that it is at least plausible that this object is related to this object in precisely this way. Are you convinced? Who is utterly not convinced? So we have established this relation. Can I remove? So properties of the Green's function or resolvent. So we have that g of z is equal to the integral of rho of lambda over z minus lambda. So the first property if you take z in absolute value very large what you have here is that downstairs this object will dominate over this. So this term will dominate over lambda so you can pull out a factor 1 over z and then we have this object here is equal to 1 because the density is normalized. So we have that the Green's function goes for large absolute value of z as 1 over z. This is sometimes useful because this property is basically equivalent to the property of normalization of rho. So when you derive your function g of z the first thing to check is that it goes as 1 over z for large models of z. So you need to impose that this coefficient is equal to 1. For example you might have parametric case where you have here a certain function of parameters and you need to impose that this function is equal to 1. Otherwise your original spectral density is not normalized. In the Gaussian case. Excellent question. So the right way to do it is to consider the spectral density where you have put your where you have rescaled your eigenvalues in such a way that they are of order 1. So for example in the Gaussian case you would rescale the eigenvalues by square root of n and your eigenvalues lambda i tilde would be of order of order 1 for n going to infinity. So the edges would not scale with n and that's the definition that I implicitly assumed for our spectral density. But that's a good point. Thanks. Good. So the second property it is the generating function of trace moments. What does this mean? You can interpret this object as a geometric series. So you take integral d lambda rho of lambda of z minus lambda. So you can pull out a factor of z downstairs. So this is 1 over z integral d lambda rho of lambda over 1 minus lambda over z. And then you can interpret this object here as the result of a geometric series. 1 over 1 minus x is the result of a geometric series. So you can write 1 over z d lambda rho of lambda and you have a summation k from 0 to infinity k I call it k lambda over z to the power k. So if you go like this, this would be summation k from 0 to infinity of mu k over z k plus 1. Because you have z to the k and there is an extra z here. And this mu k are the integral d lambda rho of lambda lambda to the k. So this is the limit n to infinity of 1 over n trace of x to the k. So the resolvent encodes all the moments of the spectral density. The integral of spectral density times lambda to the k which can be trivially mapped into this object, the average of the trace of x to the power k. So all of these are encoded in g of z. If you develop, if you make a Taylor expansion, a Loran expansion sorry in z. So sometimes this property is also useful. But the third property is the most useful of them all. Can I raise it? Yes? So the third property is how to reconstruct rho of lambda from g of z. And for this we need a very important formula which is the Sochotsky-Clemay. So this guy was Polish, half Polish, half Russian. Meaning it was Polish when Poland was under Russia. And this guy is Slovenian. So he came from not far away from here. So this formula is like this. So I first write it and then explain it. So if you have a rational function with a complex part and you take the limit of the imaginary part to zero, then the imaginary part of this object reconstructs a delta function. Which means that if you need a delta function in any of your calculation, you can always replace it with the imaginary part of a rational function. The proof, suppose that you have phi which is a complex function phi, which is a complex valued function defined and continuous on the real line. And then you take two numbers a less than zero less than b. So what you compute is integral between a and b of your function phi divided by this object here, y plus i epsilon. So you just take a test function, you divide it by this left hand side, y plus i epsilon, and you compute the integral over an integral of the real line containing zero. Now to perform this integral, you multiply up and down the integral by y minus i epsilon. So I haven't done anything, just multiplying up and down by y minus i epsilon. So the denominator becomes y square plus epsilon square, and the numerator is y minus i epsilon. And now I separate the real and imaginary part of this object. So you have y over y square plus epsilon square minus i integral psi of y epsilon divided y square plus epsilon square. And now the theorem is basically established because this object here in the limit epsilon going to zero reconstruct a one over y term. So here you have an integral of psi of y over y, and this integral will be singular in general because you have a divergence of the type one over y at zero. So this integral in the limit epsilon to zero is only defined if you take the Cauchy principle value. So the real part of this object reconstruct the Cauchy principle value. The imaginary part reconstruct a delta function because you have a representation of the delta function of this type. So delta of y delta function can be written as limit epsilon to zero plus of one over pi epsilon divided y square plus epsilon square. So this object here tends to pi delta of y. So this is the Cauchy or Lorentzian representation of the delta function. So this formula is very useful because it gives a representation of the Dirac delta of y in terms of a rational function. So we will use it a lot. Why is it useful for us for reconstructing the spectral density from the resolvent? Well, this is easy. So suppose that I compute my resolvent function not in z but in x minus i epsilon. Yeah, this bit here. So you have a here, b here and zero is in between. So when you squeeze the denominator by sending epsilon to zero, so this object here becomes a y over y square. So it becomes a one over y. So in general for a well-behaved function, test function psi of y, if you divide it by y, the integrand will acquire a singularity of the type one over y at zero, which is non-integrable. So the only way to make sense of this integral is by taking the Cauchy principle value, which means that basically you are interpreting this object as the integral up to let's say zero minus epsilon plus the integral from zero plus epsilon onwards and then you take the limit from epsilon going to zero, which means that the two divergences on the two sides of zero would cancel out. So if you now compute the resolvent not in z but in x minus i epsilon, so here you have the point x on the real axis and you compute it here. Just you move off the real line a bit. Well, this is by definition integral d lambda rho of lambda over x minus i epsilon minus lambda. This is just by definition. And then suppose that you take the imaginary part of this object, the imaginary part of g of x minus i epsilon, basically here you can play the same trick that we are playing here. So this imaginary part will be this object here. So this will be the integral over d lambda rho of lambda pi times delta of x minus lambda. So the delta function will kill this integral and this will become pi rho of x. If you send of course epsilon to zero plus. So eventually the spectral density at the point x will be given by one over pi, the limit epsilon to zero plus of the imaginary part of the greens function evaluated at x minus a bit of offset on the complex plane. So this is the inversion formula that we were after. We can compute the spectral density at x if we know the resolvent in the complex plane. So it is not enough to know the resolvent on the real numbers. We need to use the complex structure of the argument of g. But by using this formula we can reconstruct the spectral density if we know g of z in the vicinity of the point x. So this will be our main tool because if we manage to find an equation for example for g of z, we will solve this equation, find g of z, compute it in x minus i epsilon, take the imaginary part, take the limit epsilon to zero in this order, divide by one over pi and we will have our spectral density. This is what we are going to do for the Gaussian ensemble. We can do a five minutes break now. So what I wanted to do now in this half an hour, let's see if I can, is to try to use this formalism to compute the spectral density for Gaussian ensembles. In principle you don't need to do this because we have much better tools for Gaussian ensembles. But this is just for training purposes, just to show you how we can put this method to work in practice. So g of z for Gaussian ensembles. So the idea would be to establish an equation for g of z and then from here to invert this expression using Sokowski's Plemage formula and obtain the spectral density in the limit n to infinity. So we can start from the joint pdf of the eigenvalues of Gaussian matrices that I gave you before without proof. So remember we have the Gaussian weight and then this interaction all to all term and zn of beta is the normalization factor which makes sure that this is a proper joint pdf. So it should be normalized to one. So zn of beta will be then the integral of d'1 d'n exponential of minus one of summation i1 to n'i2 and then you have this product to the power beta. Okay fine so now in this integral we make a change of variable, change of variables. So I send lambda i into root beta n lambda i. So if I do this change of variable what I get is that our normalization constant is equal to a certain number that we don't care much about. This number is basically just the product of these terms and some terms coming from here. Okay times and here comes the interesting bit. I can write this object in this form the integrand. So what I'm doing here I'm basically raising this term into the exponential. So I'm rewriting this object as exponential of beta over 2 summation j different from k log of lambda j minus lambda k. So I'm only rewriting this product by writing it as in the exponential. So why am I doing this? Well I'm doing this because if you do this trick here and after the rescaling of the eigenvalues then you can rewrite this integrand in this nice exponential form. Apart from this n exponential of minus beta h. Have you ever seen exponential of minus beta h before? Where? Sorry? So let's write what is h. So h is one half summation i1 to n lambda i square minus 1 over 2n summation j different from k log of lambda, sorry, i different from j log lambda i minus lambda j. So look at this very nice formula. This formula here explains why we chose the letter beta for the Dyson index. The letter beta we could have called it gamma or delta but if we choose beta then you recognize here exponential of minus beta h. So the Gibbs Boltzmann weight of an associated thermodynamical system. And here we are integrating over the position, the degrees of freedom which are the position of your particles. So this object which originally was just a normalization constant of a joint probability density can now be interpreted as what? As the canonical partition function of an associated thermodynamical system. So one of the great advancement in the field was when people realized that the eigenvalues of random matrices, in this case Gaussian random matrices, but the concept is much more general can be viewed as an associated thermodynamical system. So we can treat eigenvalues as particles as interacting as a gas of interacting particles. And what is the corresponding Hamiltonian? Hamiltonian contains two pieces that are in competition. There is a quadratic part and the null tool interaction part. So we have our eigenvalues behave as a gas of charged particles that are basically confined inside a quadratic well so they cannot escape to infinity. But they are also repelling each other like this object is repelling this one, it is repelling this one, it is repelling this one. So every particle feels the presence of all the others. So this is a classical long range thermodynamical system. But it is a long range system so it is not the standard type of systems that we see in a standard course in statistical mechanics. Because every particle feels the presence of all the others. Now then thermodynamics tells us what we should do. So this type of gases is called two-dimensional Hulomb gas. Can someone explain to me why it is called two-dimensional? The eigenvalues are online, right? So why is this a two-dimensional system? Yeah, so the interaction between particles is logarithmic. So Hulomb interaction is logarithmic in two dimensions, not in one, not in three. So we need to think about this as a system of charged particles which are basically leaving on a plane. But are forced to stay on a line. And the physics of the system is entirely different from assuming that they are interacting with a 1D Hulomb potential, which is linear. So we have a two-dimensional system of charged particles confined to a line. And this is different from having a one-dimensional system. Is it clear? Okay. What is the one-dimensional Hulomb potential? How does it go with the distance? How does it scale with the distance? One-dimensional? Sorry? So how does the potential... Sorry? One over? In how many dimensions? In three. In three-dimensional it scales as one over R. In one dimension? R? No. So the electromagnetic potential scales as R in one dimension. It scales as log R in two dimensions. And it scales as inverse power of R in higher dimensions. Okay? So we need to have a two-dimensional potential to reconstruct this object. But your eigenvalues are real so you are forced to constrain your particles to leave on a line. So they leave on a line but they interact with the wrong potential. With the potential corresponding to another dimension. So now with this intuition we can ask our question. So question for let's say a very large number of particles. What is the most likely configuration of your particles at equilibrium? So now we can again forget completely about random matrices. We have a problem in classical teradynamics. We have a system of charged particles in two-dimensions constrained on a line. Which are subject to two competing interactions. A quadratic confinement well and a null-to-all repulsion term. So what is the most likely configurations of these charges at equilibrium? How do we find it? What is the standard rule? This is the canonical partition function. So what shall we do? Sorry? The derivative of what? You said take the derivative, right? Well, yeah, but it's occupation probabilities. You will need to have like a discretization of your real axis. Otherwise you cannot say that a certain box contains a certain number of particles. Yeah? Yeah. Yeah. Excellent. So you see that here we have that if n is very big, then the inverse temperature doesn't really matter. So whatever temperature you start with, the limit n to infinity will drive your system to the zero temperature limit. Indeed at zero temperature limits you don't have entropy. So your particles don't fluctuate at all. So the configuration that they will arrange themselves into is the configuration of lowest energy. So what is the configuration of your particles that minimize the total energy of your system? Well, what we have to do is to take the derivative of h with respect to lambda j and take this equal to zero for all particles. Or equivalently we need to force every particle to be in thermodynamic or mechanical equilibrium. So the force acting on a single particle must be balanced and must be equal to zero. Okay? This is the precise definition of equilibrium. Otherwise the particles would keep moving. Just pretend and nod. They would like to take a picture of you. You need a coffee? Smile, just do something. Okay. So I assume that it is clear that what we have to do is to take the derivative of our energy function with respect to the position of each particle and equate it to zero for all j from one to n. So if we do that we need to take the derivative of this object with respect to lambda j and take it equal to zero. Or actually I put it lambda i just the same. So if you do that then you get to this equation. So the most probable position of the eigenvalue or the particle i should satisfy this stability condition, j different from i, one over lambda i minus lambda j. So this is basically a mechanical stability condition for your particles, which is obtained by imposing that the total energy of the system is stationary. So the system will occupy the minimum amount, it will use up the minimum amount of energy that it can, if the particles are in this, satisfy this equation. So they will not move from there because any deviation from this configuration will cost energy. Picture? Good. Now the idea is to use this equation which is valid for let's say large but finite n to find an equation for the resolvent. We find an equation for the resolvent, we solve this equation and we invert it and then we get the access to the spectral density for the Gaussian ensemble. And then we compare this result with numerical diagonalization of matrices and we see that we have done the right thing. So computing the equation is quite, well it requires just a couple of manipulations of this stability condition equation. So for example in the left hand side what we do is we multiply by 1 over n z minus lambda i and then sum over i. So here on the left hand side we have lambda i divided by n z minus lambda i sum i 1 to n. Now I will use the green chalk. So I will now add and subtract z. So what you get here is lambda i minus z which cancel with z minus lambda i and produce a factor of minus 1. So this object is minus 1. And then what is left is a factor of z times the definition of your let's say finite n resolvent. So you see that simply by algebraic manipulations on the stability condition so the fact that your particles need to occupy the lowest energy level we manage to get something in here which resembles which is linked to the resolvent. This is just by purely algebraic manipulations. Now we do the same thing on the right hand side and in the end we will take the limit n to infinity and we will get basically an equation for the resolvent that we can then solve, try to solve. Now if we do the same thing on the right hand side so let's call r the right hand side of this equation. So we need to multiply by 1 over n z minus lambda i and sum over i. So there is a 1 over n here and the 1 over n here. So in total we have 1 over n square summation over i summation over j different from i. Then we have 1 over z minus lambda i 1 over lambda i minus lambda j. Now it is quite easy to show that you can basically split these two fractions here and rewrite this as 1 over z minus lambda j which multiplies 1 over z minus lambda i plus 1 over lambda i minus lambda j. Why? Because if you do this sum you get z minus lambda i lambda i minus lambda j and on top you get lambda i minus lambda j plus z minus lambda i which is then cancelled by this guy here. So I am just rewriting this product by splitting the two fractions. And then what we have basically is the sum of two terms a plus b. So we now analyze the two terms separately. So a is equal to 1 over n square summation over i 1 to n summation over j different from i. And we have 1 over z minus lambda j 1 over z minus lambda i. It is the product of these two guys here. So what we can do here is we add and subtract the term j equal to i. So we have that this is equal to 1 over n square summation over i summation over j. 1 over z minus lambda i lambda j 1 over z minus lambda i. And then we need to subtract the term where i is equal to j. So it will be summation i 1 to n 1 over z minus lambda i square. I am just adding and subtracting the term i equal to j. Just nod. So this object here we can split the two sums so summation over i times summation over j. So this object here together with 1 over n square is just g n square of z. I recall that g n of z was defined as 1 over n summation over i 1 over z minus lambda i. So here we have reconstructed the square of our would be resolvent. It is the finite n version of it. And here what we have is this object here is plus 1 over n the derivative of g n. So if you differentiate this object with respect to z you get a square in the game. But there is also minus sign that will cancel out this object. So our first term a is equal to g n square plus 1 over n g n prime. Now we work out what b is. Then we sum them up and equate the result to the left hand side. So b is equal to what? Is equal to 1 over n square summation over i 1 to n summation over j different from i. 1 over z minus lambda j 1 over lambda i minus lambda j. Now I will ask you to convince yourself that this object is a symmetric function in lambda 1 lambda n. Well it's best if I actually so suppose that you take n equal to 2. And you try to work out what this guy is. So you just write it explicitly. So this would be 1 over z minus lambda 2 1 over lambda 1 minus lambda 2 plus 1 over z minus lambda 1 times 1 over lambda 2 minus lambda 1. Just this object here for n equals to 2. And now you do the same experiment with this guy here. So we take what r was. So r is 1 over z minus lambda 1 instead of z minus lambda 2. 1 over lambda 1 minus lambda 2 plus 1 over z minus lambda 2 times 1 over lambda 2 minus lambda 1. So if you compare this object and this object you see that b is exactly equal to r with a sign difference. So there is a minus sign between the two. And this is just an example for n equals to 2 but it remains true for any n. So what you have is that b is equal to minus r. So now you have an equation for r. So you have that r is equal to a so g n square z plus 1 over n g n prime z minus so plus b which is minus r. Which means that r the right hand side is 1 over 2 g n square z plus 1 over 2 n g n prime of z. Now all you have to do is to equate the left hand side which was this one to the right hand side and solve the equation for g. Can I raise this? So the equation that we want to solve would be in principle plus 1 over 2 n g n prime of z. So in principle we could try to solve this differential equation. So let's classify this differential equation. It has a name. Excellent. So this is a Riccati equation from the name of the Venetian Nobel band who first studied it. It's first derivative equal to a second degree polynomial. This is a Riccati equation. Fortunately we don't even need to bother because in the limit n to infinity since g of n is of order 1 this object will be sub leading in n. Because it is divided by 1 over n. So what we need to do is we discard the first derivative and our differential equation becomes just an algebraic equation for the result. So what you have to solve is 1 of g n square of z minus z g n of z plus 1 equal to 0. And since we have taken the limit n to infinity we just cheat a bit and we replace this object with our. Proper result. Good. So you have a second degree equation which is solved by this object. So now we know what we have to do. So we have to compute g at x minus i epsilon and then take the imaginary part of it. And then divide by 1 over pi and then take the limit epsilon to 0. Quite a lot of things to do. I will do it in 2 minutes so that we can start afresh tomorrow. And also you get your first spectral density computed. So g of x minus i epsilon is equal to x minus i epsilon plus minus square root of x square minus epsilon square minus 2. Plus i minus 2x epsilon. I just put x minus epsilon here raised it to the power 2 and then separate real and imaginary part of the radicand. Okay. Now lemma. So the principal square root of the complex number a plus ib is p plus iq. So this lemma gives you so if you have a square root of a complex number we can write this principal square root as a complex number so p plus iq. And there is a formula that connects p and q to a and b. Okay. Maybe you know it but it's not so easy to find. So this is the real part and this is the imaginary part of a principal square root in Cartesian coordinates. So I let you work this out for this square root. So what you get is that in the limit epsilon to 0 plus the result is equal to plus minus sign of minus x over pi root 2 square root of x 2 minus 2 minus x 2 plus 2. So you see you have to take another square root in here which will become a perfect square in the limit epsilon to 0. So you have an absolute value in here. And then you have your result that if absolute value of x is larger than root 2 the density is 0. And if the absolute value of x is smaller than root 2 the density rho of x is equal to 1 over pi root 2 minus x square. Have you ever seen this low before? No? Okay. So it has a name. So it is named after Wigner again. So it's called Wigner semicircle. Why is it called semicircle? Why is it called a semicircle? Yeah. Except that it is not true, right? So this is not an equation of the semicircle. It's a semi ellipse but for some reasons it's called semicircle. Okay. So your density, your energy, your average spectral density assume this semi elliptical low for Gaussian ensembles for any beta in the large, in the larger limit. And actually you can check this on page 7. You find a comparison with numerical diagonalization of Gaussian matrices. So if you remember at the beginning we rescaled everything by root beta n. So if you rescaled your eigenvalues that you have sampled in your numerical code by root beta n then the curves collapse on this semi ellipse very nicely. Okay. Which means that on average you will have positive, half positive and half negative eigenvalues. And most likely the eigenvalues will be around zero. So there is a higher probability of getting eigenvalues around zero. Okay. Just try to, there is the code here. You can try in MATLAB to reproduce this nice picture. See you tomorrow.