text:
# Translocation of structured polynucleotides through nanopores

## General theoretical framework

Fig. 1 depicts schematically the driven translocation of a structured polynucleotide from the cis to the trans side of the pore. We seek here a convenient reduced description of this translocation process, rather than modeling the full three-dimensional polymer dynamics explicitly. Our approach is similar in spirit to the existing models for the case of unstructured polymers Sung96 ; Muthukumar99 ; Lubensky99 , where the translocation dynamics is formulated in terms of a single variable, e.g. the number of nucleotides, $`m`$, on the cis side, see Fig. 1. The dynamics, $`m(t)`$, is stochastic and can be described by ‘hopping rates’, $`k_{}(m)`$ and $`k_+(m)`$, for forward and backward motion of the nucleotide chain through the pore with a stepsize of one monomer. The external force on the molecule leads to an imbalance in the hopping rates, $`k_{}(m)>k_+(m)`$, and hence a mean drift towards the trans side. For unstructured molecules the one-dimensional description is permissible, if the relaxation of the polymer degrees of freedom on both sides of the pore is faster than the hopping process. This assumption does not hold for arbitrarily long polymers, since the relaxation time increases with the polymer length Lubensky99 ; Chuang01 , however for lengths on the order of a thousand bases, the one-dimensional description is adequate under typical experimental conditions Lubensky99 . The residual effect of the polymer ends is then only to introduce an entropic barrier for translocation, which leads to a weak $`m`$-dependence of the hopping rates.

For structured molecules, the translocation dynamics is considerably more complicated, since the dynamics of the ‘reaction coordinate’, $`m(t)`$, is then coupled to the dynamics of the basepairing patterns on both sides: the structure on the cis side, $`𝒮_{cis}(t)`$, affects the forward rate, while the structure on the trans side, $`𝒮_{trans}(t)`$, affects the backward rate,

$`m`$ $`\stackrel{k_{}(m,𝒮_{cis}(t))}{\text{}}`$ $`m1`$
$`m`$ $`\stackrel{k_+(m,𝒮_{trans}(t))}{\text{}}`$ $`m+1.`$ (1)
In two limiting cases however, the process can be modeled by a one-dimensional Brownian walk as for unstructured molecules, but with a complex sequence/structure-dependent free energy landscape $`(m)`$ along the coordinate $`m`$: (A) If the dynamics of the basepairing patterns $`𝒮_{cis}(t)`$ and $`𝒮_{trans}(t)`$ is much faster than the hopping process, the landscape is determined by the ensemble free energy of all basepairing patterns on the cis and trans side. (B) In the opposite limit, the basepairing pattern on the cis side is essentially frozen and is unzipped basepair by basepair as it is driven through the pore. The landscape is then determined by the basepairing energetics of the particular molecular structure prior to translocation, see below. In both cases, the free energy naturally decomposes into three parts,

$$(m)=_{cis}(m)+_{trans}(m)+_{\mathrm{ext}}(m),$$
(2)
where $`_{cis}(m)`$ and $`_{trans}(m)`$ denote the intrinsic binding free energies of the cis and trans parts of the molecule, while $`_{\mathrm{ext}}(m)`$ describes the effect of the external force. Given $`(m)`$, the simplest form for the hopping rates $`k_\pm (m)`$ which satisfies the detailed balance condition $`k_+(m)/k_{}(m+1)=e^{\beta [(m+1)(m)]}`$ (with $`\beta =1/k_BT`$) is

$$k_\pm (m)=k_0e^{\beta \mathrm{max}\{(m\pm 1)(m),\mathrm{\hspace{0.17em}0}\}}.$$
(3)
Here, $`k_0`$ denotes a microscopic rate constant, which can in principle be tuned by adjusting the properties of the pore. It can be interpreted as a friction coefficient and corresponds approximately to the bare hopping rate for unstructured molecules at zero external force (typical experimental estimates for $`k_0`$ are on the order of $`10^5\mathrm{s}^1`$ Meller01 ). The dynamics of the translocation process, as described by Eqs. (2) and (3) is dominated by energetic barriers due to basepairing, whereas the above-mentioned entropic barrier is completely negligible for structured molecules. These energetic barriers lead to arrests during translocation, as clearly observed already in the experiments with simple hairpins Vercoutere01  and double-stranded DNA Sauer-Budge03 .

## Pulling through a pore

Qualitative aspects. We now make use of the theoretical framework constructed above to investigate which information on structured molecules could be derived from pore translocation experiments. To this end, it is useful to compare unzipping by driven translocation through a nanopore with the more conventional way of unzipping by applying a force on the ends of a biopolymer, see e.g. Rief97 ; Liphardt01 ; Onoa03 . As illustrated in Fig. 2, the two approaches differ fundamentally: Pulling on the ends induces a spontaneous unfolding order for the individual structural elements, which is a function of their relative stabilities and the topology of the structure. In contrast, the nanopore prescribes a linear order along the sequence, and unfolds an RNA molecule much as enzymes such as the ribosome do in cells. This difference suggests that the two approaches can also yield different types of information about the molecule under study. As demonstrated by Onoa et al. Onoa03 , clever use of the pulling on the ends approach can reveal detailed information on the (un)folding pathway of an RNA molecule with known structure. However, when the structure of an RNA molecule is unknown, pulling on the ends can provide, by itself, little information beyond a count of the number of structural elements that unfold separately Onoa03 ; Gerland03 . In the following we therefore focus on the question of how much structural information may in principle be obtained with the nanopore approach.

Let us suppose that we were able to observe the trajectories $`m(t)`$ of the molecules during the translocation process. We could then assign a position within the sequence to each arrest during translocation. Since an arrest is caused by a kinetic barrier, i.e. a stem trapped at the entrance to the pore, we could thereby identify the positions of the stems in the structure. Such information can indeed be sufficient to reconstruct almost the entire basepairing pattern of a molecule, as we demonstrate explicitly using an example below. If the translocation dynamics is in the strongly driven limit (B) where the structure on the cis side is essentially frozen, then the reconstructed structure would correspond to the initial structure of the molecule before translocation. We concentrate on this limit in the following, including a discussion of its attainability. However, it may be noteworthy that in the slow translocation limit (A) one would also obtain useful structural information, namely on the average structure of the molecule (with respect to the thermodynamic ensemble of all structures McCaskill90 ). As long as the molecule is ‘well-designed’ this average structure will be dominated by the ground-state, i.e. the minimum binding free energy structure111The worst case for the purpose of structure determination corresponds to the regime where the typical timescale for the translocation of say a single hairpin is comparable to the timescale for structural rearrangements involving the formation of new stems: in this case, the structure on the cis side may relax after a stem is unzipped, so that one would oberve only the signatures of the relaxed structure rather than the original structure. This regime should be avoided by a proper choice of the driving force and the friction coefficient of the pore (R. Bundschuh and U. Gerland, to be published)..

How could one possibly observe the trajectories $`m(t)`$ during translocation? For the purpose of structure determination, we will need $`m(t)`$ with a spatial resolution below the typical length of a stem in an RNA structure (5–10 basepairs). This may be achievable through a refinement of the current nanopore technology, such that careful analysis of the ionic current allows a count (or even sequencing) of the bases that have passed the pore Akeson99 ; Meller00 . With artificial solid-state pores Li01 ; Storm03  it is also conceivable to use a tunneling current through leads within the membrane as a probe to count (or sequence) the bases as they pass through the pore. Here, we explore yet another option, namely pulling the molecule mechanically through the pore, with a device that can record force-extension curves, e.g. an atomic force microscope or optical tweezers. The explicit discussion of this case with an exemplary RNA sequence serves us to gauge the more general capability of nanopores as single-molecule tools for the study of biopolymer folding.

Quantitative aspects. Mechanical unfolding of a biopolymer yields characteristic sawtooth-shaped signatures in the force-extension curve (FEC) indicating the opening of structural elements within the molecule, see e.g. Rief97 ; Onoa03 . From the relative positions of these sawteeth one can determine length changes within the molecule with an extremely high resolution of about 1 nm. In the usual setup where the molecule is unfolded by pulling on its ends, such length changes can only be used to infer the ‘stored length’ of a structural element, but not its precise position along the backbone of the molecule, cf. Fig. 2. In contrast, for mechanical pulling through a pore, the relative positions of the resulting sawteeth will correspond directly to the relative positions of the structural elements in the sequence222The absolute position can be inferred by adding a known structural element, e.g. a strong C-G hairpin, to one end of the RNA, which can then function as a reference point.. One conceivable way to prepare the initial condition where an RNA molecule is almost entirely on the cis side, with one end threaded through the pore and attached to a pulling device on the trans side, is to start with an attached molecule on the trans side and to apply a voltage pulse across the pore that suffices to drive the molecule as far as possible to the cis side.

To apply our general model to the particular case of mechanical pulling in the strongly driven (fast pulling) limit, we need to specify the form of the three terms in the free energy landscape (2). The second term, i.e. the binding free energy on the trans side, may be set to zero,

$$_{trans}(m)=0,$$
(4)
since the reformation of structure after translocation is suppressed at high tensions in the RNA single strand333For instance, Liphardt et al. Liphardt01  observed refolding rates for a single hairpin around 1 s-1 at the unfolding force $`f_{1/2}14`$ pN. At a pulling speed of say $`1\mu `$m/s, the translocation of an RNA molecule with a thousand bases would therefore be terminated before refolding of a structural element on the trans side occurs.. The third term, $`_{\mathrm{ext}}(m)`$, describes the effect of the mechanical stress on the RNA, which stretches the single-stranded trans part of the molecule. The elastic response of this single-strand may be modeled by a freely jointed chain (FJC) polymer model. Assuming for simplicity a constant pulling speed $`v`$, the third term then takes the form

$$_{\mathrm{ext}}(m)=_{\mathrm{FJC}+\mathrm{spring}}(vt;Nm).$$
(5)
Here, the function $`_{\mathrm{FJC}+\mathrm{spring}}(R_t;n)`$ denotes the combined free energy of a single-stranded RNA of $`n`$ bases in series with a linear spring, stretched to a total extension $`R_t=vt`$ Gerland03 . (The linear spring takes into account the stiffness of the force-measuring device, see the Appendix for details.) By assumption, the first term, $`_{cis}(m)`$, represents the binding free energy of the remaining part of the initial structure on the cis side. $`_{cis}(m)`$ can be calculated for any initial structure, based on the free energy rules for RNA secondary structure Walter94  with a natural extension for pseudoknotted structures, see the Appendix. Our assumption of a frozen structure on the cis side is most likely an oversimplification for realistic pulling speeds, since small fluctuations in the secondary structure are known to occur already on timescales on the order of tens of microseconds Bonnet98 . However, since the pore pulling approach is sensitive only to stem positions, we expect that it is unaffected by small fluctuations and sensitive only to major rearrangements which significantly change the secondary structure. Such rearrangements are typically slow, sometimes even on the timescale of hours Pan98 ; Zhuang00 .

Reconstruction of secondary structures. To illustrate the problem and the method, we use an exemplary RNA, the well-studied self-splicing intron of Tetrahymena thermophila Cech93  with a sequence of 419 bases (Genbank # V01416). In its correctly folded active state, the basepairing pattern of this ribozyme contains a pseudoknot (see Fig. 3a), while its best characterized long-lived folding intermediate Pan98 ; Zhuang02  has a known alternative structure without pseudoknot (see Fig. 3b) Pan98 . We will investigate whether one can in principle use the pulling-through-a-pore approach not only to discriminate between these two different conformations in individual molecules, but also to reconstruct both structures from the FEC’s.

To obtain FEC’s for these structures, we performed Monte-Carlo simulations of the stochastic process defined by Eqs. (General theoretical framework3), and used Eq. (6) from the Appendix to calculate the force and extension time traces. We performed all calculations at the same pulling speed ($`v=0.1`$ nm/time step, which roughly corresponds to 10 $`\mu `$m/s given typical values for $`k_0`$, see above), and the same stiffness of the force-measuring device ($`\lambda =0.5`$ pN/nm). Fig. 4 displays three such FEC’s (solid lines) for the non-pseudoknotted structure of Fig. 3b, corresponding to unzipping from the 3’ end. These FEC’s show the sawtooth-like behavior which is characteristic for the sequential opening of structural elements (a very similar behavior was observed in the experiments of Onoa et al. Onoa03  where the molecule was rapidly unzipped by pulling on its ends). The rising parts of the sawteeth correspond to stretching of single strand on the trans side as a stacked region is “trapped” in front of the pore on the cis side. When a stacked region opens, some single strand is freed to pass the pore, which leads to relaxation of the tension and causes the downstrokes in the FEC’s. Note that the FEC’s do not share all of their sawteeth, which reflects the importance of thermal fluctuations for this type of single molecule experiments (this property is manifest also in the experiment of Onoa et al. Onoa03 ).

The most relevant information contained in the FEC’s are the positions of the translocation arrests, during which the required force for the opening of basepairs is built up. To extract these positions, we use FEC’s of freely jointed chains with different lengths: The dashed lines in Fig. 4 show some examples of such FEC’s where the chain length $`n`$ coincides with the length of the RNA single strand on the trans side during such an arrest. With an automated procedure described in the Appendix we obtain all of these positions (above a threshold for the duration of an arrest).

Since the bases around the position of an arrest are very likely basepaired with another segment of the sequence further to the 5’ end, we represent this information by a closing angular bracket, ‘$``$’, above that position in the RNA sequence (written from 5’ to 3’), see Fig. 5. Of course, the molecule can also be pulled through the pore in the other direction, i.e. from the 5’ end. This yields information on the positions of segments that have downstream binding partners. The same procedure then leads to the opening brackets, ‘$``$’, also shown in Fig. 5.

Bracket representations are a widely used short hand notation for RNA secondary structures. For the structures in Fig. 3, such a representation is shown in the third row of Fig. 5. Note that two types of brackets have to be used for the pseudoknotted native structure, in order to make the association between opening and closing brackets unambiguous. We observe that the angular brackets extracted from the FEC’s can be viewed as an incomplete bracket representation of the RNA secondary structure. Can we complete it using only the given sequence of the RNA molecule?

This task is a sequence alignment problem, which consists of matching each opening (closing) bracket with an associated downstream (upstream) binding sequence. Several circumstances conspire to make this, somewhat surprisingly, a nontrivial problem: (i) stems, i.e. contiguous basepaired regions, are usually short, typically 5–10 basepairs, (ii) structural elements often lead to a different number of angular brackets in the two pulling directions, i.e. not every opening bracket has a corresponding closing bracket and vice versa, and (iii) sequence segments containing several U’s have many possible binding partners, since U’s can pair with A’s and G’s.

To overcome this problem, we developed a probabilistic sequence alignment algorithm (see Appendix), which identifies the most likely set of stems that is consistent with all angular brackets and where all paired sequence segments contain at least one angular bracket on each side. The output of this algorithm is shown in the fourth rows of Fig. 5, where lower case letters indicate paired sequence segments and the alphabetic order represents the confidence level (confidence is largest for ‘a’). In Fig. 3 the bases involved in this reconstructed set of stems are colored, with green (red) indicating (in)correct basepairing. We observe that the two different basepairing patterns (for the same sequence) are clearly distinguished and the large scale secondary structure is captured in both cases. In particular, the pseudoknot in the native structure is correctly identified. The only incorrectly predicted stem is the least significant one (‘k’) in the pseudoknotted structure.

While these results seem satisfactory as a proof of principle, we stress that our reconstruction algorithm can certainly be improved upon, e.g. by allowing for mismatches in longer stems, which should help to fill in many of the missed basepairs. Also, one could make use of the known basepairing energies in the reconstruction.

## Discussion and Outlook

Our theoretical study has led us to a simple coarse-grained model, Eqs. (2-3), for the translocation of structured polynucleotides, which is applicable in the two opposite limits of very slow and very rapid translocation. This model is a useful starting point for a more detailed description that remains valid in the entire parameter regime. Here, we have applied the model to demonstrate that the physics of the translocation process can in principle be exploited to use nanopores for secondary structure determination (including pseudoknots) on the single-molecule level. Indeed, the nanopore technique would be a useful addition to the existing repertoire of structure determination methods: RNA secondary structure can be predicted computationally to some extent Zuker81 ; McCaskill90 ; Vienna94  based on experimentally determined free energy rules Walter94 , however this approach is unreliable for RNA molecules exceeding $`100`$ bases and cannot take pseudoknots properly into account. Including pseudoknots, which are often crucial to the function of RNA enzymes Pan98 ; Wadkins99 , is not only computationally expensive Rivas99 ; Isambert00 , but is also limited by a lack of experimental information on the corresponding binding free energies. Experimentally, X-ray crystallography Cate96  or NMR Colmenarejo99  provide detailed structures, but these techniques are cumbersome and limited to small molecules or isolated domains of larger RNAs. Structural information for larger RNAs can currently only be obtained from comparative sequence analysis Gutell92 , which requires large sets of homologous RNA sequences, or from indirect biochemical methods Pan98 .

Throughout this paper, we have focused on basepairing only, which is permissible under ionic conditions that disfavor tertiary interactions, e.g. low sodium and no magnesium. However, once the translocation of a molecule is well characterized under these conditions, it becomes interesting to switch to the native ionic conditions and examine the effects of tertiary interactions. Generally, one can expect more cooperativity in the presence of tertiary interactions, i.e. larger domains will open in a single step, as observed by Onoa et al. Onoa03 . This suggests a hierarchical approach to structure determination with nanopores: first unzip under low ionic conditions to obtain the secondary structure, and then repeat in the presence of magnesium to identify how the secondary structure elements are grouped into larger tertiary structure domains (such as the P4-P6 domain in the Tetrahymena ribozyme). It is worthwhile to stress the advantage of RNA as a model system to separately study the effect of secondary and tertiary structure. In contrast, the secondary structure of proteins is not stable in the absence of tertiary structure, and hence one may expect that single-domain proteins will unfold and translocate across a pore in a single step.

Nanopores could in principle also be used to probe the kinetics of large-scale secondary structure rearrangements in single-molecules. For instance, it would be useful to attach larger objects to both ends of a molecule that is already threaded through the pore, allowing the same molecule to be driven forth and back through the pore, over and over again. By varying the time interval between successive reversals of the driving force, one could then probe structural relaxation over a broad range of time scales. More generally, nanopores may emerge as a new tool to probe intra- and inter-molecular interactions in single biomolecules. For instance, one could probe the biophysics of combined binding and folding in the context of RNA-protein interactions.

Acknowledgments. We are grateful to D. Branton, S. Ling, J. Liphardt, and D. Lubensky for stimulating discussions. This research was supported by the National Science Foundation through grant No. 0211308, 0216576, and 0225630.

## Appendix

Calculation of free energy landscape. Given a secondary structure of the molecule, we obtain $`_{cis}(m)`$ by eliminating all basepairs involving the terminal $`Nm`$ bases, and calculating the binding free energy of the remaining structure according to the free energy rules for RNA secondary structure Walter94 . \[We take the free energy parameters as supplied with the Vienna RNA package (version 1.3.1) at room temperature $`T=25^oC`$. The salt concentrations at which these parameters were measured are $`[\mathrm{Na}^+]=1\mathrm{M}`$ and $`[\mathrm{Mg}^{++}]=0\mathrm{M}`$.\] For pseudoknotted structures, the free energy rules currently include no prescription, however the following extension appears reasonable: we first eliminate basepairs in stems that give rise to the pseudoknot(s) and calculate the free energy of the remaining structure according to the standard rules. We then add the free energies of the eliminated stems separately, including the free energy for the loops created by these stems, again according to the standard free energy rules (however, the bases in these loops that are involved in other stems are removed before calculating the loop free energy).

The trans part of the molecule is tethered at both ends, by the pore and the pulling device, respectively. The pulling device can be described by a linear spring, while the configurational entropy of the RNA single strand can be modeled by a freely jointed chain (FJC) with extensible segments. \[For the few bases that are inside the pore, we neglect the effect of the confinement on the entropy.\] We denote by $`R_t`$ the total extension of the trans part in series with the linear spring. The free energy (5) can then be expressed in terms of the total end-to-end distance distribution $`W_{\mathrm{FJC}+\mathrm{spring}}`$,

$$_{\mathrm{FJC}+\mathrm{spring}}(R_t;n)=k_BT\mathrm{log}W_{\mathrm{FJC}+\mathrm{spring}}(R_t;n),$$

which can in turn be written as the convolution of the individual end-to-end distance distributions of the FJC and the spring Gerland03 ,

$$W_{\mathrm{FJC}+\mathrm{spring}}(R_t;n)=\underset{0}{\overset{\mathrm{}}{}}dRW_{\mathrm{FJC}}(R;n)W_{\mathrm{spring}}(R_tR).$$

Here, $`W_{\mathrm{spring}}(R_s)=\mathrm{exp}(\beta \lambda R_{s}^{}{}_{}{}^{2}/2)/\sqrt{2\pi /\beta \lambda }`$, where $`\lambda `$ denotes the inherent stiffness of the pulling device. We calculate the end-to-end distance distribution of the freely jointed chain, $`W_{\mathrm{FJC}}(R;n)`$, as described previously Gerland01 . The polymer parameters we use were obtained from a fit Mezard01  to FEC’s of single-stranded DNA Maier00  (base-to-base length $`0.7`$nm, Kuhn length $`1.9`$nm, and stretch modulus $`815`$pN), since we are unaware of corresponding data for the chemically very similar RNA.

Calculation of FEC’s. We obtain several trajectories $`m(t)`$ with a Monte Carlo simulation of Eqs. (23) with $`m(0)=N`$, $`R_t(0)=0`$ as initial condition and incrementing $`R_t`$ at the constant rate $`v`$. The simulation is stopped when all bases have translocated ($`m=0`$). From the time trace $`m(t)`$, we calculate the force-extension curve $`f(R)`$ using

$$f=\frac{}{R_t}_{\mathrm{FJC}+\mathrm{spring}}(R_t=vt;Nm(t))$$
(6)
and $`R=vtf/\lambda `$. Here, $`f`$ and $`R`$ are both thermal averages over the polymer and spring degrees of freedom at fixed total extension $`R_t`$ and fixed basepairing pattern.

Extraction of parentheses positions from FEC. For every point on a FEC, we determine the length $`n`$ of the freely jointed chain whose FEC passes closest to the point (using the polymer parameters for single-stranded RNA as given above). We take a histogram of the resulting lengths $`n`$ over three independent FEC’s for each structure. In this histogram, the lengths $`n`$ that correspond to start positions of stably basepaired regions appear as peaks, since the length of single-stranded RNA on the trans side remains approximately constant while the force required to unzip the basepairs builds up. \[A similar procedure was applied in Ref. Koch02  to identify the positions of proteins bound to double-stranded DNA as it is being unzipped.\] We keep all $`n`$-values where the histogram exceeds a threshold of $`30`$ counts (a count is made every Monte Carlo time-step). Since thermal noise makes the molecule fluctuate back and forth by a few bases while the force is building up for the next stem to open, we pick out of each contiguous stretch in the remaining $`n`$-values only the largest. Finally, we increment the extracted $`n`$-values by one and mark the corresponding position in the sequence with a parenthesis.

Reconstruction of basepairing pattern. The FEC’s do not reveal which opening and closing parentheses are paired with each other. However, given the sequence of the RNA, we can match the parentheses by sequence complementarity. \[To keep the number of false basepair predictions to a minimum, we consider only stems where we have at least one parenthesis at each end.\] Here, we summarize the essential steps in our sequence alignment algorithm, while a detailed presentation and characterization will be given elsewhere (R. Bundschuh and U. Gerland, to be published): First, we find all possible gapless local alignments between a subsequence containing a parenthesis and subsequences to the open side of the parenthesis, using the scoring scheme 2 for GC, 1 for AU, and 0 for GU. We keep only those alignments with a score larger than 5 and where the matching sequence segment also contains a matching parenthesis. We consider the remaining alignments as possible stems in the secondary structure. To pick the most likely set of mutually consistent stems, we assign an alignment E-value to each stem Karlin90 . We then iteratively include the most likely stem into the structure prediction, and remove all other stems it excludes due to overlapping basepairs from the list of allowed stems.


format: text

text:
# Relation between Wigner functions and the source functions used in the description of Bose-Einstein correlations in multiple particle production.

## 1 Introduction

When discussing Bose-Einstein correlations among identical pions111Here and in the following we call for definiteness the identical bosons pions the results, however, are valid for any identical bosons produced in multiple particle production processes, one often uses the source function $`S(X,K)`$ (cf. e.g.  and references quoted there) related to the single particle density matrix in the momentum representation by the formula222Often an equivalent formula with the density matrix replaced by the average
$`<\widehat{a}_p^{}\widehat{a}_p^{}>=\rho (𝐩^{},𝐩)`$ is used. Also some authors include on the left hand side a factor $`\sqrt{EE^{}}`$. Since this factor does not affect our argument and is easy to introduce at any stage, we skip it.

$$\rho (\text{p},\text{p}^{})=d^4Xe^{iqX}S(X,K),$$
(1)
where $`q=pp^{}`$ and $`K=(p+p^{})/2`$ are four-vectors constructed from the on-shell single particle four momenta $`p`$ and $`p^{}`$. The four components of $`X`$ are integration variables and, therefore, there is much freedom in their interpretation. It is usual, however, to identify $`S(X,K)`$ with the space-time distribution of the sources producing the pions with momentum $`K`$ and then $`X`$ is interpreted as the space-time position of the source. A source functions combines conveniently the information and/or prejudice about the space-time ($`X`$) distribution of the sources and about the momentum distribution of the produced particles. It yields the density matrix, which can be used to get the (two- as well as more particle) correlation function. Comparing the predicted correlation functions with experiment one can improve the model used to find $`S(X,K)`$, fix its parameters etc. Note that the density matrix on the left hand side of formula (1) does not depend on time. Since after freezeout the particles propagate freely (we do not discuss the final state interactions here), this means that for times after the freeze out period the elements of the density matrix should be interpreted as the matrix elements of the time dependent (Schrödinger picture) density operator between the time dependent states $`|\text{p},t=e^{iE(p)t}|\text{p}`$, or equivalently as the matrix element of the time independent (interaction picture) density operator between the time independent states $`|\text{p}`$.

Formula (1) looks similar to the formula relating the Wigner function $`W(\text{X},\text{K})`$ (cf. e.g.  and references contained there) to the density matrix:

$$\rho (\text{p},\text{p}^{},t)=d^3Xe^{i\text{q}\text{X}}W(\text{X},\text{K},t).$$
(2)
In fact the building block of the source function was originally a Wigner function (cf.  formula (7)). Because of the similarity between formulae (1) and (2), the source function is often referred to as a Wigner function, a kind of Wigner function, a pseudo-Wigner function etc. (cf. e.g. ref. , from which all these names have been taken). A discussion of the actual relation between the source function and the Wigner functions, however, seems to be missing in the literature. Let us begin with some general remarks.

Formula (2) contains the time dependent density matrix, i.e. the time dependent density operator in the representation of the time independent states $`|\text{p}`$ . Thus

$$\rho (\text{p},\text{p}^{},t)=e^{iq_0t}\rho (\text{p},\text{p}^{}).$$
(3)
Relation (1) follows from (2), if we put

$$S(X,K)=W(𝐗,𝐊,t_0)\delta (X_0t_0).$$
(4)
This choice can be always made and has a simple physical interpretation. Choose a time $`t_0`$ after freezeout, when all the pions are already present, but none have been registered yet. The state at $`t_0`$ yields the initial condition from which any distribution at registration can be calculated. Assuming that all the pions appeared simultaneously at $`t_0`$ may be poor physics, but this does not change the fact that the predictions at registration time are correct. There is an infinite choice of other source functions, which also satisfy relation (1). Since this relation is the only link between a source function and physical reality, there can be no objection of principle against the choice (4). From the practical point of view, however, another source function, which can be calculated from some model, may be more convenient and, of course, it is just as good, if it reproduces the density matrix equally well. Note that the identity of functions in (4) does not mean that their argument have the same physical interpretation. In fact, on the left hand side $`X`$ and $`K`$ denote the position of the source in space-time and the four-momentum of the pion, while on the right hand side $`\text{X}=\frac{1}{2}(\text{x}+\text{x}^{})`$ and $`\text{K}=\frac{1}{2}(\text{p}+\text{p}^{})`$ are defined by reference to the arguments of the density matrices of the pion in the coordinate representation and in the momentum representation respectively, while $`t_0`$ is the time. It is often reasonable to assume that X and K on both sides are good approximations to the position at $`t=t_0`$ and momentum for $`tt_0`$ of the pion, but this is certainly not always the case.

Let us conclude this section with some more remarks. If one assumes that all the particles were created simultaneously at some common freezeout time $`t=t_0`$, then a ”realistic” source function should be proportional to $`\delta (X_0t_0)`$ and the proportionality coefficient is the Wigner function as in (4). The source function cannot be equal to a Wigner function, because the dimensions are different. The delta function in formula (4) brings into the dimension the necessary factor 1/time. What could be equal to the source function is what we will refer to as the differential Wigner function defined by the relation

$$\stackrel{~}{W}(\text{X},\text{K},t_0)=W_{dt_0}(\text{X},\text{K},t_0)/dt_0.$$
(5)
Here $`W_{dt_0}`$ is the Wigner function at time $`t>t_0`$ for the particles from sources created in the time interval $`dt_0`$ around time $`t_0`$. Note that we are using the Wigner functions normalized to the numbers of particles and not to one. When the amplitudes of the pions produced at different times add incoherently, formula (2) can be rewritten in terms of the differential Wigner function $`\stackrel{~}{W}`$ as:

$$\rho (p,p^{})=d^4Xe^{iqX}\stackrel{~}{W}(X,\text{K}),$$
(6)
where the fourth component of the vector $`X=(\text{X},X_0)`$ is $`t_0`$. Comparing this formula with formula (1) we find a particular solution for the source function

$$S(X,K)=\stackrel{~}{W}(X,\text{K}).$$
(7)
We stress that this is only one out of the infinity of source functions, which when substituted into formula (1) yield the correct density matrix and that for the sources incoherence in the creation time has been assumed.

## 2 GGLP models

In the seminal paper of G. Goldhaber, S. Goldhaber, W. Lee and A. Pais  the assumption was made that pion production amplitudes at different space points are incoherent. In this simple case there is no need to distinguish between the sources and the pions at $`t=t_0`$. A natural generalization (cf. e.g. ) was to assume more generally that pions production amplitudes at different space-time points are incoherent. Further we refer to such models as GGLP models.

Let us begin with the simpler case, when all the pions appear simultaneously at $`t=t_0`$. For $`t>t_0`$ they are assumed to propagate freely. The time independent density matrix for $`t>t_0`$ is:

$$\rho (\text{p},\text{p}^{})=d^3x\text{p},t|\text{x},t;t_0\rho (\text{x},t_0)\text{x},t;t_0|\text{p}^{}t.$$
(8)
Here $`|\text{x},t;t_0=\mathrm{exp}[i(tt_0)\widehat{H_0}]|𝐱`$ is the time dependent state, evolving according to the free particle Hamiltonian $`\widehat{H_0}`$ from the initial state $`|\text{x}`$, which corresponds to a pion localized at $`𝐱`$ at time $`t_0`$. One could also say that in space-time the position of the source is $`(\text{x},t_0)`$. Let us note for further use the formula

$$\text{p},t|\text{x},t;t_0=e^{iE(p)t_0}\text{p}|\text{x},$$
(9)
where $`E(p)=\sqrt{\text{p}^2+m_\pi ^2}`$. Formula (8) is perfectly acceptable from the point of view of quantum mechanics, in spite of the fact that according to the picture described in ref.  a pion with momentum $`𝐩`$ is produced at point x, which would contradict the uncertainty principle. Since the density matrix is diagonal in x (in the x, which labels the incoherent states, not in position at time $`t`$) and $`t=t^{}=t_0`$ , we can identify the components of the vector x with $`X_1,X_2,X_3`$ and the parameter $`t_0`$ with $`X_0`$. The matrix element (8) can be rewritten as

$$\rho (\text{p},\text{p}^{})=\frac{d^3X}{(2\pi )^3}\mathrm{exp}(iqX)\rho (X).$$
(10)
Using (3) and comparing with (2) this yields

$$W(X,\text{K})=\frac{\rho (X)}{(2\pi )^3}.$$
(11)
The source function is given by formula (4), because by assumption all the particle production happens at $`t=t_0`$. The components of $`X`$ are the three components of the position vector of the point where the particle was produced and its time of production (freezeout time). Since after freezeout the momentum vector of each particle is conserved separately, the momenta entering K can be taken at any time not earlier than $`t_0`$. The actual time of observation does not occur in the formula at all.

Note that for this class of models the Wigner function does not depend on momenta. With the present interpretation this would lead to disaster when comparing with experiment. The Goldhabers Lee and Pais  proceeded differently. They were interested in $`n`$-pion exclusive states for $`n=4,5,6`$ with no more than two pions of a given charge. Using matrix (8) they constructed the probability distributions $`P(p_1,\mathrm{},p_n)`$. These were interpreted as momentum distributions only after the sets of momenta contradicting energy and/or momentum conservation were projected out. In practice they interpreted the $`P`$-s as integrands of the corresponding phase space integrals. This could not be generalized to inclusive processes, led to rather complicated integrations and gave results, which at that time were considered acceptable, but which soon were shown to be in violent disagreement with experiment . Most of these difficulties were overcome by Kopylov, Podgoretsky and coworkers as described in the next section.

Let us consider now the more general case, when the pions are produced at different moments of time, but there is no interference between the amplitudes corresponding to production at different moments. Then at times larger than the largest freezeout time the density matrix is given by an integral over the freezeout time $`t_0`$ of the (differential version of the) density matrix (10) corresponding to particles produced at given time $`t_0`$:

$$\rho (𝐩,𝐩^{})=\frac{d^4X}{(2\pi )^3}\mathrm{exp}(iqX)\stackrel{~}{\rho }(X).$$
(12)
Note that in order to keep the dimensions right $`\rho (X)`$, which is a density in space, has been replaced by $`\stackrel{~}{\rho }(X)`$ which is a density in space time. In fact, $`\rho (X)`$ are diagonal elements of the time independent single pion density matrix in the $`|\text{x},t;t_0`$ representation, while $`\stackrel{~}{\rho }(X)`$ are the diagonal elements of the corresponding differential density matrix related to $`\rho (X)`$ by a formula analogous to (5). Comparing with (1) it is seen that one solution for $`S(X,K)`$ is

$$S(X,K)=\frac{\stackrel{~}{\rho }(X)}{(2\pi )^3}.$$
(13)
Since there is incoherence in time, the relation between the source function and Wigner’s function is $`S(X,K)=\stackrel{~}{W}(X,\text{K})`$ (7).

This seems to be the perfect situation: the source function reproduces exactly the distribution of (incoherent) sources in space-time and coincides with the differential Wigner function. However, (13) is only one of the infinity of source functions, which yield the same density matrix. A quite different one can be obtained e.g. by calculating the Wigner function corresponding to the density matrix $`\rho (\text{p}_1,\text{p}_2,t)`$ and substituting it into formula (4). We can conclude only that for this class of (unrealistic) models there is an acceptable source function, which can be interpreted as a differential Wigner function. Therefore, it can be argued that this particular source function is easier to find than others.

## 3 Including coherence in space

An important extension of the GGLP approach is to assume that at any given creation time there may be coherence between the amplitudes of the pions produced at different points in space, though pions produced at different times $`t_0`$ do not interfere. Models of this type, also in the more general version with coherence in time, have been introduced by Kopylov and Podgoretsky . In terms of sources this may be quite complicated. For instance, in the model described by Kopylov and Podgoretsky  all the sources are produced simultaneously and live for some time. The averaging over the energy of each source (not of the pion!), however, kills the interference between the amplitudes for the production of the pion by the same source at different moments of time. Therefore one can replace each original source by a set of sources distributed in time and producing the pions instantly. These are not quite satisfactory models, because it is difficult to explain why the creation amplitudes from neighboring time moments should not interfere, while those from neighboring space points do. Moreover, relativistically the creation processes which are at different times in one Lorentz frame may be simultaneous in another Lorentz frame. As seen from formula (4), however, even in the simplest case of simultaneous production, such models can reproduce correctly any single particle density matrix and consequently any single particle momentum distribution without doing phase space integrations. Thus they can be applied to describe inclusive processes.

Including interference in space and summing without interference over time corresponds to the following generalization of formula (8)

$$\rho (\text{p},\text{p}^{})=𝑑t_0d^3xd^3x^{}\text{p},t|\text{x},t;t_0\stackrel{~}{\rho }(\text{x},\text{x}^{},t_0)\text{x}^{},t;t_0|\text{p}^{},t.$$
(14)
There are many equivalent ways of rewriting this relation. For instance, one may use a representation $`|\alpha `$, where the (differential) density matrix $`\stackrel{~}{\rho }(\text{x},\text{x}^{},t_0)`$ is diagonal. Then (14) gets replaced by

$$\rho (\text{p},\text{p}^{})=𝑑t_0𝑑\alpha e^{iq_0t_0}\text{p}|\alpha \stackrel{~}{\rho }(\alpha ;t_0)\alpha |\text{p}^{}.$$
(15)
It is instructive to compare this model, to the GGLP models. At $`t_0`$ the localized states $`|\text{x}=|\text{x},t_0;t_0`$ have been replaced by the states $`|\alpha =|\alpha ,t_0;t_0`$. The plane waves $`\text{p}|\text{x}`$ have been replaced by the functions $`\text{p}|\alpha `$, which are known under a variety of names: as waved packets, as sources, or as currents. The source function for such models was introduced by Pratt . His statement that333This has been changed to our notation. Pratt has written $`g,\stackrel{}{p},x`$ where we have written $`S,\text{K},X`$$`S(X,\text{K})`$ can be identified as the probability of emitting a pion of momentum K from space-time point X” is, however, only an approximation.

Kopylov and Podgoretsky  assumed that $`(\alpha ,t_0)`$ is a point in space-time, so that formula (15) can be rewritten as

$$\rho (\text{p},\text{p}^{})=d^4x_0\text{p}|\psi _{x_0}\stackrel{~}{\rho }(x_0)\psi _{x_0}|\text{p}^{}e^{iq_0t_0}.$$
(16)
This approach yields an alternative method of describing the distribution of pions in phase space. In the integrand $`\text{p}|\psi _{x_0}`$ is the probability amplitude for finding a pion with momentum p produced by a source labelled $`x_0`$. When the states $`|\psi _{x_0}`$ correspond to particle well localized in space, p and $`\text{x}_0`$ give a reasonably good description of the position of the pion in phase space. One could object that this is only a rough description, but the same is true for the Wigner function: K and X give only approximately the momentum and position of the pion. An exact determination of the pion position in phase space is possible only in classical physics.

In ref.  the states $`|\alpha =|\psi _{x_0}`$ were supposed to be related by space time translations so that

$`\psi _{x_0}(x)`$ $``$ $`\text{x}|\psi _{x_0}=\psi (xx_0),`$ (17)
$`\text{p}|\psi _{x_0}`$ $`=`$ $`e^{ipx_0}\varphi (\text{p}),\varphi (\text{p})=\text{p}|\psi _0.`$ (18)
The second formula is, of course, equivalent to the first. Another choice has been made in the ”covariant current ensemble formalism” , , , . There each source is labelled by a position in space-time $`x_0`$ and a four-momentum $`p_0`$. Usually $`\text{x}_0`$ and $`\text{p}_0`$ denote the centers of the wave packet in ordinary space and in momentum space respectively. The time component of $`x_0`$ is $`t_0`$ – the freeze out time of the wave packet. The time component of $`p_0`$ is calculated from the condition $`p_0^2=m^2`$, where $`m`$ can , but does not have to  be the pion mass. The nice feature of this parametrization of the wave packets is that one can substitute a classical trajectory $`\text{x}_0(t),\text{p}_0(t)`$ for the source and remain in agreement with the Heisenberg uncertainty principle for the pions. The choice for the scalar products is

$$\text{p}|\psi _{x_0,p_0}=e^{ipx_0}j(\frac{pp_0}{m}).$$
(19)
Thus, the sources are related by space-time translations and when each of the currents $`j`$ is considered in its rest frame where $`p_0=(m,0)`$, they are identical. Formulae (17) and (19) have been applied also to relativistically covariant models cf. e.g. , .

Another way of rewriting relation (14) is

$$\rho (𝐩,𝐩^{})=d^4X\mathrm{exp}(iqX)\frac{d^3y}{(2\pi )^3}\mathrm{exp}(i𝐊\text{y})\stackrel{~}{\rho }(𝐱,𝐱^{},t_0),$$
(20)
where $`𝐗=(𝐱+𝐱^{})/2`$, $`𝐲=𝐱𝐱^{}`$ and $`X_0=t_0`$. According to formula (1) a possible choice of the source function is

$$S(X,K)=\frac{d^3y}{(2\pi )^3}\mathrm{exp}(i\text{K}\text{y})\stackrel{~}{\rho }(\text{x},\text{x}^{},t_0).$$
(21)
The differential Wigner function is the Wigner transform of the differential density matrix as it should. The relation of this particular source function to the Wigner function is again given by formula (7). We conclude that among the infinitely many source functions which give the same density matrix in the momentum representation there is one, which can be related to the Wigner function as described. Note that this source function depends on K, but does not depend on $`K_0`$.

## 4 Including coherence in space and time

Finally let us consider the case, when neither coherence in space, nor coherence in time is assumed. Then the density matrix is

$$\rho (\text{p},\text{p}^{})=d^4xd^4x^{}\text{p},t|\text{x},t;t_0\stackrel{~}{\rho }(x,x^{})\text{x}^{},t;t_0^{}|\text{p}^{},t,$$
(22)
where $`x=(\text{x},t_0)`$ and $`x^{}=(\text{x}^{},t_0^{})`$. For an application of a formula of this type see e.g. ref. , where $`_FM_F(x)M_F(x^{})`$ stand for our $`\stackrel{~}{\rho }(x,x^{})`$.

Formula (22) can be rewritten as

$$\rho (\text{p},\text{p}^{})=d^4Xd^4y\mathrm{exp}[iqX+iKy]\stackrel{~}{\rho }(x,x^{})/(2\pi )^3,$$
(23)
where $`X=\frac{1}{2}(x+x^{})`$ and $`y=xx^{}`$. Comparison with formula (1) gives as one of the solutions for the source function

$$S(X,K)=\frac{d^4y}{(2\pi )^3}\mathrm{exp}[iKy]\stackrel{~}{\rho }(x,x^{}).$$
(24)
The right hand side is rather remote from what one usually calls a Wigner function. The differential Wigner function was at least proportional to a Wigner function, though only for the particles from sources produced in the time interval $`dt_0`$ around the time $`t=t_0`$. Function (24) is proportional to the contribution to the Wigner function from the interference of the production amplitude in the time interval $`dt_0`$ around $`t=t_0`$ with the production amplitude in the time interval $`dt_0^{}`$ around $`t=t_0^{}`$, integrated over $`t_0^{}t_0`$ at fixed $`t_0^{}+t_0`$.

Formula (22) rewritten in terms of wave packets reads

$$\rho (\text{p},\text{p}^{})=𝑑t_0𝑑t_0^{}𝑑\alpha \text{p},t|\alpha ,t;t_0\stackrel{~}{\rho }(\alpha ,t_0,t_0^{})\alpha ,t;t_0^{}|\text{p}^{},t.$$
(25)
Note that the state $`|\alpha ,t;t_0`$ may, but does not have to, be connected to the state $`|\alpha ,t;t_0^{}`$ by a smooth hamiltonian evolution. The source at time $`t_0`$ can be, as well, something quite different from the source at time $`t_0^{}`$. On the other hand, the evolution of $`|\alpha ,t;t_0`$ in the time $`t`$ for $`t`$ later than the freezeout of the particle is the ordinary free particle evolution. The various models are defined by the choice of $`\alpha `$, of the states $`|\alpha ,t;t_0`$ and of the weight function $`\stackrel{~}{\rho }(\alpha ,t_0,t_0^{})`$. Examples can be found in refs. , , .

## 5 Conclusions

A source function cannot be equal to a Wigner function, because they have different dimensions. Moreover, for a given state of the system its Wigner function is well defined, while its source function is not. The problem is, therefore, to chose some special source function and try to relate it somehow to a Wigner function.

When all the particles are produced simultaneously and when this is assumed to mean that the source function is proportional to a delta function in time, the proportionality coefficient is unambiguously defined as the Wigner function of the pions at the time of freeze out (4). When it is assumed that the production amplitudes at different times add incoherently, one can use the source function proportional to the differential Wigner function as given in (13). When the production amplitudes at different times interfere, a source function can be related to a piece of the Wigner function as given by formula (24) and explained below this formula. In this case, however, the use of wave packets (or sources, or currents) may be a more natural way to analyze the phase space distribution of pions.

If, as is often done, a model is defined by postulating the source function (cf.  for examples), the question about the relation of this source function to a Wigner function cannot be answered without making an assumption about the coherence or incoherence of the amplitudes for pion production at different times. When the source function is proportional to a delta function in time, one can relate it to a Wigner function by relation (4). When it does not depend on $`K_0`$, one can use formula (7). In the later case there is no guarantee that the production process was such as the differential Wigner functions suggest.


text:
# On the (absence of a) relationship between bound and scattering states in quantum mechanics. Application to 12C + 𝛼

## I Introduction

Can bound-state properties be deduced from scattering data? The answer to this question showed interesting evolutions throughout the history of quantum scattering theory (see Foreword by R. G. Newton in Ref. Chadan and Sabatier (1977)), starting with Heisenberg’s conjecture in the early forties that the scattering (S) matrix contains all physical information, including that regarding bound states. This conjecture was later seriously weakened in the framework of the potential model: in 1949, Bargmann Bargmann (1949a, b) constructed potentials that share the same S matrix for partial wave $`\mathrm{}=0`$ (so-called phase-equivalent potentials because they share the same scattering phase shifts) but have different bound-state properties \[different energies, or identical energies but different wave functions\]. On the other hand, that same year, Levinson (see for instance Ref. Chadan and Sabatier (1977)) proved that, for a given partial wave $`\mathrm{}`$, the number of bound states, $`N_{\mathrm{}}`$, can be deduced from the phase shifts. These results were then explained by the inverse-problem theory Chadan and Sabatier (1977), which shows that for a given partial wave the potential is uniquely determined by (i) the scattering phase shifts at all energies, (ii) the bound-state energies and (iii) one additional real parameter for each bound-state wave function \[e.g., the asymptotic normalization constant (ANC), see definition (1)\]. This only applies to local central potentials displaying no singularity. Levinson’s theorem, however, was generalized by Swan in 1963 Swan (1963), who proved that $`N_{\mathrm{}}`$ cannot be deduced from the S matrix for potentials displaying an $`r^2`$ repulsive singularity at the origin. Since then, the construction of phase-equivalent potentials for a given partial wave has known important progress Baye (1987a, b); Ancarani and Baye (1992); Baye (1993); Baye and Sparenberg (1994) thanks to the use of the algebraic formalism of supersymmetric quantum mechanics Witten (1981); Sukumar (1985a, b) used below (see Fig. 2).

As far as only one partial wave is concerned, the answer to the above question is thus rather clear-cut: neither the number of bound states, nor their energies nor their ANCs can be deduced from the S matrix alone. Moreover, these bound-state properties are independent of one another, which means that knowing, for instance, both the S matrix and a bound-state energy does not constrain the value of the ANC as long as no other physical information is available on this particular bound state (e.g., its radius or its lifetime) or on the considered system in general (e.g., the range of the interaction). This impossibility proof being established in the framework of the potential model, which is a particular case of more general models (e.g., microscopic or phenomenological models), it holds for these more general models too. However, several works Iwinski et al. (1984); Blokhintsev et al. (1984, 1993); Angulo and Descouvemont (2000); Tischhauser et al. (2002) seem to ignore this argument and aim at extracting the ANC of a bound-state with known energy from the corresponding partial-wave S matrix. In Refs. Angulo and Descouvemont (2000); Tischhauser et al. (2002) for instance, the reaction- (R-)matrix formalism Wigner and Eisenbund (1947); Lane and Thomas (1958) is used to extract the ANC of an 16O bound state from 12C + $`\alpha `$ elastic-scattering data. In the R-matrix phenomenological model, it is rather natural to assume that bound states have a measurable impact on the partial-wave S matrix, like elastic-scattering resonances, because both bound and resonant states are characterized by an energy and a reduced width. However, in scattering theory, bound and resonant states have very different natures, and the above impossibility argument seriously questions the ANC extraction of Refs. Angulo and Descouvemont (2000); Tischhauser et al. (2002). In Sec. II, I show how that misuse of the R-matrix formalism actually explains the inconsistency between the ANC values obtained in these references. Then in Sec. III, I explicitly construct simple potential models which illustrate the independence of the bound-state ANC from the phase shifts for a single partial wave.

Now, considering several partial waves simultaneously leads to a less clear situation. The inverse-problem theory implies that a regular interaction potential is uniquely determined by phase shifts of a given partial wave with angular momentum high enough so that the partial wave contains no bound state. Using inversion techniques Chadan and Sabatier (1977), a unique potential can thus in principle be deduced from such phase shifts; this potential can in turn be used for lower partial waves, for which it may sustain bound states. In that case, bound state properties (of low partial waves) are totally determined by scattering phase shifts (of a high partial wave). However, this program presents three difficulties. First, it assumes that the potential is local and energy independent, which, in general, is not a realistic approximation when the interacting particles have an internal structure (see the discussion in Secs. IV and V for the particular case of 12C + $`\alpha `$). Second, it requires that the interaction be independent of the partial wave, which is not always the case: scattering phase shifts of several partial waves are not necessarily compatible with one unique potential. Third, it is based on the knowledge of the phase shifts of a high-angular-momentum partial wave at all energies, which are generally not available: at high energy, the complexity of multichannel effects increases, which makes the data difficult to analyze; and at low energy, the phase shifts of high-angular-momentum partial waves are very small because of the centrifugal barrier. Hence, considering several partial waves at the same time may help deducing bound-state properties from scattering data, but the practical applicability of the method could strongly depend on the considered system. In Sec. IV, I show that for the particular case of 12C + $`\alpha `$ encouraging results are obtained by considering several partial waves at the same time, in the framework of a simplified model in which 12C and $`\alpha `$ are considered as rigid nuclei in their $`0^+`$ ground state, interacting through a local energy- and angular-momentum-independent potential. These results are then compared with more sophisticated models available in the literature. Conclusions and perspectives for the 12C + $`\alpha `$ system are given in Sec. V.

## II  R-matrix analyses

The 12C + $`\alpha `$ system is particularly important in nuclear astrophysics since the 12C ($`\alpha `$, $`\gamma `$) 16O capture is a key reaction in the helium-burning phase of red giant stars: by competing with the triple $`\alpha `$ reaction (leading to 12C), it determines the ratio of carbon and oxygen which results from stellar nucleosynthesis. Experimentally, the Coulomb repulsion between 12C and $`\alpha `$ makes the direct measurement of the capture reaction impossible at the very low energies of astrophysical interest (around 300 keV in the center-of-mass frame). Hence, available experimental results at higher energies ($`E_{\mathrm{c}.\mathrm{m}.}>1`$ MeV) have to be extrapolated to these low energies. This extrapolation, generally performed with the R-matrix model (the only model to date able to precisely fit all the available data), is itself problematic because of the presence of two 16O bound states just below the 12C + $`\alpha `$ threshold. Whereas the influence of the $`E_{\mathrm{c}.\mathrm{m}.}=45.1\pm 0.1`$ keV state is well understood Azuma et al. (1994); Brune et al. (1999), the influence of the $`E_{\mathrm{c}.\mathrm{m}.}=244.9\pm 0.6`$ keV state is still the major source of uncertainty on the reaction rate Brune et al. (1999); Angulo and Descouvemont (2000); Tischhauser et al. (2002). The unknown key quantity is its ANC (the energies of both states are well known), which is related to its $`\alpha `$ reduced width in the R-matrix formalism (I do not use reduced widths here since they depend on the R-matrix radius, which makes comparison between different models difficult).

The $`245`$ keV subthreshold state has a positive parity and a $`J=2`$ angular momentum ($`2^+`$ state); microscopic cluster models (see Ref. Descouvemont (1987) and references therein) suggest that its structure is dominated by a 12C cluster and an $`\alpha `$ cluster, both in their fundamental $`0^+`$ state, with a relative motion of angular momentum $`\mathrm{}=2`$. It is thus on the phase shifts of this partial wave that one could naively expect to see the influence of the $`2^+`$ state, and it is from these phase shifts that one could try to extract its ANC, $`C_{12}`$ (where the second index refers to $`\mathrm{}=2`$ and the first index means that this state is the lowest one for this partial wave). This constant is defined by the asymptotic behavior of the radial part, $`u_{12}(r)`$, of the relative wave function between the 12C and $`\alpha `$ nuclei for this particular bound state:

$`u_{12}(r)\underset{r\mathrm{}}{}C_{12}W_{6.615,5/2}(0.3751r),`$ (1)
where $`r`$ is the radial coordinate (in fm) and $`W`$ is the Whittaker function Abramowitz and Stegun (1965) (a decreasing exponential function in the absence of Coulomb interaction).

In Ref. Angulo and Descouvemont (2000), the $`\mathrm{}=2`$ phase shifts from Ref. Plaga et al. (1987) are analyzed by the R-matrix model, leading to a minimum in the $`\chi ^2`$ (see Fig. 1 of Ref. Angulo and Descouvemont (2000)) for $`C_{12}=402\times 10^3`$ fm-1/2. However, this extraction is very delicate since the influence of the subthreshold state has to be disentangled from the background phase shifts in the energy region where experimental data are available; this influence is smaller (see Fig. 2 of Ref. Angulo and Descouvemont (2000)) than the main components of this background, i.e., the hard-sphere and high-energy-pole phase shifts. As a consequence, this minimum strongly depends on the range of data fitted (see Fig. 1 of Ref. Angulo and Descouvemont (2000)). Moreover, I have checked that the minimum is highly sensitive to the high-energy background-pole energy, which is fixed at 10 MeV in Ref. Angulo and Descouvemont (2000). Without this restriction, very different values of $`C_{12}`$ are obtained; hence, the error on $`C_{12}`$, which is not estimated in Ref. Angulo and Descouvemont (2000), is actually very large.

In Ref. Tischhauser et al. (2002), the 12C + $`\alpha `$ elastic scattering has been remeasured with very high precision. An R-matrix analysis with free high-energy-pole energies leads to a $`\chi ^2`$ minimum (see Fig. 2 (a) of Ref. Tischhauser et al. (2002)) for $`C_{12}=154\pm 18\times 10^3`$ fm-1/2. The error bar on $`C_{12}`$ is estimated by following the guideline Press et al. (1989); Azuma et al. (1994)

$`\chi ^2<\chi _{\mathrm{min}}^2+9\chi _{\mathrm{min}}^2/\nu ,`$ (2)
where $`\chi _{\mathrm{min}}^2`$ is the minimal $`\chi ^2`$ and $`\nu `$ is the number of degrees of freedom. In the present case, one has $`\chi _{\mathrm{min}}^2=18941`$ (the $`\chi ^2`$ of Fig. 2 in Ref. Tischhauser et al. (2002) have to be divided by the 32 reference angles), for 11392 data points and 65 parameters (32 R-matrix parameters and 33 experimental parameters), which leads to $`\nu =11327`$. Following the theory of $`\chi ^2`$ minimization Press et al. (1989), one finds that the probability of having $`\chi ^218941`$ for $`\nu =11327`$ is of the order of $`10^{391}`$, which probably indicates that some systematic errors have been underestimated in Ref. Tischhauser et al. (2002). Consequently, Eq. (2) cannot be used to define error bars on parameters since this equation only takes into account statistical errors and requires that the fit be good in the sense of $`\chi ^2`$ theory. A less rigorous way of checking the validity of a fit is to use the “chi-by-eye” approach Press et al. (1989). In the case of Ref. Tischhauser et al. (2002), the fit looks good (see Fig. 1 of this reference), which is related to the fact that $`\chi _{\mathrm{min}}^2/\nu =1.672`$ is not that far from 1. However, the fits obtained with $`C_{12}=136`$ or $`172\times 10^3`$ fm-1/2, which correspond to $`\chi ^2/\nu =1.674`$, look probably equally good. I thus consider that the error bar on $`C_{12}`$ obtained in Ref. Tischhauser et al. (2002) is underestimated.

In summary, the values of $`C_{12}`$ obtained in both Refs. Angulo and Descouvemont (2000) and Tischhauser et al. (2002) are most probably unreliable, for different reasons. This explains why, despite the fact that the $`\mathrm{}=2`$ phase shifts of Refs. Plaga et al. (1987) and Tischhauser et al. (2002); buchmann:03  are compatible with one another, their R-matrix analyses give totally incompatible results for $`C_{12}`$. One might think that this inconsistency is due to the fact that in Ref. Tischhauser et al. (2002) all relevant partial waves are taken into account whereas in Ref. Angulo and Descouvemont (2000) only $`\mathrm{}=2`$ is considered. However, since the R-matrix formalism is based on a partial-wave decomposition of the cross section, $`C_{12}`$ is actually independent of the $`\mathrm{}2`$ parameters (in Ref. Tischhauser et al. (2002), the only influence of the $`\mathrm{}2`$ phase shifts on $`\mathrm{}=2`$ comes through the choice of a common R-matrix radius for all partial waves; however, this does not seem to strongly constrain the fit of the $`\mathrm{}=2`$ phase shifts, as discussed in the next paragraph). I argue that this incompatibility rather illustrates the impossibility, well-known in quantum scattering theory, of deducing a bound-state ANC when only the bound-state energy and the elastic-scattering phase shifts of the corresponding partial wave are known.

Let us however notice that, in the case of the R-matrix phenomenological model, this impossibility is probably not as strong as in general scattering theory. As shown in the next Section (see Fig. 2), the range of a 12C + $`\alpha `$ interaction that reproduces both the $`\mathrm{}=2`$ phase shifts and the bound-state energy increases with the value of $`C_{12}`$. But the R-matrix formalism is only valid when the nuclear interaction is negligible above the R-matrix radius. In the R-matrix theory, where the number of R-matrix poles is infinite Lane and Thomas (1958), this does not cause any problem: the results are independent of the R-matrix radius, which can thus be chosen as large as necessary (in other words, it would be possible to calculate an R-matrix for any of the potentials plotted in Fig. 2 but the R-matrix radius should be larger than the potential range). In the R-matrix phenomenological model, where the number of poles is limited to the number of bound states and resonances in the energy interval considered, plus one possible high-energy pole to help simulating the background, a constraint appears on the R-matrix radius: in Ref. Vogt (1996), a value of 5.5-6 fm is found, based on the $`\mathrm{}=1`$ partial wave; in Ref. Tischhauser et al. (2002) \[see Fig. 2 (b)\] 5.5 fm leads to the best fit (I have checked that this constraint is mainly due to the $`\mathrm{}=0`$ partial wave, for which the background is maximal); in Ref. Angulo and Descouvemont (2000), the radius is fixed at 6.5 fm. Hence, very large values of $`C_{12}`$ can probably be excluded by an R-matrix fit; very small values, on the other hand, cannot be excluded. I could actually find very good fits of the $`\mathrm{}=2`$ phase shifts with $`C_{12}=0`$ for R-matrix radii in the range 4.5-6.5 fm. This confirms the doubts expressed above on the analyses of Refs. Angulo and Descouvemont (2000); Tischhauser et al. (2002) and it shows that, in this particular case, the quality of the fit does not strongly depend on the precise value of the R-matrix radius.

## III  Inversion of one partial wave

Up to now, the discussion is based on general arguments from scattering theory. Let me now construct simplified two-body potential models that illustrate these arguments in the particular case of 12C + $`\alpha `$, using the supersymmetric-quantum-mechanics Witten (1981); Sukumar (1985a, b) inversion technique of Refs. Sparenberg and Baye (1997); Sparenberg (2000). Starting from a phase-shift analysis buchmann:03  of the purely elastic data ($`E_{\mathrm{c}.\mathrm{m}.}<4.966`$ MeV) of Ref. Tischhauser et al. (2002), I have removed, using the R-matrix formalism, resonances that have a more complicated structure than 12C + $`\alpha `$ in their $`0^+`$ ground states (such resonances cannot be explained by the present simplified model). For $`\mathrm{}=2`$, this means removing the narrow resonance at 2.683 MeV and the wide resonance at 4.358 MeV Tilley et al. (1993), which leads to the phase shifts plotted in Fig. 1 (the spread of the data above 4.2 MeV is due to the difficulty of performing a reliable phase-shift analysis in the vicinity of a resonance). By solving the so-called “singular inverse problem” defined in Ref. Sparenberg and Baye (1997), I have constructed a potential that reproduces these phase shifts with very high precision. This potential has an $`r^2`$ singularity at the origin and satisfies Swan’s theorem Swan (1963); to simplify further developments, its nuclear part has been approximated by the analytical expression

$`V(r)=V\mathrm{exp}(r^2/R^2)/r^2,`$ (3)
where $`V=43.40`$ MeV fm2 and $`R=5.091`$ fm. The corresponding $`\mathrm{}=2`$ effective potential (centrifugal barrier plus Coulomb and nuclear interactions) is represented in Fig. 2 (solid line) and its phase shifts are shown in Fig. 1. Surprisingly, a potential of the form (3) is also able to reproduce the $`\mathrm{}=0`$ background phase shifts with satisfactory precision, which suggests an alternative way of parameterizing the 12C + $`\alpha `$ background in future R-matrix analyses (in the spirit of the “hybrid” R-matrix model of Refs. Johnson (1973); Koonin et al. (1974); Langanke and Koonin (1985)). The fit of the $`\mathrm{}=2`$ scattering phase shifts is very good, despite the fact that potential (3) has, by construction, no bound state (it is purely repulsive); this is a first proof that, from the point of view of scattering theory, the $`\mathrm{}=2`$ background phase shifts are totally disconnected from the $`2^+`$ subthreshold bound state.

Now, since this state is known to exist and to have a structure such that it could be well simulated by the present simplified model, it is physically meaningful to perform a further step in the inversion of the $`\mathrm{}=2`$ phase shifts, i.e., to add a bound state to the above potential without modifying its phase shifts Sparenberg (2000). This phase-equivalent bound-state addition may be carried out through a further use of supersymmetric quantum mechanics Baye (1987b); it introduces two arbitrary parameters: the bound-state energy and its ANC. In the present case, the energy is known experimentally but not the ANC; hence, there is a one-parameter family of $`\mathrm{}=2`$ effective potentials that share identical phase shifts (Fig. 1) and have a bound-state at $`245`$ keV. Such potentials are represented by dashed lines for three arbitrary values of $`C_{12}`$ in Fig. 2; they constitute generalizations (to $`\mathrm{}0`$ and in the presence of a Coulomb interaction) of Bargmann’s potentials (see Fig. 1 in Ref. Bargmann (1949b)). The corresponding bound-state wave functions are plotted in Fig. 3 (dashed lines), where their very different asymptotic behaviors are clearly seen. This proves that, from the $`\mathrm{}=2`$ phase shifts and bound-state energy only, there is no way of telling what is the best value for $`C_{12}`$, in contrast with R-matrix attempts described above (see however the discussion of the R-matrix phenomenological model in the previous Section).

## IV  Inversion of several partial waves

The general discussion from the introduction suggests to look at other partial waves simultaneously. The $`2^+`$ subthreshold state is actually a member of a rotational band, together with the $`0^+`$ bound state at $`\mathrm{1\hspace{0.17em}113}\pm 1`$ keV and the $`4^+`$ and $`6^+`$ resonances at $`\mathrm{3\hspace{0.17em}194}\pm 3`$ keV (width: $`26\pm 3`$ keV) and $`\mathrm{9\hspace{0.17em}113}\pm 7`$ keV (width: $`420\pm 20`$ keV) respectively Tilley et al. (1993). This is indicated by cluster-model calculations (see Ref. Descouvemont (1987) and references therein), which give a fair description of these four states in a single-channel calculation where both the 12C and $`\alpha `$ nuclei are in their $`0^+`$ ground state. This suggests to ask the following question: could all these states be reproduced by a unique $`\mathrm{}`$-independent potential? Of course, in principle, the internal structure of the nuclei could play a role and a satisfactory potential could be non local and energy dependent. It is however worth testing the simplest hypothesis, i.e., that of an energy- and $`\mathrm{}`$-independent local potential, because if such a potential were sufficient, it would be strongly constrained by the data of the four partial waves considered simultaneously, which could help solving the $`C_{12}`$ indefiniteness.

I could not find such a potential within the above one-parameter potential family. This is not very surprising since general arguments from the microscopic cluster model show Friedrich (1981) that $`\mathrm{}`$-independent interaction potentials between nuclei should have low-energy non-physical bound states that simulate the Pauli principle between the constituting nucleons \[because of these so-called Pauli forbidden states (PFSs), such potentials are deep, in contrast with shallow potentials that only sustain physical bound states\]. For the $`\mathrm{}=2`$ partial wave of 12C + $`\alpha `$, the expected number of PFSs is two Langanke and Koonin (1985).

In the context of inversion, such bound states complicate things a lot since each of them introduces two additional arbitrary parameters (neither their energy nor their ANC is known experimentally). Following Ref. Sparenberg (2000), I introduce these states phenomenologically by regularizing the above inversion potential with the expression $`V_{\mathrm{reg}}+ar^2+br^3`$ for $`rr_{\mathrm{reg}}`$, where $`V_{\mathrm{reg}}`$ and $`r_{\mathrm{reg}}`$ are fitting parameters ($`a`$ and $`b`$ are calculated to make the potential and its first derivative continuous in $`r_{\mathrm{reg}}`$). Since this regularization slightly alters the phase shifts, parameters $`V`$ and $`R`$ in Eq. (3) have to be adapted to maintain the fit of Fig. 1. Altogether, there are thus five parameters ($`V`$, $`R`$, $`C_{12}`$, $`V_{\mathrm{reg}}`$ and $`r_{\mathrm{reg}}`$) which can be varied to fit all the $`\mathrm{}=0,2,4,6`$ data at the same time. The best fit I have obtained is for $`V=63.81`$ MeV fm2, $`R=4.925`$ fm, $`V_{\mathrm{reg}}=107.9241`$ MeV, $`r_{\mathrm{reg}}=4.1`$ fm. By construction, this inversion potential perfectly reproduces the $`2^+`$ state energy. The wave function is represented in Fig. 3 (solid line); its two nodes are due to its orthogonality to the PFSs at $`46`$ and $`18`$ MeV, which confirms that the potential is deep. The corresponding value for $`C_{12}`$ is $`137\times 10^3`$ fm-1/2.

I consider this prediction as reliable because all the even phase shifts in the elastic region are also rather well fitted (for example, the $`4^+`$ resonance is at $`\mathrm{3\hspace{0.17em}196}`$ keV with width 22 keV). However, this potential underestimates the $`0^+`$ bound-state energy ($`\mathrm{1\hspace{0.17em}696}`$ keV), as well as the $`6^+`$ resonance energy ($`\mathrm{8\hspace{0.17em}985}`$ keV) and width (259 keV). These inaccuracies are probably due to the extreme simplicity of the present model (energy- and $`\mathrm{}`$-independent local potential) and it is worth comparing it with more sophisticated models of the literature, which better take into account many-body effects. The two non-local $`\mathrm{}`$-dependent potentials of Ref. Langanke and Koonin (1985), which also fit the physical properties of the members of the rotational band, lead to similar wave functions for the 2+ subthreshold state (2 nodes). The Gaussian potential provides $`C_{12}=144\times 10^3`$ fm-1/2 while the Woods-Saxon potential provides $`153\times 10^3`$ fm-1/2. In Ref. Buck and Rubio (1984), a simple local potential is given, which precisely reproduces the energies of the members of the rotational band, as well as the resonance widths, provided a slight $`\mathrm{}`$ dependence is allowed. Although the number of PFSs is three for this potential, which adds a node in the subthreshold-bound-state wave function (see dash-dotted line in Fig. 3), it provides a value of $`C_{12}`$ very similar to the inversion potential: $`136\times 10^3`$ fm-1/2.

## V  Conclusion

I have illustrated on 12C + $`\alpha `$ two general results from the inverse problem in quantum scattering theory: (i) for a given partial wave, the scattering phase shifts, the number of bound states, their energies and their ANCs are independent of one another; (ii) when several partial waves can be described by the same theoretical model, these quantities may be related to one another. This is reasonably the case for 12C + $`\alpha `$, for which I have deduced the ANC value of the $`2^+`$ subthreshold state from the simultaneous analysis of the $`\mathrm{}=0,2,4,6`$ partial waves, with the help of a local energy- and angular-momentum-independent potential constructed by supersymmetric inversion.

Since this potential lies on rather restrictive hypotheses and is not able to exactly fit all the data, I have compared it with results from more sophisticated potential models available in the literature. All theoretical predictions lie in a rather limited range, $`C_{12}=144.5\pm 8.5\times 10^3`$ fm-1/2, which confirms somehow their overall validity. The error bar can thus be considered as a sensible estimate of the theoretical uncertainty on $`C_{12}`$ in a 12C + $`\alpha `$ potential model. The fact that both local (present work and Ref. Buck and Rubio (1984)) and non-local Langanke and Koonin (1985) potentials give consistent results is an indication that non-local effects are not very important for this particular state. This confirms that its structure is dominated by rigid 12C and $`\alpha `$ clusters with $`\mathrm{}=2`$ relative angular momentum. A similar smallness of non-local effects has been found in Ref. Baye and Descouvemont (1985) for the 16O + $`\alpha `$ system.

All the potential models discussed above neglect multichannel effects, which means that the predicted value of $`C_{12}`$ might be overestimated. This would be consistent with the experimental value obtained from transfer reactions Brune et al. (1999), $`C_{12}=114\pm 10\times 10^3`$ fm-1/2, and with the theoretical estimate from the microscopic multiconfiguration cluster model of Ref. Descouvemont (1987), $`C_{12}=134\times 10^3`$ fm-1/2. Whereas compatible with the unreliable value from Ref. Tischhauser et al. (2002), the above theoretical estimate excludes that from Ref. Angulo and Descouvemont (2000), as well as the cascade-transition value $`C_{12}=228_{37}^{+33}\times 10^3`$ fm-1/2 from Ref. Buchmann (2001). A new measurement of the cascade transitions is planned at TRIUMF, which hopefully should clarify the situation. In the meanwhile, R-matrix fits should reflect the uncertainty on $`C_{12}`$.

Let me end with a general discussion of possible theoretical models for 12C + $`\alpha `$ and similar systems relevant to nuclear astrophysics. Ideally, ab initio calculations should be performed; however, despite the success of the microscopic cluster model, such calculations are still too primitive to provide precise predictions of all available data. Hence, phenomenological models still have an important role to play. Among them, the potential model is probably the most physical and intuitive one. However, simulating the full complexity of the many-body nuclear system requires a non-local Langanke and Koonin (1985) multichannel Baldock et al. (1984) treatment. This considerably complicates the model and a good quality fit of all available data with the potential model is not available to date. In contrast, the R-matrix model allows such a fit because it remains rather simple even in the presence of many channels. The price to pay for that is a large number of parameters and a very phenomenological description, in particular of background terms. An interesting approach which combines the advantages of both the potential and R-matrix models is that of “hybrid” models Johnson (1973); Koonin et al. (1974); Langanke and Koonin (1985). The potentials constructed in the present work renew the interest for these models. In particular, purely repulsive potentials of the form (3) could be used as a new description of the background, a possibility which will be explored in a future work.

###### Acknowledgements.

I acknowledge very useful discussions with C. Barbieri, F. Barker, D. Baye, L. Buchmann, P. Descouvemont, A. Mukhamedzhanov and B. Jennings. I also thank L. Trache and C. Brune for drawing my attention to Refs. Brune et al. (1999) and Iwinski et al. (1984) respectively.


text:
# Island Resonances

## Introduction

Resonances driven by the error field multipoles often show themselves as islands in phase plots of the particle motion. This happens when there is another non-linear field present, besides the error field multipoles, which produces a strong enough dependence of the tune on the amplitude of the particle motion. Non-linear fields that produce a dependence of the tune on the ampliude can be the quadrupole fringe field , space charge fields and octupole correctors. In some cases it may be advantageous to introduce a non-linear field that produces a dependence of the tune on the amplitude of the particle motion. In order to be able to measure the strength of the island resonance, and then to correct it, it is helpful to study how these resonances show themselves. In particular , one can study how the presence of the island resonances distort the tune dependence on the amplitude and the emittance growth caused by the island resonance. These effects are to some extent measureable, and can lead to a way to correct the resonance. Islands in 4-dimensional phase space, unlike islands in 2-dimensional phase space, are not easy to visualize. This study will show how to visualize the islands in 4-dimensional phase space by studying the tune dependence on the amplitude of the particle. Again, this effect is to some extent measureable, and can lead to a way to measure the strength of the resonance.

## Resonances in 2-dimenional phase space

This section will study an island resonance in 2-dimensional phase space, the $`3\nu _x=19`$ resonance in the storage ring of the SNS (Spallation Neutron Source at Oak Ridge). This resonance is driven by the random error sextupole, $`b_2`$, in the magnets. In the absence of space charge forces, a non-linear field that produces an appreciable tune dependence on the amplitude of the particle motion is provided by the fringe field of the quadrupoles , which causes the tune to increase with the amplitude . Starting with a zero amplitude tune of $`\nu _x=6.3267`$ , $`3\nu _x=19.02`$, one finds the phase space plot shown in Fig. 1. This plot was generated starting with $`x_0=24mm`$, $`p_{x_0}=0`$ and then increasing $`x_0`$ in steps of 2 mm keeping $`p_{x_0}=0`$ and tracking the particle for 1000 turns for each $`x_0`$. The tune starts below the resonance, and increases with amplitude. The islands appear when the tune reaches $`\nu _x=6.3333`$. The islands are reached when the initial amplitude is increased to $`x_0=30mm`$, $`p_{x_0}=0`$, which corresponds to an initial horizontal emittance of about $`ϵ_{x0}=115`$ mm.mrad.

The distortion of the tune dependence on amplitude is shown in Fig. 2 where the horizontal tune is plotted against $`ϵ_{x0}`$. Fig. 2 is generated by doing a search along the $`p_{x_0}=0`$ direction, increasing $`x_0`$ in steps of 2 mm. Fig. 2 shows a flat region where $`\nu _x=6.3333`$ which will be seen to indicate the region where $`x_0,p_{x_0}`$ is crossing the island. When $`x_0,p_{x_0}`$ are inside the island, the paticle motion is doing a slow oscilation around the fixed point at the center of the island. The particle motion then contains more than one tune, and the tune with the largest ampltude is the tune of the fixed point, 6.3333. Inside the island, the tune shown in Fig. 2, is this tune with the largest amplitude. One can verify that the beginning and the end of the flat region occur at the same $`ϵ_{x0}`$ which correpond to the borders of the islands along the $`p_{x_0}=0`$ direction. Fig. 2 can be used to find the width of the island along the $`p_{x_0}=0`$ direction which is given by the width of the flat region in Fig. 2.

The islands generated by the $`3\nu _x=19`$ resonance indicate a growth in the particle emittance. If one starts a particle inside the islands with the initial emittance $`ϵ_{x0}`$, then as the particle moves around the fixed point, its emittance will change and reach the maximun value of $`ϵ_{x,max}`$. One can use as a measure of the emittance growth the quantity

$$dϵ_{x,max}=(ϵ_{x,max}ϵ_{x0})/ϵ_{x0}$$

. $`dϵ_{x,max}`$ is plotted against $`ϵ_{x0}`$ in Fig. 3. It shows a maximun value for $`dϵ_{x,max}`$ which can be used as a measure of the strength of the resonance.

The results shown in Fig. 2 and Fig. 3 can to some extent be measured experimentally to find a way to correct this resonance with sextupole correctors. An interesting plot that may corrrespond more closely to something that might be measured is to plot the emittance growth, $`dϵ_{x,max}`$ vs the tune, $`\nu _x`$. This plot can be found by combining Fig. 2 and Fig. 3 and is shown in Fig. 4. In this plot, there is a peak in the emittance growth which occurs when $`\nu _x`$=6.3333, the resonance tune. The peak could be used as a measure of the resonance strength to correct the resonance.

Using the flat region in the tune dependence on amplitude along the $`p_{x0}=0`$ direction to indicate the strength of the resonance can sometimes lead to an error in setting the strengths of the sextupole correctors. Changing the strength of the correctors can cause the islands to move, and the width of the island as indicated by the flat region in the tune dependence on the amplitude may, for example, appear smaller because one is now crossing the island at a place where it is narrower. Also, when the islands move the the emittance growth , $`dϵ_{x,max}`$ may seem smaller. In order to avoid this error due to the movement of the islands when the corrector strengths are changed, one has to measure the tune dependence on amplitude along enough different directions in phase space so that one finds the largest width of the island and the largest emittance growth for the $`x_0,p_{x0}`$ inside the island. The same sort of argument also applies in finding the emittance growth dependence on the amplitude. In Fig. 5 and Fig. 6, results are given found by going out along many directions in phase space. One needs just enough directions to cover one of the three islands in Fig. 1. Here 6 directions were used in the region of phase space covered by one island.

In Fig. 5 that plots $`\nu _x`$ versus $`ϵ_{x0}`$, for each $`ϵ_{x0}`$ there are now 6 points plotted corresponding to the 6 directions. Points that have the tune, nux=6.3333, lie inside the island.The width of the island in $`ϵ_{x0}`$ as given by the width of the points that lie on the nux=6.3333 line and is much larger than the width found using just the $`p_{x0}`$=0 direction as it includes the direction that crosses the island close to where the island is widest. In the same way, Fig. 6 now shows the largest emittance growth for the $`x_0,p_{x0}`$ that are inside the island.

In applying the above results to develope a procedure for using the correctors to correct a resonance, one will probably use a sample of the $`x_0,p_{x0}`$ that is convenient for a particular ring and its injection system. One would then measure either the tune or the emittance growth of enough particles in the sample to be able to find the width of the island from the tune measurement or the largest emittance growth for paticles in the sample. One has to have a large enough sample that one is not misled by the movement of the islands when the excitation of the correctors are changed. To measure the largest emittance growth for particles in the sample, one can measure the frequency spectrum of the betatron oscillations. Particles inside an island will have the frequency that corresponds to the resonant tune 6.3333. The largest amplitude found in the frequency spectrum for the resonant frequency, can be used as a measure of the largest emittance growth. Something like this was done at RHIC . The measurements can be done for a weak beam (no space charge effects) as it has been shown  that the resonance correction for resonances generated by magnet errors in the absence of space charge will also work fairly well in the presence of space charge.

An interesting plot that may corrrespond more closely to something that might be measured is to plot the emittance growth, $`dϵ_{x,max}`$ vs the tune, $`\nu _x`$. This plot can be found by combining Fig. 5 and Fig. 6 and is shown in Fig. 7

The correction of the island resonances for the $`3\nu _x=19`$ example can be done with two sextupole correctors properly located around the ring.The two correctors can then be adjusted one at a time, and for each setting of the correctors, the island width can measured from the dependnce of the tune on the amplitude as shown in Fig. 5, or the dependence of the amplitude growth on the tune as shown by combining Fig. 5 and Fig. 6. Simulating this procedure gives the results shown in Fig. 8 and Fig. 9. The emittance growth ,the maximun depx,max, is reduced by almost a factor of 10 by the best setting of the correctors that was found.

The theoretical result for the width of the island indicates that the width can be made zero when the correctors are set so that $`d\nu _{30}`$=0 where $`d\nu _{30}`$ is the stopband  of the $`3\nu _x=19`$ resonance due to the field errors in the magnets that are driving this resonance. According to the theoretical result the required setting of the correctors does not depend on the field producing the dependence of the tune on amplitude. This result is not useful for setting the correctors as the errors in the magnets are not known. However, this result helps to explain why the setting of the correctors in the absence of space charge fields will also correct the resonance when the space charge fields are present. It was found for the $`3\nu _x=19`$ resonance that setting of the correctors to make $`d\nu _{30}`$=0 will only reduce the emittance growth, the maximun $`dϵ_{x,max}`$, by a factor of about 3 instead of the factor of 10 found above. It was found that for the SNS  this weaker correction that makes $`d\nu _{30}`$=0 was good enough to reduce the beam losses due to the resonance and space charge fields.

## Resonances in 4-dimenional phase space

This section will study an island resonance in 4-dimensional phase space, the $`\nu _x+2\nu _y=19`$ in the storage ring of the SNS (Spallation Neutron Source at Oak Ridge). This resonance is driven by the random error sextupole, $`b_2`$, in the magnets. Unlike, the $`3\nu _x=19`$ resonance in 2-dimensional phase space, it is difficult to visualise the islands associated with the $`\nu _x+2\nu _y=19`$ resonance by looking at the particle motion in 4-dimensional phase space. Instead, we will look at the distortion due to the resonance on the dependence of the tune and the emittance growth on the amplitude of the particle motion. In Fig. 10, the spread in total emittance ,

$$dϵ_t=(ϵ_{t,max}ϵ_{t,min})/(ϵ_{t,max}+ϵ_{t,min})$$

is plotted against the initial total emittance, $`ϵ_{t0}`$. $`ϵ_{t,max}`$ is the largest total emittance and $`ϵ_{t,min}`$ is the smallest total emittance reached by particle with the initial coordinates, $`x_0,p_{x0},y_0,p_{y0}`$, and initial total emittance, $`ϵ_{t0}`$,tracked for 1000 turns. Points were found along a particular direction in phase space given by $`p_{x0}=0`$, $`p_{y0}=0`$ and $`ϵ_x=ϵ_y`$.The emittance spread becomes large in the region $`ϵ_{t0}`$=125 to $`ϵ_{t0}`$=300, and this region may be used as a measure of the width of the island along this particular direction , $`p_{x0}=0`$, $`p_{y0}=0`$ and $`ϵ_x=ϵ_y`$. The low value of $`dϵ_t`$ near the center of region shows that this particular direction goes close to the center of the island where we expect $`dϵ_t`$=0. Fig. 11 shows the emittance growth,

$$dϵ_{t,max}=(ϵ_{t,max}ϵ_{t0})/ϵ_{t0}$$

as a function of $`ϵ_{t0}`$. The goal of a correction scheme would be to reduce the maximun value of $`dϵ_{t,max}`$ in ths plot. In the 2-dimensional case of the $`3\nu _x=19`$ resonance, the plot of $`\nu _x`$ vs. $`ϵ_{x0}`$ had a flat region where $`\nu _x`$ had the constant value of 6.3333 which gave the width of the island. It was found for the $`\nu _x+2\nu _y=19`$ resonance, a flat region would be found if we plotted $`\nu _x+2\nu _y`$ vs. $`ϵ_{t0}`$ and the constant value in the flat region is $`\nu _x+2\nu _y=19`$. This is shown in Fig. 12. The width of the flat region coincides with the width of the island as found from Fig. 10, which plots the emittance spread vs. $`ϵ_{t0}`$. For coupled motion, the particle motion contains more than one tune, and what is plotted in Fig. 12 is the main tune or the tune with the largest amplitude.

An interesting plot that may corrrespond more closely to something that might be measured is to plot the emittance growth, $`dϵ_{t,max}`$ vs $`\nu _x+2\nu _y`$. This plot can be found by combining Fig. 11 and Fig. 12 and is shown in Fig. 13. In this plot, there is a peak in the emittance growth which occurs when $`\nu _x+2\nu _y=19`$. The peak could be used as a measure of the resonance strength to correct the resonance.

As in the 2-dimesional case, to correct the $`\nu _x+2\nu _y=19`$ resonance , one has to repeat the above calculations for all directions in phase space that pass through a particular island, in order not to be misled by the movement of the islands when the correctors are changed. These results are shown in Fig. 14, Fig. 15.and Fig. 16 where for each $`ϵ_{t0}`$ runs are done for 6 directions in $`x_0,p_{x0}`$ and for 6 directions in $`y_0,p_{y0}`$. In these figures there are, then, 36 points for each value of $`ϵ_{t0}`$. The results indicated by these figures for the width of the island and the maximun emittance growth are not too different from those found above where only one direction in phase space was used, as the particular direction used appears to pass close to the center of the island.

An interesting plot that may corrrespond more closely to something that might be measured is to plot the emittance growth , $`dϵ_{t,max}`$ vs $`\nu _x+2\nu _y`$. This plot can be found by combining Fig. 15 and Fig. 16 and is shown in Fig. 17. In this plot, there is a peak in the emittance growth which occurs when $`\nu _x+2\nu _y=19`$. The peak could be used as a measure of the resonance strength to correct the resonance.

The correction of the island resonances for the $`\nu _x+2\nu _y=19`$ resonance can be done with two sextupole correctors properly located around the ring. Assuming that one has a way to measure either the width of the island or the maximun emittance growth using the above results as a guide , then the two correctors can be adjusted one at a time to reduce the emittance growth or the island width. Simulating this procedure gives the results shown in Fig. 18 and Fig. 19. The emittance growth ,the maximun $`dϵ_{t,max}`$, is reduced by almost a factor of 15 by the best setting of the correctors that was found. In this case, setting the correctors to zero the stopband, $`d\nu _{12}`$=0, reduces the maximun emittance growth by a factor of 5.

## Higher order resonances

In the above, results were presented for the 3nux=19 and the $`\nu _x+2\nu _y=19`$ resonances. Similar studies have also been done for all four of the third order resonances $`m\nu _x+n\nu _y=19`$, $`m+n=3`$, and for all five of the fourth order resonances $`m\nu _x+n\nu _y=25`$, $`m+n=4`$, with similar results. It is our expectation that similar results would be found for the higher order resonance $`m\nu _x+n\nu _y=q`$, where m,n,q are integers. In particular, if $`m\nu _x+n\nu _y`$ is plotted versus $`ϵ_{t0}`$, this plot will contain a flat region where $`m\nu _x+n\nu _y`$ is costant at $`m\nu _x+n\nu _y=q`$, and the width of this flat region will give the width of the islands associated with the $`m\nu _x+n\nu _y=q`$ resonance.

## ACKNOWLEDGEMENTS

We would like to thank A.V. Fedotov for many discussions.