text:
# Incorporation of Density Matrix Wavefunctions in Monte Carlo Simulations: Application to the Frustrated Heisenberg Model

## I Introduction

The Density Matrix Technique (DMRG) has proven to be a very efficient method to determine the groundstate properties of low dimensional systems . For a quantum chain it produces extremely accurate values for the energy and the correlation functions. In two dimensional systems the calculational effort increases rapidly with the size of the system. The most favorable geometry is that of a long small strip. In practice the width of the strip is limited to around 8 to 10 lattice sites. Greens Function Monte Carlo (GFMC) is not directly limited by the size of the system but by the efficiency of the importance sampling. When the system has a minus sign problem the statistics is ruined in the long run and accurate estimates are impossible. Many proposals  have been made to alleviate or avoid the minus sign problem with varying success, but all of them introduce uncontrollable errors in the sampling. In the DMRG calculation of the wavefunction the minus sign problem is not manifestly present. In all proposed cures of the minus sign problem the errors decrease when the guiding wavefunction approaches the groundstate.

The idea of this paper is that DMRG wavefunctions are much better, also for larger systems, than the educated guesses which usually feature as guiding wave functions. Moreover DMRG is a general technique to construct a wavefunction without knowing too much about the nature of the groundstate, with the possibility to systematically increase the accuracy. Thus DMRG wavefunctions would do very well when they could be used as guiding functions in the importance sampling of the GFMC. There is a complicating factor which prevents a straightforward implementation of this idea due to the fact that interesting systems are so large that it is impossible to use a wavefunction via a look–up table. The value of the wavefunction in a configuration has to be calculated by an in–line algorithm. This has limited the guiding wavefunctions to simple expressions which are fast to evaluate. Consequently such guiding wavefunctions are not an accurate representation of the true groundstate wavefunction, in particular if the physics of the groundstate is not well understood. In this paper we describe a method to read out the DMRG wavefunction in an efficient way by using a special representation of the DMRG wavefunction.

A second problem is that a good guiding wavefunction alleviates the minus sign problem, but cannot remove it as long as it is not exact. We resolve this dilemma by applying the method of Stochastic Reconfiguration which has recently been proposed by Sorella . The viability of our method is tested for the frustrated Heisenberg model.

The behavior of the Heisenberg antiferromagnet has been intruiging for a long time and still is in the center of research. The groundstate of the antiferromagnetic 1-dimensional chain with nearest neighbor coupling is exactly known. In higher dimensions only approximate theories or simulation results are available. The source of the complexity of the groundstate are the large quantum fluctuations which counteract the tendency of classical ordering. The unfrustrated 2–dimensional Heisenberg antiferromagnet orders in a Néel state and by numerical methods the properties of this state can be analyzed accurately . The situation is worse when the interactions are competing as in a 2-dimensional square lattice with antiferromagnetic nearest neighbor $`J_1`$ and next nearest neighbor $`J_2`$ coupling. This spin system with continuous symmetry can order in 2 dimensions at zero temperature, but it is clear that the magnetic order is frustrated by the opposing tendencies of the two types of interaction. The ratio $`J_2/J_1`$ is a convenient parameter for the frustration. For small values the system orders antiferromagnetically in a Néel type arrangement, which accomodates the nearest neighbor interaction. For large ratios a magnetic order in alternating columns of aligned spins (columnar phase) will prevail; in this regime the roles of the two couplings are reversed: the nearest neighbor interaction frustrates the order imposed by the next nearest neighbor interaction. In between, for ratios of the order of 0.5, the frustration is maximal and it is not clear which sort of groundstate results. This problem has been attacked by various methods but not yet by DMRG and only very recently by GFMC . This paper addresses the issue by studying the spin correlations.

A simple road to the answer is not possible since the behavior of the system with frustration presents some fundamental problems. The most severe obstacle is that frustration implies a sign problem which prevents the straightforward use of the GFMC simulation technique. Moreover the frustration substantially complicates the structure of the groundstate wavefunction. Generally frustration encourages the formation of local structures such as dimers and plaquettes which are at odds, but not incompatible, with long range magnetic order. These correlation patterns are the most interesting part of the intermediate phase and the main goal of this investigation.

Many attempts have been made to clarify the situation. Often simple approximations such as mean–field or spin–wave theory give useful information about the qualitative behavior of the phase diagram. A fairly sophistocated mean–field theory using the Schwinger boson representation does not give an intermediate phase . Given the complexity of the phase diagram and the subtlety of the effects it is not clear whether such approximate methods can give in this case a reliable clue to the qualitative behavior of the system.

Exact calculations have been performed on small systems up to size $`6\times 6`$ by Schulz et al. . Although this information is very accurate and unbiased to possible phases, the extrapolation to larger systems is a long way, the more so in view of indications that the anticipated finite size behavior only applies for larger systems. Another drawback of these small systems is that the groundstate is assumed to have the full symmetry of the lattice. Therefore the symmetry breaking, associated with the formation of dimers, ladders or plaquettes, which is typical for the intermediate state, can not be observed directly.

More convincing are the systematic series expansion as reported recently by Kotov et al. , and by Singh et al. , which bear on an infinite system. They start with an independent dimers (plaquettes) and study the series expansion in the coupling between the dimers (plaquettes). By the choice of the state, around which the perturbation expansion is made, the type of spatial symmetry breaking is fixed. These studies favor in the intermediate regime the dimer state over the plaquette state. Their dimer state has dimers organized in ladders in which the chains and the rungs have nearly equal strength. So the system breaks the translational invariance only in one direction. The energy differences are however small and the series is finite, so further investigation is useful. Our simulations yield correlations in good agreement with theirs, but do not confirm the picture of translational invariant ladders. Instead we find an additional weaker symmetry breaking along the ladders, such that we come closer to the plaquette picture.

Very recently Capriotti and Sorella  have carried out a GFMC simulation for $`J_2=0.5J_1`$ and have studied the susceptibilities for the orientational and translational symmetry breaking. They conclude that the groundstate is a plaquette state with full symmetry between the horizontal and vertical direction.

From the purely theoretical side the problem has been discussed by Sachdev and Read  on the basis of a large spin expansion. From their analysis a scenario emerges in which the Néel phase disappears upon increasing frustration in a continuous way. Then a gapped spatial–inhomogeneous phase with dimerlike correlations appears. For even higher frustration ratios a first order transition takes place to the columnar phase. Although this scenario is qualitative, without precise location of the phase transition points, it definitively excludes dimer formation in the magnetically ordered Néel and columnar phase. It is remarkable that two quite different order parameters (the magnetic order and the dimer order) disappear simultaneously and continuously on opposite sides of the phase transition. In this scenario, this is taken as an indication of some kind of duality of the two phases.

Given all these predictions it is of utmost interest to further study the nature of the intermediate state. Due to the smallness of the differences in energy between the various possibilities, the energy will not be the ideal test for the phase diagram. Therefore we have decided to focus directly on the spin correlations as a function of the ratio $`J_2/J_1`$. In this paper we first investigate the 2–dimensional frustrated Heisenberg model by constructing the DMRG wave function of the groundstate for long strips up to a width of 8 sites. The groundstate energy and the spin stiffness which are calculated, confirm the overal picture described above, but the results are not accurate enough to allow for a conclusive extrapolation to larger systems. Then we study an open 10x10 lattice by means of the GFMC technique using DMRG wavefunctions as guiding wavefunction for the importance sampling. The GFMC are supplemented with Stochastic Reconfiguration as proposed by Sorella  as an extension of the Fixed Node technique . This method avoids the minus sign problem by replacing the walkers regularly by a new set of positive sign with the same statistical properties. The first observation is that GFMC improves the energy of the DMRG in a substantial and systematic way as can be tested in the unfrustrated model where sufficient information is available from different sources. Secondly the spin correlations become more accurate and less dependent on the technique used for constructing the DMRG wavefunction. The DMRG technique is focussed on the energy of the system and less on the correlations. The GFMC probes mostly the local correlations of the system as all the moves are small and correspond to local changes of the configurations. With these spin correlations we investigate the phase diagram for various values of the frustration ratio $`J_2/J_1`$.

The paper begins with the definition of the model to avoid ambiguities. Then a short description of our implementation of the DMRG method is given. We go into more detail about the way how the constructed wavefunctions can be used as guiding wavefunctions in the GFMC simulation. This is a delicate problem since the full construction of a DMRG wavefunction takes several hours on a workstation. Therefore we separate off the construction of the wavefunction and cast it in a form where the configurations can be obtained from each other by matrix operations on a vector. So the length of the computation of the wavefunction in a configuration scales with the square of the number of states included in the DMRG wavefunction. But even then the actual construction of the value of the wavefunction in a given configuration is so time consuming that utmost effiency must be reached in obtaining the wavefunction for successive configurations. The remaining sections are used to outline the GFMC and the Stochastic Reconfiguration and to discuss the results. We concentrate on the correlation functions since we see them as most significant for the structure of the phases. We give first a global evaluation of the correlation function patterns for a wide set of frustration ratios and then focus on a number of points to see the dependence on the guiding wavefunction and to deduce the trends. The paper closes with a discussion and a comparison with other results in the literature.

## II The Hamiltonian

The hamiltonian of the system refers to spins on a square lattice.

$$=J_1\underset{(i,j)}{}𝐒_i𝐒_j+J_2\underset{[i,j]}{}𝐒_i𝐒_j.$$
(1)
The $`𝐒_i`$ are spin $`\frac{1}{2}`$ operators and the sum is over pairs of nearest neigbors $`(i,j)`$ and over pairs of next nearest neighbors $`[i,j]`$ on a quadratic lattice. Both coupling constants $`J_1`$ and $`J_2`$ are supposed to be positive. $`J_1`$ tries to align the nearest neigbor spin in an antiferromagnetic way and $`J_2`$ tries to do the same with the next nearest neighbors. So the spin system is frustrated, implying an intrinsic minus sign in the simulations that cannot be gauged away by a rotation of the spin operators.

In order to prepare for the representation of the hamiltonian we express the spin components in spin raising and lowering operators

$$𝐒_i𝐒_j=\frac{1}{2}(S_i^+S_j^{}+S_i^{}S_j^+)+S_i^zS_j^z.$$
(2)
We will use the $`z`$ component representation of the spins and a complete state of the spins will be represented as

$$|R=|s_1,s_2,\mathrm{},s_N,$$
(3)
where the $`s_j`$ are eigenvalues of the $`S_j^z`$ operator. The diagonal matrix elements of the hamiltonian are in the representation (3) given by

$$R||R=J_1\underset{(i,j)}{}s_is_j+J_2\underset{[i,j]}{}s_is_j.$$
(4)
The off-diagonal elements are between two nearby configurations $`R^{}`$ and $`R`$. $`R^{}`$ is the same as $`R`$ except at a pair of nearest neighbors sites $`(i,j)`$ or next nearest neighbor sites $`[i,j]`$, for which the spins $`s_i`$ and $`s_j`$ are opposite. In $`R^{}`$ the pair is turned over by the hamiltonian. Then

$$R^{}||R=\frac{1}{2}J_1\mathrm{or}R^{}||R=\frac{1}{2}J_2,$$
(5)
depending on whether a nearest or a next nearest pair is flipped.

## III The DMRG Procedure

The DMRG procedure approximates the groundstate wavefunction by searching through various representations in bases of a given dimension $`m`$ . Here we take the standard method (with two connecting sites) for granted and make the preparations for the extraction of the wavefunction.

The system is mapped onto a 1-dimensional chain (see Fig. 2) and separated into two parts: a left and right hand part. They are connected by one site. Each part is represented in a basis of at most $`m`$ states. With a representation of all the operators in the hamiltonian in these bases one can find the groundstate of the system. We thus have several representations of the groundstate depending on the way in which the system is divided up into subsystems. The point is to see how these representations are connected and how they possibly can be improved. We take a representation for the right hand parts and improve those on the left. So we assume that for a given division we have the groundstate of the whole system and we want to enlarge the left hand side at the expense of the right hand side. The first step is to include the connecting site in the left hand part. This enlarges the basis for the left hand side from $`m`$ to $`2m`$ and a selection has to be made of $`m`$ basis states. This goes with the help of the density matrix for the left hand side as induced by the wavefunction for the whole system. For later use we write out the basic equations for the density matrix in the configuration representation. Let, at a certain stage in the computation, $`|\mathrm{\Phi }`$ be the approximation to the groundstate. The configurations of the right hand part and the left hand part are denoted by $`R_r`$ and $`R_l`$. Then the density matrix for the left hand part reads

$$R_l|\rho |R_l^{}=\underset{R_r}{}R_l,R_r|\mathrm{\Phi }\mathrm{\Phi }|R_l^{},R_r.$$
(6)
In practice we do not solve the eigenvalues of the density matrix in the configuration representation, but in a projection on a smaller basis. White  has shown that the best way to represent the state $`|\mathrm{\Phi }`$ is to select the $`m`$ eigenstates $`|\alpha `$ with the largest eigenvalue

$$\underset{R_l^{}}{}R_l|\rho |R_l^{}R_l^{}|\alpha =\lambda _\alpha R_l|\alpha .$$
(7)
The next step is to break up the right hand part into a connecting site and a remainder. With the basis for this remainder and the newly acquired basis for the left hand part we can again compute the groundstate of the whole system as indicated in the lower part of the figure. Now we are in the same position as we started, with the difference that the connecting site has moved one position to the right. Thus we may repeat the cycle till the right hand part is so small that it can exactly be represented by $`m`$ states or less. Then we have constructed for the left hand part a new set of bases, all containing $`m`$ states, for system parts of variable length. Next we reverse the roles of left and right and move back in order to improve the bases for the right hand parts with the just constructed bases for the left hand part.

The process may be iterated till it converges towards a steady state. The great virtue of the method is that it is variational. In each step the energy will lower till it saturates. In 1-dimensional system the method has proven to be very accurate . So one wonders what the main trouble is in higher dimensions.

In Fig. 3 we have drawn 2 possible ways to map the system on a 1-dimensional chain. One sees that if we divide again the chain into a left hand part and a right hand part and a connecting site, quite a few sites of the left hand part are nearest or next nearest neighbors of sites of the right hand part. So the coupling between the two parts of the chain is not only through the connecting site but also through sites which are relatively far away from each other in the 1-dimensional path. The operators for the spins on these sites are not as well represented as those of the connecting site, which is fully represented by the two possible spin states. Yet the correlations between the interacting sites count as much for the energy of the system as those interacting with the connecting site. One may say that the further away two interacting sites are in the 1-dimensional chain the poorer their influence is accounted for. This consideration explains in part why open systems can be calculated more accurately than closed systems, even in 1-dimensional systems.

It is an open question which map of the 2-dimensional onto a 1-dimensional chain gives the best representation of the groundstate of the system. Also other divisions of the system than those suggested by a map on a 1–dimensional chain are possible and we have been experimenting with arrangements which reflect better the 2–dimensional character of the lattice. They are promising but the software for these is not as sophisticated as the one developed by White  for the 1-dimensional chain. We therefore have restricted our calculations to the two paths shown here. The second choice. the “meandering” path, was motivated by the fact that it has the strongest correlated sites most nearby in the chain and this choice was indeed justified by a lower energy for a given dimension $`m`$ of the representation than for the “straight” path.

The DMRG calculations as well as the corresponding GFMC simulations are carried out for both paths. The meandering path has to be preferred over the straight path as the DMRG wavefunctions generally give a better energy value and the simulations suffer less from fluctuations. Nevertheless we have also investigated the straight path, since the path chosen leaves its imprints on the resulting correlation pattern and the paths break the symmetries in different ways. Both paths have an orientational preference. In open systems the translational symmetry is broken anyway, but the meandering path has in addition a staggering in the horizontal direction. This together with the horizontal nearest neighbor sites appearing in the meandering path gives a preference for horizontal dimerlike correlations in this path. On the other hand the straight path prefers the dimers in the vertical direction. Comparing the results of the two choices, allows us to draw further conclusions on the nature of the intermediate state.

## IV Extracting configurations from the DMRG wavefunction

It is clear that the wavefunction which results from a DMRG-procedure is quite involved and it is not simple to extract its value for a given configuration. We assume now that the DMRG wavefunction has been obtained by some procedure and we will give below an algorithm to obtain efficiently the value for an arbitrary configuration (see also  for an alternative description).

The first step is the construction of a set of representations for the wavefunction in terms of two parts (without a connecting site in between). Let the left hand part contain $`l`$ sites and the other part $`Nl`$ sites. We denote the $`m`$ basis states of the left hand part by the index $`\alpha `$ and those of the right hand part by $`\overline{\alpha }`$. The eigenstates of the two parts are closely linked and related as follows

$$\{\begin{array}{ccc}\hfill R_l|\alpha & =& \frac{1}{\sqrt{\lambda _\alpha }}\underset{R_r}{}\mathrm{\Phi }|R_l,R_rR_r|\overline{\alpha },\hfill \\ \hfill R_r|\overline{\alpha }& =& \frac{1}{\sqrt{\lambda _\alpha }}\underset{R_l}{}\alpha |R_lR_l,R_r|\mathrm{\Phi }.\hfill \end{array}$$
(8)
It means that for every eigenvalue $`\lambda _\alpha `$ there is and eigenstate $`\alpha `$ for the left hand part and an $`\overline{\alpha }`$ for the right hand density matrix. The proof of (8) follows from insertion in the density matrix eigenvalue equation (7).

The second step is a relation for the groundstate wavefunction in terms of these eigenfunctions. Generally we have

$$R_l,R_r|\mathrm{\Phi }=\underset{\alpha ,\overline{\beta }}{}R_l|\alpha R_r|\overline{\beta }\alpha \overline{\beta }|\mathrm{\Phi },$$
(9)
while due to (8) we find

$$\begin{array}{ccc}\hfill \alpha \overline{\beta }|\mathrm{\Phi }& =& \underset{R_l,R_r}{}\alpha |R_l\overline{\beta }|R_rR_l,R_r|\mathrm{\Phi }\hfill \\ & =& \sqrt{\lambda _\alpha }\underset{R_r}{}\overline{\beta }|R_rR_r|\overline{\alpha }=\delta _{\alpha ,\beta }\sqrt{\lambda _\alpha }.\hfill \end{array}$$
(10)
Thus we can represent the groundstate as

$$R_l,R_r|\mathrm{\Phi }=\underset{\alpha }{}\sqrt{\lambda _\alpha ^l}R_l|\alpha _lR_r|\overline{\alpha }_{Nl}.$$
(11)
For this part of the problem we have to compute and store the set of $`m`$ eigenvalues $`\lambda _\alpha ^l`$ for each division $`l`$. We point out again that we have on the left hand side the wavefunction and on the right hand side representations for given division $`l`$, which all lead to the same wavefunction. The last step is to see the connection between these representations.

As intermediary we consider a representation of the wavefunction with one site $`s_l`$ separating the spins $`s_1\mathrm{}s_{l1}`$ on the left hand side from $`s_{l+1}\mathrm{}s_N`$ on the right hand side. Using the same basis as in (11) we have

$$\begin{array}{c}s_1\mathrm{}s_{l1},s_l,s_{l+1}\mathrm{}s_N|\mathrm{\Phi }=\hfill \\ _{\alpha ,\alpha ^{}}s_1\mathrm{}s_{l1}|\alpha \varphi _{\alpha ,\alpha ^{}}^l(s_l)s_{l+1}\mathrm{}s_N|\overline{\alpha ^{}}.\hfill \end{array}$$
(12)
We compare this representation in two ways with (11). First we combine the middle site with the left hand part. This leads to $`m`$ states which can be expressed as linear combinations of the states of the enlarged segment

$$\underset{\alpha }{}s_1\mathrm{}s_{l1}|\alpha \varphi _{\alpha ,\alpha ^{}}^l(s_l)=\underset{\alpha ^{\prime \prime }}{}s_1\mathrm{}s_l|\alpha ^{\prime \prime }T_{\alpha ^{\prime \prime },\alpha ^{}}^l.$$
(13)
In fact this relation is the very essence of the DMRG procedure. The wave function in the larger space is projected on the eigenstates of the the density matrix of that space. Since the process of zipping back forth has converged there is indeed a fixed relation (13). However when we insert (13) into (12) and compare it with (11) we conclude that the matrix $`T`$ must be diagonal

$$T_{\alpha ^{\prime \prime },\alpha ^{}}^l=\delta _{\alpha ^{\prime \prime },\alpha ^{}}\sqrt{\lambda _\alpha ^{}^l}.$$
(14)
This leads to the recursion relation

$$s_1\mathrm{}s_l|\alpha ^{}=\underset{\alpha }{}s_1\mathrm{}s_{l1}|\alpha A_{\alpha ,\alpha ^{}}^l(s_l)$$
(15)
with

$$A_{\alpha ,\alpha ^{}}^l(s_l)=\varphi _{\alpha ,\alpha ^{}}^l(s_l)/\sqrt{\lambda _\alpha ^{}^l}.$$
(16)
The second combination concerns the contraction of the middle site with the right hand part. This leads to the recursion relation

$$s_l\mathrm{}s_N|\overline{\alpha }=\underset{\alpha ^{}}{}B_{\alpha ,\alpha ^{}}^{l1}(s_l)s_{l+1}\mathrm{}s_N|\overline{\alpha }^{}$$
(17)
with

$$B_{\alpha ,\alpha ^{}}^{l1}(s)=\varphi _{\alpha ,\alpha ^{}}^l(s)/\sqrt{\lambda _\alpha ^{l1}}.$$
(18)
The $`A`$ and $`B`$ matrices are the essential ingredients of the calculation of the wavefunction. From (18) and (16) follows that they are related as

$$B_{\alpha ,\alpha ^{}}^{l1}(s)=\sqrt{\lambda _\alpha ^{}^l/\lambda _\alpha ^{l1}}A_{\alpha ,\alpha ^{}}^l(s).$$
(19)
By the recursion relations the basis states are expressed as products of $`m\times m`$ matrices. The determination of the DMRG wavefunction and the matrices $`A`$ (or $`B`$) is part of the determination of the DMRG wavefunction which is indeed lengthy but fortunately no part of the simulation. The matrices can be stored and contain the information to calculate the wavefunction for any configuration. The value of the wavefunction is now obtained as the product of matrices acting on a vector. Thus the calculational effort scales with $`m^2`$. Using relation (19) one reconfirms by direct calculation that the wavefunction is indeed independent of the division $`l`$.

When the simulation is in the configuration $`R`$, all the $`R_l|\alpha _l`$ and the $`R_r|\overline{\alpha }_{Nl}`$ are calculated and stored, with the purpose to calculate the wavefunctions more efficiently for the configurations $`R^{}`$ which are connected to $`R`$ by the hamiltonian and which are the candidates for a move. The structure of of these nearby states is $`R^{}=s_1\mathrm{}s_{l_2}\mathrm{}s_{l_1}\mathrm{}s_N(l_2>l_1)`$. So we have that for $`R^{}`$ the representation

$$R^{}|\mathrm{\Phi }=\underset{\alpha }{}\sqrt{\lambda _\alpha ^{l_2}}s_1\mathrm{}s_{l_2}\mathrm{}s_{l_1}|\alpha s_{l_2+1}\mathrm{}s_N|\overline{\alpha }$$
(20)
holds. Now we see the advantage of having the wavefunction stored for all the divisions. The second factor in (20) is already tabulated; the first factor involves a number of matrix multiplications equal to the distance in the chain of the two spins $`l_1`$ and $`l_2`$ till one reaches a tabulated function. One can use the tables for a certain number of moves but after a while it starts to pay off to make a fresh list.

## V Green Function Monte Carlo simulations

The GFMC technique employs the operator

$$𝒢=1ϵ$$
(21)
and uses the fact that the groundstate $`|\mathrm{\Psi }_0`$ results from

$$|\mathrm{\Psi }_0𝒢^n|\mathrm{\Phi },ϵ1,nϵ1$$
(22)
where in principle $`|\mathrm{\Phi }`$ may be any function which is non-orthogonal to the groundstate. In view of the possible symmetry breaking, the overlap is a point of serious concern on which we come back in the discussion. In practice we will use the best $`|\mathrm{\Phi }`$ that we can construct conveniently by the DMRG–procedure described above. The closer $`|\mathrm{\Phi }`$ is to the groundstate the smaller the number of factors $`n`$ in the product needs to be in order to find the groundstate. Evaluating (22) in the spin representation gives for the projection on the trial wavefunction the following long product

$$\mathrm{\Phi }|\mathrm{\Psi }_0\underset{𝐑}{}\mathrm{\Phi }|R_M\left[\underset{i=1}{\overset{M}{}}R_i|𝒢|R_{i1}\right]R_0|\mathrm{\Phi }.$$
(23)
Here the sum is over paths $`𝐑=(R_M,\mathrm{}R_1,R_0)`$ which will be generated by a Markov process. The Markov process involves a transition probability $`T(R_iR_{i1})`$ and the averaging process uses a weight $`m(R)`$. Its is natural to connect the transition probabilities to the matrix elements of the Greens Function $`𝒢`$. But here comes the sign problem into the game: the transition probabilities have to be positive (and normalized). So we put the transition rate proportional to the absolute value of the matrix element of the Greens Function

$$T(RR^{})=\frac{|R|𝒢|R^{}|}{_{R^{\prime \prime }}|R^{\prime \prime }|𝒢|R^{}|}.$$
(24)
This implies that we have to use a sign function $`s(R,R^{})`$

$$s(R,R^{})=\frac{R|𝒢|R^{}}{|R|𝒢|R^{}|}$$
(25)
and a weight factor

$$m(R)=\underset{R^{}}{}|R^{}|𝒢|R|.$$
(26)
All these factors together form the matrixelement of the Green Function

$$R|𝒢|R^{}=T(RR^{})s(R,R^{})m(R^{}).$$
(27)
If the matrix elements of the Greens Function were all positive, or could be made positive by a suitable transformation, we would not have to introduce the sign function. We leave its consequences to the next section. By the representation (27) we can write the contribution of the path as a product of transition probabilities, signs and local weights. The transition probabilities control the growth of the Markov chain. The signs and weights constitute the weight of a path

$$M(𝐑)=m_f(R_M)\left[\underset{i=1}{\overset{M}{}}s(R_i,R_{i1})m(R_{i1})\right]m_i(R_0).$$
(28)
The initial and final weight have to be chosen such that the weight of the paths corresponds to the expansion (23). For the innerproduct $`\mathrm{\Phi }|\mathrm{\Psi }_0`$ we get

$$m_i(R)=R|\mathrm{\Phi },m_f(R)=\mathrm{\Phi }|R.$$
(29)
With this final weight we have projected the groundstate on the trial wave. This allows us to calculate the so–called mixed averages. For that purpose we define the local estimator $`𝒪`$

$$O(R)=\frac{\mathrm{\Phi }|𝒪|R}{\mathrm{\Phi }|R},$$
(30)
which yields the mixed average

$$𝒪_m\frac{\mathrm{\Phi }|𝒪|\mathrm{\Psi }_0}{\mathrm{\Phi }|\mathrm{\Psi }_0}=\frac{_𝐑O(R_M)M(𝐑)}{_𝐑M(𝐑)}.$$
(31)
For operators not commuting with the hamiltonian the mixed average is an approximation to the groundstate average. Later on we will improve on it.

In this raw form the GFMC would hardly work because all paths are generated with equal weight. One can do better by importance sampling in which one transforms the problem to a Greens Function with matrix elements

$$R|\stackrel{~}{𝒢}|R^{}=\frac{\mathrm{\Phi }|RR|𝒢|R^{}}{\mathrm{\Phi }|R^{}}.$$
(32)
Generally this can be seen as a similarity transformation on the operators and from now on everywhere operators with a tilde are related to their counterpart without a tilde as in (32). It gives only a minor change in the formulation. The transition rates are based on the matrix elements of $`\stackrel{~}{𝒢}`$ and so are the signs and weights. Thus we have a set of definitions like (24)–(26) with everywhere a tilde on top. It leads also to a change of the initial and final weight

$$\stackrel{~}{m}_i(R)=|\mathrm{\Phi }|R|^2,\stackrel{~}{m}_f(R)=1.$$
(33)
By chosing these weights the formula (30) for the average still applies with a weight $`\stackrel{~}{M}(𝐑)`$ made up as in (28) with the weights and signs with a tilde. Using the tilde operators the local estimator (30) reads

$$O(R)=\underset{R^{}}{}R^{}|\stackrel{~}{𝒪}|R.$$
(34)
We will speak about the various paths in terms of independent walkers that sample these paths. As some walkers become more important than others in the process, it is wise to improve the variance by branching, which we will discuss later with the sign problem. Before we embark on the discussion of the sign problem we want to summarize a number of aspects of the GFMC simulation relevant to our work.

* The steps in the Markov process are small, only the ones induced by one term of the hamiltonian feature in a transition to a new state. This makes the subsequent states quite correlated. So many steps have to be performed before a statistical independent configuration is reached; on the average a number of the order of the number of sites.
* In every configuration the wave function for a number of neighboring states (the ones which are reachable by the Greens Function), has to be evaluated. This is a time consuming operation and it makes the simulation quite slow, because out of the possibilities (of order $`N`$) only one is chosen and all the information gathered on the others is virtually useless.
* The necessity to choose a small $`ϵ`$ in the Greens Function seems a further slow down of the method, but it can be avoided by the technique of continuous time steps developed by Ceperley and Trivedi . In this method the possibility of staying in the same configuration (the diagonal element of the Greens Function) is eliminated and replaced by a waiting time before a move to another state is made. (For further details in relation to the present paper we refer to )
* The average (30) can be improved by replacing it by

$$𝒪_{im}=2\frac{\mathrm{\Phi }|𝒪|\mathrm{\Psi }_0}{\mathrm{\Phi }|\mathrm{\Psi }_0}\frac{\mathrm{\Phi }|𝒪|\mathrm{\Phi }}{\mathrm{\Phi }|\mathrm{\Psi }_0}$$
(35)
of which the error with respect to the true average is of second order in the deviation of $`|\mathrm{\Phi }`$ from $`|\mathrm{\Psi }_0`$. For conserved operators, such as the energy, this correction is not needed since the the mixed average gives already the correct value.

## VI The sign problem and its remedies

In the hamiltonian (1) the $`z`$ component of the spin operator keeps the spin configuration invariant, whereas the $`x`$ and $`y`$ components change the configuration. The typical change is that a pair of nearest or next nearest neighbors is spin reversed. Inspecting the Greens Function it means that all changes to another configuration involve a minus sign! Thus the Greens Function is as far as possible from the ideal of positive matrix elements. The diagonal terms are positive, but they always are positive for sufficiently small $`ϵ`$. Importance sampling can remove minus signs in the transition rates, when the ratio of the guiding wavefunction involves also a minus sign. For $`J_20`$ no guiding wavefunction can remove the minus sign problem completely. In Fig. 3 we show a loop of two nearest neighbor spin flips followed by a flip in a next nearest neighbor pair, such that the starting configuration is restored. The product of the ratios of the guiding wavefunction drops out in this loop, but the product of the three matrix elements has a minus sign. So at least one of the transitions must involve a minus sign.

For unfrustrated systems these loops do not exist and one can remove the minus sign by a transformation of the spin operators

$$S_i^xS_i^x,S_i^yS_i^y,S_i^zS_i^z$$
(36)
which leave the commutation operators invariant. Applying this transformation on every other spin (the white fields of a checkerboard) all flips involving a pair of nearest neighbors then give a positive matrix element for the Greens Function. So when $`J_2=0`$ the appearant sign problem is transformed away. For sufficiently small $`J_2`$, Marshall  has shown that the wave function of the system has only positive components (after the “Marshall” sign flip (36)). So the minus sign problem is not due to the wave function but to the frustration. (For the Hubbard model it is the guiding wave function which must have minus signs due to the Pauli principle, while the bare transition probabilities can be taken positive).

Due to the minus sign the weight of a long path picks up a arbitrary sign. Generally the weights are also growing along a path. Thus if various paths are traced out by a number of independent walkers, the average over the paths or the walkers becomes a sum over large terms of both signs, or differently phrased: the average becomes small with respect to the variance; the signal gets lost in the noise.

Ceperley and Alder  constructed a method, Fixed Node Monte Carlo (FNMC), which avoids the minus sign problem at the expense of introducing an approximation. Their method is designed for continuum systems and handling fermion wavefunctions. They argued that the configuration space in which the wavefunction has a given sign, say positive, is sufficient for exploring the properties of the groundstate, since the other half of the configuration space contains identical information. Thus they designed a method in which the walkers remain in one domain of a given sign, essentially by forbidding to cross the nodes of the wavefunction. The approximation is that one has to take the nodal structure of the guiding wavefunction for granted and one cannot improve on that, at least not without sophistocation (nodal release). The method is variational in the sense that errors in the nodal structure always raise the groundstate energy.

It seems trivial to take over this idea to the lattice but it is not. The reason is that in continuum systems one can make smaller steps when a walker approaches a node without introducing errors. In a lattice system the configuration space is discrete; so the location of the node is not strictly defined. The important part is that, loosely speaking the nodes are between configurations and one cannot make smaller moves than displacing a particle over a lattice distance or flip a pair of spins. Van Bemmel et al.  adapted the FNMC concept to lattice systems preserving its variational character. This extension to the lattice suffers from the same shortcoming as the method of Ceperley and Alder: the “nodal” structure of the guiding wavefunction is given and cannot be improved by the Monte Carlo process. Recently Sorella  proposed a modification which overcomes this drawback. It is based on two ingredients.

Sorella noticed that the following effective hamiltonian yields also an upper bound to the energy:

$$\begin{array}{ccc}\hfill R|\stackrel{~}{}_{\mathrm{eff}}|R^{}& =& \{\begin{array}{ccc}\hfill R|\stackrel{~}{}|R^{}& \mathrm{if}& R|\stackrel{~}{}|R^{}<0\hfill \\ \hfill \gamma R|\stackrel{~}{}|R^{}& \mathrm{if}& R|\stackrel{~}{}|R^{}>0(\gamma 0)\hfill \end{array}\hfill \\ \hfill R|\stackrel{~}{}_{\mathrm{eff}}|R& =& R|\stackrel{~}{}|R+(1+\gamma )V_{\mathrm{sf}}(R)\hfill \end{array}$$
(37)
Here the “sign flip” potential is the same as that of ten Haaf et al.  and given by

$$V_{\mathrm{sf}}(R)=\underset{R_{\mathrm{na}}^{}}{}R_{\mathrm{na}}^{}|\stackrel{~}{}|R$$
(38)
where the subscript “na” (not–allowed) on $`R^{}`$ restricts the summation to the moves for which the matrix element of the hamiltonian is positive (37).

If the guiding wavefunction were to coincide with the true wavefunction, the simulation of the effective hamiltonian, which is sign free by construction, yields exact averages. So one may expect that good guiding wavefunctions lead to good upperbounds for the energy. This upperbound increases with $`\gamma `$, indicating that $`\gamma =0`$ seems the best choice, which is the effective hamiltonian of ten Haaf et al.. That hamiltonian however is a truncated version of the true hamiltonian in which all the dangerous moves are eliminated. The sign flip potential must correct this truncation by suppressing the probability that the walker will stay in a configuration with a large potential.

The second ingredient uses the fact that the hamiltonian (37) explores a larger phase space and therefore contains more information than the truncated one. Parallel to the simulation of the effective hamiltonian one can calculate the weights for the true hamiltonian. As we saw in the summary of the GFMC method, forcefully made positive transition rates still contain the correct weights when supplemented with sign functions. For the true weights of the transition probabilities as given in (37), the “sign function” must be chosen as

$$s(R,R^{})=\{\begin{array}{ccc}1\hfill & \mathrm{if}& R|\stackrel{~}{}|R^{}>0\hfill \\ 1/\gamma \hfill & \mathrm{if}& R|\stackrel{~}{}|R^{}<0\hfill \\ \frac{1ϵR|\stackrel{~}{}|R}{1ϵR|\stackrel{~}{}_{\mathrm{eff}}|R}\hfill & \mathrm{if}& R=R^{}\hfill \end{array}$$
(39)
With these “signs” in the weights a proper average can be calculated, but these averages suffer from the sign problem, the more so the smaller $`\gamma `$ is as one sees from (39). So some intermediate value of $`\gamma `$ has to be chosen. Fortunately the results are not too sensitively dependent on $`\gamma `$; the value $`\gamma =0.5`$ is a good compromise and has been taken in our simulations.

In any simulation some walkers obtain a large weight and others a small one. To lower the variance branching is regularly applied, which means a multiplication of the heavily weighted walkers in favor of the removal of those with small weight. It is not difficult to do this in an unbiased way. Sorella  proposed to use the branching much more effectively in conjunction with the signs defined in (39). The average sign is an indicator of the usefulness of the set of walkers. Start with a set of walkers with positive sign. When the average sign becomes unacceptably low, the process is stopped and a reconfiguration takes place. The walkers are replaced by another set with positive weights only, such that a number of measurable quantities gives the same average. The more observables are included the more faithful is the replacement. The construction of the equivalent set requires the solution of a set of linear equations. With the new set of walkers one continues the simulation on the basis of the effective hamiltonian and one keeps track of the true weights with signs. The reconfiguration on the basis of some observables gives at the same time a measure for these obervables. Thus measurement and reconfiguration go together. As the number of observables that can be included is limited some biases are necessarily introduced. Sorella showed that the error in the energy of the guiding wave function is easily reduced by a factor of 10, whereas reduction by the FNMC of ten Haaf et al.  rather gives only a factor 2.

## VII Results for the DMRG

In this section we give a brief summary of the results of a pure DMRG–calculation. Extensive details can be found in . The systems are strips of widths up to $`W=8`$ and of various lengths $`L`$. They are periodic in the small direction and open in the long direction. The periodicity enables us to study the spin stiffness. We have chosen open boundaries in the long direction to avoid the errors in the DMRG wavefunction due to periodic boundaries. Since we have good control of the scaling behavior in $`L`$ we extrapolate to $`L\mathrm{}`$ . In the small direction we are restricted to $`W=2,4,6`$ and 8 as odd values are not compatible with the antiferromagnetic character of the system. For wider system sizes the number of states which has to be taken into account exceeds the possiblities of the present workstations. Our criterion is that the value of the energy does not drift anymore appreciably upon the inclusion of more states. This does not mean that the wavefunction is virtually exact, since the energy is a rather insensitive probe for the wavefunction. For instance correlation functions still improve from the inclusion of more states. In Fig. 5 we present the energy as function of the ratio $`J_2/J_1`$, for strip widths 4,6 and 8 together with the best extrapolation to infinite width systems. The figure strongly suggests that the infinite system undergoes a first order phase transition around a value 0.6. This can be attributed to the transition to a columnar order (lines of opposite magnetisation). It is impossible to deduce more information from such an energy curve as other phase transitions are likely to be continuous with small differences in energy between the phases.

The spin stiffness can be calculated with the DMRG–wavefunction for systems which are periodic in at least one direction .

The result of the computation is plotted in Fig. 5. One observes a substantial decrease of $`\rho _s`$ in the frustrated region indicating the appearance of a magnetically disordered phase. In contrast to the energy the data do not allow a meaningful extrapolation to large widths. The lack of clear finite size scaling behavior in the regime of small values of $`W`$ prevents to draw firm conclusions on the disappearence of the stiffness in the middle regime.

For the correlation functions following from the DMRG wavefunction we refer to .

## VIII Results for GFMC with SR

We now come to the crux of this study: the simulations of the system with GFMC, using the DMRG wavefunctions to guide the importance sampling. All the simulations have been carried out for $`10\times 10`$ lattice with open boundaries. Standardly we have 6000 walkers and we run the simulations for about $`10^4`$ measurements. These measuring points are not fully independent and the variance is determined by chopping up the simulations into 50-100 groups, often carried out in parallel on different computers. We first give an overall assessment of the correlation function pattern and then analyze some values of the ratio $`J_2/J_1`$.

In the first series we have used the guiding wavefunction on the basis of the meandering path Fig. 3(b), because it gives a better energy than the straight option (a) . The number of basis states is $`m=75`$, which is small enough to carry out the simulations with reasonable speed and large enough that trends begin to manifest themselves. Measurements of a number of correlation functions are made in conjunction with Stochastic Reconfiguration as described in section 7. The details of these calculations are given in Table I. Note that the DMRG guiding wavefunction gives a better energy for the meandering path than for the straight path for values of $`J_2/J_1`$ up to 0.6. From 0.7 on this difference is virtually absent. This undoubtly has to do with the change to the columnar state which can equally well be realized by both paths. The value of $`ϵ`$ has been chosen as a compromise: independent measurements require a large $`ϵ`$ but the minus sign problem requires to apply often Stochastic Reconfiguration i.e. a small $`ϵ`$. One sees that in the heavily frustrated region the $`ϵ`$ must be taken small. In fact the more detailed calculations for $`J_2=0.3J_1`$ and $`J_2=0.5J_1`$ were carried out with $`ϵ=0.01`$.

In Fig. 6 and 7 we have plotted a sequence of visualizations of the correlations. From top to bottom (zig–zag) they give the correlations for the values of $`J_2/J_1`$. In order to highlight the differences a distinction is made between correlations which are above average (solid lines) and below average (broken lines). All nearest neighbor spin correlations shown are negative. In all the pictures one sees the influence of the boundaries on the spin correlations. Only 1/4 of the lattice has been pictured, the other segments follow by symmetry. The upper right corner, which corresponds to the center of the lattice, is the most significant for the behavior of the bulk. The overall trend is that spatial variations in the correlation functions occur in growing size with $`J_2/J_1`$. On the side of low $`J_1/J_2`$ (Néel phase) one sees dimer patterns in the horizontal direction, they turn over to vertical dimers (around $`J_2=0.7J_1`$) and rapidly disappear in the columnar phase. This is again support for the fact that the columnar phase is separated from the intermediate state by a first order phase transition.

Open boundary conditions have the disadvantage of boundary effects, which make it more difficult to distinguish between spontaneous and induced breaking of the translational symmetry. On the other hand for open boundaries, dimers, plaquettes or any other interruption of the translational symmetry have a natural reference frame. The correlations are not only influenced by the boundaries of the system, also the guiding DMRG–wavefunction leaves its imprint on the results. This is mainly due to the fact that we have only mixed estimators for the correlation functions, which show a mix of the guiding wavefunction and the true wavefunction. The improved estimator, used in these pictures, corrects for this effect to linear order in the deviation.<sup>*</sup><sup>*</sup>*Forward walking allows to make a pure estimate of the correlations, but requires much more calculations.  The ladder like structure in the DMRG path is reflected in a ladder like pattern in the correlations as an inspection of the correlations in the DMRG wavefunctions (not shown here) reveals. But ladders are clearly also present in the GFMC results shown in the pictures.

In order to eliminate the influence of the guiding wavefunction we scrutinize some of values of $`J_2/J_1`$ in more detail, by inspecting how the results depend on the size of the basis in the DMRG wavefunction and on the choice of the DMRG path. Since we are mostly interested in the behavior of the infinite lattice, we discuss mainly the behavior of the correlations in and around the central plaquette. So we study a sequence of DMRG wavefunctions for $`m=32`$, 75, 100, 128 and 150(200) and carry out for each of them extensive GFMC simulations. First we look to the case $`J_2=0`$, which is easy because we know that it must be Néel ordered and therefore it serves as a check on the calculations. Then we take $`J_2=0.3J_1`$ which is the most difficult case since it is likely to be close to a phase transition. Finally we inspect $`J_2=0.5J_1`$ where we are fairly sure that some dimerlike phase is realized.

### A $`J_2=0`$

For the unfrustrated Heisenberg model we have several checkpoints for our calculations. We can find to a high degree of accuracy the groundstate energy and we are sure that the Néel phase is homogeneous, i.e. that the correlations show no spatial variation other than that of the antiferromagnet. We have two ways of estimating the energy of a $`10\times 10`$ lattice. The first method is based on finite size interpolation. From DMRG calculations  we have an exact value for a $`4\times 4`$ lattice, an accurate value for the $`6\times 6`$ lattice and a good value for the $`8\times 8`$ lattice. There is also the very accurate calculation of Sandvik  for an infinitely large lattice, yielding the value of $`e_0=0.669437(5)`$. The leading finite size correction goes as $`1/L`$. Including also a $`1/L^2`$ term we have esimated the value for a $`10\times 10`$ lattice as 0.629(1) and incorporated this value in Table II(a). We stress that this is an interpolation for which the value of Sandvik is the most important input.

The second method is less well founded and uses the experience that DMRG energy estimates can be improved considerably by extrapolating to zero truncation error. When plotted as function of this truncation error the energy is often remarkably linear. In Table II(b) we give for a series of bases $`m=32,75,100,128`$ and 150, the values of the truncation error and the corresponding DMRG energy per site together with the extrapolation on the basis of linear behavior.The value for $`m=150`$ is not in line with the others. This can be explained by the fact that the construction of this DMRG wavefunction was slightly different from the others in which the basis was built up gradually. Note that the two estimates are compatible. In Table II(b) we have also listed the values of the GFMC simulations for the corresponding values of $`m`$. They do agree quite well with these estimates in particular with the one based on finite size scaling. We point out that one would have to go very far in the number of states in the DMRG calculation to obtain an accuracy that is easily obtained with GFMC. Thus the combination of GFMC and DMRG does really better than the individual components. One might wonder why there is still a drift to lower energy values in the GFMC simulations (which is also present in the tables to come). The reason is that the DMRG wavefunction is strictly zero outside a certain domain of configurations, because the truncation of the basis involves also the elimination of certain combinations of conserved quantities of the constituing parts. The domain of the wavefunction grows with the size of the basis.

Turning now to the correlations it seems that they are homogeneous in the center of the lattice for $`J_2=0`$. However a closer inspection reveals small differences. In Table III we list the asymmetries in the horizontal and vertical directions of the spin correlations in and around the central plaquette as function of the number of states. If we number the spins on the lattice as $`𝐒_{n,m}`$ with $`1n,m10`$, the central plaquette has the coordinates (5,5), (5,6), (6,5) and (6,6). We then define the asymmetry parameters $`\mathrm{\Delta }_x`$ and $`\mathrm{\Delta }_y`$ as

$$\{\begin{array}{c}\mathrm{\Delta }_x=\frac{1}{4}𝐒_{4,5}𝐒_{5,5}+𝐒_{4,6}𝐒_{5,6}+𝐒_{6,5}𝐒_{7,5}+𝐒_{6,6}𝐒_{7,6}\frac{1}{2}𝐒_{5,5}𝐒_{6,5}+𝐒_{5,6}𝐒_{6,6}\hfill \\ \mathrm{\Delta }_y=\frac{1}{4}𝐒_{5,4}𝐒_{5,5}+𝐒_{6,4}𝐒_{6,5}+𝐒_{5,6}𝐒_{5,7}+𝐒_{6,6}𝐒_{6,7}\frac{1}{2}𝐒_{5,5}𝐒_{5,6}+𝐒_{6,5}𝐒_{6,6}\hfill \end{array}$$
(40)
So $`\mathrm{\Delta }_x`$ is the average value of the correlations on the 4 horizontal bonds which are connected to the central plaquette minus the average of the values on the 2 horizontal bonds in the plaquette. Similarly $`\mathrm{\Delta }_y`$ corresponds to the vertical direction. The values for the asymmetry in Table III in the vertical direction are so small that they have no significance. Note that the anticipated decrease in $`\mathrm{\Delta }_x`$ is slow in DMRG and therefore also slow in the mixed estimator of the GFMC. The improved estimator (35) however is truely an improvement! So one sees that all the observed small deviations from the homogeneous state will disappear with the increase of the number of states in the basis of the DMRG wavefunction. (In general the accuracy of the correlations is determined by that of the GFMC simulations. We get as variance a number of the order 0.01, implying twice that value for the improved estimator) The vanishing of $`\mathrm{\Delta }_x`$ and $`\mathrm{\Delta }_y`$ also prove that finite size effects are small in the center of the $`10\times 10`$ lattice. From these data we may conclude that the GFMC can make up for the errors in the DMRG wavefunction for a relative low number of basis states. We have not carried out a similar series for the straight path since this will certainly show no dimers as will become clear from the following cases.

### B $`J_2=0.3J_1`$

This case is the most difficult to analyze since it is expected to be close to a continuous phase transition from the Néel state to a dimerlike state. As is known  the DMRG structure of the wavefunction is not very adequate to cope with the long–range correlation in the spins typical for a critical point. In Table IV we have presented the same data as in Table III but now for $`J_2=0.3`$. There is no pattern in the energy as function of the truncation error $`\delta `$. The decrease of the energy as function of the size of the basis $`m`$ is in the DMRG wavefunctions is not saturated. The GFMC simulations lead to a notably lower energy and they do hardly show a leveling off as function of the basis of the guiding wavefunction. All these points are indicators that the DMRG wavefunction is rather far from convergence and that more accurate data would require a much larger basis. As far as the staggering in the correlations is concerned the values for $`\mathrm{\Delta }_x`$ are significant, also because the simulation results generally increase the values. Those for $`\mathrm{\Delta }_y`$ are not small enough to be considered as noise. Given the fact that most authors locate the phase transition at higher values $`J_20.4J_1`$ we would expect both $`\mathrm{\Delta }`$’s to vanish. So either the dimerlike state is realized for values as low as $`J_2=0.3J_1`$ or dimer formation already starts in the Néel state.

To get more insight in the nature of the groundstate we have also carried out the same set of simulations on the straight path (a) in Fig. 3. This guiding wavefunction shows virtually no formation of dimers in any direction as can be observed from Table V. In spite of the fact that the trends indicated in the table have not come to convergence one may draw a few conclusions from the comparison of the two sets of simulations. The overal impression is that the meandering guiding wavefunction represents a groundstate of a different symmetry as compared to the straight path guiding wavefunction. The meandering wavefunction prefers dimers in the horizontal direction and the straight wavefunction leads to some dimerization in the vertical direction. The difference also shows up in the energy, it is not only large on the DMRG level but it also persists at the GFMC level. We see similar trends in the next case.

### C $`J_2=0.5J_1`$

By any estimate this value of the next nearest neighbor coupling leads to a dimerlike state if it exists at all. No accurate data are available on the energy of the $`10\times 10`$ system to compare to our results. In Table VI we list the data for a set of DMRG wavefunctions with bases $`m=32`$, 75, 100, 128, 150 and 200. The DMRG values of the energy (with exception of the value for $`m=32`$) can be extrapolated to zero truncation error with the limiting value $`E_0=48.4(1)`$, which corresponds very well with the level in the GFMC values for larger sizes of the basis. This indicates again that the GFMC simulations can make up for the shortcoming of the DMRG wavefunction. One would indeed have to enlarge the basis to $`m`$ of the order of 1000 in order to achieve the value of the energy of the simulations which use DMRG guiding wavefunctions with a basis of the order of 100.

The staggering in the correlations expressed by the quantities $`\mathrm{\Delta }_x`$ for the horizontal direction and $`\mathrm{\Delta }_y`$ for the vertical direction, has values that are significant. If one looks to the contributions of the DMRG wavefunction and the GFMC simulation separately, one observes that the overall values do agree quite well, with the tendency that the GFMC simulations lowers the staggerring in the horizontal direction and slightly increases it in the vertical direction. So we may conclude that indeed in the groundstate of the $`J_2=0.5J_1`$ system, the correlations of the spins are not translation invariant but show a staggering. However these results neither confirm the picture that the dimerstate is the lowest (as suggested by Kotov et al. ) nor that the plaquettestate is the groundstate (as concluded by Capriotti and Sorella ). We comment on these discrepancies further in the discussion.

Again it is worthwhile to compare these results with a simulation on the basis of the straight path (a) in Fig. 3. Here it is manifest that the straight path prefers to have the dimers in the vertical direction. Again the impression is that the straight path leads to a different symmetry as compared to the meandering path. It is not only the different preference in the main direction of the dimers, also the secondary dimerization in the perpendicular direction, notably in the meandering case, is not present in the straight case. The fairly large difference in energy on the DMRG level becomes quite small on the GFMC level.

## IX Discussion

We have presented a method to employ the DMRG wavefunctions as guiding wavefunctions for a GFMC simulation of the groundstate. Generally the combination is much better than the two individual methods. The GFMC simulations considerably improve the DMRG wavefunction. In the intermediate regime the properties of the GFMC simulations depend on the guiding wavefunction as the results for two different DMRG guiding wavefunctions show.

The method has been used to observe spin correlations in the frustrated Heisenberg model on a square lattice. In this discussion we focus on the intermediate region where the model is most frustrated and which is the “piece de resistance” of the present research. We see patterns of strongly correlated nearest neighbor spins, to be called dimers. To indicate what me mean by strong and weak we give the values in and around the central square of the $`10\times 10`$ lattice, for the case $`J_2=0.5J_1`$. In Fig. 7(c) we have given the values of the central square extrapolated to an infinite lattice.

The values are based on the improved estimator and it is interesting to see the trends. The horizontal strong correlation of -0.42 is the result of the DMRG value -0.44 and the GFMC value -0.43, while the weak bond -0.15 is the result of the DMRG value -0.09 and the GFMC value -0.12. Thus the GFMC weakens the order parameter $`D_x`$ associated with the staggering. For the vertical direction there is hardly a change from DMRG to GFMC. One has to go to the next decimal to see the difference. The strong bond equals -0.368 and is coming from the DMRG value -0.375 and the GFMC value -0.371, while the improved weak bond of -0.271 is the resulting value of -0.275 for DMRG and -0.273 for GFMC.

Before we comment on this result we discuss the influence of the choice of the guiding wave function. We note that for both points $`J_2=0.3J_1`$ and $`J_2=0.5J_1`$ the two choices for the DMRG wavefunction give different results. First of all the main staggering is for the meandering path (b) of Fig. 3 in the horizontal direction, while the straight path (a) of Fig. 3 prefers the dimers in the vertical direction. There is not much difference in the values of the strong and weak correlations. Secondly the straight path shows no appreciable staggering in the other direction, so one may wonder whether the observed effect for the meandering path is real. In our opinion this difference has to do with the effect that the DMRG wavefunction “locks in” on a certain symmetry. The straight path yields a groundstate which is truely dimerlike in the sense that it is translational invariant in the direction perpendicular to the dimers. The meandering path locks in on a different groundstate which holds the middle between a dimerlike and a plaquettelike state. The GFMC simulations cannot overcome this difference in symmetry, likely because the two lowest states with different symmetry are virtually orthogonal. On the DMRG level there is a large difference in energy between the two states, favoring the meandering path strongly, on the GFMC level this difference has become very small. With this observation in mind we compare our result with other findings.

The results of the series expansions ,  and  are shown in Fig. 7(a). Their correlations organize themselves in spinladders. The correlations on the rungs of the ladder are $`0.45\pm 0.5`$ which compares well with our strongest horizontal correlation and this holds also for the weak horizontal correlation (–0.12 vs –0.15). The most noticeble difference is the value of our weak correlation in the vertical direction (–0.27 vs –0.36) while the strong correlation (–0.37 vs –0.36) agrees. There is no real conflict between our result and theirs since the symmetry they find is fixed by the state around which the series expansion is made. So our claim is only that our state with different symmetry is the lower one. In fact in the paper of Singh et al. , it is noted that the susceptibility to a staggering operator in the perpendicular direction (our $`\mathrm{\Delta }_y`$) becomes very large in the dimer state for $`J_2=0.5J_1`$ which we take as an indication of the nearby lower state. The analytical calculations in  and  however do not support the existence of the state we find.

Neither do we find support for the plaquette state found in , which we have sketched in Fig. 7(b). The evidence of this investigation is based on the boundedness of the susceptibility for the operator which breaks the orientational symmetry and the divergence of the susceptibility for the order parameter breaking translational invariance (corresponding to $`\mathrm{\Delta }_x`$). They have not separately investigated the values of $`\mathrm{\Delta }_x`$ and $`\mathrm{\Delta }_y`$ since their groundstate has the symmetry of the lattice and one would find automatically the same answer. They conclude that in absence of orientational order parameter and with the presence of the translational order parameter the state must be plaquettelike. We believe that their result is influenced by the guiding wavefunction for which the one-step Lanczos approximation is taken. This wavefunction certainly has the symmetry of the square and again GFMC cannot find a groundstate with a different symmetry.

Finally we comment on the fact that we find the dimerization already for values as low as $`J_2=0.3J_1`$ at least for the meandering path. As we have mentioned earlier the results as function of the number of states have not sufficiently converged to make a firm conclusion, the more so since there is a large difference between DMRG and GFMC. Still it could be an indication that the phase transition from the Néel state to the dimer state takes place for lower values than the estimated $`J_2=0.38J_1`$ .

Thus many questions are left over, amongst others how the order parameters behave as function of the frustation ratio in the intermediate region. We feel that the combination of DMRG and GFMC is a good tool to investigate these issues since they demonstrate ad oculos the correlations in the intermediate state.

Acknowledgement The authors are indebted to Steve White for making his software available. One of us (M. S. L. duC. de J.) gratefully acknowledges the hospitality of Steve for a stay at Irvine of 3 months, where the basis of this work was laid. The authors have also benefitted from illuminating discussions with Subir Sachdev and Jan Zaanen. The authors want to acknowledge the efficient help of Michael Patra with the simulations on the cluster of PC’s of the Instituut-Lorentz.


format: text

text:
# Self-Duality in Superconductor-Insulator Quantum Phase Transitions

## ACKNOWLEDGMENTS

I’m grateful to M. Krusius for the kind hospitality at the Low Temperature Laboratory in Helsinki. I wish to acknowledge useful conversations with J. Hove, N. Kopnin, M. Krusius, K. Nguyen, M. Paalanen, A. Sudbø, and especially G. Volovik.

This work was funded in part by the EU sponsored programme Transfer and Mobility of Researchers under contract No. ERBFMGECT980122.


text:
# Memory Effects in Granular Material

## Abstract

We present a combined experimental and theoretical study of memory effects in vibration-induced compaction of granular materials. In particular, the response of the system to an abrupt change in shaking intensity is measured. At short times after the perturbation a granular analog of aging in glasses is observed. Using a simple two-state model, we are able to explain this short-time response. We also discuss the possibility for the system to obey an approximate pseudo-fluctuation-dissipation theorem relationship and relate our work to earlier experimental and theoretical studies of the problem.

PACS numbers: 45.70.-n, 61.43.Fs, 81.05.Rm

Granular materials comprise an important class of complex systems whose simple fundamental mechanics gives rise to rich macroscopic phenomenology . Recent experiments on granular compaction  suggest they are an ideal system for studying jamming, a phenomenon lying outside the domain of conventional statistical physics, yet highly reminiscent of glassiness. These studies showed that a loose packing of glass beads subjected to vertical “tapping” slowly compacts, asymptoting to a higher steady state packing fraction. This “equilibrium” packing fraction is somewhat lower than the random close packing limit, $`\rho _{\mathrm{rcp}}0.64`$, and is a decreasing function of the vibration intensity, typically parameterized by $`\mathrm{\Gamma }`$, the peak applied acceleration normalized by gravity, $`\mathrm{g}`$. The relaxation dynamics are extremely slow, taking many thousands of taps for the packing fraction, $`\rho `$, to approach its steady state value. During this evolution, $`\rho `$ increases logarithmically with the number of taps, $`t`$, which is typical for self–inhibiting processes . The average time scale $`\tau `$ of the relaxation decreases with $`\mathrm{\Gamma }`$, and in this sense the shaking intensity plays, at least qualitatively, the role of temperature. For small $`\mathrm{\Gamma }`$, the relaxation rate becomes so slow that the system cannot reach the steady state density within the experimental time scale. It was also found that compaction can be maximized through an annealing procedure. This process involves a slow “cooling” of the system starting from a high shaking intensity $`\mathrm{\Gamma }`$. These slow relaxation and annealing properties of this system are reminiscent of conventional glasses. Another qualitative similarity to glasses is observable in the density fluctuation spectrum of the granular system near equilibrium. The spectrum was found to be strongly non-Lorentzian , revealing the existence of multiple time scales in the system. The shortest and the longest relaxation timescales differ by as much as three order of magnitude, and the behavior of the spectrum for the intermediate frequencies is highly non-trivial; in certain regimes it can be fitted with a power law.

These previous experimental observations are suggestive of glassy behavior and this connection has been explored in recent models of compaction using ideas from magnetic systems . However, a more direct test of the glassy nature of granular compaction comes from measurements of the response of the system to sudden perturbations in the effective temperature, given by $`\mathrm{\Gamma }`$. This idea originates from classical experiments for the study of aging in glasses  and has recently been explored using computer simulations . In this letter we present direct experimental observations of memory effects in a vibrated granular system obtained by measuring the short-time response to an instantaneous change in tapping acceleration $`\mathrm{\Gamma }`$ and propose a simple theoretical framework.

We used the experimental set-up described in refs. : 1 mm–diameter glass beads were vertically shaken in a tall, evacuated, 19 mm-diameter glass tube, and the packing density of the beads was measured using capacitors mounted at four heights along the column.

The simplest form of this experiment consists of a single instantaneous change of vibration intensity from $`\mathrm{\Gamma }_1`$ to $`\mathrm{\Gamma }_2`$ after $`t_0`$ taps. For $`\mathrm{\Gamma }_2<\mathrm{\Gamma }_1`$ (Fig. 1a) we found that on short time scales the compaction rate increases. This is in sharp contrast to what one may expect from the long-time behavior found in previous experiments where the relaxation is slower for smaller vibration accelerations. For $`\mathrm{\Gamma }_2>\mathrm{\Gamma }_1`$ (Fig. 1b) we found that the system dilates immediately following $`t_0`$.

These results too, are opposite from the long-time behavior seen in previous experiments where the compaction rate increased: Not only does the compaction rate decrease, it becomes negative (i.e. the system dilates). Note that after several taps the “anomalous” dilation ceases and there is a crossover to the “normal” behavior, with the relaxation rate becoming the same as in constant–$`\mathrm{\Gamma }`$ mode. Thus, most of the shaking history is forgotten after a short time.

These data constitute a short-term memory effect: the future evolution of $`\rho `$ after time $`t_0`$ depends not only on $`\rho (t_0)`$, but also information about the previous tapping history, contained in other “hidden” variables. In order to demonstrate this in a more explicit manner, we modified the above experiment. In this second set of three experiments the systems was driven to the same density $`\rho _0`$ with three different accelerations $`\mathrm{\Gamma }_0`$, $`\mathrm{\Gamma }_1`$, and $`\mathrm{\Gamma }_2`$. After $`\rho _0`$ was achieved at time $`t_0`$, the system was tapped with the same intensity $`\mathrm{\Gamma }_0`$ for all three experiments. As seen in Figure 2, the evolution for $`t>t_0`$ strongly depends on the pre-history. The need for extra state variables in the problem is consistent with strongly non-Lorentzian behavior of the fluctuation spectrum, observed in earlier experiments. Indeed, if the evolution could be prescribed by a single master equation for local density, it would result in a single-relaxation-time exponential decay of the density fluctuations near equilibrium. Instead, a wide distribution of characteristic times is suggested by the spectrum.

To give a theoretical interpretation of the above results, we view the problem as an evolution in the space of discrete “microscopic” states corresponding to different realizations of the packing topology (i.e. of the contact network). For each tap there is a possibility for a transition from one microscopic state to another. Since the dynamics is dissipative and the system is under external gravity, a transition to a denser configuration is typically more probable that the reverse one. We now assume that the short–term dynamics of the system are dominated by a number of local flip–flop modes with relatively high transition rates in both directions.

This model replaces the complicated configuration space with a set of independent two-state systems, each of which is characterized by two transition rates, $`\kappa _{eg}>\kappa _{ge}`$. $`\kappa _{eg}/\kappa _{ge}`$ gives the ratio of the equilibrium probabilities of populating each state: “ground” and “excited”. As we have argued, the higher probability ground state is typically the one with higher density, i.e. the volume change $`v`$ between the ground and the excited states is normally positive (see Fig. 3 for a schematic description of the model). Our two-state approximation is close in its spirit to recent Grinev–Edwards and de Gennes models .

We now introduce the concept of a base-line density, $`\rho _b`$, which corresponds to all the elementary modes at their ground states. Obviously, the experimentally–observed density is different from $`\rho _b`$ due to a non-zero fraction of excited states:

$$\rho =\rho _b(t)\left(1\frac{1}{V}\underset{n}{}v^{(n)}\left(1+\frac{\kappa _{ge}^{(n)}}{\kappa _{eg}^{(n)}}\right)^1\right).$$
(1)
The summation here is performed over all the dominant two-state modes, $`V`$ is the total volume and $`v^{(n)}`$ is the volume difference between the excited and the ground $`n`$-th state. Assuming that the vibration intensity $`\mathrm{\Gamma }`$ is a qualitative analog of temperature, we expect the population of the excited states, $`P(\mathrm{\Gamma })=(1+\kappa _{ge}^{(n)}/\kappa _{eg}^{(n)})^1`$, to grow with $`\mathrm{\Gamma }`$, starting from zero at $`\mathrm{\Gamma }=0`$. Hence, for a given $`\rho _b`$, the total density $`\rho `$ will be lower at higher acceleration. This explains the observed effect of an abrupt change of $`\mathrm{\Gamma }`$. After a switch from $`\mathrm{\Gamma }_1`$ to $`\mathrm{\Gamma }_2`$ at time $`t_0=0`$, the flip-flop mode contribution to the total density would relax to its new equilibrium value in the following way:

$$\mathrm{\Delta }_{\mathrm{\Gamma }_1\mathrm{\Gamma }_2}(t)=\rho _bvF_{\mathrm{\Gamma }_1,\mathrm{\Gamma }_2}(v,\kappa )\left(1exp(\kappa t)\right)𝑑v𝑑\kappa $$
(2)
Here $`\kappa `$ is the relaxation rate of an individual mode, and the distribution function $`F_{\mathrm{\Gamma }_1,\mathrm{\Gamma }_2}(v,\kappa )`$ is introduced as follows:

$`F_{\mathrm{\Gamma }_1,\mathrm{\Gamma }_2}(v,\kappa ){\displaystyle \frac{1}{V}}{\displaystyle \underset{n}{}}\left(P^{(n)}(\mathrm{\Gamma }_2)P^{(n)}(\mathrm{\Gamma }_1)\right)\delta (vv^{(n)})`$ (3)
$`\delta \left(\kappa \kappa _{ge}^{(n)}(\mathrm{\Gamma }_2)\kappa _{eg}^{(n)}(\mathrm{\Gamma }_2)\right).`$ (4)
Since the observed density changes in compaction experiment are normally less than $`1\%`$ of the total density, one can estimate the typical separation between neighboring flip-flop systems as 5 particle sizes, which yields a good justification for our no-coupling approximation. According to Eq. (2), if $`F_{\mathrm{\Gamma }_1,\mathrm{\Gamma }_2}`$ does not vanish in the limit $`\kappa 0`$, the late stage of the relaxation of $`\mathrm{\Delta }_{\mathrm{\Gamma }_1\mathrm{\Gamma }_2}(t)`$ is given by the power law:

$`\delta _{\mathrm{\Gamma }_1,\mathrm{\Gamma }_2}(t)\mathrm{\Delta }_{\mathrm{\Gamma }_1\mathrm{\Gamma }_2}(t)\mathrm{\Delta }_{\mathrm{\Gamma }_1\mathrm{\Gamma }_2}(\mathrm{})=`$ (5)
$`\rho _b{\displaystyle vF_{\mathrm{\Gamma }_1,\mathrm{\Gamma }_2}(v,\kappa )e^{\kappa t}𝑑v𝑑\kappa }{\displaystyle \frac{1}{t}}.`$ (6)
Note that $`\rho _b`$ is also dependent on time: although this cannot be described within our two-state approximation, the collection of elementary modes slowly evolves. Thus, one can observe two different processes: on short time scales, a fast relaxation due to the flip–flop modes is dominant, while over the long times, the dynamics are determined by the logarithmically slow evolution of the baseline density $`\rho _b(t)`$. The crossover between the two regimes is particularly obvious in Fig. 1b, where it results in a non-monotonic evolution. For experiments performed at sufficiently late stages of the density relaxation, the dynamics of the baseline density could be neglected compared to the contribution of the flip-flop modes (note that what we call a late-stage relaxation corresponds in fact to mesoscopic time scales which are always shorter than the relaxation time for $`\rho _b`$). It has to be emphasized that the described experiments provide us with a tool for study of the response of the system, which is not limited to the nearly–equilibrium regime.

One can use our simple model to predict the response of the system to a more complicated pattern of changes of $`\mathrm{\Gamma }`$. First, we reach, using annealing dynamics, a “quasi-steady” state at amplitude $`\mathrm{\Gamma }_0`$, so that one can consider $`\rho _b`$ constant later on. Let us switch the shaking acceleration from $`\mathrm{\Gamma }_0`$ to $`\mathrm{\Gamma }_1`$ for a finite number of taps $`\delta t`$, and then switch it back to $`\mathrm{\Gamma }_0`$. During the intermediate $`\mathrm{\Gamma }_1`$–stage, the system does not have enough time to completely relax to its new equilibrium. In our two-state model, the modes whose relaxation rate (at $`\mathrm{\Gamma }_1`$) is below $`\delta t^1`$ remain unrelaxed. Assuming that the slow modes at $`\mathrm{\Gamma }_1`$ are mostly the same as at $`\mathrm{\Gamma }_0`$, we can calculate the backward density relaxation similarly to Eq. (5), with $`F(v,\kappa )`$ effectively depleted below a minimal rate, $`\kappa _0`$. This cut-off frequency, $`\kappa _0`$, is expected to decrease monotonically with increasing perturbation duration $`\delta t`$. The resulting density relaxation after returning to $`\mathrm{\Gamma }_0`$ is given by:

$$\delta _{\mathrm{\Gamma }_1,\mathrm{\Gamma }_0}(t)\frac{exp(\kappa _0t)}{t}.$$
(7)
We tested these predictions by performing this three stage experiment, varying the duration, $`\delta t`$, of the perturbation($`\mathrm{\Gamma }_1`$) stage (Fig. 4). As predicted, the time needed to recover the steady-state density increases with the number of taps $`\delta t`$ spent in the “hot” regime $`\mathrm{\Gamma }_1>\mathrm{\Gamma }_0`$. In the coordinates chosen, the relaxation curves should follow the $`\delta t=\mathrm{}`$ dynamics until the saturation at the cut-off time, $`\kappa _0^1(\delta t)`$. We approximate the distribution function $`F`$ by a constant above this low frequency cut-off at $`\kappa _0^1(\delta t)`$, up to a high-frequency cut-off, $`\kappa _h1\mathrm{t}\mathrm{a}\mathrm{p}^1`$. This eliminates the unphysical low-$`t`$ divergence in Eq. (7). Figure 4 shows fits of the data to Eq. (2), where $`\kappa _0(\delta t)`$ is determined from the fit. The best-fit is achieved at $`\kappa _h=0.4`$, and the variation of this parameter would result in a simple rescaling of the time axis.

Figure 4 demonstrates good agreement between model and experiment, with some systematic error at the earliest relaxation stage (which is an expected result of our oversimplified description of the short time dynamics). For the late stage relaxation, we conclude that (i) within our experimental precision, the $`\delta t=\mathrm{}`$ relaxation is consistent with the predicted $`1/t`$ law; (ii) finite–$`\delta t`$ relaxation curves can be parameterized by a low frequency cut-off, $`\kappa _0`$; and (iii) $`\kappa _0`$ is a decreasing function of the waiting time $`\delta t`$, shown in the insert of the Figure 4. We now relate our picture to previous experimental and theoretical results. As discussed earlier, the wide range of relaxation times reveals itself both in our response measurements and in the the fluctuation spectra of the density. It is tempting to relate these two kinds of data through an analog of a Fluctuation-Dissipation Theorem (FDT). Of course, there is no fundamental reason for FDT to be applicable to the granular system . Even though the above two-state model could be mapped onto a thermal system (in which FDT is expected to work), the thermodynamic variable conjugate to density in the context of such a mapping has no clear physical meaning. Nevertheless, below we outline the pseudo-FDT relationship expected for the granular system under rather natural approximation. Namely, we neglect the correlation between the volume change $`v`$ and the life time $`\kappa ^1`$ of an individual mode, i.e. assume $`F_{0,\mathrm{\Gamma }}(v,\kappa )=f(\kappa )g(v)`$. Then the density autocorrelation function can be written as follows:

$`\delta \rho (0)\delta \rho (t)_\mathrm{\Gamma }={\displaystyle \frac{\rho ^2}{2V}}{\displaystyle \left(v^2v^2\right)exp(\kappa t)f(\kappa )𝑑\kappa }=`$ (8)
$`\rho {\displaystyle \frac{v^2v^2}{2vV}}\delta _{0,\mathrm{\Gamma }}(t).`$ (9)
Thus, the density correlator is simply proportional to the response function corresponding to the switch between a very low acceleration (at which virtually all the modes are in their ground states) and the given one, $`\mathrm{\Gamma }`$. An experimental check of this relationship requires further high-precision studies of both the relaxation dynamics and the fluctuation spectrum.

Our model also gives a simple interpretation to the decreasing dependence of the steady-state density on $`\mathrm{\Gamma }`$: it can be attributed to the growth of the population of the excited states, $`P(\mathrm{\Gamma })`$. Indeed, the corresponding correction to the total density is about $`1\%`$, i.e. of the same order as the variation of the equilibrium packing fraction with $`\mathrm{\Gamma }`$ .

The slow dynamics associated with the evolution of the base-line density can also be addressed within our approach. For doing so we need to account for the coupling of individual modes. Namely, it would be a reasonable assumption that a relaxation of one mode to its ground state may frustrate such a transition for some of its neighbors (e.g. in 3D the most compact local cluster can be created only at the expense of less dense neighboring regions). Thus, we arrive at an effective anti-ferromagnetic (AF) coupling (of an infinite strength) between the two-state modes. This extension of our model makes it remarkably similar to the so-called reversible Parking Lot Model (PLM) , which has been successful in describing many aspects of granular compaction experiments . Recent simulations based on the “tetris model”  for compaction also find slow glassy responses to changes in $`\mathrm{\Gamma }`$, but do not capture the short-term memory effect described here. .

In the PLM, D-dimensional space (the parking lot) gets packed with finite-size objects (cars), which may arrive and depart with fixed rates and which are not allowed to overlap. Now, the transcendental relationship between PLM and the granular compaction experiment is easier to explain: both PLM and our coupled flip-flop model belong to the same generic class of AF-type systems (in the case of PLM, a local two-state mode is represented by a particle whose center of mass may or may not be placed at point x; the mode coupling is due to the hard-core interactions). The PLM is known to capture the slow dynamics of granular compaction and some features of its fluctuation spectrum . In fact, we performed numerical simulations of the PLM that display  the memory effects discussed in this work.

In conclusion, we used a sequence of abrupt switches of the shaking intensity $`\mathrm{\Gamma }`$ to study the response of a vibrated granular system. This technique can be used in the vicinity of the steady state density, as well as far from equilibrium. The major result is the direct demonstration of a memory effect in the system: the evolution is not predetermined by the local density alone, and its description requires introduction of additional “hidden” variables. Our phenomenological model for this behavior is built on minimal assumptions about the dynamics of the system and produces results which are generic and are expected to be valid for a wide class of more realistic microscopic models.

Acknowledgments: We would like to thank Sue Coppersmith, Leo Kadanoff, Sid Nagel, and Tom Witten for insightful discussions and Fernando Villarruel for help with the experiments. This work was supported by the NSF under Award CTS-9710991 and by the MRSEC Program of the NSF under Award DMR-9808595.


text:
# TAUP - 2621-2000WIS2000/3/Feb.-DPP The contribution of Σ^∗→Λ𝜋 to measured Λ polarizationSupported in part by grants from US-Israel Bi-National Science Foundation and from the Israeli Science Foundation.

## Abstract

Calculations of the polarization of $`\mathrm{\Lambda }`$ and $`\overline{\mathrm{\Lambda }}`$ particles after fragmentation of a polarized quark produced in processes like $`Z`$-decay and deep inelastic polarized lepton scattering must include $`\mathrm{\Lambda }`$ and $`\overline{\mathrm{\Lambda }}`$ produced as decay products of $`\mathrm{\Sigma }^0`$ and $`\mathrm{\Sigma }^{}`$ as well as those produced directly. These decay contributions are significant and not feasibly included in theoretical calculations based on QCD without additional input from other experimental data. Furthermore these contributions depend on the spin structure of the $`\mathrm{\Sigma }^0`$ or $`\mathrm{\Sigma }^{}`$ and are not directly related to the structure function of the $`\mathrm{\Lambda }`$

Draft: , 4:11 PM

The interpretation of studies of the nucleon spin structure functions is that quarks in the nucleon carry only $``$30% of the nucleon spin and that the strange (and non-strange) sea is polarized opposite to the polarization of the valence quarks. An attempt to shed light on this very problematic conclusion was made through measurement of the polarization of $`\mathrm{\Lambda }`$ and $`\overline{\mathrm{\Lambda }}`$ produced near the Z pole in $`e^+e^{}`$ collisions  and in polarized lepton Deep Inelastic Scattering on unpolarized targets . Several theoretical works were published on this subject in which predictions and calculations relevant to the interpretation of the experimental results were presented . A difficulty present in most experiments is that they cannot distinguish between $`\mathrm{\Lambda }`$ and $`\overline{\mathrm{\Lambda }}`$ produced directly or as decay products. The main contribution from decays is the $`\mathrm{\Sigma }^{}\mathrm{\Lambda }\pi `$. The purpose of the present note is to emphasize the importance of taking into account the contribution from this process.

To illustrate this point, if 100% polarized strange quarks hadronize directly to a $`\mathrm{\Lambda }`$ , the $`\mathrm{\Lambda }`$ polarization will be also 100% if the spin structure of the $`\mathrm{\Lambda }`$ is as expected by the naïve quark model. On the other hand, if they hadronize to a $`\mathrm{\Sigma }^{}`$ the $`\mathrm{\Lambda }`$ particles resulting from its decay will be only 55% polarized . If the spin structure of the $`\mathrm{\Lambda }`$ is as derived from SU(3) symmetry using the measured spin structure of the proton, 100% strange quarks will result in 73% polarized $`\mathrm{\Lambda }`$ if produced directly compared with the same 55% if coming from $`\mathrm{\Sigma }^{}`$ decay.

The polarization of the $`\mathrm{\Lambda }`$ particles observed in any experiment can be written

$$P(\mathrm{\Lambda })=\frac{N_{nd}P_{nd}(\mathrm{\Lambda })+N_\mathrm{\Sigma }^{}BR(\mathrm{\Sigma }^{}\mathrm{\Lambda }\pi )P_\mathrm{\Sigma }^{}(\mathrm{\Lambda })+N_{\mathrm{\Sigma }^o}P_{\mathrm{\Sigma }^o}(\mathrm{\Lambda })}{N_\mathrm{\Lambda }}$$
(1)
where $`N_\mathrm{\Lambda }`$, $`N_\mathrm{\Sigma }^{}`$ and $`N_{\mathrm{\Sigma }^o}`$ denote respectively the numbers of $`\mathrm{\Lambda }`$’s, $`\mathrm{\Sigma }^{}`$’s and $`\mathrm{\Sigma }^o`$’s produced in the experiment, $`BR(\mathrm{\Sigma }^{}\mathrm{\Lambda }\pi )`$ denotes the branching ratio for the $`\mathrm{\Sigma }^{}\mathrm{\Lambda }\pi `$ decay, $`P_\mathrm{\Sigma }^{}(\mathrm{\Lambda })`$ and $`P_{\mathrm{\Sigma }^o}(\mathrm{\Lambda })`$ denote respectively the polarizations of the $`\mathrm{\Lambda }`$’s produced via the $`\mathrm{\Sigma }^{}\mathrm{\Lambda }\pi `$ decay and the $`\mathrm{\Sigma }^o\mathrm{\Lambda }\gamma `$ decay, and $`N_{nd}`$ and $`P_{nd}(\mathrm{\Lambda })`$ denote the number and polarization of $`\mathrm{\Lambda }`$’s produced via all other ways; i.e. which which do not go via the $`\mathrm{\Sigma }^{}`$ or $`\mathrm{\Sigma }^o`$,

$$N_{nd}=N_\mathrm{\Lambda }N_\mathrm{\Sigma }^{}BR(\mathrm{\Sigma }^{}\mathrm{\Lambda }\pi )N_{\mathrm{\Sigma }^o}$$
(2)
The individual terms in the numerator of eq.(1) are all distinct and measurable. Any calculation of the polarization of the final observed $`\mathrm{\Lambda }`$ must consider all these contributions if they are not separated experimentally.

One might argue that the $`\mathrm{\Sigma }^{}\mathrm{\Lambda }\pi `$ decay is a strong interaction described in terms of quarks and gluons in QCD and should be included in the inclusive polarized fragmentation function. Clearly the $`\mathrm{\Sigma }^{}`$ intermediate state must be already included in any fragmentation function which takes into account $`all`$ strong interactions in the description of a process in which a struck quark turns into a $`\mathrm{\Lambda }`$ plus anything else. This point of view is implied in the treatments  which attempt to use the $`\mathrm{\Lambda }`$ polarization data to extract fragmentation functions. But the expression for the $`\mathrm{\Lambda }`$ polarization eq.(1) is rigorous. Thus a theoretical formulation which gives a prediction for this polarization must also include a prediction for the precise values for all parameters appearing in eq.(1)

We immediately find a crucial weak point in all attempts to obtain a theoretical estimate for the value of the $`\mathrm{\Lambda }`$ polarization. Any theoretical attempt to obtain the value of the branching ratio $`BR(\mathrm{\Sigma }^{}\mathrm{\Lambda }\pi )`$ must take into account fine threshold effects like the small SU(3) breaking produced by the $`\mathrm{\Lambda }\mathrm{\Sigma }`$ mass difference which vanishes in the SU(3) symmetry limit. Note that in the SU(3) symmetry limit the predicted ratio of the branching ratios of the two $`\mathrm{\Sigma }^{}`$ decay modes is in strong disagreement with experiment:

$$\left(\frac{BR(\mathrm{\Sigma }^{}\mathrm{\Lambda }\pi )}{BR(\mathrm{\Sigma }^{}\mathrm{\Sigma }\pi )}\right)_{theo}=1/2\left(\frac{BR(\mathrm{\Sigma }^{}\mathrm{\Lambda }\pi )}{BR(\mathrm{\Sigma }^{}\mathrm{\Sigma }\pi )}\right)_{exp}=7.3\pm 1.2$$
(3)
The disagreement is more than an order of magnitude.

There is also the problem of obtaining the values of $`N_{nd}`$, $`N_\mathrm{\Sigma }^{}`$ and $`N_{\mathrm{\Sigma }^o}`$. If one assumes a purely statistical model in which all states of two nonstrange quarks and one strange quark are equally probable the ratio

$$N_{nd}/N_\mathrm{\Sigma }^{}/N_{\mathrm{\Sigma }^o}=1:6:1$$
(4)
where the factor 6 arises from the $`(2J+1)`$ spin factor and the three charge states of the $`\mathrm{\Sigma }^{}`$ which all decay into $`\mathrm{\Lambda }\pi `$. The experimental values are very different. There is also the problem that all three charge states are equally produced by the fragmentation of a struck $`s`$ quark, while a struck $`u`$ or $`d`$ quark can only produce the two charge states containing the struck quark. All these factors complicate any attempt at this stage to predict values of $`N_{nd}`$, $`N_\mathrm{\Sigma }^{}`$ and $`N_{\mathrm{\Sigma }^o}`$ from any purely theoretical model without any external experimental input.

Thus predictions for polarization of the $`\mathrm{\Lambda }`$’s observed in any experiment must include as input the known experimental branching ratio $`BR(\mathrm{\Sigma }^{}\mathrm{\Lambda }\pi )`$ as well as the values of $`N_{nd}`$, $`N_\mathrm{\Sigma }^{}`$ and $`N_{\mathrm{\Sigma }^o}`$ obtained from other experiments or from Monte Carlo programs which rely on a number of free parameters which are adjusted to fit vast quantities of data. Note that the $`\mathrm{\Sigma }^o`$ decays electromagnetically. Its decay is never included in any strong interaction fragmentation function and the $`\mathrm{\Lambda }^{}s`$ produced via its production and decay must be considered separately in all fragmentation models.

We now note that the polarization of $`\mathrm{\Lambda }`$ produced from the decay of a $`\mathrm{\Sigma }^{}`$ or $`\mathrm{\Sigma }^o`$ is proportional to the polarization of the decaying $`\mathrm{\Sigma }^{}`$ or $`\mathrm{\Sigma }^o`$ with coefficients depending only on angular momentum Clebsch-Gordan coefficients and completely independent of the spin structure of the $`\mathrm{\Lambda }`$:

$$P_\mathrm{\Sigma }^{}(\mathrm{\Lambda })=P_\mathrm{\Sigma }^{}C(\mathrm{\Sigma }^{});P_{\mathrm{\Sigma }^o}(\mathrm{\Lambda })=P_{\mathrm{\Sigma }^o}C(\mathrm{\Sigma }^o);P_{\mathrm{\Sigma }^{}\mathrm{\Sigma }^o}(\mathrm{\Lambda })=P_\mathrm{\Sigma }^{}C(\mathrm{\Sigma }^{}\mathrm{\Sigma })$$
(5)
where $`P_\mathrm{\Sigma }^{}`$ and $`P_{\mathrm{\Sigma }^o}`$ denote the polarizations respectively of the $`\mathrm{\Sigma }^{}`$ and $`\mathrm{\Sigma }^o`$ before their decays, and $`C(\mathrm{\Sigma }^{})`$, $`C(\mathrm{\Sigma }^o)`$ and $`C(\mathrm{\Sigma }^{}\mathrm{\Sigma })`$ denote the model-independent functions of Clebsch-Gordan coefficients describing the ratio of the polarization of the final $`\mathrm{\Lambda }`$ to the polarization of the decaying baryon. The explicit values of these functions are given in ref. , where it is shown that the polarization of the final $`\mathrm{\Lambda }`$ in all models for the baryons and the dynamics of the decay process is given by the polarization of the strange quark in the simple constituent quark model for the decaying baryon. We immediately note that only the polarization of the directly-produced $`\mathrm{\Lambda }`$ can depend upon the spin-flavor structure of the $`\mathrm{\Lambda }`$. The other terms in eq.(1) depend upon the spin-flavor structure of the $`\mathrm{\Sigma }^{}`$ and the $`\mathrm{\Sigma }^o`$, but are independent of the spin-flavor structure of the $`\mathrm{\Lambda }`$.

The expression eq.(1) for the polarization of all the $`\mathrm{\Lambda }`$’s produced in a given experiment is easily generalized to obtain the polarization of $`\mathrm{\Lambda }`$’s restricted to a given domain of various kinematic variables. It is necessary to note that the momenta of $`\mathrm{\Lambda }`$’s produced from a decay of a $`\mathrm{\Sigma }^{}`$ or $`\mathrm{\Sigma }^o`$ are different from those of the parent baryon. Thus to obtain the polarization of $`\mathrm{\Lambda }`$’s produced in a given kinematic range one must integrate the expressions for $`N_\mathrm{\Sigma }^{}`$ and $`N_{\mathrm{\Sigma }^o}`$ over momenta with the appropriate weighting factors and angular distributions needed to produce the $`\mathrm{\Lambda }`$’s in the correct kinematic range. All this cannot be done reliably in present theoretical calculations and has to be taken from Monte Carlo simulations that are tuned and tested and thus reproduce well many related experimental observables.

The polarization of $`\mathrm{\Lambda }`$ and $`\overline{\mathrm{\Lambda }}`$ produced at the Z-pole  was calculated using two ingredients: the polarization of the $`s`$ and $`\overline{s}`$ quarks produced at the pole and their hadronization into $`\mathrm{\Lambda }`$ and $`\overline{\mathrm{\Lambda }}`$ . The first part was derived from weak interactions  and can be predicted directly. The hadronization of the $`s,\overline{s}`$ directly to $`\mathrm{\Lambda }`$ and $`\overline{\mathrm{\Lambda }}`$ or as decay products was calulated using Monte Carlo simulations. The authors found that about 20% of the $`\mathrm{\Lambda }`$ polarization is contributed by $`\mathrm{\Lambda }`$ particles originating from $`\mathrm{\Sigma }^{}`$ decay. These calculations used the naïve quark model wave functions for the baryons. This shows that inclusion of this contribution was essential for obtaining agreement with the data. In any attempt to go beyond the naïve quark model and include the information obtained from DIS on the spin-flavor structure of the proton, it is clearly necessary not only to have a model for the spin-flavor structure of the $`\mathrm{\Lambda }`$, but also of the $`\mathrm{\Sigma }^{}`$ and the $`\mathrm{\Sigma }^o`$ as well.

In studies of $`\mathrm{\Lambda }`$ polarization in deep inelastic scattering of polarized muon  the $`x_F`$ dependence of this contribution is studied (see fig. 1). It is found that up to $`x_F0.5`$ the contribution from $`\mathrm{\Sigma }^{}`$ decay to the $`\mathrm{\Lambda }`$ polarization is dominant. It is only for larger $`x_F`$ values that polarization from directly produced $`\mathrm{\Lambda }`$ ’s becomes dominant. It should be noted that the $`\mathrm{\Lambda }`$ yield is dropping fast for large $`x_F`$ and consequently most of the data is taken in the region where the $`\mathrm{\Lambda }`$ polarization is dominated by $`\mathrm{\Sigma }^{}`$ decay. The contribution from this process must therefore be considered very carefully. This contribution was addressed only in some of the theoretical calculations of this process .