text:
# Actions for dynamical fermion simulations: are we ready to go?

## 1 Introduction

No doubt, the lattice approach to quantum field theory has established itself as an integral part of high energy theoretical physics and has produced many results that are relevant to understand or even construct quantum field theories, as well as for phenomenology helping to interpret experimentally obtained data. Still, we are by far not at the end of the story. Serious dynamical fermion simulations with realistic number of quarks, i.e. two light and one heavy flavour, have just started and reveal already a number of unexpected problems. It is the objective of this presentation to address these problems and discuss their relevance and possible consequences. Nevertheless, we have made a lot of progress and we have understood much better to keep under control systematic errors such as discretization and finite size effects, the chiral limit and non-perturbative renormalization. However, still more patience has to be invested in our lattice computations to tell our colleagues from experiment trustworthy numbers. New computer architectures are essential and will help tremendously to step forward. Of course, we are all hoping for a miraculous invention of a new algorithm that would solve all our problems at once, although reality tells us that we will proceed in only small steps.

## 2 Conceptual questions

Certainly, for a non-lattice person (and even for a lattice practitioner) the plethora of actions that are used in todays lattice simulations is, to say it moderately, confusing. The initial Wilson fermion action  got many friends such as Symanzik improved actions  employing various gauge actions such as the Wilson , various forms of Symanzik improved gauge actions  and renormalization group inspired actions . A recent development concerning Wilson fermions is the formulation of twisted mass QCD . Of course, there is the approach of staggered fermions  which comes again in improved forms  and with different gauge actions. Newer developments are chiral invariant lattice fermions such as domain wall , overlap  and perfect actions . If exact chiral symmetry is not insisted on, but only chirally improved fermions are considered, designer actions such as FLIC , Hypercube  or truncated fixed point fermions  enter the stage.

For me, the criterion for choosing an action should be guided by the demand that

* the results obtained with these actions agree among themselves in the continuum limit and thus provide valuable cross checks
* and if using a particular action will not waste computer time.

There are a number of advantages and disadvantages of the above mentioned actions, concerning field theoretical conceptual aspects, lattice artefacts, chiral and renormalization properties and simulation cost that will be addressed below. Let me start by discussing locality properties of the actions used.

### 2.1 Rigorous actions

The list of the above actions can be ordered by the level of rigor they obey in constructing a lattice quantum field theory with the emphasis on their locality properties. The initial, standard Wilson actions (including the twisted mass formulation without clover term) have been shown to be reflection positive , they have a positive definite transfer matrix  and thus obey the reconstruction theorem (see first two refs. in ) guaranteeing thus the reconstruction of the Minkowski Green’s functions. Naive staggered fermions are on the same level  with the peculiarity that the transfer matrix has to be constructed over two lattice spacings.

### 2.2 Non-rigorous but local actions

The next class of actions are the ones that are ultra-local, i.e. they have an interaction range of only a few lattice points, but for which reflection positivity cannot be proven. To this class belong the Symanzik improved and the designer actions mentioned above. However, there cannot be a serious doubt that these actions are theoretically sound as the improvement terms added are irrelevant and vanish in the continuum limit. We might speak of “physical positivity”.

### 2.3 Exponentially localized actions

The chirally invariant lattice actions even give up ultra-locality  and only show an exponential localization , as shown for domain wall , overlap  and perfect fermions . Again, we have no reason to doubt that these actions define a local quantum field theory and describe QCD in the continuum limit.

### 2.4 Non-local actions (?)

Every time a lattice Dirac operator is constructed, where the action ranges over all points on the lattice, the question of locality has to be asked anew. An example for such an operator are improved staggered fermions for exact, say, two flavour simulations . It is not clear, how such a staggered fermion operator can be defined and also in this contribution no solution will be provided. The suggestion that is used in exact dynamical simulations, is to take the square root of the staggered operator , implemented by some polynomial approximation, which clearly defines an action that spreads over all lattice points. It should hence be tested, whether this operator is at least exponentially localized (see my plenary talk  for such a test on a Wilson operator).

There is a clear danger that through the square root operation a non-locality is introduced111There might be additional shortcomings of such a definition of a two flavour theory such as the violation of the optical theorem. However, it is not clear (at least not to me), how such a violation can be detected in practical simulations..

For a fixed set of parameters, i.e. fixed pion mass and fixed coupling, it should be expected that the localization of the operator is exponential, with couplings decaying like $`e^{r/r_{\mathrm{local}}}`$ with $`r_{\mathrm{local}}`$ the localization range. In the continuum limit, staying at a fixed value of $`r_0m_\pi `$, in principle two outcomes are possible. In scenario I, $`r_{\mathrm{local}}m_\pi =\mathrm{constant}\mathrm{for}a0.`$ In this case we would obtain a continuum field theory with a localization range of the order of the Compton wavelength of the pion. This would be clearly un-acceptable. In order to save such a theory at least to some extent, it would be necessary that the localization range is much smaller than the Compton wavelength of the heaviest physical particle of the target continuum field theory.

In scenario II, $`r_{\mathrm{local}}m_\pi 0\mathrm{for}a0.`$ In this case we would get a local continuum quantum field theory. One may suggest a “wishlist”, i.e. to

perform a check of the localization properties of the square root staggered fermion operator.

Such a work is in progress . We assume first, by lack of an alternative, that the two flavour staggered operator $`M_{\mathrm{staggered}}`$ is defined by taking the square root of the four flavor staggered operator. Second, we think of an exact simulation which will employ this square root operator. In particular, following the arguments for deriving a strict locality bound of ref. , we expect the localization range to satisfy the bound $`r_{\mathrm{local}}m_{\mathrm{quark}}=\mathrm{constant}`$. This behaviour follows directly from considering admissible gauge field configurations, the knowledge of the maximum eigenvalue of the staggered operator and the fact that the smallest eigenvalue is determined by the square of the bare quark mass. It remains to be seen, whether this theoretical bound is satisfied in practical simulations or whether it grossly over-estimates the localization range as it was found in the overlap case .

### 2.5 Perturbation theory

For most of the actions just discussed, despite their sometimes complicated nature, perturbation theory has been worked out. There are famous seminal papers for staggered fermions, see the first ref. in  and , as well as for Wilson fermions . These works build an important cornerstone to understand many aspects of lattice field theory such as renormalization. For more information on perturbation theory I refer to the recent review  and the talk by H. Trottier at this conference .

The only aspect of perturbation theory I would like to mention here is the Reisz power counting theorem . This theorem guarantees that lattice integrals that appear in lattice perturbation theory exist in the continuum limit. The theorem is satisfied for Wilson fermions. However, staggered fermions do actually not obey one of the conditions of the theorem (see ). This does certainly not imply that there is a problem for staggered fermions in perturbation theory. However, for a consistent and well founded perturbative discussion of staggered fermion it would be very helpful to

construct a “Reisz theorem” for staggered fermions.

Such a theorem for staggered fermions might also be very useful to analyze and underpin improved staggered actions, possibly leading to a disentanglement of the tastes in perturbation theory.

## 3 Algorithmic questions

Algorithms are the backbone of simulations in lattice field theory and testing and improving the algorithms cannot be appreciated enough. A review of existing algorithms for dynamical fermion simulations can be found in . Promising new ideas can be found in ,  and , see also .

The problem with our dynamical fermion simulations is the “cost wall”. We know to some extent the scaling behaviour of the algorithms and one formula for the cost $`C`$ in the case of Wilson fermions according to ref.  is

$$C=A\left(\frac{m_\pi }{m_\rho }\right)^{z_\pi }L^{z_L}a^{z_a}.$$
(1)
Here the exponents are given by $`z_\pi =6`$, $`z_L=5`$ and $`z_a=7`$; the value of the constant $`A`$ can be found in . Similar formulae are found by other groups . Although this cost formula might still be plagued by a rather large uncertainty, it provides a good measure for the costs of the algorithms employed.

In fig. 1 we show a comparison of simulations with staggered (inexact, i.e. using the $`R`$-algorithm ) and Wilson fermions. Several things are remarkable in this graph. The first is that for both kind of fermions at some value of the quark mass the costs explode. This is the famous cost wall where the simulations become essentially un-feasible. The second point is that for staggered fermions this point is reached at much lower values of the quark mass than for Wilson fermions. For the left plot we have used measured cost data from improved staggered fermion simulations  and the cost formula eq. (1) for Wilson fermions. A third observation is that the functional dependence of the cost on the quark mass seems to be the same for both actions, at least, given the accuracy of the data used here. The right plot shows again a cost comparison at a lower value of the lattice spacing, exhibiting a frightening need of teraflop years. The apparent advantage of the simulation cost using improved staggered fermions makes it even more important and urgent to test the theoretical basis of this approach as discussed in section 2.

One shortcoming of the staggered fermion simulations so far is that only a non-exact algorithm is used . In principle we know how to execute an exact odd flavour algorithm . A comparison of the in-exact algorithm and the exact one is shown in fig. 2. A striking effect of this graph is that the extrapolation to zero step size as needed in the in-exact simulation is non-monotonic. Hence such an extrapolation has to be done with very great care using many simulation points. Since exact odd flavour algorithms are on the market, it should be urged that they are indeed used in the staggered fermion simulations.

As a conclusion, it would be very useful to

use and test exact algorithms and to perform a fair cost comparison of exact algorithms including the approach to the continuum limit.

The actions discussed in section 2 are sometimes very complicated, including often fattening of links. A good news is that methods have been developed that allow for simulations of such actions either by reject/accept steps  or by computing the force explicitly  using special projection methods.

As a last topic in this algorithm section I would like to discuss a problem of principle. Random matrix theory (RMT)  predicts the eigenvalue distribution for the lowest non-vanishing eigenvalue $`\lambda `$ in topological charge sector zero to be $`P(\lambda )=\frac{z}{2}e^{\frac{1}{4}z^2}`$ with $`z=\lambda \mathrm{\Sigma }V`$ and $`\mathrm{\Sigma }`$ the infinite volume scalar condensate . This prediction is confirmed in practical simulations .

Such a distribution predicts the occurrence of very small eigenvalues to be frequent, and hence close to the critical point huge fluctuations have to be expected, rendering the simulations very expensive. Of course, a non-vanishing quark mass will regulate the problem, but it can do so only to a limited amount when the physical point at realistic values of the quark mass is reached. The conclusion is therefore that we have to expect that dynamical fermion simulations are becoming even more problematic than the cost formula of eq. (1) would imply when approaching small quark masses. As a result, the need to make contact with chiral perturbation theory is most urgent.

## 4 Scaling in the quenched approximation

Before turning to any results of unquenched simulations, let me discuss scaling in the quenched case first. Surprising results were presented by Aoki  at Lattice2000. He used different actions: staggered fermions with plaquette action, Wilson and improved Wilson fermions with plaquette and renormalization group inspired gauge actions. His results for $`m_N/m_V`$ for the different actions as a function of $`m_{\mathrm{PS}}/m_V`$ (i.e., the Edinburgh plot) did not agree in the continuum limit even for for heavy quark masses. The staggered results disagreed with the Wilson results, which agreed among themselves. Surprisingly, there was not much response to this striking and worrisome observation afterwards.

The question arises, whether something could be wrong with the continuum limit of improved staggered fermions. In Aoki’s article, a particular form of the continuum extrapolation was taken, linear in the lattice spacing for Wilson type fermions and an $`a^2`$ dependence for the staggered results. In fig. 3 I show an alternative fitting procedure, allowing for an additional $`a^2`$-dependence for the Wilson plaquette data222I leave out the improved Wilson data since they are not very consistent among themselves and it would certainly be important to generate a new set of data..

Although this way of fitting the data might not be the final truth, it provides a possibility with the conclusion that there is nothing wrong with the continuum limit of staggered fermions. The continuum value for $`m_N/m_V`$ is also consistent with results for various other improved actions as included in the plot. It would be very important, to supplement more and better data for such a plot and to

perform a precise scaling analysis for various fermion actions in the quenched approximation.

Such a test would be very re-assuring that all the calculations that are done with the many different actions lead to consistent values in the continuum limit and can hence lead to a good control over the systematic errors. For attempts for such a precise scaling tests I refer to the figure (done together with J. Zanotti) in my plenary talk  and the talks by A. and P. Hasenfratz at the conference .

## 5 Problems with dynamical fermion simulations

I now turn to somewhat surprising problems that are already encountered in practical simulations with dynamical fermions performed as of today. In the light of these problems, it may not be wise to start unthinkingly large scale and very expensive simulations with dynamical fermions but rather to concentrate on these problems for a while and try to understand and solve them.

### 5.1 Gauge actions

As mentioned before, improved gauge actions are not reflection positive. This leads to complex eigenvalues in the transfer matrix . In fig. 4 the static potential $`V(r_0)`$ is shown at a distance of the hadronic scale $`r_0`$  as function of the euclidian time. The values of $`V(r_0)`$ for the improved gauge actions approach their (in time) asymptotic value from below, indicating the effects of the complex eigenvalues . Although for asymptotic times the value of the potential can be determined, this suggests that for these actions larger time extents would have to be used. This can have consequences in particular for determinations of glue ball masses from such actions.

From a free field analysis  it also appears that different gauge actions might lead to very different lattice artefacts. In addition, the lattice spacing dependence need not to be very smooth approaching the continuum limit thus driving a reliable continuum extrapolation difficult. It might be useful therefore to perform simulations with different actions to check the continuum limit. An example, from my own experience, for such a procedure can be found in ref.  where simulations using different fermion actions were performed to extract moments of parton distribution functions. We will come back to this point towards the end of this proceedings.

Another, last difficulty of improved gauge actions is that they are very inefficient to sample topological charge sectors . Although this is not too problematic for quenched simulations since there many configurations can be generated fast, it may have serious consequences for dynamical fermions when such actions are employed.

### 5.2 Wilson fermions

The CP-PACS collaboration started $`N_f=3`$ simulations with improved Wilson fermions and the plaquette gauge action. By performing thermal cycles333In this standard procedure to detect phase transitions, one starts at some value of $`\kappa _{\mathrm{start}}`$, increases the value of $`\kappa `$ by some $`\mathrm{\Delta }\kappa `$ until a value of $`\kappa _{\mathrm{final}}`$ is reached and then runs back to $`\kappa _{\mathrm{start}}`$., clear indications of a first order phase transitions were detected as can be seen in fig. 5 where a pronounced hysteresis effect shows up. A similar effect appears also for $`N_f=2`$ flavours  if again the plaquette gauge action and non-perturbatively improved Wilson fermions  are used. It might be speculated that these phase transitions are related to the phase structure of the fundamental-adjoint pure gauge action. Working in the vicinity of such a phase transition can lead, e.g., to unexpectedly large lattice artefacts that may render the continuum extrapolation difficult.

Using the Iwasaki action  or the Symanzik gauge action the signs of the first order phase transitions seem to disappear completely . However, it is not clear, whether such effects eventually will strike back. Therefore, it would be important to

study and understand the $`T=0`$ phase diagram of QCD, to investigate different gauge actions and to understand the nature of the phase transition.

In fig. 6 we show a comparison of the lattice artefacts in the quark mass for quenched and unquenched simulations  using non-perturbatively improved Wilson fermions. The figure indicates that for dynamical fermions the artefacts are substantially larger, although they are still smaller than using naive Wilson fermions . Whether this increase of lattice artefacts is due to the possible phase transition discussed above or originates from a completely different source, is an open question and gives additional motivation to understand the phase diagram better.

During the conference it has been pointed out  that in dynamical simulations using Wilson or staggered fermions the hadronic scale $`r_0`$ depends rather strongly on the quark mass at fixed values of $`\beta `$. Clearly, this effect has to be taken into account when quantities are analyzed as functions of, say, $`r_0m_{\mathrm{PS}}`$. This becomes especially important when comparisons with chiral perturbation are performed. In addition, it was shown that in the quantity $`r_0(m)/r_0(m_{\mathrm{ref}})`$ with $`m_{\mathrm{ref}}`$ determined by $`m_{\mathrm{PS}}/m_\mathrm{V}=0.45`$, staggered fermions exhibit surprisingly large lattice artefacts. Clearly, it would be important to check and clarify this finding further with results from improved staggered fermion simulations.

### 5.3 Staggered fermions

For simulations with staggered fermions, naive or improved, there is a problem with the topological charge identification when coarse lattices at large values of the coupling are used. This can be seen quite evidently in studies that compare the predictions of the eigenvalue distributions predicted by RMT with the results of numerical staggered fermion simulations. Data obtained in topological charge sector $`1`$ are fitted with the topology $`0`$ prediction of RMT, while the curves for topology $`1`$ fail completely . These studies were mainly done with unimproved staggered fermions. For newer results with improved staggered fermions, see .

In  it was argued that this failure of staggered fermions with respect to questions related to topology has to occur when the gauge fields are not smooth enough and one is not working at small enough values of the coupling. The reason is that with such rough gauge fields, the taste breaking effects are strong and the symmetry breaking pattern is not the one of the continuum.

Indeed, in the first reference of  where simulations of the Schwinger model have been performed, there are indications that for larger values of $`\beta `$ the predictions of RMT and the numerical data get reconciled. Nevertheless, it would be very important

to check explicitly that the predictions of RMT and other properties of topology are reproduced by, preferably, improved staggered fermions.

## 6 Results for dynamical Wilson fermions

This section is devoted to some results that have been obtained recently with dynamical Wilson fermions. Results concerning staggered fermions can be found in  and I will only shortly discuss staggered fermions in the context of chiral perturbation theory here.

### 6.1 Chiral perturbation theory

The algorithms we are using at the moment will sooner or later hit a cost wall where the exponential growth of the cost will set in with full power when the quark mass is lowered below a threshold value. In order to obtain results for physical values of the pion and $`\rho `$ masses it becomes therefore mandatory to make contact to chiral perturbation theory ($`\chi `$PT) and then use $`\chi `$PT to perform the extrapolations to the “physical point”. Finding out the window where on the one hand lattice simulations can be performed and, on the other hand, $`\chi `$PT can be controlled and higher order effects can be safely neglected, is of central importance in lattice QCD and much effort will be invested to resolve this question. Fig. 1 suggests that improved staggered fermions may reach this overlap region sooner than Wilson fermions.

There are two main strategies to confront lattice results with predictions of $`\chi `$PT:

* the first strategy is to extrapolate lattice data to the continuum limit first and then compare to $`\chi `$PT;
* the second strategy is to take the lattice artefacts into account in $`\chi `$PT itself and compare then to lattice data.

Let me discuss the first strategy first. The big advantage of this procedure is that in the continuum chiral symmetry is restored and therefore a direct comparison to $`\chi `$PT becomes possible444Of course, using chiral invariant lattice fermions would also fulfill this requirement, but, they are, unfortunately, very expensive to simulate.. Clearly, the bad side of the coin is that performing the continuum limit is demanding and needs large computer resources.

This computational demand is at present “tricked out” by taking lattice data at only one value of the lattice spacing. Of course, the improved lattice actions in both the Wilson and the staggered cases substantially reduce the lattice artefacts. However, only future simulations can tell, whether the lattice data were indeed close enough to their continuum values.

In ref.  a recent analysis of chiral perturbation theory applied to lattice data are given. Older discussions are numerous and I have to refer to , the plenary discussion at Lattice 2002 and  for other examples (taken from the analysis of moments of parton distribution functions, my own hobby-horse).

In the formulae used to perform the analyzes in , the parameters of the chiral Lagrangian were fixed by taking experimental input values obtained at the physical point. In the calculation chiral symmetry was consistently kept and an “improvement term” to cancel cut-off effects was added. The numerical data are taken from improved Wilson CP-PACS simulations and are only at one value of the lattice spacing. Fig. 7 shows a comparison of the predictions from $`\chi `$PT in this setup and the data, which seems to be not so bad, even for quite large pion masses. However, in ref.  itself it is claimed: We stress again that applying the expressions to pion masses above 600 MeV is only done for illustrative purposes, for a realistic chiral extrapolation smaller pion masses are mandatory.

The reason for this becomes clear by looking at the size of corrections that can appear in a 4-loop calculation. These corrections can be very large, depending on the precise form of performing $`\chi `$PT. I do not have the space (and it is not the goal of this article) to discuss in detail which kind of $`\chi `$PT is correct, what flaws there might be and where conceptual difficulties appear. However, I believe, that finally it has to be pure $`\chi `$PT that has to be used, without any additional modeling or assumptions. Only then can we safely employ $`\chi `$PT to do for us the job of extrapolating to the physical point.

The second strategy is to take possible lattice artefacts into account in $`\chi `$PT itself. This idea was put to the fore by S. Sharpe  and was then further developed in  for Wilson and in  for staggered fermions. The advantage of this approach clearly is that chiral symmetry violating effects, taste violating effects in the case of staggered fermions and even different actions can be taken into account from the beginning. The disadvantage is that the corresponding chiral Lagrangian has twice as many low energy constants, and that the additional parameters depend on the bare coupling $`g_0`$ leading to a plethora of fit parameters.

This approach has been used already in practical simulations. For results with improved staggered fermions I refer to the contribution of S. Gottlieb at this conference . The only remark, I want to add here is that $`\chi `$PT provides a justification of taking the square root of the fermion determinant. This situation corresponds to a partially quenched setup in which two quarks are taken to be quenched . The idea can be generalized using the replica trick: take arbitrary numbers $`N_u`$, $`N_d`$, $`N_s`$ for the $`u`$, $`d`$ and $`s`$ quarks. Perform the computation in this general form and set in the end $`N_u=N_d=N_s=1/4`$ . In , although in a different context, it has been demonstrated that such a procedure is, at least to 1-loop order, completely equivalent to the super symmetric method and, although a proof to all orders is still missing, there is a good chance that the argument might work.

After this excursion to staggered fermions, let me now turn back to Wilson fermions. There are two works at the moment where simulation data are confronted to the predictions of $`\chi `$PT when lattice artefacts are included in $`\chi `$PT. The first work is ref.  and the second is ref. . We show in fig. 8 a plot of $`m_{\mathrm{PS}}^2/m_q`$ as a function of the quark mass at various values of $`\beta `$. For each value of $`\beta `$ two curves are shown. The curves showing a very strong curvature represent continuum $`\chi `$PT while the lines that fit the data represent $`\chi `$PT with lattice artefacts taken into account. Clearly, continuum $`\chi `$PT is not able to describe the simulation data while lattice $`\chi `$PT fits the data well. Another observation is, however, that very large lattice artefacts in $`m_{\mathrm{PS}}^2/m_q`$ seem to exist. Data at various values of $`\beta `$ come out to be quite different.

Another comparison is shown in fig. 9 (from ). Here a ratio of decay constants of the type, $`f_{VS}^2/(f_{VV}f_{SS})`$ is plotted with $`V`$ and $`S`$ denoting valence and sea quantities. In this double ratio, which is only one example of the many fits that can be found in ref. , a good agreement with the predictions of $`\chi `$PT is seen. However, there seem to be tremendous cancellations of cut-off effects, keeping in mind that the simulations have been done at a value of the lattice spacing of $`a=0.28`$fm.

Another, somewhat worrisome observation is that the universal $`\mathrm{\Lambda }`$-parameters come out to be very different, depending on the value of $`\beta `$ that is used (in case of ) or how the lattice spacing is defined (in case of ). Clearly, we are just in the beginning of using the new approach of $`\chi `$PT where lattice spacing effects are taken into account and the first results, using this method, are somewhat puzzling. In particular, I would like to advocate to de-double the double ratios in order to be able to disentangle the cut-off effects and show their size, furthering thus our understanding of $`\chi `$PT and the connection to lattice QCD.

## 7 Numerical results

In this section, I would like to discuss a few selected results from recent dynamical Wilson fermion simulations.

### 7.1 Finite size effects

It seems that we finally see the exponential finite size effects of hadron masses. The general form of these finite size effects are $`M(L)M=\frac{3}{16\pi ^2ML}_{\mathrm{}}^{\mathrm{}}F(iy)e^{\sqrt{M_\pi ^2+y^2}L}`$ where $`F`$ is the $`\pi \pi `$ forward scattering amplitude in infinite volume . A calculation of the integral leads to the exponential finite size correction $`M(L)ML^{3/2}e^{m_\pi L}`$. Corrections to this formula have been computed within $`\chi `$PT . An alternative calculation was performed and also tested in ref. . At this conference there were two contributions  that showed that indeed the finite size effects are exponential and not power like of the form $`M=m_{\mathrm{}}+c/L^3`$ as anticipated in ref. . In fig. 10 a comparison of numerical simulation data to various fits of the finite size effects are shown. It seems that indeed the exponential form of finite size corrections are preferred by the data confirming thus the exponential finite size effects in lattice simulations. A similar picture and conclusion can be found in . Clearly, these first results for the finite size effects need to be confirmed in future simulations.

### 7.2 Fundamental parameters of QCD

The running coupling and running quark masses are basic and important parameters of QCD that can be extracted from lattice simulations. The status of these calculations were reported at this conference . Although the running itself has been determined already to a quite high precision, the determination of the physical scale is still missing. In addition, it would be necessary to complement the present simulations with even larger lattices to keep the continuum limit under better control. Both of these goals will certainly be done in the near future. In the past, the Schrödinger functional, which was used for the scale dependent renormalization , was the domain of the ALPHA collaboration. Nowadays, there are new groups who use this renormalization scheme. In particular, a serious attempt for $`N_f=3`$ flavours of dynamical fermions has been started by the CP-PACS and JLQCD collaborations .

The CP-PACS and JLQCD collaborations seem not only to join their forces to compute the coupling and quark masses but also for other quantities. This is a convincing indication that dynamical fermions are a real challenge and a very difficult problem. Present simulations employ the Iwasaki gauge action of ref.  combined with the Symanzik improved Wilson fermion action. The reason for this choice is the occurrence of a phase transition when the standard gauge action is employed as mentioned above. The value of $`c_{\mathrm{SW}}`$ is the non-perturbatively improved one .

In the $`N_f=3`$ simulations, the physical extent of the lattice is about $`L1.6`$ fm. With such a lattice, it is not possible to perform simulations for baryonic quantities reliably since the finite size effects will be very large . Hence, at the moment one concentrates on mesons. In fig. 11 the relative difference of the K and $`\mathrm{\Phi }`$ meson masses as compared to the experimental data is shown. The open triangles represent the quenched approximation where for different physical input value a large variation is found. It is very reassuring that for dynamical fermions this intrinsic uncertainty of the quenched approximation seems to be eliminated completely.

Similar to what has been said above for the quenched approximation, it would be very interesting to

combine Wilson, staggered and other dynamical fermion action results to check the different approaches in the chiral, continuum and infinite volume limits.

### 7.3 Other topics

An article like the present one has to be selective and a number of topics had to be left out, unfortunately. In the first place, I could not discuss progress with dynamical domain wall  or even overlap  fermions. For domain wall fermions it seems that an appropriate gauge action has been determined where an additional adjoint piece is added which seems to solve some problems with unexpected large values of the residual mass. The simulations for producing physics have been started and I am sure that at the next conference we will see a lot of interesting results. I would also like to mention some works that discuss the study of violations of CP  and even CPT  without being able to go into any detail.

I could also not go into results for moments of parton distribution functions from dynamical simulations , nor could I discuss questions concerning topology in depth . Attempts to tackle the notoriously difficult problem of the $`\eta ^{}`$ meson by making use of the eigenvalue spectrum can be found in . My personal opinion here is that this is an interesting approach, but it still has a number of open questions that should be addressed in the future.

## 8 Conclusion

Obviously, there are dangerous animals on the lattice . Our present dynamical fermion simulations are still in a rather early stage but they meet already such dangerous animals: unexpected phase transitions, questions of field theoretically non-local actions, actions with unitarity violations, occurrence of large lattice artefacts, algorithms that become extremely costly when approaching the physical point, uncertainties in the continuum limit of different lattice actions, difficulties with the chiral extrapolations.

In the real world, we have learned to tame lions and tigers and in lattice field theory this will be done for sure, too. However, we have to spend work for achieving this: we have to explore and understand the phase diagram of (zero temperature) QCD, we have to study the locality properties of taking the square root, we have to explore better the phase diagrams of extended gauge actions and understand the properties of such actions. We should perform scaling tests first in the quenched approximation to demonstrate to ourselves and to the rest of the world that universality holds (with all the caveats of the notion of universality in the quenched approximation). Participate in algorithm development. It is here, where a breakthrough can happen which would solve the cost-wall problem we have with existing algorithms; it would also be very helpful to join the ILDG  initiative to combine our efforts and exchange configurations.

One recommendation could be to not put all cards on one particular action. In fact, it would be wise to first select carefully actions that do not show the flaws and problems mentioned in this contribution. Even after such identifications have been made, it may be worthwhile to perform simulations with two (good) actions to check for systematic errors.

Another recommendation could be to use actions tailored for the physics problem one is interested in. If problems with rather heavy quark masses are considered, (theoretically sound) staggered fermions or improved Wilson fermions are a good choice. Going lighter, designer actions such as FLIC, hypercube, truncated perfect actions or domain wall fermions with rather small extra dimension might enter the game. Finally, for really light quarks, chiral invariant (up to machine precision) domain wall and overlap fermions should be used.

## Acknowledgments

When preparing this talk and delivering it in Tsukuba, I felt like in front of a minefield . I thank all colleagues who helped me with discussions, suggestions, sending material and pointing out pitfalls to cross the minefield without stepping on the mines: S. Aoki, C. Bernard, M. Creutz, P. Damgaard, C. Davies, S. Dürr, M. Golterman, S. Gottlieb, A. Kronfeld, M. Lüscher, I. Montvay, R. Sommer, T. Wettig.

F. Knechtli, S. Gottlieb, M. Hasenbusch and J. Zanotti helped me directly with the talk by preparing figures or even running some simulations. I also thank K.I. Ishikawa for lending me his laptop for the time of the conference and T. Kennedy for giving me his for the presentation. Last but not least I want to thank B. Willems for the minefield picture  that I would never have been able to draw myself.


format: text

text:
# Formation of First Stars Triggered by Collisions and Shockwaves: Prospect for High Star Formation Efficiency and High Ionizing Photon Escape Fraction

## 1 Introduction

The Wilkinson Microwave Anisotropy Probe (WMAP) polarization observations measured Thomson optical depth $`\tau _e=0.17\pm 0.04`$ ($`68\%`$) (Kogut et al. 2003), requiring a full reionization to have occured first at $`z=20_9^{+10}`$, such as in the double reionization picture (Cen 2003a). This high Thomson optical depth is further confirmed by a variety of combined data sets (Tegmark et al. 2003). To ionize the universe at such a high redshift demands both a high formation efficient ($`ϵ`$) for Population III stars and a large escape fraction ($`f_{esc}`$) for the produced ionizing photons. Quantitatively, one would require the product of the two to be greater than $`0.01`$ in order to produce $`\tau _e>0.11`$ (e.g., Cen 2003b; Haiman & Holder 2003; Wyithe & Loeb 2003; Ciardi, Ferrara & White 2003; Ricotti & Ostriker 2003) to be consistent with the WMAP observations at $`1.5\sigma `$ level. Conventional wisdom may view these assumptions somewhat too optimistic at first sight.

In this Letter we suggest that star formation at high redshift ($`z>10`$) may occur in massive starbursts, triggered by supersonic collisions between constituent cold atomic clouds and blastwaves produced by explosions of first stars. The consequences are (1) high star formation efficiency and (2) high ionizing photon escape fraction.

## 2 Collisions between Cold Atomic Gas Clouds in High Redshift Galaxies

The slope of density power spectrum is $`3`$ in the cold dark matter cosmological models. As a result, galaxies over a wide range of mass, $`10^410^8\mathrm{M}_{}`$ of interest here at high redshift, are expected to form nearly simultaneously, when they form. This would imply that large galaxies of mass $`10^710^8\mathrm{M}_{}`$ (virial temperature $`T_v>10^4`$K at $`z>10`$), have most of their mass in small clouds, which are formerly minihalos. How much mass of a large galaxy is in minihalos? This is computable approximately using the Extended Press-Schechter (EPS) theory (Bond et al. 1991; Bower 1991; Lacey & Cole 1993, LC hereafter). We use the instantaneous merger rate for a halo of mass $`M_p`$ to merger with a halo of mass $`\mathrm{\Delta }M=MM_p`$ to form a halo of mass $`M`$ (LC):

$$R_N(M_pM)dM=(\frac{2}{\pi })^{1/2}\frac{1}{t}|\frac{d\mathrm{ln}\delta _c}{d\mathrm{ln}t}|\frac{1}{M}|\frac{d\mathrm{ln}\sigma }{d\mathrm{ln}M}|\frac{\delta _c(t)}{\sigma }\frac{1}{(1\sigma ^2/\sigma _p^2)^{3/2}}\mathrm{exp}[\frac{\delta _c(t)^2}{2}(\frac{1}{\sigma ^2}\frac{1}{\sigma _p^2})]dM,$$
(1)
where $`\sigma `$ is the density variance of a sphere in the linear density field extrapolated to $`z=0`$ that contains mass $`M`$; $`\sigma _p`$ is $`\sigma `$ for mass $`M_p`$; $`\delta _c(t)`$ is the critical density for collapse (equation 2.1 of LC); $`t`$ is age of the universe at the epoch examined. The total mass accretion/merger rate of halos of mass in the range $`M_1`$ to $`M_2`$ to form a halo of mass $`M`$ at $`t`$ is then

$$R_M(M,M1,M2)=_{MM_2}^{MM_1}R_N(M_pM)(MM_p)𝑑M_p.$$
(2)
Figure 1 shows the mass fraction of accreted objects in the mass range $`2\times 10^410^7\mathrm{M}_{}`$ by a large halo with $`M=1\times 10^8\mathrm{M}_{}`$ as a function of redshift. We see that $`5080\%`$ of all mass in large galaxies with virial temperature greater than $`10^4`$K at $`z>10`$ are in minihalos. The standard WMAP-normalized cold dark matter model (Spergel et al. 2003) is used: $`\mathrm{\Omega }_M=0.27`$, $`\mathrm{\Lambda }=0.73`$, $`\mathrm{\Omega }_b=0.04`$, $`H_0=70`$km/s/Mpc, $`n=0.99`$ and $`\sigma _8=0.90`$.

The constituent minihalos do not remain intact in the process of infall and virialization. First consider the effect due to gravitational tidal field of the large, host halo on minihalos. For simplicity we will assume that all dark matter halos are singular isothermal spheres truncated at the respective virial radii, $`\rho (r)=\rho _v(r/r_v)^2`$ for $`rr_v`$, where $`\rho _v=100\overline{\rho }`$, $`\overline{\rho }`$ is the mean density of the universe at that redshift and $`r_v`$ is the virial radius. For illustration we will consider a large halo of mass $`1\times 10^8\mathrm{M}_{}`$ at $`z=20`$ having a virial temperature of $`3.7\times 10^4`$K. We assume that $`2/3`$ of the total mass of the large halo were originally from minihalos (see Figure 1), and all minihalos have the same mass of $`1\times 10^6\mathrm{M}_{}`$ with a virial temperature of $`1700`$K, totaling $`67`$ such minihalos in the large halo. We denote the virial radius of the large halo as $`r_L`$ ($`=0.71`$kpc proper) and that of the minihalos as $`r_M`$ ($`=0.15`$kpc proper). Then, a minihalo at radius $`r`$ from the center of the large halo would have its matter exterior to $`r_T=r_M(r/r_L)`$ be tidally stripped off. The mean collision rate between minihalos (ignoring gravitational focus due to high speed) is

$$<n\sigma v>=\frac{_0^{r_L}vn_M(r)\sigma (r)n_M(r)4\pi r^2𝑑r}{_0^{r_L}n_M(r)4\pi r^2𝑑r},$$
(3)
where $`\sigma (r)=\pi r_T^2`$ is the cross section of a minihalo at radius $`r`$; $`n_M(r)`$ is the number density of minihalos at radius $`r`$; $`v`$ is the velocity of each minihalo. Assuming also an isothermal distribution of minihalos, $`n_M(r)=n_0(r/r_L)^2`$, we have

$$<n\sigma v>=\pi r_M^2vn_0.$$
(4)
This result is obtained only considering the gravitational tidal effect. Let us now consider the effect due to external inter-minihalo medium (IMM) on the internal gas of minihalos. Consider a minihalo located at $`r_L`$. The IMM with $`T=3.7\times 10^4`$K and a density of $`\rho _v/3`$ has a pressure higher than the minihalo gas with $`T=1700`$K and a density of $`\rho _v`$. Hence, the gas density in minihalo would be raised by adiabatic compression by the external gas by a factor of $`3.3`$ and the temperature increases by a factor of $`2.2`$, resulting in a pressure equilibrium (ignoring self-gravity of a minihalo for this simple calculation). The compression would reduce the cross section by a factor of $`2.2`$ (A better treatment will be using virial theorem with external pressure but is not given here). The percentage reduction in cross section at other radii is the same. Thus, the average collision time between minihalos is

$$t_C=\frac{2.2}{\pi r_M^2vn_0}=4.3\times 10^7\mathrm{yr},$$
(5)
where the minihalo clouds in the large halo are assumed to move at a velocity equal to $`v=\sqrt{3}\sigma _v=30\mathrm{k}\mathrm{m}\mathrm{s}^1`$ ($`\sigma _v=17.4\mathrm{km}\mathrm{s}^1`$ is 1-d the velocity dispersion of the large halo for our illustrated case). It is seen that the collision time among minihalos is about $`25\%`$ of the age of the Universe at $`z=20`$. The time it takes for a minihalo to cross the large halo is $`4.6\times 10^7`$yr. Thus, on average, a minihalo experiences a collision for every crossing time.

Gas within minihalos may be stripped off by ram pressure when traveling through the IMM of the large host halo. Applying the same physical argument of Gunn & Gott (1972), if the gravitational restoring force of a minihalo is exceeded by ram pressure force, gas and dark matter of the minihalo will be displaced, i.e., gas being stripped off. This may be expressed by the condition $`\rho _{IMM}v_0^2=\rho _Mv_M^2`$, where $`\rho _{IMM}`$ is the density of the IMM; $`v_0`$ is the traveling velocity of the minihalo; $`\rho _M`$ is the mean density of a gas parcel of the minihalo in question; $`v_M`$ is the 3-d velocity dispersion of the minihalo. We see that this condition is similar to the pressure equilibrium condition above (without considering the self-gravity of minihalos). Thus, ram pressure stripping is relevant and likely to act at the same time when the cloud is compressed by external thermal pressure. It is possible that gas in outmost layers of a minihalo may be subject to ram pressure stripping while most interior gas will not be dislodged. We will ignore the change to $`t_C`$ (equation 5) possibly caused by ram pressure stripping, because other effects such as deformation of minihalo gas clouds may easily dominate but are hard to characterize simply.

What would happen in a head-on collision between two clouds of cold gas? Smith (1980) and Struck-Marcell (1982a,b) have studied this issue of cloud collision inside galaxies for different ranges of Mach numbers using both analytic method and simulations, in somewhat different environments. Here we focus on the thermodynamic properties of post-shock gas, computed assuming that post-shock region is isobaric, while leaving detailed treatment to future 3-d hydrodynamic simulations. All cooling/heating processes for a primordial gas are incorporated, including atomic cooling, $`\mathrm{H}_2`$ formation and cooling and Compton cooling with the rates from Abel et al. (1997). Radiation background is assumed to exert no significant effect, as shown by Mac Low & Shull (1986) in the context of intergalactic shocks, which are physically similar to collision-induced shocks considered here. Figure 2 shows the evolution of the postshock gas following a head-on collision between two cold atomic gas clouds. Note that the time for the shock to propagate through an entire cloud is $`t_s=2r_M/(4v_0/3)=7.3\times 10^6`$yr. Since the shock heated cloud has already cooled to a much lower temperature and the Jeans mass has become much smaller than the total baryonic mass of the merged cloud at a time scale shorter than for the shock to transverse the entire cloud, the subsequent rarefaction wave following the thermalization (shock) process is insignificant, so are the compression and expansion phases following that. This indicates that the isobaric assumption is approximately valid and our simple calculation should have captured the essence of the process. We see that a molecular hydrogen fraction of $`8\times 10^4`$ is reached and gas is cooled to $`100`$K within a fraction of the Hubble time after the collision ($`10^610^7`$yr). The final Jeans mass is about $`100\mathrm{M}_{}`$, indicating that stars formed may be massive. One interesting difference in terms of star formation between this picture and that of Abel et al. (2001) is that multiple stars, possibly a star cluster, may form, which, unlike in a single star case of Abel et al. (2001), may provide an additional mechanism to remove the angular momentum of the collapsing gas through orbital coupling.

In summary, supersonic collisions between cold atomic gas clouds within large halos can trigger star formation and star formation efficiency may be quite high, because a large portion of the colliding clouds may be able to cool and form stars. The spatially distributed star formation will be much more favorable to the escape of ionizing photons.

## 3 Star Formation Triggered by Shockwaves

The possible spread in time among collisions may be $`10^710^8`$yr, longer than the lifetime of massive first stars of $`3\times 10^6`$yr. Thus, it is likely that collisions between minihalo clouds can not be considered coeval. Consider two possible scenarios.

In the first scenario, the first stars were so massive that they end their lives imploding and display no explosive events, for example, if the stars are more massive than $`300\mathrm{M}_{}`$ and non-rotating (Rakavy, Shaviv, & Zinamon 1967; Bond, Arnett, & Carr 1984; Glatzel, Fricke, & El Eid 1985; Woosley 1986), the subsequent star formation in collisions that follow would proceed normally. Star formation efficiency should be high and so is the ionizing photon escape fraction, since star formation will be spatially distributed and extended (see below).

In the second scenario, first stars end in explosions and launch a blastwave. Before computing the dynamics of blastwaves from the explosion of first stars, let us first consider the photoionization effect by the first stars on non-star-forming minihalos, before blastwave encompassing them. The ionizing photons from the massive stars formed in collisions will propagate outward from the minihalo. Since the ionization front (assuming an ionizing photon emissivity of $`1.6\times 10^{48}`$photons s-1 $`\mathrm{M}_{}^{}{}_{}{}^{1}`$; Bromm, Kudritzki, & Loeb 2001) travels at a velocity of $`v_{IF}=F/n=3478(d/1\mathrm{k}\mathrm{p}\mathrm{c})^2(ϵ/0.1)(n/0.1\mathrm{cm}^3)^1`$km/s (assuming no attenuation up to distance $`d`$ from the emitting source), much faster than the sound speed in photoionized gas, an R-type ionization front thus develops. Subsequently, a sound wave (not a shockwave) will develop driving the ionized gas outward at the sound speed, which would travel $`0.046`$kpc ($`r_M/3`$) in the lifetime of the massive stars. As a result, the gas density immediately surrounding the stars will be much reduced. Given that the expansion distance is comparable to the original minihalo radius one may expect that the final gas distribution surrounding the stars at the end of their lives may have a density comparable to mean density of the original minihalo. Without detailed simulations we resort to some order of magnitude estimates. Let us assume that the clumping factor of the gas in the star forming minihalo cloud is $`C`$ (in terms of pre-shock gas density), then it can be shown that the total amount of gas required to balance the ionizing photon emission from recombination is $`M_{rec}=1.0\times 10^5(e/0.1)(C/16)^1\mathrm{M}_{}`$ for an ionizing photon luminosity of $`4.2\times 10^{52}(e/0.1)`$sec-1, which should be compared to $`2.6\times 10^5\mathrm{M}_{}`$ of the total gas mass of the colliding clouds. The fact that the former is of the same order of the latter and the fact that the post-collision gas medium of the minihalos is likely to be quite irregular tell us that (1) a significant fraction of the ionizing photons is likely to escape the star-forming minihalos, and (2) the Strömgren sphere size is not much larger the minihalo itself, which would imply that the ionizing effect from the star-forming minihalos on other minihalos in the same galaxy may be relatively mild. In a fraction of the solid angle, $`f_\mathrm{\Omega }`$, from the star-forming minihalo optical depth may be small enough to allow for the ionizing photons to escape. For those directions, ignoring the recombination effect due to the IMM, it can be shown that the minimal distance from the star-forming minihalo for the recipient minihalo of ionizing photons to experience a D-type ionization front is $`d_{min}=2.2(\frac{\rho }{\rho _M})^{1/2}`$kpc, which should be compared to the virial radius of the large-host halo of $`r_L=0.71`$kpc. Thus, within the fraction of the solid angle where ionizing photons are able to escape into the IMM, it appears that ionization front is R-type initially at the virial radius of a minihalo but turns into a D-type front when penetrating deeper into the minihalo (at $`0.57r_M`$ for a minihalo at a distance of $`r_L`$). A D-type front may transmit a weak shock causing slight compression of the gas, while the regions cruised through by the R-type front may expand outward. In brief, we see that a small fraction of non-star-forming minihalos within the solid angle of escaped ionizing photons from star-forming minihalos may experience some dynamic effect due to photoionization, likely resulting in a very compact central region and evaporated outskirts (Barkana & Loeb 1999; Shapiro et al. 2003). Most of the minihalos are, however, not irradiated by ionizing photons from the first stars and would remain neutral and intact.

Now consider the effects of the blastwaves produced by the explosions of the first stars in the same large host galaxy. We will assume that the first stars each have mass of $`150\mathrm{M}_{}`$ with supernova explosion energy of $`3\times 10^{52}`$erg. A star formation efficiency of $`ϵ`$ would give a total explosive energy of $`5.2\times 10^{54}(ϵ/0.1)`$erg. It is not required that the mass of the first stars be in the range $`130300\mathrm{M}_{}`$ to induce pair-instability supernovae; our conclusions remain little changed if the mass of the first stars is the range of $`20130\mathrm{M}_{}`$. We will now estimate how far and at what speed the blastwave from such an explosion would travel. For simplicity and easy tractability we assume that the blastwave has two phases: the Sedov-Taylor phase followed by a momentum conserving phase (with an instantaneous cooling phase in-between). The assumed second phase is conservative in the sense that it ignores any pressure of possibly still hot interior gas on the shock shell. The formula for the first phase is adopted from Ostriker & McKee (1988) assuming a uniform density locally. Cooling rate for a primordial gas is taken from Sutherland & Dopita (1993) with Compton cooling added. Figure 3 shows the velocity of the blastwave after it has traveled a distance such that it would, on average, hit $`1`$ (solid curve) and $`8`$ (dashed curve) other minihalos, as a function of the ambient gas density (depending upon the radius). We see that for density in the relevant range the expected shock velocity sweeping minihalos may be in the range $`v_s=201000`$km/s. We note that, with significantly reduced density surrounding the exploding stars due to photoionization, the blastwave only starts cooling way after breaking out the minihalo gas clouds in which the stars were formed. Thus, neglecting the gas in the host minihalo does not introduce significant errors.

What would happen to the minihalos when shocks of this magnitude travel through them? A shock with a velocity of $`93`$km/s in the ambient medium would produce a shock of velocity $`30`$km/s ($`=93(\rho _{IMM}/\rho _M)^{1/2}=93/\sqrt{9.3}`$) upon entering a minihalo, producing results comparable to those shown in Figure 2. One important difference is that, in this case of encompassing blastwaves, the external thermal pressure is high and could aid in compressing the shocked and cooled minihalo gas clouds, which will likely enhance star formation in the minihalos. The reader may consult Kang et al. (1990) for a 1-d treatment in the context of globular cluster formation and Klein, McKee, & Colella (1994) for a recent calculation. In short, blastwaves from first massive stars could cause star formation in other minihalos in the same galaxy, yielding a widespread star formation. Most likely, one generation of stars may be able to trigger star formation in majority of the minihalos.

A higher Mach shock compresses the minihalo gas to a higher density, resulting in a lower Jeans mass (roughly inversely proportional to the Mach number). Thus, earlier collisions between the faster blastwave and minihalos could lead to formation of lower mass stars, even in metal-free gas. We suggest that the observed, extremely low metallicity, low mass star HE0107-5240 (Christlieb et al. 2002) may be formed during some early time, high Mach shocked minihalos, which are slightly contaminated with metals (see Hideyuki & Nomoto 2003).

One interesting consequence of sweeping blastwave is that a large bulk velocity will be imparted on the compressed minihalo clouds. The expected bulk velocity in this case is about $`3v_s/4\sqrt{9.3}=0.25v_s`$, where $`v_s`$ is the blastwave velocity in the ambient gas and the factor $`\sqrt{9.3}`$ is due to the density ratio of minihalo gas to the ambient gas. Thus, these minihalos may subsequently move at a velocity larger than the escape velocity of the large host halo ($`v_{esc}=35\mathrm{k}\mathrm{m}\mathrm{s}^1`$), if $`v_s140\mathrm{k}\mathrm{m}\mathrm{s}^1`$. Consequently, star formation may be widespread spatially and it is even possible that some stars may form after the minihalos have left the original host galaxy. One immediate implication is that ionizing photons from these stars will be much easier to escape, giving rise to escape fraction possibly well in excess of $`10\%`$, whereas in the conventional picture of central star formation where the escape fraction is estimated to be below $`1\%`$ (Wood & Loeb 2000; Ricotti & Shull 2001).

Another implication is that, if some stars have moved outside their original large host galaxy when they explode, the resulting blastwaves will travel much further thanks to the lower density in the intergalactic medium, which, in turn, could shock heat, accelerate and compress minihalos in the intergalactic medium to form stars. If this chain reaction can be triggered, star formation in minihalos, which are much more abundant (in total mass) than large halos and are otherwise unable to form stars, could produce many more ionizing photons than what can be produced by large galaxies alone. In Figure 3 the three solid dots show the velocity of the blastwave from a minihalo with an explosion energy of $`2.6\times 10^{54}`$erg (for $`ϵ=0.1`$) after it has traveled a distance such that it would have, on average, hit $`0.016`$, $`0.13`$ and $`1`$ other minihalos, respectively (the fraction of mass in minihalos of mass $`10^6\mathrm{M}_{}`$ $`f_M=10^3`$ is assumed). If minihalos are strongly clustered with a correlation length smaller than $`0.2d`$ (where $`d`$ is the mean separation equal to $`14(f/10^3)^{1/3}`$kpc proper for $`M_M=10^6\mathrm{M}_{}`$ at $`z=20`$), a chain reaction may be sustainable. Minihalos may very well be strongly clustered at high redshift (Barkana & Loeb 2003). We will study this separately.

In any event, blastwaves from stray stars propagating through the intergalactic medium would create dense shells. Figure 4 illustrate the temporal evolution of shocked gas at $`z=20`$ (assuming isobaric evolution); preshock gas is at the mean gas density at $`z=20`$. We see that cooling time is longer in this case than shown in Figure 2; here Compton cooling becomes important and its cooling time scale (i.e, e-folding time) is $`6\times 10^6`$yr. On the longer time scales shown in Figure 4 the detailed evolution would be affected by concomitant dynamic evolution driven by gravitational instability of large-scale density perturbations. Further investigation is needed to see if star formation in shells occurs. It is conceivable that a chain reaction may occur, which would result in a picture analogous to the explosion scenario proposed for larger galaxies at lower redshift (Ostriker & Cowie 1981).

## 4 Conclusions

In the cold dark matter model large, high-redshift ($`z>10`$) galaxies with virial temperature $`T_v>10^4`$K are shown to be mostly comprised of cold atomic clouds, which were originally minihalos and accreted into larger halos. It is shown that supersonic collisions of cold atomic gas clouds on time scale of $`10^7`$yr will shock heat and compress the gas, which then cools and form $`\mathrm{H}_2`$ molecules to further cool to $`100`$K. Massive star formation ensues. The resulting star formation is spatially distributed.

It is then shown that most of the star formation in minihalos may be triggered by shockwaves produced by the first stars formed in collisions. Those shockwave compressed minihalos will be more widespread spatially due to large imparted velocities and some can possibly escape into the intergalactic medium. If minihalos in the intergalactic space are strongly clustered, a chain reaction of star formation in minihalos may be triggered initiated by explosions of these stray stars. In this case, far more ionizing photons can be produced than large halos alone can. Separately, whether cooling shells of swept-up mean density gas in the intergalactic medium can form stars and possibly trigger an explosion-induced chain reaction even in the absence of strong minihalo clustering is unclear.

More importantly with respect to cosmological reionization, widespread star formation would allow an ionizing photon escape fraction much higher than in the conventional picture where star formation occurs in the central regions of galaxies and whose the escape fraction is estimated to be below $`1\%`$ for galaxies more massive than $`10^7\mathrm{M}_{}`$.

Such star formation may be the norm at $`z>10`$ and may create the physical basis to reconcile the standard cosmological model with the WMAP polarization observations. An optical depth $`\tau _e=0.100.14`$ can be achieved realistically even in the absence of chain reactions forementioned. A still higher optical depth may indicate that some further enhancement in star formation has occurred, such as through possible chain reactions. Detailed simulations of star formation triggered by collisions and blastwaves at high redshift galaxies would be required to firm up this picture and quantify key parameters, including star formation efficiency and ionizing photon escape fraction.

Finally, we note that these widespread star formation outlined above may give rise to metal enrichment of the intergalactic medium far more uniform than any other processes that have been proposed (e.g., Gnedin 1998; Aguirre et al. 2001; Cen et al. 2003).

This research is supported in part by grants AST-0206299 and NAG5-13381.


text:
# Interaction of moving discrete breathers with vacancies

## 1 Introduction

Discrete breathers (DBs) are classical, spatially localized, time-periodic, numerically exact solutions which can be sustained by many non-linear lattices . Their existence, which is proven by rigorous theorems  and a large amount of numerical results, is not restricted to special or integrable models. On the contrary, they can be found, in principle, in any discrete, nonlinear system and in any dimension . They have been observed in experiments involving different systems, as Josephson-junctions arrays , waveguide arrays , molecular crystals  and antiferromagnetic systems . They are also thought to play an important role in DNA denaturation .

Usually they are pinned to the lattice but under certain circumstances may become highly mobile . In this case, an interesting problem arises: the interaction between moving discrete breathers (MBs) and local inhomogeneities. This problem has been addressed within different frameworks: scattering of breathers or solitons by impurities , by lattice junctions  and by bending points  of a given chain.

Vacancies in a crystal are another type of local inhomogeneities. By using numerical methods, some of the authors found that a localized energy packet in the form of a moving discrete breather can put a vacancy into movement . This vacancy migration induced by localized excitations had already been suggested as an explanation for some experimental results .

The aim of this paper is to find an explanation to the somehow qualitative results established in . The main result is that vacancies mobility is highly dependent on the existence and stability of the breathers adjacent to the vacancy, hereafter called vacancy breathers (VBs).

## 2 The model

The simplest way to model a vacancy in a crystal consists in considering a one-dimensional chain of interacting particles submitted to a periodic substrate potential. It is known in other context as the Frenkel-Kontorova (FK) model . A vacancy can be represented by an empty well of the substrate potential as Fig. 1 shows.

The Hamiltonian of this system is

$$H=\underset{n=1}{\overset{N}{}}\frac{1}{2}m\dot{x}_n^2+V(x_n)+W(x_nx_{n1}),$$
(1)
being $`x_n`$ the absolute coordinate of the $`n`$-th particle. We have chosen the sine-Gordon potential

$$V(x_n)=\frac{a^2}{4\pi ^2}[1\mathrm{cos}(2\pi x_n/a)]$$
(2)
as the simplest periodic, substrate potential, with linear frequency normalized to the unity $`\omega _0=\sqrt{V^{\prime \prime }(0)}=1`$. We have chosen the Morse potential for the interaction between particles because its force weakens as the distance between particles grows and its hard part prevents the particles from crossing. It is given by

$$W(x_nx_{n1})=\frac{C}{2b^2}[e^{b(x_nx_{n1}a)}1]^2.$$
(3)
The parameter $`C=W^{\prime \prime }(a)`$ is the curvature of the Morse potential. We have taken $`C=0.5`$ in order that MBs exist in this system for a breather frequency $`\omega _b=0.9`$. The strength of the interaction potential can be modulated, without changing its curvature, by varying $`b`$, being $`b^1`$ a measure of the well width and $`C/2b^2`$ the well depth. We have also normalized the lattice period $`a`$ and the masses to the unity. The dynamical equations are given by

$$\ddot{x}_n+V^{}(x_n)+[W^{}(x_nx_{n1})W^{}(x_{n+1}x_n)]=0.$$
(4)
In this system, it is possible to generate DBs numerically using the standard methods from the anticontinuous limit . We can also induce translational motion of DBs by using a simplified form of the marginal mode method . It consists of adding a perturbation $`\stackrel{}{v}=\lambda (\mathrm{},0,1/\sqrt{2},0,1/\sqrt{2},0,\mathrm{})`$ to the velocities of the stationary breather, with the nonzero values at the neighboring sites of the initial breather center. The resulting DB kinetics is very smooth and resembles that of a classical free particle with constant velocity. Therefore, we can consider the total energy of a MB as the sum of a the internal energy, equal to the one of the stationary breather, plus the translational energy, equal to the energy of the perturbation added $`K=\lambda ^2/2`$.

In order to facilitate systematic studies, we have introduced dissipation at the boundaries.

## 3 Interaction of moving breathers with a single vacancy

In order to investigate vacancy mobility, we have generated DBs far from a vacancy and then launched this breather against it. Our numerical calculations show that the outcome of the scattering is extremely sensitive to the initial conditions . The incident breather can be reflected, trapped or transmitted (see Fig. 2), always losing energy as it occurs in the interaction between a MB and an impurity , whereas the vacancy can either move forward or backward or remain at rest. This scenario is very different to the one arising in the continuum limit where the vacancy (anti-kink) always moves backwards and the breather is always transmitted after the collision.

We have studied the dependence of the vacancy mobility on the parameter $`b`$, which controls the strength of the interaction between neighboring particles of the chain. For each value of $`b`$ we have analyzed a set of 601 numerical experiments of scattering corresponding to different initial conditions. As initial conditions we have considered a set of MBs with increasing kinetic energy, which were obtained choosing the values of the parameter $`\lambda `$ uniformly distributed in the interval $`[0.10,0.16]`$. The lower bound of this distribution is chosen so that the threshold value for vacancies movement (see below) is surpassed. The upper bound can be increased without any dramatic change in the results. However, an extremely high value of $`\lambda `$ can destroy the moving breather.

Fig. 3 shows the probability that the vacancy moves forward, backward or remains at rest (averaged over the set of initial conditions) with respect to the parameter $`b`$. Note that, as it should be expected, the backward movement is the most probable behavior because in this case the incident breather pushes forward the particle to the left of the vacancy. In fact, the forward movement of the vacancy requires a strong enough interaction, i.e. $`b`$ has to be smaller that a critical value $`b_f0.7`$.

Another interesting finding is shown in Fig. 4: the translational energy of the MB has to be higher than a minimal value $`K_{min}`$ in a certain interval $`b(0.5,0.8)`$ in order to move the vacancy. It is worth remarking that $`K_{min}=0`$ means the minimum kinetic energy needed to move a breather, which is actually different to zero. However, we have included this notation in order to have a clear picture of the vacancy movement.

## 4 Vacancy modes and vacancy breathers

We can try to understand these results by studying: a) the linear modes, which can be obtained by linearizing Eqs. 4 around the equilibrium state shown in Fig. 1; and b) the vacancy breathers, i.e., the nonlinear localized modes centered at the particles adjacent to the vacancy, with the same frequency as the incident MB.

Fig. 5a shows the dependence of the linear frequencies with respect to $`b`$. Note that the minimum linear frequency corresponds to a linear localized mode which we call linear vacancy mode (LVM). On the other hand, there are three types of VB: 1–site vacancy breathers (1VB), which consists of a single excited particle adjacent to the vacancy; and 2–site vacancy breathers, obtained by exciting the two particles adjacent to the vacancy, which can vibrate in-phase (2VBp) or anti-phase (2VBa). Fig. 5b shows the bifurcation pattern for these VBs. A letter at the beginning of the abbreviation specifies the breather stability: S stable, hU harmonically unstable, sU subharmonically unstable or oU oscillatory unstable. For example sU2VBa means a subharmonically–unstable, two–site, vacancy breather in anti–phase.

For weak interaction ($`b`$ large) there exist stable one-site VBs at both sides of the vacancy as stated by the theory . As $`b`$ decreases the S1VB disappears through an inverse pitchfork bifurcation to two unstable two-site VB in phase (hU2VBp) which becomes stable (S2VBp). Further decrease of $`b`$ brings about the annihilation of the S2VBp when the breather frequency coincides with the frequency of the LVM with the same profile.

The two–site VB in antiphase is stable (S2VBa) for the largest values of $`b`$ represented in Fig. 5b. When the interaction increases ($`b`$ decreases) this 2VB becomes oscillatory unstable (oU2VBa), and finally subharmonically unstable (sU2VBa).

From this bifurcation diagram we can deduce that the stability of the VBs is correlated to the existence of a kinetic energy threshold to move the vacancy, in the interval $`b(0.5,0.8)`$. Comparison of Fig. 4 and 5b suggests that the existence of a subharmonically or harmonically unstable, 2-site, vacancy breather is a necessary condition in order to get an optimum vacancy mobility (no kinetic energy threshold).

Generally speaking, any approaching breather with any kinetic energy is able to move the vacancy if at least one the two 2VBs exists and is either subharmonically or harmonically unstable. We cannot expect, however, an exact agreement between the numerically exact bifurcation values for a given frequency with the observed changes in behavior for the simulations. The reason for that is that a moving breather is obtained by perturbation with an asymmetric mode. Therefore, it is no longer an exact solution of the dynamical equations, its frequency is shifted and it is not unique. The bifurcation values for other frequencies will change and interaction of phonons emitted by the MB are expected to have an (unknown) influence, hence the probability analysis of the simulations outcome. However, the only important exception to the first assessment in this paragraph in the region $`b(0.73,0.8)`$ is noteworthy. The role of the unstable vacancy breather is probably to act as an intermediate structure with high oscillations for the particles nearest to the vacancy, which, being unstable, leads to one of them changing place. Clearly the stable S2VB does not play that role and its existence in the vicinity of that region seems to be the reason for that discrepancy –proximity in $`b`$ for a $`\omega _\mathrm{b}=0.9`$ also means proximity in $`\omega _\mathrm{b}`$ for the a given value of $`b`$, which can be excited by the bundle of the MB frequencies–.

Therefore our conclusion is that vacancy mobility is governed, for a particular choice of the substrate potential, by the strength of the interaction between the particles adjacent to the vacancy, which is correlated to the existence and stability properties of vacancy breathers. This correlation manifests in two different ways. Firstly, the interaction has to be strong enough to move the vacancy forward, or equivalently, there must exist no linear vacancy modes. Secondly, a threshold value of the translational energy of the incident MB in order to move the vacancy does not exist if an harmonically or subharmonically unstable VB exists.

Another choice of the interaction potential would alter the values of $`b`$ which separate the different regimes. However, the correlation between the vacancy mobility and the strength of the potential would be similar.

## 5 Double vacancy

We have studied the mobility of a double vacancy (two empty neighboring wells of the substrate potential) to test our conjecture. The outcome is similar to the single vacancy case except for the fact that the left vacancy never moves forwards. Thus, if the right vacancy moves forwards, the double vacancy splits in two (see Fig. 6). Fig. 7 shows the linear spectrum and the bifurcation diagram for the VB in this case. In these figures the lowest value of b is approximately 0.5 since there are no equilibrium states for $`b0.5`$. This result also implies the non existence of double vacancies equilibrium states for a harmonic interaction potential. As b increases, we observe two–site VBs emerging from linear localized modes, and afterwards an inverse pitchfork bifurcation similar to the one described for a single vacancy. There is not any harmonically or subharmonically unstable VB for $`b<0.75`$. Therefore, we should expect that there is a minimum value for the kinetic energy of the incident MBs able to move the double vacancy. This is confirmed by numerical computations as Fig.8 shows.

## 6 Summary

We have observed numerically the interaction of moving breathers with vacancies in the simplest, physically consistent model. We have found that the breathers can be reflected, transmitted and trapped, both by the single and the double vacancy. The single vacancy can move forward, backwards and remain and rest, but the double vacancy cannot move forward, instead it can split into two single vacancies, one at rest and the other moving forwards. This phenomenology is very different from the continuous analogue, where the breather is always transmitted and the vacancy (anti-kink) always moves backwards.

The vacancy migration is strongly correlated with the existence and stability of vacancy breathers, which we have studied in detail. In order to move the vacancy the interaction has to be strong enough and the energy of the incident MB has to be above a threshold. This threshold disappears approximately around bifurcation points where vacancy breathers become unstable.

We think that the study of the properties of breathers next to a defect may help to understand the mobility of other class of point defects such as interstitials, which might be related with recent defect migration observed experimentally in ion-irradiated silicon.

## 7 Acknowledgments

The authors acknowledges sponsorship by the Ministerio de Educación y Ciencia, Spain, project FIS2004-01183. They also acknowledge JC Eilbeck and F Palmero for useful discussions.