The Numerical Renormalization Group

Introduction

Finally here¹ it is. A simple (relatively) program implementing the Numerical Renormalization Group. I tried to implement it as simple as possible to be easy to understand. Here is the program in action, running an Anderson model at half-filling:

I already used charts generated by the program in The Kondo Effect and Renormalization Groups posts. The first one shows the spectral function where you can notice the broadening of the quantum dot energy levels because of interaction with the leads and also the Kondo resonance, the next one shows the renormalization group flow for the same model.

Links

I already gave some links about the subject in the The Kondo Effect, Quantum Dots and Renormalization Groups posts, here I’ll give some more, including with some information on extensions of NRG.

To theory

Here are two lectures of Theo Costi: Numerical Renormalization Group for Quantum Impurities² and Numerical renormalization group and multi-orbital Kondo physics³.

A paper by Oliveira: The Numerical Renormalization Group and the Problem of Impurities in Metals⁴.

A paper on DM-NRG (Density Matrix – NRG) by Hofstetter: Generalized Numerical Renormalization Group for Dynamical Quantities⁵.

A paper on TD-NRG by Frithjof B. Anders and Avraham Schiller: Spin Precession and Real Time Dynamics in the Kondo Model: A Time-Dependent Numerical Renormalization-Group Study⁶.

A paper on using non-Abelian symmetries to improve DM-NRG by A. I. Toth, C. P. Moca, O. Legeza and G. Zarand: Density matrix numerical renormalization group for non-Abelian symmetries⁷.

Sindel Michael PhD thesis: Numerical Renormalization Group studies of Quantum Impurity Models in the Strong Coupling Limit⁸.

To code

Again, here¹ is the link to the code presented in this post. It is only a toy program, just to illustrate the concepts.

Programs used in research are here: Flexible DN-NRG⁹ and here: NRG Ljubljana¹⁰. I don’t know much about the later, but by looking at the code I can tell that both are capable of using symmetries to speed up the calculations and improve accuracy, they work for non-flat density of states, NRG can be used together with Dynamical mean-field theory, they allow multiple channels, the spectral function can be calculated much better with DM-NRG (at finite temperature, too), they implement z-averaging, they can calculate out from the spectral function the Green function (with Hilbert transform) and more.

Some brief theory

Obviously it would be quite hard to present here the whole theory, I’ll try to sketch the main ideas only. For details, please check out the links.
From now on, both in the theory presented and in the code a flat density of states is assumed in the conduction band and also a constant coupling between the impurity and the electronic bath (that is, not dependent on energy).

Derivation

Approximately those are the steps of the derivation, please see the links for details:

1. Logarithmic division of the conduction band. This was the first step in the derivation and the most important one, it will lead to killing the logarithmic divergence mentioned in The Kondo Effect post. The conduction band is split into intervals that decrease towards low energies: .
2. Fourier expansion of each interval, using orthogonal basis functions that are zero outside the interval.
3. Approximating by using only the Wannier orbital centered on the impurity. It turns out that the impurity couples directly only to the fundamental harmonic, the other terms couple only indirectly through the first term and decouple in the limit .
4. Lanczos tridiagonalization: Wilson chain. The previous step ended up with the impurity coupled to a p=0 term for each interval, in order to be able to solve it iteratively, it is tridiagonalized using the Lanczos algorithm. The result is the impurity coupled to a single site which in turn couples to the next one and so on, the couplings decreasing exponentially. This ensures a energy scale separation at each step.
5. The renormalization group transform. A recursion relation is obtained: that allows to define a renormalization group transform .
6. Iteration with truncation of the eigenvalues and eigenvectors. Because the Hilbert space grows exponentially one cannot keep up all states and energies while new sites are added, a truncation scheme is applied. The high energy states that are discarded do not affect the low energy states that are added. Each step corresponds to an energy scale and also to a temperature.

The Algorithm

In the following, the states for each site are in order:
Here is a brief description of the algorithm together with relevant code:

• Start with the Hamiltonian for the ‘impurity’ site (for example a quantum dot or even several of them coupled), perhaps together with the first site of the Wilson chain. Diagonalize it and transform accordingly the operators that are carried along. In the code, startIteration = -1 means the ‘impurity’ alone, startIteration = 0 means the ‘impurity’ together with the first site from the Wilson chain. Here is the code for the Anderson model quantum dot:

B is the magnetic field, U the Coulomb interaction energy, eps is the uni-electronic energy level (for electron-hole symmetry, -U/2). delta gives the hybridization. The up operator is the spectral operator for the spectral function.

• Continue with the loop that adds up one site from the Wilson chain, iteratively. Here is the code for the loop together with the diagonalization of the Hamiltonian before the first iteration:

Each ‘step’ several things are done, they are implemented in void NRGAlgorithm::Step(int iter).

First, the Hilbert space is enlarged with the states of the newly added site, then the Hamiltonian is set by adding the for the new site.

t is easy for a flat conduction band and constant coupling:

If that’s not the case (for example when using Dynamical mean-field theory the band is adjusted self-consistently and it’s not constant), you will have to do the Lanczos tridiagonalization in the code, it’s not that simple. You’ll also end up having on site energy, so setting up the extended Hamiltonian will not be as simple as here.

Here is the code that sets up the extended Hamiltonian in one step:

Don’t forget about the anti-commutation relations. A minus sign appears for states with one electron from the added site because of that.

Then the Hamiltonian is diagonalized and truncated:

Then the unitary transformation is applied to all operators, including the f operators needed for the addition of the next site. We need them all in the same basis as the Hamiltonian.

Then prepare for the next iteration and that’s about it:

How to extract information

I deliberately not calculated all easy things that can be computed, just to avoid unnecessary complications, to keep the program simple. Again, for more information check the links. Even the spectral function could be calculated faster using Fourier Transform (FFT implementation) and convolution but I guess the code is more clear as it is.

The easiest thing to calculate is the partition function. Just stop the algoritm at a step corresponding to a particular temperature and calculate it, the spectrum is available:

From here you can already extract quite a bit of information.

By noticing that not only the spectrum is available, but also the eigenvectors, one can compute the density matrix:

At this point I must emphasize again that this is a toy program and accessing this⁵ is recommended.

Since the program can bring along the ‘static’ operators, extending them, changing the basis and truncating them, one can do even more by using an operator for an observable, O:

All the above are quite easy to calculate (not exactly if you want to use the complete set of eigenvectors and go with the reduced density matrix, for that see the link provided).

I did not calculate any of them, I chose to calculate the spectral function instead, for the limit T=0. Very shortly, the spectral function for the d operator in the T=0 limit is:

There is enough information about how to derive it and implement it in the links, the code is in the SpectralOperator class. The truncated spectrum is accumulated each even step – in SpectralOperator::PassSpectral – by using a weight for the overlapping interval, to avoid double counting. Before getting the spectral function into the chart, the discrete values are broadened in SpectralOperator::GetSpectrum() using a log Gaussian:

A spectral function calculated this way does not respect the spectral sum rule, one can do much better with DM-NRG⁵.

Having the spectral function allows one to calculate all sorts of things, like conductivity or spin susceptibility. If spectral function is not enough, the whole Green function can be computed using a Hilbert transform.

How to improve it

I already gave some hints about more that could be implemented, above, but here are some, again:

• The best thing to do is to use symmetries to speed up computation and improve precision. If you find operators (corresponding to a symmetry) that commute with the Hamiltonian (that is, the observable values are conserved), it means that you can find a common eigenbasis. The eigenvalues of such operator are good quantum numbers and you can label states with them. Obviously you may find more than one such operator. It’s easier if they are abelian, if they are not you’ll have to deal with Clebsch-Gordan coefficients, please see this for details⁷.
  Anyway, you can sort states by the quantum numbers to have the the Hamiltonian matrix block diagonal, because subspaces for different quantum numbers decouple. The Hamiltonian acting on such a state leaves it unchanged (apart from the multiplication with the eigenvalue). Also operators will be composed of non-zero blocks (sub-matrices) and the rest of the big matrix will be zero.
  For example, you may find that the number of electrons/the charge conserves. That’s the U(1) symmetry (the same goes for Sz). If you look at the initial Hamiltonian for the Anderson model, you’ll see that you have two 1×1 matrices for the 0 and 2 electrons and a 2×2 matrix for the 1 electron states. Having to diagonalize a matrix that is N/2 x N/2 is much faster because diagonalization is very expensive, being an order of .
  For spin conservation, the symmetry is SU(2) and if you use more then one symmetry you can split the Hamiltonian in even smaller sectors. I would use charge and Sz first, that is U(1)xU(1). Again, in the beginning you might want to avoid non-abelian symmetries. Here is a Wikipedia page that might be also informative regarding commuting operators: Complete set of commuting observables.
• One could make the program more general, allowing for a more general density of states for the conduction band and frequency dependent couplings of the impurity. I already mentioned the possibility of coupling an NRG program to a dynamical mean field theory one as an impurity solver.
• Make it multi-channel.
• Now both the Hamiltonian for the ‘impurity’ and the hopping operators are hardwired, a serious program allows specifying them (along with static and dynamic operators) in configuration, to be easily changed.
• Implement z-averaging. I also think that adaptive discretization schemes were considered.
• Implement DM-NRG⁶.
• And more… please check out the links for more details, also see above for some of the many things that can be calculated.

The program

Here is a short description of the program, it will be also available on GitHub¹.

Description

Everything related with NRG is either in the NRG namespace or has the class name starting with NRG. There are three kind of classes for NRG implementation, one is dealing with data passing around and adjusting and controlling the algorithm running, one is the operators, derived from the Operator abstract class and one is the NRG algorithms, derived from the abstract class NRGAlgorithm.

The rest of the program is very simple, just an interface to the NRG. It allows starting/stopping the computation, some configuration settings and it displays the charts. That’s about it.

Libraries

Besides mfc and other typical VC++ runtime libraries, the program uses GDI+ for drawing.
The program deals with matrices using Eigen¹¹ library.

Classes

The NRG Namespace:

ControllerInterface and ResultsRetrieverInterface are interfaces that allow by deriving from them classes that respectively cancel calculation and get the results from it.

The operators are derived from the Operator class. Operator::Extend() extends the operator by adding new states for the new Wilson site. Added states are in the most significant bits position. The changeSign member allows extending the operator matrix for fermionic operator type (if true) or bosonic operator type (if false). The minus sign there is due of anti-commutation. Classes derived from it are: Hamiltonian, the hopping operators FUpOperator and FDownOperator and the spectral operator, SpectralOperator. This one is a regular operator with some methods added that allow calculating the spectral function for the operator. DUpOperator is the spectral operator that is used for generating the spectral function for the Anderson and two quantum dots models.

The NRGAlgorithm class implements the NRG. From this class three examples are derived: QDAnderson, a quantum dot with the Anderson model, QDKondo, a quantum dot with the Kondo model, TwoQDRKKY, two quantum dots coupled by spin-spin interaction, only one being coupled to the leads. The later should be considered only qualitatively, to have better precision one should use symmetries for calculation. Anyway, it’s enough to show the split of the Kondo resonance due of the two stage Kondo effect.

NRGComputationThread is the class that implements the computation thread for NRG, the calculations run in a different thread to avoid locking the UI.

NRGController is derived from NRG::ControllerInterface and allows cancelling computations (the thread checks it each computation step).

NRGResultsData is derived from NRG::ResultsRetrieverInterface and allows passing the results to the main thread and allows it to check if computation is finished.

The options are implemented by Options and they are saved/loaded into/from registry. The options UI are implemented by COptionsPropertySheet, CNRGPropertyPage, CParametersPropertyPage and CChartsPropertyPage.

The charts are implemented by the Chart class. It’s pretty messy and far from perfect, I might improve it in the future. It uses GDI+ for drawing.

CAboutBox needs no explanation.

CMainFrame is the main frame window, it implements/routes commands.

CnrgApp is the application class. There aren’t many changes in there except initializing and shutting down GDI+, setting the registry key and loading the options from registry.

CnrgDoc is the ‘document’ class. It contains the computation thread, the thread controller and the computation data objects. It also contains the chart objects. The most important member is CnrgDoc::StartComputation() the others are pretty straightforward.

CnrgView is the ‘view’ class. Has some changes compared with the class generated by App Wizard, related with drawing/printing. There is a timer implemented there which allows checking for computation finish and updating the charts. There is also some handling of the cursor, making it a ‘wait’ cursor during calculations.

CNumberEdit implements an edit box for double and float values. By setting allowNegative one can control if negative numbers can be entered or not.

ComputationThread is the base class for the NRG thread. There is not much in there, just starting the thread.

Some results

Besides the simple Anderson and Kondo models I already mentioned and even supplied some results from the code already, here are some more results for a little bit more complex situations.

Two quantum dots interacting by spin-spin interaction

As a little bit more complex example I implemented the class TwoQDRKKY (RKKY comes from Ruderman-Kittel-Kasuya-Yoshida interaction). It has a Hamiltonian set up for two quantum dots coupled by antiferromagnetic spin-spin interaction, only one being connected to the leads. Here is the code that sets the Hamiltonian, the comments in the code should be enough for understanding it:

Here is a spectral function I got:

And here are two links to papers dealing with such setup: Strongly correlated regimes in a double quantum-dot device¹² and Two-stage Kondo effect in side-coupled quantum dots: Renormalized perturbative scaling theory and Numerical Renormalization Group analysis¹³.

Notice how the Kondo resonance is split (for a larger coupling between the quantum dots it is completely destroyed). As temperature goes lower, at Kondo temperature the electron from the quantum dot is screened by the electrons in the leads, the electron being ‘locked’ in a singlet state. This prevents inelastic scattering. At higher temperatures this cannot happen because the thermal energy is high enough to overcome the coupling strength. This is the usual Kondo effect already presented.

When the first quantum dot is screened by the electrons in the leads, they present together for the second quantum dot an effective fermionic bath. As it is coupled also by an antiferromagnetic interaction, as the temperature is lowered towards an energy equal with the coupling strength, the two quantum dots will lock into a singlet. Electrons from the leads need to break the coupling between the two quantum dots in order to pass through and since they do not have the thermal energy to do it anymore, the conductance drops. This is the second stage Kondo effect. Please check the links for the details.

As a warning, this setup is quite complex and it would need to use symmetries to avoid numerical errors. Without symmetries one would need to use a big matrix and numerical errors start to kick in due of diagonalization. As a consequence it’s quite hard to get a symmetrical spectral function. Here are the parameters I used for the chart: number of kept states 250 (this is very low!), lambda 2.5, 30 iterations, U=2.82, epsilon=-U/2, delta=U/16, J=0.02, no magnetic field. For spectral function, b=0.6, the step was 0.001.

Magnetic field for Anderson model

Now let’s go back to the Anderson model and apply a magnetic field. A weak magnetic field, for illustration purposes I want the Zeeman energy being quite a bit lower than the Kondo temperature.
Here is the spectral function for the operator:

It’s easy to guess how the ‘down’ spectral function looks like. Here is the spectral function for both spin up and down operators:

From this one is easy to figure out the conductance, as in the previous case the Kondo peak splits, when the temperature is low enough electrons from the leads will not be able to break the coupling with the magnetic field, the conductance drops.

For this spectral function I had to temporarily change the QDAnderson::Init() implementation, here is the relevant code:

The usage of FDownOperator might be confusing, but DUpOperator has the same implementation as FUpOperator except that it has the code for calculating the spectrum. The same would go for a DDownOperator, I did not bother to implement it since I could use FDownOperator.

For more information on this please see this Nature article (picked at random): Temperature and magnetic field dependence of a Kondo system in the weak coupling regime¹⁴.

The End

With this post I end for now the set of posts on this topic, although I intend to implement at least another renormalization group program, a DMRG one. But until then I’ll probably have some other posts on other topics, possibly on the Monte Carlo methods.

As usual, if you notice any mistakes/issues, please let me know, I’ll try to fix them.

CodeProject