{'fraction_non_alphanumeric': 0.03921096944912196, 'fraction_numerical': 0.007216742843396681, 'mean_word_length': 4.121971252566735, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Transforming Signs to Phase Distributions in Quantum Simulations

A method is developed which speeds up averaging in quantum simulations where minus signs cause difficulties. A Langevin equation method in conjunction with a replication algorithm is used enabling one to average over a continuously varying complex number. Instead of ensemble averaging this number directly, the phase of the complex number is followed over time. The method is illustrated in some simple cases where the answers obtained can be compared to exact results, and also compared to conventional averaging procedures which converge orders of magnitude slower than this method. Limitations of this method are also described.


I. INTRODUCTION
Simulation of quantum many body systems have been greatly hindered by the "sign problem". Quantum systems with many degrees of freedom are often simulated using stochastic methods. For a variety of problems, usually involving fermions, an average over many negative and positive numbers must be taken. Sometimes the sign does not effect results [1] but in many interesting cases the resulting average is very small because of cancellations, so that in order to obtain good statistics, the number of realizations that must be taken is exponential in the inverse temperature, making it infeasible to get information about ground state properties of a system.
The main idea of this work can be understood by considering the following example. It will be seen below that it describes the quantum mechanics of n bosons restricted to one site. Consider the stochastic equation for z(t) which is a complex function of time:ż = i f (t)z (1) Here f (t) is real gaussian noise with correlation function f (t)f (t ′ ) = vδ(t − t ′ ), where the angled brackets denote an average over the noise f (t). Choosing z(0) = 1 the above equation is easily solved giving where Here θ(t) describes a random walk in time, and therefore the the probability distribution of θ at time t, P (θ, t), is a gaussian distribution with zero mean and variance σ 2 = vt. Below we will see that the ground state energy of n bosons is simply obtained from the asymptotic behavior of z(t) n , which is solved for as follows: Now let us turn to the problem of how, without the aid of this analytical solution, one could obtain z(t) n numerically. Because |z(t)| = 1, to obtain the correct average would require over e n 2 vt independent runs of eqn. (1) to obtain adequate statistics. This can be very large.
The alternative to this direct averaging procedure is to obtain the probability distribution P (θ, t). In this case it is much better to keep track of the total angle θ, not θ mod (2π). Keeping track of the total angle is done as eqn. (1) is evolved. The increment in θ found in a single time step is added to the previous value of θ and a histogram of θ at a given time can be obtained by running eqn. (1) over many realizations of f .
Does analytically continuing θ and then finding P (θ, t) help in obtaining z(t) n ? In this case and in the less trivial cases considered below, it appears to reduce the computation time from something exponential to algebraic. Of course assumptions must be made in order to obtain such a reduction, the most important one being that P (θ, t) is a smooth function of θ. There are also problems encountered in deciding what function to use in fitting P , as will be discussed later. In the case of interest here, fitting ln P (θ, t) to a quadratic function of θ gives accurate results. In fact the ground state energies should be obtained to within an error N −1/2 where N is the number of independent runs.
Having motivated the method that will be considered, the Langevin equation that simulates quantum many body systems will be described and numerical results will be given.

II. SIMULATION METHOD
A system of fermions or bosons can be simulated by considering a stochastic equation which is a generalization of eqn. (1). A field φ(r, t) is governed by the equatioṅ . This describes a system of bosons of mass m interacting with potential evolving in imaginary time. The wave function at time t, is To simulate n fermions, n fields φ 1 , φ 2 , . . . , φ n are evolved starting from different initial conditions but with the same f (r, t). The wave function at time t in this case can be expressed in terms of the average of a Slater determinant: Note that the potential in eqn. (6) is the negative of the correlation function. Therefore to simulate repulsive potentials requires an f (r, t) that is complex [2]. This method is closely related to a simulation method for the density matrix [3,4] and also to projector Monte-Carlo [5].
In order to obtain the ground state energy for quantum systems such as the ones described above requires determining the asymptotic exponential decay of the wave function. In principal this can be done by averaging over many realizations of f , however there are two reasons why this is impractical.
The first is a problem of importance sampling. As an example, consider determining the ground state energy E(n) of n bosons. This is obtained by calculating < φ n (r, t) > and fitting this to an exponential exp(E(n)t) for long times.
The problem is that E(1) = E(n)/n. As a result the typical behavior of φ n (r, t) and the average behavior are very different. What dominates the average are some very rare realizations of f . Another way of describing such behavior is in terms of the notion of intermittency [6]. In order to get the correct average, one must use importance sampling.
There are two ways to do this. One is to use Monte-Carlo keeping track of all time paths [3]. The second is to use a replication algorithm [7,8]. In order to keep track of angles, it is much more convenient to use the second method.
A large number of copies of systems described by eqn. (5) are run in parallel. The number of configurations that are replicated is proportional to their weight. For example, if we wish to measure |φ| n+1 / |φ| n at some site r, we should choose a weight proportional to |φ| n . The proportionality factor is chosen so that the number of copies of the system stays almost constant. This method is similar to others described in detail in refs. [7] and [8].
The second problem is one encountered when simulating repulsive bosons and fermions. Here non-positivity leads to the problem discussed in the introduction, that an exponentially large number of runs must be taken to obtain adequate statistics. The phase method proposed here is designed to reduce the number of runs. This will be illustrated numerically for a system of repulsive bosons. One should note that although it is possible to choose a method for this problem where there are no minus signs [9,10], this is an instructive example. It will be seen that the method used generalizes to fermions.

III. OBTAINING THE PHASE
Consider n bosons on a one dimensional lattice with L sites interacting via an on-site repulsive interaction U . This has been the focus of recent investigation [11] using "world line" Monte Carlo [9,10], which is clearly better suited to this problem than the method here. We use the method above only to illustrate techniques to cope with negative signs. In obtaining ground state energies, any component of |Ψ > can be measured. As will be seen below, in order to make the phase method work it is important to make a judicious choice for this. We start by choosing to measure Ψ(x 1 , x 2 , . . . , x n ) where the x i 's are all different sites. Later we will see why using the same sites does not work as well. The prescription for computing Ψ is given by eqn. (7).
The replication method outlined above was the case where all the x i 's are the same. To generalize the discussion above, one just chooses a weight w ≡ |φ(x 1 )φ(x 2 ) . . . φ(x n )|. The ground state energy of n + 1 particles can be easily computed by knowing where the last average is done with respect to the weight w, and θ is the total phase angle of φ(x 1 ) . . . φ(x n+1 ).
To test this out numerically, the number of lattice sites was chosen to be 8, and we started out by using a small value for U = 0.5. An average of 4497 copies were replicated. From this data P (Θ, t) was obtained with good statistics. In order to obtain an improved estimate of the tails of this distribution, the replication algorithm was altered to weight large Θ, that is by letting w → w exp(const.Θ). By changing the constant in the exponent, P (Θ, t) was probed for different regions of Θ. By combining these, one obtains a good estimate for P (Θ, t) over nine orders of magnitude. An example of ln P (Θ, t) is shown in fig. 1 for n = 8. The solid line is a quadratic fit to the data.
By taking the fourier transform of this for different times, as prescribed by eqn. (9), the ground state energy can be obtained. This is shown in fig. 2   Without using P (Θ, t), one can compute the weighted average of |φ(x n+1 )| exp(iΘ) directly. For comparison the average of this, over the same copies, is shown fig. 3. It is clear that results obtained by direct averaging are far inferior to those obtained using P (Θ, t).
FIG. 3. The result of averaging the field over the same runs as was used to obtain fig. 1. An exponential decay is not discernable and many more averages would be needed to see it. The real part is plotted here.

IV. IMPROVEMENTS TO THE METHOD
The method as it stands has difficulties for large U . The histogram of Θ has too large a variance. This can be understood by considering a variant of eqn. (1) in which there are two independent random processes with different and therefore z i = e iθi(t) , i = 1, 2 . Suppose we are interested in numerically computing the sum This average will be dominated for long times by the z i with the smaller of the v i 's, say v 1 . Under what circumstances will the numerical method of finding the histogram for the phase angle of the above sum give a sensible answer? Eqn.
(10) can be thought of as the sum of two randomly rotating vectors with magnitudes |A i |. From this observation, it is easy to see that for long times, the variance of the phase angle will be the same as that for z n 1 if |A 1 | > |A 2 |. For |A 1 | < |A 2 | a broader distribution will be obtained. Unless very careful fitting to a non-quadratic expression is done, an erroneously high energy will be obtained.
This example shows that very careful fitting is necessary if the wave function has too high an amplitude of excited states. One way to reduce this effect is to choose to measure Ψ on different sites as was done in the numerical example above. Another way to improve this method is to filter φ in time. It is easily shown that applying a linear filter to φ(t), γ(t) ≡ φ(t) * g(t) does not change the ground state energy if the filtering function g(t) decays sufficiently rapidly.
For example one can evolve the equationγ (r, t) = − cγ(r, t) + φ(r, t) simultaneously with evolving eqn. (5) and measure γ n instead of φ n . This method works quite well. The energies for an eight site lattice for the case U = 4 are shown in fig. 4. The energies are in good agreement with the exact results up to n = 6 after which it becomes apparent that the higher energy states start to dominate the answer.
The method will presumably improve if in addition to filtering, different sites are used as in the example above. Lastly, the case of fermions has been considered. The determinant formulation of this problem can be used to obtain a continuous phase as a function of time. Fig. (5) shows the complex determinant d on a log polar plot, that is with radius − ln |d(t)| and angle arg(d(t)), for 3 particles on a 4 site chain, with U = 2.5. The phase randomly rotates, similar to the case of bosons.
In conclusion, a new method for dealing with cancellations in quantum systems has been developed and tried out in some simple cases. Further theoretical development of the determination of P (Θ, t) could greatly improve fitting.
Recent theoretical work [12] relating the phase used here to the Berry phase for smooth paths should be further explored in the context of this simulation method. It would also be interesting to consider more difficult problems such as the two dimensional Hubbard model with repulsive interactions using this method.

ACKNOWLEDGMENT
The author thanks Eliot Dresselhaus, Richard Scalettar, and Peter Young for useful discussions. This research was supported by the NSF under grant DMR-9112767.