Quantum impurity in a Tomonaga-Luttinger liquid: continuous-time quantum Monte Carlo approach

Quantum impurity in a Tomonaga-Luttinger liquid: continuous-time quantum Monte Carlo approach


We develop a continuous-time quantum Monte Carlo (CTQMC) method for quantum impurities coupled to interacting quantum wires described by a Tomonaga-Luttinger liquid. The method is negative-sign free for any values of the Tomonaga-Luttinger parameter, which is rigorously proved, and thus, efficient low-temperature calculations are possible. Duality between electrons and bosons in one dimensional systems allows us to construct a simple formula for the CTQMC algorithm in these systems. We show that the CTQMC for Tomonaga-Luttinger liquids can be implemented with only minor modifications of previous CTQMC codes developed for impurities coupled to non-interacting fermions. We apply this method to the Kane-Fisher model of a potential scatterer in a spin-less quantum wire and to a single spin coupled with the edge state of a two-dimensional topological insulator assuming an anisotropic XXZ coupling. Various dynamical response functions such as the electron Green’s function and spin-spin correlation functions are calculated numerically and their scaling properties are discussed.



Electronic correlations play fundamental roles in determining low-energy phenomena in one-dimensional electron systems [1]. Bosonization is a powerful technique to treat such correlations exactly. In the presence of impurities, however, it is well known that impurities in one-dimensional systems drastically influence transport properties of the systems. Such an example has been found in a classical model by Kane and Fisher in their pioneering work about a backward scattering potential problem in a spinless quantum wire [3]. There, for the case of repulsive interaction [Tomonaga-Luttinger (TL) parameter ], the conductance vanishes at zero temperature () and the potential barrier becomes infinitely strong and cut the wire into two parts, while for attractive cases with , the potential becomes zero in the low-energy limit and there remains a finite value of the conductance at .

For acquiring knowledge about thermodynamic, transport, and dynamical properties of such systems, bosonization combined with perturbative renormalization group methods [3], Bethe ansatz [6], and functional renormalization group [7] have been intensively used so far. To obtain numerically exact results for bosonized impurity problems, path integral Monte Carlo approaches have been employed [8]. A bosonic numerical renormalization group method [11] is also a powerful technique to investigate their low-energy properties. While they are very useful approaches, there is still need for even more powerful numerical approaches that allow to compute wide range of temperature properties and dynamical correlation functions even for more complex models in an exact way.

To this end, in this paper, we will develop a continuous-time quantum Monte Carlo (CTQMC) method [12] in the Tomonaga-Luttinger liquid (TLL) in one-dimensional systems coupled to an impurity. The CTQMC previously has been developed to describe quantum impurities coupled to non-interacting environments. It has mainly been used for fermionic systems and extensively used in the frame work of the dynamical mean field theory [16] as an exact numerical solver for the effective-impurity problem in it. Recent development [17] of the algorithm also enables us to treat bosonic systems and mixture of bosons and fermions. The advantage of the CTQMC is that this allows us to calculate various quantities at low temperatures in efficient ways and for some simple models there is no negative sign problem [13].

Our algorithm of CTQMC for TLL has advantages in the following points. (i) Bosonization allows us to treat correlation arising from strong interactions in the environment exactly. (ii) There is no negative sign problem for any parameters, which is, indeed, proved analytically. This enables us to carry out low-temperature analysis with high precision. (iii) There are close relations to the fermionic version of CTQMC, although the whole algorithm is written in the bosonization language. This enables ones to implement the CTQMC for the TLL easily from their fermionic CTQMC code. (iv) The method can be applicable to not only potential scattering problems but also to Kondo-type problems [19] without a negative sign problem. (v) The electron Green’s functions, the boson-boson correlations, conductance, the spin-spin correlation functions, and various local correlators are calculable. (vi) Compared with lattice QMCs, our formalism is free from finite size effects at low temperatures and phase spaces for the random walk is expected to be much smaller, and, thus, the computational cost is much lower.

This paper is organized as follows. In Section 2, we will explain the models used in this paper. Section Section 3 will be devoted to illustrate our algorithm of the CTQMC for TLL. The method will be applied to the Kane-Fisher model [3] in Section 4.1 and the XXZ Kondo problem [19] in Section 4.2. We will discuss possible extension in Section 5 and summarize the present results in Section 6.


In this section, we will introduce our model. First, we will show a one-dimensional Tomonaga-Luttinger liquid Hamiltonian and explain our notation of the bosonization we will use throughout this paper. The second part is an introduction of impurity-electron interactions. We will use a general expression that can be used in two models we will discuss in Section 4.

2.1One-dimensional bulk Hamiltonian

We consider spin-less fermions in one-dimensional (1D) systems whose non-interacting Hamiltonian is given by [2]

where is the fermion creation operator at the position and () refer to left(right)-moving component. :: indicates the normal ordering of the operator , is the system size, and is the Fermi velocity. The fermion field satisfies the anticommutation relation

Following the standard bosonization procedure [2], we define bosons as

where is the short-distance cutoff. The Bose fields satisfy

where the term is explicitly written. Two Klein factors and have been introduced in Eq. (Equation 3) to reproduce the anticommutation relation of [Eq. (Equation 2)]. Their anticommutation relation is


Note that , and the two bosons are independent fields commuting with each other and with the Klein factors, which is a physically correct description as is evident from the definition of and [4].

Electron-electron interactions are easily taken into account in the bosonized theory and then the bosonized Hamiltonian reads

with and is the TL parameter that characterizes the bosonic theory: corresponds to the noninteracting case and () describes repulsive (attractive) interactions, respectively. The velocity is now renormalized as . Throughout this paper, we will be interested in the repulsive case.

For later purposes, we introduce another representation following Delft and Schoeller as [4]

Then, the Hamiltonian (Equation 4) is rewritten as

At , a simple relation holds [4],

2.2Impurity potentials

Now, we introduce interactions between a quantum impurity located at and the interacting electrons of the TLL. We consider a coupling by single-particle scattering (generalization to more complicated interactions is straightforward but, of course, model dependent). The interactions, , are decomposed into two parts, a forward-scattering channel, , and a backward-scattering channel described by ,

The scattering of electrons can change the state of the impurity (e.g., flip a spin). This is described by the impurity operators and . They will be discussed in later sections. In the most general case, their form can be derived by first bosonizing the model and then identifying the two terms discussed above by comparing them to the bosonized version of Eqs. (Equation 6) and ( ?) discussed below. Note that has the dimension of [energy][length] and ’s are dimensionless operators.

In terms of the bosons, Eqs. (Equation 6) and ( ?) read

where , , and . The vertex operator is defined as

This normalization of the vertex operator leads to following the bare (i.e., in the absence of the impurity) two-point correlator as a function of imaginary time ,

at for and [see also the definition of multipoint correlators in Eq. (Equation 12) below].

3Continuous-time Quantum Monte Carlo Method

In this section, we will explain how continuous-time quantum Monte Carlo method can be applied to the impurity problem in the TLLs. We will demonstrate that the configuration weight for a given snapshot is easily calculated by the technique developed in fermionic CTQMCs. We will therefore omit detailed explanations about update operations, since these are essentially the same as in the fermionic CTQMCs [15].

3.1Partition function

We want to evaluate the partition function ,

within a Monte Carlo approach. Here, the “noninteracting” part is the sum of the one-dimensional TLL and the local impurity Hamiltonian, . In this paper, we will analyze models with (e.g., a magnetic impurity in the absence of magnetic fields). Via perturbative expansion of , we can express as

where and and indicates the time-ordered product. In this paper, we will discuss situations where the forward-scattering part () can be eliminated by an appropriate unitary transformation or where due to symmetry requirements. Thus, we retain only and in order to distinguish the two terms in , we define

A general th order term in the partition function is expressed as

Due to the fermion number conservation encoded by the Klein factors, the number of in Eq. (Equation 9) has to be the same as that for , and therefore only even terms with , contribute. Now, consider a fixed series of times . Then becomes

with , , or . The partition sum is obtained by averaging over all , all times, and all using a Monte Carlo procedure.

Since the product of Klein factors gives factor unity, we can write Eq. (Equation 10) as a product of a TLL correlator and a correlator involving only impurity operators

In the following subsections, we will analyze the two sectors in details.

3.2Impurity average

Here, we discuss the local part . is the time-ordered product of ’s:

Here, is the average with respect to the impurity Hamiltonian. As noted above, both and appear times and, for later convenience, we define new indices with , such that the operators () are evaluated at the times () with time-ordering in each index, .

3.3Boson average

For the bosonic part, is the time-ordered -point correlation function of :

Here, the boson average is evaluated using the Gaussian TLL Hamiltonian (Equation 5). It is well known that the correlation function of vertex operators,

are calculated as [4]


Here, , and, in order to prevent the divergence, the cutoff function is necessary and it satisfies

In actual calculations, we will use the following function throughout this paper:

For very high temperatures (not considered in this paper), it is sometimes useful to use a smooth function in order to remove the discontinuity appearing in physical quantities such as

with being a positive constant. In the limit, Eq. (Equation 12) vanishes unless . Thus, a “neutrality condition,” , has to be fulfilled. In our case, this is automatically enforced by the fermion number conservation and we obtain

with . An important observation is that Eq. (Equation 14) is positive definite, and, thus, our Monte Carlo method is negative-sign free if .

Equation (Equation 14) can be further simplified via the “generalized” Wick’s theorem [4], which is valid if and only if . We utilize this theorem, although actual numerical calculations are done with finite . The theorem might be most easily obtained by comparing the partition function for noninteracting spinless fermion in one dimension and that in the bosonization representation. The result is

The matrix is given by

and the index corresponds to . This form is particularly useful, since we can use the fast-update algorithm developed in the conventional fermionic CTQMC methods [12].


In this section, we will apply our CTQMC method to two models. One is the Kane-Fisher model describing a backward-scattering impurity potential in a (spinless) quantum wire [3]. The other is the XXZ Kondo problem [19] in helical liquids, i.e., on the edge of two-dimensional topological insulators.

4.1Kane-Fisher model

The Kane-Fisher model is defined by considering forward scattering with and in Eq. (Equation 6) and by setting in Eq. ( ?) as a potential scatterer has no internal degrees of freedom. Since contains only and only , we can separately analyze the two. The part is trivial because it can be absorbed into terms in by a unitary transformation [4]. Thus, in the following, we analyze part in details. Note that because and thus , we can use the positivity of (Equation 14) to conclude immediately that there is no negative-sign problem. Throughout this subsection, we set and fix for the unit of energy, where is the relevant microscopic unit of length, which is, for example, set by the typical width of the potential.

The model itself has been extensively analyzed by various authors [3], and now its low-energy properties are well understood. We will study this problem as a benchmark of our algorithm. We will discuss a physical quantity that has not been investigated so far: the electron Green’s function in imaginary time. This is the most natural quantity for imaginary-time algorithms like CTQMC. Though it is not a directly measurable quantity in experiments, the numerically-exact results can be used to obtain the density of states by using analytical continuation techniques [22] (not covered in this paper).

Electron Green’s function

Let us consider the local Green’s function for ,


Here, is the average over Klein factors. The correlation function for the part is trivial, leading to

Here is given by Eq. (Equation 13). Note that for our model nonlocal Green’s function can be expressed in terms of local correlators as the bosons of the TLL are noninteracting. In the following, we explain how one can calculate the and the Klein factor parts in our CTQMC.

In the Monte Carlo simulations, time-ordered averages of an operator is estimated as

where is the number of Monte Carlo samplings and we have defined an operator form of , see Eq. (Equation 10). For , this is given by

For the Green’s function , we need to calculate Eq. (Equation 17) with with . As derived in Appendix A, we need to sample the following quantity for ,

Here, is the number of vertices between and in the MC snapshot and a similar expression is obtained for . The notation represents that and are added to . Note that to derive Eq. (Equation 18), we have used the generalized Wick’s theorem mentioned before. Then we obtain

and a similar expression is applied to .

An alternative approach is, however, more efficient in regimes where high orders of perturbation theory are needed. Following Refs. [12], we can also derive an alternative expression for calculating the Green’s function . The quantity that corresponds to Eq. (Equation 18) is now given by

Note that in Eq. (Equation 20), and are chosen in a given snapshot , while in Eq. (Equation 18), and are external ones. This implies that by computing one snapshot with pairs of time variables one obtains contribution to the Green’s function for about different , which helps to reduce the statistical error. The notation represents that and are removed from .

Since the ratio of two determinants in Eq. (Equation 20) is simply , which is calculated in every MC process, this also reduces computational costs [12]. We can also derive a similar expression for , . Summing over all possible combinations for a given snapshot at th order and dividing by , we obtain

The two alternative formulas (Equation 19) and (Equation 21) can be used for checking the program code.

Bench mark for

Figure 1: (Color online) G_L^+(\tau) vs \tau/\beta for several coupling constants \lambda_B=0.1(top)-0.5(bottom) and inverse temperatures \beta. The exact result (dashed line) is compared to two numerical methods [points with error bars: Eq. (), solid lines: Eq. ()]. A comparison of panel (a) (cutoff a=1) to panel (b) (cutoff a=0.25) shows that small deviations from the exact result vanish for small a. Panel (c) shows that highly accurate results can be obtained even for very low T.
Figure 1: (Color online) vs for several coupling constants (top)-(bottom) and inverse temperatures . The exact result (dashed line) is compared to two numerical methods [points with error bars: Eq. (), solid lines: Eq. ()]. A comparison of panel (a) (cutoff ) to panel (b) (cutoff ) shows that small deviations from the exact result vanish for small . Panel (c) shows that highly accurate results can be obtained even for very low .

In this subsection, we show the results for , i.e., a system of noninteracting electrons. We compare the electron Green’s function obtained in the CTQMC and the exact results as follows:

which can be easily obtained from the equations of motion for the Green’s functions.

Figure 1 shows as a function of for and several parameters and . In each plot, the exact result (dashed lines) and the result of the two methods described above are shown. The points with error bars (indicating the statistical error arising from the Monte Carlo sampling) are obtained from Eq. (Equation 18), while the solid lines have been calculated from Eq. (Equation 20) (the statistical error can be read off from the size of the noise in the curves). As one can see, the numerical data and the exact results are consistent with each other. More precisely, the Green’s functions are only identical in the limit [we used this limit both in the derivation of the CTQMC approach and in Eq. (Equation 22)]. Figure 1 shows that tiny systematic deviations of the exact and the numerical result visible for become smaller than the noise for . In the following, we will always use as the universal properties for and discussed in the following are independent of the cutoff.

Figure 1 shows that highly accurate results are also obtained for low . Comparing the two computational methods (using the same computational time), we first note that both give reliable results. Which method is preferable depends in general both on the perturbation order and the type of binning in time used to extract data. For the parameter regime used in our calculations, we found the second approach to be more efficient. For very high orders of perturbation theory and small number of bins, however, the first approach can beat the second one in efficiency. In Section 5.1, we will discuss that in regimes, where the nonlinearities are irrelevant (attractive interactions), the second method is inefficient.

Universal scaling function for electron Green’s function

The main prediction of Kane and Fisher [3] is that for repulsive interactions, , even a weak impurity effectively cuts the chain: Electrons scatter so efficiently from the slowly decaying Friedel oscillations that for and at the Fermi energy one obtains perfect reflection. The fact that the impurity cuts the quantum wire can also be measured by tunneling spectroscopy, i.e., by considering the local Green’s function close to the impurity. Based on the assumption that the wire is perfectly cut by the impurity, one expects for ,

as has been derived by Furusaki [5]. This prediction can be checked analytically for , where an exact analytic result can be derived [4]. Equation (Equation 23) should be compared to obtained for , in the absence of the impurity.

Note that for the computation of the physical electron Green’s function one has to consider a further contribution, , in addition to . It is possible to calculate using our approach, but we do not discuss it here for simplicity.

From general scaling arguments and the analysis of Kane and Fisher [3], one expects for a weak potential scatterer (small ) a crossover from to described by a universal (but -dependent) scaling function ,

where is the Green’s function for , see Eq. ( ?). All dependence on the strength of the impurity potential is thereby encoded in the characteristic energy scale with

for small . The universal scaling form (Equation 24) is expected to be valid whenever is much smaller than the cutoff energy . For and weak , the short time dynamics is determined by the noninteracting result, , while , see Eq. (Equation 23).

In the following, we show our CTQMC results, which confirm the expected behavior and allow us to calculate the full scaling function describing the crossover from weak to strong coupling. To our knowledge, this is the first demonstration of the numerically-exact Green’s function in this model.

Figure 2 shows versus for several parameter sets and (a) , (b) , and (c) for various temperatures . Our numerical results reproduce the analytically expected behaviors: First, for wide ranges of the curves scale on top of each other (we have not used an appropriately rescaled temperature, therefore the upturns occur at different points). Second, we obtain the analytically expected asymptotic behavior with , while , see Eq. (Equation 23). Third, our result provides the full crossover function from weak to strong coupling.

Figure 2: (Color online) G_L^+(\tau)/[s(\tau)]^{-g/2} versus T^*\tau for various parameter sets (\beta,\lambda_B), a=1 and (a) g=0.3, (b) g=0.5, and (c) g=0.75. Straight lines show the \tau^{(g-1/g)/2} dependence expected from the fixed point where the impurity cuts the chain (the factor \tau^{g/2} originates from the asymptotic form of [s(\tau)]^{-g/2}).
Figure 2: (Color online) versus for various parameter sets , and (a) , (b) , and (c) . Straight lines show the dependence expected from the fixed point where the impurity cuts the chain (the factor originates from the asymptotic form of ).
Figure 3: (Color online) G_L^+(\tau)/[s(\tau)]^{-g/2} vs \tau/\beta for a=1 and various parameters (\beta,\lambda_B) keeping the ratio T/T^*\simeq 0.014 fixed.
Figure 3: (Color online) vs for and various parameters keeping the ratio fixed.

To prove that scaling works also at finite , we show in Figure 3 as a function of for a wide range of coupling constants using a fixed ratio of . The perfect collapse of the data shows that temperature only enters in the combination as predicted by Eq. (Equation 24).

Finally, we show the scaling functions for , , and in Figure 4. The curves can simply be obtained from the small data (we use ) shown in Figure 2, which are independent of within the scatter of the curve. For close to , becomes exponentially small and it becomes more difficult to extract the low and long results.

Figure 4: (Color online) {\mathcal F}_g(x,0) vs x for g=0.75, g=0.5, and g=0.3 with a=1. Data for \tau/\beta<1/6 of each of the lines in Fig.  are used in order to extract the T=0 limit.
Figure 4: (Color online) vs for , , and with . Data for of each of the lines in Fig. are used in order to extract the limit.

4.2XXZ Kondo model in helical liquids

Along edges of two-dimensional topological insulators, a special one-dimensional electron system is realized [19]. Namely right- and left-moving electrons have opposite spin polarizations, up and down, respectively. The topological protection of these edge channels is reflected by unusual scattering properties: due to time-reversal symmetry a static potential cannot scatter a right-moving spin-up electron into a left-moving spin-down electron.

The situation differs in the presence of a magnetic impurity. Using a spin-flip process, a right mover can be converted in a left mover (and vice versa) due to the exchange interaction with the quantum impurity. Therefore, it is an interesting problem to study magnetic quantum impurities at the edge of a topological insulator to investigate whether and how topological protection is affected by their presence.

In this subsection, we examine the spin-1/2 XXZ Kondo model [19]. We restrict our analysis to the case where the total spin in the direction is conserved. Although this symmetry is broken in real materials, e.g., by Rashba interactions [23], it is important to clarify also the basic properties of this simplified problem. Note that the transport properties in the presence and absence of this symmetry qualitatively differ as the current in a helical edge state (proportional to ) can only degrade by processes where spin conservation is violated.

We will use different units from those in the previous subsection, and use the energy unit for all and as a unit of length, in order to use the same high-energy cutoff as in previous studies [19]. The main results in this subsection are the phase diagram in - plane, which has been discussed perturbatively [19] and the numerically-exact time and temperature dependence of the spin-spin correlation functions for general interaction parameters.


For XXZ Kondo model, and in Eqs. (Equation 6) and ( ?) are given as

We have used a slightly different quantization axis of the impurity spin from in Refs. [20]: and . Since term is a pure potential scattering in the charge sector and equivalently sector, this does not affect the CTQMC and the following discussions, we will concentrate on the sector, , hereafter.

First, we write the Hamiltonian in the bosonization basis

where describes the coupling of the -component of the spin, while parametrizes the strength of spin-flip terms [20].

For the CTQMC, it is useful to transform via a unitary transformation [21],

This erases the term in Eq. (Equation 26), since

Thus, the Hamiltonian is transformed to

with and , where is defined as

As will be discussed in Appendix B, it is sufficient to consider cases for , i.e., .

The CTQMC algorithm for this model is similar to that for the Kane-Fisher model. Indeed, an exact relation between the partition functions of the two models is known [25]. Only even order terms remain finite due to the fact that the impurity spin is 1/2, i.e., and/or the total fermion number conservation. This also restricts configuration space for the impurity spin in . We just need to generate configurations in which and appear alternatively: . Thus, we can use algorithm similar to the “segment” representation used in the Anderson model, which accelerates the acceptance rate in the MC samplings [13].

Spin-spin correlation function

In this subsection, we explain how to calculate the dynamical spin-spin correlation functions.

First, let us discuss the transverse local spin susceptibility, , where is defined as

Noting that , we can calculate

\[\chi^{+-}(\tau)=\frac{1}{\beta} \left\langle \sum_{ij}^k {\mathcal M}_{ij}\Big[\delta(\tau-\tau_{ij})+\delta(\beta+\tau_{ij}-\tau)\Big]\right\rangle, \ \\]

by sampling the following quantity:

where and are chosen in a given snapshot at the th order as in Eq. (Equation 20) and the corresponding vertex operators at and should have and , respectively. For , the same expression holds with regarding now and . We also use symmetry properties to obtain results for . is, indeed, derived in an almost identical way as in Appendix A.2.

Second, as , the longitudinal spin susceptibility is directly evaluated as