Taming the pion condensation in QCD at finite baryon density

In the Monte Carlo study of QCD at finite baryon density based upon the phase reweighting method, the pion condensation in the phase-quenched theory and associated zero-mode prevent us to go to the low-temperature high-density region. We propose a method to circumvent them by a simple modification of the density of state method. We first argue that the standard version of the density of state method, which is invented to solve the overlapping problem, is effective only for a certain `good' class of observables. We then modify it so as to solve the overlap problem for `bad' observables as well. While, in the standard version of the density of state method, we usually constrain an observable we are interested in, we fix a different observable in our new method which has a sharp peak at some particular value characterizing the correct vacuum of the target theory. In the finite-density QCD, such an observable is the pion condensate. The average phase becomes vanishingly small as the value of the pion condensate becomes large, hence it is enough to consider configurations with small values of pion condensate, where the zero mode does not appear. We demonstrate an effectiveness of our method by using a toy model (the chiral random matrix theory) which captures the properties of finite-density QCD qualitatively. We also argue how to apply our method to other theories including finite-density QCD. Although the example we study numerically is based on the phase reweighting method, the same idea can be applied to more general reweighting methods and we show how this idea can be applied to find a possible QCD critical point.


Introduction: importance sampling and un-importance sampling
The sign problem is a severe obstacle for Monte Carlo methods based on the importance sampling, and it prevents us, for example, from studying lattice QCD at finite baryon density directly by Monte Carlo simulations, since the fermion determinant becomes complex at a finite baryon chemical potential. (For an introductory review from lattice perspective, see [1]. A review from the point of view of nuclear theory can be found in [2].) Several methods have been proposed to overcome this difficulty (for various previous attempts, see e.g. [3,4,5,6,7,8,9]), and some of them are based on the phase-reweighting technique, which, however, fail to work at high density due to the unphysical pion condensation. In this paper we propose a new method to tame the pion condensation problem of reweighting methods, whose basic idea can also be applied to some classes of sign problems. As a bonus, a zero mode associated with the unphysical pion condensation is eliminated.
Let us begin with identifying the physical origin of the sign problem. We consider a field theory on Euclidean spacetime with a complex action, Then the path-integral weight e −S is not real and positive anymore, and hence the importance sampling cannot be applied as it is. Therefore one performs the importance sampling by using a real and positive weight which 'approximates' the complex weight and take into account the effect of the non-positivity by using so-called reweighting methods. The simplest example is the phase-reweighting method, in which the phase-quenched weight e −S R is adopted; the expectation value of an operatorÔ in the full theory is obtained by using an identity Ô full = e iS I ·Ô P.Q. e iS I P.Q. , where · full and · P.Q. stand for expectation values in the full and the phase-quenched theories, respectively 4 . Then the right hand side is calculable in principle. In practice, however, both e iS I P.Q. and e iS I ·Ô P.Q. can become extremely small in some cases and then the right hand side is essentially 0/0, which is not easy to evaluate numerically. This is the sign problem. The sign problem becomes even severer when the vacua of the full and phase-quenched theories are different; this is so-called 'overlap problem'. In order to understand it, let us consider a certain observableÔ which characterizes the vacua of these two theories; the vacua are characterized by Ô full = K full and Ô P.Q. = K P.Q. , where K full = K P.Q. in general. Let the full theory is proportional to ρ P.Q. (x) · e iS I x , where e iS I x is the average phase factor with the value ofÔ fixed to x. Since ρ full (x) ∼ ρ P.Q. (x) · e iS I x peaks around K full , the phase factor e iS I x ∼ ρ full (x)/ρ P.Q. (x) is vanishingly small around x = K P.Q. = K full . (This point is clearly demonstrated in [10] by using a solvable model.) This means that, although most configurations sampled in the phase quenched simulation are around x = K P.Q. , their contribution vanishes due to huge sign fluctuation, and the true peak of the full theory appears from the tail of ρ P.Q. (x). In other words the phase-quenched simulation is the un-importance sampling, in the sense that the most of computational resources are wasted to sample un-important configurations. In fact it is even worse -the sign fluctuation becomes violent in order to erase un-important configurations, and unless one has huge amount of configurations so that vanishingly small value of the phase factor can be measured precisely, the error bar becomes large; essentially the only contribution of the un-important samples is to make the error bar larger. Such a waste of computational resources, which arises due to the lack of the overlap between vacua in full and phase-quenched theories, is the overlap problem.
In terms of the above general argument, we consider the massless two-flavor QCD with the finite baryon chemical potential (QCD B ). In this theory, two quarks (up and down) has the same value of the chemical potential µ, which coupled to the baryon number of quarks, +1/3. The partition function in Euclidean space-time is given by where D µ is the gauge covariant derivative acting on quark fields ψ = u, d, D µ (A)ψ = (∂ µ − iA µ )ψ with the gauge field A µ , and S G (A) is the action for the gauge field. The determinant factor satisfies and hence it is complex at µ = 0, so that the sign problem exists in QCD B . The phase quenched theory is described by the partition function, This theory is QCD with a finite isospin chemical potential (QCD I ), in which up and down quarks have chemical potential +µ and −µ, respectively. Hence this chemical potential couples to the isospin number, +1/2 for up and −1/2 for down. In the full theory, nothing happens until the nucleon, whose mass is about 1 GeV, condenses. On the other hand, in the phase quenched theory, the massless charged pion π + =dγ 5 u condenses as soon as µ is turned on. Therefore the overlapping problem arises due to the pion condensation.
In this paper we propose a simple way to tame the sign problem caused by the overlap problems associated with the pion condensation in the phase quenched theory. We first notice that, if one eliminates the pion condensate by hand (for example by adding delta-function like potential), two theories, QCD B and QCD I , become equivalent when the number of colors N c is sent to infinity [11,12] (N c = 3 is the usual QCD), which means that the overlap problem is just a 1/N c effect if we fix the pion condensate 5,6 . This consideration leads to our main idea that the overlap problem can be avoided by pinning down an appropriate observable, which characterizes the difference between full and phase quenched theories (in the case of the finite density QCD, the pion condensate), to the right value (zero pion condensation in QCD B ), and the sign fluctuation becomes milder there. Away from the correct vacuum, the sign fluctuation becomes severer. This is not drawback anymore, because the severe sign fluctuation is simply telling us that such configurations are not important. When the sign fluctuation becomes severer, we do not have to measure the average sign. Rather, we can safely omit such configurations. The sign fluctuation is not a problem anymore, rather it reduces numerical costs of our simulations. Furthermore, this methods automatically avoids a zero mode associated with the unphysical pion condensation, since we do not have to consider the large-π + region, where the zero mode appears.
Our method is a natural generalization of the density of state method. In order to illustrate the advantage of our method, we first review the traditional density of state method, explain what is good and what is insufficient, and then introduce our method. (Our method could be regarded as a simplified version of the multi-parameter factorization method [19], which has been applied for a supersymmetric matrix model 7 . We also explain how our method can be combined with the multi-parameter factorization.) In this paper, we demonstrate our idea taking the chiral random matrix theory (RMT) as a simpler example. Because RMT is analytically solvable and computationally much cheaper than QCD, we can test the method thoroughly. We explain basic ideas in Sec. 2 using the chiral RMT. Note that our main idea does not rely on the detail of the theory and the method can be generalized to QCD and other theories. In Sec. 3 we give simulation results of the chiral RMT to show how our method works. In Sec. 4 we briefly discuss strategies for the finite-density QCD using our method, and give more generic reweighing method in Sec. 5. Our conclusion and discussion are given in Sec.6.

Methodology
In this section we explain our method using the chiral RMT as a concrete example.

β = 2 RMT
The action of the β = 2 RMT [13,14] with chemical potential [15] is given by where 5 The equivalence at Nc = ∞ holds if one takes the massless limit after taking the large-Nc. Strictly speaking, at very large µ, other isospin-charged particles like the rho-meson would condense and lead to the overlap problem, and then their condensates must be fixed to be zero. 6 The remaining overlap problem is due to the gas of pions. The overlap problem is mild as long as the pion does not condensate, and even the phase quench is exact at large-Nc. We will comment on this point later. 7 We would like to thank J. Nishimura for a comment on this point. and Here Φ is N × N complex matrix. From now on we take the number of flavors N f to be two (up and down quarks). We assign µ 1 = µ 2 = µ for the full theory (finite baryon chemical potential) and µ 1 = +µ, µ 2 = −µ for the phase-quenched theory (isospin chemical potential). We call these matrix models as RMT B and RMT I , respectively. Hereafter we will take a massless limit (m u = m d = 0). Therefore 'chiral condensate" ūu +dd will not be discussed in this paper. The pion condensate is identically zero unless we introduce a source term. We introduce a source term to RMT I as where c is a real number, and Then the 'pion condensate' is real and satisfy As an observable we will measure the baryon density ν B , which is defined by In QCD, we havē where the charge conjugations are defined by C is the charge conjugation matrix satisfying C −1 γ µ C = −(γ µ ) T , and T stands for the transpose. This also implies Using these we see that which means ν B in the phase quenched theory can be regarded as the isospin density ν I in QCD I . It is easy to see explicitly that these properties also hold in the RMT.

Standard density of state method: when it works and when it fails
First let us explain the standard density of state method (in the context of the finite-density QCD, see e.g. [16,17,18]), in order to illustrate the essence of our method explained in the next subsection.
Suppose we want to measure a certain quantityÔ, for exampleÔ = ν B . In the density of state method, one first classifies configurations obtained from the phase-quenched simulation in terms of values ofÔ. Let the number of configurations (or equivalently, the height of the histogram) at i , and the average sign be e iS I (Ô) i . (Here we have assumed theÔ takes only real values for simplicity.) Then we have a trivial relation, In the same manner, Therefore the phase reweighting can be done as In a naive phase-quenched simulation, the configurations are generated with the weight ρ It is commonly believed that the density of state method solves the overlap problem completely, because all the values ofÔ are scanned. In fact this is not really true, because this method is not based on the idea of the importance sampling. As we have explained in the introduction, the role of the sign fluctuation is to erase contributions from the wrong vacuum (the vacuum of the phase quenched theory) and to realize the correct vacuum of the full theory. Hence the sign fluctuation should be mild around the true vacuum. Therefore, if the correct vacuum can be characterized by tuning the value ofÔ (e.g.Ô is the pion condensate in QCD B ), e iS I (Ô) i and e iS I ·Ô (Ô) i can have a sharp peak at a particular value of i. On the other hand, if the correct vacuum cannot be specified by simply tuningÔ, their distributions do not show a peak structure and hence the remaining sign problem is not under control.
In summary, the standard density of state method is effective when the quantity of interest characterizes the vacuum of the full theory. It is unlikely, however, that a quantity one takes without considering properties of the full theory correctly characterizes its vacuum. Therefore there exists a danger that one has to spend a lot of computational resources to determine a small value of the average for the remaining sign precisely. It has been sometimes reported that the remaining sign problem can spoil the density of state method, due to this inappropriate choice of observables fixed [19].

Our method
As we have seen in the previous subsection, the traditional density of state method is effective only when the quantity of interest characterizes the vacuum of the full theory, since otherwise the overlap problem still exists. We solve this problem by slightly changing the viewpoint; we do not fix the quantity we want to measure. We fix an observable which characterizes the correct vacuum, called a good observable (If more than one quantities are needed to be specified in order to characterize the vacuum, we must fix all of them.) In the case of the finite-density QCD, pion condensate is such a good observable, since the phase quenched theory becomes exact at large N c as long as the pion condensate is forbidden by hand [11,12]. Note that we have to understand the physics of the full and phase quenched theories in order to find an appropriate observable which characterizes the correct vacuum.
Let us explain the detail of our idea by using the finite-density QCD. We first classify the configurations in the phase-quenched simulation by the values of the pion condensate π + . (Here π + is defined in terms of QCD I , i.e. the chemical potentials for up and down quarks in the operator are +µ and −µ, respectively, rather than +µ and +µ.) Let the height of the histogram at x i < π + < x i+1 be ρ i . We also calculate the average sign at x i < π + < x i+1 , which we denote e iS I i . Then we have a trivial relation, where ρ i is the relative weight factor of x i < π + < x i+1 in the phase quenched simulation 9 . In the same manner, whereÔ is an arbitrary operator we are interested in other than π + . Therefore the phase reweighting can be done as We can expect that i e iS I i · ρ i and i e iS I ·Ô i · ρ i takes non-negligible values only around the vacuum of the full theory, π + = 0. Therefore we only have to study there; when i e iS I i · ρ i and i e iS I ·Ô i · ρ i become so small that the precise determination is difficult, they do not affect the results and hence we can simply omit them. (In fact, unless we have extremely large statistics with which we can determine the small average phase at nonzero π + , adding such configurations just increases the error. Therefore, by throwing away un-important configurations we can make the result more precise.) Note again that this method is not purely numerical; we know the important samples based on physics. The huge sign fluctuation then tells us that we do not have to measure them, so that it does not increase the simulation cost. Instead it reduces the cost. Therefore the sign problem turns into the sign blessing in this situation.
The actual simulation for RMT goes as follows 10 . We add a deformation term for |π + − x| ≥ . The constraint parameter γ is taken sufficiently large so that all samples lie in |π + − x| < during the simulation. Note that there are two options: • Introduce the source both for S and ∆S. In this case we have to make the zero source extrapolation in the end.
• Introduce the source only for ∆S. In this case we do not need to take the zero source extrapolation. We take this option in this paper.
It is important to stress that π + ≥ 0 as long as c > 0, so that π + 0 implies that only π + 0 configurations contribute in the full theory.
Firstly we have to determine the distribution of the histogram of the pion condensate in the phase quenched simulation precisely. By introducing the deformation and tuning x and we can sample the tail effectively. The histograms obtained are 'partial' ones restricted at [x− , x+ ]. This situation is like the leftmost panel in Fig. 2; here the simulation has been done for x = x 0 (left), x = x 0 + (center) and x = x 0 + 2 , with the common value of , the number of configurations in the partial histograms are A i for π + < x and B i for π + > x. In order to obtain the full histogram, we rescale them as in the second panel, and then glue them as in the third panel. Note that, if the difference of two theories are characterized by many observables, we must fix all of them. In the case of finite density QCD, for example, as the isospin chemical potential becomes larger, not just pion but also other fields such as the ρ-meson can condense. Then we should add deformation terms to fix them.
Our method could be regarded as an improved version of the multi-parameter factorization method [19], which has been applied for a supersymmetric matrix model. (The 'factorization method' is essentially the same as the density of state method.) For this improvement, a good understanding about the physics under consideration is crucial. In the multiple-parameter factorization method, one labels the configurations by values of a set of multiple observables, O 1 ,Ô 2 , · · · ,Ô n , and and similarly for Ô full . One can expect that e iS I (Ô 1 ,··· ,Ôn) i 1 ,i 2 ,··· ,in has a single peak by introducing sufficiently many observables. In other words, the overlap problem can be solved by fixing sufficiently many observables. Suppose the quantity in consideration, sayÔ 1 , does not characterize the vacuum. Then the overlap problem is not solved. In the case of QCD, the overlap problem can be solved by takingÔ 2 to be the pion condensate. (When necessary one should also add ρ-condensate asÔ 3 etc.) But then we do not even have to fixÔ 1 , because it is not the source of the overlap problem anyways. Then, by lettingÔ 1 take any value, we arrive at our method. The point is that we only have to fix the quantities characterizing the correct vacuum, and for that purpose we have to understand the difference between physics of full and phase quenched theories.
For that, nonperturbative arguments like the large N c equivalence [11,12] play important roles. (In the supersymmetric matrix model studied in [19], they identified the observables which characterize the vacuum by using another numerical method, and then applied the multi-parameter factorization method.) Whether one has to fix multiple observables or not is a problem-specific issue, which depends on theories and parameter regions.
A good understanding about the vacuum structure of the full and phase-quenched theories is very important for this method to work. In the case of QCD, we already know the right quantity to fix. (We could also choose other quantities, but then we would have to fix multiple quantities, which makes actual calculation more difficult.) We know that only π + ∼ 0 is important, and can safely neglect the parameter region with small average sign. We do not even have to study large π + region. If we didn't know the right quantity to fix, we would have to measure the small sign rather precisely, in order to make sure that such parameter region is not important.
If we consider other theories for which the physical interpretation of the phase quenched theory is not clear, the simplest way to find 'good' observables would be to calculate various observables in the phase-quenched theory whose counterparts in the full theory trivially vanish due to symmetries. In case that the full theory is not understood well, one has to try purely numerical method: fix various quantities, scan the parameter space and find a nice peak structure of the average phase, as suggested in [19].

A comment on 'silver-blaze' region µ < µ c
In QCD B , at zero temperature and at the 'silver blaze' region µ < µ c , ν B must be zero. Therefore, at µ < µ c , it is possible to reduce the overlapping problem by setting ν B , rather than the pion condensate, to zero. To demonstrate this in the RMT, however, we should take ν B to be some negative value, since ν B becomes negative in this region due to an RMT-artifact. In the RMT I , the observable ν B corresponds to the isospin density ν I , which is positive. In Sec. 3.1, we also consider the behavior of the average phase as a function of ν B .

Simulation results
In this section we show the simulation results of N f = 2 RMT. In order to see the effect of the pion condensate, let us start with a nonzero mass. (When mass is zero, µ = 0 is already at the border of the pion condensation.) In Fig. 3 we show the distribution of the condition number, |(minimum eigenvalue of D f )/(maximum eigenvalue of D f )|, and the pion condensate π + , for N = 4, m = 0.35, and c = 0.02 , µ = 0 and µ = 0.7. (Note that the source c is introduced only for the constraint term ∆S.) At µ = 0, pion does not condense, and hence π + takes small values. At µ = 0.7, the distribution of π + has a long tail, which is the signature of the pion condensation. We can see that the condition number becomes smaller as π + increases, as expected from the fact that the pion condensate is caused by near zero-mode in the Dirac spectrum.
In Fig. 4 we show the average phase e iS I i , relative weight ρ i , and the reweighted relative weight ρ i · e iS I i for µ = 0.7. The weight in the phase-quenched simulation has a long tail reflecting the pion condensation. The average phase becomes small as π + becomes large, so that this fat tail is removed in the reweighted relative weight ρ i · e iS I i . Therefore the large-π + region gives only a negligible contribution; in fact, as shown in Fig. 5, π + <x ρ i · e iS I i and π + <x ρ i · e iS I · ν B i calculated at limited range of π + quickly converges. We can terminate the sum at around π + ∼ 0.15, so that we will not see small condition numbers. Next let us consider the massless limit, where the sign problem is severe. At each N , the value of the baryon density ν B can be calculated analytically [20]. For example, ν B = −180µ + 1440µ 3 − 5760µ 5 + 15360µ 7 − 24960µ 9 + 24576µ 11 − 14336µ 13 + 4096µ 15 45 − 360µ 2 + 1440µ 4 − 3840µ 6 + 7680µ 8 − 9984µ 10 + 8192µ 12 − 4096µ 14 + 1024µ 16 (28) for N = 4.

Figure 5:
π + <x ρ i · e iS I i and π + <x ρ i · e iS I · ν B i calculated at limited range of π + . N = 4, m = 0.35, µ = 0.7 and c = 0.02. The normalization is the same as in Fig. 4, i.e. the peak of ρ i · e iS I i is normalized to be 1.
Let us consider N = 4, µ = 0.7 and c = 0.02 as an example. We take = 0.01, x = 0.01, 0.02, 0.03, · · · . For each bin, we collected 1, 000, 000 configurations. In Fig. 6 we show the average phase e iS I i , relative weight ρ i , and the reweighted relative weight ρ i · e iS I i . The weight in the phase-quenched simulation has a long tail reflecting the pion condensation. However the average phase is extremely small at this tail and the reweighted relative weight does not have a fat tail. Also the peak is shifted to a small-π + region. In Fig. 7, ρ i · e iS I · ν B i and e iS I · ν B i are plotted. e iS I · ν B i behaves similarly to e iS I i : it approaches zero very quickly. As a result, ρ i · e iS I · ν B i does not have a fat tail either.
From Fig. 6 and Fig. 7, we can see that π + 0.08 is negligible. It can be explicitly seen from π + <x ρ i e iS I i and π + <x ρ i · e iS I ·ν B i shown in Fig. 8. In Fig. 9, we plot ν B | π + <x as a function of x. We can see a good convergence to the analytic value. In a usual phase-reweighting method, most computational resources are wasted to evaluate a very small average sign at π + 0.08, in order to prove that this region is not important. But from the beginning, we knew it is irrelevant. Then why do we have to waste resources there ?
Let us also see the plots at µ = 0.4, which is below µ c . As shown in Fig.10, the average phase, though small, seems to remain finite. However ρ i in the phase-quenched ensemble approaches zero faster than at µ = 0.7 at large π + , and the distribution after the reweighting is similar to that at µ = 0.7. As for the baryon density, the behavior is very different from the counterpart at µ = 0.7. As shown in Fig. 11, ν B e iS I i seems to take a nonzero value at large π + . However ρ i · ν B e iS I i goes to zero rather quickly because ρ i becomes zero. (Note that the baryon density takes a negative value here because of an artifact of RMT.) In Fig. 12, we compare our results of ν B π + <x at x = 0.10 and 0.30 with exact results at   π + <x ρ i · e iS I i and π + <x ρ i · e iS I · ν B i calculated at limited range of π + . N = 4, m = 0, µ = 0.7 and c = 0.02. The normalization is the same as in Fig. 6 and Fig. 7, i.e. the peak of ρ i · e iS I i is normalized to be 1.   : ν B · e iS I i and relative wight with and without phase, ρ i and ρ i · ν B · e iS I i . N = 4, m = 0, µ = 0.4 and c = 0.02. Although ν B · e iS I i seems to take a nonzero value at large π + , ρ i · ν B · e iS I i goes to zero rather quickly because ρ i becomes zero. Note that the baryon density takes a negative value here because of an artifact of RMT. several values of µ. From this figure we conclude that our method reproduces exact results quite well, though convergences are slower at µ < µ c . Figure 12: N = 4, m = 0 and c = 0.02, exact value vs. ν B π + <x , x = 0.10 and x = 0.30 for several values of µ. Data points for x = 0.10 are shifted to x-direction slightly so that they do not overlap with those for x = 0.30. The convergence to the exact value is slower at µ < µ c , because of a fatter tail.
Next let us consider N = 8, m = 0, µ = 0.7 and c = 0.02. We take = 0.01, x = 0.01, 0.02, 0.03, · · · . At x ≥ 0.05, we collected 10, 000, 000 -13, 800, 000 configurations for each bin. In the chiral limit, the phase fluctuation becomes severer as N increases: see Fig. 13 in which the average phase for N = 8 and N = 4 are shown. Still, the sign problem can be controlled by fixing π + to be small. In Fig. 14 we show ρ i · e iS I i . We can see the dominant contribution comes from the small-π + region. It is reasonable to omit configurations with π + 0.12, and there we can evaluate the baryon density reasonably well, as shown in Fig. 15.

Fixing ν B
As we have mentioned in Sec. 2.3.1, at µ < µ c , the baryon density ν B could be used to pin down the correct vacuum. (More precisely, we fix the real part, Re[ν B ].) So let us see the correlation between the average phase and ν B at µ = 0.4, which is below µ c , and at µ = 0.7, which is above µ c . In Fig. 16, the histogram of the real part of the baryon density Re[ν B ], and average phase at µ = 0.4 are shown. The average phase is larger at small Re[ν B ] region as expected. However, the average phase remain non-negligible even at large Re[ν B ] region, which suggests the baryon density is not as good observable as the pion condensate, though it could be used to make the corresponding density of states at µ < µ c .
Moreover, as shown in Fig. 17, at µ = 0.7 the average phase oscillates around zero, a very     complicated cancellation takes place and hence one has to study whole the configurations in order to estimate ν B precisely. In order to make this point clearer, we show Re[ν B ]<x ρ i · e iS I i and Re[ν B ]<x ρ i · e iS I · ν B i in Fig. 18. We can see a large x-dependence at 0 x 4. As a result, the convergence of ν B Re[ν B ]<x is very slow, as shown in Fig. 19. ν B Re[ν B ]<x becomes close enough to the correct value of ν B only at x 4, where almost all the configurations in the phase-quenched simulation are considered. Therefore ν B is not an appropriate observable to single out the correct vacuum at µ > µ c .

Strategies for the full QCD simulations
In this section we discuss a few strategies to apply our idea to the full QCD simulations.

Low-T , large-µ region
In the low temperature and high density region, one has to overcome the pion condensate by introducing the constraint term like (26). However since the simulation cost is not small, one needs to choose the constraint term in a clever manner so that efficient algorithms e.g. the Hybrid Monte Carlo (HMC) method are applicable. For that purpose, we introduce the gaussian term  again, but this time we do not set it to zero near a. Instead we take the above Gaussian form for all values of π + (x). (And, again, we introduce the source only for ∆S.) This four-fermi term can be made fermion bi-linear by introducing an auxiliary field, which allows us to apply the HMC method. A simple method for reconstructing the histogram of π + in the phase-quenched simulation with this deformation term can be found in [21].

High-T , small-µ region
The overlap problem is not severe at high-T , small-µ region. Still, at large volume, the sign fluctuation becomes very violent and simulation cost increases.
The origin of the overlap problem in this region is the gas of charged pion. Since the pion is light, the gas of pions can easily be excited with the isospin chemical potential, and hence the isospin density ν I takes non-negligible value. On the other hand, with the baryon chemical potential, only the gas of baryons can be excited. Since baryons are heavy, the baryon density ν B must be small. By recalling ν B in QCD B corresponds to ν I in QCD I (i.e. ν B B = ν I · e iS I I / e iS I I ), it is natural to think that the overlap problem can be suppressed by taking ν I small.
Given that the overlap problem is not severe compared to the low-T , large-µ region, important configurations with small ν I would be contained to some extent in the phase-quenched ensemble.
Therefore, with the re-analysis already existing configurations by calculating ν I , classifying configurations in terms of the values of ν I , and then applying our method, it would be possible to overcome the sign and overlapping problems. Note that one does not even have to calculate the determinant at large ν I , and hence it may reduce the cost for the reweighting, while increasing the accuracy.
However this method does not solve the overlap problem, because the pion condensation at µ > µ c takes place even in such reweighting calculation. In fact, the pion condensation has been observed even in the quench simulation, in which the configurations are generated by using pure Yang-Mills action without fermions. Therefore, configurations generated at µ 0 < µ c are not necessarily important ones at µ > µ c ; What one should actually do is to calculate the pion condensate π + at (+µ, −µ) (not at (+µ 0 , −µ 0 )), classify the configurations and apply our method by using configurations with small π + . Note again that one does not even have to calculate the determinant at large π + , and hence it is possible to reduce the cost for the reweighting, while increasing the accuracy.
In the phase reweighting method, once the pion condensate is set zero the sign problem is 1/N c suppressed. For generic reweightings given in eq. (31), such a nice property does not exist, and hence more violent cancellation is expected. It should become severer as one goes deep inside the pion condensation. Still, however, it would be useful to study the chiral and deconfinement transitions by stepping a little bit inside. This can be a very important application practically, since the QCD critical point might be there.

Conclusion and future directions
In this paper, we have improved the density of state method which can make the sign problem milder. As a by-product, the problem of the zero-mode in the finite-density QCD, which is associated with the pion condensation in the phase quenched theory, can be avoided. We demonstrated our idea thoroughly by using numerically cheaper and analytically solvable toy model, the chiral random matrix theory.
There are variety of directions for future studies. For the finite-density QCD, it is important to study high-T low-µ region thoroughly. In addition to the confirmation of the effectiveness of the method, it would provide us with better numerical understanding about the nature of the QCD thermal transition. For this, we can use the simplified method described in Sec. 4.2 which does not require generation of new configurations. It is also interesting to go a little bit into the pion condensation region in order to search the QCD critical point, by using the method described in Sec. 5. Again, we do not need new configurations; re-analysis of existing configurations can provide us with better results. Needless to say, the study of the low-temperature high-density region sketched in Sec. 4.2 is the most interesting thing to do. Although the remaining sign problem would become severe at large volume, interesting phenomena would be seen already at small volume. Note that even the phase quench simulation can work up to the 1/N c -correction, once the pion condensation is erased [11,12].
Our method is quite general. We hope to report other applications, such as models in condensed matter physics and supersymmetric gauge theory, in near future.
Science Foundation under Grant No. PHYS-1066293 and the hospitality of the Aspen Center for Physics.