format: text

text:
# 1 Introduction

## 1 Introduction

The classical semi-implicit (SI) technique (Robert et al., 1972) has been widely used in NWP since it provides efficient and simple algorithms, at least for spectral models. This classical SI method requires the definition of a constant in time linear "reference" operator $`^{}`$, which usually consists in the linearisation of the original system $``$, around a stationary reference-state, noted $`𝒳^{}`$. For a given state $`𝒳`$ of the atmosphere, the evolution of the system, $`(𝒳/t)=.𝒳`$, is then time-discretised through:

$$\frac{\delta 𝒳}{\delta t}=(^{}).𝒳+^{}.\overline{[𝒳]}^t$$
(1)
where $`(\delta /\delta t)`$ is the discretised time-derivative operator, and $`\overline{[\text{ }]}^t`$ is the implicit-centred temporal average operator. The terms linked to the reference operator $`^{}`$ are thus treated in a centred-implicit way, whilst the residual non-linear terms are treated explicitly. For this scheme, there is no formal proof of the stability in real atmospheric conditions, due to the explicit treatment of non-linear residuals. This prompted the authors of pioneering NWP applications of the SI scheme to examine theoretically its stability in idealised contexts.

In a seminal study following this approach, Simmons et al., 1978 (SHB78 hereafter) analysed the stability of the SI scheme for the hydrostatic primitive equations (HPE) system with a Leap-Frog (3-TL hereafter) time-discretisation by considering the linearised equations around a stationary state $`\overline{𝒳}`$ (referred to as "atmospheric state" hereafter) when the resulting linear "atmospheric" operator $`\overline{}`$ deviates from the linear "reference" operator $`^{}`$ of the SI scheme, thus generating potentially unstable explicitly-treated residuals.

In the vertically-continuous context they performed a stability analysis valid when the eigenfunctions of $`\overline{}`$ and $`^{}`$ are identical. They showed that when the atmospheric and reference temperature profiles (respectively $`\overline{T}`$ and $`T^{}`$) are isothermal, the stability of the SI scheme requires:

$$0\overline{T}2T^{},$$
(2)
hence $`T^{}`$ cannot be chosen arbitrarily for applying the 3-TL SI scheme to the HPE system.

In the finite-difference vertically-discretised context, thew showed that a "vertically-discretised analysis" of stability following the same principle simply resulted in the solution of a standard eigenvalue problem. They found empirically that a large static-stability for the reference-state is necessary to maintain the stability of the scheme for realistic thermal atmospheric profiles. As a consequence, they recommended to use a warm isothermal state as reference-state, a rule which was then widely adopted for NWP applications using SI schemes.

Finally, they examined the effect of applying a second-order time-filter in the temporal average of linear terms, and found that an improved stability is obtained, but at the expense of an increased misrepresentation of the wave propagation.

Côté et al., 1983 (CBS83 hereafter), still in the HPE context, examined the stability of the 3-TL SI scheme for a finite-element vertical discretisation using the above vertically-discretised analysis method. They established a stability criterion for the 3-TL SI scheme in terms of the atmospheric and reference static stabilities ($`\overline{\gamma }`$ and $`\gamma ^{}`$ respectively):

$$0\overline{\gamma }2\gamma ^{},$$
(3)
therefore generalizing (2) to not necessarily isothermal thermal profiles.

Still in the HPE context, Simmons and Temperton, 1997 showed with the same method that extrapolating two-time level (2-TL) schemes have more stringent stability constraints than their 3-TL counterpart. For instance, in the isothermal framework of the SHB78 analysis, the stability of the 2-TL SI scheme requires:

$$0\overline{T}T^{}.$$
(4)
As a consequence, they recommended to use a warmer reference temperature than in the 3-TL case. ST97 also showed that the 2-TL SI scheme was intrinsically damping when stable, a characteristic which was not present in the 3-TL SI scheme. However, 3-TL schemes require a time-filter to damp the unstable computational modes, and this time-filter also damps the transient physical modes. They argue that all these aspects being considered, the effective damping of 2-TL and 3-TL schemes is of comparable overall intensity. This particular debate will be ignored in this paper, in order to focus on stability aspects only.

The relevance of the SI method for solving numerically the fully compressible Euler Equations (EE) was then advocated (Tanguay et al., 1990), and some numerical models in EE using this technique were effectively developed: Caya and Laprise (1999), presented a model in EE with a 3-TL SI scheme (with a moderate time-decentering first-order accurate in time). Semazzi et al. (1995) and Qian et al. (1998) also showed a model in EE but with a 2-TL SI scheme (with a strong time-filter however). Besides, the need of more robust schemes than the classical SI one for solving the EE system became recognized (e.g. Côté et al. 1998, BHBG95), probably motivated by some pathological behaviours with the classical SI scheme under some circumstances.

Schemes with evolution terms treated in a more centred-implicit way are usually believed to have an increased robustness, hence fulfilling the latter emerging need, and some of them were developed for fine-scale models in the EE system. Bubnová et al., 1995 used a 3-TL scheme in which the leading non-linear terms of the EE system are treated in a centred-implicit way, through a partially iterative method. For a 2-TL scheme, Côté et al, 1998, used a fully iterative method, aiming to treat all the evolution terms of the EE system in a centred-implicit way. Cullen (2000) examined the benefit of using such a fully iterative scheme for the HPE system, arguing that an improved accuracy could be obtained beside the improved stability (this latter being not strictly required however for current HPE applications). As a formal justification, he examined the stability of this iterative scheme for the 2-TL shallow-water (SW) system. The analysis was limited to a scheme called "predictor/corrector", which consists in a single additional iteration after the SI scheme. The salient result was that the additional iteration in the "predictor/corrector" scheme allows to recover a extended range of stability as in (2) instead of (4). In the following, these fully iterative schemes with a more centred-implicit treatment of the evolution terms will be referred to as "iterative centred-implicit" (ICI) schemes.

From the theoretical point of view, the current situation is that no stability analysis has been provided for the EE system with SI scheme, and for ICI schemes with more iterations, stability analyses are available only for the SW system.

Here we present a general method to carry out space-continuous stability analyses of the various time-discretisation schemes mentioned above, for any usual meteorological system of NWP interest (SW, HPE, EE), on canonical problems similar to those examined in SHB78. Some original results concerning the EE system and iterative schemes are presented.

This work may also be viewed as a first theoretical investigation about the suitability of various time-discretisation schemes for solving numerically the EE system.

## 2 General framework for analyses

The general framework for the stability analyses presented here is basically the same as in most earlier studies: The flow is assumed adiabatic inviscid and frictionless in a non-rotating dry perfect-gas atmosphere with a Cartesian coordinate system. Moreover, the flow is assumed linear around an "atmospheric" basic-state $`\overline{𝒳}`$. The actual evolution of the atmospheric flow is thus described by $`\overline{}`$, the linear-tangent operator to $``$ around $`\overline{𝒳}`$. The atmospheric basic-state $`\overline{𝒳}`$ is chosen stationary, resting, horizontally homogeneous, and hydrostatically balanced. The governing equation for the flow is then:

$$\frac{𝒳^{}}{t}=\overline{}.𝒳^{}$$
(5)
where $`𝒳^{}=𝒳\overline{𝒳}`$, and the primes are dropped henceforth for clarity.

Following the usual practice in NWP, the linear operator $`^{}`$ in (1) is taken as the tangent-linear operator to $``$ around a reference-state $`𝒳^{}`$ which is also chosen stationary, resting, horizontally homogeneous, and hydrostatically balanced. Since $`\overline{𝒳}`$ is a resting state, the linear Lagrangian time-derivative coincides with the Eulerian time-derivative, and the LHS operator of (5) hence holds for Eulerian models as well as for semi-Lagrangian models .

## 3 The class of ICI schemes in the linear framework

In the restricted resting and linear framework of section 2, classical SI schemes as well as the iterative schemes mentioned in the introduction can be gathered in a single class of ICI schemes, differing only by their number of iterations. These ICI scheme are first presented for a 2-TL discretisation. In this case, the fully implicit-centred (FIC) scheme writes:

$$\frac{𝒳^+𝒳^0}{\mathrm{\Delta }t}=\frac{\overline{}.𝒳^++\overline{}.𝒳^0}{2}$$
(6)
where, according to a standard practice for time-discretised equations, the superscripts "+" and "0" stand for time levels $`(t+\mathrm{\Delta }t)`$ and $`t`$ respectively. The principle of ICI schemes is to approach the FIC solution by starting from an initial "guess" noted $`𝒳^{+(0)}`$, then iterating the following algorithm:

$`{\displaystyle \frac{𝒳^{+(n)}𝒳^0}{\mathrm{\Delta }t}}`$ $`=`$ $`{\displaystyle \frac{\overline{}.𝒳^{+(n1)}+\overline{}.𝒳^0}{2}}+{\displaystyle \frac{^{}.𝒳^{+(n)}^{}.𝒳^{+(n1)}}{2}}`$ (7)
$``$ $`{\displaystyle \frac{\left(\overline{}^{}\right).𝒳^{+(n1)}+\left(\overline{}^{}\right).𝒳^0}{2}}+{\displaystyle \frac{^{}.𝒳^{+(n)}+^{}.𝒳^0}{2}}`$ (8)
for $`n=1,2,\mathrm{},N_{\mathrm{iter}}`$. The $`𝒳^+`$ state, valid a $`(t+\mathrm{\Delta }t)`$ is then taken as the last iterated state $`𝒳^{+(N_{\mathrm{iter}})}`$. An examination of this scheme for a model with a single prognostic variable without spatial dependency shows that it acts as a fixed-point algorithm for solving the implicit non-linear scalar equation $`f(x)=x`$, by using an estimate $`f_{}^{}{}_{}{}^{}`$ of the derivative $`f^{}`$ as a preconditioner for convergence. The method converges if $`\left|(f^{}f_{}^{}{}_{}{}^{})/f_{}^{}{}_{}{}^{}\right|<1`$. This is a weaker condition than the one for the classical (i.e. not preconditioned) fixed-point method: $`\left|f^{}\right|<1`$. The initial guess $`𝒳^{+(0)}`$ is arbitrary in ICI schemes, but choosing an appropriate initial guess may help decreasing the magnitude of the discrepancy between FIC and ICI schemes after a fixed number of iterations.

For 3-TL schemes, the ICI scheme can be defined by:

$`{\displaystyle \frac{𝒳^{+(n)}𝒳^{}}{2\mathrm{\Delta }t}}`$ $`=`$ $`{\displaystyle \frac{\overline{}.𝒳^{+(n1)}+\overline{}.𝒳^{}}{2}}+{\displaystyle \frac{^{}.𝒳^{+(n)}^{}.𝒳^{+(n1)}}{2}}`$ (9)
$``$ $`{\displaystyle \frac{\left(\overline{}^{}\right).𝒳^{+(n1)}+\left(\overline{}^{}\right).𝒳^{}}{2}}+{\displaystyle \frac{^{}.𝒳^{+(n)}+^{}.𝒳^{}}{2}}`$ (10)
where the superscript "-" denotes a variable taken at the time-level $`(t\mathrm{\Delta }t)`$.

Here follows a list of some schemes proposed in the literature, with the corresponding characteristics ($`𝒳^{+(0)}`$ ,$`N_{\mathrm{iter}}`$), in the restricted framework of section 2:

* Classical 2-TL SI extrapolating scheme: $`N_{\mathrm{iter}}=1`$ and $`𝒳^{+(0)}=(2𝒳^0𝒳^{})`$.
* 2-TL non-extrapolating SI scheme: $`N_{\mathrm{iter}}=1`$ and $`𝒳^{+(0)}=𝒳^0`$. However, this scheme is not used in practice since it is only first-order accurate in time, as mentioned in Cullen (2000).
* "Predictor/corrector" scheme of Cullen (2000): $`N_{\mathrm{iter}}=2`$ and $`𝒳^{+(0)}=𝒳^0`$.
* Iterative scheme of Côté et al., 1998: general iterative ICI scheme, but used with $`N_{\mathrm{iter}}=2`$ in practice. The choice of $`𝒳^{+(0)}`$ is not explicitly indicated.
* FIC scheme: $`N_{\mathrm{iter}}=\mathrm{}`$ (does not depend on the choice of $`𝒳^{+(0)}`$). This scheme can not be achieved in practice for numerical models, but it may be useful for theoretical examination of the asymptotic behaviour of the ICI schemes.

For 3-TL discretisations, the SI scheme, which corresponds to an ICI scheme with $`N_{\mathrm{iter}}=1`$ and $`𝒳^{+(0)}=(2𝒳^0𝒳^{})`$ is the only one to be used in practice. However, iterated 3-TL ICI schemes could be used as well, and the 3-TL FIC scheme is equivalent to a 2-TL FIC scheme with a time-step twice as long.

In the general framework, when $``$ is not linear, these various schemes cannot be gathered in the unique formalism (7) or (9).

Addition of a second-order time-filter in the above definitions for 2-TL schemes is straightforward. Two main variants are usually considered, depending on whether the filtering is applied only to the time-averages of linear terms (as in SHB78), or also to the time-averages of non-linear terms, in (8). For instance, in a 2-TL SI scheme, the first variant of the filter consists in replacing $`^{}(𝒳^{+(n)}+𝒳^0)`$ by $`^{}[(1+\kappa )𝒳^{+(n)}+(12\kappa )𝒳^0+\kappa 𝒳^{}]`$ in (8), where $`\kappa `$ is a (small) positive parameter, and for the second variant, the same modification is also applied to the first RHS term of (8). The scheme then becomes essentially a 3-TL scheme since information at level $`𝒳^{}`$ is always used. However, the use of large values of $`\kappa `$ (e.g. $`\kappa =0.5`$, which eliminates the $`𝒳^0`$ contribution) is known to deteriorate the solution through a spurious damping of transient perturbations (e.g. Hereil and Laprise, 1996).

For 3-TL schemes, second-order time-filters are uneffective, and a first-order accurate time-decentering must be used. This consists in replacing $`^{}(𝒳^{+(n)}+𝒳^{})`$ by $`^{}[(1+ϵ)𝒳^{+(n)}+(1ϵ)𝒳^{}]`$ in (10), where $`ϵ`$ is a (small) positive parameter. The scheme then ceases to be second-order accurate in time since the time-average is no longer centred in time. This results in a spurious damping of transient perturbations even for moderate values of $`ϵ`$ (due to the weaker time-selectivity of the filter $`ϵ`$ compared to $`\kappa `$).

## 4 Conditions for space-continuous analyses

### 4.1 Conditions on the upper and lower boundaries

Space-continuous analyses are much easier to carry out when the equation system is defined in the whole unbounded atmosphere, because the expression of the normal modes of the system is more general. The following space-continuous analyses will restrict to this case (although this is not strictly required). For systems in which the vertical direction is represented (e.g. HPE and EE systems), this means that the upper and lower boundary conditions must not appear explicitly in the set of governing equations. However for systems cast in mass-based coordinates \[such as HPE system in pressure-based coordinate (e.g. SHB78), and EE system in hydrostatic pressure-based coordinate (Laprise, 1992, L92 hereafter)\], the upper and lower boundary conditions actually appear inside the set of equations through vertical integral operators with definite bounds at the boundaries of the vertical domain. When they are present, it is assumed that these integral operators can be eliminated, i.e. that $`\overline{}`$, $`^{}`$ can be transformed to "unbounded" operators by application of appropriate vertical linear differential operators to the prognostic equations which originally involve integral operators, in order that e.g. (5) rewrites as:

$$\frac{}{t}\left(\begin{array}{c}l_1𝒳_1\\ \mathrm{}\\ l_P𝒳_P\end{array}\right)=\left(\begin{array}{ccc}l_1\overline{}_{11}& \mathrm{}& l_P\overline{}_{1P}\\ \mathrm{}& \mathrm{}& \mathrm{}\\ l_1\overline{}_{P1}& \mathrm{}& l_P\overline{}_{PP}\end{array}\right).\left(\begin{array}{c}𝒳_1\\ \mathrm{}\\ 𝒳_P\end{array}\right)$$
(11)
where $`P`$ is the number of prognostic variables of the unbounded system, $`(l_1,\mathrm{},l_P)`$ are linear vertical operators, and $`(l_1\overline{}_{11},\mathrm{},l_P\overline{}_{PP})`$ are linear spatial operators which no longer contain any reference to the upper and lower boundaries. The transformed system obtained for $`\overline{}`$ can then be written as:

$$\frac{l.𝒳}{t}=l.\overline{}.𝒳$$
(12)
where $`l`$ is the diagonal matrix $`(l_1,\mathrm{},l_P)`$. A similar condition must be true for $`^{}`$ as well: it is assumed that applying the same operator $`l`$ to $`^{}`$ leads to an operator $`l^{}`$ for which $`l_i_{ij}^{}`$ does not contain any reference to the upper and lower boundaries for $`(i,j)(1,\mathrm{},P)`$. The first condition for the following analyses is:

* \[C1\]: There exists a linear operator $`l`$ such as $`l\overline{}`$ and $`l^{}`$ have no reference to the upper and lower boundaries.

The system (12) is henceforth referred to as the "unbounded" system.

### 4.2 Conditions on the stability of the $`\overline{𝒳}`$ state

The aim of these analyses is to determine in which conditions a stationary state $`\overline{𝒳}`$ for $`l\overline{}`$ will remain a stable equilibrium-state in the time-discretised context, provided it is a stable equilibrium-state in the time-continuous context. Hence the analyses will be restricted to stationary states $`\overline{𝒳}`$ which are in stable equilibrium. Given the linear context used here, a physical transposition of this condition is:

* \[C2\]: For any perturbation $`𝒳(t=0)`$ around $`\overline{𝒳}`$ with a bounded energy-density, the time-evolution $`𝒳(t)`$ resulting from (12) must have a bounded energy-density.

The condition is formulated with energy-density instead of total energy because the domain is unbounded in space. The complex eigenmodes of the unbounded system (12) are the complex functions of space $`𝒳(𝐫)`$ which satisfy:

$$l\overline{}𝒳(𝐫)=\overline{\lambda }l𝒳(𝐫)$$
(13)
where $`\overline{\lambda }\text{ }\mathrm{C}`$ and $`𝐫`$ denotes the spatial dependency. The time-evolution of the mode is bounded in time if $`\overline{\lambda }`$ is a pure imaginary number, and $`𝒳(𝐫)`$ is then a "normal mode" in the usual terminology. Theoretical arguments out of the scope of this paper show that a mathematical transposition of \[C2\] writes:

* \[C2’\]: For any complex eigenmode of $`l\overline{}`$, $`\overline{\lambda }i\mathrm{I}\mathrm{R}𝒳(𝐫)`$ has a bounded energy-density.

### 4.3 Conditions on the normal modes of the linear unbounded system

The time-continuous normal modes of the original system (5) are the complex functions of space $`𝒳(𝐫)`$ which have an oscillatory temporal evolution. Hence they satisfy $`\overline{}𝒳(𝐫)=i\overline{\omega }𝒳(𝐫)`$. Similarly, the time-continuous normal modes of the linear unbounded system (12) are the complex functions of space $`𝒳(𝐫)`$ for which $`l𝒳(𝐫)`$ has an oscillatory temporal evolution. Hence they satisfy $`l.\overline{}𝒳(𝐫)=i\overline{\omega }l.𝒳(𝐫)`$.

Any normal mode of the original system is also necessarily a normal mode of the unbounded system, with the same frequency $`\overline{\omega }`$. For any normal mode of the unbounded system , we can choose the origin of space o in such a way that this mode writes:

$$𝒳(𝐫)=𝒳(𝐨)f(𝐫)=\widehat{𝒳}f(𝐫)$$
(14)
with $`\widehat{𝒳}=(\widehat{𝒳}_1,\mathrm{},\widehat{𝒳}_P)\text{ }\mathrm{C}^P`$, and $`f=(f_1,\mathrm{},f_P)`$ is a vector of space-dependent functions (the product $`\widehat{𝒳}f`$ is understood "component by component"). For the indices $`i`$ such as $`𝒳_i`$ uniformly vanishes, $`\widehat{𝒳}_i=0`$ and $`f_i(𝐫)=1`$ are assumed by convention. The function $`f`$ and the vector $`\widehat{𝒳}`$ are respectively termed the "structure" and the "polarisation vector" of the mode. The two following conditions are required for the proposed analyses:

* \[C3\]: For any normal mode $`𝒳`$ of the unbounded linear atmospheric system with a structure $`f(𝐫)`$, $`l_i.f_i(𝐫)`$ must be proportional to $`f_i(𝐫)`$:

$$i(1,\mathrm{},P),l_if_i(𝐫)=\xi _if_i(𝐫)\mathrm{with}\xi _i\text{ }\mathrm{C}^{}.$$
(15)
and:

* \[C4\]: For any normal mode $`𝒳`$ of the unbounded linear atmospheric system with a structure $`f(𝐫)`$, $`\overline{}_{ij}.f_j(𝐫)`$ \[resp. $`_{ij}^{}.f_j(𝐫)`$\] must be proportional to $`f_i(𝐫)`$:

$$(i,j),l_i\overline{}_{ij}.f_j(𝐫)=\overline{\mu }_{ij}f_i(𝐫)\mathrm{and}l_i_{ij}^{}.f_j(𝐫)=\mu _{ij}^{}f_i(𝐫),\mathrm{with}(\overline{\mu }_{ij},\mu _{ij}^{})\text{ }\mathrm{C}.$$
(16)
As will be seen below, these latter two conditions have the important consequence that for any normal mode of the unbounded system, each individual time-discretised prognostic equation for $`𝒳_i(𝐫)`$ becomes a scalar equation. This key ingredient makes the analysis straightforward for every member of the ICI class.

### 4.4 Comments

Since the set of normal modes for the unbounded system encompasses the set of normal mode of the bounded system, the transformation from the bounded system to the unbounded system is not likely to "mask" some instabilities of the original system, unless the causes of the instability lie in the boundary conditions themselves. However, discretised-analyses allow to clarify this point by showing that the stability of the bounded and unbounded systems are actually found to be similar in practice. Besides, discretised-analyses of the unbounded system are by nature impossible to perform, hence the continuous analysis is the only way to estimate the intrinsic stability of the unbounded system and to demonstrate that instabilities, when they occur in a practical application, are not due to a weakness in the spatial discretisation or even to the boundary conditions, but actually to the time-discretised propagation of free modes inside the atmosphere.

In spite of their apparently abstract and constraining form, conditions \[C1\]–\[C4\] are easy to verify with routine normal mode analysis techniques when examining a particular concrete meteorological system. The condition \[C2’\] restricts the set of stationary states $`\overline{𝒳}`$ around which the analysis is meaningful, and conditions \[C1\], \[C3\], \[C4\] restrict the spectrum of meteorological contexts accessible to the analysis, since they require a qualitative similarity between $`l`$, $`\overline{}`$ and $`^{}`$ operators. As stated in SHB78, analyses performed under this type of conditions "…grossly exaggerate the stability of the scheme…" since in more realistic meteorological contexts, the atmospheric and reference operators can be qualitatively much more different than imposed by this condition.

## 5 Time-Discretised Space-Continuous Analysis

The analysis examine the stability of the time-discretised system for perturbations which have a time-continuous normal mode structure. Hence we consider a given function $`f=(f_1,\mathrm{},f_P)`$ which is a normal mode structure for the time-continuous system, and we determine the normal modes of the time-discretised system which have the same structure $`f`$, by solving the equation:

$$\widehat{𝒳}_{(t=\mathrm{\Delta }t)}f(𝐫)=\lambda \widehat{𝒳}_{(t=0)}f(𝐫)$$
(17)
where $`\widehat{𝒳}_{(t=0)}`$ and $`\lambda `$ are the unknowns, and $`\widehat{𝒳}_{(t=\mathrm{\Delta }t)}`$ is determined using the time-discretisation scheme (7) or (9). For schemes using three time levels (as Leap-Frog or extrapolating 2-TL) a similar relationship $`\widehat{𝒳}_{(t=\mathrm{\Delta }t)}f(𝐫)=\lambda ^1\widehat{𝒳}_{(t=0)}f(𝐫)`$ must be added. If for some solution, $`\left|\lambda \right|>1`$ (resp. $`<1`$) the scheme is unstable (resp. damping) for this particular mode. The ratio $`\mathrm{Arg}(\lambda )/(i\overline{\omega }\mathrm{\Delta }t)`$ gives the relative phase-speed error of the scheme for this mode.

The analysis is described here for a 2-TL discretisation (7), but the transformation to a 3-TL scheme as well as the addition of time-filters such as $`\kappa `$ or $`ϵ`$ are straightforward. In the following of this section, the notation $`𝒳(𝐫,t)`$ is replaced by the usual superscript notation for time-discretised variables $`𝒳^t(𝐫)`$ as in section 3.

As a consequence of the discussion in section 3 and applying (17) , $`𝒳^{+(0)}`$ can be written as $`𝒳^{+(0)}(𝐫)=\mu (\lambda )𝒳^0(𝐫)`$, where $`\mu (\lambda )`$ depends on the choice of the initial guess $`𝒳^{+(0)}`$ (e.g. $`\mu =1`$ for a 2-TL non-extrapolating scheme, $`\mu =21/\lambda `$ for a 2-TL extrapolating scheme, etc…). The original unbounded system (12) is thus time-discretised following (7):

$`l𝒳^{+(0)}(𝐫)`$ $`=`$ $`\mu (\lambda )l𝒳^0(𝐫)`$
$`{\displaystyle \frac{l𝒳^{+(n)}(𝐫)l𝒳^0(𝐫)}{\mathrm{\Delta }t}}`$ $`=`$ $`{\displaystyle \frac{l\overline{}.𝒳^{+(n1)}(𝐫)+l\overline{}.𝒳^0(𝐫)}{2}}`$
$`+{\displaystyle \frac{l^{}.𝒳^{+(n)}(𝐫)l^{}.𝒳^{+(n1)}(𝐫)}{2}},n(1,\mathrm{},N_{\mathrm{iter}})`$
$`l𝒳^+(𝐫)l\lambda 𝒳^0`$ $`=`$ $`l𝒳^{+(N_{\mathrm{iter}})}(𝐫)`$
For schemes with $`N_{\mathrm{iter}}2`$, it is assumed that the intermediate states $`l𝒳^{+(n)}(𝐫)`$ for $`i(1,\mathrm{},N_{\mathrm{iter}}1)`$ have the same structure $`f`$ as the one currently examined, which allows to define the polarisation vector $`\widehat{𝒳}_i^{+(n)}`$ by:

$$l_i𝒳_i^{+(n)}(𝐫)=\widehat{𝒳}_i^{+(n)}l_if_i(𝐫).$$
(18)
Hence, using (15) and (16), the space dependency $`f_i(𝐫)`$ eliminates. For the generalized state-vector $`𝒵=(\widehat{𝒳}^0,\widehat{𝒳}^{+(0)},\widehat{𝒳}^{+(1)},\mathrm{},\widehat{𝒳}^{+(N_{\mathrm{iter}})})`$ of length $`(N_{\mathrm{iter}}+2)\times P`$, the above system writes:

$$\left(\begin{array}{cccccc}\mu (\lambda )I_P& I_P& 0_P& \mathrm{}& \mathrm{}& 0_P\\ M_1& M_2& M_3& \mathrm{}& & \mathrm{}\\ \mathrm{}& 0_P& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& 0_P\\ M_1& \mathrm{}& & \mathrm{}& M_2& M_3\\ \lambda I_P& 0_P& \mathrm{}& \mathrm{}& 0_P& I_P\end{array}\right).𝒵=𝐌.𝒵=0.$$
(19)
where $`I_P`$ (resp. $`0_P`$) is the unit (resp. null) $`P`$-order matrix, and:

$`\left(M_1\right)_{ij}`$ $`=`$ $`\delta _{ij}{\displaystyle \frac{\mathrm{\Delta }t}{2}}{\displaystyle \frac{\overline{\mu }_{ij}}{\xi _i}}`$ (20)
$`\left(M_2\right)_{ij}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }t}{2}}{\displaystyle \frac{1}{\xi _i}}\left(\overline{\mu }_{ij}\mu _{ij}^{}\right)`$ (21)
$`\left(M_3\right)_{ij}`$ $`=`$ $`+\delta _{ij}{\displaystyle \frac{\mathrm{\Delta }t}{2}}{\displaystyle \frac{\mu _{ij}^{}}{\xi _i}}`$ (22)
where $`\delta _{ij}`$ is the $`(i,j)`$ Kronecker symbol. The possible values of $`\lambda `$ for the normal mode structure $`f(𝐫)`$ that we examine, are thus given by the roots of the following polynomial equation in $`\lambda `$:

$$\mathrm{Det}(𝐌)=0$$
(23)
The dependencies to $`\lambda `$ are limited to the top- and bottom-left blocks. For a non-extrapolating 2-TL scheme the degree of the polynomial is $`P`$, and there are $`P`$ physical modes associated with this structure $`f(𝐫)`$. For time-schemes making use of three time-levels (i.e. 3-TL schemes, or extrapolating 2-TL schemes, or 2-TL schemes with a time-filter), the degree becomes $`2P`$, and there are $`P`$ additional computational modes. The growth-rate for any of these modes is given by the modulus of the corresponding complex root of (23). The growth rate of the time scheme for the considered structure $`f`$ is then defined by the maximum value of the modulus of these $`P`$ or $`2P`$ roots:

$$\mathrm{\Gamma }(f)=\mathrm{Max}\left(\left|\lambda _i(f)\right|\right),i(1,\mathrm{},P)\left[\mathrm{or}(1,\mathrm{},2P)\right]$$
(24)
Analytical solution of (23) is not possible for large values of $`P`$, and a numerical solution is often needed.

In this paper we will call "asymptotic growth-rate" the growth-rate for the matrix $`𝐌`$ in the limit of large time-steps $`\mathrm{\Delta }t\mathrm{}`$. The analysis of the asymptotic growth-rate is easier than for finite time-steps, since the matrix $`𝐌`$ of (23) degenerates to a matrix $`𝐌^{}`$ in which $`\delta _{ij}`$ vanishes and $`\mathrm{\Delta }t/2`$ eliminates. Another advantage of asymptotic growth-rates is that they appear to be independent of the structure $`f`$ in most of the cases examined below, thus simplifying considerably the interpretation of the results. When the asymptotic growth-rate is independent of the structure $`f`$ the growth-rate of the scheme can be defined by the growth-rate obtained for this scheme with any structure $`f`$.

When the growth rate for a given scheme is one (or less) for any mode of any normal mode structure $`f`$, the scheme is then said to be "unconditionally stable" (being understood "in $`\mathrm{\Delta }t`$"). The criterion for unconditional stability obtained through $`\mathrm{Det}(𝐌^{})=0`$ is not only of academic interest since the considered time-schemes are actually used with large time-steps in NWP: the practice shows that a scheme which is not unconditionally stable in the simplified context of these analyses has few chance to be robust enough for use in real conditions.

## 6 Simple examples: 1D systems

### 6.1 1D Shallow-water system

The linearised 1D shallow water system in an horizontal direction $`x`$ can be classically written in terms of the wind $`u`$ along $`x`$, and the geopotential $`\varphi `$:

$`{\displaystyle \frac{u}{t}}`$ $`=`$ $`{\displaystyle \frac{\varphi }{x}}`$ (25)
$`{\displaystyle \frac{\varphi }{t}}`$ $`=`$ $`\overline{\varphi }{\displaystyle \frac{u}{x}}`$ (26)
This system is also valid for the external mode of an isothermal atmosphere in HPE and EE system, replacing $`\overline{\varphi }`$ by $`4(R^2/C_p)\overline{T}`$ (which is done in the following). The reference system is obtained by replacement of $`\overline{T}`$ by $`T^{}`$, and a "non-linearity" factor is defined through: $`\alpha =(\overline{T}T^{})/T^{}`$. Solution of (13) implies that if $`\overline{T}<0`$, $`\overline{\lambda }i\mathrm{I}\mathrm{R}u=\widehat{u}\mathrm{exp}(rx)`$, which has not a bounded energy-density for $`rx+\mathrm{}`$. Hence \[C2’\] requires $`\overline{T}0`$ (i.e. $`\alpha 1`$). The boundary conditions do not appear explicitly in the system, hence $`l`$ can be taken as the identity operator to satisfy \[C1\], and \[C3\]. In the notations of section 4 and 5 we have $`P=2`$ , $`𝒳_1=u`$ and $`𝒳_2=\varphi `$. The normal modes of the system write $`\psi (x)=\widehat{\psi }\mathrm{exp}(ikx)`$ with $`k\mathrm{I}\mathrm{R}`$ and $`\psi =(u,\varphi )`$. Conditions \[C1\] – \[C4\] are easily checked to be satisfied.

For a 3-TL SI scheme, (23) writes:

$$\left(\frac{\lambda ^21}{2\mathrm{\Delta }t}\right)^2=\frac{k^2c^2}{4}\left(\lambda ^2+1\right)\left(\lambda ^2+1+2\alpha \lambda \right),$$
(27)
where $`c^{}=2\sqrt{(R/C_p)RT^{}}`$. In the limit of long time-steps, the LHS term disappears, and the four roots of the RHS give the "asymptotic" numerical growth-rate for the two physical and two computational modes of the system. The two roots of the first factor have a neutral stability, while those of the second factor have a modulus equal to 1 if $`1\alpha 1`$. The criterion (on $`\overline{𝒳}`$, $`𝒳^{}`$) for unconditional stability (in $`\mathrm{\Delta }t`$) of the 3-TL SI scheme is thus : $`0\overline{T}2T^{}`$. Some further algebraic manipulations from (23) with $`\mathrm{\Delta }t=\mathrm{}`$ show that this criterion remains unchanged when increasing $`N_{\mathrm{iter}}`$.

For a 2-TL SI non-extrapolating ($`\mu =1`$) scheme, (23) becomes:

$$\left(\frac{\lambda 1}{\mathrm{\Delta }t}\right)^2=\frac{k^2c^2}{4}\left(\lambda +1\right)\left(\lambda +1+2\alpha \right),$$
(28)
and the criterion for unconditional stability becomes more constraining than for the 3-TL scheme: $`1\alpha 0`$ (i.e. $`0\overline{T}T^{}`$). The asymptotic growth-rate of a 2-TL non-extrapolating ICI scheme with $`N_{\mathrm{iter}}`$ iterations is given by:

$$\mathrm{\Gamma }=\mathrm{Max}(1,\left|2(\alpha )^{N_{\mathrm{iter}}}1\right|)$$
(29)
The domain for unconditional stability is thus $`1\alpha 0`$ for odd values of $`N_{\mathrm{iter}}`$, and $`1\alpha 1`$ for even values of $`N_{\mathrm{iter}}`$.

### 6.2 1D vertical acoustic system in mass-based coordinate

We consider a vertical 1D compressible atmospheric column satisfying the conditions listed in section 2. A regular mass-based coordinate $`\sigma =(\pi /\pi _s)`$ is chosen following L92, by making $`\eta =\sigma `$, $`A(\sigma )=0`$ and $`B(\sigma )=\sigma `$ in equation (31) of L92. The variable $`\pi `$ denotes the hydrostatic pressure, and $`\pi _s`$, the surface hydrostatic pressure. The system is readily obtained from equations (36) – (45) of L92 by removing all horizontal dependencies. The surface hydrostatic-pressure does not evolve in time (see equation (45) in L92).

The equations are linearised around a resting atmospheric-state $`\overline{𝒳}`$ and a resting reference-state $`𝒳^{}`$, both satisfying the conditions of section 2. The temperatures $`\overline{T}`$ and $`T^{}`$ are taken uniform, and we still define the "non-linearity" factor by: $`\alpha =(\overline{T}T^{})/T^{}`$. The pressure values $`\overline{p}`$ and $`p^{}`$ are assumed to be equal to a common value $`\pi _0`$ at the origin ($`\sigma =1`$). Since $`\overline{𝒳}`$ and $`𝒳^{}`$ are hydrostatically balanced, $`\overline{p}=\overline{\pi }=\sigma \pi _0`$ and $`p^{}=\pi ^{}=\sigma \pi _0`$ at any level. The thermodynamics equation decouples and the linear system around $`\overline{𝒳}`$ for the vertical velocity $`w`$ , and the pressure deviation $`p^{}=p\overline{p}`$ writes in standard notations:

$`{\displaystyle \frac{w}{t}}`$ $`=`$ $`{\displaystyle \frac{g}{\pi _0}}{\displaystyle \frac{p^{}}{\sigma }}`$ (30)
$`{\displaystyle \frac{p^{}}{t}}`$ $`=`$ $`{\displaystyle \frac{C_p}{C_v}}{\displaystyle \frac{g\pi _0}{R\overline{T}}}\sigma ^2{\displaystyle \frac{w}{\sigma }}`$ (31)
The same derivation holds for $`^{}`$ and leads to an operator formally identical to the RHS of (30)-(31), still acting on $`(w,p^{})`$, but with $`\overline{T}`$ replaced by $`T^{}`$. The solution of (13), implies that if $`\overline{T}<0`$, $`\overline{\lambda }i\mathrm{I}\mathrm{R}w=\widehat{w}\sigma ^r`$ with $`r<1`$ or $`r>0`$. For the mode with $`r<1`$, the energy-density is not bounded when $`\sigma 0`$. If $`\overline{T}0`$, the structure of the normal modes of (30)-(31) is given by:

$`w(\sigma )`$ $`=`$ $`\widehat{w}\sigma ^{(i\nu 1/2)}=\widehat{w}f_1(\sigma )`$ (32)
$`p^{}(\sigma )`$ $`=`$ $`\widehat{p^{}}\sigma ^{(i\nu +1/2)}=\widehat{p^{}}f_2(\sigma )`$ (33)
where $`\nu `$ is a real number, and they have a bounded energy-density. The condition \[C2’\] therefore requires $`\overline{T}0`$ (i.e. $`\alpha 1`$). Finally, \[C4\] is trivially checked to be satisfied.

For a 3-TL SI scheme, (23) writes:

$$\left(\frac{\lambda ^21}{2\mathrm{\Delta }t}\right)^2=\frac{(\nu ^2+1/4)c^2}{4H^2}\left(\lambda ^2+1\right)\left(\lambda ^2+1\frac{2\alpha \lambda }{1+\alpha }\right),$$
(34)
where $`c^{}=\sqrt{(C_p/C_v)RT^{}}`$ and $`H^{}=RT^{}/g`$. Comparison of (34) and (27) shows that the stability of the 1D vertical system for $`\alpha `$ is the same as the one of the previous shallow-water system for $`\alpha ^{}=\alpha /(1+\alpha )`$. Hence the criteria for unconditional stability directly follows from those of the previous case, by similarity arguments. The criterion for unconditional stability of the 3-TL SI scheme is $`\alpha (1/2)`$, i.e. $`\overline{T}(1/2)T^{}`$, and this criterion remains unchanged when increasing $`N_{\mathrm{iter}}`$.

For a 2-TL SI non-extrapolating ($`\mu =1`$) scheme, (23) becomes:

$$\left(\frac{\lambda 1}{\mathrm{\Delta }t}\right)^2=\frac{(\nu ^2+1/4)c^2}{4H^2}\left(\lambda +1\right)\left(\lambda +1\frac{2\alpha }{1+\alpha }\right),$$
(35)
and the criterion for unconditional stability ($`\alpha 0`$, i.e. $`\overline{T}T^{}`$), which is more constraining than for the 3-TL SI scheme. For iterated 2-TL schemes, the criterion for unconditional stability is $`\alpha (1/2)`$ for even values of $`N_{\mathrm{iter}}`$, and $`\alpha 0`$ for odd values.

### 6.3 1D vertical acoustic system in height-based coordinate

In this example we show that the stability properties of the 1D vertical system may depend on the coordinate. The framework is taken as in the previous example except that the vertical coordinate is the height $`z`$. The linearised system $`\overline{}`$ (cf. e.g. Caya and laprise, 1999) writes:

$`{\displaystyle \frac{w}{t}}`$ $`=`$ $`R\overline{T}{\displaystyle \frac{q^{}}{z}}+g{\displaystyle \frac{T^{}}{\overline{T}}}`$ (36)
$`{\displaystyle \frac{T^{}}{t}}`$ $`=`$ $`{\displaystyle \frac{R\overline{T}}{C_v}}{\displaystyle \frac{}{z}}w`$ (37)
$`{\displaystyle \frac{q^{}}{t}}`$ $`=`$ $`\left({\displaystyle \frac{g}{R\overline{T}}}{\displaystyle \frac{C_p}{C_v}}{\displaystyle \frac{}{z}}\right)w,`$ (38)
where $`q^{}=q\overline{q}`$, $`q=\mathrm{ln}(p/p_0)`$, $`\overline{q}=gz/R\overline{T}`$, $`p_0`$ is a reference pressure, and $`p`$ is the true pressure. The normal modes of $`\overline{}`$ have the following form, for $`\psi =(w,T^{},q^{})`$:

$$\psi =\widehat{\psi }\mathrm{exp}\left[\left(i\nu +\frac{1}{2\overline{H}}\right)z\right]$$
(39)
where $`\overline{H}=(R\overline{T}/g)`$. The reference system $`^{}`$ is defined in a similar way replacing $`\overline{T}`$ by $`T^{}`$ and $`\overline{q}`$ by $`q^{}=gz/RT^{}`$. It should be noted that the structure of the normal modes of $`^{}`$ is not the same as for $`\overline{}`$ since the characteristic height for $`^{}`$ is $`H^{}=(RT^{}/g)`$

For a 2-TL SI non-extrapolating ($`\mu =1`$) scheme, (23) becomes:

$$\left(\frac{\lambda 1}{\mathrm{\Delta }t}\right)^2=\frac{c^2}{4}\left(i\nu +\frac{1}{2\overline{H}}\right)\left(\lambda +1+2\alpha \right)\left[i\nu \left(\lambda +1\right)\frac{1}{H^{}}\left(\frac{1+2\alpha }{1+\alpha }\right)\left(\lambda +1\frac{4\alpha }{1+2\alpha }\right)\right]$$
(40)
where $`c^{}=\sqrt{(C_p/C_v)RT^{}}`$. In the height-coordinate framework, the asymptotic growth-rate depends on the structure since $`\nu `$ appears in one of the factors which become dominant at large time-steps. The interpretation of the results is thus slightly complicated in comparison to the case of a mass-based coordinate. For the most external structure ($`\nu =0`$), the asymptotic growth-rate is given by the roots of:

$$\left(\lambda +1+2\alpha \right)\left(\lambda +1\frac{4\alpha }{1+2\alpha }\right)=0.$$
(41)
This polynomial consists in a combination of two factors similar to those obtained in the two previous examples (through a formal replacement of $`2\alpha `$ by $`\alpha `$ for the second factor). As a consequence, the unconditional stability domains can be readily deduced from these previous examples: the external structure $`\nu =0`$ is unstable for any value $`\alpha 0`$ when $`\mathrm{\Delta }t\mathrm{}`$. The instability is thus much more severe than in the case of a mass-based vertical coordinate, for which $`\alpha 0`$ was sufficient to ensure unconditional stability. Moreover, slightly shorter structures with vertical wavelengths of the order of $`(1/\overline{H})`$ are found to be more unstable than the external one, and for these structures, the unconditional stability criterion ($`\alpha =0`$) remains unchanged when $`N_{\mathrm{iter}}`$ is increased. Fig. 1 shows the asymptotic growth-rates for the 2-TL SI $`(N_{\mathrm{iter}}=1)`$ scheme and the the 2-TL ICI scheme with $`N_{\mathrm{iter}}=2`$ for $`\nu =0.0001\mathrm{m}^1`$, a structure for which the instability is close to its maximum. The severe instability of the 2-TL SI scheme is only alleviated but not suppressed by choosing $`N_{\mathrm{iter}}=2`$.

For the 3-TL SI scheme, the external structure $`\nu =0`$ is unconditionally stable for $`0.25\alpha 1`$, but slightly shorter structures as above are found unstable at large time-steps as soon as $`\alpha 0`$ (very short modes are stable however). Fig. 1 depicts the asymptotic growth-rates for two structures: the external structure $`\nu =0`$, and a long structure $`\nu =0.0001\mathrm{m}^1`$. The growth-rate of the long structure for a moderate time-step $`\mathrm{\Delta }t=30`$ s with a time-decentering $`ϵ=0.1`$ (as in Caya and Laprise, 1999) is also depicted: the practical instability becomes small in these conditions, and the 3-TL scheme cannot be positively rejected, especially considering the fact that dissipative processes could act in a way to stabilize the scheme. The practical impact of this predicted weak instability for NWP applications could be easily evaluated with a $`z`$-coordinate model, using an experimental set-up similar to the one used here, and then progressively extending the set-up to approach real-case experimental conditions.

### 6.4 Comments

In the three simple examples examined above, the criterion for unconditional stability is seen to be more constraining in the 2-TL non-extrapolating SI scheme than in the 3-TL SI scheme. The 2-TL extrapolating SI scheme is found to have similar domains of unconditional stability than its non-extrapolating counterpart (not shown). For mass-based coordinates, if both vertically propagating acoustic waves and external gravity waves are simultaneously allowed by a given equation system, the above analyses suggest that 2-TL SI schemes are so constraining that there is no domain for unconditional stability. For height-based coordinates, the long vertically propagating acoustic waves are always unstable in the 2-TL SI scheme.

This leads to suspect that, in opposition to 3-TL SI schemes, classical 2-TL SI schemes are not suitable for the EE system with any vertical coordinate. The HPE system with 2-TL SI scheme did not suffer from this problem since vertically propagating waves are not allowed in HPE (i.e. the 1D column atmosphere is stationary). The intrinsic instability of the 2-TL SI scheme for EE system is confirmed in section 7 for mass-based coordinates.

When a second-order time-filter with parameter $`\kappa `$ is applied to the system examined in the first example for a 2-TL SI non-extrapolating scheme, (28) becomes:

$$\left(\frac{\lambda 1}{\mathrm{\Delta }t}\right)^2=\frac{k^2c^2}{4}\left[\left(\lambda +1\right)+\kappa \left(\lambda 2+\frac{1}{\lambda }\right)\right]\left[\left(\lambda +1+2\alpha \right)+\kappa \left(\lambda 2+\frac{1}{\lambda }\right)\right],$$
(42)
and the criterion for unconditional stability becomes $`1\alpha 2\kappa `$. For the second example, the similarity argument shows that the criterion for unconditional stability becomes: $`\alpha 2\kappa /(1+2\kappa )`$. The domains of stability of 3-TL and 2-TL time-filtered ICI schemes for the two first examples are summarized in Table 1.

The application of a time-filter thus allows to alleviate the stability constraints for 2-TL SI schemes, and a non-vanishing domain for unconditional stability is recovered. The width of the unconditional stability domain increases with $`\kappa `$. This is found to hold for the 1D vertical system in height-based coordinates as well, which is consistent with the results of Semazzi et al. (1995) and Qian et al (1998): they succeeded to solve numerically the EE system at low-resolution with a 2-TL SI scheme, however, the use of a large value $`\kappa =0.5`$ was required to stabilize the model. As a consequence the forecasts suffered from a dramatic loss of energy with increasing forecast-range, and ceased to be of meteorological interest after 2-3 days. Moreover, at high resolutions (and consequently steep orography) the use of a time-filter $`\kappa `$ is found experimentally to be an insufficient solution for eliminating the intrinsic instability of the scheme (not shown).

If a high level of accuracy is desired for the EE system with a 2-TL classical SI time-discretization and high resolution, a more robust scheme (e.g. with a larger value of $`N_{\mathrm{iter}}`$) must be used. The above analyses show that the unconditional stability domain is dramatically reduced for odd values of $`N_{\mathrm{iter}}`$, hence ICI schemes with even values of $`N_{\mathrm{iter}}`$ are preferable for solving the EE system with a 2-TL scheme.

The 1D vertical system in mass-based coordinates has been found to be more stable than its counterpart in height-based coordinates in a general way. For mass-coordinates, the 3-TL SI and the 2-TL ICI schemes with even values of $`N_{\mathrm{iter}}`$ have an extended domain of unconditional stability, whilst for height-based coordinates, they are unstable as soon as $`\alpha 0`$. For 3-TL SI schemes, the necessity to have recourse to a first-order time-decentering $`ϵ>0`$ to overcome this instability is a significant drawback since it results in a spurious damping of transient phenomena, similarly to $`\kappa `$ but in an even less selective way, as mentioned above. We think these differences give a substantial theoretical advantage to mass-based coordinates for solving the EE system with classical SI or ICI schemes.

## 7 Analysis of the EE system for isothermal atmospheres in mass-coordinate

The analysis of the isothermal HPE system for ICI schemes does not substantially modify the general conclusions drawn for the shallow-water case (not shown), hence the case of the EE system is directly examined.

In this section, the EE system is cast in the pure unstretched terrain-following coordinate $`\sigma `$ which can be classically derived from the hydrostatic-pressure coordinate $`\pi `$ of L92, through $`\sigma =(\pi /\pi _s)[0,1]`$, where $`\pi _s`$ is the hydrostatic surface-pressure. The nonhydrostatic prognostic variables are the non-dimensionalised nonhydrostatic pressure departure $`𝒫=(p\pi )/\pi `$ (where $`p`$ is the true pressure), and the vertical divergence $`𝖽`$ which writes in $`\sigma `$ coordinate:

$$𝖽=\frac{g}{RT}(1+𝒫)\sigma \frac{w}{\eta }$$
(43)
The adiabatic system writes:

$`{\displaystyle \frac{d𝐕}{dt}}`$ $`=`$ $`RTq{\displaystyle \frac{RT}{(1+𝒫)}}𝒫\left(1+𝒫+\sigma {\displaystyle \frac{𝒫}{\sigma }}\right)\varphi `$ (44)
$`{\displaystyle \frac{d𝖽}{dt}}`$ $`=`$ $`{\displaystyle \frac{g^2(1+𝒫)}{RT}}\left(\sigma {\displaystyle \frac{}{\sigma }}\right)\left(1+\sigma {\displaystyle \frac{}{\sigma }}\right)𝒫`$ (45)
$`+`$ $`𝖽(𝐕D_3)+{\displaystyle \frac{g(1+𝒫)}{RT}}\left[w\left(\sigma {\displaystyle \frac{𝐕}{\sigma }}\right)\right]`$
$`{\displaystyle \frac{dT}{dt}}`$ $`=`$ $`{\displaystyle \frac{RT}{C_v}}D_3`$ (46)
$`{\displaystyle \frac{d𝒫}{dt}}`$ $`=`$ $`\left(1+𝒫\right)\left({\displaystyle \frac{C_p}{C_v}}D_3+{\displaystyle \frac{\dot{\pi }}{\pi }}\right)`$ (47)
$`{\displaystyle \frac{q}{t}}`$ $`=`$ $`{\displaystyle _0^1}\left(𝐕+𝐕q\right)𝑑\sigma ^{}`$ (48)
where:

$`D_3`$ $`=`$ $`𝐕+𝖽+{\displaystyle \frac{(1+𝒫)}{RT}}\varphi .\left(\sigma {\displaystyle \frac{𝐕}{\sigma }}\right)`$ (49)
$`\varphi `$ $`=`$ $`R{\displaystyle _\sigma ^1}\left({\displaystyle \frac{T}{1+𝒫}}\right){\displaystyle \frac{d\sigma ^{}}{\sigma ^{}}}`$ (50)
$`{\displaystyle \frac{\dot{\pi }}{\pi }}`$ $`=`$ $`𝐕q{\displaystyle \frac{1}{\sigma }}{\displaystyle _0^\sigma }\left(𝐕+𝐕q\right)𝑑\sigma ^{},`$ (51)
V is the horizontal wind, and $``$ is the horizontal derivative operator. The domain is restricted to a vertical plane along $`(x,\sigma )`$ directions for clarity. The system is linearized around a resting isothermal and hydrostatically-balanced state $`\overline{𝒳}`$:

$`{\displaystyle \frac{D}{t}}`$ $`=`$ $`R𝒢^2T+R\overline{T}(𝒢)^2𝒫R\overline{T}^2q`$ (52)
$`{\displaystyle \frac{𝖽}{t}}`$ $`=`$ $`{\displaystyle \frac{g^2}{R\overline{T}}}\left(1+\sigma {\displaystyle \frac{}{\sigma }}\right)\left(\sigma {\displaystyle \frac{}{\sigma }}\right)𝒫`$ (53)
$`{\displaystyle \frac{T}{t}}`$ $`=`$ $`{\displaystyle \frac{R\overline{T}}{C_v}}(D+𝖽)`$ (54)
$`{\displaystyle \frac{𝒫}{t}}`$ $`=`$ $`𝒮D{\displaystyle \frac{C_p}{C_v}}(D+𝖽)`$ (55)
$`{\displaystyle \frac{q}{t}}`$ $`=`$ $`𝒩D,`$ (56)
where the vertical integral operators $`𝒢`$, $`𝒮`$ and $`𝒩`$ are defined by:

$`𝒢X`$ $`=`$ $`{\displaystyle _\sigma ^1}(X/\sigma ^{})𝑑\sigma ^{}`$ (57)
$`𝒮X`$ $`=`$ $`(1/\sigma ){\displaystyle _0^\sigma }X𝑑\sigma ^{}`$ (58)
$`𝒩X`$ $`=`$ $`{\displaystyle _0^1}X𝑑\sigma ^{}`$ (59)
The $`^{}`$ operator is similar to the RHS of this system, simply replacing $`\overline{T}`$ by $`T^{}`$.

### 7.1 Verification of conditions \[C1\] – \[C4\]

The linear operator $`l_1=\sigma (/\sigma )`$ is applied to (52), and $`l_4=[+\sigma (/\sigma )]`$ to (55). The $`q`$ equation (56) decouples and we obtain a linear unbounded system, in which (52) and (55) are replaced by:

$`\left(\sigma {\displaystyle \frac{}{\sigma }}\right){\displaystyle \frac{D}{t}}`$ $`=`$ $`R^2TR\overline{T}\left(\sigma {\displaystyle \frac{}{\sigma }}+\right)^2𝒫`$ (60)
$`\left(\sigma {\displaystyle \frac{}{\sigma }}+\right){\displaystyle \frac{𝒫}{t}}`$ $`=`$ $`D{\displaystyle \frac{C_p}{C_v}}\left(\sigma {\displaystyle \frac{}{\sigma }}+\right)\left(D+𝖽\right)`$ (61)
Hence we have $`P=4`$, $`𝒳=(D,𝖽,T,𝒫)`$. Using of the same operators ($`l_1`$, $`l_4`$), the reference operator is also made free of any reference to the upper and lower boundary conditions, which shows that the condition \[C1\] is satisfied. Solution of (13) shows that \[C2\] requires $`\overline{T}0`$ (i.e. $`\alpha 1`$). The normal modes of the system are then:

$$\psi (x,\sigma )=\widehat{\psi }\mathrm{exp}(ikx)\sigma ^{(i\nu 1/2)}$$
(62)
where $`(k,\nu )\mathrm{I}\mathrm{R}`$ and $`\psi `$ represents $`D`$, $`𝖽`$, $`T`$ or $`𝒫`$. In this particular case, the $`f_1,\mathrm{},f_4`$ functions are all identical. The verification of \[C3\], \[C4\] proceeds easily, as in previous sections.

### 7.2 Results

As a first illustration of the results, the growth-rates of 2-TL non-extrapolating ICI schemes are shown in Fig. 2 as function of $`\alpha =(\overline{T}T^{})/T^{}`$ with a moderate time-step $`\mathrm{\Delta }t=20s`$, for three particular mode structures:

(i) an external mode ($`k=0.0005\mathrm{m}^1`$, $`\nu =0`$)

(ii) a vertical very internal mode ($`k=0`$, $`\nu =100`$)

(ii) an intermediate slantwise mode ($`k=0.0005\mathrm{m}^1`$, $`\nu =3`$)

The suspicions raised in the simple 1D examples are confirmed: the internal vertically-propagating mode is unstable for $`\alpha <0`$ while the external gravity mode is unstable for $`\alpha >0`$. Moreover, intermediate, slantwise-propagating modes are unstable in the whole domain, and the acoustic external mode (Lamb-wave) appears to be unstable for $`\alpha <0`$ as well. The domain of stability vanishes, which confirms that 2-TL SI scheme is not relevant for solving the EE system. The effect of introducing a time-filter to save the situation is discussed below.

The asymptotic growth-rates resulting from the EE system for $`\mathrm{\Delta }t=\mathrm{}`$ are now examined. Similarly to most previous cases, they are independent of the geometry ($`k`$, $`\nu `$) of the mode. Fig. 3 shows the asymptotic growth-rates as a function of $`\alpha `$ for 2-TL non-extrapolating ICI schemes with $`N_{\mathrm{iter}}=(1,2,3,4)`$. As stated above, the SI scheme ($`N_{\mathrm{iter}}=1`$) is unstable for any value of $`\alpha `$. For even values of $`N_{\mathrm{iter}}`$, the scheme has an "optimal" domain of unconditional stability $`1/2\alpha 1`$, while for odd values, the scheme is unstable for all values of $`\alpha `$.

For 3-TL ICI schemes the domain of unconditional stability is $`1/2\alpha 1`$ independently of the values of $`N_{\mathrm{iter}}`$ and $`\kappa `$. The curves (not shown) are similar to those obtained for even values of $`N_{\mathrm{iter}}`$ for 2-TL ICI schemes.

The impact of applying a time-filter $`\kappa =0.1`$ to 2-TL ICI schemes is depicted on Fig 4 for the second variant (the first variant behave qualitatively in the same way). The global impact is to "lower" the curves of the asymptotic growth-rates, and consequently, to increase the width of the unconditional-stability domain. However, large values of $`\kappa `$ (e.g. $`\kappa 0.5`$) are required in order to obtain a wide stability domain, especially for the 2-TL SI scheme, and this strategy is known to be irrelevant for NWP purpose.

Finally, it is worth noting that the results obtained for the EE system are fully compatible with the conclusions that can be drawn from the intersection of the domains of unconditional stability in Table 1 for the two simple frameworks examined above in mass-coordinate. The ability of these very refined frameworks to capture the essence of the behaviour of the time-discretised EE system in the limit of long time-steps makes them very useful tools to fully understand the underlying causes of its stability or instability.

## 8 Conclusion

A general method for investigating the stability of the ICI class of time-discretisations on canonical problems with various space-continuous equation systems has been presented. These ICI schemes are based on a separation of evolution terms between a simple linear operator and "non-linear" residuals. The method has been validated by confirming earlier results, then the application to new frameworks (equation systems or time-discretisation schemes) allowed to extend these results. The main conclusions drawn from this study are:

* Even on very simple (1D) examples, the stability properties of time-discretisations for a given equation system are very dependent on fundamental choices (as e.g. the choice of the vertical coordinate). Hence, the import of conclusions drawn from a given analysis must be carefully limited to the examined framework.
* For the EE system, height-based coordinates have a theoretical disavantage compared to mass-based coordinates since they exhibit an intrinsic instability for (long) vertically propagating waves.
* The 2-TL SI scheme is found not to be appropriate for the EE set of equations, whatever coordinate is employed (using a time-filter results in an unacceptable degradation of the solution).
* For the EE system, the 2-TL scheme with $`N_{\mathrm{iter}}=2`$ brings a dramatic increase of the stability compared to the 2-TL SI scheme ($`N_{\mathrm{iter}}=1`$). This statement holds for even values of $`N_{\mathrm{iter}}`$, while odd values leads to a significantly weaker stability.
* As a consequence of the latter point, the 2-TL ICI scheme with $`N_{\mathrm{iter}}=2`$ seems worth to be considered for the EE system.

However, as mentioned in SHB78, the stability inferred from this type of analyses is overestimated, and flows in which the non-linearity comes from other sources than the discrepancy between the atmospheric and reference temperature profiles could reveal new instabilities in practice. For instance, in spite of its apparent "optimal" stability in the simplified context of this paper, the 3-TL SI scheme has proved to be not stable enough for solving numerically the EE system in realistic highly non-linear conditions at high resolutions, due to other terms treated explicitly. This point clearly demonstrates the limitations of this type of academic exercise.

Nevertheless, in spite of its necessary limitations, this study can serve to distinguish schemes which are definitely not relevant for practical use from the others, and give a first theoretical justification for those which are worth considering. In agreement with Cote et al (1998) and Cullen (2000) we think that ICI schemes with $`N_{\mathrm{iter}}2`$ are among the most approriate alternatives for integrating the EE system in highly non-linear conditions at fine-scale, including from the point of view of efficiency.

## References

Bubnová, R., G. Hello, P. Bénard, and J.F. Geleyn, 1995: Integration of the fully elastic equations cast in the hydrostatic pressure terrain-following coordinate in the framework of the ARPEGE/Aladin NWP system. Mon. Wea. Rev., 123, 515-535.

Caya, D., and R. Laprise, 1999: A semi-implicit semi-Lagrangian regional climate model: the Canadian RCM. Mon. Wea. Rev., 127, 341-362.

Côté, J., M. Béland, and A. Staniforth, 1983: Stability of vertical discretization schemes for semi-implicit primitive equation models: theory and application. Mon. Wea. Rev., 111, 1189-1207.

Côté, J., S. Gravel, A. Méthot, A. Patoine, M. Roch, and A. Staniforth, 1998: The Operational CMC-MRB Global Environmental Multiscale (GEM) Model. Part I: Design Considerations and Formulation. Mon. Wea. Rev., 126, 1373-1395.

Cullen, M. J. P., 2000: Alternative implementations of the semi-Lagrangian semi-implicit schemes in the ECMWF model. Q. J. R. Meteorol. Soc., 127, 2787-2802.

Hereil, P., and R. Laprise, 1996: Sensitivity of Internal Gravity Waves Solutions to the Time Step of a Semi-Implicit Semi-Lagrangian Nonhydrostatic Model. Mon. Wea. Rev., 124, 972-999.

Laprise, R., 1992: The Euler equations of motion with hydrostatic pressure as an independent variable. Mon. Wea. Rev., 120, 197-207.

Qian, J.-H., F. H. M. Semazzi, and J. S. Scroggs, 1998: A global nonhydrostatic semi-Lagrangian atmospheric model with orography. Mon. Wea. Rev., 126, 747-771.

Robert, A. J., J. Henderson, and C. Turnbull, 1972: An implicit time integration scheme for baroclinic models of the atmosphere . Mon. Wea. Rev., 100, 329-335.

Semazzi, F. H. M., J. H. Qian, and J. S. Scroggs, 1995: A global nonhydrostatic semi-Lagrangian atmospheric model without orography. Mon. Wea. Rev., 123, 2534-2550.

Simmons, A. J., B. Hoskins, and D. Burridge, 1978: Stability of the semi-implicit method of time integration. Mon. Wea. Rev., 106, 405-412.

Simmons, A. J., C. Temperton, 1997: Stability of a two-time-level semi-implicit integration scheme for gravity wave motion. Mon. Wea. Rev., 125, 600-615.

Tanguay, M., A. Robert, and R. Laprise, 1990: A Semi-Implicit Semi-Larangian Fully Compressible Regional Forecast Model. Mon. Wea. Rev., 118, 1970-1980.

## List of Tables

Table 1: Domains of unconditional stability for time-filtered schemes for the two first 1D examples.

|  | Shallow-water | 1D vertical (mass) |
| --- | --- | --- |
| 3-TL ICI | $`1\alpha 1`$ | $`1/2\alpha `$ |
| 2-TL ICI ($`N_{\mathrm{iter}}`$ even) | $`1\alpha 1`$ | $`1/2\alpha `$ |
| 2-TL ICI ($`N_{\mathrm{iter}}`$ odd) | $`1\alpha (2\kappa )^{1/N_{\mathrm{iter}}}`$ | $`\frac{2\kappa ^{\left(1/N_{\mathrm{iter}}\right)}}{1+2\kappa ^{\left(1/N_{\mathrm{iter}}\right)}}\alpha `$ |

Table 1: Domains of unconditional stability for time-filtered schemes for the two first 1D examples.

## List of Figures

Fig. 1: Asymptotic growth-rates $`\mathrm{\Gamma }`$ for 1D vertical system in $`z`$ coordinate as a function of the nonlinearity parameter $`\alpha `$. thin line: long mode ($`\nu =0.0001\mathrm{m}^1`$) with 2-TL SI scheme; thick line: long mode with 2-TL ICI scheme $`N_{\mathrm{iter}}=2`$; dotted line: external mode ($`\nu =0`$)with 3-TL SI scheme; dashed line: long mode with 3-TL SI scheme. Circles: practical growth-rate of 3-TL SI scheme for the long mode with $`\mathrm{\Delta }t=30`$ s and $`ϵ=0.1`$.

Fig. 2: Growth-rate $`\mathrm{\Gamma }`$ with $`\mathrm{\Delta }t=20s`$ for the EE system with a 2-TL SI scheme as a function of the nonlinearity parameter $`\alpha `$. solid line: external mode (i); dashed line: slantwise mode (ii); dot-dashed line: internal mode (iii). The left-part of the solid line represents an acoustic external mode.

Fig. 3: Asymptotic growth-rates $`\mathrm{\Gamma }`$ for EE system with 2-TL ICI scheme as a function of the nonlinearity parameter $`\alpha `$. solid line: $`N_{\mathrm{iter}}=1`$; dashed line: $`N_{\mathrm{iter}}=2`$; dot-dashed line: $`N_{\mathrm{iter}}=3`$; dotted line: $`N_{\mathrm{iter}}=4`$.

Fig. 4: Same as Fig 3, but with a time-filter $`\kappa =0.1`$. solid line: $`N_{\mathrm{iter}}=1`$; dashed line: $`N_{\mathrm{iter}}=2`$; dot-dashed line: $`N_{\mathrm{iter}}=3`$; dotted line: $`N_{\mathrm{iter}}=4`$.


text:
# Spin alignment of the high 𝑝_𝑇 vector mesons in polarized 𝑝𝑝 collisions at high energies

## I Introduction

Spin effects in high energy fragmentation processes have attracted much attention recently \[See, e.g., Refs. \[1-12\] and the references given therein\]. One of the important issues is the study of spin transfer from the fragmenting quark to the produced hadrons. It can be studied by measuring the polarizations of the hadrons produced in the fragmentation of a polarized quark. Hyperon polarization has been frequently used in this connection. On the other hand, it has been pointed out that the vector meson polarization can also be used for this purpose since the polarization of a vector meson can also be studied easily in experiments. Furthermore, compared with hyperon, the production rate of a vector meson is usually higher and the origin is simpler due to the much smaller contribution from the decay of heavier hadrons. This implies less theoretical uncertainties in calculations and better statistics in experiments. Therefore, the study of polarization of the vector mesons may provide us important information to understand the role of spin in hadronic reactions.

The polarization of a vector meson is described by the spin density matrix $`\rho `$ or its element $`\rho _{m,m^{}}`$, where $`m`$ and $`m^{}`$ label the spin component along the quantization axis, usually the $`z`$-axis of the frame. The diagonal elements $`\rho _{11}`$, $`\rho _{00}`$ and $`\rho _{11}`$ for the unit-trace matrix are the relative intensities of meson spin component $`m`$ to take the values $`1`$, $`0`$, and $`1`$ respectively. In experiment, $`\rho _{00}`$ can be measured from the angular distribution of its decay products and a deviation of $`\rho _{00}`$ from $`1/3`$ indicates spin alignment. In the helicity basis, i.e., in the case that the $`z`$-axis is chosen as the moving direction of the vector meson, the matrix is usually called helicity density matrix and $`\rho _{00}`$ represents the probability of the vector meson to be in the helicity zero state. Measurements have been carried out in different reactions in particular in $`e^+e^{}`$ annihilation at LEP recently . The data show that $`\rho _{00}`$ is significantly larger than $`1/3`$ for vector mesons with large $`z`$ and increases with increasing $`z`$ (here $`z2p_V/\sqrt{s}`$, where $`p_V`$ is the momentum of the vector meson, $`\sqrt{s}`$ is the total $`e^+e^{}`$ center of mass energy). This means that there exist a significant spin alignment for vector mesons in helicity frame in $`e^+e^{}`$ annihilation at the $`Z^0`$ pole, and the effect is more significant for the large momentum fraction region. The simple statistic model based on spin counting gives the result of $`\rho _{00}=1/3`$. A QCD-inspired model is given in Ref. \[References\], where a fast quark combines a soft anti-quark to form a vector meson and the soft anti-quark preferentially has the same helicity as the fast one. This model leads to $`\rho _{00}=0`$. Another model describes the production of vector meson through the channel $`qqV`$, with the vector meson coupling to the quark like a vector current. In this case, the helicity-conserving vector current gives rise to $`\rho _{00}=1`$. None of these predictions is consistent with the data.

In a recent paper, we calculated the helicity density matrix of vector mesons inclusively produced in $`e^+e^{}`$ annihilation at $`Z^0`$ pole by taking the spin of a vector meson which contains a fragmenting quark as the sum of the spin of the polarized fragmenting quark and that of the anti-quark created in the fragmentation process. By comparing the obtained results with the data, we showed that the experimental results for $`\rho _{00}`$ imply a significant polarization for the anti-quark which is created in the fragmentation process and combines with the fragmenting quark to form the vector meson. It should be polarized in the opposite direction as that of the fragmenting quark and the polarization can approximately be written as,

$$P_z=\alpha P_f,$$
(1)
where $`\alpha 0.5`$ is a constant; $`P_z`$ is the polarization of the anti-quark in the moving direction of the fragmenting quark and $`P_f`$ is the longitudinal polarization of the fragmenting quark of flavor $`f`$. The results of $`\rho _{00}`$ obtained from above relation are in good agreement with the data for different mesons. Since the calculations are rather straight-forward, the relation given by Eq. (1) can be considered as a direct implication of the data in $`e^+e^{}`$ annihilation. It implies that there exists a spin interaction between the fragmenting quark and the anti-quark during the hadronization process. It should be interesting to extend the calculations to other processes and make further tests to see whether the relation is universal, in the sense that it is true for quark fragmentation in all different reactions. Hence we made similar calculations for polarized lepton-nucleon deep-inelastic scattering (DIS) process and polarized $`pp`$ collisions. The results for DIS have been given in Ref.\[References\]. In this paper, we present the calculations and the results for high $`p_T`$ vector mesons in polarized $`pp`$ collisions, where factorization theorem can be used so that fragmentation effects and other effects such as those from the structure functions can be studied separately.

Compared with those in the lepton-induced reactions such as $`e^+e^{}`$ annihilations and DIS, the study of high $`p_T`$ vector mesons in $`pp`$ collisions has the following properties. (1) The involved hard scatterings are strong interaction processes rather than the electroweak processes in the lepton-induced reactions, so the corresponding cross section should in general be larger. Furthermore, the luminosity of the incoming proton beams can in general be made higher than that for leptons. Hence, the statistics in experiment can be better. (2) Not only longitudinally but also transversely polarized quarks can be produced, so that we can study the properties in both cases. (3) There are many different hard subprocesses which contribute to the high $`p_T`$ hadron production. This makes the study more interesting and, at the same time, more complicated than those in lepton-induced reactions. For example, the gluon fragmentation is also involved here, which is unclear even for the unpolarized case. It is very important to know whether its contribution is significant in different kinematic regions. We will first make an estimation of it in next section using a Monte Carlo event generator and then present the calculation method of spin alignment of vector mesons in $`pp`$ collisions. In section III, we calculate the spin alignment of vector mesons with high $`p_T`$ in the helicity frame in $`pp`$ collisions with one longitudinally polarized proton beam and present the results. In Section IV, we extend the study to the the spin alignment of high $`p_T`$ vector mesons for the case that one proton beam is transversely polarized. Finally, a short summary is given in Section V.

## II The calculation method

In this section, we first summarize the calculation method of spin alignment of vector mesons that was used in $`e^+e^{}`$ annihilation, and then extend it to the inclusive high $`p_T`$ vector meson production in longitudinally polarized $`pp`$ collisions.

### II.1 For $`e^+e^{}`$ annihilation process

The calculation method for vector mesons produced in $`e^+e^{}`$ annihilation has been presented in Ref.\[References\]. The main points are summarized as follows. For $`e^+e^{}`$ annihilation at a given energy, we need to consider vector mesons produced in the fragmentation of a polarized initial quark $`q_f^0`$ of flavor $`f`$ and those of polarized initial anti-quark $`\overline{q}_f^0`$. The polarization of $`q_f^0`$ or $`\overline{q}_f^0`$ is a given constant which can be calculated using the standard model for electroweak interaction. To calculate the spin density matrices of these vector mesons, we divide them into the following two groups and consider them separately: (A) those which contain the fragmenting quark $`q_f^0`$; (B) those which don’t contain the fragmenting quark. The spin density matrix $`\rho ^V(z)`$ for the vector meson $`V`$ is given by:

$$\rho ^V(z)=\underset{f}{}\frac{n(z|A,f)}{n(z)}\rho ^V(A,f)+\frac{n(z|B)}{n(z)}\rho ^V(B),$$
(2)
where $`n(z|A,f)`$ and $`\rho ^V(A,f)`$ are respectively the average number and spin density matrix of the vector mesons from group (A); $`n(z|B)`$ and $`\rho ^V(B)`$ are those from (B); $`n(z)`$= $`_fn(z|A,f)`$\+ $`n(z|B)`$ is the total number of vector mesons and $`z`$ is the momentum fraction of the initial quark’s momentum carried off by the vector meson. The average numbers $`n(z|A,f)`$ and $`n(z|B)`$ are independent of the spin properties and can be calculated using an event generator based on a fragmentation model which gives a good description of the unpolarized data. We used generator $`\mathrm{pythia}`$ in our calculations.

In terms of the usually used cross section and fragmentation functions, Eq.(2) has the following form,

$$\rho ^V(z)=\frac{_fD_{V/f}^A(z)\rho ^V(A,f)\sigma _{e^+e^{}q_f^0\overline{q}_f^0}}{_fD_{V/f}(z)\sigma _{e^+e^{}q_f^0\overline{q}_f^0}}+\frac{_fD_{V/f}^B(z)\rho ^V(B)\sigma _{e^+e^{}q_f^0\overline{q}_f^0}}{_fD_{V/f}(z)\sigma _{e^+e^{}q_f^0\overline{q}_f^0}},$$
(3)
where $`D_{V/f}^A(z)`$ and $`D_{V/f}^B(z)`$ are the fragmentation functions of the initial quark $`q_f^0`$ (anti-quark) into the vector mesons of group (A) and (B) with momentum fraction $`z`$ respectively; $`D_{V/f}(z)=D_{V/f}^A(z)+D_{V/f}^B(z)`$ is the usual fragmentation function; $`\sigma _{e^+e^{}q_f^0\overline{q}_f^0}`$ is the cross section of the process $`e^+e^{}q_f^0\overline{q}_f^0`$. We see that, $`D_{V/f}^A(z)\sigma _{e^+e^{}q_f^0\overline{q}_f^0}`$ is just the cross section of the production of the meson $`V`$ which contains the initial quark via the fragmentation of $`f`$-flavor quark in $`e^+e^{}`$ annihilation, i.e., $`e^+e^{}q_f^0\overline{q}_f^0V(q_f^0\overline{q})X`$. It is proportional to the average number $`n(z|A,f)`$ in Eq. (2). The denominator is the total cross section of inclusive vector meson production, which corresponds to $`n(z)`$ of Eq. (2).

There are many different possibilities for the production of of the vector mesons in group (B), we take them as unpolarized, i.e., $`\rho ^V(B)`$=1/3. For those vector mesons from group (A), i.e., those which contain $`q_f^0`$ and an anti-quark $`\overline{q}`$ created in the fragmentation, the spin is taken as the sum of the spins of $`q_f^0`$ and $`\overline{q}`$. The polarization of $`q_f^0`$ is taken as the same as that before the fragmentation. Then the spin density matrix $`\rho ^V(A,f)`$ can be calculated from the direct product of the spin density matrix $`\rho ^{q_f^0}`$ for $`q_f^0`$ and $`\rho ^{\overline{q}}`$ for $`\overline{q}`$. Transforming the direct product, $`\rho ^{q_f^0\overline{q}}`$=$`\rho ^{q_f^0}`$$``$$`\rho ^{\overline{q}}`$, to the coupled basis $`|s,s_z`$ (where $`\stackrel{}{s}`$=$`\stackrel{}{s}^q`$+$`\stackrel{}{s}^{\overline{q}}`$), we can obtain the spin density matrix of the vector meson $`\rho ^V(A,f)`$. In the helicity frame of $`q_f^0`$, i.e., $`z`$-axis is taken as the moving direction of $`q_f^0`$, the density matrix of the vector meson of type (A) is given by,

$$\rho ^V(A,f)=\frac{1}{3+P_fP_z}\left(\begin{array}{ccc}\left(1+P_f\right)\left(1+P_z\right)& \frac{1+P_f}{\sqrt{2}}\left(P_xiP_y\right)& 0\\ \frac{\left(1+P_f\right)}{2}\left(P_x+iP_y\right)& \left(1P_fP_z\right)& \frac{1P_f}{\sqrt{2}}\left(P_xiP_y\right)\\ 0& \frac{1P_f}{\sqrt{2}}\left(P_x+iP_y\right)& \left(1P_f\right)\left(1P_z\right)\end{array}\right),$$
(4)
where $`P_f`$ is the longitudinal polarizations of $`q_f^0`$ and $`\stackrel{}{P}=(P_x,P_y,P_z)`$ is the polarization vector of the anti-quark $`\overline{q}`$. Hence, the 00-component of the density matrix takes the following simple form,

$$\rho _{00}^V(A,f)=(1P_fP_z)/(3+P_fP_z),$$
(5)
in which there is only one unknown variable $`P_z`$, the polarization of $`\overline{q}`$ along the moving or polarization direction of $`q_f^0`$.

Using Eqs. (2) and (5), we can determine $`P_z`$ in different cases by fitting the data in the $`e^+e^{}`$ annihilation at $`Z^0`$ pole for the production of different vector mesons. As mentioned in the introduction, we found that, the data for different vector mesons can reasonably be fitted if we take the $`P_z`$ satisfying the relation given by Eq.(1). Now, we insert Eq.(1) into Eq. (5) and obtain that,

$$\rho _{00}^V(A,f)=(1+\alpha P_f^2)/(3\alpha P_f^2).$$
(6)
Finally, from Eqs. (2) and (6), we have

$$\rho _{00}^V(z)=\underset{f}{}\frac{1+\alpha P_f^2}{3\alpha P_f^2}\frac{n(z|A,f)}{n(z)}+\frac{1}{3}\frac{n(z|B)}{n(z)}.$$
(7)
### II.2 Polarization of the outgoing quark in the hard subprocess in $`pp`$ collisions

Now we come to the high $`p_T`$ vector meson production in polarized $`pp`$ collisions. Such vector mesons come mainly from the fragmentation of outgoing partons in the hard scattering subprocesses. The outgoing parton can be a gluon or a quark, and the cross section of gluon-involved subprocess is even larger than others. However, the gluon distribution and fragmentation functions are both poorly known yet, especially for the polarized case. Fortunately, the momentum fractions carried by the gluons in a proton are usually very small, so the gluon contributions to very high $`p_T`$ hard scattering subprocesses are suppressed. In addition, it is known that the gluon fragmentation function is softer than the quark fragmentation function. Consequently, for the final hadron production with high $`p_T`$, the contribution from gluon fragmentation is much smaller than that from quark fragmentation. To see it numerically, we calculate these contributions using $`\mathrm{pythia}`$. In Fig. 1, we show the results for inclusive $`\rho ^0`$ production with $`p_T`$$`>`$13 GeV at $`\sqrt{s}`$=500 GeV in $`pp`$ collisions. We see that, at such high $`p_T`$ cutoff, the contribution of gluon fragmentation is indeed much smaller than that of quark fragmentation, especially for the large $`\eta `$ region ($`\eta `$ is the pseudo-rapidity of the vector meson). Thus, as we did in studying hyperon polarization in high $`p_T`$ jets, we neglect the contribution of polarized gluon fragmentation to the spin alignment of vector mesons, and consider only the contribution of polarized quark’s fragmentation.

For quark fragmentation, we can use the same method of calculating spin alignment as that used in $`e^+e^{}`$ annihilation. To do it, we need to know the polarization of the outgoing quark in the hard subprocess. Such quark is also longitudinally polarized when one proton beam is longitudinally polarized and its polarization can be calculated in pQCD. The obtained polarization can be written as,

$$P_q(x_a,Q^2,y)=D_L(y)P_q^{in}(x_a,Q^2)$$
(8)
where $`D_L(y)`$ is the longitudinal polarization transfer factor in the scattering from the incoming parton to the outgoing parton; $`D_L(y)`$ is only a function of $`y`$ which is defined as $`yk_b(k_ak_c)/k_ak_b`$ for the subprocess $`abcd`$, where $`k_a`$, $`k_b`$, $`k_c`$ and $`k_d`$ are the four-momenta of the partons $`a`$, $`b`$, $`c`$ and $`d`$ respectively. $`D_L(y)`$ has been calculated using pQCD and the leading order results can be found in many publications, for example, in Ref.\[References\]. We take the parton $`a`$ as the parton from the polarized proton beam, and the polarization of the incoming parton $`a`$ is given by $`P_q^{in}(x_a,Q^2)=\mathrm{\Delta }q(x_a,Q^2)/q(x_a,Q^2)`$, where $`\mathrm{\Delta }q(x_a,Q^2)`$ and $`q(x_a,Q^2)`$ are the helicity and unpolarized distribution functions at momentum fraction $`x_a`$ and scale $`Q^2`$. Therefore, the longitudinal polarization $`P_q`$ is the function of variables $`x_a`$, $`Q^2`$ and $`y`$. Since we are now discussing $`pp`$ collisions with one polarized beam, it is clear that the outgoing quarks with positive momenta in $`z`$ direction have a much larger probability of being polarized than those with negative momenta when the positive axis of $`z`$ is taken as the direction of the incoming polarized proton.

### II.3 Calculation method for spin alignment in longitudinally polarized $`pp`$ collisions

Similar to that in $`e^+e^{}`$ annihilation, to calculate the spin density matrix of the high $`p_T`$ vector mesons in $`pp`$ collisions, we also divide them into two groups (A) and (B) according to whether they contain an outgoing quark in the hard scattering and consider them separately. The difference is that, now the fragmenting quark is an outgoing quark in the hard scattering subprocesses and its polarization is not a constant. It depends on the origin of this quark and is given by Eq.(8). Obviously, it is a function of $`x_a`$, $`Q^2`$ and $`y`$. Hence, $`\rho _{00}`$ for vector mesons of group (A) should be also a function of $`x_a`$, $`Q^2`$ and $`y`$ and it is given by,

$$\rho _{00}^V(x_a,Q^2,y|A)=(1+\alpha P_q^2)/(3\alpha P_q^2).$$
(9)
To obtain $`\rho _{00}`$ for vector mesons at a given pseudo-rapidity $`\eta `$, we need to replace the products in Eq.(3) with the corresponding convolutions, i.e.,

$$\rho ^V(\eta )=\frac{_f[D_{V/f}^A(z_c)\rho ^V(A,f)+D_{V/f}^B(z_c)\rho ^V(B)]\sigma _{ppq_fX}}{_fD_{V/f}(z_c)\sigma _{ppq_fX}}$$
(10)
To write it more precisely, we have,

$$\rho _{00}^V(\eta )=\frac{d^2p_T_{abcd}dx_adx_bf_a(x_a,\mu ^2)f_b(x_b,\mu ^2)\frac{d\widehat{\sigma }}{d\eta }\left[D_{V/c}^A(z_c,\mu ^2)\rho _{00}^V(A,c)+D_{V/c}^B(z_c,\mu ^2)\rho _{00}^V\left(B\right)\right]}{d^2p_T_{abcd}𝑑x_a𝑑x_bf_a(x_a,\mu ^2)f_b(x_b,\mu ^2)\frac{d\widehat{\sigma }}{d\eta }\left(abcd\right)D_{V/c}(z_c,\mu ^2)},$$
(11)
where $`p_T`$ is the transverse momentum of the vector meson; $`f_a(x_a,\mu ^2)`$ and $`f_b(x_b,\mu ^2)`$ are the unpolarized distribution functions of partons $`a`$ and $`b`$ in proton at the scale $`\mu `$, $`x_a`$ and $`x_b`$ are the corresponding momentum fractions carried by $`a`$ and $`b`$; $`D_{V/c}(z_c,\mu ^2)`$ is the usual fragmentation function of parton $`c`$ into vector meson $`V`$, $`D_{V/c}^A(z_c,\mu ^2)`$ and $`D_{V/c}^B(z_c,\mu ^2)`$ are those of parton $`c`$ into vector mesons of group (A) and (B) respectively, $`z_c`$ is the momentum fraction of parton $`c`$ carried by the produced $`V`$; $`d\widehat{\sigma }/d\eta `$ is the cross section at the parton level. The cross section of the hard subprocess can be calculated using perturbative QCD. The summation in Eq. (11) runs over all possible subprocesses.

The fragmentation functions are independent of the spin properties and can be calculated using an available Monte Carlo event generator. In practice, we generate a $`pp`$ collision event and then search for the produced vector meson $`V`$. We then calculate its contribution to $`\rho _{00}^V(\eta )`$ by tracing back its origin: If the vector meson $`V`$ belongs to group (B), then its contribution to $`\rho _{00}^V(\eta )`$ is given by $`\rho _{00}(B)=1/3`$. If it belongs to group (A), we then need to further trace back the origin of the fragmenting quark to calculate its polarization $`P_q`$ using Eq.(8). We insert the $`P_q`$ into Eq.(9) to get the contribution $`\rho _{00}(A)`$ of such vector meson to $`\rho _{00}^V(\eta )`$. After running the program for sufficiently large number of events, we obtain the final $`\rho _{00}^V(\eta )`$. In terms of a mathematic formula, the calculation in this procedure is expressed as,

$$\rho _{00}^V(\eta )=[\underset{i=1}{\overset{N(\eta |A)}{}}\rho _{00}^i(x_a,Q^2,y|A)+\frac{1}{3}N(\eta |B)]/N(\eta ),$$
(12)
where $`N(\eta |A)`$ and $`N(\eta |B)`$ are the numbers of vector mesons of group (A) and (B) as a function of $`\eta `$, and $`N(\eta )=N(\eta |A)+N(\eta |B)`$ is the total number of vector mesons. The results for $`\rho _{00}^V(\eta )`$ in longitudinally polarized $`pp`$ collisions are presented in next section.

## III Results of longitudinally polarized case

Using the method described in last section, we calculate the $`\rho _{00}`$’s for different vector mesons in $`pp`$ collisions with one proton beam longitudinally polarized. We now summarize the obtained results in the following. Using the generator $`\mathrm{pythia}`$, we first calculate the different contributions to vector meson production, i.e., those of group (A) and (B). As an example, we show the results for $`K^+`$ at $`\sqrt{s}`$=500 GeV and $`p_T`$$`>`$13 GeV in Fig. 2. From these results, we can see clearly that, the decay contribution is indeed very small. The contribution of those containing an outgoing quark in the hard scattering subprocess is very high, even larger than $`90\%`$ for large $`\eta `$ region. The results of $`\rho _{00}`$’s for different vector mesons are shown in Fig. 3. We see that, the spin alignment for $`K^+`$, $`\rho ^\pm `$ and $`\rho ^0`$ are significantly high. The obtained $`\rho _{00}`$’s for these mesons increase from 1/3 to about 0.45 with increasing $`\eta `$. There is nearly no spin alignment in the negative $`\eta `$ region since the scattered quark moving in the negative $`z`$ direction is almost unpolarized. The spin alignment for $`K^0`$ is smaller in both positive and negative $`\eta `$ region, this is because $`d`$ quark fragmentation dominate $`K^0`$ production at high $`p_T`$ and the polarization of $`d`$ quark $`|\mathrm{\Delta }d(x,Q^2)/d(x,Q^2)|`$ is smaller than that of $`u`$ quark $`|\mathrm{\Delta }u(x,Q^2)/u(x,Q^2)|`$.

We can also make calculations for other vector mesons, such as $`K^{}`$, $`\overline{K}^0`$, and $`\varphi `$. However, since these mesons do not contain $`u`$ or $`d`$ valence quark, which dominate the high $`p_T`$ hadron production in $`pp`$ collisions, their production rates at high $`p_T`$ are small. This makes the statistics much worse than those for $`K^+`$, $`K^0`$, $`\rho ^\pm `$ and $`\rho ^0`$. For the same reason, the $`s`$ quark distribution is much more important. Hence, the spin alignment for these mesons are much more sensitive to the helicity distribution of the strange quark in nucleon, which is less precisely determined yet. It is thus not a good choice to study the spin transfer in fragmentation process by measuring the spin alignment of these mesons. But, it is possible to use them to study the helicity distribution of the strange sea in nucleon if they can be measured with good accuracy.

## IV Transversely polarized case

When one of the proton beams is transversely polarized, the incoming quark $`a`$ for the elementary hard scattering process can also be transversely polarized in the same direction. Its polarization $`P_{aT}`$ is determined by the transversity distribution function $`\delta q(x_a,Q^2)`$ and the unpolarized distribution function $`q(x_a,Q^2)`$, i.e., $`P_{aT}=\delta q(x_a,Q^2)/q(x_a,Q^2)`$. There is no gluon transversity distribution at leading twist, thus the incoming gluon in the hard scattering is not transversely polarized, which is different from the longitudinally polarized case. The transverse polarization can also be transferred from the incoming quark to the outgoing quark in the hard scattering. The polarization direction of the outgoing quark is also transverse to the moving direction of the quark but is in general different from that of the incoming quark. Both the magnitude and the direction of the polarization of the outgoing quark can be calculated using pQCD for the hard elementary process. The results are summarized in the following.

We recall that, for a quark $`a`$ with transverse polarization $`P_{aT}`$, the spin density matrix $`\rho _a`$ in the helicity basis is given by,

$$\rho _a^{(in)}=\frac{1}{2}\left(\begin{array}{cc}1& P_{aT}e^{i\varphi }\\ P_{aT}e^{i\varphi }& 1\end{array}\right),$$
(13)
where $`\varphi `$ is the angle between the polarization direction and the $`x`$ axis. The $`xy`$-plane is perpendicular to the moving direction of the quark, which is taken as the $`z`$-direction; and the $`x`$-axis is taken as the normal direction of the scattering plane. The spin density matrix of the outgoing quark in the helicity basis can be obtained from $`\rho _a^{(in)}`$ given in Eq. (13) and the helicity amplitudes of the scattering. The result is given by,

$$\rho _q=\frac{1}{2}\left(\begin{array}{cc}1& P_{aT}D_T(y)e^{i\varphi }\\ P_{aT}D_T(y)e^{i\varphi }& 1\end{array}\right),$$
(14)
where $`D_T(y)`$ is a real function of $`y`$ defined in Section 2B for the scattering and is called the polarization transfer factor. $`D_T(y)`$ for different subprocesses have been calculated using pQCD and can be found, for example, in Ref.\[References\].

From Eq. (14), we see that the outgoing quark is also transversely polarized. The magnitude of the the polarization is given by $`D_T(y)P_{aT}`$, i.e,

$$P_{qT}(x_a,Q^2,y)=D_T(y)\delta q(x_a,Q^2)/q(x_a,Q^2).$$
(15)
We see also that, the angle between the polarization direction and the $`x`$-axis remains the same before and after the scattering. More precisely, the polarization vector of the outgoing quark is in the plane transverse to the moving direction of the quark and the angle between the polarization direction and the normal of the scattering plane is the same as that for the incoming quark. In another word, the direction of transverse polarization of the incoming and that of the outgoing quark are related to each other by a rotation around the normal of the scattering plane, which changes the moving direction of the quark from the incoming to the outgoing direction. (C.f. Fig. 2 of Ref.\[References\]). We note that the direction of the polarization of the outgoing quark depends not only on the scattering angle but also on the azimuthal direction. In contrast, the magnitude of the polarization depends only on the scattering angle. The dependence is given by the corresponding spin transfer function $`D_T(y)`$, which is only a function of $`y`$ defined above.

Having seen that the outgoing quarks are also transversely polarized if the incoming protons are transversely polarized, we now discuss the spin alignment of the vector mesons produced in the fragmentation of such quarks. We recall that, the relation given in Eq. (1) is obtained for the fragmentation of the longitudinally polarized quarks by fitting the $`e^+e^{}`$ annihilation data. It shows that the anti-quark which combines the fragmenting quark to form the vector meson is polarized in the opposite direction as that of the fragmenting quark. To test whether the relation is in general true for the fragmentation of a longitudinally polarized quark in any high energy reactions, we extended the calculations of the spin alignment of vector mesons to other lepton-induced reactions in \[References\] and longitudinally polarized $`pp`$ collisions in last sections. Now we further assume that the relation is also true for the fragmentation of a quark polarized in any direction, i.e, we extend the relation to,

$$\stackrel{}{P}_{\overline{q}}=\alpha \stackrel{}{P}_q.$$
(16)
We use this relation to calculate the spin alignment of the vector mesons in transversely polarized $`pp`$ collisions so that this extension can be tested in future experiments. The calculations are carried out in a way similar to that in the longitudinally polarized case. We now present the calculations and the results in the following two different frames.

(1) In the helicity frame: In the helicity frame of the outgoing quark of the hard scattering, its spin density matrix is off-diagonal and is given by Eq. (14). From the relation given by Eq. (16), the spin density matrix of the anti-quark which combines with the fragmenting quark to form the vector meson is,

$$\rho _{\overline{q}}=\frac{1}{2}\left(\begin{array}{cc}1& \alpha P_{qT}e^{i\varphi }\\ \alpha P_{qT}e^{i\varphi }& 1\end{array}\right).$$
(17)
We build the direct product of $`\rho _q`$ given by Eq.(14) with $`\rho _{\overline{q}}`$ given by Eq.(17), transform it to the coupled basis, and obtain the spin density matrix of the vector meson of group (A) as,

$$\rho ^{(h)}(A)=\frac{1}{3\alpha P_{qT}^2}\left(\begin{array}{ccc}1& \frac{1}{\sqrt{2}}P_{qT}\left(1\alpha \right)e^{i\varphi }& \alpha P_{qT}^2e^{2i\varphi }\\ \frac{1}{\sqrt{2}}P_{qT}\left(1\alpha \right)e^{i\varphi }& 1\alpha P_{qT}^2& \frac{1}{\sqrt{2}}P_{qT}\left(1\alpha \right)e^{i\varphi }\\ \alpha P_{qT}^2e^{2i\varphi }& \frac{1}{\sqrt{2}}P_{qT}\left(1\alpha \right)e^{i\varphi }& 1\end{array}\right),$$
(18)
where the superscript “$`(h)`$" denotes the helicity frame. Hence, the $`\rho _{00}^{(h)}(A)`$ for such vector mesons in the helicity basis is given by,

$$\rho _{00}^{(h)}(x_a,Q^2,y|A)=(1\alpha P_{qT}^2)/(3\alpha P_{qT}^2).$$
(19)
We can see that, as can be expected, if we look only at the $`00`$-component of the spin density matrix, the result is independent of the azimuthal direction of the transverse polarization, i.e. independent of the angle $`\varphi `$. It depends only on the magnitude of the transverse polarization $`P_{qT}`$, thus on the variables $`x_a,Q^2`$ and $`y`$. We see also that $`\rho _{00}^{(h)}(A)`$ is smaller than $`1/3`$, which is different from the results in the longitudinally polarized case. This shows that, if the fragmenting quark is transversely polarized, the produced vector meson should have a larger probability to be in the helicity $`\pm 1`$ states. As long as there exists the contribution of the vector mesons of group (A), the final results for $`\rho _{00}^{(h)}(\eta )`$ in this frame should also be smaller than $`1/3`$. We know from Fig. 2 that the contribution of the group (A) increase with increasing $`\eta `$, so we expect a decreasing $`\rho _{00}^{(h)}(\eta )`$ as $`\eta `$ increases for vector mesons such as $`K^+`$ and $`\rho ^\pm `$. It is significantly smaller than $`1/3`$ for large $`\eta `$. Hence, the measurements of $`\rho _{00}^{(h)}(\eta )`$ of the vector mesons in transversely polarized $`pp`$ collisions should give a good check of the extension of the relation of Eq. (1) to the transversely polarized case.

(2) In the transversity frame: To compare with the results in the helicity frame for the longitudinally polarized case, we also calculate the spin alignment of the vector mesons in the transversely polarized case in the “transversity frame”. In this frame, similar to the helicity frame for longitudinally polarized case, we take the (transverse) polarization direction of the fragmenting quark as the quantization axis. The spin density matrix of the outgoing quark and that of the anti-quark are both diagonal in this frame, and they are given by

$$\rho _q^{(t)}=\frac{1}{2}\left(\begin{array}{cc}1+P_{qT}& 0\\ 0& 1P_{qT}\end{array}\right),$$
(20)
$$\rho _{\overline{q}}^{(t)}=\frac{1}{2}\left(\begin{array}{cc}1\alpha P_{qT}& 0\\ 0& 1+\alpha P_{qT}\end{array}\right),$$
(21)
where the superscript “$`(t)`$" denotes the transversity frame. We see that they have completely the same form as those for the corresponding quarks and anti-quarks in longitudinally polarized case in the helicity frame. The corresponding spin density matrix $`\rho ^{(t)}(A)`$, which can be obtained from the direct product $`\rho _q^{(t)}\rho _{\overline{q}}^{(t)}`$, is also diagonal, and the result for the $`00`$-component $`\rho _{00}^{(t)}(A)`$ is given by,

$$\rho _{00}^{(t)}(x_a,Q^2,y|A)=(1+\alpha P_{qT}^2)/(3\alpha P_{qT}^2).$$
(22)
We note that, since $`\rho ^{(t)}(A)`$ is diagonal in the transversity frame, there is a simple relation between $`\rho _{00}^{(h)}(A)`$ and $`\rho _{00}^{(t)}(A)`$. The relation can easily be obtained by using $`\rho _{00}^{(h)}(A)=h=0|\rho (A)|h=0=_{m_T,m_T^{}}h=0|m_Tm_T|\rho (A)|m_T^{}m_T^{}|h=0`$ and $`m_T|\rho (A)|m_T^{}=\rho _{m_T,m_T^{}}^{(t)}=\rho _{m_T,m_T}^{(t)}\delta _{m_T,m_T^{}}`$. (Here, $`h`$ denotes the helicity and $`m_T,m_T^{}`$ denote the spin component along the quantization axis in the transversity frame.) It is given by,

$$2\rho _{00}^{(h)}(x_a,Q^2,y|A)+\rho _{00}^{(t)}(x_a,Q^2,y|A)=1.$$
(23)
We see that, $`\rho _{00}^{(t)}(x_a,Q^2,y|A)`$ can also be obtained from this relation and the result for $`\rho _{00}^{(h)}(x_a,Q^2,y|A)`$ given by Eq. (19). We also see that, the obtained $`\rho _{00}^{(t)}(x_a,Q^2,y|A)`$ given by Eq.(22) is in general larger than $`1/3`$. This means that, the probability for the produced vector mesons to be in the spin states which are transverse to the moving direction are relatively smaller. We see that, although the values of $`\rho _{00}^{(t)}(x_a,Q^2,y|A)`$ and $`\rho _{00}^{(h)}(x_a,Q^2,y|A)`$ are different, they represent the same physical meaning.

We should note that, to measure $`\rho _{00}^{(t)}`$ in experiments, we need to determine the quantization axis in this frame. Since quantization axis is chosen as the polarization direction of the quark before the fragmentation, it can be different for different vector mesons. The polarization direction is determined by Eq.(14). In practice, we first find out the normal direction of the scattering plane and then determine the polarization direction according to the rules described below Eq. (15). (C.f., also Fig.2 of Ref.\[References\]). The moving direction of the fragmenting quark can be replaced by the thrust axis of the quark jet.

Inserting Eqs. (19) and (22) into Eq.(12), we can calculate $`\rho _{00}(\eta )`$ for different vector mesons in the above two frames in $`pp`$ collisions with one beam transversely polarized. For these calculations, we first need to calculate $`P_{qT}`$ using Eq.(15). So far, the transversity distribution $`\delta q(x,Q^2)`$ is not known yet. In order to get a feeling of how large the spin alignment can be, we make an estimation using the simple form for $`\delta q(x)`$ in the light-cone SU(6) quark-spectator model as input. In this model, the $`\delta q(x)`$ is given by,

$`\delta u_v(x)=[u_v(x){\displaystyle \frac{1}{2}}d_v(x)]\stackrel{~}{M}_S(x){\displaystyle \frac{1}{6}}d_v(x)\stackrel{~}{M}_V(x);`$
$$\delta d_v(x)=\frac{1}{3}d_v(x)\stackrel{~}{M}_V(x),$$
(24)
where $`\stackrel{~}{M}_S(x)`$ and $`\stackrel{~}{M}_V(x)`$ are related to the Melsh-Wigner rotation, which can be calculated according to the procedure given in \[References\]. The obtained results of $`\rho _{00}(\eta )`$ for different vector mesons at $`\sqrt{s}=500`$ GeV and $`p_T>13`$ GeV in the helicity frame are shown in Fig. 4 and those in the transversity frame are shown in Fig. 5. For comparison, we also show the results obtained by taking the upper limit of $`\delta q(x,Q^2)`$ in Soffer’s inequality,

$$|\delta q(x,Q^2)|\frac{1}{2}[\mathrm{\Delta }q(x,Q^2)+q(x,Q^2)].$$
(25)
From these figures, we see that, in the helicity frame, the $`\rho _{00}`$’s of vector mesons are smaller than $`1/3`$ in the large $`\eta `$ region, but the effect is not very significant. In the transversity frame, however, the $`\rho _{00}`$’s of $`K^+`$, $`\rho ^\pm `$ and $`\rho ^0`$ are larger than $`1/3`$ and increase to about 0.5 with increasing $`\eta `$, which are similar to the results of longitudinally polarized case. However, the results of $`K^0`$’s $`\rho _{00}`$ in both frames are much near the unpolarized case of 1/3. This is because, the $`d`$ quark fragmentation dominates the $`K^0`$’s production at high $`p_T`$, and according to the Soffer’s inequality, $`|\delta d(x,Q^2)|`$ should be small, since $`\mathrm{\Delta }d(x,Q^2)`$ is negative for valence quark.

As we mentioned above, the gluon is not transversely polarized, so it is much safer to assume that there is no contribution of gluon fragmentation to the spin alignment of vector mesons in transversely polarized $`pp`$ collisions. Therefore, the high $`p_T`$ cut which we invoked in longitudinally polarized case to ensure a small contribution of gluon fragmentation is unnecessary now. Here, only a relative high $`p_T`$ cut is required to ensure the validity of pQCD. We therefore make also the calculations for $`p_T>8`$ GeV at $`\sqrt{s}=500`$ GeV and the results are given in Fig. 6 and Fig. 7. The results are much similar to those obtained in the case of higher $`p_T`$ cut $`p_T>13`$ GeV.

## V Summary

In summary, we calculate the spin alignment of vector mesons with high $`p_T`$ in high energy $`pp`$ collisions with one longitudinally polarized beam by taking the spin of vector meson as the sum of the spin of the scattered quark and that of the secondary anti-quark produced in the fragmentation process. The results for different vector mesons at RHIC energy are presented. These results show that there exist quite significant spin alignment for $`K^+`$, $`\rho ^\pm `$ and $`\rho ^0`$. We also study the spin alignment of vector mesons in the case that one of the proton beams is transversely polarized. Using the same method as that in the longitudinally polarized case, we obtain the results in the helicity and the transversity frame. The results show that, in the large $`\eta `$ region, the $`\rho _{00}`$’s in the helicity frame are smaller than $`1/3`$, while the $`\rho _{00}`$’s in the transversity frame are larger than $`1/3`$. The magnitudes of the effect depend on the transversity distribution of the quark in proton. They are expected to be more significant in the large $`\eta `$ region. Measurements of them can provide useful information for spin effects in strong interaction and the spin structure of nucleon.

We thank Li Shi-yuan, Si Zong-guo, Xie Qu-bing and other members in the theoretical particle physics group of Shandong University for helpful discussions. This work was supported in part by the National Science Foundation of China (NSFC) under the approval number 10175037 and the Education Ministry of China under Huo Ying-dong Foundation.

Figure captions

Fig.1: (a) Contributions of quark and gluon fragmentation to $`\rho ^0`$ production, and (b) the ratio of the quark fragmentation’s contribution to the total rate, in $`pp`$$``$$`\rho ^0X`$ at $`\sqrt{s}`$=500 GeV and $`p_T`$$`>`$13 GeV.

Fig.2: Different contributions to $`K^+`$ production in $`pp`$$``$$`K^+X`$ at $`\sqrt{s}`$=500 GeV and $`p_T`$$`>`$13 GeV as functions of $`\eta `$. Here ‘containing hard u" and similar denote the contributions of those containing the outgoing quark from the hard scattering subprocess.

Fig.3: Spin alignment of different vector mesons with $`p_T`$$`>`$13 GeV in $`pp`$$``$$`VX`$ at $`\sqrt{s}`$=500 GeV in helicity frame when one proton beam is longitudinally polarized. The standard LO set of GRSV2000 and LO set of GRV98 are used for helicity and unpolarized distribution functions respectively. The scale $`\mu `$ is taken as the transverse momentum of the hard scattering subprocess.

Fig.4: Spin alignment of different vector mesons with $`p_T`$$`>`$13 GeV in $`pp`$$``$$`VX`$ at $`\sqrt{s}`$=500 GeV in the helicity frame when one of the proton beams is transversely polarized. The solid lines denote the results obtained using the light-cone quark-spectator model for $`\delta q(x)`$; the dotted lines correspond to the results obtained using the upper limit in Soffer’s inequality.

Fig.5: Spin alignment of different vector mesons with $`p_T`$$`>`$13 GeV in $`pp`$$``$$`VX`$ at $`\sqrt{s}`$=500 GeV in the transversity frame when one of the proton beams is transversely polarized. The solid lines denote the results obtained using the light-cone quark-spectator model for $`\delta q(x)`$; the dotted lines correspond to the results obtained using the upper limit in Soffer’s inequality.

Fig.6: The same as Fig. 4, but for $`p_T`$$`>`$8 GeV.

Fig.7: The same as Fig. 5, but for $`p_T`$$`>`$8 GeV.