format: text

text:
# 1 The shifts 𝒎̃_𝒃-𝒎_𝒃 in the quark mass (a) and 𝝁_𝝅^𝟐-𝝁̃_𝝅^𝟐 in the kinetic operator (b) introduced by imposing a lower cut in the photon energy in 𝑩→𝑿_𝒔+𝜸. Blue and maroon curves correspond to two different ansätze for the heavy quark distribution function, 𝑭_𝟏 and 𝑭_𝟐, respectively.

Bicocca-FT-03-21

UND-HEP-03-BIG 03

hep-ph/0308165

On Extracting Heavy Quark Parameters
from Moments with Cuts

Contributed to the Lepton-Photon Conference, August 11-15, 2003
FNAL, Batavia IL

I.I. Bigi$`^𝒂`$ and N. Uraltsev$`^𝒃\mathbf{}`$

a Department of Physics, University of Notre Dame du Lac

Notre Dame, IN 46556, USA
b INFN, Sezione di Milano, Milan, Italy

Abstract

We point out that applying the photon energy cut significantly modifies the moments of energy spectrum in $`𝑩\mathbf{}𝑿_𝒔\mathbf{+}𝜸`$ decays, with a certain class of effects not accounted for in the mostly used OPE expressions. This leads to a systematic bias in the extracted values of the $`𝒃`$ quark mass and other heavy quark parameters. The apparent $`𝒃`$ quark mass increases typically by $`\mathrm{𝟕𝟎}\text{MeV}`$ or more, together with an even more dramatic downward shift in the kinetic expectation value. Accounting for these cut-related shifts brings different measurements into a good agreement, when the OPE-based theory employs the robust approach. These nonperturbative effects are exponential in the effective hardness severely lowered by high cuts, and do not signify a breakdown of the $`\mathbf{\hspace{0.17em}1}\mathbf{/}𝒎_𝒃`$ expansion itself. Similar effects in semileptonic $`𝒃\mathbf{}𝒄`$ decays are briefly addressed. We stress the utility of the second moment of $`𝑬_𝜸`$ once these effects are incorporated.

On leave of absence from Department of Physics, University of Notre Dame, Notre Dame, IN 46556, USA
and St. Petersburg Nuclear Physics Institute, Gatchina, St. Petersburg 188300, Russia

Treating effects of strong interactions in decays of beauty particles is of primary importance nowadays when many precision experimental studies are being carried out. Inclusive distributions in radiative and semileptonic decays are the portal to accurately determining the nonperturbative heavy quark parameters controlling short-distance observables in $`𝑩`$ decays. Their utility rests on an consistent expansion in $`\mathrm{𝟏}\mathbf{/}𝒎_𝒃`$, the inverse b quark mass . Recent experimental data generally show a nontrivial agreement between quite different – and a priori unrelated – measurements at the nonperturbative level, on one hand, and consistency with the QCD-based OPE theory.

In order to enjoy in full the potential of a small expansion parameter provided by the heavy quark mass, the observable in question must be sufficiently inclusive. However experimental cuts imposed for practical reasons – to suppress backgrounds etc. – often essentially degrade the effective hardness $`𝓠`$ of the process. This brings in another expansion parameter $`\mathrm{𝟏}\mathbf{/}𝓠`$ effectively replacing $`\mathrm{𝟏}\mathbf{/}𝒎_𝒃`$ in certain QCD effects. The reliability of the expansion greatly deteriorates for $`𝓠\mathbf{}𝒎_𝒃`$. This phenomenon is particularly important in $`𝒃\mathbf{}𝒔\mathbf{+}𝜸`$ decays where experiments so far have imposed cut $`𝑬_𝜸\mathbf{>}\mathrm{𝟐}\text{GeV}`$ or even higher.

The theoretical aspects of such limitations have been discussed during the last couple of years . In particular, the effective ‘hardness’ of the inclusive $`𝒃\mathbf{}𝒔\mathbf{+}𝜸`$ decays with $`𝑬_𝜸\mathbf{>}\mathrm{𝟐}\text{GeV}`$ amounts to only about $`\mathbf{1.25}\text{GeV}`$, which casts doubts on the precision of the routinely used expressions incorporated into the fits of heavy quark parameters.

In the present letter we point out that these effects are numerically significant, lead to a systematic bias that often exceeds naive error estimates and therefore cannot be ignored. Evaluating them in the most straightforward (although somewhat simplified) way we find, for instance for $`𝒃\mathbf{}𝒔\mathbf{+}𝜸`$ decays

$`\stackrel{\mathbf{~}}{𝒎}_𝒃`$ $`\text{ }\mathbf{}`$ $`𝒎_𝒃\mathbf{+}\mathrm{𝟕𝟎}\text{MeV}`$
$`\stackrel{\mathbf{~}}{𝝁}_𝝅^\mathrm{𝟐}`$ $`\text{ }\mathbf{}`$ $`𝝁_𝝅^\mathrm{𝟐}\mathbf{}\mathbf{(}\mathbf{0.15}\mathbf{÷}\mathbf{0.2}\mathbf{)}\text{GeV}^\mathrm{𝟐}`$ (1)
where $`\stackrel{\mathbf{~}}{𝒎}_𝒃`$ and $`\stackrel{\mathbf{~}}{𝝁}_𝝅^\mathrm{𝟐}`$ are the apparent values of the $`𝒃`$ quark mass and of the kinetic expectation value, respectively, as extracted from the $`𝒃\mathbf{}𝒔\mathbf{+}𝜸`$ spectrum with $`𝑬_𝜸\mathbf{>}\mathrm{𝟐}\text{GeV}`$ in a usual way. Correcting for these effects eliminates alleged problems for the OPE in describing different data and rather leads to a too good agreement between the data on different types of inclusive decays.

Moreover, this resolves the controversy noted previously: while the values of $`\overline{𝚲}`$ and $`𝝁_𝝅^\mathrm{𝟐}`$ reportedly extracted from the CLEO $`𝒃\mathbf{}𝒔\mathbf{+}𝜸`$ spectrum were found to be significantly below the theoretical expectations, the theoretically obtained spectrum itself turned out to yield a good description of the observed spectrum when we evaluated it based on these theoretically preferred values of parameters .

The bias in Eqs. (1) depends on the position of the cut (more precisely, on the gap $`\frac{𝒎_𝒃}{\mathrm{𝟐}}\mathbf{}𝑬_{\mathrm{𝐜𝐮𝐭}}`$) and the actual values of other heavy quark parameters. The quoted estimates assume moderate values, $`𝒎_𝒃\mathbf{(}\mathrm{𝟏}\text{GeV}\mathbf{)}\mathbf{}\mathbf{4.6}\text{GeV}`$ and $`𝝁_𝝅^\mathrm{𝟐}\mathbf{(}\mathrm{𝟏}\text{GeV}\mathbf{)}\mathbf{}\mathbf{0.43}\text{GeV}^\mathrm{𝟐}`$. If $`𝒎_𝒃`$ becomes lower and/or the true $`𝝁_𝝅^\mathrm{𝟐}`$ increases as may be indicated by the most recent data, the bias increases further.

OPE and cuts The origin of these effects and why they are missed in the standard application of the OPE and in estimates of the theoretical accuracy, have been discussed elsewhere . In brief, considering a constrained fraction of the $`𝑩\mathbf{}𝑿_𝒔\mathbf{+}𝜸`$ events

$$\mathrm{𝟏}\mathbf{}𝚽_𝜸\mathbf{(}𝑬\mathbf{)}\mathbf{=}\frac{\mathrm{𝟏}}{𝚪_{𝒃𝒔𝜸}}\mathbf{}_𝑬^{\frac{𝑴_𝑩}{\mathrm{𝟐}}}𝐝𝑬_𝜸\frac{𝐝𝚪_{𝒃𝒔𝜸}}{𝐝𝑬_𝜸}\mathbf{.}$$
(2)
(or similarly truncated photon energy moments), the simple-minded approach routinely expands the spectrum in powers of $`\mathrm{𝟏}\mathbf{/}𝒎_𝒃`$. Ignoring perturbative bremsstrahlung one obtains a $`𝜹`$-like spectrum peculiar for two-body decays, and the expansion around the free-quark kinematics does not change this – it only generates higher derivatives of $`𝜹\mathbf{(}𝑬_𝜸\mathbf{}\frac{𝒎_𝒃}{\mathrm{𝟐}}\mathbf{)}`$:

$$\frac{\mathrm{𝟏}}{𝚪_{𝒃𝒔𝜸}^\mathrm{𝟎}}\frac{𝐝𝚪_{𝒃𝒔𝜸}}{𝐝𝑬_𝜸}\mathbf{=}𝒂𝜹\mathbf{(}𝑬_𝜸\mathbf{}\frac{𝒎_𝒃}{\mathrm{𝟐}}\mathbf{)}\mathbf{+}𝒃𝜹^{\mathbf{}}\mathbf{(}𝑬_𝜸\mathbf{}\frac{𝒎_𝒃}{\mathrm{𝟐}}\mathbf{)}\mathbf{+}𝒄𝜹^{\mathbf{\prime \prime }}\mathbf{(}𝑬_𝜸\mathbf{}\frac{𝒎_𝒃}{\mathrm{𝟐}}\mathbf{)}\mathbf{+}\mathbf{}$$
(3)
where $`𝒂`$, $`𝒃`$, … are given by the $`𝑩`$ meson expectation values of local $`𝒃`$-quark operators. Naively computing $`\mathrm{𝟏}\mathbf{}𝚽_𝜸\mathbf{(}𝑬\mathbf{)}`$, or spectral moments over the restricted domain in this way would yield unity in Eq. (2) for any $`𝑬\mathbf{>}\frac{𝒎_𝒃}{\mathrm{𝟐}}`$ – a result clearly meaningless on physical grounds. The actual behavior of the spectrum and the moments is described by the heavy quark distribution function. Its tail is indeed exponentially suppressed by a typical factor $`𝒆^{\mathbf{}𝒄𝓠\mathbf{/}𝝁_{\mathrm{𝐡𝐚𝐝𝐫}}}`$ at $`𝓠\mathbf{(}𝑬_𝜸\mathbf{)}\mathbf{}𝝁_{\mathrm{𝐡𝐚𝐝𝐫}}`$. However, for $`𝓠\mathbf{(}𝑬_𝜸\mathbf{)}\mathbf{}𝝁_{\mathrm{𝐡𝐚𝐝𝐫}}`$ there is little suppression of the missed tail contribution; the error becomes of order one.

This obvious point is missed in the naive application of the OPE and in the way to gauge the theoretical uncertainty based on it. Conceptually this is related to the limited range of convergence of the OPE for the width, determined in this case by the support of the heavy quark distribution function .

Here we rather concentrate on the numerical consequences for $`𝑩\mathbf{}𝑿_𝒔\mathbf{+}𝜸`$ decays. To this end we first turn off perturbative effects altogether. The spectrum then is described by the nonperturbative light-cone distribution function $`𝑭\mathbf{(}𝒌_\mathbf{+}\mathbf{)}`$:

$$\frac{\mathrm{𝟏}}{𝚪}\frac{𝐝𝚪}{𝐝𝑬_𝜸}\mathbf{=}\mathrm{𝟐}𝑭\mathbf{(}\mathrm{𝟐}𝑬_𝜸\mathbf{}𝒎_𝒃\mathbf{)}\mathbf{.}$$
(4)
Although not necessary for our purpose, one can imagine a theoretical heavy quark limit with fixed hardness $`𝓠`$:

$$𝓠\mathbf{}𝒎_𝒃\mathbf{}\mathrm{𝟐}𝑬_{\mathrm{𝐜𝐮𝐭}}\mathbf{=}\text{fixed}\mathbf{,}\text{ }𝒎_𝒃\mathbf{}\mathbf{}\mathbf{.}$$
(5)
The fully integrated moments of $`𝑭\mathbf{(}𝒙\mathbf{)}`$ and therefore of the spectrum give then directly the underlying heavy quark parameters:

$$\mathbf{}_\mathrm{𝟎}^{\mathbf{}}𝒌𝑭\mathbf{(}\overline{𝚲}\mathbf{}𝒌\mathbf{)}𝐝𝒌\mathbf{=}\overline{𝚲}\mathbf{,}\text{ }\mathbf{}_\mathrm{𝟎}^{\mathbf{}}\mathbf{(}𝒌\mathbf{}\overline{𝚲}\mathbf{)}^\mathrm{𝟐}𝑭\mathbf{(}\overline{𝚲}\mathbf{}𝒌\mathbf{)}𝐝𝒌\mathbf{=}\frac{𝝁_𝝅^\mathrm{𝟐}}{\mathrm{𝟑}}\mathbf{,}\text{ etc.}$$
(6)
As mentioned above, in the standard practical-OPE–based formulae these relations remain the same for the moments with the cut as well (provided $`𝑬_{\mathrm{𝐜𝐮𝐭}}\mathbf{<}\frac{𝒎_𝒃}{\mathrm{𝟐}}`$ which is always assumed) – yet not in reality. Paralleling the routinely used way we therefore introduce

$$\stackrel{\mathbf{~}}{𝚲}\mathbf{(}𝑬_{\mathrm{𝐜𝐮𝐭}}\mathbf{)}\mathbf{=}\frac{\mathbf{}_{𝑬_{\mathrm{𝐜𝐮𝐭}}}\mathbf{\left(}𝑴_𝑩\mathbf{}\mathrm{𝟐}𝑬_𝜸\mathbf{\right)}\frac{𝐝𝚪}{𝐝𝑬_𝜸}𝐝𝑬_𝜸}{\mathbf{}_{𝑬_{\mathrm{𝐜𝐮𝐭}}}\frac{𝐝𝚪}{𝐝𝑬_𝜸}𝐝𝑬_𝜸}\mathbf{,}\text{ }\stackrel{\mathbf{~}}{𝝁}_𝝅^\mathrm{𝟐}\mathbf{(}𝑬_{\mathrm{𝐜𝐮𝐭}}\mathbf{)}\mathbf{=}\mathrm{𝟑}\mathbf{[}\frac{\mathbf{}_{𝑬_{\mathrm{𝐜𝐮𝐭}}}\mathbf{\left(}𝑴_𝑩\mathbf{}\mathrm{𝟐}𝑬_𝜸\mathbf{\right)}^\mathrm{𝟐}\frac{𝐝𝚪}{𝐝𝑬_𝜸}𝐝𝑬_𝜸}{\mathbf{}_{𝑬_{\mathrm{𝐜𝐮𝐭}}}\frac{𝐝𝚪}{𝐝𝑬_𝜸}𝐝𝑬_𝜸}\mathbf{}\stackrel{\mathbf{~}}{𝚲}^\mathrm{𝟐}\mathbf{(}𝑬_{\mathrm{𝐜𝐮𝐭}}\mathbf{)}\mathbf{]}\mathbf{.}$$
(7)
Clearly the apparent value $`\stackrel{\mathbf{~}}{𝚲}\mathbf{(}𝑬_{\mathrm{𝐜𝐮𝐭}}\mathbf{)}`$ is always below the actual $`\overline{𝚲}`$. This is illustrated by Fig. 1a where we plot the cut-related ‘bias’ – the difference between $`\overline{𝚲}`$ and $`\stackrel{\mathbf{~}}{𝚲}\mathbf{(}𝑬_{\mathrm{𝐜𝐮𝐭}}\mathbf{)}`$ as a function of energy $`𝑬_{\mathrm{𝐜𝐮𝐭}}`$. It turns out to be quite significant.

The naive extraction of the kinetic expectation value through the variance of the truncated distribution undercounts it even more dramatically as Fig. 1b illustrates, since higher moments are more sensitive to the tail of the distribution.

To gauge the sensitivity to the choice of the heavy quark distribution function we follow Ref.  evaluating the effect for two ansätze – one exponential in $`𝒌_\mathbf{+}`$, and the other in $`𝒌_\mathbf{+}^\mathrm{𝟐}`$ yielding an even faster decreasing tail:

$$𝑭_\mathrm{𝟏}\mathbf{(}𝒌_\mathbf{+}\mathbf{)}\mathbf{=}𝑵_\mathrm{𝟏}\mathbf{(}\overline{𝚲}\mathbf{}𝒌_\mathbf{+}\mathbf{)}^𝜶𝒆^{𝒄𝒌_\mathbf{+}}𝜽\mathbf{(}\overline{𝚲}\mathbf{}𝒌_\mathbf{+}\mathbf{)}\mathbf{,}\text{ }𝑭_\mathrm{𝟐}\mathbf{(}𝒌_\mathbf{+}\mathbf{)}\mathbf{=}𝑵_\mathrm{𝟐}\mathbf{(}\overline{𝚲}\mathbf{}𝒌_\mathbf{+}\mathbf{)}^𝜷𝒆^{\mathbf{}𝒅\mathbf{(}\overline{𝚲}\mathbf{}𝒌_\mathbf{+}\mathbf{)}^\mathrm{𝟐}}𝜽\mathbf{(}\overline{𝚲}\mathbf{}𝒌_\mathbf{+}\mathbf{)}\mathbf{;}$$
(8)
the parameters are adjusted in such a way as to yield the same actual $`𝒎_𝒃`$ ($`\overline{𝚲}`$) and $`𝝁_𝝅^\mathrm{𝟐}`$. (For $`𝒎_𝒃\mathbf{=}\mathbf{4.6}\text{GeV}`$ and $`𝝁_𝝅^\mathrm{𝟐}\mathbf{=}\mathbf{0.43}\text{GeV}^\mathrm{𝟐}`$ we have $`𝜶\mathbf{=}\mathbf{2.22}`$ and $`𝜷\mathbf{=}\mathbf{0.773}`$.) Curiously, the ‘deficit’ $`\overline{𝚲}\mathbf{}\stackrel{\mathbf{~}}{𝚲}\mathbf{(}𝑬_{\mathrm{𝐜𝐮𝐭}}\mathbf{)}`$ practically does not depend on the choice at $`𝑬_{\mathrm{𝐜𝐮𝐭}}`$ around $`\mathrm{𝟐}\text{GeV}`$, and even the deficit in $`𝝁_𝝅^\mathrm{𝟐}`$ is reasonably stable.

Why are the effects of the cut so significant? They are exponential in the inverse hadronic scale $`𝝁_{\mathrm{𝐡𝐚𝐝𝐫}}`$, but are governed by the hardness $`𝓠\mathbf{}𝒎_𝒃\mathbf{}\mathrm{𝟐}𝑬_{\mathrm{𝐜𝐮𝐭}}`$ rather than by $`𝒎_𝒃`$:

$$\overline{𝚲}\mathbf{}\stackrel{\mathbf{~}}{𝚲}\mathbf{(}𝑬_{\mathrm{𝐜𝐮𝐭}}\mathbf{)}\mathbf{}𝝁_{\mathrm{𝐡𝐚𝐝𝐫}}𝒆^{\mathbf{}\frac{𝓠}{𝝁_{\mathrm{𝐡𝐚𝐝𝐫}}}}\mathbf{,}\text{ }𝝁_𝝅^\mathrm{𝟐}\mathbf{}\stackrel{\mathbf{~}}{𝝁}_𝝅^\mathrm{𝟐}\mathbf{}𝝁_{\mathrm{𝐡𝐚𝐝𝐫}}^\mathrm{𝟐}𝒆^{\mathbf{}\frac{𝓠}{𝝁_{\mathrm{𝐡𝐚𝐝𝐫}}}}$$
(9)
(the exponent may have a form of a power of $`𝓠\mathbf{/}𝝁_{\mathrm{𝐡𝐚𝐝𝐫}}`$). Even at $`𝒎_𝒃\mathbf{}\mathbf{}`$ these effects survive unless $`𝓠`$ is made large as well!

As explained in Ref.  the bias terms (9) are associated from a theoretical viewpoint with the limited (in fact, zero) convergence radius of the OPE. This becomes practically relevant due to the presence of a subseries in powers of $`\mathrm{𝟏}\mathbf{/}𝓠`$ rather than $`\mathrm{𝟏}\mathbf{/}𝒎_𝒃`$. The limitations on convergence appear due to a factorial growth of the power coefficients, a rather universal property of the OPE. In this respect one may associate this effect with quark-hadron duality . Yet it has no features peculiar to local quark hadron duality violation intrinsic to inclusive decay widths in actual Minkowski world. (For a discussion of the notorious subtleties in the notion of quark-hadron duality, see reviews ). For instance, these effects are truly exponential and do not oscillate.111Peculiarities of real local duality would appear here only at the next-to-leading order in $`1/m_b`$ and, therefore are not of much interest.

The validity of the routinely applied expressions for the moments with cuts for actual decays is additionally complicated by perturbative corrections. Incorporated into fits are naive sums of pure perturbative and pure nonperturbative terms:

$$𝑴_𝒏^{\mathrm{𝐧𝐩}}\mathbf{}𝑴_𝒏^{\mathrm{𝐧𝐩}}\mathbf{+}𝑴_𝒏^{\mathrm{𝐩𝐞𝐫𝐭}}\mathbf{(}𝜶_𝒔\mathbf{,}𝒎_𝒃\mathbf{,}𝑬_{\mathrm{𝐜𝐮𝐭}}\mathbf{)}$$
(10)
where nonperturbative corrections to the moments $`𝑴_𝒏`$ still do not depend on $`𝑬_{\mathrm{𝐜𝐮𝐭}}`$. This is not true in general, but would hold if the actual spectrum were exactly a convolution of the perturbative and nonperturbative spectra,

$$\frac{𝐝𝚪}{𝐝𝑬_𝜸}\mathbf{=}\mathbf{}𝐝𝒌\frac{𝐝𝚪^{\mathrm{𝐩𝐞𝐫𝐭}}\mathbf{(}𝑬_𝜸\mathbf{}𝒌\mathbf{)}}{𝐝𝑬}\frac{𝐝𝚪^{\mathrm{𝐧𝐩}}\mathbf{(}\frac{𝒎_𝒃}{\mathrm{𝟐}}\mathbf{+}𝒌\mathbf{)}}{𝐝𝑬}$$
(11)
provided no cut is introduced (or its effect on the pure nonperturbative distribution is negligible). We hasten to add, though that the effects we consider are unrelated to this complication and rather represent an independent phenomenon – they are significant even in the complete absence of perturbative corrections.

Since the perturbative effects can potentially modify the effect, we have evaluated the cut-induced deficit in $`\stackrel{\mathbf{~}}{𝚲}`$ and $`\stackrel{\mathbf{~}}{𝝁}_𝝅^\mathrm{𝟐}`$ including short-distance corrections. Namely, we considered $`\stackrel{\mathbf{~}}{𝚲}`$, $`\stackrel{\mathbf{~}}{𝝁}_𝝅^\mathrm{𝟐}`$ in Eqs. (7) for the complete spectrum obtained by the convolution (11) of the perturbative and primordial (nonperturbative) ones, and compared them to the naive sum Eq. (10) which would indeed hold for a sufficiently low cut. Including the perturbative spectrum as detailed in Ref.  we found no appreciable change in $`\overline{𝚲}\mathbf{}\stackrel{\mathbf{~}}{𝚲}`$ or $`𝝁_𝝅^\mathrm{𝟐}\mathbf{}\stackrel{\mathbf{~}}{𝝁}_𝝅^\mathrm{𝟐}`$ at realistic cuts (a small increase emerged only at $`𝑬_{\mathrm{𝐜𝐮𝐭}}\mathbf{}\mathbf{<}\mathrm{𝟏}\text{GeV}`$).

Based on these results we conclude:
$`\mathbf{}`$ the value of $`𝒎_𝒃`$ as routinely extracted from the $`𝒃\mathbf{}𝒔\mathbf{+}𝜸`$ spectrum is to be decreased by an amount of order $`\mathrm{𝟕𝟎}\text{MeV}`$;
$`\mathbf{}`$ relative corrections to $`𝝁_𝝅^\mathrm{𝟐}`$ are even more significant and can naturally constitute a shift of $`\mathbf{\hspace{0.17em}0.2}\text{GeV}^\mathrm{𝟐}`$. This arises on top of other potential effects.

Practical implications Accepting the above shifts at face value and using the rather arbitrary choice $`𝒎_𝒃\mathbf{=}\mathbf{4.595}\text{GeV}`$, $`𝒎_𝒄\mathbf{=}\mathbf{1.15}\text{GeV}`$, $`𝝁_𝝅^\mathrm{𝟐}\mathbf{=}\mathbf{0.45}\text{GeV}^\mathrm{𝟐}`$, $`\stackrel{\mathbf{~}}{𝝆}_𝑫^\mathrm{𝟑}\mathbf{=}\mathbf{0.06}\text{GeV}^\mathrm{𝟑}`$ and $`𝝆_{𝑳𝑺}^\mathrm{𝟑}\mathbf{=}\mathbf{}\mathbf{0.15}\text{GeV}^\mathrm{𝟑}`$ adjusted to accommodate $`\mathbf{}𝑴_𝑿^\mathrm{𝟐}\mathbf{}_{𝑬_{\mathbf{}}\mathbf{>}\mathrm{𝟏}\mathrm{𝐆𝐞𝐕}}`$, we obtain222The two values for the second $`E_\gamma `$-moment correspond to the two ansätze; they are obtained discarding higher-order power corrections to the light-cone distribution function.

$`\text{}\mathbf{}𝑴_𝑿^\mathrm{𝟐}\mathbf{}`$ $`\text{ }\mathbf{}`$ $`\mathbf{4.434}\text{GeV}^\mathrm{𝟐}\text{ }\mathbf{[}\text{cf. }\mathbf{(}\mathbf{4.538}\mathbf{\pm }\mathbf{0.093}\mathbf{)}\text{GeV}\text{ (New DELPHI)}\mathbf{]}`$
$`\text{}\mathbf{}𝑴_𝑿^\mathrm{𝟐}\mathbf{}_{𝑬_{\mathbf{}}\mathbf{>}\mathbf{1.5}\mathrm{𝐆𝐞𝐕}}`$ $`\text{ }\mathbf{}`$ $`\mathbf{4.177}\text{GeV}^\mathrm{𝟐}\text{ }\mathbf{[}\text{cf. }\mathbf{4.180}\text{GeV}^\mathrm{𝟐}\text{ (BaBar), }\mathbf{4.189}\text{GeV}^\mathrm{𝟐}\text{ (CLEO)}\mathbf{]}`$
$`\text{}\mathbf{}𝑬_{\mathbf{}}\mathbf{}`$ $`\text{ }\mathbf{}`$ $`\mathbf{1.389}\text{GeV}\text{ }\mathbf{[}\text{cf. }\mathbf{(}\mathbf{1.383}\mathbf{\pm }\mathbf{0.015}\mathbf{)}\text{GeV}\text{ (DELPHI)}\mathbf{]}`$
$`\text{}\mathbf{}𝑬_𝜸\mathbf{}_{𝑬_𝜸\mathbf{>}\mathrm{𝟐}\mathrm{𝐆𝐞𝐕}}`$ $`\text{ }\mathbf{}`$ $`\mathbf{2.329}\text{GeV}\text{ }\mathbf{[}\text{cf. }\mathbf{(}\mathbf{2.346}\mathbf{\pm }\mathbf{0.034}\mathbf{)}\text{GeV}\text{ (CLEO)}\mathbf{]}`$
$`\text{}\mathbf{}𝑬_𝜸^\mathrm{𝟐}\mathbf{}\overline{𝑬}_𝜸^\mathrm{𝟐}\mathbf{}_{𝑬_𝜸\mathbf{>}\mathrm{𝟐}\mathrm{𝐆𝐞𝐕}}`$ $`\text{ }\mathbf{}`$ $`\text{ }_{\mathbf{0.0233}}^{\mathbf{0.0202}}\text{GeV}^\mathrm{𝟐}\text{ }\mathbf{[}\text{cf. }\mathbf{(}\mathbf{0.0226}\mathbf{\pm }\mathbf{0.0066}\mathbf{\pm }\mathbf{0.0020}\mathbf{)}\text{GeV}^\mathrm{𝟐}\text{ (CLEO)}\mathbf{]}`$ (12)
(experimental data are from Refs. ). $`𝑬_{\mathrm{𝐜𝐮𝐭}}^{\mathbf{}}`$-dependence of $`\mathbf{}𝑴_𝑿^\mathrm{𝟐}\mathbf{}`$ is also reproduced, see Fig. 2. The counterpart of the above exponential cut-related effects for the semileptonic decays has not been incorporated here, however.

It should be noted that the cut-induced shifts $`\stackrel{\mathbf{~}}{𝒎}_𝒃\mathbf{}𝒎_𝒃`$, $`\stackrel{\mathbf{~}}{𝝁}_𝝅^\mathrm{𝟐}\mathbf{}𝝁_𝝅^\mathrm{𝟐}`$ are not unambiguously determined in terms of a few known HQ parameters, but rather depend on the actual shape of the heavy quark distribution function; in particular it is driven by the asymptotics of certain high-dimension expectation values. The evaluation presented above relies on the most natural assumptions about the function. Strictly speaking, they can vary for more contrived ansätze. Therefore, the estimates in Eqs. (1) can be conservatively viewed as the minimal possible theoretical inaccuracy of the usual naive evaluations. The safest way to tackle possible ambiguity is to resort to more inclusive moments without too severe cuts.

Cuts in semileptonic moments. As argued in Ref.  similar cut-related ‘exponential’ biases missed in the naive OPE applications affect the truncated moments in the semileptonic decays as well. Their description, even simplified is less transparent and would be more involved, though. In particular, the light-cone distribution functions for $`𝒃\mathbf{}𝒔\mathbf{+}𝜸`$ is replaced by a different function, the form of which actually is not universal. Yet the qualitative trend is expected the same – it should smoothly interpolate the case of literal OPE expression at low $`𝑬_{\mathbf{}}`$ cuts and the values at a high cut dictated simply by the actual hadrons kinematics. This would replace step- or $`𝜹`$-like behavior in the formal OPE expressions.

Numerical aspects are less certain. Keeping in mind that for $`𝑩\mathbf{}𝑿_𝒄\mathbf{}𝝂`$ the effective hardness

$$𝓠_{\mathrm{𝐬𝐥}}\mathbf{}𝒎_𝒃\mathbf{}𝑬_{\mathrm{𝐜𝐮𝐭}}\mathbf{}\sqrt{𝑬_{\mathrm{𝐜𝐮𝐭}}^\mathrm{𝟐}\mathbf{+}𝒎_𝒄^\mathrm{𝟐}}$$
(13)
at $`𝑬_{\mathrm{𝐜𝐮𝐭}}\mathbf{=}\mathbf{1.5}\text{GeV}`$ is about $`\mathbf{1.25}\text{GeV}`$ , nearly the same as for $`𝑩\mathbf{}𝑿_𝒔\mathbf{+}𝜸`$ with $`𝑬_𝜸\text{ }\mathbf{}\mathbf{>}\mathrm{𝟐}\text{GeV}`$, we may expect quite significant effects. To get a rough idea of the possible magnitude of the bias we can use a simplified rule of thumb – assume that $`𝒎_𝒃`$ in the semileptonic decay can be just replaced by an effective larger value $`\stackrel{\mathbf{~}}{𝒎}_𝒃`$ apparent in $`𝒃\mathbf{}𝒔\mathbf{+}𝜸`$ at the commensurate cut yielding the same hardness. In other words, we mimic the effect of decreasing hardness by an additional effective nonperturbative running of the heavy quark mass at low scales. As illustrated above, this is to increase $`𝒎_𝒃`$ (or, equivalently decrease $`\overline{𝚲}`$) by about $`\mathrm{𝟕𝟎}\text{MeV}`$, a significant change.

Yet the semileptonic decay characteristics strongly depend on both $`𝒎_𝒃`$ and $`𝒎_𝒄`$. To stay on the conservative side we assume the apparent shift in $`𝒎_𝒄`$ as high as in $`𝒎_𝒃`$ (which would reflect just the heavy flavor symmetry):

$$𝒎_𝒄\mathbf{}\stackrel{\mathbf{~}}{𝒎}_𝒄\mathbf{}𝒎_𝒄\mathbf{+}\mathbf{(}\stackrel{\mathbf{~}}{𝒎}_𝒃\mathbf{}𝒎_𝒃\mathbf{)}\mathbf{.}$$
(14)
In actuality the corrections to $`𝒎_𝒄`$ are typically somewhat smaller due to $`\mathrm{𝟏}\mathbf{/}𝒎_𝒄^𝒌`$ terms.

As discussed elsewhere , for actual $`𝑩`$ decays both lepton moments and $`\mathbf{}𝑴_𝑿^\mathrm{𝟐}\mathbf{}`$ depend on more or less the same combination of masses $`𝒎_𝒃\mathbf{}\mathbf{0.65}𝒎_𝒄`$. This means that the literal ansatz (14) would suppress the effect by a factor of $`\mathrm{𝟑}`$ to $`\mathrm{𝟒}`$, yet $`\mathrm{𝟏}\mathbf{/}𝒎_𝒄`$ corrections may be thought to softening this suppression. We then expect the exponential terms in semileptonic decays with $`𝑬_{\mathrm{𝐜𝐮𝐭}}\mathbf{}\mathbf{1.5}\text{GeV}`$ to introduce effects on the same scale as shifting $`𝒎_𝒃`$ upward by up to $`\mathrm{𝟐𝟓}`$ to $`\mathrm{𝟑𝟎}\text{MeV}`$ (assuming fixed $`𝒎_𝒄`$ and other heavy quark parameters). This rule of thumb is useful to get an idea of the ultimate theoretical accuracy one can count on.

For example, the CLEO’s cut moment $`𝑹_\mathrm{𝟏}\mathbf{=}\mathbf{}𝑬_{\mathbf{}}\mathbf{}𝑬_{\mathbf{}}\mathbf{>}\mathbf{1.5}\text{GeV }`$ is approximately given by

$$𝑹_\mathrm{𝟏}\mathbf{=}\mathbf{1.776}\text{GeV}\mathbf{+}\mathbf{0.27}\mathbf{(}𝒎_𝒃\mathbf{}\mathbf{4.595}\text{GeV}\mathbf{)}\mathbf{}\mathbf{0.17}\mathbf{(}𝒎_𝒄\mathbf{}\mathbf{1.15}\text{GeV}\mathbf{)}\text{ }\text{at }\mathbf{|}𝑽_{𝒖𝒃}\mathbf{/}𝑽_{𝒄𝒃}\mathbf{|}\mathbf{=}\mathbf{0.08}$$
(15)
(the above mentioned values of the nonperturbative parameters are assumed). An increase in $`𝒎_𝒃`$ by only $`\mathrm{𝟐𝟎}\text{MeV}`$ would then change

$$𝑹_\mathrm{𝟏}\mathbf{}𝑹_\mathrm{𝟏}\mathbf{+}\mathbf{0.0055}\text{GeV}$$
(16)
and would perfectly fit the central CLEO’s value $`\mathbf{1.7810}\text{GeV}`$. It is worth noting that Eqs. (15)-(16) make it explicit that the imposed cut on $`𝑬_{\mathbf{}}`$ degrades the theoretical calculability of $`𝑹_\mathrm{𝟏}`$ far beyond its experimental error bars, the fact repeatedly emphasized over the last year. Unfortunately, this was not reflected in the fits of parameters which placed much weight on the values of $`𝑹_\mathrm{𝟎}`$$`𝑹_\mathrm{𝟐}`$ just owing to their small experimental uncertainties, whilst paying less attention to actual theoretical errors.333N.U. acknowledges that a similar in spirit criticism of the theory error treatment in  was expressed by D. Hitlin at the BaBar Workshop, SLAC December 2002.

Similar reservations apply to the quality of theory for CLEO’s $`𝑹_\mathrm{𝟐}`$ representing the second moment with the cut, with the only difference that the effective hardness only deteriorates for higher moments.

The ratio $`𝑹_\mathrm{𝟎}`$ is the normalized decay rate with the cut on $`𝑬_{\mathbf{}}`$ as high as $`\mathbf{1.7}\text{GeV}`$, and for it hardness $`𝓠`$ is below $`\mathrm{𝟏}\text{GeV}`$. A precision – beyond just semiquantitative – treatment of nonperturbative effects is then questionable, and far more significant corrections should be allowed for.

We note that there is a good agreement of most data referring to sufficiently ‘hard’ decay distributions with the theory based on the OPE in QCD, in the ‘robust’ OPE approach. The latter was put forward  to get rid of vulnerable and unnecessary assumptions of the usually employed fits of the data. This consistency likewise refers to the absolute values of the heavy quark parameters necessary to accommodate the data. Theory anticipates, however that the expansion becomes deceptive with increase of the experimental cuts. Here we have addressed the most obvious effects, those from the variety of ‘exponential’ terms in the effective hardness. While presently not amenable to precise theoretical treatment, they can be estimated using quite natural assumptions and are found to be far too significant for $`𝑬_{\mathrm{𝐜𝐮𝐭}}^{\mathbf{}}\text{ }\mathbf{}\mathbf{>}\mathbf{1.5}\text{GeV}`$ and $`𝑬_{\mathrm{𝐜𝐮𝐭}}^𝜸\text{ }\mathbf{}\mathbf{>}\mathrm{𝟐}\text{GeV}`$ often employed in experiment. Taking these corrections at face value and incorporating in our predictions, we get a good, more than qualitative agreement with “less short-distance” inclusive decays as well.

To summarize, the cuts essentially decreasing the hardness in $`𝑩`$ decays introduce terms some of which are exponentially suppressed though only in the effective hardness, but not suppressed by powers of the heavy quark mass. They significantly change the extracted values of the heavy quark parameters leading to an apparent suppression of the magnitude of the nonperturbative ones in $`𝑩`$ mesons. Accounting for such effects appears necessary in $`𝑩\mathbf{}𝑿_𝒔\mathbf{+}𝜸`$ decays unless the cut on photon energy is pushed well below present $`\mathrm{𝟐}\text{GeV}`$. This brings the practical predictions for various inclusive characteristics into a good agreement. Moreover, it seems that the second moment of the photon spectrum discarded for the fits so far, does impose informative constraints provided these effects are properly incorporated.

Acknowledgments: N.U. is indebted to P. Gambino for joint work heavily used in the present analysis. He is also thankful to O. Buchmueller and U. Langenegger for many discussions which initiated this study, and for useful comments. This work was supported in part by the NSF under grant number PHY-0087419.


text:
# Stellar abundance patterns. What is the possible level of completeness today ?

## 1. Introduction

If we will try to overview the results of determinations of stellar abundances, we will find that a lot of elements with Z$`>`$30 are omitted. Elements heavier than iron group are synthesized primarily by means of neutron capture processes. According to the work by Burbidge et al. (1957, B2FH), in order to reconstruct the Galactic evolutionary history of heavy elements , one has to consider two major mechanisms of neutron addition: the s-process and r-process. Small amount of nuclei can be created in p-process. New process of heavy elements creation was proposed by Woosley & Hoffman (1992) - this is $`\alpha `$-process for heavy elements. A review of the developments of the theory of elements creation was made by Wallerstein et al. (1997). This paper was devoted to the 40th anniversary of B2FH work. The theory of element creation needs a detailed abundance pattern for comparison of theoretical predictions with observations. New data on stellar abundance patterns for stars of different types can significantly influence the theory.

The full isotopic pattern is available only for meteoritic matter. The precision of abundance determinations in stars is not so high. The best abundance samples for the stars are:

Sun – 73 elements (Grevesse & Sauval, 1998)

Procyon A – 55 elements (Yushchenko & Gopka, 1996a,b),

Przybylski’s star – 54 elements (Cowley et al., 2000),

$`\chi `$ Lupi – 51 elements (Leckrone et al., 1999),

Gopka (2000) investigated the abundances of several heavy elements if the atmosphere of Sirius A. Now we are preparing a review of chemical composition of Sirius A. The abundance pattern of this star consists of 50 elements We can see a significant progress. In 1988 Reinolds et al. pointed, that the abundance pattern for Canopus (38 elements) is the third after the Sun and Przybylski’s star (Wegner & Petford, 1974 – 51 elements).

We realize, that this review is not full, but it is sufficient to make some conclusions. First of all - detailed abundance determination is a very long process. In many cases it cost a decades of efforts of many scientists. The majority of above cited papers are the review papers or the papers which are the final paper in a long series of articles on abundances in this star. What are the necessary conditions for future progress ? The base of the improvement are:

1) the increasing of level/noise ratio and spectral resolution stellar spectroscopy, investigation of the chemical abundance using UV spectra;

2) the new atomic data;

3) the using of the spectrum synthesis method not only for limited number of lines in the spectrum, but for majority of lines, taking in account hyperfine and isotopic splitting, magnetic fields, spots, detailed analysis of spectral binaries, individual atmosphere models, etc.

Let us show a short overview for these three items. We will show only some aspects which, we hope, will be interesting for the community of this conference. It is difficult to point poor better observations or poor better atomic data or software influence on final results. Usually it is the combination of all three items.

## 2. High quality observations

The majority of cited papers, which give us the most complete stellar abundance samples, deals with the spectra obtained at the best telescopes and spectrographs. But it should be noted, that we have 2.5 meter telescope since 1917, 5 meter - since 1948, 6 meter - since 1975, 10 meters - since 1991. But while this telescopes were the single in it’s class, no significant progress were available. Only during last two decades, when the number of telescopes over 2 meters became near 50, we can see a lot of new results.

During the last decade observations with coude-echelle spectrographs were started at the 1 meter telescope of SAO RAS and 2 meter telescope of ICAMER. The last telescope is located at 3100 meters elevation. Now it is possible to obtain spectra with resolutions 45000, 80000, 125000, 190000. Now we are testing the new spectrograph which permit us to obtain the spectra of bright stars with the resolution comparable with solar spectral atlases.

In the last section of this paper we review the preliminary results of one of the first detailed investigations of stellar abundances with this spectrograph. The mild barium star $`\zeta `$ Cyg was observed in 2000 with relatively low resolution - near 80000. This investigation will be the base for future observation of this and other barium stars with highest spectral resolution in visual and near UV wavelengths regions.

## 3. New atomic data – line identification in Przybylski’s star

The more detailed review of used observations and methodics one can find near this poster in the poster of Shavrina et al. (2002). Here we will point only important items:

1) wavelength coverage 6123-6175 ÅÅ and 6676-6732 ÅÅ – 108 Å;

2) spectral resolution R=100000, signal to noise ratio S/N$`>`$100;

3) new lanthanides lines from DREAM database (Biemont et al., 2002);

4) identification of lines and calculation of abundances were made using spectrum synthesis method. We used Kurucz (1995) SYNTHE program for calculation synthetic spectra and Yushchenko (1998) URAN program for line identification and spectrum synthesis in automatic mode. Abundances for all identified lines were found with spectrum synthesis method.

In tables 1,2 we give the summary of results. In table 1 one can find the total number of lines of different elements and ions in DREAM database. The numbers in the last column of this table \- the part of these lines, which were present in our previous line list. The total number of lanthanides and thorium lines in DREAM database is 56150.

We tried to use these lines to find the abundances in Przybylski’s star. We used 108 Å wavelength interval in red part of the spectrum. In table 2 we show our mean abundances, errors and number of lines from 108 Å interval in red region of the spectrum and the corresponding values from 2693 Å interval in blue-red parts of the spectrum according to Cowley et al. (2000).

The brief inspection of the abundances show, that the results are similar in the range of uncertainties. In the last column of this table we give the ratio of lines, identified in the work of Cowley et al (2000) and in this work.

From the ratio of wavelengths coverage in two studies we can expect that the ratio of numbers of identified lines will be near 25 or more, if we take in account, that the majority of lanthanides lines are located in the blue part of the spectrum. But only one value in the last column of the table exceeds 10.

It means that the lanthanides lines from DREAM database give us a powerful instrument for examination of the atmospheres of peculiar stars.

It should be noted that not only lanthanides show the overabundance in the atmosphere of Przybylski’s star, but all heavy elements. Unidentified lines can be the lines of other heavy elements.

## 4. Chemical composition of $`\zeta `$ Cyg

$`\zeta `$ Cyg is one of the brightest middle barium stars. The spectrum was obtained at ICAMER 2 meter telescope with resolving power R=80000 in the wavelength range 3495-10000 Å with signal to noise ratio in visual and infrared region more than 100. As a first step we tried to find atmosphere parameters. We tested the parameters used by Zacs (1994) and Boyarchyk et al. (2001). We analyzed iron lines in the spectrum and found that following atmosphere parameters are valid:

Teff=5050 K, $`lgg`$=2.8, $`v_{micro}`$=1.45 km/s, $`v_{macro}`$=3-3.5 km/s,

We interpolated Kurucz (1995) atmosphere model with these parameters and calculated synthetic spectrum for all observed region. This spectrum was used for identification of spectral lines. Abundance calculations for all elements except iron were made with spectrum synthesis method. We used Kurucz (1995) SYNTHE program for calculation of synthetic spectra and Yushchenko (1998) URAN program for approximation of the observed program by calculated one in automatic mode.

We tried to obtain abundances by direct comparison with the solar spectrum. We found solar oscillator strengths for majority of the used lines. We used Liege Solar atlas (Delbouille et al., 1974), Holweger-Muller atmosphere model, microturbulent velocity 1.0 km/s, macroturbulent velocity 1.8 km/s. SYNTHE and URAN codes were used for approximation of observed solar spectrum by synthetic one.

In tables 3,4 we show the results of abundance determinations. The results of Gratton (1985), Zacs (1994), Boyarchuk et al.(2001) are displayed for comparison. The chemical elements with lines, that have no counterparts in the solar spectrum are marked by asterisk.

We found the abundances of 48 elements in the atmosphere of $`\zeta `$ Cyg. The abundances of Li, N, Ba are known from previous investigations of this star. The total abundance sample consist of 51 elements. Detailed analysis of these results will be made later.

## 5. Conclusion

We show that the abundances of 50-55 elements are the best results for stellar abundance samples in sharp-lined stars. We can point some ways how to increase these number for several types of the stars. First of all F stars of main sequence. The detailed chemical composition of Procyon is a summary of many papers, and significant part of this papers were used Griffin’s (1979) atlas of this star. New atlas of the star observed with higher signal to noise ratio and spectral resolution will help to find several new elements.

The other way is the observations of UV spectrum of Procyon or Procyon type stars. We have no determinations of chemical composition of Procyon from UV spectral data.

For barium stars additional abundances can be obtained from observations in the near UV wavelengths region - 3100-3500 Å.

And let we will not forget about the Sun. There are a dozen of chemical elements with unknown abundances in the solar atmosphere. For example Gopka et al. (2001) found the abundance of arsenic in the Sun. The spectral atlases of solar spots can be used for abundance determinations, if the good model of the spot will be made.

It should be noted that all these ways will lead us to very crowded spectral regions. Usual model atmosphere method is not available in this case. It is necessary to use spectrum synthesis and automatic spectrum synthesis. The programs for automatic spectrum synthesis were described by Cowley (1995), Valenty and Piskunov (1995), Tsymbal & Cowley (2000), Erspamer & North (2002), Yushchenko (1998). These programs are very different, but all of them can help us to obtain detailed abundance samples for different type stars.

Here before we have showed the significance of new atomic data for lanthanides. But it should be noted that Corliss & Bozman (1961) oscillator strengths are still in use for several elements. We have zero or very limited information about the abundances of elements with Z=50-55, 33-36. One of the peaks of solar system abundances is near Z=50-55. And we have zero information about this peak in the Sun and in other stars. One of the first attempts to find the abundance of tellurium (Z=52) (Yushchenko & Gopka, 1996a,b) lead us to possible overabundance of this element in Procyon.

Maybe a decade later, when 50-70 elements will be a common number of elements in every paper on chemical abundance in sharp-lined stars, we will be able to resolve a lot of modern problems in this field and to ask new questions.

### Acknowledgments.

This work was supported by the grant of Post-Doc Program, Chonbuk National University, (South) Korea (2002).

## References

Biemont J., Palmeri P., Quinet P. Database of rare earths at Mons University http://www.umh.ac.be/ãstro/dream.html

Boyarchuk A., Antipova L.,Boyarchuk M, Savanov I. 2001, Astronomy reports, 45, 301

Burbidge E.M., Burbidge G.R., Fawler W.A., Hoyle F. 1957, Reviews Modern Physics, 29, 547 (B2FH)

Corliss C.H., Bozman W.R. 1962, Experimental Transition probabilities for spectral lines of seventy elements. NBS Monograph 53, Moscow, 1968 ( Russian edition )

Cowley C. 1996, ASP Conf. Ser., 108, 170

Cowley C., Ryabchikova T., Kupka F., Bord D., Mathys G., Bidelman W. 2000, MNRAS, 317, 299

Delbouille L., Roland G., Neven L. 1973, Photometric Atlas of the Solar Spectrum from $`\lambda `$ 3000 to $`\lambda `$ 10000. Liege: Institut d’Astrophysique de l’Universite’ de Liege.. (Digitized version from $`\lambda `$ 3601 Å to $`\lambda `$ 8400 Å)

Erspamer D., North P. 2002, A&A, 383, 227

Gopka V.F. 2000, Kinematics and physics of Selectial Bodies. Supplement Series, 3, 364

Gopka V.F., Yushchenko A.V., Shavrina A.V., Andrievsky, S.M., Chernysheva L.V. 2001, Kinematics and Physics of selectial bodies, 17, 37

Gratton R.G. 1985, A&A, 148, 105

Grevesse N., Sauval A.J. 1998, Space Science Reviews, 85, 161

Griffin R., Griffin R. 1979, A photometric atlas of the spectrum of Procyon $`\lambda `$$`\lambda `$ 3140-7475 Å, Cambridge Philosophical Society

Holweger H., Muller E.A. 1974, Solar Phys. 39, 19

Kurucz R.L. 1995, ASP Conf. Ser. 81, 583

Leckrone D.S., Proffitt C.R.. Wahlgren G.M.. Johansson S.G., Brage T. 1999, AJ, 117, 1454

Reynolds S.E., Hearnshaw J.B., Cottrell P.L. 1988, MNRAS, 235, 1423

Shavrina A., Polosukhina N., Pavlenko Ya., Yushchenko A., Quinet P., Hack J., North P., Gopka V., Zverko J., Ziznovsky J. 2002, This conference

Tsymbal V., Cowley C. 2000, private communication

Valenty J.A., Piskunov N.E. 1996, ASP Conf. Ser. 108, 175

Wallerstein G., Iben I., Parker P., Boesgaard A.M., Hale G.M., Champagne A.E., Barnes C.A., Kappeler F., Smith V.V., Hoffman R.D., Timmes F.X., Sneden C., Boyd R.N., Meyer B.S., Lambert D.L. 1997, Reviews Modern Physics, 69, 995

Wegner G., Petford A.D. 1974, MNRAS, 168, 557

Woosley S.E., Hoofman R.D. 1992, ApJ, 395, 202

Yushchenko A.V. 1998, Proc. of the 20th conf. of variable star research, Brno, Czech Republic, November 5-9, 1997, 201, Ed. Dusek J.

Yushchenko A.V., Gopka V.F. 1996a, Astronomy Letters, 22, 412

Yushchenko A.V., Gopka V.F. 1996b, Odessa Astronomical Publications, 9, 86

Zacs L. 1994, A&A, 283, 937


text:
# Interstitial Mn in (Ga,Mn)As: Binding energy and exchange coupling

## I Introduction

Diluted magnetic semiconductors, such as (Ga,Mn)As, are important materials combining ferromagnetic behavior with a sensitivity to doping characteristic for semiconductors. This leads to their interesting physical properties and makes possible their applications e.g. in spin-electronics Dietl00 ; Ohno03 . Although (Ga,Mn)As has been extensively studied in the last years, some aspects of the incorporation of Mn into the crystal lattice still remain unclear. It was assumed that in the well defined samples Mn simply substitutes for the host cation. Only recently it was pointed out Masek01 ; Maca02  that the presence of interstitial Mn may explain some peculiar properties of (Ga,Mn)As such as the low doping efficiency of Mn acceptors Beschoten99 . The main reason for considering the Mn atoms on the interstitial positions was that they act as double donors and partly compensate the Mn acceptors in the substitutional positions.

At the same time, channelling Rutherford backscattering experiments proved that a large fraction of Mn atoms indeed occupies the interstitial positions (Mnint) in the as-grown samples Yu02 . A correlation between the removal of Mn interstitials and increase of the conductivity, the Curie temperature, and saturation magnetization has been found.

In addition to the self-compensation effect, the interstitial Mn atom also reduce the number of local moments that participate in the ferromagnetic state. This was explained by pairing of Mnint with the Mn atoms in the substitutional positions (MnGa) due to their Coulomb attraction Yu02 . At the bonding distance, the antiferromagnetic superexchange within the pair is assumed to outweigh the hole mediated ferromagnetic exchange. As a result, the moments of the paired Mn atoms have opposite directions and the pair as a whole has no magnetic moment.

The ab initio studies of the interstitial Mn showed the differences in the electronic structure of the Mn atom in the interstitial and substitutional positions Maca02 ; Sanvito02  and also a possible reaction path for the incorporation of Mn into the GaAs lattice Erwin02 . It was also found that the increase of the lattice constant of Ga1-xMnxAs with increasing content of Mn is partly due to the presence of Mn in the interstitial positions Masek03 .

Recently, Blinowski and Kacman Blinowski02 ; Blinowski03  investigated the coupling of the local moment on the interstitial Mn to the spin of the carriers in the valence band. Using a simplified tight-binding model they showed that the local moment on Mn atom in the tetrahedral position with four Ga neighbors, T(Ga4), is effectively decoupled from the spin of the holes in the valence band.

They claimed that the coupling is weak because, in addition to the reduced number of the holes, also the coupling constant $`J_{pd}`$ is small for the interstitial Mn in the T(Ga4) position. As a result, the local moment is not subject to the ferromagnetic coupling with the moments around it. This opens the way for the AF exchange to be important in the Mnint-MnGa pair, as anticipated in Ref. Yu02 .

There are, however, several open questions concerning the magnetic interactions of Mn in the interstitial positions. First of all, the spin-polarized band structures Maca02  did not show any indication of different values of $`J_{pd}`$ for substitutional and interstitial Mn, at least for the T(As4) positions. In addition, the Mnint-MnGa pair as a whole is expected Blinowski02 ; Blinowski03  to have only a small magnetic moment, but the exchange coupling with MnGa – being uncompensated by the contribution of Mnint – should strongly polarize the holes.

That is why we performed a more detailed study of the interstitial Mn and its spin interactions. We use the density-functional, full-potential linearized augmented plane wave calculations (FPLAPW WIEN ) to obtain the electronic structure of (Ga,Mn)As with Mn atoms in various crystallographic positions. The calculated total energies are used to compare different positions of the interstitial Mn and to estimate the strength of the MnGa-Mnint pair interactions. The splitting of the valence band for the majority- and minority-spin electrons is used to compare the corresponding values of the exchange parameters $`J_{pd}`$.

## II Mn in various interstitial positions

We compare the total energies for three interstitial positions of Mn in GaAs. There are two inequivalent tetrahedral positions in the zinc-blende structure of GaAs, T(As4) and T(Ga4). They are surrounded by four As and Ga atoms, respectively. The (unrelaxed) distances of these nearest neighbors are equal to the length $`d_1`$ of the Ga-As bond, i.e. to the Mn-As distance for the substitutional MnGa. The local arrangement around the interstitial Mn is, in contrast to MnGa, characterized by another six close neighbors at the distance $`d_2=\sqrt{9/8}d_11.155d_1`$. In the hexagonal position, the interstitial Mn has three Ga and three As atoms at the same distance $`d_{hex}=\sqrt{11/12}d_10.957d_1`$, and no other close neighbors.

We represent GaAs with Mnint by hexagonal Ga12As12Mn supercells. The c-axis coincides with the body diagonal of the conventional cubic unit cell. The symmetry of the cell does not change if we shift the Mn interstitial along the c-axis from T(As4) to hexagonal and T(Ga4) positions. This makes possible to perform all calculations under the same conditions.

The results are summarized in Table I. Surprisingly, the binding energy of the interstitial Mn does not depend much on its nearest neighbors. The interstitial Mn has minimum energy in the T(As4) position. However, the difference of the total energies obtained for Mn in T(Ga4) and T(As4) positions is of order of a few meV and can be neglected in practice. This means that, without intervention of other charged defects, the Mn interstitials can be found with an almost equal probability in either T(As4) or T(Ga4) position. The total energy corresponding to the hexagonal interstitial position of Mn is remarkably higher and represents a barrier $``$ 0.5 eV separating the tetrahedral positions. Calculations with similar results have been recently done by Boguslawski Edmonds03 .

The weak influence of the nearest neighbors on the interstitial Mn in T(As4) and T(Ga4) positions can be observed also in the densities of states in Fig. 1. The total DOS and the distribution of the Mn d-states are almost identical. Also the DOS of Ga12As12Mn with Mn in the hexagonal position is quite similar and differs mainly by the overlap of the valence band for the majority-spin electrons ($``$) with the conduction band for the minority-spin ($``$) electrons. This tendency to close the gap in the electron spectrum correlates well with the increase of the total energy.

The lattice relaxation around the Mn impurity is not very important for the substitutional Mn as shown by Mirbt et al. Mirbt02 . This is not the case of the interstitial Mn which tends to expand the crystal lattice Masek03 . The addition of Mn into either T(As4) or T(Ga4) positions results in a remarkable repulsion of the nearest and next-nearest neighbors. The positions of the more distant atoms in the supercell are changed much less. The energy gain due to the relaxation is in both cases approximately 20 meV and it does not change the ordering of the total energies for the T(As4) or T(Ga4) positions. The lattice relaxation around Mn in the hexagonal position has not been considered because of its minor importance.

In the case of the T(As4) position, the distance of the four nearest As neighbors was found to increase by $``$ 1.5 % from 2.45 Å to 2.47 Å. The distance of the six next-nearest Ga neighbors increases by $``$ 0.7 % from 2.835 Å to 2.88 Å.

The relaxation around Mn in the T(Ga4) position is larger. The nearest Ga neighbors are pushed to the distance 2.515 Å and the relaxed distance of the next-nearest neighbors increases by $``$ 0.5 %. The enhanced relaxation of the nearest neighbors in this case is due to the Coulomb repulsion between Mn and Ga atoms which are both positively charged.

In all cases, there are two electrons in the conduction band, i.e., the Mnint always acts as a double donor. The electrons in the conduction band are almost completely spin-polarized. They accumulate in the minority-spin conduction band, so that the total spin of the cell (i.e., the spin per Mn) is reduced to $`3/2`$ in accordance with our previous calculations Maca02 .

For comparison, we performed the same set of calculations also for hypothetical Ga9Zn3As12Mn crystals with Zn atoms substituted at the sites most distant from the interstitial Mn. The presence of Zn has only a little effect on the density of states, but the material is converted into the p-type with one hole in the valence band. In this case, the T(Ga4) turns to be the stable interstitial position of Mn. The energy difference between the T(As4) and T(Ga4) positions, approximately 60 meV, is high enough for a preferential occupation of the T(Ga4) position at typical growth conditions. The remarkable increase of the energy of the hexagonal position indicates that the mobility of the interstitial Mn may depend on the type and degree of the doping.

## III Binding energy and exchange coupling of a Mn – Mn pair

The hexagonal unit cell used in Sect. II is well suited also for the study of the Mnint-MnGa pair interactions. We consider a hypothetical Ga11MnAs12Mn crystal with one substitutional and one interstitial Mn in the unit cell. We consider three positions of the interstitial Mn shown in Fig. 2. The positions T(Ga3Mn) and T(As4), denoted (a) and (b) in Fig. 2, are located on the c-axis of our unit cell and correspond to the initial stage of dissociation of the MnGa-Mnint pair. On the other hand, the T(As4) position denoted (c) is very close to MnGa and is - together with T(Ga3Mn) - a candidate for the ground state of Mnint.

In the T(Ga3Mn) position, Mnint and MnGa are the nearest neighbors at the distance $`d_1`$. We performed the density-functional calculations with two initial conditions, corresponding to parallel and antiparallel alignment of their local moments. In both cases, the self-consistent procedure converges to a locally stable solution without changing the initial alignment of the local moments. The resulting total energies are given in Table II. The antiparallel alignment is energetically more favorable than the parallel alignment, in a good correspondence with the expectations Yu02 ; Blinowski03 . The coupling is strong enough so that the AF state of the Mnint-MnGa pair is stable with respect to the thermal fluctuations.

The Table II shows also the local spins on MnGa and Mnint atoms defined as integrals of the spin density over the corresponding atomic spheres. Although these quantities are not directly related to the size of the observable local moments, we can see that the local moments of MnGa and Mnint are comparable and that the total magnetic moment of the pair in the ground state is much smaller than the magnetic moment of a single Mn.

Assuming that the exchange coupling between more distant Mn atoms from different unit cells is much smaller than the exchange interaction $`J(d_1)`$ within the closest pair we can estimate its strength from the difference $`\mathrm{\Delta }E(d_1)=E_{}(d_1)E_{}(d_1)2J(d_1)S^2`$ of the total energies. Using saturated values $`S=5/2`$ for both local moments we obtain the lower estimate for $`J(d_1)`$, namely $`J(d_1)`$ -26 meV.

The coupling between Mnint and MnGa remains antiferromagnetic also if Mnint moves to any of the adjacent T(As4) positions (b) and (c) in Fig. 2. The energy difference between the parallel and antiparallel alignment of the local moments decreases with the increasing distance of the Mn atoms, as expected for the superexchange. It remains almost the same for the nearest and close next-nearest pairs and it is reduced approximately to one half of $`\mathrm{\Delta }E(d_1)`$ for the doubled distance corresponding to the configuration (b).

It should be pointed out, however, that the density-functional calculations tend to overestimate the strength of the exchange coupling. The reason for this is that the exchange splitting of the Mn d-states, i.e. the separation of occupied majority-spin and empty minority-spin states on the energy scale levels is systematically underestimated. According to our calculations, the exchange splitting $`\epsilon _d()\epsilon _d()`$ deduced from the spin-polarized spectral distribution of Mn d-states ranges from 2 eV to 3 eV. This is roughly one half of the realistic estimate for the exchange splitting in Mn (cf. e.g. Ref. Drchal91 ). Correspondingly, the above given value of $`J(d_1)`$ should be divided by four. In this way, we end with a value close to the result of Ref. Blinowski03 .

Finally, we estimate the binding energy of the Mnint-MnGa pair. We compare the energies corresponding to the (a) and (b) configurations from Fig. 2, both in the magnetic ground state with the antiparallel alignment of the local moments. The partially dissociated pair has a higher energy and the energy difference $`E_{}(d_2)E_{}(d_1)`$ is approximately 0.33 eV. It is, however, only the lower estimate for the binding energy of the Mnint-MnGa pair because the dissociation of the pair is far from being complete in our periodic model. Nevertheless, even the value obtained here indicates that the Mn interstitials are strongly attracted by the MnGa atoms. As long as the concentration of Mnint is lower than the concentration of MnGa we can assume that most of Mnint atoms are involved in the pairs and that the blocking mechanism proposed in Ref. Yu02  works.

## IV Exchange interaction of interstitial Mn and holes

Due to the hybridization with the spin-polarized Mn d-states, also the distribution of the states derived from the GaAs valence band depends on the spin. In particular, the valence-band states for the majority-spin electrons hybridize with the occupied d-states and are pushed to higher energies. The minority-spin states are, on the other hand, pushed down due to their hybridization with unfilled d-states.

This effect is formally described by the Kondo exchange interaction between the local spins $`𝐒_i`$ at sites $`R_i`$ occupied by Mn and the spin density $`𝐬(r)`$ due to the itinerant holes Moriya85 ,

$$H_{int}=J_{pd}\underset{i}{}𝐒_i𝐬(R_i).$$
(1)
The exchange parameter $`J_{pd}`$ characterizes the strength of the coupling. Within the mean-field theory, the Kondo exchange interaction results in the splitting of the valence band edge $`E_v`$. The splitting $`\mathrm{\Delta }E_v`$ is proportional to the size $`S`$ of the local spins and to the concentration $`x`$ of magnetic ions,

$$\mathrm{\Delta }E_vE_v()E_v()=\frac{4x}{a^3}J_{pd}S$$
(2)
assuming spin $`\frac{1}{2}`$ for the holes, $`a`$ is the lattice constant. We use Eq. (2) to determine the exchange parameter $`J_{pd}`$ from our spin-polarized band structures. The results of the calculations for the unrelaxed geometries presented here and in Ref. Maca02  are summarized in Table III.

The resulting values of $`J_{pd}`$ are overestimated from the same reasons as discussed above in Sect. III. Being divided by a factor of two (the reduction factor for the energy nominator), they approach the realistic values Ohno99 . In this work, however, we concentrate on the comparison of $`J_{pd}`$ for different geometries. The overestimate of $`J_{pd}`$ due to the reduced band gap is expected to be similar in all cases and not very important in this respect.

In Fig. 3, we plot $`\mathrm{\Delta }E_v`$ against $`x`$ to visualize the linear increase of the band splitting with the concentration of Mn. The most noticeable result is that all points in Fig. 3, obtained for substitutional Mn as well for Mn in the interstitial T(As4) and T(Ga4) positions lie around the same line. This means that the value of $`J_{pd}`$ does not depend much on the position of the Mn atom in the lattice.

This result is not very surprising in the case of the substitutional Mn and T(As4) position. In either case, Mn is surrounded by four As atoms at the same distance $`d_1`$ and the same degree of the hybridization of the valence band states with Mn d-states can be expected.

The results for the T(Ga4) position are not simple to interpret. The value of $`J_{pd}`$ for both Ga12MnAs12 and Ga16MnAs16 is smaller than but comparable with the result for the substitutional MnGa site. Assuming for simplicity that only the hybridization of the Mn d-states with p-states of the nearest neighbors is relevant for the spin splitting $`\mathrm{\Delta }E_v`$, the values of $`J_{pd}`$ simply reflect the orbital composition of the valence band. It is well known that the top of the valence band in GaAs, as well as in other III-V semiconductors, is composed of both anion and cation p-states, in the proportion approximately $`c_p^{\mathrm{As}}3/4`$, $`c_p^{\mathrm{Ga}}1/4`$ Harrison80 . This means that the same proportion (roughly 3:1) should be expected also for the exchange parameters $`J_{pd}`$ parameters corresponding to Mn atoms in substitutional and T(Ga4) positions. The strong deviation of the actual density-functional results from this model expectation may indicate that the hybridization of the Mn d-states with more distant neighbors has a remarkable influence on the magnetic interactions.

The fact that the magnetic behavior of the interstitial Mn is rather insensitive to its position in the crystal has two important implications concerning the MnGa-Mnint pairs. First of all, the pair as a whole does not interact much with the spin of the holes because the effects due to MnGa and Mnint compensate one another. In this respect, our ab initio results overcome the troubles of the simplified tight-binding models mentioned in the Introduction.

In addition, the effective ’annihilation’ of MnGa due to pairing with Mnint is not restricted to the closest pairs with Mnint in the T(Ga3Mn) position, but it works as long as the Mn atoms are close one to another and the exchange coupling in the pair remains antiferromagnetic.

## V Conclusions

We used the FPLAPW method to obtain total energies of supercells simulating various geometric and magnetic configurations of Mn atoms in (Ga,Mn)As. In absence of other defects, the ground state of the Mn interstitials is the tetrahedral T(As4) position. The energy of the T(Ga4) position, however, is almost the same. The situation changes in the p-type material where the T(Ga4) position has a lower energy. The hexagonal interstitial position has much higher energy and represents a barrier for diffusion of Mn from one to another interstitial site. The barrier, and consequently also the mobility of the interstitial Mn, depends on the doping.

The exchange coupling $`J_{pd}`$ of Mn interstitials with the holes in the valence band is, for both T(As4) and T(Ga4), close to the value of $`J_{pd}`$ obtained for the substitutional Mn. This is not consistent with the simplest tight-binding picture of Mn d-states that hybridize only with the nearest neighbors. In this way, our result indicate that the hybridization with more distant neighbors may be also important for the magnetic interactions.

The Mn interstitials are attracted to the substitutional Mn and form stable and magnetically inactive pairs. The density-functional estimate for both binding energy of the pair and for the energy of the antiferromagnetic coupling is of order of 0.3 eV. This fits well with the present day notion of the interstitial Mn in (Ga,Mn)As. Yu02 ; Blinowski03 . In contrary to the general opinion, however, we found that the efficient pairing is not restricted to Mnint in the T(Ga3Mn) position and we showed the importance of the close next-nearest neighbors for the properties of the interstitial Mn.

### Acknowledgement

This work has been done within the project AVOZ1-010-914 of the ASCR. The financial support provided by the Grant Agency of the ASCR (Grant No. A1010214) and by RTN project No. HPRN-CT-2000-00143 the EC is acknowledged.


text:
# Collapse of a semiflexible polymer in poor solvent

## I Introduction and motivation

The conformation of individual polymer chains depends on the solvent properties de Gennes (1979); Bird et al. (1987); Doi and Edwards (1986). In good solvent the monomers prefer being surrounded by solvent molecules, effectively repelling each other. This effect leads to a swollen coil conformation for flexible polymers in good solvent. In poor solvent, conversely, the monomers try to exclude the solvent and effectively attract one another, and a flexible chain forms a compact globule of roughly spherical shape to minimize the contacts between monomers and solvent. The dynamics of the coil-globule transition in flexible chains are well known deGennes (1985); Ostrovsky and Baryam (1994); Dawson et al. (1997); Halperin and Goldbart (2000) and involve the formation of a pearl necklace and the gradual diffusion of large pearls from the chain ends Chu et al. (1995); Abrams et al. (2002).

Many polymers, however, exhibit substantial bending stiffness, i.e., they are semiflexible. Such polymers can be described by the worm-like chain model; examples include biopolymers (F-actin, DNA) as well as synthetic polymers (xanthan, PBO, PPTA, kevlar). Semiflexible polymers in good solvents form open, extended structures, while in poor solvent the effective self-affinity of the chain contrast this tendency and leads to chain collapse. The collapsed state configurations and the pathways to their formation are the result of the interplay between two opposing potentials: the bending potential related to the chain stiffness and the attractive potential due to the poor solvency of the environment.

A compact globule is energetically unfavorable for semiflexible polymers because it involves highly bent states. The collapsed ground state is instead a torus, which reduces the monomer-solvent contacts without causing excessive bending penalty. Theoretical analysis Grosberg (1979); Bloomfield (1997), Monte Carlo simulations Noguchi et al. (1996) and experimental evidences Baeza et al. (1987); Martin et al. (2000) have detailed the stability, the features and the packing of the collapsed tori. The dynamics of the collapse of worm-like chains in poor solvent have been investigated more recently Schnurr et al. (2000); Noguchi and Yoshikawa (2000); Sakaue and Yoshikawa (2002), strongly suggesting a possible generic pathway for the collapse of semiflexible polymers, featuring a series of long-lived, partially collapsed, racquet-shaped intermediate states, before the eventual collapse to a torus. Such intermediate states form an energetically driven cascade of increasingly compact conformations with sharp transitions between their energy levels.

In the present work, we verify this hypothesis and further explore the effect of the relevant length scales on the phenomenon. The natural length scale for the polymer stiffness is its persistence length $`L_p`$, defined as the ratio of the bending stiffness $`\kappa `$ to the thermal energy $`k_BT`$, which represents the length along the chain for which tangent vectors remain correlated in good solvent. The balance between attractive forces and bending forces can be considered in term of the so-called “condensation length” $`L_0=\sqrt{kTL_p/u_0}`$, where $`u_0`$ is the value of the attractive potential per unit length. We choose as time scale for the collapse event the rotational diffusion time of a perfectly rigid rod with $`L=L_p`$, $`\tau _{\text{diff},L_p}=\zeta L_p^3/(72k_BT)`$, where $`\zeta `$ is the transverse friction coefficient fri . Therefore we use $`\tau _{\text{diff},L_p}`$ as unit of time and $`k_BT`$ as unit of energy in presenting all the results of the present work. We investigate the sequences and the likelihood of different collapsed structures formed at various solvent quality (described by $`u_0`$) and chain stiffness.

## II Simulation details

The dynamics of the continuous wormlike chain are described by a inertialess Langevin equation:

$$𝜻\left\{\frac{𝐫}{t}\dot{𝜸}𝐫\right\}=\frac{\delta U}{\delta 𝐫}+𝜼$$
(1)
Equation 1 balances the hydrodynamic force exerted on the chain by the surrounding fluid with the force due to the global potential acting on the chain $`U=U_{\text{bend}}+U_{\text{conn}}+U_{\text{solv}}`$, and the Brownian force per unit length $`𝜼`$, with correlations $`𝜼(s,t)𝜼(s^{},t^{})=2k_BT𝜻\delta (tt^{})\delta (ss^{})`$. $`U_{\text{bend}}=kTL_p_0^L𝑑s𝐮(s)/s`$ is the bending potential related to the chain stiffness, $`U_{\text{conn}}`$ is the connector potential, $`U_{\text{solv}}`$ describes the monomer-monomer interaction, and the Brownian term accounts for the rate of fluctuating momentum transfer due to the random collisions of small solvent particles with the chain.

The potential $`U_{\text{solv}}`$ chosen to account for the effect of poor solvent is an exponential attractive potential with an excluded volume repulsive term:

$`U_{\text{solv}}={\displaystyle 𝑑s𝑑s^{}\stackrel{~}{u}_0\mathrm{exp}\left(\frac{|𝐫(s)𝐫(s^{})|}{R_{\text{attr}}}\right)}`$
$`+{\displaystyle 𝑑s𝑑s^{}\stackrel{~}{u}_0\left(\frac{\sigma }{|𝐫(s)𝐫(s^{})|}\right)^{12}},`$ (2)
where $`\stackrel{~}{u}_0`$ defines the depth of the potential well, $`R_{\text{attr}}`$ is the range of the attractive forces and $`\sigma `$ is the radius of the excluded volume region. With this potential, $`u_0=R_{\text{attr}}\stackrel{~}{u}_0`$. Different choices for the shape of the potential should not change the stable and metastable configurations for the wormlike chain in poor solvent, while they may modify the speed and the dynamics of the collapse. Our choice of the exponential potential is related to the physics of the condensation phenomenon. DNA and polyelectrolytes collapse because of charge inversion and counterion induced attraction Grosberg et al. (2002). The attractive forces are therefore due to highly screened electrostatic charges, which typically show an exponential decay, mimicked by the potential selected in this work.

The Brownian Dynamics code used in this work solves in dimensionless, discretized form the Langevin equation of motion of a linear wormlike chain of $`N`$ beads with positions $`𝐑_1,\mathrm{},𝐑_N`$, connected by $`N1`$ connectors of equilibrium length $`aL/(N1)`$ with unit tangent vectors $`𝐮_i(𝐑_{i+1}𝐑_i)/a`$. For most of the simulations reported in this work, we use quadratic springs as connectors, with a force constant equal to 100 $`kT`$, which ensures a variation of the connectors’ length within 30 $`\%`$. This choice allows for longer time steps without interfering with the spring relaxation time. We use a midstep algorithm Grassia and Hinch (1996) to compute the bead trajectories generated by the equation of motion

$$\zeta _b\left\{\frac{d𝐑_i}{dt}\dot{𝜸}𝐑_i\right\}=𝐅_i=𝐅_i^{\text{bend}}+𝐅_i^{\text{elast}}+𝐅_i^{\text{solv}}+𝐅_i^{\text{rand}}.$$
(3)
Here, $`\zeta _b=\zeta a`$ is a bead friction coefficient, the bending force is:

$$𝐅_i^{\mathrm{𝖻𝖾𝗇𝖽}}=\kappa /a\underset{j=2}{\overset{N1}{}}(𝐮_j𝐮_{j1})/𝐑_i,$$
(4)
the elastic force due to the spring is

$$𝐅_i^{\text{elast}}=H\left[(|𝐮_{i1}|a)𝐮_{i1}(|𝐮_i|a)𝐮_i\right],$$
(5)
the potential force related to the effective interactions between the beads of the chain — repulsive within the hard core, attractive at short distance is

$$𝐅_i^{\text{solv}}=\stackrel{~}{u}_0\underset{ji}{}\left[\frac{\mathrm{exp}(\frac{|𝚫𝐑_{\mathrm{𝐢𝐣}}|}{R_{\text{attr}}})}{R_{\text{attr}}}\frac{12\sigma ^{12}}{|𝚫𝐑_{\mathrm{𝐢𝐣}}|^{\mathrm{𝟏𝟑}}}\right]\frac{𝚫𝐑_{\mathrm{𝐢𝐣}}}{|𝚫𝐑_{\mathrm{𝐢𝐣}}|}$$
(6)
where $`𝚫𝐑_{\mathrm{𝐢𝐣}}=𝐑_𝐢𝐑_𝐣`$, and $`𝐅_i^{\text{rand}}`$ is a random force.

Another relevant issue in the simulations is the accurate description of the high curvature of the collapsed states. The discretized computational model for the wormlike chain retains its validity when subsequent connectors are very nearly parallel: to satisfy this condition, a large number of connectors within the characteristic length scale of the collapsed structure is needed. The typical size of the torus and of the racquet heads are of the same order of magnitude of $`L_p`$. A good accuracy in the description of their curvature requires a large number of connectors per persistence length; for all the simulations presented here, we use 20 connectors per $`L_p`$.

We also performed simulations with a bead-rod model Pasquali et al. (2001); Pasquali and Morse (2002), which ensures local inextensibility of the worm-like chain. The computational cost for this model is however significantly higher ($``$ 10-50 times) than that for a bead-spring chain model, mainly because a smaller time step is required for convergence. While performing most of the simulations using the bead-spring model, we compared the results for a small number of conditions with the bead-rod model, as discussed in Section III.4. In that Section we also compare some of the results with that obtained using a bead-spring phantom chain, i.e., without excluded volume interactions.

For simplicity, we use isotropic drag in the simulations, and we neglect hydrodynamic interactions between segments of the chain. Including such interactions would likely give minor corrections to the evaluation of the time scales of the phenomena, but — because at the beginning of the simulations the chains are locally straight, and because after the collapse the hydrodynamic interactions would be dominated by the attractive forces — it would not change the kinetic pathways and the relative stability of the intermediate states.

## III Results

This section is organized as follows: Section III.1 reports qualitatively the results of our simulations, describing the various possible intermediate configurations and the typical pathways. Section III.2 presents a more quantitative description of the collapsed states in term of the correlation function between tangent vectors along the chain and of the energy levels associated with each different collapsed state, which confirm that the toroid is the stable configuration, and show how the collapse proceeds through a cascade of subsequently more energetic favorable metastable states — the various multiple headed racquets shapes, as predicted from theory Schnurr et al. (2002). The decay rate of uncollapsed states is investigated as a function of a proper combination of $`L`$, $`L_0`$ and $`L_p`$ in Section III.3, together with a statistical analysis of the preferred collapsed configurations. Section III.4 provides a comparison between the results obtained with the different computational models (phantom chain and bead-rod chain). Some preliminary results of the effect of shear flow on the kinetics and pathways of collapse are presented in Section III.5.

### III.1 Molecule configurations and collapse pathways

Figure 1 shows a typical snapshot of an ensemble of semiflexible chains during collapse. The 3-D structures of the molecules are shown here as 2-D projections on a convenient plane. At a given time molecules coexist in different configurations, from completely uncollapsed (molecule (a)), to partially collapsed (molecules (b), (f) and (g)), to higher order racquet (molecules (c) and (e)), to toroid (molecule (d)).

The fraction of collapsed molecules increases monotonically with time, confirming that collapsed states are energetically favorable in a bad solvent. Many other ensembles with different values of $`L/L_p`$ and $`L_0/L_p`$ show a similar behavior. However, if the attractive forces are not strong enough or the bending stiffness is too high there is no evidence of collapse and the open conformations persist in time. As already observed in Schnurr et al. (2000), a sequence of collapse can evolve following two different generic patterns. The first pattern involves the direct formation of a torus from the open conformation, when the two ends of the chain meet while the tangent vectors are locally parallel. In this instance, the chain collapses without formation of any intermediate, the end to end distance drops quickly to a value of the order of the torus dimension, and then fluctuates as the chain keeps folding into the final structure. Meanwhile the total energy of the chain keeps decreasing monotonically till reaching the equilibrium value. Conversely, the second pattern involves the formation of intermediate metastable states: initially a single racquet-headed shape which folds quickly into a short-lived, 3-heads racquet configuration, and then rearranges as a toroid via a subsequent folding. Figure 2 shows the two different pathways. We notice here that the direct formation of a torus is a more unlikely pattern in 3-D — and in fact observed with scarce frequency — than in 2-D. Whereas in 2-D the relative orientation of the two ends is described by a single angle, two angles are required in 3-D and both need to be close to $`2\pi `$ for the direct formation of a torus to take place. Most of the chain will therefore initially form a single headed racquet, as described below mov .

The time evolution of the total energy of an ensemble of collapsing chains is shown in Figure 3. The total energy $`U`$ of the chain during collapse is defined as the sum of the bending energy $`U_{\text{bend}}`$ and the interaction energy along the chain $`U_{\text{solv}}`$. In the discretized form,

$$U=kTL_p\underset{i=1}{\overset{N2}{}}\mathrm{cos}^2\theta _i+\underset{j=1}{\overset{N}{}}\underset{k=j+1}{\overset{N}{}}u_{\text{solv},jk},$$
(7)
where $`\theta _i`$ is the bond angle and $`u_{\text{solv},jk}`$ is the interaction energy evaluated between the beads $`j`$ and $`k`$ along the chain.

The time series in Figure 3 reports most of the possible evolutions observed in our simulations. Few molecules remain in the uncollapsed states for the whole duration of the run, equivalent to $`250\tau _{\text{diff},L_p}`$. For these molecules, the thermal fluctuations never lead two branches of the chain close enough together to initiate the collapse. The other molecules typically reach lower energy states by collapsing in sequentially higher order metastable conformations. The energy levels for each racquet-head shape of increasing order are well separated, and there is a rather large gap from the highest order racquet observed ($`6^{\mathrm{th}}`$ order) and the stable tori, as expected from theory Schnurr et al. (2002). Most of the chains initially form a single head racquet: after the first monomer-monomer contact and the formation of the racquet head, the two sides of the chain tend to become parallel and form a neck region, which ensure the maximum number of monomer-monomer contacts without increasing the bending penalty. The molecules then fold again on themselves, forming higher order shapes. Depending on the way in which the molecule folds back on itself, the next configuration can be a 2-heads or a 3-heads racquet shape. Figure 3 shows examples of both pathways. Successive folding leads to higher order racquets. The transition between each single metastable conformation is very rapid and driven by the deterministic attractive and bending forces. In contrast, the initiating event for each successive collapse is related to the Brownian motion of the molecule, and therefore the time interval between folding events is extremely variable, and completely random for the simulations that we have performed. In the simulations of longer chains ( $`L=10L_p`$ ) we observe combinations of the conformational elements described above, e.g., chains partially folded in a racquet shape at one end while still uncollapsed at the other. Those conformations have intermediate energies between the well-defined levels corresponding to racquets and tori and quickly disappear in favor of more compact and ordered structures. In contrast, the metastable high order racquets are extremely long-lived, and the energetic barrier for reaching the equilibrium shape is of the order of several $`k_BT`$. While the partial unfolding of a torus in a racquet is never observed, simulations show the formation of a torus via metastable configurations, and this, together with the lower total energy shown by the torus (see Figure 3), confirms once more that the latter is the equilibrium configuration.

### III.2 Quantitative description of collapsed states

The configurations and the shapes of semiflexible molecules in bad solvent can be studied through direct visualization of the collapsing chain. While this method is a valuable way to gather information on the kinetics and pathways of collapse, it is time consuming and not quantitative. We therefore define a spatial correlation matrix $`𝐌`$ to analyze in a more systematic way the collapsed shapes. The elements of $`𝐌`$ are defined as $`M_{ij}=𝐮_i𝐮_j`$. With this definition, $`𝐌`$ is symmetric; a perfectly straight rod has $`M_{ij}1i,j`$, as all the connecting normalized vectors have the same direction and orientation. An uncollapsed semiflexible chain with persistence length $`L_p`$ will have $`M_{ii}1`$, while the off-diagonal terms $`M_{ij}`$ would be smaller and decaying exponentially with the distance along the chain between the connectors $`i`$ and $`j`$.

When the chain collapses, the actual shape assumed by the molecule can be inferred from the spatial correlation between connectors. The spatial correlation matrix for a single head racquet, for example, will be a block matrix: one block will show values close to $`1`$, corresponding to the almost straight filament in the neck of the racquet, followed by a region with rapidly decreasing values from $`1`$ to $`1`$ — the racquet head — and by another block with values close to $`1`$, indicating the other filament of the neck , with same direction but opposite orientation.

The spatial correlation matrix can be simplified into a spatial covariogram for the molecule. For standard statistical analysis, the covariogram is defined as

$$C(h)=1/N(h)\underset{j=1}{\overset{N}{}}\underset{i=1}{\overset{N}{}}\{Z_i\mu \}\{Z_j\mu \}$$
(8)
where $`h`$ is the distance (in time or in space) between the observations $`Z_i`$ and $`Z_j`$, $`N(h)`$ is the number of observations at a distance $`h`$ and $`\mu `$ is the expected average value of the observation. In the present analysis, for each molecule we construct the spatial covariogram $`C_{\text{sp}}`$ assuming $`Z_i=𝐮_i`$ and $`\mu =\mathrm{𝟎}`$. We therefore obtain a scalar function of the distance along the chain $`s`$, bounded between $`+1`$ and $`1`$, which contains all the relevant information about the shape of the molecule.

Figure 4 shows three relevant examples of this type of analysis: the actual shape of the semiflexible molecules is shown in comparison with the spatial correlation matrix and the spatial covariogram. For a uncollapsed chain (Figure 4-a), the matrix shows randomly alternating regions of positive and negative correlation along the length of the chain, corresponding to the typical open configuration of a semiflexible chain in good solvent. The corresponding covariogram shows how the connectors correlation decays from 1 to values around 0 along the length of the chain; the decay is exponential at short distances and related to the stiffness of the chain, i.e., $`C_{\text{sp}}=\mathrm{exp}(s/L_p)`$. At longer distances, the correlation randomly oscillates between negative and positive values.

For a racquet shape with multiple heads, $`𝐌`$ has a structure of alternating blocks with strong positive and strong negative correlation, as shown in Figure 4-b. The 5 straight filaments in the neck of the racquet correspond to these 5 blocks in the matrix; the size of the blocks is approximately the same and indicates the length of the racquet neck. The covariogram shows 4 zero-crossing, each corresponding to a head of the racquet, and can be fitted quite accurately with a $`4^{\text{th}}`$ order harmonic; the parameters of this fit are related to the curvature of the racquet heads.

For a torus, the spatial correlation matrix appears completely different (Figure 4-c), showing a structure with diagonal bands. In fact, the correlation between pairs of connectors at the same distance is the same, regardless of their position along the chain; this translates into bands along the diagonals of the matrix. The distance between a diagonal with values close to $`1`$ and one with values close to $`1`$ corresponds to the diameter of the collapsed torus. The same information can be extracted from the covariogram, which is a perfect sinusoid: the period of the sinusoidal function corresponds to $`\pi `$ times the diameter of the torus, and the number of complete periods in along the chain indicates the number of loops formed by the collapsed molecule.

### III.3 Kinetics of collapse

The fraction of collapsed semiflexible molecules in a bad solvent increases monotonically with time, as already stated. We have studied the effect of the relevant length scales of the phenomenon, $`L_p`$, $`L_0`$ and $`L`$, on the kinetics of this collapse. By systematically varying the total length of the chain ($`L=3L_p`$, $`5L_p`$, $`8L_p`$ and $`10L_p`$), as well as the strength ($`\stackrel{~}{u}_0=0.25kT/a^2`$, $`0.5kT/a^2`$ and $`1kT/a^2`$) and the range ($`R_{\text{attr}}=5/40L_p`$ and $`1/6L_p`$) of the attractive forces, we have explored the effect of these parameters on the decay rate of uncollapsed molecules. We define the time to collapse $`t_{\text{coll}}`$ as the time required by the molecule to form the first metastable collapsed structure. We obtain from simulations the average value of $`t_{\text{coll}}`$ in an ensemble of collapsing chain under each condition, normalized with $`\tau _{\text{diff},L_p}`$. The inverse of $`t_{\text{coll}}`$ is taken as a measure of the decay rate $`\text{D}_{\text{coll}}`$ of uncollapsed molecules. The distribution of the times to collapse follows the expected lognormal distribution, typical of such dynamical processes, with a peak close to the average value and a long tail.

At a given value of $`L/L_p`$, the decay rate increases, i.e., the molecules collapse faster, for decreasing values of $`L_0`$, i.e., for increasing values of $`u_0`$. Similarly, the decay rate increases for longer chains at a given value of the attractive potential and persistence length. We can collapse the values of $`\text{D}_{\text{coll}}`$ on a single curve as a function of a proper combination of the two dimensionless ratios, $`L/L_p`$ and $`L_0/L_p`$, $`\left(L/L_p\right)^{1/3}\left(L_p/L_0\right)`$. Figure 5 shows the mastercurve for the decay rate.

Table 1 reports the statistics of chains conformations for different values of $`L/L_p`$ at two different times, i.e., t1 $`=100\tau _{\text{diff},L_p}`$ and t2 $`=200\tau _{\text{diff},L_p}`$. At t1, the the single head racquet is a conformation much more likely than the torus for all the explored values of $`L/L_p`$. This confirms that the direct formation of a torus is an unlikely event, as discussed in Section III.1, and most of the chains collapse through a sequence of folding events. At t2, the percentage of molecules in the 1-head racquet configuration does not vary substantially with respect to t1, because — while uncollapsed molecules fold to form new 1-head racquets — some of these configurations further collapse in higher order metastable racquets and stable tori. The number of chains collapsed in tori or multiple head racquets monotonically increases within the observed span of time. For shorter chains the likelihood of tori and multiple racquets is similar, while the metastable states seem to be initially more favorable for longer chains ($`L/L_p=8`$ and $`L/Lp=10`$). Overall, the results show that the intermediates conformations in the collapse phenomenon are extremely long-lived: only $`8.4\%`$ of the simulated chain has reached the equilibrium, torus state after $`200\tau _{\text{diff},L_p}`$.

### III.4 Comparison of different computational models

We have performed simulation of collapsing chains with different computational models, in order to verify the results obtained and the appropriateness of the model chosen. In particular, we have performed simulations with a bead-rod model, where the chain connectors are rigid rods rather than stiff elastic spring. As noted in Section II, the bead-rod model is a constrained model which ensures local inextensibility of the chain, but is computationally more expensive, due to more stringent time step requirements for convergence. While local inextensibility is fundamental when analyzing the polymer contribution to the stress tensor Pasquali et al. (2001), its effect on the pathways and dynamics of collapse is less relevant, as confirmed by our results. In fact, we observe qualitatively similar pathways to collapse in the simulations performed with the bead-rod model. The preferred pathway includes the formation of metastable racquet shapes, similarly to what observed in the bead-spring simulations. The kinetics of the collapse are slower for the locally inextensible model; we believe that this is due to how the constraints are imposed, i.e., by projecting the random forces onto the constraints, and therefore not allowing any fluctuations along the chain. A more detailed comparison between the collapse kinetics two model is prevented by the high computational cost of simulating the constrained system.

We also perform simulations with a phantom bead-spring chain, i.e., a chain with no repulsive hard core. Once more, we find that the dynamics of collapse are qualitatively the same, although the actual shapes of the torii and of the racquets formed by the phantom chain are slightly different. The absence of hard core repulsion allows in fact for tighter structures, as shown in Figure 6. In particular, in the metastable multiple racquets the neck for phantom chains appears much narrower than for the chain where excluded volume interactions are taken into account, while it is noticeable that the shape and the curvature of the racquet head appears unchanged. The traces in time of the energy level for the phantom chain show the features already discussed in Section III.1: well defined energy states corresponding to metastable structures, with rapid transitions between them (see Figure 6, inset). However, the average life of metastable states appear slightly shorter than in the complete model: intuitively, this is due to the wider range of motion of the phantom chain, which permits folding paths otherwise prevented by the excluded volume interactions.

### III.5 Effect of shear flow

We discuss here the results of some preliminary work on the effects of steady shear flow on the collapse dynamics of semiflexible chains, In particular we monitor the decay rate $`\text{D}_{\text{coll}}`$ at different flow strengths, and we compare the pathways of collapse with those at equilibrium. Flexible chains expand in shear flow, i.e., the average end to end distance increases while the molecules tend to orient with the flow; the effect of the flow is therefore counteracting that of the attractive forces, which tend to form compact globula. Stiff semiflexible chains, with $`L_pL`$ shrink in shear flow, due to a buckling instability Montesi et al. (2003); such rigid chain however would not collapse to form tori and multiple racquets in bad solvent, as the open, rod-like shape is stable at equilibrium Schnurr et al. (2002). In this work, we consider the collapse dynamics of less stiff semiflexible molecules: is not clear a priori what to expect for the collapse dynamics under shear flow.

We show here that strong enough shear flow speeds up sensibly the collapse kinetics, increasing the likelihood for the two ends of the chain to meet. Figure 7 show the trend of the decay rate as a function of Wi, where Wi $`=\dot{\gamma }\tau _{\text{diff},L_p}`$, for 3 different values of $`L/L_p`$. At Wi $`=0.1`$, the decay rate does not change much with respect to the equilibrium value, and there is not a clear trend among the different $`L/L_p`$. At Wi $`=1`$, the increase in decay rate becomes instead significant, as the average time to first contact halves; for Wi $`=10`$, the effect of shear flow is completely dominating and the decay rate under the 3 different conditions is almost equivalent. Under shear flow, the semiflexible molecules undergo a sequence of compressions and extensions while tumbling in the plane of shear; in the compression region, the molecule tends to shrink, therefore raising the probability that two sections far apart along the chain backbone come in contact. The attractive forces then trap the molecule in the partially collapsed configuration, preventing it from straightening back in the extension region. However, this continuous tumbling of the chain and alternation of extension and compression does not allow the formation of the higher order, compact structures observed at equilibrium. In particular, the shear flow seems to inhibit the formation of the torus, which was never observed in the high shear flow simulations. Figure 8 shows two typical configurations observed in strong shear flow; while they show the basic features of the metastable racquets, the structures are less compact and do not reach the same level of energy of the equilibrium conditions.

## IV Discussion

The dynamics of collapse of semiflexible molecules in bad solvent has been investigated in the literature. Direct, visual experimental evidence of the formation of torii and high order racquets has been obtained studying biopolymers with sufficiently large persistence lengths Baeza et al. (1987); Martin et al. (2000). Theoretical work and Monte Carlo simulations Bloomfield (1997); Noguchi et al. (1996); Schnurr et al. (2002) have confirmed that the torus is the stable conformation, but intermediate, metastable racquet shaped conformations exist and are long lived. Previous Brownian Dynamics simulations in 2-D Schnurr et al. (2000) have firstly suggested the possible general pathways to the dynamics of collapse.

The present work provides new insights on the collapse of semiflexible chains, confirming the suggested pathways through 3-D simulations, and showing the effect of the two relevant dimensionless ratios $`L/L_p`$ and $`L_0/L_p`$ on the kinetics of collapse. The time series of the total energy of each chain show clearly the fast transition between well defined energy levels, corresponding to different metastable intermediate shapes, and the final, stable torus. We also showed that two possible pathways towards the stable conformation are possible: direct formation of a torus and successive folding of the chain in progressively higher order racquets. The latter one appears to be more probable for the range of parameters investigated. The decay rate of uncollapsed states for all the different conditions of the simulation can be plotted as a function of a proper combination of $`L/L_p`$ and $`L_0/L_p`$. Further work is needed to confirm the validity of this empirical scaling and verify it in a wider range of parameters.

We also performed simulations with different computational models, therefore confirming that the observed metastable conformations and pathways are general and do not depend on the choice of the model. However, the kinetics, i.e., the actual value of the decay rate can be different: in particular, for a constrained bead-rod model the kinetics of collapse appear slower, possibly due to the way the Brownian noise is projected onto the constraints.

We show preliminary results on the kinetics and pathways of collapse of semiflexible chains in shear flow. Strong enough shear flow increases the decay rate of uncollapsed states, but prevents the formation of the compact, well defined configuration observed at rest, and seems to inhibit the formation of the torus. While these results have been obtained in unbounded shear flow, future work will consider the interaction of the collapsing molecules with a solid surface, which could be attractive or repulsive for the chain. Also, future possible direction of research include the effect of local defects along the chain, which make the molecule locally more flexible or stiffer.

The authors wish to thank Bernhard Schnurr and F. Gittes for useful discussions, and R. Adam Horch for preparing the visualizations of the collapsing DNA molecules. This work was partially supported by the National Science Foundation under grant PHY99-07949, through award CTS-CAREER-0134389 and through the Nanoscale Science and Engineering Initiative award EEC-0118007. Computational resources were provided by the Rice Parallel Computational Engineering Cluster (NSF-MRI-0116289) and by the SARA Computing and Networking Services in Amsterdam through Vrije Universiteit.