stricted by a radius R$_{0}$ (c.f. Eq. (2.37)). In the rest of the sample the particles do not interact with each other and the field ˜ φ$_{r}$ ( τ ) vanishes there. The interacting part of the free energy functional should be proportional to the number of the sites N$_{d}$ located inside the sphere with the radius R$_{0}$ . In order to obtain a good precision one needs to consider a large number of sites in the bath, such that the total number of the sites N total d should consider- ably exceed N$_{d}$ . At the same time, considering the bath should not demand a large computation time because the HS field is identically zero in that region. We write the exact partition function Z of the original model in the form The operator M r,r ′ and the function ˜ Φ$_{r,r}$ ′ ( τ ) are defined as where a = 1 , 2 and where β = N Δ, φ$_{r}$ ( τ ) = φ$_{r,l}$ for ( l − 1) Δ ≤ τ < l Δ (see Fig. 6), and V$_{0}$ is the on-site repulsion. The values φ$_{r,l}$ satisfy the periodicity conditions φ$_{r,l}$ = φ$_{r,l}$$_{+}$$_{N}$ . The weight W$_{Δ}$ [ ˜ φ ] has the form Z = lim Δ → 0 ∫ Z$_{b}$ [ ˜ φ ] W$_{Δ}$ [ ˜ φ ] ∏ N l =1 dφ$_{r,l}$ ∫ W$_{Δ}$ [ ˜ φ ] ∏ N l =1 dφ$_{r,l}$ (4.44) = N Δ, ˜ φ$_{r}$ ( τ ) = φ$_{r,l}$ for ( l 1) Δ τ < l Δ W$_{Δ}$ [ ˜ φ ] = exp [ − Δ 2 V$_{0}$ ∑ r N ∑ l =1 φ 2 r,l ] . (4.45) M r,r ′ ( z ) = ˆ ε$_{r}$ − ˆ ε$_{r}$ ′ − uσ ˜ Φ$_{r,r}$ ′ ( τ ) , ˜ Φ$_{r,r}$ ′ ( τ ) = ˜ φ$_{r}$ ( τ ) − ˜ φ$_{r}$ ′ ( τ ) , (4.51) A$_{1;}$$_{r,r}$ ′ ( z ) = ( A ′ 1; r,r ′ ( z ) A ′′ 1; r,r ′ ( z ) ) , B$_{1}$ = ( 0 1 ) , A$_{2;}$$_{r,r}$ ′ ( z ) = ( A ′′ 2; r,r ′ ( z ) A ′ 2; r,r ′ ( z ) ) , B$_{2}$ = ( 1 0 ) . n$_{r,r}$ ′ = ∑ k n$_{k}$e ik ( r − r $^{$^{′}$)}$, (4.52) n$_{k}$ = [exp { βϵ$_{k}$ } + 1] − 1 . A$_{a}$$_{;}$$_{r,r}$ ′ ( τ, σ, u ) = A$_{a}$$_{;}$$_{r,r}$ ′ ( τ + β, σ, u ) (4.53) µ ′ = µ − V$_{0}$/ 2 (4.46) Z$_{b}$ [ ˜ φ ] = Z$_{0}$ exp [ − 1 2 ∑ r,σ,a ∫ β 0 ∫ 1 0 σ ˜ φ$_{r}$ ( τ ) A ′ a ; r,r ( z ) dudτ ] . (4.47) Z$_{0}$ = ∏ k,σ [ 1 + exp ( − βϵ$_{k}$ )] (4.48) ϵ$_{k}$ = − t ( k ) − µ + V$_{0}$/ 2 . (4.49) ( γ M r,r ′ + ∂ ∂τ M r,r ′ − ∂ ∂τ − γ ) A$_{a}$$_{;}$$_{r,r}$ ′ ( z ) = − un$_{r,r}$ ′ σ ˜ Φ$_{r,r}$ ′ ( τ ) B$_{a}$, (4.50) in the region under the radius R$_{0}$ . Outside this region the HS field is equal to zero identically. In order to write the correct functional Z$_{b}$ [ ˜ φ ] in the presence of the bath, the equations of Section II have to be adjusted as concerns the shift of the chemical potential by the amount − V$_{0}$/ 2, Eq. (2.14). We assume that main system is located in a potential well and that the height of the walls equals V$_{0}$/ 2. In other words, we shift the chemical potential µ everywhere in the space by V$_{0}$/ 2, although the electron-electron interaction is present only in the main system, using the value This homogeneous shift is convenient because at µ ′ = 0 the presence of the bath does not violate the particle-hole symmetry. Then, the functional Z$_{b}$ [ ˜ φ ] equals where and In Eq. (4.49), µ is the chemical potential of the original fermion model, Eq. (2.2). The functions A ′ a ; r,r ( z ) should be found from the linear equations These arguments are based on the fact that the solu- tions A$_{r,r}$ ′ ( z ) are real and non-singular, such that inte- grating over u does not generate imaginary parts. This property of the solution follows from the spectral expan- sion, Eqs. (2.52, 2.54) supplemented by the proof that the real eigenvalues E entering the sum in Eq. (2.52) in pairs E and − E cannot be equal to zero (see Appendix A). Moreover, due to the fact that the eigenstates en- ter in pairs with the eigenvalues E and − E, the entire exponent in Eqs. (2.54, 4.47) does not become singular even for very small E due to the mutual compensation of the divergencies in the pairs. This follows also from the symmetry property of the Green function, Eq. (2.53). As all the parameters of the equations are real, there is no danger to obtain an imaginary part in the exponent in Eq. (2.54) and therefore the weight Z$_{b}$ [ ˜ φ ] is always positive. In order to calculate the partition function Z , one should solve Eq. (4.50-4.53) for a set of functions φ$_{r,l}$ defined on the slices and a small parameter γ , substi- tute the solution A ′ a ; r,r ′ ( z$_{l}$ ) into Eq. (4.47) and find the function Z$_{b}$ [ ˜ φ ]. This function is positive because the so- lution A$_{a}$$_{;}$$_{r,r}$ ′ ( z$_{l}$ ) must be real. Then, one can calculate Z in Eq. (4.44) using Z$_{b}$ [ ˜ φ ] as the probability in the MC sampling for finite Δ diminishing this length as much as necessary. We emphasize here that there should not be any principle limitation of the computational accuracy provided a sufficiently large system is considered and suf- ficiently small length of the time slices is used. The real parameter γ should be taken as small as neces- sary to be sure that the final result is insensitive to its value. Both the functions ˜ φ$_{r}$ ( τ ) and A$_{r,r}$ ′ ( z ), z = ( τ, σ, u ), should be periodic in time: ˜ φ$_{r}$ ( τ ) = ˜ φ$_{r}$ ( τ + β ) and Finally, n$_{r,r}$ ′ is the Fermi distribution 23