{'fraction_non_alphanumeric': 0.03066452550279116, 'fraction_numerical': 0.011011699931176875, 'mean_word_length': 4.293462861768873, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantum Monte Carlo Study of Disordered Fermions

We study a strongly correlated fermionic model with attractive interactions in the presence of disorder in two spatial dimensions. Our model has been designed so that it can be solved using the recently discovered meron-cluster approach. Although the model is unconventional it has the same symmetries of the Hubbard model. Since the naive algorithm is inefficient, we develop a new algorithm by combining the meron-cluster technique with the directed-loop update. This combination allows us to compute the pair susceptibility and the winding number susceptibility accurately. We find that the s-wave superconductivity, present in the clean model, does not disappear until the disorder reaches a temperature dependent critical strength. The critical behavior as a function of disorder close to the phase transition belongs to the Berezinsky-Kosterlitz-Thouless universality class as expected. The fermionic degrees of freedom, although present, do not appear to play an important role near the phase transition.


I. INTRODUCTION
A variety of systems show quantum coherence over large distances at low temperatures. Superfluidity and superconductivity are two striking physical phenomena showing such behavior, which have been extensively studied over the years. However, when correlations between the microscopic degrees of freedom become strong it is difficult to study these phenomena theoretically from first principles. The calculations must take into account strong fluctuations over many length scales which is only possible numerically. When the microscopic degrees of freedom involve bosonic variables one can usually devise efficient quantum Monte Carlo methods to solve the problem. 1 On the other hand, it is still difficult to study a variety of models from first principles when the microscopic theory is fermionic. For example, the critical temperature below which superconductivity is seen in the well-known attractive Hubbard model was only determined recently, 2 using the determinantal Monte Carlo method 3, 4 and on lattices only as large as 18 × 18. The main approaches to dealing with fermionic systems can be viewed as arguments that universality allows one to replace the microscopic theory with an effective low energy theory. The resulting effective theory is usually either a Fermi-liquid theory, a BCS-type mean-field theory, or some bosonic theory. 5 A key element in furthering microscopic understanding, then, is to validate the universality arguments and determine the low-energy effective theory; in practice, this has proved very difficult for systems with strong correlations.
Real systems usually contain impurities. Thus, in addition to understanding superconductivity in clean systems, the effects of impurities in the form of disorder need to be incorporated in the studies. In certain systems like two-dimensional superconducting films and Josephson-junction arrays, it has been discovered that superconductivity can be destroyed by tuning parameters such as the film thickness. 6,7 Since these tuning parameters change the effective strength of the disorder, it is believed that the superconductor-to-insulator (SI) phase transitions in these systems can be understood as being driven by disorder. Among the models used to explain the experiments, the attractive Hubbard model with disordered chemical potential is one well-known starting point. 8,9 The relevance of disorder for superconductivity was first addressed by Anderson, 10 where he argued that superconductivity is insensitive to perturbations that do not destroy time reversal invariance. Using a BCS type trial wave function Ma and Lee 11 showed that superconductivity can persist even below the mobility edge. Clearly, these studies suggest that an SI transition is an effect of strong disorder which makes it a difficult subject for analytic work. Fisher et al. 12 have argued that the effective theory describing the transition is bosonic, and then developed a deeper understanding of the purely bosonic superfluid-insulator transition using renormalization group arguments along with scaling. A variety of quantum Monte Carlo work has been done over the years on these purely bosonic microscopic theories. 13,14,15,16,17,18 If fermions do not play an important role near the transition, it is likely that these studies will also be useful in understanding the universality of the fermionic SI transition. Recently determinantal quantum Monte Carlo studies of the attractive (fermionic) Hubbard model with disorder have been performed, 9 which show that it is indeed possible to drive an SI transition by increasing the disorder and, as expected, the critical disorder is large. However, the system sizes explored were quite small, 8 × 8. Other studies of disorder effects also involved only very small systems. 19,20 Motivated by the physics of the SI transition, in this article we study the effects of disorder in a strongly correlated fermionic model. Our model is unconventional and has been built so that it can be studied using the recently discovered meron cluster algorithms. 21,22 These novel algorithms are so efficient that lattices as large as L = 128 were studied recently, and it was shown with great precision both that the low temperature phase of the clean model is indeed superconducting and that the finite temperature phase transition belongs to the Berezinski-Kosterlitz-Thouless (BKT) universality class. 23,24 Here we use the same model to explore the effects of disorder on superconductivity and, in particular, focus on the role of fermions. Unfortunately, the naive extension of the earlier algorithm becomes inefficient in the presence of disorder; hence, we also develop a new algorithm by combining the meron-cluster formulation with the directed loop algorithm. 25 This new algorithm allows us to measure the relevant observables very accurately.
Our paper is organized as follows: In Section II, we introduce our model and define the observables that we use later. In Section III, we rewrite the model in a cluster representation. Section IV explains the new directed-loop algorithm we have developed. Section V contains our results, and Section VI contains our conclusions and directions for the future.

II. THE MODEL
The model we consider in this article was motivated by the ability to solve the fermion sign problem using the the meron-cluster algorithm. 22 The Hamiltonian of the model can be written as where H (2) ij consists of all the nearest neighbor interactions between sites i and j on an L × L square lattice and H where U represents the Hubbard interaction between spin-up and spin-down electrons and µ i is the local chemical potential, through which disorder is introduced in the model. The term H (2) ij is unconventional and is given by Here c † iσ and c iσ are the usual creation and destruction operators of spin σ at site i, n = c † c, and n ij = σ n iσ + n jσ . S i is the spin operator on site i defined by and J i is the pseudo-spin operator defined by J z i = 1 2 (n i↑ + n i↓ − 1).
J + and J − are related to pair creation and annihilation operators. In our notation i x(y) refers to the x(y) component of the site i. Although the Hamiltonian we study is unconventional, it has all the relevant symmetries of the Hubbard model when J 3 = 1. In particular when µ = 0 the Hamiltonian is invariant under the SU (2) spin and SU (2) pseudo-spin transformations. When µ = 0, the pseudo-spin symmetry is broken to the U (1) fermion number symmetry. One can introduce repulsion or attraction by making U sufficiently positive or negative respectively. The important difference with the Hubbard model is that when U = 0 the model is still strongly interacting and by setting J 3 = 1 we can break the pseudo-spin symmetry. Further, the model simplifies in the U → −∞ limit; in this limit the model can be mapped to the simple Hamiltonian Clearly, when J 3 = 1, one obtains the antiferromagnetic Heisenberg model involving pseudo-spins, while J 3 = 0 leads to the XY model. When J 3 = 1, U → ∞ and µ i = 0 one gets the antiferromagnetic spin model An interesting aspect of this model is that the fermion sign problem can be solved using the meron cluster approach for 0 ≤ J 3 ≤ 1 when U < 0 at any value of µ i . Also when J 3 = 1 the sign problem can be solved when U > 0 when µ < U/2. Thus, we think the model offers a rich phase diagram and deserves to be investigated. In this paper we will consider −∞ < U ≤ 0 and investigate the physics when J 3 = 0 and J 3 = 1. We introduce disorder through µ i = µ + δµ i where δµ i is a random number distributed uniformly from −∆ to ∆.
In order to probe superconductivity in this system we will focus on two observables. The first is the pair susceptibility defined as where V is the spatial volume, T is the temperature, and . The second is the winding number susceptibility defined as is the current-current correlation function where J µ (i) is the conserved fermion current, at the site i in theμ =x,ŷ direction. In order to estimate the importance of fermions we will also look at the density of singly occupied sites defined as and at the total density of electrons defined by A comparison of n s and n will tell us how many sites have formed local pairs.

III. CLUSTER REPRESENTATION
It is possible to rewrite the partition function of our model in discrete time in terms of a statistical mechanics of closed loops on a space-time lattice. 22 We first divide β, the length in the Euclidean time direction, into M equal steps such that ǫ = β/M . Interactions between nearest neighbor sites are introduced in a checker-board type manner, so that on each time slice every site interacts with a unique neighbor. This then introduces 4 extra time slices for every ǫ time step. In the cluster representation the nearest neighbor interactions occur in the form of three types of bond configurations on spacetime plaquettes as shown in Fig. 1. Their weights are given by, Given a configuration {C} of bonds, one can connect them together to form many closed loops; we denote them by C α , α = 1, 2, ..., N c . The partition function can then be written as where ω P is the weight of the bond associated with the plaquette P and takes one of the values given in Eq. (15).
is the weight corresponding to each loop C α that arises due to the fermionic degrees of freedom associated with the loop. Here where the sum is over all all space-time points that belong to the cluster C α . Thus S Cα is just the size of the cluster. On the other hand where w iτ is +1 when the cluster is going forward in time and −1 when it is going backward in time at the site iτ . If the cluster moves horizontally, our convention is that the temporal direction is reversed. In order to determine w iτ one can start from any point and traverse the cluster in either direction. Note that is the temporal winding of the cluster C α . Finally, the factor σ(C α ) in Eq. (17) is a sign factor associated with the cluster topology, that arises due to the fermion permutation signs. 22 Following Ref. 22 we call the cluster a meron if σ(C α ) = −1. If N h (C α ) is the number of horizontal hops in the cluster C α , then the cluster is a meron if and only if N h (C α )/2 + W t (C α ) is even. An example of a bond configuration in two dimensions is shown in Fig. 2.
In the cluster representation it is easy to show that where δ (iτ ;jτ ′ )∈Cα imposes the restriction that both the spacetime sites (iτ ) and (jτ ′ ) belong to the cluster C α and where s iτ,jτ ′ kτ ′′ is +1 while going from (iτ ) to (jτ ′ ) and −1 while continuing from (jτ ′ ) to (iτ ). The winding number susceptibility is given by where W µα refers to the spatial winding of the loop clusters along the spatial direction µ. The density of single occupation turns out to be and is a measure of the number of unpaired fermions.

IV. DIRECTED-LOOP ALGORITHM
A simple Monte Carlo algorithm for the current problem involves visiting every interaction plaquette and updating the bond configuration on it by replacing it with one of the three choices shown in Fig. 1. Since the Boltzmann weight also depends on the structure of the loops formed by these bonds, the decision involves figuring out the connectivity of the sites of the plaquette (referred to as P, Q, R, S) due to the bonds on other plaquettes. This connectivity can be one of three types as shown in Fig. 3. Thus, choosing a new interaction involves finding weights of nine configurations (three bond configurations for three types of connectivity) and the new bond configuration can be found by using a heat bath or a Metropolis step.
We have found that in the presence of a chemical potential µ, this simple algorithm is inefficient. This behavior can be understood by noting that a chemical potential is similar to a magnetic field in a a quantum spin model and in the context of quantum spin models there is evidence that this type of naive algorithm becomes inefficient in the presence of magnetic fields. 26,27 Today it is well known that quantum spin systems in the presence of a magnetic field can be solved efficiently using the directed-loop algorithm. 25 However, until now this algorithm has been formulated only in the spin representation and not in the cluster representation. Unfortunately, the sign problem in the fermionic model can only be solved in the cluster representation. In this article we show how one can extend the directed-loop algorithm to the cluster representation which then leads to an efficient algorithm for the fermionic model even in the presence of a chemical potential.
The basic idea behind the directed loop algorithm is to extend the configuration space so that configurations that contribute to certain two point correlation functions (denoted {C (2) iτ,jτ ′ }) are sampled along with the configurations that contribute to the partition function (denoted {C}). The configurations C (2) iτ,jτ ′ have two reference space-time points, iτ and jτ ′ ; during the directed-loop update one of these points, say iτ , is held fixed while the other point jτ ′ is moved around. The directed loop update begins with a configuration in the set {C} and chooses a site iτ at random and probabilistically creates a configuration in the set {C iτ,jτ ′ } is created, the point jτ ′ is moved around while satisfying detailed balance and thus sampling other configurations in the set C (2) iτ,jτ ′ with the correct Boltzmann weight. Finally, when the two points meet again, i.e.., when jτ ′ = iτ , the two points may be removed to obtain a configuration in the set {C} in accord with detailed balance. Thus, since at every step detailed balance is satisfied, it is easy to show that the directed-loop update, which starts from a configuration in {C} and ends on another configuration in {C}, satisfies detailed balance. During the loop update, all the sites encountered contribute to the two point correlation function.
In the current work we have used the pair correlations to develop the directed loop algorithm. Thus, the configurations in {C (2) iτ,jτ ′ } are the ones that contribute to the pair correlations in Eq. (21). Thus, the weight of such a configuration is taken to be where the sites iτ and jτ ′ are forced to remain on the same cluster α 0 . The weight of a configuration in the set {C} is The essential steps of the directed loop algorithm are as follows: (i) Start with the initial configuration which belongs to the set {C}.
(ii) Select a space-time site iτ at random and propose to create a configuration C iτ,jτ ′ assuming jτ ′ = iτ . Ac- cept the proposal with probability where C α0 is the cluster which contains the site iτ .
(iii) If the proposal is not accepted then the update is complete. Otherwise we go on.
(iv) The site jτ ′ is moved to the next site by picking the plaquette it is connected to which is not the one just visited. Since each site is connected to two plaquettes the plaquette is unique unless one is at the beginning. In that case one chooses one of the two randomly.
(v) Each plaquette update involves choosing one of ten possible configurations. These configurations depend on the connectivity of the plaquette and an example is given in Fig. 4. One of these ten is chosen using a heat bath.
(vi) Steps (iv)-(v) are repeated until jτ ′ reaches iτ , at which stage the transition to the {C} sector is made with probability Min If the transition is made then the directed-loop algorithm ends, otherwise one goes back to step (iv) assuming one is in at the beginning of the loop.
In the above algorithm, the pair susceptibility can be computed using the formula where S is the number of sites visited during the directed loop update. Other observables such as χ w and n s can be computed using the formula of Eq. (23) and (24). We have tested the efficiency of the directed loop algorithm by comparing the results with spin model results in the limit U = −∞ which can be obtained using the usual directed loop algorithm. Since each update of the plaquette in the fermion algorithm requires knowing the connectivity of the loop, the algorithm is indeed much slower than the directed loop algorithm of a spin model which does not require this step. Unfortunately, this is a price one has to pay for being able to compute quantities in a fermionic theory. In Ref. 23 a trick was used to reduce the time to determine the connectivity of the clusters. The trick was to use a "tree" structure to store the information about the cluster connectivity which allowed one to obtain the relevant information in a time that grew like the logarithm of the cluster size. This was not implemented in the current work but could be implemented if necessary.

V. RESULTS
In this Section we discuss the results obtained from extensive simulations for lattice sizes up to L = 32. In our work we have fixed ǫ = 0.25 in order to avoid changes in the time discretization errors. We have found that this value of ǫ is reasonably small and the results at smaller values appear to join within our error bars. Further, since our desire is to understand universal physics of disorder, we believe that fixing ǫ should not be a major concern since it only changes the transfer matrix by a small amount. For a given disorder realization we typically discard the first 1000 directed loop updates for equilibration and then average over 20000 directed loop updates in blocks of 1000 to generate each of our statistics. All quantities plotted at a given value of ∆ have been averaged over 20 disorder realizations.

A. Superconductivity with Disorder
It is known from earlier studies 23 that our model has a low temperature superconducting phase when U = 0, J 3 = 0 and µ = ∆ = 0. The superconductivity disappears at a finite temperature and the phase transition belongs to the expected BKT universality class. Here we study the effects of disorder on this system by keeping µ = 0 but ∆ = 0. In the clean model the BKT predictions for the leading finite size scaling form of the pair susceptibility and the winding number susceptibility are known: and where A, B and L 0 are constants which depend on the temperature. We further expect 0 ≤ η ≤ 0. 25  B = 2 at the phase transition. 23,28 Here, instead of varying the temperature we fix the temperature at T = 0.25 and study the effects of disorder by increasing its strength through the parameter ∆. If the BKT universality holds we then expect the same finite size behavior to be true where the constants A, B and L 0 now depend on ∆.
In Fig. 5 we plot our results for χ p and χ w as a function of L. The fits are shown in Table I. The data appear to be consistent with a BKT transition around a critical disorder of ∼ 1.5. Note in particular the excellent power law behavior of χ p for ∆ ≤ 1.55. The expected form for χ w is seen for ∆ ≤ 1.35, but it does not fit very well for ∆ ≈ 1.5. From the rapid decay of χ w (confirmed by saturation in χ p , not shown), I: Fitting results for the pair susceptibility (χp) and winding number susceptibility (χw) for J3 = 0, U = 0, and µ = 0 (Fig. 5). The fitting formulas are Eqs. (30)-(31); the χ 2 given is per degree of freedom from the fit. the system is clearly no longer superconducting when ∆ = 3.
To obtain a good estimate of the critical disorder, ∆ c , we assume that the forms (30)-(31) hold close to the transition and that the deviation of the exponent η and constant B is linear in ∆: A joint fit of η and B to this form for ∆ in the interval [1.35, 1.65] yields ∆ c = 1.53 (4). When U = −∞ our model reduces to the XY model in the pseudo-spin variables. As seen in Fig. 5, ∆ = 1.65 is still in the superconducting phase in the XY limit. 29 Thus, the effect of fermions is to disorder the superconductor more quickly, as can be intuitively expected because of the increased entropy.
As discussed earlier, when J 3 = 1 and µ = 0 the model has an additional SU (2) pseudo-spin symmetry, as in the attractive Hubbard model. Thus, due to the Mermin-Wagner theorem, superconductivity is only possible when µ = 0. The J 3 = 1 model also has been studied earlier in the absence of disorder, 24 and a BKT transition was established using universal finite size scaling. We have extended these calculations to the disordered regime. We again fix the temperature at T = 0.25 and study the effects of disorder with µ = 1. Fig. 6 gives our results and Table II shows the corresponding fits. The expected BKT transition behavior is indeed seen for ∆ ≤ 0.5. Note that the values obtained for η are constant, within statistical error, in this range. Thus to extract the critical disorder ∆ c we use only the χ w data. A fit of the ∆ ≤ 0.5 values to the form in Eq. (33) yields ∆ c = 0.13 (5); the data for η is consistent with this value. The critical disorder found here is much smaller than in the J 3 = 0 case above, indicating as expected that superconductivity is weaker when J 3 = 1.

B. Fermionic degrees of freedom
Since we are studying a strongly correlated electronic system with on-site attraction between spin-up and spin-down electrons, one might worry that the electrons have formed local pairs on each site and the model is essentially bosonic. Hence, it is important to look at observables that extract fermionic information in the model. One such observable is the density of singly occupied sites defined in Eq. (13). In Fig. 7 we plot n s as a function of temperature when J 3 = 0 and µ = 0. Since there is particle-hole symmetry at µ = 0 we have n = 1. For this calculation we have fixed ǫ = β/M = 0.25 for β ≥ 4 and M = 64 for β < 4. We show two different values of ∆ -one for weak disorder and the other for strong disorder -at different values of U . The data shown was obtained at L = 16; we have observed that the density does not vary much as the lattice size increases. First we note that the density of singly occupied sites decreases significantly as the temperature is lowered for U = 0. When the disorder is small the pairs begin to break close to the II: Fitting results for the pair susceptibility (χp) and the winding number susceptibility (χw) for J3 = 1, U = −0.1, and µ = 1 (Fig. 6). The fitting formulas are Eqs. (30)-(31); the χ 2 given is per degree of freedom from the fit.  figure). This leads us to conclude that the properties of the phase transition are most likely governed by a bosonic model such as the disordered quantum XY model, especially in the strongly disordered case. Thus, our work gives credence to the "dirty boson" scenario. The dotted line indicates the approximate critical value of β = 1/T where the system undergoes a BKT transition to a superconductor. The total particle density is n = 1 due to particle-hole symmetry.

VI. DIRECTIONS FOR THE FUTURE
In this article we have studied the effects of disorder on a strongly correlated electronic system. Our model was known to be superconducting in the clean limit and our motivation in this work was to study the effects of disorder through a position dependent chemical potential. Unfortunately, a naive extension of the meron cluster algorithm was found to be inefficient in the presence of a disordered chemical potential. Hence in this work we constructed a new and efficient algorithm by combining the the meron-cluster approach with the directed-loop algorithm. Earlier work on the directedloop algorithm involved quantum spin systems and always was constructed in the spin representation. In this work we have shown that the algorithm can be constructed even in the cluster representation, which is essential in the fermionic system in order to solve the fermion sign problem. Our new algorithm was quite successful and allowed us to compute the pair and winding number susceptibilities accurately.
We found that disorder significantly suppresses superconductivity and the system undergoes a phase transition which appears consistent with the BKT universality class. Although this scenario has been expected, 15 our work is the first, as far as we know, to study the universal scaling predictions of the BKT transition in a fermionic system with disorder. We could go to lattices as large as L = 32 thanks to our new algorithm.
We found that when J 3 = 0 superconductivity is stronger than when J 3 = 1.
We also found that there indeed are fermionic excitations in the system, but they are not affected by the disorder. The role of these fermions remains an interesting open question. For example do the background fermions form a Fermi-liquid in the weak disorder regime? If this is the case, then it would be interesting to ask whether the fermions become localized or do they remain extended. Is the phase transition between a superconductor and an insulator or whether it is a transition between a superconductor and a metal. Finally, although we have focused on attractive interactions in this work, we can study the repulsive model by setting U positive. In that case it is possible to add a chemical potential such that µ ≤ U/2 when J 3 = 1 without introducing a sign problem. 22 These studies have the potential to increase the fermionic effects further.