text:
# 1 Introduction

## 1 Introduction

The (Euclidean) action of many interesting systems in fields ranging from condensed matter physics to high-energy physics has an imaginary part. Some examples in high-energy physics are QCD at finite baryon density, Chern-Simons theories, systems with topological terms (like the $`\theta `$-term in QCD) and systems with chiral fermions. While this is not a conceptual problem per se, it severely limits the application of Monte Carlo methods, which otherwise might provide a powerful tool to understand the properties of these systems from first principles. So far there is no general solution to this ‘complex action problem’.

In Ref.  a new Monte Carlo approach to systems with a complex action was proposed. This method utilizes a simple factorization property of distribution functions of observables. Since the property holds quite generally, the approach can be applied to any system with a complex action. Most notably, the method eliminates the so-called overlap problem, which occurs when one applies the standard re-weighting technique to include the effect of the imaginary part. Ultimately we hope that this method will enable us, among other things, to explore the phase diagram of QCD at finite baryon density, where interesting phases such as a superconducting phase have been conjectured to appear .

As a first step toward achieving this goal we test the new approach in a Random Matrix Theory for finite density QCD . Random Matrix Theory was originally introduced to describe the spectrum of the Dirac operator at zero chemical potential  and has been studied intensively in the literature. (See Ref.  for a review). The particular extended model we study can be regarded as a schematic model for QCD at finite baryon density. As one increases the ‘chemical potential’, the model undergoes a first order phase transition, where the imaginary part of the action plays a crucial role. Since it is solvable even for finite matrix size $`N`$ , it serves as a useful testing ground for simulation techniques for QCD at finite density. For instance, the problem with quenched simulations  has been clarified in Ref. . It was also used to test the so-called Glasgow method , and the source of the problems was identified .

In this article we apply the new method to both phases of this model (below and above the critical point) and obtain the expectation value of the ‘quark number density’. The results nicely reproduce the exact results known for finite $`N`$. The values of $`N`$ that are accessible by the new method turn out to be large enough to extract the large-$`N`$ limit. Moreover, our results provide a clear understanding of how the first order phase transition is induced by the imaginary part of the action.

The method  proposed for simulating $`\theta `$-vacuum like systems can be regarded as a special case of the factorization method. A simplified version of the method was sufficient because the observable was identical to the imaginary part of the action. The essence of the factorization method is that it avoids the overlap problem by the use of constrained simulations. Results for 2d CP3 and other models are very promising.

The remainder of this article is organized as follows. In Section 2, we define a Random Matrix Theory for QCD at finite baryon density and review the known exact results. In Section 3, we explain the complex action problem associated with the standard re-weighting technique. The application of the factorization method to a Monte Carlo study of the Random Matrix Model is discussed in Section 4. Section 5 shows the results which nicely reproduce the exact values for the quark number density. Section 6 is devoted to summary and discussions.

## 2 Random Matrix Theory for finite density QCD

The Random Matrix Model we study in this article is defined by the partition function

$$Z=\text{d}W\text{e}^{N\text{tr}(W^{}W)}detD,$$
(2.1)
where $`W`$ is a $`N\times N`$ complex matrix, and the $`D`$ is a $`2N\times 2N`$ matrix given by

$$D=\left(\begin{array}{cc}m& iW+\mu \\ iW^{}+\mu & m\end{array}\right).$$
(2.2)
The parameters $`m`$ and $`\mu `$ correspond to the ‘quark mass’ and the ‘chemical potential’, respectively. The size of the matrix can be thought of as the total number of ‘low-lying’ modes for given total volume. Since the density of these modes is taken to be unity, $`N`$ can be interpreted as the volume of space time. This model has the global symmetries of the Dirac operator of QCD at nonzero baryon chemical potential, where the matrix $`W`$ in (2.2) is a complicated function of the background gauge field. Thus the above model can be thought of as a schematic model for QCD at finite baryon density, where the path integral over the gauge field is simply replaced by the Gaussian integral over $`W`$. Interesting observables are the ‘chiral condensate’ and the ‘quark number density’ defined by

$`\mathrm{\Sigma }`$ $`=`$ $`{\displaystyle \frac{1}{2N}}\mathrm{tr}(D^1)`$ (2.3)
$`\nu `$ $`=`$ $`{\displaystyle \frac{1}{2N}}\mathrm{tr}(\gamma _4D^1),\text{ }\gamma _4=\left(\begin{array}{cc}0& 11\\ 11& 0\end{array}\right).`$ (2.6)
The model was first solved in the large-$`N`$ limit . Later it was noticed that the model can be solved even for finite $`N`$ . Throughout this paper, we consider the massless case ($`m=0`$) for simplicity. Then the partition function can be expressed as

$$Z(\mu )=\pi \text{e}^\kappa N^{(N+1)}N!\left[1+\frac{(1)^{N+1}}{N!}\gamma (N+1,\kappa )\right],$$
(2.7)
where $`\kappa =N\mu ^2`$ and $`\gamma (n,x)`$ is the incomplete $`\gamma `$-function defined by

$$\gamma (n,x)=_0^x\text{e}^tt^{n1}𝑑t.$$
(2.8)
From this one obtains the vacuum expectation value (VEV) of the quark number density as

$`\nu `$ $`=`$ $`{\displaystyle \frac{1}{2N}}{\displaystyle \frac{}{\mu }}\mathrm{ln}Z(\mu )`$ (2.9)
$`=`$ $`\mu \left[1+{\displaystyle \frac{\kappa ^N\text{e}^\kappa }{(1)^{N+1}N+\gamma (N+1,\kappa )}}\right].`$ (2.10)
In Fig. 1 we plot $`\nu `$ as a function of the chemical potential $`\mu `$ for $`N=8,16,32,64,128`$. The large-$`N`$ limit of this formula can be taken easily by applying the saddle-point analysis to the incomplete $`\gamma `$-function. We obtain

$$\underset{N\mathrm{}}{lim}\nu =\{\begin{array}{cc}\mu \hfill & \text{for }\mu <\mu _\mathrm{c}\hfill \\ 1/\mu \hfill & \text{for }\mu >\mu _\mathrm{c},\hfill \end{array}$$
(2.11)
where $`\mu _\mathrm{c}`$ is the solution to the equation $`1+\mu ^2+\mathrm{ln}(\mu ^2)=0`$, and its numerical value is given by $`\mu _\mathrm{c}=0.527\mathrm{}`$. We find that the quark number density $`\nu `$ has a discontinuity at $`\mu =\mu _\mathrm{c}`$. Thus the schematic model reproduces qualitatively the first order phase transition expected to occur in ‘real’ QCD at nonzero baryon density.

## 3 The complex action problem

In this section we describe the complex action problem that appears in standard Monte Carlo studies of the model (2.1). Let us first rewrite (2.1) as

$$Z=\text{d}W\text{e}^{S_0+i\mathrm{\Gamma }},$$
(3.1)
where we have introduced $`S_0`$ and $`\mathrm{\Gamma }`$ by

$`S_0`$ $`=`$ $`N\mathrm{tr}(W^{}W)\mathrm{ln}|detD|`$ (3.2)
$`detD`$ $`=`$ $`\text{e}^{i\mathrm{\Gamma }}|detD|.`$ (3.3)
In this form it becomes manifest that the system has a complex action, where the problematic imaginary part $`\mathrm{\Gamma }`$ is given by the phase of the fermion determinant. Since the weight $`\text{e}^{S_0+i\mathrm{\Gamma }}`$ in (3.1) is not positive definite, we cannot regard it as a probability density. Hence it seems difficult to apply the idea of standard Monte Carlo simulations, which reduces the problem of obtaining VEVs to that of taking an average over an ensemble generated by the probability density. One way to proceed is to apply the reweighting method and rewrite the VEV $`\nu `$ as

$$\nu =\frac{\nu \text{e}^{i\mathrm{\Gamma }}_0}{\text{e}^{i\mathrm{\Gamma }}_0},$$
(3.4)
where the symbol $`_0`$ denotes a VEV with respect to the so-called phase quenched partition function

$$Z_0=\text{d}W\text{e}^{N\text{tr}(W^{}W)}|detD|=\text{d}W\text{e}^{S_0}.$$
(3.5)
Since the system (3.5) has a positive definite weight, the VEV $`_0`$ can be evaluated by standard Monte Carlo simulations. However, the fluctuations of the phase $`\mathrm{\Gamma }`$ in (3.4) grows linearly with the size of the matrix $`D`$, which is of O($`N`$). Due to huge cancellations, both the denominator and the numerator of the r.h.s. of (3.4) vanish as $`\text{e}^{\mathrm{const}.N}`$ as $`N`$ increases, while the ‘observables’ $`\text{e}^{i\mathrm{\Gamma }}`$ and $`\nu \text{e}^{i\mathrm{\Gamma }}`$ are of O(1) for each configuration. As a result, the number of configurations required to obtain the VEVs with some fixed accuracy grows as $`\text{e}^{\mathrm{const}.N}`$ (We remind the reader that $`N`$ can be considered as the volume of space time). This is the notorious ‘complex action problem’ (or rather the ‘sign problem’, as we see below). See Refs.  for simulation results for ‘real’ finite density QCD obtained by this reweighting technique.

In fact we may simplify the expression (3.4) slightly by using a symmetry. We note that the fermion determinant $`detD`$, as well as the observable $`\nu `$, becomes complex conjugate under the transformation

$$WW,$$
(3.6)
while the Gaussian action remains invariant. From this we find that

$`\nu =\nu _\mathrm{R}+i\nu _\mathrm{I}`$ (3.7)
$`\nu _\mathrm{R}={\displaystyle \frac{\nu _\mathrm{R}\mathrm{cos}\mathrm{\Gamma }_0}{\mathrm{cos}\mathrm{\Gamma }_0}};\nu _\mathrm{I}=i{\displaystyle \frac{\nu _\mathrm{I}\mathrm{sin}\mathrm{\Gamma }_0}{\mathrm{cos}\mathrm{\Gamma }_0}},`$ (3.8)
where $`\nu _\mathrm{R}`$ and $`\nu _\mathrm{I}`$ denote the real part and the imaginary part of $`\nu `$, respectively. In Eq. (3.8) the problem takes the form of the ‘sign problem’, since $`\mathrm{cos}\mathrm{\Gamma }`$ and $`\mathrm{sin}\mathrm{\Gamma }`$ flip their sign violently as a function of the configuration $`W`$. Note that both terms in the r.h.s. of (3.7) are real, meaning in particular that their sum $`\nu `$ is also real.

The model (3.5) is solvable in the large-$`N`$ limit . For $`m=0`$ one obtains

$$\underset{N\mathrm{}}{lim}\nu _0=\{\begin{array}{cc}\mu \hfill & \text{for }\mu <1\hfill \\ 1/\mu \hfill & \text{for }\mu >1.\hfill \end{array}$$
(3.9)
In this case the VEV of the quark number density is a continuous function of the chemical potential $`\mu `$ unlike in (2.11). This remarkable difference between (2.11) and (3.9) is precisely due to the imaginary part $`\mathrm{\Gamma }`$ of the action. The two results agree trivially at $`\mu =0`$, but interestingly they also agree at $`\mu >1`$. This is because the eigenvalues of $`W`$ are located inside the complex unit circle . If $`\mu >1`$, the effect of the fermion determinant is $`1/N`$ suppressed, and the quenched approximation, as well as the phase quenched approximation, becomes exact in the thermodynamic limit. Note also that the symmetry under (3.6) implies

$$\nu _\mathrm{I}_0=0;\nu _\mathrm{R}_0=\nu _0.$$
(3.10)
We can also predict $`\nu _\mathrm{R}`$ and $`\nu _\mathrm{I}`$ separately for the unquenched model in the large $`N`$ limit. For that we note that $`\nu (\mu )=\nu (\mu )^{}`$. Therefore we have

$`\nu _\mathrm{R}={\displaystyle \frac{1}{2}}\left\{\nu (\mu )\nu (\mu )\right\},\nu _\mathrm{I}={\displaystyle \frac{1}{2i}}\left\{\nu (\mu )+\nu (\mu )\right\},`$ (3.11)
where in the calculation of $`\nu (\mu )`$ the sign of the chemical potential in the fermion determinant is not changed. In calculating $`\nu (\mu )`$, the fermion determinant $`detD(\mu )`$ cancels the poles of $`\nu (\mu )`$. This is responsible for the difference between quenched and unquenched results for $`\mu <1`$. However, the cancellation of the poles does not occur in calculating $`\nu (\mu )`$. Therefore we expect that $`\nu (\mu )`$ coincides with the corresponding quenched result $`\nu (\mu )_0`$ even for $`\mu <1`$. This implies in particular that for $`\mu <\mu _\mathrm{c}`$ we obtain the results $`\nu _\mathrm{R}=0`$ and $`i\nu _\mathrm{I}=\mu `$ in sharp contrast to the quenched result (3.10). For $`\mu >1`$, on the other hand, both $`\nu _\mathrm{R}`$ and $`\nu _\mathrm{I}`$ agree with the corresponding quenched results. These results are indeed observed in our Monte Carlo simulations, which shall be discussed in Section 5.

## 4 The factorization method

### 4.1 The basic formulae

In this section, we explain how the factorization method  can be used to obtain the VEVs $`\nu _\mathrm{R}`$ and $`\nu _\mathrm{I}`$. The fundamental objects of the method are the distribution functions

$$\rho _\mathrm{R}(x)\stackrel{\mathrm{def}}{=}\delta (x\mathrm{Re}(\nu ))\mathrm{and}\rho _\mathrm{I}(y)\stackrel{\mathrm{def}}{=}\delta (y\mathrm{Im}(\nu ))$$
(4.1)
associated with the complex valued $`\nu `$ . In a unified notation that will be used below these equations can be rewritten as

$$\rho _i(x)\stackrel{\mathrm{def}}{=}\delta (x\nu _i)i=\mathrm{R},\mathrm{I}.$$
(4.2)
In terms of these functions, the VEVs of $`\nu _i`$ can be expressed as

$$\nu _i=_{\mathrm{}}^{\mathrm{}}\text{d}xx\rho _i(x).$$
(4.3)
What is essential for the method are the constrained partition functions

$$Z_i(x)=\text{d}W\text{e}^{S_0}\delta (x\nu _i)$$
(4.4)
and the average of the phase $`\text{e}^{i\mathrm{\Gamma }}`$ with respect to these partition functions

$$\phi _i(x)\stackrel{\mathrm{def}}{=}\text{e}^{i\mathrm{\Gamma }}_{i,x}.$$
(4.5)
If we also introduce the distribution functions for the phase quenched model by

$$\rho _i^{(0)}(x)\stackrel{\mathrm{def}}{=}\delta (x\nu _i)_0,$$
(4.6)
the distribution functions for the unquenched model trivially factorize as

$$\rho _i(x)=\frac{1}{C}\rho _i^{(0)}(x)\phi _i(x)i=\mathrm{R},\mathrm{I},$$
(4.7)
where the normalization constant $`C`$ is given by

$$C\stackrel{\mathrm{def}}{=}\text{e}^{i\mathrm{\Gamma }}_0.$$
(4.8)
Plugging (4.7) into (4.3), we get

$$\nu _i=\frac{1}{C}_{\mathrm{}}^{\mathrm{}}\text{d}xx\rho _i^{(0)}(x)\phi _i(x).$$
(4.9)
Note also that Eq. (4.7) implies that the constant $`C`$ can be written as

$$C=_{\mathrm{}}^{\mathrm{}}\text{d}x\rho _\mathrm{R}^{(0)}(x)\phi _\mathrm{R}(x).$$
(4.10)
Thus the VEVs $`\nu _i`$ ($`i=\mathrm{R},\mathrm{I}`$) can be expressed solely in terms of the functions $`\rho _i^{(0)}(x)`$ and $`\phi _i(x)`$ ($`i=\mathrm{R},\mathrm{I}`$), which can be calculated by standard Monte Carlo simulations. Clearly the derivation presented here is quite general, and we can obtain similar formulae in any complex action system.

In the present case, we can slightly simplify the formulae due to the symmetry under (3.6). Using the properties

$`\phi _\mathrm{R}(x)^{}`$ $`=`$ $`\phi _\mathrm{R}(x),`$ (4.11)
$`\phi _\mathrm{I}(x)^{}`$ $`=`$ $`\phi _\mathrm{I}(x),`$ (4.12)
$`\rho _\mathrm{I}^{(0)}(x)`$ $`=`$ $`\rho _\mathrm{I}^{(0)}(x),`$ (4.13)
which follow from the symmetry, we arrive at

$`\nu _\mathrm{R}`$ $`=`$ $`{\displaystyle \frac{1}{C}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}\text{d}xx\rho _\mathrm{R}^{(0)}(x)w_\mathrm{R}(x),`$ (4.14)
$`\nu _\mathrm{I}`$ $`=`$ $`{\displaystyle \frac{2i}{C}}{\displaystyle _0^{\mathrm{}}}\text{d}xx\rho _\mathrm{I}^{(0)}(x)w_\mathrm{I}(x),`$ (4.15)
$`C`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}\text{d}x\rho _\mathrm{R}^{(0)}(x)w_\mathrm{R}(x),`$ (4.16)
where the weight functions $`w_i(x)`$ are defined by

$`w_\mathrm{R}(x)`$ $`\stackrel{\mathrm{def}}{=}`$ $`\mathrm{cos}\mathrm{\Gamma }_{\text{R},x},`$ (4.17)
$`w_\mathrm{I}(x)`$ $`\stackrel{\mathrm{def}}{=}`$ $`\mathrm{sin}\mathrm{\Gamma }_{\text{I},x}.`$ (4.18)
Thus the problem reduces to the calculation of the four real functions $`\rho _i^{(0)}(x)`$, $`w_i(x)`$ ($`i=\mathrm{R},\mathrm{I}`$), which we will now discuss.

### 4.2 Monte Carlo evaluation of $`\rho _i^{(0)}(x)`$ and $`w_i(x)`$

In order to obtain $`w_i(x)`$ ($`i=\mathrm{R},\mathrm{I}`$), we need to simulate (4.4). In practice, we simulate a partition function where the $`\delta `$-function is replaced by a sharply peaked potential

$$Z_{i,V}=\text{d}W\text{e}^{S_0}\text{e}^{V(\nu _i)}.$$
(4.19)
In this study we use a Gaussian potential

$$V(x)=\frac{1}{2}\gamma (x\xi )^2,$$
(4.20)
where $`\gamma `$ and $`\xi `$ are real parameters. The results are insensitive to the choice of $`\gamma `$ as far as they are big enough (we used $`\gamma =1000.0`$). Let us denote the VEV associated with the partition function (4.19) as $`𝒪_{i,V}`$. The expectation value $`\mathrm{cos}\mathrm{\Gamma }_{\mathrm{R},V}`$ represents the value of $`w_\mathrm{R}(x)`$ at $`x=\nu _\mathrm{R}_{\mathrm{R},V}`$, while the expectation value $`\mathrm{sin}\mathrm{\Gamma }_{\mathrm{I},V}`$ represents the value of $`w_\mathrm{I}(x)`$ at $`x=\nu _\mathrm{I}_{\mathrm{I},V}`$.

In fact we can also obtain the functions $`\rho _i^{(0)}(x)`$ from the same simulation (4.19). For that we note that the distribution of $`\nu _i`$ in the system (4.19) is given by

$$\rho _{i,V}(x)\stackrel{\mathrm{def}}{=}\delta (x\nu _i)_{i,V}\rho _i^{(0)}(x)\text{e}^{V(x)},$$
(4.21)
which typically has a peak. The position of the peak, which we denote as $`x=\stackrel{~}{x}`$, is given by the solution to

$$0=\frac{d}{dx}\mathrm{ln}\rho _i(x)=f_i^{(0)}(x)V^{}(x),$$
(4.22)
where we have introduced

$$f_i^{(0)}(x)\stackrel{\mathrm{def}}{=}\frac{d}{dx}\mathrm{ln}\rho _i^{(0)}(x).$$
(4.23)
Therefore, the quantity $`V^{}(\stackrel{~}{x})`$ represents the value of $`f_i^{(0)}(x)`$ at $`x=\stackrel{~}{x}`$. Since the value of $`\gamma `$ is taken to be large, the peak is sharp and we can safely approximate $`\stackrel{~}{x}`$ by the VEV $`\nu _i_{i,V}`$. Once we obtain $`f_i^{(0)}(x)`$ for various $`x`$, we can calculate $`\rho _i^{(0)}(x)`$ by integrating (4.23) and exponentiating, where the integration constant can be determined by the normalization of the distribution $`\rho _i^{(0)}(x)`$. If the distribution $`\rho _i^{(0)}(x)`$ would be Gaussian (as it is the case near its peak; see Appendix), the above procedure is exact even for finite $`\gamma `$.

Monte Carlo simulation of (4.19) can be performed by using the Hybrid Monte Carlo algorithm in much the same way as in Refs. . The required computational effort for the present model is O($`N^3`$).

### 4.3 The virtues of the method

Comparing (4.14), (4.15) and (4.16) with Eq. (3.8), we notice that

$`\nu _\mathrm{R}\mathrm{cos}\mathrm{\Gamma }_0`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}\text{d}xx\rho _\mathrm{R}^{(0)}(x)w_\mathrm{R}(x),`$ (4.24)
$`\nu _\mathrm{I}\mathrm{sin}\mathrm{\Gamma }_0`$ $`=`$ $`2{\displaystyle _0^{\mathrm{}}}\text{d}xx\rho _\mathrm{I}^{(0)}(x)w_\mathrm{I}(x),`$ (4.25)
$`\mathrm{cos}\mathrm{\Gamma }_0`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}\text{d}x\rho _\mathrm{R}^{(0)}(x)w_\mathrm{R}(x).`$ (4.26)
Thus the new method as it stands simply amounts to using the standard reweighting formula (3.8), but calculating each VEV by using (4.24), (4.26) and (4.26). We now explain the virtues of the present method.

If we are to obtain the VEVs by directly simulating the system (3.5), for most of the time we sample configurations whose $`\nu _i`$ takes a value close to the peak of $`\rho _i^{(0)}(x)`$. However, from the r.h.s. of (4.24)$``$(4.26) it is clear that we have to sample configurations whose $`\nu _i`$ takes a value where $`\rho _i^{(0)}(x)|w_i(x)|`$ becomes large, in order to obtain the VEVs accurately. In general these two regions of configuration space have little overlap, which makes the sampling ineffective. This is the overlap problem. Since the overlap becomes exponentially small as the system size increases, this composes some portion of the complex action problem. The use of (4.24)$``$(4.26) and calculating the relevant functions as explained in the previous section avoids this problem by ‘forcing’ the simulation to sample the important region.

The knowledge of the weight factor $`w_i(x)`$, which is provided by the present approach, allows us to probe directly the effect of the imaginary part $`\mathrm{\Gamma }`$ on the observable of our concern. We can understand which values of the observable are enhanced or suppressed by the effect of $`\mathrm{\Gamma }`$. On the other hand, the standard reweighting technique simply gives the integrals on the l.h.s. of (4.24)$``$(4.26), from which we can hardly imagine how they resulted from the effect of $`\mathrm{\Gamma }`$. In the particular model we are studying, the imaginary part $`\mathrm{\Gamma }`$ has a dramatic effect on the VEV of the quark number as we discussed in Section 3. We will see in the next section that the weight factor $`w_i(x)`$ indeed provides a clear understanding of this phenomenon.

As we also see in the next section, a Monte Carlo calculation of the weight factor $`w_i(x)`$ becomes increasingly difficult with the system size. In that sense the complex action problem is still there. This should be contrasted with the meron-cluster algorithm , which has been applied to a special class of complex action systems with computational efforts increasing at most by some power of the system size. However, as mentioned, the factorization method eliminates the overlap problem and this is a substantial step forward which allows us to get closer to the thermodynamic limit. Indeed, we are able to obtain the thermodynamic limit of the ‘quark number density’ in this random matrix model with modest computer resources. See also Ref.  for an idea to ameliorate the overlap problem in ‘real’ QCD at finite baryon density by interpolations in the ($`\mu `$, $`T`$) plane (notice however ).

Using the generic scaling properties of the weight factor $`w_i(x)`$, one may extrapolate the results obtained by direct Monte Carlo evaluations to larger system size. Such an extrapolation is expected to be particularly useful in cases where the distribution function turns out to be positive definite. In those cases we can actually even avoid using the reweighting formula (3.8) by reducing the question of obtaining the expectation value to that of finding the minimum of the free energy, which is (minus) the log of the distribution function. Here, the error in obtaining the scaling function propagates to the final result without significant magnifications. Therefore, the extrapolation can be a powerful tool to probe the thermodynamic limit from the accessible system size. Indeed such a technique has been used in Ref.  to discuss the dynamical generation of four-dimensional space-time in a nonperturbative formulation  of type IIB superstring theory in ten dimensions.

## 5 Reproducing exact results by the new method

In this section we apply the factorization method to the model (2.1) as described in the previous section and show that it reproduces the known results $`\nu `$ given by (2.10) for arbitrary $`N`$. Let us recall the dramatic difference between (2.11) and (3.9). It is particularly interesting to see how this occurs due to the effects of $`\mathrm{\Gamma }`$ in the present approach. First we focus on two values of $`\mu `$, $`\mu =0.2`$ and $`\mu =1.0`$, which are on opposite sides of the first order phase transition point $`\mu =\mu _\mathrm{c}=0.527\mathrm{}`$.

### 5.1 $`\mu <\mu _\mathrm{c}`$

We start with the results for $`\mu =0.2`$. In Figs. 2 and 3, we plot $`w_\mathrm{R}(x)`$ and $`w_\mathrm{I}(x)`$ respectively. In Figs. 4 and 5, we show the results for the functions $`f_i^{(0)}(x)`$. (The behavior of these functions can be understood theoretically as discussed in the Appendix.) By integrating these functions and exponentiating, we obtain the $`\rho _i^{(0)}(x)`$, which is shown in Figs. 6 and 7. In Fig. 8 and 9, we show $`\rho _\mathrm{R}^{(0)}(x)w_\mathrm{R}(x)`$ and $`\rho _\mathrm{I}^{(0)}(x)w_\mathrm{I}(x)`$, respectively. Using them we obtain $`\nu _\mathrm{R}`$ and $`i\nu _\mathrm{I}`$. Summing these values, we get $`\nu `$, which should be compared with the exact result obtained from Eq. (2.10). The results are shown in Table 1.

Note that the sign change of $`w_\mathrm{R}(x)`$ occurs near the peak of $`\rho _\mathrm{R}^{(0)}(x)`$, so that the product $`\rho _\mathrm{R}^{(0)}(x)w_\mathrm{R}(x)`$ has a positive regime and a negative regime, which cancel each other resulting in $`\nu _\mathrm{R}0`$. Thus the main contribution to $`\nu `$ comes from the imaginary part $`\nu _\mathrm{I}`$ in contrast to the results (3.10) for the phase quenched system. This is consistent with the theoretical argument at the end of Section 3 that the contribution to $`\nu `$ comes solely from the imaginary part in the large $`N`$ limit.

### 5.2 $`\mu >\mu _\mathrm{c}`$

We perform a similar analysis at $`\mu =1.0`$. The results are presented in Figs. 10, …, 17. Again the results nicely reproduce the exact results as can be seen in Table 1. Note in particular that the finite $`N`$ effect is 1.6% at $`N=32`$, meaning that the accessible values of $`N`$ are large enough to extract the large $`N`$ limit.

Here, $`w_\mathrm{R}(x)`$ is approximately constant in the region where $`\rho _\mathrm{R}^{(0)}(x)`$ is peaked, so the shape of the product $`\rho _\mathrm{R}^{(0)}(x)w_\mathrm{R}(x)`$ is similar to $`\rho _\mathrm{R}^{(0)}(x)`$. On the other hand, the peak of $`\rho _\mathrm{I}^{(0)}(x)`$ at $`x=0`$ is slightly shifted by multiplying $`w_\mathrm{I}(x)`$, but the first moment of the product $`\rho _\mathrm{I}^{(0)}(x)w_\mathrm{I}(x)`$ is still small. This is also the case for $`\mu =0.2`$, but the difference comes from the normalization constant $`C`$ in the formula (4.15) for $`\nu _\mathrm{I}`$. As in (4.16), the constant $`C`$ is obtained by integrating $`\rho _\mathrm{R}^{(0)}(x)w_\mathrm{R}(x)`$, where cancellations occur at $`\mu =0.2`$, but not at $`\mu =1.0`$. As a result we obtain $`\nu _\mathrm{I}0`$ at $`\mu =1.0`$ as $`N`$ increases. Thus the main contribution to $`\nu `$ comes from the real part $`\nu _\mathrm{R}`$, and moreover, it is close to $`\nu _\mathrm{R}_0`$. Again this is consistent with the theoretical argument given at the end of Section 3.

It is interesting that the $`w_\mathrm{R}(x)`$ changes from positive to negative for $`\mu =0.2`$, but it changes from negative to positive for $`\mu =1.0`$. Similarly $`w_\mathrm{I}(x)`$ is positive at $`x>0`$ for $`\mu =0.2`$, but it is negative at $`x>0`$ for $`\mu =1.0`$. Thus the behavior of $`w_i(x)`$ changes drastically as the chemical potential $`\mu `$ crosses its critical value $`\mu _\mathrm{c}`$.

### 5.3 $`\mu \mu _\mathrm{c}`$

Let us see in more detail what is going on in the critical regime. In Figs. 18 and 19 we plot $`w_\mathrm{R}(x)`$ and $`w_\mathrm{I}(x)`$ respectively for $`N=8`$ at various $`\mu `$. The final results for $`\nu _i`$ are plotted in Fig. 20 and the corresponding data are listed in Table 2. Note that we were able to reproduce the exact result even at the turning point ($`\mu =0.6139`$). These results provide a clear understanding of how the first order phase transition occurs due to the effects of $`\mathrm{\Gamma }`$.

The fact that $`|w_i(x)|`$ becomes small near the critical regime reveals an increasing difficulty in approaching the critical point. In Figs. 21, 22, we plot $`\mathrm{ln}(\mathrm{max}_x|w_\mathrm{I}(x)|)`$ against $`|\mu \mu _\mathrm{c}|`$ at $`N=16,32`$ for $`\mu <\mu _\mathrm{c}`$ and $`\mu >\mu _\mathrm{c}`$ respectively. Our data can be nicely fitted to

$$\mathrm{ln}(\underset{x}{\mathrm{max}}|w_\mathrm{I}(x)|)=a\mathrm{exp}(b|\mu \mu _\mathrm{c}|)+\mathrm{const}..$$
(5.1)
For fixed $`\mu `$, on the other hand, we observe that $`\mathrm{max}_x|w_\mathrm{I}(x)|`$ decreases exponentially as $`\mathrm{exp}(cN)`$. We plot the result for $`\mu =0.2`$ in Fig. 23. This is consistent with the scaling behavior of the weight function, which plays a crucial role in the extrapolation  as we mentioned at the end of Section 4.3.

## 6 Summary

In this article we have clarified the properties of the factorization method for systems with a complex action. This method circumvents the overlap problem, and we hope that this feature alone will make various interesting questions accessible by the present day computer resources. Indeed, we were able to reproduce the exact results for the quark number density in a schematic model for QCD at finite baryon density. The achieved system size was already large enough to obtain the thermodynamic limit. We therefore expect that the factorization method is useful to explore the phase diagram of the ‘real’ finite density QCD. The method itself is quite general, and it can be applied to any system with a complex action, although the actual gains may depend on the system.

We also emphasize that in the case where the distribution functions turn out to be positive definite, the method can be even more powerful by utilizing the scaling property of the weight factor. This extrapolation appeared particularly valuable in the study of spontaneous symmetry breaking in the type IIB matrix model. Although this is a very interesting problem which would tell us a great deal about nonperturbative string dynamics and the dynamical origin of the space-time dimensionality, we expect that this method will be useful for other systems as well.

## Acknowledgments

We would like to thank W. Bietenholz, Z. Burda, S. Iso, Y. Kikukawa, E. Kiritsis, T. Onogi and B. Peterssen for discussions. The authors are also grateful to I. Kostov, P. Di Francesco and other organizers of the Sixth Claude Itzykson Meeting (Matrix Models 2001), which triggered the current collaboration. K.N.A.’s research was partially supported by RTN grants HPRN-CT-2000-00122, HPRN-CT-2000-00131 and HPRN-CT-1999-00161 and the INTAS contract N 99 0590. The work of J.N. was supported in part by Grant-in-Aid for Scientific Research (No. 14740163) from the Ministry of Education, Culture, Sports, Science and Technology. The work of J.J.M.V. was partially supported by the US DOE grant DE-FG-88ER40388.

## Appendix: Large $`N`$ behavior of $`\rho _i^{(0)}(x)`$

In this Appendix, we discuss the large $`N`$ behavior of functions $`\rho _i^{(0)}(x)`$. From Figures 4, 5, 12, 13, we find that $`\rho _i^{(0)}(x)`$ is well approximated by the Gaussian distribution near the peak, but there is a transition to a power-like tail ($`\rho _i^{(0)}(x)x^4`$) at large $`|x|`$. The function $`f_i^{(0)}(x)`$ scales in this power regime. In the Gaussian regime, on the other hand, the function $`\frac{1}{N}f_i^{(0)}(x)`$, with the normalization factor $`1/N`$, scales as is shown in Figs. 24 and 25 for $`\mu =0.2`$. As a result the extent of the Gaussian regime shrinks as $`1/\sqrt{N}`$. Such a behavior is obtained if $`N`$ independently distributed eigenvalues of $`W`$ contribute to $`\nu `$. Indeed, the correlations of the eigenvalues of the matrix $`W`$ in our model decay exponentially on the scale of the average level spacing.

The observed power tail can be understood as follows. The point is that large values of the baryon number result from eigenvalues of the matrix $`W`$ close to $`\pm i\mu `$. The probability of finding one eigenvalue $`\lambda `$ inside a radius $`\alpha `$ can be easily obtained for $`\mu =0`$ from generalizing a calculation in the book of Mehta  (the chapter on RMT’s without hermiticity conditions). One finds that

$$P(\lambda <\alpha )=\frac{N^2}{2}\alpha ^4\mathrm{for}\alpha 0.$$
(A.1)
If we assume stationarity of this result we obtain

$$P(|\nu |>x)=P\left(\lambda ,|\lambda i\mu |<\frac{1}{2Nx}\right)=\frac{1}{32N^2x^4}.$$
(A.2)
For the probability density we thus find

$$P(|\nu |=x)=\frac{1}{8N^2x^5},$$
(A.3)
so that the real and imaginary parts of the quark number density behave as $`1/x^4`$ for large $`|x|`$. Since the level spacing distribution is a universal feature of Random Matrix Theories, we expect that such power-like tails are generic and also occur in QCD.

In fact our results in Figures 4, 5, 12, 13 suggest that $`f_\mathrm{R}^{(0)}(x)`$ and $`f_\mathrm{I}^{(0)}(x)`$ are given at large $`|x|`$ by

$`f_\mathrm{R}^{(0)}(x)`$ $``$ $`{\displaystyle \frac{4}{x\mu }}`$ (A.4)
$`f_\mathrm{I}^{(0)}(x)`$ $``$ $`{\displaystyle \frac{4}{x}},`$ (A.5)
which agrees with the above argument.