Describing systems of interacting fermions by boson models:exact mapping in arbitrary dimension and applications

Describing systems of interacting fermions by boson models:
exact mapping in arbitrary dimension and applications

K. B. Efetov    C. Pépin    H. Meier Theoretische Physik III, Ruhr-Universität Bochum, 44780 Bochum, Germany
IPhT, CEA-Saclay, L’Orme des Merisiers, 91191 Gif-sur-Yvette, France
July 5, 2019

We develop a new method that allows us to map models of interacting fermions onto bosonic models describing collective excitations in an arbitrary dimension. This mapping becomes exact in the thermodynamic limit in the presence of a bath. The boson models can be written either in the form of a model of non-interacting bosons in a fluctuating auxiliary field or in the form of a superfield theory of interacting bosons. We show how one can study the latter version using perturbation theory. Using the developed diagrammatic technique we compared the first two orders of perturbation theory with the corresponding results for the original fermion model and found a perfect agreement. As concerns the former representation, we suggest a scheme that may be suitable for Monte Carlo simulations and demonstrate that it is free of the fermionic sign problem. We discuss in details the properties of the bosonic representation and argue that there should not be any obstacles preventing from an efficient computation.

71.10.Ay, 71.10.Pm, 75.40.Cx

I Introduction

Many body electron systems with interaction traditionally attract a strong attention. Numerous interesting phenomena of the condensed matter physics originate from the electron-electron, electron-phonon and other types of the interaction and do not exist in the ideal Fermi gas of the electrons.

Generally, models with interactions cannot be solved exactly in the space of dimensionality and one has to use approximate schemes of calculation. The most developed method is the perturbation theory with respect to the interaction. Within this theory one starts with the ideal Fermi gas and calculates corrections to physical quantities assuming that the interaction is weak. Sometimes it is sufficient to compute only the first several orders of this perturbation theory, but very often in order to capture important physics, one has to sum certain series.

The formalism of the perturbation theory for the interacting electron gas is well developed (see, e.g., Ref. agd, ). Each term of the perturbation theory is represented diagrammatically, which enables one to carry out quite complicated calculations. Many interesting phenomena like, e.g., superconductivity are not seen in the lowest orders of the diagrammatic expansion and show up after summation of certain ladder diagrams. Consideration of chain diagrams is necessary in the microscopic theory of the Landau Fermi liquid and there are many examples of this type.

At the same time, such sequences of diagrams correspond to physically well-defined collective bosonic excitations. In the Fermi liquid theory, the chain diagrams describe particle density and spin excitations that can be spoken of as bosonic quasiparticles. From this and many other examples, where a well-defined physical quasiparticle is represented by an infinite sequence of diagrams, one can conclude that the conventional diagrammatic technique is not always the most convenient technique for the description of interacting many body systems.

The inconveniences of the conventional perturbative methods are especially evident in situations when the interaction between the quasiparticles is important. This happens in one-dimensional systems, but models in higher dimensions, , used for describing high superconducting cuprateslee (), heavy fermionsnorman (), etc., can hardly be successfully studied by conventional diagrammatic approaches either. These difficulties call for constructing different approaches dealing with collective excitations rather than with single particles.

Such an approach called bosonization is well known for one dimensional () systems. A huge number of publications are devoted to it (see, for a review, Refs. tsvelik, ; giamarchi, ). Attempts to bosonize higher dimensional fermionic systems have been undertaken in the past starting from the work by Luther luther () where the Fermi surface of a special form (square or cubic) was considered. This idea was further developed by Haldanehaldane () who suggested to patch the Fermi surface by a large number of segments and continued in a number of publicationshoughton1 (); neto (); kopietz (); kopietzb (); khveshchenko (); khveshchenko1 (); castellani (). All these bosonization schemes are based on the assumption that only small momenta are transferred by the electron-electron interaction such that, after any scattering event, the electron remains in the same segment of the Fermi surface. This case is, of course, relevant to a long range interaction but a generalization to an arbitrary interaction and arbitrary momentum transfer is hardly possible within this scheme. At the end one can reproduce results of the random phase approximation (RPA) that can be interpreted in terms of quasiparticles but the interaction between these quasiparticles cannot be taken into account properly.

Another method based on quasiclassical equations has been developed recently for an arbitrary interactionaleiner (). Using this approach anomalous contributions to the specific heat and magnetic spin susceptibilityschwiete () have been calculated in all dimensions and new logarithmic contributions were found. In , the results obtained are in agreement with those known from renormalization group considerationsdl () and exact solutions for spin chainsaffleck (); lukyanov (). They have also been confirmed in by the conventional diagrammatic expansioncms () where the first non-vanishing logarithmic contribution was calculated. However, the same diagrammatic method did not exactly reproduce chubukov () the logarithmic contributions to for , where is the anomalous correction to the specific heat. This discrepancy has been attributed to the observation that not all the effects of Fermi surface curvature have been taken into account in Ref. aleiner, .

One can conclude from this discussion that the problem of the bosonization in has not been completely solved yet but such a method is very desirable because it would be a new tool for analytical calculations. Remarkably, fermionic models are shown to be problematic for numerical computations as well. The powerful Monte Carlo (MC) method suffers from the famous fermionic sign problemblankenbecler (); hirsch (); linden (); dosantos (); troyer (). The MC method implies that the partition function can be represented as a sum of terms with positive probabilities and in this case, the computation time is proportional to a power of the size of the system. This makes the MC method very robust compared to other approaches. However, the fermion determinant is negative for some paths and, actually, its average sign becomes zero in the sampling process. As a result, the computation time grows exponentially with the size of the system or inverse temperature and the advantages of MC get lost.

We are not able to list here all publications where the sign problem was discussed but (except special cases like electron systems with attraction, systems with repulsion but with half filling and some others) the solution has not been found. Moreover, according to Ref. troyer, , the sign problem should be NP-hard cook (), which means that its resolution is almost impossible. We are not aware about any mathematically rigorous proof of the NP-hardness of the fermionic sign problem but it is anyway clear that it is the main problem for MC simulations.

This problem does not arise in (non-frustrated) bosonic systems, though. Therefore, it is quite natural to think about a possibility of circumventing the sign problem by mapping the initial fermionic model onto a bosonic one. Of course, it is impossible to convert real fermions into real bosons but one can think of recasting the interacting part of the fermionic action in the language of collective bosonic excitations. This leads us again to the idea of the bosonization.

In this paper, we suggest a new bosonization scheme that allows one to map in the fermionic model to a model describing bosonic excitations exactly in the sense that the mapping works in any dimension at any temperature and for arbitrary interactions. It can be written in a form of a model of non-interacting bosons in an effective Hubbard-Stratonovich (HS) field. The transformation from the fermionic system to the bosonic one corresponds to using in quantum mechanics the density matrix instead of the wave functions. We argue that the representation of the model in terms of bosons in the external field should be free of the negative sign problem and can be used for MC simulations.

It should be emphasized, though, that for finite systems used in the MC simulations the new bosonized model is not exactly equivalent to the initial fermionic one. In order to derive the boson model, we need a certain regularization and a bath attached to the system. This is necessary to avoid some uncertainties in the bosonized expressions. Only after the regularization is carried out we obtain a bosonic model free of the sign problem. Completely exact mapping of the fermionic model onto the bosonic one is impossible because, technically, several electron Green functions with coinciding Matsubara frequencies and momenta cannot be converted into bosonic excitations. If we did only identical transformations we would not be able to get rid of the sign problem. At the same time, one can see from the derivation that considering a sufficiently large system with the bath any desired precision can be achieved.

Following an alternative route one can average over the HS field before making any approximation. This goal is achieved using the famous BRST (Becci-Rouet-Stora-Tyutin) transformation based on the introduction of superfieldsbrst (). This allows us to represent the partition function in such a form that the average over the HS field can be done immediately. After the averaging one comes to a model of bosons with quadratic, cubic and quartic interactions. We show how this model can be studied with the help of perturbation theory in the interaction terms and compare the first orders with the corresponding terms of the conventional perturbation theory demonstrating a full agreement.

The paper is organized as follows:

In Section II, we formulate the model of interacting fermions, decouple the interaction with the help of the HS transformation and map the model onto the model of bosons in the external field. We introduce a regularization that allows us to avoid singularities in the equation of motion.

In Section III, we introduce superfields, average over the HS field and derive the model of interacting bosons. After that we develop a diagrammatic technique for the superfield theory describing the bosonic excitations and demonstrate in the first orders in the interaction its agreement with the conventional perturbation theory.

Section IV is devoted to the discussion of the bosonization from the point of view of MC simulations. We discuss in detail how the negative sign contributions originate in the fermionic formulation and why they do not exist in the bosonic regularized model. Properties of the partition functions are considered in the complex plane of the electron-electron interaction, which helps to understand in details what happens in the procedure of the bosonization. We derive formulas directly suitable for numerical simulations and give an example of numerical convergence in a case of static fields.

We discuss the obtained results in Section V.

The main results of this paper have already been presented in a shorter publication epm ().

Ii Fermionic model and bosonization

ii.1 Reduction to fermions in the auxiliary field

Our method of the bosonization is quite general and is applicable to a broad class of models. Therefore we start our considerations with a general model of interacting fermions on a lattice of an arbitrary dimension. We write the Hamiltonian of the system as follows


where is the bare part describing non-interacting electrons,


with a chemical potential , a hopping amplitude between the sites and and the spin .

The term describes the electron-electron interaction,


The partition function can obtained from the Hamiltonian using the standard relation


where is temperature.

In order to proceed with the mapping of the fermion model onto a boson on, we will decouple the interaction term by integration over an auxiliary HS field. In principle, the decoupling can be performed in different ways and we should choose one of them. The choice can be different depending on the problem studied. In this work we are mostly interested in systems with an on-site repulsion. We also keep in mind a possibility of using the method for MC computations. The latter desire enforces us to use real HS fields depending only on one coordinate .

At the same time, it would be advantageous to keep an arbitrary sign of the off-diagonal interaction The on-site attraction is less interesting for us because MC schemes can be free of sign in this case dosantos () anyway. The bosonization procedure can be carried out also in this case with minimal changes but we do not do this here.

In order to carry out the HS transformation with the real fields we split the matrix  into its diagonal and off-diagonal parts,


Assuming that  is bounded from above, we can add a sufficiently large constant  to , ensuring the matrix


to be positive definite. (Its Fourier transform should be positive for all momenta .) Introducing the coupling constant


assuming it is site independent, and using the fermionic commutation relations we rewrite the interaction term as


while the chemical potential should be replaced by .

Due a somewhat arbitrary choice of the constant the coupling constants and are also not uniquely defined. However, this does not create any problems.

As the interaction term does not commute with the HS transformation can be performed first subdividing the interval into segments of the length with and then integrating over the auxiliary field for each .

In this Section we are interested in analytical calculations and therefore we write all formulas in the continuous limit () with respect to the “imaginary time” . Transformations for the discrete time (finite ) will be performed in Section IV

We rewrite the partition function , Eq. (4), as


where is the time ordering operator (the time grows from the right to the left) and is the operator in the interaction representation,


The term can be decoupled as follows


where and are creation and annihilation operators in the interaction representation [introduced analogously to Eq. (10)], is a real field with the bosonic periodicity, .

Before decoupling the term Eq. (8), we rewrite it in the form


where , is the fermion density per spin direction, and is the number of sites in the system.

The second term in Eq. (13) renormalizes once more the chemical potential, such that now it equals


The last term in Eq. (13), , is a trivial contribution to the thermodynamic potential .

The fermion density should be calculated using the shifted chemical potential , Eq. (14), and, thus solving a self-consistency equation. We assume that this has been done and consider the first term in in Eq. (13) decoupling it by integration over another field ,



and is also a real periodic field, . The matrix is a matrix inverse to the matrix . The integral over in Eqs. (II.1, 16) converges because is positive definite.

Introducing the field


we write finally the partition function , Eq. (9), in the form


with equal to


and the Hamiltonian now given by


with the modified chemical potential from Eq. (14).

We calculate the trace over the fermionic operators introducing an additional variable , replacing by in Eq. (19) and writing in the form

In Eq. (II.1), is the partition function of the ideal Fermi gas,


where are the eigenvalues of the Hamiltonian , Eq. (20), and


is the electron Green function in the field . In Eq. (22), the operators and are written in the interaction representation using , Eq. (20), as the bare Hamiltonian and is the Green function of the ideal Fermi gas described by the Hamiltonian . The symbol in Eq. (23) stands for the Gibbs averaging over the states of the Hamiltonian , Eq. (20).

The Green function Eq. (23), satisfies the following equation


where for an arbitrary function .

A conjugated equation can be written as


The electron Green function , Eq. (23), satisfies the fermionic boundary conditions

Eqs. (II.1, 22, 24, 25) can serve as the starting point of our bosonization scheme. In the next subsection, we derive equations for the bosonic excitations in an external time dependent HS field.

ii.2 Equations for bosons in the fluctuating field

ii.2.1 General equations.

According to the results of the previous subsection one can calculate the partition function solving Eq. (24) or (25) for a given configuration of the field , substitute the solution into Eq. (II.1), thus obtaining the functional , and then calculate using Eq. (18). Actually, solving Eqs. (24, 25) one obtains more information than necessary because the Green function that should be obtained from this equations depends on two times, and whereas these time should be taken equal in Eq. (II.1).

It is not difficult to derive closed equations for the Green function entering Eq. (II.1). Subtracting Eq. (24) from Eq. (25) and putting , we obtain


We write here equal times in the functions because we assume that they are continuous. This is not so for the Green functions because they have a jump at equal times. For discontinuous functions one would have to take slightly different times and , too.

It follows from Eq. (26) that


which shows that the total number of particles in the system of the non-interacting fermions in the auxiliary HS field does not depend on time. To avoid confusion, we remind the reader that we started from the grand canonical formulation for the model of the interacting particles, Eqs. (1, 4), and writing about the particle conservation now we do not mean that we change the formulation of our initial model.

For subsequent manipulations, it is convenient to introduce a new function


where .

The function has the property


is periodic,


and, hence, describes bosons.

With this function, the functional takes a simple form


Eq. (31) is a reformulation of (II.1) in terms of , Eq. (28).

Having written Eq. (26) we can derive a closed equation for the function Eq. (28). Our procedure is very similar to deriving the kinetic equation starting from equations for Green functions keldysh ().

We rewrite Eq. (26) at , and subtract it from Eq. (26). Using the definition of , we come to the following equation for this function


where is the Fermi distribution function of the ideal gas in the coordinate representation


The solution of Eq. (32), as it stands, is not unique. One can easily see that there is a non-physical solution for any . Moreover, for static fields , one can add to the solution of Eq. (32) a combination of the form


where are eigenfunctions of the operator and are time independent coefficients.

There can be less trivial solutions of the homogeneous equation corresponding to Eq. (32) when is a discontinuous function of time with large jumps. Such solutions correspond to poles in the Green functions and should be treated very carefully. The sign problem in the MC simulations appears as a result of the poles in the integrals over in Eqs. (31) which arise when the fields are discontinuous in time. We discuss this question in Section IV.

It is clear that the ambiguity in solving Eq. (32) following from the existence of solutions of the homogeneous equation creates problems in both numerical and analytical treatment. This means that a procedure fixing the proper solutions of Eq. (32) is necessary and we introduce it below.

ii.2.2 Regularization of the bosonized model.

Figure 1: Simple electronic diagrams.

Seeking for the unique correct solution of Eq. (32) we are guided by the conventional diagrammatic technique for fermions which is well defined and all corresponding diagrams can be calculated at least in principle. In order to visualize how we come to solving Eq. (32) instead of summing the conventional diagrams we consider the simplest loop represented in Fig. 1 (a). The expression corresponding to this diagram can be written as


where and and are fermionic and bosonic frequencies, respectively.

The expression in the second line of Eq. (35) corresponds to the solution of Eq. (32) in the lowest order in , which demonstrates how the bosonic modes are obtained from the fermionic lines.

However, the first and the second lines in Eq. (35) are different at The first line gives the density of states at the Fermi surface


but the expression in the second line is not defined. Of course, one could speak about the limit instead of just putting but this is not justified for any finite physical system with discrete energy levels.

The uncertainty of the expression in the second line of Eq. (35) is reflected in the existence of solutions [like those given by Eq. (34)] of the homogeneous equation corresponding to Eq. (32).

We see that our bosonization scheme is not exactly equivalent to the initial fermionic model and needs a regularization in order to avoid the uncertainty. This uncertainty appears not only in the situation represented in Fig. 1 (a), where the field enters with . The same problem is encountered when calculating the contribution of, e.g., Fig. 1 (b), using the bosonization scheme because this graph contains fermionic Green functions at coinciding momenta and frequencies but the interaction lines may carry any momenta and frequencies. We will come back to this point when discussing in detail the diagrammatics of the bosonized theory.

The hint at a possible regularization is given by the form of Eq. (35). If the momenta were continuous rather then discrete and we could simply neglect the contribution of the state with , the bosonization scheme would work in this order. This is not sufficient, though, because the momenta and frequencies of the two horizontal fermion lines in Fig. 1 (b) obtained after integration over the HS field coincide for any translationally invariant systems and one should slightly violate the translational invariance in order to split the momenta.

This goal can be achieved if we assume that the system of the interacting electrons considered here is imbedded in a “bath”. We introduce this bath considering a model of the interacting electrons on a -dimensional lattice with the total number of sites . However, we assume that the interaction , Eq. (3), entering the Hamiltonian , Eq. (1), vanishes outside a subsystem consisting of a considerably smaller number of sites . As an example, we can suggest the following form of the interaction


which means that the interaction is finite inside the sphere of the radius containing sites and vanishes outside this sphere. In other words, we attach metallic leads to the system of the interacting electrons and assume that the inter-electron interaction vanishes in the leads.

It is clear that the leads cannot change physics of the system of the interacting electrons if the latter is sufficiently large. However, this model is reasonable even for a system with a small number sites (quantum dot), although properties of the systems with and without leads can be different.

The form of the interaction , Eq. (37), formally violates the translational invariance and the matrix elements


are finite for .

Of course, the matrix element with remains finite and, for such momenta, the two Green functions in Fig. 1 (b) still have the same momenta and frequencies. However, we can neglect them because they give a small contribution provided the number of the contributing momenta is large. The latter is achieved for a large number of the sites in the bath.

We emphasize, that neglecting the contribution of graphs containing several electron Green functions with coinciding momenta and frequencies is possible because their contribution is not singular. One can estimate the relative error of this approximation as being proportional to , where is the total number of sites in the entire system including the bath.

Following this idea we could simply exclude in the bosonization approach all propagators with “by hand”. However, this is possible only when doing a perturbation theory and such an approach is not sufficient for non-perturbative and numerical studies.

The bosonic modes with correspond to solutions of the homogeneous equation in Eq. (32). In order to discard these modes by a regular procedure we regularize Eq. (32) slightly changing it. This is not a trivial task because this slight modification should guarantee the complete absence of the solutions of the homogeneous equation corresponding to Eq. (32).

We regularize the bosonized model replacing Eq. (32) by a system of two equations


where , is a matrix



The functions satisfy the bosonic boundary conditions


and Pauli matrices are used


The real parameter should be put to zero, , at the end of calculations.

With this modification we write the function in the form


The exponent in Eq. (43) is an even function of . It is not difficult to see that putting in Eqs. (39-43) one returns to Eqs. (31, 32). However, keeping in Eqs. (39-43) the parameter finite allows us to discard all the solutions of the homogeneous equation in Eq. (32).

This can be understood rather easily because the operator is hermitian and, therefore, its eigenvalues are real. Moreover, as will be shown in the next subsection, they appear in pairs: if an eigenvalue exists, then the eigenvalue exists, too. If, when changing parameters of the operator , an eigenvalue with its counterpart turned to zero at some point, this would mean that after crossing this point they become imaginary. However, the latter is forbidden by the hermiticity of the operator . A general proof of the absence of the solutions of the homogeneous equation in Eq. (39) or, in other words, absence of zero eigenvalues of the operator , Eq. (40), is given in Appendix A.

The absence of zero eigenvalues makes the operator invertible. Since all the parameters of Eqs. (39-40) are real, the solutions are real, too. This means that the exponent in Eq. (43) is real and is real and positive. This property of the can be very important for numerical computations using the MC method.

Using the regularization with the parameter , Eqs. (39-43), we can prove a stronger than Eq. (29) relation for ,


Eq. (44) is fulfilled automatically for any solution of Eqs. (39-40). This property can be proven putting in Eqs. (39, 40) and summing over . Then, one obtains the equations


where , .

Only the trivial solution satifies the boundary condition, Eq. (41), leading to Eq. (44).

ii.2.3 Spectral expansion.

In this subsection we discuss important properties of the solutions of Eqs. (39, 40) using spectral expansions in eigenfunctions of this equation. In order to better understand the difference between solutions of Eq. (32) and Eqs. (39, 40), let us consider both the equations.

The operator in the L.H.S. of Eq. (32) is not hermitian, which results in a double set of eigenfunctions and satisfying the equations


with the boundary conditions , . Two other equations can be written taking the complex conjugate of Eqs. (46),


The orthogonality conditions for the functions can be written as

(the same for the complex conjugates).

The eigenvalues are generally complex and is not excluded. For example, the function , Eq. (34), corresponds to the zero eigenvalue for static .

All non-zero eigenvalues appear in pairs. If an eigenvalue with an eigenfunction exists, then the eigenvalue with the eigenfunction exists, too. We see that nothing prevents the eigenvalues from turning at some points to zero and being complex, which is an unpleasant feature of the theory.

In contrast, the operator is hermitian. The structure of the two-component vectors can be represented as

and they can be found from the equation


supplemented by the boundary condition , where the eigenvalues must be real.

The eigenstates with eigenfunctions and eigenvalues have their counterparts with the eigenfunctions and eigenvalues .

The orthogonality condition for the eigenvectors follows from Eq. (48) and can be written as


where the symbol stands for the complex conjugation and transposition .

The completeness of the system of the eigenvectors expressed by the relation


allows us to write the spectral expansion for the matrix Green function satisfying the equation


in the form


Since all eigenvalues are not equal to zero, the Green function , Eq. (52), is uniquely defined and does not have singularities. Using the properties of the eigenvectors explained after Eq. (48) we obtain a symmetry property for these functions


With the Green function Eqs. (51, 52), we write the function , Eq. (43), in the form