text:
# Fixing the conformal window in QCD11footnote 1Extended version of a talk given at the QCDNET2000 network meeting, Paris, September 11-14 2000 (http://qcd.th.u-psud.fr/QCDNET/).

## 1 Introduction

The notion of an infrared (IR) finite coupling has been used extensively in recent years, especially in connection  with the phenomenology of power corrections in QCD. The present investigation is motivated by the desire to understand better the theoretical background behind such an assumption. In particular, given an IR finite coupling $`\alpha `$, does it remain finite within perturbation theory itself (such as the two-loop coupling with opposite signs one and two loop beta function coefficients), or does one need a non-perturbative contribution $`\delta \alpha `$ to cancell ($`\alpha =\alpha _{PT}+\delta \alpha `$) the Landau singularities present in its perturbative part $`\alpha _{PT}`$? The answer I shall suggest is a mixed one: the perturbative part of the QCD coupling is always IR finite but, below the so called “conformal window” (the range of $`N_f`$ values where the theory is scale invariant at large distances and flows to a non-trivial IR fixed point), one still needs a $`\delta \alpha `$ term since the perturbative coupling is no more causal there, despite being IR finite. As the main outcome, one obtains an equation to determine the lower boundary of the conformal window in QCD. The plan of the paper is as follows. In section 2 I review the evidence and present a formal argument for the existence of Landau singularities in the perturbative coupling. A more physical argument, relating Landau singularities to the very existence of the conformal window and a two-phase structure of QCD is given in section 3, which also suggests a clean way to disentangle “perturbative” from “non-perturbative” below the conformal window. In section 4, two scenarios for causality breaking are described. In section 5, an equation to determine the bottom of the conformal window in QCD is suggested, and is solved through the Banks-Zaks expansion in section 6. Section 7 gives evidence, through a modified Banks-Zaks expansion, for the existence of a negative ultraviolet (UV) fixed point in QCD, necessary for the consistency of the present approach. Section 8 contains the conclusions.

## 2 Evidence for Landau singularities in the perturbative coupling

The only present evidence for a Landau singularity in the perturbative renormalized222In QED, the well established “triviality” property gives only direct evidence  for a singularity in the bare coupling constant. coupling is still the old Landau-Pomeranchuk leading log QED calculation, now reformulated in QCD as a $`N_f\mathrm{}`$ (“large $`\beta _0`$”) limit. In this limit, the perturbative coupling is one-loop: $`\alpha _{PT}(k^2)=1/\beta _0\mathrm{log}(k^2/\mathrm{\Lambda }^2)`$, where $`\mathrm{\Lambda }`$ is the Landau pole. The question is whether there is a singularity at finite $`N_f`$. Some light on this problem can be shed by considering further the $`N_f`$ dependence. Indeed, another (conflicting) piece of information is available at the other end of the spectrum, around the value $`N_f=N_f^0=16.5`$ (I consider $`N_c=3`$) where the one loop coefficient $`\beta _0=\frac{1}{4}(11\frac{2}{3}N_f)`$ of the beta function vanishes (“small $`\beta _0`$” limit). For $`N_f`$ slightly below $`16.5`$ a weak coupling (Banks-Zaks) IR fixed point develops  and the perturbative coupling is causal, i.e. there are no Landau singularities in the whole first sheet of the complex $`k^2`$ plane. Can then the perturbative coupling remain causal down to $`N_f=\mathrm{}`$? I shall assume that a Landau singularity cannot arise “spontaneously” in the limit, i.e. that “the limit of a sequence of causal couplings must itself be causal”. In such a case, the existence of a Landau pole at $`N_f\mathrm{}`$ implies the existence of a finite value $`N_f^{}`$ below which Landau singularities appear on the first sheet of the complex $`k^2`$ plane and perturbative causality is lost, which is the common wisdom (at $`N_f^{}`$ itself, according to the above philosophy, the coupling must still be causal). The range $`N_f^{}<N_f<N_f^0`$ where the perturbative coupling is causal and flows to a finite IR fixed point is taken as the definition of the “conformal window” for the sake of the present discussion. Let us now refine this definition, and give a more physical argument for the existence of Landau singularities, which illuminates their physical meaning.

## 3 Landau singularities and conformal window

I assume the existence of a two-phase structure in QCD as the number of flavors $`N_f`$ is varied.

i) For $`N_f^{}<N_f<N_f^0`$ (the conformal window) the theory is scale invariant at large distances, and the vacuum is “perturbative”, in the sense there is no confinement nor chiral symmetry breaking. Conformal window amplitudes (generically noted as $`D_{\overline{PT}}(Q^2)`$, where $`Q`$ stands for an external scale) are in this generalized sense “perturbative”, i.e. could in principle be determined from information contained in perturbation theory to all orders (although they should also include contributions from all instanton sectors): this motivates the subscript $`\overline{PT}`$.

ii) For $`0<N_f<N_f^{}`$ there is a phase transition to a non-trivial vacuum, with confinement and chiral symmetry breaking, as expected in standard QCD.

A direct, physical motivation for Landau singularities can now be given: they trace the lower boundary $`N_f=N_f^{}`$ of the conformal window. This statement is implied from the following two postulates:

1) Conformal window amplitudes $`D_{\overline{PT}}(Q^2)`$ can be analytically continued in $`N_f`$ below the bottom $`N_f^{}`$ of the conformal window.

2) For $`N_f<N_f^{}`$, the (analytically continued) conformal window amplitudes $`D_{\overline{PT}}(Q^2)`$ must differ from the full QCD amplitude $`D(Q^2)`$, since one enters a new phase, i.e. we have

$$D(Q^2)=D_{\overline{PT}}(Q^2)+D_{\overline{NP}}(Q^2)$$
(1)
(whereas $`D(Q^2)D_{\overline{PT}}(Q^2)`$ within the conformal window). Assuming QCD to be a unique theory at given $`N_f`$, $`D_{\overline{PT}}(Q^2)`$ cannot provide a consistent solution if $`N_f<N_f^{}`$: this must be signalled by the appearance of unphysical Landau singularities in $`D_{\overline{PT}}(Q^2)`$. $`N_f^{}`$ should thus coincide with the value of $`N_f`$ below which (first sheet) Landau singularities appear in $`D_{\overline{PT}}(Q^2)`$. The occurence of a “genuine” non-perturbative component $`D_{\overline{NP}}(Q^2)`$ is then necessary below $`N_f^{}`$ in order to cancell the Landau singularities present in $`D_{\overline{PT}}(Q^2)`$. If these assumptions are correct, they provide an interesting connection between information contained in principle in “perturbation theory” (over all instanton sectors), which fix the structure of the conformal window amplitudes and “genuine” non-perturbative phenomena, which fix the bottom of the conformal window. In addition, eq.(1) provide a neat way to disentangle the “perturbative” from the genuine “non-perturbative” part of an amplitude, for instance the part of the gluon condensate related to renormalons from the one reflecting the presence of the non-trivial vacuum. Note also $`D_{\overline{PT}}(Q^2)`$ and $`D_{\overline{NP}}(Q^2)`$ are separately free of renormalons ambiguities, but contain Landau singularities below $`N_f^{}`$, so the renormalon and Landau singularity problems are also disentangled. In order to get a precise condition to determine $`N_f^{}`$, one needs to look in more details how causality can be broken in the perturbative coupling.

## 4 Scenarios for causality breaking

There are two main scenarios:

i) The “standard” one where the IR fixed point present within the conformal window just disappears when $`N_f<N_f^{}`$ while a real, space-like Landau singularity is generated in the perturbative coupling. For instance, two simple zeroes of the beta function can merge into a double zero when $`N_f=N_f^{}`$ before moving to the complex plane (a plausible scenario  in supersymmetric QCD (SQCD)).

ii) Alternatively, it is possible for the fixed point to be still present333This assumption is consistent with the suggestion  that the perturbative coupling has a non-trivial IR fixed point down to $`N_f=2`$ in QCD. However the full non-perturbative coupling must still differ by a $`\delta \alpha `$ term, since the perturbative coupling is non-causal below $`N_f=N_f^{}`$. in the perturbative part of the coupling in some range of $`N_f`$ below $`N_f^{}`$. The motivation behind this assumption is the observation  that for QCD effective charges associated to Euclidean correlators (the only ones for which the notion of $`k^2`$ plane analyticity makes sense) the Banks-Zaks expansion in QCD (as opposed to SQCD ) seems to converge down to fairly small values of $`N_f`$. In this case there can be no space-like Landau singularity, and causality must be violated by the appearance of complex Landau singularities on the first sheet of the $`k^2`$ plane. The example below suggests that they arise as the result of the continuous migration to the first sheet, through the time-like cut, of some second sheet singularities already present when $`N_f>N_f^{}`$. I shall assume that this is the scenario which prevails in QCD. As the simplest example, consider the two-loop coupling, which satisfies the renormalization group (RG) equation $`d\alpha /d\mathrm{log}k^2=\beta _0\alpha ^2\beta _1\alpha ^3`$. If $`\beta _0>0`$ but $`\beta _1<0`$, there is an IR fixed point at $`\alpha _{IR}=\beta _0/\beta _1`$. It has been shown  that this coupling has a pair of complex conjugate Landau singularities on the second (or higher) sheet if

$$0<\gamma _{2loop}=\frac{\beta _0^2}{\beta _1}<1$$
(2)
where $`\gamma _{2loop}`$ is the 2-loop critical exponent (see below). For $`\gamma _{2loop}>1`$, the second sheet singularities move to the first sheet through the time-like cut, which is reached when $`\gamma _{2loop}=1`$. The latter condition thus determines the bottom of the conformal window in this model. Note that in the limit $`\beta _10^{}`$ where $`\gamma _{2loop}=+\mathrm{}`$, one gets the one loop coupling and the complex conjugate singularities collapse to a space-like Landau pole. This limit is thus the analogue of the $`N_f\mathrm{}`$ limit in QCD.

## 5 An equation to determine the bottom of the conformal window in QCD

Assuming from now on that the second scenario described above applies, i.e. that the IR fixed point persists in some range of $`N_f`$ below $`N_f^{}`$, one needs some information on the location of Landau singularities in coupling constant space to derive an equation for $`N_f^{}`$. I shall assume throughout that there are no complex (in $`\alpha `$ space) Landau singularities (such as complex poles in the beta function), i.e. that the Landau singularities originate only from the $`\alpha <0`$ or from the $`\alpha >\alpha _{IR}`$ regions, and argue that the condition

$$0<\gamma <1$$
(3)
is both necessary and sufficient for causality in QCD. Consequently, the lower boundary of the conformal window is obtained from the equation

$$\gamma (N_f=N_f^{})=1$$
(4)
where $`\gamma `$ is the critical exponent defined as the derivative of the beta function at the fixed point

$$\gamma =\frac{d\beta (\alpha )}{d\alpha }|_{\alpha =\alpha _{IR}}$$
(5)
As is well known, this is a universal quantity, independant of the definition of the coupling, and eq.(4) is a renormalization scheme invariant condition, as it should. The argument proceeds in two steps.

1) Let us first assume (this will be justified below) there is an $`\alpha >\alpha _{IR}`$ UV Landau singularity (for instance a pole in the beta function at $`\alpha _P>\alpha _{IR}`$), in the domain of attraction of $`\alpha _{IR}`$. One can then show  that eq.(3) is a necessary condition for causality. I give an improved version of the argument of . Solving the RG equation $`d\alpha /d\mathrm{log}k^2=\beta (\alpha )`$ around $`\alpha =\alpha _{IR}`$, one gets

$$\alpha (k^2)=\alpha _{IR}\left(\frac{k^2}{\mathrm{\Lambda }^2}\right)^\gamma +\mathrm{}$$
(6)
There are thus rays

$$k^2=|k^2|\mathrm{exp}\left(\pm \frac{i\pi }{\gamma }\right)$$
(7)
in the complex $`k^2`$ plane, which in the infrared limit $`|k^2|0`$ are mapped by eq.(6) to positive real values of the coupling larger then $`\alpha _{IR}`$. Assuming an expansion

$$\beta (\alpha )=\gamma (\alpha \alpha _{IR})+\gamma _1(\alpha \alpha _{IR})^2+\mathrm{}$$
(8)
the corrections to eq.(6) are given by a series

$$\mathrm{log}(k^2/\mathrm{\Lambda }^2)=\frac{1}{\gamma }\mathrm{log}(\alpha _{IR}\alpha )+\frac{\gamma _1}{\gamma ^2}(\alpha _{IR}\alpha )+\mathrm{}$$
(9)
with real coefficients, showing that the only contribution to the phase for $`\alpha >\alpha _{IR}`$ comes from the logarithm on the right hand side of eq.(9). The trajectories in the $`k^2`$ plane which map to the $`\alpha >\alpha _{IR}`$ region are thus straight lines to all orders of perturbation theory around $`\alpha _{IR}`$. This fact suggests that even away from the infrared limit, these trajectories are given by the rays eq.(7). As $`|k^2|`$ is increased along these rays, the coupling will flow to the assumed UV Landau singularities , reached at some finite value of $`|k^2|`$. If $`\gamma >1`$ the rays, hence also the singularities, are located on the first sheet of the $`k^2`$ plane, showing that eq.(3) is a necessary condition for causality.

2) To assert whether eq.(3) is also sufficient for causality, one has to make sure that no other sources of (first sheet) Landau singularities are present, but the one arising from the $`\alpha >\alpha _{IR}`$ region. Barring Landau singularities at complex $`\alpha `$, a potential problem can still arise from an eventual UV Landau singularity at $`\alpha <0`$, in the domain of attraction of the trivial IR fixed point $`\alpha =0^{}`$. Indeed at weak coupling the solution of the RG equation is controlled by the 2-loop beta function

$$\mathrm{log}(k^2/\mathrm{\Lambda }^2)=\frac{1}{\beta _0\alpha }+\frac{\beta _1}{\beta _0^2}\mathrm{log}\alpha +const+\mathrm{}.$$
(10)
where the $`const`$ is real. For $`\alpha <0`$ the right hand side of eq.(10) acquires a $`\pm i\pi \frac{\beta _1}{\beta _0^2}`$ imaginary part, which implies the rays

$$k^2=|k^2|\mathrm{exp}\left(\pm i\pi \frac{\beta _1}{\beta _0^2}\right)$$
(11)
map to the $`\alpha <0`$ region. Along the rays eq.(11), we are effectively in a QED like situation: increasing $`|k^2|`$, the coupling is either attracted to a non-trivial UV fixed point, or reaches an UV Landau singularity at some finite $`|k^2|`$. In the latter case, one must require that the condition:

$$|\frac{\beta _0^2}{\beta _1}|<1$$
(12)
is satisfied in the whole $`N_f`$ range where eq.(3) is valid, which will confine the rays, hence the singularities to the second (or higher) sheet.

However in QCD condition eq.(12) can be satisfied only if $`\beta _1<0`$, and coincides444Eq.(12) is however a necessary condition for causality for any beta function which admits an UV Landau singularity at negative $`\alpha `$ (in the domain of attraction of the $`0^{}`$ trivial IR fixed point), and applies also if $`\beta _1`$ is positive in a general theory! with the 2-loop causality condition eq.(2), which requires $`N_f>9.7`$. Therefore, to preserve causality within the conformal window as determined555It is a priori possible in QCD to have an $`\alpha <0`$ UV Landau singularity rather then an $`\alpha <0`$ UV fixed point. In such a case the bottom of the conformal window would be given as in the two-loop model by the condition $`\beta _0^2/\beta _1=1`$, yielding $`N_f^{}9.7`$. Indeed, at such large $`N_f`$, any eventual UV Landau singularity from the $`\alpha >\alpha _{IR}`$ region has not yet reached the first sheet since $`\gamma <0.4`$ as shown by Fig. 1. This possibility is however disfavored as explained in the text. by eq.(3), $`\gamma `$ should reach $`1`$ in the region $`N_f>9.7`$, which is clearly excluded (see Fig. 1 below). In order that eq.(3) be also a sufficient condition for causality one must thus check that a non-trivial (finite or infinite) UV fixed point $`\alpha _{UV}`$ is present at negative $`\alpha `$! A minimal example satisfying this requirement is the 3-loop beta function $`\beta (\alpha )=\beta _0\alpha ^2\beta _1\alpha ^3\beta _2\alpha ^4`$ with $`\beta _0>0`$ and $`\beta _2<0`$ ($`\beta _1`$ can have any sign).

It is worth mentionning eq.(12) is always violated  in the lower part of the conformal window in SQCD as determined by duality , and the previous argument thus implies the existence of a negative UV fixed point in this theory. In fact the “exact” NSVZ  beta function for $`N_f=0`$ does exhibit an (infinite) UV fixed point as $`\alpha \mathrm{}`$, which might be the parent of a similar one present within Seiberg conformal window.

There is some evidence that a negative UV fixed point is indeed also present in QCD. At the three-loop level, QCD effective charges associated to Euclidean correlators have typically $`\beta _2<0`$, and appear  to be causal and admit a negative UV fixed point, even somewhat below the two-loop causality boundary $`N_f=9.7`$. This evidence could be washed out in yet higher orders (e.g. the Pade improved three-loop beta functions ). However, a different and more sytematic derivation is given in section 7. The plausible existence of a nearby negative UV fixed point666The non-trivial UV fixed point is actually not relevant to the proper analytic continuation of the coupling at complex $`k^2`$, which must be consistent  with (UV) asymptotic freedom. This means that in presence of this fixed point, the correct analytic continuation must involve complex rather then negative values of $`\alpha `$ along the rays eq.(11), and one should approach the non-trivial (rather then the trivial) IR fixed point as $`|k^2|0`$, and the trivial (rather then the non-trivial) UV fixed point as $`|k^2|\mathrm{}`$. This is possible since the solution of eq.(10) is not unique for a given (complex) $`k^2`$. For the same reason, any eventual IR Landau singularity arising from the region $`\alpha <\alpha _{UV}`$, in the domain of attraction of the non-trivial UV fixed point, is not relevant to the correct analytic continuation. On the other hand, any UV Landau singularity in the domain of attraction of either the trivial or the non-trivial IR fixed points as considered above is relevant to the proper analytic continuation, since the coupling will flow to the only UV fixed point available, namely the trivial one, once the UV Landau singularity is passed.  in QCD, hence the absence of an UV Landau singularity at negative $`\alpha `$, justifies a posteriori (barring complex $`\alpha `$ singularities) the above assumption that there must be an $`\alpha >\alpha _{IR}`$ singularity, to provide the necessary causality violation below the conformal window.

## 6 Computing $`N_f^{}`$ through the Banks-Zaks expansion

One can try to use the Banks-Zaks expansion to compute $`\gamma `$ and determine $`N_f^{}`$, the lower boundary of the conformal window. The Banks-Zaks expansion  is an expansion of the fixed point in powers of the distance $`N_f^0N_f=16.5N_f`$ from the top of the conformal window, which is proportional to $`\beta _0`$. The solution of the equation

$$\beta (\alpha )=\beta _0\alpha ^2\beta _1\alpha ^3\beta _2\alpha ^4+\mathrm{}=0$$
(13)
in the limit $`\beta _00`$, with $`\beta _i`$ ($`i1`$) finite is obtained as a power series

$$\alpha _{IR}=a+𝒪(a^2)$$
(14)
where the expansion parameter $`a\frac{8}{321}(16.5N_f)=\frac{16}{107}\beta _0`$. The Banks-Zaks expansion for the critical exponent eq.(5) is presently known  up to next-to-next-to leading order :

$$\gamma =\frac{107}{16}a^2(1+4.75a8.89a^2+\mathrm{})$$
(15)
Using the trunkated expansion eq.(15), one finds that $`\gamma <1`$ for $`N_f5`$, with $`\gamma =1`$ reached for $`N_f=N_f^{}4`$. To assess whether it is reasonable to use perturbation theory down to $`N_f=N_f^{}`$, let us look at the magnitude of the successive terms within the parenthesis in eq.(15). They are given by: $`1,1.44,0.82`$. Although the next to leading term gives a very large correction, and the series seem at best poorly converging at $`N_f=N_f^{}`$, one can observe that the next-to next to leading term still gives a moderate correction to the sum of the first two terms, which might be considered together  as building the “leading” contribution, since they are both derived from information contained  in the minimal 2-loop beta function necessary to get a non-trivial fixed point. Indeed, keeping only the first two terms in eq.(15), one finds that $`\gamma =1`$ is reached for $`N_f=N_f^{}6`$. On the other hand, using a $`[1,1]`$ Pade approximant as a model777The alternative $`[0,2]`$ Pade yields a result ($`\gamma <0.26`$ for any $`N_f`$) inconsistent with the present framework. It also predicts a not very plausible $`𝒪(a^5)`$ coefficient of $`192`$. for extrapolation of the perturbative series, one gets

$$\gamma =\frac{107}{16}a^2\frac{1+6.62a}{1+1.87a}$$
(16)
which yields $`\gamma =1`$ for $`N_f=N_f^{}5`$. Fig. 1 shows $`\gamma `$ as a function of $`N_f`$:

Note that in the obtained range of $`N_f^{}`$ values ($`4<N_f^{}<6`$), $`\beta _1`$ is still positive ($`\beta _1`$ changes sign for $`N_f=8.05`$) and of the same sign as $`\beta _0`$. The fixed point must thus arise from the contributions of higher then 2 loop beta function corrections, although I am assuming the Banks-Zaks expansion is still converging there. This is consistent with the previously mentionned fact that QCD effective charges have negative 3-loop beta function coefficients in the above range of $`N_f`$ values.

## 7 Further evidence for a negative UV fixed point in QCD

Additional evidence for the existence of a couple of (negative-positive) UV-IR fixed points is provided by the following modified Banks-Zaks argument. Assume $`\beta _1=0`$, i.e. $`N_f=8.05`$. Then a real fixed point can still exist at the three-loop level if $`\beta _2<0`$, and actually one gets a pair of opposite signs zeroes at $`\stackrel{~}{\alpha }=\pm (\beta _0/\beta _2)^{1/2}`$, an IR and an UV fixed point. If $`\beta _0`$ is small enough, they are weakly coupled, and calculable through a modified Banks-Zaks expansion around $`N_f=16.5`$, applied to the auxiliary function $`\stackrel{~}{\beta }(\alpha )\beta (\alpha )+\beta _1\alpha ^3`$ with the two-loop term removed. One gets

$$\stackrel{~}{\alpha }=\pm \stackrel{~}{a}\left(1\frac{1}{2}\frac{\beta _{3,0}}{\beta _{2,0}}\stackrel{~}{a}+\mathrm{}\right)$$
(17)
where the expansion parameter $`\stackrel{~}{a}(\beta _0/\beta _{2,0})^{1/2}`$, and $`\beta _{i,0}`$, the $`N_f=16.5`$ values of $`\beta _i`$ ($`i=2,3`$), are scheme dependent and can be obtained  from the knowledge of the coefficients in eq.(15). For effective charges associated to Euclidean correlators, the correction in eq.(17) ranges from $`0.1`$ to $`0.7`$ at $`N_f=8.05`$ (where $`\stackrel{~}{\beta }(\alpha )`$ coincides with $`\beta (\alpha )`$), which is encouraging evidence for a couple of UV and IR fixed points around this value of $`N_f`$ (which is within the alleged conformal window, but below the 2-loop causality region). Additional support is given by the calculation of the auxiliary critical exponent $`\stackrel{~}{\gamma }d\stackrel{~}{\beta }(\alpha )/d\alpha |_{\alpha =\stackrel{~}{\alpha }_{IR}}`$, whose modified Banks-Zaks expansion is

$$\stackrel{~}{\gamma }=2\beta _0\stackrel{~}{a}(1+𝒪(\stackrel{~}{a}^2))$$
(18)
(no $`𝒪(\stackrel{~}{a})`$ correction!), which yields at $`N_f=8.05`$ (where it coincides with $`\gamma `$) $`0.6<\stackrel{~}{\gamma }<0.7`$ for effective charges associated to Euclidean correlators, in reasonable agreement with the standard Banks-Zaks result (Fig. 1) $`0.5<\gamma <0.6`$.

## 8 Conclusions

1) I have suggested a direct, physical motivation for Landau singularities, assuming a two-phase structure of QCD: they should trace the lower boundary $`N_f^{}`$ of the conformal window. This approach avoids the notoriously tricky disentangling of the “perturbative” from the “non-perturbative” part of the QCD amplitudes within the conformal window, since they are by definition entirely “perturbative” there. On the other hand, such a disentangling is naturally achieved below the conformal window, by introducing the analytic continuation of the conformal window amplitudes to the $`N_f<N_f^{}`$ region.

2) Assuming that the perturbative QCD coupling has a non-trivial IR fixed point $`\alpha _{IR}`$ even in some range below $`N_f^{}`$ leads to the equation $`\gamma (N_f=N_f^{})=1`$ to determine $`N_f^{}`$ from the critical exponent $`\gamma `$ at the IR fixed point. Using the available terms in the Banks-Zaks expansion, this equation yields $`4N_f^{}6`$. It would clearly be desirable to have more terms to better control the accuracy of the Banks-Zaks expansion. Note that this condition is inconsistent with the existence of an analogue of Seiberg duality in QCD, which would rather imply $`\gamma (N_f=N_f^{})=0`$.

3) Some conditions on the QCD beta function are required: the only source of (UV) Landau singularities must come from the $`\alpha >\alpha _{IR}`$ region. One needs in particular a negative UV fixed point $`\alpha _{UV}`$. There is indeed some evidence for such a fixed point in QCD.

4) A negative UV fixed point is also required in the SQCD case, where duality fixes the conformal window.

5) It is possible the finite IR fixed point persists in the perturbative QCD coupling down to the $`N_f\mathrm{}`$ one-loop limit. A simple example is provided by a beta function with one positive pole $`\alpha _P`$ (the required Landau singularity) and two opposite sign zeroes $`\alpha _{IR}`$ and $`\alpha _{UV}`$: $`\beta (\alpha )=\beta _0\alpha ^2(1\alpha /\alpha _{IR})(1\alpha /\alpha _{UV})/(1\alpha /\alpha _P)`$ where $`\alpha _{UV}<0`$ and $`0<\alpha _{IR}<\alpha _P`$. The one-loop limit is achieved for $`\alpha _{IR}=\alpha _P`$ and $`\alpha _{UV}=\mathrm{}`$. A paper developing further these issues is under preparation.

###### Acknowledgments.

I thank A. Armoni and A.H. Mueller for discussions.


format: text

text:
# Intermittent Switching between Potential Flow and Turbulence in Superfluid Helium at mK Temperatures

## Abstract

Superfluid flow around an oscillating microsphere is investigated at temperatures down to 25 mK. Stable laminar flow below a critical velocity and turbulence at large drives are found to be separated below 0.5 K by an intermediate range of driving forces where the flow is unstable, intermittently switching between laminar and turbulent phases. We have recorded time series of this switching phenomenon and have made a statistical analysis of the switching probability. The mean lifetime of the turbulent phases grows with increasing drive and becomes infinite at a critical value. Stability of the laminar phases above the critical velocity is limited by natural background radioactivity or cosmic rays.

Turbulence in a classical liquid is based on the nonlinearity of the Navier-Stokes equation of motion  while turbulence in superfluid helium, which also has been a very active area of research for several decades, is described by two-fluid hydrodynamics of the normal and superfluid components of the liquid . Because of the quantum mechanical properties of the liquid, turbulence in superfluid helium consists of a tangle of individual vortex lines having the same quantized circulation $`h/m_4`$, where $`h`$ is Planck’s constant and $`m_4`$ is the mass of a helium-4 atom. In this respect, superfluid turbulence is different from classical turbulence but it has also striking similarities like a Kolmogorov energy spectrum in the inertial range . Both cases represent different branches of the general phenomenon of fluid turbulence. A completely new aspect arises when superfluid turbulence is considered in the limit of very low temperatures . In this case all what is left of the normal fluid component is a very dilute gas of ballistically propagating phonons and therefore turbulence in the pure superfluid can be investigated. While a number of experiments have been performed which demonstrate the production of individual vortex lines both in superfluid 4He and in superfluid 3He  and the decay of the vortex tangle , the transition to turbulence in the flow of the superfluid around a macroscopic body so far has not been investigated at these very low temperatures.

In our present work, we have studied superfluid flow around an oscillating microsphere at temperatures down to 25 mK. Laminar and turbulent flow can be identified by the drag force on the sphere, which in case of laminar flow is linear and determined by ballistic phonon scattering, while for turbulent flow, which exists above a critical velocity, the drag force is nonlinear and much larger . Here, we report on the observation of an intermediate regime located between stable laminar flow below the critical velocity and stable turbulence above a critical driving force, where the flow is found to be unstable, switching intermittently between laminar and turbulent phases. We have measured time series of this switching phenomenon at different constant temperatures and driving forces. We find, firstly, the mean lifetime of the turbulent phases to grow with increasing drive and to diverge at a critical value. Secondly, the stability of the laminar phases above the critical velocity is analyzed and the probability for switching to turbulence is obtained. Surprisingly, we find the long term stability of the laminar phases close to the critical velocity to be limited ultimately by natural radioactivity or cosmic rays.

Our experimental technique has been described in earlier work  and is only briefly summarized here as follows. The microsphere (radius $`R=0.1\text{mm}`$, mass $`m=27\mu \text{g}`$) is made of ferromagnetic $`\mathrm{SmCo}_5`$. Therefore, we can use superconducting levitation, i.e., the repulsion of a permanent magnet from a superconducting surface, to let it float in the middle between the superconducting electrodes of a parallel plate capacitor made of niobium (spacing $`d=1\text{mm}`$, diameter $`4\text{mm}`$) without any mechanical suspension. Lateral stability is provided by flux lines trapped in the superconductor. The sphere carries an electric charge $`q`$ of about 1 pC when it levitates, because we apply a voltage of several hundred volts to the bottom electrode when cooling the capacitor through the critical temperature $`T_c=9.2\text{K}`$ of niobium. Vertical oscillations of the sphere about its equilibrium position can now be excited by applying a small ac voltage $`U_{ac}`$ ranging from $`0.1\text{mV}`$ up to $`20\text{V}`$ at the resonance frequency of the sphere $`120\text{Hz}`$. These oscillations induce an ac current $`qv/d`$ in the electrodes, where $`v`$ is the velocity amplitude, which is detected by an electrometer. The damping of the oscillations is very small, $`Q`$-values in excess of one million are achieved if the capacitor is evacuated . We then fill the capacitor with pure superfluid helium and measure the velocity amplitude at resonance as a function of the driving force $`F=qU_{ac}/d`$ at constant temperature between 25 mK and 1 K.

As reported earlier  there is a linear regime $`v(F)`$ at small amplitudes, see Fig. 1. This indicates a linear drag force $`\lambda (T)v`$ where the drag coefficient $`\lambda T^4`$ is determined by ballistic scattering of thermal phonons which can be considered as a dilute gas of excitations whose density rapidly decreases as $`T^4`$. The flow of the superfluid around the sphere is frictionless potential flow. At large driving forces, however, the sphere creates vorticity in the superfluid. Now the drag force is large and nonlinear, see Fig. 1. We find the drag force to be given by $`F_D=\gamma v^2F_0`$, where $`\gamma =c_D\rho \pi R^2/2`$ is the classical turbulent drag coefficient ($`\rho `$ is the density of the liquid and $`c_D0.4`$ for a sphere) and $`F_0`$ is a constant. In contrast to a classical liquid where $`F_0=0`$ in our case turbulent drag ($`F_D>0`$) exists only above a critical velocity $`v_c=19.4\text{mm/s}`$. A similar drag force has recently been calculated for a cylinder in a two-dimensional dilute Bose-Einstein condensate, although the critical velocity is very much different in that case .

Most interesting is our new observation of an intermediate regime where neither potential flow nor turbulence are stable but instead the flow switches between both patterns intermittently . This regime of instability is observed only below $`0.5\text{K}`$ and extends from the critical velocity $`v_c`$ up to a critical driving force where enough power is being delivered to the superfluid turbulence to become stable. We have recorded time series (up to 36 hours per run) of this switching phenomenon at various driving forces and temperatures. An example is shown in Fig. 2 where a section of 17 minutes is shown for three different driving forces at $`300\text{mK}`$. The amplitude switches between a low level corresponding to turbulent drag and an exponential recovery of the laminar level corresponding to phonon scattering. If the laminar phase lasts long enough the stationary value for laminar flow is reached (see upper time series). When the laminar flow breaks down a rapid drop to the turbulent level occurs. It is obvious from Fig. 2 that with increasing driving force the lifetimes of the laminar phases shrink while those of the turbulent phases grow.

In the following we discuss the results of a statistical analysis of these time series based on reliability theory . A typical distribution of the lifetimes of the turbulent phases is shown in Fig. 3a where the number $`P`$ of turbulent phases having a lifetime longer than a given threshold $`t`$ is plotted on a logarithmic scale versus $`t`$. We always find an exponential distribution $`P(t)\mathrm{exp}(\nu t)`$ with $`1/\nu `$ being the mean lifetime. This implies a probability for the breakdown of the turbulent phase (”failure rate”) $`d\mathrm{ln}P/dt=\nu `$ being constant in time (i.e no history dependence). The values of $`\nu `$ are independent of temperature and collapse to an universal drive dependence if the strongly temperature dependent laminar drag $`\lambda v`$ (where $`v=v_t`$ is the velocity amplitude of the turbulent phase) is subtracted from the external driving force. Obviously, $`\nu `$ depends only on the turbulent drag force $`F\lambda v_t`$ and vanishes at a critical drive approximately as the fourth power of the distance from a critical value, see Fig. 3b. At the critical drive the power delivered to the turbulent superfluid is $`0.6\text{pW}`$ which corresponds approximately to the production of one vortex ring of the size of the sphere per half period of the oscillation.

Stability of the laminar phases can be analyzed in two ways, either the distribution of the signal height $`\mathrm{\Delta }v`$ (velocity at breakdown minus turbulent velocity $`v_t`$) is considered or the distribution of the lifetimes. Both variables are related by the exponential recovery $`\mathrm{\Delta }v(t)=\mathrm{\Delta }v_{max}(1\mathrm{exp}(t/\tau ))`$, where $`\mathrm{\Delta }v_{max}`$ is the difference between the stationary laminar velocity amplitude $`F/\lambda `$ and $`v_t`$, and $`\tau =2m/\lambda `$ is the time constant. In Fig. 4a the distribution of the signal height $`\mathrm{\Delta }v`$ is shown to be of the Gaussian form $`P(\mathrm{\Delta }v)\mathrm{exp}((\mathrm{\Delta }v/v_W)^2)`$ where $`v_W`$ is a fitting parameter. This implies a failure rate (per unit velocity increment) $`d\mathrm{ln}P/d\mathrm{\Delta }v=2\mathrm{\Delta }v/v_W^2`$ to be proportional to the signal height $`\mathrm{\Delta }v`$. The corresponding probability density of $`\mathrm{\Delta }v`$ ($`dP/d\mathrm{\Delta }v`$) is called a Weibull distribution. The fitting parameter $`v_W4.8\text{mm/s}`$ is found to be completely independent of both temperature and driving force, see Fig. 4b. We expect it to depend on the size or the surface roughness of the sphere but experiments with different spheres remain to be performed in the future. Because $`\mathrm{\Delta }v(t)`$ saturates for $`t\tau `$, it is more informative in this case to analyze the lifetime distribution, in particular near $`v_c`$ where the laminar phase becomes stable. From $`P(t)=P(\mathrm{\Delta }v(t))\mathrm{exp}((\mathrm{\Delta }v(t)/v_W)^2)`$ we conclude that for $`t\tau `$ $`P(t)`$ must be constant and therefore the failure rate (per unit time) is $`d\mathrm{ln}P/dt=v_W^2d(\mathrm{\Delta }v)^2/dt=0`$. This simply means that if the laminar phase has survived for a time $`t\tau `$ it will be stable indefinitely although the velocity is clearly above $`v_c`$. Our experimental results, however, demonstrate that the lifetime is finite: $`\mathrm{ln}P(t)`$ slowly decreases at large $`t`$ and from a straight line fit to the data in Fig. 5 we obtain a mean lifetime of 25 minutes. We can rule out mechanical vibrations of the cryostat as the source of this instability: jumping on the floor, slamming the door, even transferring liquid helium to the cryostat had no effect. Therefore, another mechanism must be responsible. In fact, placing a small radioactive source $`{}_{}{}^{60}\mathrm{Co}`$ ($`74\text{kBq}`$) near the cryostat has a dramatic effect: the mean lifetime of the laminar phase is reduced to 3.0 minutes. We have measured the dose rate of the source at the position of the measuring cell inside the cryostat (taking into account a measured 20% loss in the dewar walls) to be 440 nGy/h ($`\pm `$ 5%). Comparing this value with a measured dose rate due to natural background radiation in our laboratory of 50 nGy/h ($`\pm `$ 10%), which is typical for our area, we obtain an increase of the dose rate due to the source by a factor of $`(440+50)/50=9.8`$. Within the error bars this compares well with the measured reduction of the mean lifetime of the laminar phases by a factor of $`25/3.0=8.3`$. This effect may be attributed to the creation of ions in the superfluid, which produce local vorticity either during the creation and recombination processes or when the ions are accelerated by the electric field which exists in our measuring cell . This vorticity then initiates the breakdown of the potential flow if the velocity amplitude of the sphere is larger than the critical velocity $`v_c`$. We therefore conclude that the lifetime of the metastable laminar phases which we observe above $`v_c`$ is only limited by natural background radioactivity.

In summary, we have found that the transition from potential flow to turbulence around a sphere in superfluid helium below 0.5 K occurs by intermittent switching between both flow patterns. We understand now why above 0.5 K the intermittent switching changes into the hysteretic behavior observed earlier : the velocity increases $`\mathrm{\Delta }v`$ become very small at higher phonon drag which implies a very low failure rate of the laminar phase. But if it fails (i.e. when the driving force $`F`$ is largely increased) the following turbulent phase is stable because the critical drive is exceeded.

At present we have no theoretical model which could describe our experimental results for the failure rates of the laminar and turbulent phases. But because the transition to turbulence in superfluid helium at mK temperatures is a new phenomenon, we hope that our results will stimulate further theoretical progress.

We are most grateful to H. von Philipsborn for measuring the dose rates and for illuminating discussions. We thank M. Creuzburg for clarifying comments on dosimetry and J. Reisinger for the loan of the radioactive source. We had helpful discussions on turbulence with H.R. Brand. This work is supported by the Deutsche Forschungsgemeinschaft.


text:
# Quantum Monte Carlo calculation of the finite temperature Mott-Hubbard transition

\[

## Abstract

We present clear numerical evidence for the coexistence of metallic and insulating dynamical mean field theory(DMFT) solutions in a half-filled single-band Hubbard model with bare semicircular density of states at finite temperatures. Quantum Monte Carlo(QMC) method is used to solve the DMFT equations. We discuss important technical aspects of the DMFT-QMC which need to be taken into account in order to obtain the reliable results near the coexistence region. Among them are the critical slowing down of the iterative solutions near phase boundaries, the convergence criteria for the DMFT iterations, the interpolation of the discretized Green’s function and the reduction of QMC statistical and systematic errors. Comparison of our results with those of other numerical methods is presented in a phase diagram.

\]

The Mott-Hubbard Transition  is an outstanding open problem in condensed matter physics. The theoretical progress in this field has been impeded by non-perturbative nature of the strongly correlated system.

In recent years considerable progress has been made with a single band Hubbard model by using DMFT approximation which becomes exact in the limit of infinite coordinations ($`z=\mathrm{}`$. The most striking point of these works is a discontinuous metal-insulator transition(MIT) at finite temperature $`T`$ which takes place at the critical interaction strength line $`U_c(T)`$. This line lies inside a coexistence region where two phases, one metallic and one insulating coexist. The coexistence region is bounded by two lines $`U_{c1}(T)`$ and $`U_{c2}(T)`$ and $`U_{c1}(T)<U_c(T)<U_{c2}(T)`$ .

While this picture is now generally accepted  by all groups working on the problem, the origin of the controversial results over the coexistence region in a single-band Hubbard model at finite temperature within QMC method, and more generally the difficulties in carrying out the DMFT-QMC algorithm , have not been clarified. This is an important problem, since DMFT-QMC is currently being applied to many similar models with more complicated orbital and spin structure.

In this paper we provide clear evidence for the coexistence of metallic and insulating solutions in a half-filled single-band Hubbard model with bare semicircular density of states at finite temperature within DMFT-QMC calculation . We discuss several aspects of DMFT-QMC method which have to be taken into account in order to get the reliable results near a coexistence region. First, we identify the source of the main difficulties in running DMFT-QMC near the coexistence region: uncontrolled errors in DMFT-QMC in the region where critical slowing down is dominant. The critical slowing down of the DMFT iteration follows from the Landau analysis  and can only be avoided by choosing a different rearrangement of the DMFT equation such as the one suggested by the Newton method. In this paper, however we use this critical slowing down as an indicator which allows us to locate boundaries of the coexistence region, and to fully confirm ,using QMC methods, the predictions of the Landau analysis. We discuss how to minimize errors in DMFT-QMC calculation by using correct interpolation of the discretized Green’s functions in direct Fourier Transformation($`FT`$) as well as a proper treatment of discontinuities of Green’s functions in inverse Fourier Transformation($`iFT`$). Also extrapolation of data to the limit of $`\delta \tau (=\beta /L)0`$ is used to reduce Trotter errors. Finally we obtain reliable coexistent insulating and metallic solutions within DMFT-QMC method and estimate the location of the first order MIT line $`U_c(T)`$.

In the limit of infinite dimensions a single-band Hubbard model is mapped onto an Anderson impurity model supplemented by a self-consistency condition . DMFT equation of this model,

$$t^2G(i\omega )[\mathrm{\Delta }]=\mathrm{\Delta }(i\omega ),$$
(1)
is solved iteratively within QMC calculation, where $`\mathrm{\Delta }(i\omega )`$ and $`G(i\omega )[\mathrm{\Delta }]`$ are a hybridization function and a local Green’s function of the impurity model. Our energy unit is set to be $`t=1/2`$ which comes from a scaled hopping amplitude, $`t/\sqrt{z}`$, in an infinite dimensional single-band Hubbard model.

Within DMFT-QMC calculation, the coexistence region should be explored with a great care due to the following two reasons: critical slowing down and errors in DMFT-QMC. The critical slowing down is the critical phenomena which manifests itself by slowing down of the convergence rate of a numerical method. The critical slowing down at a phase boundary can be proved with the argument based on Landau free energy : no convergence at an inflection point of Landau free energy within iteration method. At a boundary of the coexistence region one solution never converges while another one converges fast. However around Mott endpoint (where $`U_{c1}`$ and $`U_{c2}`$ merge) both solutions converge slowly because the narrow coexistence region around Mott endpoint brings two boundaries close enough so as to make a whole coexistence region numerically unstable as was pointed out in ref. . Unless carefully taken care of, fluctuation of a solution originating from errors in DMFT-QMC may result in undesirable transition of one solution to another one in DMFT iteration around Mott endpoint where critical slowing down is dominant.

In this paper we show that one can use the existence of critical slowing down as an indicator of boundaries of the coexistence region. Near a coexistence boundary, $`U_{c1}(T)`$ or $`U_{c2}(T)`$, metallic and insulating solutions show different convergence patterns within DMFT iterations: one converges very slowly while the other converges in a few iterations. A convergence criterion is applied to determine required number of iterations to get a convergent solution. We locate the coexistence boundaries at which such iteration numbers diverge.

The convergence of a solution is determined by monitoring of the most fluctuating component $`Im\mathrm{\Sigma }(i\pi T)`$, imaginary part of self energy at the first Matsubara frequency, as a function of number of iterations. QMC statistical error of $`Im\mathrm{\Sigma }(i\pi T)`$ at each iteration is controlled by the number of QMC sweeps while the change in $`Im\mathrm{\Sigma }(i\pi T)`$ between iterations is controlled by the number of iterations. The iteration procedure is stopped when the change in $`Im\mathrm{\Sigma }(i\pi T)`$ from iteration to iteration is less than the QMC statistical fluctuation of $`Im\mathrm{\Sigma }(i\pi T)`$, $`\delta _{QMC}Im\mathrm{\Sigma }\delta ImG/(ImG)^2`$. Notice that the QMC fluctuation of insulating solution is five to twenty times larger than QMC fluctuation of metallic solution as in Fig. 1 mostly because $`ImG_{ins}ImG_{met}`$. This convergence criterion is more strict than the conventional one, which requires to stop iteration procedure when the norm of $`Im\mathrm{\Sigma }^{(i+1)}Im\mathrm{\Sigma }^{(i)}`$ is smaller than a tolerance. For the conventional criterion underestimates the contribution of the most fluctuating component to the convergence by taking average of self energy differences between iterations over whole frequency range. The convergence criteria are illustrated in Fig. 1.

With this convergence criterion we performed the DMFT-QMC calculation to find solutions of a half-filled single-band Hubbard model at temperatures $`T`$=1/57, 1/64, 1/100, which turn out to be well below the critical temperature $`T_c`$. A careful selection of initial seeding to feed DMFT iteration is an essential part of finding of the coexistence at finite temperature. At a given temperature, the best educated seeding is used to start DMFT iteration at initial value of interaction strength $`U`$: small interaction strength $`U`$ for metallic solution and large $`U`$ for insulating solution. After DMFT-QMC run we obtain convergent metallic and insulating solutions within a few iterations. We increase interaction strength $`U`$ by small increment(0.01) for metallic solution and decrease $`U`$ for insulating one. To reduce computational time, the convergent solution at the previous value of interaction strength $`U`$ is used as initial seeding for DMFT-QMC run at the next value of interaction strength $`U`$. We observe

convergence pattern of solutions and repeat the procedures until critical slowing down is noticeably developed. By plotting inverse number of required iterations for the convergence of insulating solution in Fig. 2(a) and of metallic solution in Fig. 2(b), we identify upper bound of $`U_{c1}(T)`$ and lower bound of $`U_{c2}(T)`$. To find lower bound of $`U_{c1}(T)`$ and upper bound of $`U_{c2}(T)`$ we approach coexistence boundaries, $`U_{c1}(T)`$ and $`U_{c2}(T)`$, from outside the coexistence region. Outside the coexistence region only one solution exists. For interaction strength $`U(T)<U_{c1}(T)`$ DMFT iteration which starts with insulating seeding would produce convergent metallic solution. As $`U`$ gets closer to $`U_{c1}`$ from outside the coexistence region, convergence rate become slower. Such critical slowing down is clearly captured with open symbols in Fig. 2(a). By approaching coexistence boundaries from both inside and outside the coexistence region we obtain upper and lower bounds of $`U_{c1}`$ and $`U_{c2}`$ at several temperatures: for $`\beta =57`$, $`2.36U_{c1}(T)2.4`$ and $`2.42U_{c2}(T)2.53`$; for $`\beta =64`$, $`2.36U_{c1}(T)2.41`$ and $`2.48U_{c2}(T)2.52`$; for $`\beta =100`$, $`2.35U_{c1}(T)2.45`$ and $`2.58U_{c2}(T)2.62`$.

At $`\beta =57`$ the coexistence region is narrow enough so as to bring the critical slowing down into play over the whole coexistence region. Thus at temperatures above $`\beta =57`$ the coexistence region could exist as well but will be inapproachable numerically for it is hard to achieve convergent solutions.

All known errors in DMFT-QMC should be reduced as much as possible in order to obtain irrefutable evidence of the coexistence. Those errors are the QMC statistical errors and the systematic errors. The latter originate from two sources: the finite sampling size errors in both direct $`FT`$ and $`iFT`$ and the Trotter finite size errors.

To reduce the finite sampling size errors in direct $`FT`$ we use cubic spline interpolation between grid points of Green’s functions. The main idea of cubic spline interp-

olation is to get an interpolating function that has continuous first and second derivatives at grid points. Boundary conditions are necessary to close the set of liner equations to determine interpolation function. The usual way, when boundary conditions are unknown, is to put the second derivatives at end points to zero (natural spline). In the system we treat we can obtain boundary conditions from Green’s function. To do so, we expand Green’s function $`G(i\omega )`$ via k-th derivatives of $`G(\tau )`$ at end points: $`G(i\omega )`$=$`_{k=0}^{\mathrm{}}(1)^{k+1}(G^{(k)}(\beta ^{})+G^{(k)}(0^+))/(i\omega _n)^{k+1}`$. One can also obtain correct high frequency behavior of Green’s functions using the moments of the spectral density . Making correspondence between these two expansions we get the following formula for the second derivatives of Green’s functions in the case of half filled Hubbard model within DMFT: $`G^{(2)}(0^+)+G^{(2)}(\beta ^{})=(U^2/4+t^2)`$ for local Green’s function and $`G_o^{(2)}(0^+)+G_o^{(2)}(\beta ^{})=t^2`$ for Weiss function. The symmetry property of Green’s function at half filling is taken into account: $`G^{(k)}(0^+)=(1)^kG^{(k)}(\beta ^1)`$.

Despite of a lot of frequency points used to define Green’s function $`G(i\omega )`$, the finite frequency cutoff in $`iFT`$ results in removal of the discontinuity of $`G(\tau )`$ at end points. To correct this situation we should extract $`1/i\omega `$ tail from Green’s function before making $`iFT`$. To get correct $`iFT`$ of $`G(i\omega )`$ one should make numerical $`iFT`$ for $`G(i\omega )1/i\omega `$ function and analytical $`iFT`$ for the tail $`1/i\omega `$. In the case of half-filling the correct $`iFT`$ is: $`G(\tau )=T_ne^{i\omega _n\tau }[G(i\omega _n)\frac{1}{i\omega }]1/2`$ for $`\tau >0`$.

The corrections in direct $`FT`$ and $`iFT`$ are implemented in DMFT-QMC: correct boundary conditions in spline interpolation and tail correction in $`iFT`$. All artificial oscillations at both Im$`\mathrm{\Sigma }(i\omega )`$ and $`G(\tau )`$ presented in DMFT-QMC without corrections disappears in DMFT-QMC with corrections.

Having implemented the corrections in DMFT-QMC, we perform DMFT-QMC run to get the coexistence at numerically most stable points, i.e., at mid-points between the coexistence boundaries. The statistical errors are reduced by increase of QMC sweeps per site upto $`3\times 10^5`$. Several DMFT-QMC calculations are performed with different values of $`\delta \tau (=\frac{\beta }{L})`$. Trotter errors are removed by extrapolation of data to the limit of $`\delta \tau 0`$. In Fig. 3 it is shown convergence of local Green’s function $`G(\tau )`$ in the limit of $`\delta \tau 0`$. Insulating solutions are all within the QMC fluctuation $`(10^4)`$. However metallic solutions show monotonous increases with decreasing of $`\delta \tau `$. In Fig. 4 it is shown convergence of double occupancy, $`d=n_{}n_{}`$, as a function of $`\delta \tau ^2`$. The convergence of metallic and insulating double occupancies in the limit of $`\delta \tau 0`$ confirms the coexistence within DMFT-QMC. In the limit of $`\delta \tau 0`$ we obtain $`d_{met}0.039`$ and $`d_{ins}0.026`$ for $`\beta =64`$ and $`U=2.44`$; $`d_{met}0.036`$ and $`d_{ins}0.023`$ for $`\beta =100`$ and $`U=2.53`$. With reducing all other finite size errors except ones originating from Trotter approximation, the effect of Trotter errors on the convergence of solution is found to be small compared to errors in direct $`FT`$ and $`iFT`$.

We now use our QMC results together with previous numerical results obtained with other methods to draw the phase diagram in Fig. 5. To supplement our QMC data, we use the estimates of the Mott endpoint within QMC . We also include estimates for the coexistence boundaries  from Exact Diagonalization(ED). At zero temperature, the value of $`U_{c2}`$ was obtained with great accuracy in ref  using the projective self consistent method. Other exact methods such as the numerical renormalization group method(NRG)  give a very similar value for this quantity. The region near the Mott endpoint cannot be studied with numerically exact methods for the reasons discussed in this paper. However we can use the Landau theory to scale the results of the IPT near the critical region as described in ref . This approach was used to determine the $`U_{c1}(T)`$ and $`U_{c2}(T)`$

lines near the Mott endpoint as well as the location of the first order MIT line $`U_c(T)`$. Two crossover lines above the Mott endpoint  were scaled in a similar manner. The physical meaning of the crossovers were discussed in ref  in connection with experimental observations in $`V_2O_3`$. Landau MIT line and crossover lines are mapped onto QMC phase diagram by shifting of $`U_c(T_c)`$ and using the temperature scaling discussed in ref .

Near zero temperature, we determine the location of the finite temperature MIT line, by following ref  and assuming an entropy difference of $`ln(2)`$ between the metallic and the insulatining phase. The MIT lines obtained from those two methods agree well with MIT line obtained by ED method  which we also included in Fig. 5. Given the systematic finite size errors of the various numerical approximations, it is clear that there is quantitative consensus ( at the level of ten to twenty percent) between numerical calculations on where the coexistence region of fully frustrated Hubbard model in infinite dimensions lies. We note in passing, that the line of second order phase transitions reported in ref. , is actually crossover line at high temperature and the $`U_{c2}(T)`$ line at low temperature.

In conclusion, we find that QMC can be used to establish rigorously the existence of a coexistence region in the frustrated Hubbard model in infinite dimensions. Furthermore the numerical results are in full agreement with the predictions of the Landau theory of ref.  below $`T_c`$. Finally we notice that the methodology described in this paper should be useful when applying DMFT to other phase transitions in strongly correlated systems.

We thank G. Kotliar for many useful and fruitful discussions. This research was motivated by the stimulating discussions at the workshop on Theoretical Methods for Strongly Correlated Fermions at the CRM in Montreal. We have used the supercomputer at the Center for Advanced Information Processing in Rutgers.


text:
# Analysis of isoplanatic high resolution stellar fields by the StarFinder code

## 1  Introduction

StarFinder is a code developed within the frame of the ESO PAPAO (Currie et al. 2000a , 2000b ) program for the reduction of Adaptive Optics data. It was originally designed (Diolaiti et al. dio1 (1999)) to analyze AO images of very crowded fields, like the PUEO frame of the Galactic Center shown in this paper, which contains about 1000 detectable stars in a field of view of 13$`\times `$13 arcsec2. This image is an example of an AO observation: the PSF is Nyquist-sampled and characterized by a complex shape, showing a sharp peak, one or more fragmented diffraction rings and an extended irregular halo. Moreover, due to the small field of view, the imaged region is approximately isoplanatic and the PSF may be considered space invariant. Under the assumptions of isoplanatism and Nyquist-sampling, StarFinder models the observed stellar field as a superposition of shifted scaled replicas of the PSF lying on a smooth background due to faint undetected stars, possible faint diffuse objects and noise.

The procedure derives first a PSF digital template from the brightest isolated field stars; then a catalogue of presumed objects is formed, searching for the relative intensity maxima in the image frame. In the following step the images of the suspected stars are analyzed in order of decreasing luminosity. In this phase a catalogue including the accepted objects is formed and a synthetic image of the observed field is constructed, placing an intensity-scaled PSF template in the position of each identified star. Each suspected object in the original list is accepted on the basis of its correlation coefficient with the PSF template; the relative astrometry and photometry of the source are determined by means of a fit, taking into account the contribution of the local non-uniform background and of the already detected stars, known from the synthetic field. As the analysis proceeds, fainter and fainter sources can be successfully analyzed, discriminating their central peaks from the secondary bumps of light in the distorted diffraction rings of the neighboring more luminous already identified stars: in this way the synthetic field becomes more and more similar to the observed image. Residual unexplained features of the image may be further analyzed checking for indications of blended groups, at separations smaller than the PSF FWHM.

StarFinder should be regarded as a tool for high precision astrometry and photometry of crowded fields acquired under the above assumptions of accurate PSF knowledge, isoplanatism and correct sampling. An application of StarFinder to HST undersampled data handled by dithering techniques (Fruchter et al. fruchter (1997)) has been shown in Aloisi et al. (aloisi (2000)). In this paper it can be seen that the results obtained by our method are comparable to those obtained by DAOPHOT (Stetson stetson (1987)).

As far as the PSF is concerned, the accurate knowledge of its features outside the central peak is fundamental to perform a deep study of a stellar field, achieving accurate photometry of faint stars and avoiding either false detections or star loss (Esslinger & Edmunds esslinger (1998)). If the PSF template cannot be extracted directly from the field, due to extreme crowding or lack of bright isolated stars, StarFinder can still be applied using a PSF estimated by means of other methods, as the reconstruction technique proposed by V$`\stackrel{´}{e}`$ran et al. (vera1 (1997)) for AO observations or the TINY TIM simulation software for HST (Krist & Hook krist (1999)).

Much more intriguing and difficult to solve is the case of a field with space variant PSF, either due to anisoplanatic effects as in AO observations or to design and control in HST images. In general the analysis of an anisoplanatic field requires the knowledge of the local PSF. A method to reconstruct the off-axis PSF in AO imaging has been proposed by Fusco et al. (fusco (2000)). StarFinder, in its present version, can analyze frames with space invariant PSF or sub- frames in which the isoplanicity condition is nearly satisfyied, as will be shown in Section 4; the extension to the space variant case is in progress, with preliminary results presented by Diolaiti et al. (2000a , 2000b ).

The paper is organized as follows: the general features of the algorithm are described in Sect. 2; Sect. 3 deals with more technical aspects and might be skipped on a first reading; the method is validated on simulated and experimental data in Sect. 4; details on the IDL code are presented in Sect. 5; Sect. 6 includes our conclusions and future plans.

## 2  Code description

### 2.1  PSF estimation

The accuracy of the PSF estimate is primary in StarFinder, since the PSF array is used as a template for all the stars in the image to be analyzed. The user selects the most suitable stars, which are background subtracted, cleaned from the most contaminating sources around, centered with sub-pixel accuracy, normalized and superposed with a median operation. The centering is performed by an iterative shift of the stellar image in order to cancel the sub-pixel offset of its centroid (see Christou & Bonaccini christou (1996)). The median operation, which is performed pixel-by-pixel, is preferred to the mean because it is less sensitive to anomalous pixels (outliers) which might appear in one or more stellar images among the selected ones. The retrieved PSF is post-processed in order to reject any residual spurious feature and to smooth the noise in the extended halo.

The PSF extraction procedure also reconstructs the core of saturated stars by replacing the corrupted pixels with the central part of the PSF estimate. Accurate positioning is achieved by means of a cross-correlation technique, while the scaling factor is determined with a least squares fit to the wings of the star to repair. For a detailed description of the procedure, see Sect. 3.5.

### 2.2 Standard analysis of a stellar field

At first we build a list of objects, the presumed stars, which satisfy the condition

$$i(x,y)>b(x,y)+t$$
(1)
where $`i(x,y)`$ is the observed intensity, $`b(x,y)`$ the background emission and $`t`$ a fair detection threshold. The presumed stars are analyzed one by one by decreasing intensity. To illustrate a generic step of the analysis, we consider the $`(n+1)th`$ object in the list, after the examination of the first $`n`$. A small sub-image of fixed size is extracted around the object (Fig. 1, left). This sub-image may contain brighter stars formerly analyzed, fainter objects neglected in the current step and features of other stars lying outside the sub-image. The information on the brighter sources is recorded in a synthetic stellar field (Fig. 1, right), defined as the sum of two terms: a superposition of PSF replicas, one for each star detected up to this point, and an estimate of the background, assumed to be non uniform in general.

The local contribution due to the brighter stars and the background, derived from the synthetic field, is subtracted from the sub-image. If a statistically significant residual remains, it is compared to the central core of the PSF by means of a correlation check. If the correlation coefficient is higher than a pre-fixed threshold then the object of interest is rated similar to the PSF and accepted. The accurate determination of its position and relative flux is attained by means of a local fit, in which the observed sub-image is approximated with the multi-component model described in Sect. 3.6. The actual size of the fitting region is comparable to the diameter of the first diffraction ring of the PSF. This choice ensures that the information represented by the shape of a high-Strehl PSF is included in the fitting process to achieve better accuracy and prevents the growth of the number of sources to be fitted together. For the central object of our example a single component fit is performed and the contribution due to the brighter stars is considered as a fixed additive term. A multi-component fit is performed when the star is in a very compact group, at separations comparable to the PSF FWHM. If the fit is acceptable, the parameters of the new detected star are saved and those of the already known sources, which have been possibly re-fitted, are upgraded. The new star and an upgrade of the re-fitted sources are added to the synthetic field.

This analysis is performed for each object in the initial list (a flow-chart illustrating the operations of StarFinder is in Figure 2). To achieve a better astrometric and photometric accuracy, at the end of the examination of all the objects, the detected stars are fitted again, this time considering all the known sources. This step may be iterated a pre-fixed number of times or until a convergence condition is met.

At the end of the analysis, it is possible to stop the algorithm or instead perform a new search for lost objects removing the detected stars and using an upgraded background estimate. It should be stressed that this image subtraction is just a tool to highlight significant residuals. Any further analysis is performed on the original frame, in order to take into account the effects arising from the superposition of the PSFs of neighboring sources. Generally, after 2-3 iterations of the main loop, the number of detected stars approaches a stable value.

### 2.3 Crowding and blending effects

A binary star with different separation values (Fig. 3) has been simulated to show how the code works with crowded sources. With a separation of 2 PSF FWHM the two components are well separated and the code analyzes them with the standard procedure described in the previous sub-section. In the other cases (separation from 1 to 0.5 PSF FWHM) the secondary component is not detected as a separate relative maximum and it is lost. However, if the separation is not as small as 0.5 FWHM, a further iteration of the main loop enables the algorithm to detect the fainter component by subtracting the brighter one. This strategy forces the two stars to pass the correlation test, the principal component as a presumed single object and the secondary in a subsequent iteration of the loop. In a way the iteration of the main loop is a de-blending strategy, because it enables the algorithm to detect stars whose intensity peak is not directly visible in the observed data.

This strategy fails when:

* the binary is very noisy and the two components have similar flux. In a similar situation the principal component may not pass the correlation test, preventing also the detection of the secondary in a further iteration. These noisy blends may be recovered at the end of the overall analysis.
* The two components are almost equally bright and have a separation close to the lower limit (1/2 PSF FWHM). The residual corresponding to the secondary after subtracting the principal component may have a distorted shape and might not pass the correlation check.

The latter case may be handled by a method based on a thresholding technique. The object is cut at a prefixed level, about 20% below the central peak, and transformed to a binary array, setting to 1 all the pixel above the threshold and to 0 the pixels below. If the area of the pixels with value 1 is more extended than the PSF, the object is considered a blend and the secondary star may be detected by subtracting the brighter one; then a two-components fit allows accurate astrometry and photometry of the two sources. This strategy can be iteratively applied to multiple blends. It should be stressed that the area measurement is not reliable when the value of the cutting threshold is comparable to the noise level: for this reason the de-blending procedure is applied only to objects with a suitable signal-to-noise ratio. Moreover the area measurement is reliable when the data are adequately sampled. This de-blending procedure is applied at the end of the last iteration of the main loop, when all the resolved sources have been detected: in this way the probability that a single object may appear artificially blurred because of the contamination of still unknown sources is largely reduced.

In a normal case, like the simulated field of Sect. 4.1, two or three iterations of the main loop find almost all ($``$99%) the stars that StarFinder may detect. The de-blending procedure described above adds $``$1% more stars, without additional false detection. The number of lost objects belonging to the first category described above is negligible ($`<<`$1%).

Normally we perform two or three iterations of the main loop and apply the de-blending strategy only in very crowded fields.

## 3 Technical aspects

### 3.1 Bad pixels

StarFinder includes an automatic procedure to repair known bad pixels, which are replaced with the median of the good data points in a suitable neighborhood. The use of this procedure is particularly recommended in the extraction of the PSF. Indeed the stellar images, which will be combined to form the PSF, are centered with sub-pixel accuracy by means of an interpolation technique and the presence of bad pixels may produce unpredictable interpolation errors. Even if replaced, the bad pixels are excluded from any further computation concerning the correlation coefficient and the fitting process.

### 3.2 Background estimation

A reliable estimate of the background is necessary to define the detection condition and to compute the correlation coefficient of the presumed stars with the PSF.

A straightforward technique to estimate the background is to smooth the image by median filtering, replacing each pixel with the median computed over a suitable neighborhood, of size comparable to the characteristic width of a stellar image. This method tends to over-estimate the background underneath strong peaks. A more accurate approximation, even in crowded fields, is described in Bertin & Arnouts (bertin (1996)). The image is partitioned in sub-regions arranged in a regular grid and a local estimate is calculated for each sub-patch by means of an IDL implementation of the DAOPHOT SKY routine (Stetson stetson (1987)), due to Landsman (landsman (1995)). This array of sky measurements is smoothed by median filtering and interpolated onto the same grid of the input image.

It should be stressed that the background computation is unavoidably affected by the presence of bright sources. In general a more accurate estimate can be obtained after the analysis of the stellar field, when the most contaminating sources are known and can be subtracted.

### 3.3  Noise estimation

The estimate of the noise is useful to define the detection threshold and to compute the formal errors on the retrieved astrometry and photometry.

The overall effect of photon and instrumental noise can be computed if the required parameters (detector gain, read-out-noise, dark current, etc.) are known. Otherwise an estimate of the mean background noise can still be obtained by means of histogram fitting techniques (Almoznino et al. almoznino (1993), Bijaoui bijaoui (1980)). Assuming that the intensity of the sky radiation is distributed normally around a typical value, the histogram of the observed intensity levels should be quite similar to a gaussian distribution, whose mode and standard deviation represent respectively the sky value and the associated noise. Actually the contamination due to stellar sources produces a high-intensity tail and an artificial broadening of the histogram, which prevents an accurate estimate of the background noise. This problem can be partially overcome by removing the signal from the image, leaving only the pixel-to-pixel variations associated to pure noise: a reasonable estimate of the signal to subtract for this purpose can be obtained by smoothing the data with a median filtering technique. After this operation the histogram is symmetric around its mode and the background noise standard deviation can be estimated by means of a gaussian fit to the histogram itself.

### 3.4  Correlation coefficient

False detections, associated to noise spikes or residual PSF features of bright stars, are recognized on the basis of their low correlation coefficient with the PSF, which represents a template for each true star in the image. The correlation coefficient (Gonzalez & Woods gonzales (1992)) is defined as

$`c(a,b)=`$ (2)
$`{\displaystyle \frac{_{x,y}\left[i(x,y)\overline{ı}\right]\left[p(xa,yb)\overline{p}\right]}{\sqrt{_{x,y}\left[i(x,y)\overline{ı}\right]^2}\sqrt{_{x,y}\left[p(xa,yb)\overline{p}\right]^2}}}`$
where $`i(x,y)`$ and $`p(x,y)`$ are the object and the PSF respectively, $`\overline{ı}`$ and $`\overline{p}`$ are the corresponding mean values. Maximizing the correlation coefficient as a function of the offset $`(a,b)`$ yields an objective measure of similarity. After maximizing $`c(a,b)`$ for integral offsets, it is possible to repeat the procedure for sub-pixel shifts, improving the positioning accuracy.

The correlation coefficient is computed on the core of the PSF: typically the central spike of the diffraction pattern is considered, out to the first dark ring. A fair correlation threshold must be fixed in order to discriminate and reject unlikely detections, without losing faint stars contaminated by the background noise; a value of 0.7 or 0.8 is acceptable in most cases.

The correlation coefficient represents also an effective tool to select the stars with the highest photometric reliability, since generally a very high correlation value is associated to resolved single sources.

### 3.5  Saturated stars

Saturated stars provide precious information on the PSF halo, so it may be useful to include them, appropriately reconstructed, in the PSF extraction process. In addition, the repaired saturated stars can be recognized by the star detection algorithm and their contribution be taken into account during the analysis of near and fainter sources.

The core of a saturated star is replaced with a shifted scaled replica of a preliminary estimate of the PSF. The repaired star is defined as

$$i_{rep.}(x,y)=\{\begin{array}{cc}i(x,y)\hfill & \text{ if }i(x,y)<T\hfill \\ fp(xx_0,yy_0)\hfill & \text{otherwise}\hfill \end{array}$$
(3)
where $`T`$ is the upper linearity threshold of the detector. The position of the center $`(x_0,y_0)`$ is estimated by means of the maximization of the correlation coefficient, which is not sensitive to the intensity levels. The scaling factor $`f`$ is calculated with a least squares fit to the wings of the saturated star. Of course the saturated pixels are excluded from all the computations. The background is temporarily subtracted in order to prevent affecting the repair process.

This procedure has been applied to the brightest star (IRS 7) of the Galactic Center image shown in Sect. 4.2. This source is not saturated in the original data, but it has been artificially corrupted by an upper cut at half maximum. The repair procedure is able to reconstruct the original peak with an error of $``$5%, imposing a positioning accuracy of 1/2 pixel.

### 3.6 Fitting procedure

A sub-image centered on the star of interest is extracted and approximated with the model

$`h(x,y)`$ $`=`$ $`s_0(x,y)+{\displaystyle \underset{n=1}{\overset{N_s}{}}}f_np(xx_n,yy_n)`$ (4)
$`+b_0+b_1x+b_2y`$
where $`s_0(x,y)`$ is the fixed contribution of known stars outside the sub- image support, $`N_s`$ is the number of point sources within the sub-image, $`x_n`$, $`y_n`$, $`f_n`$ are the position and flux of the $`nth`$ source, $`p(x,y)`$ is the PSF and $`b_0`$, $`b_1`$, $`b_2`$ are the coefficients of a slanting plane representing the local background. It should be stressed that the retrieved astrometry and photometry are referred to the absolute centering and normalization of the PSF array. The optimization of the parameters is performed by minimizing the least squares error between the data and the model. If the noise is known, it is possible to weigh the data by their inverse standard deviation and obtain a statistically optimal fit (Beck & Arnold beck (1977)). In this case, formal error estimates on the parameters can be obtained (Bevington & Robinson bevington (1994)). The optimization is performed by means of an iterative Newton-Gauss technique with linearized Hessian (Beck & Arnold beck (1977), Luenberger luen (1984) ). The Moore-Penrose generalized inverse concept (Rust & Burrus rust (1972), Bendinelli et al. bendinelli (1987), Lorenzutta lorenzutta (1967)) is applied to invert the Hessian matrix at every step. The iterations are stopped when the parameters approach a stable value. The major difficulty of a Newton-like method is represented by the computation of the model derivatives with respect to the parameters, some of which (stellar positions) yield non-linear dependence. For this purpose the Fourier shift theorem is applied:

$`p(xx_n,yy_n)=`$ (5)
$`FT^1\left[FT\left[p(x,y)\right]e^{i2\pi (ux_n+vy_n)/N}\right]`$
where $`FT`$ represents the discrete Fourier transform operation, $`N`$ is the sub-image size and $`u`$, $`v`$ are spatial frequencies. The derivative with respect to $`x_n`$ is then

$`{\displaystyle \frac{p(xx_n,yy_n)}{x_n}}=`$ (6)
$`FT^1\left[FT\left[p(x,y)\right]e^{i2\pi \left(ux_n+vy_n\right)/N}\left(i{\displaystyle \frac{2\pi u}{N}}\right)\right]=`$
$`FT^1\left[\left(i{\displaystyle \frac{2\pi u}{N}}\right)FT\left[p(xx_n,yy_n)\right]\right]`$
and requires in practice an interpolation of the PSF to compute $`p(xx_n,yy_n)`$ . In principle this interpolation-based method can only be applied to Nyquist-sampled images and this is currently the main limit of the algorithm. A similar technique has been described by V$`\stackrel{´}{e}`$ran & Rigaut (vera2 (1998)).

### 3.7  Sampling and interpolation

A band-limited function in one dimension is Nyquist-sampled if the step size fulfills the condition

$$\mathrm{\Delta }x\frac{1}{2}f_c$$
(7)
where $`f_c`$ is the so-called cut-off frequency of the spectrum. In this case it is possible to reconstruct the original continuous function from a set of equally-spaced samples by means of the so-called sinc interpolation (Mariotti mariotti (1988)). For a nearly diffraction-limited PSF, the critical sampling condition is generally stated by saying that the its FWHM must contain at least two pixels. Many accurate and efficient interpolation schemes exist to perform the PSF shift required by the fitting procedure. Probably the most simple is a straightforward application of the Fourier shift theorem, but one may also use bicubic splines or the fast sinc interpolation algorithm described in Yaroslavksy (yaro (1997)). In StarFinder we have applied the cubic convolution interpolation method (Park & Schowengerdt park (1983)), which approaches very closely the optimal sinc interpolation for Nyquist-sampled data; this algorithm is implemented in the IDL function INTERPOLATE.

Even though the interpolation of the PSF array is allowed only on Nyquist-sampled data, the cubic convolution technique seems to produce acceptable results even with marginally under-sampled images. Tests performed on Airy diffraction patterns, with a sampling step twice as large as the critical sampling step size, indicate that the interpolation-induced oscillations amount to a few percent of the image peak, as opposed to 10-20% of other interpolation techniques like Fourier shift or bicubic splines.

## 4  Applications

### 4.1 Simulated field

The code was first applied to a simulated image (Fig. 4) including 1000 stars, placed randomly in the frame and with a given magnitude distribution (Fig. 6). Each star is a scaled copy of a long exposure high-Strehl PSF, obtained with the ADONIS AO system at the ESO 3.6 m telescope. The image is 368 $`\times `$ 368 pixels large ($`13\mathrm{}\times 13\mathrm{}`$) and has a stellar density of 6 stars arcsec-2. Photon, readout noise and a background nebulosity, normalized to the same flux of the stellar sources, were added to the image. The faintest stars have a peak signal-to-noise ratio of 5.

We performed a standard reduction of the artificial image using the ”default” parameters of the method.

The PSF was estimated by superposing the images of the four brightest stars in the field. The retrieved PSF is very similar to the true one (Fig. 5).

Figure 6 shows the good agreement between the true and the observed luminosity function. The sole statistically significant discrepancy is in the bin from magnitude 7 to 8, where $``$25% of the stars were lost; the other small differences are due to photometric errors which shift some objects to a neighboring magnitude interval. The lost stars are $``$10% of the total number of sources and are generally faint: about 90% have magnitude between 7 and 8, the rest is in the bin between magnitude 6 and 7. The only bright lost star has magnitude $``$4.5 and is the secondary component of a very close binary, with a separation of just 1/2 pixel.

Roughly 70% of the lost stars are located at a distance $``$ 1 PSF FWHM from the nearest object, $``$20% are on the first diffraction ring of a brighter source and only $``$10% are isolated. It should be stressed however that $``$15% of the lost stars can be recognized by visual inspection as faint objects in the halo of the four brightest stars in the field, independently of their separation from the nearest source. It is apparent that the blending effect and the contamination by the halo of very bright stars do account for the lost stars. Note that the number of false detections in this simulated field is negligible (1 case out of 1000).

The plots in Figures 7, 8 show the astrometric and photometric accuracy of StarFinder. About 80% of the detected stars have both astrometric error smaller than 0.1 PSF FWHM and photometric error smaller than 0.1 magnitudes. The stars with less accurate astrometry or photometry are generally faint: $``$75% belong to the magnitude interval from 7 to 8, $``$20% to the bin between magnitude 6 and 7 and only $``$5% are distributed in lower bins of the luminosity function.

These stars have, in $``$45% of the cases, a lost source in their immediate neighborhood within the first diffraction ring of the PSF. The others are faint stars dispersed in the halo of the brightest sources or in the noisy background nebulosity.

After discussing the performance of the code with a standard analysis, it is interesting to examine how the results are affected by the main parameters of the method. Applying the de-blending strategy we detected $``$+10% more binaries in the separation range between 1/2 and 1 PSF FWHM, even though the overall detection gain is less than 1% referred to the total number of sources. With a higher number of iterations of the main loop we detected $``$+30% more binaries in the range from 1 to 2 PSF FWHM. Decreasing the detection threshold from 3 to 1 times the noise standard deviation, we found $``$+25% more binaries in the range between 1/2 and 1 PSF FWHM, but with 10 faint (mag $`>8`$) false detections instead of 1. Increasing the threshold on the correlation coefficient, from 0.7 to 0.8, we reported no false detection, but the number of lost stars increased by about 60%; the additional lost sources were generally fainter than magnitude 7, but not necessarily in crowded groups. Lowering the correlation threshold to 0.6 we detected more faint isolated stars and binaries, at separations between 1 and 2 PSF FWHM, but with a higher probability of false detections (2 instead of 1). Finally the astrometric and photometric accuracy approaches a stable level after a few ($``$2) re-fitting iterations.

### 4.2  Galactic Center

The code was run on a 15 min exposure time K band (2.2$`\mu `$m) image of the Galactic Center (Fig. 9), taken with the PUEO AO system on the 3.6m CFH telescope (Rigaut et al. rigaut (1998)). The Strehl ratio in the image is $``$45%. The PSF FWHM is $`0.13\mathrm{}`$, with a sampling of $`0.034\mathrm{}`$ /pixel. The Adaptive Optics guide star was a $`m_R`$=14.5 star (called star 2 in Biretta et al. biretta (1982)) located about (to the upper left) $`20\mathrm{}`$ from the center of the image, out of the field of view of the figure. There is therefore a slightly elongation of the PSF towards the direction of the guide star. However, due to the fact that the isoplanatic patch was much larger than the $`13\mathrm{}\times 13\mathrm{}`$ shown in the figure, a space-invarient PSF fits the data very well, as we will show.

A standard analysis was performed, analogous to the one described in Sect. 4.1 for the synthetic stellar field. About 1000 stars were detected, with a correlation coefficient of at least 0.7; the reconstructed image is shown in Figure 10.

We evaluated the accuracy of the algorithm by means of an experiment with synthetic stars. We created 10 frames adding to the original image 10% of synthetic stars at random positions for each magnitude bin of the estimated luminosity function. The 10 frames were analyzed separately. As in the simulated case, a distance tolerance of 1 PSF FWHM was adopted to find coincidences between the detected stars and their true counterparts. The lists of detected artificial stars were merged together and the astrometric and photometric errors were computed and plotted as a function of the true magnitude (Figs. 11, 12).

The plots show no apparent photometric bias and high astrometric and photometric accuracy: the stars brighter than magnitude 5, for instance, have a mean astrometric error $`<0.5`$ mas and a mean absolute photometric error $`<0.01`$ mag. It should be stressed however that the artificial sources are contaminated by the background noise present in the observed data and by the photon noise due to neighboring stars, but no additional noise was added. Figure 13 shows a comparison between the mean luminosity function retrieved in the 10 experiments and the truth. Assuming an expected error for each bin equal to the square root of the corresponding number of counts, according to the Poisson statistic, the only significant differences occur for magnitudes fainter than 9.

It is also interesting to consider the magnitude distribution of the false detection cases (dashed-dotted line in Figure 13), i.e. the detected stars for which we found no counterpart in the original list, within a distance of 1 PSF FWHM. The mean percentage of false detections in the 10 experiments is 2% of the total number of stars. The false detections are almost always very faint (mag $`>8`$); their number is comparable to the square root of the total counts only in the last magnitude bin, for magnitudes fainter than 10. The percentage of false detections reported in these experiments seems to confirm the analysis performed by visual inspection on the stars detected in the original frame.

## 5  The StarFinder code

The StarFinder code has been provided with a collection of auxiliary routines for data visualization and basic image processing, in order to allow the user to analyze a stellar field, produce an output list of objects and compare different lists, e.g. referred to different observations of the same target. The input image is supposed to be just calibrated.

The code is entirely written in the IDL language and has been tested on Windows and Unix platforms supporting IDL v. 5.0 or later. A widget-based graphical user interface has been created. The main widget appearing on the computer screen is nothing more than an interface to call secondary widget-based applications, in order to perform various operations on the image. The basic documentation about the code can be found in the on-line help pages. IDL users might wish to run interactively the StarFinder routines, without the widget facilities: complete documentation on each module is available for this purpose.

## 6  Conclusions and future developments

The elaboration of real and simulated data seems to prove the effectiveness of StarFinder in analyzing crowded isoplanatic stellar fields characterized by high Strehl ratio PSFs and correct sampling, reaching in this case the full utilization of the data information content. The code can be applied also to low Strehl data with results comparable to those attainable by other methods. Moreover we are evaluating its performance on undersampled images.

StarFinder is also reasonably fast: the analysis of the Galactic Center image ($`368\times 368`$ pixels, $``$1000 stars) requires between 5 and 10 minutes on a normal PC (Pentium Pro - 64Mb RAM - 350MHz). The graphic interface makes it accessible to users unfamiliar with IDL.

In the near future we plan to provide the code with tools to handle images with spatially variable PSF and new elements which might be helpful in the analysis of a stellar field. An interesting application, still in progress, aims at cleaning mixed fields from the contamination of foreground stellar images, leaving out true or suspected diffuse sources.

The StarFinder package and its technical documentation are directly available at: http://www.bo.astro.it or can be obtained on request to E. Diolaiti.

###### Acknowledgements.

François Rigaut is acknowledged for kindly providing the PUEO image of the Galactic Center and for supporting the initial development of this method. The ESO La Silla 3.6m telescope team is acknowledged for the support during the technical observations with Adonis, from which we have retrieved ’real’ PSFs used to simulate stellar fields. This work was partially supported by the ESO Instrumentation Division Adaptive Optics Dept. and by the Italian Ministry for University and Research (MURST) under grant Cofin 98-02-32.