text:
# Geometric Numerical Integration Applied to The Elastic Pendulum at Higher Order Resonance

## 1 Introduction

Higher order resonances are known to have a long time-scale behaviour. From an asymptotic point of view, a first order approximation (such as first order averaging) would not be able to clarify the interesting dynamics in such a system. Numerically, this means that the integration times needed to capture such behaviour are significantly increased. In this paper we present a reasonably cheap method to achieve a qualitatively good result even after long integration times.

Geometric numerical integration methods for (ordinary) differential equations () have emerged in the last decade as alternatives to traditional methods (e.g. Runge-Kutta methods). Geometric methods are designed to preserve certain properties of a given ODE exactly (i.e. without truncation error). The use of geometric methods is particularly important for long integration times. Examples of geometric integration methods include symplectic integrators, volume-preserving integrators, symmetry-preserving integrators, integrators that preserve first integrals (e.g. energy), Lie-group integrators, etc. A survey is given in .

It is well known that resonances play an important role in determining the dynamics of a given system. In practice, higher order resonances occur more often than lower order ones, but their analysis is more complicated. In , Sanders was the first to give an upper bound on the size of the resonance domain (the region where interesting dynamics takes place) in two degrees of freedom Hamiltonian systems. Numerical studies by van den Broek , however, provided evidence that the resonance domain is actually much smaller. In , Tuwankotta and Verhulst derived improved estimates for the size of the resonance domain, and provided numerical evidence that for the $`4:1`$ and the $`6:1`$ resonances of the elastic pendulum, their estimates are sharp. The numerical method they used in their analysis 333A Runge-Kutta method of order $`78`$, however, was not powerful enough to be applied to higher order resonances. In this paper we construct a symplectic integration method, and use it to show numerically that the estimates of the size of the resonance domain in  are also sharp for the $`4:3`$ and the $`3:1`$ resonances.

Another subtle problem regarding to this resonance manifold is the bifurcation of this manifold as the energy increases. To study this problem numerically one would need a numerical method which is reasonably cheap and accurate after a long integration times.

In this paper we will use the elastic pendulum as an example. The elastic pendulum is a well known (classical) mechanical problem which has been studied by many authors. One of the reasons is that the elastic pendulum can serve as a model for many problems in different fields. See the references in . In itself, the elastic pendulum is a very rich dynamical system. For different resonances, it can serve as an example of a chaotic system, an auto-parametric excitation system (), or even a linearizable system. The system also has (discrete) symmetries which turn out to cause degeneracy in the normal form.

We will first give a brief introduction to the splitting method which is the main ingredient for the symplectic integrator in this paper. We will then collect the analytical results on the elastic pendulum that have been found by various authors. Mostly, in this paper we will be concerned with the higher order resonances in the system. All of this will be done in the next two sections of the paper. In the fourth section we will compare our symplectic integrator with the standard 4-th order Runge-Kutta method and also with an order $`78`$ Runge-Kutta method. We end the fourth section by calculating the size of the resonance domain of the elastic pendulum at higher order resonance.

## 2 Symplectic Integration

Consider a symplectic space $`\mathrm{\Omega }=^{2n},n`$ where each element $`𝝃`$ in $`\mathrm{\Omega }`$ has coordinate $`(𝒒,𝒑)`$ and the symplectic form is $`\text{d}𝒒\text{d}𝒑`$. For any two functions $`F,G𝒞^{\mathrm{}}(\mathrm{\Omega })`$ define

$$\{F,G\}=\underset{1}{\overset{n}{}}\left(\frac{F}{q_j}\frac{G}{p_j}\frac{G}{q_j}\frac{F}{p_j}\right)𝒞^{\mathrm{}}(\mathrm{\Omega }),$$

which is called the Poisson bracket of $`F`$ and $`G`$. Every function $`H𝒞^{\mathrm{}}(\mathrm{\Omega })`$ generates a (Hamiltonian) vector field defined by $`\{q_i,H\},\{p_i,H\},i=1,\mathrm{},n`$. The dynamics of $`H`$ is then governed by the equations of motion of the form

$$\begin{array}{cc}\hfill \dot{q_i}=& \{q_i,H\}\hfill \\ \hfill \dot{p_i}=& \{p_i,H\},i=1,\mathrm{},n.\hfill \end{array}$$

Let $`X`$ and $`Y`$ be two Hamiltonian vector fields, defined in $`\mathrm{\Omega }`$, associated with Hamiltonians $`H_X`$ and $`H_Y`$ in $`𝒞^{\mathrm{}}(\mathrm{\Omega })`$ respectively. Consider another vector field $`[X,Y]`$ which is just the commutator of the vector fields $`X`$ and $`Y`$ . Then $`[X,Y]`$ is also a Hamiltonian vector field with Hamiltonian $`H_{[X,Y]}=\{H_X,H_Y\}`$. See for example  for details.

We can write the flow of the Hamiltonian vector fields $`X`$ as

$$\phi _{X;t}=\text{ exp}(tX)I+tX+\frac{1}{2!}(tX)^2+\frac{1}{3!}(tX)^3+\mathrm{}$$

(and so does the flow of $`Y`$). By the BCH formula, there exists a (formal) Hamiltonian vector field $`Z`$ such that

$$Z=(X+Y)+\frac{t}{2}[X,Y]+\frac{t^2}{12}\left([X,X,Y]+[Y,Y,X]\right)+O(t^3)$$
(1)
and $`\text{exp}(tZ)=\text{exp}(tX)\text{exp}(tY)`$, where $`[X,X,Y]=[X,[X,Y]]`$ , and so on. Moreover, Yoshida (in ) shows that $`\text{exp}(tX)\text{exp}(tY)\text{exp}(tX)=\text{exp}(tZ),`$ where

$$Z=(2X+Y)+\frac{t^2}{6}\left([Y,Y,X][X,X,Y]\right)+O(t^4).$$
(2)
We note that in terms of the flow, the multiplication of the exponentials above means composition of the corresponding flow, i.e. $`\phi _{Y;t}\phi _{X;t}`$.

Let $`\tau `$ be a small positive number and consider a Hamiltonian system with Hamiltonian $`H(𝝃)=H_X(𝝃)+H_Y(𝝃)`$, where $`𝝃\mathrm{\Omega }`$, and $`\dot{𝝃}=X+Y`$. Using (1) we have that $`\phi _{Y;\tau }\phi _{X;\tau }`$ is (approximately) the flow of a Hamiltonian system

$$\dot{𝝃}=(X+Y)+\frac{\tau }{2}[X,Y]+\frac{\tau ^2}{12}\left([X,X,Y]+[Y,Y,X]\right)+O(\tau ^3),$$

with Hamiltonian

$$H_\tau =H_X+H_Y+\frac{\tau }{2}\{H_X,H_Y\}+\frac{\tau ^2}{12}\left(\{H_X,H_X,H_Y\}+\{H_Y,H_Y,H_X\}\right)+O(\tau ^3).$$

Note that $`\{H,K,F\}=\{H,\{K,F\}\}`$. This mean that $`HH_\tau =O(\tau )`$ or, in other words

$$\phi _{Y;\tau }\phi _{X;\tau }=\phi _{X+Y}(\tau )+O(\tau ^2).$$
(3)
As before and using (2), we conclude that

$$\phi _{X;\frac{\tau }{2}}\phi _{Y;\tau }\phi _{X;\frac{\tau }{2}}=\phi _{X+Y}(\tau )+O(\tau ^3).$$
(4)
Suppose that $`\psi _{X;\tau }`$ and $`\psi _{Y;\tau }`$ are numerical integrators of system $`\dot{𝝃}=X`$ and $`\dot{𝝃}=Y`$ (respectively). We can use symmetric composition (see ) to improve the accuracy of $`\psi _{X+Y;\tau }`$. If $`\psi _{Y;\tau }`$ and $`\psi _{X;\tau }`$ are symplectic, then the composition forms a symplectic numerical integrator for $`X+Y`$. See  for more discussion; also  for references. If we can split $`H`$ into two (or more) parts which Poisson commute with each other (i.e. the Poisson brackets between each pair vanish), then we have $`H=H_\tau `$. This implies that in this case the accuracy of the approximation depends only on the accuracy of the integrators for $`X`$ and $`Y`$. An example of this case is when we are integrating the Birkhoff normal form of a Hamiltonian system.

## 3 The Elastic Pendulum

Consider a spring with spring constant $`s`$ and length $`l_{}`$ to which a mass $`m`$ is attached. Let $`g`$ be the gravitational constant and $`l`$ the length of the spring under load in the vertical position, and let $`r`$ be the distance between the mass $`m`$ and the suspension point. The spring can both oscillate in the radial direction and swing like a pendulum. This is called the elastic pendulum. See Figure (1) for illustration and  (or ) for references.

The phase space is $`^4`$ with canonical coordinate $`𝝃=(z,\phi ,p_z,p_\phi )`$, where $`z=(rl_{})/l_{}`$. Writing the linear frequencies of the Hamiltonian as $`\omega _z=\sqrt{s/m}`$ and $`\omega _\phi =\sqrt{g/l}`$, the Hamiltonian of the elastic pendulum becomes

$$H=\frac{1}{2\sigma }\left(p_z^2+\frac{p_\phi ^2}{(z+1)^2}\right)+\frac{\sigma }{2}\omega _{z}^{}{}_{}{}^{2}\left(z+\left(\frac{\omega _\phi }{\omega _z}\right)^2\right)^2\sigma \omega _{\phi }^{}{}_{}{}^{2}(z+1)\mathrm{cos}\phi ,$$
(5)
where $`\sigma =ml^2`$. By choosing the right physical dimensions, we can scale out $`\sigma `$. We remark that for the elastic pendulum as illustrated in Figure 1, we have $`\omega _z\omega _\phi `$. See  for details. It is clear that this system possesses symmetry

$$T:(z,\phi ,p_z,p_\phi ,t)(z,\phi ,p_z,p_\phi ,t)$$
(6)
and the reversing symmetries

$$\begin{array}{c}R_1:(z,\phi ,p_z,p_\phi ,t)(z,\phi ,p_z,p_\phi ,t),\\ R_2:(z,\phi ,p_z,p_\phi ,t)(z,\phi ,p_z,p_\phi ,t).\end{array}$$
(7)
If there exist two integers $`k_1`$ and $`k_2`$ such that $`k_1\omega _z+k_2\omega _\phi =0`$, then we say $`\omega _z`$ and $`\omega _\phi `$ are in resonance. Assuming $`(|k_1|,|k_2|)=1`$, we can divide the resonances in two types, e.g. lower order resonance if $`|k_1|+|k_2|<5`$ and higher order resonance if $`|k_1|+|k_2|5`$. In the theory of normal forms, the type of normal form of the Hamiltonian is highly dependent on the type of resonance in the system. See .

In general, the elastic pendulum has at least one fixed point which is the origin of phase space. This fixed point is elliptic. For some of the resonances, there is also another fixed point which is of the saddle type, i.e. $`(z,\phi ,p_z,p_\phi )=(2(\omega _\phi /\omega _z)^2,\pi ,0,0)`$. From the definition of $`z`$, it is clear that the latter fixed point only exists for $`\omega _z/\omega _\phi >\sqrt{2}`$. The elastic pendulum also has a special periodic solution in which $`\phi =p_\phi =0`$ (the normal mode). This normal mode is an exact solution of the system derived from (5). We note that there is no nontrivial solution of the form $`(0,\phi (t),0,p_\phi (t))`$.

Now we turn our attention to the neighborhood of the origin. We refer to  for the complete derivation of the following Taylor expansion of the Hamiltonian (we have dropped the bar)

$$H=H_2+\epsilon H_3+\epsilon ^2H_4+\epsilon ^3H_5+\mathrm{},$$
(8)
with

$$\begin{array}{cc}\hfill H_2=& \frac{1}{2}\omega _z\left(z^2+p_z^2\right)+\frac{1}{2}\omega _\phi \left(\phi ^2+p_\phi ^2\right)\hfill \\ \hfill H_3=& \frac{\omega _\phi }{\sqrt{\omega _z}}\left(\frac{1}{2}z\phi ^2zp_\phi ^2\right)\hfill \\ \hfill H_4=& \left(\frac{3}{2}\frac{\omega _\phi }{\omega _z}z^2p_\phi ^2\frac{1}{24}\phi ^4\right)\hfill \\ \hfill H_5=& \frac{1}{\sqrt{\omega _z}}\left(\frac{1}{24}z\phi ^4+2\frac{\omega _\phi }{\omega _z}z^3p_\phi ^2\right)\hfill \\ & \mathrm{}\hfill \end{array}$$

In  the $`2:1`$-resonance of the elastic pendulum has been studied intensively. At this specific resonance, the system exhibits an interesting phenomenon called auto-parametric excitation, e.g. if we start at any initial condition arbitrarily close to the normal mode, then we will see energy interchanging between the oscillating and swinging motion. In , the author shows that the normal mode solution (which is the vertical oscillation) is unstable and therefore, gives an explanation of the auto-parametric behavior.

Next we consider two limiting cases of the resonances, i.e. when $`\omega _z/\omega _\phi \mathrm{}`$ and $`\omega _z/\omega _\phi 1`$. The first limiting case can be interpreted as a case with a very large spring constant so that the vertical oscillation can be neglected. The spring pendulum then becomes an ordinary pendulum; thus the system is integrable. The other limiting case is interpreted as the case where $`l_{}=0`$ (or very weak spring) 444This case is unrealistic for the model illustrated in Figure 1. A more realistic mechanical model with the same Hamiltonian (5) can be constructed by only allowing some part of the spring to swing. Using the transformation $`r=l(z+1)`$, $`x=r\mathrm{cos}\phi `$ and $`y=r\mathrm{sin}\phi ,`$ we transform the Hamiltonian (5) to the Hamiltonian of the harmonic oscillator. Thus this case is also integrable. Furthermore, in this case all solutions are periodic with the same period which is known as isochronism. This means that we can remove the dependence of the period of oscillation of the mathematical pendulum on the amplitude, using this specific spring. We note that this isochronism is not derived from the normal form (as in ) but exact.

All other resonances are higher order resonances. In two degrees of freedom (which is the case we consider), for fixed small energy the phase space of the system near the origin looks like the phase space of decoupled harmonic oscillator. A consequence of this fact is that in the neighborhood of the origin, there is no interaction between the two degrees of freedom. The normal mode (if it exists), then becomes elliptic (thus stable).

Another possible feature of this type of resonance is the existence of a resonance manifold containing periodic solutions (see  paragraph 4.8). We remark that the existence of this resonance manifold does not depend on whether the system is integrable or not. In the resonance domain (i.e. the neighborhood of the resonant manifold), interesting dynamics (in the sense of energy interchanging between the two degrees of freedom) takes place (see ). Both the size of the domain where the dynamics takes place and the time-scale of interaction are characterized by $`\epsilon `$ and the order of the resonance, i.e. the estimate of the size of the domain is

$$d_\epsilon =O\left(\epsilon ^{\frac{m+n4}{2}}\right)$$
(9)
and the time-scale of interactions is $`O(\epsilon ^{\frac{m+n}{2}})`$ for $`\omega _z:\omega _\phi =m:n`$ with $`(m,n)=1`$. 555Due to a particular symmetry, some of the lower order resonances become higher order resonances (). In those cases, $`(m,n)=1`$ need not hold. We note that for some of the higher order resonances where $`\omega _z/\omega _\phi 1`$ the resonance manifold fails to exist. See  for details.

## 4 Numerical Studies on the Elastic Pendulum

One of the aims of this study is to construct a numerical Poincaré map ($`𝒫`$) for the elastic pendulum in higher order resonance. As is explained in the previous section, interesting dynamics of the higher order resonances takes place in a rather small part of phase space. Moreover, the interaction time-scale is also rather long. For these two reasons, we need a numerical method which preserves qualitative behavior after a long time of integration. Obviously by decreasing the time step of any standard integrator (e.g. Runge-Kutta method), we would get a better result. As a consequence however, the actual computation time would become prohibitively long. Under these constraints, we would like to propose by means of an example that symplectic integrators offer reliable and reasonably cheap methods to obtain qualitatively good phase portraits.

We have selected four of the most prominent higher order resonances in the elastic pendulum. For each of the chosen resonances, we derive its corresponding equations of motion from (8). This is done because the dependence on the small parameter $`\epsilon `$ is more visible there than in (5). Also from the asymptotic analysis point of view, we know that (8) truncated to a sufficient degree has enough ingredients of the dynamics of (5).

The map $`𝒫`$ is constructed as follows. We choose the initial values $`𝝃_{\mathbf{}}`$ in such a way that they all lie in the approximate energy manifold $`H_2=E_{}`$ and in the section $`\mathrm{\Sigma }=\{𝝃=(z,\phi ,p_z,p_\phi )|z=0,p_z>0\}`$. We follow the numerically constructed trajectory corresponding to $`𝝃_{\mathbf{}}`$ and take the intersection of the trajectory with section $`\mathrm{\Sigma }`$. The intersection point is defined as $`𝒫(𝝃_{\mathbf{}})`$. Starting from $`𝒫(𝝃_{\mathbf{}})`$ as an initial value, we go on integrating and in the same way we find $`𝒫^2(𝝃_{\mathbf{}})`$, and so on.

The best way of measuring the performance of a numerical integrator is by comparing with an exact solution. Due to the presence of the normal mode solution (as an exact solution), we can check the performance of the numerical integrator in this way (we will do this in section 4.2) . Nevertheless, we should remark that none of the nonlinear terms play a part in this normal mode solution. Recall that the normal mode is found in the invariant manifold $`\{(z,\phi ,p_z,p_\phi |\phi =p_\phi =0\}`$ and in this manifold the equations of motion of (8) are linear.

Another way of measuring the performance of an integrator is to compare it with other methods. One of the best known methods for time integration are the Runge-Kutta methods (see ). We will compare our integrator with a higher order ($`7`$-$`8`$ order) Runge-Kutta method (RK78). The RK78 is based on the method of Runge-Kutta-Felbergh (). The advantage of this method is that it provides step-size control. As is indicated by the name of the method, to choose the optimal step size it compares the discretizations using $`7`$-th order and $`8`$-th order Runge-Kutta methods. A nice discussion on lower order methods of this type, can be found in  pp. 448-454. The coefficients in this method are not uniquely determined. For RK78 that we used in this paper, the coefficients were calculated by C. Simo from the University of Barcelona. We will also compare the symplectic integrator (SI) to the standard $`4`$-th order Runge-Kutta method.

We will first describe the splitting of the Hamiltonian which is at the core of the symplectic integration method in this paper. By combining the flow of each part of the Hamiltonian, we construct a $`4`$-th order symplectic integrator. The symplecticity is obvious since it is the composition of exact Hamiltonian flows. Next we will show the numerical comparison between the three integrators, RK78, SI and RK4. We compare them to an exact solution. We will also show the performance of the numerical integrators with respect to energy preservation. We note that SI are not designed to preserve energy (see ). Because RK78 is a higher order method (thus more accurate), we will also compare the orbit of RK4 and SI. We will end this section with results on the size of the resonance domain calculated by the SI method.

### 4.1 The Splitting of the Hamiltonian

Consider again the expanded Hamiltonian of the elastic pendulum (8). We split this Hamiltonian into integrable parts: $`H=H^1+H^2+H^3`$, where

$$\begin{array}{cc}\hfill H^1=& \epsilon \frac{\omega _\phi }{2\sqrt{\omega _z}}z\phi ^2\epsilon ^2\frac{1}{24}\phi ^4\epsilon ^3\frac{1}{24\sqrt{\omega _z}}z\phi ^4+\mathrm{}\hfill \\ \hfill H^2=& \epsilon \frac{\omega _\phi }{\sqrt{\omega _z}}zp_\phi ^2+\epsilon ^2\frac{3}{2}\frac{\omega _\phi }{\omega _z}z^2p_\phi ^2\epsilon ^3\frac{2\omega _\phi }{\omega _z\sqrt{\omega _z}}z^3p_\phi ^2+\mathrm{}\hfill \\ \hfill H^3=& \frac{1}{2}\omega _z\left(z^2+p_z^2\right)+\frac{1}{2}\omega _\phi \left(\phi ^2+p_\phi ^2\right).\hfill \end{array}$$
(10)
Note that the equations of motion derived from each part of the Hamiltonian can be integrated exactly; thus we know the exact flow $`\phi _{1;\tau }`$, $`\phi _{2:\tau }`$, and $`\phi _{3;\tau }`$ corresponding to $`H^1`$,$`H^2`$, and $`H^3`$ respectively. This splitting has the following advantages.

* It preserves the Hamiltonian structure of the system.
* It preserves the symmetry (6) and reversing symmetries (7) of $`H`$.
* $`H^1`$ and $`H^2`$ are of $`O(\epsilon )`$ compared with $`H`$ (or $`H^3`$).

Note that, for each resonance we will truncate (10) up to and including the degree where the resonant terms of the lowest order occur.

We define

$$\phi _\tau =\phi _{1;\tau /2}\phi _{2;\tau /2}\phi _{3;\tau }\phi _{2;\tau /2}\phi _{1;\tau /2}.$$
(11)
From section 2 we know that this is a second order method. Next we define $`\gamma =1/(2\sqrt[3]{2})`$ and $`\psi _\tau =\phi _{\gamma \tau }\phi _{(12\gamma )\tau }\phi _{\gamma \tau }`$ to get a fourth order method. This is known as the generalized Yoshida method (see ). By, Symplectic Integrator (SI) we will mean this fourth order method. This composition preserves the symplectic structure of the system, as well as the symmetry (6) and the reversing symmetries (7). This is in contrast with the Runge-Kutta methods which only preserves the symmetry (6), but not the symplectic structure, nor the reversing symmetries (7). As a consequence the Runge-Kutta methods do not preserve the KAM tori caused by symplecticity or reversibility.

### 4.2 Numerical Comparison between RK4, RK78 and SI

We start by comparing the three numerical methods, i.e. RK4, RK78, and SI. We choose the $`4:1`$-resonance, which is the most prominent higher order resonance, as a test problem. We fix the value of the energy ($`H_2`$) to be $`5`$ and take $`\epsilon =0.05`$. Starting at the initial condition $`z(0)=0,\phi (0)=0,p_z=\sqrt{5/2},`$ and $`p_\phi (0)=0`$, we know that the exact solution we are approximating is given by $`(\sqrt{5/2}\mathrm{sin}(4t),0,\sqrt{5/2}\mathrm{cos}(4t),0)`$. We integrate the equations of motion up to $`t=10^5`$ seconds and keep the result of the last $`10`$ seconds to have time series $`\overline{z}(t_n)`$ and $`\overline{p_z}(t_n)`$ produced by each integrator. Then we define a sequence $`s_n=99990+5n/100,n=0,1,\mathrm{},200`$. Using an interpolation method, for each of the time series we calculate the numerical $`\overline{z}(s_n)`$. In figure 2 we plot the error function $`\overline{z}(s_n)z(s_n)`$ for each integrator.

The plots in Figure in 2 clearly indicate the superiority of RK78 compared with the other methods (due to the higher order method). The error generated by RK78 is of order $`10^7`$ for an integration time of $`10^5`$ seconds. The minimum time step taken by RK78 is $`0.0228`$ and the maximum is $`0.0238`$. The error generated by SI on the other hand, is of order $`10^5`$. The CPU time of RK78 during this integration is $`667.75`$ seconds. SI completes the computation after $`446.72`$ seconds while RK4 only needs $`149.83`$ seconds.

We will now measure how well these integrators preserve energy. We start integrating from an initial condition $`z(0)=0`$, $`\phi (0)=1.55`$, $`p_\phi (0)=0`$ and $`p_z(0)`$ is determined from $`H^3=5`$ (in other word we integrate on the energy manifold $`H=5+O(\epsilon )`$). The small parameter is $`\epsilon =0.05`$ and we integrate for $`t=10^5`$ seconds.

For RK78, the integration takes $`667.42`$ second of CPU time. For RK4 and SI we used the same time step, that is $`10^2`$. RK4 takes $`377.35`$ seconds while SI takes $`807.01`$ seconds of CPU time. It is clear that SI, for this size of time step, is inefficient with regard to CPU time. This is due to the fact that to construct a higher order method we have to compose the flow several times. We plot the results of the last $`10`$ seconds of the integrations in Figure 3. We note that in these $`10`$ seconds, the largest time step used by RK78 is 0.02421…while the smallest is 0.02310…. It is clear from this, that even though the CPU time of RK4 is very good, the result in the sense of conservation of energy is rather poor relative to the other methods.

We increase the time step to $`0.05`$ and integrate the equations of motion starting at the same initial condition and for the same time. The CPU time of SI is now $`149.74`$ while for the RK4 it is $`76.07`$. Again, in Figure 3 (the right hand plots) we plot the energy against time. A significant difference between RK4 and SI then appears in the energy plots. The results of symplectic integration are still good compared with the higher order method RK78. On the other hand, the results from RK4 are far below the other two.

### 4.3 Computation of the Size of the Resonance Domain

Finally, we calculate the resonance domain for some of the most prominent higher order resonances for the elastic pendulum. In Figure 4 we give an example of the resonance domain for the $`4:1`$ resonance. We note that RK4 fails to produce the section.

On the other hand, the results from SI are still accurate. We compare the results from SI and RK78 in Figure 4. After $`5\times 10^3`$ seconds, one loop in the plot is completed. For that time of integration, RK78 takes $`34.92`$ seconds of CPU time, while SI takes only $`16.35`$ seconds. This is very useful since to calculate for smaller values of $`\epsilon `$ and higher resonance cases, the integration time is a lot longer which makes it impractical to use RK78.

In Table 1 we list the four most prominent higher order resonances for the elastic pendulum. This table is adopted from  where the authors list six of them.

The numerical size of the domain in table 1 is computed as follows. We first draw several orbits of the Poincaré maps $`𝒫`$. Using a twist map argument, we can locate the resonance domain. By adjusting the initial condition manually, we then approximate the heteroclinic cycle of $`𝒫`$. See figure 4 for illustration. Using interpolation we construct the function $`r_o(\theta )`$ which represent the distance of a point in the outer cycle to the origin and $`\theta `$ is the angle with respect to the positive horizontal axis. We do the same for the inner cycle and then calculate $`\mathrm{max}_\theta |r_o(\theta )r_i(\theta )|`$. The higher the resonance is, the more difficult to compute the size of the domain in this way.

For resonances with very high order, manually approximating the heteroclinic cycles would become impractical, and one could do the following. First we have to calculate the location of the fixed points of the iterated Poincaré maps numerically. Then we can construct approximations of the stable and unstable manifolds of one of the saddle points. By shooting to the next saddle point, we can make corrections to the approximate stable and unstable manifold of the fixed point.

## 5 Discussion

In this section we summarize the previous sections. First the performance of the integrators is summarized in table 2.

As indicated in table 2, for the $`4:3`$ and the $`3:1`$ resonances, the higher order Runge-Kutta method fails to produce the section. This is caused by the dissipation term, artificially introduced by this numerical method, which after a long time of integration starts to be more significant. On the other hand, we conclude that the results of our symplectic integrator are reliable. This conclusion is also supported by the numerical calculations of the size of the resonance domain (listed in Table 2).

In order to force the higher order Runge-Kutta method to be able to produce the section, one could also do the following. Keeping in mind that RK78 has automatic step size control based on the smoothness of the vector field, one could manually set the maximum time step for RK78 to be smaller than $`0.02310`$. This would make the integration times extremely long however.

We should remark that in this paper we have made a number of simplifications. One is that we have not used the original Hamiltonian. The truncated Taylor expansion of (5) is polynomial. Somehow this may have a smoothing effect on the Hamiltonian system. It would be interesting to see the effect of this simplification on the dynamics of the full system. Another simplification is that, instead of choosing our initial conditions in the energy manifold $`H=C`$, we are choosing them in $`H^3=C`$. By using the full Hamiltonian instead of the truncated Taylor expansion of the Hamiltonian, it would become easy to choose the initial conditions in the original energy manifold. Nevertheless, since in this paper we always start in the section $`\mathrm{\Sigma }`$, we know that we are actually approximating the original energy manifold up to order $`\epsilon ^2`$.

We also have not used the presence of the small parameter $`\epsilon `$ in the system. As noted in , it may be possible to improve our symplectic integrator using this small parameter. Still related to this small parameter, one also might ask whether it would be possible to go to even smaller values of $`\epsilon `$. In this paper we took $`e^3<\epsilon <e^{0.5}`$. As noted in the previous section, the method that we apply in this paper can not be used for computing the size of the resonance domain for very high order resonances. This is due to the fact that the resonance domain then becomes exceedingly small. This is more or less the same difficulty we might encounter if we decrease the value of $`\epsilon `$.

Another interesting posibility is to numerically follow the resonance manifold as the energy increases. As noted in the introduction, this is numerically difficult problem. Since this symplectic integration method offers a cheap and accurate way of producing the resonance domain, it might be posible to numerically study the bifurcation of the resonance manifold as the energy increases. Again, we note that to do so we would have to use the full Hamiltonian.

Acknowledgements

J.M. Tuwankotta thanks the School of Mathematical and Statistical Sciences, La Trobe University, Australia for their hospitality when he was visiting the university. Thanks to David McLaren of La Trobe University, and Ferdinand Verhulst, Menno Verbeek and Michiel Hochstenbach of Universiteit Utrecht, the Netherlands for their support and help during the execution of this research. Many thanks also to Santi Goenarso for every support she has given.

We are grateful to the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO) and to the Australian Research Council (ARC) for financial support.


format: text

text:
# References

DFTT 32/2000

INFNCA-TH0012

hep-ph/0008186

$`𝚲`$ polarization from unpolarized quark fragmentation
M. Anselmino1, D. Boer2, U. D’Alesio3 and F. Murgia3
1 Dipartimento di Fisica Teorica, Università di Torino and

INFN, Sezione di Torino, Via P. Giuria 1, I-10125 Torino, Italy

2 RIKEN-BNL Research Center

Brookhaven National Laboratory, Upton, NY 11973, USA

3 Istituto Nazionale di Fisica Nucleare, Sezione di Cagliari

and Dipartimento di Fisica, Università di Cagliari

C.P. 170, I-09042 Monserrato (CA), Italy

Abstract:
The longstanding problem of explaining the observed polarization of $`\mathrm{\Lambda }`$ hyperons inclusively produced in the high energy collisions of unpolarized hadrons is tackled by considering spin and $`𝒌_{}`$ dependent quark fragmentation functions. The data on $`\mathrm{\Lambda }`$’s and $`\overline{\mathrm{\Lambda }}`$’s produced in $`pN`$ processes are used to determine simple phenomenological expressions for these new “polarizing fragmentation functions”, which describe the experiments remarkably well.

1. Introduction It is well known since a long time that $`\mathrm{\Lambda }`$ hyperons produced with $`x_F\stackrel{>}{}0.2`$ and $`p_T\stackrel{>}{}`$ 1 GeV/$`c`$ in the collision of two unpolarized hadrons, $`AB\mathrm{\Lambda }^{}X`$, are polarized perpendicularly to the production plane, as allowed by parity invariance; a huge amount of experimental information, for a wide energy range of the unpolarized beams, is available on such single spin asymmetries :

$$P_\mathrm{\Lambda }=\frac{d\sigma ^{AB\mathrm{\Lambda }^{}X}d\sigma ^{AB\mathrm{\Lambda }^{}X}}{d\sigma ^{AB\mathrm{\Lambda }^{}X}+d\sigma ^{AB\mathrm{\Lambda }^{}X}}$$
(1)
Similar effects have been observed for several other hyperons, but we shall consider here only $`\mathrm{\Lambda }`$’s and $`\overline{\mathrm{\Lambda }}`$’s.

Despite the wealth of data and the many years they have been known, no convincing theoretical explanation or understanding of the phenomenon exist . The perturbative QCD dynamics forbids any sizeable single spin asymmetry at the partonic level ; the polarization of hyperons resulting from the strong interaction of unpolarized hadrons must then originate from nonperturbative features, presumably in the hadronization process. A number of models attempting some understanding of the data in this perspective - only achieve partial explanations.

In the last years other single spin asymmetries observed in $`p^{}p\pi X`$ reactions  have attracted a lot of theoretical activity -; a phenomenological description of such asymmetries appears now possible with the introduction of new distribution  and/or fragmentation  functions which are spin and $`𝒌_{}`$ dependent; $`𝒌_{}`$ denotes either the transverse momentum of a quark inside a nucleon or of a hadron with respect to the fragmenting quark.

In particular the effect first discussed by Collins  – that is, the azimuthal angle dependence of the number of hadrons produced in the fragmentation of a transversely polarized quark – has been recently observed ; were such results confirmed the role of these new fragmentation functions would be of great phenomenological importance.

We consider here an effect similar to that suggested by Collins, namely a spin and $`𝒌_{}`$ dependence in the fragmentation of an unpolarized quark into a polarized hadron: a function describing this mechanism was first introduced in Ref.  and denoted by $`D_{1T}^{}`$. This function is introduced in a frame defined by two light-like four-vectors $`n_+`$ and $`n_{}`$, satisfying $`n_+n_{}`$ = 1, and by the plane transverse to them. The four-momentum $`P`$ of the outgoing hadron – a $`\mathrm{\Lambda }`$ hyperon in the present investigation – is in the $`n_{}`$ direction (up to a mass term correction). The function $`D_{1T}^{}`$ is then defined as (displayed in the $`n_+A=0`$ gauge)

$`{\displaystyle \frac{ϵ_T^{ij}k_{Ti}S_{Tj}}{M_h}}D_{1T}^{}(z,k_{}){\displaystyle \underset{X}{}}{\displaystyle \frac{dy^+d^2𝒚_T}{4z(2\pi )^3}e^{iky}}`$ (2)
$`\times \mathrm{Tr}0|\psi (y)|P,S_T;XP,S_T;X|\overline{\psi }(0)\gamma ^{}|0|_{y^{}=0},`$
where the final state depends on the transverse part ($`S_T`$) of the spin vector $`S`$ of the produced $`\mathrm{\Lambda }`$ only, i.e. one should interpret it as $`|P,S_T;X(|P,S=S_T;X|P,S=S_T;X)/2`$, such that the contribution from unpolarized fragmentation cancels out. Furthermore, $`k_{}=|𝒌_{}|`$ is the modulus of the transverse momentum of the hadron in a frame where the fragmenting quark has no transverse momentum. More details on this type of definition of fragmentation (or decay) functions can be found in Refs. .

In the notations of Ref.  a similar function is defined as:

$`\mathrm{\Delta }^ND_{h^{}/a}(z,𝒌_{})`$ $``$ $`\widehat{D}_{h^{}/a}(z,𝒌_{})\widehat{D}_{h^{}/a}(z,𝒌_{})`$
$`=`$ $`\widehat{D}_{h^{}/a}(z,𝒌_{})\widehat{D}_{h^{}/a}(z,𝒌_{}),`$
and denotes the difference between the density numbers $`\widehat{D}_{h^{}/a}(z,𝒌_{})`$ and $`\widehat{D}_{h^{}/a}(z,`$ $`𝒌_{})`$ of spin 1/2 hadrons $`h`$, with longitudinal momentum fraction $`z`$, transverse momentum $`𝒌_{}`$ and transverse polarization $``$ or $``$, inside a jet originated by the fragmentation of an unpolarized parton $`a`$. From the above definition it is clear that the $`𝒌_{}`$ integral of the function vanishes.

The exact relation between $`D_{1T}^{}`$ and $`\mathrm{\Delta }^ND_{h^{}/a}`$ is given by (notice that also $`D_{1T}^{}`$ should have labels $`h`$ and $`a`$ which are often omitted):

$$\mathrm{\Delta }^ND_{h^{}/a}(z,𝒌_{})=2\frac{k_{}}{zM_h}\mathrm{sin}\varphi D_{1T}^{}(z,k_{}),$$
(4)
where $`\varphi `$ is the angle between $`𝒌_{}`$ and the transverse polarization vector of the hadron, which shows that the function $`\mathrm{\Delta }^ND_{h^{}/a}(z,𝒌_{})`$ vanishes in case the transverse momentum and transverse spin are parallel.

In the sequel we shall refer to $`\mathrm{\Delta }^ND_{h^{}/a}`$ and $`D_{1T}^{}`$ as “polarizing fragmentation functions”.

In analogy to Collins’ suggestion for the fragmentation of a transversely polarized quark we write :

$$\widehat{D}_{h^{}/q}(z,𝒌_{})=\frac{1}{2}\widehat{D}_{h/q}(z,k_{})+\frac{1}{2}\mathrm{\Delta }^ND_{h^{}/q}(z,k_{})\frac{\widehat{𝑷}_h(𝒑_q\times 𝒌_{})}{|𝒑_q\times 𝒌_{}|}$$
(5)
for an unpolarized quark with momentum $`𝒑_q`$ which fragments into a spin 1/2 hadron $`h`$ with momentum $`𝒑_h=z𝒑_q+𝒌_{}`$ and polarization vector along the $`=\widehat{𝑷}_h`$ direction; $`\widehat{D}_{h/q}(z,k_{})=D_{h^{}/q}(z,𝒌_{})+D_{h^{}/q}(z,𝒌_{})`$ is the $`k_{}`$ dependent unpolarized fragmentation function. Notice that $`\widehat{𝑷}_h(𝒑_q\times 𝒌_{})=𝒑_q(𝒌_{}\times \widehat{𝑷}_h)\mathrm{sin}\varphi `$.

A QCD factorization theorem gives for the high energy and large $`p_T`$ process $`AB\mathrm{\Lambda }^{}X`$, at leading twist with collinear parton configurations:

$$d\sigma ^{AB\mathrm{\Lambda }^{}X}=\underset{a,b,c,d}{}f_{a/A}(x_a)f_{b/B}(x_b)d\widehat{\sigma }^{abcd}D_{\mathrm{\Lambda }^{}/c}(z)$$
(6)
and

$$d\sigma ^{AB\mathrm{\Lambda }^{}X}=\underset{a,b,c,d}{}f_{a/A}(x_a)f_{b/B}(x_b)d\widehat{\sigma }^{abcd}D_{\mathrm{\Lambda }^{}/c}(z).$$
(7)
Here and in the sequel we shall fix the scattering plane as the $`x`$$`z`$ plane, with hadron $`A`$ moving along $`+\widehat{𝒛}`$ and the detected $`\mathrm{\Lambda }`$ produced in the first $`x`$$`z`$ quadrant; the $``$, $``$ directions are then respectively $`+\widehat{𝒚}`$ and $`\widehat{𝒚}`$.

In the absence of intrinsic $`𝒌_{}`$ (or rather when integrated over) the fragmentation functions $`D_{\mathrm{\Lambda }^{}/c}(z)=d^2𝒌_{}\widehat{D}_{\mathrm{\Lambda }^{}/q}(z,𝒌_{})`$ (or $`D_{\mathrm{\Lambda }^{}/c}(z)`$) cannot depend on the hadron polarization, so that one has $`d\sigma ^{}=d\sigma ^{}`$, which implies $`P_\mathrm{\Lambda }=0`$.

However, by taking into account intrinsic $`𝒌_{}`$ in the hadronization process, and assuming that the factorization theorem holds also when $`𝒌_{}`$’s are included , one has, using Eq. (5) instead of $`D_{\mathrm{\Lambda }/c}(z)`$ in Eqs. (1), (6) and (7):

$`{\displaystyle \frac{E_\mathrm{\Lambda }d^3\sigma ^{AB\mathrm{\Lambda }X}}{d^3𝒑_\mathrm{\Lambda }}}P_\mathrm{\Lambda }`$ $`=`$ $`{\displaystyle \underset{a,b,c,d}{}}{\displaystyle \frac{dx_adx_bdz}{\pi z^2}d^2𝒌_{}f_{a/A}(x_a)f_{b/B}(x_b)}`$ (8)
$`\times `$ $`\widehat{s}\delta (\widehat{s}+\widehat{t}+\widehat{u}){\displaystyle \frac{d\widehat{\sigma }^{abcd}}{d\widehat{t}}}(x_a,x_b,𝒌_{})\mathrm{\Delta }^ND_{\mathrm{\Lambda }^{}/c}(z,𝒌_{})`$
where $`\widehat{s},\widehat{t}`$ and $`\widehat{u}`$ are the Mandelstam variables for the elementary process: for collinear configurations $`\widehat{s}=x_ax_bs,\widehat{t}=x_at/z`$ and $`\widehat{u}=x_bu/z`$ and the modifications due to intrinsic $`𝒌_{}`$ will be taken into account in the numerical evaluations. $`E_\mathrm{\Lambda }d^3\sigma ^{AB\mathrm{\Lambda }X}/d^3𝒑_\mathrm{\Lambda }`$ is the unpolarized cross-section:

$`{\displaystyle \frac{E_\mathrm{\Lambda }d^3\sigma ^{AB\mathrm{\Lambda }X}}{d^3𝒑_\mathrm{\Lambda }}}`$ $`=`$ $`{\displaystyle \underset{a,b,c,d}{}}{\displaystyle \frac{dx_adx_bdz}{\pi z^2}d^2𝒌_{}f_{a/A}(x_a)f_{b/B}(x_b)}`$ (9)
$`\times `$ $`\widehat{s}\delta (\widehat{s}+\widehat{t}+\widehat{u}){\displaystyle \frac{d\widehat{\sigma }^{abcd}}{d\widehat{t}}}(x_a,x_b,𝒌_{})\widehat{D}_{\mathrm{\Lambda }/c}(z,k_{}).`$
In Eq. (8) $`𝒌_{}`$ is considered only where its absence would lead to zero polarization: that is, leading collinear configurations are assumed for partons $`a`$ and $`b`$ inside unpolarized hadrons $`A`$ and $`B`$, while transverse motion is considered in the fragmentation process. The final hadron, detected with momentum $`𝒑_\mathrm{\Lambda }`$, is generated by the hadronization of a parton $`c`$ whose momentum, $`𝒑_c=(𝒑_\mathrm{\Lambda }𝒌_{})/z`$, varies with $`𝒌_{}`$; also the corresponding elementary process, $`abcd`$, depends on $`𝒌_{}`$.

$`P_\mathrm{\Lambda }`$ is a function of the hyperon momentum $`𝒑_\mathrm{\Lambda }=𝒑_L+𝒑_T`$ and is normally measured in the $`AB`$ c.m. frame as a function of $`x_F2p_L/\sqrt{s}`$ and $`p_T`$.

Notice that, in principle, there might be another contribution to the polarization of a final hadron produced at large $`p_T`$ in the high energy collision of two unpolarized hadrons; in analogy to Sivers’ effect  one might introduce a new spin and $`𝒌_{}`$ dependent distribution function:

$`\mathrm{\Delta }^Nf_{a^{}/A}(x_a,𝒌_a)`$ $``$ $`\widehat{f}_{a^{}/A}(x_a,𝒌_a)\widehat{f}_{a^{}/A}(x_a,𝒌_a)`$
$`=`$ $`\widehat{f}_{a^{}/A}(x_a,𝒌_a)\widehat{f}_{a^{}/A}(x_a,𝒌_a),`$
i.e. the difference between the density numbers $`\widehat{f}_{a^{}/A}(x_a,𝒌_a)`$ and $`\widehat{f}_{a^{}/A}(x_a,`$ $`𝒌_a)`$ of partons $`a`$, with longitudinal momentum fraction $`x_a`$, transverse momentum $`𝒌_a`$ and transverse polarization $``$ or $``$, inside an unpolarized hadron $`A`$.

This idea was first applied to unpolarized lepto-production  and to single spin asymmetries in $`pp^{}`$ scattering ; the corresponding function, related to $`\mathrm{\Delta }^Nf_{a^{}/A}(x_a,𝒌_a)`$, was denoted by $`h_1^{}`$. In the present case of transversely polarized $`\mathrm{\Lambda }`$ production this function would enter the cross-section accompanied by the transversity fragmentation function $`H_1(z)`$ (or $`\mathrm{\Delta }D_{h^{}/a^{}}`$). We shall not consider this contribution here; not only because of the problems concerning $`\mathrm{\Delta }^Nf_{a^{}/A}(x_a,𝒌_a)`$ discussed below, but also because the experimental evidence that the hyperon polarization is somewhat independent of the nature of the hadronic target suggests that the mechanism responsible for the polarization is in the hadronization process111This is not in contradiction with the observed spin transfer ($`D_{NN}`$) in $`pp^{}\mathrm{\Lambda }^{}X`$ , which in the factorized approach can be described in terms of the ordinary transversity distribution and fragmentation functions, rather than in terms of $`\mathrm{\Delta }^ND_{h^{}/a}`$.. A clean test of this should come from a measurement of $`P_\mathrm{\Lambda }`$ in unpolarized DIS processes, $`\mathrm{}p\mathrm{\Lambda }^{}X`$ .

The $`𝒌_{}`$ dependent functions considered in Refs.  ($`\mathrm{\Delta }^Nf_{a/A^{}}`$, $`\mathrm{\Delta }^Nf_{a^{}/A}`$, $`\mathrm{\Delta }^ND_{h/a^{}}`$ and $`\mathrm{\Delta }^ND_{h^{}/a}`$, or, respectively, $`f_{1T}^{}`$, $`h_1^{}`$, $`H_1^{}`$ and $`D_{1T}^{}`$) have the common feature that the transverse momentum direction is correlated with the direction of the transverse spin of either a quark or a hadron, via a $`\mathrm{sin}\varphi `$ dependence, as can be seen from Eq. (2) for example. The reason for considering these functions is that this “handedness” of the transverse momentum compared to the transverse spin is displayed by the single spin asymmetry data in, for instance, $`pp^{}\pi X`$ and $`pp\mathrm{\Lambda }^{}X`$. However, the problem of using such functions is that naively they appear to be absent due to time reversal invariance. This conclusion would be valid if the hadronic state appearing in the definition of such functions is treated as a plane wave state. One can then show that the functions are odd under the application of time reversal invariance, whereas hermiticity requires them to be even. If, however, initial or final state interactions are present, then time reversal symmetry will not prevent the appearance of nonzero (naively) T-odd functions. In the case of a state like $`|P_h,S_h;X`$ final state interactions are certainly present and nonzero fragmentation functions $`\mathrm{\Delta }^ND_{h/a^{}}`$ and $`\mathrm{\Delta }^ND_{h^{}/a}`$ are expected. However, for distribution functions this issue poses severe problems and since we will only consider fragmentation functions here, we refer to Refs.  for more detailed discussions on this topic.

The main difference between the function $`\mathrm{\Delta }^ND_{h/a^{}}`$ as originally proposed by Collins, and the function under present investigation $`\mathrm{\Delta }^ND_{h^{}/a}`$, is that the former is a so-called chiral-odd function, which means that it couples quarks with left- and right-handed chiralities, whereas the latter function is chiral-even. Since the pQCD interactions conserve chirality, chiral-odd functions must always be accompanied by a mass term or appear in pairs. Both options restrict the accessibility of such functions and make them harder to be determined separately. On the other hand, the chiral-even fragmentation function can simply occur accompanied by the unpolarized (chiral-even) distribution functions, which are the best determined quantities, allowing for a much cleaner extraction of the fragmentation function itself. Moreover, since chiral-even functions can appear in charged current mediated processes (as opposed to chiral-odd functions), more methods of extraction are available .

As it was studied in Ref.  the Collins function $`H_1^{}`$ (or $`\mathrm{\Delta }^ND_{h/a^{}}`$) satisfies a sum rule arising from momentum conservation in the transverse directions. The same holds for the other $`𝒌_{}`$-odd, T-odd fragmentation function $`D_{1T}^{}`$ ,

$$\underset{h}{}𝑑z𝑑k_{}^2\frac{k_{}^2}{zM_h}D_{1T}^{}(z,k_{})=0,$$
(11)
or, in terms of the function $`\mathrm{\Delta }^ND_{h^{}/a}`$,

$$\underset{h}{}C_hM_h\underset{h}{}𝑑zd^2𝒌_{}k_{}\mathrm{sin}\varphi \mathrm{\Delta }^ND_{h^{}/a}(z,𝒌_{})=0,$$
(12)
which is equivalent to Eq. (11) via Eq. (4).

Notice that the above sum rule can be written as

$$\underset{h}{}𝑑zd^2𝒌_{}𝒌_{}\widehat{D}_{h^{}/a}(z,𝒌_{})=0,$$
(13)
and the same holds independently for $`\widehat{D}_{h^{}/a}(z,𝒌_{})`$. Eq. (13) has a clear nontrivial physical meaning: for each polarization direction ($``$ or $``$) the total transverse momentum carried by spin 1/2 hadrons222Strictly speaking, the sum over $`h`$ is over all hadrons that carry transverse polarization, which might be true also for higher spin hadrons. is zero.

The sum over hadrons prevents a straightforward application of the sum rule to the case of $`\mathrm{\Lambda }`$ production alone. It can not be used as a constraint on the parameterization of the function to be fitted to the data. However, we note that the integral $`C_hM_h`$ for each hadron type $`h`$ is the same function of the energy scale (implicit in all expressions), apart from as yet unknown normalization. In this sense it closely resembles the tensor charge. In other words, the running of the functions are the same for any type of hadron and there is no mixing with other functions. The ratios for different types of hadrons are constants, which allow for checks of consistency between sets of data obtained at different energies, without the need for evolution. These constants are universal, if indeed the T-odd fragmentation functions are universal. This universality is the main point of interest here: one wants to see if $`\mathrm{\Lambda }`$ polarization from unpolarized quark fragmentation is independent of the initial state, as it is implicit when writing down the factorized cross-sections Eqs. (8) and (9). At the present time, this can not be verified due to lack of data, but some predictions can be given  that will allow tests of this universality.

We only consider the quark fragmenting into a $`\mathrm{\Lambda }`$ and not into other hyperons, like the $`\mathrm{\Sigma }^0`$. The latter actually produces a significant amount (20-30%) of the $`\mathrm{\Lambda }`$’s via the decay $`\mathrm{\Sigma }^0\mathrm{\Lambda }\gamma `$. The reason we do not introduce a separate $`\mathrm{\Sigma }^0`$ fragmentation function at this stage is that the factorized approach by itself does not address such a separation (it is about a generic spin-1/2 hadron of type $`h`$), unless one introduces some additional input, like a model based on $`SU(3)`$, or unless one applies it to separate sets of data for each hyperon (which are not available yet). Apart from that, for each quark flavour such a $`\mathrm{\Sigma }^0`$ fragmentation function would evolve in the same way as the $`\mathrm{\Lambda }`$ fragmentation function, which implies that their relative fraction stays constant. In this way we can view the $`\mathrm{\Lambda }`$ fragmentation function as an effective fragmentation function that includes the contamination due to $`\mathrm{\Sigma }^0`$’s. Indeed, the fragmentation functions we shall use in next Section have been obtained from fits to inclusive $`\mathrm{\Lambda }`$ productions in $`e^+e^{}`$ processes, independently of their origin.

At a later stage one might make the distinction that the $`\mathrm{\Sigma }^0`$ fragmentation is a different fraction of the effective $`\mathrm{\Lambda }`$ fragmentation function for different quark flavours. One can insert all this information with hindsight and correct for it, but the present approach cannot be used to acquire this information unless the data would clearly distinguish the $`\mathrm{\Lambda }`$’s coming from $`\mathrm{\Sigma }^0`$’s. Our approach of using an effective $`\mathrm{\Lambda }`$ fragmentation function would be more problematic if the $`\mathrm{\Sigma }^0`$ would decay into other final states than $`\mathrm{\Lambda }\gamma `$ (which branching ratio happens to be 100%): then only a part of the total $`\mathrm{\Sigma }^0`$ fragmentation function would be included into the effective $`\mathrm{\Lambda }`$ fragmentation function and this would be energy dependent.

In the case of longitudinally polarized $`\mathrm{\Lambda }`$ production the $`\mathrm{\Sigma }^0\mathrm{\Lambda }\gamma `$ background gives rise to a depolarizing effect of about 10% , but in the present situation of transversely polarized $`\mathrm{\Lambda }`$ production this is not the case, since the photon does not carry away any of the transverse polarization and it hardly affects the definition of the plane compared to which the transverse polarization is measured. Therefore, the $`\mathrm{\Sigma }^0\mathrm{\Lambda }\gamma `$ background does not produce a significant depolarizing effect for the transverse $`\mathrm{\Lambda }`$ polarization and an effective $`\mathrm{\Lambda }`$ polarizing fragmentation function can be used also.

We shall now consider both $`\mathrm{\Lambda }`$ and $`\overline{\mathrm{\Lambda }}`$ production and attempt a determination of the polarizing fragmentation functions $`\mathrm{\Delta }^ND_{\mathrm{\Lambda }^{}/q}`$ by comparing results for $`P_\mathrm{\Lambda }`$ and $`P_{\overline{\mathrm{\Lambda }}}`$ from Eqs. (8) and (9) with data -.

2. Numerical fits of data on $`P_𝚲`$ from $`pN\mathbf{}𝚲X`$ processes Eq. (8), for proton-nucleon processes, can be rewritten as:

$`{\displaystyle \frac{E_\mathrm{\Lambda }d^3\sigma ^{pN\mathrm{\Lambda }X}}{d^3𝒑_\mathrm{\Lambda }}}P_\mathrm{\Lambda }={\displaystyle \underset{a,b,c,d}{}}{\displaystyle _{(+k_{})}}d^2𝒌_{}[{\displaystyle _{z_{\mathrm{min}}}^1}dz{\displaystyle _{x_{a\mathrm{min}}}^1}dx_a{\displaystyle \frac{1}{\pi z}}{\displaystyle \frac{\overline{x}_b^2s}{(t\mathrm{\Phi }_t)}}`$ (14)
$`\times f_{a/p}(x_a)f_{b/N}(\overline{x}_b){\displaystyle \frac{d\widehat{\sigma }^{abcd}}{d\widehat{t}}}(x_a,\overline{x}_b,𝒌_{})\{𝒌_{}𝒌_{}\}]\mathrm{\Delta }^ND_{\mathrm{\Lambda }^{}/c}(z,𝒌_{}),`$
which deserves several comments and some explanation of notations.

* In deriving Eq. (14) from Eq. (8) we have used the fact that $`\mathrm{\Delta }^ND_{\mathrm{\Lambda }^{}/c}(z,𝒌_{})`$, Eq. (S0.Ex2), is an odd function of $`𝒌_{}`$; the $`(+k_{})`$ integration region of $`𝒌_{}`$ runs over one half-plane of its components.
* The $`x_b`$ integration has been performed by exploiting the $`\delta (\widehat{s}+\widehat{t}+\widehat{u})`$ function; the resulting value of $`x_b`$ is given by

$$\overline{x}_b=\frac{x_at\mathrm{\Phi }_t}{x_azs+u\mathrm{\Phi }_u},$$
(15)
where $`\mathrm{\Phi }_t`$ and $`\mathrm{\Phi }_u`$ are defined below.
* $`𝒌_{}`$ could have any direction in the plane perpendicular to $`𝒑_c`$; however, due to parity conservation in the hadronization process – Eq. (5) – the only $`𝒌_{}`$ component contributing to the polarizing fragmentation function is that perpendicular to $`𝑷_\mathrm{\Lambda }`$, i.e. the component lying in the production plane, the $`xz`$ plane in our conventions. To simplify the kinematics we shall then consider only the leading contributions of $`𝒌_{}`$ vectors in the $`xz`$ plane.
* $`s,t`$ and $`u`$ are the Mandelstam variables for the $`pN\mathrm{\Lambda }X`$ process; in the simple planar configuration discussed above ($`𝒑_c`$ and $`𝒌_{}`$ both lying in the $`xz`$ production plane) they are related to the corresponding variables for the elementary processes by:

$`\widehat{s}`$ $`=`$ $`2p_ap_b=x_ax_bs`$
$`\widehat{t}`$ $`=`$ $`2p_ap_c=(x_a/z)t\mathrm{\Phi }_t(\pm k_{})`$ (16)
$`\widehat{u}`$ $`=`$ $`2p_bp_c=(x_b/z)u\mathrm{\Phi }_u(\pm k_{})`$
with

$`t\mathrm{\Phi }_t(\pm k_{})=g(k_{})\left\{t2k_{}{\displaystyle \frac{\sqrt{stu}}{t+u}}[1g(k_{})]{\displaystyle \frac{tu}{2}}\right\}`$ (17)
$`u\mathrm{\Phi }_u(\pm k_{})=g(k_{})\left\{u\pm 2k_{}{\displaystyle \frac{\sqrt{stu}}{t+u}}+[1g(k_{})]{\displaystyle \frac{tu}{2}}\right\}`$ (18)
where $`g(k_{})=\sqrt{1(k_{}/p_\mathrm{\Lambda })^2}`$ and where $`\pm k_{}`$ refers respectively to the configuration in which $`𝒌_{}`$ points to the left or to the right of $`𝒑_c`$. At leading order in $`k_{}/p_\mathrm{\Lambda }`$ one has:

$$\mathrm{\Phi }_t(k_{})=1\frac{k_{}}{p_\mathrm{\Lambda }}\sqrt{\frac{u}{t}}\mathrm{\Phi }_u(k_{})=1+\frac{k_{}}{p_\mathrm{\Lambda }}\sqrt{\frac{t}{u}}$$
(19)
* The lower integration limits in Eq. (14) are given by:

$$x_{a\mathrm{min}}=\frac{u\mathrm{\Phi }_u(\pm k_{})}{zs+t\mathrm{\Phi }_t(\pm k_{})}z\frac{t\mathrm{\Phi }_t(\pm k_{})+u\mathrm{\Phi }_u(\pm k_{})}{s}$$
(20)
Notice that the integration limits are slightly different for $`+k_{}`$ and $`k_{}`$; when replacing $`𝒌_{}`$ with $`𝒌_{}`$ inside the square bracket of Eq. (14), one should not forget to change accordingly also the $`z`$ and $`x_a`$ integration limits, and the value of $`\overline{x}_b`$, Eq. (15), although the results are only marginally affected by this.
* Eq. (14) can be schematically written as

$`d\sigma ^{pN\mathrm{\Lambda }X}P_\mathrm{\Lambda }=d\sigma ^{pN\mathrm{\Lambda }^{}X}d\sigma ^{pN\mathrm{\Lambda }^{}X}={\displaystyle \underset{a,b,c,d}{}}f_{a/p}(x_a)f_{b/N}(x_b)`$ (21)
$``$ $`[d\widehat{\sigma }^{abcd}(x_a,x_b,𝒌_{})d\widehat{\sigma }^{abcd}(x_a,x_b,𝒌_{})]\mathrm{\Delta }^ND_{\mathrm{\Lambda }^{}/c}(z,𝒌_{})`$
which shows clearly that $`P_\mathrm{\Lambda }`$ is a higher twist effect, despite the fact that the polarizing fragmentation function $`\mathrm{\Delta }^ND_{h^{}/a}`$ is a leading twist function: this is due to the difference in the square brackets, $`[d\widehat{\sigma }(+𝒌_{})d\widehat{\sigma }(𝒌_{})]k_{}/p_T`$, similarly to what happens for the single spin asymmetries in $`p^{}p\pi X`$ .
* In the computation of the unpolarized cross-section $`E_\mathrm{\Lambda }d^3\sigma ^{pN\mathrm{\Lambda }X}/d^3𝒑_\mathrm{\Lambda }`$ intrinsic transverse motion is significant only in limited phase space regions: we have checked that the values obtained in most of the kinematical regions of available data do not vary whether or not taking into account $`𝒌_{}`$. Notice that when taking into account $`𝒌_{}`$ the expression for $`E_\mathrm{\Lambda }d^3\sigma ^{pN\mathrm{\Lambda }X}/d^3𝒑_\mathrm{\Lambda }`$ is the same as Eq. (14) with the $``$ inside the square brackets replaced by a $`+`$ and $`\mathrm{\Delta }^ND_{\mathrm{\Lambda }^{}/c}(z,𝒌_{})`$ replaced by $`\widehat{D}_{\mathrm{\Lambda }/c}(z,k_{})`$.
* When computing the cross-sections for scattering off nuclei, $`pA\mathrm{\Lambda }^{}X`$, for which plenty of data are available, we have adopted the most simple incoherent sum, neglecting nuclear effects. That is, for the scattering off a nucleus with A nucleons and Z protons we use:

$$d\sigma ^{pA\mathrm{\Lambda }X}=Zd\sigma ^{pp\mathrm{\Lambda }X}+(AZ)d\sigma ^{pn\mathrm{\Lambda }X}.$$
(22)
We have checked that different ways of taking into account nuclear effects leave results for $`P_\mathrm{\Lambda }`$ – a ratio of cross-sections – almost unchanged. The partonic distribution functions in a neutron are obtained from the usual proton distribution functions by applying isospin invariance.
* Eq. (14) holds for any spin 1/2 baryon; we shall use it also for $`\overline{\mathrm{\Lambda }}`$’s, using charge conjugation invariance to obtain $`\overline{q}\overline{\mathrm{\Lambda }}`$ fragmentation properties from $`q\mathrm{\Lambda }`$ ones, which implies $`D_{\overline{\mathrm{\Lambda }}/\overline{q}}=D_{\mathrm{\Lambda }/q}`$ and $`\mathrm{\Delta }^ND_{\overline{\mathrm{\Lambda }}^{}/\overline{q}}=\mathrm{\Delta }^ND_{\mathrm{\Lambda }^{}/q}`$.

We now use Eq. (14) in order to see whether or not it can reproduce the data and in order to obtain information on the new polarizing fragmentation functions. To do so we introduce a simple parameterization for these functions and fix the parameters by fitting the existing data on $`P_\mathrm{\Lambda }`$ and $`P_{\overline{\mathrm{\Lambda }}}`$ -.

We assume that $`\mathrm{\Delta }^ND_{\mathrm{\Lambda }^{}/c}(z,𝒌_{})`$ is strongly peaked around an average value $`𝒌_{}^0`$ lying in the production plane, so that we can expect:

$$_{(+k_{})}d^2𝒌_{}\mathrm{\Delta }^ND_{\mathrm{\Lambda }^{}/c}(z,𝒌_{})F(𝒌_{})\mathrm{\Delta }_0^ND_{\mathrm{\Lambda }^{}/c}(z,k_{}^0)F(𝒌_{}^0).$$
(23)
Consistently, since in this case $`F(𝒌_{})`$ depends weakly on $`𝒌_{}`$ when $`k_{}p_T`$, in the computation of the unpolarized cross-section we use:

$$_{(+k_{})}d^2𝒌_{}\widehat{D}_{\mathrm{\Lambda }/c}(z,k_{})F(𝒌_{})\frac{1}{2}D_{\mathrm{\Lambda }/c}(z)F(𝒌_{}^0).$$
(24)
The average $`k_{}^0`$ value depends on $`z`$ and we parameterize this dependence in a most natural way:

$$\frac{k_{}^0(z)}{M}=Kz^a(1z)^b,$$
(25)
where $`M`$ is a momentum scale ($`M=1`$ GeV/$`c`$); in performing the fit we demand that $`K`$, $`a`$ and $`b`$ are constrained so that they satisfy the kinematical bound $`p_q^2=(p_\mathrm{\Lambda }^2k_{}^2)/z^2p_\mathrm{\Lambda }^2`$, from which $`k_{}^2(1z^2)p_\mathrm{\Lambda }^2`$ and

$$k_{}^0(z)\stackrel{<}{}(p_\mathrm{\Lambda })_{\mathrm{min}}\sqrt{1z}(1\mathrm{GeV}/c)\sqrt{1z},$$
(26)
which implies $`a0`$ and $`b0.5`$. The values of $`K`$, $`a`$ and $`b`$ resulting from our best fit will have a clear physical meaning.

We parameterize $`\mathrm{\Delta }_0^ND_{\mathrm{\Lambda }^{}/c}(z,k_{}^0)`$ in a similar simple form: we know that it has to be zero when $`k_{}=0`$ and $`z=1`$; in addition, the positivity of the fragmentation functions – Eq. (5) – requires, at any $`k_{}`$ value, $`|\mathrm{\Delta }^ND_{h^{}/q}(z,k_{})|\widehat{D}_{h/q}(z,k_{})`$. However, according to Eqs. (23), (24) and to take into account the mentioned \[see comment after Eq. (20)\] difference of the integration regions $`(+k_{})`$ and $`(k_{})`$, which is significant at the boundaries of the kinematical ranges (when $`p_Tk_{}`$ and when $`p_Tp_{T\mathrm{max}}`$) we prefer to impose the more stringent bound $`|\mathrm{\Delta }_0^ND_{\mathrm{\Lambda }^{}/c}(z,k_{}^0)|D_{\mathrm{\Lambda }/c}(z)/2`$. Following Ref. , this is automatically satisfied by taking:

$`\mathrm{\Delta }_0^ND_{\mathrm{\Lambda }^{}/q}(z,k_{}^0)`$ $`=`$ $`N_q^{}{\displaystyle \frac{k_{}^0(z)}{M}}\left[{\displaystyle \frac{z^{\alpha _q}(1z)^{\beta _q}}{c_q^{c_q}d_q^{d_q}/(c_q+d_q)^{c_q+d_q}}}\right]{\displaystyle \frac{D_{\mathrm{\Lambda }/q}(z)}{2}}`$
$`=`$ $`N_q^{}K{\displaystyle \frac{z^{c_q}(1z)^{d_q}}{c_q^{c_q}d_q^{d_q}/(c_q+d_q)^{c_q+d_q}}}{\displaystyle \frac{D_{\mathrm{\Lambda }/q}(z)}{2}}`$
$``$ $`N_qz^{c_q}(1z)^{d_q}{\displaystyle \frac{D_{\mathrm{\Lambda }/q}(z)}{2}},`$
where we have used Eq. (25), and we require $`c_q=a+\alpha _q>0`$, $`d_q=b+\beta _q>0`$, and $`|N_q^{}|K1`$.

We are now almost fully equipped to compute $`P_\mathrm{\Lambda }`$ and $`P_{\overline{\mathrm{\Lambda }}}`$; let us briefly discuss the remaining quantities appearing in Eq. (14).

1. We sum over all possible elementary interactions, computed at lowest order in pQCD. The polarizing fragmentation functions – parameterized as in Eq. (S0.Ex11) – are supposed to be non vanishing only for $`\mathrm{\Lambda }`$ valence quarks, $`u`$, $`d`$ and $`s`$; similarly for $`\overline{\mathrm{\Lambda }}`$ valence antiquarks $`\overline{u}`$, $`\overline{d}`$ and $`\overline{s}`$. All contributions to the unpolarized fragmentation functions – from quarks, antiquarks and gluons – are taken into account.
2. We adopt the unpolarized distribution functions of MRST . We have explicitly checked that a different choice makes no difference in our conclusions. We fix the QCD hard scale of distribution (and fragmentation) functions at 2 (GeV$`/c)^2`$, corresponding to an average value $`p_T1.5`$ GeV$`/c`$. Since the range of $`p_T`$ values of the data is rather limited, no evolution effect would be visible anyway.
3. We use the set of unpolarized fragmentation functions of Ref. , which allows a separate determination of $`D_{\mathrm{\Lambda }/q}`$ and $`D_{\overline{\mathrm{\Lambda }}/q}`$, and which includes $`\mathrm{\Lambda }`$’s both directly and indirectly produced; it also differentiates between the $`s`$ quark contribution and the $`u`$ and $`d`$ ones: the non strange fragmentation functions $`D_{\mathrm{\Lambda }/u}=D_{\mathrm{\Lambda }/d}`$ are suppressed by an $`SU(3)`$ symmetry breaking factor $`\lambda =0.07`$ as compared to $`D_{\mathrm{\Lambda }/s}`$. In our parameterization of $`\mathrm{\Delta }_0^ND_{\mathrm{\Lambda }^{}/q}(z,k_{}^0)`$, Eq. (S0.Ex11), we use the same $`D_{\mathrm{\Lambda }/q}`$ as given in Ref. , keeping the same parameters $`\alpha _q`$ and $`\beta _q`$ ($`c_q`$ and $`d_q`$) for all quark flavours, but allowing for different values of $`N_u=N_d`$ and $`N_s`$.

We can now use pQCD dynamics, together with the chosen distribution and fragmentation functions, and the parameterized expressions of the polarizing fragmentation functions, in Eq. (14) to compute $`P_\mathrm{\Lambda }`$ and $`P_{\overline{\mathrm{\Lambda }}}`$; by comparing with data we obtain the best fit values of the parameters $`K`$, $`a`$, $`b`$, $`N_u=N_d`$, $`N_s`$, $`c_q`$ and $`d_q`$ introduced in Eqs. (25) and (S0.Ex11). Notice that we remain with 7 free parameters, $`c_q`$ and $`d_q`$ being the same for all flavours.

Our best fit results ($`\chi ^2`$/d.o.f. = 1.57), taking into account all data -, are shown in Figs. 1-5. They correspond to the best fit parameter values:

$$K=0.69a=0.36b=0.53$$
(28)
$$N_u=4.30N_s=1.13c_q=6.58d_q=0.67.$$
(29)
Let us comment in some details on our results.

In Fig. 1 and 2 we present our best fits to $`P_\mathrm{\Lambda }`$ as a function of $`p_T`$ for different $`x_F`$ values, as indicated in the figures: the famous approximately flat $`p_T`$ dependence, for $`p_T`$ greater than 1 GeV/$`c`$, is well reproduced. Such a behaviour, as expected, does not continue indefinitely with $`p_T`$ and we have explicitly checked that at larger values of $`p_T`$ the values of $`P_\mathrm{\Lambda }`$ drop to zero: the shape of such a decrease, contrary to what happens in the $`p_T`$ region of the data shown here, strongly depends on the assumptions about the nuclear corrections. It may be interesting to note that this fall-off has not yet been observed experimentally, but is expected to be first seen in the near-future BNL-RHIC data.

Also the increase of $`|P_\mathrm{\Lambda }|`$ with $`x_F`$ at fixed $`p_T`$ values can be well described, as shown in Fig. 3; the two curves correspond to $`p_T=1.5`$ GeV/$`c`$ (solid) and $`p_T=3`$ GeV/$`c`$ (dashed-dotted).

The best fits of Figs. 1 and 2 are compared to $`p`$$`Be`$ data -; these are collected at two different energies, in the $`p`$$`Be`$ c.m. frame, $`\sqrt{s}82`$ GeV - and $`\sqrt{s}116`$ GeV . Our calculations are performed at $`\sqrt{s}=80`$ GeV; we have explicitly checked that by varying the energy between 80 and 120 GeV, our results for $`P_\mathrm{\Lambda }`$ vary, in the kinematical range considered here, at most by 10%, in agreement with the observed energy independence of the data.

Some data from $`p`$$`p`$ collisions are also available; in Ref.  a linear fit to $`P_\mathrm{\Lambda }(x_F)`$ is performed, collecting all data with $`p_T0.96`$ GeV/$`c`$, for an average value $`p_T=1.1`$ GeV/$`c`$. In Fig. 4 we show these data and the linear fit (central line); the upper and lower lines show the fit error band. The solid line is our computation, at $`p_T=1.1`$ GeV/$`c`$, with all parameters fixed as in Eqs. (28) and (29): our fit reproduces the data with good accuracy.

In Fig. 5 we show our best fit results for $`P_{\overline{\mathrm{\Lambda }}}`$ as a function of $`p_T`$ for different $`x_F`$ values, as indicated in the figure: in this case all data  are compatible with zero, with large errors, and the measured $`x_F`$ range is not as wide as for $`P_\mathrm{\Lambda }`$, as expected from the lack of overlapping between $`\overline{\mathrm{\Lambda }}`$ and nucleon valence quarks.

The resulting values of the parameters, Eqs. (28) and (29), are very realistic; notice in particular that $`b`$ essentially reaches its kinematical limit 0.5 and the whole function (25) giving the average $`k_{}`$ value of a $`\mathrm{\Lambda }`$ inside a jet turns out to be very reasonable.

Mostly $`u`$ and $`d`$ quarks contribute to $`P_\mathrm{\Lambda }`$, resulting in a negative value of $`N_u`$; instead, $`u`$, $`d`$ and $`s`$ quarks all contribute significantly to $`P_{\overline{\mathrm{\Lambda }}}`$ and opposite signs for $`N_u`$ and $`N_s`$ are found, inducing cancellations.

In Fig. 6 we plot $`|\mathrm{\Delta }_0^ND_{\mathrm{\Lambda }^{}/u,d}|`$ and $`\mathrm{\Delta }_0^ND_{\mathrm{\Lambda }^{}/s}`$ as given by Eq. (S0.Ex11) with the best fit parameters (28) and (29). We show, for comparison, also $`D_{\mathrm{\Lambda }/u,d}`$ and $`D_{\mathrm{\Lambda }/s}`$: notice how a tiny value of the polarizing fragmentation functions, in a limited large $`z`$ region, is enough to allow a good description of the data. This also shows that taking into account only valence quark contributions to $`\mathrm{\Delta }^ND_{\mathrm{\Lambda }^{}/q}`$ – as we have done – is a justified assumption.

A different set of unpolarized fragmentation functions can be found in the literature : it holds for the quark fragmentation into $`\mathrm{\Lambda }+\overline{\mathrm{\Lambda }}`$ and gives $`D_{(\mathrm{\Lambda }+\overline{\mathrm{\Lambda }})/q}`$ rather than a separate expression of $`D_{\mathrm{\Lambda }/q}`$ and $`D_{\overline{\mathrm{\Lambda }}/q}`$; it would be appropriate to compute the $`\mathrm{\Lambda }+\overline{\mathrm{\Lambda }}`$ single spin asymmetry:

$`P_{\mathrm{\Lambda }+\overline{\mathrm{\Lambda }}}`$ $`=`$ $`{\displaystyle \frac{d\sigma ^\mathrm{\Lambda }^{}+d\sigma ^{\overline{\mathrm{\Lambda }}^{}}d\sigma ^\mathrm{\Lambda }^{}d\sigma ^{\overline{\mathrm{\Lambda }}^{}}}{d\sigma ^\mathrm{\Lambda }^{}+d\sigma ^{\overline{\mathrm{\Lambda }}^{}}+d\sigma ^\mathrm{\Lambda }^{}+d\sigma ^{\overline{\mathrm{\Lambda }}^{}}}}`$ (30)
$`=`$ $`(P_\mathrm{\Lambda }+{\displaystyle \frac{d\sigma ^{\overline{\mathrm{\Lambda }}}}{d\sigma ^\mathrm{\Lambda }}}P_{\overline{\mathrm{\Lambda }}})(1+{\displaystyle \frac{d\sigma ^{\overline{\mathrm{\Lambda }}}}{d\sigma ^\mathrm{\Lambda }}})^1`$
However, since one knows from experiments on $`pN`$ reactions that in the kinematical range of interest $`d\sigma ^{\overline{\mathrm{\Lambda }}}d\sigma ^\mathrm{\Lambda }`$ (and this is confirmed in our scheme, simply due to the dominance of $`q`$ over $`\overline{q}`$ in a nucleon), one can safely assume

$$P_{\mathrm{\Lambda }+\overline{\mathrm{\Lambda }}}P_\mathrm{\Lambda }.$$
(31)
This set – differently from the one of Ref.  – assumes $`SU(3)`$ symmetry and takes $`D_{\mathrm{\Lambda }/u}=D_{\mathrm{\Lambda }/d}=D_{\mathrm{\Lambda }/s}`$.

We have also computed $`P_\mathrm{\Lambda }`$ with this second set of fragmentation functions; as in the previous case we have parameterized $`\mathrm{\Delta }^ND_{\mathrm{\Lambda }^{}/q}`$ according to Eq. (S0.Ex11), with the same values of $`c_q`$ and $`d_q`$ for each flavour, but different values of $`N_{u,d}`$ and $`N_s`$. We can equally well ($`\chi ^2`$/d.o.f. = 1.85) fit the data on $`P_\mathrm{\Lambda }`$, obtaining fits almost indistinguishable from those of Figs. 1-4, with the best fit values:

$$K=0.66a=0.37b=0.50$$
(32)
$$N_u=28.13N_s=57.53c_q=11.64d_q=1.23$$
(33)
Notice that, also in this case of $`SU(3)`$ symmetric fragmentation functions $`D_{\mathrm{\Lambda }/q}`$, and using only data on $`P_\mathrm{\Lambda }`$, one reaches similar conclusions about the polarizing fragmentation functions $`\mathrm{\Delta }^ND_{\mathrm{\Lambda }^{}/q}`$: $`N_{u,d}N_s`$ and not only is there a difference in magnitude, but once more one finds negative values for $`\mathrm{\Delta }_0^ND_{\mathrm{\Lambda }^{}/u,d}`$ and positive ones for $`\mathrm{\Delta }_0^ND_{\mathrm{\Lambda }^{}/s}`$. This seems to be a well established general trend. Plots analogous to those of Fig. 6, would show also in this case, Eqs. (S0.Ex11) and (33), $`\mathrm{\Delta }_0^ND_{\mathrm{\Lambda }^{}/s}>|\mathrm{\Delta }_0^ND_{\mathrm{\Lambda }^{}/u,d}|`$.

Very recently, new sets of quark and antiquark fragmentation functions into a $`\mathrm{\Lambda }`$, based on a bag model calculation and a fit to $`e^+e^{}`$ data, have been published . Both a $`SU(3)`$ flavour symmetric and a $`SU(3)`$ symmetry breaking set are available, although at a rather too low energy scale ($`\mu =0.25`$ GeV). Nevertheless, we have also tried using these sets in our scheme to fit the data on $`P_\mathrm{\Lambda }`$ and $`P_{\overline{\mathrm{\Lambda }}}`$: once more we can fit the data, better with the symmetric than the asymmetric set, and, again, negative $`\mathrm{\Delta }_0^ND_{\mathrm{\Lambda }^{}/u,d}`$ and positive $`\mathrm{\Delta }_0^ND_{\mathrm{\Lambda }^{}/s}`$ are found, with $`\mathrm{\Delta }_0^ND_{\mathrm{\Lambda }^{}/s}>|\mathrm{\Delta }_0^ND_{\mathrm{\Lambda }^{}/u,d}|`$.

3. Conclusions

A phenomenological approach has been developed towards a consistent explanation and predictions of transverse single spin effects in processes with inclusively produced hadrons; we work in a kinematical region where pQCD and the factorization scheme can be used, but $`p_T`$ is not much larger than intrinsic $`k_{}`$, so that higher twist contributions are still important. This applies to several processes for which data are available, like $`p^{}p\pi X`$  and $`pN\mathrm{\Lambda }^{}X`$. The single spin effect originates in the fragmentation process, either of a transversely polarized quark into an unpolarized hadron – Collins’ effect  – or of an unpolarized quark into a transversely polarized hadron – the polarizing fragmentation functions . Single spin effects in quark distribution functions  have also been discussed .

We have considered here the well known and longstanding problem of the polarization of $`\mathrm{\Lambda }`$ hyperons, produced at large $`p_T`$ in the collision of two unpolarized hadrons. We have assumed a generalized factorization scheme – with the inclusion of intrinsic transverse motion – with pQCD dynamics; the new, spin and $`𝒌_{}`$ dependent, polarizing fragmentation functions $`\mathrm{\Delta }^ND_{\mathrm{\Lambda }^{}/q}`$ have been parameterized in a simple way and data on $`pBe\mathrm{\Lambda }^{}X`$, $`pBe\overline{\mathrm{\Lambda }}^{}X`$ and $`pp\mathrm{\Lambda }^{}X`$ have been used to determine the values of the parameters which give a best fit to the experimental measurements.

The data can be described with remarkable accuracy in all their features: the large negative values of the $`\mathrm{\Lambda }`$ polarization, the increase of its magnitude with $`x_F`$, the puzzling flat $`p_T\stackrel{>}{}1`$ GeV/$`c`$ dependence and the $`\sqrt{s}`$ independence; data from $`p`$-$`p`$ processes are in agreement with data from $`p`$-$`Be`$ interactions and also the tiny or zero values of $`\overline{\mathrm{\Lambda }}`$ polarization are well reproduced. The resulting functions $`\mathrm{\Delta }^ND_{\mathrm{\Lambda }^{}/q}`$ are very reasonable and realistic.

Different sets of unpolarized fragmentation functions – either $`SU(3)`$ symmetric or not – lead to very similar conclusions about these new polarizing fragmentation functions describing the hadronization process of an unpolarized quark into a polarized $`\mathrm{\Lambda }`$: they have opposite signs for $`u`$ and $`d`$ quarks as compared with $`s`$ quarks and their magnitudes are larger for $`s`$ quarks. They are sizeable – with respect to the unpolarized fragmentation functions – only in limited $`z`$ regions, yet they can describe remarkably well the experiments.

These polarizing fragmentation functions have a partonic interpretation and a formal definition, Eq. (2); they are free from the ambiguities related to initial state interactions which might affect analogous distribution functions and we expect them to be universal, process independent functions. Our parameterization of $`\mathrm{\Delta }^ND_{\mathrm{\Lambda }^{}/q}`$ should allow us to give prediction for $`\mathrm{\Lambda }`$ polarization in other processes; a study of $`\mathrm{}p\mathrm{\Lambda }^{}X`$ and $`e^+e^{}\mathrm{\Lambda }^{}X`$ is in progress .

Acknowledgements D.B. thanks RIKEN, Brookhaven National Laboratory and the U.S. Department of Energy for support under contract # DE-AC02-98CH10886; M.A. thanks RIKEN BNL Research Center for support during a Workshop last year, when this work was initiated.


text:
# Weak ferromagnetism and field-induced spin reorientation in K2V3O8

\[

## Abstract

Magnetization and neutron diffraction measurements indicate long-range antiferromagnetic ordering below $`T_N`$=4 K in the 2D, S=1/2 Heisenberg antiferromagnet K2V3O8. The ordered state exhibits “weak ferromagnetism” and novel, field-induced spin reorientations. These experimental observations are well described by a classical, two-spin Heisenberg model incorporating Dzyaloshinskii-Moriya interactions and an additional c-axis anisotropy. This additional anisotropy can be accounted for by inclusion of the symmetric anisotropy term recently described by Kaplan, Shekhtman, Entin-Wohlman, and Aharony. This suggests that K2V3O8 may be a very unique system where the $`qualitative`$ behavior relies on the presence of this symmetric anisotropy.

\]

One of the most exotic and intriguing anisotropies which can occur in magnetic solid state systems is the antisymmetric Dzyaloshinskii-Moriya (DM) interaction. These interactions, originally proposed by Dzyaloshinskii  and later derived by Moriya  through extension of Anderson’s superexchange model , arise from spin-orbit coupling in the presence of sufficiently low local lattice symmetry. The form of this interaction is typically written $`𝐃(𝐒_\mathrm{𝟏}\times 𝐒_\mathrm{𝟐}`$) where $`𝐃`$ is a vector determined by the symmetry of the lattice. Competition between the vector product in the DM interaction and the scaler product in the Heisenberg interaction often leads to novel phenomena such as spin-canting  and spiral spin structures .

Relatively recently, Kaplan , and later Shekhtman, Entin-Wohlman, and Aharony (KSEA)  (motivated by relevance to parent compounds of high-Tc superconductors), re-examined the DM interaction with particular emphasis on an often ignored symmetric component. This symmetric anisotropy always accompanies the typical antisymmetric DM term but is typically ignored as its strength is an order of magnitude weaker. KSEA showed that, in the general case where there is no $`frustration`$, this symmetric term cannot be ignored as its inclusion restores the original O(3) symmetry of the Heisenberg Hamiltonian. We report in this letter the novel properties of a S=1/2, 2d Heisenberg system, K2V3O8, and suggest that they can be explained by a model including both symmetric and antisymmetric DM interactions. This implies that K2V3O8 is a novel system where the $`qualitative`$ behavior relies heavily on inclusion of the symmetric component of the DM interaction.

K2V3O8 crystallizes in a tetragonal unit cell with space group P4bm and lattice constants a=8.870 Å and c=5.215 Å. The structure is shown in Fig. 1(a) and consists of magnetic V4+-O5 pyramids and non-magnetic V5+-O4 tetrahedra with interstitial K+ ions. Previous magnetic measurements consisted of powder magnetization from 5-300 K  which were best described by a 2D Heisenberg model with coupling constant J=12.6 K and a g-value of 1.89. EPR measurements on single crystals suggest very small anisotropy with g-values of gc=1.922 and gab=1.972. Liu et al.  suggested from observation of field-dependent magnetization that K2V3O8 may order at lower temperatures.

Single crystal plates of K2V3O8 (typical dimensions: 1$`\times `$1$`\times `$0.1 cm3) were grown in a platinum crucible by cooling appropriate amounts of VO2 in a molten KVO3 flux. Additional details of the crystal growth will be reported elsewhere. One of these single crystals was mounted in a Quantum Design SQUID magnetometer equipped with a sample rotator and the magnetization is shown in Fig. 2 as a function of both temperature (a) and applied magnetic field (b) for fields applied along both the c-axis and within the tetragonal basal plane. The M(T) data (Fig. 2(a)) were measured by cooling in the presence of an applied magnetic field while the M(H) data (Fig. 2(b)) were taken by cooling to base temperature in zero-field. From Fig. 2(a), the magnetization is seen to be very isotropic for temperatures in excess of about 8 K confirming the Heisenberg nature of the interactions. At high temperatures, we observe a Curie-Weiss susceptibility with an effective moment ($`\mu `$=g$`\mu _B`$\[S(S+1)\]1/2) of 1.7$`\mu _B`$ per formula unit (FU). This result is consistent with the presence of two non-magnetic V-ions per FU and one magnetic ion with S=1/2.

The temperature dependence in the presence of a weak magnetic field shows a clear ordering phase transition below a temperature of about 4 K. The rapid decrease of the magnetization below this temperature with field applied along the c-axis is indicative of antiferromagnetic ordering state with an easy axis along c. However, measurements within the basal plane show a clear ferromagnetic enhancement upon entering the ordered state. The ferromagnetic ordered moment is small ($``$10-2 $`\mu _B`$) and thus we conclude that the system is antiferromagnetic with spins primarily along the c-axis and a small canted moment in the basal plane. This represents the first observation of an ordering phase transition in this material.

The behavior of the magnetization as a function of field for temperatures well below the Néel temperature is shown in Fig. 2(b). For field applied along c, we observe an abrupt increase in magnetization for fields in excess of about 0.85 T consistent with a discontinuous spin-flop phase transition, further confirming the easy nature of the c-axis. For fields within the basal plane, we observe an intriguing field-induced phase transition which occurs at a magnetic field of HSR=0.65 T. There was no evidence of irreversibility in M(H) for either field orientation or in angular rotations of the sample in the presence of an applied field and no evidence of thermal hysteresis was observed.

To investigate the detailed microscopic spin arrangement, neutron diffraction measurements were performed on the HB1A and HB1 triple axis spectrometers at the High Flux Isotope Reactor in Oak Ridge National Laboratory. For the measurements in zero-field, the sample was mounted in the ($`h`$0$`l`$) scattering plane and comparison of scattering at 1.6 K and 10 K indicated the presence of magnetic Bragg peaks (in the ($`h`$0$`l`$) plane) described by integer indices with $`h`$ odd. These peaks suggest a magnetic unit cell in which the corner and face-center moments within the basal plane are aligned antiparallel with parallel alignment along the c-axis. Examination of the integrated intensities of 6 magnetic reflections indicated that the moments lie along the c-axis consistent with conclusions drawn from the magnetization measurements. The value of the ordered moment ($`\mu `$=g$`\mu _B`$S) was found to be 0.72(4) $`\mu _B`$ at T=1.6 K, reduced from the expected value of 1 $`\mu _B`$. This can likely be attributed to quantum fluctuations in the low-dimensional S=1/2 system. The zero-field magnetic structure is shown schematically by the red arrows in Fig. 1(b).

To examine the behavior in the presence of a basal plane magnetic field, the crystal was mounted in the ($`h`$0$`l`$) scattering plane and a vertical field applied along the (010) direction. The field dependence of the (100) and (101) magnetic Bragg peaks are shown in Figs. 3(a) and 3(b). As the intensity of magnetic neutron scattering is proportional to the component of magnetic moment normal to the wavevector transfer, Q, the observed disappearance of the (100) magnetic Bragg peak and the concomitant enhancement of the (101) reflection are consistent with a continuous rotation of the spins from the (001) to the (100) direction with increasing magnetic field strength. The rotation is complete at a field strength of HSR=0.65 T in excellent agreement with the anomaly in the magnetization. This is shown schematically by the blue arrows in Fig. 1(b). Thus, we conclude that this transition is a peculiar spin-reorientation phase transition whereby the spins rotate from the c-axis to the basal plane while remaining normal to the applied field direction.

A second crystal was mounted in the (hk0) scattering plane to allow measurements with field applied along the c-axis. The results of these measurements are shown in Fig 3(c) and (d). As a function of applied magnetic field, the (100) Bragg peak intensity drops abruptly upon passing through a field of 0.85 T suggesting a discontinuous reorientation of the spins. The temperature dependence in an applied field of 1 T shows a strikingly different behavior for the (100) and (010) Bragg peaks and the weaker intensity for (010) suggests that the moments have “flopped” from the c-axis into the basal plane selecting a specific direction within the basal plane in close proximity to the (010) direction. No evidence for a structural distortion has been observed and the selection of a specific direction likely results from a slight misalignment of the magnetic field with the c-axis. This misalignment was less than 1 for the measurements performed but even this level may induce a net field in the basal plane of sufficient strength to overcome domain energies thus selecting a specific magnetic domain. Such a domain selection has been observed in the tetragonal antiferromagnet MnF2 .

As mentioned in the introduction, weak ferromagnets often have at their origin the competition between the antisymmetric DM interaction and the symmetric Heisenberg interaction. To consider the DM interactions in K2V3O8, we will consider the projection of the crystal structure onto the tetragonal basal plane as shown in Fig. 4(a). The plane of inversion symmetry between near-neighbor V4+ spins results in DM vectors located in that inversion plane  giving components along (110) and (001) directions, denoted by $`𝐃_{\mathrm{𝐱𝐲}}`$ and $`𝐃_𝐳`$ respectively (these directional vectors are indicted in the figure relative to the central site). The four-fold rotation symmetry through the central site gives rise to the configuration of $`𝐃`$ vectors shown in the figure. From the neutron diffraction measurements, we know that the magnetic unit cell contains two spins in both zero and non-zero fields and, hence, we will consider a two-spin, classical Hamiltonian. Noting that in averaging over the the near-neighbor spins within a unit cell, the components of the DM vectors in the basal plane (i.e. the $`𝐃_{\mathrm{𝐱𝐲}}`$ terms) cancel, we write the ground state energy as

$$F=J𝐌_\mathrm{𝟏}𝐌_\mathrm{𝟐}+D_z[𝐌_\mathrm{𝟏}\times 𝐌_\mathrm{𝟐}]_z$$
(1)
where $`𝐌_\mathrm{𝟏}`$ and $`𝐌_\mathrm{𝟐}`$ are the sublattice magnetizations and $`[𝐌_\mathrm{𝟏}\times 𝐌_\mathrm{𝟐}]_z`$ is the z-component of the vector cross-product.

This ground state energy will be minimized when the spins are in the basal plane and, consequently, an additional c-axis anisotropy is required to account for the easy c-axis observed experimentally. Thus, we add a c-axis anisotropy of the form $`E_zM_{1z}M_{2z}`$ and, with addition of the Zeeman energy for nonzero field, we arrive at the ground state energy,

$`F`$ $`=`$ $`J𝐌_\mathrm{𝟏}𝐌_\mathrm{𝟐}+D_z[𝐌_\mathrm{𝟏}\times 𝐌_\mathrm{𝟐}]_z+E_zM_{1z}M_{2z}`$ (2)
$`+`$ $`g\mu _B𝐇[𝐌_\mathrm{𝟏}+𝐌_\mathrm{𝟐}].`$ (3)
This rather simple, classical ground state energy describes all the low temperature properties of the system. For appropriate values of $`D_z`$ and $`E_z`$, the spins will point along the c-axis with canting observed in small fields due to the normal competition between DM and Heisenberg interactions. For fields applied in the basal plane, the Zeeman energy and the DM interaction favor arrangement of the spins in the basal plane. This competes with the c-axis anisotropy causing a continuous rotation of the spins to the basal plane with increasing field strength. The spins remain primarily normal to the applied field direction to minimize the Zeeman energy. Finally, with field applied along the c-axis, a spin-flop transition occurs due to competition between Zeeman energy and the anisotropies in the system causing an abrupt rotation of the spins into the basal plane.

As we have two field-induced phase transitions, we can estimate the values for Dz/J and Ez/J required to produce phase transitions at the measured field values (using J=12.6K and g=1.89 as determined from Liu et al. ). This results in a value of Dz/J=0.22 and Ez/J=3.8$`\times `$10-2 (this value of D/J seems large considering the measured g-values  but compares favorably to the DM spiral system BaCuGe2O7 which has a similar basal plane structure and for which D/J was found to be $``$0.18 ). Using these values, we can calculate the induced magnetic moment along the field direction for field applied along both the a- and c-axes and the results are plotted in Fig. 4(d) as a function of applied magnetic field in Tesla. These results can be directly compared to the M(H) data plotted in Fig. 2(b) and the qualitative agreement is remarkable suggesting that this relatively simple model adequately describes the low temperature properties of the system. To further emphasize this, we have superimposed on Fig. 3(a) and 3(b) the calculated magnetic Bragg peak intensity for the (100) and (101) reflections in the presence of an applied magnetic field along the (010) direction. These results are represented by the solid lines in Figs. 3(a) and 3(b) and the agreement between the calculations and the data is excellent, particularly for the (100) reflection. The agreement for the (101) reflection is very good for field strengths up to HSR but seems to deviate for higher fields. The reason for this disagreement is unclear but could result from a field-dependent $`<`$S$`>`$. Such a scenario can occur via field-dependent suppression of quantum fluctuations, as has been observed in some ABX3 antiferromagnets .

To shed some light on the nature of the c-axis anisotropy, we calculate the symmetric anisotropy term of KSEA. The addition of this term causes a rescaling of J to (J-D2/4J) and the inclusion of a symmetric anisotropy term (1/2J)$`𝐌_\mathrm{𝟏}`$$``$$`𝐀`$$``$$`𝐌_\mathrm{𝟐}`$ where $`𝐀`$ is a 3$`\times `$3 matrix with elements Auv=DuDv . Making these changes to the Hamiltonian and averaging over near-neighbor spins produces the ground state energy,

$`F`$ $`=`$ $`(J{\displaystyle \frac{D_z^2}{4J}})𝐌_\mathrm{𝟏}𝐌_\mathrm{𝟐}+D_z[𝐌_\mathrm{𝟏}\times 𝐌_\mathrm{𝟐}]_z`$ (4)
$`+`$ $`({\displaystyle \frac{D_z^2}{2J}}{\displaystyle \frac{D_{xy}^2}{4J}})M_{1z}M_{2z}+g\mu _B𝐇[𝐌_\mathrm{𝟏}+𝐌_\mathrm{𝟐}].`$ (5)
This energy has precisely the same form as Eq. 2 and, consequently, the c-axis anisotropy falls out naturally upon inclusion of both the symmetric and antisymmetric anisotropies. The maximum value for the c-axis anisotropy in Eq. 3 is D$`{}_{}{}^{2}{}_{z}{}^{}`$/2J which, for Dz/J of 0.22, gives a maximum anisotropy (Ez/J) of 2.4$`\times `$10-2. Within the level of approximation involved in this model, this is in reasonably good agreement with the above value of 3.8$`\times `$10-2. Consequently, we conclude that K2V3O8 represents a unique system where the $`qualitative`$ features rely heavily on the inclusion of the KSEA symmetric anisotropy term. The only known, clear evidence of this symmetric anisotropy occurs in the DM spiral system Ba2CuGe2O7 where the addition of this term is necessary to produce $`quantitative`$ agreement between experiment and theory . As a caveat, it is important to note that we cannot rule out other possible mechanisms for this additional anisotropy. Rather, we simply note that inclusion of this symmetric anisotropy alone, which is necessary to properly consider DM interactions, seems to adequately describe the properties of the system.

In summary, we have observed an ordering phase transition in K2V3O8 with a primarily antiferromagnetic ordered state accompanied by weak ferromagnetism. Field-induced spin reorientation phase transitions have been observed with field applied along both the c-axis and in the tetragonal basal plane. The spin arrangement has been determined by single-crystal neutron diffraction measurements. Remarkably, this rich magnetic behavior can be well described by a simple, classical, two-spin Heisenberg model with inclusion of DM interactions and an additional c-axis anisotropy. This additional anisotropy can, at least partially, be accounted for by inclusion of the KSEA symmetric anisotropy term making K2V3O8 a unique system where introduction of this interaction is necessary to describe the $`qualitative`$ behavior of the system.

We would like to acknowledge stimulating discussions with G. Murthy, M.W. Meisel, R. Coldea, and D.A. Tennant. Oak Ridge National Laboratory is managed by UT-Battelle, LLC, for the U.S. Dept. of Energy under contract DE-AC05-00OR22725.