The Sampling Procedure

3.3. Continuous-Time Hybridization Expansion Quantum Monte-Carlo Algorithm 59

With the definition ofg_C, it is immediately clear from this expression, that g¯σC(τ) =

Z

dτ⁰dτ⁰⁰∆(τ⁰−τ⁰⁰−τ) 1 w_C

∂w_C

∂Γσ¯(τ⁰⁰−τ⁰)

= Z

dτ⁰dτ⁰⁰∆(τ⁰−τ⁰⁰−τ) det(Mσ¯({τ_l^c†¯σ, τ_l^cσ¯}))∂det(M_σ¯⁻¹({τ^c

†

¯ σ

l , τ_l^c¯σ})

∂Γ_σ¯(τ⁰⁰−τ⁰). (3.90) Combining the rules for functional differentiation, i.e.

∂Γσ(τi−τl)

∂Γ_σ¯(τ⁰⁰−τ⁰⁰)=δ(τi−τl−τ⁰⁰+τ⁰)δσ¯σ, (3.91) and matrix inversion

Mij =det( ˜M⁻¹)(ji)

det(M⁻¹) = 1 det(M⁻¹)

∂det(M⁻¹)

∂(M⁻¹)_ji , (3.92)

where we used the cofactors( ˜M⁻¹)(ji)for the matrix inversion, which has a possible representation as the derivative of the determinant with respect to the corresponding matrix element(M⁻¹)_ji, we finally obtain the formula for the observables:

g¯σC(τ) = ∆(τ_i^cσ¯ −τ^c

†

¯ σ

j −τ)(M¯σ({τ_l^c†¯σ, τ_l^c¯σ}))ji. (3.93) The delta-function in this formula must either be measured on a fine grid inτ-space, i.e. on a discretized time interval[0, β]or, as we will do, the observables have to be performed in another representation which is much more suitable for continuous, smooth functions (see chapter ...). The average occupation hnσifor theσ-component and the average double occupancyhn↑n_↓iare other observables, which can be directly measured in the Monte-Carlo sampling, which becomes clear when taking the derivatives

βhnσi= 1 Z

∂Z

∂µσ

and βhn_↑n_↓i=−1 Z

∂Z

∂U. (3.94)

As above, this leads to the observables nσC = 1

βw_C

∂wC

∂µσ

= τ_σ^occ

β and (n↑n↓)C =− 1 βw_C

∂wC

∂U =τ^double

β , (3.95)

which can be accumulated without any effort during the Monte-Carlo sampling and are determined more accurately than the Green’s functions since they are static variables.

With all the necessary ingredients, namely the configurations, the weights and the observables, we can skip and jump to the next section and set up the Markov-Chain and determine the acceptance probabilities for the sampling process.

Figure 3.9: Illustration of four subsequent possible steps in the Markov-Chain sampling. First, in a), the chain starts with a configurationCconsisting of two segments for up and three segments for spin-down. In the next step of the Markov-Chain, b), a specific spin-down segment is removed. In step c), an anti-segment for the spin-up component is removed, leading to the merging of the two present segments into a single one, and finally, in d), a distinct spin-up segment is inserted again.

Suppose the current configuration C consists ofn_σ segments, occupying theσ-state of the impurity.

Then with probabilityp = 0.5we pick a spin-componentσ = {↑,↓}to update. Possible updates are either to insert a new segment into a free time interval or to remove one of then_σ existing segments from the impurity. For the insertion, which will be chosen with probabilityp= 0.5, a timeτ^c†σhas to be selected at which the creation operatorc^†_σ(τ^c†σ)shall be inserted into the impurity. This time is chosen from the interval[0, β]with uniform distributed probability, i.e. p = ^dτ_β. After the timeτ^c†σ has been determined, the algorithm has to check if the impurity is empty or occupied at this time. If the impurity is occupied, this move is rejected and the sampling then proceeds with the next step, otherwise when the impurity is empty, the lengthlmaxof the empty time interval from timeτ^c†σ to the starting time of the next segment is determined. The annihilation operatorc_σ(τ^cσ)is then inserted at a uniformly chosen timeτ^cσ in the interval[τ^c†σ, τ^c†σ+l_max], i.e. with probabilityp=_l^dτ

max. The proposal probability for the insertion of this segment is then determined by the product of the two probabilitiesp^prop(nσ →nσ+ 1) = _βl^dτ2

max. On

the other hand, fornσ+ 1segments occupying the impurity, the probability to uniformly pick a certain segment to remove is simplyp^prop(nσ+ 1→nσ) =_n¹

σ+1. With this, the detailed balance condition states p^acc(nσ→nσ+ 1)

p^acc(n_σ+ 1→n_σ) = p^prop(nσ+ 1→nσ) p^prop(n_σ→n_σ+ 1)

wn_σ+1

w_nσ

= (nσ+ 1)

βl_max e^{µσδτσocc−U δτdouble}

det(M_σ⁻¹(nσ+ 1)) det(Mσ⁻¹(nσ))

, (3.96)

where in the last step, the formula for the weights (3.84) was inserted and we introduced the quantities

δτ_σ^occ andδτ^doublewhich are the change in the occupation time and double occupation time between the

two configurations. Note that we took the absolute value of the determinant ratio since we have shown previously that the prefactors(−1)^scoming from the determinants cancel with the prefactors from the trace-factor, see (3.81). The Metropolis condition fixes these probabilities to

p^acc(n_σ→n_σ+ 1) = min

1,(nσ+ 1)

βl_max e^{µσδτσocc−U δτdouble}

det(M_σ⁻¹(nσ+ 1)) det(Mσ⁻¹(nσ))

, (3.97) p^acc(nσ+ 1→nσ) = min

1, βlmax

(nσ+ 1)e^{−µσδτσocc+U δτdouble}

det(M_σ⁻¹(nσ)) det(Mσ⁻¹(n_σ+ 1))

. (3.98) By including updates that insert and remove segments, ergodicity has been obtained on a formal level since any configurationC = (n_↑, n_↓, ...)can be reached from any other configurationC˜= (˜n_↑,˜n_↓, ...)

9Please note that the insertion and removal of segments alone is not sufficient to obtain ergodicity, i.e. the insertion and removal of antisegments is required for the Markov chain and not a tool that simply improves the sampling procedure.

3.3. Continuous-Time Hybridization Expansion Quantum Monte-Carlo Algorithm 61

inN = n_↑+n_↓+ ˜n_↑+ ˜n_↓ steps. However, it turns out that the insertion and removal of segments does not guarantee ergodicity in finite sampling processes. This fact can also be conjectured from the fact that a completely empty impurity and a completely filled impurity are not treated on equal footing with these kind of updates as it is very unlikely to insert a segment spanning the whole time interval[0, β].

This can be adjusted by considering a complementary set of updates, namely the insertion and removal of so-called anti-segments. While a segment is nothing but a time interval (τ^c†σ, τ^cσ) for which the impurity is occupied, an anti-segment is just a contrary time interval(τ^cσ, τ^c†σ)for which the impurity is unoccupied. To express this in the segment picture, the insertion of an anti-segment(˜τ^cσ,τ˜^c†σ)is nothing but the splitting of a single segment(τ^c†σ, τ^cσ)into two distinct segments(τ^c†σ,τ˜^cσ) (˜τ^c†σ, τ^cσ), while the removal of an anti-segment(˜τ^cσ,τ˜^c†σ)is the recombination of two subsequent segments into a single one. The update probabilities are then determined by the formulas (3.97) and (3.98) since the anti-segment moves follow completely the same procedure as the segment moves and the weights are identical¹⁰. As for the CT-AUX solver, further updates can also be implemented in the CT-HYB solver, which lead to a decrease of the auto-correlation times but are not required for ergodicity. A set of such updates would be the possibility to not only insert and remove (anti-) segments but also to allow for shifting segments in time or to shift only the starting- or end point of a segment in time. These updates are implemented straight-forwardly by using the detailed balance condition and are not derived here. However, in some situations a complete exchange of the segments for the two different components is required, especially when the system is approaching magnetically ordered phases [67].

Steps of the Sampling Process

The Monte-Carlo flow diagram for the CT-HYB solver is very similar to the diagram for the CT-AUX solver, Fig. 3.4. The steps are as follows:

1. Initialization: At the beginning of the sampling process one possible configuration must be prede-termined as the starting configuration to initialize the sampling procedure.

2. For the current configurationCa possible update has to be chosen, either the insertion/removal of a segment(τ^c†σ, τ^cσ)or insertion/removal of an anti-segment(τ^cσ, τ^c†σ). Where all four of those possibilities have the same probability to not influence the detailed balance formulas above.

3. A new configuration must be proposed, either with one segment(τ^c†σ, τ^cσ)less or more (or one anti-segment less or more), depending on the previous step.

4. The acceptance probability for this new configuration is determined according to the Metropolis probabilities (3.97) and (3.98).

5. Either the new configuration is accepted and the old configuration replaced, or the new configura-tion is rejected, and the old configuraconfigura-tion stays unchanged (i.e. the new configuraconfigura-tion is discarded).

6. A measurement of the observables, i.e. theg_C’s, is performed (these may be many observables).

7. Go back to step 2, until the specified number of iterations has been achieved.

Finally, after deriving the configurations, weights and observables in the previous chapter and setting up the sampling procedure, we will also compare the CT-HYB solver with exact diagonalization to check whether it is working correctly or not. In Fig. 3.10, we show a comparison between ED and CT-HYB Green’s functions for two distinct hybridization functions, determined by the bath energies_σand hy-bridization parametersV_σ, which are chosen to be the same for both spin components. Together with these results, the performance of the CT-HYB algorithm can be discussed in the same way as the CT-AUX solver from Fig. 3.6, which shows the average perturbation order in dependence ofβtandβU. Since CT-HYB is an expansion in the hybridization, it does not scale withU but rather with the factorβt, which is

10In fact they are only identical after a rearrangement of the matrix M has been performed, since, in contrast to the segment insertion/removal, for the anti-segment insertion/removal the inserted/removed columns do not belong to the same segments.

Although this is a trivial modification of the matrixMafter the insertion/removal of an anti-segment it must not be forgotten or neglected since it totally changes the measurements.

Figure 3.10: Comparison of the interacting Green’s functiong(τ) =Gσ(τ)of an AIM with finite number of bath orbitalsS = 3obtained from a CT-HYB algorithm and the exact results from diagonalization of the Hamiltonian. The accuracy does not depend on the interaction parameterU, since CT-HYB is an expansion in the hybridization. The Green’s functions have been obtained for (left) bath energies σ = (1,−1,1)and hybridization parametersVσ = (1,−1,1)and (right) bath energiesσ = (3,2,1) and hybridization parametersVσ= (1,2,3). Although a very small number of sampling steps,N = 10⁶, was used, the agreement is already very good and can be improved by taking the usual number of steps N= 10⁸, which turned out to be optimal for most simulations that we performed.

the Hubbard hopping parameter and indirectly proportional to the hybridization.

In this chapter, we have found the basic formalism for the continuous-time Monte-Carlo impurity solvers CT-HYB and CT-AUX, for which we will discuss several extensions and improvements in the following chapter and which we will use later in combination with RDMFT to investigate the properties of fermionic lattice models with artificial gauge fields.

4. Extensions and Improvements to the CT-QMC Methods

The continuous-time Monte-Carlo impurity solvers introduced in the previous chapter, namely the CT-AUX and the CT-HYB solver, were so far only discussed for the case of a single impurity and with spin a good quantum number, i.e. without terms that mix the spins, as they would occur for systems which show superfluid ordering or spin-changing hopping parameters. In this chapter, we will discuss both, the extension of the CT-AUX and the CT-HYB solvers to systems with superfluid pairing and spin-changing hopping processes. We will also extend the CT-AUX solver to cluster impurities, containing more than just a single lattice site. Additionally, improvements of the measurements in the CT-HYB solver are discussed, which allow for a direct measurement of the self-energy [74] and which perform the measurements in a different basis compared to the previous chapter, i.e. the measurements will be performed in a polynomial basis instead of the discretized imaginary-time interval [16].

Im Dokument Topological phases of interacting fermions in optical lattices with artificial gauge fields (Seite 63-67)